function [peak_threshold, extent_threshold] = ...
   stat_threshold(search_volume, voxel_volume, fwhm,  ... 
   df, p_val_peak, cluster_threshold, p_val_extent, D)
%STAT_THRESHOLD 
%
% [PEAK_THRESHOLD, EXTENT_THRESHOLD, EXTENT_THRESHOLD_1] = STAT_THRESHOLD(  
%   [ SEARCH_VOLUME [, VOXEL_VOLUME [, FWHM  [, DF [, P_VAL_PEAK 
%   [, CLUSTER_THRESHOLD [, P_VAL_EXTENT [, D ]]]]]]]] )
%
% If P-values are supplied, returns the thresholds for local maxima 
% or peaks (PEAK_THRESHOLD) and spatial extent of contiguous voxels 
% above CLUSTER_THRESHOLD (EXTENT_THRESHOLD) for t or F statistic  
% maps from fmrilm or multistat. If thresholds are supplied then it
% returns P-values. For peaks (local maxima), two methods are used:  
% random field theory (Worsley et al. 1996, Human Brain Mapping, 4:58-73),
% and a simple Bonferroni correction based on the number of voxels
% in the search region. The final threshold is the minimum of the two.
% EXTENT_THRESHOLD_1 is the extent of a single cluster chosen in advance
% of looking at the data, e.g. the nearest cluster to a pre-chosen voxel
% or ROI - see Friston KJ. Testing for anatomically specified regional 
% effects. Human Brain Mapping. 5(2):133-6, 1997.
%
% The random field theory threshold is based on the assumption that
% the search region is a sphere, which is a very tight lower bound for any
% non-spherical region, and that the data is isotropic, i.e. equal
% fwhm in each direction. If the data is not isotropic, then replacing
% fwhm by the geometric mean of the three x,y,z fwhm's is a very good
% approximation. 
%
% For clusters, the method of Cao, 1999, Advances in Applied Probability,
% 31:577-593 is used. The cluster size is only accurate for large 
% CLUSTER_THRESHOLD (say >3) and large SEARCH_VOLUME. 
%
% SEARCH_VOLUME is the volume of the search region in mm^3. Default is
% 1000000, i.e. 1000 cc. The method for finding PEAK_THRESHOLD works well 
% for any value of SEARCH_VOLUME, even zero, which simply gives the  
% threshold for the image at a single voxel. 
%
% VOXEL_VOLUME is the volume of one voxel in mm^3. Default is 1*1*1=1.
%
% FWHM is the fwhm of a smoothing kernel applied to the data. Default is
% 0.0, i.e. no smoothing, which is roughly the case for raw fMRI data.
% For motion corrected fMRI data, use at least 6mm;
% for PET data, this would be the scanner fwhm of about 6mm.
%
% 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.
%
% P_VAL_PEAK is the row vector of desired P-values for peaks in the range 
% from 0 to 0.2. If the first element is greater than 1, then they are
% treated as peak values and P-values are returned. Default is 0.05.
%
% CLUSTER_THRESHOLD is the scalar threshold of the image for clusters.
% If it is less than 0.2, it is taken as a probability p, and the 
% cluter_thresh is chosen so that P( T > cluster_thresh ) = p. 
% Default is 0.001, i.e. P( T > cluster_thresh ) = 0.001. 
%
% P_VAL_EXTENT is the row vector of desired P-values for spatial extents
% of clusters of contiguous voxels above CLUSTER_THRESHOLD. 
% If the first element is greater than 1, then they are treated as 
% extent values and P-values are returned. Default is 0.05.
%
% D is the number of dimensions; default is 3.
%
% Examples:
%
% T statistic, 1000000mm^3 search volume, 1mm voxel volume, 20mm smoothing,
% 30 degrees of freedom, P=0.05 and 0.01 for peaks, cluster threshold chosen
% so that P=0.001 at a single voxel, P=0.05 and 0.01 for extent:
%
% stat_threshold(1000000,1,20,30,[0.05 0.01],0.001,[0.05 0.01]);
%
% peak_threshold = 5.0589    5.7803
% Cluster_threshold = 3.3852
% extent_threshold = 3340.2    6315.5
% extent_threshold_1 = 2872.9    5709.5
%
% Check: Suppose we are given peak values of 5.0589, 5.7803 and
% spatial extents of 3841.9, 6944.4 above a threshold of 3.3852,
% then we should get back our original P-values:
% 
% stat_threshold(1000000,1,20,30,[5.0589 5.7803],3.3852,[3340.2 6315.5]);
%
% P_val_peak = 0.0500    0.0100
% P_val_extent = 0.0500    0.0100
% P_val_extent_1 = 0.0372    0.0074
%
% Another way of doing this is to use the fact that T^2 has an F 
% distribution with 1 and 30 degrees of freedom, and double the P values: 
% 
% stat_threshold(1000000,1,20,[1 30],[0.1 0.02],0.002,[0.1 0.02]);
%
% peak_threshold = 25.5927   33.4115
% Cluster_threshold = 11.4596
% extent_threshold = 3297.8    6305.1
% extent_threshold_1 = 1869.0    4402.3
%
% Note that sqrt([25.5927 33.4115 11.4596])=[5.0589 5.7803 3.3852]
% in agreement with the T statistic thresholds, but the spatial extent
% thresholds are very close but not exactly the same.

