Weighted triangulations don't work

[DanielVandH]

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

DelaunayTriangulation.jl/test/helper_functions.jl

Lines 1461 to 1571 in 41ccf23

	#=
	These weighted Delaunay examples come from the following MATLAB code, which exploits
	the relationship between the convex hull of the lifting map and the Delaunay triangulation.
	1:78 are random triangulations, and 79:155 are triangulations of convex polygons.
	For computing these triangulations, we first define the MATLAB functions:
	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
	Then, for the random triangulations, we use:
	digits(16)
	rng(123)
	ctr = 1;
	for pass = 1:3
	if pass == 1
	range = 4:20;
	elseif pass == 2
	range = 21:25:1500;
	else
	range = 15000;
	end
	for n = range
	points = randn(n, 2);
	weights = randn(n, 1);
	s = 10*rand;
	u = rand;
	if u < 1/250
	weights = 0 * weights;
	elseif (1/250 <= u) && (u < 1/5)
	weights = weights.^2;
	elseif (1/5 <= u) && (u < 0.6)
	weights = s * weights;
	elseif (u <= 0.6) && (u < 0.65)
	sparse_prob = rand(n);
	weights(sparse_prob >= 0.95) = 0;
	end
	save_to_file(points, weights, ctr)
	ctr = ctr + 1;
	end
	end
	The triangulations of convex polygons are computed with the code:
	digits(16)
	rng(123)
	ctr = 79;
	for pass = 1:3
	if pass == 1
	range = 6:20;
	elseif pass == 2
	range = 21:25:1500;
	else
	range = [100, 15000];
	end
	for n = range
	if pass < 3
	n0 = 0;
	while n0 <= 3
	points = randn(n, 2);
	convex_poly = points(convhull(points), :);
	points = convex_poly(1:end-1, :);
	n0 = size(points, 1);
	end
	weights = randn(n0, 1);
	else
	theta = linspace(0, 2*pi, n);
	x = cos(theta);
	y = sin(theta);
	points = [x(:), y(:)];
	points = points(1:end-1, :);
	n0 = size(points, 1);
	weights = randn(n0, 1);
	end
	s = 10*rand;
	u = rand;
	if u < 1/250
	weights = 0 * weights;
	elseif (1/250 <= u) && (u < 1/5)
	weights = weights.^2;
	elseif (1/5 <= u) && (u < 0.6)
	weights = s * weights;
	elseif (u <= 0.6) && (u < 0.65)
	sparse_prob = rand(n0);
	weights(sparse_prob >= 0.95) = 0;
	end
	save_to_file(points, weights, ctr)
	ctr = ctr + 1;
	end
	end
	=#

) 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:

1. Find one triangle whose open circumdisk contains $u$. (Actually, we find a triangle containing $u$ instead of just an open circumdisk, but anyway.)
2. Find all the others by a depth-first search.
3. Delete all the triangles containing $u$ in their open circumdisk.
4. 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

DelaunayTriangulation.jl/test/helper_functions.jl

