spline-analyzer-torch

PyTorch reimplementation of a PhD dissertation algorithm for unsupervised aerial image segmentation and anomaly detection.
Deterministic, device-agnostic, fully tested against the original NumPy research codebase.

---

Highlights

Task	Unsupervised pixel-level segmentation of UAV/aerial RGB images
Density model	Non-parametric S₂,₀ biquadratic spline over (mean, σ) feature space
Segmentation	GA (gradient-ascent convergence potentials) · DL (valley-scan + TLS discriminant lines)
Anomaly detection	Quantile gates + optional GMM post-selection
Validation	Strict numerical golden tests against original numpy implementation
Device	CPU · CUDA · MPS (device-agnostic API)
Python	≥ 3.10 · PyTorch ≥ 2.2

---

Scientific Background

This project is a clean-room PyTorch port of the core algorithm from a dissertation on non-parametric texture-based segmentation of drone imagery.

Pipeline:

RGB image
  └─ FeatureExtractor ──────────→ (mean_map, sigma_map)   [local window stats]
       └─ SplineS20.fit ──────→ S₂,₀ density surface     [histogram → 9 coeff maps]
            └─ find_modes ────→ local maxima               [Newton-refined, Hessian-verified]
                 └─ DLClusterer / GAClusterer ─→ label map [valley-scan TLS / gradient-ascent]
                      └─ QuantileAnomalyDetector ─→ mask   [quantile gates + GMM]

Feature parameterization:

• Intensity: I = (R + G + B) / 3
• Local mean μ and std σ via sliding window (OpenCV boxFilter, BORDER_REFLECT)

S₂,₀ spline:

• 9 coefficient maps C₀₀…C₂₂ per histogram cell
• Analytical gradient and Hessian (used for mode detection and GA clustering)
• Float64 coefficients, float32 evaluation

---

Project Layout

spline-analyzer-torch/
├── src/spline_torch/
│   ├── __init__.py          # Public API: AnalysisConfig, SplineAnalyzer
│   ├── config.py            # AnalysisConfig dataclass
│   ├── features.py          # FeatureExtractor (OpenCV exact-parity backend)
│   ├── spline.py            # SplineS20: fit / evaluate / gradient / hessian / find_modes
│   ├── clustering.py        # GAClusterer, DLClusterer (valley-scan + TLS)
│   ├── anomaly.py           # QuantileAnomalyDetector, apply_post_selection_gmm
│   ├── pipeline.py          # SplineAnalyzer orchestration API
│   └── _geometry.py         # Rhombus geometry, TLS line fit
├── tests/
│   ├── conftest.py
│   ├── test_features.py
│   ├── test_spline.py
│   ├── test_clustering.py
│   ├── test_anomaly.py
│   ├── test_pipeline_golden.py
│   ├── test_golden_real_scenes.py   # Strict golden vs original implementation
│   ├── test_config.py
│   ├── test_properties.py
│   └── golden/
│       ├── generate_reference_golden.py
│       └── generate_real_scene_goldens.py
├── benchmarks/
│   ├── bench_features.py
│   ├── bench_spline.py
│   └── gpu_vs_cpu_table.py
├── notebooks/
│   ├── 01_quickstart.ipynb   # Minimal CPU smoke-run on synthetic input
│   └── 02_aerial_demo.ipynb  # Aerial scene walkthrough
├── pyproject.toml
├── LICENSE                   # MIT
└── README.md

---

Installation

# Basic install
pip install -e .

# With dev dependencies (pytest, ruff, black, hypothesis)
pip install -e .[dev]

---

Quick Start

import torch
from spline_torch import AnalysisConfig, SplineAnalyzer

config = AnalysisConfig(
    window_radius=8,          # local statistics window radius
    bins_mean=64,             # spline histogram grid
    bins_sigma=64,
    clustering_method="dl",   # "dl" (valley-scan) or "ga" (gradient-ascent)
    p_auto_confidence=0.0,    # GMM post-selection off
)

analyzer = SplineAnalyzer(config=config, device="cpu")  # or "cuda", "mps"

# image: (3, H, W) RGB float32 tensor, values in [0, 255]
image = torch.rand(3, 720, 1280) * 255

result = analyzer.run(image)
print(result.labels.shape)        # (H, W) int64 cluster map
print(result.anomaly_mask.shape)  # (H, W) bool
print(result.n_clusters)          # number of detected modes
print(f"anomaly rate: {result.anomaly_rate:.3%}")

Stage-by-stage API

features     = analyzer.extract_features(image)
spline_state = analyzer.fit_density(features)
clusters     = analyzer.segment(features, spline_state)
anomalies    = analyzer.detect_anomalies(features, clusters.labels)

---

