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).0001 function avw = avw_homogeneity_correction(avw,threshold) 0002 0003 % avw_homogeneity_correction - 2D correction of MRI RF nonuniformity 0004 % 0005 % avw = avw_homogeneity_correction(avw,[threshold]) 0006 % 0007 % threshold - optional integer, lower intensity of volume 0008 % for high-pass cutoff (default = 5), 0009 % assumes T1 weighted volume 0010 % 0011 % a) fills in holes with the average intensity 0012 % b) generates a low frequency background image 0013 % c) calculates a correction factor 0014 % 0015 % The algorithm was designed for correction of 2D slices. 0016 % It was adapted to iterate over slices in a 3D volume. 0017 % This is not an ideal solution. Most 3D cerebral MRI 0018 % volumes have RF inhomogeneity in the sagittal 0019 % direction (patient inferior to superior). This function 0020 % has been adapted here to work in this plane, given an 0021 % input avw data struct from the mri_toolbox. However, 0022 % a better solution should work in 3D. 0023 % This function is also severely limited by only handling 0024 % volumes with even numbered slices. If possible, it will 0025 % be replaced with a better 3D function at some stage. I 0026 % advise, from limited experience, to use some other tools 0027 % for this purpose (see the FSL tools, ie, FAST). 0028 % 0029 0030 % $Revision: 1.1 $ $Date: 2004/11/12 01:30:25 $ 0031 0032 % Licence: GNU GPL, no express or implied warranties 0033 % History: homocor.m rev 0 4/8/00 Gary Glover 0034 % @(#)vol_homocor.m 1.10 Kalina Christoff 2000-05-29 0035 % - see Kalina's comments at the end of this function 0036 % 07/2003, Darren.Weber_at_radiology.ucsf.edu 0037 % - adapting to mri_toolbox 0038 % 0039 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0040 0041 0042 version = '[$Revision: 1.1 $]'; 0043 fprintf('\nAVW_HOMOGENEITY_CORRECTION [v%s]\n',version(12:16)); tic; 0044 0045 0046 0047 % This algorithm generates a low freq guassian matrix 0048 % in such a way that it requires a square matrix (see below), 0049 % so this function assumes there is a square in-plane slice 0050 % matrix in the input avw.img file (most MRI is square in 0051 % the in-plane acquisition slices). Hence, it first 0052 % tries to identify the in-plane acquisition orientation. 0053 % There may be better ways to generate a low freq guassian 0054 % matrix (?). [DLWeber] 0055 0056 0057 xdim = double(avw.hdr.dime.dim(2)); 0058 ydim = double(avw.hdr.dime.dim(3)); 0059 zdim = double(avw.hdr.dime.dim(4)); 0060 0061 % guess the slice acquisition orientation 0062 % eg. avw.hdr.dime.dim=[4 256 256 124 1 0 0 0] 0063 % the inplane dims are the same in the axial plane 0064 % -------------------------------------------------- 0065 0066 if xdim == ydim, 0067 orientation='axial'; 0068 inplane_index=[1 2]; slice_index=3; 0069 0070 elseif ydim == zdim, 0071 orientation='sagittal'; 0072 inplane_index=[2 3]; slice_index=1; 0073 0074 elseif xdim == zdim, 0075 orientation='coronal'; 0076 inplane_index=[1 3]; slice_index=2; 0077 end 0078 0079 fprintf('...assuming the slices were acquired in %s orientation\n', orientation); 0080 0081 0082 0083 if ~exist('threshold','var'), 0084 threshold = 5; 0085 elseif isempty(threshold), 0086 threshold = 5; 0087 end; 0088 threshold = threshold * .01; 0089 0090 0091 npix = size(avw.img,inplane_index(1)); 0092 npixHalf = npix/2; 0093 0094 if mod(npix,2) == 0, 0095 %OK 0096 else, 0097 slices = npix 0098 error('sorry, this simple function can only handle even numbered slices.') 0099 end 0100 0101 fprintf('...please wait - processing slice '); 0102 0103 % begin looping through slices 0104 % ---------------------------- 0105 for slice = 1:size(avw.img,slice_index); 0106 0107 if slice<10, fprintf('\b%d',slice); 0108 elseif slice<100, fprintf('\b\b%d',slice); 0109 elseif slice<1000, fprintf('\b\b\b%d',slice); 0110 end 0111 0112 % squeeze out the singleton dimensions for coronal and sagittal 0113 if slice_index==1; img=squeeze(avw.img(slice,:,:)); end; 0114 if slice_index==2; img=squeeze(avw.img(:,slice,:)); end; 0115 if slice_index==3; img=squeeze(avw.img(:,:,slice)); end; 0116 0117 intensityMax = max(max(img)); 0118 intensityThreshold = intensityMax*threshold; 0119 0120 % fill in holes with image average intensity 0121 0122 imgThresholdMask = (img > intensityThreshold); 0123 imgThresholdMaskSum = sum(sum(imgThresholdMask)); 0124 0125 imgMasked = img .* imgThresholdMask; 0126 imgAverage = sum(sum(imgMasked))/imgThresholdMaskSum; 0127 imgFilled = imgMasked + (1 - imgThresholdMask) .* imgAverage; 0128 0129 % make low freq image 0130 0131 z = fftshift(fft2(imgFilled)); 0132 0133 0134 a2 = 1/(npixHalf/8); 0135 y2 = (1:npixHalf).^2; % use 15 as default win for guassian 0136 0137 for x = 1:npixHalf, 0138 r2 = y2 + x*x; 0139 win1(x,:) = exp(-a2*r2); 0140 end 0141 0142 % allocate win1 to 4 quadrants of win 0143 win(npixHalf+1:npix,npixHalf+1:npix) = win1; 0144 win(npixHalf+1:npix,npixHalf:-1:1) = win1; 0145 win(1:npixHalf,:) = win(npix:-1:npixHalf+1,:); 0146 0147 zl = win.*z; 0148 0149 imgLowFreq = abs(ifft2(fftshift(zl))); 0150 0151 % image correction 0152 0153 imgLowFreqAverage = sum(sum(imgLowFreq.*imgThresholdMask))/imgThresholdMaskSum; 0154 imgCorrectionFactor = (imgLowFreqAverage./imgLowFreq).*imgThresholdMask; 0155 imgCorrected = img .* imgCorrectionFactor; 0156 0157 % update the volume with the corrected values 0158 % ------------------------------------------- 0159 if slice_index==1; avw.img(slice,:,:) = imgCorrected; end 0160 if slice_index==2; avw.img(:,slice,:) = imgCorrected; end 0161 if slice_index==3; avw.img(:,:,slice) = imgCorrected; end 0162 0163 % end slice loop 0164 end 0165 0166 t=toc; fprintf('...done (%5.2f sec).\n',t); 0167 0168 return 0169 0170 0171 % *Corecting for inhomogeneities in anatomical images* 0172 % 0173 % * 0174 % ------------------------------------------------------------------------ 0175 % * 0176 % 0177 % Here is a small but rather useful program developed by Gary Glover for 0178 % correcting inhomogeneities in anatomical images. 0179 % 0180 % Inhomogeneities in signal intensity can be somewhat problematic at high 0181 % magentic fields (e.g. 3 Tesla) and may lead to reduced quality of the 0182 % image due to slow spatial modulation in overall intensity. 0183 % 0184 % Gary's program identifies and removes such low-frequency intensity 0185 % modulations. The transformations are computed for each slice separately. 0186 % 0187 % I have modified the program so that it can handle SPM analyze format 0188 % images. From within our lab it can be used directly by typing 0189 % "vol_homocor" at the matlab5 prompt, or otherwise it can be downloaded 0190 % from here 0191 % 0192 % vol_homocor.m 0193 % 0194 % and saved as a file named "vol_homocor.m" in your matlab directory. 0195 % 0196 % Below is an example of an anatomical image that benefitted a lot from 0197 % correction. Of course, as it goes, some anatomical images are less 0198 % problematic than others, so the usefullness of correction may vary from 0199 % subject to subject. 0200 % 0201 % ------------------------------------------------------------------------ 0202 % 0203 % Kalina Christoff / 2000-05-27 0204 % kalina@psych.stanford.edu