Intersection points between a line and a polygon.
P = intersectLinePolygon(LINE, POLY)
Returns the intersection points of the lines LINE with polygon POLY.
LINE is a 1-by-4 row vector containing parametric representation of the
line (in the format [x0 y0 dx dy], see the function 'createLine' for
details).
POLY is a NV-by-2 array containing coordinates of the polygon vertices
P is a K-by-2 array containing the coordinates of the K intersection
points.
P = intersectLinePolygon(LINE, POLY, TOL)
Specifies the tolerance for geometric tests. Default is 1e-14.
[P, INDS] = intersectLinePolygon(...)
Also returns the indices of edges involved in intersections. INDS is a
K-by-1 column vector, such that P(i,:) corresponds to intersection of
the line with the i-th edge of the polygon. If the intersection occurs
at a polygon vertex, the index of only one of the two neighbor edges is
returned.
Note that due to numerical approximations, the use of function
'isPointOnEdge' may give results not consistent with this function.
[P, INDS, POS] = intersectLinePolygon(...)
Also returns the relative position of each intersection point along the
line. The position can be used to sort the points.
Examples
% compute intersections between a square and an horizontal line
poly = [0 0;10 0;10 10;0 10];
line = [5 5 1 0];
intersectLinePolygon(line, poly)
ans =
10 5
0 5
% also return indices of edges
[inters, inds] = intersectLinePolygon(line, poly)
inters =
10 5
0 5
inds =
4
2
% Potentially prolematic case
% create a polygon with various configurations at y=50
poly = [10 30;30 50;45 30; 50 50; 60 70; 70 50; ...
90 30; 80 80; 50 80; 40 50; 30 80; 20 80];
figure; axis([0 100 0 100]); hold on;
drawPolygon(poly, 'b'); drawPoint(poly, 'b.');
% Computes intersection with horizontal line at y=50
line = [10 50 2 0]; drawLine(line, 'm');
points = intersectLinePolygon(line, poly);
% result is a 6-by-2 numeric array, with a double intersection
% (indices 2 and 3), resulting in five displayed intersections.
drawPoint(points, 'ko');
% sort intersection points according to x-coordinate
points2 = sortrows(points, 1);
% display pairs of successive intersection points as colored lines
for i = 1:2:size(points2, 1)
drawEdge(points2(i,:), points2(i+1,:), 'linewidth', 2);
end
References
https://web.cs.ucdavis.edu/~ma/ECS175_S00/Notes/0411_b.pdf
https://alienryderflex.com/polygon_fill/
See Also
lines2d, polygons2d, intersectLines, intersectRayPolygon, polygonEdges0001 function [points, edgeInds, linePositions] = intersectLinePolygon(line, poly, varargin) 0002 % Intersection points between a line and a polygon. 0003 % 0004 % P = intersectLinePolygon(LINE, POLY) 0005 % Returns the intersection points of the lines LINE with polygon POLY. 0006 % LINE is a 1-by-4 row vector containing parametric representation of the 0007 % line (in the format [x0 y0 dx dy], see the function 'createLine' for 0008 % details). 0009 % POLY is a NV-by-2 array containing coordinates of the polygon vertices 0010 % P is a K-by-2 array containing the coordinates of the K intersection 0011 % points. 0012 % 0013 % P = intersectLinePolygon(LINE, POLY, TOL) 0014 % Specifies the tolerance for geometric tests. Default is 1e-14. 0015 % 0016 % [P, INDS] = intersectLinePolygon(...) 0017 % Also returns the indices of edges involved in intersections. INDS is a 0018 % K-by-1 column vector, such that P(i,:) corresponds to intersection of 0019 % the line with the i-th edge of the polygon. If the intersection occurs 0020 % at a polygon vertex, the index of only one of the two neighbor edges is 0021 % returned. 0022 % Note that due to numerical approximations, the use of function 0023 % 'isPointOnEdge' may give results not consistent with this function. 0024 % 0025 % [P, INDS, POS] = intersectLinePolygon(...) 0026 % Also returns the relative position of each intersection point along the 0027 % line. The position can be used to sort the points. 0028 % 0029 % Examples 0030 % % compute intersections between a square and an horizontal line 0031 % poly = [0 0;10 0;10 10;0 10]; 0032 % line = [5 5 1 0]; 0033 % intersectLinePolygon(line, poly) 0034 % ans = 0035 % 10 5 0036 % 0 5 0037 % % also return indices of edges 0038 % [inters, inds] = intersectLinePolygon(line, poly) 0039 % inters = 0040 % 10 5 0041 % 0 5 0042 % inds = 0043 % 4 0044 % 2 0045 % 0046 % % Potentially prolematic case 0047 % % create a polygon with various configurations at y=50 0048 % poly = [10 30;30 50;45 30; 50 50; 60 70; 70 50; ... 0049 % 90 30; 80 80; 50 80; 40 50; 30 80; 20 80]; 0050 % figure; axis([0 100 0 100]); hold on; 0051 % drawPolygon(poly, 'b'); drawPoint(poly, 'b.'); 0052 % % Computes intersection with horizontal line at y=50 0053 % line = [10 50 2 0]; drawLine(line, 'm'); 0054 % points = intersectLinePolygon(line, poly); 0055 % % result is a 6-by-2 numeric array, with a double intersection 0056 % % (indices 2 and 3), resulting in five displayed intersections. 0057 % drawPoint(points, 'ko'); 0058 % % sort intersection points according to x-coordinate 0059 % points2 = sortrows(points, 1); 0060 % % display pairs of successive intersection points as colored lines 0061 % for i = 1:2:size(points2, 1) 0062 % drawEdge(points2(i,:), points2(i+1,:), 'linewidth', 2); 0063 % end 0064 % 0065 % References 0066 % https://web.cs.ucdavis.edu/~ma/ECS175_S00/Notes/0411_b.pdf 0067 % https://alienryderflex.com/polygon_fill/ 0068 % 0069 % See Also 0070 % lines2d, polygons2d, intersectLines, intersectRayPolygon, polygonEdges 0071 % 0072 0073 % ------ 0074 % Author: David Legland 0075 % e-mail: david.legland@inrae.fr 0076 % Created: 2003-10-31, using Matlab 7.9.0.529 (R2009b) 0077 % Copyright 2012 INRA - Cepia Software Platform. 0078 0079 % 2022-02-14: new algorithm, more robust to numerical issues, but that can 0080 % change behaviour compared to previous implementation 0081 0082 % line origin and angle 0083 ox = line(1); 0084 oy = line(2); 0085 dx = line(3); 0086 dy = line(4); 0087 0088 % create transform matrix that project line onto the horizontal axis 0089 % (then, computation of intersections rely only on the y-coordinate) 0090 theta = atan2(dy, dx); 0091 s = hypot(dx, dy); 0092 cot = cos(theta) / s; 0093 sit = sin(theta) / s; 0094 transfo = [cot sit 0; -sit cot 0; 0 0 1] * [1 0 -ox; 0 1 -oy; 0 0 1]; 0095 0096 % number of vertices in polygon 0097 nVertices = size(poly, 1); 0098 0099 % create arrays for storing x-coordinates and edge inds of intersections 0100 linePositions = []; 0101 edgeInds = []; 0102 0103 % retrieve previous vertex and its y-coordinate 0104 ivp = nVertices; 0105 previousVertex = transformPoint(poly(ivp,:), transfo); 0106 yvp = previousVertex(2); 0107 0108 % iterate over indices of first edge vertex 0109 for iv = 1:nVertices 0110 % current vertex and its y-coordinate 0111 currentVertex = transformPoint(poly(iv,:), transfo); 0112 yv = currentVertex(2); 0113 0114 % check conditions for intersection 0115 % either if: 0116 % 1) previous vertex is above or on, and current vertex is strictly below 0117 % 2) previous vertex is strictly below, and current vertex is above or on 0118 if yvp >= 0 && yv < 0 || yvp < 0 && yv >= 0 0119 % slope of current edge (dy cannot be zero due to above condition) 0120 edgeDx = currentVertex(1) - previousVertex(1); 0121 edgeDy = currentVertex(2) - previousVertex(2); 0122 % position of intersection on the horizontal line 0123 currentPos = currentVertex(1) - yv * edgeDx / edgeDy ; 0124 % add to list of intersections 0125 linePositions = [linePositions ; currentPos]; %#ok<AGROW> 0126 % keep list if edge indices 0127 edgeInds = [edgeInds ; ivp]; %#ok<AGROW> 0128 end 0129 0130 % switch current vertex to previous vertex 0131 previousVertex = currentVertex; 0132 ivp = iv; % keep edge index for optional output 0133 yvp = yv; 0134 end 0135 0136 % format output result into a N-by-2 array of points 0137 if ~isempty(linePositions) 0138 points = [linePositions zeros(size(linePositions))]; 0139 points = transformPoint(points, inv(transfo)); 0140 linePositions = linePositions; 0141 else 0142 points = []; 0143 end