📄 avw_shrinkwrap.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_shrinkwrap</title> <meta name="keywords" content="avw_shrinkwrap"> <meta name="description" content="avw_shrinkwrap - Tesselate the surface of a 3D Analyze 7.5 avw struct"> <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_shrinkwrap.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_shrinkwrap</h1><h2><a name="_name"></a>PURPOSE <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2><div class="box"><strong>avw_shrinkwrap - Tesselate the surface of a 3D Analyze 7.5 avw struct</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 [FV, Edges] = avw_shrinkwrap(avw,FV,interpVal,fitval,fittol,fititer,fitchange,fitvattr) </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_shrinkwrap - Tesselate the surface of a 3D Analyze 7.5 avw struct
[FV, Edges] = avw_shrinkwrap(avw,FV,interpVal,...
fitval,fittol,fititer,fitchange,fitvattr)
avw - an Analyze 7.5 data struct (see avw_read)
FV - input tesselation; if empty, sphere tesselation
is created. FV has fields FV.vertices, FV.faces
interpVal - use radial interpolation (faster) or incremental
radial shrink (slower), 0|1 (default 1, faster)
fitval - image intensity to shrink wrap (default 20)
fittol - image intensity tolerance (default 5)
fititer - max number of iterations to fit (default 200)
fitchange - least sig. change in intensity values
between fit iterations (default 2)
fitvattr - vertex attraction (constraint), 0:1, smaller
values are less constraint; close to 0 for
no constraint is useful when dealing with
binary volumes, otherwise 0.4 (40%) seems OK
FV - a struct with 2562 vertices and 5120 faces
Edges - a [2562,2562] matrix of edge connectivity for FV
This function has been developed to find the scalp surface
for MRI of the human head. It is not a sophisticated, robust
algorithm!
It starts with a sphere tesselation (large radius) and moves
each vertex point toward the center of the volume until it
lies at or near the fitval. The vertices are constrained to
move only along the radial projection from the origin and they
are also required to stay within a small distance of their
neighbours. The function is not optimised for speed, but
it should produce reasonable results.
Example of creating a scalp tesselation for SPM T1 MRI template:
avw = avw_read('T1');
FV = avw_shrinkwrap(avw,[],0,0,[],intensity,5.0,50,0.5,0.4);
patch('vertices',FV.vertices,'faces',FV.faces,'facecolor',[.6 .6 .6]);
The output vertex coordinates are in meters with an origin at (0,0,0),
which lies at the center of the MRI volume. This function uses the
avw.hdr.dime.pixdim values to convert from voxel coordinates into
meters.
See also: ISOSURFACE, SPHERE_TRI, MESH_REFINE_TRI4, MESH_BEM_SHELLS_FUNC, MESH_BEM_SHELLS_SCRIPT</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)"><li><a href="avw_center.html" class="code" title="function center = avw_center(avw)">avw_center</a> avw_center - find center of a volume</li></ul>This function is called by:<ul style="list-style-image:url(../matlabicon.gif)"><li><a href="mesh_3shell_script.html" class="code" title="">mesh_3shell_script</a> mesh_3shell_script</li></ul><!-- crossreference --><h2><a name="_subfunctions"></a>SUBFUNCTIONS <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2><ul style="list-style-image:url(../matlabicon.gif)"><li><a href="#_sub1" class="code">function [fit] = estimate_scalp(FV,vol,origin,fit),</a></li><li><a href="#_sub2" class="code">function [FV, intensityAtR, R] = locate_val(FV,vol,origin,fit),</a></li><li><a href="#_sub3" class="code">function [FV, intensityAtR, R] = shrink_wrap(FV,vol,origin,fit),</a></li><li><a href="#_sub4" class="code">function [FV] = constrain_vertex(FV,index,origin),</a></li><li><a href="#_sub5" class="code">function [val] = vol_val(vol,x,y,z),</a></li><li><a href="#_sub6" class="code">function [FV] = vertex_outliers(FV, R, origin),</a></li></ul><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 [FV, Edges] = avw_shrinkwrap(avw,FV,interpVal,</a><span class="keyword">...</span>0002 fitval,fittol,fititer,fitchange,fitvattr)0003 0004 <span class="comment">% avw_shrinkwrap - Tesselate the surface of a 3D Analyze 7.5 avw struct</span>0005 <span class="comment">%</span>0006 <span class="comment">% [FV, Edges] = avw_shrinkwrap(avw,FV,interpVal,...</span>0007 <span class="comment">% fitval,fittol,fititer,fitchange,fitvattr)</span>0008 <span class="comment">%</span>0009 <span class="comment">% avw - an Analyze 7.5 data struct (see avw_read)</span>0010 <span class="comment">% FV - input tesselation; if empty, sphere tesselation</span>0011 <span class="comment">% is created. FV has fields FV.vertices, FV.faces</span>0012 <span class="comment">% interpVal - use radial interpolation (faster) or incremental</span>0013 <span class="comment">% radial shrink (slower), 0|1 (default 1, faster)</span>0014 <span class="comment">% fitval - image intensity to shrink wrap (default 20)</span>0015 <span class="comment">% fittol - image intensity tolerance (default 5)</span>0016 <span class="comment">% fititer - max number of iterations to fit (default 200)</span>0017 <span class="comment">% fitchange - least sig. change in intensity values</span>0018 <span class="comment">% between fit iterations (default 2)</span>0019 <span class="comment">% fitvattr - vertex attraction (constraint), 0:1, smaller</span>0020 <span class="comment">% values are less constraint; close to 0 for</span>0021 <span class="comment">% no constraint is useful when dealing with</span>0022 <span class="comment">% binary volumes, otherwise 0.4 (40%) seems OK</span>0023 <span class="comment">%</span>0024 <span class="comment">% FV - a struct with 2562 vertices and 5120 faces</span>0025 <span class="comment">% Edges - a [2562,2562] matrix of edge connectivity for FV</span>0026 <span class="comment">%</span>0027 <span class="comment">% This function has been developed to find the scalp surface</span>0028 <span class="comment">% for MRI of the human head. It is not a sophisticated, robust</span>0029 <span class="comment">% algorithm!</span>0030 <span class="comment">%</span>0031 <span class="comment">% It starts with a sphere tesselation (large radius) and moves</span>0032 <span class="comment">% each vertex point toward the center of the volume until it</span>0033 <span class="comment">% lies at or near the fitval. The vertices are constrained to</span>0034 <span class="comment">% move only along the radial projection from the origin and they</span>0035 <span class="comment">% are also required to stay within a small distance of their</span>0036 <span class="comment">% neighbours. The function is not optimised for speed, but</span>0037 <span class="comment">% it should produce reasonable results.</span>0038 <span class="comment">%</span>0039 <span class="comment">% Example of creating a scalp tesselation for SPM T1 MRI template:</span>0040 <span class="comment">%</span>0041 <span class="comment">% avw = avw_read('T1');</span>0042 <span class="comment">% FV = avw_shrinkwrap(avw,[],0,0,[],intensity,5.0,50,0.5,0.4);</span>0043 <span class="comment">% patch('vertices',FV.vertices,'faces',FV.faces,'facecolor',[.6 .6 .6]);</span>0044 <span class="comment">%</span>0045 <span class="comment">% The output vertex coordinates are in meters with an origin at (0,0,0),</span>0046 <span class="comment">% which lies at the center of the MRI volume. This function uses the</span>0047 <span class="comment">% avw.hdr.dime.pixdim values to convert from voxel coordinates into</span>0048 <span class="comment">% meters.</span>0049 <span class="comment">%</span>0050 <span class="comment">% See also: ISOSURFACE, SPHERE_TRI, MESH_REFINE_TRI4,</span>0051 <span class="comment">% MESH_BEM_SHELLS_FUNC, MESH_BEM_SHELLS_SCRIPT</span>0052 <span class="comment">%</span>0053 0054 <span class="comment">% $Revision: 1.2 $ $Date: 2004/02/07 01:41:51 $</span>0055 0056 <span class="comment">% Licence: GNU GPL, no implied or express warranties</span>0057 <span class="comment">% History: 08/2003, Darren.Weber_at_radiology.ucsf.edu</span>0058 <span class="comment">% - adapted from mesh_shrinkwrap</span>0059 <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>0060 0061 <span class="comment">% Parse arguments</span>0062 0063 <span class="keyword">if</span> ~exist(<span class="string">'avw'</span>,<span class="string">'var'</span>), error(<span class="string">'...no input volume\n'</span>);0064 <span class="keyword">elseif</span> isempty(avw), error(<span class="string">'...empty input volume\n'</span>);0065 <span class="keyword">end</span>0066 0067 <span class="comment">% Find approximate center of volume</span>0068 center = <a href="avw_center.html" class="code" title="function center = avw_center(avw)">avw_center</a>(avw);0069 origin = center.abs.voxels;0070 0071 version = <span class="string">'[$Revision: 1.2 $]'</span>;0072 fprintf(<span class="string">'AVW_SHRINKWRAP [v%s]\n'</span>,version(12:16)); tic;0073 0074 <span class="keyword">if</span> ~exist(<span class="string">'interpVal'</span>,<span class="string">'var'</span>), interpVal = 0;0075 <span class="keyword">elseif</span> isempty(interpVal), interpVal = 0;0076 <span class="keyword">end</span>0077 0078 <span class="keyword">if</span> ~exist(<span class="string">'fitval'</span>,<span class="string">'var'</span>), fit.val = 20;0079 <span class="keyword">elseif</span> isempty(fitval), fit.val = 20;0080 <span class="keyword">else</span> fit.val = fitval;0081 <span class="keyword">end</span>0082 0083 <span class="keyword">if</span> ~exist(<span class="string">'fittol'</span>,<span class="string">'var'</span>), fit.tol = 5;0084 <span class="keyword">elseif</span> isempty(fittol), fit.tol = 5;0085 <span class="keyword">else</span> fit.tol = fittol;0086 <span class="keyword">end</span>0087 0088 <span class="keyword">if</span> fit.val <= fit.tol,0089 error(<span class="string">'...must use fit tolerance < fit value\n'</span>);0090 <span class="keyword">end</span>0091 0092 <span class="keyword">if</span> ~exist(<span class="string">'fititer'</span>,<span class="string">'var'</span>), fit.iter = 200;0093 <span class="keyword">elseif</span> isempty(fititer), fit.iter = 200;0094 <span class="keyword">else</span> fit.iter = fititer;0095 <span class="keyword">end</span>0096 0097 <span class="keyword">if</span> ~exist(<span class="string">'fitchange'</span>,<span class="string">'var'</span>),fit.change = 2;0098 <span class="keyword">elseif</span> isempty(fitchange), fit.change = 2;0099 <span class="keyword">else</span> fit.change = fitchange;0100 <span class="keyword">end</span>0101 0102 <span class="keyword">if</span> ~exist(<span class="string">'fitvattr'</span>,<span class="string">'var'</span>), fit.vattr = 0.4;0103 <span class="keyword">elseif</span> isempty(fitvattr), fit.vattr = 0.4;0104 <span class="keyword">else</span> fit.vattr = fitvattr;0105 <span class="keyword">end</span>0106 <span class="keyword">if</span> fit.vattr > 1,0107 fprintf(<span class="string">'...fit vertattr (v) must be 0 <= v <= 1, setting v = 1\n'</span>);0108 fit.vattr = 1;0109 <span class="keyword">end</span>0110 <span class="keyword">if</span> fit.vattr < 0,0111 fprintf(<span class="string">'...fit vertattr (v) must be 0 <= v <= 1, setting v = 0.\n'</span>);0112 fit.vattr = 0;0113 <span class="keyword">end</span>0114 0115 0116 <span class="comment">% get size of volume, in voxels</span>0117 [xdim,ydim,zdim] = size(avw.img);0118 0119 0120 <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>0121 <span class="comment">% Check whether to create a sphere tesselation</span>0122 <span class="comment">% or use an input tesselation as the start point</span>0123 0124 sphere = 0;0125 <span class="keyword">if</span> ~exist(<span class="string">'FV'</span>,<span class="string">'var'</span>),0126 sphere = 1;0127 <span class="keyword">elseif</span> ~isfield(FV,<span class="string">'vertices'</span>),0128 sphere = 1;0129 <span class="keyword">elseif</span> ~isfield(FV,<span class="string">'faces'</span>),0130 sphere = 1;0131 <span class="keyword">elseif</span> isempty(FV.vertices),0132 sphere = 1;0133 <span class="keyword">elseif</span> isempty(FV.faces),0134 sphere = 1;0135 <span class="keyword">end</span>0136 <span class="keyword">if</span> sphere,0137 <span class="comment">% Create a sphere tesselation to encompass the volume</span>0138 diameter = max([xdim ydim zdim]);0139 radius = diameter / 1.5;0140 FV = sphere_tri(<span class="string">'ico'</span>,4,radius); <span class="comment">% 2562 vertices</span>0141 <span class="comment">% Shift the center of the sphere to the center of the volume</span>0142 FV.vertices = FV.vertices + repmat(origin,size(FV.vertices,1),1);0143 <span class="keyword">else</span>0144 fprintf(<span class="string">'...using input FV tesselation...\n'</span>);0145 <span class="keyword">end</span>0146 0147 <span class="comment">% the 'edge' matrix is the connectivity of all vertices,</span>0148 <span class="comment">% used to find neighbours during movement of vertices,</span>0149 <span class="comment">% including smoothing the tesselation</span>0150 FV.edge = mesh_edges(FV);0151 0152 0153 0154 0155 0156 <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>0157 <span class="comment">% Estimate scalp intensity</span>0158 fit = <a href="#_sub1" class="code" title="subfunction [fit] = estimate_scalp(FV,vol,origin,fit),">estimate_scalp</a>(FV,avw.img,origin,fit);0159 0160 0161 0162 <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>0163 <span class="comment">% Now begin recursion</span>0164 fprintf(<span class="string">'...fitting...\n'</span>); tic;0165 0166 i = 1;0167 Fminima = 0;0168 intensityAtRMean = [0 0];0169 0170 <span class="keyword">while</span> i <= fit.iter,0171 0172 <span class="keyword">if</span> interpVal,0173 <span class="comment">% use radial interpolation method, moving directly</span>0174 <span class="comment">% to the intensity value nearest correct intensity</span>0175 [FV, intensityAtR, R] = <a href="#_sub2" class="code" title="subfunction [FV, intensityAtR, R] = locate_val(FV,vol,origin,fit),">locate_val</a>(FV,avw.img,origin,fit);0176 <span class="keyword">else</span>0177 <span class="comment">% use incremental method, moving along radial line</span>0178 <span class="comment">% gradually until finding correct intensity</span>0179 [FV, intensityAtR, R] = <a href="#_sub3" class="code" title="subfunction [FV, intensityAtR, R] = shrink_wrap(FV,vol,origin,fit),">shrink_wrap</a>(FV,avw.img,origin,fit);0180 <span class="keyword">end</span>0181 0182 intensityAtRMean = [ intensityAtRMean(2), mean(intensityAtR) ] ;0183 0184 fprintf(<span class="string">'...distance: mean = %8.4f voxels, std = %8.4f voxels\n'</span>,mean(R),std(R));0185 fprintf(<span class="string">'...intensity: mean = %8.4f, std = %8.4f\n'</span>,<span class="keyword">...</span>0186 mean(intensityAtR),std(intensityAtR));0187 fprintf(<span class="string">'...real iteration: %3d\n'</span>,i);0188 ⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -