Some modifications based on this FEX
x = cumsum(rand(1,10));
% Generate random P1 and P2
x = x/max(x);
y = rand(size(x));
P1 = [x; y];
y = cumsum(rand(1,10));
y = x/max(x);
x = rand(size(y));
P2 = [x; y];
X = poly2poly(P1, P2);
% Check with graphical output
figure
hold on
plot(P1(1,:),P1(2,:))
plot(P2(1,:),P2(2,:))
for k=1:size(X,2)
plot(X(1,k),X(2,k),'o');
text(X(1,k),X(2,k),num2str(k));
end
function X = poly2poly(P1, P2)
% function X = poly2poly(P1, P2)
% Intersection of two 2D polygons P1 and P2.
%
% INPUTS:
% P1 and P2 are two-row arrays, each column is a vertice
% They might or might not be wrapped around
% OUTPUT:
% X is two-row array, each column is an intersecting point
%
% Author: Bruno Luong <brunoluong@yahoo.com>
% History:
% Original 20-May-2010
% Wrap around: Pad the first point to the end if necessary
% if ~isequal(P1(:,1),P1(:,end))
% P1 = [P1 P1(:,1)];
% end
% if ~isequal(P2(:,1),P2(:,end))
% P2 = [P2 P2(:,1)];
% end
% swap P1 P2 so that we loop on a smaller one
if size(P1,2) > size(P2,2)
[P1 P2] = deal(P2, P1);
end
% We increment the intermediate results by this amount
increment = 10;
% Empty buffer
X = zeros(2,0);
filled = 0;
sizec = 0;
% Loop over segments of P1
for n=2:size(P1,2)
cn = seg2poly(P1(:,n-1:n), P2);
m = size(cn,2);
filled = filled+m;
% Buffer too small
if sizec < filled
sizec = filled+increment;
X(2,sizec) = 0;
end
% Store the result
X(:,filled+(-m+1:0)) = cn;
end
% remove the tail
X(:,filled+1:end) = [];
end % poly2poly
%%
function X = seg2poly(s1, P)
% function X = seg2poly(s1, P)
% Check if a line segment s1 intersects with a polygon P.
% INPUTS:
% s is (2 x 2) where
% s(:,1) is the first point
% s(:,2) is the the second point of the segment.
% P is (2 x n) array, each column is a vertices
% OUTPUT
% X is (2 x m) array, each column is an intersecting point
%
% Author: Bruno Luong <brunoluong@yahoo.com>
% History:
% Original 20-May-2010
% Translate so that first point is origin
a = s1(:,1);
M = bsxfun(@minus, P, a);
b = s1(:,2)-a;
% Check if the points are on the left/right side
x = [b(2) -b(1)]*M;
sx = sign(x);
% x -coordinates has opposite signs
ind = sx(1:end-1).*sx(2:end) <= 0;
if any(ind)
ind = find(ind);
% cross point to the y-axis (along the segment)
x1 = x(ind);
x2 = x(ind+1);
d = b.'/(b(1)^2+b(2)^2);
y1 = d*M(:,ind);
y2 = d*M(:,ind+1);
dx = x2-x1;
% We won't bother with the degenerate case of dx=0 and x1=0
y = (y1.*x2-y2.*x1)./dx;
% Check if the cross point is inside the segment
ind = y>=0 & y<1;
if any(ind)
y = sort(y(ind));
X = bsxfun(@plus, a, b*y);
else
X = zeros(2,0);
end
else
X = zeros(2,0);
end
end % seg2poly
