mesh_shrinkwrap_old

0001 0001 function [FV, Edges] = mesh_shrinkwrap(vol,FV,smooth,vthresh,interpVal,...
0002 0002     fitval,fittol,fititer,fitchange,fitvattr)
0003 0003 
0004 0004 % mesh_shrinkwrap - Tesselate the surface of a 3D volume
0005 0005 %
0006 0006 % [FV, Edges] = mesh_shrinkwrap(vol,FV,smooth,vthresh,interpVal,...
0007 0007 %               fitval,fittol,fititer,fitchange,fitvattr)
0008 0008 %
0009 0009 % vol       - a 3D image volume
0010 0010 % FV        - input tesselation; if empty, sphere tesselation
0011 0011 %             is created.  FV has fields FV.vertices, FV.faces
0012 0012 % smooth    - Gaussian smoothing of vol (5x5x5 kernel), 0|1 (default 0)
0013 0013 % vthresh   - binarise vol at fitval, 0|1 (default 0)
0014 0014 % interpVal - use radial interpolation (faster) or incremental
0015 0015 %             radial shrink (slower), 0|1 (default 1, faster)
0016 0016 % fitval    - image intensity to shrink wrap (default 20)
0017 0017 % fittol    - image intensity tolerance (default 5)
0018 0018 % fititer   - max number of iterations to fit (default 200)
0019 0019 % fitchange - least sig. change in intensity values
0020 0020 %             between fit iterations (default 2)
0021 0021 % fitvattr  - vertex attraction (constraint), 0:1, smaller
0022 0022 %             values are less constraint; close to 0 for
0023 0023 %             no constraint is useful when dealing with
0024 0024 %             binary volumes, otherwise 0.4 (40%) seems OK
0025 0025 %
0026 0026 % FV        - a struct with 2562 vertices and 5120 faces
0027 0027 % Edges     - a [2562,2562] matrix of edge connectivity for FV
0028 0028 %
0029 0029 % An alternative to isosurface for volumes with an external
0030 0030 % surface that can be shrink-wrapped. It has been developed to
0031 0031 % find the scalp surface for MRI of the human head.
0032 0032 %
0033 0033 % It starts with a sphere tesselation (large radius) and moves
0034 0034 % each vertex point toward the centre of the volume until it
0035 0035 % lies at or near the fitval.  The function is not optimised
0036 0036 % for speed, but it should produce reasonable results.
0037 0037 %
0038 0038 % Example of creating a scalp tesselation for SPM T1 MRI template:
0039 0039 %
0040 0040 %   avw = avw_read('T1');
0041 0041 %   FV = mesh_shrinkwrap(avw.img,[],0,0,intensity,5.0,50,0.5,0.4);
0042 0042 %   patch('vertices',FV.vertices,'faces',FV.faces,'facecolor',[.6 .6 .6]);
0043 0043 %
0044 0044 % Example of creating a skull from FSL BET skull volume:
0045 0045 %
0046 0046 %   avw = avw_read('T1_skull');
0047 0047 %   FV = mesh_shrinkwrap(avw.img,[],1,0,intensity,0.2,10,0.005,0.1);
0048 0048 %   patch('vertices',FV.vertices,'faces',FV.faces,'facecolor',[.6 .6 .6]);
0049 0049 %
0050 0050 % ***** NOTE *****
0051 0051 % An important limitation at present is that it doesn't
0052 0052 % read the image dimensions, but assumes they are 1mm^3.  Further
0053 0053 % versions might take a standard analyze volume and check the
0054 0054 % header to scale the result according to the image dimensions. At
0055 0055 % present, this must be done with the results of this function.  The
0056 0056 % output vertex coordinates are in mm with an origin at (0,0,0),
0057 0057 % which lies at the center of the MRI volume.
0058 0058 %
0059 0059 % See also: ISOSURFACE, SPHERE_TRI, MESH_REFINE_TRI4,
0060 0060 %           MESH_BEM_SHELLS_FUNC, MESH_BEM_SHELLS_SCRIPT
0061 0061 %
0062 0062 
0063 0063 % $Revision: 1.1 $ $Date: 2004/11/12 01:32:35 $
0064 0064 
0065 0065 % Licence:  GNU GPL, no implied or express warranties
0066 0066 % History:  07/2002, Darren.Weber_at_radiology.ucsf.edu
0067 0067 %                    - created to provide alternative to
0068 0068 %                      isosurface in matlab R12
0069 0069 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0070 0070 
0071 0071 version = '[$Revision: 1.1 $]';
0072 0072 fprintf('MESH_SHRINKWRAP [v%s]\n',version(12:16));  tic;
0073 0073 
0074 0074 % Parse arguments
0075 0075 
0076 0076 if ~exist('vol','var'), error('...no input volume\n');
0077 0077 elseif isempty(vol),    error('...empty input volume\n');
0078 0078 end
0079 0079 if isstruct(vol),
0080 0080     if isfield(vol,'img'), vol = vol.img;
0081 0081     else
0082 0082         error('...input volume is a struct, it must be a 3D image volume only\n');
0083 0083     end
0084 0084 end
0085 0085 
0086 0086 if ~exist('smooth','var'),   smooth = 0;
0087 0087 elseif isempty(smooth),      smooth = 0;
0088 0088 end
0089 0089 
0090 0090 if ~exist('vthresh','var'),  vthresh = 0;
0091 0091 elseif isempty(vthresh),     vthresh = 0;
0092 0092 end
0093 0093 
0094 0094 if ~exist('interpVal','var'), interpVal = 1;
0095 0095 elseif isempty(interpVal),    interpVal = 1;
0096 0096 end
0097 0097 
0098 0098 if ~exist('fitval','var'),   fit.val = 20;
0099 0099 elseif isempty(fitval),      fit.val = 20;
0100 0100 else                         fit.val = fitval;
0101 0101 end
0102 0102 
0103 0103 if ~exist('fittol','var'),   fit.tol = 5;
0104 0104 elseif isempty(fittol),      fit.tol = 5;
0105 0105 else                         fit.tol = fittol;
0106 0106 end
0107 0107 
0108 0108 if fit.val <= fit.tol,
0109 0109     error('...must use fit tolerance < fit value\n');
0110 0110 end
0111 0111 
0112 0112 if ~exist('fititer','var'),  fit.iter = 200;
0113 0113 elseif isempty(fititer),     fit.iter = 200;
0114 0114 else                         fit.iter = fititer;
0115 0115 end
0116 0116 
0117 0117 if ~exist('fitchange','var'),fit.change = 2;
0118 0118 elseif isempty(fitchange),   fit.change = 2;
0119 0119 else                         fit.change = fitchange;
0120 0120 end
0121 0121 
0122 0122 if ~exist('fitvattr','var'), fit.vattr = 0.4;
0123 0123 elseif isempty(fitvattr),    fit.vattr = 0.4;
0124 0124 else                         fit.vattr = fitvattr;
0125 0125 end
0126 0126 if fit.vattr > 1,
0127 0127     fprintf('...fit vertattr (v) must be 0 <= v <= 1, setting v = 1\n');
0128 0128     fit.vattr = 1;
0129 0129 end
0130 0130 if fit.vattr < 0,
0131 0131     fprintf('...fit vertattr (v) must be 0 <= v <= 1, setting v = 0.\n');
0132 0132     fit.vattr = 0;
0133 0133 end
0134 0134 
0135 0135 % MAIN
0136 0136 
0137 0137 % Find approximate centre of volume
0138 0138 xdim = size(vol,1);
0139 0139 ydim = size(vol,2);
0140 0140 zdim = size(vol,3);
0141 0141 
0142 0142 origin(1) = floor(xdim/2);
0143 0143 origin(2) = floor(ydim/2);
0144 0144 origin(3) = floor(zdim/2);
0145 0145 
0146 0146 % Check whether to create a sphere tesselation
0147 0147 % or use an input tesselation as the start point
0148 0148 sphere = 0;
0149 0149 if ~exist('FV','var'),
0150 0150     sphere = 1;
0151 0151 elseif ~isfield(FV,'vertices'),
0152 0152     sphere = 1;
0153 0153 elseif ~isfield(FV,'faces'),
0154 0154     sphere = 1;
0155 0155 elseif isempty(FV.vertices),
0156 0156     sphere = 1;
0157 0157 elseif isempty(FV.faces),
0158 0158     sphere = 1;
0159 0159 end
0160 0160 if sphere,
0161 0161     % Create a sphere tesselation to encompass the volume
0162 0162     radius = max([xdim ydim zdim]) / 1.5;
0163 0163     FV = sphere_tri('ico',4,radius); % 2562 vertices
0164 0164 else
0165 0165     fprintf('...using input FV tesselation...\n');
0166 0166 end
0167 0167 
0168 0168 
0169 0169 % the 'edge' matrix is the connectivity of all vertices,
0170 0170 % used to find neighbours during movement of vertices,
0171 0171 % including smoothing the tesselation
0172 0172 FV.edge = mesh_edges(FV);
0173 0173 
0174 0174 
0175 0175 % Shift the centre of the sphere to the centre of the volume
0176 0176 centre = repmat(origin,size(FV.vertices,1),1);
0177 0177 FV.vertices = FV.vertices + centre;
0178 0178 
0179 0179 
0180 0180 % 'Binarise' the volume, removing all values below
0181 0181 % a threshold, setting all others to threshold
0182 0182 if vthresh,
0183 0183     fprintf('...thresholding volume...'); tic;
0184 0184     Vindex = find(vol < fit.val);
0185 0185     vol(Vindex) = 0;
0186 0186     Vindex = find(vol >= fit.val);
0187 0187     vol(Vindex) = fit.val;
0188 0188     t = toc; fprintf('done (%5.2f sec)\n',t);
0189 0189 end
0190 0190 binvol = find(vol > 1);
0191 0191 
0192 0192 % smooth the volume
0193 0193 if smooth,
0194 0194     fprintf('...gaussian smoothing (5-10 minutes)...'); tic;
0195 0195     vol = smooth3(vol,'gaussian',5,.8);
0196 0196     t = toc; fprintf('done (%5.2f sec)\n',t);
0197 0197 end
0198 0198 
0199 0199 
0200 0200 % Now begin recursion
0201 0201 fprintf('...fitting...\n');    tic;
0202 0202 
0203 0203 i = 1;
0204 0204 Fminima = 0;
0205 0205 intensityAtDMean = [0 0];
0206 0206 
0207 0207 while i <= fit.iter,
0208 0208     
0209 0209     if interpVal,
0210 0210         % use radial interpolation method, moving directly
0211 0211         % to the intensity value nearest correct intensity
0212 0212         [FV, intensityAtD, D] = locate_val(FV,vol,origin,fit);
0213 0213     else
0214 0214         % use incremental method, moving along radial line
0215 0215         % gradually until finding correct intensity
0216 0216         [FV, intensityAtD, D] = shrink_wrap(FV,vol,origin,fit);
0217 0217     end
0218 0218     
0219 0219     intensityAtDMean(1) = intensityAtDMean(2);
0220 0220     intensityAtDMean(2) = mean(intensityAtD);
0221 0221     
0222 0222     fprintf('...distance:  mean = %8.4f mm, std = %8.4f mm\n',mean(D),std(D));
0223 0223     fprintf('...intensity: mean = %8.4f,    std = %8.4f\n',...
0224 0224         mean(intensityAtD),std(intensityAtD));
0225 0225     fprintf('...real iteration: %3d\n',i);
0226 0226     
0227 0227     % Is the mean distance reasonable?
0228 0228     if mean(D) < 0.5,
0229 0229         error('...mean distance < 0.5 mm!\n');
0230 0230     end
0231 0231     
0232 0232     % MDifVal is the mean of the absolute difference
0233 0233     % between the vertex intensity and the fit intensity
0234 0234     MDifVal = abs(intensityAtDMean(2) - fit.val);
0235 0235     
0236 0236     % Is the mean difference within the tolerance range?
0237 0237     if MDifVal < fit.tol,
0238 0238         fprintf('...mean intensity difference < tolerance (%5.2f +/- %5.2f)\n',...
0239 0239             fit.val,fit.tol);
0240 0240         break;
0241 0241     else
0242 0242         fprintf('...mean intensity difference > tolerance (%5.2f +/- %5.2f)\n',...
0243 0243             fit.val,fit.tol);
0244 0244     end
0245 0245     
0246 0246     % How much has the intensity values changed?
0247 0247     if (i > 1) & intensityAtDMean(2) > 0,
0248 0248         if intensityAtDMean(2) - intensityAtDMean(1) < fit.change,
0249 0249             fprintf('...no significant intensity change (< %5.2f) in this iteration\n',...
0250 0250                 fit.change);
0251 0251             Fminima = Fminima + 1;
0252 0252             if Fminima >= 5,
0253 0253                 fprintf('...no significant intensity change in last 5 iterations\n');
0254 0254                 break;
0255 0255             end
0256 0256         else
0257 0257             Fminima = 0;
0258 0258         end
0259 0259     end
0260 0260     
0261 0261     % Ensure that iterations begin when MDifVal is truly initialised
0262 0262     if isnan(MDifVal),
0263 0263         i = 1;
0264 0264     else,
0265 0265         i = i + 1;
0266 0266     end
0267 0267     
0268 0268 end
0269 0269 
0270 0270 FV = mesh_smooth(FV,origin,fit.vattr);
0271 0271 
0272 0272 % Remove large edges matrix from FV
0273 0273 Edges = FV.edge;
0274 0274 FV = struct('vertices',FV.vertices,'faces',FV.faces);
0275 0275 
0276 0276 % Now center the output vertices at 0,0,0 by subtracting
0277 0277 % the volume centroid
0278 0278 FV.vertices(:,1) = FV.vertices(:,1) - origin(1);
0279 0279 FV.vertices(:,2) = FV.vertices(:,2) - origin(2);
0280 0280 FV.vertices(:,3) = FV.vertices(:,3) - origin(3);
0281 0281 
0282 0282 t=toc; fprintf('...done (%5.2f sec).\n\n',t);
0283 0283 
0284 0284 return
0285 0285 
0286 0286 
0287 0287 
0288 0288 
0289 0289 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0290 0290 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0291 0291 function [FV, intensityAtD, D] = locate_val(FV,vol,origin,fit),
0292 0292     
0293 0293     xo = origin(1); yo = origin(2); zo = origin(3);
0294 0294     
0295 0295     Nvert = size(FV.vertices,1);
0296 0296     progress = round(.1 * Nvert);
0297 0297     
0298 0298     % Initialise difference intensity & distance arrays
0299 0299     intensityAtD = zeros(Nvert,1);
0300 0300     D    = intensityAtD;
0301 0301     
0302 0302     % Find distance and direction cosines for line from
0303 0303     % origin to all vertices
0304 0304     XV = FV.vertices(:,1);
0305 0305     YV = FV.vertices(:,2);
0306 0306     ZV = FV.vertices(:,3);
0307 0307     DV = sqrt( (XV-xo).^2 + (YV-yo).^2 + (ZV-zo).^2 );
0308 0308     LV = (XV-xo)./DV; % cos alpha
0309 0309     MV = (YV-yo)./DV; % cos beta
0310 0310     NV = (ZV-zo)./DV; % cos gamma
0311 0311     
0312 0312     % Check for binary volume data, if empty, binary
0313 0313     binvol = find(vol > 1);
0314 0314     
0315 0315     % Locate each vertex at a given fit value
0316 0316     tic
0317 0317     for v = 1:Nvert,
0318 0318         
0319 0319         if v > progress,
0320 0320             fprintf('...interp3 processed %4d of %4d vertices',progress,Nvert);
0321 0321             t = toc; fprintf(' (%5.2f sec)\n',t);
0322 0322             progress = progress + progress;
0323 0323         end
0324 0324         
0325 0325         % Find direction cosines for line from origin to vertex
0326 0326         x = XV(v);
0327 0327         y = YV(v);
0328 0328         z = ZV(v);
0329 0329         d = DV(v);
0330 0330         l = LV(v); % cos alpha
0331 0331         m = MV(v); % cos beta
0332 0332         n = NV(v); % cos gamma
0333 0333         
0334 0334         % find discrete points between origin
0335 0335         % and vertex + 50% of vertex distance
0336 0336         points = 250;
0337 0337         
0338 0338         Darray = linspace(0,(d + .2 * d),points);
0339 0339         
0340 0340         L = repmat(l,1,points);
0341 0341         M = repmat(m,1,points);
0342 0342         N = repmat(n,1,points);
0343 0343         
0344 0344         XI = (L .* Darray) + xo;
0345 0345         YI = (M .* Darray) + yo;
0346 0346         ZI = (N .* Darray) + zo;
0347 0347         
0348 0348         % interpolate volume values at these points
0349 0349         VI = interp3(vol,YI,XI,ZI,'*linear');
0350 0350         
0351 0351         % do we have a binary volume (no values > 1)
0352 0352         if isempty(binvol),
0353 0353             maxindex = max(find(VI>0));
0354 0354             if maxindex,
0355 0355                 D(v) = Darray(maxindex);
0356 0356             end
0357 0357         else
0358 0358             % find the finite values of VI
0359 0359             index = max(find(VI(isfinite(VI))));
0360 0360             if index,
0361 0361                 
0362 0362                 % Find nearest volume value to the required fit value
0363 0363                 nearest = max(find(VI >= fit.val));
0364 0364                 
0365 0365                 %[ nearest, value ] = NearestArrayPoint( VI, fit.val );
0366 0366                 
0367 0367                 % Check this nearest index against a differential
0368 0368                 % negative peak value
0369 0369                 %diffVI = diff(VI);
0370 0370                 %if max(VI) > 1,
0371 0371                 %    diffindex = find(diffVI < -20);
0372 0372                 %else
0373 0373                 % probably a binary volume
0374 0374                 %    diffindex = find(diffVI < 0);
0375 0375                 %end
0376 0376                 
0377 0377                 % now set d
0378 0378                 if nearest,
0379 0379                     D(v) = Darray(nearest);
0380 0380                 end
0381 0381             end
0382 0382         end
0383 0383         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0384 0384         % Constrain relocation by fit.vattr,
0385 0385         % some % of distance from neighbours
0386 0386         
0387 0387         vi = find(FV.edge(v,:));  % the neighbours' indices
0388 0388         X = FV.vertices(vi,1);    % the neighbours' vertices
0389 0389         Y = FV.vertices(vi,2);
0390 0390         Z = FV.vertices(vi,3);
0391 0391         
0392 0392         % Find neighbour distances
0393 0393         DN = sqrt( (X-xo).^2 + (Y-yo).^2 + (Z-zo).^2 );
0394 0394         % Find mean distance of neighbours
0395 0395         DNmean = mean(DN);
0396 0396         
0397 0397         minattr = fit.vattr;
0398 0398         maxattr = 1 + fit.vattr;
0399 0399         
0400 0400         if D(v) < (minattr * DNmean),
0401 0401             D(v) = minattr * DNmean;
0402 0402         end
0403 0403         if D(v) > (maxattr * DNmean),
0404 0404             D(v) = maxattr * DNmean;
0405 0405         end
0406 0406         if D(v) == 0, D(v) = DNmean; end
0407 0407         
0408 0408         % relocate vertex to new distance
0409 0409         x = (l * D(v)) + xo;
0410 0410         y = (m * D(v)) + yo;
0411 0411         z = (n * D(v)) + zo;
0412 0412         
0413 0413         FV.vertices(v,:) = [ x y z ];
0414 0414         
0415 0415         % Find intensity value at this distance
0416 0416         intensityAtD(v) = interp1(Darray,VI,D(v),'linear');
0417 0417         
0418 0418     end
0419 0419     
0420 0420     if isempty(binvol),
0421 0421         % check outliers and smooth twice for binary volumes
0422 0422         FV = vertex_outliers(FV, D, origin);
0423 0423         FV = mesh_smooth(FV,origin,fit.vattr);
0424 0424     end
0425 0425     FV = mesh_smooth(FV,origin,fit.vattr);
0426 0426     
0427 0427 return
0428 0428 
0429 0429 
0430 0430 
0431 0431 
0432 0432 
0433 0433 
0434 0434 
0435 0435 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0436 0436 function [FV, intensityAtD, D] = shrink_wrap(FV,vol,origin,fit),
0437 0437     
0438 0438     xo = origin(1); yo = origin(2); zo = origin(3);
0439 0439     
0440 0440     Nvert = size(FV.vertices,1);
0441 0441     
0442 0442     intensityAtD = zeros(Nvert,1);
0443 0443     D    = intensityAtD;
0444 0444     
0445 0445     for v = 1:Nvert,
0446 0446         
0447 0447         x = FV.vertices(v,1);
0448 0448         y = FV.vertices(v,2);
0449 0449         z = FV.vertices(v,3);
0450 0450         
0451 0451         % Find distance of vertex from origin
0452 0452         d = sqrt( (x-xo)^2 + (y-yo)^2 + (z-zo)^2 );
0453 0453         
0454 0454         % check whether vertex is already at fit.val in vol
0455 0455         
0456 0456         volval = vol_val(vol,x,y,z);
0457 0457         
0458 0458         if isnan(volval) | (volval < fit.val - fit.tol) | (volval > fit.val + fit.tol),
0459 0459             
0460 0460             % Find direction cosines for line from centre to vertex
0461 0461             l = (x-xo)/d; % cos alpha
0462 0462             m = (y-yo)/d; % cos beta
0463 0463             n = (z-zo)/d; % cos gamma
0464 0464             
0465 0465             % now modify d by fit.dist
0466 0466             if isnan(volval) | (volval < fit.val - fit.tol),
0467 0467                 d = d - (d * fit.vattr);
0468 0468             else
0469 0469                 d = d + (d * fit.vattr);
0470 0470             end
0471 0471             
0472 0472             
0473 0473             %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0474 0474             % Constrain relocation by fit.vattr,
0475 0475             % some % of distance from neighbours
0476 0476             
0477 0477             vi = find(FV.edge(v,:));  % the neighbours' indices
0478 0478             X = FV.vertices(vi,1);    % the neighbours' vertices
0479 0479             Y = FV.vertices(vi,2);
0480 0480             Z = FV.vertices(vi,3);
0481 0481             
0482 0482             % Find neighbour distances
0483 0483             DN = sqrt( (X-xo).^2 + (Y-yo).^2 + (Z-zo).^2 );
0484 0484             % Find mean distance of neighbours
0485 0485             DNmean = mean(DN);
0486 0486             
0487 0487             minattr = fit.vattr;
0488 0488             maxattr = 1 + fit.vattr;
0489 0489             
0490 0490             if d < (minattr * DNmean),
0491 0491                 d = minattr * DNmean;
0492 0492             elseif d > (maxattr * DNmean),
0493 0493                 d = maxattr * DNmean;
0494 0494             end
0495 0495             
0496 0496             % locate vertex at this new distance
0497 0497             x = (l * d) + xo;
0498 0498             y = (m * d) + yo;
0499 0499             z = (n * d) + zo;
0500 0500             
0501 0501             FV.vertices(v,:) = [ x y z ];
0502 0502         end
0503 0503         
0504 0504         intensityAtD(v) = vol_val(vol,x,y,z);
0505 0505         D(v) = d;
0506 0506     end
0507 0507     
0508 0508     FV = mesh_smooth(FV,origin,fit.vattr);
0509 0509     
0510 0510 return
0511 0511 
0512 0512 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0513 0513 % Currently not calling this function (Oct 02)
0514 0514 function [val] = vol_val(vol,x,y,z),
0515 0515     
0516 0516     % This function just ensures that xyz are
0517 0517     % actually within the volume before trying
0518 0518     % to get a volume value
0519 0519     
0520 0520     val = nan; % assume zero value
0521 0521     
0522 0522     x = round(x);
0523 0523     y = round(y);
0524 0524     z = round(z);
0525 0525     
0526 0526     if x > 0 & y > 0 & z > 0,
0527 0527         
0528 0528         % get bounds of vol
0529 0529         Xv = size(vol,1);
0530 0530         Yv = size(vol,2);
0531 0531         Zv = size(vol,3);
0532 0532         
0533 0533         if x <= Xv & y <= Yv & z <= Zv,
0534 0534             % OK return volume value at xyz
0535 0535             val = vol(x,y,z);
0536 0536         end
0537 0537     end
0538 0538     
0539 0539 return
0540 0540 
0541 0541 
0542 0542 
0543 0543 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0544 0544 function [FV] = vertex_outliers(FV, D, origin),
0545 0545     
0546 0546     xo = origin(1); yo = origin(2); zo = origin(3);
0547 0547     
0548 0548     % Screen FV for outlying vertices, using
0549 0549     % mean +/- 2 * stdev of distance from origin
0550 0550     DistMean = mean(D);
0551 0551     DistStDev = std(D);
0552 0552     DistMax = DistMean + 2 * DistStDev;
0553 0553     DistMin = DistMean - 2 * DistStDev;
0554 0554     
0555 0555     for v = 1:size(FV.vertices,1),
0556 0556         
0557 0557         if D(v) >= DistMax,
0558 0558             D(v) = DistMean;
0559 0559             relocate = 1;
0560 0560         elseif D(v) <= DistMin,
0561 0561             D(v) = DistMean;
0562 0562             relocate = 1;
0563 0563         else
0564 0564             relocate = 0;
0565 0565         end
0566 0566         
0567 0567         if relocate,
0568 0568             x = FV.vertices(v,1);
0569 0569             y = FV.vertices(v,2);
0570 0570             z = FV.vertices(v,3);
0571 0571             
0572 0572             % Find direction cosines for line from centre to vertex
0573 0573             d = sqrt( (x-xo)^2 + (y-yo)^2 + (z-zo)^2 );
0574 0574             l = (x-xo)/d; % cos alpha
0575 0575             m = (y-yo)/d; % cos beta
0576 0576             n = (z-zo)/d; % cos gamma
0577 0577             
0578 0578             % relocate vertex to this new distance
0579 0579             x = (l * D(v)) + xo;
0580 0580             y = (m * D(v)) + yo;
0581 0581             z = (n * D(v)) + zo;
0582 0582             
0583 0583             FV.vertices(v,:) = [ x y z ];
0584 0584         end
0585 0585     end
0586 0586 return