TESSELLATION_STATS - Calculate statistics of the tesselation and hunt for suspicious faces and vertices
function [FaceNormal, FaceArea, FaceCenter, VertexNormal, VertexArea, SuspectFace, NumFacesEachVertex] = tessellation_stats(FV,VERBOSE);
Optional slightly faster calls:
Return the statistics of the planar triangles:
function [FaceNormal, FaceArea, FaceCenter] = tessellation_stats(FV,VERBOSE);
Return also the statistics of the vertices
function [FaceNormal, FaceArea, FaceCenter, VertexNormal, VertexArea] = tessellation_stats(FV,VERBOSE);
For INPUT of FV.vertices of size numVert x 3 and FV.faces of size numTri x 3,
then OUTPUTs are:
FaceNormal is 3 x numTri, the normal of each face
FaceArea is 1 x numTri, the area of each face
FaceCenter is 3 x numTri, the location of the center of each face
Additional calculation of the vertex statistics:
VertexNormal is 3 x numVert, the average normal assigned to each vertex,
see description below.
VertexArea is 1 x numVert, the average area assigned to each vertex, see
description below.
Additional calculation of properly ordered and arranged triangles
SuspectFace is 1 x numTri, each element giving the number of times the triangle
was identified as suspicious, see calculation below.
NumFacesEachVertex is 1 x numVert, each element giving the number of faces
attached to a vertex. Vertices with 0 faces are unassigned, 1 and 2 faces
are most likely at edges.
VERTEX NORMALS RELATIVE TO MATLAB'S
If the same FV is fed into Matlab's "h = patch(FV);" function, then vnorm =
get(h,'vertexnormals'); returns exactly the same as VertexNormal calculated
here, EXCEPT: vnorm is the reverse direction (CW ordering is positive in
Matlab), and the norm of vnorm is twice that of VertexNormal (inconsequential
scaling difference).
SURFACE TEST
This routine also performs a basic test of the ordering of the
vertices. Each face should be entered in the same CCW or CW direction. If a
triangle is ordered differently from it's neighbors in a smooth region, then
its normal vector will point in the opposite direction. A list of possible
problem faces will be displayed, where the figure number is the face number
in question. In very irregular surfaces (such as the cortex), the problem
face may be simply tucked in too tightly to the local surface for the
algorithm to catch.
SuspectFace is 1 x numTri, gives the number of tests that detected a possibly
bad face ordered in a direction opposite of the adjacent faces. A "patch" is
formed at a vertex by finding all faces that are attached to a given vertex.
The unweighted average of the face normals is computed, then dotted with
each of the face normals. A face is marked as suspect if it's normal points
in the opposite direction of the unweighted average normal. The process is
repeated for all vertices, so that each triangular face is included in three
tests. The SuspectFace is the cumulative detection score, such that
SuspectFace(i) == 3 indicates that all three vertices marked the ith face as
suspect. Detections of 1 and 2 indicate that only one or two of the vertices
marked the face as suspect.
If VERBOSE is true, then this routine will offer to display all triangular
faces that had SuspectFace == 3, i.e. all three vertices triggered the
detection.
NumFacesEachVertex, since we are generally expecting closed surfaces in brain
analysis, representing cortical surfaces or boundary elements, then vertices
with 0, 1, or 2 faces only are worthy of further attention.
For constant approximations across the triangle, then use FaceNormal and FaceCenter.
For linear approximations, then use VertexNormals and vertices.
VertexArea is the effective area spanned by a vertex, analogous to the FaceArea
spanned by triangle.
CALCULATION of VERTEX NORMALS and AREAS
Vertices in a tesselation are technically point discontinuities in the surface
description, Edges are line discontinuities. In a closed surface, then there
are precisely 2(N-2) faces for N vertices, so there are nearly half as many
vertex parameters as triangle centers representing the same surface. Hence
vertex representations of linear variation across the triangles are a popular
alternative to triangle center representation of a constant variation.
The average area assigned to each vertex is found by dividing the area of a
face by three and assigning this equally to all three vertices. Each vertex
therefore has as its area the sum of 1/3 of the areas of the attached to the
vertex.
The average normal is found by weighting the unit length normals of each face
attached to a vertex by the area of that face, then averaging. The length of
the average normal reflects this averaging.
Simple Example of 2-d triangles to check units
(1,5)
(1,0) ----------- (1,1)
| |
| |
(0,0) ------------(1,0)
number vertices as
5
4 ----------- 3
1 ----------- 2
FV.vertices = [0 0 0; 1 0 0; 1 1 0; 0 1 0; 1 5 0];% in the plane of the screen
FV.faces = [1 2 3; 1 3 4; 3 5 4]; % outward from screen0001 function [FaceNormal, FaceArea, FaceCenter, VertexNormal, VertexArea, SuspectFace, NumFacesEachVertex] = tessellation_stats(FV,VERBOSE); 0002 %TESSELLATION_STATS - Calculate statistics of the tesselation and hunt for suspicious faces and vertices 0003 % function [FaceNormal, FaceArea, FaceCenter, VertexNormal, VertexArea, SuspectFace, NumFacesEachVertex] = tessellation_stats(FV,VERBOSE); 0004 % 0005 % Optional slightly faster calls: 0006 % 0007 % Return the statistics of the planar triangles: 0008 % function [FaceNormal, FaceArea, FaceCenter] = tessellation_stats(FV,VERBOSE); 0009 % 0010 % Return also the statistics of the vertices 0011 % function [FaceNormal, FaceArea, FaceCenter, VertexNormal, VertexArea] = tessellation_stats(FV,VERBOSE); 0012 % 0013 % For INPUT of FV.vertices of size numVert x 3 and FV.faces of size numTri x 3, 0014 % then OUTPUTs are: 0015 % 0016 % FaceNormal is 3 x numTri, the normal of each face 0017 % FaceArea is 1 x numTri, the area of each face 0018 % FaceCenter is 3 x numTri, the location of the center of each face 0019 % 0020 % Additional calculation of the vertex statistics: 0021 % 0022 % VertexNormal is 3 x numVert, the average normal assigned to each vertex, 0023 % see description below. 0024 % VertexArea is 1 x numVert, the average area assigned to each vertex, see 0025 % description below. 0026 % 0027 % Additional calculation of properly ordered and arranged triangles 0028 % 0029 % SuspectFace is 1 x numTri, each element giving the number of times the triangle 0030 % was identified as suspicious, see calculation below. 0031 % NumFacesEachVertex is 1 x numVert, each element giving the number of faces 0032 % attached to a vertex. Vertices with 0 faces are unassigned, 1 and 2 faces 0033 % are most likely at edges. 0034 % 0035 % VERTEX NORMALS RELATIVE TO MATLAB'S 0036 % 0037 % If the same FV is fed into Matlab's "h = patch(FV);" function, then vnorm = 0038 % get(h,'vertexnormals'); returns exactly the same as VertexNormal calculated 0039 % here, EXCEPT: vnorm is the reverse direction (CW ordering is positive in 0040 % Matlab), and the norm of vnorm is twice that of VertexNormal (inconsequential 0041 % scaling difference). 0042 % 0043 % SURFACE TEST 0044 % 0045 % This routine also performs a basic test of the ordering of the 0046 % vertices. Each face should be entered in the same CCW or CW direction. If a 0047 % triangle is ordered differently from it's neighbors in a smooth region, then 0048 % its normal vector will point in the opposite direction. A list of possible 0049 % problem faces will be displayed, where the figure number is the face number 0050 % in question. In very irregular surfaces (such as the cortex), the problem 0051 % face may be simply tucked in too tightly to the local surface for the 0052 % algorithm to catch. 0053 % 0054 % SuspectFace is 1 x numTri, gives the number of tests that detected a possibly 0055 % bad face ordered in a direction opposite of the adjacent faces. A "patch" is 0056 % formed at a vertex by finding all faces that are attached to a given vertex. 0057 % The unweighted average of the face normals is computed, then dotted with 0058 % each of the face normals. A face is marked as suspect if it's normal points 0059 % in the opposite direction of the unweighted average normal. The process is 0060 % repeated for all vertices, so that each triangular face is included in three 0061 % tests. The SuspectFace is the cumulative detection score, such that 0062 % SuspectFace(i) == 3 indicates that all three vertices marked the ith face as 0063 % suspect. Detections of 1 and 2 indicate that only one or two of the vertices 0064 % marked the face as suspect. 0065 % 0066 % If VERBOSE is true, then this routine will offer to display all triangular 0067 % faces that had SuspectFace == 3, i.e. all three vertices triggered the 0068 % detection. 0069 % 0070 % NumFacesEachVertex, since we are generally expecting closed surfaces in brain 0071 % analysis, representing cortical surfaces or boundary elements, then vertices 0072 % with 0, 1, or 2 faces only are worthy of further attention. 0073 % 0074 % For constant approximations across the triangle, then use FaceNormal and FaceCenter. 0075 % For linear approximations, then use VertexNormals and vertices. 0076 % VertexArea is the effective area spanned by a vertex, analogous to the FaceArea 0077 % spanned by triangle. 0078 % 0079 % CALCULATION of VERTEX NORMALS and AREAS 0080 % 0081 % Vertices in a tesselation are technically point discontinuities in the surface 0082 % description, Edges are line discontinuities. In a closed surface, then there 0083 % are precisely 2(N-2) faces for N vertices, so there are nearly half as many 0084 % vertex parameters as triangle centers representing the same surface. Hence 0085 % vertex representations of linear variation across the triangles are a popular 0086 % alternative to triangle center representation of a constant variation. 0087 % 0088 % The average area assigned to each vertex is found by dividing the area of a 0089 % face by three and assigning this equally to all three vertices. Each vertex 0090 % therefore has as its area the sum of 1/3 of the areas of the attached to the 0091 % vertex. 0092 % The average normal is found by weighting the unit length normals of each face 0093 % attached to a vertex by the area of that face, then averaging. The length of 0094 % the average normal reflects this averaging. 0095 % 0096 % Simple Example of 2-d triangles to check units 0097 % (1,5) 0098 % 0099 % (1,0) ----------- (1,1) 0100 % | | 0101 % | | 0102 % (0,0) ------------(1,0) 0103 % 0104 % 0105 % number vertices as 0106 % 5 0107 % 0108 % 4 ----------- 3 0109 % 0110 % 1 ----------- 2 0111 % 0112 % FV.vertices = [0 0 0; 1 0 0; 1 1 0; 0 1 0; 1 5 0];% in the plane of the screen 0113 % FV.faces = [1 2 3; 1 3 4; 3 5 4]; % outward from screen 0114 0115 %<autobegin> ---------------------- 21-May-2004 15:14:52 ----------------------- 0116 % --------- Automatically Generated Comments Block Using AUTO_COMMENTS --------- 0117 % 0118 % CATEGORY: Utility - Numeric 0119 % 0120 % Subfunctions in this file, in order of occurrence in file: 0121 % c = cross(a,b); 0122 % 0123 % At Check-in: $Author: psdlw $ $Revision: 1.1 $ $Date: 2004/11/12 01:32:36 $ 0124 % 0125 % This software is part of BrainStorm Toolbox Version 2.0 (Alpha) 19-May-2004 0126 % 0127 % Principal Investigators and Developers: 0128 % ** Richard M. Leahy, PhD, Signal & Image Processing Institute, 0129 % University of Southern California, Los Angeles, CA 0130 % ** John C. Mosher, PhD, Biophysics Group, 0131 % Los Alamos National Laboratory, Los Alamos, NM 0132 % ** Sylvain Baillet, PhD, Cognitive Neuroscience & Brain Imaging Laboratory, 0133 % CNRS, Hopital de la Salpetriere, Paris, France 0134 % 0135 % See BrainStorm website at http://neuroimage.usc.edu for further information. 0136 % 0137 % Copyright (c) 2004 BrainStorm by the University of Southern California 0138 % This software distributed under the terms of the GNU General Public License 0139 % as published by the Free Software Foundation. Further details on the GPL 0140 % license can be found at http://www.gnu.org/copyleft/gpl.html . 0141 % 0142 % FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE 0143 % UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY 0144 % WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF 0145 % MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY 0146 % LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE. 0147 %<autoend> ------------------------ 21-May-2004 15:14:52 ----------------------- 0148 0149 0150 % ----------------------------- Script History --------------------------------- 0151 % JCM 13-May-2004 Creation 0152 % JCM 13-May-2004 Further efficiencies for fast retrievals of constant or linear 0153 % tesselation statistics, or longer surface checks of 0154 % neighborhoods. 0155 % JCM 14-May-2004 Surface checks dramatically improved in speed using no do 0156 % loops. 0157 % JCM 15-May-2004 Comments updating, code cleaning, removal or renaming of 0158 % parameters. 0159 % ----------------------------- Script History --------------------------------- 0160 0161 0162 if ~exist('VERBOSE','var'), 0163 VERBOSE = 1; % default talkative case 0164 end 0165 0166 % Output decisions 0167 VERTEX = 0; % default flag, don't compute vertices 0168 SUSPECT = 0; % default flag, don't test ordering 0169 if nargout > 3, 0170 VERTEX = 1; % compute vertices information 0171 end 0172 if nargout > 5, 0173 SUSPECT = 1; % test ordering 0174 end 0175 0176 numTri = size(FV.faces,1); % number of faces 0177 numVert = size(FV.vertices,1); % number of vertices 0178 0179 % --------------------------- Triangle Statistics ------------------------------ 0180 % calculate the area and normals for each triangle 0181 0182 if VERBOSE, 0183 disp(sprintf('Generating statistics for %.0f faces. . .',numTri)); 0184 end 0185 0186 Vertices = FV.vertices(FV.faces',:); 0187 % each set of three rows of Vertices is one triangle 0188 0189 % now difference them to get the vectors on two sides 0190 dVertices = diff(Vertices); 0191 dVertices(3:3:end,:) = []; % remove the transition between triangles 0192 0193 % now each pair of rows in dVertices represents each triangle 0194 % row 1 is vector from 1 to 2 0195 % row 2 is vector from 2 to 3 0196 0197 % v1 = dVertices(1:2:end,:)'; % each column is side one of a triangle 0198 % v2 = dVertices(2:2:end,:)'; % side two 0199 0200 % right-hand rule, counter-clockwise ordering of the triangle yields a positive 0201 % upward area and normal. 0202 % Call a fast subfunction of this function: 0203 WeightedNormals = cross(dVertices(1:2:end,:)',dVertices(2:2:end,:)')/2; 0204 % each column is the normal for each triangle 0205 % the length the vector gives the area 0206 0207 FaceArea = sqrt(sum(WeightedNormals .* WeightedNormals)); % the area 0208 FaceNormal = WeightedNormals ./ (FaceArea([1 1 1],:)); % normalize them 0209 0210 % now calculate the centers of each triangle 0211 FaceCenter = cumsum(Vertices); 0212 FaceCenter = FaceCenter(3:3:end,:); % every third summation for every triangle 0213 FaceCenter = diff([0 0 0;FaceCenter])'/3; % average of each summation 0214 % each column is the mean of the vertices of the triangles 0215 % so now we know the center, area, and the normal vector of each triangle 0216 0217 if ~VERTEX 0218 % only wanted tesselation statistics 0219 return 0220 end 0221 0222 0223 % --------------------------- Faces Connectivity ------------------------------- 0224 % first calculate what faces go with which vertex 0225 0226 [VertexNumbering,I] = sort(FV.faces(:)); % sorted Vertex numbers 0227 0228 FacesNumbering = rem(I-1,numTri)+1; % triangle number for each Vertex 0229 0230 % For the ith vertex, then FacesNumbering(VertexNumbering == i) returns the indices of the 0231 % polygons attached to the ith vertex. 0232 % 0233 % For the set of vertices in the row vector sv (e.g sv = [3 5 115 121]), then use 0234 % [i,ignore] = find(VertexNumbering(:,ones(1,length(sv))) == (ones(size(VertexNumbering,1),1)*sv)); 0235 % (compares the Vertex numbers to the indices, extracts the row indices into i) 0236 % then FacesNumbering(i) returns the indices. Apply unique to clean up. 0237 0238 % So now we know what faces are connected to each vertex 0239 0240 0241 % --------------------------- Vertices Statistics ------------------------------ 0242 % For each vertex, there is a neighborhood of triangles 0243 % Find the mean of the centers of these triangles, and see if all normals are in 0244 % the same direction away from this center. 0245 0246 % fast analysis, no do-loops 0247 0248 if VERBOSE, 0249 disp(sprintf('Generating statistics for %.0f vertices. . .',numVert)); 0250 end 0251 0252 % First, sort and replicate weighted norms for each vertex using FacesNumbering 0253 AllWeightedNorms = cumsum(WeightedNormals(:,FacesNumbering),2); 0254 AllAverages = cumsum(FaceArea(FacesNumbering)); 0255 0256 % now extract each sum 0257 VertNdx = find(diff([VertexNumbering;Inf])); % each column where a new vertex starts 0258 % pull out the sum for each vertex, then difference from previous vertex to get 0259 % the sum for just that vertex 0260 SortedWeightedNorms = diff([[0;0;0] AllWeightedNorms(:,VertNdx)],1,2); 0261 SortedAverages = diff([0 AllAverages(:,VertNdx)]); 0262 0263 % divide by the number of faces used in the sum to get a mean 0264 NumFacesEachPatch = diff([0;VertNdx]); % the number of faces for each vertex 0265 NumFacesEachPatch = NumFacesEachPatch(:)'; % ensure row vector 0266 % the average weighted norm 0267 SortedWeightedNorms = SortedWeightedNorms ./ NumFacesEachPatch([1 1 1],:); 0268 SortedAverages = SortedAverages/3; % 1/3 assignment 0269 % the average area assigned to each vertex. 1/3 of the area of each triangle is 0270 % assigned equally to it's vertices 0271 0272 % now make sure it' assigned to the right vertex numbers 0273 VertexNormal = zeros(3,numVert); % each column is the average surface normal for each vertex 0274 VertexNormal(:,VertexNumbering(VertNdx)) = SortedWeightedNorms; 0275 VertexArea = zeros(1,numVert); 0276 VertexArea(VertexNumbering(VertNdx)) = SortedAverages; 0277 0278 0279 if ~SUSPECT, 0280 % don't want ordering tests, we're done, 0281 return 0282 end 0283 0284 % ------------------------------ Surface Check --------------------------------- 0285 % fast analysis, no do-loops 0286 0287 if VERBOSE, 0288 disp(sprintf('Examining the patch around %.0f vertices. . .',numVert)); 0289 end 0290 0291 % Want to detect if a few of the normals in a vertex patch are pointed in the 0292 % opposite direction. Rather than use the weighted vertex normal, we will form 0293 % the unweighted normal 0294 % First, sort and replicate unweighted norms for each vertex using FacesNumbering 0295 AllNorms = cumsum(FaceNormal(:,FacesNumbering),2); 0296 SortedNorms = diff([[0;0;0] AllNorms(:,VertNdx)],1,2); 0297 % the average unweighted norm 0298 SortedNorms = SortedNorms ./ NumFacesEachPatch([1 1 1],:); 0299 % now make sure it' assigned to the right vertex numbers 0300 VertexUnweightedNormal = zeros(3,numVert); % each column is the average surface normal for each vertex 0301 VertexUnweightedNormal(:,VertexNumbering(VertNdx)) = SortedNorms; 0302 0303 % form the dot product of each normal with it's average normal 0304 0305 InOut = sign(sum(FaceNormal(:,FacesNumbering) .* VertexUnweightedNormal(:,VertexNumbering))); 0306 RevDir = find(InOut < 0); % normals in the reverse direction. 0307 0308 % FacesNumbering(RevDir) gives me the Triangle numbers for ones that reversed 0309 BadTri = sort(FacesNumbering(RevDir)); 0310 BadTriNdx = find(diff([BadTri;Inf])); % changes in the BadTri numbers 0311 NumBad = diff([0;BadTriNdx]); % how many times was the triangle marked as bad 0312 0313 SuspectFace = zeros(1,numTri); % number of times a face is detected as bad. 0314 SuspectFace(BadTri(BadTriNdx)) = NumBad; 0315 0316 0317 % ---------------------------- Vertex Statistics -------------------------------- 0318 % How many faces are attached to each vertex, including unassigned ones 0319 NumFacesEachVertex = zeros(1,numVert); % 0320 NumFacesEachVertex(VertexNumbering(VertNdx)) = NumFacesEachPatch(:)'; 0321 0322 0323 % ----------------------------- VERBOSE Display -------------------------------- 0324 0325 if VERBOSE, 0326 0327 % ------------------------ Closed Surface Check ----------------------------- 0328 0329 % there may be more vertices in the description than actually used 0330 if length(VertNdx) ~= numVert, 0331 disp(sprintf('\nThere are %.0f vertices that are unused in the faces description.',... 0332 numVert - length(VertNdx))); 0333 end 0334 0335 disp(' ') 0336 disp('In a truly closed surface of triangles, then the number of triangles is') 0337 disp(sprintf(' number of triangles == 2 * (the number of vertices - 2).')); 0338 disp(sprintf('With %.0f triangles comprising %.0f vertices,',... 0339 numTri, length(VertNdx))) 0340 disp(sprintf('the numbers %.0f == %.0f here suggest:',numTri,2*(length(VertNdx) - 2))) 0341 disp(' ') 0342 if numTri == 2*(length(VertNdx) - 2), 0343 disp('CLOSED SURFACE') 0344 else 0345 disp('OPEN SURFACE') 0346 end 0347 0348 % ------------------------- Suspicious Vertices ----------------------------- 0349 0350 maxFace = max(NumFacesEachVertex); 0351 0352 disp(' '); 0353 0354 for i = 0:maxFace, 0355 disp(sprintf('There are %6.0f vertices with %3.0f faces attached.',... 0356 sum(NumFacesEachVertex == i),i)); 0357 end 0358 0359 % one and two faces are not in the interior regions 0360 ndx = find(NumFacesEachVertex == 1 | NumFacesEachVertex == 2); 0361 0362 0363 if length(ndx) > 0, 0364 PLOTLEN = input(sprintf('Plot how many of these %.0f isolated (1 or 2) vertices: ',length(ndx))); 0365 else 0366 % none found 0367 PLOTLEN = 0; % don't bother asking 0368 end 0369 0370 for i = 1:PLOTLEN, 0371 0372 % get the vertices for the bad face 0373 iVert = ndx(i); 0374 0375 % get the faces attached to this vertex 0376 % For the ith vertex, then FacesNumbering(VertexNumbering == i) returns the indices of the 0377 % polygons attached to the ith vertex. 0378 0379 fndx = FacesNumbering(VertexNumbering == iVert); % the faces attached to these vertices 0380 0381 figure 0382 set(gcf,'Name',sprintf('Vertex %.0f with %.0f faces attached',iVert,NumFacesEachVertex(iVert))); 0383 h = patch('vertices',FV.vertices,'faces',FV.faces(fndx,:),'facecolor','r','edgecolor','k'); 0384 axis equal 0385 axis vis3d 0386 hold on 0387 plot3(FaceCenter(1,fndx),FaceCenter(2,fndx),FaceCenter(3,fndx),'*'); 0388 plot3(FV.vertices(iVert,1),FV.vertices(iVert,2),FV.vertices(iVert,3),'g+'); 0389 ma = mean(FaceArea(fndx)); % mean area 0390 quiver3(FaceCenter(1,fndx),FaceCenter(2,fndx),FaceCenter(3,fndx),... 0391 FaceNormal(1,fndx),FaceNormal(2,fndx),FaceNormal(3,fndx),.25); 0392 set(h,'FaceAlpha',.8) 0393 hold off 0394 cameratoolbar('Show'); % activate the camera toolbar 0395 ret = cameratoolbar; % for strange reason, this activates the default orbit mode. 0396 drawnow 0397 end 0398 0399 0400 % -------------------------- Suspicious Faces ------------------------------- 0401 0402 disp(' '); 0403 0404 for i = 0:2, 0405 disp(sprintf('There are %6.0f faces with %1.0f suspect detects.',... 0406 sum(SuspectFace == i),i)); 0407 end 0408 0409 ndx = find(SuspectFace > 2); 0410 disp(sprintf(... 0411 '\nThere are %.0f faces with more than two suspicious direction tests.\n',... 0412 length(ndx))); 0413 0414 if length(ndx) > 0, 0415 PLOTLEN = input(sprintf('Plot how many of these %.0f suspect faces: ',length(ndx))); 0416 else 0417 % none found 0418 PLOTLEN = 0; % don't bother asking 0419 end 0420 0421 if PLOTLEN > 0 0422 disp(sprintf('Suspect face is painted green.')); 0423 end 0424 0425 for i = 1:PLOTLEN, 0426 0427 % get the vertices for the bad face 0428 iVert = FV.faces(ndx(i),:); 0429 iVert = iVert(:)'; % ensure row vector 0430 0431 % get the faces attached to these vertices 0432 % For the set of vertices in the row vector sv (e.g sv = [3 5 115 121]), then 0433 % use 0434 % [i,ignore] = find(VertexNumbering(:,ones(1,length(sv))) == (ones(size(VertexNumbering,1),1)*sv)); 0435 % (compares the Vertex numbers to the indices, extracts the row indices into i) 0436 % then FacesNumbering(i) returns the indices. 0437 0438 [fndx, ignore] = find(VertexNumbering(:,ones(1,length(iVert))) == (ones(size(VertexNumbering,1),1)*iVert)); 0439 fndx = FacesNumbering(fndx); % the faces attached to these vertices 0440 fndx = unique(fndx); % ensure unique 0441 nifndx = fndx; 0442 nifndx(find(nifndx == ndx(i))) = []; % remove the ith face 0443 0444 figure 0445 set(gcf,'Name',sprintf('Face %.0f with %.0f suspect detects',ndx(i),SuspectFace(ndx(i)))); 0446 % not the ith face 0447 h = patch('vertices',FV.vertices,'faces',FV.faces(nifndx,:),'facecolor','r','edgecolor','k'); 0448 % the ith face 0449 hi = patch('vertices',FV.vertices,'faces',FV.faces(ndx(i),:),'facecolor','g'); 0450 axis equal 0451 axis vis3d 0452 hold on 0453 plot3(FaceCenter(1,fndx),FaceCenter(2,fndx),FaceCenter(3,fndx),'*') 0454 ma = mean(FaceArea(fndx)); % mean area 0455 quiver3(FaceCenter(1,fndx),FaceCenter(2,fndx),FaceCenter(3,fndx),... 0456 FaceNormal(1,fndx),FaceNormal(2,fndx),FaceNormal(3,fndx),0.25); 0457 set(h,'FaceAlpha',.8) 0458 hold off 0459 cameratoolbar('Show'); % activate the camera toolbar 0460 ret = cameratoolbar; % for strange reason, this activates the default orbit mode. 0461 drawnow 0462 end 0463 end 0464 0465 % done 0466 0467 0468 % ------------------------------ CROSS PRODUCT --------------------------------- 0469 function c = cross(a,b); 0470 % fast computation, no permutes 0471 0472 % Calculate cross product 0473 c = [a(2,:).*b(3,:)-a(3,:).*b(2,:) 0474 a(3,:).*b(1,:)-a(1,:).*b(3,:) 0475 a(1,:).*b(2,:)-a(2,:).*b(1,:)];