Lines 1461 to 1663 in 41ccf23

	#=
	These weighted Delaunay examples come from the following MATLAB code, which exploits
	the relationship between the convex hull of the lifting map and the Delaunay triangulation.
	1:78 are random triangulations, and 79:155 are triangulations of convex polygons.
	For computing these triangulations, we first define the MATLAB functions:
	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
	Then, for the random triangulations, we use:
	digits(16)
	rng(123)
	ctr = 1;
	for pass = 1:3
	if pass == 1
	range = 4:20;
	elseif pass == 2
	range = 21:25:1500;
	else
	range = 15000;
	end
	for n = range
	points = randn(n, 2);
	weights = randn(n, 1);
	s = 10*rand;
	u = rand;
	if u < 1/250
	weights = 0 * weights;
	elseif (1/250 <= u) && (u < 1/5)
	weights = weights.^2;
	elseif (1/5 <= u) && (u < 0.6)
	weights = s * weights;
	elseif (u <= 0.6) && (u < 0.65)
	sparse_prob = rand(n);
	weights(sparse_prob >= 0.95) = 0;
	end
	save_to_file(points, weights, ctr)
	ctr = ctr + 1;
	end
	end
	The triangulations of convex polygons are computed with the code:
	digits(16)
	rng(123)
	ctr = 79;
	for pass = 1:3
	if pass == 1
	range = 6:20;
	elseif pass == 2
	range = 21:25:1500;
	else
	range = [100, 15000];
	end
	for n = range
	if pass < 3
	n0 = 0;
	while n0 <= 3
	points = randn(n, 2);
	convex_poly = points(convhull(points), :);
	points = convex_poly(1:end-1, :);
	n0 = size(points, 1);
	end
	weights = randn(n0, 1);
	else
	theta = linspace(0, 2*pi, n);
	x = cos(theta);
	y = sin(theta);
	points = [x(:), y(:)];
	points = points(1:end-1, :);
	n0 = size(points, 1);
	weights = randn(n0, 1);
	end
	s = 10*rand;
	u = rand;
	if u < 1/250
	weights = 0 * weights;
	elseif (1/250 <= u) && (u < 1/5)
	weights = weights.^2;
	elseif (1/5 <= u) && (u < 0.6)
	weights = s * weights;
	elseif (u <= 0.6) && (u < 0.65)
	sparse_prob = rand(n0);
	weights(sparse_prob >= 0.95) = 0;
	end
	save_to_file(points, weights, ctr)
	ctr = ctr + 1;
	end
	end
	=#
	const WGT_DIR = joinpath(dirname(dirname(pathof(DelaunayTriangulation))), "test", "triangulation", "weighted_delaunay_mats")
	const NUM_WEGT = length(readdir(WGT_DIR))
	const NUM_CWEGT = NUM_WEGT - 78 # number of convex polygon weighted examples
	function get_weighted_example(i)
	# get the points, weights, and triangles
	file = "$i.txt"
	mat = readdlm(joinpath(WGT_DIR, file), ',')
	num_pts = Int(mat[1, 1])
	points = mat[2:num_pts, 1:2]
	weights = mat[2:num_pts, 3]
	triangles = Int.(mat[num_pts+1:end, :])
	triangles = Set([Tuple(tri) for tri in eachrow(triangles)])
	_triangles = empty(triangles)
	for T in triangles
	i, j, k = T
	p, q, r = eachrow(points[[i, j, k], :])
	or = DT.triangle_orientation(p, q, r)
	if !DT.is_positively_oriented(or)
	push!(_triangles, (j, i, k))
	else
	push!(_triangles, (i, j, k))
	end
	end
	triangles = _triangles
	points = ElasticMatrix(points')
	# since some of the vertices might be submerged, we can't just get the
	# convex hull of points to find the boundary (in cases of collinear points). instead, let's find
	# all the edges which have only one adjoining triangle, and then sort them.
	d = Dict{NTuple{2,Int},Int}()
	for T in triangles
	i, j, k = T
	d[(i, j)] = k
	d[(j, k)] = i
	d[(k, i)] = j
	end
	boundary_nodes = Set{Int}()
	for (i, j) in keys(d)
	if !haskey(d, (j, i))
	push!(boundary_nodes, i, j)
	end
	end
	boundary_nodes = collect(boundary_nodes)
	DT.sort_convex_polygon!(boundary_nodes, points)
	push!(boundary_nodes, boundary_nodes[1])
	# now get all the submerged vertices
	submerged_vertices = Set{Int}()
	nonsubmerged_vertices = Set{Int}()
	for T in triangles
	push!(nonsubmerged_vertices, T...)
	end
	for i in axes(points, 2)
	if i ∉ nonsubmerged_vertices
	push!(submerged_vertices, i)
	end
	end
	# return
	tri = Triangulation(points, triangles, boundary_nodes; weights)
	unlock_convex_hull!(tri)
	return (tri=tri, submerged_vertices=submerged_vertices, nonsubmerged_vertices=nonsubmerged_vertices, weights=weights)
	end
	get_unstructured_weighted_example(i) = get_weighted_example(i)
	function get_convex_polygon_weighted_example(i)
	tri, submerged_vertices, nonsubmerged_vertices, weights = get_weighted_example(78 + i)
	S = eachindex(weights)
	representative_point_list = DT.get_representative_point_list(tri)
	cx, cy = DT.mean_points(get_points(tri), S)
	representative_point_list[1] = DT.RepresentativeCoordinates(cx, cy, length(S))
	return (tri=tri, submerged_vertices=submerged_vertices, nonsubmerged_vertices=nonsubmerged_vertices, weights=weights, S=S)
	end
	function get_all_distances_to_witness_planes(tri, i)
	distances = Dict{NTuple{3,Int},Float64}()
	for T in each_solid_triangle(tri)
	T = DT.sort_triangle(T)
	δ = DT.get_distance_to_witness_plane(tri, i, T)
	distances[T] = δ
	end
	return distances
	end
	function get_nearest_power_point(tri, i)
	all_dists = Dict{Int,Float64}()
	for j in each_solid_vertex(tri)
	p = DT.get_point(tri, i)
	q = DT.get_point(tri, j)
	all_dists[j] = norm(p .- q)^2 - get_weight(tri, i) - get_weight(tri, j)
	end
	return argmin(all_dists)
	end

. 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)

