[FaceNormal, FaceArea, FaceCenter, VertexNormal, VertexArea, SuspectFace, NumFacesEachVertex, Duplicated_Faces, Not_Twice_Faces] = tessellation_stats(FV,VERBOSE);

TESSELLATION_STATS - Calculate statistics of the tesselation and hunt for suspicious faces and vertices function [FaceNormal, FaceArea, FaceCenter, VertexNormal, VertexArea, SuspectFace, NumFacesEachVertex, Duplicated_Faces, Not_Twice_Faces] = tessellation_stats(FV,VERBOSE); 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. Duplicated_Faces, Not_Twice_Faces Each is a cell array, each element contains the one or more faces that are adjacent a bad edge. Duplicated_Faces are adjacent an edge that was specified more than once in the same direction. Not_Twice_Faces are adjacent an edge that was not specified exactly once in each direction, but are not in the set of duplicated 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.

- sort_key C:\BrainStorm_2001\Toolbox\sort_key.m

- tessellation_outwards C:\BrainStorm_2001\Toolbox\tessellation_outwards.m

function [FaceNormal, FaceArea, FaceCenter, VertexNormal, VertexArea, SuspectFace, NumFacesEachVertex, Duplicated_Faces, Not_Twice_Faces] = tessellation_stats(FV,VERBOSE); %TESSELLATION_STATS - Calculate statistics of the tesselation and hunt for suspicious faces and vertices % function [FaceNormal, FaceArea, FaceCenter, VertexNormal, VertexArea, SuspectFace, NumFacesEachVertex, Duplicated_Faces, Not_Twice_Faces] = tessellation_stats(FV,VERBOSE); % 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. % % Duplicated_Faces, Not_Twice_Faces % Each is a cell array, each element contains the one or more faces that are % adjacent a bad edge. Duplicated_Faces are adjacent an edge that was % specified more than once in the same direction. Not_Twice_Faces are adjacent % an edge that was not specified exactly once in each direction, but are not % in the set of duplicated 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. %<autobegin> ---------------------- 26-May-2004 11:34:37 ----------------------- % --------- Automatically Generated Comments Block Using AUTO_COMMENTS --------- % % CATEGORY: Utility - Numeric % % Alphabetical list of external functions (non-Matlab): % toolbox\sort_key.m % % Subfunctions in this file, in order of occurrence in file: % c = cross(a,b); % % At Check-in: $Author: Mosher $ $Revision: 9 $ $Date: 5/26/04 10:02a $ % % This software is part of BrainStorm Toolbox Version 2.0 (Alpha) 24-May-2004 % % Principal Investigators and Developers: % ** Richard M. Leahy, PhD, Signal & Image Processing Institute, % University of Southern California, Los Angeles, CA % ** John C. Mosher, PhD, Biophysics Group, % Los Alamos National Laboratory, Los Alamos, NM % ** Sylvain Baillet, PhD, Cognitive Neuroscience & Brain Imaging Laboratory, % CNRS, Hopital de la Salpetriere, Paris, France % % See BrainStorm website at http://neuroimage.usc.edu for further information. % % Copyright (c) 2004 BrainStorm by the University of Southern California % This software distributed under the terms of the GNU General Public License % as published by the Free Software Foundation. Further details on the GPL % license can be found at http://www.gnu.org/copyleft/gpl.html . % % FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE % UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY % WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF % MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY % LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE. %<autoend> ------------------------ 26-May-2004 11:34:37 ----------------------- % ----------------------------- Script History --------------------------------- % JCM 13-May-2004 Creation % JCM 13-May-2004 Further efficiencies for fast retrievals of constant or linear % tesselation statistics, or longer surface checks of % neighborhoods. % JCM 14-May-2004 Surface checks dramatically improved in speed using no do % loops. % JCM 15-May-2004 Comments updating, code cleaning, removal or renaming of % parameters. % ----------------------------- Script History --------------------------------- % 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 screen if ~exist('VERBOSE','var'), VERBOSE = 1; % default talkative case end % Output decisions VERTEX = 0; % default flag, don't compute vertices SUSPECT = 0; % default flag, don't test ordering if nargout > 3, VERTEX = 1; % compute vertices information end if nargout > 5, SUSPECT = 1; % test ordering end numTri = size(FV.faces,1); % number of faces numVert = size(FV.vertices,1); % number of vertices % --------------------------- Triangle Statistics ------------------------------ % calculate the area and normals for each triangle if VERBOSE, disp(sprintf('Generating statistics for %.0f faces. . .',numTri)); end Vertices = FV.vertices(FV.faces',:); % each set of three rows of Vertices is one triangle % now difference them to get the vectors on two sides dVertices = diff(Vertices); dVertices(3:3:end,:) = []; % remove the transition between triangles % now each pair of rows in dVertices represents each triangle % row 1 is vector from 1 to 2 % row 2 is vector from 2 to 3 % v1 = dVertices(1:2:end,:)'; % each column is side one of a triangle % v2 = dVertices(2:2:end,:)'; % side two % right-hand rule, counter-clockwise ordering of the triangle yields a positive % upward area and normal. % Call a fast subfunction of this function: WeightedNormals = cross(dVertices(1:2:end,:)',dVertices(2:2:end,:)')/2; % each column is the normal for each triangle % the length the vector gives the area FaceArea = sqrt(sum(WeightedNormals .* WeightedNormals)); % the area FaceNormal = WeightedNormals ./ (FaceArea([1 1 1],:)); % normalize them % now calculate the centers of each triangle FaceCenter = cumsum(Vertices); FaceCenter = FaceCenter(3:3:end,:); % every third summation for every triangle FaceCenter = diff([0 0 0;FaceCenter])'/3; % average of each summation % each column is the mean of the vertices of the triangles % so now we know the center, area, and the normal vector of each triangle if ~VERTEX % only wanted tesselation statistics return end % --------------------------- Faces Connectivity ------------------------------- % first calculate what faces go with which vertex [VertexNumbering,I] = sort(FV.faces(:)); % sorted Vertex numbers FacesNumbering = rem(I-1,numTri)+1; % triangle number for each Vertex % For the ith vertex, then FacesNumbering(VertexNumbering == i) returns the indices of the % polygons attached to the ith vertex. % % For the set of vertices in the row vector sv (e.g sv = [3 5 115 121]), then use % [i,ignore] = find(VertexNumbering(:,ones(1,length(sv))) == (ones(size(VertexNumbering,1),1)*sv)); % (compares the Vertex numbers to the indices, extracts the row indices into i) % then FacesNumbering(i) returns the indices. Apply unique to clean up. % So now we know what faces are connected to each vertex % --------------------------- Vertices Statistics ------------------------------ % For each vertex, there is a neighborhood of triangles % Find the mean of the centers of these triangles, and see if all normals are in % the same direction away from this center. % fast analysis, no do-loops if VERBOSE, disp(sprintf('Generating statistics for %.0f vertices. . .',numVert)); end % First, sort and replicate weighted norms for each vertex using FacesNumbering AllWeightedNorms = cumsum(WeightedNormals(:,FacesNumbering),2); AllAverages = cumsum(FaceArea(FacesNumbering)); % now extract each sum VertNdx = find(diff([VertexNumbering;Inf])); % each column where a new vertex starts % pull out the sum for each vertex, then difference from previous vertex to get % the sum for just that vertex SortedWeightedNorms = diff([[0;0;0] AllWeightedNorms(:,VertNdx)],1,2); SortedAverages = diff([0 AllAverages(:,VertNdx)]); % divide by the number of faces used in the sum to get a mean NumFacesEachPatch = diff([0;VertNdx]); % the number of faces for each vertex NumFacesEachPatch = NumFacesEachPatch(:)'; % ensure row vector % the average weighted norm SortedWeightedNorms = SortedWeightedNorms ./ NumFacesEachPatch([1 1 1],:); SortedAverages = SortedAverages/3; % 1/3 assignment % the average area assigned to each vertex. 1/3 of the area of each triangle is % assigned equally to it's vertices % now make sure it' assigned to the right vertex numbers VertexNormal = zeros(3,numVert); % each column is the average surface normal for each vertex VertexNormal(:,VertexNumbering(VertNdx)) = SortedWeightedNorms; VertexArea = zeros(1,numVert); VertexArea(VertexNumbering(VertNdx)) = SortedAverages; if ~SUSPECT, % don't want ordering tests, we're done, return end % ------------------------------ Surface Check --------------------------------- % fast analysis, no do-loops if VERBOSE, disp(sprintf('Examining the patch around %.0f vertices. . .',numVert)); end % Want to detect if a few of the normals in a vertex patch are pointed in the % opposite direction. Rather than use the weighted vertex normal, we will form % the unweighted normal % First, sort and replicate unweighted norms for each vertex using FacesNumbering AllNorms = cumsum(FaceNormal(:,FacesNumbering),2); SortedNorms = diff([[0;0;0] AllNorms(:,VertNdx)],1,2); % the average unweighted norm SortedNorms = SortedNorms ./ NumFacesEachPatch([1 1 1],:); % now make sure it' assigned to the right vertex numbers VertexUnweightedNormal = zeros(3,numVert); % each column is the average surface normal for each vertex VertexUnweightedNormal(:,VertexNumbering(VertNdx)) = SortedNorms; % form the dot product of each normal with it's average normal InOut = sign(sum(FaceNormal(:,FacesNumbering) .* VertexUnweightedNormal(:,VertexNumbering))); RevDir = find(InOut < 0); % normals in the reverse direction. % FacesNumbering(RevDir) gives me the Triangle numbers for ones that reversed BadTri = sort(FacesNumbering(RevDir)); BadTriNdx = find(diff([BadTri;Inf])); % changes in the BadTri numbers NumBad = diff([0;BadTriNdx]); % how many times was the triangle marked as bad SuspectFace = zeros(1,numTri); % number of times a face is detected as bad. SuspectFace(BadTri(BadTriNdx)) = NumBad; % ---------------------------- Vertex Statistics -------------------------------- % How many faces are attached to each vertex, including unassigned ones NumFacesEachVertex = zeros(1,numVert); % NumFacesEachVertex(VertexNumbering(VertNdx)) = NumFacesEachPatch(:)'; % ----------------------------- Edge Statistics -------------------------------- % schematic representation of edges [odd edge vertex, even edge vertex]; % Edges = FV.faces(:,[1 2 2 3 3 1]); % each pair of columns is an edge % The row number is the face number % Order all of the odd columns adjacent the even columns % First column is the odd edge number, then second column is the even edge number Edges = reshape(FV.faces(:,[1 2 3 2 3 1]),[],2); % columns reordered for reshaping if VERBOSE, disp(sprintf('Sorting %.0f edge descriptions . . .',size(Edges,1))); end [Edges,EdgeDirection] = sort(Edges,2); % sort each row, retaining direction % EdgeDirection(i,:) gives the ordering [1 2] or [2 1] for each Edge(i,:) [Edges_and_Directions, EdgeFaceNumber] = sort_key([Edges EdgeDirection]); EdgeFaceNumber = rem(EdgeFaceNumber-1,size(FV.faces,1))+1; % adjust the block lex ordering % EdgeFaceNumber(i) gives us the corresponding face number for the Edge(i,:) % so the vector [Edges_and_Directions(i,:) EdgeFaceNnumber(i)] % gives use information about the ith edge % In a triangle manifold, each edge should have been used twice, once in each % direction. By the keyed sorting, the directions are also sorted for each edge % Diff the edge numberings with their directions diff_Edges_and_Directions = diff([0 0 0 0; Edges_and_Directions]); % Each time an edge changes, then diff catches it numEdges = sum(any(diff_Edges_and_Directions(:,1:2),2)); if VERBOSE, disp(sprintf('Examining %.0f distinct edges . . .',numEdges)); end % if an edge in the same direction was specified more than once, then we can % detect as ndx_Duplicated_Edge = find(~any(diff_Edges_and_Directions,2)); % if an edge is repeated properly then the difference is zero in the edge % numbering but is different in the edge direction ndx_Repeated_Edge = find(~any(diff_Edges_and_Directions(:,1:2),2)); % all edges should have been specified twice, find those that were not diff_ndx = diff([0;ndx_Repeated_Edge]); ndx_NotTwice = find(diff_ndx ~= 2); % something wrong with this edge ndx_NotTwice = ndx_Repeated_Edge(ndx_NotTwice); % in the original Edges_and_Directions % so Edges_and_Directions([ndx_NotTwice;ndx_Duplicated_Edge]) are problem edges, % maybe duplicated between the sets. An edge specified only once would be caught % in NotTwice only. Duplicated_Edges = unique(Edges_and_Directions(ndx_Duplicated_Edge,1:2),'rows'); Not_Twice_Edges = unique(Edges_and_Directions(ndx_NotTwice,1:2),'rows'); % remove the duplicated edges from the not Twice Not_Twice_Edges = setdiff(Not_Twice_Edges,Duplicated_Edges,'rows'); % now get all faces for these edges, for visualization purposes and tracking % there should not be that many, so don't worry about do loop Duplicated_Faces = cell(1,size(Duplicated_Edges,1)); for i = 1:size(Duplicated_Edges,1), Duplicated_Faces{i} = EdgeFaceNumber(find(Edges_and_Directions(:,1) == Duplicated_Edges(i,1) & ... Edges_and_Directions(:,2) == Duplicated_Edges(i,2)))'; end Not_Twice_Faces = cell(1,size(Not_Twice_Edges,1)); for i = 1:size(Not_Twice_Edges,1), Not_Twice_Faces{i} = EdgeFaceNumber(find(Edges_and_Directions(:,1) == Not_Twice_Edges(i,1) & ... Edges_and_Directions(:,2) == Not_Twice_Edges(i,2)))'; end % ----------------------------- VERBOSE Display -------------------------------- if VERBOSE, % ------------------------ Closed Surface Check ----------------------------- % there may be more vertices in the description than actually used if length(VertNdx) ~= numVert, disp(sprintf('\nThere are %.0f vertices that are unused in the faces description.',... numVert - length(VertNdx))); end EulerCharacteristic = numTri + length(VertNdx) - numEdges; disp(' ') disp('Poincaré Formula Check for Triangles:') disp('In a truly closed surface of triangles, then the number of triangles is') disp(sprintf(' number of triangles == 2 * (the number of vertices - 2).')); disp(sprintf('With %.0f triangles comprising %.0f vertices,',... numTri, length(VertNdx))) disp(sprintf('the numbers %.0f == %.0f here suggest:',numTri,2*(length(VertNdx) - 2))) disp(' ') if numTri == 2*(length(VertNdx) - 2), disp('CLOSED SURFACE') else disp('OPEN SURFACE') end disp(' ') disp('General Poincaré Formula Check:') disp(sprintf(... 'There are %.0f faces, %.0f vertices, %.0f edges,',numTri,length(VertNdx),numEdges)); disp(sprintf('such that the Poincaré Formula is %.0f + %.0f - %.0f = %.0f',... numTri,length(VertNdx),numEdges,EulerCharacteristic)) if EulerCharacteristic > 0 & ~rem(EulerCharacteristic,2), % positive even number disp(sprintf('\nThe Euler Characteristic of %.0f suggests a surface of genus %.0f',... EulerCharacteristic,round((EulerCharacteristic - 2)/2))); end disp(' ') if EulerCharacteristic == 2, disp('CLOSED SURFACE') else disp('Not surface of genus 0 (not a closed surface)'); end % Edges at the boundary or adjacent faces that are reverse wound disp(' '); if ~isempty(Duplicated_Faces), disp(sprintf('BAD, there are %.0f edges that were duplicated in the face descriptions.',... length(Duplicated_Faces))); disp(sprintf('The faces adjacent these edges are in Duplicated_Faces')); else disp('Good, no edges were duplicated in the face descriptions.'); end disp(' '); if ~isempty(Not_Twice_Faces), disp(sprintf('BAD, there are %.0f edges that were not specified twice, once in each direction.',... length(Not_Twice_Faces))); disp(sprintf('The faces adjacent these edges are in Not_Twice_Faces')); else disp('Good, all edges were properly used twice, once in each direction.'); end % ------------------------- Suspicious Vertices ----------------------------- maxFace = max(NumFacesEachVertex); disp(' '); for i = 0:maxFace, disp(sprintf('There are %6.0f vertices with %3.0f faces attached.',... sum(NumFacesEachVertex == i),i)); end % one and two faces are not in the interior regions ndx = find(NumFacesEachVertex == 1 | NumFacesEachVertex == 2); if length(ndx) > 0, PLOTLEN = input(sprintf('Plot how many of these %.0f isolated (1 or 2) vertices: ',length(ndx))); else % none found PLOTLEN = 0; % don't bother asking end for i = 1:PLOTLEN, % get the vertices for the bad face iVert = ndx(i); % get the faces attached to this vertex % For the ith vertex, then FacesNumbering(VertexNumbering == i) returns the indices of the % polygons attached to the ith vertex. fndx = FacesNumbering(VertexNumbering == iVert); % the faces attached to these vertices figure set(gcf,'Name',sprintf('Vertex %.0f with %.0f faces attached',iVert,NumFacesEachVertex(iVert))); h = patch('vertices',FV.vertices,'faces',FV.faces(fndx,:),'facecolor','r','edgecolor','k'); axis equal axis vis3d hold on plot3(FaceCenter(1,fndx),FaceCenter(2,fndx),FaceCenter(3,fndx),'*'); plot3(FV.vertices(iVert,1),FV.vertices(iVert,2),FV.vertices(iVert,3),'g+'); ma = mean(FaceArea(fndx)); % mean area quiver3(FaceCenter(1,fndx),FaceCenter(2,fndx),FaceCenter(3,fndx),... FaceNormal(1,fndx),FaceNormal(2,fndx),FaceNormal(3,fndx),.25); set(h,'FaceAlpha',.8) hold off cameratoolbar('Show'); % activate the camera toolbar ret = cameratoolbar; % for strange reason, this activates the default orbit mode. drawnow end % -------------------------- Suspicious Faces ------------------------------- disp(' '); for i = 0:2, disp(sprintf('There are %6.0f faces with %1.0f suspect detects.',... sum(SuspectFace == i),i)); end ndx = find(SuspectFace > 2); disp(sprintf(... '\nThere are %.0f faces with more than two suspicious direction tests.\n',... length(ndx))); if length(ndx) > 0, PLOTLEN = input(sprintf('Plot how many of these %.0f suspect faces: ',length(ndx))); else % none found PLOTLEN = 0; % don't bother asking end if PLOTLEN > 0 disp(sprintf('Suspect face is painted green.')); end for i = 1:PLOTLEN, % get the vertices for the bad face iVert = FV.faces(ndx(i),:); iVert = iVert(:)'; % ensure row vector % get the faces attached to these vertices % For the set of vertices in the row vector sv (e.g sv = [3 5 115 121]), then % use % [i,ignore] = find(VertexNumbering(:,ones(1,length(sv))) == (ones(size(VertexNumbering,1),1)*sv)); % (compares the Vertex numbers to the indices, extracts the row indices into i) % then FacesNumbering(i) returns the indices. [fndx, ignore] = find(VertexNumbering(:,ones(1,length(iVert))) == (ones(size(VertexNumbering,1),1)*iVert)); fndx = FacesNumbering(fndx); % the faces attached to these vertices fndx = unique(fndx); % ensure unique nifndx = fndx; nifndx(find(nifndx == ndx(i))) = []; % remove the ith face figure set(gcf,'Name',sprintf('Face %.0f with %.0f suspect detects',ndx(i),SuspectFace(ndx(i)))); % not the ith face h = patch('vertices',FV.vertices,'faces',FV.faces(nifndx,:),'facecolor','r','edgecolor','k'); % the ith face hi = patch('vertices',FV.vertices,'faces',FV.faces(ndx(i),:),'facecolor','g'); axis equal axis vis3d hold on plot3(FaceCenter(1,fndx),FaceCenter(2,fndx),FaceCenter(3,fndx),'*') ma = mean(FaceArea(fndx)); % mean area quiver3(FaceCenter(1,fndx),FaceCenter(2,fndx),FaceCenter(3,fndx),... FaceNormal(1,fndx),FaceNormal(2,fndx),FaceNormal(3,fndx),0.25); set(h,'FaceAlpha',.8) hold off cameratoolbar('Show'); % activate the camera toolbar ret = cameratoolbar; % for strange reason, this activates the default orbit mode. drawnow end end % done % ------------------------------ CROSS PRODUCT --------------------------------- function c = cross(a,b); % fast computation, no permutes % Calculate cross product c = [a(2,:).*b(3,:)-a(3,:).*b(2,:) a(3,:).*b(1,:)-a(1,:).*b(3,:) a(1,:).*b(2,:)-a(2,:).*b(1,:)];

Produced by color_mat2html, a customized BrainStorm 2.0 (Alpha) version of mat2html on Tue Oct 12 12:05:14 2004

Cross-Directory links are: ON