%   FWPolygons.m - test program for 
%   area, centroid, perimeter and area moments of polygons 
% -------------------------------------------------------------------
% originally by H.J. Sommer III - 02.05.14 - tested under MATLAB v5.2
% modified for Matching Food and Wine polygons by
% Stefano De Marchi - Sept. - Dic. 2009 
%--------------------------------------------------------------------
clear all
close all

rt=1;   %Reflexion (rt=0) or Roto-translation (rt=1) 

% constants
d2r = pi / 180;
scrsz = get(0,'ScreenSize');
figure('Position',[0 0 scrsz(4)-150 scrsz(4)-150]);

% plot of 10 concentric circles and lines
% of wine and food characteristics 
axis square
[coeff_w, coeff_f]=cerchi;

%------------------
% first polygon
%------------------
% get wine characteristics from GUI
wine=myInput(1);
for j=1:length(wine)
    cw(j)=str2num(wine{j});
end
ncw=length(coeff_w);
for i=1:ncw
   if cw(i)==0 & i>1
     x(i)=x(i-1); 
     y(i)=y(i-1);  
   else
     nw=coeff_w(i);
     x(i)=cw(i)*cos(nw);
     y(i)=cw(i)*sin(nw);
   end
end
if cw(1)==0,
    x=x(2:ncw);
    y=y(2:ncw);
end
% get geometry
[ geom1, iner1, cpmo1 ] = polygeom( x, y );

% show results
area = geom1(1);
x_cen1 = geom1(2);
y_cen1 = geom1(3);
perimeter = geom1(4);
disp( [ ' ' ] )
disp( [ 'polygon of wine characteristics ' ] )
disp( [ ' ' ] )
disp( [ '      area     x_cen     y_cen     perim' ] )
disp( [ area  x_cen1  y_cen1  perimeter ] )

I1 = cpmo1(1);
angle1 = cpmo1(2);
I2 = cpmo1(3);
angle2 = cpmo1(4);
disp( [ ' ' ] )
disp( [ '        I1        I2' ] )
disp( [ I1 I2 ] )
disp( [ '    angle1    angle2' ] )
disp( [ angle1/d2r angle2/d2r ] )

