[Master Index]
[Index for Toolbox]
tessellation_stats
(Toolbox/tessellation_stats.m in BrainStorm 2.0 (Alpha))
Function Synopsis
[FaceNormal, FaceArea, FaceCenter, VertexNormal, VertexArea, SuspectFace, NumFacesEachVertex, Duplicated_Faces, Not_Twice_Faces] = tessellation_stats(FV,VERBOSE);
Help Text
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.
Cross-Reference Information
This function calls
- sort_key C:\BrainStorm_2001\Toolbox\sort_key.m
This function is called by
Listing of function C:\BrainStorm_2001\Toolbox\tessellation_stats.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