## How can I quickly find the intersections between many individual line segments and a polygon?

### Rose Palermo (view profile)

on 22 Mar 2018
Latest activity Commented on by Rose Palermo

on 2 Apr 2018

### Kelly Kearney (view profile)

I've been using polyxpoly and looping through each line segment and the polygon, but when looping through thousands of line segments, this takes hours.
The code below uses a circle as the polygon. For each point around the polygon, a ray is extended in all directions. I check if the ray is in or on the polygon. If it is, a line segment is sent in that direction. I then find the intersection between all of those line segments and the polygon. This loops through for each point on the circle.
It's too slow for my applications now. I don't think I can make the polyxpoly part a vector. If anyone knows how to make this part much faster, I'd appreciate your ideas.
Code:
% polygon. in this example, a circle
theta = 0 : 0.1 : 2*pi;
% add the first point to the end to close the polygon
x = [x x(1)];
y = [y y(1)];
plot(x,y,'k','LineWidth',2)
hold on
% find 1/2 the smallest distance from one point to the next
eps = 0.5*sqrt(min((x(1:end-1)-x(2:end)).^2 + (y(1:end-1)-y(2:end)).^2));
% find 2 * the maximum distance between 2 points
D = 2*max(sqrt(x.^2 + y.^2));
% angles around from point
theta = linspace(0,2*pi,1000);
for j=1:length(x)
dx = D*cos(theta);
dy = D*sin(theta);
% check if a small ray is in or on the polygon
dx_eps = eps * cos(theta);
dy_eps = eps * sin(theta);
%inside of the bounding polygon ( in this case a circle)
[in, on] = inpolygon(x(j)+dx_eps,y(j)+dy_eps,x,y);
ind_inon = find(in|on);
% for the points that are in or on the polygon, find the intersection
% of a line segment in that direction with the polygon
for i=2:length(ind_inon)
% line segment to find intersections for each angle from point j
x1 = [x(j), x(j) + dx(ind_inon(i))];
y1 = [y(j), y(j) + dy(ind_inon(i))];
[xi,yi] = polyxpoly(x1,y1,x,y); % find intersection between ray and each polygon
% find shortest intersection that isn't 0
[m,k] = min((x(j)-xi).^2 + (y(j)-yi).^2);
while m < eps^2
xi(k) = [];
yi(k) = [];
[m,k] = min((x(j)-xi).^2 + (y(j)-yi).^2);
end
%if there is no intersection (lines are parallel and on top of each
%other), then use point j because I want a filler in there
if isempty(xi)
x_int_save(ind_inon(i)) = x(j);
y_int_save(ind_inon(i)) = y(j);
continue
end
% save intersecting points
% point int = points that make up fetch area polygon
x_int = xi(k);
y_int = yi(k);
x_int_save(ind_inon(i)) = x_int;
y_int_save(ind_inon(i)) = y_int;
plot([x(j) x_int],[y(j) y_int]); drawnow % plot all intersecting rays
end
end
end

Jan

on 22 Mar 2018
Rose Palermo

### Rose Palermo (view profile)

on 22 Mar 2018
Thank you. Posted it in the body of my question.
Image Analyst

### Image Analyst (view profile)

on 23 Mar 2018
Note:
To use 'polyxpoly', the following product must be licensed, installed, and enabled:
Mapping Toolbox

### Kelly Kearney (view profile)

on 23 Mar 2018
Edited by Kelly Kearney

### Kelly Kearney (view profile)

on 23 Mar 2018

Is your polygon always convex? If so, you can do the ray extension from each polygon vertex in all directions at once, and assume that all intersection points lead to the desired line segments. This would reduce the number of calls to polyxpoly from #vertices * #rays to just #vertices.
th = linspace(0,2*pi,60);
r = 2;
x = r.*cos(th);
y = r.*sin(th);
nvert = length(x) - 1; % # vertices in polygon, assuming closed
nray = 1000;
theta = linspace(0,2*pi,nray);
D = 2*max(sqrt(x.^2 + y.^2));
xseg = nan(3,nray,nvert);
yseg = nan(3,nray,nvert);
for iv = 1:nvert
xend = x(iv) + D.*cos(theta);
yend = y(iv) + D.*sin(theta);
xseg(1:2,:,iv) = [ones(size(theta))*x(iv); xend];
yseg(1:2,:,iv) = [ones(size(theta))*y(iv); yend];
end
tic;
[xin, yin] = deal(cell(nvert,1));
for iv = 1:nvert
[xi,yi] = polyxpoly(reshape(xseg(:,:,iv),[],1), reshape(yseg(:,:,iv),[],1), x, y);
xin{iv} = [ones(size(xi))*x(iv) xi]';
yin{iv} = [ones(size(xi))*y(iv) yi]';
end
toc;
On my computer, this takes ~1 sec.
The down side is that this method doesn't save the info regarding which ray/angle each intersection point was associated with, so if that info is important to your application, you'll have to do a bit of back-calculation of angles to match them up.
Likewise, if your polygons can be concave, things get more difficult, since you'll have to sort which bits of the resulting line segments are in and outside the polygon. But again, you should be able to calculate that via some post-polyxpoly processing.

