mesh_vertex_normals - Calculate vertex surface normals [normals,unit_normals] = mesh_vertex_normals(FV,index) FV.vertices - vertices of mesh, Nx3 Cartesian XYZ FV.faces - triangulation of vertices (Mx3 matrix) index is an array of vertex indices (default is all vertices) normals - vertex surface normals (Nx3 matrix) unit_normals - normalized normals! This routine first calculates the surface normals of each face. It then finds all the faces that contain a specific vertex and sums across the face normals (weighted by the face area). If the faces are defined according to the right-hand rule, all their normals will be "outward" normals and the average should be a sensible value for the vertex normal. See also the VertexNormals property of the patch and surface commands.
0001 function [vertNormals,vertNormalsUnit] = mesh_vertex_normals(FV,index) 0002 0003 % mesh_vertex_normals - Calculate vertex surface normals 0004 % 0005 % [normals,unit_normals] = mesh_vertex_normals(FV,index) 0006 % 0007 % FV.vertices - vertices of mesh, Nx3 Cartesian XYZ 0008 % FV.faces - triangulation of vertices (Mx3 matrix) 0009 % 0010 % index is an array of vertex indices (default is all vertices) 0011 % 0012 % normals - vertex surface normals (Nx3 matrix) 0013 % unit_normals - normalized normals! 0014 % 0015 % This routine first calculates the surface normals of each face. It then 0016 % finds all the faces that contain a specific vertex and sums across 0017 % the face normals (weighted by the face area). 0018 % 0019 % If the faces are defined 0020 % according to the right-hand rule, all their normals will be "outward" 0021 % normals and the average should be a sensible value for the vertex normal. 0022 % 0023 % See also the VertexNormals property of the patch and surface commands. 0024 % 0025 0026 % $Revision: 1.2 $ $Date: 2005/08/14 21:18:25 $ 0027 0028 % 0029 % Licence: GNU GPL, no implied or express warranties 0030 % History: 04/2004, Darren.Weber_at_radiology.ucsf.edu 0031 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0032 0033 0034 tic; 0035 fprintf('...calculating vertex surface normals...'); 0036 0037 Nvert = size(FV.vertices,1); 0038 0039 if ~exist('index','var'), index = 1:Nvert; end 0040 if isempty(index), index = 1:Nvert; end 0041 0042 vertNormals = zeros(length(index),3); 0043 vertNormalsUnit = vertNormals; 0044 0045 for i = 1:length(index), 0046 0047 vi = index(i); 0048 0049 % To calculate a vertex normal you need to perform an iteration 0050 % beginning with a triangle that holds the vertex and perform a cross 0051 % product on the two vectors of that triangle that extend from that 0052 % vertex. This will result in a normal for that triangle. Then go 0053 % around to each triangle containing that vertex and perform a cross 0054 % product for that triangle. You will end up with a set of normals for 0055 % all the triangles. To compute the vertex normal, sum all the face 0056 % normal vectors. 0057 0058 % get all the faces that contain this vertex 0059 [faceIndexI,faceIndexJ] = find(FV.faces == vi); 0060 0061 Nface = length(faceIndexI); 0062 0063 faceNormals = zeros(Nface,3); 0064 0065 for face = 1:Nface, 0066 0067 f = faceIndexI(face); 0068 0069 vi1 = FV.faces(f,1); 0070 vi2 = FV.faces(f,2); 0071 vi3 = FV.faces(f,3); 0072 0073 v1 = FV.vertices(vi1,:); 0074 v2 = FV.vertices(vi2,:); 0075 v3 = FV.vertices(vi3,:); 0076 0077 % If the vertices are given in clockwise order, when viewed from the 0078 % outside, then following calculates the "outward" surface normals. 0079 0080 edge1 = v2 - v1; 0081 edge2 = v3 - v1; 0082 0083 faceNormals(face,:) = cross( edge2, edge1 ); 0084 0085 end 0086 0087 %faceNormalsMagnitude = vector_magnitude(faceNormals); 0088 %faceNormalsUnit = faceNormals ./ repmat(faceNormalsMagnitude,1,3); 0089 0090 % Area of Triangle = || AB x AC|| / 2 ; ie, the absolute value of 0091 % the length of the cross product of AB and AC divided by 2 0092 %faceArea = abs(faceNormalsMagnitude) / 2; 0093 0094 % Weight the faceNormals by the faceArea 0095 %vertNormals(v,:) = sum( faceNormals .* repmat(faceArea,1,3) ) / sum( faceArea ); 0096 0097 vertNormals(i,:) = sum(faceNormals) / Nface; 0098 0099 vertNormalsMagnitude = norm(vertNormals(i,:)); 0100 vertNormalsUnit(i,:) = vertNormals(i,:) / vertNormalsMagnitude; 0101 0102 end 0103 0104 t=toc; 0105 fprintf('done (%5.2f sec).\n',t); 0106 0107 return