API Reference

abstract type AbstractSurface{M<:Manifold,C<:CRS} end

A surface of a PointSurface.

struct PointSurface{M,C,S,T} <: AbstractSurface{M,C}

This is a typical representation of a surface via points.

Type Parameters

• M<:Manifold - manifold type
• C<:CRS - coordinate reference system
• S - shadow type
• T<:AbstractTopology - topology type for surface-local connectivity

struct SurfaceElement{M,C,N,A}

Representation of a point on a <:PointSurface.

struct PointBoundary{M,C} <: Domain{M,C}

A boundary of points.

Fields

• surfaces: Named surfaces forming the boundary

Type Parameters

• M <: Manifold: The manifold type
• C <: CRS: The coordinate reference system

struct PointVolume{M,C,T,V} <: Domain{M,C}

Interior volume points with optional topology.

Type Parameters

• M<:Manifold - manifold type
• C<:CRS - coordinate reference system
• T<:AbstractTopology - topology type for volume-local connectivity
• V<:AbstractVector{Point{M,C}} - storage type (allows GPU arrays)

struct PointCloud{M,C,T} <: Domain{M,C}

A point cloud with optional topology (connectivity).

Type Parameters

• M<:Manifold - manifold type
• C<:CRS - coordinate reference system
• T<:AbstractTopology - topology type for cloud-level connectivity

points(surf::PointSurface)

Return vector of points from surface. Alias for point(surf).

points(vol::PointVolume)

Return vector of points from volume.

points(boundary::PointBoundary)

Return vector of all points from all surfaces in the boundary.

points(cloud::PointCloud)

Return vector of all points (boundary + volume).

point(surf::PointSurface)

Return the vector of point coordinates for all surface elements.

normal(surf::PointSurface)

Return the vector of outward unit normal vectors for all surface elements.

area(surf::PointSurface)

Return the vector of surface areas for all surface elements.

topology(surf::PointSurface)

Return the topology of the surface.

topology(vol::PointVolume)

Return the topology of the volume.

topology(cloud::PointCloud)

Return the topology of the point cloud.

abstract type AbstractTopology{S}

Abstract base type for point cloud topology (connectivity). Type parameter S is the storage format for neighbor indices.

struct NoTopology <: AbstractTopology{Nothing}

Singleton type representing no topology. Default for PointCloud.

mutable struct KNNTopology{S} <: AbstractTopology{S}

k-nearest neighbors topology.

Fields

• neighbors::S - neighbor indices storage
• k::Int - number of neighbors per point

mutable struct RadiusTopology{S,R} <: AbstractTopology{S}

Radius-based topology where neighbors are all points within a given radius.

Fields

• neighbors::S - neighbor indices storage
• radius::R - search radius (scalar or function of position)

set_topology(surf::PointSurface, ::Type{KNNTopology}, k::Int)

Build and return new surface with k-nearest neighbor topology.

set_topology(surf::PointSurface, ::Type{RadiusTopology}, radius)

Build and return new surface with radius-based topology.

set_topology(vol::PointVolume, ::Type{KNNTopology}, k::Int)

Build and return new volume with k-nearest neighbor topology.

set_topology(vol::PointVolume, ::Type{RadiusTopology}, radius)

Build and return new volume with radius-based topology.

set_topology(cloud::PointCloud, ::Type{KNNTopology}, k::Int)

Build and return new cloud with k-nearest neighbor topology.

set_topology(cloud::PointCloud, ::Type{RadiusTopology}, radius)

Build and return new cloud with radius-based topology.

rebuild_topology!(topo::NoTopology, points)

No-op for NoTopology (nothing to rebuild).

rebuild_topology!(topo::KNNTopology, points)

Rebuild k-nearest neighbor topology in place.

rebuild_topology!(topo::RadiusTopology, points)

Rebuild radius-based topology in place.

rebuild_topology!(surf::PointSurface)

Rebuild topology in place using same parameters. No-op if NoTopology.

rebuild_topology!(vol::PointVolume)

Rebuild topology in place using same parameters. No-op if NoTopology.

rebuild_topology!(cloud::PointCloud)

Rebuild topology in place using same parameters. No-op if NoTopology.

neighbors(t::AbstractTopology)

Return the neighbor storage from a topology.

neighbors(t::AbstractTopology, i::Int)

Return neighbors of point i.

neighbors(surf::PointSurface)

Return all neighbor lists from the surface topology. Throws error if no topology.

neighbors(surf::PointSurface, i::Int)

Return neighbors of point i in surface-local indices. Throws error if no topology.

neighbors(vol::PointVolume)

Return all neighbor lists from the volume topology. Throws error if no topology.

neighbors(vol::PointVolume, i::Int)

Return neighbors of point i in volume-local indices. Throws error if no topology.

neighbors(cloud::PointCloud)

Return all neighbor lists from the topology. Throws error if no topology or invalid.

neighbors(cloud::PointCloud, i::Int)

Return neighbors of point i. Throws error if no topology or invalid.

hastopology(surf::PointSurface)

