Skip to content

Latest commit

 

History

14 Commits

Folders and files

NameName
Last commit message
Last commit date
 
 
 
 
 
 
 
 
 
 
 
 
 
 

Repository files navigation

clearmap-ng

This crate is not ready for use. Some new name will be picked in the future

This is a rewrite of ClearMap, but

  • It aims to reproduce some of the results for now, primarily as a test that it fulfills the same role
  • The long-term goal is efficiency and precision. Any compatibility is accidental. This is a rewrite, not a translation
  • It builds upon blockflow, which is a from-scratch implementation of all primitives
  • Because the only thing we carry over is the equations we aim to solve, we consider the code code sufficiently separate to warrant a separate license (moving away from the mix of licences in our translation, clearmap-rs)

Finally, this is used to drive development of blockflow by piloting it on a serious use case

Below is random text for the LLM. Cleaning to be done later; it might not be accurate

The three parts

params — two regimes, one implementation

Every stage runs in either ClearMap's parameters or better ones, from the same code. The rule that makes that enforceable rather than aspirational:

ClearMap's choices must be a point in the parameter space, never a constant in the code.

It fails in two directions, and one sweep catches both:

failure what it looks like how check_parameter_space sees it
hard-coded an arbitrary choice the field is there, the stage ignores it a step off that axis changes nothing → dead_axes
cannot reproduce ClearMap the point is not in the space the run at ClearMap's point differs from the oracle → reproduces_clearmap

Tuned<T> is the carrier: both the value in use and ClearMap's own, so Regime is derived rather than stored and parity mode is a restore rather than a rebuild. Parameters::axes() lists one departure per arbitrary choice, and a choice missing from that list is a choice nothing checks.

The second failure is the one that matters for the claim this crate makes about itself, and it is the one a parity test cannot find on its own: a stage with a decorative parameter passes parity perfectly while having ClearMap's choice baked in — it reproduces the original instead of being able to.

oracle — which comparison, and what it is allowed to claim

tier what is compared where it runs
1, routine this crate against clearmap-rs in process, here
2, periodic clearmap-rs against ClearMap out of band; recorded in oracle::registry

A tier-1 comparison establishes parity only where that stage's own parity against ClearMap is digest-recorded, and agreement otherwise. The failure to design against:

the new pipeline agreeing with clearmap-rs on a stage neither has ever had checked against ClearMap. Two implementations sharing an assumption agree perfectly.

So the label is structural, not a flag:

  • ClearMapDigest's fields are private to oracle::registry, so no other module — not just no other crate — can mint one. Every digest in the program is a pub const there, and every one names the run it came from.
  • Stage::CLEARMAP is an associated const pointing at one of those, or None.
  • Comparison::run derives the claim from it. There is no parity: bool, no tier argument and no expect_parity.
  • Verdict's fields are private and it has no constructor, so an agreement result cannot be edited into a parity one.
  • A digest names its stage, and a comparison is refused when the stage it names is not the stage being compared.

Two compile_fail doctests pin the first and last of those; tests/tier_labelling.rs walks the paths that are expressible and shows each landing on Agreement or an error.

stages

stage what it is for how it is checked
clip the pipeline's first step, and the proof the machinery works end to end tier 1 against clearmap-rs, licensed as parity by a recorded ClearMap digest
crop run a subset, without a probe binary's --crop-start hacks the sub-array identity, which is exact and needs no oracle

clip was chosen because it is the simplest stage carrying a recorded ClearMap digest, and without a digest a stage cannot exercise the parity path at all. Its four arbitrary choices — clip range, bin count, whether the top of the range is bins or bins - 1, and truncate versus round — are fields, not constants.

crop has no ClearMap counterpart and does not need one. crop(V, r) must be V[r], which is exact, holds for every region rather than for one fixture, and needs nothing outside this crate. A crop is a translated block grid — configuration of machinery blockflow already has, not a new framework feature.

Two limitations of crop, both real

Reach 0 at the crop's own boundaries. BlockGeometry::derive treats a read clamped at the phase's volume edge as trustworthy, which is right at a real edge of the array — and a cropping phase's edges are not edges of the image below. Until reach carries its coordinate space (change 6 of BLOCK_OPS.md), a cropping phase must use reach 0, so anything with a reach runs as a later phase over the cropped image and sees the crop edge as a real edge. tests/crop_identity.rs asserts the interior region explicitly and separately asserts that the boundary shell really does differ.

