API Reference

Core Types

The fundamental data structures for representing point clouds.

WhatsThePoint.AbstractSurface Type

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

A surface of a PointSurface.

WhatsThePoint.PointSurface Type

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

This is a typical representation of a surface via points.

Type Parameters

• M<:Manifold - manifold type

• C<:CRS - coordinate reference system

• T<:AbstractTopology - topology type for surface-local connectivity

• G<:StructVector - storage type for surface elements

WhatsThePoint.SurfaceElement Type

struct SurfaceElement{M,C,N,A}

Representation of a point on a <:PointSurface.

WhatsThePoint.PointBoundary Type

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

WhatsThePoint.PointVolume Type

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)

WhatsThePoint.PointCloud Type

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

Accessors

Common accessor functions for point cloud types.

WhatsThePoint.points Function

points(surf::PointSurface)

Return vector of point coordinates for all surface elements.

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

Meshes.normal Function

normal(surf::PointSurface)

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

Meshes.area Function

area(surf::PointSurface)

Return the vector of surface areas for all surface elements.

WhatsThePoint.topology Function

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.

WhatsThePoint.boundary Function

boundary(cloud::PointCloud)

Return the cloud's PointBoundary.

WhatsThePoint.volume Function

volume(cloud::PointCloud)

Return the cloud's PointVolume of interior points.

WhatsThePoint.surfaces Function

surfaces(boundary::PointBoundary)

Return an iterator over the boundary's PointSurfaces (without their names).

WhatsThePoint.namedsurfaces Function

namedsurfaces(boundary::PointBoundary)

Return the ordered dictionary mapping surface names (Symbol) to their PointSurfaces.

Base.names Method

names(boundary::PointBoundary) -> Vector{Symbol}

Return the names of all surfaces in the boundary.

WhatsThePoint.hassurface Function

hassurface(boundary::PointBoundary, name) -> Bool

Return true if the boundary has a surface named name.

Topology

Point connectivity for meshless stencils.

WhatsThePoint.AbstractTopology Type

abstract type AbstractTopology{S}

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

WhatsThePoint.NoTopology Type

struct NoTopology <: AbstractTopology{Nothing}

Singleton type representing no topology. Default for PointCloud.

WhatsThePoint.KNNTopology Type

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

k-nearest neighbors topology.

Fields

• neighbors::S - neighbor indices storage

• k::Int - number of neighbors per point

WhatsThePoint.RadiusTopology Type

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)

WhatsThePoint.set_topology Function

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.

WhatsThePoint.rebuild_topology! Function

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.

WhatsThePoint.neighbors Function

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.

WhatsThePoint.hastopology Function

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

Discretization

Volume point generation algorithms and spacing types.

Meshes.discretize Function

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

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() in 3D; in 2D, an Orthtree built from the boundary loops — pass alg = FornbergFlyer() for the older height-field fill, which requires ConstantSpacing)

• max_points: Maximum number of volume points to generate. For the Orthtree algorithm, defaults to an automatic estimate from the spacing integral (∫ 1/h(x)ᴺ dx, N the spatial dimension) when nothing; other algorithms default to 10_000_000.

Example

mesh = import_mesh("model.stl", u"m")
boundary = PointBoundary(mesh)
cloud = discretize(boundary, 3.0m; alg=Orthtree(mesh))

Note

WhatsThePoint's discretize generates volume fill points from a boundary. This differs from Meshes.jl's discretize which converts continuous geometry into a mesh. No dispatch collision exists — argument types are distinct.

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

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

For the Orthtree algorithm, max_points defaults to an automatic estimate from the spacing integral (∫ 1/h(x)ᴺ dx, N the spatial dimension) when nothing. Other algorithms default to 10_000_000.

WhatsThePoint.AbstractNodeGenerationAlgorithm Type

AbstractNodeGenerationAlgorithm

Abstract supertype for volume discretization algorithms (SlakKosec, VanDerSandeFornberg, FornbergFlyer, Orthtree).

WhatsThePoint.AbstractSpacing Type

AbstractSpacing

Interface for spacing functions that control node density during discretization.

Subtypes must be callable with a single Point or Vec argument and return a Unitful.Length representing the desired node spacing at that location.

(s::MySpacing)(p::Union{Point, Vec}) -> Unitful.Length

See ConstantSpacing, LogLike, and BoundaryLayerSpacing for concrete implementations.

WhatsThePoint.SlakKosec Type

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", u"m")
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 the mesh once; boundary and octree share it
mesh = import_mesh("model.stl", u"m")
boundary = PointBoundary(mesh)
cloud = PointCloud(boundary)
octree = TriangleOctree(mesh; min_ratio=1e-6, 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)

WhatsThePoint.VanDerSandeFornberg Type

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

WhatsThePoint.FornbergFlyer Type

FornbergFlyer <: AbstractNodeGenerationAlgorithm

2D volume discretization algorithm using a height-field approach projected onto the x-axis. Requires ConstantSpacing and offers no Poisson-disk guarantee, so it is no longer the 2D default — Orthtree is. Pass alg = FornbergFlyer() explicitly to use it.

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

WhatsThePoint.Orthtree Type

Orthtree <: AbstractNodeGenerationAlgorithm

Spacing-driven volume discretization algorithm.

Note: This is not solution-adaptive (AMR). Refinement is determined a priori by a prescribed spacing function, not by computed solution features.

Note

Orthtree is a discretization algorithm that generates volume fill points. TriangleOctree is a separate spatial data structure used internally for mesh geometry queries. They serve different purposes.

Algorithm

Uses two trees internally (octrees in 3D, quadtrees in 2D):

• Geometry index: TriangleOctree (3D) / SegmentQuadtree (2D) — captures geometry (surfaces, curvature)

• Node tree: Spacing-driven subdivision where h_box ≤ alpha * h_spacing(x)

Points are generated by weighted allocation across tree leaves based on local spacing requirements.

Key Parameters

• min_ratio: Geometry-index resolution (default: auto from geometry complexity)

• node_min_ratio: Node tree resolution (default: auto from spacing if provided)

• alpha: Subdivision aggressiveness, h_box ≤ alpha * h_spacing (default: 2.0, use 1.0 for fine boundary layers)

• placement: :random, :jittered, :lattice, or :bridson (default: :bridson — global graded Poisson-disk, the recommended production sampler)

• bridson_factor: Poisson-disk radius relative to h(x) for :bridson (default: 0.75)

• boundary_oversampling: Oversampling near boundaries (default: 2.0)

• max_growth: Lipschitz cap on the spacing gradient |∇h| (default: 0.0 = off). When > 0, the prescribed spacing is gradient-limited so neighbouring points differ in spacing by no more than this rate — steep variations stay sharp where the geometry forces them but transition smoothly, which RBF-FD stencils need. 0.1–0.2 matches CFD boundary-layer growth ratios of 1.1–1.2. See "Gradient-limited spacing" below.

Gradient-limited spacing (max_growth)

A raw spacing function may vary faster than a meshless stencil can tolerate: two adjacent points with very different target spacings give an asymmetric, ill-conditioned neighbourhood. With max_growth = g > 0 the algorithm replaces the prescribed field h₀(x) with its g-Lipschitz envelope h(x) = minᵧ (h₀(y) + g·‖x − y‖) — the steepest field that is everywhere ≤ h₀ and grows no faster than g. The limiter runs on the node-tree leaves (a multi-source min-plus relaxation over a k-NN graph of leaf centres) and, because limiting can make the field finer than h₀ in a transition band, it then refines any leaf the new field out-resolves and re-limits, to a fixpoint. The sampler, the grid resolution, and the point-count estimate all then read the limited field, so the delivered point distribution grades smoothly. The limiter is a no-op (no refinement) when h₀ is already g-smooth.

Placement modes

