avw_shrinkwrap.html

mri_toolbox是一个工具用来MRI. 来自于SourceForge, 我上传这个软件,希望能结识对医疗软件感兴趣的兄弟.

0469     <span class="keyword">end</span>
0470     <span class="keyword">if</span> R(v) == 0, R(v) = neighbour_distances_mean; <span class="keyword">end</span>
0471     
0472     <span class="comment">% relocate vertex to new distance</span>
0473     x = (l * R(v)) + xo;
0474     y = (m * R(v)) + yo;
0475     z = (n * R(v)) + zo;
0476     
0477     FV.vertices(v,:) = [ x y z ];
0478     
0479     <span class="comment">% Find intensity value at this distance</span>
0480     intensityAtR(v) = interp1(Rarray,VI,R(v),<span class="string">'linear'</span>);
0481     
0482 <span class="keyword">end</span>
0483 
0484 <span class="keyword">if</span> isempty(binvol),
0485     <span class="comment">% check outliers and smooth twice for binary volumes</span>
0486     FV = <a href="#_sub6" class="code" title="subfunction [FV] = vertex_outliers(FV, R, origin),">vertex_outliers</a>(FV, R, origin);
0487     FV = mesh_smooth(FV,origin,fit.vattr);
0488 <span class="keyword">end</span>
0489 FV = mesh_smooth(FV,origin,fit.vattr);
0490 
0491 <span class="keyword">return</span>
0492 
0493 
0494 
0495 
0496 
0497 
0498 
0499 
0500 
0501 <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0502 <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0503 <a name="_sub3" href="#_subfunctions" class="code">function [FV, intensityAtR, R] = shrink_wrap(FV,vol,origin,fit),</a>
0504 
0505 xo = origin(1); yo = origin(2); zo = origin(3);
0506 
0507 Nvert = size(FV.vertices,1);
0508 progress = round(.1 * Nvert);
0509 
0510 <span class="comment">% Initialise difference intensity &amp; distance arrays</span>
0511 R = zeros(Nvert,1);
0512 intensityAtR = R;
0513 
0514 <span class="comment">% Check for binary volume data, if empty, binary</span>
0515 binvol = find(vol &gt; 1);
0516 
0517 Imin = 0;
0518 Imax = max(max(max(vol)));
0519 
0520 <span class="comment">% Manipulate each vertex</span>
0521 tic
0522 <span class="keyword">for</span> v = 1:Nvert,
0523   
0524   <span class="keyword">if</span> v &gt; progress,
0525     fprintf(<span class="string">'...shrink-wrap processed %4d of %4d vertices'</span>,progress,Nvert);
0526     t = toc; fprintf(<span class="string">' (%5.2f sec)\n'</span>,t);
0527     progress = progress + progress;
0528   <span class="keyword">end</span>
0529   
0530   <span class="comment">% Find direction cosines for line from origin to vertex</span>
0531   x = FV.vertices(v,1);
0532   y = FV.vertices(v,2);
0533   z = FV.vertices(v,3);
0534   r = sqrt( (x-xo).^2 + (y-yo).^2 + (z-zo).^2 );
0535   l = (x-xo)/r; <span class="comment">% cos alpha</span>
0536   m = (y-yo)/r; <span class="comment">% cos beta</span>
0537   n = (z-zo)/r; <span class="comment">% cos gamma</span>
0538   
0539   <span class="comment">% interpolate volume values at this point ( x,y are swapped here</span>
0540   <span class="comment">% because the Analyze volume is a left handed coordinate system )</span>
0541   intensity_old = interp3(vol,y,x,z,<span class="string">'*nearest'</span>);
0542   
0543   <span class="comment">% move vertex closer to the origin, in a radial direction</span>
0544   r_change = r * 0.05;
0545   r_new = r - r_change;
0546   
0547   <span class="comment">% Calculate new vertex coordinates</span>
0548   x = (l * r_new) + xo; <span class="comment">% cos alpha</span>
0549   y = (m * r_new) + yo; <span class="comment">% cos beta</span>
0550   z = (n * r_new) + zo; <span class="comment">% cos gamma</span>
0551   
0552   <span class="comment">% interpolate volume values at this point ( x,y are swapped here</span>
0553   <span class="comment">% because the Analyze volume is a left handed coordinate system )</span>
0554   intensity_new = interp3(vol,y,x,z,<span class="string">'*nearest'</span>);
0555   
0556   intensity_dif = intensity_new - intensity_old;
0557   
0558   <span class="keyword">if</span> intensity_dif == 0,
0559     <span class="comment">% relocate vertex to new distance</span>
0560     FV.vertices(v,1) = x;
0561     FV.vertices(v,2) = y;
0562     FV.vertices(v,3) = z;
0563     intensityAtR(v,1) = intensity_new;
0564     R(v) = r_new;
0565   <span class="keyword">elseif</span> (fit.val - intensity_new) &lt; (fit.val - intensity_old),
0566     <span class="comment">% relocate vertex to new distance</span>
0567     FV.vertices(v,1) = x;
0568     FV.vertices(v,2) = y;
0569     FV.vertices(v,3) = z;
0570     intensityAtR(v,1) = intensity_new;
0571     R(v) = r_new;
0572   <span class="keyword">else</span>
0573     intensityAtR(v,1) = intensity_old;
0574     R(v) = r;
0575   <span class="keyword">end</span>
0576   
0577   FV = <a href="#_sub4" class="code" title="subfunction [FV] = constrain_vertex(FV,index,origin),">constrain_vertex</a>(FV,v,origin);
0578   
0579 <span class="keyword">end</span>
0580 
0581 FV = mesh_smooth(FV,origin,fit.vattr);
0582 
0583 <span class="keyword">return</span>
0584 
0585 
0586 
0587 <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0588 <a name="_sub4" href="#_subfunctions" class="code">function [FV] = constrain_vertex(FV,index,origin),</a>
0589 
0590 <span class="comment">% This function adapts Smith, S. (2002), Fast robust automated brain</span>
0591 <span class="comment">% extraction.  Human Brain Mapping, 17: 143-155.  It corresponds to update</span>
0592 <span class="comment">% component 2: surface smoothness control.  It assumes that vertices are</span>
0593 <span class="comment">% moved along a radial line from an origin, given by direction</span>
0594 <span class="comment">% cosines, rather than calculating the surface normal vector.</span>
0595 
0596 xo = origin(1); yo = origin(2); zo = origin(3);
0597 
0598 v = FV.vertices(index,:);
0599 x = FV.vertices(index,1);
0600 y = FV.vertices(index,2);
0601 z = FV.vertices(index,3);
0602 
0603 <span class="comment">% Find radial distance of vertex from origin</span>
0604 r = sqrt( (x-xo)^2 + (y-yo)^2 + (z-zo)^2 );
0605 
0606 <span class="comment">% Calculate unit vector</span>
0607 v_unit_vector = ( v - origin ) / r;
0608 
0609 <span class="comment">% Find direction cosines for line from center to vertex</span>
0610 l = (x-xo)/r; <span class="comment">% cos alpha</span>
0611 m = (y-yo)/r; <span class="comment">% cos beta</span>
0612 n = (z-zo)/r; <span class="comment">% cos gamma</span>
0613 
0614 <span class="comment">% Find neighbouring vertex coordinates</span>
0615 vi = find(FV.edge(index,:));  <span class="comment">% the indices</span>
0616 neighbour_vertices = FV.vertices(vi,:);
0617 X = neighbour_vertices(:,1);
0618 Y = neighbour_vertices(:,2);
0619 Z = neighbour_vertices(:,3);
0620 
0621 <span class="comment">% Find neighbour radial distances</span>
0622 r_neighbours = sqrt( (X-xo).^2 + (Y-yo).^2 + (Z-zo).^2 );
0623 r_neighbours_mean = mean(r_neighbours);
0624 
0625 <span class="comment">% Find difference in radial distance between the vertex of interest and its</span>
0626 <span class="comment">% neighbours; this value approximates the magnitude of sn in</span>
0627 <span class="comment">% Smith (2002, eq. 1 to 4)</span>
0628 r_diff = r - r_neighbours_mean;
0629 
0630 <span class="comment">% Find the vector sn, in the direction of the vertex of interest, given the</span>
0631 <span class="comment">% difference in radial distance between vertex and mean of neighbours</span>
0632 sn = r_diff * v_unit_vector;
0633 
0634 <span class="comment">% Find distances between vertex and neighbours, using edge lengths.</span>
0635 <span class="comment">% The mean value is l in Smith (2002, eq. 4)</span>
0636 edge_distance = FV.edge(index,vi);
0637 edge_distance_mean = mean(edge_distance);
0638 
0639 <span class="comment">% Calculate radius of local curvature, solve Smith (2002, eq. 4)</span>
0640 <span class="keyword">if</span> r_diff,
0641   radius_of_curvature = (edge_distance_mean ^ 2) / (2 * r_diff);
0642 <span class="keyword">else</span>
0643   radius_of_curvature = 10000;
0644 <span class="keyword">end</span>
0645 
0646 <span class="comment">% Define limits for radius of curvature</span>
0647 radius_min =  3.33; <span class="comment">% mm</span>
0648 radius_max = 10.00; <span class="comment">% mm</span>
0649 
0650 <span class="comment">% Sigmoid function parameters,</span>
0651 <span class="comment">% &quot;where E and F control the scale and offset of the sigmoid&quot;</span>
0652 E = mean([(1 / radius_min),  (1 / radius_max)]);
0653 F = 6 * ( (1 / radius_min) - (1 / radius_max) );
0654 
0655 Fsigmoid = (1 + tanh( F * (1 / radius_of_curvature - E))) / 2;
0656 
0657 <span class="comment">% multiply sigmoid function by sn</span>
0658 move_vector = Fsigmoid * sn;
0659 
0660 FV.vertices(index,:) = v + move_vector;
0661 
0662 <span class="keyword">return</span>
0663 
0664 
0665 
0666 
0667 
0668 
0669 <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0670 <a name="_sub5" href="#_subfunctions" class="code">function [val] = vol_val(vol,x,y,z),</a>
0671 
0672 <span class="comment">% This function just ensures that xyz are</span>
0673 <span class="comment">% actually within the volume before trying</span>
0674 <span class="comment">% to get a volume value</span>
0675 
0676 val = nan; <span class="comment">% assume zero value</span>
0677 
0678 x = round(x);
0679 y = round(y);
0680 z = round(z);
0681 
0682 <span class="keyword">if</span> x &gt; 0 &amp; y &gt; 0 &amp; z &gt; 0,
0683     
0684     <span class="comment">% get bounds of vol</span>
0685     xdim = size(vol,1);
0686     ydim = size(vol,2);
0687     zdim = size(vol,3);
0688     
0689     <span class="keyword">if</span> x &lt;= xdim &amp; y &lt;= ydim &amp; z &lt;= zdim,
0690         <span class="comment">% OK return volume value at xyz</span>
0691         val = vol(x,y,z);
0692     <span class="keyword">end</span>
0693 <span class="keyword">end</span>
0694 
0695 <span class="keyword">return</span>
0696 
0697 
0698 
0699 
0700 
0701 
0702 <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0703 <a name="_sub6" href="#_subfunctions" class="code">function [FV] = vertex_outliers(FV, R, origin),</a>
0704 
0705 xo = origin(1); yo = origin(2); zo = origin(3);
0706 
0707 <span class="comment">% Screen FV for outlying vertices, using</span>
0708 <span class="comment">% mean +/- 2 * stdev of distance from origin</span>
0709 DistMean = mean(R);
0710 DistStDev = std(R);
0711 DistMax = DistMean + 2 * DistStDev;
0712 DistMin = DistMean - 2 * DistStDev;
0713 
0714 <span class="keyword">for</span> v = 1:size(FV.vertices,1),
0715     
0716     <span class="keyword">if</span> R(v) &gt;= DistMax,
0717         R(v) = DistMean;
0718         relocate = 1;
0719     <span class="keyword">elseif</span> R(v) &lt;= DistMin,
0720         R(v) = DistMean;
0721         relocate = 1;
0722     <span class="keyword">else</span>
0723         relocate = 0;
0724     <span class="keyword">end</span>
0725     
0726     <span class="keyword">if</span> relocate,
0727         x = FV.vertices(v,1);
0728         y = FV.vertices(v,2);
0729         z = FV.vertices(v,3);
0730         
0731         <span class="comment">% Find direction cosines for line from center to vertex</span>
0732         d = sqrt( (x-xo)^2 + (y-yo)^2 + (z-zo)^2 );
0733         l = (x-xo)/d; <span class="comment">% cos alpha</span>
0734         m = (y-yo)/d; <span class="comment">% cos beta</span>
0735         n = (z-zo)/d; <span class="comment">% cos gamma</span>
0736         
0737         <span class="comment">% relocate vertex to this new distance</span>
0738         x = (l * R(v)) + xo;
0739         y = (m * R(v)) + yo;
0740         z = (n * R(v)) + zo;
0741         
0742         FV.vertices(v,:) = [ x y z ];
0743     <span class="keyword">end</span>
0744 <span class="keyword">end</span>
0745 <span class="keyword">return</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> &copy; 2003</address>
</body>
</html>