Bridge returns NaNs for PolyArea with intersecting Rings

[Ethan-Russell]

I was trying to do a geojoin involving United States census block level data, and ran into an error. I tracked it down to a specific geometry, where calling simplexify would introduce NaN into the vertices. I have created a minimum working example below.

julia> ring1 = Ring(
               Point(0., 0.),
               Point(0., 3.),
               Point(2., 3.),
               Point(2., 2.),
               Point(3., 2.),
               Point(3., 0.)
           );

julia> ring2 = Ring(
               Point(1., 1.),
               Point(1., 2.),
               Point(2., 2.),
               Point(2., 1.)
           );

julia> pa = PolyArea(ring1, ring2);

julia> simplexify(pa)
10 SimpleMesh{2,Float64}
  10 vertices
  ├─ Point(0.0, 0.0)
  ├─ Point(0.0, 3.0)
  ├─ Point(2.0, 3.0)
  ├─ Point(NaN, NaN)
  ├─ Point(NaN, NaN)
  ├─ Point(1.0, 2.0)
  ├─ Point(1.0, 1.0)
  ├─ Point(2.0, 1.0)
  ├─ Point(3.0, 2.0)
  └─ Point(3.0, 0.0)
  10 elements
  ├─ Triangle(7, 1, 2)
  ├─ Triangle(6, 2, 3)
  ├─ Triangle(5, 6, 3)
  ├─ Triangle(3, 4, 5)
  ├─ Triangle(6, 7, 2)
  ├─ Triangle(10, 1, 7)
  ├─ Triangle(10, 7, 8)
  ├─ Triangle(10, 8, 5)
  ├─ Triangle(9, 10, 5)
  └─ Triangle(5, 4, 9)

Note the NaN's in the list of vertices above.

Here is a visualization of the two rings composing the PolyArea:

[juliohm]

Thank you Ethan for reporting the issue. Can you please apply the Bridge transform to this geometry? The issue is probably there. Em qui., 9 de nov. de 2023 14:02, Ethan Russell ***@***.***> escreveu:

…

I was trying to do a geojoin involving United States census block level data, and ran into an error. I tracked it down to a specific geometry, where calling simplexify would introduce NaN into the vertices. I have created a minimum working example below. julia> ring1 = Ring( Point(0., 0.), Point(0., 3.), Point(2., 3.), Point(2., 2.), Point(3., 2.), Point(3., 0.) ); julia> ring2 = Ring( Point(1., 1.), Point(1., 2.), Point(2., 2.), Point(2., 1.) ); julia> pa = PolyArea(ring1, ring2); julia> simplexify(pa)10 SimpleMesh{2,Float64} 10 vertices ├─ Point(0.0, 0.0) ├─ Point(0.0, 3.0) ├─ Point(2.0, 3.0) ├─ Point(NaN, NaN) ├─ Point(NaN, NaN) ├─ Point(1.0, 2.0) ├─ Point(1.0, 1.0) ├─ Point(2.0, 1.0) ├─ Point(3.0, 2.0) └─ Point(3.0, 0.0) 10 elements ├─ Triangle(7, 1, 2) ├─ Triangle(6, 2, 3) ├─ Triangle(5, 6, 3) ├─ Triangle(3, 4, 5) ├─ Triangle(6, 7, 2) ├─ Triangle(10, 1, 7) ├─ Triangle(10, 7, 8) ├─ Triangle(10, 8, 5) ├─ Triangle(9, 10, 5) └─ Triangle(5, 4, 9) Note the NaN's in the list of vertices above. Here is a visualization of the two rings composing the PolyArea: [image: mwe] <https://user-images.githubusercontent.com/44298283/281810464-28f1210d-94e2-4f4a-9d1d-1cb18f68878a.png> — Reply to this email directly, view it on GitHub <#629>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAZQW3PWBYWCNQLXHTHMMRLYDUED5AVCNFSM6AAAAAA7E4QP5CVHI2DSMVQWIX3LMV43ASLTON2WKOZRHE4DMMBVGI4TMNY> . You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

[juliohm]

Also, could you please confirm the predicate you are using in the geojoin?

[juliohm]

To clarify what is happening, the simplexify algorithm will discretize the polygon with holes in three steps:

1. Connect the holes with bridges using the Bridge transform
2. Run discretization algorithm such as Dehn1898 or FIST
3. Undo bridges to reintroduce holes

The Bridge transform is currently implemented assuming that the inner rings don't touch the outer ring. Perhaps we need a fix there to accommodate this corner case.

[Ethan-Russell]

From bridge.jl

function apply(transform::Bridge, poly::Polygon{Dim,T}) where {Dim,T}
  # sort rings lexicographically
  rpoly, rinds = apply(Repair{9}(), poly)

  # retrieve bridge width
  δ = isnothing(transform.δ) ? zero(T) : transform.δ

  ring, dups = if hasholes(rpoly)
    bridge(rings(rpoly), rinds, δ)
  else
    first(rings(rpoly)), []
  end

  PolyArea(ring), dups
end

So you were right to suspect the bridge. The NaNs are introduced in the bridge call (since rpoly has holes)

[juliohm]

One possible fix would be to "enlarge" the outer ring with a fixed epsilon value, perform the same operations already implemented in the bridge, and then undo the "enlarge". This could be implemented inside the function you copied/pasted above. Would you like to give it a try?

This fix could be implemented in the function above or inside the bridge function itself.

[Ethan-Russell]

I'm looking in to a fix for the corner case that would include the intersecting point in the bridge - do you think that could work? I will create a fork if I am able to get it working.

[juliohm]

I think the easiest solution will be a variation of what I described where we first "enlarge" the outer ring to make sure that it doesn't touch/intersect the inner rings. That way we can run the algorithm safely, and then undo the "enlargement" in a post-processing step.

This enlargement can be as easy as a Translate to the origin, followed by a Stretch, then undo the translation. Or you could even just move the touching point outwards using the normal vectors nearby.

[juliohm]

@Ethan-Russell I've just added the two new transforms as explained above.

The Expand transform is similar to the Stretch transform, but "enlarges" the geometry in place: e1c1ad3

The Repair{10} expands the outer ring of the polygon to handle these edge cases where the outer ring is touching the inner rings: 889b565

Now we just need to figure out how to apply this repair before and after the discretization to make sure that the points of the original polygon remain present in the mesh.