Problematic polygon for `FIST`

[juliohm]

We fixed the random number generator in our tests for reproducibility and uncovered a bug in our FIST discretization.

MWE:

The discretize call below is hanging:

using Meshes
using Random

# helper function to read *.line files containing polygons
# generated with RPG (https://github.com/cgalab/genpoly-rpg)
function readpoly(T, fname)
  open(fname, "r") do f
    # read outer chain
    n = parse(Int, readline(f))
    outer = map(1:n) do _
      coords = readline(f)
      x, y = parse.(T, split(coords))
      Point(x, y)
    end

    # read inner chains
    inners = []
    while !eof(f)
      n = parse(Int, readline(f))
      inner = map(1:n) do _
        coords = readline(f)
        x, y = parse.(T, split(coords))
        Point(x, y)
      end
      push!(inners, inner)
    end

    # return polygonal area
    @assert first(outer) == last(outer)
    @assert all(first(i) == last(i) for i in inners)
    rings = [outer, inners...]
    PolyArea([r[begin:(end - 1)] for r in rings])
  end
end

rng = MersenneTwister(123)

poly = readpoly(Float32, "test/data/hole1.line")

mesh = discretize(poly, FIST(rng))

[juliohm]

Reduced the MWE to a small ring that doesn't satisfy the ear condition in the FIST paper:

julia> v
5-element Vector{Point2f}:
 Point(-0.5f0, 0.3296139f0)
 Point(-0.19128194f0, -0.5f0)
 Point(-0.37872985f0, 0.29592824f0)
 Point(0.21377224f0, -0.0076110554f0)
 Point(-0.20127837f0, 0.24671146f0)

julia> r = Ring(v)
Ring{2,Float32}
├─ Point(-0.5f0, 0.3296139f0)
├─ Point(-0.19128194f0, -0.5f0)
├─ Point(-0.37872985f0, 0.29592824f0)
├─ Point(0.21377224f0, -0.0076110554f0)
└─ Point(-0.20127837f0, 0.24671146f0)

julia> Meshes.earsccw(r)
Int64[]

In that case, the ring enters in the recovery process of

Meshes.jl/src/discretization/fist.jl

Lines 86 to 105 in 2e8148f

	else # recovery process
	# check if consecutive edges vᵢ-1 -- vᵢ and vᵢ+1 -- vᵢ+2
	# intersect and fix the issue by clipping ear (vᵢ, vᵢ+1, vᵢ+2)
	v = vertices(𝒫)
	for i in 1:n
	s1 = Segment(v[i - 1], v[i])
	s2 = Segment(v[i + 1], v[i + 2])
	λ(I) = type(I) == Crossing
	if intersection(λ, s1, s2)
	# 1. push a new triangle to 𝒯
	push!(𝒯, connect((inds[i], inds[i + 1], inds[i + 2])))
	# 2. remove the vertex from 𝒫
	inds = inds[setdiff(1:n, mod1(i + 1, n))]
	𝒫 = Ring(stdpts[inds])
	n = nvertices(𝒫)
	clipped = true
	break
	end
	end
	end

and never exits it. We need to implement the other recovery steps in the paper besides the one already implemented.

[juliohm]

The new recovery process has been added, and the ring above is successfully triangulated. More tests are needed before we reintroduce the hole1 benchmark.