function [coordNORMALS,varargout] = COMPUTE_mesh_normals(meshdataIN,invertYN)
% COMPUTE_mesh_normals  Calculate the normals for each facet of a triangular mesh
%==========================================================================
% AUTHOR        Adam H. Aitkenhead
% CONTACT       adam.aitkenhead@physics.cr.man.ac.uk
% INSTITUTION   The Christie NHS Foundation Trust
% DATE          March 2010
% PURPOSE       Calculate the normal vectors for each facet of a triangular
%               mesh.  The ordering of the vertices
%               (clockwise/anticlockwise) is also checked for all facets if
%               this is requested as one of the outputs.
%
% USAGE         [coordNORMALS] = COMPUTE_mesh_normals(meshdataIN)
%       ..or..  [coordNORMALS,meshdataOUT] = COMPUTE_mesh_normals(meshdataIN,invertYN)
%
% INPUTS
%
%    meshdataIN   - (structure)  Structure containing the faces and
%                   vertices of the mesh, in the same format as that
%                   produced by the isosurface command.
%         ..or..  - (Nx3x3 array)  The vertex coordinates for each facet,
%                   with:  1 row for each facet
%                          3 columns for the x,y,z coordinates
%                          3 pages for the three vertices
%    invertYN     - (optional)  A flag to say whether the mesh is to be
%                   inverted or not.  Should be 'y' or 'n'.
%
% OUTPUTS
%
%    coordNORMALS - Nx3 array   - The normal vectors for each facet, with:
%                          1 row for each facet
%                          3 columns for the x,y,z components
%
%    meshdataOUT  - (optional)  - The mesh data with the ordering of the
%                   vertices (clockwise/anticlockwise) checked.  Uses the
%                   same format as <meshdataIN>.
%
% NOTES       - Computing <meshdataOUT> to check the ordering of the
%               vertices in each facet may be slow for large meshes.
%             - It may not be possible to compute <meshdataOUT> for
%               non-manifold meshes.
%==========================================================================

%==========================================================================
% VERSION  USER  CHANGES
% -------  ----  -------
% 100331   AHA   Original version
% 101129   AHA   Can now check the ordering of the facet vertices.
% 101130   AHA   <meshdataIN> can now be in either of two formats.
% 101201   AHA   Only check the vertex ordering if that is required as one
%                of the outputs, as it can be slow for large meshes.
% 101201   AHA   Add the flag invertYN and make it possible to invert the
%                mesh
% 111004   AHA   Housekeeping tidy-up
%
%----------------------------------------------------------------------------
%Copyright © 2013. Adam H. Aitkenhead. All rights reserved.
%Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
%1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
%2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in
%   the documentation and/or other materials provided with the distribution.
%3. Neither the name of the Christie NHS Foundation Trust nor the names of its contributors may be used to endorse or promote products
%   derived from this software without specific prior written permission.
%THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ''AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, 
%THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS 
%BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 
%SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER 
%IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED 
%OF THE POSSIBILITY OF SUCH DAMAGE.
%==========================================================================

%======================================================
% Read the input parameters
%======================================================

if isstruct(meshdataIN)==1
  faces         = meshdataIN.faces;
  vertex        = meshdataIN.vertices;
  coordVERTICES = zeros(size(faces,1),3,3);
  for loopa = 1:size(faces,1)
    coordVERTICES(loopa,:,1) = vertex(faces(loopa,1),:);
    coordVERTICES(loopa,:,2) = vertex(faces(loopa,2),:);
    coordVERTICES(loopa,:,3) = vertex(faces(loopa,3),:);
  end
else
  coordVERTICES = meshdataIN;
end

%======================================================
% Invert the mesh if required
%======================================================

if exist('invertYN','var')==1 && isempty(invertYN)==0 && ischar(invertYN)==1 && ( strncmpi(invertYN,'y',1)==1 || strncmpi(invertYN,'i',1)==1 )
  coV           = zeros(size(coordVERTICES));
  coV(:,:,1)    = coordVERTICES(:,:,1);
  coV(:,:,2)    = coordVERTICES(:,:,3);
  coV(:,:,3)    = coordVERTICES(:,:,2);
  coordVERTICES = coV;
end

