Panic in bool-ops intersection

[RobWalt]

I get a panic for some boolean operations.

thread 'tests' panicked at 'called `Option::unwrap()` on a `None` value', /home/robw/.cargo/registry/src/github.com-1ecc6299db9ec823/geo-0.21.0/src/algorithm/sweep/active.rs:50:37
stack backtrace:
   0: rust_begin_unwind
             at /rustc/fe5b13d681f25ee6474be29d748c65adcd91f69e/library/std/src/panicking.rs:584:5
   1: core::panicking::panic_fmt
             at /rustc/fe5b13d681f25ee6474be29d748c65adcd91f69e/library/core/src/panicking.rs:143:14
   2: core::panicking::panic
             at /rustc/fe5b13d681f25ee6474be29d748c65adcd91f69e/library/core/src/panicking.rs:48:5
   3: core::option::Option<T>::unwrap
             at /rustc/fe5b13d681f25ee6474be29d748c65adcd91f69e/library/core/src/option.rs:755:21
   4: <geo::algorithm::sweep::active::Active<T> as core::cmp::Ord>::cmp
             at /home/robw/.cargo/registry/src/github.com-1ecc6299db9ec823/geo-0.21.0/src/algorithm/sweep/active.rs:50:9
   5: alloc::collections::btree::search::<impl alloc::collections::btree::node::NodeRef<BorrowType,K,V,Type>>::find_key_index
             at /rustc/fe5b13d681f25ee6474be29d748c65adcd91f69e/library/alloc/src/collections/btree/search.rs:211:19
   6: alloc::collections::btree::search::<impl alloc::collections::btree::node::NodeRef<BorrowType,K,V,Type>>::find_upper_bound_index
             at /rustc/fe5b13d681f25ee6474be29d748c65adcd91f69e/library/alloc/src/collections/btree/search.rs:266:45
   7: alloc::collections::btree::search::<impl alloc::collections::btree::node::NodeRef<BorrowType,K,V,alloc::collections::btree::node::marker::LeafOrInternal>>::search_tree_for_bifurcation
             at /rustc/fe5b13d681f25ee6474be29d748c65adcd91f69e/library/alloc/src/collections/btree/search.rs:119:26
   8: alloc::collections::btree::navigate::<impl alloc::collections::btree::node::NodeRef<BorrowType,K,V,alloc::collections::btree::node::marker::LeafOrInternal>>::find_leaf_edges_spanning_range
             at /rustc/fe5b13d681f25ee6474be29d748c65adcd91f69e/library/alloc/src/collections/btree/navigate.rs:258:15
   9: alloc::collections::btree::navigate::<impl alloc::collections::btree::node::NodeRef<alloc::collections::btree::node::marker::Immut,K,V,alloc::collections::btree::node::marker::LeafOrInternal>>::range_search
             at /rustc/fe5b13d681f25ee6474be29d748c65adcd91f69e/library/alloc/src/collections/btree/navigate.rs:308:18
  10: alloc::collections::btree::map::BTreeMap<K,V>::range
             at /rustc/fe5b13d681f25ee6474be29d748c65adcd91f69e/library/alloc/src/collections/btree/map.rs:1077:28
  11: alloc::collections::btree::set::BTreeSet<T>::range
             at /rustc/fe5b13d681f25ee6474be29d748c65adcd91f69e/library/alloc/src/collections/btree/set.rs:298:23
  12: <alloc::collections::btree::set::BTreeSet<geo::algorithm::sweep::active::Active<T>> as geo::algorithm::sweep::active::ActiveSet>::previous
             at /home/robw/.cargo/registry/src/github.com-1ecc6299db9ec823/geo-0.21.0/src/algorithm/sweep/active.rs:73:9
  13: geo::algorithm::sweep::proc::Sweep<C>::handle_event
             at /home/robw/.cargo/registry/src/github.com-1ecc6299db9ec823/geo-0.21.0/src/algorithm/sweep/proc.rs:67:20
  14: geo::algorithm::sweep::proc::Sweep<C>::handle_event
             at /home/robw/.cargo/registry/src/github.com-1ecc6299db9ec823/geo-0.21.0/src/algorithm/sweep/proc.rs:95:40
  15: geo::algorithm::sweep::proc::Sweep<C>::next_event::{{closure}}
             at /home/robw/.cargo/registry/src/github.com-1ecc6299db9ec823/geo-0.21.0/src/algorithm/sweep/proc.rs:45:13
  16: core::option::Option<T>::map
             at /rustc/fe5b13d681f25ee6474be29d748c65adcd91f69e/library/core/src/option.rs:912:29
  17: geo::algorithm::sweep::proc::Sweep<C>::next_event
             at /home/robw/.cargo/registry/src/github.com-1ecc6299db9ec823/geo-0.21.0/src/algorithm/sweep/proc.rs:43:9
  18: <geo::algorithm::sweep::iter::CrossingsIter<C> as core::iter::traits::iterator::Iterator>::next
             at /home/robw/.cargo/registry/src/github.com-1ecc6299db9ec823/geo-0.21.0/src/algorithm/sweep/iter.rs:149:26
  19: geo::algorithm::bool_ops::op::Op<T>::sweep
             at /home/robw/.cargo/registry/src/github.com-1ecc6299db9ec823/geo-0.21.0/src/algorithm/bool_ops/op.rs:65:30
  20: <geo_types::polygon::Polygon<T> as geo::algorithm::bool_ops::BooleanOps>::boolean_op
             at /home/robw/.cargo/registry/src/github.com-1ecc6299db9ec823/geo-0.21.0/src/algorithm/bool_ops/mod.rs:57:21
  21: geo::algorithm::bool_ops::BooleanOps::intersection
             at /home/robw/.cargo/registry/src/github.com-1ecc6299db9ec823/geo-0.21.0/src/algorithm/bool_ops/mod.rs:29:9
  ...

