avw_homogeneity_correction.m

function avw = avw_homogeneity_correction(avw,threshold)

% avw_homogeneity_correction - 2D correction of MRI RF nonuniformity
% 
% avw = avw_homogeneity_correction(avw,[threshold])
%
% threshold - optional integer, lower intensity of volume
%             for high-pass cutoff (default = 5),
%             assumes T1 weighted volume
% 
%	a) fills in holes with the average intensity
%   b) generates a low frequency background image 
%   c) calculates a correction factor
%
% The algorithm was designed for correction of 2D slices.  
% It was adapted to iterate over slices in a 3D volume.  
% This is not an ideal solution. Most 3D cerebral MRI 
% volumes have RF inhomogeneity in the sagittal 
% direction (patient inferior to superior).  This function
% has been adapted here to work in this plane, given an
% input avw data struct from the mri_toolbox.  However,
% a better solution should work in 3D.
% This function is also severely limited by only handling
% volumes with even numbered slices.  If possible, it will
% be replaced with a better 3D function at some stage.  I
% advise, from limited experience, to use some other tools
% for this purpose (see the FSL tools, ie, FAST).
% 

% $Revision: 1.1 $ $Date: 2004/11/12 01:30:25 $

% Licence:  GNU GPL, no express or implied warranties
% History:  homocor.m  rev 0    4/8/00	          Gary Glover
%           @(#)vol_homocor.m    1.10  Kalina Christoff    2000-05-29
%                    - see Kalina's comments at the end of this function
%           07/2003, Darren.Weber_at_radiology.ucsf.edu
%                    - adapting to mri_toolbox
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


version = '[$Revision: 1.1 $]';
fprintf('\nAVW_HOMOGENEITY_CORRECTION [v%s]\n',version(12:16));  tic;



% This algorithm generates a low freq guassian matrix
% in such a way that it requires a square matrix (see below), 
% so this function assumes there is a square in-plane slice 
% matrix in the input avw.img file (most MRI is square in
% the in-plane acquisition slices).  Hence, it first 
% tries to identify the in-plane acquisition orientation.  
% There may be better ways to generate a low freq guassian 
% matrix (?).  [DLWeber]


xdim = double(avw.hdr.dime.dim(2));
ydim = double(avw.hdr.dime.dim(3));
zdim = double(avw.hdr.dime.dim(4));

% guess the slice acquisition orientation
% eg. avw.hdr.dime.dim=[4 256 256 124 1 0 0 0] 
% the inplane dims are the same in the axial plane
% --------------------------------------------------

if     xdim == ydim,
    orientation='axial';
    inplane_index=[1 2]; slice_index=3;
    
elseif ydim == zdim,
    orientation='sagittal';
    inplane_index=[2 3]; slice_index=1;
    
elseif xdim == zdim,
    orientation='coronal'; 
    inplane_index=[1 3]; slice_index=2;
end

fprintf('...assuming the slices were acquired in %s orientation\n', orientation);



if ~exist('threshold','var'),
    threshold = 5;
elseif isempty(threshold),
    threshold = 5;
end;
threshold = threshold * .01;


npix = size(avw.img,inplane_index(1));
npixHalf = npix/2;

if mod(npix,2) == 0,
    %OK
else,
    slices = npix
    error('sorry, this simple function can only handle even numbered slices.')
end

fprintf('...please wait - processing slice    ');

% begin looping through slices
% ----------------------------
for slice = 1:size(avw.img,slice_index);
    
    if     slice<10,   fprintf('\b%d',slice); 
    elseif slice<100,  fprintf('\b\b%d',slice); 
    elseif slice<1000, fprintf('\b\b\b%d',slice);
    end
    
    % squeeze out the singleton dimensions for coronal and sagittal
    if slice_index==1; img=squeeze(avw.img(slice,:,:)); end;  
    if slice_index==2; img=squeeze(avw.img(:,slice,:)); end;
    if slice_index==3; img=squeeze(avw.img(:,:,slice)); end;
    
    intensityMax = max(max(img));
    intensityThreshold = intensityMax*threshold;
    
    % fill in holes with image average intensity
    
    imgThresholdMask = (img > intensityThreshold);
    imgThresholdMaskSum = sum(sum(imgThresholdMask));
    
    imgMasked = img .* imgThresholdMask;
    imgAverage = sum(sum(imgMasked))/imgThresholdMaskSum;
    imgFilled = imgMasked + (1 - imgThresholdMask) .* imgAverage;
    
    %  make low freq image
    
    z = fftshift(fft2(imgFilled));
    
    
    a2 = 1/(npixHalf/8);
    y2 = (1:npixHalf).^2;  % use 15 as default win for guassian
    
    for x = 1:npixHalf,
        r2 = y2 + x*x;
        win1(x,:) = exp(-a2*r2);
    end
    
    % allocate win1 to 4 quadrants of win
    win(npixHalf+1:npix,npixHalf+1:npix) = win1;
    win(npixHalf+1:npix,npixHalf:-1:1) = win1;
    win(1:npixHalf,:) = win(npix:-1:npixHalf+1,:);
    
    zl = win.*z;
    
    imgLowFreq = abs(ifft2(fftshift(zl)));
    
    % image correction
    
    imgLowFreqAverage = sum(sum(imgLowFreq.*imgThresholdMask))/imgThresholdMaskSum;
    imgCorrectionFactor = (imgLowFreqAverage./imgLowFreq).*imgThresholdMask;
    imgCorrected = img .* imgCorrectionFactor;
    
    % update the volume with the corrected values
    % -------------------------------------------
    if slice_index==1; avw.img(slice,:,:) = imgCorrected; end
    if slice_index==2; avw.img(:,slice,:) = imgCorrected; end
    if slice_index==3; avw.img(:,:,slice) = imgCorrected; end
    
    % end slice loop
end

t=toc; fprintf('...done (%5.2f sec).\n',t);

return


% *Corecting for inhomogeneities in anatomical images*
% 
% *
% ------------------------------------------------------------------------
% *
% 
% Here is a small but rather useful program developed by Gary Glover for
% correcting inhomogeneities in anatomical images.
% 
% Inhomogeneities in signal intensity can be somewhat problematic at high
% magentic fields (e.g. 3 Tesla) and may lead to reduced quality of the
% image due to slow spatial modulation in overall intensity.
% 
% Gary's program identifies and removes such low-frequency intensity
% modulations. The transformations are computed for each slice separately.
% 
% I have modified the program so that it can handle SPM analyze format
% images. From within our lab it can be used directly by typing
% "vol_homocor" at the matlab5 prompt, or otherwise it can be downloaded
% from here
% 
% vol_homocor.m
% 
% and saved as a file named "vol_homocor.m" in your matlab directory.
% 
% Below is an example of an anatomical image that benefitted a lot from
% correction. Of course, as it goes, some anatomical images are less
% problematic than others, so the usefullness of correction may vary from
% subject to subject.
% 
% ------------------------------------------------------------------------
% 
% Kalina Christoff / 2000-05-27
% kalina@psych.stanford.edu