Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

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

Open
pablosanjose opened this issue Jul 30, 2024 · 7 comments · May be fixed by #195
Open

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

pablosanjose opened this issue Jul 30, 2024 · 7 comments · May be fixed by #195
Labels
clipping About polygon clipping or processing documentation Improvements or additions to documentation

Comments

@pablosanjose
Copy link

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.

Screenshot 2024-07-30 at 15 21 45
@asinghvi17 asinghvi17 added bug Something isn't working clipping About polygon clipping or processing labels Jul 31, 2024
@asinghvi17
Copy link
Member

asinghvi17 commented Jul 31, 2024 •
edited
Loading

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

iTerm2 Jj28s7

@asinghvi17
Copy link
Member

asinghvi17 commented Jul 31, 2024 •
edited
Loading

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

iTerm2 b55qob

@asinghvi17 asinghvi17 added documentation Improvements or additions to documentation and removed bug Something isn't working labels Jul 31, 2024
@asinghvi17 asinghvi17 changed the title Intersection bug: "Clipping tracing hit every point" Intersection bug: "Clipping tracing hit every point" when polygons are not closed Jul 31, 2024
@rafaqz
Copy link
Member

rafaqz commented Aug 1, 2024

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
Copy link
Author

pablosanjose commented Aug 1, 2024 •
edited
Loading

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
Copy link
Member

rafaqz commented Aug 1, 2024 •
edited
Loading

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
Copy link
Author

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
Copy link
Member

rafaqz commented Aug 1, 2024

It also uses exact arithmetic :)

This was referenced Aug 3, 2024
Open
Open
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
clipping About polygon clipping or processing documentation Improvements or additions to documentation
Projects
None yet
Milestone
No milestone
Development

Successfully merging a pull request may close this issue.

Use lazy iterator to ensure all rings in create_[a/b]_list are closed JuliaGeo/GeometryOps.jl
3 participants
@rafaqz @pablosanjose @asinghvi17