For some very basic shapes:

Shape_A:
   Line { start: Coordinate { x: 17.724912058920285, y: -16.37118892052372 }, end: Coordinate { x: 18.06452454246989, y: -17.693907532504 } }
    Line { start: Coordinate { x: 18.06452454246989, y: -17.693907532504 }, end: Coordinate { x: 19.09389292605319, y: -17.924001641855178 } }
    Line { start: Coordinate { x: 19.09389292605319, y: -17.924001641855178 }, end: Coordinate { x: 17.724912058920285, y: -16.37118892052372 } }
Shape_B:
   Line { start: Coordinate { x: 17.576085274796423, y: -15.791540153598898 }, end: Coordinate { x: 17.19432983818328, y: -17.499393422066746 } }
    Line { start: Coordinate { x: 17.19432983818328, y: -17.499393422066746 }, end: Coordinate { x: 18.06452454246989, y: -17.693907532504 } }
    Line { start: Coordinate { x: 18.06452454246989, y: -17.693907532504 }, end: Coordinate { x: 17.576085274796423, y: -15.791540153598898 } }

I use T=f64 and the error is probably related to some overflow/underflow. A very similar thing happens in #531 at the same place in the code although there some other algorithm uses the function that panics.

Here is a short snippet of code to reproduce the panic with the values from above:

use geo::{BooleanOps, Coordinate as GeoCoord, LineString, Polygon};

type Coordinate = GeoCoord<f64>;

fn main() {
    let ls1_coords = vec![
        Coordinate {
            x: 17.724912058920285,
            y: -16.37118892052372,
        },
        Coordinate {
            x: 18.06452454246989,
            y: -17.693907532504,
        },
        Coordinate {
            x: 19.09389292605319,
            y: -17.924001641855178,
        },
        Coordinate {
            x: 17.724912058920285,
            y: -16.37118892052372,
        },
    ];
    let p1 = Polygon::new(LineString::new(ls1_coords), vec![]);

    let ls2_coords = vec![
        Coordinate {
            x: 17.576085274796423,
            y: -15.791540153598898,
        },
        Coordinate {
            x: 17.19432983818328,
            y: -17.499393422066746,
        },
        Coordinate {
            x: 18.06452454246989,
            y: -17.693907532504,
        },
        Coordinate {
            x: 17.576085274796423,
            y: -15.791540153598898,
        },
    ];
    let p2 = Polygon::new(LineString::new(ls2_coords), vec![]);

    p1.intersection(&p2);
}

[RobWalt]

The error seems to occur due to the edge case of overlap only on one line which is shared by both polygons and the individual lines having different WindingOrders.

