API Reference

cluster(smld::SMLMData.BasicSMLD, cfg::AbstractClusterConfig) -> (smld_out, ClusterInfo)

Cluster the localizations in smld using the backend selected by the concrete type of cfg. The input smld is not modified — each backend deep-copies the input emitters and writes cluster labels onto the copy's emitter.id (0 marks noise, 1..K mark distinct clusters). If cfg.remove_unclustered is true the returned SMLD contains only clustered emitters.

Dispatch on AbstractClusterConfig alone has no implementation — concrete backends supply a method specialized on their own config type. Calling cluster with an unsupported config raises a clear error.

See also: AbstractClusterConfig, ClusterInfo.

cluster(smld::BasicSMLD, cfg::PointHysteresisConfig;
        seed::AbstractVector{Bool}, support::AbstractVector{Bool})
    -> (smld_out, ClusterInfo)

Run point-graph hysteresis on smld using cfg's knobs and the per-cell seed and support boolean masks.

Required keyword arguments

• seed: high-confidence foreground points (length must equal length(smld.emitters)). Components must contain at least one to be kept.
• support: candidate foreground points (same length). BFS only crosses through support points. seed must imply support (every seed point is also a support point); violations raise ArgumentError.

Algorithm

1. Build a KDTree over the (per-dataset or pooled) emitter coordinates.
2. Iterate emitters in input order. Whenever an unvisited support point is found, BFS through the support set following kNN edges (graph_k neighbors per node, excluding self). Mark every reached node visited.
3. After the BFS, if the component contains at least one seed point AND has at least cfg.min_points members, assign it a fresh cluster id (local within the per-dataset group, V3). Otherwise leave its members at id = 0 (the "noise" contract from AbstractClusterConfig — input ids are zeroed on the deepcopy before BFS so prior labels never leak).

Convention

• The config is non-mutating: emitters are deep-copied (V9) and ids are zeroed on the copy before BFS, so unclustered emitters come out as id = 0 regardless of input id values.
• Asymmetric kNN BFS: the BFS only crosses edges where the current node has the neighbor in its k-NN. For tight clusters with graph_k << N, outlier points may be unreachable; this is correct behavior matching the reference implementation, but at small N the user should prefer larger graph_k (or use full-mesh min(N-1, ...)).

See also: PointHysteresisConfig, LocalContrastFeature.

AbstractClusterConfig <: SMLMData.AbstractSMLMConfig

Abstract supertype for SMLMClustering backend configurations.

Each backend defines a concrete subtype (e.g. DBSCANConfig, HDBSCANConfig, VoronoiConfig, HierarchicalConfig) and adds a cluster(smld, cfg) method specialized on it.

Shared fields (expected on every concrete subtype)

• min_points::Int: minimum points for a valid cluster
• use_3d::Bool: whether the z-coordinate participates in clustering
• per_dataset::Bool: if true, cluster within each dataset independently so that (dataset, id) uniquely identifies a cluster across a multi-dataset SMLD
• remove_unclustered::Bool: if true, drop emitters with id == 0 from the returned SMLD (noise/rejected localizations)

Algorithm-specific fields (eps_nm, min_cluster_size, density_factor, cut_nm, ...) live on the concrete subtype that needs them.

ClusterInfo <: SMLMData.AbstractSMLMInfo

Secondary output from cluster() — summary of a clustering run.

Fields

• n_locs_in::Int: number of input localizations
• n_clustered::Int: number of localizations assigned to a cluster (id > 0)
• n_noise::Int: number of localizations tagged as noise (id == 0)
• n_clusters::Int: number of distinct clusters formed
• cluster_sizes::Vector{Int}: size of each cluster, indexed by cluster id (cluster_sizes[k] is the size of cluster k); length equals n_clusters
• algorithm::Symbol: backend identifier (:dbscan, :hdbscan, :voronoi, :hierarchical)
• elapsed_s::Float64: wall-clock time spent in the cluster call, in seconds

Example

(smld_out, info) = cluster(smld, DBSCANConfig(eps_nm=50.0, min_points=5))
println("$(info.n_clustered)/$(info.n_locs_in) clustered into $(info.n_clusters) clusters")

DBSCANConfig(; eps_nm, min_points=5, use_3d=false, per_dataset=true, remove_unclustered=false)

Configuration for DBSCAN clustering of SMLM localizations.

Fields

• eps_nm::Float64: neighborhood radius in nanometers. Coordinates on AbstractEmitter subtypes are in microns; the backend converts internally.
• min_points::Int = 5: minimum number of points in an ε-neighborhood for a point to be a core point (classical DBSCAN minPts). Also used as the minimum cluster size.
• use_3d::Bool = false: if true, cluster in (x, y, z); requires Emitter3DFit emitters. Otherwise cluster in (x, y).
• per_dataset::Bool = true: if true, cluster within each dataset index independently so that (dataset, id) uniquely identifies a cluster in a multi-dataset SMLD. If false, all emitters are clustered together and id alone identifies the cluster.
• remove_unclustered::Bool = false: if true, emitters tagged as noise (id == 0) are dropped from the returned SMLD.

Example

cfg = DBSCANConfig(eps_nm=50.0, min_points=5)
(smld_out, info) = cluster(smld, cfg)

See also: AbstractClusterConfig, ClusterInfo, cluster.

PrecisionDBSCANConfig(; nsigma, min_points=5, use_3d=false, per_dataset=true,
                      remove_unclustered=false)

Configuration for precision-weighted DBSCAN of SMLM localizations: the neighbor test is ‖pᵢ − pⱼ‖ < nsigma · (σ_effᵢ + σ_effⱼ), where each emitter's σ_eff is the geometric mean of its per-axis localization precisions (√(σ_x·σ_y) in 2D, ∛(σ_x·σ_y·σ_z) in 3D). Unlike DBSCANConfig's fixed eps_nm, the neighborhood adapts to each localization's uncertainty.

Fields

