function [wam_pts, err]=polygon_wam(polygon_vertices,N)

%--------------------------------------------------------------------------
% OBJECT.
%----------
%
% THIS ROUTINE COMPUTES THE NODES "wam_pts" AND OF A POLYGON. THE POLYGON 
% IS SUBDIVIDED IN TRIANGLES BY A MINIMAL TRIANGULATION (see [1]) AND THEN, 
% IN EACH TRIANGLE "T" A WAM OF THE UNIT SIMPLEX IS MAPPED ON "T".
%
% PS. WE USED A BASIC TRIANGULATION ALGORITHM FOR THIS PURPOSE. A FASTER
% ONE HAS BEEN STUDIED BY CHAZELLE AND SEIDEL. OUR ALGORITHM WORKS IN
% POLYGONS THAT ARE "SIMPLE", I.E. WITHOUT SELF-INTERSECTIONS.
%
%--------------------------------------------------------------------------
% INPUT.
%----------
%
% "polygon_vertices": VERTICES OF THE POLYGON (WITHOUT REPEATING THE 1ST
%    INDEX AS LAST INDEX (TO CLOSE THE POLYGON)). IT IS AN "N x 2" MATRIX.
%
% "N": WAM DEGREE.
%
%--------------------------------------------------------------------------
% OUTPUT.
%----------
%
% "wam_pts": COLUMN VECTORS OF WAM POINTS.
%
% "err": IF "err=1" THEN THE ROUTINE IS NOT BEEN COMPLETED AND A WARNING IS
%    WRITTEN ON THE DISPLAY.
%
%--------------------------------------------------------------------------
% FUNCTIONS USED IN THIS ROUTINE.
%---------------------------------
%
% 1. minimal_triangulation
% 2. ref_wam
% 3. triangle_map
% 4. find_ear (USED BY "minimal_triangulation")
% 5. polyarea (MATLAB BUILT_IN FUNCTION, USED BY "minimal_triangulation").
% 6. setdiff (MATLAB BUILT_IN FUNCTION, USED BY "minimal_triangulation").
% 7. verLessThan (MATLAB BUILT_IN FUNCTION, USED BY "find_ear").
% 8. inpolygon (MATLAB BUILT_IN FUNCTION, USED BY "find_ear").
%
% THE FUNCTIONS "minimal_triangulation", "ref_wam", "triangle_map",
% "find_ear" ARE ATTACHED TO THIS FILE.
%
% THE FUNCTIONS "polyarea", "setdiff", "verLessThan", "inpolygon" ARE
% MATLAB BUILT-IN.
%
%--------------------------------------------------------------------------
% TEST.
%-------
%
% THIS ROUTINE HAS BEEN TESTED IN MATLAB 7.6.0.324 (R2008a).
%
%--------------------------------------------------------------------------
% REFERENCES.
%-------------
%
% [1] Federica Basaglia, Un nuovo metodo di cubatura su poligoni (in
%     italian)
%
% [2] Chazelle and J. Incerpi, Triangulation and shape complexity,
%     ACM Trans. on Graphics, vol. 3, pp. 135-152, 1984.
%
% [3] B. Chazelle, Triangulating a simple polygon in linear time,
%     Discrete Comput. Geom., vol. 6, pp. 485-524, 1991.
%
% [4] David Eberly, Triangulation by Ear Clipping,
%           http://www.geometrictools.com/
%
% [5] R. Seidel, A simple and fast incremental randomized algorithm for
%     computing trapezoidal decompositions and for triangulating polygons,
%     Computational Geometry: Theory and Applications, vol. 1, no.1,
%     pp. 51-64, 1991.
%
% [6] Subhash Suri, Polygon Triangulation
%
% [7] http://cgm.cs.mcgill.ca/~godfried/teaching/cg-projects/97/
%            Ian/cutting_ears.html
%
%--------------------------------------------------------------------------