DelaunayTriangulation.jl/src/algorithms/triangulation/basic_operations/add_point.jl

Lines 111 to 120 in 41ccf23

	"""
	add_point!(tri::Triangulation, x, y, w; kwargs...)
	Adds the point `(x, y)` into `tri` with weight `w`. This function requires that [`add_weight!`](@ref)
	is defined on the weights stored in `tri`. The `kwargs` match those from [`add_point!(tri::Triangulation, ::Any)`](@ref).
	"""
	function add_point!(tri::Triangulation, x, y, w; kwargs...)
	add_weight!(tri, w)
	return add_point!(tri, x, y; kwargs...)
	end

point_position_relative_to_circumcircle

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

DelaunayTriangulation.jl/src/predicates/predicates.jl

Lines 476 to 520 in 41ccf23

	@doc """
	point_position_relative_to_circumcircle(tri::Triangulation, i, j, k, ℓ) -> Certificate
	point_position_relative_to_circumcircle(tri::Triangulation, T, ℓ) -> Certificate
	Tests the position of the vertex `ℓ` of `tri` relative to the circumcircle of the triangle `T = (i, j, k)`. The returned value is a [`Certificate`](@ref), which is one of:
	- `Outside`: `ℓ` is outside of the circumcircle of `T`.
	- `On`: `ℓ` is on the circumcircle of `T`.
	- `Inside`: `ℓ` is inside the circumcircle of `T`.
	!!! note "Ghost triangles"
	The circumcircle of a ghost triangle is defined as the oriented outer halfplane of the solid edge of the triangle. See [`point_position_relative_to_oriented_outer_halfplane`](@ref).
	!!! note "Weighted triangulations"
	If `tri` is a weighted triangulation, then an orientation test is instead applied, testing the orientation of the lifted companions of each point to determine if
	`ℓ` is above or below the witness plane relative to `(i, j, k)`. For ghost triangles, the same rule applies, although if the vertex is on the solid edge of the ghost triangle then,
	in addition to checking [`point_position_relative_to_oriented_outer_halfplane`](@ref), we must check if the new vertex is not submerged by the adjoining solid triangle.
	"""
	point_position_relative_to_circumcircle
	function point_position_relative_to_circumcircle(tri::Triangulation, i, j, k, ℓ)
	u, v, w = sort_triangle(i, j, k)
	a, b, c, p = get_point(tri, u, v, w, ℓ)
	if is_ghost_vertex(w)
	cert = point_position_relative_to_oriented_outer_halfplane(a, b, p)
	if is_on(cert) && is_weighted(tri)
	u′, v′, w′ = replace_ghost_triangle_with_boundary_triangle(tri, (u, v, w))
	sub_cert = point_position_relative_to_witness_plane(tri, u′, v′, w′, ℓ)
	is_above(sub_cert) && return Cert.Outside
	return cert
	else
	return cert
	end
	elseif !is_weighted(tri)
	return point_position_relative_to_circle(a, b, c, p)
	else
	cert = point_position_relative_to_witness_plane(tri, i, j, k, ℓ)
	return is_above(cert) ? Cert.Outside :
	is_below(cert) ? Cert.Inside :
	Cert.On
	end
	end
	point_position_relative_to_circumcircle(tri::Triangulation, T, ℓ) = point_position_relative_to_circumcircle(tri, geti(T), getj(T), getk(T), ℓ)

