function fwhm_file=fwhm(input_file_wresid, df, df_resid, output_file_base, df_limit)

% FWHM
%
% Estimates the effective FWHM in mm of the residuals from a linear model,
% as if the residuals were white noise smoothed with a Gaussian filter 
% whose fwhm was FWHM. Applies a first-order correction for small degrees 
% of freedom and small FWHM relative to the voxel size. Great care has been
% taken to make sure that FWHM is unbiased, particularly for low df, so that
% if it is smoothed spatially then it remains unbiased. The bias in FWHM
% is less than 1/10 of the voxel size for FWHM > 2*voxel size. However the 
% continuity correction for small FWHM breaks down if FWHM < voxel size,
% in which case FWHM is generally too large. If FWHM > 50, FWHM = 50.
%
%   FWHM_FILE = FWHM( INPUT_FILE_WRESID, DF1 [, DF2  ...
%                    [, OUTPUT_FILE_BASE [, DF_LIMIT]]] )
%
% INPUT_FILE_WRESID is a .mnc or .img file of (whitened) standardised
% residuals from FMRILM or MULTISTAT, obtained by setting WHICH_STATS(7)=1.
%
% DF is the df of the estimated standard deviation, equal to the first output
% parameter from FMRILM or MULTISTAT. Can be Inf.
%
% DF_RESID is the residual df, equal to the number of observations (scans/
% runs/subjects) - rank of linear model. For residuals from FMRILM, this is 
% the first output parameter DF; for residuals from MULTISTAT, this is
% the second output parameter DF_RESID. Default is DF.
%
% OUTPUT_FILE_BASE: Default is INPUT_FILE_WRESID minus its extension.
% 
% DF_LIMIT controls which method is used. If DF_RESID > DF_LIMIT, then
% the FWHM is calculated assuming the Gaussian filter is arbitrary.
% However if DF_RESID is small, this gives inaccurate results, so
% if DF_RESID <= DF_LIMIT, the FWHM is calculated assuming that the axes of
% the Gaussian filter are aligned with the x, y and z axes of the data. 
% Technically, the methods are right only if DF=Inf or DF=DF_RESID,
% and approximate inbetween. Default is 4. 
%
% FWHM_FILE = OUTPUT_FILE_BASE_FWHM.mnc or .img is the output FWHM file.

%############################################################################
% COPYRIGHT:   Copyright 2001 K.J. Worsley, 
%              Department of Mathematics and Statistics,
%              McConnell Brain Imaging Center, 
%              Montreal Neurological Institute,
%              McGill University, Montreal, Quebec, Canada. 
%              worsley@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 the above copyright
%              notice appear 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 < 3
   df_resid=df
end
if nargin < 4 
   first_file=deblank(input_file_wresid);
   sinf=size(first_file,2);
   if input_file_wresid(sinf) == 'z'
      output_file_base=first_file(1,1:(sinf-7))
   else
      output_file_base=first_file(1,1:(sinf-4)) 
   end
end
if nargin < 5
   df_limit=4
end

d=fmris_read_image(input_file_wresid,0,0);
d.dim
n=d.dim(4);
numslices=d.dim(3);
J=d.dim(2);
I=d.dim(1);
Steps=d.vox
file_path=d.file_path;
ext=input_file_wresid((length(input_file_wresid)-2):length(input_file_wresid));