%--------------------------------------------------------------------------
%% Copyright (C) 2007-2009 Alvise Sommariva, Marco Vianello.
%%
%% This program is free software; you can redistribute it and/or modify
%% it under the terms of the GNU General Public License as published by
%% the Free Software Foundation; either version 2 of the License, or
%% (at your option) any later version.
%%
%% This program is distributed in the hope that it will be useful,
%% but WITHOUT ANY WARRANTY; without even the implied warranty of
%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%% GNU General Public License for more details.
%%
%% You should have received a copy of the GNU General Public License
%% along with this program; if not, write to the Free Software
%% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
%%
%% Author:  Federica Basaglia
%%          Alvise Sommariva <alvise@euler.math.unipd.it>
%%          Marco Vianello   <marcov@euler.math.unipd.it>
%%
%% Date: November 17, 2009
%--------------------------------------------------------------------------

% COMPUTE A MINIMAL TRIANGULATION. THE MATRIX tri HAS POSITIVE
% INTEGERS AS COMPONENTS.
[tri,err]=minimal_triangulation(polygon_vertices(:,1),polygon_vertices(:,2));

ref_pts=ref_wam_pts(N);
wam_pts=[];

for index_triangulation=1:size(tri,1)

    % COMPUTING FOR ANY TRIANGLE A CUBATURE RULE.
    % loc_triangle_pts IS THE LOCAL TRIANGLE THAT IS CONSIDERED.
    % A CUBATURE FORMULA IS COMPUTED ON THE UNIT TRIANGLE AND ADAPTED BY
    % triangle_map AND polyarea TO THE LOCAL TRIANGLE.

    loc_tri_indices=tri(index_triangulation,:);
    loc_polygon_vertices=polygon_vertices(loc_tri_indices,:);
    loc_triangle_pts=triangle_map(loc_polygon_vertices,ref_pts);

    wam_pts=[wam_pts; loc_triangle_pts];

end


%---------------------------------------
% minimal_triangulation.
%---------------------------------------

function [tri,err]=minimal_triangulation(x,y)

%--------------------------------------------------------------------------
% OBJECT.
%----------
%
% COMPUTE MINIMAL TRIANGULATION (BASIC ALGORITHM) ON SIMPLE DOMAINS,
% I.E. WITHOUT SELF-INTERSECTIONS.
%--------------------------------------------------------------------------
% INPUT.
%----------
%
% x,y: VERTICES OF THE POLYGONS. THE FIRST VERTEX IS NOT REPEATED.
%      THEY MUST BE COLUMN VECTORS.
%
%--------------------------------------------------------------------------
% OUTPUT.
%----------
%
% tri: MATRIX IN WHICH EACH ROW IS FORMED BY 3 POSITIVE INTEGERS IN WHICH
%      EACH ONE REPRESENTS A POINTER TO A VERTEX. EXAMPLE: [1 2 3; 1 3 4]
%      DESCRIBE TWO TRIANGLES T1=[x(1) y(1); x(2) y(2); x(3) y(3)] AND
%      T2=[x(1) y(1); x(3) y(3); x(4) y(4)].
%
%--------------------------------------------------------------------------

err=0;

number_pts=length(x);
pts=[x y];
pts_remaining=pts;
pointer_pts_remaining=(1:number_pts);
tri=[];

while length(pts_remaining(:,1)) > 3

    [tri_loc,err]=find_ear(pts_remaining);

    if (err == 1)
        tri=[];
        fprintf('\n \t [WARNING]: THE FUNCTION IS NOT ABLE');
        fprintf(' TO COMPUTE THE TRIANGULATION');
        return;
    end

    tri_add=pointer_pts_remaining(tri_loc);
    tri=[tri; tri_add];
    pointer_pt_remove=tri_add(2);
    
pointer_pts_remaining=setdiff(pointer_pts_remaining,pointer_pt_remove);
    pts_remaining=pts(pointer_pts_remaining,:);

end

tri=[tri; pointer_pts_remaining];




%--------------------------------
% find ear.
%--------------------------------
function [tri,err]=find_ear(vertices_pts)