%############################################################################
% 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 < 1
   search_volume=1000000
end
if nargin < 2
   voxel_volume=1.0
end
if nargin < 3
   fwhm=0.0
end
if nargin < 4
   df = Inf
end
if nargin < 5
   p_val_peak = 0.05
end
if nargin < 6
   cluster_threshold = 0.001
end
if nargin < 7
   p_val_extent = 0.05
end
if nargin < 8
   D=3   
end

if D>3
   fprintf('D>3 not implemented yet');
   thresh=NaN
   return
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))];

% Random field theory
if fwhm>0 & D>0 & search_volume>0
   resels=zeros(1,3);
   if D==2
      radius=(search_volume/pi)^(1/2);
      resels(1)=pi*radius/fwhm;
   end
   if D==3
      radius=(search_volume/(4/3*pi))^(1/3);
      resels(1)=4*radius/fwhm;
      resels(2)=2*pi*(radius/fwhm)^2;
   end 
   resels(D)=search_volume/fwhm^D;
   f1=sqrt(2*log(2)/pi);
   rho=zeros(3,length(t));
   if df2==Inf
      u=df1*t;
      b1=exp(-(df1-2)/2*log(2)-gammaln(df1/2));
      c1=exp((df1-1)/2*log(u)-u/2);
      rho(1,:)=f1*b1*c1;
      
      c2=exp((df1-2)/2*log(u)-u/2);
      rho(2,:)=f1^2*b1*c2.*(u-(df1-1));
      
      c3=exp((df1-3)/2*log(u)-u/2);
      rho(3,:)=f1^3*b1*c3.*(u.^2-(2*df1-1)*u+(df1-1)*(df1-2));
   else   
      u=df1*t/df2;
      b1=exp(gammaln((df1+df2-1)/2)-gammaln(df1/2)-gammaln(df2/2));
      c1=exp((df1-1)/2*log(u)-(df1+df2-2)/2*log(1+u));
      rho(1,:)=f1*b1*c1*sqrt(2);
      
      b2=exp(gammaln((df1+df2-2)/2)-gammaln(df1/2)-gammaln(df2/2));
      c2=exp((df1-2)/2*log(u)-(df1+df2-2)/2*log(1+u));
      rho(2,:)=f1^2*b2*c2.*((df2-1)*u-(df1-1));
      
      b3=exp(gammaln((df1+df2-3)/2)-gammaln(df1/2)-gammaln(df2/2));
      c3=exp((df1-3)/2*log(u)-(df1+df2-2)/2*log(1+u));
      p2=(df2-1)*(df2-2);
      p1=-(2*df1*df2-df1-df2-1);
      p0=(df1-1)*(df1-2);
      rho(3,:)=f1^3*b3*c3/sqrt(2).*(p2*u.^2+p1*u+p0);
   end
   pval_rf=1-exp(-(pt+resels*rho));
else
   if D==0 | search_volume==0
      pval_rf=pt;
   else
      pval_rf=1.0;
   end
end

% Bonferroni 
pval_bon=min((search_volume/voxel_volume+1)*pt,1);

% Minimum of the two:
pval=min(pval_rf,pval_bon);

if p_val_peak(1) < 0.2 
   % Make sure it is monotonic:
   monotonic=find(diff(pval)>0);
   peak_threshold=interp1(pval(monotonic),t(monotonic),p_val_peak ...
      *(1+is_tstat)).^(1/(1+is_tstat))
end
if p_val_peak(1) > 1
   % p_val_peak is treated as a peak value:
   P_val_peak=interp1(t,pval,p_val_peak.^(1+is_tstat))/(1+is_tstat)
   peak_threshold=P_val_peak;
end

% Cluster size:

if fwhm<=0 | D==0 | search_volume==0
   return
end

