Currently, weighted triangulations do not work, and so we error if the user tries to create one (actually, weighted triangulations of convex polygons do work with triangulate_convex, but anyway). The point of this issue is to (1) explain the theory behind weighted triangulations, (2) the algorithm, and (3) our implementation. I'll finish by giving an actually example of it breaking and explaining all the pieces that are in the code and tests thus far.

Weighted Triangulations

The following explanation is based on Delaunay Mesh Generation by Cheng, Dey, and Shewchuk (2013).

Suppose we have a set of points $\mathcal P \subset \mathbb R$. To each point $p_i \in \mathcal P$ we associate a weight $w_i$ and define the lifted point (or lifted companion) $p_i^+ = (x_i, y_i, x_i^2 + y_i^2 - w_i)$; this $z$ coordinate $x_i^2 + y_i^2 - w_i$ is the distance of the lifted point from the paraboloid $z = x^2 + y^2$, in particular:

$w_i > 0$: The point is below the paraboloid.
$w_i < 0$: The point is above the paraboloid.
$w_i = 0$: The point is on the paraboloid.

The weighted Delaunay triangulation, which I'll denote $\mathcal W\mathcal T(\mathcal P)$, is the projection of the underside of the convex hull of $\mathcal P^+$, the set of all lifted points, onto the plane. For example, here is a MATLAB script that can compute a weighted Delaunay triangulation (I also give this script here https://github.com/JuliaGeometry/DelaunayTriangulation.jl/blob/41ccf2319bd637f92bcda3df817e91c776b83216/test/helper_functions.jl#L1461-L1571) The save_to_file function computes the triangulation and saves it to a .txt file.

function [cells, lifted] = LiftingMap(points,weights)
num = size(points, 1);
lifted = zeros(num, 3);
for i = 1:num
    p = points(i, :);
    lifted(i, :) = [p sum(p.^2) - weights(i)];
end
ch = convhull(lifted);
sep_ch = reshape(lifted(ch, :), [size(ch) 3]);
vecs = diff(sep_ch, 1, 2);
n = cross(vecs(:, 1, :), vecs(:, 2, :), 3);
downward_faces = n(:, :, 3) <= 0;
cells = ch(downward_faces, :);
end

function save_to_file(points, weights, n)
[lift_cells, ~] = LiftingMap(points,weights);
res = [(1+size(points, 1)) * ones(1, 3); [points weights]; lift_cells];
% so the points are in res(2:res(1), 1:2), the weights are in
% res(2:res(1), 3), and the triangles are in res((res(1)+1):end, :).
writematrix(res, fullfile(pwd, 'tempmats', strcat(string(n), '.txt')));
end

One issue with weighted triangulations, and indeed a big part of the reason they don't yet work in this package, is the problem of submerged vertices. Given a vertex $v$, if the associated with $w_v$ is so small that its lifted companion $p_v^+$ is not on the underside of the convex hull of $\mathcal P^+$, then it will not appear in $\mathcal W\mathcal T(\mathcal P)$ at all. In the example below, the fourth vertex appears with a weight of $w_4 = -0.35$ but changing the weight to $w_4=-0.4$ makes it no longer appear in the triangulation, i.e. it becomes a submerged vertex.

points = [
    0.0 0.0
    1.0 0.0
    0.0 1.0
    0.25 0.25];
weights = [
    0.0
    0.0
    0.0
    -0.35
];
tri = LiftingMap(points, weights);
weights(4) = -0.4; % from -0.35 to -0.4
tri2 = LiftingMap(points, weights);
figure;
subplot(1, 2, 1);
triplot(tri, points(:, 1), points(:, 2));
subplot(1, 2, 2);
triplot(tri2, points(:, 1), points(:, 2));

Lastly, note that a weighted triangulation where $w_v \equiv w$ for all $v$, i.e. the weights are all constants, is just the Delaunay triangulation.

Computing Weighted Delaunay Triangulations (In Theory)

In this package, we use the Bowyer-Watson algorithm for computing triangulations. To understand how weighted Delaunay triangulations can be computed, the same algorithm can be used with some slight modifications. It might be helpful to read https://juliageometry.github.io/DelaunayTriangulation.jl/dev/math/delaunay/#Bowyer-Watson-Algorithm for the unweighted case before proceeding. To summarise, the Bowyer-Watson algorithm for inserting a point $u$ into a triangulation is:

Find one triangle whose open circumdisk contains $u$. (Actually, we find a triangle containing $u$ instead of just an open circumdisk, but anyway.)
Find all the others by a depth-first search.
Delete all the triangles containing $u$ in their open circumdisk.
For each edge of the evacuated polyhedral cavity, adjoin the vertices to $u$ to create a new triangle.

Note that this assumes the point $u$ is inside the triangulation. For points outside of the triangulation, we use ghost triangles as described here https://juliageometry.github.io/DelaunayTriangulation.jl/dev/math/delaunay/#Exterior-Points. For weighted triangulations, the rules for working with ghost triangles change slightly: A ghost triangle is deleted if a vertex is inserted in its outer halfplane, except perhaps if the vertex lies on the solid edge of the ghost triangle, in which case the ghost triangle is deleted if the new vertex is not submerged.

For the next modification for weighted triangulations, we need to replace how the incircle predicate is computed. In particular, instead of checking incircle(u, v, w, x) > 0 (see the DigCavity function in the mentioned textbook on p.63), we need to look at orient3(u⁺, v⁺, w⁺, x⁺), which tests whether $x^+$ is below the witness plane of the triangle $uvw$; the witness plane to a non-submerged point is a non-vertical plane in space that contains the triangle on the convex hull of $\mathcal P^+$ such that no other vertex in the convex hull is below the witness plane (see Def 2.5 in the book). It is also important that, when we put an input triangle whose open circumdisk contains $u$ into our algorithm, it must be a triangle whose witness plane is above the inserted vertice's lifted point. If $u$ is submerged, no such triangle exists and so it should not be inserted.

Finally, while the book does mention point location using the jump and march algorithm (as we do in this package), in two dimensions they also mention a conflict graph method (Section 2.6) which seems to have better handling for weighted triangulations. When I emailed Shewchuk previously about this, he also referred back to this conflict graph approach. I would really prefer to be able to avoid implementing this since I'm not sure how well it works with dynamic updates and I think we can get around it.

What's Implemented in the Package

Before I list what's actually not working with my implementation of weighted triangulations, let me list out what I've implemented in the package already (whether there are any bugs or not is to be determined I guess). Most of these functions have corresponding tests you can search for. Some of these functions were implemented so long ago that I don't really remember their purpose. A lot of tests are here https://github.com/JuliaGeometry/DelaunayTriangulation.jl/blob/41ccf2319bd637f92bcda3df817e91c776b83216/test/triangulation/weighted.jl, and some test triangulations are generated in the code at https://github.com/JuliaGeometry/DelaunayTriangulation.jl/blob/41ccf2319bd637f92bcda3df817e91c776b83216/test/helper_functions.jl#L1461-L1663. Some tests in test/triangulation/weighted.jl are broken because triangulate errors when weights are provided, but you can Revise the check_config function to skip over that if you want.

`add_point!(tri, x, y, w)`

https://github.com/JuliaGeometry/DelaunayTriangulation.jl/blob/41ccf2319bd637f92bcda3df817e91c776b83216/src/algorithms/triangulation/basic_operations/add_point.jl#L111-L120

`point_position_relative_to_circumcircle`

This function already uses an is_weighted check to convert to the test of the witness plane orientation.

https://github.com/JuliaGeometry/DelaunayTriangulation.jl/blob/41ccf2319bd637f92bcda3df817e91c776b83216/src/predicates/predicates.jl#L476-L520

`point_position_relative_to_witness_plane`

As it says on the tin.

https://github.com/JuliaGeometry/DelaunayTriangulation.jl/blob/41ccf2319bd637f92bcda3df817e91c776b83216/src/predicates/predicates.jl#L521-L541

`is_submerged`

https://github.com/JuliaGeometry/DelaunayTriangulation.jl/blob/41ccf2319bd637f92bcda3df817e91c776b83216/src/data_structures/triangulation/methods/weights.jl#L156-L177

`get_weighted_nearest_neighbour`

I have no idea why this is needed. Maybe it was for power diagrams later?

https://github.com/JuliaGeometry/DelaunayTriangulation.jl/blob/41ccf2319bd637f92bcda3df817e91c776b83216/src/data_structures/triangulation/methods/weights.jl#L123-L154

`get_distance_to_witness_plane`

As above, I don't know why this is needed, except maybe it's for power diagrams or something.

https://github.com/JuliaGeometry/DelaunayTriangulation.jl/blob/41ccf2319bd637f92bcda3df817e91c776b83216/src/data_structures/triangulation/methods/weights.jl#L90-L121

`get_power_distance`

Must be for power diagrams.

https://github.com/JuliaGeometry/DelaunayTriangulation.jl/blob/41ccf2319bd637f92bcda3df817e91c776b83216/src/data_structures/triangulation/methods/weights.jl#L77-L88

`triangulate_convex`

Triangulating convex polygons with weights actually already works. You can just pass triangulate_convex a set of weights and it'll work. This functions uses a slightly specialised version of the Bowyer-Watson algorithm though, so might not be too much to take from that. This is tested at https://github.com/JuliaGeometry/DelaunayTriangulation.jl/blob/41ccf2319bd637f92bcda3df817e91c776b83216/test/triangulation/weighted.jl#L256-L263. Something that's new.. the tests get stuck at fi = 77 (the last iteration) when using validate_triangulation. I don't know why.. it used to work.

Inside Bowyer-Watson

Also see

https://github.com/JuliaGeometry/DelaunayTriangulation.jl/blob/ee4fa67f169657df665430807c0d31621f632e00/src/algorithms/triangulation/unconstrained_triangulation.jl#L147-L154

A Broken Example