point_position_relative_to_witness_plane

As it says on the tin.

DelaunayTriangulation.jl/src/predicates/predicates.jl

Lines 521 to 541 in 41ccf23

	"""
	point_position_relative_to_witness_plane(tri::Triangulation, i, j, k, ℓ) -> Certificate
	Given a positively oriented triangle `T = (i, j, k)` of `tri` and a vertex `ℓ` of `tri`, returns the position of `ℓ` relative to the witness plane of `T`. The returned value is a [`Certificate`](@ref), which is one of:
	- `Above`: `ℓ` is above the witness plane of `T`.
	- `On`: `ℓ` is on the witness plane of `T`.
	- `Below`: `ℓ` is below the witness plane of `T`.
	# Extended help
	The witness plane of a triangle is defined as the plane through the triangle `(i⁺, j⁺, k⁺)` where, for example, `pᵢ⁺ = (x, y, x^2 + y^2 - wᵢ)` is the lifted companion of `i`
	and `(x, y)` are the coordinates of the `i`th vertex. Moreover, by the orientation of `ℓ` relative to this witness plane we are referring to `ℓ⁺`'s position, not the plane point `ℓ`.
	"""
	function point_position_relative_to_witness_plane(tri::Triangulation, i, j, k, ℓ)
	p⁺ = get_lifted_point(tri, i)
	q⁺ = get_lifted_point(tri, j)
	r⁺ = get_lifted_point(tri, k)
	a⁺ = get_lifted_point(tri, ℓ)
	cert = orient_predicate(p⁺, q⁺, r⁺, a⁺)
	return convert_certificate(cert, Cert.Above, Cert.On, Cert.Below)
	end

is_submerged

DelaunayTriangulation.jl/src/data_structures/triangulation/methods/weights.jl

Lines 156 to 177 in 41ccf23

	@doc """
	is_submerged(tri::Triangulation, i) -> Bool
	is_submerged(tri::Triangulation, i, V) -> Bool
	Returns `true` if the vertex `i` is submerged in `tri` and `false` otherwise. In the
	second method, `V` is a triangle containing `tri`.
	"""
	is_submerged
	function is_submerged(tri::Triangulation, i)
	# A source that mentions that testing if `i` is submerged only needs to consider the triangle that contains it
	# is given in https://otik.uk.zcu.cz/bitstream/11025/21574/1/Zemek.pdf on p.17.
	# (If the link dies, it is the PhD thesis of Michal Zemek, "Regular Triangulation in 3D and Its Applications".)
	is_ghost_vertex(i) && return false
	q = get_point(tri, i)
	V = find_triangle(tri, q)
	return is_submerged(tri, i, V)
	end
	function is_submerged(tri::Triangulation, i, V)
	is_ghost_vertex(i) && return false
	cert = point_position_relative_to_circumcircle(tri, V, i)
	return is_outside(cert)
	end

get_weighted_nearest_neighbour

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

DelaunayTriangulation.jl/src/data_structures/triangulation/methods/weights.jl

Lines 123 to 154 in 41ccf23

	"""
	get_weighted_nearest_neighbour(tri::Triangulation, i, j = rand(each_solid_vertex(tri))) -> Vertex
	Using a greedy search, finds the closest vertex in `tri` to the vertex `i` (which might not already be in `tri`),
	measuring distance in lifted space (i.e., using the power distance - see [`get_power_distance`](@ref)).
	The search starts from the vertex `j` which should be in `tri`.
	"""
	function get_weighted_nearest_neighbour(tri::Triangulation, i, j=rand(each_solid_vertex(tri)))
	if has_vertex(tri, i)
	return i
	else
	return _get_weighted_nearest_neighbour(tri, i, j)
	end
	end
	function _get_weighted_nearest_neighbour(tri, i, j)
	min_δ² = get_power_distance(tri, i, j)
	min_j = j
	for k in get_neighbours(tri, j)
	if !is_ghost_vertex(k)
	δ² = get_power_distance(tri, i, k)
	if δ² < min_δ²
	min_δ² = δ²
	min_j = k
	end
	end
	end
	if min_j == j
	return j
	else
	return _get_weighted_nearest_neighbour(tri, i, min_j)
	end
	end

