%% LOAD DATA AND DRAW FIRST FIGURE

% clear both figures
clf(figure(1))
clf(figure(2))

x = [581.04 475.18 872.73 670.43 964.48 978.59 1072.7 1145.6 1503.2 ...
    1439.7 1376.1 1164.4 1150.3 886.85 893.9 868.03 578.68].';
y = [1069.7 956.75 455.69 286.32 323.95 265.14 248.68 248.68 269.85 ...
    293.37 302.78 495.68 490.98 730.92 770.91 829.72 1067.3].';

% first figure: draw countour
figure(1)
plot(x,y,'ro','Markersize',8,'Linewidth',2)
hold on
plot(x, y, 'b-','Linew',1.5);
grid

%% ANALIZE EACH TRIANGLE

% Here each triangle is analized on its own. The analysis consists in
% setting two points on each side of the triangle. The points are set very
% close to the corners, depending on a value "delta". Once the points are
% set it is checked wether every corner lies inside the polygon, using the
% function "inpolygon()". It is necessary to set the points at a certain
% "delta" from the corner, because every corner lies inside the polygon,
% even if the side does not.

% create set of triangles and permute them to plot one at the time
tri = delaunay(x,y);
tri = permute(tri,[3,2,1]);

for i = 1:length(tri)
    
    % plot one triangle at the time
    htr(i) = triplot(tri(1,:,i), x, y,'g');
    
    % delta: distance of points from the corners
    delta = 0.001;
    
    % slope of the sides -- for vertical slopes: Inf
    m = (htr(i).YData(2:3:end) - htr(i).YData(1:3:end))...
        ./(htr(i).XData(2:3:end) - htr(i).XData(1:3:end));
    
    % x-sign of the points: sais wether point has to be on the positive or on
    % the negative side of the corner on the x axis
    signx = sign(htr(i).XData(2:3:end) - htr(i).XData(1:3:end));
   
    % x-values of the points
    xq = [htr(i).XData(1:3:end) + signx(1:end)*delta...
        htr(i).XData(2:3:end) - signx(1:end)*delta];
    
    % y-sign of the points -- for vertical slopes: 0
    signy = sign(xq - [htr(i).XData(1:3:end) htr(i).XData(2:3:end)]);
    % sign for vertical slopes is computed here
    signy_temp = [sign(htr(i).YData(2:3:end) - htr(i).YData(1:3:end))...
        sign(htr(i).YData(1:3:end) - htr(i).YData(2:3:end))].*(signy==0);
    signy = signy + signy_temp;
    
    % y-values of the points -- for vertical slopes: Nan
    yq = [htr(i).YData(1:3:end) htr(i).YData(2:3:end)]...
        + [m m]*delta.*signy;
    % y-values for vertical slopes is computed here
    yq_temp = (isnan(abs(yq))|(abs(yq) == inf)).*...
        ([htr(i).YData(1:3:end) htr(i).YData(2:3:end)] + ...
        signy*delta);
    % infinite or nan values are set to 0 to be substituted
    yq(isnan(yq)|(abs(yq) == inf)) = 0;
    yq = yq + yq_temp;
    
    % logical array: says for every triangle wether all points are inside
    % or outside
    in(i) = logical(prod((inpolygon(xq,yq,x,y))));
    
end

%% PLOT ONLY INSIDE TRIANGLES

% outside triangles are removed using logical array "in"
tri = permute(tri(:,:,in),[3,2,1]);

% second figure: no outside triangles
figure(2)
hold on
htr = triplot(tri, x, y,'g');
plot(x,y,'LineWidth',2) % polygon
axis equal