Robust Winding Order

[rmanoka]

This PR introduces a few traits to abstract the predicates logic for different scalar types.

• Introduces Kernel, and HasKernel traits to allow customizing the predicates used for different scalars (integer vs float)
• Defines RobustKernel meant for float types (f32, f64)
• Defines SimpleKernel meant for exact arithmetic types (i32, i64; can be extended to support arbitrary precision types)

Also port a few key, but simple algos to use the predicates. This improves the correctness of these algos.

• Winding: shift from area logic into using the orient2d logic.
• Convex Hull: provide an alternative graham_hull (Graham Scan implementation) using the predicates. This has the advantage of being more correct, and also works on geoms,with int. coords.

[rmanoka]

Here is another option to allow the user to switch between different implementations of the predicate. It is similar to the spade crate, and also CGAL. We abstract the predicates into a trait Kernel as follows:

pub trait Kernel {
    type Scalar: CoordinateType;

    fn orient2d(
        p: Coordinate<Self::Scalar>,
        q: Coordinate<Self::Scalar>,
        r: Coordinate<Self::Scalar>,
    ) -> Option<WindingOrder>;
}

Then we implement different kernels: struct FloatKernel, struct SimpleKernel, struct UnsignedIntegerKernel etc. For instance, a FloatKernel could be like this:

#[derive(Default)]
pub struct FloatKernel<T>(PhantomData<T>);

impl<T: Float> Kernel for FloatKernel<T> {
    type Scalar = T;

    fn orient2d(
        p: Coordinate<Self::Scalar>,
        q: Coordinate<Self::Scalar>,
        r: Coordinate<Self::Scalar>,
    ) -> Option<WindingOrder> {
        //  implementation using robust::orient2d predicate
    }
}

Then, the Winding trait becomes:

pub trait RobustWinding<K>
where
    K: Kernel
{
    /// Return the winding order of this object
    fn robust_winding_order(&self) -> Option<WindingOrder>;
}

and implementation for LineString:

impl<T: CoordinateType, K: Kernel<Scalar = T>> RobustWinding<K> for LineString<T> {
    fn robust_winding_order(&self) -> Option<WindingOrder> {
        // ....
    }
}

The usage becomes a bit verbose, but we allow the user to choose the Kernel they want to use:

#[test]
fn robust_winding_order() {
    use crate::Point;
    // 3 points forming a triangle
    let a = Point::new(0., 0.);
    let b = Point::new(2., 0.);
    let c = Point::new(1., 2.);

    use super::kernels::FloatKernel;

    // That triangle, but in clockwise ordering
    let cw_line = LineString::from(vec![a.0, c.0, b.0, a.0]);

    let orientation = RobustWinding::<FloatKernel<_>>::robust_winding_order(&cw_line);
}

The added verbosity is unfortunate, but it makes the usage more flexible. If we could go ahead and modify the geometry structures, then we can reduce the verbosity:

pub struct Point<K: Kernel<Scalar = T>, T>(pub Coordinate<T>, pub PhantomData<K>)

Unfortunately, that's too many changes to the existing types.

Please let me know what you think about this.

Edit: slight improvement to usability.

[rmanoka]

I'm outta my league here, but based your earlier comments, I understand that this Kernel would only be appropriate for integer types. Is that right?

It is meant for signed integer, and exact arithmetic types (arbitrary precision arithmetic types like eg. gmp float or (signed) integers).

[rmanoka]

@urschrei Could you opine on the traits design for using different predicates for different types.

Please see HasKernel assignments and the usage syntax for signed integers and for floats; they're actually using different predicate implementations underneath. Should we go this way? I like this, but there are caveats.

A more flexible, but much more verbose approach is explored here.

[rmanoka]

Added tests to verify robustness, and also to visualize the non-robustness of the naive implementation. The tests are ported from the demos in from https://github.com/danshapero/predicates . They take 3 points, modify the middle point at small increments over a grid of values, and output a png map. The output is quite intriguing; had to see it to believe the non-robustness.

The naive orientation map; the output jumps between cw (black), colinear (gray), and ccw(white).

The robust orientation produces a clean output!

To run tests:

cargo test test_naive -- --ignored
cargo test test_robust -- --ignored

[michaelkirk]

The naive orientation map; the output jumps between cw (black), colinear (gray), and ccw(white).
The robust orientation produces a clean output!