Examples

Scenes from the public LoveDA dataset (Wang et al., 2021 — see Running on Public Datasets below), run through SplineAnalyzer.run() with the default AnalysisConfig (window_radius=8, bins_mean=64, bins_sigma=64, clustering_method="dl", p_auto_confidence=0.5). LoveDA is released for research use, so these scenes are shown directly.

Source	Input	Segmentation	Anomaly mask
LoveDA U002 (urban)
LoveDA R002 (rural)
LoveDA U004 (urban)

Cluster colors are arbitrary per-run IDs, not semantic land-cover classes. The anomaly mask (red overlay) is the raw output of the default quantile gate — at default settings it flags roughly half of all pixels as anomalous on these scenes (50.0–51.0%), which is a property of running the 5%/95% quantile thresholds unmodified on real (non-synthetic) imagery, not a rendering error.

The real UAV scenes referenced in Testing (DSC_1.jpg, DSC_5.jpg) are a separate, non-redistributable source used only for strict golden-test validation against the original dissertation code — they are not shown here and are not included in this repository.

---

Benchmark Results

Evaluated on 150 real aerial scenes from three independent sources (UAV footage, LoveDA, OpenEarthMap); 1050 measurements total (150 scenes × 7 algorithms). Ranked by ARI:

Method	ARI	NMI	mIoU	mDice	Silhouette
GA + p_auto=0.5 (this method)	0.134	0.151	0.409	0.267	0.635
DL + p_auto=0.5 (this method)	0.121	0.136	0.402	0.255	0.611
GMM (Oracle K)	0.117	0.134	0.248	0.300	0.395
Agglomerative (Oracle K)	0.115	0.122	0.244	0.288	0.447
K-Means (Oracle K)	0.114	0.123	0.244	0.292	0.457

The proposed method ranks first among seven compared algorithms on ARI, NMI, mIoU and Silhouette, while inferring the number of clusters automatically; the classical baselines are given the true cluster count (Oracle K). It does not lead on mDice (0.267 vs 0.300 for GMM). The mIoU margin is the largest: 0.409 vs 0.244 for K-Means on the pooled set; on the UAV subset the method exceeds K-Means on mIoU by 68.9%.

Effect of post-selection (p_auto)

Metric	Without	With p_auto=0.5	Relative gain
ARI	0.100	0.134	+33.3%
NMI	0.116	0.151	+30.4%
Silhouette	0.405	0.635	+57.0%

Robustness and determinism

• Under additive Gaussian noise from 0% to 50%, ARI degrades by less than 0.2%; GMM degrades by 24.3% on the same data.
• Identical results across 10 independent runs (σ=0); GMM: σ=0.083.
• ≈0.4 s per 3,000-point patch.

Anomaly detection

Spline Tail Probability vs specialised detectors (synthetic data):

Method	ROC-AUC	F1
Isolation Forest	0.940	0.736
Spline Tail Probability	0.932	0.635
LOF	0.840	0.696
One-Class SVM	0.745	0.535

The spline detector is competitive with, not superior to, Isolation Forest. For fine-grained anomaly classification a hybrid pipeline is recommended: spline method for structural scene segmentation, Isolation Forest for anomaly classification. The method also tends to under-segment, slightly underestimating cluster count.

---

Configuration

Parameter	Default	Range	Description
window_radius	8	[1, 64]	Local statistics window radius (pixels)
stride	1	[1, 8]	Feature map downsampling stride
bins_mean	64	[16, 256]	Histogram bins — mean axis
bins_sigma	64	[16, 256]	Histogram bins — sigma axis
min_mode_distance	6	[1, 20]	Minimum mode separation (grid units)
mode_threshold_rel	0.05	[0.0, 1.0]	Mode detection relative threshold
lower_quantile	0.05	(0, 1)	Lower anomaly quantile
upper_quantile	0.95	(0, 1)	Upper anomaly quantile
sigma_gate_enabled	True	bool	Robust sigma gate for anomaly stabilization
sigma_gate_mad_scale	1.0	[0.0, 10.0]	MAD scale for sigma gate
p_auto_confidence	0.50	[0.0, 1.0]	GMM post-selection confidence
clustering_method	"dl"	"ga", "dl"	Segmentation strategy
eval_dtype	"float32"	"float32", "float64"	Runtime evaluation dtype

---

Numerical Precision Policy

Stage	Dtype	Location	Rationale
Histogram build	float64	CPU/NumPy	Deterministic accumulation
Coefficient maps C₀₀…C₂₂	float64	CPU/NumPy	Maximum coefficient accuracy
evaluate() / gradient()	float64 input → float32 output	PyTorch	Scientific fidelity with deployment efficiency
Feature extraction	float32	CPU/NumPy	OpenCV-native precision

