function [fdr_thresh, bon_thresh] = fdr_threshold( input_file, input_thresh, ... 
   mask_file, mask_thresh, df, fdr)
%FALSE_DISCOVERY_RATE 
%
% FDR_THRESH = FDR_THRESHOLD( INPUT_FILE, [INPUT_THRESH, 
%   [MASK_FILE, [MASK_THRESH, [DF, [FDR ]]]]] )
%
% 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 in minc format.
%
% 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 MASK_FILE > MASK_THRESH.
% Default is 0.
%
% 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. Default is 0.05.

%############################################################################
% COPYRIGHT:   Copyright 2000 K.J. Worsley and C. Liao, 
%              Department of Mathematics and Statistics,
%              McConnell Brain Imaging Center, 
%              Montreal Neurological Institute,
%              McGill University, Montreal, Quebec, Canada. 
%              worsley@math.mcgill.ca, liao@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=0
end
if nargin < 5
   df = Inf
end
if nargin < 6
   fdr=0.05   
end

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

% Find lower limit for F statistic values:
if df1<=20
   % Upper 0.2 quantiles of sqrt(F_(nu,Inf)) for nu=1:20:
   q=[1.28 1.27 1.24 1.22 1.21 1.19 1.18 1.17 1.17 1.16 1.15 ...
         1.15 1.14 1.14 1.13 1.13 1.13 1.12 1.12 1.12];
   tlim=q(df1);
else
   % Wilson-Hilferty approximation for nu>20: 
   tlim=(1-2/9/df1+0.84162*sqrt(2/9/df1))^(3/2);
   tlim=round(tlim*100)/100;
end
% Values of the F statistic based on squares of t values:
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))];

if isempty(input_thresh)
   monotonic=find(diff(pt)>0);
   input_thresh=interp1(pt(monotonic),t(monotonic),fdr ...
      *(1+is_tstat)).^(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);
   if ~isempty(mask_file)
      dm=fmris_read_image(mask_file,slice,1);
      mask=reshape(dm.data,numpix,1);
      Y=[Y; img(find(img>input_thresh & mask>mask_thresh))];
      num_voxels=num_voxels+sum(mask>mask_thresh);
   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=min(num_voxels*pt,1);
monotonic=find(diff(pval_bon)>0);
bon_thresh=interp1(pval_bon(monotonic),t(monotonic),fdr ...
   *(1+is_tstat)).^(1/(1+is_tstat))

% False discovery rate:
P_val=interp1(t,pt,Y.^(1+is_tstat),'linear',0)/(1+is_tstat);
[P_sort, index]=sort(P_val);
num_vox_above_thr=length(P_sort);
c=1;
% c=sum(1./(1:num_voxels));
r=max(find(P_sort./(1:num_vox_above_thr)'*num_voxels*c <= fdr));
fdr_thresh=Y(index(r))

% end