function [Geo] = Polygons_intersection(S,Display_result,Accuracy)
% This function computes n-times intersection region of shapes collection
% and allows to identify every intersection region in which shapes
% intersect.
%
% The function takes one argument as input, a structure S containing
% geometrical descriptions of shapes, and delivers one output argument,
% a structure containing the different shape intersections, with the
% shape indexes involded in the intersection and the associated area.
% The second argument is optional. Display_result is a binary value which
% indicates whether the result should be displayed (1) or not (0).
%
% Input: S : Structure containing geometrical description of polygons.
% S(i) contains all information relative to the i-th shape
% S(i).P(j) gives access to the geometrical description of the j-th
% element of the i-th shape.
% XData : S(i).P(j).x : Vector
% YData : S(i).P(j).y : Vector
% Hole : S(i).P(j).hole : Binary value (1= hole, 0= fill).
% This binary variable indicates whether the consider polygon
% S(i).P(j) is a hole or not. Hole are represented with
% dotted lines on figures.
% Display_result: Binary value used to display or not result
% Accuracy: Relative accuracy. The accuracy of the algorithm
% will the product of this parameter and the area of the largest
% polygon
%
% Output: Geo : Structure containing geometrical description of the
% intersection polygons, their area, and shape indexes used
% to compute the intersection.
%
% Geo(i).index contains the polygon indexes
% Geo(i).P contains the geometry of the polygon
% P is a structure containing XData &
% YData of polygons.
% XData : Geo(i).P(j).x : Vector
% YData : Geo(i).P(j).y : Vector
% Hole : Geo(i).P(j).hole : Binary value (1
% = hole, 0= fill)
% Geo(i).area contains the area of the i-th shape.
%
% This function starts to evaluate the intersection between all shapes.
% Results are stored in an intersection matrix M, where M(i,j) contains the
% intersection between shape i & j.
% Then, for each shape which intersects other shapes, one lists all
% possible polygons involved in the intersection with the function
% Polygons_intersection_Create_combination. The results is a list of
% possible intersection, each element of this list is tested and each time
% there is intersection, a new shape is created with the attribute of the
% polygons involved in the intersection.
%
% Then, new shapes that are the results of several intersection
% computations are substracted to polygons lower shape indexes.
%
% This function uses Polygon Clipper, a function posted by Sebastian Hlz:
%
% "The Polygon Clipper (based on the gpc-library) is used to performe
% algebraic operations on two polygons.
% Given two arbitrary polygons (which may self intersect, may contain
% holes, may be constructed of several contours) the Polygon Clipper
% is used to calculate the resulting polygon for the operations diff,
% union, AND, XOR.
%
% Credit for the gpc-library goes to Alan Murta."
%
%
%
% Guillaume JACQUENOT
% guillaume at jacquenot at gmail dot com
% 2007_10_08
% 2009_06_16
if nargin==0
% Creates initial geometry for an example.
S(1).P(1).x = [1.5 1.5 0.5 1.5+cos(pi:pi/30:3*pi/2)];
S(1).P(1).y = [0.5 1.5 1.5 1.5+sin(pi:pi/30:3*pi/2)];
S(1).P(1).hole = 0;
S(1).P(2).x = [0.7 0.9 0.9 0.7 0.7]+0.5;
S(1).P(2).y = [0.9 0.9 1.1 1.1 0.9]+0.27;
S(1).P(2).hole = 1;
S(2).P(1).x = [0.4 0.4 1.6 1.6 1.25 0.4];
S(2).P(1).y = [1 0.9 0.725 0.9 1.25 1];
S(2).P(1).hole = 0;
S(3).P(1).x = [0.5 0.75 1.125 1.5 1.125 1.125 0.75 1 0.8 0.5];
S(3).P(1).y = [0.5 0.125 0.125 0.375 1 0.3 0.43 1.125 1.4 0.5]+0.05;
S(3).P(1).hole = 0;
Display_result = 1;
Accuracy = 1e-6;
Polygons_intersection(S,Display_result,Accuracy);
clear S
S(1).P.x = [0 -10.9180 -9.8212 1.0968 0];
S(1).P.y = [0 -1.3406 -10.2735 -8.9329 0];
S(2).P.x = [6 -4.9180 -3.8212 7.0968 6];
S(2).P.y = [0 -1.3406 -10.2735 -8.9329 0];
S(3).P.x = [6 -4.9180 -3.8212 7.0968 6];
S(3).P.y = [5 3.6594 -5.2735 -3.9329 5];
S(4).P.x = [0 -10.9180 -9.8212 1.0968 0];
S(4).P.y = [5 3.6594 -5.2735 -3.9329 5];
S(1).P.hole = 0;
S(2).P.hole = 0;
S(3).P.hole = 0;
S(4).P.hole = 0;
Display_result = 1;
Accuracy = 1e-6;
Polygons_intersection(S,Display_result,Accuracy);
return;
elseif nargin == 1
Display_result = 0;
Accuracy = 1e-6;
elseif nargin == 2
Accuracy = 1e-6;
end
% Number of polygons stored in N_polygons
N_polygons = numel(S);
for i=1:N_polygons
S(i).index = i;
end
% Compute the area of all elements
S = Polygons_intersection_Compute_area(S);
% Compute total accuracy
Accuracy = Accuracy * max([S(:).area]);
% Preallocating memory to store intersection area between polygons.
Res = zeros(N_polygons);
% Initiate index, which will used to store domains
index = 1;
% Storing_index is used to know where the information of the intersection
% polygons i & j is stored in the output structure Geo
Storing_index = zeros(N_polygons);
% Performs a double sweep on all polygons to perform a first polygon
% intersection
% index is used to count the number of intersections
% Geo is used to the output structure containing all intersection results
% Compute a first level intersection
for i = 1:N_polygons-1
for j = i+1:N_polygons
P = PolygonClip(S(i).P,S(j).P,1);
% if the intersection is not empty
if size(P,2)~=0
% for all elements of P
for p = 1:size(P,2)
Res(i,j) =Res(i,j) + (1-2*P(p).hole)*polyarea(P(p).x,P(p).y);
end
Storing_index(i,j)=index;
% Initiates Geo, structure containing all results.
Geo(index).index = [i,j];
Geo(index).P = P;
Geo(index).area = Res(i,j);
index = index +1;
end
end
end
index_first = index -1;
% This function creates a list of all possible overlapping elements
% computed from an intersection area matrix
Index2 = Polygons_intersection_Create_combination(Res);
V = zeros(1,numel(Index2));
for i=1:numel(Index2)
Geo(index_first+i).index = Index2(i).data;
j=1;
str2compare = Index2(i).data(1:end-1);
N_str2compare = numel(str2compare);
bool=0;
while (j<=index_first+i-1) && bool==0
if length(str2compare)==length(Geo(j).index);
if sum(str2compare == Geo(j).index)==N_str2compare
bool=1;
else
j=j+1;
end
else
j=j+1;
end
end
V(i)=j;
end
% One computes possible polygon intersections.
To_delete = [];
for i=1:numel(Index2)
if ~isempty(Geo(V(i)).P)
P=PolygonClip(Geo(V(i)).P,S(Index2(i).data(end)).P,1);
if size(P,2)~=0
area = 0;
for p = 1:size(P,2)
area = area + (1-2*P(p).hole)*polyarea(P(p).x,P(p).y);
Geo(index_first+i).P = P;
Geo(index_first+i).area = area;
end
else
To_delete = [To_delete,index_first+i];
end
else
To_delete = [To_delete,index_first+i];
end
end
Geo(To_delete)=[];
% Posttreatment
N_Geo = numel(Geo);
for i=1:N_Geo
for j=1:numel(Geo(i).P)
% Polygons are closed
if ((Geo(i).P(j).x(end)~=Geo(i).P(j).x(1)) ||...
(Geo(i).P(j).y(end)~=Geo(i).P(j).y(1)))
Geo(i).P(j).x = [Geo(i).P(j).x;Geo(i).P(j).x(1)];
Geo(i).P(j).y = [Geo(i).P(j).y;Geo(i).P(j).y(1)];
end
end
end
% N_Geo = numel(Geo);
% colour = rand(N_Geo,3);
% figure, hold on, grid on
% for i=1:N_Geo
% for j=1:numel(Geo(i).P)
% h = plot(Geo(i).P(j).x,Geo(i).P(j).y,'color',colour(i,:),'LineWidth',3);
% if Geo(i).P(j).hole==1
% set(h,'LineStyle','--');
% end
% % text(S(n).x(1)+0.2,S(n).y(1)+0.2,num2str(n),...
% % 'color',colour(n,:),'FontSize',14);
% end
% end
% Now that the all intersection polygons have been computed, one calls the
% following function that will splits the original polygons into smaller
% ones
Geo = Polygons_intersection_Posttreatment(S,Geo,Display_result,Accuracy);
