volara-metrics 0.1.1

Blockwise tasks for computing simple dataset metrics

volara-metrics

A volara plugin for computing and visualizing various dataset metrics.

Metrics

Segment Persistence metric

This metric is a useful proxy for seeing how well segments persist across the z dimension, without ground truth data.

It works by using a sliding window (defined by k) and comparing unique segments on the left slab with the right slab of the window. Using some overlap score (i.e IoU), we can generally see how well segments are persisting in a comparative sense (i.e some section vs other sections). This should not be considered in an absolute way - as there are valid ways for segments to end across z (true splits, spines, etc). But it allows for seeing how well a target section (e.g. a registered boundary) segments compared to the rest of the sections.

If you have a segmentation, it can be run like below. Here we just assume the full roi is 200 sections, and we thus tile over it in z columns by setting the block size in z to the roi size in z. This will return a heatmap of the segment persistence that shows region dependent segmentation quality.

from funlib.geometry import Coordinate
from pathlib import Path
from volara.datasets import Labels, Raw
from volara-metrics import SegmentPersistenceTask

f = Path("test.zarr")

labels_data = Labels(store=f / "segments")
heatmap_data = Raw(store=f / "heatmap")

# assuming full z roi is 200 sections, can just make that the block size
block_size = Coordinate(200, 128, 128)

persistence_task = SegmentPersistenceTask(
    labels_data=labels_data,
    heatmap_data=heatmap_data,
    block_size=block_size,
    num_workers=10,
    k=3
)

persistence_task.run_blockwise()

Since segments can be tedious to get for larger volumes (i.e requiring full agglomeration), a quicker option is to simply use the affinities and compute a pseudo segmentation prior to computing the persistence. This is likely less accurate than using agglomerated segments, but seems to correlate well enough and can be run faster.

from funlib.geometry import Coordinate
from pathlib import Path
from volara.datasets import Affs, Raw
from volara-metrics import SegmentPersistenceTask

f = Path("test_jm200_crop_updated.zarr")

affs_data = Affs(
    store=f / "affs",
    neighborhood=[
        Coordinate(1, 0, 0),
        Coordinate(0, 1, 0),
        Coordinate(0, 0, 1),
    ])

heatmap_data = Raw(store=f / "heatmap")

# assuming full z roi is 200 sections, can just make that the block size
block_size = Coordinate(200, 128, 128)

persistence_task = SegmentPersistenceTask(
    affs_data=affs_data,
    heatmap_data=heatmap_data,
    block_size=block_size,
    num_workers=10,
    k=3,
    fg_thresh=0.35,
    seed_thresh=0.75,
    smooth_sigma=1.0,
    z_link_thresh=0.5,
)

persistence_task.run_blockwise()

We can also return the per block metrics (as block .npz files) and then aggregate across boundaries later. This allows us to compute statistics for a large volume (without requiring in memory processing) while still retaining region dependent quality.

from funlib.geometry import Coordinate
from pathlib import Path
from volara.datasets import Labels, Raw
from volara-metrics import SegmentPersistenceTask
from volara-metrics.utils import aggregate_block_metrics

f = Path("test.zarr")

labels_data = Labels(store=f / "segments")
heatmap_data = Raw(store=f / "heatmap")

# assuming full z roi is 200 sections, can just make that the block size
block_size = Coordinate(200, 128, 128)

persistence_task = SegmentPersistenceTask(
    labels_data=labels_data,
    heatmap_data=heatmap_data,
    block_size=block_size,
    num_workers=10,
    k=3,
    return_metrics=True,
    save_raw_scores=True,
    metrics_section_boundary=100,
)

persistence_task.run_blockwise()

metrics_dir = f / "heatmap" / "metrics"

summary = aggregate_block_metrics(
    metrics_dir,
    stitched_boundary_idx=100,
)

print(summary)

Slice raggedness metric

This metric scores the raggedness of a cut surface -- a proxy for how much tissue a physical slice disrupted, and thus how poorly it is expected to register/persist across the cut. It runs on a 2D (y, x) heightmap z(y, x) (the cut-surface height field), with the lateral pixel pitch given in the same physical unit as the heights so slopes are dimensionless.