%======================
% Initialise array to hold the normal vectors
%======================

facetCOUNT   = size(coordVERTICES,1);
coordNORMALS = zeros(facetCOUNT,3);

%======================
% Check the vertex ordering for each facet
%======================

if nargout==2
  startfacet  = 1;
  edgepointA  = 1;
  checkedlist = false(facetCOUNT,1);
  waitinglist = false(facetCOUNT,1);

  while min(checkedlist)==0
    
    checkedlist(startfacet) = 1;

    edgepointB = edgepointA + 1;
    if edgepointB==4
      edgepointB = 1;
    end
    
    %Find points which match edgepointA
    sameX = coordVERTICES(:,1,:)==coordVERTICES(startfacet,1,edgepointA);
    sameY = coordVERTICES(:,2,:)==coordVERTICES(startfacet,2,edgepointA);
    sameZ = coordVERTICES(:,3,:)==coordVERTICES(startfacet,3,edgepointA);
    [tempa,tempb] = find(sameX & sameY & sameZ);
    matchpointA = [tempa,tempb];
    matchpointA = matchpointA(matchpointA(:,1)~=startfacet,:);
  
    %Find points which match edgepointB
    sameX = coordVERTICES(:,1,:)==coordVERTICES(startfacet,1,edgepointB);
    sameY = coordVERTICES(:,2,:)==coordVERTICES(startfacet,2,edgepointB);
    sameZ = coordVERTICES(:,3,:)==coordVERTICES(startfacet,3,edgepointB);
    [tempa,tempb] = find(sameX & sameY & sameZ);
    matchpointB = [tempa,tempb];
    matchpointB = matchpointB(matchpointB(:,1)~=startfacet,:);
  
    %Find edges which match both edgepointA and edgepointB -> giving the adjacent edge
    [memberA,memberB] = ismember(matchpointA(:,1),matchpointB(:,1));
    matchfacet = matchpointA(memberA,1);
  
    if numel(matchfacet)~=1
      if exist('warningdone','var')==0
        warning('Mesh is non-manifold.')
        warningdone = 1;
      end
    else
      matchpointA = matchpointA(memberA,2);
      matchpointB = matchpointB(memberB(memberA),2);
      
      if checkedlist(matchfacet)==0 && waitinglist(matchfacet)==0
        %Ensure the adjacent edge is traveled in the opposite direction to the original edge  
        if matchpointB-matchpointA==1 || matchpointB-matchpointA==-2
          %Direction needs to be flipped
          [ coordVERTICES(matchfacet,:,matchpointA) , coordVERTICES(matchfacet,:,matchpointB) ] = deal( coordVERTICES(matchfacet,:,matchpointB) , coordVERTICES(matchfacet,:,matchpointA) );
        end
      end
    end
  
    waitinglist(matchfacet) = 1;
    
    if edgepointA<3
      edgepointA = edgepointA + 1;
    elseif edgepointA==3
      edgepointA = 1;
      checkedlist(startfacet) = 1;
      startfacet = find(waitinglist==1 & checkedlist==0,1,'first');
    end
  
  end
end

%======================
% Compute the normal vector for each facet
%======================

for loopFACE = 1:facetCOUNT
  
  %Find the coordinates for each vertex.
  cornerA = coordVERTICES(loopFACE,1:3,1);
  cornerB = coordVERTICES(loopFACE,1:3,2);
  cornerC = coordVERTICES(loopFACE,1:3,3);
  
  %Compute the vectors AB and AC
  AB = cornerB-cornerA;
  AC = cornerC-cornerA;
    
  %Determine the cross product AB x AC
  ABxAC = cross(AB,AC);
    
  %Normalise to give a unit vector
  ABxAC = ABxAC / norm(ABxAC);
  coordNORMALS(loopFACE,1:3) = ABxAC;
  
end %loopFACE

%======================================================
% Prepare the output parameters
%======================================================

if nargout==2
  if isstruct(meshdataIN)==1
    [faces,vertices] = CONVERT_meshformat(coordVERTICES);
    meshdataOUT = struct('vertices',vertices,'faces',faces);
  else
    meshdataOUT = coordVERTICES;
  end
  varargout(1) = {meshdataOUT};
end


%======================================================
end %function