Chunk alignment is the fast path, and it is visible rather than enforced. Crop::alignment predicts from the plan; ZarrEnvironment::unaligned_reads() counts what happened; tests/crop_alignment_zarr.rs asserts the two agree in both directions, and that the output is byte-identical either way. Alignment is performance, never correctness.

Building

cargo test                  # 22 tests, no storage backend
cargo test --features zarr  # 25, adding the alignment counters

zarr is off by default for the same reason it is off in blockflow: zarrs is 118 extra packages and a C toolchain.

Note on build location. /data is at 100 %. A debug build of this crate plus its clearmap-rs dev-dependency is around 10 GB, so set CARGO_TARGET_DIR somewhere on /big rather than letting target/ land here.

Benchmark: this crate against ClearMap Python, materialized

One crop, both sides, same machine, minutes apart. Every figure is a minimum over repeats, and both sides reproduce the recorded surfaces — the harness on this side refuses to print a timing at all if the answer moved.

  • cropvasc-pipeline-axial, 96 x 96 x 32 = 294 912 voxels, crop_start = [16, 640, 1792], --resolution 1 1 1
  • ClearMap Python — 1 process, OMP_NUM_THREADS=1, --unblocked so every stage runs on the whole volume, best of 3
  • clearmap-ng — one block (block = 96 over a 96-wide crop, so the grid is the volume), 40 rayon threads, best of 5
op ClearMap materialized ns/voxel clearmap-ng ns/voxel py / ng vs ClearMap blocked
binarize 1.5526 s 5 265 1.6643 s 5 643 0.93x — slower 0.94x
smooth 0.3385 s 1 148 0.0443 s 150 7.6x 17.5x
binary_fill * 0.1073 s 364 0.0140 s 48 7.7x 8.4x
combine_binary 0.7203 s 2 442 0.0454 s 154 15.9x 18.0x
skeletonize 0.1539 s 522 0.0729 s 247 2.1x 5.9x
build_graph 4.1384 s 14 033 0.0054 s 18 766x 793x
clean_graph 0.0206 s 70 0.0004 s 1.2 52x 56x
reduce_graph 0.0320 s 109 0.0003 s 1.1 107x 119x
matched total 7.064 s 1.847 s 3.8x 4.4x

The one number that matters

binarize is 90.1% of this crate's entire chain — 1.664 s of 1.847 s — and it is the only op that loses. Every other op is already under 73 ms: making all seven of them infinitely fast would take the chain from 1.847 s to 1.664 s, a 10% win. In ClearMap the same stage is 22% of its chain, whose bulk is build_graph at 59%.

So the aggregate 3.8x is bought entirely by ops that were already cheap, and the one op that costs this crate anything is one it loses — to a single-threaded Python program, on forty threads. An aggregate that hides that is worse than no aggregate.

ClearMap's own blocking costs it 2-3x

The last column is the same reference run without --unblocked, at its shipped --block 40,6,5. Blocking is not free for it: smooth goes 0.339 -> 0.774 s and skeletonize 0.154 -> 0.434 s, while binarize does not move because this crop is already one block for it. So "materialized" is true of ClearMap's ops — they are whole-array numpy/scipy — but not of its volume handling, which blocks by default and at mp.cpu_count() processes rather than the 1 used here.

Three things that would mislead you

build_graph, clean_graph and reduce_graph are not per-voxel speedups and must not be scaled. They run at 5.4 / 0.4 / 0.3 ms on this side — at or near the timer's resolution — while the reference's figures are dominated by graph_tool object construction and Python loop overhead that does not grow with the volume. A 766x ratio on a 1 684-vertex graph says nothing about a million-vertex one.

* binary_fill is inert on this crop. filled.npy is smoothed.npy byte for byte, so both sides are timed searching for nothing to fill. the_filled_surface_is_inert_on_this_crop asserts the identity, so this row cannot quietly become a real comparison.

† The reference's combine_binary produces combined and final, and its separate final stage is an idempotent re-run of the second half. combine::run here produces both in one call, so it maps to combine_binary alone; adding final would double-charge it.

The crop is smaller than binarize's own windows, and it matters