• nsigma::Float64: neighbor radius in units of the summed precision σ_effᵢ + σ_effⱼ.
• min_points::Int = 5: classical DBSCAN minPts — the minimum neighborhood size (the point itself plus its active neighbors) for a core point — and the minimum cluster size (clusters smaller than min_points become noise). Identical in meaning and default to DBSCANConfig.
• use_3d::Bool = false: cluster in (x, y, z) using σ_z; requires Emitter3DFit.
• per_dataset::Bool = true: cluster within each dataset index independently.
• remove_unclustered::Bool = false: drop noise emitters (id == 0) from the output.

Cluster ids are written to emitter.id (0 = noise, 1..K = clusters); returns (smld_out, ClusterInfo) with algorithm = :precision_dbscan.

For the lower-level reuse-the-graph primitive (build once, relabel many times with varying σ_eff/nsigma), use build_precision_neighbor_graph + precision_dbscan_labels directly.

See also: DBSCANConfig, AbstractClusterConfig, cluster.

HDBSCANConfig(; min_points=5, min_cluster_size=nothing, knn_graph_k=30,
                cluster_selection_method=:eom, allow_single_cluster=false,
                halo_trim_frac=0.10, use_3d=false, per_dataset=true,
                remove_unclustered=false)

Configuration for HDBSCAN clustering of SMLM localizations. Pure-Julia implementation of Campello/Moulavi/Sander 2013.

Fields

