Community hub
Recent from talks
Contribute something
Nothing was collected or created yet.
Point in polygon
View on Wikipedia
In computational geometry, the point-in-polygon (PIP) problem asks whether a given point in the plane lies inside, outside, or on the boundary of a polygon. It is a special case of point location problems and finds applications in areas that deal with processing geometrical data, such as computer graphics, computer vision, geographic information systems (GIS), motion planning, and computer-aided design (CAD).
An early description of the problem in computer graphics shows two common approaches (ray casting and angle summation) in use as early as 1974.[1]
An attempt of computer graphics veterans to trace the history of the problem and some tricks for its solution can be found in an issue of the Ray Tracing News.[2]
Ray casting algorithm
[edit]
One simple way of finding whether the point is inside or outside a simple polygon is to test how many times a ray, starting from the point and going in any fixed direction, intersects the edges of the polygon. If the point is on the outside of the polygon the ray will intersect its edge an even number of times. If the point is on the inside of the polygon then it will intersect the edge an odd number of times. The status of a point on the edge of the polygon depends on the details of the ray intersection algorithm.
This algorithm is sometimes also known as the crossing number algorithm or the even–odd rule algorithm, and was known as early as 1962.[3] The algorithm is based on a simple observation that if a point moves along a ray from infinity to the probe point and if it crosses the boundary of a polygon, possibly several times, then it alternately goes from the outside to inside, then from the inside to the outside, etc. As a result, after every two "border crossings" the moving point goes outside. This observation may be mathematically proved using the Jordan curve theorem.
Limited precision
[edit]If implemented on a computer with finite precision arithmetics, the results may be incorrect if the point lies very close to that boundary, because of rounding errors. For some applications, like video games or other entertainment products, this is not a large concern since they often favor speed over precision. However, for a formally correct computer program, one would have to introduce a numerical tolerance ε and test in line whether P (the point) lies within ε of L (the Line), in which case the algorithm should stop and report "P lies very close to the boundary."
Most implementations of the ray casting algorithm consecutively check intersections of a ray with all sides of the polygon in turn. In this case the following problem must be addressed. If the ray passes exactly through a vertex of a polygon, then it will intersect 2 segments at their endpoints. While it is OK for the case of the topmost vertex in the example or the vertex between crossing 4 and 5, the case of the rightmost vertex (in the example) requires that we count one intersection for the algorithm to work correctly. A similar problem arises with horizontal segments that happen to fall on the ray. The issue is solved as follows: If the intersection point is a vertex of a tested polygon side, then the intersection counts only if the other vertex of the side lies below the ray. This is effectively equivalent to considering vertices on the ray as lying slightly above the ray.
Once again, the case of the ray passing through a vertex may pose numerical problems in finite precision arithmetics: for two sides adjacent to the same vertex the straightforward computation of the intersection with a ray may not give the vertex in both cases. If the polygon is specified by its vertices, then this problem is eliminated by checking the y-coordinates of the ray and the ends of the tested polygon side before actual computation of the intersection. In other cases, when polygon sides are computed from other types of data, other tricks must be applied for the numerical robustness of the algorithm.
Winding number algorithm
[edit]Another technique used to check if a point is inside a polygon is to compute the given point's winding number with respect to the polygon. If the winding number is non-zero, the point lies inside the polygon. This algorithm is sometimes also known as the nonzero-rule algorithm.
One way to compute the winding number is to sum up the angles subtended by each side of the polygon.[4] However, this involves costly inverse trigonometric functions, which generally makes this algorithm performance-inefficient (slower) compared to the ray casting algorithm. Luckily, these inverse trigonometric functions do not need to be computed. Since the result, the sum of all angles, can add up to 0 or (or multiples of ) only, it is sufficient to track through which quadrants the polygon winds,[5] as it turns around the test point, which makes the winding number algorithm comparable in speed to counting the boundary crossings.
An improved algorithm to calculate the winding number was developed by Dan Sunday in 2001.[6] It does not use angles in calculations, nor any trigonometry, and functions exactly the same as the ray casting algorithms described above. Sunday's algorithm works by considering an infinite horizontal ray cast from the point being checked. Whenever that ray crosses an edge of the polygon, Juan Pineda's edge crossing algorithm (1988)[7] is used to determine how the crossing will affect the winding number. As Sunday describes it, if the edge crosses the ray going "upwards", the winding number is incremented; if it crosses the ray "downwards", the number is decremented. Sunday's algorithm gives the correct answer for nonsimple polygons, whereas the boundary crossing algorithm fails in this case.[6]
Implementations
[edit]SVG
[edit]Similar methods are used in SVG for defining a way of filling with color various shapes (such as path, polyline, polygon, text etc.).[8] The algorithm of filling is influenced by 'fill-rule' attribute. The value may be either nonzero or evenodd. For example, in a pentagram, there is a central "hole" (visible background) with evenodd, and none with nonzero attribute.[9]
For simple polygons, the algorithms will give the same result. However, for complex polygons, the algorithms may give different results for points in the regions where the polygon intersects itself, where the polygon does not have a clearly defined inside and outside. One solution using the even-odd rule is to transform (complex) polygons into simpler ones that are even-odd-equivalent before the intersection check.[10] This, however, is computationally expensive. It is less expensive to use the fast non-zero winding number algorithm, which gives the correct result even when the polygon overlaps itself.
Point in polygon queries
[edit]The point in polygon problem may be considered in the general repeated geometric query setting: given a single polygon and a sequence of query points, quickly find the answers for each query point. Clearly, any of the general approaches for planar point location may be used. Simpler solutions are available for some special polygons.
Special cases
[edit]![[icon]](https://en.wikipedia.org/wiki/File:Wiki_letter_w_cropped.svg){kind=small} | This section needs expansion. You can help by adding to it. (August 2013) |
Simpler algorithms are possible for monotone polygons, star-shaped polygons, convex polygons and triangles.
The triangle case can be solved easily by use of a barycentric coordinate system, parametric equation or dot product.[11] The dot product method extends naturally to any convex polygon.
References
[edit]- ^ Ivan Sutherland et al.,"A Characterization of Ten Hidden-Surface Algorithms" 1974, ACM Computing Surveys vol. 6 no. 1.
- ^ Mark Vandewettering; Eric Haines; Edward John Kalenda; et al. (October 1, 1990), "Point in Polygon, One More Time...", Ray Tracing News, 3 (4)
- ^ Shimrat, M., "Algorithm 112: Position of point relative to polygon" 1962, Communications of the ACM Volume 5 Issue 8, Aug. 1962. https://dl.acm.org/doi/10.1145/368637.368653
- ^ Hormann, K.; Agathos, A. (2001). "The point in polygon problem for arbitrary polygons". Computational Geometry. 20 (3): 131. doi:10.1016/S0925-7721(01)00012-8.
- ^ Weiler, Kevin (1994), "An Incremental Angle Point in Polygon Test", in Heckbert, Paul S. (ed.), Graphics Gems IV, San Diego, CA, USA: Academic Press Professional, Inc., pp. 16–23, ISBN 0-12-336155-9
- ^ a b Sunday, Dan (2001). "Inclusion of a Point in a Polygon". Archived from the original on 26 January 2013.
- ^ Pineda, Juan (August 1988). A Parallel Algorithm for Polygon Rasterization (PDF). SIGGRAPH'88. Computer Graphics. Vol. 22, no. 4. Atlanta. Retrieved 8 August 2021.
- ^ "Painting: Filling, Stroking, Colors and Paint Servers – SVG Tiny 1.2". www.w3.org. Retrieved 2021-07-24.
- ^ "Painting: Filling, Stroking, Colors and Paint Servers – SVG Tiny 1.2". www.w3.org. Retrieved 2021-07-24.
- ^ Michael Galetzka, Patrick Glauner (2017). A Simple and Correct Even-Odd Algorithm for the Point-in-Polygon Problem for Complex Polygons. Proceedings of the 12th International Joint Conference on Computer Vision, Imaging and Computer Graphics Theory and Applications (VISIGRAPP 2017), Volume 1: GRAPP.
- ^ Accurate point in triangle test "...the most famous methods to solve it"
See also
[edit]- Java Topology Suite (JTS)
- Discussion: http://www.ics.uci.edu/~eppstein/161/960307.html
- Winding number versus crossing number methods: https://web.archive.org/web/20131210180851/http://geomalgorithms.com/a03-_inclusion.html, also available at Scribd https://www.scribd.com/document/735206906/Inclusion-of-a-Point-in-a-Polygon (subscription required)
Point in polygon
View on GrokipediaFundamentals
Problem Definition
The point-in-polygon (PIP) problem is a fundamental decision problem in computational geometry that involves determining whether a given point in the plane lies in the interior, exterior, or on the boundary of a polygon.[6] This problem arises frequently in applications requiring spatial queries, such as identifying regions enclosed by boundaries defined by vertex coordinates. The polygon is typically represented as a closed chain of line segments connecting an ordered sequence of vertices, and the task is to classify the point's location relative to this boundary.[7] Polygons in the PIP problem are often assumed to be simple, meaning their edges do not intersect except at vertices, forming a single closed loop without self-intersections.[8] The vertices are commonly listed in counterclockwise order for the outer boundary to ensure consistent orientation, though some algorithms accommodate clockwise ordering or require preprocessing to standardize it.[9] The problem can extend to polygons with holes, where interior holes are represented as additional simple polygons with opposite orientation (e.g., clockwise), and a point is considered inside the overall region if it lies within the outer polygon but outside all holes.[10] Self-intersecting polygons, while more complex, can also be handled by certain methods that account for multiple boundary crossings.[6] The classification of the point's location includes three cases: interior (strictly inside the polygon, not on any edge), exterior (strictly outside), and boundary (lying exactly on an edge or at a vertex).[11] Points on the boundary require special handling in algorithms to avoid ambiguity, often using criteria like whether the point coincides with a vertex or lies midway on an edge segment. For polygons with holes, boundary classification applies similarly to both outer and inner boundaries. To illustrate, consider a simple triangular polygon with vertices at (0,0), (4,0), and (2,3). A test point at (2,1) lies in the interior, as it is enclosed by the three edges, while a point at (5,1) is exterior. Algorithms like ray casting address this classification by analyzing intersections between a ray from the test point and the polygon edges.[7]Mathematical Prerequisites
The Jordan curve theorem provides the foundational topological basis for point-in-polygon problems, asserting that every simple closed curve in the plane divides the Euclidean plane into exactly two regions: an interior (bounded) region and an exterior (unbounded) region, with the curve itself forming the boundary between them.[12] This theorem guarantees that for a simple closed polygonal chain, points can be unambiguously classified as interior or exterior, excluding the boundary, and underpins the binary nature of containment tests in computational geometry.[13] A simple polygon is defined as a closed polygonal chain in the plane that does not intersect itself, forming a Jordan curve whose vertices connect via straight line segments.[13] Convex polygons represent a special case where the line segment joining any two points in the polygon lies entirely within it, equivalent to all interior angles being less than or equal to 180 degrees and the polygon containing its convex hull.[14] Polygon orientation refers to the direction of traversal along the boundary, typically counterclockwise for the exterior boundary (positive orientation) or clockwise for interior boundaries, which can be determined by the sign of the signed area computed via the shoelace formula.[15] Complex polygons, such as those with holes, consist of an outer simple boundary enclosing one or more disjoint inner simple boundaries (holes), where the region between the outer and inner boundaries defines the polygon's interior, requiring consistent orientation conventions (counterclockwise for outer, clockwise for holes) to maintain topological integrity.[13] In coordinate geometry, determining point-line segment relations is essential, particularly tests for whether a point lies on a line segment. To check collinearity of a point with segment endpoints and , compute the cross product of vectors and ; if zero, the points are collinear, and further verify that lies within the bounding box of and (i.e., and similarly for -coordinates).[16] Vector concepts, specifically the two-dimensional cross product, enable orientation tests for triples of points. For points , , and , the orientation is given by the sign of the scalar where a positive value indicates a counterclockwise (left) turn from to , negative indicates clockwise (right), and zero indicates collinearity.[16] This cross product magnitude also equals twice the signed area of the triangle formed by the points, reinforcing its role in geometric predicates.[15]Core Algorithms
Ray Casting Algorithm
The ray casting algorithm, also known as the crossing number method, determines whether a given point lies inside a polygon by casting a ray from the point and counting the number of intersections with the polygon's edges. A horizontal ray is typically extended from the test point to positive infinity along the x-axis, and the parity of the intersection count decides the location: an odd count indicates the point is inside the polygon, while an even count indicates it is outside. This approach relies on the even-odd rule, which assumes the polygon is simple and closed, ensuring that each boundary crossing toggles the interior/exterior status.[17][18] The step-by-step process begins by representing the polygon as a sequence of vertices for to , with the edge from vertex connecting back to vertex 0 to close the shape. For a test point , initialize a counter to zero. Iterate over each edge defined by endpoints and . An intersection occurs if the ray crosses the edge, which requires the y-coordinate of the point to lie between the edge's y-coordinates, including the lower endpoint but excluding the upper endpoint to handle vertices consistently, and the computed x-intersection to be at or beyond the test point's x-coordinate. This ensures vertices are not double-counted, as the ray only "crosses" when entering or exiting the polygon interior.[17][10] The key condition for intersection with a non-horizontal edge is that , allowing the ray to count crossings that include the lower vertex but not the upper to avoid double-counting at vertices. The x-coordinate of the intersection point is then calculated via linear interpolation along the edge: If , increment the counter, as the crossing lies to the right of the test point (note the strict inequality here aligns with common implementations treating boundary points separately). This formula assumes ; horizontal edges (where ) are skipped, as they do not cross the horizontal ray.[17][18] Vertical edges, where but , are handled by the same y-straddle condition; the intersection x is simply , and the counter increments if , treating the edge as a crossing if the ray reaches it. The choice of a horizontal ray to the right simplifies computations by aligning with the x-axis, avoiding trigonometric functions, though other directions can be used with adjusted intersection tests. Care is taken with the inequality directions to consistently classify boundary points, often considering points on edges as inside or outside based on application needs.[17][10] For clarity, the following pseudocode outlines the algorithm for a polygon with vertices stored in arrays and :function isInside(x, y, n, px[], py[]):
count = 0
for i = 0 to n-1:
j = (i + 1) % n
if ((py[i] > y) != (py[j] > y)) and (x < px[i] + (px[j] - px[i]) * (y - py[i]) / (py[j] - py[i])):
count = 1 - count // Toggle parity
return count == 1 // Odd parity means inside
function isInside(x, y, n, px[], py[]):
count = 0
for i = 0 to n-1:
j = (i + 1) % n
if ((py[i] > y) != (py[j] > y)) and (x < px[i] + (px[j] - px[i]) * (y - py[i]) / (py[j] - py[i])):
count = 1 - count // Toggle parity
return count == 1 // Odd parity means inside
Winding Number Algorithm
The winding number algorithm assesses whether a test point lies inside a polygon by calculating the winding number, which quantifies the net number of complete revolutions the polygon's boundary makes around the point when traversed in a consistent direction, typically counterclockwise. This integer value represents the topological winding of the curve around the point; a non-zero winding number indicates that the point is enclosed by the polygon, while zero signifies it is outside. The method originates from complex analysis and topology but is adapted for computational geometry to handle planar polygons robustly.[17] To compute the winding number, the algorithm sums the signed angular differences subtended by each polygon edge at the test point. For each edge connecting vertices and , the signed angles from the test point to these vertices are determined using the two-argument arctangent function, atan2, which accounts for the quadrant and provides angles in the range . The difference between consecutive angles is normalized to lie within to avoid accumulation errors from wrapping around , and these differences are accumulated over all edges. The total sum is then divided by and rounded to the nearest integer to yield the winding number . The core formula is where is the number of edges, and is the normalized signed angle difference for the -th edge, with as the test point coordinates. The point is inside if , often using a small epsilon threshold in floating-point implementations for numerical stability.[19] This non-zero winding rule contrasts with the even-odd rule, which classifies a point as inside based on the parity of ray-boundary intersections rather than directional enclosure. The winding number excels in handling self-intersecting polygons, such as star shapes, where it distinguishes true interiors from intersection artifacts by tracking net enclosure, and in multiply-connected regions like polygons with holes, where separate boundary components can contribute positively or negatively to the total winding. These properties make it suitable for complex topological scenarios without requiring decomposition.[17] A standard pseudocode implementation for the atan2-based winding number computation is as follows:function winding_number(point p, polygon vertices):
sum = 0.0
n = length(vertices)
for i from 0 to n-1:
v1 = vertices[i]
v2 = vertices[(i + 1) mod n]
theta1 = atan2(v1.y - p.y, v1.x - p.x)
theta2 = atan2(v2.y - p.y, v2.x - p.x)
dtheta = theta2 - theta1
// Normalize delta to [-pi, pi]
if dtheta > pi:
dtheta -= 2 * pi
else if dtheta < -pi:
dtheta += 2 * pi
sum += dtheta
w = round(sum / (2 * pi))
return w // Non-zero indicates inside
function winding_number(point p, polygon vertices):
sum = 0.0
n = length(vertices)
for i from 0 to n-1:
v1 = vertices[i]
v2 = vertices[(i + 1) mod n]
theta1 = atan2(v1.y - p.y, v1.x - p.x)
theta2 = atan2(v2.y - p.y, v2.x - p.x)
dtheta = theta2 - theta1
// Normalize delta to [-pi, pi]
if dtheta > pi:
dtheta -= 2 * pi
else if dtheta < -pi:
dtheta += 2 * pi
sum += dtheta
w = round(sum / (2 * pi))
return w // Non-zero indicates inside
Implementation Challenges
Precision and Numerical Issues
Floating-point arithmetic introduces significant challenges in point-in-polygon (PIP) implementations, primarily due to rounding errors that accumulate during computations such as intersection tests in ray casting and angle summations in winding number algorithms. These errors often manifest as misclassifications of points near polygon boundaries, where small perturbations in coordinate representations can alter the outcome of critical comparisons or calculations. For instance, in intersection tests, the evaluation of determinants for orientation predicates may yield incorrect signs due to roundoff, leading to erroneous counts of edge crossings.[20] Similarly, angle calculations in winding number methods suffer from imprecise summation of floating-point values, potentially resulting in totals that deviate from the expected integer multiples of .[21] In the ray casting algorithm, specific vulnerabilities arise from exact equality checks, such as comparing the query point's y-coordinate to an edge's endpoints to determine crossings, where floating-point inexactness can cause failures to detect or properly handle horizontal edges. This may lead to under- or over-counting intersections, particularly when the ray aligns closely with edge orientations. Historical developments in computational geometry during the 1970s and 1980s exacerbated these issues, as varying floating-point formats across hardware led to inconsistent behaviors in PIP tests; the standardization of IEEE 754 in 1985 provided a uniform framework but did not eliminate the need for careful error management in geometric predicates.[20][22] Mitigation strategies focus on enhancing numerical robustness without sacrificing performance. A common approach involves introducing epsilon tolerances in comparisons, such as treating as a near-equality to avoid strict floating-point matches, where is typically on the order of machine epsilon (e.g., for double precision). More advanced techniques employ robust predicates via adaptive precision arithmetic, which refines computations iteratively until the result's sign is reliably determined, or exact arithmetic libraries that simulate arbitrary-precision operations for critical tests. For example, in degenerate cases like a point exactly on an edge, unmitigated implementations may trigger division by zero in parametric intersection formulas or produce infinite angles in winding computations, but robust methods bound errors to ensure consistent classification. These approaches, rooted in seminal work on geometric predicates, have been formally verified in modern PIP implementations to guarantee correctness under IEEE 754 semantics.[20][21]Edge and Vertex Cases
In point-in-polygon tests, edge and vertex cases arise when the query point lies exactly on a polygon edge or at a vertex, requiring careful conventions to resolve ambiguities in classification as inside, outside, or on the boundary. These cases are critical because they can lead to inconsistencies if not handled uniformly, depending on the application—such as treating boundaries as inside for closed regions in geographic systems or as outside for strict containment in graphics. Common conventions include classifying boundary points as inside to ensure continuity in filled areas, though some implementations opt for a distinct "on boundary" category to avoid misclassification.[19] For the ray casting algorithm, vertex handling prevents double-counting intersections that could erroneously toggle the inside/outside parity. A standard rule counts a crossing only if the ray intersects the upper vertex of an edge (where the y-coordinate matches the vertex but the x-coordinate is to the right), treating the ray as passing infinitesimally above horizontal edges and ignoring those below to maintain odd/even consistency. This approach, often called "simulation of simplicity," avoids ambiguity by conceptually perturbing the ray slightly without actual computation. For tangent rays grazing a vertex, the rule ensures the intersection is counted once, as in cases where the ray aligns collinearly with an edge endpoint; diagramatically, imagine a horizontal ray from the point hitting a polygon's vertex—only the leftward edge's upper endpoint triggers a count if the vertex is the higher y-extremum.[23][17] In the winding number algorithm, vertex cases are managed through continuous angle accumulation to prevent discontinuous jumps at cusps or reflex angles. The winding number sums signed angular contributions from each edge, classifying the point as inside if nonzero; at a vertex, the angle is computed incrementally using atan2 differences, ensuring smooth wrapping by checking quadrant transitions and determinants to resolve collinear alignments (e.g., a 180-degree turn at a concave vertex contributes ±1 to the total revolutions without reset). This handles special polygons like concave shapes, where reflex vertices (interior angles >180°) alter the local winding but preserve global enclosure, and cusps, where edges meet at 0° or 360°, by treating them as degenerate edges with zero angular contribution. For polygons with holes, consistent counterclockwise orientation for outer boundaries and clockwise for holes ensures the winding number subtracts inner contributions, maintaining nonzero for enclosed regions. Diagramatically, for a collinear point on an edge, the winding might accumulate to ±0.5 per side, but boundary detection via zero determinant flags it separately.[17][19] Standard libraries adopt specific inclusion rules to resolve these cases, but implementations vary, leading to inconsistent behavior. The CGAL library'sbounded_side_2 function uses the even-odd rule and returns a distinct ON_BOUNDARY for points on edges or vertices, handling degeneracies like ray-vertex intersections robustly without misclassification. In contrast, other codes like PNPOLY classify boundary points consistently as inside or outside but may fail on vertices due to roundoff, while CSG and barycentric methods show higher error rates (up to 53% false negatives on vertices in benchmarks). Studies on real-world polygons (e.g., 304 vertices) reveal false positives/negatives near edges ranging from 0.4%/0% in CGAL to 23.9%/23.7% in CSG, highlighting the need for degeneracy-aware designs.[24][25][23]