Strict numerical equivalence with the original NumPy research code is validated by golden tests on real UAV scenes.

---

Testing

pytest                      # 39 tests (unit + golden + property)
pytest -q --tb=short        # compact output
pytest tests/test_golden_real_scenes.py   # strict real-scene golden comparison

Golden Workflow

The project ships with two golden test tiers:

1. Self-consistency golden (no original code required):

python tests/golden/generate_reference_golden.py
pytest tests/test_pipeline_golden.py

2. Strict real-scene golden (requires original dissertation codebase in ../src/):

python tests/golden/generate_real_scene_goldens.py
pytest tests/test_golden_real_scenes.py

Validated with real UAV scenes (DSC_1.jpg, DSC_5.jpg):

Quantity	Tolerance	Status
mean_map, sigma_map	atol 1e-4	✓ strict
histogram	atol 1e-10	✓ strict
C₀₀…C₂₂ coefficient maps	atol 1e-10	✓ strict
evaluate() probe grid	atol 1e-10	✓ strict
gradient() probe grid	atol 1e-10	✓ strict
mode positions	atol 1e-6	✓ strict
label map (pixel agreement)	≥ 98 %	✓ structural

These specific source scenes are not redistributed with this repository. test_golden_real_scenes.py skips automatically when the referenced images/goldens are absent, so the rest of the suite runs standalone.

Running on Public Datasets

The pipeline itself is dataset-agnostic — it accepts any (3, H, W) RGB tensor. To try it on real aerial/UAV imagery without private data, use one of these public, freely available datasets:

• LoveDA — Wang, J., Zheng, Z., Ma, A., Lu, X., & Zhong, Y. (2021). "LoveDA: A Remote Sensing Land-Cover Dataset for Domain Adaptive Semantic Segmentation." NeurIPS 2021 Track on Datasets and Benchmarks. arXiv:2110.08733 · Zenodo DOI: 10.5281/zenodo.5706578
• OpenEarthMap — Xia, J., Yokoya, N., Adriano, B., & Broni-Bediako, C. (2023). "OpenEarthMap: A Benchmark Dataset for Global High-Resolution Land Cover Mapping." WACV 2023, 6254–6264. arXiv:2210.10732

Download a scene from either source, load it as an RGB tensor, and pass it to SplineAnalyzer.run() as shown in Quick Start — no other setup required.

---

Benchmarks

python benchmarks/bench_features.py --device cpu --height 720 --width 1280
python benchmarks/bench_spline.py   --device cpu --height 720 --width 1280
python benchmarks/gpu_vs_cpu_table.py --height 720 --width 1280 --out artifacts/benchmark_table.md

Measured on CPU, shape 720×1280:

Stage	Fit (ms)	Eval (ms)	Gradient (ms)	Features (ms)
pipeline-core	61.8	123.2	134.4	33.3

CUDA/MPS unavailable in the test environment; accelerator columns omitted.

---

Key Design Decisions

• No JAX in the main pipeline. Optional jax extra exists only for research side-experiments.
• OpenCV for feature extraction. The exact BORDER_REFLECT padding of cv2.boxFilter is required to match the dissertation's numerical outputs. A pure-PyTorch fallback is retained but not guaranteed bit-exact.
• sklearn GMM retained. apply_post_selection_gmm wraps sklearn.GaussianMixture deliberately — it is a post-processing step, not a core tensor operation.
• Histogram stays on CPU float64. Moving it to GPU/float32 would break strict numerical reproducibility with the research code.
• Deterministic mode. torch.use_deterministic_algorithms(True, warn_only=True) is set on SplineAnalyzer construction.

---

Future Work

• Optional CUDA kernel for spline coefficient lookup (hot path for GPU workloads).
• Cross-dataset validation benchmarks (LoveDA, custom UAV sets).
• Deployment packaging examples (ONNX export, TorchScript).

---

Intended Use

This repository provides a research method for unsupervised image segmentation and anomaly detection on high-resolution imagery, developed for aerial and remote-sensing data. Applicable domains include remote sensing, aerial-surveillance imagery analysis, precision agriculture, environmental and land-cover monitoring, and industrial quality inspection.

No real-world classified imagery is included in this repository or its test fixtures.

---

License

MIT — see LICENSE.

Status

Core pipeline stages are wired and executable:

• features
• spline density
• segmentation (GA + baseline DL)
• anomaly detection
• end-to-end pipeline API

Planned upgrades (golden equivalence depth, full valley-TLS DL, richer benchmarks) are tracked as incremental work.

Measured benchmark numbers above are from a 2026-07-11 CPU run at shape 720×1280; CUDA/MPS were unavailable in that environment.