:random, :jittered, and :lattice sample each tree leaf independently (point counts allocated by leaf volume / local spacing). :bridson instead runs a single global advancing-front Poisson-disk pass (Bridson 2007, graded to h(x)) seeded from the boundary points: every generated point keeps a distance of at least min(rᵢ, rⱼ) with r = bridson_factor · h(x) from every other point — including the boundary — by construction. The front saturates at the disk-packing density, so max_points acts as a cap, not a target; a warning is emitted if the cap truncates the front before saturation (which would leave unfilled regions).

bridson_factor (default 0.75) sets the disk radius r = bridson_factor·h relative to the local spacing. A saturated graded front packs ≈ 0.46–0.52 points per r³ (measured on this implementation, k = 30 attempts; geometry-dependent — ≈ 1.09/h³ on a convex box, ≈ 1.22/h³ on the non-convex Stanford bunny). At the 0.75 default that is ≈ 1.1–1.2× the prescribed 1/h³ density, i.e. the front slightly over-fills the nominal budget; the automatic max_points estimate carries matching headroom (see _BRIDSON_CAP_HEADROOM) so the inward-advancing front saturates rather than truncating — a truncated front leaves the deep interior (filled last) empty. Use 1.0 for strict d_NN ≥ h Poisson-disk sampling — ≈ 50% fewer points than 1/h³, which starts repel's spacing-equilibrium force in its attractive branch and degrades rather than polishes the seeding (measured on the cavity gate).

Examples

# Automatic (recommended)
alg = Orthtree(mesh; spacing, alpha=1.0)
cloud = discretize(boundary, spacing; alg, max_points=100_000)

# Manual geometry resolution
alg = Orthtree(mesh; min_ratio=1e-3, spacing, alpha=1.0)

# 2D: closed loop(s) of ordered boundary points (SegmentQuadtree inside)
bnd = PointBoundary(loop_points)
alg = Orthtree(bnd; spacing)
cloud = discretize(bnd, spacing; alg)

WhatsThePoint.ConstantSpacing Type

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

Constant node spacing.

WhatsThePoint.BoundaryLayerSpacing Type

BoundaryLayerSpacing <: VariableSpacing

Smooth spacing transition from fine spacing at the boundary to coarse spacing in the bulk.

Uses physical boundary layer intuition with clear parameters:

• at_wall: Spacing at the boundary surface (fine)

• bulk: Spacing far from boundaries (coarse)

• layer_thickness: Distance over which transition occurs

Example

# Fine 0.5m spacing at walls, coarse 10m in bulk, 8m boundary layer
spacing = BoundaryLayerSpacing(boundary, at_wall=0.5m, bulk=10m, layer_thickness=8m)

Internally uses sigmoid: h(d) = at_wall + (bulk - at_wall) * σ(d) where σ(d) = 1 / (1 + exp(-(d - δ/2) / (δ/6))) and δ = layer_thickness.

WhatsThePoint.LogLike Type

LogLike <: VariableSpacing

Node spacing based on a log-like function of the distance to nearest boundary     where is the distance to the nearest boundary, is base_size, and    is the characteristic length controlling the growth rate, with the conventional growth rate parameter.

Boundary Operations

Normal computation and surface manipulation.

WhatsThePoint.compute_normals Function

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.

WhatsThePoint.orient_normals! Function

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.

WhatsThePoint.update_normals! Function

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.

WhatsThePoint.split_surface! Function

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.

WhatsThePoint.combine_surfaces! Function

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.

Shadow Points

Virtual points offset inward from the boundary for Hermite-type boundary condition enforcement.

WhatsThePoint.ShadowPoints Type

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.

WhatsThePoint.generate_shadows Function

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.

Geometry and Queries

Point-in-volume testing, orthtree acceleration, and spatial utilities.

WhatsThePoint.isinside Function

Fast interior/exterior test using a geometry index.

isinside(testpoint::Point{𝔼{2}}, pts::AbstractVector{<:Point{𝔼{2}}}) -> Bool
isinside(testpoint::Point{𝔼{N}}, cloud::Union{PointCloud, PointBoundary}) -> Bool

Test whether testpoint lies inside the closed domain defined by the boundary points.

For 2D, uses the winding number algorithm — pts must be ordered sequentially around the polygon boundary (clockwise or counter-clockwise). An ArgumentError is thrown if the points do not form a valid ordered polygon.

For 3D, uses a Green's function approach over the boundary surfaces.

Note

WhatsThePoint's isinside tests point-in-polygon/volume membership for meshless point clouds. This is distinct from Meshes.jl's isinside which operates on geometric domain objects. No dispatch collision exists — argument types differ.

WhatsThePoint.TriangleOctree Type

Octree spatial index for triangle mesh queries. Accelerates isinside(), signed distance, etc. Carries the mesh as a TriangleIndex{T} — no SimpleMesh reference, no Meshes.jl in the runtime.

WhatsThePoint.SegmentQuadtree Type

Quadtree spatial index for 2D boundary queries — the 𝔼{2} implementation of AbstractGeometryIndex, mirroring TriangleOctree. Accelerates isinside, signed distance, and leaf classification over the segments of a SegmentIndex.

WhatsThePoint.num_leaves Function

num_leaves(octree::TriangleOctree) -> Int

Return the number of leaf nodes in the octree's spatial subdivision.

num_leaves(quadtree::SegmentQuadtree) -> Int

Return the number of leaf nodes in the quadtree's spatial subdivision.

WhatsThePoint.num_triangles Function

num_triangles(octree::TriangleOctree) -> Int

Return the number of triangles indexed by the octree.

WhatsThePoint.num_segments Function

num_segments(quadtree::SegmentQuadtree) -> Int

Return the number of boundary segments indexed by the quadtree.

WhatsThePoint.has_consistent_normals Function

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

WhatsThePoint.emptyspace Function

emptyspace(testpoint, points)

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

Meshes.boundingbox Function

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

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

Meshes.centroid Function

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

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

Node Repulsion

Point distribution optimization.

WhatsThePoint.repel Function

repel(cloud::PointCloud, spacing; kwargs...) -> PointCloud

Optimize the point distribution via node repulsion. Only volume points move; boundary points form a fixed wall. Points pushed outside the domain are discarded (isinside filter). Returns a new cloud with NoTopology.

Each iteration rebuilds the k-NN graph from the current positions, computes the repel force on every point, and moves it by the adaptive step α_i = clamp(1/|F_i|, α_min, α_max) scaled by the local spacing and capped at one spacing unit. Convergence is the force norm max_i(|F_i|·s_i), which vanishes at equilibrium.

Keywords

• force_model = ClippedSpacingForce(β): force law, any RepelForceModel; the default is repulsive below r = s and zero beyond, so a cloud that already satisfies the Poisson-disk criterion is preserved rather than re-packed. β = 0.2 feeds the default and is ignored when force_model is passed explicitly.

• α, α_min: step-size bounds. Defaults: α = 0.05·min(spacing), α_min = α/100.

• k = 21: neighborhood size.

• max_iters = 1000, tol = 1e-6: iteration and convergence limits.

• cv_target = 0.0: quality-based stop — end the relaxation once the movable points' d_NN/s coefficient of variation drops to this value (read off the sweep's nearest-neighbor data, no extra cost). The natural setting is the raw quality of the direct generation pipeline (≈ 0.07 on the cavity): relaxing past the quality a re-seed would give is wasted budget. The stop returns the configuration the measurement describes (the pre-sweep snapshot), so a cloud already at target comes back unchanged. Off when 0.

• stall_after = 50: stop when that same CV has not improved by ≥0.1 % for this many consecutive iterations. The force residual of a saturated repulsion-only packing plateaus at a nonzero value instead of reaching tol, so cv_target/stall_after are the practical stops for the default force; CV keeps creeping down for hundreds of iterations, making stall_after the backstop (on by default so default runs terminate instead of burning max_iters) and cv_target the primary. Pass 0 to disable and rely on tol/max_iters alone.

