Asked by Rose Palermo
on 22 Mar 2018

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;

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)];

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

Answer by Kelly Kearney
on 23 Mar 2018

Edited by Kelly Kearney
on 23 Mar 2018

Accepted Answer

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
on 26 Mar 2018

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;

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
on 2 Apr 2018

Thank you so much! This is perfect.

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 3 Comments

## Jan (view profile)

## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/390016-how-can-i-quickly-find-the-intersections-between-many-individual-line-segments-and-a-polygon#comment_548463

## Rose Palermo (view profile)

## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/390016-how-can-i-quickly-find-the-intersections-between-many-individual-line-segments-and-a-polygon#comment_548487

## Image Analyst (view profile)

## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/390016-how-can-i-quickly-find-the-intersections-between-many-individual-line-segments-and-a-polygon#comment_548858

Sign in to comment.