density-clustering

Header-only C++20 density-based clustering. A small, dependency-free toolkit of the most useful density-based clustering algorithms: find clusters of arbitrary shape, at differing densities, with noise — without knowing the number of clusters in advance. Built for time-critical use on large clouds (millions of points).

It offers four algorithms behind one consistent API:

• OPTICS (Ankerst et al., 1999) — the namesake. Produces a reachability ordering (a reachability plot) you cut at one or many density scales. Background: paper, Wikipedia, YouTube.
• HDBSCAN* (Campello et al., 2013) — a parameter-light cluster hierarchy condensed by stability; no epsilon, just min_cluster_size. Background: John Healy's talk HDBSCAN, Fast Density Based Clustering, the How and the Why.
• sOPTICS / sHDBSCAN — scalable, approximate cosine variants via random projections (sDBSCAN/sOPTICS, Xu & Pham 2024) for the high-dimensional or very-large-n regime.

The only vendored dependency is nanoflann; everything else is the C++ standard library. (The C++ namespace is optics::, and the public header is <optics/optics.hpp> — a historical name from when the library began as an OPTICS-only implementation; it is kept stable for compatibility.)

When to use density-based clustering (vs k-means & DBSCAN)

These methods do not need the number of clusters up front, they find arbitrarily-shaped clusters, and they label low-density points as noise. OPTICS and HDBSCAN* go one step beyond DBSCAN: they do not commit to a single global density, so they recover clusters that sit at different densities — the case a single DBSCAN eps cannot capture.

Generated by python tools/compare_algorithms.py — OPTICS labels come from this library; k-means and DBSCAN from scikit-learn. On moons, k-means cuts across the crescents because it assumes convex, comparably-sized blobs, while OPTICS and DBSCAN recover the true shapes. On varied, k-means has no notion of noise and splits the elongated cluster, while the density methods isolate the structure and leave the sparse background as noise. On density — two tight clusters beside one sparse cluster — no single DBSCAN eps works (the radius that connects the sparse cluster also merges the dense pair, so DBSCAN collapses to two), and k-means ignores density entirely; only OPTICS, reading the reachability hierarchy with the ξ method, recovers all three. This is the case density-based clustering is built for.

	needs k?	arbitrary shapes	noise	varying density	hierarchy	speed	deterministic
k-means	yes	no	no	no	no	very fast	no (random init)
DBSCAN	no	yes	yes	no (one global eps)	no	fast	yes
OPTICS	no	yes	yes	yes (ordering across scales)	yes (ξ)	slower (query-bound)	yes
HDBSCAN*	no	yes	yes	yes (stability hierarchy)	yes (condensed tree)	slower (query-bound)	yes

Be honest about the trade-offs. If your clusters are roughly convex and you know k, k-means is far faster and perfectly adequate. If a single density threshold separates your clusters, DBSCAN is simpler and quick. Reach for OPTICS or HDBSCAN* when you don't know k, when clusters sit at different densities (where one DBSCAN eps can't win), or when you want to see the cluster structure before committing to a cut. These are the most general of the methods — and the most expensive: runtime is dominated by neighbor queries and cannot match k-means. That said, we are much faster than other OPTICS implementations — even faster than scikit-learn's DBSCAN — with scalable approximate variants (sOPTICS / sHDBSCAN) for the hardest high-dimensional / very-large-n regimes. → Quantitative runtime comparisons (this library's backends vs scikit-learn OPTICS, DBSCAN, and k-means, across sample sizes and dimensions) are in perf/README.md.

The algorithms — which one, and when

The four split along two axes — what structure they build (a 1-D reachability ordering vs. a cluster hierarchy) and how they find neighbors (an exact KD-tree search vs. approximate CEOs random projections):

	output	neighbors	metric	exact?	reach for it when
OPTICS	reachability ordering	exact KD-tree	Euclidean	yes	the default — low/medium-D, you want the exact ordering
HDBSCAN*	cluster hierarchy → labels	exact k-NN	Euclidean	yes	you want a parameter-light hierarchy (min_cluster_size only, no epsilon) that tolerates varying density
sOPTICS	reachability ordering	CEOs random projections	cosine (L2/L1 via kernel features)	approx	high dimension (≳ 16-D) or very large n, and your similarity is angular
sHDBSCAN	cluster hierarchy → labels	CEOs random projections	cosine	approx	HDBSCAN* clustering at n ≳ 1e5, where the exact O(n²) MST is infeasible

