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.
% ---------------------- 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.
% ------------------------ 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,:)];