INTERSECTLINEPOLYGON 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 %INTERSECTLINEPOLYGON 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 2003-2024 INRA - Cepia Software Platform 0078 0079 % line origin and angle 0080 ox = line(1); 0081 oy = line(2); 0082 dx = line(3); 0083 dy = line(4); 0084 0085 % create transform matrix that project line onto the horizontal axis 0086 % (then, computation of intersections rely only on the y-coordinate) 0087 theta = atan2(dy, dx); 0088 s = hypot(dx, dy); 0089 cot = cos(theta) / s; 0090 sit = sin(theta) / s; 0091 transfo = [cot sit 0; -sit cot 0; 0 0 1] * [1 0 -ox; 0 1 -oy; 0 0 1]; 0092 0093 % number of vertices in polygon 0094 nVertices = size(poly, 1); 0095 0096 % create arrays for storing x-coordinates and edge inds of intersections 0097 linePositions = []; 0098 edgeInds = []; 0099 0100 % retrieve previous vertex and its y-coordinate 0101 ivp = nVertices; 0102 previousVertex = transformPoint(poly(ivp,:), transfo); 0103 yvp = previousVertex(2); 0104 0105 % iterate over indices of first edge vertex 0106 for iv = 1:nVertices 0107 % current vertex and its y-coordinate 0108 currentVertex = transformPoint(poly(iv,:), transfo); 0109 yv = currentVertex(2); 0110 0111 % check conditions for intersection 0112 % either if: 0113 % 1) previous vertex is above or on, and current vertex is strictly below 0114 % 2) previous vertex is strictly below, and current vertex is above or on 0115 if yvp >= 0 && yv < 0 || yvp < 0 && yv >= 0 0116 % slope of current edge (dy cannot be zero due to above condition) 0117 edgeDx = currentVertex(1) - previousVertex(1); 0118 edgeDy = currentVertex(2) - previousVertex(2); 0119 % position of intersection on the horizontal line 0120 currentPos = currentVertex(1) - yv * edgeDx / edgeDy ; 0121 % add to list of intersections 0122 linePositions = [linePositions ; currentPos]; %#ok<AGROW> 0123 % keep list if edge indices 0124 edgeInds = [edgeInds ; ivp]; %#ok<AGROW> 0125 end 0126 0127 % switch current vertex to previous vertex 0128 previousVertex = currentVertex; 0129 ivp = iv; % keep edge index for optional output 0130 yvp = yv; 0131 end 0132 0133 % format output result into a N-by-2 array of points 0134 if ~isempty(linePositions) 0135 points = [linePositions zeros(size(linePositions))]; 0136 points = transformPoint(points, inv(transfo)); 0137 else 0138 points = []; 0139 end