This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

So that e.g. this doesn't hang and instead throws:

a = 0.0
b = 1.0
c = 0.0
d = 1.0
nx = 23
ny = 28
T, adj, adj2v, DG, pts = triangulate_structured(a, b, c, d, nx, ny)
triangulate_bowyer(pts)

Documentation should clearly describe what -1 means e.g. in sets returned by get_neighbours(tri).

See Chap. 6 of Cheng, Dey and Shewchuk.

using DelaunayTriangulation, StableRNGs
j = 96
rng = StableRNG(j)
points = rand(rng, 2, 500)
tri = triangulate(points; rng)
k = 19
rng = StableRNG(k)
i = 183
delete_point!(tri, i, rng=rng)

The resulting triangulation is not valid:

julia> get_adjacent(tri, 92, 9)
33

julia> get_adjacent(tri, 9, 92)
0

I'm not sure what a simple fix is for this. The problem region is here:

The red dots show the link (the surrounding polygon). I believe the problem is that those three red dots in the first blue triangle, which form a concave indent shown in the isolated second image, are trying to connect to each other somehow. The edge (92, 9) is the outer-most edge of that blue triangle (i.e. the one opposite the indented vertex).

One fix might be that we just need to avoid reusing the same triangulation. I'm not a fan of this, but it might help avoid this issue. In particular,

function _delete_interior_point!(tri, S, rng)
    triangulate_convex!(tri, S; rng, update_ghost_edges=false) # Works even for boundary nodes, due to the replacement we make above
    return nothing
end

should become

function _delete_interior_point!(tri::Triangulation{P,Ts,I,E,Es,BN,BNM,B,BIR}, S, rng) where {P,Ts,I,E,Es,BN,BNM,B,BIR}
	temp_tri = triangulate_convex(get_points(tri), S;
        IntegerType=I,
        EdgeType=E,
        TriangleType=triangle_type(Ts),
        EdgesType=Es,
        TrianglesType=Ts,
        representative_point_list = get_representative_point_list(tri),
        rng,
        add_ghost_triangles=false,
        add_convex_hull=false,
        compute_centers=false,
        delete_empty_features=false)
	for T in each_triangle(temp_tri)
        add_triangle!(tri, T; protect_boundary=true)
    end
	return nothing
end

This also happens on my development branch for v1.0, which has slightly improved the logic in this function. There I've implemented this fix and it works. v1.0 also uses smarter caching, so that this temporary triangulation doesn't lead to as many allocations. I'll think about this a bit more before going with this fix.

I'm encountering something which sounds weird.
I noticed your package yesterday (9/3/23) and use julia> ]add DelaunayTriangulation to install your package.
It turns out that the version installed is 0.7.2, according to the Project.toml inside the package.
At that version, the function triplot works properly, as your documention "Unconstrained Triangulation" suggests "fig, ax, sc = triplot(tri)". But I've been confused about how this function is included in the DelaunayTriangulation.jl, it seems to be inside this macro:

function triplot end # line 213, DelaunayTriangulation.jl

some lines here

@static if !isdefined(Base, :get_extension)
include("../ext/DelaunayTriangulationMakieCoreExt.jl")
end

some lines here

export triplot

However, It seems that get_extension is indeed defined in the Base module (I'm using win11+julia 1.9.3).
Thus it appears that the "include(***)" enclosed in the @static-end block should not be reached.
But this is not true right?, since you really wrote a line (352) that define this:

@doc (@doc _triplot) DelaunayTriangulation.triplot(args...; kwargs...) = _triplot(args...; kwargs...)

So I was confused from then on.
Until today I found some updates is available within my env, thus I use julia> ]up to update my packages.
After doing that I got the new version of your package (v 0.8.7).
But this time I get the error as mentioned in the title of this comment.
I find that you've removed this @static macro which made things involved.
But I clearly found that this triplot example is still in your documentation(v 0.8.7).

