📄 avw_homogeneity_correction.html
字号:
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN" "http://www.w3.org/TR/REC-html40/loose.dtd"><html><head> <title>Description of avw_homogeneity_correction</title> <meta name="keywords" content="avw_homogeneity_correction"> <meta name="description" content="avw_homogeneity_correction - 2D correction of MRI RF nonuniformity"> <meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1"> <meta name="generator" content="m2html © 2003 Guillaume Flandin"> <meta name="robots" content="index, follow"> <link type="text/css" rel="stylesheet" href="../m2html.css"></head><body><a name="_top"></a><div><a href="../index.html">Home</a> > <a href="index.html">mri_toolbox</a> > avw_homogeneity_correction.m</div><!--<table width="100%"><tr><td align="left"><a href="../index.html"><img alt="<" border="0" src="../left.png"> Master index</a></td><td align="right"><a href="index.html">Index for mri_toolbox <img alt=">" border="0" src="../right.png"></a></td></tr></table>--><h1>avw_homogeneity_correction</h1><h2><a name="_name"></a>PURPOSE <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2><div class="box"><strong>avw_homogeneity_correction - 2D correction of MRI RF nonuniformity</strong></div><h2><a name="_synopsis"></a>SYNOPSIS <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2><div class="box"><strong>function avw = avw_homogeneity_correction(avw,threshold) </strong></div><h2><a name="_description"></a>DESCRIPTION <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2><div class="fragment"><pre class="comment"> 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).</pre></div><!-- crossreference --><h2><a name="_cross"></a>CROSS-REFERENCE INFORMATION <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>This function calls:<ul style="list-style-image:url(../matlabicon.gif)"></ul>This function is called by:<ul style="list-style-image:url(../matlabicon.gif)"></ul><!-- crossreference --><h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2><div class="fragment"><pre>0001 <a name="_sub0" href="#_subfunctions" class="code">function avw = avw_homogeneity_correction(avw,threshold)</a>0002 0003 <span class="comment">% avw_homogeneity_correction - 2D correction of MRI RF nonuniformity</span>0004 <span class="comment">%</span>0005 <span class="comment">% avw = avw_homogeneity_correction(avw,[threshold])</span>0006 <span class="comment">%</span>0007 <span class="comment">% threshold - optional integer, lower intensity of volume</span>0008 <span class="comment">% for high-pass cutoff (default = 5),</span>0009 <span class="comment">% assumes T1 weighted volume</span>0010 <span class="comment">%</span>0011 <span class="comment">% a) fills in holes with the average intensity</span>0012 <span class="comment">% b) generates a low frequency background image</span>0013 <span class="comment">% c) calculates a correction factor</span>0014 <span class="comment">%</span>0015 <span class="comment">% The algorithm was designed for correction of 2D slices.</span>0016 <span class="comment">% It was adapted to iterate over slices in a 3D volume.</span>0017 <span class="comment">% This is not an ideal solution. Most 3D cerebral MRI</span>0018 <span class="comment">% volumes have RF inhomogeneity in the sagittal</span>0019 <span class="comment">% direction (patient inferior to superior). This function</span>0020 <span class="comment">% has been adapted here to work in this plane, given an</span>0021 <span class="comment">% input avw data struct from the mri_toolbox. However,</span>0022 <span class="comment">% a better solution should work in 3D.</span>0023 <span class="comment">% This function is also severely limited by only handling</span>0024 <span class="comment">% volumes with even numbered slices. If possible, it will</span>0025 <span class="comment">% be replaced with a better 3D function at some stage. I</span>0026 <span class="comment">% advise, from limited experience, to use some other tools</span>0027 <span class="comment">% for this purpose (see the FSL tools, ie, FAST).</span>0028 <span class="comment">%</span>0029 0030 <span class="comment">% $Revision: 1.2 $ $Date: 2004/02/07 01:41:51 $</span>0031 0032 <span class="comment">% Licence: GNU GPL, no express or implied warranties</span>0033 <span class="comment">% History: homocor.m rev 0 4/8/00 Gary Glover</span>0034 <span class="comment">% @(#)vol_homocor.m 1.10 Kalina Christoff 2000-05-29</span>0035 <span class="comment">% - see Kalina's comments at the end of this function</span>0036 <span class="comment">% 07/2003, Darren.Weber_at_radiology.ucsf.edu</span>0037 <span class="comment">% - adapting to mri_toolbox</span>0038 <span class="comment">%</span>0039 <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>0040 0041 0042 version = <span class="string">'[$Revision: 1.2 $]'</span>;0043 fprintf(<span class="string">'\nAVW_HOMOGENEITY_CORRECTION [v%s]\n'</span>,version(12:16)); tic;0044 0045 0046 0047 <span class="comment">% This algorithm generates a low freq guassian matrix</span>0048 <span class="comment">% in such a way that it requires a square matrix (see below),</span>0049 <span class="comment">% so this function assumes there is a square in-plane slice</span>0050 <span class="comment">% matrix in the input avw.img file (most MRI is square in</span>0051 <span class="comment">% the in-plane acquisition slices). Hence, it first</span>0052 <span class="comment">% tries to identify the in-plane acquisition orientation.</span>0053 <span class="comment">% There may be better ways to generate a low freq guassian</span>0054 <span class="comment">% matrix (?). [DLWeber]</span>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 <span class="comment">% guess the slice acquisition orientation</span>0062 <span class="comment">% eg. avw.hdr.dime.dim=[4 256 256 124 1 0 0 0]</span>0063 <span class="comment">% the inplane dims are the same in the axial plane</span>0064 <span class="comment">% --------------------------------------------------</span>0065 0066 <span class="keyword">if</span> xdim == ydim,0067 orientation=<span class="string">'axial'</span>;0068 inplane_index=[1 2]; slice_index=3;0069 0070 <span class="keyword">elseif</span> ydim == zdim,0071 orientation=<span class="string">'sagittal'</span>;0072 inplane_index=[2 3]; slice_index=1;0073 0074 <span class="keyword">elseif</span> xdim == zdim,0075 orientation=<span class="string">'coronal'</span>; 0076 inplane_index=[1 3]; slice_index=2;0077 <span class="keyword">end</span>0078 0079 fprintf(<span class="string">'...assuming the slices were acquired in %s orientation\n'</span>, orientation);0080 0081 0082 0083 <span class="keyword">if</span> ~exist(<span class="string">'threshold'</span>,<span class="string">'var'</span>),0084 threshold = 5;0085 <span class="keyword">elseif</span> isempty(threshold),0086 threshold = 5;0087 <span class="keyword">end</span>;0088 threshold = threshold * .01;0089 0090 0091 npix = size(avw.img,inplane_index(1));0092 npixHalf = npix/2;0093 0094 <span class="keyword">if</span> mod(npix,2) == 0,0095 <span class="comment">%OK</span>0096 <span class="keyword">else</span>,0097 slices = npix0098 error(<span class="string">'sorry, this simple function can only handle even numbered slices.'</span>)0099 <span class="keyword">end</span>0100 0101 fprintf(<span class="string">'...please wait - processing slice '</span>);0102 0103 <span class="comment">% begin looping through slices</span>0104 <span class="comment">% ----------------------------</span>0105 <span class="keyword">for</span> slice = 1:size(avw.img,slice_index);0106 0107 <span class="keyword">if</span> slice<10, fprintf(<span class="string">'\b%d'</span>,slice); 0108 <span class="keyword">elseif</span> slice<100, fprintf(<span class="string">'\b\b%d'</span>,slice); 0109 <span class="keyword">elseif</span> slice<1000, fprintf(<span class="string">'\b\b\b%d'</span>,slice);0110 <span class="keyword">end</span>0111 0112 <span class="comment">% squeeze out the singleton dimensions for coronal and sagittal</span>0113 <span class="keyword">if</span> slice_index==1; img=squeeze(avw.img(slice,:,:)); <span class="keyword">end</span>; 0114 <span class="keyword">if</span> slice_index==2; img=squeeze(avw.img(:,slice,:)); <span class="keyword">end</span>;0115 <span class="keyword">if</span> slice_index==3; img=squeeze(avw.img(:,:,slice)); <span class="keyword">end</span>;0116 0117 intensityMax = max(max(img));0118 intensityThreshold = intensityMax*threshold;0119 0120 <span class="comment">% fill in holes with image average intensity</span>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 <span class="comment">% make low freq image</span>0130 0131 z = fftshift(fft2(imgFilled));0132 0133 0134 a2 = 1/(npixHalf/8);0135 y2 = (1:npixHalf).^2; <span class="comment">% use 15 as default win for guassian</span>0136 0137 <span class="keyword">for</span> x = 1:npixHalf,0138 r2 = y2 + x*x;0139 win1(x,:) = exp(-a2*r2);0140 <span class="keyword">end</span>0141 0142 <span class="comment">% allocate win1 to 4 quadrants of win</span>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 <span class="comment">% image correction</span>0152 0153 imgLowFreqAverage = sum(sum(imgLowFreq.*imgThresholdMask))/imgThresholdMaskSum;0154 imgCorrectionFactor = (imgLowFreqAverage./imgLowFreq).*imgThresholdMask;0155 imgCorrected = img .* imgCorrectionFactor;0156 0157 <span class="comment">% update the volume with the corrected values</span>0158 <span class="comment">% -------------------------------------------</span>0159 <span class="keyword">if</span> slice_index==1; avw.img(slice,:,:) = imgCorrected; <span class="keyword">end</span>0160 <span class="keyword">if</span> slice_index==2; avw.img(:,slice,:) = imgCorrected; <span class="keyword">end</span>0161 <span class="keyword">if</span> slice_index==3; avw.img(:,:,slice) = imgCorrected; <span class="keyword">end</span>0162 0163 <span class="comment">% end slice loop</span>0164 <span class="keyword">end</span>0165 0166 t=toc; fprintf(<span class="string">'...done (%5.2f sec).\n'</span>,t);0167 0168 <span class="keyword">return</span>0169 0170 0171 <span class="comment">% *Corecting for inhomogeneities in anatomical images*</span>0172 <span class="comment">%</span>0173 <span class="comment">% *</span>0174 <span class="comment">% ------------------------------------------------------------------------</span>0175 <span class="comment">% *</span>0176 <span class="comment">%</span>0177 <span class="comment">% Here is a small but rather useful program developed by Gary Glover for</span>0178 <span class="comment">% correcting inhomogeneities in anatomical images.</span>0179 <span class="comment">%</span>0180 <span class="comment">% Inhomogeneities in signal intensity can be somewhat problematic at high</span>0181 <span class="comment">% magentic fields (e.g. 3 Tesla) and may lead to reduced quality of the</span>0182 <span class="comment">% image due to slow spatial modulation in overall intensity.</span>0183 <span class="comment">%</span>0184 <span class="comment">% Gary's program identifies and removes such low-frequency intensity</span>0185 <span class="comment">% modulations. The transformations are computed for each slice separately.</span>0186 <span class="comment">%</span>0187 <span class="comment">% I have modified the program so that it can handle SPM analyze format</span>0188 <span class="comment">% images. From within our lab it can be used directly by typing</span>0189 <span class="comment">% "vol_homocor" at the matlab5 prompt, or otherwise it can be downloaded</span>0190 <span class="comment">% from here</span>0191 <span class="comment">%</span>0192 <span class="comment">% vol_homocor.m</span>0193 <span class="comment">%</span>0194 <span class="comment">% and saved as a file named "vol_homocor.m" in your matlab directory.</span>0195 <span class="comment">%</span>0196 <span class="comment">% Below is an example of an anatomical image that benefitted a lot from</span>0197 <span class="comment">% correction. Of course, as it goes, some anatomical images are less</span>0198 <span class="comment">% problematic than others, so the usefullness of correction may vary from</span>0199 <span class="comment">% subject to subject.</span>0200 <span class="comment">%</span>0201 <span class="comment">% ------------------------------------------------------------------------</span>0202 <span class="comment">%</span>0203 <span class="comment">% Kalina Christoff / 2000-05-27</span>0204 <span class="comment">% kalina@psych.stanford.edu</span></pre></div><hr><address>Generated on Fri 21-May-2004 12:38:21 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/">m2html</a></strong> © 2003</address></body></html>⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -