Winding algorithm and floating point error

Question

I have written a Java implementation of the winding algorithm for solving the point in polygon problem.

This is the code:

public boolean isPointInsidePolyline(final Vec3D pt) {
    boolean isCw = this.isPolylineClockwise();
    float count = 0;
    for (Line line : this.getLines()) { // gets all the line segments of the polyline
        if (line.contains(pt)) { // checks if the point is located on the line
            return false;
        }
        Vec3D pt1 = line.getStart();
        Vec3D pt2 = line.getEnd();
        Vec3D dir1 = pt1.sub(pt).normalize();
        Vec3D dir2 = pt2.sub(pt).normalize();
        float angle = dir1.angleBetween(dir2);
        if (dir1.dot(dir2) < 1) {
            if (angle != 0) {
                if (Points.isPointSequenceClockwise(
                    Arrays.asList(pt,pt1,pt2)) == isCw
                ) {
                    count += angle;
                } else {
                    count -= angle;
                }
            }
        }
    }
    return count > eps; // eps = 0.001f;
}

The vector library used is from toxiclibs.

The algorithm overall seems to work well, but I have encountered edge cases where it fails due to floating-point errors in the calculations. Just changing the eps value doesn't seem like a robust solution as the epsilon value is not based on something.

I would like to ask if anyone else who has encountered a similar problem has managed to find a better solution to mitigate floating-point errors.

Sample data to demonstrate an example that is failing below.

Polyline (described by its points with x,y,z coordinates):

-199.21376f, -2.1003783E-5f, 0,
-204.51433f, -1.302136E-4f, 0,
-204.79602f, -6.259075E-6f, 0,
-212.74092f, 7.808674E-8f, 0,
-223.72076f, 1.8734062E-6f, 0,
-224.89467f, -4.392134E-6f, 0,
-225.87766f, 1.9165287E-7f, 0,
-327.65207f, 1.855702E-5f, 0,
-355.3976f, 3.3298825E-6f, 0,
-366.83472f, 3.8695987E-5f, 0,
-367.7069f, -2.677933E-5f, 0,
-368.92123f, 1.2206072E-4f, 0,
-370.22495f, 2.6854828E-5f, 0,
-371.65765f, 3.9196784E-5f, 0,
-399.56082f, 3.5489844E-5f, 0,
-400.69135f, 0.013866073f, 0,
-401.04166f, 0.0067931935f, 0,
-399.67267f, 48.833225f, 0,
-398.56805f, 84.82277f, 0,
-397.50436f, 92.38736f, 0,
-397.18918f, 93.921135f, 0,
-397.0501f, 94.37826f, 0,
-395.92285f, 97.33734f, 0,
-393.09155f, 104.53852f, 0,
-392.62076f, 105.47283f, 0,
-387.8274f, 113.93879f, 0,
-387.5292f, 114.346664f, 0,
-386.7448f, 115.26469f, 0,
-379.80237f, 122.86422f, 0,
-374.12714f, 127.701904f, 0,
-372.0301f, 129.37302f, 0,
-370.6804f, 130.21136f, 0,
-365.25098f, 133.43129f, 0,
-363.58215f, 134.21275f, 0,
-357.24515f, 136.94856f, 0,
-352.32858f, 138.49944f, 0,
-348.04324f, 139.74101f, 0,
-344.63644f, 140.3836f, 0,
-340.69135f, 141.17427f, 0,
-339.6991f, 142.09425f, 0,
-337.85547f, 143.91138f, 0,
-336.91998f, 150.82086f, 0,
-336.6831f, 151.84744f, 0,
-335.75406f, 151.8479f, 0,
-335.27155f, 151.84789f, 0,
-327.8963f, 151.84781f, 0,
-290.45605f, 151.84741f, 0,
-286.98816f, 138.05928f, 0,
-285.6575f, 132.89183f, 0,
-282.18118f, 120.361824f, 0,
-280.883f, 116.26616f, 0,
-281.07944f, 113.260635f, 0,
-281.1565f, 105.80201f, 0,
-281.10834f, 103.23026f, 0,
-280.75845f, 92.780136f, 0,
-237.04051f, 92.7801f, 0,
-215.375f, 92.77999f, 0,
-199.18573f, 92.78001f, 0,
-199.20132f, 65.52554f, 0,
-199.23079f, 8.813348f, 0,
-199.21376f, -2.1003798E-5f, 0

And the point I checked for inclusion has coordinates (201.28201f, 2.0f, 0)

Just by looking at the min and max x coordinates we can easily see that the point is not located inside the polyline, but the unit test I wrote still returned true.

Answer 1

The point-in-polygon problem has always been tricky, due to the floating-point representation. In particular, testing on which side of an edge a point lies is ambiguous when the point is very close. Using a tolerance is indeed not a solution as the choice of the value is arbitrary (you don't know in advance how accurate are the vertices themselves), and that just displaces the problem: points close to the tolerance limit are ambiguous as well. Testing point alignment is difficult.

On another hand, I don't like the winding number algorithm for two reasons:

• it uses a costly inverse trigonometric function,
• error analysis is difficult.

I prefer the simple half-line method:

• consider an horizontal line by the test point;
• find all polygon edges that cross this line. The test is simple and fast, based on ordinate comparisons;
• for all crossing edges, determine if the test point lies to the left or the right of the intersection with the line. This costs a 2x2 determinant computation (there is no need to explicitly compute the intersection);
• the point is inside if the number of intersections is odd.

I don't claim that this approach avoids the ambiguity problem (which is there forever). Anyway, it provides a better level of safety, is easier to interpret, and has a low cost. If you detect crossing by assigning every vertex an "above" or "below" status (but no "on" status), then it is guaranteed that the number of crossing edges is even, whatever the configuration. The algorithm will never fail (whereas the winding number cannot conclude when the test point is a vertex).

---

Update:

Sorry, I did not notice that you are addressing the problem in 3D. (In 3D, polygons are in fact not perfectly flat and this can cause additional ambiguities: a point cannot be classified as inside or outside a skew polygon). You can deal with this by projecting all points to a plane (such as a coordinate plane, or a best fit plane).