Rose Palermo

### Rose Palermo (view profile)

on 26 Mar 2018
My polygons are most often concave, which is why I've been using the minimum (finding the shortest intersection that isn't zero) for each angle. For each intersection I need which ray it was (which theta) and the first intersection (similar to the line of sight of each ray). This is much faster than what I've been doing, so if there's a way to get out that information then this would work really well for me.
Kelly Kearney

### Kelly Kearney (view profile)

on 26 Mar 2018
% An example polygon with lots of concavity
th = linspace(0,2*pi,60);
r = 2 + rand(size(th))-0.5;
x = r.*cos(th);
y = r.*sin(th);
x(end) = x(1); % close it
y(end) = y(1);
nvert = length(x) - 1; % # vertices in polygon, assuming closed
% Set up rays
nray = 1000;
theta = linspace(0,2*pi,nray+1);
theta = theta(1:end-1); % eliminate duplicate at 2pi
D = 2*max(sqrt(x.^2 + y.^2));
% Set up output: light-of-sight intersection point for each vertex/angle
% pair
xlos = nan(nray,nvert);
ylos = nan(nray,nvert);
% Loop over polygon vertices...
for iv = 1:nvert
xend = x(iv) + D.*cos(theta);
yend = y(iv) + D.*sin(theta);
xseg = [ones(size(theta))*x(iv); xend; nan(size(theta))];
yseg = [ones(size(theta))*y(iv); yend; nan(size(theta))];
[xi,yi] = polyxpoly(xseg(:), yseg(:), x, y);
% Unique intersection points
xyi = setdiff(unique([xi yi], 'rows'), [x(iv) y(iv)], 'rows');
xi = xyi(:,1);
yi = xyi(:,2);
% Figure out which points are the first intersection points along each
% ray by grouping by angle
ang = atan2(yi - y(iv), xi - x(iv));
ang = interp1(wrapToPi(theta), theta, ang, 'nearest'); % remove rounding error
d = sqrt((yi - y(iv)).^2 + (xi - x(iv)).^2);
[G, unqang] = findgroups(ang);
dmin = splitapply(@min, d, G);
[~,loc] = ismember([unqang dmin], [ang d], 'rows');
ximin = xi(loc);
yimin = yi(loc);
% Now figure out whether each ray was inside or outside the polygon
% before intersecting. Can check simply by seeing if bisector is
% inside or outside.
dx = ximin - x(iv);
dy = yimin - y(iv);
isin = inpolygon(x(iv) + dx/2, y(iv) + dy/2, x, y);
ximin = ximin(isin);
yimin = yimin(isin);
unqang = unqang(isin);
% Replace in proper slot
[~,loc] = ismember(unqang, theta);
xlos(loc,iv) = ximin;
ylos(loc,iv) = yimin;
end
% For plotting, define each resulting line segment from vertex to
% line-of-sight point
xplt = permute(cat(3, x(1:end-1).*ones(nray,1), xlos), [3 1 2]);
yplt = permute(cat(3, y(1:end-1).*ones(nray,1), ylos), [3 1 2]);
% To check, let's plot just the rays associated with angle idx=50
plot(x,y,'k');
hold on;
for iv = 1:nvert
plot(xplt(:,50,iv), yplt(:,50,iv));
end
title(sprintf('theta = %.2f from each vertex', theta(50)));
On my computer, the additional calculations bring things up to about 3 seconds total.
Rose Palermo

### Rose Palermo (view profile)

on 2 Apr 2018
Thank you so much! This is perfect.