
function [xv,yv,zv,R10,RFM,abcd]=maptopolygon2(XV,YV,ZV)

%--------------------------------------------------------------------------
% Object:
% This routine takes a polygon in the 3D space whose vertices are (XV,YV,ZV)
% and maps it to a flat polygon with vertices (xv,yv,0) by a "rigid
% motion".
% The vector "zv" and the rotation matrix "R10" allow to pull back the
% points.
%--------------------------------------------------------------------------
% Input:
%
% XV,YV,ZV: coordinates of the vertices of a polygon in the three
% dimensional space (e.g. the vertices of a facet of a tetrahedra); last
% vertex may not be repeated.
%--------------------------------------------------------------------------
% Output:
%
% xv,yv: coordinates of the vertices of a polygon in the two
% dimensional space that is obtained by a "rigid motion" of a polygon with
% vertices (XV,YV,ZV); last vertex may not be repeated.
%
% zv, R10: the polygon P with coordinates (xv,yv,0) is
%          (R10*[XV'; YV'; ZV'])'-[0V 0V ZV]
%     where 0V is the null vector with the same dimension of ZV.
%
% R10,RFM: rotation and reflection matrices involved in the process.
%     It is "P=R10'*(RFM*P2+PS)" where "PS" is a certain shift, "P2" a 
%     the plane xy and "P" is a point in the 3D plane.
%
% abcd: vector whose components [a b c d] determine the equation of the
%     plane "a*x+b*y+c*z+d=0" containing the vertices of the polygon.
%     
%--------------------------------------------------------------------------

% Determine the needed rigid motion. We map A in 0, B in a point with
% coordinates (dist(A,B),0,0) and the versor of the normal into (0,0,1).

A=[XV(1) YV(1) ZV(1)]; A=A';
B=[XV(2) YV(2) ZV(2)]; B=B';
C=[XV(3) YV(3) ZV(3)]; C=C';

% Equation of the plane containing the facets.
% See e.g.:
% https://math.stackexchange.com/questions/2686606/equation-of-a-plane-passing-through-3-points
n=cross(B-A,C-A);
a=n(1); b=n(2); c=n(3);
d=-A'*n;

abcd=[a b c d];

% Rotation to the xy plane.
% https://math.stackexchange.com/questions/1167717/transform-a-plane-to-the-xy-plane
% Consider correction due to user "paterry".

if (a^2+b^2) > 0 % rotation is required
    den=sqrt(a^2+b^2+c^2);
    cos_theta=c/den;
    sin_theta=sqrt(a^2+b^2)/den;
    
    den2=sqrt(a^2+b^2);
    u1=b/den2;
    u2=-a/den2;
    
    term_cos=(1-cos_theta);
    
    R_row1=[cos_theta+u1^2*term_cos  u1*u2*term_cos u2*sin_theta];
    R_row2=[u1*u2*term_cos cos_theta+u2^2*term_cos  -u1*sin_theta];
    R_row3=[-u2*sin_theta u1*sin_theta cos_theta];
    R10=[R_row1; R_row2; R_row3];
else % no rotation is required
    R10=eye(3);
end

xyv=R10*[XV'; YV'; ZV']; xyv=xyv';

xv=xyv(:,1);
yv=xyv(:,2);
zv=xyv(:,3);

% Check orientation polygon
[orient,signed_area] = polyorient(xv,yv);

if orient == 0
    % if orientation is clockwise, we reflect the 2D polygon w.r.t. the
    % x-axis so to map it into an equivalent one with the right orientation
    RFM=diag([1 -1 1]);
    yv=-yv;
    [orient,signed_area] = polyorient(xv,yv);
else
    RFM=eye(3);
end

if orient == 0
    warning('Some issues with orientation');
end













function [orient,signed_area] = polyorient(x,y)

%POLYORIENT Orientation of polygon
%   Returns the orientation and signed area of a 2D polygon
%
%   Syntax:
%      [ORIENT,SAREA] = POLYORIENT(X,Y)
%
%   Inputs:
%      X, Y   Vectors with the polygon vertices
%
%   Outputs:
%      ORIENT   Polygon orientation. 1 if the orientation is
%               counter-clockwise (direct), 0 otherwise
%      SAREA    Signed area of the polygon, negative if orientation is
%               not direct
%
%   Examples:
%      x1 = [0 0 1 1]; y1 = [1 2 2 1];
%      x2 = [0 0 1 1]; y2 = [1 0 0 1];
%      x3 = [x1 x2];   y3 = [y1 y2];
%
%      [o1,a1] = polyorient(x1,y1) % 0, -1
%      [o2,a2] = polyorient(x2,y2) % 1,  1
%      [o3,a3] = polyorient(x3,y3) % 0,  0
%
%   MMA 21-4-2006, mma@odyle.net

signed_area=0.5* sum(x.*y([2:end 1]) - y.*x([2:end 1]));
orient = signed_area > 0;