Two facts worth knowing up front, both explained in the secondary doc:

• Why the "s" variants are cosine. CEOs finds neighbors by ranking inner products (random projections preserve inner-product order), which is intrinsically angular — so points are L2-normalized onto the unit sphere where Euclidean distance is monotone in cosine distance. L2/L1 are supported by first embedding into random kernel features, then running the same cosine pipeline.
• Why the "s" variants are slower on small data. They pay a large fixed O(n·D·Dim) projection cost up front. That only pays off once dimension and/or n make the exact neighbor search expensive; below the crossover, exact OPTICS/HDBSCAN* win.

→ Full comparison — the two design seams, the cosine derivation, the small-n crossover, the HDBSCAN* MST backbones, and the efficiency roadmap — in docs/algorithms.md.

Quickstart: cluster your own data

C++ integration

This library is header-only. Simply add include/ to your include path and add #include <optics/optics.hpp>. Then, one call of optics::extract_xi() is all you need - see below for more details.

Python integration

Build + usage in python/README.md.

Testing on CSV data & visualization

You don't have to write any C++ to try it — build the bundled example and point it at a point cloud in CSV format (one numeric row per point, with a x0,x1,... header; extra trailing columns are ignored).

# 1. Build the example (MSVC shown; use linux-gcc / linux-clang on Linux/macOS)
cmake --preset msvc
cmake --build --preset msvc --target cluster_csv

# 2. Use our tool datasets.py to generate the 'moons' dataset (see example above) in csv format.
#    Alternatively: Bring your own point cloud in csv format and use it in step 3. 
python tools/datasets.py --name moons --n 1500 --out data/moons.csv

# 3. Cluster it  ->  writes <out>_points.csv (labels) and <out>_reach.csv (the plot)
build/examples/Release/cluster_csv data/moons.csv data/moons 10        # 10 = min_pts

# 4. See it
python tools/visualize.py --points data/moons_points.csv --reach data/moons_reach.csv

On Linux/macOS the binary is build/examples/cluster_csv. Python deps for the tools: pip install -r requirements.txt. See examples/cluster_csv/README.md for all options, and Reading the reachability plot / Choosing parameters below.

Regular usage from a C++ application

A point is a std::array<T, Dim> of Cartesian coordinates (T is float or double); a cloud is a std::vector of those.

#include <optics/optics.hpp>
#include <array>
#include <vector>

using Point = std::array<double, 2>;

int main() {
    std::vector<Point> points = /* your data */;

    // The OPTICS cluster-ordering + reachability distances.
    auto reach = optics::compute_reachability_dists(points, /*min_pts=*/10);

    // Flat clusters by a reachability threshold ...
    auto clusters = optics::get_cluster_indices(reach, /*threshold=*/2.0);

    // ... or the hierarchical Xi method (nested cluster trees).
    auto trees = optics::get_chi_clusters(reach, /*chi=*/0.05, /*min_pts=*/10);
}

If you need more dials and knobs: The full signature of the compute_reachability_dists lets you choose the backend, neighbor-acquisition strategy, and thread count:

optics::compute_reachability_dists<T, Dim, Backend>(
    points, min_pts,
    epsilon  = -1.0,                       // auto-estimated when <= 0
    mode     = optics::NeighborMode::OnDemand,  // default; Precompute is the parallel, opt-in cache
    n_threads = 0);                        // threads for Precompute only (0 => hardware concurrency)

Convenience helpers

Two one-call helpers, both (points, min_pts, [param]) → one point-index list per cluster:

• optics::cluster_threshold(points, min_pts, threshold = auto) — compute the ordering and cut the reachability plot at a flat threshold. This is the paper's ExtractDBSCAN: the same clustering DBSCAN gives at eps = threshold — we do not run DBSCAN. Omit threshold for an educated default (a high percentile of the reachabilities; inspect the plot to tune). (Was cluster_dbscan, now a deprecated alias.)
• optics::extract_xi(points, min_pts, chi = 0.05) — hierarchical ξ (steep-area) extraction, flattened to a list of clusters. For the nested cluster tree, use optics::get_chi_clusters(reach, chi, min_pts): each node holds a (begin,end) range into the ordering, and the cluster's points are reach[i].point_index for i ∈ [begin,end).
• optics::convert_cloud<float>(int_points) — convert an integer/byte cloud (e.g. uint8 color data) to a floating-point cloud, since T must be float/double.

