Verifying that a point belongs to a nonconvex polygon

Verifying that a point belongs to a nonconvex polygon in linear time is not at all difficult. One of the most common methods is to fire a beam and count the number of intersection points. However, it is necessary to carefully consider the cases when the points of the polygon fall on the ray. This naturally raises the question of how to deal with these cases the easiest way?

Task statement

We will solve a variant of the problem in which the vertices of the polygon and the point being checked are given by integer coordinates up to 10 ⁸ . For the lack of simplicity, we will not impose any restrictions on the vertices of the polygon. Solve the coinciding vertices, vertices on the sides, self-tilting and zero area. We prohibit self-intersections, otherwise it is not possible to determine the inside of the resulting broken line. The program should return 0 if the point lies on the border of the polygon, otherwise it is necessary to output 1 if the point lies outside the polygon and -1 if the point lies inside the polygon.

General solution scheme

Let x ₁ , x ₂ , ..., x _n be the vertices of a polygon (they can be given either in a counterclockwise direction), y is the point that we need to check. If n = 1, then it suffices to check the equality of y and x ₁ . Otherwise, we calculate the expression F ( y ) = f ( x ₁ , x ₂ , y ) · f ( x ₂ , x ₃ , y ) ··· f ( x _{n -1} , x _n , y ) · f ( x _n , x ₁ , y ). The function f ( a , b , c ) takes the values ​​1, 0, -1. It determines whether the ray intercepted from point c along the x axis in the direction of increasing x intersects the segment ( a , b ) taking into account all extreme cases. If it crosses, it will return -1. If point c lies on the segment ( a , b ) it returns 0, otherwise it will return 1.
')

Optimize the function f

The entire complexity of checking whether a point belongs to a polygon is concentrated in the function f . Therefore, I had the idea to sort through its implementation and find the shortest among them. It is clear that the first step in the function f is the calculation of the coordinates a and b relative to the point c . Put these coordinates in the local variables ax , ay , bx and by . To check the implementations, I built all possible segments ( a , b ) with coordinates no more than 5 and found the right answers for them. In addition, I restricted the expressions that can be used in f to signum(ax * by - ay * bx) , signum(ax) , signum(bx) , signum(ay) , signum(by) . Thus, the function f should return a response based on the value of these 5 values, which can be -1, 0 or 1.
In the beginning, I tried to express the answer as a sign of a linear combination of these quantities. Unfortunately, nothing came of it. The next idea is to use a chain of nested ifs. As the complexity of the program, select the number of ifs + number of returns. As if conditions, we solve expressions of the form v _i in S , where v _i is one of the 5 signum expressions above, and S is the subset {-1, 0, 1}.
Now you can run brute force. As an input, the enumeration function takes a set of allowed values ​​for each of 5 variables and a set of test numbers that satisfy these conditions. As a result, the brute force function produces a tree of functions with minimal complexity, which passes all tests. To speed up the search, all results are memorized and not re-calculated.
The number of options turned out to be small enough to complete a brute force program completion in a few seconds. As a result, the program with complexity 27 turned out to be the optimal one. After I allowed to return expressions of the form v _i and - v _i besides constants, the complexity dropped to 18. This is still too much. Careful study of the program allowed us to detect the xor-like structure when checking the ax and bx signs. Therefore, I added the signum(ax * bx) variable signum(ax * bx) and the complexity dropped to 13. The experiment with the addition of the variable signum(ay * by) reduced the complexity to 11. The result was the following function:

 public static int check2(Point a, Point b, Point middle) { long ax = ax - middle.x; long ay = ay - middle.y; long bx = bx - middle.x; long by = by - middle.y; if (ay * by > 0) return 1; int s = Long.signum(ax * by - ay * bx); if (s == 0) { if (ax * bx <= 0) return 0; return 1; } if (ay < 0) return -s; if (by < 0) return s; return 1; } 

It turned out to be a little shorter in the number of lines and noticeably shorter in the number of characters of the initial version that was used to generate tests. The short length of this program makes it easy to remember:

1. At the beginning, we check whether the segment ( a , b ) does not lie strictly on one side of the ray. In this case, the intersection is not exactly.
2. Check the degenerate case where point c lies on the line ( a , b ). The answer zero can be avoided only if the points lie on one side of the y- axis. If both of them lie on the y axis, then the case when c does not lie on the segment ( a , b ) has already been considered in the previous paragraph.
3. Xor-like code to check the intersection of the beam with the segment. The sign of s depends on which point is below the ray. Here you can put a minus of s in either of two cases; sign less, replace with more or less or equal.
4. We return 1 to the case when one point lies on the ray, and the other is on some fixed side of it. Such intersections are not counted, so as not to count one intersection two times if the vertex of the polygon hit the ray.

But this is not the limit! In this implementation, ay is multiplied by by . Why not get rid of multiplication by replacing it with something like ay < 0 ^ by < 0 ? As a result, the result is not exactly the same, but you can add this variable to the search and see what happens. The next start of the algorithm ... and it gives a solution with complexity 9! The new version looks like this:

 public static int check3(Point a, Point b, Point middle) { long ax = ax - middle.x; long ay = ay - middle.y; long bx = bx - middle.x; long by = by - middle.y; int s = Long.signum(ax * by - ay * bx); if (ay < 0 ^ by < 0) { if (by < 0) return s; return -s; } if (s == 0 && ax * bx <= 0) return Long.signum(ay * by); return 1; } 

It can be improved a little more. After checking for ay < 0 ^ by < 0 Long.signum(ay * by) value Long.signum(ay * by) cannot be -1. Therefore, the last lines can be rewritten as:

 public static int check3(Point a, Point b, Point middle) { long ax = ax - middle.x; long ay = ay - middle.y; long bx = bx - middle.x; long by = by - middle.y; int s = Long.signum(ax * by - ay * bx); if (ay < 0 ^ by < 0) { if (by < 0) return s; return -s; } if (s == 0 && (ay == 0 || by == 0) && ax * bx <= 0) return 0; return 1; } 

In this version, the two top level if'a can be rearranged. Measurements show that it gives a small gain in operating time on all tests:

 public static int check3(Point a, Point b, Point middle) { long ax = ax - middle.x; long ay = ay - middle.y; long bx = bx - middle.x; long by = by - middle.y; int s = Long.signum(ax * by - ay * bx); if (s == 0 && (ay == 0 || by == 0) && ax * bx <= 0) return 0; if (ay < 0 ^ by < 0) { if (by < 0) return s; return -s; } return 1; } 

In the new algorithm, again, there is a large number of degrees of freedom. In the expression ay < 0 ^ by < 0 , instead of the less sign, you can use the sign more or lax sign. The minus sign of s again be put in either of two cases. The by variable can be changed in the nested if'e to ay .
The hardest thing in the new algorithm is to understand the case of zero. It simultaneously checks two independent cases.

In the first case, one of the points of the segment ( a , b ) coincides with the point c , in the second case the segment is parallel to the x axis and passes through the point c . These checks can be performed in other ways, but the method above is the shortest.

Testing for performance at worst shows that on a 32-bit JVM, check2 is faster than check3. On a 64-bit JVM, the result is not clear. However, the difference in both cases is quite small - 5-10%. Attempts to replace ax * bx <= 0 inequalities for some reason do not give noticeable improvements in either the 32-bit or 64-bit versions of the java machine.

Sources for those who want to experiment: PolyCheck.zip