绕数算法（Winding Algorithm）浮点误差问题及解决方案咨询

Great question—floating-point inaccuracies are a common pain point in computational geometry, and tweaking a global epsilon is rarely the most robust fix. Let’s break down why your test case is failing and explore several better solutions.

First, let’s recap your setup: you’ve implemented a winding algorithm using toxiclibs, and a point clearly outside your polygon is incorrectly returning true due to accumulated floating-point errors in angle calculations.

Your Current Implementation #

public boolean isPointInsidePolyline(final Vec3D pt) {
    boolean isCw = this.isPolylineClockwise();
    float count = 0;
    for (Line line : this.getLines()) { // 获取折线的所有线段
        if (line.contains(pt)) { // 检查点是否在线段上
            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;
}

Failing Test Case #

Polygon coordinates (x,y,z per vertex):

-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

Test point: (201.28201f, 2.0f, 0)

Solutions to Fix the Issue #

1. Switch to Cross-Product-Based Winding Count (Most Robust)

Angle calculations are inherently prone to floating-point drift. A better approach uses cross products to track the winding number directly, avoiding angle math entirely. Here’s a revised implementation:

public boolean isPointInsidePolyline(final Vec3D pt) {
    int windingNumber = 0;
    List<Line> lines = this.getLines();
    
    for (Line line : lines) {
        Vec3D p1 = line.getStart();
        Vec3D p2 = line.getEnd();

        // Check if point is on the segment (with scaled epsilon)
        if (isPointOnSegment(pt, p1, p2)) {
            return false;
        }

        // Determine if the edge crosses the horizontal ray from the point to the right
        boolean p1Below = p1.y <= pt.y;
        boolean p2Below = p2.y <= pt.y;

        if (p1Below != p2Below) {
            // 2D cross product (using z-component since all points are on z=0)
            float cross = (p1.x - pt.x) * (p2.y - pt.y) - (p1.y - pt.y) * (p2.x - pt.x);
            
            // Increment/decrement winding number based on edge direction
            if ((p2Below && cross > 1e-6) || (!p2Below && cross < -1e-6)) {
                windingNumber += p2Below ? 1 : -1;
            }
        }
    }

    // Non-zero winding number means the point is inside
    return windingNumber != 0;
}

// Helper for point-on-segment check with relative epsilon
private boolean isPointOnSegment(Vec3D pt, Vec3D p1, Vec3D p2) {
    // Check collinearity with scaled epsilon
    float cross = (p1.x - pt.x) * (p2.y - pt.y) - (p1.y - pt.y) * (p2.x - pt.x);
    if (Math.abs(cross) > 1e-6 * Math.max(Math.abs(p1.x), Math.abs(p2.x))) {
        return false;
    }

    // Check if point is within the segment's bounding box
    float minX = Math.min(p1.x, p2.x);
    float maxX = Math.max(p1.x, p2.x);
    float minY = Math.min(p1.y, p2.y);
    float maxY = Math.max(p1.y, p2.y);

    float epsilon = 1e-6 * Math.max(maxX - minX, maxY - minY);
    return pt.x >= minX - epsilon && pt.x <= maxX + epsilon &&
           pt.y >= minY - epsilon && pt.y <= maxY + epsilon;
}

This version uses cross products (numerically stable) to count how many times the polygon winds around the point. The epsilon values are scaled to the segment's size, making them more reliable than a global threshold.

2. Adjust Your Angle-Based Logic to Respect Winding Number's Integer Nature

If you want to keep your existing angle approach, remember that the total accumulated angle should be a multiple of 2π (full rotations). Instead of checking count > eps, check if the absolute value of the winding number (total angle divided by 2π) is greater than 0.5:

float windingNumber = count / (2 * (float) Math.PI);
return Math.abs(windingNumber) > 0.5f;

This ignores tiny accumulated errors that don't add up to half a full rotation, which is exactly what's happening in your test case.

3. Normalize Coordinates to Eliminate Floating-Point Errors

If your coordinates have a fixed number of decimal places, scale them to integers (e.g., multiply by 1e6 and cast to long). Integer operations are exact, so this eliminates floating-point drift entirely. You’ll need to implement integer-based vector operations, but it’s a permanent fix for precision issues.

4. Use Relative Epsilon Instead of Absolute

Replace your global eps with a value scaled to the magnitude of your coordinates. For example:

float maxCoord = getMaxCoordinateMagnitude();
float eps = 1e-6 * maxCoord;

This ensures the threshold adapts to large or small coordinate systems, avoiding false positives/negatives due to mismatched scales.

Why Your Test Case Fails #

In your example, the point is far outside the polygon, but tiny floating-point errors in angle calculations are adding up to just over your 0.001f threshold. Switching to the cross-product method or adjusting the winding number check will fix this immediately.

内容的提问来源于stack exchange，提问作者leonidasarch