if cluster_threshold >= 1
   tt=cluster_threshold.^(1+is_tstat);
else
   % cluster_threshold is treated as a probability:
   % Make sure pt is monotonic:
   monotonic=find(diff(pt)>0);
   tt=interp1(pt(monotonic),t(monotonic),cluster_threshold*(1+is_tstat));
   Cluster_threshold=tt.^(1/(1+is_tstat))
end

rhoD=interp1(t,rho(D,:),tt);
p=interp1(t,pt,tt);
% Expected number of clusters:
EL=resels(D)*rhoD/(1+is_tstat);

if df2==Inf
   if p_val_extent(1) < 1 
      pS=-log(1-p_val_extent)/EL;
      extent_threshold=p/rhoD/gamma(D/2+1)*fwhm^D*(-log(pS))^(D/2)
   else
      % p_val_extent is now treated as a spatial extent:
      pS=exp(-(p_val_extent/p*rhoD*gamma(D/2+1)/fwhm^D).^(2/D));
      P_val_extent=1-exp(-pS*EL)
      extent_threshold=P_val_extent;
   end
else
   % Find dbn of S by taking logs then using fft for convolution:
   ny=2^12;
   f=zeros(ny,5);
   a=D/2;
   b2=a*5*sqrt(2/(df1+df2));
   b1=a*log(1-(1-0.000001)^(2/(df2-D)));
   dy=(b2-b1)/ny;
   b1=round(b1/dy)*dy;
   y=((1:ny)'-1)*dy+b1;
   
   yy=exp(y.*(y<0)/a);   
   % Density of log(Beta(1,(df2-D)/2)^(D/2)):
   f(:,1)=(1-yy).^((df2-D)/2-1).*((df2-D)/2).*yy/a.*(y<0);
   
   for i=2:(D+2) 
      nu=(df1+df2-D)*(i==2)+(df2+4-i)*(i>2);
      aa=a*(i==2)-1/2*(i>2);
      yy=y/aa+log(df1+df2);
      % Density of log((chi^2_nu)^aa):
      f(:,i)=exp(nu/2*yy-exp(yy)/2-(nu/2)*log(2)-gammaln(nu/2))/abs(aa);
   end
   % Check: sum(f*dy,1) should be 1
   
   omega=2*pi*((1:ny)'-1)/ny/dy;
   shift=complex(cos(-b1*omega),sin(-b1*omega))*dy;
   prodfft=prod(fft(f),2).*shift.^(D+1);
   % Density of Y=log(B^(D/2)*U^(D/2)/sqrt(det(Q))):
   ff=real(ifft(prodfft));
   % Check: sum(ff*dy) should be 1
   
   mu0=exp(D*log(2)+gammaln(df2-D+1)-gammaln(df2+1) ...
      +gammaln(D/2+1)+gammaln((df1+df2)/2)-gammaln((df1+df2-D)/2));
   % Check: E(exp(Y))=sum(ff.*exp(y)*dy) should equal mu0   
   alpha=p/rhoD/mu0*fwhm^D;
   
   % Integrate the density to get the p-value for one cluster: 
   pS=cumsum(ff(ny:-1:1))*dy;
   pS=pS(ny:-1:1);
   % The number of clusters is Poisson with mean EL:
   pSmax=1-exp(-pS*EL);
   
   if p_val_extent(1) < 1 
      % Make sure pSmax is monotonic:
      monotonic=find(diff(pSmax)<0.000001*min(diff(pSmax)));
      yval=interp1(pSmax(monotonic),y(monotonic),p_val_extent);
      % Spatial extent is alpha*exp(Y) -dy/2 correction for mid-point rule:
      extent_threshold=alpha*exp(yval-dy/2)
   else
      % p_val_extent is now treated as a spatial extent:
      P_val_extent=interp1(y,pSmax,log(p_val_extent/alpha)+dy/2)
      extent_threshold=P_val_extent;
   end
   
   % For a single cluster:
   pSmax=1-exp(-pS);
   if p_val_extent(1) < 1 
      % Make sure pSmax is monotonic:
      monotonic=find(diff(pSmax)<0.000001*min(diff(pSmax)));
      yval=interp1(pSmax(monotonic),y(monotonic),p_val_extent);
      % Spatial extent is alpha*exp(Y) -dy/2 correction for mid-point rule:
      extent_threshold_1=alpha*exp(yval-dy/2)
   else
      % p_val_extent is now treated as a spatial extent:
      P_val_extent_1=interp1(y,pSmax,log(p_val_extent/alpha)+dy/2)
      extent_threshold_1=P_val_extent_1;
   end
end
   
% End