It works by removing the smooth form of the surface (a robust polynomial: order 1 = tilt, order 2 = tilt + gentle bow), then measuring the residual with derivative-based statistics: a 90th-percentile local slope (slope_p90) and the developed-area ratio (sdr = mean(sqrt(1+|grad|^2)) - 1). Both are mapped onto a single physical axis -- surface tilt as a fraction of vertical -- and combined into raggedness_index in [0, 1). Because the axis has a fixed ceiling (vertical), the index is absolute and batch-independent: it needs no fitted constant, no calibration, and no cohort. The task writes a per-pixel slope- magnitude roughness heatmap (so you can see where a face is ragged) and saves the scalar summary as per-block .npz files.

The blockwise task needs the volara/daisy stack: pip install volara-metrics[blockwise]. The pure metric kernel (volara_metrics.slice_metric) imports with only numpy/scipy/scikit-image, so downstream libraries can reuse the exact same math without the blockwise dependencies.

Run it on an existing heightmap. "whole" mode (the default) computes the exact metric in a single block -- appropriate when a face fits in memory:

from funlib.geometry import Coordinate
from pathlib import Path
from volara.datasets import Raw
from volara_metrics import SliceMetricTask
from volara_metrics.utils import aggregate_slice_metrics

f = Path("test.zarr")

heightmap_data = Raw(store=f / "gel2_top_heightmap")   # 2D (y, x), heights in same unit as pitch
heatmap_data = Raw(store=f / "gel2_top_roughness")     # 2D (y, x) float32 output

task = SliceMetricTask(
    heightmap_data=heightmap_data,
    heatmap_data=heatmap_data,
    block_size=Coordinate(512, 512),   # zarr chunking; the whole face is one daisy block
    mode="whole",
    detrend_order=1,                   # 1 = remove tilt, 2 = tilt + bow
    num_workers=1,
)
task.run_blockwise()

summary = aggregate_slice_metrics(f / "gel2_top_roughness" / "metrics")["summary"]
print(summary)   # {'slope_p90': ..., 'sdr': ..., 'raggedness_index': ..., ...}

The upstream heightmap is usually produced by reducing a 3D surface-prediction mask over z. That reducer, ComputeHeightmapTask, lives in volara-registration; compose the two:

from volara_registration import ComputeHeightmapTask   # reducer="max" -> top, "min" -> bot

ComputeHeightmapTask(
    seg=Raw(store=f / "gel2_top_surface_labels"),
    heightmap=heightmap_data,
    surface_label=1, reducer="max",
    block_size=Coordinate(512, 512),
).run_blockwise()

SliceMetricTask(heightmap_data=heightmap_data, heatmap_data=heatmap_data,
                block_size=Coordinate(512, 512)).run_blockwise()

For a face too large to hold in memory, use mode="tiled": the global detrend is fit once (subsampled) and each block subtracts it analytically, saving mergeable partials. slope_mean/slope_rms/sdr then aggregate exactly and slope_p90 is recovered from a merged histogram.

task = SliceMetricTask(
    heightmap_data=heightmap_data,
    heatmap_data=heatmap_data,
    block_size=Coordinate(1024, 1024),
    mode="tiled",
    num_workers=10,
)
task.run_blockwise()

summary = aggregate_slice_metrics(f / "gel2_top_roughness" / "metrics")["summary"]
print(summary)

Save the aggregated scores and render the diagnostics (mirrors the persistence "save scores + plots" workflow; the plotting needs the [plot] extra for matplotlib):

from volara_metrics.utils import (
    save_aggregated_slice_metrics,
    load_aggregated_slice_metrics,
)
from volara_metrics.plotting import (
    PlotConfig,
    SliceDiagnosticsConfig,
    plot_all_slice_diagnostics,
)
from funlib.persistence import open_ds

metrics_dir = f / "gel2_top_roughness" / "metrics"

save_aggregated_slice_metrics(metrics_dir, return_only_summary=False)
data = load_aggregated_slice_metrics(str(metrics_dir) + "/aggregated_slice_full.npz")

field = open_ds(f / "gel2_top_roughness")[:]   # the 2D roughness heatmap

plot_all_slice_diagnostics(
    field,
    data=data,
    config=SliceDiagnosticsConfig(
        plot=PlotConfig(view=False, save=True, save_dir="all_plots", dpi=200),
    ),
)

This writes roughness_heatmap.png (the local slope-magnitude map), slope_histogram.png (the |grad z| distribution with slope_mean/slope_p90 marked), and slice_summary.png (heatmap + histogram + the scalar raggedness table). The individual plot_roughness_heatmap / plot_slope_histogram / plot_slice_summary functions are available too.