The shipped parameters use an adaptive selem of 250 x 250 x 3 and an equalize selem of 200 x 200 x 5 against a 96 x 96 cross-section, so both are clipped by the volume. Repeating the reference at 256 x 256 x 32 = 2 097 152 voxels, where they are not, moves the per-voxel figures a long way:

op ns/voxel @ 96 x 96 x 32 ns/voxel @ 256 x 256 x 32
binarize 5 265 6 801 (+29%)
smooth 1 148 282
combine_binary 2 442 488
skeletonize 522 299
build_graph 14 033 1 992 (-86%)

So the small crop understates binarize by about a third and overstates every other op, and build_graph's collapse is the direct evidence for the "do not scale these" warning above — a stage dominated by a per-process constant gets cheaper per voxel as the volume grows.

At that larger crop binarize is the only op timed on both sides so far:

seconds ns/voxel mask set
ClearMap Python, materialized, 1 process 14.2628 6 801 265 919
clearmap-ng, one block, 40 threads 13.6370 6 503 265 907

Read that as parity on speed, not as a 5% win. The two runs are separate processes at different loads, and 5% is inside the drift this machine shows between nominally identical runs. What survives is the finding above: forty threads against one Python process is roughly a draw, and binarize is the whole of this crate's cost.

An open parity question at crops larger than the recorded one

The two masks above are not the same answer. They agree to 12 in the count — and that is a coincidence of cancellation, not agreement: 2 618 voxels differ, 1 315 that ClearMap sets and this crate does not and 1 303 the other way, 0.125% of the volume, spread through it rather than at a boundary.

At the recorded 96 x 96 x 32 crop the two are byte-identical — that is what registry::STEP_BINARIZE is — so this divergence is invisible to every fixture in the tree. The likely site is the sampling lattice: at 96 x 96 both the adaptive and equalize windows are clipped to a single sample, and at 256 x 256 they are not, so LocalStatistics' grid interpolation is exercised for the first time. This is not diagnosed, and the timing pair above is a comparison between two implementations that disagree, which is why it is reported as a bound rather than folded into the table.

tests/binarize_gap.rs refused to print a timing for this crop for exactly that reason, which is the harness working as designed.

The loads were not equal: the reference ran at load 10-19 and this crate at 5.5-6.4, so this crate had the easier window. That biases in its favour, which makes the binarize result stronger rather than weaker — it loses anyway.

Where that cost lives is already localised, though against clearmap-rs rather than against Python: docs/performance.md §8-9 puts it inside vesselize, and within vesselize in the percentile rather than the tubeness.

The cell pipeline, measured separately

The table above is the tube/vasculature chain. The cell chain is measured on its own, because it is dominated by ops none of the above exercises: a grey opening by a 10 x 10 disk and a local-maximum test over a 5 x 5 x 5 box.

At a 384 x 384 x 256 crop, active chain only — background correction plus maxima detection, medians of three fresh-process runs a side:

seconds ns/voxel
ClearMap Python, 1 process 3.056 81.0
clearmap-ng, best cut, 40 threads 2.816 74.6
1.09x faster

Every mask is asserted against ClearMap's recorded maxima_raw.npy at every cut measured, so this is a like-for-like time on an identical answer. This started at roughly 16x slower. Four things carry the distance.

The plan is built at the input's width. Program::assemble takes the source Dtype, so a uint16 volume runs the rank filters at uint16 and widens only where the subtraction needs f64 — worth about 1.2x. The fan-in's identity arm has to be a WidenOp rather than a VoxelwiseMapOp for this to plan at all; the latter refuses uint16.

A rank filter that asks for an end uses running extrema. Van Herk's algorithm costs three operations a voxel whatever the window is, against one per tap. A box is better still — an extremum over a box factorises, so three separable passes answer a 5 x 5 x 5 window that the gather pays 125 taps for. Getting the general case to apply took three tries, and the first is worth recording: decomposing along the fastest axis made the kernel 1.7x faster in isolation and moved the chain by nothing, because a disk stated as 10 x 10 x 1 is flat on that axis and decomposes into eighty runs of one.

A phase that copies a volume is not free. With the difference of Gaussians off — ClearMap's own point — that phase was an identity, kept so the image numbering would not move with a parameter, and described as one voxelwise pass. A materialisation is a bandwidth operation: it measured 0.31 s, 12% of the run, saturating memory at 1.4 cores' worth of the executor. The plan now points the handle at the phase that already holds those values. PlanBuilder::pixels makes the argument itself when it refuses an empty chain.

