JuliaGeometry / DelaunayTriangulation.jl

DelaunayTriangulation.jl: A Julia package for Delaunay triangulations and Voronoi tessellations in the plane
https://juliageometry.github.io/DelaunayTriangulation.jl/
MIT License
65 stars 6 forks source link

Bad lift to a surface of the Delaunay triangulation of its domain of definition #43

Closed empet closed 1 year ago

empet commented 1 year ago

Meanwhile I solved this issue

DanielVandH commented 1 year ago

What ended up being the problem @empet? If there's anything that was maybe unclear in the documentation for example that could be good for me to fix up.

empet commented 1 year ago

The docs are very good, but in an attempt to reproduce a surface triangulation I worked on a few years ago, with Delaunay.jl I made the following mistake: The surface (x,y) ->f(x,y) was defined on the planar disk, D(0,2). Hence I started with:

r = collect(0.0:0.1:2)
θ = collect(0:π/30:2π)
x = vec(r * cos.(θ)')
y = vec(r * sin.(θ)')
z = f.(x, y)
pts=[[xe, ye] for (xe, ye) in zip(x,y)]
tri=triangulate(pts)

but got:

AssertionError: Duplicate points are not allowed.

Delaunay.jl as well as scipy.spatial.Delaunay work very well with non-unique points. I replaced pts by unique(pts) and continued to transcribe the triangulation code for the former package, in terms of DelaunayTraingulation. But I plotted the triangulated surface using tri.triangles, and the points onto the surface corresponding to initial, non-unique pts, and obviuosly I got an ugly, non-smooth surface. When I clicked to send the formulated issue, I realized what was wrong. So I started with r=collect(0.05:0.1, 2), that led to unique pts and the triangulated surface was perfect. :) Final conclusion: DelaunayTrisngulation.jl is great!

DanielVandH commented 1 year ago

I see - yes duplicate points can be a bit of a hassle to handle, so I don't try and resolve it internally as is done in scipy.spatial.Delaunay (and hence Delaunay.jl). unique! might've saved you some trouble here too since it would just update pts inplace.

I had initially implemented it so that if duplicated points were detected, triangulate was re-called with skip_points set to be the points that were detected as duplicates, but I ended up deciding to leave it to the user to resolve these degeneracies rather than handling it within the package. Proceeding with them in the algorithm is also not good since it can cause some infinite loops.

Thanks for letting me know about the solution! Appreciate the feedback too :).