function [fdr_thresh, bon_thresh] = fdr_threshold( input_file, input_thresh, ... 
   mask_file, mask_thresh, df, fdr, flip, k)

%FALSE_DISCOVERY_RATE 
%
% [FDR_THRESH, BON_THRESH] = FDR_THRESHOLD( INPUT_FILE, [INPUT_THRESH, 
%   [MASK_FILE, [MASK_THRESH, [DF, [FDR [, FLIP [, K ]]]]]]] )
%
% Calculates the threshold (FDR_THRESH) for a t or F image for 
% controlling the false discovery rate (FDR). The FDR is the expected 
% proportion of false positives among the voxels above FDR_THRESH. 
% This threshold is higher than that for controlling the expected proportion 
% of false positives in the whole search volume, and usually lower than the 
% Bonferroni threshold (printed out as BON_THRESH) which controls the 
% probability that *any* voxels are above the threshold. References: 
% Benjamini, Y. and Hochberg, Y. (1995). Controlling the false discovery 
% rate: a practical and powerful approach to multiple testing. Journal of 
% the Royal Statistical Society, Series B, 57:289-300, and
% Genovese et al. (2001). Thresholding statistical maps in functional
% neuroimaging using the false discovery rate. Neuroimage, in press.
%
% INPUT_FILE is the input t or F statistic image.
%
% INPUT_THRESH is the lower limit for searching the t or F statistic image.
% If empty, it is the threshold corresponding to a P-value of FDR at a single 
% point. Default is []. 
%
% MASK_FILE is a mask file. If empty, it is ignored. Default is [].
%
% MASK_THRESH defines the search volume as the first frame of MASK_FILE 
% > MASK_THRESH. If MASK_THRESH is a vector [a b], a<=b, then mask 
% is a < MASK_FILE <= b. If empty (default), calls fmri_mask_threshold. 
%
% DF: If a scalar, then it is the df of the t statistic image;
% if DF >=1000 then DF is set it Inf, so that it calculates 
% thresholds for a Gaussian image (if DF is very large the t-dbn 
% is almost identical to the Gaussian dbn).
% If DF=[DF1, DF2] then these are the df's of the F statistic image.
% If DF2 >= 1000 then DF2 is set to Inf. Default is Inf.
%
% FDR is the desired false discovery rate. 
% If the first element is greater than 0.1, then they are treated as 
% thresholds and Q-values are returned. Default is 0.05.
%
% FLIP: INPUT_FILE is multiplied by FLIP before processing. For T statistic
% images, FLIP = -1 will look for negative peaks and clusters. Default is 1.
%
% K is the number of conjunctions. If K > 1, calculates P-values for peaks 
% (not clusters) of the minimum of K independent t or F maps. Default is 1.

%############################################################################
% COPYRIGHT:   Copyright 2002 K.J. Worsley, 
%              Department of Mathematics and Statistics,
%              McConnell Brain Imaging Center, 
%              Montreal Neurological Institute,
%              McGill University, Montreal, Quebec, Canada. 
%              worsley@math.mcgill.ca
%
%              Permission to use, copy, modify, and distribute this
%              software and its documentation for any purpose and without
%              fee is hereby granted, provided that this copyright
%              notice appears in all copies. The author and McGill University
%              make no representations about the suitability of this
%              software for any purpose.  It is provided "as is" without
%              express or implied warranty.
%############################################################################

% Defaults:

if nargin < 2
   input_thresh=[]
end
if nargin < 3
   mask_file=[]
end
if nargin < 4
   mask_thresh=[]
end
if nargin < 5
   df = Inf
end
if nargin < 6
   fdr=0.05   
end
if nargin < 7
   flip=1   
end
if nargin < 8
   k = 1
end

if ~isempty(mask_file)
   if isempty(mask_thresh)
      mask_thresh=fmri_mask_thresh(mask_file);
   end
   mask_thresh1=mask_thresh(1);
    if length(mask_thresh)>=2
        mask_thresh2=mask_thresh(2);
    else
        mask_thresh2=Inf;
    end
