Intersection bug: "Clipping tracing hit every point" when polygons are not closed

[pablosanjose]

Interection of two overlapping polygons with no particular degeneracies fails for some reason

MWE:

julia> import GeometryOps as GO, GeoInterface as GI

julia> using StaticArrays

julia> hex = GI.Polygon([SVector{2}.([[0.012661261919576613, 0.0], [0.006330630959788308, 0.010964974466321872], [-0.006330630959788304, 0.010964974466321874], [-0.012661261919576613, 1.5505573882975772e-18], [-0.006330630959788312, -0.01096497446632187], [0.006330630959788308, -0.010964974466321872]])]);

julia> hex´ = GI.Polygon([SVector{2}.([[0.012382790701850534, 0.0012522257170548476], [0.00548697322789982, 0.011870853442051934], [-0.007156936872252775, 0.011208214192248522], [-0.012905029498454658, -7.305278255197601e-5], [-0.006009212024503947, -0.010691680507549062], [0.006634698075648655, -0.01002904125774565]])]);

julia> int = GO.intersection(hex, hex´; target = GI.PolygonTrait())
ERROR: AssertionError: Clipping tracing hit every point - clipping error. Please open an issue with polygons: [[[0.012661261919576613, 0.0], [0.006330630959788308, 0.010964974466321872], [-0.006330630959788304, 0.010964974466321874], [-0.012661261919576613, 1.5505573882975772e-18], [-0.006330630959788312, -0.01096497446632187], [0.006330630959788308, -0.010964974466321872]]] and [[[0.012382790701850534, 0.0012522257170548476], [0.00548697322789982, 0.011870853442051934], [-0.007156936872252775, 0.011208214192248522], [-0.012905029498454658, -7.305278255197601e-5], [-0.006009212024503947, -0.010691680507549062], [0.006634698075648655, -0.01002904125774565]]].
Stacktrace:
 [1] _trace_polynodes(::Type{…}, a_list::Vector{…}, b_list::Vector{…}, a_idx_list::Vector{…}, f_step::typeof(GeometryOps._inter_step), poly_a::GeoInterface.Wrappers.Polygon{…}, poly_b::GeoInterface.Wrappers.Polygon{…})
   @ GeometryOps ~/.julia/packages/GeometryOps/FMWuP/src/methods/clipping/clipping_processor.jl:572
 [2] _intersection(::GeometryOps.TraitTarget{…}, ::Type{…}, ::GeoInterface.PolygonTrait, poly_a::GeoInterface.Wrappers.Polygon{…}, ::GeoInterface.PolygonTrait, poly_b::GeoInterface.Wrappers.Polygon{…}; exact::GeometryOps._True, kwargs::@Kwargs{})
   @ GeometryOps ~/.julia/packages/GeometryOps/FMWuP/src/methods/clipping/intersection.jl:74
 [3] _intersection
   @ ~/.julia/packages/GeometryOps/FMWuP/src/methods/clipping/intersection.jl:63 [inlined]
 [4] #intersection#152
   @ ~/.julia/packages/GeometryOps/FMWuP/src/methods/clipping/intersection.jl:45 [inlined]
 [5] intersection
   @ ~/.julia/packages/GeometryOps/FMWuP/src/methods/clipping/intersection.jl:42 [inlined]
 [6] top-level scope
   @ REPL[7]:1
Some type information was truncated. Use `show(err)` to see complete types.

[asinghvi17]

I can't see float64 issues being a problem here. union has the same issue and difference seems to give an incorrect answer...

~cc @skygering ~

import GeometryOps as GO, GeoInterface as GI
using GeoMakie, CairoMakie, StaticArrays
hex = GI.Polygon([SVector{2}.([[0.012661261919576613, 0.0], [0.006330630959788308, 0.010964974466321872], [-0.006330630959788304, 0.010964974466321874], [-0.012661261919576613, 1.5505573882975772e-18], [-0.006330630959788312, -0.01096497446632187], [0.006330630959788308, -0.010964974466321872]])]);

hex´ = GI.Polygon([SVector{2}.([[0.012382790701850534, 0.0012522257170548476], [0.00548697322789982, 0.011870853442051934], [-0.007156936872252775, 0.011208214192248522], [-0.012905029498454658, -7.305278255197601e-5], [-0.006009212024503947, -0.010691680507549062], [0.006634698075648655, -0.01002904125774565]])]);

dif = GO.difference(hex, hex´; target = GI.PolygonTrait())

f, a, p = poly(hex; color = Cycled(1), label = "hex")
poly!(hex´; color = Cycled(2), label = "hex´")
poly!(dif; color = Cycled(3), label = "difference")
axislegend(a)
f

[asinghvi17]

Ah, the polygons aren't closed (i.e., poly[end] == poly[begin]) @pablosanjose. You can run hex = GO.fix(hex); hex´ = GO.fix(hex´) and this will work fine.

We do need a better error message if the rings are not closed, though. Maybe that can go in trace_polynodes? Or we could fix it internally...

hex = GI.Polygon([SVector{2}.([[0.012661261919576613, 0.0], [0.006330630959788308, 0.010964974466321872], [-0.006330630959788304, 0.010964974466321874], [-0.012661261919576613, 1.5505573882975772e-18], [-0.006330630959788312, -0.01096497446632187], [0.006330630959788308, -0.010964974466321872]])]) |> GO.fix;

hex´ = GI.Polygon([SVector{2}.([[0.012382790701850534, 0.0012522257170548476], [0.00548697322789982, 0.011870853442051934], [-0.007156936872252775, 0.011208214192248522], [-0.012905029498454658, -7.305278255197601e-5], [-0.006009212024503947, -0.010691680507549062], [0.006634698075648655, -0.01002904125774565]])]) |> GO.fix;

int = GO.intersection(hex, hex´; target = GI.PolygonTrait())

f, a, p = poly(hex; color = Cycled(1), label = "hex")
poly!(hex´; color = Cycled(2), label = "hex´")
poly!(int; color = Cycled(3), label = "intersection")
axislegend(a)
f

[rafaqz]

Yes we can probably fix it on the fly if the begin and end points don't match in any ring. It would be good if unclosed rings just worked everywhere without manual fixes.

[pablosanjose]

So, this was the problem! Thanks.

Is it explained anywhere that polygons "should" be closed (identical start and end points)? Note that intersection in Meshes.jl does not have this problem, but in GeometryOps it very much makes the assumption that polys are closed without checking. If this is an intrinsic assumption in the overall design, it should be guaranteed by the contructors themselves, no? It's not just a documentation issue.

[rafaqz]

This package is very new, just a bug not intentional.

We accept geoms from basically any package so we can't know if they're closed. We just have to check.

[pablosanjose]

Of course, sorry! Actually, GeometryOps is surprisingly good already. In this particular intersection problem it is faster and allocates way less than Meshes.jl. Keep up the great work!

[rafaqz]

It also uses exact arithmetic :)