[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

    tessellation_outwards              C:\BrainStorm_2001\Toolbox\tessellation_outwards.m

### 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 Nov 2 16:17:59 2004
Cross-Directory links are: ON