function summary=stat_summary(input_file, fwhm, df, mask_file, mask_thresh, ...
   input_thresh, flip);

%STAT_SUMMARY
%
% STAT_SUMMARY( INPUT_FILE , [FWHM , [DF , [MASK_FILE, MASK_THRESH ...
%   [, INPUT_THRESH [, FLIP]]]]] )
%
% Finds local maxima and clusters of INPUT_FILE above a INPUT_THRESH
% and MASK_FILE above MASK_THRESH (ignored if empty; empty by default).
% If INPUT_THRESH < 1 then it is taken as a probability and threshold is
% chosen so that the uncorrected P-value is INPUT_THRESH (0.001 by default),
% If MASK_THRESH is a vector [a b], a<=b, then mask is a<MASK_FILE<=b. 
% Clusters are voxels connected in any of the 6 diresctions. A local maximum 
% is a voxel which is greater than or equal to all its 6 neighbours,
% and strictly greater than at least one of them. 
%
% FLIP: INPUT_FILE is multiplied by FLIP before processing. For T statistic
% images, FLIP = -1 will look for negative peaks nd clusters. Default is 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.
% 
% SUMMARY is a matrix with 12 columns. 
% Col 1: index of cluster, in descending order of cluster volume. 
% Col 2: volume of cluster in mm^3.
% Col 3: resels of cluster.
% Col 4: P-value of cluster extent.
% Col 5: P-value if the cluster was chosen in advance, e.g. nearest to an ROI.
% Col 6: values of local maxima, sorted in descending order.
% Col 7: P-value of local maxima.
% Col 8: P-value if the peak was chosen in advance, e.g. nearest to an ROI.
% Col 9: Q-value or false discovery rate ~ probability that voxel is not signal.
% Cols 10-12: x,y,z coords of local maxima in C notation i.e. starting at 0.
%        If x or y step sizes are negative, then their voxel indices 
%        are reversed so that they agree with 'register'.

if nargin < 2
   fwhm=0
end
if nargin < 3
   df=Inf
end
if nargin < 4
   mask_file=[];
end
if nargin < 6
   input_thresh=0.001
end
if nargin < 7
   flip=1
end

if input_thresh < 1
   input_thresh=stat_threshold(0,1,0,df,input_thresh)
end

[search_volume, voxel_volume]= ...
   glass_brain(input_file,input_thresh,mask_file,mask_thresh,flip);

lm=locmax(input_file,input_thresh,mask_file,mask_thresh,fwhm,flip);

[p_peak, p_cluster, p_peak1, p_cluster1]= ...
   stat_threshold(search_volume, voxel_volume, fwhm, ... 
   df, lm(:,1), input_thresh, lm(:,6), 3);

q_value = fdr_threshold( input_file, input_thresh, ... 
   mask_file, mask_thresh, df, lm(:,1),flip);

n=size(lm,1);
['clus    vol  resel  p-val  (one)   peak  p-val  (one)   q-val   x   y   z']
summary=[repmat('  ',n,1) num2str(round(lm(:,5))) repmat('  ',n,1) ...
   num2str(round(lm(:,6))) repmat('  ',n,1) ...
   num2str(round(lm(:,7))) repmat('  ',n,1) ...
   num2str(round(p_cluster*1000)/1000) repmat(' (',n,1)  ... 
   num2str(round(p_cluster1*1000)/1000) repmat(')  ',n,1)  ... 
   num2str(round(lm(:,1)*100)/100) repmat('  ',n,1)  ...  
   num2str(round(p_peak*1000)/1000) repmat(' (',n,1)  ...  
   num2str(round(p_peak1*1000)/1000) repmat(')  ',n,1)  ... 
   num2str(round(q_value'*1000)/1000) repmat('  ',n,1)  ... 
   num2str(round(lm(:,2:4)))]

return