%S=[10 3 6 1 11 20 25 3;2 1 5 6 7 4 8 11]
%S=[22 10 18 -7 -3 20 2 9;-1 5 27 11 26 7 27 13]
%S=[0 0 10 4 4 6 6 10;0 10 0 4 6 4 6 10]
S=[0 0 10 4 4 5 5 10;0 11 0 4 8 4 8 11]

%Given an ordered set of S=(Sx;Sy) (2 by N matrix) of N points in R2 from which we may
%choose any W points, then there are a total of (N choose W) combinations forming
%verticies of a W-gon.

Sx=S(1,:);Sy=S(2,:);   %S(1,:) defines the 1st row vector
N=length(Sx);
W=5;                   %Number of vertices. W=5 for a pentagon.
index=combnk(1:N,W);   %Indexed set of all combinations of W-gon vertices [1:W]
[L,W]=size(index);     %L=(N choose W)

%The next part of the algorithm orders the points defined by S(:,index)
%so that they can be SEQUENTIALLY connected in a pentagon structure.

K=[0 1 0 0 0;
   0 0 1 0 0;
   0 0 0 1 0;
   0 0 0 0 1;
   1 0 0 0 0];         % indexk Transformation matrix [1 2 3 4 5] -> [2 3 4 5 1]
k=1:W;
indexk=K*k';
indexk=indexk';
count1=0;
count2=0;
hold on

for i=1:L
    indexi=index(i,:);
    So=S(:,indexi);    %Choose W points of S according to the index(i,:)
                       %arrangement (2 by W matrix) with first row
                       %x=So(1,:) and y=So(2,:). The colon means the set of
                       %all column elements.
    ro=mean(So');      %takes the mean value of the xy-coordinates. The prime'
                       %is a matrix transpose of So with x-coor=So(:,1),y-coor=So(2,:)
    Ravg=[ro(1)*ones(1,W);ro(2)*ones(1,W)]; %(2 by W matrix)
    dR=So-Ravg;        %set of vector differences wrt ro
    theta=zeros(1,W);

    for j=1:W                       %Find angle theta(j) of dR wrt positive x-direction
        dX=dR(1,j);dY=dR(2,j);
        theta(j)=acos(dX/norm(dR(:,j)))*180/pi;
        if dY<0
            theta(j)=360-theta(j);  %since inverse cosine bijective on [0,180]
        end
    end
    [thetaSort,ithetaSort]=sort(theta);
    SO=So(:,ithetaSort);            %This is our sorted set of W CCW ordered points
    dRR=dR(:,ithetaSort);           %Set of sorted vectors defining the verticies wrt the distribution center Ravg

    %Construct set of CCW orriented vectors defining the edge of the pentagon
    phi=zeros(1,W);

    dr=SO(:,indexk)-SO(:,k);   %Vectors defining pentagons edges. Note SO has already been ordered in CCW direction.
    for m=1:W
        %Acute angle between dr(:,indexm(m)) and dr(:,m) vectors
        DOT1=dot(dr(:,indexk(m)),dr(:,m));
        CROSS1=cross([dr(:,indexk(m));0],[dr(:,m);0]);
        phim=acos(abs(DOT1)/(norm(dr(:,indexk(m)))*norm(dr(:,m))))*180/pi;
        if DOT1>0              %Defines obtuse angles if they exist
            phim=180-phim;
        end
        if CROSS1==0
            Disp('A set three collinear points was found. Please redefine \nthese points so they are not collinear with each other, \nor any other points.')
            count=count+1;
            break
        end
        phi(m)=phim;
    end
    plot(SO(1,:),SO(2,:),'b-o',S(1,:),S(2,:),'o')
    SM=sum(phi);
    if count1==1
        break
    end
    if SM<180*(W-2)
        continue               %Implies one of the interior angles is larger than 180 deg.
    else
        fprintf(1, 'Unfortunately a convex pentagon can be formed by connedcting CCW orriented lines using the points (%2f),(%2f),(%2f),(%2f),(%2f). \n',ithetaSort)
        count2=1;
        break
    end
end
hold off

%To plot a convex W-gon.
if count2==1
    plot(SO(1,:),SO(2,:),'b-o',S(1,:),S(2,:),'o')
end

%If no convex W-gons are found... done!
if i==L
    fprintf('CONGRATULATIONS, your selected set of points has NO convex %2.0f -gons!! \n',W)
end