Check if surface has a topology (not NoTopology).

hastopology(vol::PointVolume)

Check if volume has a topology (not NoTopology).

hastopology(cloud::PointCloud)

Check if point cloud has a topology (not NoTopology).

discretize(bnd::PointBoundary, spacing; alg=auto, max_points=10_000_000)

Generate volume points for the given boundary and return a new PointCloud.

spacing can be either an AbstractSpacing object or a bare Unitful.Length value (which will be wrapped in ConstantSpacing).

Keyword Arguments

• alg: Discretization algorithm (default: SlakKosec() for 3D)
• max_points: Maximum number of volume points to generate

Example

mesh = GeoIO.load("model.stl").geometry
boundary = PointBoundary(mesh)
octree = TriangleOctree(mesh; h_min=0.1)
cloud = discretize(boundary, 3.0m; alg=SlakKosec(octree), max_points=100_000)

discretize(cloud::PointCloud, spacing; alg=auto, max_points=10_000_000)

Generate volume points for an existing cloud and return a new PointCloud with the volume populated.

discretize(bnd::PointBoundary{𝔼{3}}, alg::OctreeRandom; max_points=10_000_000)

Generate volume points using OctreeRandom without requiring a spacing parameter.

OctreeRandom generates uniformly random points and does not use spacing, so this overload removes the need for a dummy spacing value.

Example

mesh = GeoIO.load("bunny.stl").geometry
boundary = PointBoundary(mesh)
cloud = discretize(boundary, OctreeRandom(mesh); max_points=100_000)

SlakKosec <: AbstractNodeGenerationAlgorithm

Slak-Kosec algorithm for volume point generation with optional octree acceleration.

The algorithm generates candidate points on spheres around existing points and accepts them if they are inside the domain and sufficiently far from existing points.

Fields

• n::Int - Number of candidate points per sphere (default: 10)
• octree::Union{Nothing,TriangleOctree} - Optional octree for fast isinside queries

Constructors

SlakKosec()                          # Default: n=10, no octree
SlakKosec(20)                        # Custom n, no octree
SlakKosec(octree::TriangleOctree)    # Use octree acceleration with n=10
SlakKosec(20, octree)                # Custom n with octree acceleration

Performance

• Without octree: Uses Green's function for isinside (~50ms per query)
• With octree: Uses spatial indexing (~0.05ms per query, 1000× faster!)

Usage Examples