% plot outline
xplot1 = [ x' ; x(1) ];
yplot1 = [ y' ; y(1) ];

 %area(xplot1,yplot1)
 hold on

% Axes of centroid
rad = 10;
xx1_1 = [ x_cen1-rad*cos(angle1)  x_cen1+rad*cos(angle1) ];
yy1_1 = [ y_cen1-rad*sin(angle1)  y_cen1+rad*sin(angle1) ];
xx2_1 = [ x_cen1-rad*cos(angle2)  x_cen1+rad*cos(angle2) ];
yy2_1 = [ y_cen1-rad*sin(angle2)  y_cen1+rad*sin(angle2) ];

plot(xplot1,yplot1,'b-o', x_cen1,y_cen1,'bo',xx1_1,yy1_1,'k:', ...
     xx2_1,yy2_1,'k:', 'LineWidth',2)

%------------------
% second polygon
%------------------
% get food characteristics from GUI
food=myInput(2);
ncf=length(coeff_f);
for j=1:length(food)
    fc(j)=str2num(food{j});
end
for i=1:ncf
 if fc==0 & i>1
    x1(i)=x1(i-1); 
    y1(i)=y1(i-1);
 else
  nf=coeff_f(i);
  x1(i)=fc(i)*cos(nf);
  y1(i)=fc(i)*sin(nf);
 end
end
if fc(1)==0,
    x1=x1(2:ncf);
    y1=y1(2:ncf);
end
% get geometry
[ geom2,iner2,cpmo2] = polygeom(x1,y1);

% show results
area = geom2(1);
x_cen2 = geom2(2);
y_cen2 = geom2(3);
perimeter = geom2(4);
disp( [ ' ' ] )
disp( [ 'polygon of food characteristics ' ] )
disp( [ ' ' ] )
disp( [ '      area     x_cen     y_cen     perim' ] )
disp( [ area  x_cen2  y_cen2  perimeter ] )

I1 = cpmo2(1);
angle3 = cpmo2(2);
I2 = cpmo2(3);
angle4 = cpmo2(4);
disp( [ ' ' ] )
disp( [ '        I1        I2' ] )
disp( [ I1 I2 ] )
disp( [ '    angle1    angle2' ] )
disp( [ angle3/d2r angle4/d2r ] )

% plot outline
xplot2 = [ x1' ; x1(1) ];
yplot2 = [ y1' ; y1(1) ];

% Axes of centroid
rad = 10;
xx1_2 = [ x_cen2-rad*cos(angle3)  x_cen2+rad*cos(angle3) ];
yy1_2 = [ y_cen2-rad*sin(angle3)  y_cen2+rad*sin(angle3) ];
xx2_2 = [ x_cen2-rad*cos(angle4)  x_cen2+rad*cos(angle4) ];
yy2_2 = [ y_cen2-rad*sin(angle4)  y_cen2+rad*sin(angle4) ];
plot(xplot2,yplot2,'g-.o',x_cen2,y_cen2,'go',...
     xx1_2,yy1_2,'k:',xx2_2,yy2_2,'k:', 'LineWidth',2)
%----------------------------------
if ~ispolycw(xplot1,yplot1) & ~ispolycw(xplot2,yplot2)
 [xplot1,yplot1]=poly2cw(xplot1,yplot1);
 [xplot2,yplot2]=poly2cw(xplot2,yplot2);
end
[xi,yi] = polybool('intersection',xplot1,yplot1,xplot2,yplot2);
disp('area intersezione con polybool prima rotazione')
interArea=polyarea(xi,yi);

if rt==0
   % reflection 
   xy1=[xplot2 -yplot2];
   % translation on the centroid of the first polygon
   [xy1]=rotoTranslation([xy1(:,1) xy1(:,2)], [x_cen1 y_cen1], 0);  
    
else
%----------   roto translation of the second polygon ---------------

delta=pi/192;
% Subdivision angles for wine characteristics
angles=[pi,pi+2*delta,-2*delta,-delta,pi/2-delta/2,pi/2+delta/2];

P1.x=x; 
P1.y=y;
P1.hole=0;

for s=1:length(angles)
[xy1]=rotoTranslation([xplot2 yplot2], [x_cen1 y_cen1], angles(s));

% second polygon
P2.x=xy1(:,1); 
P2.y=xy1(:,2);
P2.hole=0;

P3=PolygonClip(P1,P2,1);
disp('Area of intersection = ') 
interArea1=polyarea(P3(1).x,P3(1).y);
if interArea1>interArea,
    interArea=interArea1
    optimalAngle=angles(s)
end
end

[xy1]=rotoTranslation([xplot2 yplot2], [x_cen1 y_cen1], optimalAngle);

end
% get geometry
[ geom, iner, cpmo ] = polygeom( xy1(:,1), xy1(:,2) );

% show results
 area = geom(1);
 x_cen = geom(2);
 y_cen = geom(3);
 perimeter = geom(4);
% disp( [ ' ' ] )
% disp( [ '3x5 test rectangle with long axis at 30 degrees' ] )
% disp( [ ' ' ] )
% disp( [ '      area     x_cen     y_cen     perim' ] )
% disp( [ area  x_cen  y_cen  perimeter ] )
% 
% I1 = cpmo(1);
 angle1 = cpmo(2);
% I2 = cpmo(3);
 angle2 = cpmo(4);
% disp( [ ' ' ] )
% disp( [ '        I1        I2' ] )
% disp( [ I1 I2 ] )
% disp( [ '    angle1    angle2' ] )
% disp( [ angle1/d2r angle2/d2r ] )

%-----------------------------------------------
% plot outline
xplot = [ xy1(:,1) ; xy1(1,1) ];
yplot = [ xy1(:,2) ; xy1(1,2) ];
rad = 10;
x1 = [ x_cen-rad*cos(angle1)  x_cen+rad*cos(angle1) ];
y1 = [ y_cen-rad*sin(angle1)  y_cen+rad*sin(angle1) ];
x2 = [ x_cen-rad*cos(angle2)  x_cen+rad*cos(angle2) ];
y2 = [ y_cen-rad*sin(angle2)  y_cen+rad*sin(angle2) ];

figure(2)
% plot of the polygons after rotation
plot(xplot1,yplot1,'b-o', x_cen1,y_cen1,'co',...
    xx1_1,yy1_1,'k:',xx2_1,yy2_1,'k:', 'LineWidth',2)
hold on
plot(xplot,yplot,'g-.o', x_cen,y_cen,'co',...
     x1,y1,'k:', x2,y2,'k:', 'LineWidth',2)

grid
axis equal
grid
hold off



%------------------------
figure(3)
% first polygon

P1.x=x; 
P1.y=y;
P1.hole=0;

% second polygon
P2.x=xy1(:,1); 
P2.y=xy1(:,2);
P2.hole=0;

P3=PolygonClip(P1,P2,1);
disp('Area of intersection = ') 
interArea=polyarea(P3(1).x,P3(1).y)

    for i=1:3
            eval(['p=P' num2str(i) ';'])
            tcolor = [0 0 1;  0 0 1; 1 0 0];
            for np=1:length(p)
                obj=patch(p(np).x,p(np).y,i);
                if p(np).hole==1; set(obj,'facecolor','w'); end
            end
    end


