function [t,ts] = multFDR(P,df,q)
% Find common T threshold FDR over several images
% FORMAT [t,ts] = multFDR(P,df,q)
%
% P     Array of filenames of T images
% df    Common degrees-of-freedom of T images
% q     Level at which to control FDR
%
% t     Common level-q FDR T threshold
% ts    level-q FDR T threshold for each image (for comparison)
%______________________________________________________________________________
% $Id: multFDR.m,v 1.2 2004/10/15 02:56:20 nichols Exp $

if nargin<1, P = spm_get(Inf,'IMAGE','Select T images'); end
if nargin<2, df = spm_input('Enter df of T images','+1','r'); end
if nargin<3, FDRa = 0.05; end

V = spm_vol(P);

Pval = [];
Ths = [];
for i = 1:size(P,1)

    tmp = spm_read_vols(V(i));
    tmp = tmp(:); tmp(isnan(tmp)) = []; tmp(tmp==0) = [];
    tmp = 1-spm_Tcdf(tmp,df);

    Ths = [Ths; myFDR(tmp,FDRa)];
    Pval = [Pval; tmp];

end

Ths0 = myFDR(Pval,FDRa);

t = spm_invTcdf(1-Ths0,df);
ts = spm_invTcdf(1-Ths,df);

return



function [pID,pN] = myFDR(p,q)
%
% p   - vector of p-values
% q   - False Discovery Rate level
%
% pID - p-value threshold based on independence or positive dependence
% pN  - Nonparametric p-value threshold
%______________________________________________________________________________
% Based on FDR.m     1.4 Tom Nichols 02/07/02


p = sort(p(:));
V = length(p);
I = (1:V)';

cVID = 1;
cVN = sum(1./(1:V));

pID = p(max(find(p<=I/V*q/cVID)));
pN = p(max(find(p<=I/V*q/cVN)));

return