Merge pull request #146 from khanlab/copilot/update-multiotsu-options

Replace fixed histogram range/bins in multiotsu with percentile-based range and bin-width; add threshold sweep QC report

Author: akhanf

spimquant/config/snakebids.yml

@@ -142,16 +142,18 @@ parse_args:
       - otsu+k3i2
     nargs: '+'
 
-  --seg_hist_range:
-    help: "Range of intensities to use for histogram calculation in multiotsu segmentation. Only applicable when seg_method is otsu+k{}i{}. Specify 2 numbers, for min and max values. (default: %(default)s)"
-    default: 
-      - 0
-      - 1000
+  --seg_hist_percentile_range:
+    help: "Percentile range to use for histogram calculation in multiotsu segmentation. Only applicable when seg_method is otsu+k{}i{}. Specify 2 numbers for the low and high percentiles (default: %(default)s)"
+    default:
+      - 1
+      - 99
     nargs: 2
+    type: float
 
-  --seg_hist_bins:
-    help: "Number of bins to use for histogram calculation in multiotsu segmentation. Only applicable when seg_method is otsu+k{}i{}. (default: %(default)s)"
-    default: 1000
+  --seg_hist_bin_width:
+    help: "Bin width to use for histogram calculation in multiotsu segmentation. Only applicable when seg_method is otsu+k{}i{}. (default: %(default)s)"
+    default: 1
+    type: float
 
 
   --register_to_mri:

spimquant/workflow/Snakefile

@@ -82,6 +82,9 @@ if stain_for_reg is None:
 
 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+")]
+
 if len(stains_for_seg) == 0 or config["no_segmentation"]:
     do_seg = False
     do_coloc = False
@@ -223,6 +226,32 @@ rule all_fieldfrac_tune:
         ),
 
 
