Implementing weighted restricted delaunay triangulation

[manas-avi]

Hi Bruno,

I am having a lot of fun using this library to try and implement different permutations and combinations of Voronoi diagrams in different settings. It is super cool!

I was wondering if it is possible to use the existing codebase and data structures to implement weighted Restricted Delaunay triangulation (rdt). If I understand correctly, computation of rdt requires the pointer to Delaunay (Delaunay) to be able to re-castable to (GEO::Delaunay_NearestNeighbors) but if I have (GEO::PeriodicDelaunay3d) which implements weighted Delaunay triangulation I cannot recast it into (GEO::Delaunay_NearestNeighbors).

Here: https://github.com/BrunoLevy/geogram/blob/45ca0a1ac7b5bddab4b94ae20175c5b3cbefa741/src/lib/geogram/voronoi/generic_RVD.h#L129

Minimal sample code:

GEO::mesh_load(filename, M_in); GEO::SmartPointer<GEO::PeriodicDelaunay3d> delaunay_ = new GEO::PeriodicDelaunay3d(false, 0.0); GEO::RestrictedVoronoiDiagram_var RVD; RVD = GEO::RestrictedVoronoiDiagram::create(delaunay_,&M_in); GEO::Mesh surface;

// assuming points and weights are loaded in their respective arrays. delaunay_->set_vertices(points_.size()/3, points_.data()); delaunay_->set_weights(weights_.data());

RVD->compute_RDT(surface, mode); // crashes as Delaunay does not have the required neighborhood information as it fails to typecast (GEO::PeriodicDelaunay3d) to (GEO::Delaunay_NearestNeighbors).

Let me know if I can provide some other details. Thanks in advance!

Best, Manas

[BrunoLevy]

Hi, You can use RegularWeightedDelaunay3d, it has several drawbacks:

• it does not work in parallel
• it does not support periodic boundaries
• it has a weird API: both positions and weights must be provided as an array of 4D points, where the 4th coordinates correspond to sqrt(W_max - w_i) where W_max is the max of all the weights and where w_i is the ith-weight

but it is compatible with RVD.

An example is available in the Optimal Transport solver of exploragram, see for instance optimal_transport_on_surface.h and optimal_transport_on_surface.cpp.

For installing and compiling exploragram, see instructions here

[manas-avi]

Hi Bruno,

Thanks for the pointers! I got it working.

I ended up using Delaunay::create(dimension_, "NN");` where the dimension is 3 + (1 for weights) which was compatible with RVD instead of RegularWeightedDelaunay3d (it was not compatible with RVD).

https://github.com/BrunoLevy/exploragram/blob/151764f0764069c1277193e210b62e9c1596f17f/optimal_transport/optimal_transport_on_surface.cpp#L412

The weights were set up as was done here: https://github.com/BrunoLevy/exploragram/blob/151764f0764069c1277193e210b62e9c1596f17f/optimal_transport/optimal_transport.cpp#L607

A follow-up question: Just out of curiosity, is this implementation based on this paper? If not can you point me to the correct reference :)

Thanks in advance! Best, Manas

[BrunoLevy]

Hi Manas,

Yes, the one that you use is similar to this reference (but this reference it is very specific to the GPU), but the one that you use is not the fastest one. Use Delaunay::create("BPOW") to create an instance of WeightedRegularTriangulation, that implements Bowyer-Watson's algorithm (fastest). The one that you use (NN) uses nearest neighbors and convex clipping with the sqrt(Wmax - W) embedding and my "security radius" theorem. It works, but when the weights are very different, it can take up to O(N^2) to compute a single Laguerre diagram.

The algorithm is the result of three incremental steps, described in these articles:

• Voronoi Parallel Linear Enumeration, Lévy and Bonnel, International Meshing Roundtable
• A numerical algorithm for L2 optimal transport, B. Lévy, 2015, ESAIM Math. Modeling and Analysis
• Notions of optimal transport theory and how to implement them on a computer, B. Lévy and E. Schwindt, 2018, Computer and Graphics J
• A fast semi-discrete optimal transport algorithm for a unique reconstruction of the early Universe , B. Lévy, Roya Mohayaee and Sebastian von Hausegger, 2021, Monthly Notices of the Royal Astronomy Society

(but if you read only the third one you will have most of the idea)

Note: the reason why it has this weird API (weighted 3D = 4D, with 4th coordinate as sqrt(weight_max - weight)) is because when I first got interested in transport, I did not have a code for Laguerre diagrams, but I noticed this relation between (d+1) clipped Voronoi and Laguerre, so I used it (because I had a 4D Voronoi code). Now it is not a super good idea to have this API, because it is (probably) a bit confusing for the users, so I will redesign it sometime (but this will break compatibility). The class "PeriodicDelaunay3d" has the new API, standard, with 3D points and weights in a separate vector.

[manas-avi]

Thanks a lot for the advice! I will check the papers out :)