function mesh = calcMeshNormals(mesh, refNormal)
% mesh = calcMeshNormals(mesh, refNormal)
%
% Author : Ajmal Saeed Mian {ajmal@csse.uwa.edu.au}
% Computer Science. Univ of Western Australia
%
% This function takes 3D points and polygon information and creates a new mesh
% from it.
%
% Arguments: mesh - a data structure containing only mesh.vertices and
% mesh.triangles
% refNormal - (optional) vector towards the sensor
% default is [0 0 1]
%
% Return: mesh - a complete mesh data structure with the following fields
% vertices: Nx3 vertex coordinates
% triangles: M triangles using index numbers of the vertices
% resolution: the mean edge length of triangles
% stdeviation: the standard deviation o edge lengths
% triangleNormals: Mx3 normal vectors of each triangle
% vertexNormals: Nx3 normal vectors of each vertex
% vertexNtriangles: Nx1 cell of neighboring triangles
% of each vertex
% triangleNtriangles: Mx1 cell of nieghboring triangles
% of each triangle
%
% Copyright : This code is written by Ajmal Saeed Mian {ajmal@csse.uwa.edu.au}
% Computer Science, The University of Western Australia. The code
% may be used, modified and distributed for research purposes with
% acknowledgement of the author and inclusion this copyright information.
%
% Disclaimer : This code is provided as is without any warrantly.
noOfpolygons = size(mesh.triangles,1);
noOfpoints = size(mesh.vertices,1);
mesh.vertexNormals = zeros(noOfpolygons,3); % innitialize a matrix to store polygon normals
mesh.vertexNormals = zeros(noOfpoints,3); % innitialize a matrix to store point normals
mesh.vertexNtriangles = cell(noOfpoints, 1); %a cell array to store neighbouring polygons for the current point
mesh.triangleNtriangles = cell(noOfpolygons, 1); % to store neighbors of current polygon
for ii = 1:noOfpolygons %find normals of all polygons
%indices of the points from which the polygon is made
pointIndex1 = mesh.triangles(ii,1);
pointIndex2 = mesh.triangles(ii,2);
pointIndex3 = mesh.triangles(ii,3);
%coordinates of the points
point1 = mesh.vertices(pointIndex1,:);
point2 = mesh.vertices(pointIndex2,:);
point3 = mesh.vertices(pointIndex3,:);
vector1 = point2 - point1;
vector2 = point3 - point2;
normal = cross(vector1,vector2);
normal = normal / norm(normal);
if nargin == 3
theta = acos(dot(refNormal, normal));
if theta > pi/2
normal = normal * (-1);
a = mesh.triangles(ii,2);
mesh.triangles(ii,2) = mesh.triangles(ii,1);
mesh.triangles(ii,1) = a;
end
end
mesh.vertexNormals(ii,:)=normal;
%make entry of this polygon as the neighbouring polygon of the three
%vertex points
mesh.vertexNtriangles(pointIndex1,1)={[mesh.vertexNtriangles{pointIndex1,1} ii]};
mesh.vertexNtriangles(pointIndex2,1)={[mesh.vertexNtriangles{pointIndex2,1} ii]};
mesh.vertexNtriangles(pointIndex3,1)={[mesh.vertexNtriangles{pointIndex3,1} ii]};
end
% find neighbouring polygons of all polygons
for ii = 1 : noOfpolygons
polNeighbor = [];
for jj = 1 : 3
polNeighbor = [polNeighbor mesh.vertexNtriangles{mesh.triangles(ii,jj)}];
end
polNeighbor = unique(polNeighbor);
polNeighbor = setdiff(polNeighbor, [ii]);
mesh.triangleNtriangles(ii,1)={[polNeighbor]};
end
%find normals of all points
for ii = 1 : noOfpoints
polys = mesh.vertexNtriangles{ii}; % get neighboring polygons to this point
% get neighboring polygons to all these polygons
%for jj = 1 : length(polys)
% polys = [polys mesh.triangleNtriangles{polys(jj)}];
%end
%polys = unique(polys); % remove redundancy
normal2 = zeros(1,3);
for jj = 1 : length(polys)
normal2 = normal2 + mesh.vertexNormals(polys(jj),:)*...
polyArea([mesh.vertices(mesh.triangles(polys(jj),1),:); ...
mesh.vertices(mesh.triangles(polys(jj),2),:); ...
mesh.vertices(mesh.triangles(polys(jj),3),:)], mesh.vertexNormals(polys(jj),:));
end
normal2 = normal2 / norm(normal2);
mesh.vertexNormals(ii,:) = normal2;
end
edgeLength = [sqrt(sum((mesh.vertices(mesh.triangles(:,1),:) - mesh.vertices(mesh.triangles(:,2),:)).^2,2)),...
sqrt(sum((mesh.vertices(mesh.triangles(:,2),:) - mesh.vertices(mesh.triangles(:,3),:)).^2,2)),...
sqrt(sum((mesh.vertices(mesh.triangles(:,3),:) - mesh.vertices(mesh.triangles(:,1),:)).^2,2))];
resolution = mean(edgeLength(:));
stdeviation = std(edgeLength(:));
mesh.resolution = resolution;
mesh.stdeviation = stdeviation;