Block shape is a cache question, not a halo question. A slab cut on the axis the element does not reach along was expected to buy blocks nearly free; [16,384,256] measured 6.30 s where [128,128,64] measured 3.99 s.

The measurement that redirected all of this is tests/cells_op_profile.rs, which sums the executor's own OpApplied events per op and per phase. A one-block ladder had attributed 99.4% of the chain to the phase the opening sits in; per op that phase is about half, and per phase the identity copy and the flood were a third of the wall between them — none of which a phase-level ladder can show.

Reproduce with:

CLEARMAP_NG_FIXTURES=... CLEARMAP_NG_CELLS_BIG=1 \
  cargo test --release --test cells_block_sweep -- --ignored --nocapture

# one cut, for a figure comparable with a single reference run
CLEARMAP_NG_FIXTURES=... CLEARMAP_NG_CELLS_BIG=1 CLEARMAP_NG_CELLS_CUT=128,128,64 \
  cargo test --release --test cells_block_sweep -- --ignored --nocapture

# where the time goes, per op and per phase
CLEARMAP_NG_FIXTURES=... CLEARMAP_NG_CELLS_BIG=1 \
  cargo test --release --test cells_op_profile -- --ignored --nocapture

Reproducing it

Both sides need the machine-local recordings and the vendored ClearMap tree; see docs/performance.md §1 for the environment.

export CLEARMAP_NG_FIXTURES=/husky/henriksson/clearmap-img/outputs
IN="$CLEARMAP_NG_FIXTURES/../inputs/osf-sa3x8/14-16-41_tricocktail_UltraII[05 x 06]_C00_UltraII Filter0000.ome.tif"

# the reference, materialized. Drop --unblocked for the last column.
PYTHONPATH=<clearmap-rs>/ClearMap-py OMP_NUM_THREADS=1 \
  ~/miniconda3/envs/clearmap-gt2/bin/python \
  <clearmap-rs>/tools/vasculature_pipeline_probe.py \
  --input "$IN" --crop-start 16,640,1792 --crop-size 96,96,32 \
  --resolution 1,1,1 --block 40,6,5 --n-processes 1 --unblocked \
  --output-dir <dir>/py
# per-stage seconds and ns/voxel land in <dir>/py/progress.json

# this side
CLEARMAP_NG_PERF_REPEATS=5 CLEARMAP_NG_PERF_OUT=<dir>/ng \
  cargo test --release --test reference_timing -- --ignored --nocapture --test-threads=1

# the join
<clearmap-rs>/tools/pipeline_timing_join.py --python <dir>/py --rust <dir>/ng

For the larger crop, pass --crop-size 256,256,32 to the probe. This side has no harness at that crop — reference_timing's CROP is a constant tied to the recording — so binarize was taken through tests/binarize_gap.rs, pointed at the probe's output by symlinking it under a fixture root:

mkdir -p <dir>/fixtures && ln -sfn <dir>/py <dir>/fixtures/vasc-pipeline-axial
CLEARMAP_NG_FIXTURES=<dir>/fixtures CLEARMAP_NG_GAP_REPEATS=3 \
  cargo test --release --test binarize_gap \
  the_gap_between_the_two_implementations -- --ignored --nocapture --test-threads=1

It will refuse to print a timing and fail on the mask count, which is the divergence above rather than a broken invocation.

Run them one at a time and record the load. This machine is shared; a concurrent --release build in another checkout moved these absolutes by more than the differences being measured. progress.json carries machine.loadavg on both sides for exactly that reason, and a join across incomparable loads is not evidence.

Registration Phase Benchmark

Date: 2026-08-28

This records a phase-separated benchmark for the real-texture registration trial. The fixture uses local OSF-style ClearMap autofluorescence OME-TIFF planes from:

/husky/henriksson/clearmap-img/inputs/osf-sa3x8/autofluorescence

The benchmark crops planes Z0420..Z0451 to a 96 x 96 x 32 volume at crop_xyz = [520, 920, 420]. The moving volume is a shifted copy of the real crop (shift_xyz = [2, -1, 1]) and both inputs are normalized to [0, 1] before original ClearMap/elastix registration. This is a real microscope texture benchmark for iteration; it is not a full microscope-to-atlas run.