So could you help me? Thanks a lot!

This is already done for geometric predicates, and Float32 is causing a ton of problems with other computations. This is the cause of Issue #72, namely triangle areas are so hard to compute for needle triangles. I propose that internal calls to getxy, getx, and gety be replaced with _getxy, _getx, and _gety where

_getxy(p) = (_getx(p), _gety(p))
_getx(p) = Float64(getx(p))
_gety(p) = Float64(gety(p))

This first function _getxy is already implemented as f64_getxy:

"""
    f64_getxy(p)

Returns the coordinates of the point `p` as `Float64`s.
"""
function f64_getxy(p)
    x, y = getxy(p)
    return Float64(x), Float64(y)
end

Maybe instead of having definitions e.g. for AbstractPoint, AbstractTriangle, AbstractEdge, Point, Triangle, and Edge, I should just make all methods generic and instead require a specific set of methods for these types.

I overcome the following problem when installing, while it works during installation of other packages.
Thanks for your advice !

The CI will always fail since it does not have access to the Gmsh installation. I'm not sure how to add a conditional to skip over these tests when being performed by the CI (I still want them in there when using ] test manually), or if there's a way to somehow install it on each CI run (if that's even allowed with licensing?). This will also bump the code coverage to be accurate if I got it to actually run the code.

Would be important to address with sometime once I learn how.

I believe I've been a bit careless in some places with how I use the terms edges / segments and points / vertices. I should clean this up at some point, so that:

• edges are simply pairs (u, v) in the triangulation, but
• segments are the constrained edges.
  In this way, boundary edge might just refer to an edge that happens to be on the boundary, while a boundary segment should refer to a constrained edge appearing on the boundary. Similarly,
• points are the actual coordinates p[i] = (x[i], y[i]),
• while vertices are the references to the point, i.e. i
  Similarly, wherever I use BoundaryIndex it should probably just be GhostVertex or GhostIndex.

We can easily obtain the convex hull from the Bowyer-Watson code, noting that we can very easily traverse the boundary using the ghost triangles. For de Berg's method, we can use the DelaunayGraph (after #15) and then sort the points based on the polar angle.

@DanielVandH thanks for putting this package together. In the below I attempt to create a crossed square mesh (4 corners and a central point). I then wanted to test the removal of points from this set using delete_point!. Perhaps I am using it wrong, but this function does not seem to work properly. In the below example I use delete_point!(TR,5), which seems to work. However, the 5th point is the only point I can remove, any other point, e.g. 1:4 will produce this error: ERROR: BoundsError: attempt to access 5-element Vector{Point3{Float64}} at index [0].

Any help would be appreciated.

using GeometryBasics
using DelaunayTriangulation

# Create a crossed square grid (4 corners and 1 central point)
V = reshape([ Point{3,Float64}(i,j,0.0) for i in -1:2:1, j in -1:2:1],4)
push!(V,Point{3,Float64}(0.0,0.0,0.0))

# Delaunay triangulation
TR = triangulate(V)

# Get faces 
F1 =  [TriangleFace{Int64}(tr) for tr in TR.triangles]
println(F1)

# Remove a selected point
delete_point!(TR,5)

F2 =  [TriangleFace{Int64}(tr) for tr in TR.triangles]
println(F2)

For context, I am removing these red points (which are 3 or 4 connected) in this triangulation:

To obtain the following:

Currently I just repeat your Delaunay method on the set of reduced points. However I was wondering if the delete_point! approach is more efficient.

Thanks.

It would be nice if the bounding triangle were treated symbolically rather than explicitly through some large coordinates. I tried doing this initially, but just couldn't get it to work.

Need to implement constrained Delaunay triangulation at some point. The method here https://www.sciencedirect.com/science/article/pii/S0925772115000322?via%3Dihub would be good. Once this is implemented, de Berg's method will probably be removed.

Some thoughts:

• The interface for this will be important. It is important to be able to support:
  • General vertices.
  • Specified edges. These are just edges that will appear in the triangulation, and could be e.g. a matrix where each column defines an edge. For example, E = [1 2 3; 4 5 6] defines the edges [1, 4], [2, 5], and [3, 6] that have to appear in the triangulation. Could also allow E = [[1, 4], [2, 5], [3, 6]] for the same thing.
  • Specified holes. Holes will be a list of vertices v such that v[begin] == v[end], listed in clockwise order. Holes will be special as all edges that are inside the hole are removed. How this is done is to be decided, but the adjacent map should make this relatively simple. Note that the clockwise order is important to be consistent with all triangles being counterclockwise.
  • A specification for the boundary. Need to somehow handle a boundary loop as being different to a hole, and to allow the boundary to simply be the convex hull. Will need to think about this a bit more. Maybe we simply have a rule that boundaries are given counterclockwise and holes clockwise, so that if no counterclockwise boundary is given then we default to the convex hull. This could be extended to allow for the triangulation of multiply-connected domains, where there are multiple exterior boundaries.

With this interface, we could define three methods:

triangulate(pts) # Unconstrained Delaunay
triangulate(pts, edges) # Constrained Delaunay, no holes (this is a planar straight line graph)
triangulate(pts, edges, loops) # Constrained Delaunay, with holes (this is a piecewise linear complex)

(Maybe a PSLG or PLC struct needs to be defined?)

The returned value for all of these could still be a Triangulation. Triangulation would be updated to be something like

struct Triangulation
    points
    adjacent
    adjacent2vertex
    graph
    edges # better named needed probably - this will just be the list of constrained edges, not all edges in the triangulation - edges(graph) serves that purpose
    loops
    boundary_nodes
end

(The boundary_nodes field is just something extra to store boundaries for both the exterior and for the interiors.)

• With these new features, assigning a value to edges that have no associated triangle is more complicated. e.g. if (i, j) is a boundary edge, currently adj(i, j) = 0 so that (i, j) is a boundary edge and (i, j, 0) is a ghost triangle. When we have interior holes, this is more delicate. Perhaps if we have a set of loops loops = [v1, v2, v3] (say), then adj(i, j) = -1 if the edge is a part of v1, -2 if it is a part of v2, and -3 if it is a part of v3, and so on. (Note that this drops support for non-one-based-indexing.) This is probably the simplest choice, though we would need to adjust DefaultAdjacentValue to something other than -4 (de Berg being removed will free the indices -1, -2, -3.

• Point location will also be more complicated. This is where I struggled the most with the development of the unconstrained code (so many corner cases...), so this will most likely require the most thought again.

• My Gmsh code could be extended as a good way to test all this, providing the interfaces I give above for triangulate except for generate_mesh (except entirely in terms of loops, since the main point there is dealing with boundaries and holes).

Will probably have to sit on these ideas for a while before getting the time to complete it, and to just think about the interface in general. The point location work needs a lot of thought as mentioned.

i  = 265 
rng = StableRNG(i)
points = [(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0), (0.5, 0.5), (0.2, 0.8), (0.1, 0.785)]
tri = triangulate(points; rng)

hits an infinite loop, since

julia> p
(0.5, 0.5)

julia> q
(0.2, 0.8)

julia> r
(0.0, 1.0)

julia> orient(p,q,r)
1

julia> DelaunayTriangulation.triangle_area(p,q,r)
0.0

It would be nice to have a method for sorting points in space, such as with a Hilbert sort. Obviously this could be provided by the user (with the point_order kwarg), but it would be nice to have a function in the code for ease of use. Will need to find some implementation for it that works for general points without any domain restrictions.

The main objective would be to implement a biased randomised insertion order (see Sec 5.3 of the Delaunay book) based on this spatial sorting.

We are advancing geometric processing in Julia with the Meshes.jl umbrella. Would it make sense to port this code over there? We can always integrate code using dependencies, but many things like exact predicates and other low-level functions are already available there, which makes ExactPredicates.jl and other dependencies not necessary.

Currently, the method for constrained triangulation assumes that segments do not intersect, although collinear segments are automatically detected (most of the time, floating point arithmetic can only do so much). It would be nice to have a way to automatically detect this in a user's provided edge input. e.g. some routine like

1. Merge the provided segments and boundary edges into a single set (merge_constrained_edges does this).
2. Compute all intersections between the segments, recording for each intersection a tuple (p, segments), where p is the intersection point and segments are the segments involved in the intersection.
3. For each (p, segments), where segments is say (c1, c2, ...), split each ci in segments at p. A method like split_constrained_edge can do this, similar to how we handle it in process_collinear_segments!. Make sure boundary edges appropriately handled.
4. Proceed with the triangulation.

Maybe add_edge! could have an initial seeking step that looks through all inserted segments (all_constrained_edges) and just finds intersections there, in case users add edges after the initial construction of a triangulation. An option could be toggled to allow this, so that we disable it while the triangulation is being built but when using add_edge! after the fact it is enabled by default.

Currently I define default methods for specific types, e.g. getx(::AbstractVector), but perhaps its smarter to just define getx(p) = p[1] as a default and trust that users can redefine if needed, rather than throwing a MethodError. If this is done, I should provide a function that allows a user to test their interface matches the package's assumptions correctly.

Should also check if any methods really need to exist or can be better defined. Really just need to put a lot more thought into this overall for the longevity of the package.

Found while looking into MakieOrg/Makie.jl#3102.

points = [
    Float32[0.32965052, 0.7966664],
    Float32[0.015137732, 0.31555605],
    Float32[0.54775107, 0.7222514],
    Float32[0.687552, 0.6982844],
    Float32[0.65762305, 0.5177773],
    Float32[0.9649102, 0.8300816],
    Float32[0.12174326, 0.82220316],
    Float32[0.007668495, 0.7747718],
    Float32[0.3144052, 0.5493178],
    Float32[0.6848137, 0.12493092],
    Float32[0.39197737, 0.6912688],
    Float32[0.41400427, 0.025964081],
    Float32[0.35919905, 0.7255166],
    Float32[0.80712754, 0.3415957],
    Float32[0.542679, 0.51094216],
    Float32[0.092720866, 0.90151125],
    Float32[0.90992355, 0.8814645],
    Float32[0.02194357, 0.00064593554],
    Float32[0.9616154, 0.10633117],
    Float32[0.0044495463, 0.97074896],
    Float32[0.4309939, 0.5323847],
    Float32[0.90867966, 0.55974066],
    Float32[0.580766, 0.7668439],
    Float32[0.8563475, 0.88143903],
    Float32[0.18311942, 0.8367877]
]
voronoi(triangulate(points), true)

┌ Warning: Invalid input, empty interval is returned
└ @ IntervalArithmetic C:\Users\User\.julia\packages\IntervalArithmetic\EquAX\src\intervals\intervals.jl:106
┌ Warning: Invalid input, empty interval is returned
└ @ IntervalArithmetic C:\Users\User\.julia\packages\IntervalArithmetic\EquAX\src\intervals\intervals.jl:106
ERROR: DivideError: integer division error
Stacktrace:
  [1] sub!(z::Rational{BigInt}, x::Rational{BigInt}, y::Rational{BigInt})
    @ Base.GMP.MPQ .\gmp.jl:1012
  [2] -
    @ ExactPredicates .\gmp.jl:1061 [inlined]
  [3] orient_reference(u::Tuple{Float64, Float64}, v::Tuple{Float64, Float64}, w::Tuple{Float64, Float64})
    @ ExactPredicates C:\Users\User\.julia\packages\ExactPredicates\GSx4w\src\Codegen.jl:171
  [4] orient_slow(u::Tuple{Float64, Float64}, v::Tuple{Float64, Float64}, w::Tuple{Float64, Float64})
    @ ExactPredicates C:\Users\User\.julia\packages\ExactPredicates\GSx4w\src\Codegen.jl:449
  [5] orient
    @ DelaunayTriangulation C:\Users\User\.julia\packages\ExactPredicates\GSx4w\src\Codegen.jl:459 [inlined]
  [6] orient_predicate(p::Tuple{Float32, Float32}, q::Tuple{Float32, Float32}, r::Tuple{Float32, Float32})
    @ DelaunayTriangulation c:\Users\User\.julia\dev\DelaunayTriangulation\src\predicates\general.jl:18
  [7] point_position_relative_to_line
    @ DelaunayTriangulation c:\Users\User\.julia\dev\DelaunayTriangulation\src\predicates\general.jl:159 [inlined]
  [8] intersection_of_edge_and_bisector_ray(a::Tuple{Float32, Float32}, b::Tuple{Float32, Float32}, c::Tuple{Float32, Float32})
    @ DelaunayTriangulation c:\Users\User\.julia\dev\DelaunayTriangulation\src\geometry_utils\intersections.jl:158
  [9] process_ray_intersection!(vorn::VoronoiTessellation{…}, u::Int64, v::Int64, incident_polygon::Int64, intersected_edge_cache::Vector{…}, segment_intersections::Vector{…}, boundary_sites::Dict{…}, exterior_circumcenters::Set{…}, equal_circumcenter_mapping::Dict{…})
    @ DelaunayTriangulation c:\Users\User\.julia\dev\DelaunayTriangulation\src\voronoi\clipped_construction.jl:77
 [10] process_polygon!(vorn::VoronoiTessellation{…}, e::Tuple{…}, incident_polygon::Int64, boundary_sites::Dict{…}, segment_intersections::Vector{…}, intersected_edge_cache::Vector{…}, exterior_circumcenters::Set{…}, equal_circumcenter_mapping::Dict{…})
    @ DelaunayTriangulation c:\Users\User\.julia\dev\DelaunayTriangulation\src\voronoi\clipped_construction.jl:313
 [11] dequeue_and_process!(vorn::VoronoiTessellation{…}, polygon_edge_queue::Queue{…}, edges_to_process::Set{…}, intersected_edge_cache::Vector{…}, left_edge_intersectors::Set{…}, right_edge_intersectors::Set{…}, current_edge_intersectors::Set{…}, processed_pairs::Set{…}, boundary_sites::Dict{…}, segment_intersections::Vector{…}, exterior_circumcenters::Set{…}, equal_circumcenter_mapping::Dict{…})
    @ DelaunayTriangulation c:\Users\User\.julia\dev\DelaunayTriangulation\src\voronoi\clipped_construction.jl:437
 [12] find_all_intersections(vorn::VoronoiTessellation{Triangulation{…}, Tuple{…}, Int64, Tuple{…}, Set{…}, Tuple{…}})
    @ DelaunayTriangulation c:\Users\User\.julia\dev\DelaunayTriangulation\src\voronoi\clipped_construction.jl:476
 [13] clip_voronoi_tessellation!(vorn::VoronoiTessellation{Triangulation{…}, Tuple{…}, Int64, Tuple{…}, Set{…}, Tuple{…}}, is_convex::Bool)
    @ DelaunayTriangulation c:\Users\User\.julia\dev\DelaunayTriangulation\src\voronoi\main.jl:42
 [14] voronoi(tri::Triangulation{Vector{…}, Set{…}, Int64, Tuple{…}, Set{…}, Vector{…}, Dict{…}, OrderedDict{…}, OrderedDict{…}, Dict{…}}, clip::Bool)
    @ DelaunayTriangulation c:\Users\User\.julia\dev\DelaunayTriangulation\src\voronoi\main.jl:30
 [15] top-level scope
    @ Untitled-1:28
Some type information was truncated. Use `show(err)` to see complete types.

I should get this package to 1.0 when I eventually find some more time. Some goals before doing so:

• Better default methods for interface methods (#62)
• Provide supporting for clipping Voronoi tessellations on more general geometries (#48)
• I want to eventually support curved segments. I don't think this will be in 1.0, but I should at least double check that the interface I currently use could be easily extended to support a more general definition of a boundary (e.g. I provide a function defining a boundary which I can then sub-sample to get boundary notes until I have "enough", as is required for some curved boundary methods).
• Spatial sorting is a big need (#34). I just haven't found a good reference for it.
• Better support should be provided for cleaning up a user's input if needed, or at least making it possible. There is currently check_args that has to be used e.g. for disjoint domains - I wonder if that can be modified somehow (it's currently like that because I don't know of an efficient way for testing if a domain is in disjoint components). This could be something separate to triangulate, since it is nice to assume that everything is working (if a user provides duplicate points, should I remove them? error? continue? etc. This is hard to think about since it may make the user use e.g. each_point which will include the duplicate points but won't match the vertices in the triangulation). A related issue is #42. This can probably be done before or after 1.0, but I should think about it first.
• Of course, it would be nice to resolve #49 and #44, but those can be done anytime. Not a priority right now, though.
• Simplify the documentation. Currently there's too much of a mix between example, tutorial, and background explanation. Docs are being revamped currently in the new-docs branch https://github.com/DanielVandH/DelaunayTriangulation.jl/tree/new-docs.

Meanwhile I solved this issue

This will not be a breaking change, just needs to adjust the Graph data structure to rely on Graphs.jl instead.

Or I can just create the own graphs myself without relying on any Graphs dependency, and then just show an example of how to convert to Graphs.jl?

I have recently faced this compilation error.

Reinstalling both ExactPredicates and DelaunayTriangulation failed to yield any results.

julia> using DelaunayTriangulation
[ Info: Precompiling DelaunayTriangulation [927a84f5-c5f4-47a5-9785-b46e178433df]
ERROR: LoadError: MethodError: no method matching (StaticArraysCore.SVector)(::ExactPredicates.Codegen.Formula, ::ExactPredicates.Codegen.Formula)
Stacktrace:
 [1] var"@genpredicate"(__source__::LineNumberNode, __module__::Module, args::Vararg{Any})
   @ ExactPredicates.Codegen ~/.julia/packages/ExactPredicates/GSx4w/src/Codegen.jl:391
 [2] include(mod::Module, _path::String)
   @ Base ./Base.jl:419
 [3] include(x::String)
   @ ExactPredicates ~/.julia/packages/ExactPredicates/GSx4w/src/ExactPredicates.jl:4
 [4] top-level scope
   @ ~/.julia/packages/ExactPredicates/GSx4w/src/ExactPredicates.jl:54
 [5] include
   @ ./Base.jl:419 [inlined]
 [6] include_package_for_output(pkg::Base.PkgId, input::String, depot_path::Vector{String}, dl_load_path::Vector{String}, load_path::Vector{String}, concrete_deps::Vector{Pair{Base.PkgId, UInt64}}, source::String)
   @ Base ./loading.jl:1554
 [7] top-level scope
   @ stdin:1
in expression starting at /home/pedro/.julia/packages/ExactPredicates/GSx4w/src/plane.jl:6
in expression starting at /home/pedro/.julia/packages/ExactPredicates/GSx4w/src/plane.jl:6
in expression starting at /home/pedro/.julia/packages/ExactPredicates/GSx4w/src/ExactPredicates.jl:4
in expression starting at stdin:1
ERROR: LoadError: Failed to precompile ExactPredicates [429591f6-91af-11e9-00e2-59fbe8cec110] to /home/pedro/.julia/compiled/v1.8/ExactPredicates/jl_cK3mYs.
Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:35
  [2] compilecache(pkg::Base.PkgId, path::String, internal_stderr::IO, internal_stdout::IO, keep_loaded_modules::Bool)
    @ Base ./loading.jl:1707
  [3] compilecache
    @ ./loading.jl:1651 [inlined]
  [4] _require(pkg::Base.PkgId)
    @ Base ./loading.jl:1337
  [5] _require_prelocked(uuidkey::Base.PkgId)
    @ Base ./loading.jl:1200
  [6] macro expansion
    @ ./loading.jl:1180 [inlined]
  [7] macro expansion
    @ ./lock.jl:223 [inlined]
  [8] require(into::Module, mod::Symbol)
    @ Base ./loading.jl:1144
  [9] include
    @ ./Base.jl:419 [inlined]
 [10] include_package_for_output(pkg::Base.PkgId, input::String, depot_path::Vector{String}, dl_load_path::Vector{String}, load_path::Vector{String}, concrete_deps::Vector{Pair{Base.PkgId, UInt64}}, source::Nothing)
    @ Base ./loading.jl:1554
 [11] top-level scope
    @ stdin:1
in expression starting at /home/pedro/.julia/packages/DelaunayTriangulation/PegdX/src/DelaunayTriangulation.jl:1
in expression starting at stdin:1
ERROR: Failed to precompile DelaunayTriangulation [927a84f5-c5f4-47a5-9785-b46e178433df] to /home/pedro/.julia/compiled/v1.8/DelaunayTriangulation/jl_rQNYFo.
Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:35
 [2] compilecache(pkg::Base.PkgId, path::String, internal_stderr::IO, internal_stdout::IO, keep_loaded_modules::Bool)
   @ Base ./loading.jl:1707
 [3] compilecache
   @ ./loading.jl:1651 [inlined]
 [4] _require(pkg::Base.PkgId)
   @ Base ./loading.jl:1337
 [5] _require_prelocked(uuidkey::Base.PkgId)
   @ Base ./loading.jl:1200
 [6] macro expansion
   @ ./loading.jl:1180 [inlined]
 [7] macro expansion
   @ ./lock.jl:223 [inlined]
 [8] require(into::Module, mod::Symbol)
   @ Base ./loading.jl:1144

The tests currently take about an hour to complete. It should be at most 20 minutes I reckon - a lot of time is spent on the same triangulations, using randomised insertion orders to check for issues (it is not impossible that something works 99 times and fails once, I've had it previously). Perhaps we can eliminate this in some places, or somehow use multithreading in some places.

A big issue is the function validate_triangulation in tests/helper_functions.jl, where we check that a constrained triangulation is valid, in particular test_delaunay_criterion:

function test_delaunay_criterion(tri; check_ghost_triangle_delaunay=true)
    for T in each_triangle(tri)
        if DT.is_ghost_triangle(T) && !check_ghost_triangle_delaunay
            continue
        end
        for r in each_solid_vertex(tri)
            cert = DT.point_position_relative_to_circumcircle(tri, T, r)
            if DT.is_inside(cert)
                ace = get_all_constrained_edges(tri)
                if DT.is_empty(ace)
                    flag = !DT.is_inside(cert)
                    if !flag
                        println("Delaunay criterion test failed for the triangle-vertex pair ($T, $r).")
                        return false
                    end
                else # This is extremely slow. Should probably get around to cleaning this up sometime.
                    i, j, k = DT.indices(T)
                    ## We need to see if a subsegment separates T from r. Let us walk from each vertex of T to r.
                    i_history = DT.PointLocationHistory{DT.triangle_type(tri),DT.edge_type(tri),DT.integer_type(tri)}()
                    j_history = DT.PointLocationHistory{DT.triangle_type(tri),DT.edge_type(tri),DT.integer_type(tri)}()
                    k_history = DT.PointLocationHistory{DT.triangle_type(tri),DT.edge_type(tri),DT.integer_type(tri)}()
                    DT.is_boundary_index(i) || jump_and_march(tri, get_point(tri, i); point_indices=nothing, m=nothing, try_points=nothing, k=r, store_history=Val(true), history=i_history)
                    DT.is_boundary_index(j) || jump_and_march(tri, get_point(tri, j); point_indices=nothing, m=nothing, try_points=nothing, k=r, store_history=Val(true), history=j_history)
                    DT.is_boundary_index(k) || jump_and_march(tri, get_point(tri, k); point_indices=nothing, m=nothing, try_points=nothing, k=r, store_history=Val(true), history=k_history)
                    i_walk = DT.is_boundary_index(i) ? Set{DT.triangle_type(tri)}() : i_history.triangles
                    j_walk = DT.is_boundary_index(j) ? Set{DT.triangle_type(tri)}() : j_history.triangles
                    k_walk = DT.is_boundary_index(k) ? Set{DT.triangle_type(tri)}() : k_history.triangles
                    all_edges = Set{NTuple{2,DT.integer_type(tri)}}()
                    E = DT.edge_type(tri)
                    for T in Iterators.flatten((each_triangle(i_walk), each_triangle(j_walk), each_triangle(k_walk)))
                        for e in DT.triangle_edges(T)
                            u, v = DT.edge_indices(e)
                            ee = DT.construct_edge(E, min(u, v), max(u, v))
                            push!(all_edges, (min(u, v), max(u, v)))
                        end
                    end
                    all_constrained_edges = Set{NTuple{2,DT.integer_type(tri)}}()
                    for e in each_edge(ace)
                        u, v = DT.edge_indices(e)
                        ee = DT.construct_edge(E, min(u, v), max(u, v))
                        push!(all_constrained_edges, (min(u, v), max(u, v)))
                    end
                    intersect!(all_edges, all_constrained_edges)
                    flags = [DT.line_segment_intersection_type(tri, initial(e), terminal(e), i, r) for e in each_edge(all_edges) for i in filter(!DT.is_boundary_index, (i, j, k))]
                    flag = !all(DT.is_none, flags) || isempty(flags)
                    if !flag
                        println("Delaunay criterion test failed for the triangle-vertex pair ($T, $r).")
                        return false
                    end
                end
            end
        end
    end
    return true
end

I was rushing at the time that I wrote this. I reckon we need only checking the surrounding_polygon for the points, checking if the line from a point to another intersects the polygon at all, and then testing more carefully inside the polygon. triangle_line_segment_intersection would be good for this also - that's a big reason why I initially wrote that predicate...

any plans for this?

https://github.com/Kcnarf/d3-voronoi-treemap

e: Discard some of the discussion below. The development is in https://github.com/DanielVandH/DelaunayTriangulation.jl/tree/voronoi.

Doesn't always work for some reason. Some tests are commented out while this remains (sometimes) broken.

@DanielVandH I am working on a minimal example (only have a clunky one currently), but I've noticed that when one has added interior points that happen to be exactly on a boundary edges, they are then included/dded to the constrained segments provided. This seems to me to be undesirable behaviour. It would be better to exclude both "in" and "on" points from the triangulation.

The example in the README with cells is slower than with standard points. It shouldn't.

MakieOrg/Makie.jl#3079

Could be good to use some proper examples focused on applications, especially as a way to show how the interfaces should be used in a clearer way. For example:

• Using my other package FiniteVolumeMethod.jl, demonstrate adaptive mesh refinement, using a posteriori error estimates to impose area constraints on each triangular element iteratively.
• Demonstrate a cell simulation model which uses Delaunay triangulations and Voronoi tessellations to model cell connectivities and Voronoi polygons to model geometric properties of a cell (i.e. area).
• Image compression with Lloyd's algorithm.

This would accelerate de Berg's method since we don't have to compute so many properties so many times...