get_distance_to_witness_plane

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

DelaunayTriangulation.jl/src/data_structures/triangulation/methods/weights.jl

Lines 90 to 121 in 41ccf23

	""""
	get_distance_to_witness_plane(tri::Triangulation, i, V) -> Float64
	Computes the distance between the lifted companion of the vertex `i` and the witness plane to the triangle `V`. If `V` is a ghost triangle
	and `i` is not on its solid edge, then the distance is `-Inf` if it is below the ghost triangle's witness plane and `Inf` if it is above. If `V` is a ghost triangle and `i`
	is on its solid edge, then the distance returned is the distance associated with the solid triangle adjoining `V`.
	In general, the distance is positive if the lifted vertex is above the witness plane, negative if it is below,
	and zero if it is on the plane.
	See also [`point_position_relative_to_witness_plane`](@ref) and [`get_distance_to_plane`](@ref).
	"""
	function get_distance_to_witness_plane(tri::Triangulation, i, V)
	if !is_ghost_triangle(V)
	u, v, w = triangle_vertices(V)
	p⁺ = get_lifted_point(tri, u)
	q⁺ = get_lifted_point(tri, v)
	r⁺ = get_lifted_point(tri, w)
	s⁺ = get_lifted_point(tri, i)
	return get_distance_to_plane(p⁺, q⁺, r⁺, s⁺)
	else
	cert = point_position_relative_to_circumcircle(tri, V, i)
	if is_inside(cert)
	return -Inf
	elseif is_outside(cert)
	return Inf
	else # is_on(cert)
	V′ = replace_ghost_triangle_with_boundary_triangle(tri, V)
	return get_distance_to_witness_plane(tri, i, V′)
	end
	end
	end

get_power_distance

Must be for power diagrams.

DelaunayTriangulation.jl/src/data_structures/triangulation/methods/weights.jl

Lines 77 to 88 in 41ccf23

	"""
	get_power_distance(tri::Triangulation, i, j) -> Float64
	Returns the power distance between vertices `i` and `j`, defined by
	`||pᵢ - pⱼ||^2 - wᵢ - wⱼ`, where `wᵢ` and `wⱼ` are the respective weights.
	"""
	function get_power_distance(tri::Triangulation, i, j)
	dij² = edge_length_sqr(tri, i, j)
	wᵢ = get_weight(tri, i)
	wⱼ = get_weight(tri, j)
	return dij² - wᵢ - wⱼ
	end

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

DelaunayTriangulation.jl/test/triangulation/weighted.jl

Lines 256 to 263 in 41ccf23

	end
	end
	end
	@testset "Convex polygons" begin
	for fi in 1:NUM_CWEGT
	tri, submerged, nonsubmerged, weights, S = get_convex_polygon_weighted_example(fi)
	ctri = triangulate_convex(get_points(tri), S; weights)

. 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

DelaunayTriangulation.jl/src/algorithms/triangulation/unconstrained_triangulation.jl

Lines 147 to 154 in ee4fa67

	if is_weighted(tri)
	#=
	This part here is why the Bowyer-Watson algorithm is not implemented for weighted triangulations yet. I don't understand why it sometimes
	fails here for submerged points. Need to look into it some more.
	=#
	cert = point_position_relative_to_circumcircle(tri, V, _new_point) # redirects to point_position_relative_to_witness_plane
	is_outside(cert) && return V # If the point is submerged, then we don't add it
	end

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.

[DanielVandH]

This, and the development of power diagrams, is also why #148 has seen no development