Let's now give an example of a weighted triangulation not being computed correctly. For this, I have locally disabled check_config to get past the errors. Consider

using Random
Random.seed!(123)
points = [0.0 1.0 0.0 0.25 -0.5; 0.0 0.0 1.0 0.25 0.5]
weights = [0.0, 0.0, 0.0, -1.0, 1.0]
tri = triangulate(points; weights, insertion_order = [1, 2, 4, 5, 3])
# force insertion order since there are other issues like KeyErrors or infinite loops sometimes, don't know why. 
# Feel free to run it without forcing the order to see what else can happen.
# It gets it correct for insertion_order = [1, 2, 3, 4, 5]!

The fourth vertex is supposed to be submerged, so that the triangulation we want is

but instead we get

triplot(tri)

The fourth vertex is not only not submerged, but it has a dangling edge attached to it and doesn't form an actual triangle. Look at all the problems this triangulation has:

julia> DelaunayTriangulation.validate_triangulation(tri) # on main
The adjacent2vertex map is inconsistent with the adjacent map. The vertex 5 and edge (2, 3) are inconsistent with the associated edge set Set([(3, -1), (1, 2), (-1, 1), (2, 3)]) from the adjacent2vertex map, as (2, 3) maps to 4 instead of 5.
The triangle (3, 5, 2) is incorrectly stored in the adjacent map, with the edge (2, 3) instead mapping to 4.
The Delaunay criterion does not hold for the triangle-vertex pair ((3, 4, 2), 1).
The triangle (3, 2, 4) is negatively oriented.
The edge iterator has duplicate edges.
The solid_edge iterator has duplicate edges.

What's going wrong here? Why is 4 not being submerged? Let's look at it step by step. Here is code that gets the triangulation one step at a time, using exactly what triangulate does.

Random.seed!(123)
points = [0.0 1.0 0.0 0.25 -0.5; 0.0 0.0 1.0 0.25 0.5]
weights = [0.0, 0.0, 0.0, -1.0, 1.0]
tri = Triangulation(points; weights)
_tri = deepcopy(tri)

DelaunayTriangulation.initialise_bowyer_watson!(tri, [1, 2, 4, 5, 3])
_tri2 = deepcopy(tri)

remaining_points = [1, 2, 4, 5, 3][(begin+3):end]
initial_search_point = DelaunayTriangulation.get_initial_search_point(tri, 1,
    remaining_points[1], [1, 2, 4, 5, 3], DelaunayTriangulation.default_num_samples,
    Random.default_rng(), true)
DelaunayTriangulation.add_point_bowyer_watson!(tri, remaining_points[1], initial_search_point, Random.default_rng()) # new point is in (1, 4, -1)
_tri3 = deepcopy(tri) 

initial_search_point = DelaunayTriangulation.get_initial_search_point(tri, 2,
    remaining_points[2], [1, 2, 4, 5, 3], DelaunayTriangulation.default_num_samples,
    Random.default_rng(), true)
DelaunayTriangulation.add_point_bowyer_watson!(tri, remaining_points[2], initial_search_point, Random.default_rng()) # new point is in (5, 4, -1)
_tri4 = deepcopy(tri)

Let's now validate at each stage.

julia> DelaunayTriangulation.validate_triangulation(_tri)
true

julia> DelaunayTriangulation.validate_triangulation(_tri2)
true

julia> DelaunayTriangulation.validate_triangulation(_tri3)
The triangle (5, 2, 4) is degenerately oriented.

false

julia> DelaunayTriangulation.validate_triangulation(_tri4)
The adjacent2vertex map is inconsistent with the adjacent map. The vertex 5 and edge (2, 3) are inconsistent with the associated edge set Set([(3, -1), (1, 2), (-1, 1), (2, 3)]) from the adjacent2vertex map, as (2, 3) maps to 4 instead of 5.
The triangle (3, 5, 2) is incorrectly stored in the adjacent map, with the edge (2, 3) instead mapping to 4.
The Delaunay criterion does not hold for the triangle-vertex pair ((3, 4, 2), 1).
The triangle (3, 2, 4) is negatively oriented.
The edge iterator has duplicate edges.
The solid_edge iterator has duplicate edges.

false

julia> DelaunayTriangulation.get_point(_tri3, 5, 2, 4)
((-0.5, 0.5), (1.0, 0.0), (0.25, 0.25))

There's signs of trouble starting at _tri3. Why is it trying to make a triangle from a degenerate arrangement?

We can probably go even further in debugging (believe me, I have tried many times), but that might be enough for someone to have some ideas. It probably would've been better to use an example that breaks without any collinearity issues to avoid that distraction, but alas.

I think a big issue is that somehow the submerged vertices, once in the triangulation (a submerged vertex might not be submerged for a given $\mathcal P$, but it could be for $\mathcal P \cup {p}$), are not currently being correctly identified as being submerged later on, and so they stay in the triangulation and cause issues. I don't know how to fix this.

JuliaGeometry / DelaunayTriangulation.jl