%--------------------------------------------------------------------------
% OBJECT.
%----------
%
% FIND ONE EAR IN A POLYGON DESCRIBED BY THE ORDER SEQUENCE vertices_pts,
% BY A BASIC ALGORITHM, TO COMPUTE A MINIMAL TRIANGULATION (SEE CHAZELLE -
% SEIDEL ALGORITHM FOR A FASTER IMPLEMENTATION). IF THE ROUTINE IS NOT
% COMPLETED, THEN err=1.
%
% PS. WE USED A BASIC TRIANGULATION ALGORITHM FOR THIS PURPOSE. A FASTER
% ONE HAS BEEN STUDIED BY CHAZELLE AND SEIDEL. OUR ALGORITHM WORKS IN
% POLYGONS THAT ARE "SIMPLE", I.E. WITHOUT SELF-INTERSECTIONS. SEE THE
% REFERENCES ABOVE FOR MORE DETAILS.
%
%--------------------------------------------------------------------------


number_pts=size(vertices_pts,1);
vertices_index=[ (1:number_pts)'; 1];
pts_remaining=vertices_pts;
curr_starting_index=1;
ear_found=0;
err=0;

while (ear_found == 0) & (curr_starting_index < number_pts)

    trial_triangle_vertices_pointer=...
        vertices_index(curr_starting_index:curr_starting_index+2);

    % TRIANGLE TO TEST.
    pts_trial_triangle=vertices_pts(trial_triangle_vertices_pointer,:);

    pts_remaining_pointer=...
        setdiff(1:number_pts,curr_starting_index:curr_starting_index+2);

    pts_remaining=vertices_pts(pts_remaining_pointer,:);

    xv=[vertices_pts(trial_triangle_vertices_pointer,1); ...
        vertices_pts(curr_starting_index,1)];
    yv=[vertices_pts(trial_triangle_vertices_pointer,2); ...
        vertices_pts(curr_starting_index,2)];

    % CHECKING POSITION OF BARYCENTER (TO SEE IF IT IS AN EAR).
    barycenter_triangle=...
        [sum(vertices_pts(trial_triangle_vertices_pointer,1)) ...
        sum(vertices_pts(trial_triangle_vertices_pointer,2))]/3;

    xt=pts_remaining(:,1);
    yt=pts_remaining(:,2);

    in = inpolygon(xt,yt,xv,yv);

    barycenter_triangle=[sum(vertices_pts...
        (trial_triangle_vertices_pointer,1)) ...
        sum(vertices_pts(trial_triangle_vertices_pointer,2))]/3;


    % inpolygon has been modified by MATLAB. We try to patch this problem.

    local_version=verLessThan('matlab', '7.6.1');

    if local_version == 0

        
%******************************************************************
        %
        % IN = inpolygon(X,Y,xv,yv) returns a matrix IN the
        %      same size as X and Y.
        %      Each element of IN is assigned the value 1 or 0 depending on
        %      whether the point (X(p,q),Y(p,q)) is inside the polygonal
        %     region whose vertices are specified by the vectors xv and yv.
        %
        % In particular:
        %
        % IN(p,q) = 1: if (X(p,q),Y(p,q)) is inside the polygonal region
        %              or on the polygon boundary
        % IN(p,q) = 0: if (X(p,q),Y(p,q)) is outside the polygonal region
        %
        % [IN ON] = inpolygon(X,Y,xv,yv) returns a second matrix ON the
        %           same size as X and Y.
        %           Each element of ON is assigned the value 1 or 0
        %           depending on whether the point (X(p,q),Y(p,q)) is on
        %           the boundary of the polygonal region whose vertices are
        %           specified by the vectors xv and yv.
        %
        % In particular:
        %
        % ON(p,q) = 1: if (X(p,q),Y(p,q)) is on the polygon boundary
        % ON(p,q) = 0: if (X(p,q),Y(p,q)) is inside or outside the polygon
        %              boundary.
        %
        
%******************************************************************

        % Working in Matlab 7.6.
        [bar_in,bar_on]= inpolygon(barycenter_triangle(1),...
            barycenter_triangle(2),vertices_pts(vertices_index,1),...
            vertices_pts(vertices_index,2));

        if (sum(in(1:end-1)) == 0) & (bar_in == 1) & (bar_on == 0)
            ear_found=1;
            tri=curr_starting_index:curr_starting_index+2;
            return;
        end

    else

        
%******************************************************************
        %
        % IN = INPOLYGON(X, Y, XV, YV) returns a matrix IN the size of
        % X and Y.  IN(p,q) = 1 if the point (X(p,q), Y(p,q)) is
        % strictly inside the polygonal region whose vertices are
        % specified by the vectors XV and YV;  IN(p,q) is 0.5 if
        % the point is on the polygon; otherwise IN(p,q) = 0.
        %
        
%******************************************************************

        % Working in Matlab 6.1.
        [bar_in]= inpolygon(barycenter_triangle(1),...
            barycenter_triangle(2),vertices_pts(vertices_index,1),...
            vertices_pts(vertices_index,2));

        if (sum(in(1:end-1)) == 0) & (bar_in > 0.5)
            ear_found=1;
            tri=curr_starting_index:curr_starting_index+2;
            return;
        end

    end

    curr_starting_index=curr_starting_index+1;
end

err=1;




%---------------------------------------
% triangle_map.
%---------------------------------------

function pts_out=triangle_map(triangle_vertices,pts_in)

%--------------------------------------------------------------------------
% OBJECT.
%----------
%
% THIS FUNCTION MAPS POINTS pts_in OF THE UNIT TRIANGLE IN THE TRIANGLE
% WITH (ORDERED) VERTICES triangle_vertices, OBTAINING THE POINTS pts_out.
%
%--------------------------------------------------------------------------
% INPUTS.
%----------
%
% triangle_vertices: vertices of a triangle.
%
% pts_in: points on the reference triangle (0,0),(1,0),(0,1).
%
%--------------------------------------------------------------------------
% OUTPUT.
%----------
%
% pts_out: the points "pts_in" are mapped on the triangle
%          "triangle_vertices" by a classical map.
%
%--------------------------------------------------------------------------

Q1=triangle_vertices(1,:);
Q2=triangle_vertices(2,:);
Q3=triangle_vertices(3,:);

b=Q1;
P=[(Q2-Q1)' (Q3-Q1)'];

pts_out_pre=(P*pts_in')';
pts_out(:,1)=b(1)+pts_out_pre(:,1);
pts_out(:,2)=b(2)+pts_out_pre(:,2);



%---------------------------------------
% ref_wam_pts.
%---------------------------------------
function pts=ref_wam_pts(M)

%--------------------------------------------------------------------------
% OBJECT.
%----------
%
% COMPUTE FOR A FIXED DEGREE M, A WAM ON THE UNIT SIMPLEX.
%
%--------------------------------------------------------------------------
% INPUTS.
%----------
%
% M: WAM degree.
%
%--------------------------------------------------------------------------
% OUTPUT.
%----------
%
% pts_in: points on the reference triangle (0,0),(1,0),(0,1).
%
%--------------------------------------------------------------------------

N=M;
deg=M;

triangle=[0 1; 1 0; 0 0];

n1=2*deg;
n2=n1;
j1=(0:1:deg-1);
j2=(0:1:deg);
[rho,theta]=meshgrid(cos(j1*pi/n1),j2*pi/n2);

% Mapping to the simplex.

B=[rho(:).*cos(theta(:)) rho(:).*sin(theta(:))];
meshS=[B(:,1).^2 B(:,2).^2];
meshT(:,1)=triangle(3,1)*(1-meshS(:,1)-meshS(:,2))+...
    triangle(2,1)*meshS(:,2)+triangle(1,1)*meshS(:,1);
meshT(:,2)=triangle(3,2)*(1-meshS(:,1)-meshS(:,2))+...
    triangle(2,2)*meshS(:,2)+triangle(1,2)*meshS(:,1);
nodes_x=meshT(:,1);
nodes_y=meshT(:,2);

nodes_x=[nodes_x; 0];
nodes_y=[nodes_y; 0];

pts=[nodes_x nodes_y];