+rule all_otsu_hist_qc:
+    """Target rule to generate otsu threshold sweep QC HTML reports.
+
+Produces one HTML report per subject/stain/method combination for every
+otsu+k{}i{} segmentation method.  Each report shows the multi-Otsu
+histogram alongside 2D crops at a sweep of threshold values so that the
+optimal threshold can be identified visually.  Only meaningful when the
+active ``seg_method`` contains at least one ``otsu+k{}i{}`` entry.
+"""
+    input:
+        inputs["spim"].expand(
+            bids(
+                root=root,
+                datatype="qc",
+                stain="{stain}",
+                desc="{desc}",
+                suffix="otsuthreshqc.html",
+                **inputs["spim"].wildcards,
+            ),
+            stain=stains_for_seg,
+            desc=otsu_seg_methods,
+        )
+        if (do_seg and len(otsu_seg_methods) > 0)
+        else [],
+
+
 rule all_register:
     input:
         inputs["spim"].expand(

spimquant/workflow/rules/qc.smk

@@ -403,6 +403,62 @@ specified by the ``{stain}`` wildcard.
         "../scripts/qc_objectstats.py"
 
 
+rule qc_otsu_threshold_sweep:
+    """Threshold sweep QC HTML report for multiotsu segmentation.
+
+Sweeps over a range of threshold values (spanning the configurable
+percentile range of the bias-field corrected image) and generates 2D
+crops at evenly-spaced positions in the image, one figure per threshold
+value.  The otsu histogram PNG produced by the ``multiotsu`` rule is
+embedded at the top of the report.  The resulting HTML report can be
+visually assessed to select the optimal threshold before running the full
+segmentation pipeline.
+
+Only applicable when ``seg_method`` uses the ``otsu+k{}i{}`` pattern.
+"""
+    input:
+        corrected=bids(
+            root=work,
+            datatype="seg",
+            stain="{stain}",
+            level=str(config["segmentation_level"]),
+            desc="corrected{method}".format(method=config["correction_method"]),
+            suffix="SPIM.ome.zarr",
+            **inputs["spim"].wildcards,
+        ),
+        thresholds_png=bids(
+            root=root,
+            datatype="seg",
+            stain="{stain}",
+            level=str(config["segmentation_level"]),
+            desc="{desc}",
+            suffix="thresholds.png",
+            **inputs["spim"].wildcards,
+        ),
+    output:
+        html=bids(
+            root=root,
+            datatype="qc",
+            stain="{stain}",
+            desc="{desc}",
+            suffix="otsuthreshqc.html",
+            **inputs["spim"].wildcards,
+        ),
+    threads: 4
+    resources:
+        mem_mb=32000,
+        runtime=30,
+    params:
+        n_thresholds=10,
+        n_crops=5,
+        patch_size=300,
+        level=config["segmentation_level"],
+        hist_percentile_range=[float(x) for x in config["seg_hist_percentile_range"]],
+        zarrnii_kwargs={"orientation": config["orientation"]},
+    script:
+        "../scripts/qc_otsu_threshold_sweep.py"
+
+
 rule qc_roi_summary:
     """Per-ROI summary QC: top-region bar plots for a single subject.
 

spimquant/workflow/rules/segmentation.smk

@@ -119,8 +119,8 @@ rule multiotsu:
             **inputs["spim"].wildcards,
         ),
     params:
-        hist_bins=int(config["seg_hist_bins"]),
-        hist_range=[int(x) for x in config["seg_hist_range"]],
+        hist_bin_width=float(config["seg_hist_bin_width"]),
+        hist_percentile_range=[float(x) for x in config["seg_hist_percentile_range"]],
         otsu_k=lambda wildcards: int(wildcards.k),
         otsu_threshold_index=lambda wildcards: int(wildcards.i),
         zarrnii_kwargs={"orientation": config["orientation"]},

spimquant/workflow/scripts/multiotsu.py

@@ -1,3 +1,5 @@
+import numpy as np
+
 from dask_setup import get_dask_client
 from zarrnii import ZarrNii
 from zarrnii.analysis import compute_otsu_thresholds
@@ -8,15 +10,45 @@
 if __name__ == "__main__":
     with get_dask_client(snakemake.config["dask_scheduler"], snakemake.threads):
 
-        # we use the default level=0, since we are reading in the n4 output, which is already downsampled if level was >0
-        znimg = ZarrNii.from_ome_zarr(
-            snakemake.input.corrected, **snakemake.params.zarrnii_kwargs
+        zarrnii_kwargs = snakemake.params.zarrnii_kwargs
+        pct_lo, pct_hi = snakemake.params.hist_percentile_range
+        bin_width = snakemake.params.hist_bin_width
+
+        # load a downsampled version to estimate the percentile-based range
+        print("estimating intensity range from downsampled image...")
+        znimg_ds = None
+        for ds_level in [5, 4, 3, 2, 1]:
+            try:
+                candidate = ZarrNii.from_ome_zarr(
+                    snakemake.input.corrected, level=ds_level, **zarrnii_kwargs
+                )
+                znimg_ds = candidate
+                break
+            except Exception:
+                pass
+
+        if znimg_ds is None:
+            znimg_ds = ZarrNii.from_ome_zarr(
+                snakemake.input.corrected, **zarrnii_kwargs
+            )
+
+        data_ds = znimg_ds.data.compute().ravel().astype(np.float32)
+        range_lo = float(np.percentile(data_ds, pct_lo))
+        range_hi = float(np.percentile(data_ds, pct_hi))
+        print(
+            f"  📊 percentile range [{pct_lo}%, {pct_hi}%]: [{range_lo:.3f}, {range_hi:.3f}]"
         )
 
-        # first calculate histogram - using preset bins to avoid issues where bins are too large
-        # because of high intensity outliers
+        # compute number of bins from bin width
+        n_bins = max(2, int(np.ceil((range_hi - range_lo) / bin_width)))
+        print(f"  📊 bins: {n_bins} (bin width: {bin_width})")
+
+        # we use the default level=0, since we are reading in the n4 output, which is already downsampled if level was >0
+        znimg = ZarrNii.from_ome_zarr(snakemake.input.corrected, **zarrnii_kwargs)
+
+        # calculate histogram using percentile-based range and bin-width-derived bin count
         (hist_counts, bin_edges) = znimg.compute_histogram(
-            bins=snakemake.params.hist_bins, range=snakemake.params.hist_range
+            bins=n_bins, range=[range_lo, range_hi]
         )
 
         # get otsu thresholds (uses histogram)