Command:

<clearmap-rs>/tools/benchmark_registration_phases.py \
  --repeats 3 \
  --work-root /tmp/clearmap-ng-registration-phase-bench \
  --json-out /tmp/clearmap-ng-registration-phase-bench.json

The full raw JSON is stored at /tmp/clearmap-ng-registration-phase-bench.json for the run below.

All timing/RSS rows below are per-command /usr/bin/time -v measurements except rows marked internal, which come from the existing ClearMap/wrapper or clearmap-ng report payloads. Values are mean over 3 runs with min/max in parentheses.

phase wall s CPU s peak RSS
prepare subset 0.364 (0.338..0.383) 5.35 (5.21..5.46) 56.0 MB
original script total 6.522 (6.286..6.752) 88.29 (86.96..89.25) 126.7 MB
original ClearMap internal 5.085 (4.789..5.322) 76.80 (75.81..77.63) 126.7 MB
direct elastix 5.291 (5.087..5.527) 75.28 (74.73..75.77) 56.8 MB
direct transformix 0.136 (0.128..0.151) 0.12 (0.11..0.13) 39.2 MB
wrapper total 6.506 (6.455..6.552) 79.13 (78.67..79.55) 56.7 MB
wrapper internal 5.321 (5.267..5.398) 77.23 (76.83..77.53) 56.7 MB
native report total 2.277 (2.221..2.363) 2.26 (2.20..2.35) 21.0 MB
native internal 1.997 (1.949..2.076) 2.00 (1.95..2.08) 21.0 MB

Mean ratios:

ratio wall CPU RSS
direct elastix / native internal 2.650x 37.701x 2.697x
original ClearMap internal / native internal 2.547x 38.466x 6.018x
wrapper internal / native internal 2.665x 38.679x 2.695x
original script total / native report total 2.864x 39.065x 6.018x
wrapper total / direct elastix 1.230x 1.051x 0.999x
original ClearMap internal / direct elastix 0.961x 1.020x 2.232x

Output agreement:

  • Original ClearMap/elastix point: [48.929583, 47.178869, 16.763063]
  • Wrapper point: [48.929583, 47.178869, 16.763063]
  • Native point: [48.757321807536, 47.24119067785, 16.73328319771]
  • Wrapper-vs-original max error: 0.0
  • Native-vs-original max error: 0.172261192464 voxels under tolerance 0.2
  • Fixed input FNV-1a64: 76376a4a79f1be89
  • Moving input FNV-1a64: 9facb4b5d400704a
  • Original transform FNV-1a64: ec92bf47c9839c9d
  • Wrapper transform FNV-1a64: ec92bf47c9839c9d

Interpretation:

  • Original ClearMap and the wrapper are doing the same elastix work on this fixture: same inputs, same parameter file, same transformed point, and byte-identical transform output.
  • Original ClearMap internal time and direct elastix time are close, so most of the original wall time is elastix rather than Python.
  • Wrapper total is slower than direct elastix because it includes wrapper process startup, version probes, fixture parsing, .npy to .mhd preparation outside the internal timer, elastix, and transformix.
  • Native clearmap-ng is faster and uses less RSS on this trial, but this is an alternative deterministic optimizer, not the same elastix algorithm.
  • CPU time is not a fair single-core comparison: although elastix was invoked with -threads 1, the ClearMap-style parameter file omits UseMultiThreadingForMetrics, and elastix 4.900 defaults that setting to true. That is why direct elastix reports about 75 CPU seconds for about 5 wall seconds. Native CPU time is close to native wall time.

Licence

MIT — see LICENSE, and the same as blockflow.

This is still the domain crate: it holds ClearMap's arbitrary choices as named points in a parameter space, where blockflow's ops know none of them and nothing here may be moved into it. That separation is about generality and the direction of the dependency — blockflow is depended on by things that have never heard of a lightsheet, and a domain term filed there is a term all of them carry — not about licence terms, which are now the same on both sides.

clearmap-rs is GPL-3.0 and stays so. It is a dev-dependency only: it is the oracle, nothing under src/ references it, and nothing shipped links it.

About

No description, website, or topics provided.

Resources

Stars

0 stars

Watchers

0 watching

Forks

Releases

Packages

Contributors

Languages