// Software    : Scilab >= 5.5
// Function    : polygonmetrics()
// License     : CeCILL-C : http://www.cecill.info/licences
// Copyright(C): Samuel GOUGEON - France
// Version     : 1.0
// Date        : 2015-08-18
// Release     : http://fileexchange.scilab.org/toolboxes/445000
// History     :
//   2015-08-18 : First release on FileExchange
//
// DISCLAIMER: This macro is provided as is, with NO WARRANTY of any kind.
//  Use it to your own risks, awareness, and responsability.

function pm = polygonmetrics(XYZ)
//
//Computes the length, area, CoM, and unit vector of a plane or 3D simple polygon
//
// SYNTAX
// ------
//      polygonmetrics     // Displays this help
// pm = polygonmetrics(XY)
// pm = polygonmetrics(XYZ)
//
// PARAMETERS
// ----------
//  XY : Matrix of size (2,N) or (N,2) with N>2:
//       Cartesian coordinates of linked vertices of a simple plane polygon. 
//  XYZ: Matrix of size (3,N) or (N,3) with N>2
//       Cartesian coordinates of linked vertices of a simple polygon. 
//    NOTES: 
//       * If the last point is not equal to the first one, polygonmetrics()
//         closes the polygon.
//       * The polygon must be simple, without self-intersections.
//         In case of a self-crossing polygon, the area value is unreliable.
//       * Concave as well as convex polygons are accepted. 
//       * In 3D, the 'pseudo-polygon' is not necessarily plane.
//
//  pm : Geometrical parameters of the polygon. A structure with the following
//        attributes is returned:
//         pm.nbvert  : Number of vertices
//         pm.length  : Perimeter
//         pm.area    : Positive area
//         pm.com     : 3D (x,y,z) coordinates of the Center of Mass
//         pm.dir     : Resulting unit vector (a,b,c) orthogonal to the polygon
//
// DESCRIPTION
// -----------
// The same algorithm is used in 2D and in 3D, with plane and non-plane polygons:
//   - The position of the Center of Mass CoM of the polygon is computed.
//     In 2D, it is set to the vertices centroid. 
//     In 3D, its coordinates are the average of coordinates of its vertices.
//   - The vectorial area of each triangle going from the CoM to an edge is
//     computed. Its norm is the triangle's area. Its direction is orthogonal 
//     to the triangle. Its sign follows the corkscrew right hand rule.
//   - The total oriented area is the sum of vectorial areas for all triangles
//     set on all edges.
//   - The polygon's area is the norm of the total oriented area.
//   - The polygon direction is the unit vector of the total oriented area.
//
// In 2D, the resulting area matches with results from the Surveyor's formula:
//  area = |sum_(i=1=>N){ x(i+1)*y(i) - x(i)*y(i+1) }| / 2
// The algorithm somewhat generalizes the Surveyor's formula in 3D.
//
// SEE ALSO
// --------
//  mesharea : http://fileexchange.scilab.org/toolboxes/370000
//
// REFERENCE
// ---------
//  Comments, scoring and bug reports are welcome on
//  http://fileexchange.scilab.org/toolboxes/445000
//
// EXAMPLES
// --------
// da = gda(); da.isoview = "on"; da.axes_visible = "on"; clf
// yellow = color("yellow");
//
//// Rectangle dx=2, dY = 2
// x = [0 2 2 0 0]; y = [0 0 1 1 0];
// subplot(4,2,1); xfpoly(x, y, yellow)
// p = polygonmetrics([x ; y])
//
//// Self-intersecting rectangle
// x = [0 2 2 0 0]; y = [0 1 0 1 0];
// subplot(4,2,3); xfpoly(x, y, yellow)
// p = polygonmetrics([x ; y])
//
//// Plane concave polygon: 
// x=[0 3 3 2 2 3 3 0 0 0.5 0.5 1 1 0 0]; y=[0 0 1 1 2 2 3 3 1.5 1.5 2 2 1 1 0];
// subplot(2,2,2); xfpoly(x, y, yellow)
// p = polygonmetrics([x ; y])
//
//// Other plane polygon
// r = grand(1,13,"unf",1,5);
// t = linspace(0, 2*%pi, 13);
// x = r.*cos(t); y = r.*sin(t);
// subplot(2,2,3); xfpoly(x, y, yellow)
// p = polygonmetrics([x ; y ])
//
//// Equilateral triangle in a cube
// x = [0 0 1 0]; y = [0 1 1 0]; z = [0 1 0 0];
// subplot(2,2,4); param3d(x, y, z)
// e = gce(); e.fill_mode = "on"; e.background = yellow;
// ax = gca(); ax.isoview = "on"; ax.rotation_angles = [60 -50];
// p = polygonmetrics([x ; y ; z])
//
// da.isoview = "off"; da.axes_visible="off";
//

    fname = "polygonmetrics"
    if argn(2)==0
        head_comments(fname)
        pm = struct()
        return
    end
    
    // CHECKING PARAMETERS
    // -------------------
    if type(XYZ)~=1
        msg = _("%s: Wrong type for argument %d: Real matrix expected.\n")
        error(msprintf(msg, fname, 1))
    end
    if max(size(XYZ))<3
        msg = _("%s: Wrong size for input argument #%d: (N>2,2) | (N>2,3) matrix expected.\n")
        error(msprintf(msg, fname, 1))
    end
    
    nd = min(size(XYZ))
    // Formating
    if size(XYZ,1)==nd
        XYZ = XYZ'
    end
    
    // PROCESSING
    // ----------
    // Closure
    if or(XYZ($,:)~=XYZ(1,:))
        XYZ = [XYZ ; XYZ(1,:)]
    end
    // Number of vertices
    pm.nbvert = size(XYZ,1)-1
    
    // Perimeter
    dXYZ = XYZ(2:$,:) - XYZ(1:$-1,:)
    pm.length = sum(sqrt(sum(dXYZ.^2, "c")))
    
    // Center of Mass
    if nd==2
        // CoM = centroid of vertices
        xyp = prod(XYZ(2:$,:),"c")
        COM = [sum(xyp)/sum(XYZ(2:$,2)) sum(xyp)/sum(XYZ(2:$,1)) 0]
    else
        // CoM = average position of vertices
        COM = mean(XYZ(2:$,:), "r")
    end

    // Area
    if nd==2
        // pm.area = abs(sum(XYZ(2:$,1).*XYZ(1:$-1,2) - XYZ(1:$-1,1).*XYZ(2:$,2)))/2
        XYZ(:,$+1) = 0
    end
    d = XYZ - ones(XYZ(:,1))*COM
    va = d(1:$-1,:)
    vb = d(2:$,:)
    dS = cross(va', vb')
    oriented_area = sum(dS,"c")
    pm.area = norm(oriented_area)/2
    // absolute_area = sum(sqrt(sum(dS.^2,"r")))
    // pm.aarea = absolute_area/2
    pm.com = COM
    tmp = sum(dS,'c')'
    pm.dir = tmp/norm(tmp)
endfunction

