Add group-level Otsu thresholding (groupotsu+k{}i{} seg method)

docs/howto/segmentation.md

@@ -10,10 +10,12 @@ flowchart LR
     B --> C{seg_method}
     C -->|threshold| D[Threshold]
     C -->|otsu+k3i2| E[Multi-Otsu + k-means]
-    D --> F[Binary Mask]
-    E --> F
-    F --> G[Clean: remove edge / small objects]
-    G --> H[Segmentation Output]
+    C -->|groupotsu+k3i2| F[Group Multi-Otsu]
+    D --> G[Binary Mask]
+    E --> G
+    F --> G
+    G --> H[Clean: remove edge / small objects]
+    H --> I[Segmentation Output]
 ```
 
 After segmentation the binary mask is used to compute [field fraction](../reference/outputs.md#field-fraction-map-seg), [object count](../reference/outputs.md#object-count-map-seg), and [per-object region properties](../reference/outputs.md#region-properties-statistics-table-tabular).
@@ -81,6 +83,45 @@ stain_defaults:
 
 **Limitations:** Can fail on images with unusual histograms (e.g. very sparse pathology that does not form a distinct peak) or when the background is very noisy.
 
+### Group Multi-Otsu (`seg_method: groupotsu+k3i2`)
+
+A variant of Multi-Otsu that derives a **single shared threshold from the aggregate histogram of all subjects** rather than computing a threshold independently per image.  This is preferred when subjects were acquired with common acquisition settings and you want to ensure consistent, comparable quantification across the cohort.
+
+The workflow is a two-step process:
+
+**Step 1 — compute group threshold** (run once for the whole cohort):
+
+```bash
+spimquant /bids /output participant \
+    --targets all_group_otsu \
+    --seg_method groupotsu+k3i2
+```
+
+This triggers:
+
+1. For each subject: compute a percentile-clipped intensity histogram from the bias-field corrected image and save it as an NPZ file.
+2. Aggregate all subject histograms onto a common intensity grid, apply multi-level Otsu thresholding, and save the resulting thresholds as a JSON file in `{output}/group/`.
+
+**Step 2 — segment each subject using the group threshold**:
+
+```bash
+spimquant /bids /output participant \
+    --seg_method groupotsu+k3i2
+```
+
+Each subject's binary mask is produced by applying the group-level threshold from the JSON file.  A per-subject PNG is also generated showing the group threshold overlaid on the individual histogram, useful for visual quality control.
+
+**Config key:**
+
+```yaml
+seg_method:
+  - groupotsu+k3i2
+```
+
+**When to use:** Preferred when a batch of subjects shares the same acquisition protocol and you want consistent thresholding across subjects.  Reduces subject-to-subject variability in the segmentation boundary that can occur with per-subject Otsu.
+
+**Limitations:** Less adaptive than per-subject Otsu — if staining intensity varies substantially across subjects (e.g. due to different batches of antibody or tissue preparation), a single group threshold may over- or under-segment some subjects.
+
 ---
 
 ## Post-Segmentation Cleaning

spimquant/workflow/Snakefile

@@ -85,6 +85,9 @@ stains_for_seg = list(set(config["stains_for_seg"]).intersection(set(stains)))
 # seg methods that use multi-Otsu thresholding (otsu+k{}i{} pattern)
 otsu_seg_methods = [m for m in config["seg_method"] if m.startswith("otsu+")]
 
+# seg methods that use group-level multi-Otsu thresholding (groupotsu+k{}i{} pattern)
+groupotsu_seg_methods = [m for m in config["seg_method"] if m.startswith("groupotsu+")]
+
 if len(stains_for_seg) == 0 or config["no_segmentation"]:
     do_seg = False
     do_coloc = False
@@ -908,6 +911,43 @@ rule all_group:
         rules.all_group_stats_coloc.input if do_coloc else [],
 
 
+rule all_group_otsu:
+    """Target rule for group-level Otsu threshold computation.
+
+    Aggregates intensity histograms across all subjects for each stain and
+    computes a single set of Otsu thresholds that can be applied consistently
+    to every subject.  Run this rule before participant-level segmentation
+    when ``groupotsu+k{}i{}`` is included in ``seg_method``.
+
+    Example::
+
+        spimquant /bids /output participant \\
+            --targets all_group_otsu \\
+            --seg_method groupotsu+k3i2
+
+    Then run segmentation using the computed group thresholds::
+
+        spimquant /bids /output participant \\
+            --seg_method groupotsu+k3i2
+    """
+    input:
+        expand(
+            bids(
+                root=root,
+                datatype="group",
+                stain="{stain}",
+                level="{level}",
+                desc="{desc}",
+                suffix="thresholds.json",
+            ),
+            stain=stains_for_seg,
+            level=config["segmentation_level"],
+            desc=groupotsu_seg_methods,
+        )
+        if (do_seg and groupotsu_seg_methods)
+        else [],
+
+
 include: "rules/import.smk"
 include: "rules/masking.smk"
 include: "rules/templatereg.smk"

spimquant/workflow/rules/groupstats.smk

@@ -4,9 +4,67 @@ Group-level statistical analysis rules for SPIMquant.
 This module performs group-based statistical tests on segmentation statistics
 (e.g., fieldfrac, density, volume) across participants, using metadata from
 participants.tsv to define contrasts.
+
+It also provides the ``group_otsu`` rule which aggregates per-subject intensity
+histograms (produced by ``compute_subject_histogram``) to compute a single set
+of Otsu thresholds shared across all subjects.
 """
 
 
+rule group_otsu:
+    """Compute group-level Otsu thresholds from aggregated per-subject histograms.
+
+    Collects the intensity histogram NPZ files computed by
+    ``compute_subject_histogram`` for all subjects, merges them onto a
+    common intensity grid, and applies multi-level Otsu thresholding to the
+    aggregate histogram.  The resulting thresholds are saved as a JSON file
+    (consumed by ``multiotsu_group`` during participant-level segmentation)
+    and as a PNG figure for visual inspection.
+
+    This rule is the target of ``all_group_otsu`` and should be run before
+    participant-level segmentation when ``groupotsu+k{}i{}`` is used as the
+    segmentation method.
+    """
+    input:
+        histogram_npz=lambda wildcards: inputs["spim"].expand(
+            bids(
+                root=work,
+                datatype="seg",
+                stain=wildcards.stain,
+                level=wildcards.level,
+                desc="groupotsu+k{k}i{i}".format(k=wildcards.k, i=wildcards.i),
+                suffix="histogram.npz",
+                **inputs["spim"].wildcards,
+            )
+        ),
+    params:
+        otsu_k=lambda wildcards: int(wildcards.k),
+        otsu_threshold_index=lambda wildcards: int(wildcards.i),
+    output:
+        thresholds_json=bids(
+            root=root,
+            datatype="group",
+            stain="{stain}",
+            level="{level}",
+            desc="groupotsu+k{k,[0-9]+}i{i,[0-9]+}",
+            suffix="thresholds.json",
+        ),
+        thresholds_png=bids(
+            root=root,
+            datatype="group",
+            stain="{stain}",
+            level="{level}",
+            desc="groupotsu+k{k,[0-9]+}i{i,[0-9]+}",
+            suffix="thresholds.png",
+        ),
+    threads: 4
+    resources:
+        mem_mb=8000,
+        runtime=10,
+    script:
+        "../scripts/group_otsu.py"
+
+
 rule perform_group_stats:
     """Perform group-based statistical tests on segmentation statistics.
 

spimquant/workflow/rules/segmentation.smk

@@ -78,19 +78,14 @@ rule n4_biasfield:
         shrink_factor=16 if config["sloppy"] else 1,
         target_chunk_size=512,  #this sets the chunk size for this and downstream masks
     output:
-        corrected=temp(
-            directory(
-                bids(
-                    root=work,
+        corrected=bids(
+                    root=root,
                     datatype="seg",
                     stain="{stain}",
                     level="{level}",
                     desc="correctedn4",
-                    suffix="SPIM.ome.zarr",
+                    suffix="SPIM.ozx",
                     **inputs["spim"].wildcards,
-                )
-            ),
-            group_jobs=True,
         ),
     threads: 128 if config["dask_scheduler"] == "distributed" else 32
     resources:
@@ -110,12 +105,12 @@ rule multiotsu:
     """
     input:
         corrected=bids(
-            root=work,
+            root=root,
             datatype="seg",
             stain="{stain}",
             level="{level}",
             desc="corrected{method}".format(method=config["correction_method"]),
-            suffix="SPIM.ome.zarr",
+            suffix="SPIM.ozx",
             **inputs["spim"].wildcards,
         ),
     params:
@@ -152,6 +147,109 @@ rule multiotsu:
         "../scripts/multiotsu.py"
 
 
+rule compute_subject_histogram:
+    """Compute intensity histogram for a single subject for group-level Otsu thresholding.
+
+    Calculates a percentile-clipped histogram of the bias-field corrected image
+    and saves it as an NPZ file (hist_counts + bin_edges).  These per-subject
+    histogram files are later aggregated by the ``group_otsu`` rule to derive
+    a single set of thresholds that can be applied consistently across the whole
+    cohort with the ``multiotsu_group`` rule.
+    """
+    input:
+        corrected=bids(
+            root=root,
+            datatype="seg",
+            stain="{stain}",
+            level="{level}",
+            desc="corrected{method}".format(method=config["correction_method"]),
+            suffix="SPIM.ozx",
+            **inputs["spim"].wildcards,
+        ),
+    params:
+        hist_bin_width=float(config["seg_hist_bin_width"]),
+        hist_percentile_range=[float(x) for x in config["seg_hist_percentile_range"]],
+        zarrnii_kwargs={"orientation": config["orientation"]},
+    output:
+        histogram_npz=bids(
+            root=work,
+            datatype="seg",
+            stain="{stain}",
+            level="{level}",
+            desc="groupotsu+k{k,[0-9]+}i{i,[0-9]+}",
+            suffix="histogram.npz",
+            **inputs["spim"].wildcards,
+        ),
+    threads: 128 if config["dask_scheduler"] == "distributed" else 32
+    resources:
+        mem_mb=500000 if config["dask_scheduler"] == "distributed" else 250000,
+        runtime=90,
+    script:
+        "../scripts/compute_subject_histogram.py"
+
+
+rule multiotsu_group:
+    """Apply group-level Otsu threshold for segmentation.
+
+    Uses a pre-computed group-level Otsu threshold (derived from an aggregate
+    histogram across all subjects) to create a binary mask for the current
+    subject.  This ensures consistent thresholding across the whole cohort.
+    A per-subject PNG is also produced showing the group threshold overlaid
+    on this subject's own intensity histogram for visual quality control.
+
+    Run ``all_group_otsu`` before using this rule so that the group threshold
+    JSON file is available.
+    """
+    input:
+        corrected=bids(
+            root=root,
+            datatype="seg",
+            stain="{stain}",
+            level="{level}",
+            desc="corrected{method}".format(method=config["correction_method"]),
+            suffix="SPIM.ozx",
+            **inputs["spim"].wildcards,
+        ),
+        thresholds_json=bids(
+            root=root,
+            datatype="group",
+            stain="{stain}",
+            level="{level}",
+            desc="groupotsu+k{k,[0-9]+}i{i,[0-9]+}",
+            suffix="thresholds.json",
+        ),
+    params:
+        hist_bin_width=float(config["seg_hist_bin_width"]),
+        hist_percentile_range=[float(x) for x in config["seg_hist_percentile_range"]],
+        zarrnii_kwargs={"orientation": config["orientation"]},
+    output:
+        mask=bids(
+            root=root,
+            datatype="seg",
+            stain="{stain}",
+            level="{level}",
+            desc="groupotsu+k{k,[0-9]+}i{i,[0-9]+}",
+            suffix="mask.ozx",
+            **inputs["spim"].wildcards,
+        ),
+        thresholds_png=bids(
+            root=root,
+            datatype="seg",
+            stain="{stain}",
+            level="{level}",
+            desc="groupotsu+k{k,[0-9]+}i{i,[0-9]+}",
+            suffix="thresholds.png",
+            **inputs["spim"].wildcards,
+        ),
+    threads: 128 if config["dask_scheduler"] == "distributed" else 32
+    resources:
+        mem_mb=500000 if config["dask_scheduler"] == "distributed" else 250000,
+        disk_mb=2097152,
+        runtime=180,
+    script:
+        "../scripts/multiotsu_group.py"
+
+
 rule threshold:
     """Apply simple intensity threshold for segmentation.
 
@@ -160,12 +258,12 @@ rule threshold:
     """
     input:
         corrected=bids(
-            root=work,
+            root=root,
             datatype="seg",
             stain="{stain}",
             level="{level}",
             desc="corrected{method}".format(method=config["correction_method"]),
-            suffix="SPIM.ome.zarr",
+            suffix="SPIM.ozx",
             **inputs["spim"].wildcards,
         ),
     params: