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

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

[Jan]

Please post your code.

[Rose Palermo]

Thank you. Posted it in the body of my question.

[Image Analyst]

Note:

To use 'polyxpoly', the following product must be licensed, installed, and enabled:

Mapping Toolbox

Answer 1

[Kelly Kearney]

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]

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]

% 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]

Thank you so much! This is perfect.