if strcmp(computer,'PCWIN')
   if strcmp(output_file_base(2:3),':\');file_path='';end
else
   if strcmp(output_file_base(1),'/');file_path='';end
end

fwhm_file=[deblank(output_file_base(1,:)) '_fwhm.' ext]
out_handle_fwhm.file_name=fwhm_file;
out_handle_fwhm.file_path=file_path;
out_handle_fwhm.dim=[I J numslices 0];
out_handle_fwhm.parent_file=input_file_wresid;
out_handle_fwhm.precision='float';  
out_handle_fwhm.descrip='';

% Technically, the methods are right only if df=Inf or df=df_resid,
% i.e. d=df_resid or d=df_resid-1, and approximate inbetween. 

if df_resid>df_limit
   % constants for arbitrary filter method:
   if df<Inf
      d=df_resid-df_resid/df;
      bias=exp((gammaln((d)/2-1/6)-gammaln((d)/2))+ ...
         (gammaln((d-1)/2-1/6)-gammaln((d-1)/2))+ ...
         (gammaln((d-2)/2-1/6)-gammaln((d-2)/2))+ ...
          gammaln((df+1)/2)-gammaln(df/2))*df^(-3/6);
   else
      bias=exp((gammaln((df_resid)/2-1/6)-gammaln((df_resid)/2))+ ...
         (gammaln((df_resid-1)/2-1/6)-gammaln((df_resid-1)/2))+ ...
         (gammaln((df_resid-2)/2-1/6)-gammaln((df_resid-2)/2)))*2^(-3/6);
   end
   correct=1/4;
else
   % constants filter aligned with axes method:
   if df<Inf
      d=df_resid-df_resid/df;
      bias=exp((gammaln((d)/2-1/6)-gammaln((d)/2))*3+ ...
         gammaln((df+1)/2)-gammaln(df/2))*df^(-3/6);
      correct=(df_resid-14/3*df_resid/df)/(d-1/3)/4;
   else
      bias=exp((gammaln((df_resid)/2-1/6)-gammaln((df_resid)/2))*3)*2^(-3/6);
      correct=(df_resid)/(d-1/3)/4;
   end
end
Step=abs(prod(Steps))^(1/3);
cons=sqrt(4*log(2))*Step/bias;

IJ=I*J;
Im=I-1;
Jm=J-1;
nxy=conv2(ones(Im,Jm),ones(2));

% Finds an estimate of the fwhm for each of the 8 cube corners surrounding
% a voxel, then averages. 

d=fmris_read_image(input_file_wresid,1,1:n);
u=squeeze(d.data);
ux=diff(u,1,1);
uy=diff(u,1,2);
Axx=sum(ux.^2,3);
Ayy=sum(uy.^2,3);
f=zeros(I,J);
ip=[0 1 0 1];
jp=[0 0 1 1];
is=[1 -1  1 -1];
js=[1  1 -1 -1];

for slice=1:(numslices-1)
   slice
   d=fmris_read_image(input_file_wresid,slice+1,1:n);
   up=squeeze(d.data);
   uz=up-u;
   Azz=sum(uz.^2,3);
   % The 4 upper cube corners:
   for index=1:4
      i=(1:Im)+ip(index);
      j=(1:Jm)+jp(index);
      axx=Axx(:,j);
      ayy=Ayy(i,:);
      azz=Azz(i,j);
      if df_resid>df_limit
         axy=sum(ux(:,j,:).*uy(i,:,:),3)*is(index)*js(index);
         axz=sum(ux(:,j,:).*uz(i,j,:),3)*is(index);
         ayz=sum(uy(i,:,:).*uz(i,j,:),3)*js(index);
         detlam=(axx.*ayy-axy.^2).*azz-(axz.*ayy-2*axy.*ayz).*axz-axx.*ayz.^2;
         s1=axx+ayy+azz;
         s2=axy+ayz+axz;
         % Continuity correction:
         detlamtrlaminvlam2=(s1+2*s2).*axx.*ayy.*azz+ ...
            4*(2*s1-s2).*axz.*axy.*ayz- ...
            (axx+4*(ayy-ayz+azz)).*axx.*ayz.^2- ...
            (ayy+4*(axx-axz+azz)).*ayy.*axz.^2- ...
            (azz+4*(axx-axy+ayy)).*azz.*axy.^2- ...
            2*(axx.*ayz.*(axz.*ayy+axy.*azz)+ayy.*azz.*axy.*axz);
         detlam=detlam+correct*detlamtrlaminvlam2;
      else
         detlam=axx.*ayy.*azz.*(1+correct*(axx+ayy+azz));
      end
      f(i,j)=f(i,j)+(detlam>0).*(detlam+(detlam<=0)).^(-1/6);
   end
   f=cons/((slice>1)+1)*f./nxy;
   out_handle_fwhm.data=f.*(f<50)+50*(f>=50);
   fmris_write_image(out_handle_fwhm,slice,1);
   
   f=zeros(I,J);
   u=up;
   ux=diff(u,1,1);
   uy=diff(u,1,2);
   Axx=sum(ux.^2,3);
   Ayy=sum(uy.^2,3);
   % The 4 lower cube corners:
   for index=1:4
      i=(1:Im)+ip(index);
      j=(1:Jm)+jp(index);
      axx=Axx(:,j);
      ayy=Ayy(i,:);
      azz=Azz(i,j);
      if df_resid>df_limit
         axy=sum(ux(:,j,:).*uy(i,:,:),3)*is(index)*js(index);
         axz=-sum(ux(:,j,:).*uz(i,j,:),3)*is(index);
         ayz=-sum(uy(i,:,:).*uz(i,j,:),3)*js(index);
         detlam=(axx.*ayy-axy.^2).*azz-(axz.*ayy-2*axy.*ayz).*axz-axx.*ayz.^2;
         s1=axx+ayy+azz;
         s2=axy+ayz+axz;
         % Continuity correction:
         detlamtrlaminvlam2=(s1+2*s2).*axx.*ayy.*azz+ ...
            4*(2*s1-s2).*axz.*axy.*ayz- ...
            (axx+4*(ayy-ayz+azz)).*axx.*ayz.^2- ...
            (ayy+4*(axx-axz+azz)).*ayy.*axz.^2- ...
            (azz+4*(axx-axy+ayy)).*azz.*axy.^2- ...
            2*(axx.*ayz.*(axz.*ayy+axy.*azz)+ayy.*azz.*axy.*axz);
         detlam=detlam+correct*detlamtrlaminvlam2;
      else
         detlam=axx.*ayy.*azz.*(1+correct*(axx+ayy+azz));
      end
      f(i,j)=f(i,j)+(detlam>0).*(detlam+(detlam<=0)).^(-1/6);
   end
end
slice=numslices
f=cons*f./nxy;
out_handle_fwhm.data=f.*(f<50)+50*(f>=50);
fmris_write_image(out_handle_fwhm,slice,1);

return