HDBSCAN* (a separate clusterer)

If you'd rather not pick an epsilon or read a reachability plot, optics::hdbscan builds a cluster hierarchy and condenses it by stability — you supply only min_cluster_size:

#include <optics/hdbscan.hpp>

auto result = optics::hdbscan(points, /*min_cluster_size=*/15);
// result.labels[i] == -1 for noise, else the cluster id; result.probabilities[i] is the membership strength.

It scales the same way the suite does: the exact O(n²) MST is the default up to ~1e4 points; pass optics::MstAlgorithm::Auto (or Boruvka / KnnGraph) for the sub-quadratic backbones at larger n, and optics::shdbscan(...) for the approximate cosine variant in the high-dimensional / 1e5+ regime. See docs/algorithms.md.

Python (optional binding)

An optional pybind11 binding exposes OPTICS for 1/2/3/4-D NumPy clouds (off by default; the C++ library stays dependency-free):

import numpy as np, optics_py
labels = optics_py.cluster_threshold(pts, min_pts=10, threshold=2.0)   # (N, Dim) -> per-point labels
labels_xi = optics_py.extract_xi(pts, min_pts=10, chi=0.05)

Build + usage in python/README.md. For data already on disk, the cluster_csv example + tools/visualize.py need no binding at all.

Visualizing results

The core writes no images; export CSV and render with the bundled script (matplotlib). After processing a point cloud points into reachability distances reach and labels as shown above, you can export these via the optics::io:: functions like this:

#include <optics/io.hpp>
auto labels = optics::io::cluster_labels(points.size(), clusters);
optics::io::export_points_csv("points.csv", points, labels);   // any dimension
optics::io::export_reachability_csv("reach.csv", reach);

And have them plotted through the shipped python scripts:

python tools/visualize.py --points points.csv --reach reach.csv --out plot.png

visualize.py handles 2D and 3D scatter (e.g. color spaces) and falls back to a PCA projection for higher dimensions.

Reading the reachability plot

OPTICS doesn't hand you clusters directly — it produces a reachability plot: the points laid out in cluster-order, each bar its reachability distance. Read it like a landscape.

• Valleys (runs of low bars) are clusters — points packed densely together.
• Peaks (tall bars) are the jumps between clusters, or sparse/noise points.
• A horizontal cut at some height is exactly the flat get_cluster_indices(reach, threshold) extraction: every valley dipping below the line becomes a cluster. The hierarchical ξ method (get_chi_clusters / extract_xi) instead follows the valley walls, so it can pull out clusters that sit at different depths — the case a single horizontal cut (or a single DBSCAN eps) cannot capture.

Choosing parameters

• min_pts (required) — how many neighbors make a point "core". Higher values smooth the plot and ignore small/noisy groups; lower values are more sensitive. A common starting point is 2 × dim to ~20; raise it if the plot is jagged, lower it if real small clusters get swallowed.
• epsilon (optional, auto-estimated when ≤ 0) — the largest neighborhood radius considered. It mainly bounds cost and memory: large enough captures all structure; too small truncates the plot (UNDEFINED reachability, shown as full-height bars). Leave it auto unless you need to cap work on huge clouds.
• threshold (flat cut) — the reachability height at which valleys become clusters. Pick it by looking at the plot: just under the valley floors you care about. Smaller → tighter, more clusters; larger → fewer, looser clusters.
• chi (ξ method) — relative steepness (e.g. 0.05) that delimits a cluster's walls. Use ξ instead of a flat threshold when clusters live at different densities. It outputs a hierarchical ordering (a tree) of clusters that can be nested within one another (via get_chi_clusters).

For domain-specific tuning and gotchas on image/color data (flat-color regions, eps in RGB units, separating anti-aliasing/JPEG "bridge" colors), see examples/color_clustering/README.md.

Performance

OPTICS comfortably handles millions of low-dimensional points (a few seconds for 1e6 3-D points; see perf/README.md for hardware and scaling). Computation time is dominated by neighbor queries, so it is query-bound and high dimensionality is expensive (the curse of dimensionality). On dense data (e.g. flat-color images) neighborhoods grow with n and the ordering tends toward O(n²) in both time and memory — keep epsilon modest, use OnDemand mode, or downsample.

