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

19 次查看(过去 30 天)
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
  3 个评论

请先登录,再进行评论。

采纳的回答

Kelly Kearney
Kelly Kearney 2018-3-23
编辑:Kelly Kearney 2018-3-23
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.
  3 个评论
Kelly Kearney
Kelly Kearney 2018-3-26
% 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.

请先登录,再进行评论。

更多回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Elementary Polygons 的更多信息

Community Treasure Hunt

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

Start Hunting!

Translated by