
function timing=demo_incurvpolygon(example,Npts,spline_parms,Nsub,...
    plot_sets)

% OBJECT:
% Demo on "incurvpolygon" routine.
%
% INPUT:
% example: parameter that describes the example.
%          1 unit disk
%          2 cardioid case 
%          3 astroid 
%          4 deltoid 
%          5 epitrocoid
% Npts: number of points to test.
% spline_parms: this variable is useful only when X,Y are vectors and can
%      be dropped otherwise, for instance setting "spline_parms=[]" or
%      without passing its parameter to the function "incurvpolygon_v8c";
%      spline_parms is a vector whose i-th row is of the form
%                        [order(i) index(i)]
%      In this case the i-th arc will be defined by the splines
%      "Sx", "Sy" of  order "order(i)" so that the curve
%                    "((Sx(i))(t),((Sy(i))(t))"
%      with "t" in "[t(i),t(i+1)]" will interpolate all the points
%                         (X(t(k)),Y(t((k)))
%      where "k" ranges in the interval [index(i-1), index(i)] (where as
%      default index(0)=1).
% Nsub: number of subdvisions of each subinterval determined by two
%       successive breaks and x-turning points. If not sure of this
%       parameter, set it as "Nsub=[]".
% plot_sets: 1 for plotting domain, pointsets and singular points (black
%     circle)
% 
% OUTPUT
% timing: vector of indomain cputimes (each component is the cputime of an
%     experiment).
%
% ADDITIONAL DATA:
%
% built: april 2019
% modified: May 10, 2019
%
%
% CODE VERSION:
%
% v: 11.0.0 (minor issues)

% ............................. initialize ................................

if nargin < 1, example=5; end
if nargin < 2, Npts=1000; end
if nargin < 3, spline_parms=[4 20]; end
% if nargin < 4, Nsub=ceil(210/spline_parms(end,2)); end 
if nargin < 4, Nsub=4; end 
if nargin < 5, plot_sets=1; end % plot point sets.

if isempty(Nsub) == 1, Nsub=ceil(210/spline_parms(end,2)); end

write_stats=1;
SPLtypestring='not-a-knot';
Ntest=10;
a=0; b=2*pi; % spline full parametric interval

% ......................... define problem ................................

[xt,yt]=define_functions(example);

% spline analysis
N=spline_parms(end,2);

t=linspace(a,b,N); % data to interpolate
t=t';
X=feval(xt,t);
Y=feval(yt,t);

% ......................... defining pointset .............................

xm=min(X); xM=max(X); ym=min(Y); yM=max(Y);
xm=min(0.8*xm,1.2*xm); xM=max(0.8*xM,1.2*xM);
ym=min(0.8*ym,1.2*ym); yM=max(0.8*yM,1.2*yM);
Px=(xm+xM)/2+(xM-xm)*(2*rand(Npts,1)-1)/2;
Py=(ym+yM)/2+(yM-ym)*(2*rand(Npts,1)-1)/2;
P=[Px Py];

% ......................... incurvpolygon .................................

for ii=1:Ntest
    tic
    [in,on,Sx,Sy,boxVx,turning_points_Sx,turning_points_Sy,...
        singular_junctions_Sx]=incurvpolygon(P,X,Y,Nsub,spline_parms,...
        SPLtypestring);
    timing(ii)=toc;
end

% ......................... define singular points ........................

if isempty(turning_points_Sx) == 0
    singular_points=[turning_points_Sx(:,1:2)];
else
    singular_points=[];
end

if isempty(singular_junctions_Sx) == 0
    singular_points=[singular_points; singular_junctions_Sx(:,1:2)];
end

% ......................... statistics ....................................

if write_stats == 1
    fprintf('\n \t POINTS ANALYSED: %6.0f',size(P,1));
    for ii=1:size(spline_parms,1)
        fprintf('\n \t ORDER: %4.0f LAST PT: %7.0f',spline_parms(ii,1),...
            spline_parms(ii,2));
    end
    fprintf('\n \t AVER. CPUTIME: %1.3e \n \n',timing/Ntest);
end

% ............................ plots ......................................

if plot_sets == 1
    plot_spl_boxes(Sx,Sy,boxVx,P,in,singular_points);
end





%--------------------------------------------------------------------------
% Attached functions.
%--------------------------------------------------------------------------

function plot_spl_boxes(Sx,Sy,boxV,P,res,singular_points)

% ........................ plot preferences ...............................
plot_boxes=1;
plot_singpts=0;
plot_axis=1;

% ........................... plot: start .................................
clf;
hold on;

if plot_axis == 0
    axis off;
end

% ........................... plot: boxes .................................
if plot_boxes == 1
    for ii=1:size(boxV,1)
        
        boxL=boxV(ii,:);
        
        xminL=boxL(1); xmaxL=boxL(2);
        yminL=boxL(3); ymaxL=boxL(4);
        
        if boxV(ii,5) == 0
            plot([xminL xmaxL xmaxL xminL xminL],...
                [yminL yminL ymaxL ymaxL yminL],'k-');
        else
            fill([xminL xmaxL xmaxL xminL xminL],...
                [yminL yminL ymaxL ymaxL yminL],[0.85 0.85 0.85]);
            plot([xminL xmaxL xmaxL xminL xminL],...
                [yminL yminL ymaxL ymaxL yminL],'k-',...
                'LineWidth',0.5,'Color',[0.5 0.5 0.5]);
        end
        
    end
end

% ........................... plot: all points ............................

plot(P(:,1),P(:,2),'ro','LineWidth',2,...
    'MarkerEdgeColor','r',...
    'MarkerFaceColor','r',...
    'MarkerSize',2);

% ........................ plot: points in domain .........................

inpts=find(res == 1);

plot(P(inpts,1),P(inpts,2),'go','LineWidth',4,...
    'MarkerEdgeColor','c',...
    'MarkerFaceColor','c',...
    'MarkerSize',3);

% ........................ plot: singular points ..........................

if plot_singpts == 1
    if size(singular_points,1) > 0
        plot(singular_points(:,1),singular_points(:,2),'ko',...
            'MarkerSize',8);
    end
end

% ....................... plot: parametric splines ........................

L=length(Sx);
for ii=1:L
    
    SxL=Sx(ii);
    SyL=Sy(ii);
    
    SxL_breaks=SxL.breaks;
    
    
    Nbreaks=length(SxL_breaks);
    
    for kk=1:Nbreaks-1
        
        t0=SxL_breaks(kk); t1=SxL_breaks(kk+1);
        
        N=10^3;
        tt=linspace(t0,t1,N);
        
        xx=ppval(SxL,tt);
        yy=ppval(SyL,tt);
        
        
        plot(xx,yy,'k-','LineWidth',3);
    end
    
end

hold off;




function [xt,yt]=define_functions(example)

switch example
    case 1 % unit disk
        xt=@(t) cos(t);
        yt=@(t) sin(t);
    case 2 % cardioid
        xt=@(t) 1-2*cos(t)+cos(2*t);
        yt=@(t) 2*sin(t)-sin(2*t);
    case 3 % astroid
        xt=@(t) 3*cos(t)+cos(-3*t);
        yt=@(t) 3*sin(t)+sin(-3*t);
    case 4 % deltoid
        xt=@(t) (2/3)*cos(t)-(1/3)*cos(2*t);
        yt=@(t) (2/3)*sin(t)+(1/3)*sin(2*t);
    case 5 % epitrocoid
        r=2;
        R=16;
        d=1;
        xt=@(t) (R+r)*cos(t)-d*cos((R+r)*t/r);
        yt=@(t) (R+r)*sin(t)-d*sin((R+r)*t/r);
    case 6
        xt=@(x) 1+2*x+3*x.^2+4*x.^3;
        yt=@(x) 1+2*x+3*x.^2+4*x.^3;
end