• min_points::Int = 5: k for the core-distance computation. Larger values produce more conservative clusters.
• min_cluster_size::Union{Int,Nothing} = nothing: minimum cluster size in the condensed tree. When nothing, defaults to min_points.
• knn_graph_k::Int = 30: width k' of the sparse k'-NN graph used as the MST scaffold. Must be large enough that the resulting MST is connected; the backend errors out (with a clear message) if knn_graph_k is too small.
• cluster_selection_method::Symbol = :eom: how to pick clusters from the condensed tree — :eom (excess of mass; the canonical HDBSCAN choice) or :leaf (return all leaves of the condensed tree, like flat DBSCAN at varying ε).
• allow_single_cluster::Bool = false: when true, the root cluster (the whole connected mass) is a candidate for EOM selection. Useful for data that is essentially one tight blob with no internal real splits — the canonical EOM rule otherwise returns zero clusters in that case. Default matches the Python hdbscan package default.
• halo_trim_frac::Float64 = 0.10: halo-trim fraction. A point that fell out of its cluster individually (a transient / halo point) is kept only if it survived past this fraction of the cluster's λ-life; weakly-attached points that peeled off near the cluster's birth (density-connected background the MST routed onto the cluster's side) are dropped to noise. 0 disables trimming (raw condensed-tree labels — every point under a selected branch is a member); larger values trim more aggressively. Calibrate against the physical edge (the radius where cluster density crosses background); for the validated default the diffuse-cluster member radius matches that edge.
• use_3d::Bool = false: cluster in (x, y, z); requires 3D emitters.
• per_dataset::Bool = true: cluster within each dataset index independently.
• remove_unclustered::Bool = false: drop noise emitters from the output SMLD.

Outputs

cluster(smld, ::HDBSCANConfig) returns (smld_out, info) matching the shared cluster interface. Per-cluster persistence (stability) and birth lambda are written to smld_out.metadata under the keys "hdbscan_cluster_persistence" and "hdbscan_cluster_lambda_birth" — flat vectors of length info.n_clusters in cluster-id order (and concatenated across datasets when per_dataset=true).

Example

cfg = HDBSCANConfig(min_points=10, min_cluster_size=20, knn_graph_k=50)
(smld_out, info) = cluster(smld, cfg)
persistence = smld_out.metadata["hdbscan_cluster_persistence"]

See also: AbstractClusterConfig, ClusterInfo, cluster.

HierarchicalConfig(; cut_threshold=nothing, n_clusters=nothing, linkage=:ward,
                    min_points=5, use_3d=false, per_dataset=true,
                    remove_unclustered=false)

Configuration for agglomerative hierarchical clustering of SMLM localizations.

Exactly one of cut_threshold or n_clusters must be supplied.

Fields

• cut_threshold::Union{Float64,Nothing} = nothing: dendrogram cut height. Units depend on linkage: for distance-based linkages (:single, :complete, :average) the value is in nanometers and is converted to μm internally (h = cut_threshold / 1000.0). For :ward the dendrogram height is a variance-increase cost (roughly μm²) and is passed through without conversion — there is no meaningful "nm" interpretation under Ward.
• n_clusters::Union{Int,Nothing} = nothing: cut the dendrogram to produce exactly this many clusters (before min_points filtering). Natural for Ward, where cut_threshold has no intuitive unit. Mutually exclusive with cut_threshold.
• linkage::Symbol = :ward: linkage criterion — :single, :complete, :average, or :ward. Ward minimizes within-cluster variance.
• min_points::Int = 5: clusters with fewer than min_points members after cutting are relabeled noise (id = 0). Remaining clusters are renumbered compactly 1..K.
• use_3d::Bool = false: include z-coordinate in distance calculation.
• per_dataset::Bool = true: cluster within each dataset independently so that (dataset, id) uniquely identifies a cluster in a multi-dataset SMLD.
• remove_unclustered::Bool = false: drop noise emitters (id == 0) from the output.

Example

# Distance-based: cut at 200 nm under single linkage.
cfg = HierarchicalConfig(cut_threshold=200.0, linkage=:single)
(smld_out, info) = cluster(smld, cfg)

# Ward linkage: specify number of clusters directly.
cfg2 = HierarchicalConfig(n_clusters=3, linkage=:ward)
(smld_out, info) = cluster(smld, cfg2)

See also: AbstractClusterConfig, ClusterInfo, cluster.

VoronoiConfig(; density_factor=2.0, min_points=5, use_3d=false,
                per_dataset=true, remove_unclustered=false)

Configuration for Voronoi-tessellation-based (SR-Tesseler) clustering of SMLM localizations.

Algorithm

1. Build the Voronoi tessellation of the localization coordinates (per dataset when per_dataset=true), clipped to the convex hull so every generator has a finite cell.
2. A localization is dense when its cell area is smaller than mean_cell_area / density_factor (equivalently, its local density exceeds density_factor × mean_density).
3. Dense points that are Delaunay-adjacent are merged into clusters via connected components.
4. Clusters with fewer than min_points members are relabeled noise (id = 0); the rest are renumbered compactly 1..K within the group.

Fields

• density_factor::Float64 = 2.0: density threshold multiplier. Higher values require stronger local density to qualify as a cluster member (smaller area threshold = fewer dense points).
• min_points::Int = 5: minimum cluster size; smaller connected components become noise.
• use_3d::Bool = false: must be false. 3D Voronoi clustering is not supported — DelaunayTriangulation.jl is 2D only.
• per_dataset::Bool = true: cluster within each dataset independently so that (dataset, id) uniquely identifies a cluster.
• remove_unclustered::Bool = false: drop noise emitters (id == 0) from the returned SMLD.

Example

cfg = VoronoiConfig(density_factor=2.0, min_points=5)
(smld_out, info) = cluster(smld, cfg)

See also: AbstractClusterConfig, ClusterInfo, cluster.

MRFDensityClusterConfig(; n_regimes=2, regime_thresholds=nothing,
                          regime_gaussians=nothing,
                          density_estimator=:voronoi, density_k=20,
                          smoothness_lambda=nothing,
                          graph_kind=:delaunay, graph_k=8,
                          inference=:icm, icm_iters=50,
                          min_points=5, use_3d=false,
                          per_dataset=true, remove_unclustered=false)

Adaptive-density clustering via per-emitter density → multi-component GMM regime assignment → multi-class Potts MRF refinement → connected components on the foreground.

Designed for SMLM data with multiple density regimes (e.g. one low-density background + one or more higher-density structures) where a single global ε would either over-merge or over-split. The MRF smoothness term enforces spatial coherence — points surrounded by foreground neighbors stay foreground even if their individual density is borderline (fixes "missing middles"); isolated tight knots in a low-density sea get pulled back to background (fixes spurious-small-cluster artifacts).

Density estimator

Voronoi density (1/cell area) is the simplest local estimator but its noise scale is large (σlog ≈ 1) because it averages over a single Voronoi cell. For thin elongated structures (width comparable to nearest-neighbor distance, e.g. dSTORM fibers narrower than ~100 nm at high density) Voronoi cells span the patch boundary and inflate, dragging interior densities toward background and producing a false-negative band at patch interiors. The kNN density estimator `ρk = k / (π · rk²)` (Loftsgaarden-Quesenberry 1965) integrates over the k nearest neighbors and reduces noise to σlog ≈ 1/√k, recovering patch-interior detection at the cost of some boundary blur. Default density_estimator=:voronoi preserves backward compatibility; switch to :knn (with density_k=20-30) for thin-fiber-bearing datasets.

Fields

• n_regimes::Int = 2: number of density regimes. Lowest is treated as background/noise; higher regimes form the foreground.
• regime_thresholds::Union{Nothing, Vector{Float64}} = nothing: optional explicit log-density thresholds, length n_regimes - 1, sorted ascending. When provided, GMM is bypassed and points are binned directly with hard unaries. This is useful as a conservative low-contrast stabilizer, but borderline points cannot be rescued by spatial neighbors.
• regime_gaussians::Union{Nothing, NamedTuple} = nothing: optional calibrated Gaussian emission model (means, vars, weights) in log-density space. When provided, GMM is bypassed but soft unaries are preserved, so the Potts prior can still flip borderline interior points.
• density_estimator::Symbol = :voronoi: per-emitter density estimator. :voronoi uses 1/Voronoi-cell-area; :knn uses k / (π · rk²) where `rkis thedensity_k`-th nearest-neighbor distance.
• density_k::Int = 20: k for the kNN density estimator when density_estimator=:knn. Ignored for :voronoi. Larger k → lower variance, more boundary blur; recommended range 10-50.
• smoothness_lambda::Union{Nothing, Float64} = nothing: MRF smoothness weight. When nothing, auto-set per group to max(1e-6, MAD(U_max - U_min)) where MAD is the median absolute deviation of the per-emitter unary range.
• graph_kind::Symbol = :delaunay: neighbor graph for both the MRF and the CC step. :delaunay reuses the tessellation from step 1 (free, requires density_estimator=:voronoi); :knn builds a symmetrized k-NN graph (uses graph_k). When density_estimator=:knn AND graph_kind=:delaunay, a Delaunay triangulation is still computed to produce the neighbor graph.
• graph_k::Int = 8: k for the kNN graph when graph_kind=:knn. Ignored for :delaunay.
• inference::Symbol = :icm: MRF inference algorithm. Only :icm (Iterated Conditional Modes) is supported in v1; future graph-cut inference would land here.
• icm_iters::Int = 50: maximum ICM passes; terminates early when no point changes label in a full pass.
• min_points::Int = 5: minimum cluster size after CC; smaller components are demoted to noise.
• use_3d::Bool = false: must be false. 3D Voronoi tessellation is not supported (DelaunayTriangulation.jl is 2D only); passing true raises ArgumentError.
• per_dataset::Bool = true: when true, run the full pipeline per dataset so each cell's density distribution is fit independently.
• remove_unclustered::Bool = false: drop noise emitters from the output.

Outputs

Standard (smld_out, info) tuple. Per-cluster cluster ids in emitter.id (0 = noise / lowest regime / sub-min component, 1..K = cluster). HDBSCAN-style metadata stamped onto smld_out.metadata:

• metadata["mrf_regime_per_emitter"]::Vector{Int} of length n_locs_in, in original emitter order. Values: 0 (ungroupable / group too small), 1 (lowest density / background), ..., n_regimes (highest density).
• metadata["mrf_lambda_used"]::Vector{Float64}: per-group λ used.
• metadata["mrf_regime_means"]::Vector{Vector{Float64}}: per-group Gaussian component means (sorted ascending) in log-density space; vector of NaNs when hard regime_thresholds were provided.

Example

# 2-regime auto: GMM finds the foreground/background split per dataset
cfg = MRFDensityClusterConfig()
(smld_out, info) = cluster(smld, cfg)
regimes = smld_out.metadata["mrf_regime_per_emitter"]

# 3-regime with manual thresholds (e.g. learned from training data)
cfg2 = MRFDensityClusterConfig(n_regimes = 3,
                               regime_thresholds = [3.5, 5.0],
                               min_points = 10)

See also: AbstractClusterConfig, ClusterInfo, cluster, VoronoiDensityConfig (the read-only sibling that only exposes per-emitter densities without clustering).

PointHysteresisConfig(; graph_k=12, min_points=100, use_3d=false,
                        per_dataset=false, remove_unclustered=false)

Configuration for point-graph hysteresis seed-and-grow clustering.

Holds only knobs — seed and support boolean masks are passed as required keyword arguments to cluster, not on the config. The same PointHysteresisConfig instance is reusable across many SMLDs.

Fields

• graph_k::Int = 12: degree of the kNN graph used for the BFS.
• min_points::Int = 100: minimum component size for a cluster to be kept. Components with <min_points members or no seed presence are dropped to noise (id = 0).
• use_3d::Bool = false: build the kNN graph in (x, y, z); requires Emitter3DFit emitters.
• per_dataset::Bool = false: when true, cluster within each dataset independently. Cluster ids are local to each dataset's namespace so (dataset, id) is the unique identifier (V3).
• remove_unclustered::Bool = false: if true, drop emitters with id = 0 from the returned SMLD.

Example

# Build the per-cell feature, then thresholds, then hysteresis.
(_, info_f) = cluster_statistics(smld, LocalContrastFeature(density_k=200,
                                                            background_k=2000))
contrast = info_f.extras[:contrast_per_emitter]
fine = info_f.extras[:log_density_per_emitter]
fine_floor = quantile(filter(isfinite, fine), 0.35)
seed = isfinite.(contrast) .& isfinite.(fine) .& (contrast .> 0.25) .& (fine .> fine_floor)
support = isfinite.(contrast) .& isfinite.(fine) .& (contrast .> -0.05) .& (fine .> fine_floor)

cfg = PointHysteresisConfig(graph_k = 12, min_points = 150)
(smld_out, info) = cluster(smld, cfg; seed = seed, support = support)

See also: AbstractClusterConfig, LocalContrastFeature, MRFDensityClusterConfig (the GMM+Potts alternative), ClusterInfo.

calibrate_regime_gaussians(smld; n_regimes=2, density_estimator=:knn,
                            density_k=20) -> NamedTuple

Fit a calibrated soft Gaussian emission model in log-density space. The returned named tuple has fields (means, vars, weights) and can be passed directly to MRFDensityClusterConfig(regime_gaussians = ...).

Unlike regime_thresholds, this preserves soft unary costs: U[i, k] = -log(w_k * N(log_rho[i] | μ_k, σ_k²)). That makes the MRF behave more like a spatial hidden-state model: an emitter slightly below the foreground/background decision boundary can still become foreground when its neighbors provide enough Potts support.

See also: calibrate_regime_thresholds, MRFDensityClusterConfig.

calibrate_regime_thresholds(smld; n_regimes=2, density_estimator=:knn,
                             density_k=20) -> Vector{Float64}

Derive log-density regime-boundary thresholds from a calibration SMLD by fitting a 1D Gaussian mixture on its per-emitter log densities. The returned Vector{Float64} of length n_regimes - 1 is sorted ascending and is in the exact shape MRFDensityClusterConfig.regime_thresholds accepts — pass it directly into the config to bypass GMM on subsequent query ROIs.

Use case

For datasets with stable density structure across ROIs (same imaging conditions, same labeling protocol), calibrating once on a trusted ROI and reusing the thresholds across query ROIs is faster (skips per-ROI EM) and more reproducible than letting the GMM auto-fit on each ROI independently. The recommendation is especially strong at the low-contrast end of the density-ratio band (Round 013 / KB V13): when the high/low log-density modes are close together, GMM EM on a single ROI can pick a degenerate split that the Potts smoothness prior then amplifies into a uniform-label collapse — calibration-derived thresholds avoid that failure mode by fixing the regime boundary at the value learned from a higher-quality fit.

Arguments

• smld::SMLMData.BasicSMLD: the calibration ROI. The helper computes per-emitter log densities globally (no per-dataset split — the returned thresholds are a single global vector that applies uniformly across all groups when passed to MRFDensityClusterConfig.regime_thresholds).
• n_regimes::Int = 2: number of density regimes (≥ 2). Returns n_regimes - 1 thresholds.
• density_estimator::Symbol = :knn: density estimator for log-ρ. :knn uses density_k nearest neighbors (recommended for thin-fiber-bearing datasets per V12); :voronoi uses 1/cell-area.
• density_k::Int = 20: k for the kNN estimator. Ignored for :voronoi.

Boundary computation

Each threshold is the analytic Bayes decision boundary between consecutive GMM components (sorted ascending by mean), accounting for unequal variances and weights. Falls back to the midpoint between consecutive means when the quadratic discriminant is negative (one component dominates everywhere).

Returns

Vector{Float64} of length n_regimes - 1, sorted ascending. Pass directly to MRFDensityClusterConfig(regime_thresholds = ...).

Errors

• ArgumentError if n_regimes < 2, density_estimator is unknown, density_k < 1, or the calibration SMLD has fewer valid emitters than n_regimes (after filtering NaN densities).
• ErrorException if the GMM fit fails (variance collapse, all-equal densities). Try a different calibration ROI or fewer regimes.

Example

# Fit thresholds on a high-quality calibration ROI.
cal_smld = load_calibration_data()
thr = calibrate_regime_thresholds(cal_smld; n_regimes = 2,
                                   density_estimator = :knn,
                                   density_k = 20)

# Apply on query ROIs without re-fitting.
cfg = MRFDensityClusterConfig(n_regimes = 2,
                              density_estimator = :knn,
                              density_k = 20,
                              regime_thresholds = thr)
(smld_out, info) = cluster(query_smld, cfg)

See also: MRFDensityClusterConfig, cluster.

PrecisionNeighborGraph

Geometry-only neighbor cache for precision-weighted DBSCAN, built once by build_precision_neighbor_graph and reused across many precision_dbscan_labels calls with varying σ_eff / nsigma.

Treat it as an opaque handle: build it with build_precision_neighbor_graph and pass it back to precision_dbscan_labels. Its fields (n, offsets, neighbors, dists, max_radius, dims) are an implementation detail — a CSR store of every candidate pair with raw Euclidean d ≤ max_radius (full adjacency, both directions) — and may change between releases; you should not need to read them.

The cache is a valid superset for a label pass iff nsigma·(σ_eff_i + σ_eff_j) ≤ max_radius for every pair; the sufficient rule is max_radius ≥ nsigma · 2 · maximum(σ_eff) (checked at label time when check_superset = true).

build_precision_neighbor_graph(coords::AbstractMatrix{<:Real}, max_radius::Real)
    -> PrecisionNeighborGraph

Build the geometry-only neighbor cache. coords is a D×N matrix (columns are points, D ∈ {2, 3}) in whatever length unit you will also use for max_radius and σ_eff. Every pair within max_radius (raw Euclidean) is cached with its distance; the σ-weighted threshold is applied later, per call, by precision_dbscan_labels.

The per-point range queries are the expensive part and are run in parallel (@threads); each thread writes only its own point's neighbor list, so no per-thread scratch is needed and the result is independent of thread scheduling. Neighbor lists are sorted by index, making the CSR — and hence the labels — reproducible.

precision_dbscan_labels(g::PrecisionNeighborGraph, σ_eff, nsigma;
                        min_points = 1, check_superset = true) -> Vector{Int}
precision_dbscan_labels!(labels, g, σ_eff, nsigma; min_points = 1, check_superset = true)

Label the points of a prebuilt PrecisionNeighborGraph by re-thresholding the cached pairs: an edge (i, j) is active iff g.dists < nsigma·(σ_eff[i] + σ_eff[j]). σ_eff (length g.n) and nsigma may change on every call without rebuilding g — that reuse is the point of the primitive.

• min_points == 0: connected components of the active graph, via union-find. The partition is order-free / bit-exact — independent of thread scheduling and edge order. Every point is labeled 1..K (a singleton is its own component; no noise).
• min_points ≥ 1: core-point DBSCAN with the classical (self-inclusive) minPts, identical to DBSCANConfig. A point is core iff its neighborhood — itself plus its active neighbors — has ≥ min_points members (i.e. active degree ≥ min_points − 1); core points sharing an active edge merge; a non-core border point joins the lowest-id adjacent core cluster (deterministic tie-break); a non-core point with no active core neighbor is noise (0).

Cluster ids are canonical: 1..K in ascending order of first appearance by point index, so the integer labels themselves are reproducible.

check_superset (default true) asserts nsigma·2·maximum(σ_eff) ≤ g.max_radius and errors if the cache is too tight to contain every active pair (see PrecisionNeighborGraph). The mutating form fills a caller-supplied labels::Vector{Int} of length g.n.

precision_dbscan_labels!(labels, g, σ_eff, nsigma; min_points=1, check_superset=true) -> labels

In-place form of precision_dbscan_labels: fills the caller-supplied labels::Vector{Int} (which must have length g.n) and returns it. Use this in a hot loop that reuses one PrecisionNeighborGraph across many calls to avoid reallocating the label vector each time.

cluster_statistics(smld::SMLMData.BasicSMLD, cfg::AbstractStatisticsConfig)
    -> (smld, ClusterStatisticsInfo)

Compute a spatial / clustering statistic on smld using the backend selected by the concrete type of cfg. Returns a tuple (smld, info).

Pass-through SMLD semantic (NOT non-mutating copy)

Unlike cluster() — which deep-copies emitters so it can write cluster labels without touching the input — cluster_statistics() writes nothing onto the SMLD and returns the same reference as the input. The two-element tuple shape is preserved for ecosystem symmetry (every SMLM step returns (smld, info)), but callers should know the SMLD is the unmodified input, not a fresh copy. This asymmetry is intentional: copying for a read-only operation is wasted work, and statistic backends never have a label to write back.

Dispatch on AbstractStatisticsConfig alone has no implementation — concrete backends supply a method specialized on their own config type. Calling cluster_statistics with an unsupported config raises a clear error naming the available concrete backends.

See also: AbstractStatisticsConfig, ClusterStatisticsInfo, cluster (the labeling sibling).

AbstractStatisticsConfig <: SMLMData.AbstractSMLMConfig

Abstract supertype for cluster_statistics backend configurations.

Sibling to AbstractClusterConfig. Concrete subtypes (e.g. HopkinsConfig) configure spatial-statistic computations that do not modify the SMLD — they read the coordinate set, compute a scalar (and optional vector) summary, and return it alongside the input SMLD reference.

Shared fields (expected on every concrete subtype)

• use_3d::Bool: whether the z-coordinate participates in the statistic
• per_dataset::Bool: if true, compute per dataset and aggregate

Algorithm-specific fields (n_samples, seed, random_repeats, ...) live on the concrete subtype that needs them.

ClusterStatisticsInfo <: SMLMData.AbstractSMLMInfo

Secondary output from cluster_statistics() — summary of a spatial-statistic computation.

Fields

• n_locs_in::Int: number of input localizations
• statistic::Float64: primary scalar result (Hopkins H, median density, ...)
• statistic_name::Symbol: identifier for statistic (:hopkins, :median_density, ...)
• algorithm::Symbol: backend identifier (:hopkins, :voronoi_density, ...)
• elapsed_s::Float64: wall-clock time spent in the cluster_statistics call (seconds)
• extras::Dict{Symbol,Any}: per-backend supplementary outputs — vector-valued results (e.g. per-dataset Hopkins scores under :hopkins_per_dataset, per-emitter densities under :density_per_emitter) live here

Convention for vector-valued backends

Backends that produce a natural vector output (per-dataset Hopkins, per-emitter density, per-cluster silhouette) place the vector in extras under a descriptive key and a meaningful summary scalar (mean, median, ...) in statistic. This keeps the simple info.statistic access ergonomic for one-number consumers while preserving the full result for callers that need it.

Example

(_, info) = cluster_statistics(smld, HopkinsConfig(n_samples=50, seed=1))
println("Hopkins H = $(round(info.statistic, digits=3))")
per_ds = info.extras[:hopkins_per_dataset]  # Vector{Float64} when per_dataset=true

HopkinsConfig(; n_samples=20, random_repeats=1, seed=nothing,
                use_3d=false, per_dataset=true)

Configuration for the Hopkins-statistic spatial-randomness test.

Algorithm

For each repeat:

1. Draw n_samples reference points uniformly from the per-group bounding box; let u_i = NN distance from each reference point to the data.
2. Draw n_samples real points without replacement from the data; let w_i = NN distance from each sampled point to the OTHER real points (excluding itself).
3. With d = 2 (or 3 if use_3d=true), H = Σuᵢ^d / (Σuᵢ^d + Σwᵢ^d).

Repeats are averaged.

Fields

• n_samples::Int = 20: number of reference / sampled points per repeat. Must satisfy n_samples ≤ n_points per group.
• random_repeats::Int = 1: number of independent repeats to average. Higher values reduce variance at linear cost.
• seed::Union{Int,Nothing} = nothing: when set, seeds an internal Xoshiro for reproducibility. When nothing, uses the global RNG.
• use_3d::Bool = false: include the z-coordinate (and use d=3 in the formula).
• per_dataset::Bool = true: when true, compute Hopkins per dataset; the reported statistic is the mean across datasets and the full per-dataset vector is placed in extras[:hopkins_per_dataset]. When false, all emitters are pooled and a single H is returned.
• region = nothing: observation window for the uniform reference points (2D only). Hopkins is window-sensitive — sampling references over the data bounding box makes data that is uniform inside a non-convex boundary read as falsely clustered. Options: nothing → the data bounding box (default); a polygon Vector{NTuple{2,Float64}} → references are rejection-sampled inside it; :metadata → use the polygon at smld.metadata["edge_outer_polygon"] (written by classify_emitters — the pipeline channel); Dict(dataset_id => polygon) → one polygon per dataset (per_dataset = true). Incompatible with use_3d = true.

Interpretation

• H ≈ 0.5: data is consistent with uniform spatial randomness (Poisson)
• H → 1.0: strong clustering tendency
• H → 0.0: anti-clustering / regular spacing

Example

cfg = HopkinsConfig(n_samples = 50, random_repeats = 5, seed = 1)
(_, info) = cluster_statistics(smld, cfg)
println("Hopkins H = $(round(info.statistic, digits=3))")

See also: AbstractStatisticsConfig, ClusterStatisticsInfo, cluster_statistics.

VoronoiDensityConfig(; use_3d=false, per_dataset=true)

Configuration for the per-emitter Voronoi-density spatial-statistic backend.

Algorithm

1. For each group (dataset when per_dataset=true, all emitters otherwise), build the Voronoi tessellation of the emitter coordinates clipped to the convex hull (so every generator has a finite cell area).
2. For each emitter i, compute Aᵢ = its Voronoi cell area (μm²) and ρᵢ = 1/Aᵢ (μm⁻²).
3. Stitch the per-group vectors back into original emitter order so extras[:density_per_emitter][i] corresponds to smld.emitters[i].

Fields

• use_3d::Bool = false: must be false. 3D Voronoi tessellation is not supported (DelaunayTriangulation.jl is 2D only); passing true raises ArgumentError.
• per_dataset::Bool = true: when true, tessellate each dataset independently. When false, all emitters are pooled into one tessellation. Per-emitter outputs remain flat in original emitter order in either case.

Returned info

• statistic = median(non-NaN densities) — single number summarizing overall density. NaN if no group produced any valid density.
• statistic_name = :median_density
• algorithm = :voronoi_density
• extras[:density_per_emitter] — Vector{Float64} of length n_locs_in, in original emitter order. Emitters in groups smaller than 3 (untessellatable) receive NaN. Units: μm⁻².
• extras[:area_per_emitter] — Vector{Float64} of length n_locs_in, in original emitter order. Same NaN semantics. Units: μm².

Edge case handling (per V10 NaN-vs-throw rule)

• Groups with fewer than 3 points: those emitters receive NaN density and NaN area. Other groups proceed normally.
• Empty SMLD: empty per-emitter vectors, statistic = NaN.
• Group with exact-duplicate (x, y) coordinates: ArgumentError raised before triangulation (mirrors VoronoiConfig's guard — duplicate coordinates are a boundary-input issue, not a data-shape edge case).

Example

cfg = VoronoiDensityConfig()
(_, info) = cluster_statistics(smld, cfg)
ρ = info.extras[:density_per_emitter]   # Vector{Float64}, length == n_locs_in
A = info.extras[:area_per_emitter]
println("median density = $(round(info.statistic, digits=2)) μm⁻²")

See also: AbstractStatisticsConfig, ClusterStatisticsInfo, cluster_statistics, VoronoiConfig (the clustering sibling).

LocalContrastFeature(; density_k=200, background_k=2000,
                      use_3d=false, per_dataset=false)

Per-emitter local-density-contrast feature.

Algorithm

1. Build a KDTree over the (per-dataset or pooled) emitter coordinates.
2. For each emitter i, compute the kNN log-density f_i = log(density_k / (π · r_k²)) where r_k is the distance to the density_k-th nearest neighbor (excluding self). The kNN ball area normalization gives a density estimate in coordinate-units⁻²; the log keeps the downstream median/threshold comparisons in a comparable scale across a cell.
3. For each emitter i, compute the median of f_j over its background_k nearest neighbors (excluding self) — the local baseline at a coarser spatial scale. Median is used (not mean) because it tolerates a small number of locally-elevated neighbors without dragging the baseline up.
4. Contrast c_i = f_i − median(f_j over j in k_bg-NN of i). Positive contrast means point i is denser than its local surroundings.

Fields

• density_k::Int = 200: fine-scale neighborhood for the per-point log-density. Sets the spatial scale of the signal; for SMLM cells, ~80–200 nm typically.
• background_k::Int = 2000: coarse-scale neighborhood for the local baseline (must be > density_k). Sets the spatial scale of the baseline; should be large enough to span the structures you want to detect (typically 1–2 μm).
• use_3d::Bool = false: if true, build the KDTree over (x, y, z).
• per_dataset::Bool = false: if true, compute per dataset independently (each dataset gets its own KDTree). Per-emitter outputs stitch back into original emitter order in either case.

Returned info

• statistic = median(non-NaN contrast) — single scalar summary. NaN when no group has enough points to compute the feature.
• statistic_name = :median_local_contrast
• algorithm = :local_contrast
• extras[:contrast_per_emitter] — Vector{Float64} of length n_locs_in, in original emitter order. Units: nat (natural-log scale). The caller thresholds this directly.
• extras[:log_density_per_emitter] — Vector{Float64} of length n_locs_in, the fine kNN log-density f_i. Useful for absolute-density gates that complement the local-contrast gate (e.g. require f_i > q35 floor).

Edge case handling

• Groups with n ≤ density_k: those emitters receive NaN for both contrast and log-density. Avoids a degenerate kNN query against a too-small set.
• background_k ≥ n in a group: clamped to n − 1 for that group; the feature becomes "log-density minus group median," which is still well-defined.
• background_k ≤ density_k: raises ArgumentError at config use; the baseline must be coarser than the signal.
• Coincident coordinates (r_k = 0): the affected emitter receives NaN log-density and NaN contrast.

Composition

This feature is the missing primitive for hysteresis seed-and-grow on point clouds with non-stationary baseline density:

(_, info) = cluster_statistics(smld, LocalContrastFeature())
contrast = info.extras[:contrast_per_emitter]
fine = info.extras[:log_density_per_emitter]
fine_floor = quantile(filter(isfinite, fine), 0.35)
seed = isfinite.(contrast) .& (contrast .> 0.25) .& (fine .> fine_floor)
support = isfinite.(contrast) .& (contrast .> -0.05) .& (fine .> fine_floor)
(smld_out, _) = cluster(smld, PointHysteresisConfig(graph_k=12, min_points=150);
                        seed=seed, support=support)

See also: AbstractStatisticsConfig, ClusterStatisticsInfo.

classify_emitters(smld::BasicSMLD, cfg::AbstractEdgeClassifyConfig) -> (smld, info)
classify_emitters(x_um, y_um, cfg::AbstractEdgeClassifyConfig; fov_um) -> info

Classify each emitter as :outside, :membrane, or :interior. The concrete config type selects the cell gate by dispatch (OuterPolygonConfig, KdeValleyConfig); the rest of the pipeline is shared.

The SMLD method follows the package (out, Info) convention: it returns the smld with the published multi-cell mask threaded into smld.metadata["edge_cells"] (and the dominant cell's outer ring into smld.metadata["edge_outer_polygon"] for back-compat / Hopkins region=:metadata), plus an EdgeClassifyInfo. The per-emitter class lives only in info.class (with the in_cell / interior_mask accessors) — it is deliberately not mirrored into the metadata, because a per-emitter side-list desyncs the moment a downstream step subsets emitters. The coordinate method is the computational core and returns the info directly; fov_um = (xmin, xmax, ymin, ymax) in µm.

Config types for classify_emitters.

AbstractEdgeClassifyConfig is a sibling of the package's AbstractClusterConfig / AbstractStatisticsConfig (all <: SMLMData.AbstractSMLMConfig); each concrete config is dispatched as a method of classify_emitters:

• OuterPolygonConfig — multi-K density gate → multi-cell alpha-shape mask → point-in-region + membrane band.
• KdeValleyConfig — adaptive dSTORM density-valley gate (Gaussian-KDE + valley + footprint); gates on the original cloud, then builds the multi-cell mask on the footprint subset.

Struct fields are lowercase (idiomatic Julia); the UPPERCASE params.json keys are produced only at the serialization boundary by to_dict.

OuterPolygonConfig(; alpha_nm=300, membrane_nm=100,
                   fov_trunc_tol_nm=150, k_list=(16,128), rho_k_thresh=200)

Point-in-region vs the multi-cell alpha-shape mask on the multi-K-density-gated set, plus a membrane_nm band.

KdeValleyConfig(; alpha_nm=600, membrane_nm=100,
                fov_trunc_tol_nm=150, sigma_nm=150, ...)

Adaptive density-valley gate for dSTORM data. Gates on the per-FOV KDE density valley (threshold-free, no per-cell tuning). The defaults are tuned for dSTORM membrane data — notably alpha_nm = 600 (vs. the polygon default of 300) — so a bare KdeValleyConfig() is the intended entry point.

EdgeClassifyInfo{C} <: SMLMData.AbstractSMLMInfo

Result of classify_emitters. class is the authoritative per-emitter answer (:outside, :membrane, :interior), a partition of the input set.

cells::MultiCellMask is the published mask: one CellPolygon (an outer ring plus optional internal holes) per distinct cell in the FOV, ordered largest-first. It is the drawn boundary and the Hopkins region = :metadata observation window. outer_polygon is cells[1].outer (the dominant cell's outer ring), retained for back-compat.

inside_outer is geometric membership in the mask (in_region(cells)), and dist_to_outer_um is the distance to the nearest real boundary segment (NaN when not inside). Boundary segments lying on a truncated FOV edge are excluded from that distance, so a field-of-view cut is never labeled membrane — membrane is the band within membrane_nm of a real cell edge. in_cell(info) (class .!= :outside) equals inside_outer. Read class for the answer.

loops are the raw alpha-shape boundary loops (serialized to polygon_loops.tsv); config holds the concrete config that ran (honest provenance).

in_cell(info) -> BitVector

Topological cell membership, == (info.class .!= :outside). Equals inside_outer (both are mask membership via in_region).

interior_mask(info) -> BitVector

Per-emitter interior mask, == (info.class .== :interior) — the strictly-interior emitters (membrane and outside both excluded). Unlike in_cell (which keeps membrane), this is the subset a downstream analysis usually carries forward, and is the boolean to AND with any other per-emitter mask (e.g. a separate structure mask) before subsetting. classify_emitters stays non-destructive — it never drops emitters — so the caller composes masks and subsets once.

interior_fraction(info) -> Float64

CellPolygon(outer, holes = [])

One cell's footprint: an outer boundary ring (a Vector of (x, y) µm tuples, stored open / not repeating the first vertex) plus optional internal holes (each a ring). holes is empty unless the mask was built with keep_internal = true.

MultiCellMask

Vector{CellPolygon} — the edge mask of a field of view: one CellPolygon per distinct cell. Returned by classify_emitters (info.cells / metadata["edge_cells"]) and accepted as a Hopkins observation window. Cells are ordered largest-first, so mask[1] is the dominant cell.

build_mask(loops; keep_internal = false, min_cell_frac = 1/3, min_hole_frac = 0) -> MultiCellMask

Assemble a MultiCellMask from raw alpha-shape boundary loops (closed vertex rings):

1. Split each loop into simple sub-rings at self-touch / shared vertices.
2. Nest by depth (how many other rings contain it): even depth ⇒ a cell outer, odd depth ⇒ an internal region (a real void or a pinch-enclosed region).
3. Group each cell outer with the internal regions directly inside it — only when keep_internal = true; otherwise cells are solid (holes empty). With min_hole_frac > 0, holes smaller than min_hole_frac × (that cell's outer area) are dropped (filled back into the interior), so only real voids above that scale survive — sub-cell texture such as inter-cluster gaps no longer fragments the cell.
4. Drop debris: remove any cell whose outer area is below min_cell_frac × (largest cell area); min_cell_frac = 0 keeps everything.

Cells are returned largest-first.

in_region(x, y, mask::MultiCellMask) -> Bool

True when (x, y) lies in any cell of mask: inside a cell's outer ring and not inside one of that cell's holes.

region_area(mask::MultiCellMask) -> Float64

Total occupied area in µm²: Σ cell outer areas − Σ hole areas.

EdgeReport

Figure-data derivative of an EdgeClassifyInfo: the classified SMLD plus per-emitter coordinates/classes and per-class fractions needed to write diagnostics and render the standard edge figures. Produced by compute_edge_report; consumed by write_edge_report and plot_edge_report.

compute_edge_report(smld, info::EdgeClassifyInfo) -> EdgeReport

Bundle the classified SMLD and info into the figure-data report. Pass the SMLD returned by classify_emitters — its emitters are unchanged, so info.class is 1:1 with smld.emitters.

write_edge_report(report; output_dir, condition="", cell="") -> String

Write the edge-classify text/TSV diagnostics (classified.tsv, polygonloops.tsv, loopdiagnostics.csv, params.json, manifest.json) for report into output_dir (folds write_edge_artifacts). Returns output_dir.

class_codes(info::EdgeClassifyInfo) -> Vector{Int}

Per-emitter integer class code for categorical rendering: outside = 0, membrane = 1, interior = 2. The 0 for outside matches SMLMRender's reserved-gray id.

method_name(cfg) -> String

Short identifier for the edge-classification strategy selected by cfg's concrete type: "outer_polygon" for OuterPolygonConfig and "kde_valley" for KdeValleyConfig. Used in diagnostics, logging, and output manifests.

write_edge_artifacts(leaf, info::EdgeClassifyInfo, x_um, y_um; condition, cell,
                     smld_input_meta=nothing)

Write the diagnostic artifact set (classified.tsv, polygon_loops.tsv, loop_diagnostics.csv, params.json, manifest.json) into leaf. x_um/y_um are the original coordinates (the info does not retain them). Compute and IO are separate: classify_emitters does not write artifacts.

compute_concavity_metric(info::EdgeClassifyInfo, x_um, y_um;
                         buffer_um = 2.0, asym_R_nm = 1000.0, asym_gate = 0.20,
                         rho_lo = 200.0, intracellular_dense_threshold = 0.8)
    -> ConcavityMetricReport

Boundary-proximal concavity-error metric for the outer-polygon classifier. Identifies emitters classified :interior that live in deep concave bays the alpha-shape bridged across. Measured against the dominant cell (cells[1].outer); on a multi-cell FOV it characterizes that cell only.

A :interior emitter is a suspect iff: (1) it lies within buffer_um of the outer polygon boundary; (2) it is not inside an "intracellular void" loop (non-outer loop with frac_in_fov == 1 and frac_dense >= intracellular_dense_threshold); (3) its directional asymmetry at asym_R_nm is >= asym_gate; (4) its local density ρK(K=128) is `<= rholo`.

buffer_um is a diagnostic parameter (default 2.0 µm), not part of the classifier config.

ConcavityMetricReport

Boundary-proximal concavity-error metric for the outer-polygon classifier (see compute_concavity_metric).

LoopDiagnostic

Per-loop diagnostic record (one row of loop_diagnostics.csv).