MATLAB Answers

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

88 views (last 30 days)
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.
% polygon. in this example, a circle
theta = 0 : 0.1 : 2*pi;
radius = 2;
x = radius * cos(theta);
y = radius * sin(theta);
% add the first point to the end to close the polygon
x = [x x(1)];
y = [y y(1)];
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);
%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);
% 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


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

Sign in to comment.

Accepted Answer

Kelly Kearney
Kelly Kearney on 23 Mar 2018
Edited: Kelly Kearney 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];
[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]';
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 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 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;
% 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
hold on;
for iv = 1:nvert
plot(xplt(:,50,iv), yplt(:,50,iv));
title(sprintf('theta = %.2f from each vertex', theta(50)));
On my computer, the additional calculations bring things up to about 3 seconds total.

Sign in to comment.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!