end

if length(df)==1
    df=[df 0];
end
if size(df,1)==1
    df=[df; Inf Inf]
end
if size(df,2)==1
    df=[df [0; df(2,1)]]
end

% is_tstat=1 if it is a t statistic, 0 if it is an F statistic:
is_tstat=(df(1,2)==0);
if is_tstat
   df1=1;
   df2=df(1,1);
else
   df1=df(1,1);
   df2=df(1,2);
end
if df2 >= 1000
   df2=Inf;
end

% Values of the F statistic based on squares of t values:
tlim=0.1;
t=[100:-0.1:10.1 10:-0.01:tlim].^2;

% Find the upper tail probs of the F distribution by 
% cumulating the F density using the mid-point rule:
n=size(t,2);
tt=(t(1:(n-1))+t(2:n))/2;
if df2==Inf
   u=df1*tt;
   b=exp(-u/2-gammaln(df1/2)-df1/2*log(2)+(df1/2-1)*log(u));
else  
   u=df1*tt/df2;
   b=(1+u).^(-(df1+df2)/2).*u.^(df1/2-1)*df1/df2/beta(df1/2,df2/2);
end
pt=[0 -cumsum(b.*diff(t))].^k;

if isempty(input_thresh)
   monotonic=find(diff(pt)>0);
   input_thresh=interp1(pt(monotonic),t(monotonic),fdr ...
      *(1+is_tstat)^k).^(1/(1+is_tstat));
end
input_thresh=max([input_thresh tlim.^(2-is_tstat)])

d=fmris_read_image(input_file,0,0);
d.dim
%numframes=d.dim(4);
numslices=d.dim(3);
numys=d.dim(2);
numxs=d.dim(1);
numpix=numxs*numys;
Steps=d.vox
voxel_volume=abs(prod(Steps))

Y=[];
num_voxels=0;
for slice=1:numslices
   d=fmris_read_image(input_file,slice,1);
   img=reshape(d.data,numpix,1)*flip;
   if ~isempty(mask_file)
      dm=fmris_read_image(mask_file,slice,1);
      mask=reshape(dm.data,numpix,1);
      mask= (mask>mask_thresh1 & mask<=mask_thresh2);
      Y=[Y; img(find(img>input_thresh & mask))];
      num_voxels=num_voxels+sum(mask);
   else
      Y=[Y; img(find(img>input_thresh))];
      num_voxels=num_voxels+length(img);
   end
end

% Bonferroni: 
search_volume=num_voxels*voxel_volume
pval_bon=num_voxels*pt;
monotonic=find(diff(pval_bon)>0);

% False discovery rate
P_val=interp1(t,pt,Y.^(1+is_tstat),'linear',0)/(1+is_tstat)^k;
[P_sort, index]=sort(P_val);

nfdr=length(fdr);
fdr_thresh=zeros(1,nfdr);
if fdr(1)<1
   bon_thresh=interp1(pval_bon(monotonic),t(monotonic),fdr ...
      *(1+is_tstat).^k).^(1/(1+is_tstat));
   % False discovery rate threshold:
   q=P_sort./(1:length(P_sort))'*num_voxels;
   for i=1:nfdr
      r=max(find(q <= fdr(i)));
      fdr_thresh(i)=Y(index(r));
   end
   if nfdr<=1
      bon_thresh
      fdr_thresh
   end
else
   P0=interp1(t,pt,fdr.^(1+is_tstat),'linear',0)/(1+is_tstat)^k;
   bon_thresh=P0*num_voxels;
   % Q-value:
   for i=1:nfdr
      K=find(P_sort>=P0(i));
      fdr_thresh(i)=min(P_sort(K)./K*num_voxels);
   end
   if nfdr<=1
      bon_thresh
      q_value=fdr_thresh
   end
end

return