Performance bottleneck: 'is_boundary_node' and 'has_ghost_vertices' check in point location

[rkube]

Hi,
I'm using DelaunayTriangulation for numerical simulation data and ran into a performance issue. I've been profiling and interpreted the results with Claude. Below is the generated report.
Thanks for providing this excellent package. I'm happy to help with more profiling data or testing proposed fixes.

Summary

When triangulating ~30,000 points, I'm seeing unexpectedly slow performance (~15 seconds for constrained, ~12 seconds unconstrained). Profiling reveals that a significant portion of time is spent in is_boundary_node and has_ghost_vertices checks during point location, which appear to iterate over Sets rather than using O(1) lookups.

Environment

• Julia version: 1.12.3
• DelaunayTriangulation.jl version: 1.6.6
• OS: macOS

Reproduction

using DelaunayTriangulation
using BenchmarkTools

# Generate test points (or use your own ~30k point dataset)
n = 30000
rz = rand(2, n)  # 2×30000 matrix

# Benchmark unconstrained triangulation
@benchmark DelaunayTriangulation.triangulate($rz)

Profile Analysis

I profiled triangulate() on a 29,841 point dataset. The profile shows:

Unconstrained triangulation breakdown (72,296 total samples):

Function	Samples	% of Total
has_ghost_vertices	31,901	44%
is_boundary_node	20,693	29%
find_triangle (total)	22,328	31%

Call stack showing the bottleneck:

add_point_bowyer_watson!
└── find_triangle
    └── _find_triangle
        ├── is_boundary_node → iterate (over Set) — 20,693 samples
        └── num_solid_vertices
            └── num_ghost_vertices
                └── has_ghost_vertices → any (iterates Set) — 10,511 samples

Constrained triangulation adds:

Function	Samples	% of Total
delete_holes! / distance_to_polygon	5,100	9%

Root Cause Analysis

Looking at the source code:

1. has_ghost_vertices (src/data_structures/triangulation/graph.jl:296):

   has_ghost_vertices(G::Graph) = any(is_ghost_vertex, G.vertices)

   This iterates over all vertices with any() instead of maintaining a boolean flag or count.

2. is_boundary_node (src/data_structures/triangulation/methods/boundaries_and_ghosts.jl:119-121):

   is_boundary_node(tri, u) = u ∈ get_all_boundary_nodes(tri)

   This performs a Set membership test, but get_all_boundary_nodes may be reconstructing the set or the iteration pattern is inefficient.

3. num_ghost_vertices calls has_ghost_vertices which iterates the entire vertex set.

These functions are called for every step of the jump-and-march point location, and for every point insertion (~30,000 times).

Suggested Improvements

1. Cache has_ghost_vertices: Maintain a boolean flag or counter that's updated when ghost vertices are added/removed, rather than iterating on every query.

2. Optimize is_boundary_node: If this is called frequently during point location, consider:

   • Storing boundary nodes in a more cache-friendly structure
   • Avoiding the check entirely during internal Bowyer-Watson operations where boundary status may not be needed

3. Cache num_ghost_vertices / num_solid_vertices: These could be maintained incrementally rather than computed via iteration.

Profile Data

Flat profile (top functions by sample count)

 Count  Function
 31901  has_ghost_vertices
 21158  _find_triangle
 20693  is_boundary_node (via iterate)
 10757  num_ghost_vertices
 10511  num_solid_vertices
  8242  getindex (array access)
  8178  _any (from has_ghost_vertices)
  4234  Dict iterate

Tree profile excerpt

add_point_bowyer_watson!
 └── find_triangle: 32,867 samples
     └── _find_triangle: 22,328 samples
         ├── is_boundary_node: 10,432 + 10,261 samples
         │   └── iterate (Set iteration)
         └── num_solid_vertices: 10,511 samples
             └── num_ghost_vertices
                 └── has_ghost_vertices
                     └── any: 31,901 samples

Impact

For my use case (NIMROD plasma simulation mesh with ~30k points), the triangulation takes ~15 seconds when it should likely be closer to 1-2 seconds based on algorithmic complexity. The O(n) ghost/boundary checks happening O(n√n) times (once per point location step) effectively makes the algorithm O(n²√n) instead of the expected O(n log n).

Additional Context

• I also noticed the comment in delete_holes.jl:

  ## I'm really not convinced that I need to be relying so heavily on
  ## point-in-polygon querying in these functions. There should be something
  ## much smarter and faster that I could do - not pleased with any of this at all.

  This suggests the distance_to_polygon approach in delete_holes! is a known area for improvement.

[DanielVandH]

Thanks for the issue! In addition to the results showing the number of times these functions are called, do you have any results concerning the amount of time spent in each function - and code to reproduce that analysis?

[rkube]

Sure.

The profiling results and analysis reported above can be reproduced with this code

using DelaunayTriangulation
using Profile

# Generate test points (or use your own ~30k point dataset)
n = 30000
rz = rand(2, n)  # 2×30000 matrix

triangulate(rz)
@profile 

triangulate(rz)

This is the profile trace I got.

profile_rand.txt

[DanielVandH]

Those results don't seem to match up with the analysis provided in your first message, right? In the past when I've profiled, these functions have never come up as problematic, so any sort of finer-detail timing results would be useful

[rkube]

Thank you for looking into this. Attached

profile_flat_1e7.txt

is a profile using n=10^7, delay=0.0005.

[DanielVandH]

Thanks. I have some ideas how to improve these functions. I'll see what I can do and, once I've got a working implementation, ping you in a PR to try out your actual test cases that were slow.

[DanielVandH]

What results do you see from your same analysis using the branch in #240?