Metadata-Version: 2.4
Name: volara-metrics
Version: 0.1.0
Summary: Blockwise tasks for computing simple dataset metrics 
Requires-Python: >=3.11
Requires-Dist: numpy
Requires-Dist: scikit-image>=0.25.2
Requires-Dist: scipy>=1.15.3
Provides-Extra: blockwise
Requires-Dist: daisy>=1.2.2; extra == 'blockwise'
Requires-Dist: funlib-geometry>=0.3; extra == 'blockwise'
Requires-Dist: funlib-persistence>=0.7.0; extra == 'blockwise'
Requires-Dist: iohub>=0.3.0a5; extra == 'blockwise'
Requires-Dist: s3fs; extra == 'blockwise'
Requires-Dist: volara>=1.0.2; extra == 'blockwise'
Provides-Extra: dev
Requires-Dist: daisy>=1.2.2; extra == 'dev'
Requires-Dist: funlib-geometry>=0.3; extra == 'dev'
Requires-Dist: funlib-persistence>=0.7.0; extra == 'dev'
Requires-Dist: iohub>=0.3.0a5; extra == 'dev'
Requires-Dist: pytest-cov>=5.0; extra == 'dev'
Requires-Dist: pytest>=8.0; extra == 'dev'
Requires-Dist: ruff>=0.6; extra == 'dev'
Requires-Dist: s3fs; extra == 'dev'
Requires-Dist: ty>=0.0.24; extra == 'dev'
Requires-Dist: volara>=1.0.2; extra == 'dev'
Provides-Extra: docs
Requires-Dist: autodoc-pydantic>=2.0; extra == 'docs'
Requires-Dist: myst-parser>=4.0; extra == 'docs'
Requires-Dist: sphinx-autodoc-typehints>=2.0; extra == 'docs'
Requires-Dist: sphinx-material>=0.0.36; extra == 'docs'
Requires-Dist: sphinx-rtd-theme>=2.0; extra == 'docs'
Description-Content-Type: text/markdown

# volara-metrics

A [volara](https://github.com/e11bio/volara) plugin for computing and
visualizing various dataset metrics.

## Metrics

<details>
  <summary><b>Segment Persistence metric</b></summary>

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.

```py
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. 

```py
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.

```py
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)
```

</details>

<details>
  <summary><b>Slice raggedness metric</b></summary>

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:

```py
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](https://github.com/e11bio/volara-registration); compose the
two:

```py
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.

```py
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)
```

</details>
