Methods for determining whether a point belongs to a polygon
Recently on Habré there was an article in which it was described how one can determine where the point is in relation to a polygon: inside or outside. A similar problem is encountered in geometric modeling and in computer graphics quite often. And since the method described in the article was not optimal, and there was a slight chaos in the comments, the idea arose to write this article. So, what algorithms exist in modern computer graphics to determine whether a given point belongs to a polygon or not.
Before starting, I want to immediately describe the problem. Although the problem itself is simple: we have a set of points and an order in which these points connect. And given the point that we test for affiliation. It is implied that we have a closed polygon, and in the general case the edges of the polygon do not intersect with each other, that is, it is spared of self-intersections. There are no restrictions on the number of vertices, that is, a polygon with a million vertices can be easily defined. We hope that the user does not ask us what is not clear, but we cannot guarantee it either. And another nuance: since we work with computer graphics, which means that we use floating-point arithmetic, which, although it allows us to operate with numbers quite accurately, is still not free from errors.
Well, sort of decided on the problem, let's now see what methods of solution exist.
I will start with the one that is considered the most popular in the world of graphics and games: ray tracing. In short, the algorithm can be described as follows:
')
Sounds easy, doesn't it? Indeed, this is the easiest method. It ultimately boils down to a certain number of intersections of a segment (a polygon face) and a ray, that is, in essence, the intersection of two lines and testing the resulting point for belonging to a ray and a segment. In the simplest case, it is necessary to cross the beam with all segments, using trees (quadratic, binary or BVH), you can reduce this number. In general, it is said that the algorithmic complexity of this algorithm is O (log n), where n is the number of edges in the polygon.
The method is simple, but, unfortunately, in general, it is better not to use it. The reason for this is the case when we intersect with a ray a vertex of a polygon or edge that partially coincides with a ray. I illustrate this with an example.
Suppose we have a polygon, and there is a point. At the very beginning we agreed that the direction would be along the x axis. Let us release a ray from a point in the positive direction of the x axis and the ray successfully crossed the polygon at the vertex. Then the question arises, how exactly do we check this situation? Do not forget that we work with floating point numbers, and small errors are possible. Let us turn into the world of analytic geometry, so that you can operate not just with geometric concepts, but with numbers.
The equation of each segment (a polygon face): a + t ( b - a ), where a is the coordinates of one end of the segment, b is the coordinates of the other end of the segment. Obviously, if the beam intersects the segment, t must be in the interval [0,1]. If we cross the vertex by a ray, then t = 0 or t = 1. This is in the ideal world of mathematics. In the real world of computer algebra, you might have something like t = 1e-10. Or t = -1e-10. Or 0.99999999999998. Or 1.0000000000001. Since the intersection of two lines involves a division procedure, this can easily happen. And what, then, is to be done? Trust the computer and assume that if we are strictly greater than or equal to zero or strictly less than or equal to one, then we consider this as an intersection, otherwise we do not consider it? We trusted and got a situation when with one edge the parameter t = -1e-20, with another edge t = 1.000000000000000000001. As a result, we do not consider this as intersections, and in our country the beam has safely passed by and the result is wrong.
Let's look in another direction. Sent the beam in a negative direction. There is also not very good - the ray intersects the vertex inside the polygon. It can also be anything. Instead of the horizontal direction, take the vertical one? No one guarantees that you will not cross the top again. In the particular example chosen by me at the top, the point is chosen in such a way that its intersection with the ray parallel to the y axis and going from top to bottom also intersects the polygon at the vertex.
And if you think that crossing with a vertex is bad, see what else can happen:
Here we intersect the ray with a segment that coincides with this ray. How to be in this case? And if not the same, but almost the same? And imagine that in a polygon there are many almost degenerate edges, how to cross with such?
The saddest thing in this whole situation is that it seems to us: “I need something very simple for my simple purposes, this situation will not affect me.” According to Murphy's law, unfortunately, this is the situation that arises whenever you do not expect it at all. And so I smoothly turn to the second method.
In general, this method has the terrible name angle weighted pseudo normals and it is connected in the concept of so-called distance fields with a sign (signed distance fields). But once again, I do not want to frighten anyone, so just let the near point and its normal (that is, the perpendicular vector) be.
The algorithm in this case is:
Moreover, the angle itself is not necessarily considered, it is enough to check the cosine sign of this angle. If positive - inside, if negative - outside. And since we are only interested in the cosine sign, then in essence we calculate the sign of the scalar product between two vectors.
Consider an example. Point A1, the nearest point for it is on the edge. If we do everything correctly, the normal to the edge is parallel to the vector from the test point to the closest one. In the case of point A1, the angle between the vectors = 0. Or almost zero, since everything is possible because of floating point operations. Less than 90 degrees, test point A1 is inside. Test point A2. Its nearest point is a vertex, the normal to which is the average normal of the edges adjacent to this vertex. We consider the scalar product of two vectors, should be negative. We are outside.
So, like with the essence of the method figured out. What about performance and floating point issues?
As in the case of point tracing, performance is O (log n), if trees are used to store information about edges. From a computational point of view, the method, although it has a similar complexity, will be somewhat slower than tracing. First of all, because the distance between the point and the edge is a slightly more expensive operation than the intersection of two lines. Floating troubles arise in this method, usually near the edges of the polygon. And the closer we are to the edge, the greater the likelihood of incorrectly determining the sign. Fortunately, the closer we are to the edge, the smaller the distance. That is, if we, for example, say that if the resulting distance is less than the predetermined minimum (it can be a constant like DBL_EPSILON or FLT_EPSILON), then the point belongs to the edge. And if it belongs to an edge, then we already decide whether a part of a polygon is an edge or not (as a rule, a part).
Describing the previous method, a lot has been said about the shortcomings. It's time to name a few flaws and this method. First of all, this method requires that all normals to the edges are directed in the correct direction. That is, before determining whether we are outside or inside, we must do some work on the calculation of these normals and their proper orientation. Very often, especially when there is a large dump of peaks and edges at the entrance, this process is not always simple. If it is necessary to determine only one point, the process of orientation of the normals may take most of the time that could be spent on something else. Also, this method does not like when the input is a polygon with self-intersections. At the beginning, I said that in our problem such a case is not considered, but if it were considered, then this method could produce completely unobvious results.
But in general, the method is not bad, especially if we have a polygon with a large number of vertices and edges at the entrance, and we need to test a lot of points for belonging. If there are few points, ray tracing is unstable, and you want something more or less reliable, then there is a third way.
This method has been known for quite some time, but mostly remains theoretical, for the most part because it is not as effective as the previous two. Here is its essence "on the fingers." Take the unit circle with the center at the test point. Then we project each vertex of the polygon onto this circle by rays that pass through the vertex and the test point. Something like this:
In the figure, points P1, P2 and so on - the vertices of the polygon, and the points P1 ', P2' and so on - their projections onto the circle. Now, when we consider the edge of a polygon, it is possible to determine from the projections whether the rotation is counterclockwise or clockwise when moving from one vertex to another. Counterclockwise rotation will be considered a positive rotation, and clockwise rotation - negative. The angle that corresponds to each edge is the angle between the segments of the circle through the projections of the vertices of this edge. Since the rotation can be positive or negative, the angle can be positive or negative.
If you then add up all these angles, the sum should be +360 or -360 degrees for the point inside and 0 for the point outside. Plus or minus is no accident here. If the polygon is counterclockwise when specifying the traversal order, it should be +360. If the edge traversal in the polygon is counterclockwise, then it turns out to be -360. To avoid numerical errors, they usually do this: divide the resulting sum by 360 and lead to the nearest integer. The resulting number, by the way, is the index of the point. The result can be one of three: -1 means that we are inside and the polygon goes around clockwise, 0 means that we are outside, +1 means that we are outside. Everything else means that we are somewhere wrong, or that the input data is incorrect.
The algorithm in this case is as follows:
Before starting, I want to immediately describe the problem. Although the problem itself is simple: we have a set of points and an order in which these points connect. And given the point that we test for affiliation. It is implied that we have a closed polygon, and in the general case the edges of the polygon do not intersect with each other, that is, it is spared of self-intersections. There are no restrictions on the number of vertices, that is, a polygon with a million vertices can be easily defined. We hope that the user does not ask us what is not clear, but we cannot guarantee it either. And another nuance: since we work with computer graphics, which means that we use floating-point arithmetic, which, although it allows us to operate with numbers quite accurately, is still not free from errors.
Well, sort of decided on the problem, let's now see what methods of solution exist.
Method 1. Ray Tracing
I will start with the one that is considered the most popular in the world of graphics and games: ray tracing. In short, the algorithm can be described as follows:
')
- From the test point we release the beam either in a predetermined or in an arbitrary direction.
- We count the number of intersections with a polygon.
- If the number of intersections is even, we are outside. If the number of intersections is odd, we are inside.
Sounds easy, doesn't it? Indeed, this is the easiest method. It ultimately boils down to a certain number of intersections of a segment (a polygon face) and a ray, that is, in essence, the intersection of two lines and testing the resulting point for belonging to a ray and a segment. In the simplest case, it is necessary to cross the beam with all segments, using trees (quadratic, binary or BVH), you can reduce this number. In general, it is said that the algorithmic complexity of this algorithm is O (log n), where n is the number of edges in the polygon.
The method is simple, but, unfortunately, in general, it is better not to use it. The reason for this is the case when we intersect with a ray a vertex of a polygon or edge that partially coincides with a ray. I illustrate this with an example.
Suppose we have a polygon, and there is a point. At the very beginning we agreed that the direction would be along the x axis. Let us release a ray from a point in the positive direction of the x axis and the ray successfully crossed the polygon at the vertex. Then the question arises, how exactly do we check this situation? Do not forget that we work with floating point numbers, and small errors are possible. Let us turn into the world of analytic geometry, so that you can operate not just with geometric concepts, but with numbers.
The equation of each segment (a polygon face): a + t ( b - a ), where a is the coordinates of one end of the segment, b is the coordinates of the other end of the segment. Obviously, if the beam intersects the segment, t must be in the interval [0,1]. If we cross the vertex by a ray, then t = 0 or t = 1. This is in the ideal world of mathematics. In the real world of computer algebra, you might have something like t = 1e-10. Or t = -1e-10. Or 0.99999999999998. Or 1.0000000000001. Since the intersection of two lines involves a division procedure, this can easily happen. And what, then, is to be done? Trust the computer and assume that if we are strictly greater than or equal to zero or strictly less than or equal to one, then we consider this as an intersection, otherwise we do not consider it? We trusted and got a situation when with one edge the parameter t = -1e-20, with another edge t = 1.000000000000000000001. As a result, we do not consider this as intersections, and in our country the beam has safely passed by and the result is wrong.
Let's look in another direction. Sent the beam in a negative direction. There is also not very good - the ray intersects the vertex inside the polygon. It can also be anything. Instead of the horizontal direction, take the vertical one? No one guarantees that you will not cross the top again. In the particular example chosen by me at the top, the point is chosen in such a way that its intersection with the ray parallel to the y axis and going from top to bottom also intersects the polygon at the vertex.
And if you think that crossing with a vertex is bad, see what else can happen:
Here we intersect the ray with a segment that coincides with this ray. How to be in this case? And if not the same, but almost the same? And imagine that in a polygon there are many almost degenerate edges, how to cross with such?
The saddest thing in this whole situation is that it seems to us: “I need something very simple for my simple purposes, this situation will not affect me.” According to Murphy's law, unfortunately, this is the situation that arises whenever you do not expect it at all. And so I smoothly turn to the second method.
Method 2. The near point and its normal
In general, this method has the terrible name angle weighted pseudo normals and it is connected in the concept of so-called distance fields with a sign (signed distance fields). But once again, I do not want to frighten anyone, so just let the near point and its normal (that is, the perpendicular vector) be.
The algorithm in this case is:
- For the test point, look for the nearest point on the polygon. At the same time, we remember that the nearest point can be not only on the segment, but also at the vertex.
- We are looking for the normal of the nearest point. If the near point lies on the edge, then the normal is a vector perpendicular to the edge and looking outward of the polygon. If the near point is one of the vertices, then the normal is the average vector of edges adjacent to the vertex.
- We calculate the angle between the normal of the nearest point and the vector from the test point to the closest one. If the angle is less than 90 degrees, then we - inside, otherwise - outside.
Moreover, the angle itself is not necessarily considered, it is enough to check the cosine sign of this angle. If positive - inside, if negative - outside. And since we are only interested in the cosine sign, then in essence we calculate the sign of the scalar product between two vectors.
Consider an example. Point A1, the nearest point for it is on the edge. If we do everything correctly, the normal to the edge is parallel to the vector from the test point to the closest one. In the case of point A1, the angle between the vectors = 0. Or almost zero, since everything is possible because of floating point operations. Less than 90 degrees, test point A1 is inside. Test point A2. Its nearest point is a vertex, the normal to which is the average normal of the edges adjacent to this vertex. We consider the scalar product of two vectors, should be negative. We are outside.
So, like with the essence of the method figured out. What about performance and floating point issues?
As in the case of point tracing, performance is O (log n), if trees are used to store information about edges. From a computational point of view, the method, although it has a similar complexity, will be somewhat slower than tracing. First of all, because the distance between the point and the edge is a slightly more expensive operation than the intersection of two lines. Floating troubles arise in this method, usually near the edges of the polygon. And the closer we are to the edge, the greater the likelihood of incorrectly determining the sign. Fortunately, the closer we are to the edge, the smaller the distance. That is, if we, for example, say that if the resulting distance is less than the predetermined minimum (it can be a constant like DBL_EPSILON or FLT_EPSILON), then the point belongs to the edge. And if it belongs to an edge, then we already decide whether a part of a polygon is an edge or not (as a rule, a part).
Describing the previous method, a lot has been said about the shortcomings. It's time to name a few flaws and this method. First of all, this method requires that all normals to the edges are directed in the correct direction. That is, before determining whether we are outside or inside, we must do some work on the calculation of these normals and their proper orientation. Very often, especially when there is a large dump of peaks and edges at the entrance, this process is not always simple. If it is necessary to determine only one point, the process of orientation of the normals may take most of the time that could be spent on something else. Also, this method does not like when the input is a polygon with self-intersections. At the beginning, I said that in our problem such a case is not considered, but if it were considered, then this method could produce completely unobvious results.
But in general, the method is not bad, especially if we have a polygon with a large number of vertices and edges at the entrance, and we need to test a lot of points for belonging. If there are few points, ray tracing is unstable, and you want something more or less reliable, then there is a third way.
Method 3. The index of the point relative to the polygon
This method has been known for quite some time, but mostly remains theoretical, for the most part because it is not as effective as the previous two. Here is its essence "on the fingers." Take the unit circle with the center at the test point. Then we project each vertex of the polygon onto this circle by rays that pass through the vertex and the test point. Something like this:
In the figure, points P1, P2 and so on - the vertices of the polygon, and the points P1 ', P2' and so on - their projections onto the circle. Now, when we consider the edge of a polygon, it is possible to determine from the projections whether the rotation is counterclockwise or clockwise when moving from one vertex to another. Counterclockwise rotation will be considered a positive rotation, and clockwise rotation - negative. The angle that corresponds to each edge is the angle between the segments of the circle through the projections of the vertices of this edge. Since the rotation can be positive or negative, the angle can be positive or negative.
If you then add up all these angles, the sum should be +360 or -360 degrees for the point inside and 0 for the point outside. Plus or minus is no accident here. If the polygon is counterclockwise when specifying the traversal order, it should be +360. If the edge traversal in the polygon is counterclockwise, then it turns out to be -360. To avoid numerical errors, they usually do this: divide the resulting sum by 360 and lead to the nearest integer. The resulting number, by the way, is the index of the point. The result can be one of three: -1 means that we are inside and the polygon goes around clockwise, 0 means that we are outside, +1 means that we are outside. Everything else means that we are somewhere wrong, or that the input data is incorrect.
The algorithm in this case is as follows:
- Set the sum of the angles to 0.
- For all edges of a polygon:
- Using the dot product, we calculate the angle formed by the vectors starting at the test point and ending at the ends of the edge.
- We calculate the vector product of these vectors. If its z-component is positive, we add an angle to the sum of the angles; otherwise, we subtract from the same sum.
We divide the sum by 2 if we count in radians or 360 if we count in degrees. The latter is unlikely, but suddenly.
Round the result to the nearest integer. +1 or -1 means we are inside. 0 we mean outside.
Consider an example. There is a polygon whose order is counterclockwise. There is a point A, which we test. For testing, we first calculate the angle between the vectors AP1 and AP2. The vector product of the same vectors looks at us, then we add to the sum. Go ahead and consider the angle between AP2 and AP3. The vector product looks at us, the resulting angle is subtracted. And so on.
For this particular picture, I counted everything, and this is what happened:
Point A.
(AP1, AP2) = 74.13, (AP2, AP3) = 51.58, (AP3, AP4) = 89.99, (AP4, AP5) = 126.47, (AP5, AP1) = 120.99.
sum = 74.13-51.58 + 89.99 + 126.47 + 120.99 = 360. 360/360 = 1 Point - inside.
Point B.
(BP1, BP2) = 44.78, (BP2, BP3) = 89.11, (BP3, BP4) = 130.93, (BP4, BP5) = 52.97, (BP5, BP1) = 33.63.
sum = -44.78 + 89.11-130.93 + 52.97 + 33.63 = 0. The point is outside.
And traditionally describe the pros and cons of this approach. Let's start with the cons. The method is simple mathematically, but not so efficient in terms of performance. First, its algorithmic complexity is O (n) and, whatever one may say, and all edges of the polygon will have to be searched. Secondly, to calculate the angle, one will have to use the arc cosine and two root operations (the scalar product formula and its connection with the angle to help those who do not understand why). These operations are very expensive in terms of speed, and, moreover, the errors associated with them can be significant. And third, the algorithm does not directly determine the point lying on the edge. Either - outside, or - inside. There is no third. However, the last drawback is easily determined: if at least one of the angles is equal to (or almost equal to) 180 degrees, this automatically means an edge.
The disadvantages of the method are somewhat compensated by its merits. First, it is the most stable method. If the input polygon is valid, then the result is correct for all points, except perhaps the points on the edges, but see above for details. Moreover, the method allows to partially deal with incorrect input data. Is the polygon self-intersecting? It does not matter, the method is likely to determine most points correctly. Is the polygon not closed or not a polygon at all, but a little meaningful set of segments? The method will determine the points correctly in a large number of cases. In general, all the method is good, but slow and requires the calculation of arc cosines.
What would you like to finish this review? And the fact that there are more than one and even not two methods for solving the problem of determining whether a point belongs to a polygon. They serve different purposes and some are more suitable in cases where speed is important, others when quality is important. Well, do not forget that we have unpredictable input data and we work with a computer whose floating point arithmetic is prone to errors. If you need speed and quality is completely unimportant - ray tracing to help. In most real-world applications, the near-point and normal method will most likely help. If in the first place - the accuracy of determination with incomprehensible input data, the point index method should help.
If you have any questions, ask. As a person engaged in geometry and similar problems related to graphics, I will be happy to help as much as I can.
Source: https://habr.com/ru/post/301102/
All Articles