In the 1.0.0 reference benchmark (a reproducible matrix over n = 316…1e6, d = 1…128, several layouts — see docs/benchmarking.md), on identical clouds this library's OPTICS is a median ~109× faster than scikit-learn's OPTICS (range ~6×–700×), and runs in regimes (n ≥ 1e5) where scikit-learn OPTICS is too slow to finish at all. Our HDBSCAN* matches scikit-learn's HDBSCAN clustering essentially exactly (mean ARI difference 0.0002). k-means (no neighbor graph) remains the cheapest per run.

Tuning knobs — and the matrix confirms the defaults are the right ones for the common case: OnDemand (the default — lean memory, never OOMs) vs Precompute (opt-in parallel cache, faster wherever it fits in RAM but O(n × neighbors) memory); the thread count (Precompute only); the knee ε estimator (default; best where cluster densities vary); and the backend — exact NanoflannBackend (default; fastest up to ~16-D and never loses quality), ApproxNanoflannBackend (a modest win only in high dimensions, where search dominates), or the optional Boost R*-tree / HNSW. For high-dimensional or very-large-n data, switch to the approximate sOPTICS / sHDBSCAN (cosine) — they overtake the exact methods exactly in that regime.

→ Full performance analysis — the 1.0.0 reference matrix and decisions in docs/benchmarking.md; committed baselines, scaling by sample size, per-backend and scikit-learn/DBSCAN/k-means comparisons (incl. real images), and why the approximate backend helps only in high dimensions in perf/README.md.

Documentation

• docs/algorithms.md — the four algorithms compared (OPTICS / sOPTICS / HDBSCAN* / sHDBSCAN): when to use which, why the statistical variants are cosine, why they cost more on small data, and the efficiency roadmap.
• docs/benchmarking.md — the 1.0.0 reference benchmark matrix (speed/quality vs scikit-learn, decisions D1–D5), datasets, and the fairness rules for a comparison.
• perf/README.md — performance: benchmarks, scaling, backend & cross-library comparisons, the approximate-backend analysis.
• examples/cluster_csv/README.md — cluster your own CSV (2/3/4/16-D), all options.
• examples/color_clustering/README.md — the color-space guide: image pipeline, output modes, color parameter tips, and gotchas.
• python/README.md — the optional NumPy (pybind11) binding.
• tools/README.md — visualization, dataset generators, and comparison/validation scripts.
• docs/API-STABILITY.md — the frozen 1.0.0 public-API contract: what is stable, what is internal, and the semver rules.
• docs/ROADMAP-0.9.1.md — the current milestone plan; docs/ROADMAP-post-0.9.1.md — the lookahead toward 1.0.0.

Stability & versioning

From 1.0.0 the public optics:: API is frozen and follows Semantic Versioning: signatures, default arguments, struct layouts, and enum values in the stable surface won't break without a major-version bump. Everything in optics::detail::, the vendored nanoflann/hnswlib, and the test/benchmark helpers (optics::testdata::, stopwatch::) is internal and may change in any release. The approximate algorithms (compute_soptics_reachability_dists, shdbscan) keep a stable signature and seed-determinism, but their exact labels may improve across minor versions — depend on cluster quality, not specific label values. Full contract: docs/API-STABILITY.md.

Dependencies

None required: nanoflann (BSD 2-Clause) is vendored under include/optics/, and everything else is the C++ standard library. Boost is an optional alternative neighbor-search backend, enabled with -DOPTICS_ENABLE_BOOST_RTREE=ON. Bundled third-party licenses are listed in THIRD-PARTY-LICENSES.md.

Building & testing

Header-only — just add include/ to your include path and #include <optics/optics.hpp>. To build the tests/benchmark (CMake ≥ 3.21, MSVC 2022 / GCC 10+ / Clang 13+):

cmake --preset linux-gcc      # or: linux-clang, msvc
cmake --build --preset linux-gcc
ctest --preset linux-gcc

License

Distributed under the MIT Software License (X11 license). (See accompanying file LICENSE.) Bundled third-party components and their licenses are listed in THIRD-PARTY-LICENSES.md.