
function F = gaussblur(f,FWHM);
%function F = gaussblur(f,FWHM);
% f : input 2D image 
% FWHM: filtersize in mm.
% WARNING: it works for 2D only so far.
% (c) Moo K. Chung, July 2003. March 2004.
% See Chung et al., 2001. A Unified Statistical Appraoch to 
% Deformation-based morphometry, NeuroImage for implementation detail.
% to do list - implement this as 1D version.


%FWHM=12;

%ksize is the number of neighboring voxels we are averaging. it depends on FWHM.
%ksize is allowd to be an odd interger.

ksize=round(FWHM/2);
if ~mod(ksize,2)
    ksize=ksize+1;
end;
% see Chung et al.(2001) for the relationship between t and FWHM.
% Given FWHM, we need sigma=FWHM/(2*(log4)^0.5)
t= FWHM^2/(16*log(2));

d=round(ksize/2);

% we truncate and normalize the Gaussian kernel for fast computation
% 2D and 3D version of kerenl is slightly different
if ndims(f)==2
    domain_z=kron(ones(ksize,1),[1:ksize])-d;
    domain_y=kron(ones(1,ksize),[1:ksize]')-d;
    ker=exp(-(domain_y.^2 +domain_z.^2)/(4*t))/t;
    sum_ker =sum(sum(ker));
else
    %implment domain_x for 3D version here.
    domain_x=kron(ones(ksize,1),[1:ksize])-d;
    domain_z=kron(ones(ksize,1),[1:ksize])-d;
    domain_y=kron(ones(1,ksize),[1:ksize]')-d;
    ker=exp(-(domain.x.^2 + domain_y.^2 +domain_z.^2)/(4*t))/t^(3/2);
    sum_ker =sum(sum(sum(sum(ker))));
end;
ker=ker/sum_ker;

% smoothing is done by 2D convolution
F=conv2(f,ker,'same');