That's really satisfying to see.

[rmanoka]

Requesting review / discussion and merging of this PR. Currently, the new code is fairly non-intrusive (separate traits, that could later replace existing ones). Here are the key additions:

1. Kernel trait to abstract the predicate logic. We currently have one naive implementation, used for i32, i64; and the robust one used for f32, f64. Includes fairly comprehensive test of robustness.
2. RobustWinding trait: same interface as Winding trait, but uses the new goodies, and so robust. If all is good, we can replace the existing Winding trait with this.
3. graham_hull algorithm for computing convex_hull (it can be used for integer coords too). Currently not yet a trait, but we may consider using this algo for ConvexHull trait so that it works for all scalar types.

@michaelkirk @frewsxcv @urschrei Please provide feedback.

[michaelkirk]

Nice work!

In the interest of having the best possible defaults, and not imposing unnecessary choice on geo users...

1. Currently, the new code is fairly non-intrusive (separate traits, that could later replace existing ones).

It seems like there are some correctness wins with the new implementations. Is there a benefit to keeping the old ones at all? If not, I'd be in favor of replacing them.

If there is some legit benefit to keeping the old ones, is there a reason we shouldn't default to the new, more accurate, implementations?

(similar concern for graham_hull vs. the existing convex hull impl)

And re: this comment from the PR description:

This PR side-steps the whole trait issue, and just creates a new trait, only meant to be implemented for T: Float. The rationale is to get the logic right, and hope for some expert advice on the trait magic.

Is this comment outdated? Is there still some question around traits outstanding, and if so, could you restate the question?

[rmanoka]

It seems like there are some correctness wins with the new implementations. Is there a benefit to keeping the old ones at all? If not, I'd be in favor of replacing them.

Makes sense; will replace the existing Winding trait as this is a more correct impl.

If there is some legit benefit to keeping the old ones, is there a reason we shouldn't default to the new, more accurate, implementations?

(similar concern for graham_hull vs. the existing convex hull impl)

Currently graham_hull is provided as a free-standing algo, and not as a trait. I didn't want to introduce too many changes in a single PR. Also, I'm not completely sure if they have the same performance in practice, so we should do some benchmarks before switching. It is in this PR to basically show-case the Kernel ideas and to provide the functionality for integer geometries (which wasn't previously possible).

Is this comment outdated? Is there still some question around traits outstanding, and if so, could you restate the question?

Yeah :). Will edit to top comment and update the current status.

[urschrei]

This is really great work, @rmanoka. Apart from my two minor comments, I don't think there's anything else that needs to be added here. There are typos here and there, but let's defer those to a separate PR. I have a strong preference to retain the QuickHull implementation, but it's based purely on sentimentality and I wouldn't want to do so at the expense of the Graham implementation and they're both O(n log n) on average, so ¯_(ツ)_/¯

[rmanoka]

This is really great work, @rmanoka. Apart from my two minor comments, I don't think there's anything else that needs to be added here. There are typos here and there, but let's defer those to a separate PR. I have a strong preference to retain the QuickHull implementation, but it's based purely on sentimentality and I wouldn't want to do so at the expense of the Graham implementation and they're both O(n log n) on average, so ¯_(ツ)_/¯

@urschrei The graham scan algorithm has two added benefits: it uses robust predicates, and also works on integer coords. The quick hull impl. could be made robust (by fixing point_location), but do not see how to get it working on integer coords. That said, I'm okay with leaving graham scan as a stand-alone algo. that users can import and use if needed.

[urschrei]

Noo the point of all this is to integrate robust predicates to any existing traits that can use them, so the convex hull algorithm should be using Graham by default! We can fix up / revisit Quickhull later on if need be.

[rmanoka]

Noo the point of all this is to integrate robust predicates to any existing traits that can use them, so the convex hull algorithm should be using Graham by default! We can fix up / revisit Quickhull later on if need be.

Sounds good! Let's do that in a separate issue / PR; see a few simple improvements, and bug-fixes to the trait.

[rmanoka]

Should we merge this? @michaelkirk @urschrei @frewsxcv

[michaelkirk]

I believe so - thanks for your work @rmanoka !

bors r+

[bors]

Build succeeded:

• coverage
• geo
• geo-postgis
• geo-types