function gsmooth = stat_smooth(g, FWHM);
%function gsmooth = stat_smooth(g, FWHM);
%
% Peforms Gaussian kernel smoothing on collection of 2D images. It modifies
% gaussblur2D.m function.
%
% g:            Collection of input images in 3D matrix (image dimension x image
%                dimension x # of images)  
% FWHM :  filtersize in full width at half maximum (FWHM) mm.
%
%gsmooth: Collection of output images in 3D matrix  
%
% (c) 2012  Moo K. Chung     mkchung@wisc.edu
%      Waisman Laboratory for Brain Imaging and Behavior
%      University of Wisconsin-Madison
%    
%
%Reference:
% [1] Chung, M.K., Dalton, K.M., Alexander, A.L., Davidson, R.J. 2004.
%      Less white matter concentration in autism: 2D Voxel-Based 
%      Morphometry. NeuroImage 23:242-251
%      http://www.stat.wisc.edu/~mchung/papers/ni_autism.pdf
%
% [2] Chung, M.K., Worsley, K.J. , Paus, T., Cherif, C.,  Giedd, J.N.,
%      Rapoport, J.L,  Evans, A.C. 2001. A unified statistical approach to 
%      deformation-based morphometry, NeuroImage 14:595-606
%      http://www.stat.wisc.edu/~mchung/papers/ni_deformation.pdf


ksize=round(FWHM);
if ~mod(ksize,2) %ksize ihas to be an odd interger
    ksize=ksize+1;
end;

% The relationship between FWHM and kernel size sig is given in [2]
% sigma=FWHM/(2*(log4)^0.5)
sig= FWHM^2/(16*log(2));
d=round(ksize/2);

% constructing grid
domain_z=kron(ones(ksize,1),[1:ksize])-d;
domain_y=kron(ones(1,ksize),[1:ksize]')-d;

% constructing Gaussian kernel
ker=exp(-(domain_y.^2 +domain_z.^2)/(4*sig))/sig;
sum_ker =sum(sum(ker));

%normalizing kernel wights. sum(sum(ker) =1
ker=ker/sum_ker; 

% smoothing is done by 2D convolution
n=size(g);
gsmooth=zeros(n(1),n(2),n(3));
for i=1:n(3)
    gsmooth(:,:,i)=conv2(g(:,:,i),ker,'same');
end;
