Description of avw

0001 function [FV, Edges] = avw_shrinkwrap(avw,FV,interpVal,...
0002     fitval,fittol,fititer,fitchange,fitvattr)
0003 
0004 % avw_shrinkwrap - Tesselate the surface of a 3D Analyze 7.5 avw struct
0005 %
0006 % [FV, Edges] = avw_shrinkwrap(avw,FV,interpVal,...
0007 %               fitval,fittol,fititer,fitchange,fitvattr)
0008 %
0009 % avw       - an Analyze 7.5 data struct (see avw_read)
0010 % FV        - input tesselation; if empty, sphere tesselation
0011 %             is created.  FV has fields FV.vertices, FV.faces
0012 % interpVal - use radial interpolation (faster) or incremental
0013 %             radial shrink (slower), 0|1 (default 1, faster)
0014 % fitval    - image intensity to shrink wrap (default 20)
0015 % fittol    - image intensity tolerance (default 5)
0016 % fititer   - max number of iterations to fit (default 200)
0017 % fitchange - least sig. change in intensity values
0018 %             between fit iterations (default 2)
0019 % fitvattr  - vertex attraction (constraint), 0:1, smaller
0020 %             values are less constraint; close to 0 for
0021 %             no constraint is useful when dealing with
0022 %             binary volumes, otherwise 0.4 (40%) seems OK
0023 %
0024 % FV        - a struct with 2562 vertices and 5120 faces
0025 % Edges     - a [2562,2562] matrix of edge connectivity for FV
0026 %
0027 % This function has been developed to find the scalp surface
0028 % for MRI of the human head.  It is not a sophisticated, robust
0029 % algorithm!
0030 %
0031 % It starts with a sphere tesselation (large radius) and moves
0032 % each vertex point toward the center of the volume until it
0033 % lies at or near the fitval.  The vertices are constrained to
0034 % move only along the radial projection from the origin and they
0035 % are also required to stay within a small distance of their
0036 % neighbours.  The function is not optimised for speed, but
0037 % it should produce reasonable results.
0038 %
0039 % Example of creating a scalp tesselation for SPM T1 MRI template:
0040 %
0041 %   avw = avw_read('T1');
0042 %   FV = avw_shrinkwrap(avw,[],0,0,[],intensity,5.0,50,0.5,0.4);
0043 %   patch('vertices',FV.vertices,'faces',FV.faces,'facecolor',[.6 .6 .6]);
0044 %
0045 % The output vertex coordinates are in meters with an origin at (0,0,0),
0046 % which lies at the center of the MRI volume.  This function uses the
0047 % avw.hdr.dime.pixdim values to convert from voxel coordinates into
0048 % meters.
0049 %
0050 % See also: ISOSURFACE, SPHERE_TRI, MESH_REFINE_TRI4,
0051 %           MESH_BEM_SHELLS_FUNC, MESH_BEM_SHELLS_SCRIPT
0052 %
0053 
0054 % $Revision: 1.1 $ $Date: 2004/11/12 01:30:25 $
0055 
0056 % Licence:  GNU GPL, no implied or express warranties
0057 % History:  08/2003, Darren.Weber_at_radiology.ucsf.edu
0058 %                    - adapted from mesh_shrinkwrap
0059 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0060 
0061 % Parse arguments
0062 
0063 if ~exist('avw','var'), error('...no input volume\n');
0064 elseif isempty(avw),    error('...empty input volume\n');
0065 end
0066 
0067 % Find approximate center of volume
0068 center = avw_center(avw);
0069 origin = center.abs.voxels;
0070 
0071 version = '[$Revision: 1.1 $]';
0072 fprintf('AVW_SHRINKWRAP [v%s]\n',version(12:16));  tic;
0073 
0074 if ~exist('interpVal','var'), interpVal = 0;
0075 elseif isempty(interpVal),    interpVal = 0;
0076 end
0077 
0078 if ~exist('fitval','var'),   fit.val = 20;
0079 elseif isempty(fitval),      fit.val = 20;
0080 else                         fit.val = fitval;
0081 end
0082 
0083 if ~exist('fittol','var'),   fit.tol = 5;
0084 elseif isempty(fittol),      fit.tol = 5;
0085 else                         fit.tol = fittol;
0086 end
0087 
0088 if fit.val <= fit.tol,
0089     error('...must use fit tolerance < fit value\n');
0090 end
0091 
0092 if ~exist('fititer','var'),  fit.iter = 200;
0093 elseif isempty(fititer),     fit.iter = 200;
0094 else                         fit.iter = fititer;
0095 end
0096 
0097 if ~exist('fitchange','var'),fit.change = 2;
0098 elseif isempty(fitchange),   fit.change = 2;
0099 else                         fit.change = fitchange;
0100 end
0101 
0102 if ~exist('fitvattr','var'), fit.vattr = 0.4;
0103 elseif isempty(fitvattr),    fit.vattr = 0.4;
0104 else                         fit.vattr = fitvattr;
0105 end
0106 if fit.vattr > 1,
0107     fprintf('...fit vertattr (v) must be 0 <= v <= 1, setting v = 1\n');
0108     fit.vattr = 1;
0109 end
0110 if fit.vattr < 0,
0111     fprintf('...fit vertattr (v) must be 0 <= v <= 1, setting v = 0.\n');
0112     fit.vattr = 0;
0113 end
0114 
0115 
0116 % get size of volume, in voxels
0117 [xdim,ydim,zdim] = size(avw.img);
0118 
0119 
0120 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0121 % Check whether to create a sphere tesselation
0122 % or use an input tesselation as the start point
0123 
0124 sphere = 0;
0125 if ~exist('FV','var'),
0126     sphere = 1;
0127 elseif ~isfield(FV,'vertices'),
0128     sphere = 1;
0129 elseif ~isfield(FV,'faces'),
0130     sphere = 1;
0131 elseif isempty(FV.vertices),
0132     sphere = 1;
0133 elseif isempty(FV.faces),
0134     sphere = 1;
0135 end
0136 if sphere,
0137     % Create a sphere tesselation to encompass the volume
0138     diameter = max([xdim ydim zdim]);
0139     radius = diameter / 1.5;
0140     FV = sphere_tri('ico',4,radius); % 2562 vertices
0141     % Shift the center of the sphere to the center of the volume
0142     FV.vertices = FV.vertices + repmat(origin,size(FV.vertices,1),1);
0143 else
0144     fprintf('...using input FV tesselation...\n');
0145 end
0146 
0147 % the 'edge' matrix is the connectivity of all vertices,
0148 % used to find neighbours during movement of vertices,
0149 % including smoothing the tesselation
0150 FV.edge = mesh_edges(FV);
0151 
0152 
0153 
0154 
0155 
0156 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0157 % Estimate scalp intensity
0158 fit = estimate_scalp(FV,avw.img,origin,fit);
0159 
0160 
0161 
0162 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0163 % Now begin recursion
0164 fprintf('...fitting...\n');    tic;
0165 
0166 i = 1;
0167 Fminima = 0;
0168 intensityAtRMean = [0 0];
0169 
0170 while i <= fit.iter,
0171     
0172     if interpVal,
0173         % use radial interpolation method, moving directly
0174         % to the intensity value nearest correct intensity
0175         [FV, intensityAtR, R] = locate_val(FV,avw.img,origin,fit);
0176     else
0177         % use incremental method, moving along radial line
0178         % gradually until finding correct intensity
0179         [FV, intensityAtR, R] = shrink_wrap(FV,avw.img,origin,fit);
0180     end
0181     
0182     intensityAtRMean = [ intensityAtRMean(2), mean(intensityAtR) ] ;
0183     
0184     fprintf('...distance:  mean = %8.4f voxels, std = %8.4f voxels\n',mean(R),std(R));
0185     fprintf('...intensity: mean = %8.4f,        std = %8.4f\n',...
0186         mean(intensityAtR),std(intensityAtR));
0187     fprintf('...real iteration: %3d\n',i);
0188     
0189     % Is the mean distance reasonable?
0190     if mean(R) < 0.5,
0191         error('...mean distance < 0.5 voxel!\n');
0192     end
0193     
0194     % MDifVal is the mean of the absolute difference
0195     % between the vertex intensity and the fit intensity
0196     MDifVal = abs(intensityAtRMean(2) - fit.val);
0197     
0198     % Is the mean difference within the tolerance range?
0199     if MDifVal < fit.tol,
0200         fprintf('...mean intensity difference < tolerance (%5.2f +/- %5.2f)\n',...
0201             fit.val,fit.tol);
0202         break;
0203     else
0204         fprintf('...mean intensity difference > tolerance (%5.2f +/- %5.2f)\n',...
0205             fit.val,fit.tol);
0206     end
0207     
0208     % How much has the intensity values changed?
0209     if (i > 1) & intensityAtRMean(2) > 0,
0210         if intensityAtRMean(2) - intensityAtRMean(1) < fit.change,
0211             fprintf('...no significant intensity change (< %5.2f) in this iteration\n',...
0212                 fit.change);
0213             Fminima = Fminima + 1;
0214             if Fminima >= 5,
0215                 fprintf('...no significant intensity change in last 5 iterations\n');
0216                 break;
0217             end
0218         else
0219             Fminima = 0;
0220         end
0221     end
0222     
0223     % Ensure that iterations begin when MDifVal is truly initialised
0224     if isnan(MDifVal),
0225         i = 1;
0226     else,
0227         i = i + 1;
0228     end
0229     
0230 end
0231 
0232 FV = mesh_smooth(FV,origin,fit.vattr);
0233 
0234 % Remove large edges matrix from FV
0235 Edges = FV.edge;
0236 FV = struct('vertices',FV.vertices,'faces',FV.faces);
0237 
0238 % Now center the output vertices at 0,0,0 by subtracting
0239 % the volume centroid
0240 FV.vertices(:,1) = FV.vertices(:,1) - origin(1);
0241 FV.vertices(:,2) = FV.vertices(:,2) - origin(2);
0242 FV.vertices(:,3) = FV.vertices(:,3) - origin(3);
0243 
0244 % Scale the vertices by the pixdim values, after
0245 % converting them from mm to meters
0246 xpixdim = double(avw.hdr.dime.pixdim(2)) / 1000;
0247 ypixdim = double(avw.hdr.dime.pixdim(3)) / 1000;
0248 zpixdim = double(avw.hdr.dime.pixdim(4)) / 1000;
0249 
0250 FV.vertices(:,1) = FV.vertices(:,1) .* xpixdim;
0251 FV.vertices(:,2) = FV.vertices(:,2) .* ypixdim;
0252 FV.vertices(:,3) = FV.vertices(:,3) .* zpixdim;
0253 
0254 
0255 t=toc; fprintf('...done (%5.2f sec).\n\n',t);
0256 
0257 return
0258 
0259 
0260 
0261 
0262 
0263 
0264 
0265 
0266 
0267 
0268 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0269 % Subfunctions...
0270 
0271 
0272 
0273 
0274 
0275 
0276 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0277 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0278 function [fit] = estimate_scalp(FV,vol,origin,fit),
0279 
0280 xo = origin(1); yo = origin(2); zo = origin(3);
0281 
0282 Nvert = size(FV.vertices,1);
0283 
0284 % Estimate the scalp intensity, using 10% of vertices
0285 N = round(0.05 * Nvert);
0286 scalp_intensity = zeros(N,1);
0287 
0288 fprintf('...estimating scalp intensity from %d vertices\n', N);
0289 fprintf('...starting scalp intensity = %d\n', fit.val);
0290 
0291 indices = round(rand(1,N) * Nvert);
0292 
0293 i = 0;
0294 for v = indices,
0295   
0296   x = FV.vertices(v,1);
0297   y = FV.vertices(v,2);
0298   z = FV.vertices(v,3);
0299   r = sqrt( (x-xo).^2 + (y-yo).^2 + (z-zo).^2 );
0300   l = (x-xo)/r; % cos alpha
0301   m = (y-yo)/r; % cos beta
0302   n = (z-zo)/r; % cos gamma
0303   
0304   % find discrete points from zero to the vertex
0305   POINTS = 250;
0306   radial_distances = linspace(0,r,POINTS);
0307   
0308   L = repmat(l,1,POINTS);
0309   M = repmat(m,1,POINTS);
0310   N = repmat(n,1,POINTS);
0311   
0312   XI = (L .* radial_distances) + xo;
0313   YI = (M .* radial_distances) + yo;
0314   ZI = (N .* radial_distances) + zo;
0315   
0316   % interpolate volume values at these points
0317   % ( not sure why have to swap XI,YI here )
0318   VI = interp3(vol,YI,XI,ZI,'*nearest');
0319   
0320   % find location in VI where intensity gradient is steep
0321   half_max = 0.5 * max(VI);
0322   index_maxima = find(VI > half_max);
0323   % use the largest index value to locate the maxima intensity that lie
0324   % furthest toward the outside of the head
0325   index_max = index_maxima(end);
0326   
0327   i = i + 1;
0328   scalp_intensity(i,1) = VI(index_max);
0329   
0330   plot(radial_distances,VI)
0331   
0332 end
0333 
0334 fit.val = mean(scalp_intensity);
0335 fit.tol = std(scalp_intensity);
0336 
0337 fprintf('...estimated scalp intensity = %g\n', fit.val);
0338 fprintf('...estimated tolerance intensity = %g\n', fit.tol);
0339 
0340 return
0341 
0342 
0343 
0344 
0345 
0346 
0347 
0348 
0349 
0350 
0351 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0352 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0353 function [FV, intensityAtR, R] = locate_val(FV,vol,origin,fit),
0354 
0355 xo = origin(1); yo = origin(2); zo = origin(3);
0356 
0357 Nvert = size(FV.vertices,1);
0358 progress = round(.1 * Nvert);
0359 
0360 % Initialise difference intensity & distance arrays
0361 intensityAtR = zeros(Nvert,1);
0362 R = intensityAtR;
0363 
0364 % Find distance and direction cosines for line from
0365 % origin to all vertices
0366 XV = FV.vertices(:,1);
0367 YV = FV.vertices(:,2);
0368 ZV = FV.vertices(:,3);
0369 RV = sqrt( (XV-xo).^2 + (YV-yo).^2 + (ZV-zo).^2 );
0370 LV = (XV-xo)./RV; % cos alpha
0371 MV = (YV-yo)./RV; % cos beta
0372 NV = (ZV-zo)./RV; % cos gamma
0373 
0374 % Check for binary volume data, if empty, binary
0375 binvol = find(vol > 1);
0376 
0377 % Locate each vertex at a given fit value
0378 tic
0379 for v = 1:Nvert,
0380     
0381     if v > progress,
0382         fprintf('...interp3 processed %4d of %4d vertices',progress,Nvert);
0383         t = toc; fprintf(' (%5.2f sec)\n',t);
0384         progress = progress + progress;
0385     end
0386     
0387     % Find direction cosines for line from origin to vertex
0388     x = XV(v);
0389     y = YV(v);
0390     z = ZV(v);
0391     d = RV(v);
0392     l = LV(v); % cos alpha
0393     m = MV(v); % cos beta
0394     n = NV(v); % cos gamma
0395     
0396     % find discrete points between origin
0397     % and vertex + 20% of vertex distance
0398     POINTS = 250;
0399     
0400     Rarray = linspace(0,(d + .2 * d),POINTS);
0401     
0402     L = repmat(l,1,POINTS);
0403     M = repmat(m,1,POINTS);
0404     N = repmat(n,1,POINTS);
0405     
0406     XI = (L .* Rarray) + xo;
0407     YI = (M .* Rarray) + yo;
0408     ZI = (N .* Rarray) + zo;
0409     
0410     % interpolate volume values at these points
0411     % ( not sure why have to swap XI,YI here )
0412     VI = interp3(vol,YI,XI,ZI,'*linear');
0413     
0414     % do we have a binary volume (no values > 1)
0415     if isempty(binvol),
0416         maxindex = max(find(VI>0));
0417         if maxindex,
0418             R(v) = Rarray(maxindex);
0419         end
0420     else
0421         % find the finite values of VI
0422         index = max(find(VI(isfinite(VI))));
0423         if index,
0424             
0425             % Find nearest volume value to the required fit value
0426             nearest = max(find(VI >= fit.val));
0427             
0428             %[ nearest, value ] = NearestArrayPoint( VI, fit.val );
0429             
0430             % Check this nearest index against a differential
0431             % negative peak value
0432             %diffVI = diff(VI);
0433             %if max(VI) > 1,
0434             %    diffindex = find(diffVI < -20);
0435             %else
0436             % probably a binary volume
0437             %    diffindex = find(diffVI < 0);
0438             %end
0439             
0440             % now set d
0441             if nearest,
0442                 R(v) = Rarray(nearest);
0443             end
0444         end
0445     end
0446     
0447     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0448     % Constrain relocation by fit.vattr,
0449     % some % of distance from neighbours
0450     
0451     vi = find(FV.edge(v,:));  % the neighbours' indices
0452     X = FV.vertices(vi,1);    % the neighbours' vertices
0453     Y = FV.vertices(vi,2);
0454     Z = FV.vertices(vi,3);
0455     
0456     % Find neighbour distances
0457     RN = sqrt( (X-xo).^2 + (Y-yo).^2 + (Z-zo).^2 );
0458     % Find mean distance of neighbours
0459     neighbour_distances_mean = mean(RN);
0460     
0461     minattr = fit.vattr;
0462     maxattr = 1 + fit.vattr;
0463     
0464     if R(v) < (minattr * neighbour_distances_mean),
0465         R(v) = minattr * neighbour_distances_mean;
0466     end
0467     if R(v) > (maxattr * neighbour_distances_mean),
0468         R(v) = maxattr * neighbour_distances_mean;
0469     end
0470     if R(v) == 0, R(v) = neighbour_distances_mean; end
0471     
0472     % relocate vertex to new distance
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     % Find intensity value at this distance
0480     intensityAtR(v) = interp1(Rarray,VI,R(v),'linear');
0481     
0482 end
0483 
0484 if isempty(binvol),
0485     % check outliers and smooth twice for binary volumes
0486     FV = vertex_outliers(FV, R, origin);
0487     FV = mesh_smooth(FV,origin,fit.vattr);
0488 end
0489 FV = mesh_smooth(FV,origin,fit.vattr);
0490 
0491 return
0492 
0493 
0494 
0495 
0496 
0497 
0498 
0499 
0500 
0501 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0502 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0503 function [FV, intensityAtR, R] = shrink_wrap(FV,vol,origin,fit),
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 % Initialise difference intensity & distance arrays
0511 R = zeros(Nvert,1);
0512 intensityAtR = R;
0513 
0514 % Check for binary volume data, if empty, binary
0515 binvol = find(vol > 1);
0516 
0517 Imin = 0;
0518 Imax = max(max(max(vol)));
0519 
0520 % Manipulate each vertex
0521 tic
0522 for v = 1:Nvert,
0523   
0524   if v > progress,
0525     fprintf('...shrink-wrap processed %4d of %4d vertices',progress,Nvert);
0526     t = toc; fprintf(' (%5.2f sec)\n',t);
0527     progress = progress + progress;
0528   end
0529   
0530   % Find direction cosines for line from origin to vertex
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; % cos alpha
0536   m = (y-yo)/r; % cos beta
0537   n = (z-zo)/r; % cos gamma
0538   
0539   % interpolate volume values at this point ( x,y are swapped here
0540   % because the Analyze volume is a left handed coordinate system )
0541   intensity_old = interp3(vol,y,x,z,'*nearest');
0542   
0543   % move vertex closer to the origin, in a radial direction
0544   r_change = r * 0.05;
0545   r_new = r - r_change;
0546   
0547   % Calculate new vertex coordinates
0548   x = (l * r_new) + xo; % cos alpha
0549   y = (m * r_new) + yo; % cos beta
0550   z = (n * r_new) + zo; % cos gamma
0551   
0552   % interpolate volume values at this point ( x,y are swapped here
0553   % because the Analyze volume is a left handed coordinate system )
0554   intensity_new = interp3(vol,y,x,z,'*nearest');
0555   
0556   intensity_dif = intensity_new - intensity_old;
0557   
0558   if intensity_dif == 0,
0559     % relocate vertex to new distance
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   elseif (fit.val - intensity_new) < (fit.val - intensity_old),
0566     % relocate vertex to new distance
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   else
0573     intensityAtR(v,1) = intensity_old;
0574     R(v) = r;
0575   end
0576   
0577   FV = constrain_vertex(FV,v,origin);
0578   
0579 end
0580 
0581 FV = mesh_smooth(FV,origin,fit.vattr);
0582 
0583 return
0584 
0585 
0586 
0587 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0588 function [FV] = constrain_vertex(FV,index,origin),
0589 
0590 % This function adapts Smith, S. (2002), Fast robust automated brain
0591 % extraction.  Human Brain Mapping, 17: 143-155.  It corresponds to update
0592 % component 2: surface smoothness control.  It assumes that vertices are
0593 % moved along a radial line from an origin, given by direction
0594 % cosines, rather than calculating the surface normal vector.
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 % Find radial distance of vertex from origin
0604 r = sqrt( (x-xo)^2 + (y-yo)^2 + (z-zo)^2 );
0605 
0606 % Calculate unit vector
0607 v_unit_vector = ( v - origin ) / r;
0608 
0609 % Find direction cosines for line from center to vertex
0610 l = (x-xo)/r; % cos alpha
0611 m = (y-yo)/r; % cos beta
0612 n = (z-zo)/r; % cos gamma
0613 
0614 % Find neighbouring vertex coordinates
0615 vi = find(FV.edge(index,:));  % the indices
0616 neighbour_vertices = FV.vertices(vi,:);
0617 X = neighbour_vertices(:,1);
0618 Y = neighbour_vertices(:,2);
0619 Z = neighbour_vertices(:,3);
0620 
0621 % Find neighbour radial distances
0622 r_neighbours = sqrt( (X-xo).^2 + (Y-yo).^2 + (Z-zo).^2 );
0623 r_neighbours_mean = mean(r_neighbours);
0624 
0625 % Find difference in radial distance between the vertex of interest and its
0626 % neighbours; this value approximates the magnitude of sn in
0627 % Smith (2002, eq. 1 to 4)
0628 r_diff = r - r_neighbours_mean;
0629 
0630 % Find the vector sn, in the direction of the vertex of interest, given the
0631 % difference in radial distance between vertex and mean of neighbours
0632 sn = r_diff * v_unit_vector;
0633 
0634 % Find distances between vertex and neighbours, using edge lengths.
0635 % The mean value is l in Smith (2002, eq. 4)
0636 edge_distance = FV.edge(index,vi);
0637 edge_distance_mean = mean(edge_distance);
0638 
0639 % Calculate radius of local curvature, solve Smith (2002, eq. 4)
0640 if r_diff,
0641   radius_of_curvature = (edge_distance_mean ^ 2) / (2 * r_diff);
0642 else
0643   radius_of_curvature = 10000;
0644 end
0645 
0646 % Define limits for radius of curvature
0647 radius_min =  3.33; % mm
0648 radius_max = 10.00; % mm
0649 
0650 % Sigmoid function parameters,
0651 % "where E and F control the scale and offset of the sigmoid"
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 % multiply sigmoid function by sn
0658 move_vector = Fsigmoid * sn;
0659 
0660 FV.vertices(index,:) = v + move_vector;
0661 
0662 return
0663 
0664 
0665 
0666 
0667 
0668 
0669 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0670 function [val] = vol_val(vol,x,y,z),
0671 
0672 % This function just ensures that xyz are
0673 % actually within the volume before trying
0674 % to get a volume value
0675 
0676 val = nan; % assume zero value
0677 
0678 x = round(x);
0679 y = round(y);
0680 z = round(z);
0681 
0682 if x > 0 & y > 0 & z > 0,
0683     
0684     % get bounds of vol
0685     xdim = size(vol,1);
0686     ydim = size(vol,2);
0687     zdim = size(vol,3);
0688     
0689     if x <= xdim & y <= ydim & z <= zdim,
0690         % OK return volume value at xyz
0691         val = vol(x,y,z);
0692     end
0693 end
0694 
0695 return
0696 
0697 
0698 
0699 
0700 
0701 
0702 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0703 function [FV] = vertex_outliers(FV, R, origin),
0704 
0705 xo = origin(1); yo = origin(2); zo = origin(3);
0706 
0707 % Screen FV for outlying vertices, using
0708 % mean +/- 2 * stdev of distance from origin
0709 DistMean = mean(R);
0710 DistStDev = std(R);
0711 DistMax = DistMean + 2 * DistStDev;
0712 DistMin = DistMean - 2 * DistStDev;
0713 
0714 for v = 1:size(FV.vertices,1),
0715     
0716     if R(v) >= DistMax,
0717         R(v) = DistMean;
0718         relocate = 1;
0719     elseif R(v) <= DistMin,
0720         R(v) = DistMean;
0721         relocate = 1;
0722     else
0723         relocate = 0;
0724     end
0725     
0726     if relocate,
0727         x = FV.vertices(v,1);
0728         y = FV.vertices(v,2);
0729         z = FV.vertices(v,3);
0730         
0731         % Find direction cosines for line from center to vertex
0732         d = sqrt( (x-xo)^2 + (y-yo)^2 + (z-zo)^2 );
0733         l = (x-xo)/d; % cos alpha
0734         m = (y-yo)/d; % cos beta
0735         n = (z-zo)/d; % cos gamma
0736         
0737         % relocate vertex to this new distance
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     end
0744 end
0745 return
avw_shrinkwrap

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