• rebuild_every = 1: iterations between k-NN graph rebuilds (larger = cheaper, staler).

• kick_after = 0: if the closest pair freezes at the same r/s for this many iterations (a balanced standoff), kick one point by 0.1·s in a random direction to break the symmetry. Off when 0; 10–20 is reasonable.

• cull_ratio = 0.0: after relaxation, drop near-duplicates closer than cull_ratio·spacing to a kept point. A safety net — a healthy relaxation leaves nothing to cull, so a @warn is emitted whenever it fires.

• convergence: pass an empty float vector (e.g. Float64[]) to collect the per-iteration force norm; entries are computed in the cloud's machine type and convert on insertion.

• trace: pass a NamedTuple[] to record the closest pair each iteration (global boundary-then-volume indices, measured on that iteration's snapshot).

repel(cloud::PointCloud, spacing, octree::TriangleOctree; kwargs...) -> PointCloud

Node repulsion with boundary projection: all points move, boundary points are re-projected onto the mesh surface every iteration, and volume points that escape the domain bounce back — or stick to the surface (see deposit_ratio). Step size, convergence, and the shared keywords are as in the method without octree.

The returned boundary is a single surface named :boundary (use split_surface! to re-establish surface distinctions); topology is NoTopology.

Additional keywords

• deposit_ratio = 0.0: when > 0, an escaped volume point is deposited — projected onto the nearest triangle and converted into a boundary point, accepted only if no boundary point already lies within deposit_ratio·spacing of the landing site. Surface sampling then emerges from volume containment instead of the mesh tessellation; the acceptance test keeps the deposited density self-limiting (conversion is one-way). Deposited points carry the landing triangle's normal and spacing² as area. 0.5–0.7 is reasonable.

• cull_ratio = 0.0: as in the volume-only method; here it also targets boundary pairs that the deterministic projection parks on a shared edge/vertex.

WhatsThePoint.RepelForceModel Type

RepelForceModel

Abstract type for node-repulsion force laws used by repel.

A concrete subtype M <: RepelForceModel must implement compute_force(m::M, u::Real).

WhatsThePoint.ClippedSpacingForce Type

ClippedSpacingForce(β=0.2, u0=1.0)

Repulsion-only force law: F(u) = (u0² − u²) / (u² + β)² for u < u0, zero beyond. The compact support makes any configuration whose pairwise distances all exceed u0·s an exact equilibrium — the Poisson-disk property — so an already-blue-noise cloud is preserved rather than pulled toward a different packing.

This is SpacingEquilibriumForce with the attractive branch removed (identical for u < u0 when u0 = 1). The attractive branch acts as a cohesion force whose preferred bond length s is unreachable at the prescribed density 1/s³: the cloud condenses into locally denser clusters plus voids (measured on the cavity gate: spacing CV and coordination rise, separation falls, starting from an already-good cloud, at a rate proportional to the step size — an instability, not noise). Clipping the force at its root removes the mechanism; the same 300 iterations then improve a constructed blue-noise cloud (CV 0.072 → 0.044) instead of degrading it.

u0 sets the support radius (and root) in units of the local spacing s.

The force residual of a saturated repulsion-only packing plateaus at a small nonzero value instead of vanishing (a frustrated glass), so combine with the stall_after stopping criterion of repel rather than relying on tol alone.

WhatsThePoint.InverseDistanceForce Type

InverseDistanceForce(β=0.2)

Original Miotti (2023) force law F(u) = 1 / (u² + β)². Purely repulsive and monotonically decreasing. β > 0 softens the force near u = 0 and keeps it finite. Has no root, so equilibrium is reached only through damping (α).

WhatsThePoint.SpacingEquilibriumForce Type

SpacingEquilibriumForce(β=0.2)

Force law F(u) = (1 − u²) / (u² + β)² with a zero at u = 1. Repulsive for u < 1 (points closer than the target spacing push apart), attractive for u > 1 (points farther than the target pull together), and zero at the target spacing itself. β > 0 softens the amplitude near u = 0.

Shares the β-softened denominator with InverseDistanceForce; the numerator (1 − u²) introduces the equilibrium at u = 1 without changing the small-u behavior.

WhatsThePoint.StrongSpacingForce Type

StrongSpacingForce(β=0.2, γ=3)

Force law F(u) = (1 − u²) / (u² + β)^γ with a zero at u = 1 and a configurable singularity strength γ. Like SpacingEquilibriumForce but with a stronger repulsive core: at small u the force scales as u^(-2γ) instead of u^(-4). This breaks balanced standoffs where neighbor forces cancel the weaker default core.

γ = 2 recovers SpacingEquilibriumForce. γ = 3 (default) is strong enough to break typical standoffs in a few iterations. Higher values increase the repulsive kick at close range but may require more iterations to settle; the displacement cap in repel prevents runaway.

WhatsThePoint.compute_force Function

compute_force(model::RepelForceModel, u::Real) -> Real

Evaluate the force magnitude at normalized separation u = r / s, where r is the distance between two points and s is the local target spacing.

The returned scalar multiplies the unit vector (xᵢ − xⱼ) / r in the repulsion step, so positive values push xᵢ away from xⱼ and negative values pull it toward xⱼ. Concrete subtypes of RepelForceModel must implement a method for this function.

Diagnostics

WhatsThePoint.metrics Function

metrics(cloud::PointCloud; k=20)

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

Also reports the global separation and fill distances and their ratio, the quasi-uniformity quality measure most relevant to meshless stencil conditioning:

• separation — the smallest nearest-neighbor distance anywhere in the cloud. Small values signal near-coincident points (the source of singular RBF-FD stencils).

• fill — the largest nearest-neighbor distance (a proxy for the worst covering gap).

• mesh_ratio — fill / separation (≥ 1). Closer to 1 means a more uniform, blue-noise-like cloud; large values indicate clustering and voids coexisting.

Returns a NamedTuple with fields avg, std, max, min, separation, fill, mesh_ratio, and k.

WhatsThePoint.spacing_metrics Function

spacing_metrics(cloud::PointCloud, spacing::AbstractSpacing; k=20)

Measure how closely the point distribution matches the target spacing function.

For each point xᵢ, the local actual spacing is estimated as the mean distance to its k nearest neighbors (self excluded). The per-point relative error is

errorᵢ = |r̄ᵢ − s(xᵢ)| / s(xᵢ)

Returns a NamedTuple (max_error, mean_error, std_error, k). Use before and after repel (or any placement step) to quantify spacing preservation.

WhatsThePoint.spacing_fidelity_metrics Function

spacing_fidelity_metrics(cloud::PointCloud, spacing::AbstractSpacing; k=30, coord_radius=1.4)

Per-point spacing fidelity: how well each point's nearest-neighbor distance matches the prescribed spacing h(x).

Computes d_NN(i) / h(x_i) for every point and returns:

• mean_dnn_h — mean of the distribution (ideal ≈ 0.74 for 3D blue-noise)

• cv — coefficient of variation std / mean (lower = more uniform)

• p05, p50, p95 — percentiles of d_NN/h (tight spread = good)

• coordination — mean count of neighbors within coord_radius · h (ideal ≈ 12–14 for 3D blue-noise packing)

• k, coord_radius

Spacing Guidance

WhatsThePoint.suggest_spacing Function

suggest_spacing(mesh; n_points=nothing, bridson_factor=0.75, verbose=true)
suggest_spacing(boundary; ...)
suggest_spacing("model.stl", u"mm"; ...)

Quick geometry probe that recommends a baseline node spacing — the "step 0" before discretize. Reports the domain extent, enclosed volume, and three spacing landmarks, and (with verbose=true) prints a short summary.

The recommendation is driven by the shortest bounding-box axis L_min, because that axis sets how coarse a Poisson-disk fill can be before its interior collapses. With the bridson disk radius r = bridson_factor·h, the interior along L_min has width L_min − 2r, so it is empty once h ≥ L_min/(2·bridson_factor) — the reported h_ceiling. Returned landmarks:

• h_ceiling — coarsest spacing that still hosts any interior; spacings at or above it yield an empty bridson cloud. Stay well below this.

• h_baseline — recommended starting point: ≈10 points across the shortest axis (or, when n_points is given, (volume / n_points)^(1/N) capped to stay fillable). Good enough to run a first simulation, then refine where needed.

• h_fine — h_baseline/2, a second rung for resolving features.

n_baseline/n_fine are rough volume-point counts (volume / hᴺ). In 2D the reported volume is an area and N = 2 throughout.

Returns a NamedTuple with extent, min_extent, max_extent, diagonal, volume, n_triangles, bridson_factor, h_ceiling, h_baseline, h_fine, n_baseline, and n_fine (all spacings/lengths carry units).

Example

mesh = import_mesh("bunny.stl", u"m")
g = suggest_spacing(mesh)
cloud = discretize(PointBoundary(mesh), g.h_baseline; alg=Orthtree(mesh))

Geometry Inspection

WhatsThePoint.geometry_info Function

geometry_info(filepath, filepaths...; verbose=true) -> Vector{NamedTuple}

Inspect the raw bounding box of one or more mesh files before constructing anything — the "what do these numbers mean?" probe for files without unit metadata (STL coordinates are just numbers; GeoIO.jl assigns meters by default). Returns one (file, min, max, extent) named tuple per file, with values as raw (unitless) coordinate tuples. With verbose=true prints each bounding box and, for multiple files, their union — useful when parts must fit together.

Once you know the unit, pass it to import_mesh or PointBoundary — the raw numbers are reinterpreted in that unit.

Example

geometry_info("intake.stl", "exhaust.stl")
# ─── intake.stl ───
#   min:    (0.0, 0.0, 0.0)
#   max:    (120.5, 87.3, 42.0)
#   extent: (120.5, 87.3, 42.0)
# ─── exhaust.stl ───
#   ...
# ─── Union ───
#   ...

Surface Sampling

WhatsThePoint.sample_surface Function

sample_surface(mesh::SimpleMesh, spacing; factor=0.75, max_points=10_000_000,
               stall_limit=2000) -> PointSurface

Sample the surface of mesh with a graded Poisson-disk distribution: blue-noise points with pairwise separation at least min(rᵢ, rⱼ) where r = factor·spacing(x), by construction. An alternative to the face-center sampling of PointBoundary(mesh) whose point density follows the prescribed spacing instead of the mesh tessellation.

Each sample carries its parent triangle's normal; point areas preserve the total mesh surface area, distributed proportionally to r² (equal shares for constant spacing).

Dart throwing runs until saturation — stall_limit consecutive rejections — or until max_points is reached (with a warning, since a truncated pass leaves the surface under-sampled).

Visualization

WhatsThePoint.visualize Function

visualize(x; kwargs...)

Plot a PointCloud, PointBoundary, or PointSurface with Makie. Requires a Makie backend to be loaded (e.g. using GLMakie); the implementation lives in the WhatsThePointMakieExt package extension.

I/O

WhatsThePoint.import_mesh Function

import_mesh(filepath, unit::Unitful.Units) -> SimpleMesh

Load a surface mesh from a file (STL, OBJ, or any format supported by GeoIO.jl) and reinterpret its raw coordinates in unit: a stored 46 becomes 46 mm with unit = u"mm" — no conversion happens, because mesh files carry no unit metadata (GeoIO's default of meters is discarded). Topology and coordinate machine type (Float32 stays Float32) are preserved.

This is the single gateway for file geometry: the returned mesh feeds PointBoundary, TriangleOctree, and Orthtree. Use geometry_info to probe the raw extents when unsure of the unit.

WhatsThePoint.import_surface Function

import_surface(filepath, unit::Unitful.Units)

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. Coordinates are reinterpreted in unit (see import_mesh); areas therefore carry unit^2.

FileIO.save Function

save(filename::String, cloud::PointCloud; format=:jld2)

Save a point cloud to a file.

• format=:jld2 (default): Serialize via FileIO.jl.

• format=:vtk: Export to VTK format with boundary and volume points, normals, and a point type indicator (1 = boundary, 2 = volume).

save(filename::String, boundary::PointBoundary; format=:jld2)

Save a boundary to a file.

• format=:jld2 (default): Serialize via FileIO.jl.

• format=:vtk: Export to VTK format with boundary points and normals.

save(filename::String, surf::PointSurface; format=:jld2)

Save a surface to a file.

• format=:jld2 (default): Serialize via FileIO.jl.

• format=:vtk: Export to VTK format with surface points, normals, and areas.

WhatsThePoint.export_vtk Function

export_vtk(filename, cloud::PointCloud; fields=(), verbose=false)

Write cloud to a ParaView-ready .vtu (one VTK_VERTEX cell per point). Open it in ParaView and set Representation to Point Gaussian (or Points).

Always-attached point data:

• point_type — 1 = boundary, 2 = volume (colour to separate wall from bulk).

• surface_id — 1..N in names(cloud) order, 0 for volume (colour by named surface). Pass verbose=true to print the integer→name legend.

• normals — boundary normals (zero on volume points).

fields attaches solution data, so a .vtu can be re-exported after solving and viewed like any CAE result. It is an iterable of name => values pairs (e.g. a Dict or a tuple of pairs); values may be scalars or per-point vectors and must be ordered like points(cloud) — boundary points first, then volume (the natural global DOF order). Units are stripped automatically.

Examples

export_vtk("cloud", cloud)                                 # geometry only
export_vtk("sol", cloud; fields = ("T" => temp, "U" => velocity))

Unexported API

WhatsThePoint.AbstractGeometryIndex Type

AbstractGeometryIndex{M<:Manifold}

The query interface that node-generation code (build_node_octree, repel) uses to consult boundary geometry, parameterized by the manifold M. Because every seam method dispatches on M, discretization is generic over Manifold. Implementations: TriangleOctree{T} <: AbstractGeometryIndex{𝔼{3}} (3D) and SegmentQuadtree{T} <: AbstractGeometryIndex{𝔼{2}} (2D).

Seam contract — what an implementation provides (the accessors have field-backed defaults covering the standard tree/leaf_classification layout):

• geometry_tree(g) → the index's spatial subdivision tree

• leaf_classes(g) → per-leaf classification cache, or nothing

• _signed_distance(g, p) → signed distance to the boundary (negative inside)

• domain_bounds(g) → (min, max)::NTuple{2,SVector{N,T}}

• project_to_boundary(g, p, offset) → (SVector{N,T}, element_id::Int)

Provided generically on top of that contract (do not reimplement per index — tolerance and fast-path fixes must not fork between dimensions):

• classify_point(g, p, tol) → LEAF_INTERIOR / LEAF_BOUNDARY / LEAF_EXTERIOR

• isinside(p, g) → Bool

WhatsThePoint.NearestElementState Type

NearestElementState{N,T}

Mutable traversal state for _nearest_element_tree!: the best squared distance found so far, the winning element id, the closest point on it, and the feature code of that closest point (element interior / edge / vertex — the codes are geometry-specific, see geometric_utils.jl).

WhatsThePoint.SegmentIndex Type

SegmentIndex{T}

The package's runtime representation of a 2D boundary: one or more closed polyline loops, indexed, unit-stripped and machine-typed — the 2D counterpart of TriangleIndex.

Loops are oriented at construction so that segment normals point out of the domain: outer loops counter-clockwise, hole loops clockwise (nesting parity decides which is which — a loop contained in an odd number of other loops is a hole). Callers may pass loops in any orientation.

Fields:

• vertices/segments: indexed representation (unique coords + per-segment vertex indices; segment i runs vertices[segments[i][1]] → vertices[segments[i][2]]).

• normal: precomputed unit outward segment normals.

• vertex_normal: per-vertex pseudonormal (sum of the two adjacent segment normals) — the sign-exact feature normal for signed-distance queries.

• bbox_min/bbox_max: boundary bounding box.

• len_unit: the unit stripped from coordinates, re-attached at exits.

WhatsThePoint.SegmentIndex Method

SegmentIndex(T, loops, len_unit)

Build a SegmentIndex from closed loops of stripped 2D vertices. Duplicated consecutive vertices are dropped (so the explicitly-closed convention [p₁, …, pₙ, p₁] is accepted), each loop is validated (≥ 3 distinct vertices, non-degenerate signed area) and loop orientation is normalized by nesting parity so normals point out of the domain.

WhatsThePoint.SpacingCriterion Type

SpacingCriterion{T<:Real, S} <: SubdivisionCriterion

Tree subdivision criterion based on local spacing requirements.

Subdivides boxes where h_box > alpha * h_spacing(center), ensuring the tree resolution is fine enough to properly represent the spacing function.

Fields

• spacing::S: Spacing function object

• alpha::T: Subdivision aggressiveness factor

• absolute_min::T: Absolute minimum box size (prevents infinite subdivision)

Algorithm

For each box: 2. Query h_local = spacing(box_center)

3. If h_box > alpha * h_local, subdivide

4. Stop if h_box ≤ absolute_min

Smaller alpha values create finer trees (more aggressive subdivision).

WhatsThePoint.SpatialOctree Type

3D octree — the historical name and the default 3D specialization.

WhatsThePoint.SpatialTree Type

SpatialTree{N,E,T}

Adaptive 2^N-tree over N-d Cartesian space (N=3 octree, N=2 quadtree). Uses a (cell, level) integer coordinate system:

• cell are integer coordinates at refinement level

• box center = origin + (2*cell + 1) * (root_size / level) / 2

• O(1) neighbour calculation from coordinates

Children of a subdivided box are the contiguous block first_child .+ (0:2^N-1) (first_child == 0 marks a leaf).

INVARIANT: subdivide! always allocates the 2^N children contiguously, and children are never re-mutated (there is no coarsening). A future merge/coarsen feature MUST preserve this or replace first_child with an explicit child list.

Type Parameters

• N: spatial dimension

• E: element type stored per box (e.g. Int for triangle IDs)

• T: coordinate numeric type (Float64, Float32)

WhatsThePoint.SpatialTree Method

SpatialTree{N,E,T}(origin::SVector{N,T}, size::T; initial_capacity=1000)

Create an empty tree with a single root box (cubic, edge length size).

WhatsThePoint.SubdivisionCriterion Type

SubdivisionCriterion

Abstract type for tree subdivision decision logic. Subtypes implement should_subdivide(criterion, tree, box_idx) (subdivide if content/size demands it) and can_subdivide(criterion, tree, box_idx) (physical limits only — used during balancing).

WhatsThePoint.TriangleIndex Type

TriangleIndex{T}

The package's runtime representation of a triangle mesh: an indexed, unit-stripped, machine-typed cache built once from a SimpleMesh at the package boundary. After construction, query paths read contiguous SVector{3,T} arrays — they never touch the Meshes.jl object model.

Fields (all canonical — no derived cache):

• vertices/triangles: indexed representation (unique coords + per-triangle vertex indices). Also the serialization format.

• face: precomputed unit face normals.

• edge/vertex: angle-weighted pseudonormals (Bærentzen & Aanæs 2005) keyed by exact coordinates — the sign-exact feature normals for signed-distance queries.

• bbox_min/bbox_max: mesh bounding box (computed once at construction).

• len_unit: the unit stripped from mesh coordinates, re-attached at exits.

This struct is self-sufficient: nothing in the package carries a SimpleMesh reference after TriangleIndex(T, mesh) runs.

WhatsThePoint.VertexResolutionCriterion Type

Subdivide based on vertex density within box bounds.

WhatsThePoint._ClampedSpacing Type

_ClampedSpacing(inner, hmax) <: AbstractSpacing

Spacing that caps inner at hmax: min(inner(p), hmax). Used by _guard_coarse_spacing to make an everywhere-too-coarse spacing fillable without discarding the user's variation where it is already fine.

WhatsThePoint._LeafSpacing Type

_LeafSpacing(node_tree, field, fallback, len_unit) <: AbstractSpacing

Spacing that reads a precomputed per-leaf field (the gradient-limited envelope) by find_leaf lookup, falling back to the original spacing outside the tree. Lets the sampler, grid sizing, and point-count estimate all see the limited field through the usual _spacing_value / call interface.

Base.filter Method

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

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

Base.isvalid Method

isvalid(t::AbstractTopology)

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

WhatsThePoint._allocate_counts_by_volume Method

_allocate_counts_by_volume(volumes, total_count; ensure_one=false)

Proportionally allocate counts to volumes using probabilistic rounding. Fractional parts become selection probabilities, reducing clustering.

WhatsThePoint._apply_gradient_limit Method

_apply_gradient_limit(node_tree, classification, spacing, alg, geometry)
    -> (node_tree, classification, spacing_used)

Gradient-limit the spacing on the node-tree leaves with alg.max_growth, and refine any leaf the limited field now out-resolves (box_size > alpha·h), re-limiting to a fixpoint. Returns the (possibly further-subdivided) tree, the refreshed classification, and a _LeafSpacing that serves the limited field. A no-op tree-wise when h₀ is already g-smooth.

WhatsThePoint._auto_min_ratio Method

_auto_min_ratio(::Type{T}, n, dim) where {T}

Default geometry tree resolution from the element count: 1 / (4 * n^(1/dim)) — n^(1/3) for triangle meshes in 3D, n^(1/2) for segment loops in 2D. The 3D exponent applied to a many-segment 2D boundary floors subdivision far too coarsely, degrading leaf queries toward per-leaf linear scans.

Factor of 4 (vs. 2) ensures accurate geometry in high-curvature regions. Override with explicit min_ratio parameter if needed.

WhatsThePoint._box_corners Method

_box_corners(lo, hi) -> NTuple{2^N, SVector{N}}

The 2^N corners of a box, in bit order (bit d-1 selects hi[d]).

WhatsThePoint._box_probe_points Method

_box_probe_points(bbox_min, bbox_max)

Deterministic point probes covering a box: center, corners, and midpoints of faces/edges — 27 points in 3D, 9 in 2D (center + 4 corners + 4 edge midpoints).

WhatsThePoint._bridson_h_min Method

_bridson_h_min(node_tree, classification, spacing) -> T

Minimum spacing over non-exterior leaf centers — sets the background grid resolution for graded spacing. T is the node tree's coordinate type.

WhatsThePoint._bridson_inside Method

_bridson_inside(c, node_tree, classification, geometry) -> Bool

Domain test for a Bridson candidate, mirroring the trust rules of the per-leaf path: interior node-tree leaves are accepted outright, boundary leaves get an exact isinside check, exterior leaves (and anything outside the geometry bbox) are rejected.

WhatsThePoint._bridson_separated Method

_bridson_separated(grid, pts, rs, c, r_c) -> Bool

true when candidate c (local disk radius r_c) is at least min(r_c, r_q) away from every accepted point q. Any rejecting point lies within r_c of c, so scanning buckets within r_c suffices.

WhatsThePoint._build_knn_neighbors Method

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

Build k-nearest neighbor adjacency list from points.

WhatsThePoint._build_radius_neighbors Method

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

Build radius-based adjacency list from points.

WhatsThePoint._classify_leaf_conservative Method

_classify_leaf_conservative(tree, leaf_idx, geometry_query, tol_rel, tol_abs) -> Int8

WhatsThePoint._closest_pair Method

_closest_pair(nn_dist, nn_id, spacings, n_fixed) -> NamedTuple

Closest pair read off the sweep's per-point nearest-neighbor data (no extra search). Returns (; r, s, r_over_s, idx_a, idx_b) in snapshot-global indices. A frozen r_over_s across iterations indicates a balanced standoff; an oscillating one, an overshoot limit-cycle.

WhatsThePoint._compute_signed_distance_octree Method

Signed distance to the mesh: distance to the closest point, signed by the angle-weighted pseudonormal of the closest feature (exact for watertight, consistently outward-oriented meshes — Bærentzen & Aanæs 2005). Returns 0 only for points exactly on the surface (or at a degenerate fold where the feature pseudonormal vanishes).

WhatsThePoint._compute_signed_distance_quadtree Method

Signed distance to the boundary loops: distance to the closest point, signed by the pseudonormal of the closest feature (negative inside the domain) — the 2D specialization of the 3D triangle query.

WhatsThePoint._constrain_octree Method

_constrain_octree(id, xi, x_proposed, is_bnd, escaped, tri_indices,
                  octree, offset_dist, len_unit) -> Point

Wall rule for one point: boundary points are re-projected onto the mesh (falling back to projecting their previous position), volume points keep the proposed position while it stays inside, and escapees revert — flagged in escaped as deposition candidates.

WhatsThePoint._cull Method

_cull(pts, spacing, ratio) -> BitVector

Near-duplicate keep-mask plus the defect warning: the cull is a safety net, so any non-zero removal is surfaced.

WhatsThePoint._dedup_loop Method

Drop vertices coincident (within len_tol) with their predecessor, including the seam duplicate of the explicitly-closed convention [p₁, …, pₙ, p₁].

A zero-length segment has no direction: it contributes a zeroed normal, and the two coincident seam vertices each accumulate only one of their two adjacent segment normals, so the pseudonormal at the seam is half-built and points the wrong way over an angular sector — stray exterior points classify interior at convex corners, interior candidates get rejected at reflex ones.

WhatsThePoint._deposit_escaped! Method

_deposit_escaped!(p, tree, kq, escaped, is_bnd, tri_indices, octree, spacing,
                  deposit_ratio, offset_dist, len_unit) -> n_deposited

One deposition pass: each escaped volume point is projected onto its nearest triangle and converted to a boundary point, unless another boundary point already sits within deposit_ratio·spacing of the landing site. tree is the sweep's snapshot kd-tree and kq the neighbor count to inspect. Serial on purpose — earlier deposits must be visible to later candidates, because conversion is one-way and a parallel pass would over-deposit when a whole layer escapes in one iteration.

WhatsThePoint._dnn_cv Method

_dnn_cv(nn_dist, spacings, n_fixed) -> Real

Coefficient of variation of d_NN/s over the movable points, read off the sweep's per-point nearest-neighbor data. The quality monitor behind stall_after: it tracks the gate's binding spacing-CV metric for free. Computed in the promoted distance/spacing type.

WhatsThePoint._estimate_volume_points Method

_estimate_volume_points(node_tree, classification, spacing) -> Int

Estimate the auto max_points cap for the Orthtree algorithm (used when the caller leaves max_points unset) as ⌈1.5 · ∑ box_measure/h(x)ᴺ⌉ over non-exterior leaves — the discrete spacing integral ∫ 1/hᴺ dx with headroom for the super-1/hᴺ saturated Poisson-disk packing density (see _estimate_volume_points). It is a non-truncating ceiling, not a target.

WhatsThePoint._extract_min_spacing Method

_extract_min_spacing(spacing)

Extract minimum spacing value from spacing object (field access, no sampling).

WhatsThePoint._feature_pseudonormal Method

Pseudonormal of the feature (face / edge / vertex) the closest point lies on. Falls back to the face normal if the feature key is missing — cannot happen for features built from this mesh's own triangles, but keeps the query total.

WhatsThePoint._force_occupied_boundary! Method

_force_occupied_boundary!(classification, tree)

Demote every leaf that holds boundary elements to LEAF_BOUNDARY, whatever the point-sample classification said — downstream sampling trusts LEAF_INTERIOR leaves to need no exact inside test, so a leaf overlapping geometry must never keep that label. Shared by the triangle and segment indexes.

WhatsThePoint._generate_bridson Method

_generate_bridson(node_tree, classification, geometry, spacing, seeds,
                  max_points; factor=0.75, k_attempts=30)
    -> Vector{SVector{N,T}}

Graded Bridson Poisson-disk sampling of the domain volume with disk radius factor · h(x) (see the Orthtree docstring for the choice of factor). seeds (boundary points) initialize the advancing front and occupy the background grid so volume points keep their distance from the wall, but are not returned. The front runs until saturation or until max_points volume points exist; truncation warns because it leaves the far side of the front unfilled.

WhatsThePoint._generate_points_in_box Method

_generate_points_in_box(bbox_min, bbox_max, n, placement)

Generate n points: :random, :jittered (stratified), or :lattice (grid).

WhatsThePoint._gradient_limit_field Method

_gradient_limit_field(node_tree, leaves, h0_field, g; k, tol, max_sweeps)
    -> Vector{T}

Box-indexed g-Lipschitz envelope of h0_field, restricted to leaves. Multi-source min-plus relaxation over a k-NN graph of leaf centres: every leaf is a source carrying its own h₀, and h[i] ← min(h[i], min_j h[j] + g·dᵢⱼ) is swept to a fixpoint. The k-NN edges approximate Euclidean distance through multi-hop paths, so the result converges to the Euclidean envelope as the leaf sampling refines. Propagation crosses only the leaf graph, so thin exterior gaps between interior regions are bridged only if their centres are k-NN close (acceptable; a geodesic limiter would be exact but is not needed here).

WhatsThePoint._grid_insert! Method

_grid_insert!(grid, p)

Insert the point with the next sequential index into its bucket. Points must be inserted in index order (1, 2, …) — the link vector grows by one per call.

WhatsThePoint._guard_coarse_spacing Method

_guard_coarse_spacing(spacing, geometry, bridson_factor) -> spacing

Bridson safety net. When the finest prescribed spacing over the domain is at or above the Poisson-disk ceiling L_min/(2·bridson_factor) — i.e. a saturated front would leave the interior empty — emit a loud @warn and return a _ClampedSpacing capped at a usable baseline (L_min/10) so generation still yields a cloud. Otherwise returns spacing unchanged (the request is viable and is respected).

WhatsThePoint._leaf_class_from_signed_distance Method

_leaf_class_from_signed_distance(sd, tol) -> Int8

WhatsThePoint._leaf_spacing_field Method

_leaf_spacing_field(::Type{T}, node_tree, leaves, spacing) -> Vector{T}

Box-indexed vector of the prescribed spacing evaluated at each leaf centre (zero for non-leaf / exterior boxes). Parallel over leaves.

WhatsThePoint._loop_diagonal Method

Bounding-box diagonal of a loop — its intrinsic length scale.

WhatsThePoint._loop_signed_area Method

Shoelace signed area of a closed loop (positive = counter-clockwise).

The sum is centered on loop[1] and accumulated over differences: on absolute coordinates, a small loop far from the origin loses the whole area to cancellation (each term is O(|x|²) while the result is O(area)), and the sign — which decides the loop orientation, hence which side of the boundary is inside — becomes round-off noise.

WhatsThePoint._manifold_dim Method

Spatial dimension of a Euclidean manifold type (𝔼{2} → 2, 𝔼{3} → 3).

WhatsThePoint._maybe_kick! Method

_maybe_kick!(p, pair, state, kick_after, spacings, n_fixed, n_protected, len_unit)
    -> (state, kicked)

Kick one point of the closest pair by 0.1·s in a random direction once the pair has stayed frozen (same indices and r/s) for kick_after consecutive iterations. Prefers an index past n_protected (a volume point) and always picks a movable one (past n_fixed). state is the (pair, rs, count) tuple the caller threads through.

WhatsThePoint._near_duplicate_keep_mask Method

_near_duplicate_keep_mask(pts, spacings, ratio) -> BitVector

Greedy, order-preserving keep-mask: drop any point closer than ratio·spacing to a kept, lower-indexed point (so boundary points, indexed first, survive over volume points). The ball search at the largest cull radius sees every point inside the threshold, so the guarantee holds for clusters of any size.

WhatsThePoint._nearest_element_tree! Method

_nearest_element_tree!(point, tree, box_idx, state, update!)

Generic branch-and-bound nearest-element descent: prunes boxes farther than state.best_dist_sq, visits children nearest-first, and calls update!(element_id) for every element of each surviving leaf. update! is a callable (typically a small struct, not a closure, to stay allocation-free) that tightens state when the element beats the current best. Shared by the 3D triangle and 2D segment indexes.

WhatsThePoint._neighbor_offset Method

Axis-aligned neighbour offset for a direction code in 1:2N (odd = negative face).

WhatsThePoint._node_tree_like Method

Empty node tree spanning the same root box as the geometry tree.

WhatsThePoint._non_exterior_leaves Method

_non_exterior_leaves(node_tree, classification) -> Vector{Int}

Indices of every leaf that is not classified LEAF_EXTERIOR (interior + boundary). The shared leaf set behind the spacing-integral helpers below.

WhatsThePoint._point_box_distance_sq Method

Squared distance from point to the axis-aligned box [bbox_min, bbox_max].

WhatsThePoint._point_in_loop Method

Even-odd crossing test of p against a closed loop (loop vertices ordered).

WhatsThePoint._points_with_unit Method

_points_with_unit(raw, len_unit)

Rebuild Points from unit-stripped magnitudes, re-attaching the unit that was stripped on entry. Point(pt...) on bare numbers is not unit-agnostic — Cartesian attaches metres unconditionally — so a boundary in any other unit would exit numerically-in-its-own-unit but typed as metres, and PointCloud assembly would either throw on the CRS mismatch or mix coordinates 1000× apart.

WhatsThePoint._project_to_boundary Method

_project_to_boundary(sv, octree, offset_dist) -> (SVector, tri_idx)

Nearest point on the mesh surface, nudged offset_dist inward along the triangle normal. tri_idx == 0 means no triangle was found.

WhatsThePoint._rand_point_in_box Method

_rand_point_in_box(bbox_min, bbox_max)

Generate a random point uniformly distributed in a bounding box.

WhatsThePoint._raw_point Method

_raw_point(pt) -> SVector

Unitless coordinates of a point in its native machine type (dimension-generic).

WhatsThePoint._reconstruct_cloud Method

_reconstruct_cloud(cloud, p, tri_indices, is_bnd, n_boundary, octree, spacing, keep)

Rebuild a PointCloud after octree repel: kept points are partitioned by is_bnd into a single :boundary surface and the volume. Projected boundary points take the landing triangle's normal; imported ones (id ≤ n_boundary) keep their original area, deposited ones get spacing².

WhatsThePoint._relax! Method

_relax!(p, p_old, snap, spacing, force_model, constrain; kwargs...) -> Vector{T}

Shared relaxation loop behind both repel methods. p holds the movable points (updated in place); snap is the search snapshot whose first n_fixed entries are static and whose tail mirrors p, refreshed every rebuild_every iterations. constrain(id, xi, x_proposed) maps a proposed position to the final one (identity, or the octree wall rule). deposit!(p, method, i), when given, runs serially after each sweep. Kick targets prefer indices past n_protected. cv_target > 0 and stall_after > 0 add the quality-based stops: end when the movable points' d_NN/s CV (read off the sweep's nearest-neighbor data, no extra search) reaches the target, or has not improved for that many iterations. Returns the force-norm convergence history max_i(|F_i|·s_i)in the points' machine typeT.

WhatsThePoint._resolve_alg Method

_resolve_alg(bnd, spacing, alg) -> AbstractNodeGenerationAlgorithm

Pick (or validate) the algorithm for a boundary, dispatching on its manifold. Single place for the decision, so the PointBoundary and PointCloud entry points cannot drift, and so every rejected combination throws an instructive ArgumentError rather than deferring to a raw MethodError inside _discretize_volume dispatch.

WhatsThePoint._safe_direction Method

_safe_direction(xi, xj, r) -> Vec

Unit direction from xj to xi; a random unit vector when r == 0, avoiding the 0/0 NaN that traps coincident points.

WhatsThePoint._segment_feature_pseudonormal Method

Pseudonormal of the feature (segment / vertex) the closest point lies on.

WhatsThePoint._segment_len_tol Method

_segment_len_tol(diagonal)

Absolute length below which a segment counts as degenerate. b - a carries an absolute round-off of order eps(T) times the coordinate magnitudes, so the floor must scale with a bbox diagonal: an absolute eps(T) would declare every segment of a finely-sampled Float32 loop degenerate (zeroed normals, endpoint snapping), and would miss true duplicates far from the origin.

Callers pass the diagonal of the individual loop, not of the whole boundary: one global tolerance lets a large outer boundary silently decimate — or collapse outright — a finely-sampled small hole (multi-scale domains).

WhatsThePoint._signed_distance Method

Signed distance from p to the boundary loops (negative inside).

WhatsThePoint._signed_distance Method

Signed distance from p to the mesh (negative inside).

WhatsThePoint._signed_volume Method

Signed volume of a closed triangle mesh (divergence theorem, Σ dot(v1, v2 × v3) / 6). Positive iff the winding orients normals outward. Catches globally inside-out meshes, which has_consistent_normals cannot (a perfectly consistent but inverted mesh classifies its complement as interior — exactly the 2026-06-11 cavity corruption #2). Only meaningful for closed surfaces.

WhatsThePoint._triangle_axis_test Method

_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

WhatsThePoint.add_box! Method

add_box!(tree, cell::SVector{N,Int}, level, parent_idx) -> Int

Append a new leaf box at (cell, level); grows arrays as needed. Returns index.

WhatsThePoint.all_boxes Method

all_boxes(tree) -> Vector{Int}

WhatsThePoint.all_leaves Method

all_leaves(tree) -> Vector{Int}

WhatsThePoint.any_leaf_overlapping Method

any_leaf_overlapping(tree, bbox_min, bbox_max, predicate) -> Bool

true if any leaf overlapping [bbox_min, bbox_max] satisfies predicate. Prunes non-overlapping subtrees for O(log L) expected cost.

WhatsThePoint.balance_octree! Method

balance_octree!(tree, criterion::SubdivisionCriterion; redistribute! = nothing)

Enforce the 2:1 balance constraint across the tree (dimension-agnostic). Uses can_subdivide (physical limits only), not should_subdivide.

subdivide! does not move a leaf's element_lists entries into the new children — queries only scan leaves, so elements of a balance-subdivided box would vanish from the queryable tree. Callers whose queries read element lists (e.g. TriangleOctree's nearest-triangle search) must pass redistribute!(tree, box_idx), invoked right after each forced subdivision to push the parent's elements into the intersecting children.

WhatsThePoint.bounding_box Method

bounding_box(tree) -> (min_corner, max_corner)

Root box bounds.

WhatsThePoint.box_bounds Method

box_bounds(tree, box_idx) -> (min_corner, max_corner)

WhatsThePoint.box_center Method

box_center(tree, box_idx) -> SVector{N}

center = origin + (2*cell + 1) * box_size / 2, box_size = root_size / level.

WhatsThePoint.box_size Method

box_size(tree, box_idx) -> Real

Edge length of the box: root_size / level.

WhatsThePoint.build_node_octree Method

build_node_octree(geometry, spacing, alpha, node_min_ratio)

Build a spacing-driven node tree from an existing geometry index.

Creates a new SpatialTree that subdivides based on a spacing function, enabling spacing-aware point distribution. The node tree is: 2. Recursively subdivided using SpacingCriterion

3. Balanced to maintain 2:1 refinement ratio

4. Independent of the geometry-index resolution

Arguments

• geometry: Geometry index (TriangleOctree in 3D, SegmentQuadtree in 2D)

• spacing: Spacing function (e.g., ConstantSpacing, BoundaryLayerSpacing)

• alpha: Subdivision aggressiveness (h_box ≤ alpha * h_spacing)

• node_min_ratio: Minimum box size ratio relative to domain

Returns

SpatialTree{N, Int, T} with spacing-driven subdivision, where T is the geometry index's coordinate type (the source CRS machine type)

Example

tri_octree = TriangleOctree(mesh; classify_leaves=true)
spacing = BoundaryLayerSpacing(points; at_wall=0.5m, bulk=5m, layer_thickness=2m)
node_tree = build_node_octree(tri_octree, spacing, 1.0, 1e-6)

WhatsThePoint.children Method

Range of child indices of box (empty range if box is a leaf).

WhatsThePoint.classify_leaves! Method

classify_leaves!(tree, geometry_query; tolerance_relative, tolerance_absolute) -> Vector{Int8}

Conservative (2^N + 1)-point classification (center + corners) per leaf, using geometry_query(point::SVector{N,T}, tol::T) -> Int8.

WhatsThePoint.classify_node_octree Method

classify_node_octree(node_tree, geometry)

Classify node tree leaves as interior, boundary, or exterior using the geometry index's classify_point seam method. Correctness of downstream sampling (skipping isinside on LEAF_INTERIOR points) relies on the geometry-bbox early return inside classify_point, which prevents sign-vote flips from promoting far-exterior leaves into LEAF_INTERIOR.

Returns

Vector of Int8 classifications indexed by node-tree box index.

WhatsThePoint.classify_point Method

Shared classification of a point against a geometry index: bbox fast path, cached leaf classification for INTERIOR/EXTERIOR dispatch, exact signed distance only for BOUNDARY leaves. tol expands the bounding box for the exterior fast path (0 for exact bbox checks, positive for conservative classification). Implementations plug in through _signed_distance.

WhatsThePoint.closest_point_on_segment_feature Method

closest_point_on_segment_feature(p, a, b, len_tol_sq=zero(T)) -> (closest_point, feature)

Closest point on segment a → b, plus the feature it lies on: FEATURE_FACE (segment interior), FEATURE_VERTEX_1 (a), or FEATURE_VERTEX_2 (b).

len_tol_sq is the squared length below which the segment is degenerate and a is returned. The default is an exact-zero guard against 0/0 only — a tiny nonzero denom is harmless (t just clamps to an endpoint), while any geometric threshold here would need to track the loop scale, which a single per-index value cannot do for multi-scale domains (a fine hole's short segments would snap to endpoints under the outer boundary's tolerance).

WhatsThePoint.closest_point_on_triangle Method

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: 2. Projecting P onto the triangle plane

3. Computing barycentric coordinates

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

WhatsThePoint.closest_point_on_triangle_feature Method

closest_point_on_triangle_feature(P, v1, v2, v3) -> (SVector{3,T}, Int8)

Same as closest_point_on_triangle but also reports which triangle feature the closest point lies on (FEATURE_FACE, FEATURE_VERTEX_k, FEATURE_EDGE_jk). The feature is a free by-product of the Ericson region classification and selects the correct angle-weighted pseudonormal for exact signed-distance computation (Bærentzen & Aanæs 2005).

WhatsThePoint.domain_bounds Method

Axis-aligned bounds of the boundary geometry.

WhatsThePoint.domain_bounds Method

Axis-aligned bounds of the boundary geometry.

WhatsThePoint.emptyspace Method

emptyspace(testpoint, points)

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

WhatsThePoint.find_boxes_at_coords Method

3D convenience overload preserving the historical (i,j,k,N) signature.

WhatsThePoint.find_boxes_at_coords Method

find_boxes_at_coords(tree, target_cell::SVector{N,Int}, target_level) -> Vector{Int}

Box(es) covering (target_cell, target_level): exact match, the coarser box covering it, or Int[] if absent.

WhatsThePoint.find_leaf Method

find_leaf(tree, point::SVector{N}) -> Int

Leaf box containing point, O(log n) from root.

WhatsThePoint.find_neighbor Method

find_neighbor(tree, box_idx, direction) -> Vector{Int}

Neighbour(s) across the face direction (1:2N). Handles the 2:1 level difference: returns [same-level], or the finer boxes covering the face, or Int[] at a boundary.

WhatsThePoint.geometry_tree Method

Spatial subdivision tree of a geometry index.

WhatsThePoint.global_to_local Method

global_to_local(boundary::PointBoundary, global_idx::Int) -> (Symbol, Int)

Convert a boundary-global index to a (surface_name, local_index) tuple.

WhatsThePoint.global_to_local Method

global_to_local(cloud::PointCloud, global_idx::Int) -> (Symbol, Int)

Convert a cloud-global index to a (component, local_index) tuple. Returns (:volume, local_idx) for volume indices, or (surface_name, local_idx) for boundary indices.

WhatsThePoint.has_children Method

has_children(tree, box_idx) -> Bool

WhatsThePoint.is_leaf Method

is_leaf(tree, box_idx) -> Bool

WhatsThePoint.leaf_classes Method

Per-leaf classification cache of a geometry index (nothing if not built).

WhatsThePoint.local_to_global Method

local_to_global(boundary::PointBoundary, name::Symbol, local_idx::Int) -> Int

Convert a surface-local index to a boundary-global index.

WhatsThePoint.local_to_global Method

local_to_global(cloud::PointCloud, name::Symbol, local_idx::Int) -> Int

Convert a surface-local index to a cloud-global index.

WhatsThePoint.needs_balancing Method

needs_balancing(tree, box_idx) -> Bool

true if subdividing this leaf would leave a 2-level jump with any face neighbour (i.e. a neighbour that already has grandchildren).

WhatsThePoint.num_elements Method

num_elements(tree) -> Int

WhatsThePoint.project_to_boundary Method

Project p onto the boundary, returning (projected point, boundary element id). Delegates to the repel projection kernel (defined later in the module).

WhatsThePoint.segment_box_intersection Method

Segment–AABB overlap via Liang–Barsky parameter clipping.

WhatsThePoint.subdivide! Method

subdivide!(tree, box_idx) -> UnitRange{Int}

Subdivide a leaf into 2^N children (contiguous block). Returns the child range. Child m (0-based) sits at 2*cell + bits(m), where bit d-1 of m is the offset along axis d — generalizing the standard octant ordering.

WhatsThePoint.surface_offset Method

surface_offset(boundary::PointBoundary, name::Symbol) -> Int

Return the global offset for the start of surface name within the boundary. The first point of surface name has boundary-global index surface_offset(...) + 1.

WhatsThePoint.triangle_box_intersection Method

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 2. Translate triangle and box so box is centered at origin

3. Test each potential separating axis

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

WhatsThePoint.volume_to_global Method

volume_to_global(cloud::PointCloud, local_idx::Int) -> Int

Convert a volume-local index to a cloud-global index.