georust / geo

Geospatial primitives and algorithms for Rust
https://crates.io/crates/geo
Other
1.51k stars 195 forks source link

Incorrect convex hull for nearly collinear points #1028

Closed YikChingTsui closed 1 year ago

YikChingTsui commented 1 year ago

Consider these points that roughly form a triangle-like polygon: image

Points A, B, C are the 3 points on the top right, from right to left. Point A is the topmost point.

// using geo v0.25.1
let points = geo::polygon![
    geo::coord! {
        x: 813580.3_f32,
        y: 817846.6,
    },
    // A
    geo::coord! {
        x: 813544.4,
        y: 817952.9,
    },
    // B
    geo::coord! {
        x: 813529.3,
        y: 817947.8,
    },
    // C
    geo::coord! {
        x: 813517.7,
        y: 817943.9,
    },
    geo::coord! {
        x: 813462.3,
        y: 817923.8,
    },
    geo::coord! {
        x: 813290.44,
        y: 817857.06,
    },
];

let expected = points.clone();
let ch = points.convex_hull();
for p in ch.exterior() {
    eprintln!("{},{}", p.x, p.y)
}
assert_eq!(ch, expected);

The convex hull has all points of the input polygon. But Point B is actually slightly lower than the line from C to A. Therefore, the convex hull should go directly from C to A, skipping point B.

I am indeed using a particular CRS, hence why these numbers look unnatural. But plotting it on a "normal" 2D plane, it still appears that point B is unneeded for the convex hull. You can paste the printed output into desmos:

813580.3,817846.6
813544.4,817952.9
813529.3,817947.8
813517.7,817943.9
813462.3,817923.8
813290.44,817857.06
813580.3,817846.6

Furthermore, JTS correctly returns the answer as POLYGON ((813580.3 817846.6, 813290.44 817857.06, 813462.3 817923.8, 813517.7 817943.9, 813544.4 817952.9, 813580.3 817846.6)). So this shouldn't be a CRS issue.

ConvexHull c = new ConvexHull(new Coordinate[] {
                new Coordinate(813580.3, 817846.6),
                new Coordinate(813580.3, 817846.6),
                new Coordinate(813544.4, 817952.9),
                new Coordinate(813529.3, 817947.8),
                new Coordinate(813517.7, 817943.9),
                new Coordinate(813462.3, 817923.8),
                new Coordinate(813290.44, 817857.06),
}, new GeometryFactory());
var geom = c.getConvexHull();
System.out.println(geom);

The last two points (leftmost two points) can be removed without affecting the result and this bug

rmanoka commented 1 year ago

Can you try with f64 precision? The points may be too large for f32 to correctly work. You can also try subtracting out the first point, before calling the algo to see if that fixes it.

YikChingTsui commented 1 year ago

Thanks for the suggestion, f64 works correctly. Interestingly, 32-bit precision still fails even after subtracting the points so that the leftmost point is at the origin. Not the first time I've encountered floating point issues. Hopefully I won't get an edge case where 64-bit fails too, otherwise I could consider a rational/fraction type instead of floats

rmanoka commented 1 year ago

Yes, generally, floating-point numbers has some caveats when converting between display and internal representation.

For instance, 0.1_f64 + 0.2 == 0.3 evaluates to false. This is because, 0.1 is not representible in f64 (or f32), and it's approximated. Now, the degree of this approximation depends on the float type, and in your input I believe the extra point is approximated to fall outside the "true conv hull" when using f32.

If you want to use fixed point/precision, please scale and convert to integers and run the hull (the crate supports hull on i32 or i64 points).

rmanoka commented 1 year ago

Closing; the output is correct for the given input (and data type).