[rmanoka]

The issue seems to be due to finite-precision arithmetic. The line-segments that seemingly overlap, are apparently not overlapping. Refer l1 x l2 in snippet below (they intersect at a single point as opposed to an overlap).

Now, this should imply that l3 and l1 do not intersect at the starting point of l3 (because l3 and l2 are adjacent segments of input polygon A), but unfortunately, the loss in precision while calculating intersection ends up disagreeing on this. This leads to topological inconsistencies in the algorithm.

fn snip() {
    let l1 = Line::new(
        coord!(x: 17.576085274796423, y: -15.791540153598898),
        coord!(x: 18.06452454246989, y: -17.693907532504),
    );

    let l2 = Line::new(
        coord!(x: 17.724912058920285, y: -16.37118892052372),
        coord!(x: 18.06452454246989, y: -17.693907532504),
    );
    let l3 = Line::new(
        coord!(x: 17.724912058920285, y: -16.37118892052372),
        coord!(x: 19.09389292605319, y: -17.924001641855178),
    );

    let i12 = line_intersection(l1, l2);
    eprintln!("l1 x l2 = {i12:?}");
    let i13 = line_intersection(l1, l3);
    eprintln!("l1 x l3 = {i13:?}");
}

The output is:

l1 x l2 = Some(SinglePoint { intersection: Coordinate { x: 18.06452454246989, y: -17.693907532504 }, is_proper: false })
l1 x l3 = Some(SinglePoint { intersection: Coordinate { x: 17.724912058920285, y: -16.37118892052372 }, is_proper: true })

We have been trying to fix these issues using f64::nextafter. Essentially, we act like the intersection of l1 x l3 is at the start of l3 + a small epsilon. Adding this fix does solve this issue, but I am yet to wrap my head around the correctness of the fix.

[RobWalt]

I investigated the issue a bit further and found the following. I don't how closely this is related to the problem you were talking about.

So I was curious how the line_intersection function works. There are some orientation checks in there so I looked at them. There we differentiate between the cases

• CounterClockwise
• Clockwise
• Collinear

In your example above we would expect, that we get Collinear for any combination of 3 points of l1 and l2. To be more specific: Let p = l1 and q=l2 in the following snippet (copied from the source code of geo directly)

In line_intersection.rs:

    let p_q1 = RobustKernel::orient2d(p.start, p.end, q.start);
    let p_q2 = RobustKernel::orient2d(p.start, p.end, q.end);
    // [...]
    // omitted
    // [...]

    let q_p1 = RobustKernel::orient2d(q.start, q.end, p.start);
    let q_p2 = RobustKernel::orient2d(q.start, q.end, p.end);
    // [...]
    // omitted
    // [...]

    if matches!(
        (p_q1, p_q2, q_p1, q_p2),
        (Collinear, Collinear, Collinear, Collinear)
    ) {
        return collinear_intersection(p, q);
    }

Then we would expect the conditional at the bottom of the snippet to exectue. However, it doesn't. This comes from floating point precision issues in the RobustKernel::orient2d function. Let's take a look at the critical part:

fn orient2d(p: Coordinate<T>, q: Coordinate<T>, r: Coordinate<T>) -> Orientation {
        // [...]
        // omitted
        // [...]

        if orientation < 0. {
            Orientation::Clockwise
        } else if orientation > 0. {
            Orientation::CounterClockwise
        } else {
            Orientation::Collinear
        }
    }

We get collinear only in the == 0.0 case. This is never a good idea. The values of orientation are respectively:

p_q1 - orientation: -1.2490009027033011e-15
p_q2 - orientation: -0.0
q_p1 - orientation: 1.2212453270876722e-15
q_p2 - orientation: -0.0

I suggest not comparing to 0.0 at this point and using some EPSILON constant instead. I don't know a good heuristic to choose such a value though.

[RobWalt]

Nevermind, scratch that. I tried changing it for my use-case and it only made things worse.

[rmanoka]

@RobWalt FP arithmetic throws many suprises. For instance, the pt (0.1_f64, 1.0) is not collinear with segment (0, 0) - (1, 10) . This is because "0.1" can't be represented by f64 and it actually stores a slightly different value. The orient2d is based on a classical work on Robust Predicates by Shewchuk; this is used by most other geom. packages.