% SURF_TSTAT_THRESHOLD Returns the threshold for tstat-maps after mulpliple comparisions
% THRESHOLD = SURF_TSTAT_THRESHOLD(AREA, FWHM, DF, SIGNIFICANCE)
%
% AREA is the total surface area of the cortex
% FWHM is the fwhm of a smoothing kernel applied to the data. See [1] and [3] for surface smoothing.
% DF is the degrees of freedom of the tstat image.
% SIGNIFICANCE is the desired significance in the range from 0 to 0.2.
%
%
% (C) Moo K. Chung 
% Last update May 8, 2005.
% mchung@stat.wisc.edu
% http://www.stat.wisc.edu/softwares/hk/hk.html
% This is modified from Keith K. Worsley's original FSTAT code
% The details on the mathematical basis of of the algorithm can be found in [1] and [2].
% Surface smoothing and FWHM can be found in [1] and [3].
% [1] Chung, M.K., Robbins,S., Dalton, K.M., Davidson, R.J., Evans, A.C. 2005 
%     Cortical thickness analysis in autism via heat kernel smoothing. NeuroImage, in press 
%     http://www.stat.wisc.edu/~mchung/papers/ni_heatkernel.pdf
% [2] Worsley et al. 1996, Human Brain Mapping, 4:58-73
% [3] Chung, M.K., Robbins,S., Dalton, K.M., Davidson, Alexander, A.L., R.J., Evans, A.C. 2005. 
%      Cortical thickness analysis in autism via heat kernel smoothing. NeuroImage 25:1256-1265.
%      http://www.stat.wisc.edu/~mchung/papers/ni_heatkernel.pdf
%
% 
% Note: There is a mistake in NeuroImage 2003 example [1] due to typos in formula
% area = 275800; fwhm=20; df = 28-1, significance = 0.025
% resulting threshold = 5.1 should be corrected to 5.38
% The correct t-stat threshold formula is given in [2] and [3].

function threshold = surf_tstat_threshold(area, fwhm, df, significance)


t=[100:-0.1:10.1 10:-0.01:1.28];
pt=1-tcdf(t,df);

% Different way to compute the cumulative distribution of t.
% pt=0.5*betainc(df./(df+t.^2),df/2,0.5);
% Cumulating the t density using the mid-point rule is much faster:
%n=size(t,2);
%tt=(t(1:(n-1))+t(2:n))/2;
%b=(1+tt.^2/df).^(-(df+1)/2)/(sqrt(df)*beta(df/2,1/2));
%pt=[0 -cumsum(b.*diff(t))];

% Random field theory

resels_0=2;
resels_2 = area/(2*fwhm^2);
rho_0=pt;
rho_2= 4*log(2)/(2*pi)^(3/2)*gamma((df+1)/2)/gamma(df/2)/(df/2)^(1/2)*(1+t.^2/df).^(-(df-1)/2).*t;
pval_rf=resels_0*rho_0+resels_2*rho_2;

index=find(abs(pval_rf-significance)==min(abs(pval_rf-significance)));
threshold=t(index)