Standard Usage (Green's function)

using WhatsThePoint

# Load boundary
boundary = PointBoundary("model.stl")
cloud = PointCloud(boundary)

# Discretize without octree (slow for large domains)
spacing = ConstantSpacing(1.0u"m")
result = discretize(cloud, spacing; alg=SlakKosec(), max_points=10_000)

Octree-Accelerated Usage (Recommended for large domains)

using WhatsThePoint

# Load boundary points
boundary = PointBoundary("model.stl")
cloud = PointCloud(boundary)

# Build octree from STL file (Option 1: simplest)
octree = TriangleOctree("model.stl"; h_min=0.01, classify_leaves=true)

# Or from SimpleMesh (Option 2)
# mesh = GeoIO.load("model.stl").geometry
# octree = TriangleOctree(mesh; h_min=0.01, classify_leaves=true)

# Use octree-accelerated discretization (100-1000× faster!)
spacing = ConstantSpacing(1.0u"m")
alg = SlakKosec(octree)  # Pass octree to algorithm
result = discretize(cloud, spacing; alg=alg, max_points=100_000)

References

Šlak J, Kosec G. "On generation of node distributions for meshless PDE discretizations" (2019)

VanDerSandeFornberg <: AbstractNodeGenerationAlgorithm

3D volume discretization algorithm that projects a 2D grid onto the shadow plane and fills the volume layer by layer using sphere packing heights. Requires ConstantSpacing.

See: Van der Sande, K. & Fornberg, B. (2021). SIAM J. Sci. Comput., 43(1).

FornbergFlyer <: AbstractNodeGenerationAlgorithm

2D volume discretization algorithm using a height-field approach projected onto the x-axis. This is the default and only algorithm for 2D boundaries. Requires ConstantSpacing.

See: Fornberg, B. & Flyer, N. (2015). Comput. Math. Appl., 69(7).

OctreeRandom <: AbstractNodeGenerationAlgorithm

Octree-guided random point generation for volume discretization.

This algorithm uses a TriangleOctree with leaf classification to efficiently generate random points inside the mesh. It's much faster than rejection sampling because:

• Only samples in interior/boundary regions (skips exterior entirely)
• Interior boxes need no filtering (100% acceptance rate)
• Only boundary boxes require isinside() checks

Fields

• octree::TriangleOctree - Octree with leaf classification
• boundary_oversampling::Float64 - Oversampling factor for boundary boxes (default: 2.0)
• verify_interior::Bool - Per-point signed distance verification for interior leaves (default: false)

Constructors

# From mesh (recommended) — builds octree automatically with classified leaves
OctreeRandom(mesh::SimpleMesh)                          # Auto h_min
OctreeRandom(mesh; h_min=0.01, boundary_oversampling=2.0)  # Custom h_min

# From file path — loads mesh then builds octree
OctreeRandom("bunny.stl")
OctreeRandom("bunny.stl"; h_min=0.01)

# From pre-built octree (advanced)
OctreeRandom(octree::TriangleOctree)
OctreeRandom(octree, oversampling; verify_interior=false)

Performance

• Interior leaves: O(1) per point (no filtering needed)
• Boundary leaves: O(k) per point where k ≈ 10-50 triangles
• Much faster than bounding box rejection: Typical acceptance 5-20% vs 90%+ here

Usage Examples

Recommended Usage

using WhatsThePoint

mesh = GeoIO.load("bunny.stl").geometry
boundary = PointBoundary(mesh)
cloud = discretize(boundary, OctreeRandom(mesh); max_points=100_000)

With Spacing (backward compatible)

octree = TriangleOctree(mesh; h_min=0.01, classify_leaves=true)
cloud = discretize(boundary, 1.0u"m"; alg=OctreeRandom(octree), max_points=10_000)

With Custom Oversampling

alg = OctreeRandom(mesh; boundary_oversampling=2.5)
cloud = discretize(boundary, alg; max_points=10_000)

Notes

• The mesh convenience constructor always sets classify_leaves=true
• When h_min is omitted, it is computed as bbox_diagonal / (2 * cbrt(n_triangles))
• The spacing parameter is not used (random uniform distribution)
• For spacing-based discretization, use SlakKosec or VanDerSandeFornberg
• Actual point count may be slightly less than max_points due to boundary filtering
• Oversampling 1.5-3.0 is recommended; higher values waste computation

Algorithm Details

1. Identify regions: Separate octree leaves into interior (2) and boundary (1)
2. Allocate points: Distribute target count proportionally to leaf volumes
3. Generate interior points: Random sampling in interior boxes (no filtering)
4. Generate boundary points: Oversample and filter with isinside()
5. Return: Combined set of validated interior points

When to Use

Good for:

• Quick initial discretization
• Uniform random distributions
• Maximum point count targets
• Testing and prototyping

Not ideal for:

• Spacing-controlled discretization (use SlakKosec instead)
• Adaptive refinement
• Smooth point distributions

ConstantSpacing{L<:Unitful.Length} <: AbstractSpacing

Constant node spacing.

LogLike <: VariableSpacing

Node spacing based on a log-like function of the distance to nearest boundary $x/(x+a)$ where $x$ is the distance to the nearest boundary and $a$ is a parameter to control the growth rate as $a = 1 - (g - 1)$ where $g$ is the conventional growth rate parameter.

Power <: VariableSpacing

Node spacing based on a power of the distance to nearest boundary $x^{g}$ where $x$ is the distance to the nearest boundary and $g$ is the growth_rate.

compute_normals(surf::PointSurface{𝔼{N},C}; k::Int=5) where {N,C<:CRS}

Estimate the normals of a set of points that form a surface. Uses the PCA approach from "Surface Reconstruction from Unorganized Points" - Hoppe (1992).

Requires Euclidean manifold (𝔼{2} or 𝔼{3}). This function assumes flat space geometry.

compute_normals(search_method::KNearestSearch, surf::PointSurface{𝔼{N},C}) where {N,C<:CRS}

Estimate the normals of a set of points that form a surface. Uses the PCA approach from "Surface Reconstruction from Unorganized Points" - Hoppe (1992).

Requires Euclidean manifold (𝔼{2} or 𝔼{3}). This function assumes flat space geometry.

orient_normals!(search_method::KNearestSearch, normals::AbstractVector{<:AbstractVector}, points)

Correct the orientation of normals on a surface as the compute_normals function does not guarantee if the normal is inward or outward facing. Uses the approach from "Surface Reconstruction from Unorganized Points" - Hoppe (1992).

orient_normals!(normals::AbstractVector{<:AbstractVector}, points::AbstractVector{<:Point{𝔼{N}}}; k::Int=5) where {N}

Correct the orientation of normals on a surface as the compute_normals function does not guarantee if the normal is inward or outward facing. Uses the approach from "Surface Reconstruction from Unorganized Points" - Hoppe (1992).

Requires Euclidean manifold (𝔼{2} or 𝔼{3}). This function uses Euclidean dot products for orientation consistency.

update_normals!(surf::PointSurface{𝔼{N},C}; k::Int=5) where {N,C<:CRS}

Update the normals of the boundary of a surf. This is necessary whenever the points change for any reason.

Requires Euclidean manifold (𝔼{2} or 𝔼{3}). This function assumes flat space geometry.

split_surface!(cloud, angle; k=10)
split_surface!(cloud, target_surf, angle; k=10)

Split a surface into sub-surfaces based on normal angle discontinuities. Builds a k-nearest neighbor graph, removes edges where adjacent normals differ by more than angle, and labels each connected component as a separate named surface.

When called on a cloud/boundary with a single surface, that surface is split automatically. When multiple surfaces exist, specify target_surf by name.

combine_surfaces!(boundary::PointBoundary, surfs...)

Merge multiple named surfaces into one. The first name is kept and subsequent surfaces are merged into it. All original surfaces are removed and replaced by the combined surface.

ShadowPoints(Δ, order=1)
ShadowPoints(Δ::Number, order)

Shadow point configuration for generating virtual points offset inward from the boundary. Δ is the offset distance (constant or a function of position). order is the derivative order for Hermite-type boundary condition enforcement.

generate_shadows(points, normals, shadow::ShadowPoints)
generate_shadows(surf::PointSurface, shadow::ShadowPoints)
generate_shadows(cloud::PointCloud, shadow::ShadowPoints)

Generate shadow points offset inward from the boundary along the normal direction by the distance specified in shadow. Returns a vector of Point objects.

isinside(point::SVector{3,T}, octree::TriangleOctree) -> Bool

Fast interior/exterior test using octree spatial index.

Returns true if point is inside the closed surface defined by the mesh.

Performance

• Complexity: O(log M + k) where M = number of triangles, k ≈ 10-50
• Speedup: 100-1000× faster than brute-force O(M) approach

Example

using WhatsThePoint, StaticArrays

octree = TriangleOctree("model.stl"; h_min=0.01, classify_leaves=true)

point = SVector(0.5, 0.5, 0.5)
is_inside = isinside(point, octree)

isinside(points::Vector{SVector{3,T}}, octree::TriangleOctree) -> Vector{Bool}

Batch interior test for multiple points.

Example

test_points = [SVector(randn(3)...) for _ in 1:1000]
results = isinside(test_points, octree)  # Returns Vector{Bool}

isinside(point::Point{𝔼{3}}, octree::TriangleOctree{𝔼{3},C,T}) -> Bool

Bridge from Meshes.jl Point to SVector-based isinside.

Extracts coordinates from Point, converts to SVector{3,T}, and delegates to the SVector method.

isinside(points::AbstractVector{<:Point{𝔼{3}}}, octree::TriangleOctree{𝔼{3}}) -> Vector{Bool}

Batch bridge from Meshes.jl Points to SVector-based isinside.

TriangleOctree{M<:Manifold, C<:CRS, T<:Real}

An octree-based spatial index for efficient triangle mesh queries.

Fields

• tree::SpatialOctree{Int,T}: Underlying octree storing triangle indices
• mesh::SimpleMesh{M,C}: Original Meshes.jl mesh (single source of truth for triangle data)
• leaf_classification::Union{Nothing, Vector{Int8}}: Classification of empty leaves
  • 0: Exterior (outside mesh)
  • 1: Boundary (near surface but empty)
  • 2: Interior (inside mesh)

Performance

For a mesh with M triangles:

• Construction: O(M log M) - distribute triangles to boxes
• Query: O(log M + k) where k ≈ 10-50 triangles per leaf
• Memory: O(N) for octree nodes (mesh data accessed on-the-fly)

Example

# Option 1: From STL file (simplest)
octree = TriangleOctree("model.stl"; h_min=0.01, classify_leaves=true)

# Option 2: From SimpleMesh
mesh = GeoIO.load("model.stl").geometry
octree = TriangleOctree(mesh; h_min=0.01, classify_leaves=true)

# Fast queries (100-1000x speedup over brute force)
point = SVector(0.5, 0.5, 0.5)
is_inside = isinside(point, octree)

Workflow with PointBoundary

Since PointBoundary discards mesh topology, build both from the same mesh:

mesh = GeoIO.load("model.stl").geometry
boundary = PointBoundary(mesh)
octree = TriangleOctree(mesh; h_min=0.01, classify_leaves=true)
cloud = discretize(boundary, spacing; alg=OctreeRandom(octree))

has_consistent_normals(mesh::SimpleMesh) -> Bool

Check if triangle faces are consistently oriented (manifold orientation test).

This function verifies that all shared edges between adjacent triangles are traversed in OPPOSITE directions, which is the geometric requirement for a properly oriented manifold surface. This test is independent of surface curvature.

Algorithm

For each shared edge between two triangles:

• Triangle A has edge v1→v2
• Triangle B (adjacent) must have edge v2→v1 (opposite direction)
• If both traverse in the same direction → faces are incorrectly oriented

Returns

• true if all triangles are correctly oriented (manifold surface)
• false if any triangles have flipped faces (orientation errors)

Performance

O(n) complexity using edge hash map instead of O(n²) pairwise comparison.

Examples

octree = TriangleOctree("model.stl"; h_min=0.01)
if !has_consistent_normals(octree.mesh)
    @warn "Mesh has flipped triangles - will cause incorrect isinside() results"
end

Note

This test uses GEOMETRIC edge orientation, not algebraic normal dot products. It will correctly validate meshes with high curvature (sharp edges, creases).

emptyspace(testpoint, points)

Check if a point occupies empty space within a certain tolerance.

boundingbox(pts::AbstractVector{<:Point})

Compute the axis-aligned bounding box of a collection of points.

centroid(pts::AbstractVector{<:Point})

Compute the centroid (geometric center) of a collection of points.

repel(cloud::PointCloud, spacing; β=0.2, α=auto, k=21, max_iters=1000, tol=1e-6)

Optimize point distribution via node repulsion (Miotti 2023). Returns (new_cloud, convergence_vector) tuple.

The returned cloud has NoTopology since points have moved.

metrics(cloud::PointCloud; k=20)

Print distance statistics (mean, std, max, min) to the k nearest neighbors for all points in the cloud. Useful for assessing point distribution quality before and after repulsion.

import_surface(filepath::String)

Load a surface mesh from a file (STL, OBJ, or any format supported by GeoIO.jl). Returns a tuple of (points, normals, areas, mesh) where points are face centers.

export_cloud(filename::String, cloud::PointCloud)

Export a point cloud to VTK format. The output file contains boundary point coordinates and normal vectors.

save(filename::String, cloud::PointCloud)

Save a point cloud to a file using FileIO.jl serialization.

AbstractOctree{E,T} = AbstractSpatialTree{3,E,T}

Convenience alias for 3D spatial trees (octrees).

AbstractSpatialTree{N,E,T}

Abstract type for N-dimensional spatial trees.

Type Parameters

• N::Int: Spatial dimensionality (2 for quadtree, 3 for octree)
• E: Element type stored in tree (e.g., Int for indices)
• T<:Real: Coordinate numeric type (Float64, Float32, etc.)

Interface Requirements

Trees must implement:

• find_leaf(tree, point) - Locate leaf containing point
• bounding_box(tree) - Get overall bounds
• num_elements(tree) - Total elements stored

Optional:

• find_neighbors(tree, box_idx, direction) - Neighbor queries
• balance!(tree) - Enforce refinement constraints

AndCriterion{T<:Tuple} <: SubdivisionCriterion

Combine multiple criteria - all must be satisfied for subdivision.

Fields

• criteria::T: Tuple of criteria to combine (parametrized for type stability)

MaxElementsCriterion <: SubdivisionCriterion

Subdivide box if number of elements exceeds threshold.

Fields

• max_elements::Int: Maximum elements before subdivision

NearestTriangleState{T<:Real}

Mutable state for nearest-triangle search, replacing 3 separate Ref allocations with a single struct (1 heap allocation instead of 3, enables potential SROA).

SizeCriterion{T<:Real} <: SubdivisionCriterion

Subdivide box if size exceeds threshold.

Fields

• h_min::T: Minimum box size (stop subdividing when reached)

SpatialOctree{E,T<:Real} <: AbstractOctree{E,T}

Concrete octree implementation using integer coordinate system for efficient neighbor finding.

Uses (i,j,k,N) coordinate system where:

• (i,j,k) are integer coordinates at refinement level N
• Box center = origin + (2*[i,j,k] + 1) * (root_size / N) / 2
• Enables O(1) neighbor calculation from coordinates

Type Parameters

• E: Type of elements stored (e.g., Int for triangle IDs)
• T: Coordinate numeric type (e.g., Float64, Float32)

Fields

• parent::Vector{Int}: Parent box index for each node (0 = root has no parent)
• children::Vector{SVector{8,Int}}: 8 child indices per box (0 = no child)
• coords::Vector{SVector{4,Int}}: (i,j,k,N) coordinates per box where N is refinement level
• origin::SVector{3,T}: Spatial origin of root box
• root_size::T: Size of root box (assumed cubic)
• element_lists::Vector{Vector{E}}: Elements in each box
• num_boxes::Ref{Int}: Current number of boxes (mutable counter)

Interface Implementation

Implements AbstractSpatialTree interface:

• find_leaf(tree, point) - O(log n) point location
• bounding_box(tree) - Root box bounds
• num_elements(tree) - Total stored elements
• find_neighbor(tree, box, dir) - 6-directional neighbor finding
• balance_octree!(tree) - 2:1 refinement constraint

Example

using StaticArrays

origin = SVector(0.0, 0.0, 0.0)
octree = SpatialOctree{Int,Float64}(origin, 10.0)

# Subdivide root
subdivide!(octree, 1)

# Find leaf containing point
point = SVector(2.0, 2.0, 2.0)
leaf = find_leaf(octree, point)

SpatialOctree{E,T}(origin::SVector{3,T}, size::T; initial_capacity=1000)

Create empty octree with root node.

Arguments

• origin: Minimum corner of root box
• size: Edge length of root box (assumed cubic)
• initial_capacity: Initial array capacity (will grow as needed)

Example

origin = SVector(0.0, 0.0, 0.0)
octree = SpatialOctree{Int,Float64}(origin, 10.0)

SubdivisionCriterion

Abstract type for subdivision decision logic.

Allows pluggable subdivision criteria via dispatch.

Built-in Criteria

• MaxElementsCriterion(max_elements) - Subdivide if too many elements
• SizeCriterion(h_min) - Subdivide until box small enough
• AndCriterion(criteria...) - All criteria must be satisfied

Example

criterion = AndCriterion((
    MaxElementsCriterion(50),
    SizeCriterion(0.1)
))

filter(f::Function, vol::PointVolume)

Return new PointVolume with only points satisfying predicate f. Topology is stripped since point indices change.

isvalid(t::AbstractTopology)

Check if topology is valid. With immutable design, topology is always valid if it exists.

_allocate_counts_by_volume(volumes, total_count; ensure_one=false)

Allocate integer counts across leaves proportionally to volumes using largest-remainder rounding, preserving exact total count.

If ensure_one=true and total_count >= length(volumes), each leaf gets at least one allocation before distributing the remainder by volume.

_auto_h_min(::Type{T}, mesh::SimpleMesh) where {T}

Compute a default h_min for octree construction based on mesh geometry.

Heuristic: bbox_diagonal / (2 * cbrt(n_triangles)). Scales leaf size to the mesh's characteristic triangle spacing so the octree resolves surface detail without excessive subdivision.

_build_knn_neighbors(points, k::Int) -> Vector{Vector{Int}}

Build k-nearest neighbor adjacency list from points.

_build_radius_neighbors(points, radius) -> Vector{Vector{Int}}

Build radius-based adjacency list from points.

_classify_empty_leaf_conservative(tree, mesh, leaf_idx) -> Int8

Conservative classification for empty leaves.

• 2 (interior): center and all 8 corners classify inside
• 0 (exterior): center and all 8 corners classify outside
• 1 (boundary): mixed results

This avoids falsely labeling a partially-outside leaf as interior.

Implementation uses octree-accelerated signed distance, not ray casting.

_classify_leaves(tree::SpatialOctree{Int,T}, mesh::SimpleMesh) -> Vector{Int8}

Classify octree leaves as exterior (0), boundary (1), or interior (2).

Uses octree-accelerated signed distance for robust classification.

Returns

• classification::Vector{Int8}: Classification for each box
  • -1: Internal node (not a leaf, unclassified)
  • 0: Exterior leaf
  • 1: Boundary leaf
  • 2: Interior leaf

_collect_classified_leaves(octree::TriangleOctree)

Collect interior and boundary leaves with their volumes from a classified octree.

Returns (interior_leaves, interior_volumes, boundary_leaves, boundary_volumes).

_compute_bbox(::Type{T}, mesh::SimpleMesh) where {T} -> (bbox_min, bbox_max)

Compute bounding box from mesh triangle data.

Performance: Single-pass min/max accumulation over all triangles.

_compute_local_signed_distance(
    point::SVector{3,T},
    mesh::SimpleMesh,
    tri_indices::Vector{Int}
) -> T

Compute signed distance from point to nearest triangle in the local set.

This is the key performance optimization: instead of checking all M triangles, we only check the k≈10-50 triangles in the point's octree leaf.

Performance: Fused computation avoids redundant closest_point calculation.

_compute_signed_distance_octree(point, mesh, tree) -> T

Compute signed distance to the closest triangle using octree branch-and-bound.

This avoids global O(M) triangle scans by searching only relevant octree regions.

_create_root_octree(::Type{T}, mesh::SimpleMesh, n_triangles::Int) where {T} -> SpatialOctree{Int,T}

Create root octree covering mesh bounding box with all triangles in root.

_edge_key(v1::SVector{3,T}, v2::SVector{3,T}) where {T} -> Tuple

Create a canonical edge key from two vertices (order-independent). Uses component-wise comparison for consistent ordering.

_extract_vertex(::Type{T}, vert) where {T} -> SVector{3,T}

Extract coordinates from a Meshes.jl vertex to an SVector{3,T}.

_get_triangle_normal(::Type{T}, mesh::SimpleMesh, tri_idx::Int) where {T} -> SVector{3,T}

Extract and normalize triangle normal from mesh. Accesses mesh data on-the-fly.

_get_triangle_vertices(::Type{T}, mesh::SimpleMesh, tri_idx::Int) where {T} -> (SVector{3,T}, SVector{3,T}, SVector{3,T})

Extract triangle vertices from mesh as SVectors. Accesses mesh data on-the-fly.

_normalize_normal(::Type{T}, n_vec) where {T} -> SVector{3,T}

Extract and normalize a Meshes.jl Vec normal to a unit SVector{3,T}.

_subdivide_triangle_octree!(tree, mesh, box_idx, criterion)

Recursively subdivide a box in the triangle octree.

For each box that needs subdivision:

1. Create 8 child boxes
2. For each triangle in parent, find which children it intersects
3. Distribute triangle indices to children
4. Recursively subdivide children if they meet criteria

_triangle_axis_test(
    axis::SVector{3,T},
    v0::SVector{3,T},
    v1::SVector{3,T},
    v2::SVector{3,T},
    half::SVector{3,T}
) where {T<:Real} -> Bool

Internal helper for triangle-box intersection separating axis test.

Tests if the projection intervals of the triangle vertices and box overlap along the given axis.

Arguments

• axis: Separating axis direction
• v0, v1, v2: Triangle vertices in box-centered coordinates
• half: Box half-extents

Returns

true if intervals overlap (potential intersection), false if separated

add_box!(octree::SpatialOctree, i::Int, j::Int, k::Int, N::Int, parent_idx::Int) -> Int

Add new box to octree. Returns box index.

Automatically grows arrays if capacity exceeded.

Arguments

• i, j, k: Integer coordinates at level N
• N: Refinement level (N=1 is root, N=2 is first subdivision, etc.)
• parent_idx: Index of parent box

Returns

Index of newly created box

all_boxes(octree::SpatialOctree) -> Vector{Int}

Return indices of all boxes (leaves and internal nodes).

all_leaves(octree::SpatialOctree) -> Vector{Int}

Return indices of all leaf boxes.

balance_octree!(octree::SpatialOctree, criterion::SubdivisionCriterion)

Enforce 2:1 balance constraint across entire octree.

Iteratively subdivides boxes that violate the 2:1 constraint until all adjacent boxes differ by at most one refinement level.

Arguments

• criterion: Subdivision criterion (only size constraints are enforced)

Algorithm

1. Collect all leaves
2. Check each leaf for balance violations
3. Subdivide violating neighbors (only respecting physical limits like h_min)
4. Repeat until no violations

Note

• Uses can_subdivide (not should_subdivide) to ignore element count criteria
• Balancing is a geometric constraint, not an optimization decision
• Maximum iterations limit prevents infinite loops. If hit, tree may not be fully balanced.

bounding_box(octree::SpatialOctree) -> (SVector{3}, SVector{3})

Get overall bounding box of tree (root box).

box_bounds(octree::SpatialOctree, box_idx::Int) -> (SVector{3}, SVector{3})

Get (mincorner, maxcorner) of box.

Returns

Tuple of (mincorner, maxcorner) as SVector{3,C}

box_center(octree::SpatialOctree, box_idx::Int) -> SVector{3}

Compute spatial center of box using (i,j,k,N) coordinates.

Formula

center = origin + (2*[i,j,k] + 1) * box_size / 2

where box_size = root_size / N

box_size(tree::AbstractSpatialTree, box_idx::Int) -> Real

Get edge length of box.

Required

Trees must implement this function.

box_size(octree::SpatialOctree, box_idx::Int) -> Real

Get edge length of box.

Box size at refinement level N is: root_size / N

can_subdivide(criterion::SubdivisionCriterion, tree, box_idx) -> Bool

Check if box CAN be subdivided based on physical constraints only.

Unlike should_subdivide, this ignores content-based criteria (like element count) and only checks physical limits (like minimum size). Used during balancing where subdivision is required for geometric correctness, not optimization.

Arguments

• criterion: Subdivision criterion (only size constraints are checked)
• tree: Spatial tree
• box_idx: Index of box to check

Returns

true if box can physically be subdivided, false if at minimum size limit.

Example

# For balancing, we only respect size limits
if needs_balancing(leaf)
    if can_subdivide(criterion, tree, leaf)
        subdivide!(tree, leaf)
    end
end

closest_point_on_triangle(
    P::SVector{3,T},
    v1::SVector{3,T},
    v2::SVector{3,T},
    v3::SVector{3,T}
) where {T<:Real} -> SVector{3,T}

Compute the closest point on triangle (v1, v2, v3) to point P.

Uses barycentric coordinate method from Ericson's "Real-Time Collision Detection". The closest point is computed by:

1. Projecting P onto the triangle plane
2. Computing barycentric coordinates
3. Clamping to triangle if outside

Algorithm

The triangle can be parameterized as: T(u,v) = v1 + u(v2-v1) + v(v3-v1) for u,v ≥ 0, u+v ≤ 1

We find the closest point by solving a constrained minimization problem.

Returns

Point on triangle surface closest to P (may be on edge or vertex).

References

Ericson, "Real-Time Collision Detection", Chapter 5.1.5

distance_point_triangle(
    P::SVector{3,T},
    v1::SVector{3,T},
    v2::SVector{3,T},
    v3::SVector{3,T}
) where {T<:Real} -> T

Compute unsigned distance from point P to triangle (v1, v2, v3).

Returns

Unsigned distance (always positive or zero)

distance_point_triangle(
    P::SVector{3,T},
    v1::SVector{3,T},
    v2::SVector{3,T},
    v3::SVector{3,T},
    normal::SVector{3,T}
) where {T<:Real} -> T

Compute signed distance from point P to triangle (v1, v2, v3).

The distance is:

• Positive if P is on the side of the triangle that the normal points to
• Negative if P is on the opposite side
• Zero if P is on the triangle plane

Algorithm

1. Find closest point Q on triangle to P
2. Compute distance ||P - Q||
3. Determine sign based on which side of triangle P is on

Arguments

• P: Query point
• v1, v2, v3: Triangle vertices in counterclockwise order
• normal: Outward-pointing unit normal vector

Returns

Signed distance (positive = outside, negative = inside for closed surface)

Example

using StaticArrays

# Triangle in xy-plane
v1 = SVector(0.0, 0.0, 0.0)
v2 = SVector(1.0, 0.0, 0.0)
v3 = SVector(0.0, 1.0, 0.0)
normal = SVector(0.0, 0.0, 1.0)

# Point above triangle
P = SVector(0.25, 0.25, 1.0)
d = distance_point_triangle(P, v1, v2, v3, normal)
# d ≈ 1.0 (positive, on normal side)

find_boxes_at_coords(octree::SpatialOctree, i_target::Int, j_target::Int, k_target::Int, N_target::Int) -> Vector{Int}

Find box(es) at given (i,j,k,N) coordinates.

• If exact match found at level Ntarget, returns [boxidx]
• If location is covered by coarser box, returns [coarserboxidx]
• If location is subdivided finer, returns all descendants at that location

Returns

Vector of box indices covering the target location

find_leaf(octree::SpatialOctree, point::SVector{3}) -> Int

Find leaf box containing point. Returns box index.

Traverses tree from root to leaf in O(log n) time.

Arguments

• point: Query point in same coordinate system as octree

Returns

Index of leaf box containing point

Throws

AssertionError if point is outside octree bounds

find_neighbor(octree::SpatialOctree, box_idx::Int, direction::Int) -> Vector{Int}

Find neighbor(s) in given direction. Returns vector of neighbor indices.

Handles 2:1 refinement level difference:

• If neighbor exists at same level: returns [neighbor_idx]
• If neighbor is subdivided (finer): returns children on shared face (up to 4)
• If neighbor doesn't exist (boundary): returns empty vector

Arguments

• box_idx: Box to find neighbor of
• direction: Direction code (1-6, see neighbor_direction)

Returns

Vector of neighbor box indices (may be empty if at boundary)

has_children(tree::AbstractSpatialTree, box_idx::Int) -> Bool

Check if box has been subdivided.

Required

Trees must implement this function.

has_children(octree::SpatialOctree, box_idx::Int) -> Bool

Check if box has been subdivided.

is_leaf(tree::AbstractSpatialTree, box_idx::Int) -> Bool

Check if box is a leaf (has no children).

Required

Trees must implement this function.

is_leaf(octree::SpatialOctree, box_idx::Int) -> Bool

Check if box has no children.

needs_balancing(octree::SpatialOctree, box_idx::Int) -> Bool

Check if subdividing this box would violate 2:1 balance with any neighbor.

Returns true if any neighbor has grandchildren (2-level refinement difference).

neighbor_direction(direction::Int) -> (Int, Int, Int)

Convert direction code to (di, dj, dk) offset.

Direction Codes

• 1: -x (left)
• 2: +x (right)
• 3: -y (bottom)
• 4: +y (top)
• 5: -z (front)
• 6: +z (back)

num_elements(octree::SpatialOctree) -> Int

Total number of elements stored in tree.

should_subdivide(criterion::SubdivisionCriterion, tree, box_idx) -> Bool

Determine if box should be subdivided based on criterion.

Arguments

• criterion: Subdivision criterion to evaluate
• tree: Spatial tree
• box_idx: Index of box to check

Returns

true if box should be subdivided, false otherwise.

subdivide!(octree::SpatialOctree, box_idx::Int) -> SVector{8,Int}

Subdivide box into 8 children. Returns child indices [1:8].

Child Ordering (Standard Octree Convention)

1: (0,0,0) - bottom-left-front (x-, y-, z-) 2: (1,0,0) - bottom-right-front (x+, y-, z-) 3: (0,1,0) - top-left-front (x-, y+, z-) 4: (1,1,0) - top-right-front (x+, y+, z-) 5: (0,0,1) - bottom-left-back (x-, y-, z+) 6: (1,0,1) - bottom-right-back (x+, y-, z+) 7: (0,1,1) - top-left-back (x-, y+, z+) 8: (1,1,1) - top-right-back (x+, y+, z+)

Arguments

• box_idx: Index of box to subdivide (must be leaf)

Returns

SVector{8,Int} of child indices in standard order

triangle_box_intersection(
    v1::SVector{3,T},
    v2::SVector{3,T},
    v3::SVector{3,T},
    box_min::SVector{3,T},
    box_max::SVector{3,T}
) where {T<:Real} -> Bool

Test if triangle (v1, v2, v3) intersects axis-aligned box.

Uses the Separating Axis Theorem (SAT) with 13 potential separating axes:

• 3 box face normals (x, y, z axes)
• 1 triangle normal
• 9 edge-edge cross products

If any axis separates the triangle and box, they don't intersect.

Algorithm

1. Translate triangle and box so box is centered at origin
2. Test each potential separating axis
3. Return false if any axis separates, true otherwise

References

• Akenine-Möller, "Fast 3D Triangle-Box Overlap Testing" (2001)
• Ericson, "Real-Time Collision Detection", Chapter 5.2.9

Performance

Optimized with early-out tests. Average case is much faster than worst case.