From a8692fd288fd079982b5efae4ad1c2ddd4d5205b Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 10 Apr 2026 23:52:17 +0000 Subject: [PATCH 1/3] Initial plan From 47f1c2b2766b5d0a99529fbe1fe8b0a1e12d33a9 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sat, 11 Apr 2026 00:03:25 +0000 Subject: [PATCH 2/3] Add batch histogram computation, aggregation, Otsu threshold, and QC rules Agent-Logs-Url: https://github.com/khanlab/SPIMquant/sessions/1911c21b-2dbc-4019-96d4-4c15c24c890a Co-authored-by: akhanf <11492701+akhanf@users.noreply.github.com> --- spimquant/workflow/Snakefile | 26 ++ spimquant/workflow/rules/qc.smk | 56 ++++ spimquant/workflow/rules/segmentation.smk | 191 +++++++++--- .../workflow/scripts/aggregate_histograms.py | 59 ++++ .../scripts/compute_spim_histogram.py | 51 ++++ .../scripts/multiotsu_from_batch_histogram.py | 76 +++++ .../scripts/qc_batch_otsu_threshold_sweep.py | 274 ++++++++++++++++++ 7 files changed, 696 insertions(+), 37 deletions(-) create mode 100644 spimquant/workflow/scripts/aggregate_histograms.py create mode 100644 spimquant/workflow/scripts/compute_spim_histogram.py create mode 100644 spimquant/workflow/scripts/multiotsu_from_batch_histogram.py create mode 100644 spimquant/workflow/scripts/qc_batch_otsu_threshold_sweep.py diff --git a/spimquant/workflow/Snakefile b/spimquant/workflow/Snakefile index 0a123fa..e19f13f 100644 --- a/spimquant/workflow/Snakefile +++ b/spimquant/workflow/Snakefile @@ -252,6 +252,32 @@ active ``seg_method`` contains at least one ``otsu+k{}i{}`` entry. else [], +rule all_batch_otsu_hist_qc: + """Target rule to generate batch-wide Otsu threshold sweep QC reports. + +Produces one HTML report per stain/method combination (across all subjects) +for every ``otsu+k{}i{}`` segmentation method. Each report shows the +aggregated histogram with Multi-Otsu threshold positions and one mid-volume +crop per subject at each threshold, enabling visual selection of a threshold +that generalises across the entire acquisition batch. Only meaningful when +the active ``seg_method`` contains at least one ``otsu+k{}i{}`` entry. +""" + input: + expand( + bids( + root=root, + datatype="group", + stain="{stain}", + desc="{desc}", + suffix="batchotsuthreshqc.html", + ), + 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( diff --git a/spimquant/workflow/rules/qc.smk b/spimquant/workflow/rules/qc.smk index 92ca884..1d668c1 100644 --- a/spimquant/workflow/rules/qc.smk +++ b/spimquant/workflow/rules/qc.smk @@ -459,6 +459,62 @@ Only applicable when ``seg_method`` uses the ``otsu+k{}i{}`` pattern. "../scripts/qc_otsu_threshold_sweep.py" +rule qc_batch_otsu_threshold_sweep: + """Batch threshold sweep QC HTML report for multiotsu segmentation. + +Uses the batch-wide Multi-Otsu thresholds (computed from the aggregated +histogram across all subjects in the dataset) to sweep over a range of +threshold values and display one mid-volume 2D crop per subject at each +threshold. Subjects are shown as columns so that the user can immediately +see whether a given threshold generalises across the entire acquisition batch. + +The aggregated histogram PNG is embedded at the top of the report. + +Only applicable when ``seg_method`` uses the ``otsu+k{}i{}`` pattern. +""" + input: + corrected_zarrs=inputs["spim"].expand( + 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, + ), + allow_missing=True, + ), + thresholds_png=bids( + root=root, + datatype="group", + stain="{stain}", + level=str(config["segmentation_level"]), + desc="{desc}", + suffix="thresholds.png", + ), + output: + html=bids( + root=root, + datatype="group", + stain="{stain}", + desc="{desc}", + suffix="batchotsuthreshqc.html", + ), + threads: 4 + resources: + mem_mb=32000, + runtime=60, + params: + n_thresholds=10, + 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_batch_otsu_threshold_sweep.py" + + rule qc_roi_summary: """Per-ROI summary QC: top-region bar plots for a single subject. diff --git a/spimquant/workflow/rules/segmentation.smk b/spimquant/workflow/rules/segmentation.smk index 70904b6..1cfff63 100644 --- a/spimquant/workflow/rules/segmentation.smk +++ b/spimquant/workflow/rules/segmentation.smk @@ -1,7 +1,7 @@ """ Segmentation and quantification workflow for SPIMquant. -This module performs intensity correction, segmentation, and quantitative analysis of +This module performs intensity correction, segmentation, and quantitative analysis of pathology markers from SPIM microscopy data. It handles multi-channel data with different stains (e.g., beta-amyloid, alpha-synuclein, Iba1) and produces per-region statistics. @@ -27,9 +27,6 @@ rule gaussian_biasfield: """simple bias field correction with gaussian""" input: spim=inputs["spim"].path, - params: - proc_level=5, - zarrnii_kwargs={"orientation": config["orientation"]}, output: corrected=temp( directory( @@ -64,6 +61,9 @@ rule gaussian_biasfield: mem_mb=256000, disk_mb=2097152, runtime=15, + params: + proc_level=5, + zarrnii_kwargs={"orientation": config["orientation"]}, script: "../scripts/gaussian_biasfield.py" @@ -72,11 +72,6 @@ rule n4_biasfield: """N4 bias field correction with antspyx""" input: spim=inputs["spim"].path, - params: - proc_level=5, - zarrnii_kwargs={"orientation": config["orientation"]}, - 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( @@ -96,18 +91,23 @@ rule n4_biasfield: resources: mem_mb=500000 if config["dask_scheduler"] == "distributed" else 250000, runtime=180, + params: + proc_level=5, + zarrnii_kwargs={"orientation": config["orientation"]}, + shrink_factor=16 if config["sloppy"] else 1, + target_chunk_size=512, #this sets the chunk size for this and downstream masks script: "../scripts/n4_biasfield.py" rule multiotsu: """Perform multi-Otsu thresholding for segmentation. - - Applies multi-level Otsu thresholding to identify intensity classes in the - corrected image. The k parameter determines number of classes, and the i parameter - selects which threshold index to use for creating the binary mask. Outputs a - histogram visualization of the threshold selection. - """ + +Applies multi-level Otsu thresholding to identify intensity classes in the +corrected image. The k parameter determines number of classes, and the i parameter +selects which threshold index to use for creating the binary mask. Outputs a +histogram visualization of the threshold selection. +""" input: corrected=bids( root=work, @@ -118,12 +118,6 @@ rule multiotsu: suffix="SPIM.ome.zarr", **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"]], - otsu_k=lambda wildcards: int(wildcards.k), - otsu_threshold_index=lambda wildcards: int(wildcards.i), - zarrnii_kwargs={"orientation": config["orientation"]}, output: mask=bids( root=root, @@ -148,16 +142,23 @@ rule multiotsu: mem_mb=500000 if config["dask_scheduler"] == "distributed" else 250000, disk_mb=2097152, runtime=180, + params: + 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"]}, script: "../scripts/multiotsu.py" rule threshold: """Apply simple intensity threshold for segmentation. - - Creates binary mask by thresholding the corrected image at a specified intensity value. - Simpler alternative to multi-Otsu for cases where the threshold is known a priori. - """ + +Creates binary mask by thresholding the corrected image at a specified intensity value. +Simpler alternative to multi-Otsu for cases where the threshold is known a priori. + +""" input: corrected=bids( root=work, @@ -168,9 +169,6 @@ rule threshold: suffix="SPIM.ome.zarr", **inputs["spim"].wildcards, ), - params: - threshold=lambda wildcards: int(wildcards.threshold), - zarrnii_kwargs={"orientation": config["orientation"]}, output: mask=bids( root=root, @@ -185,17 +183,136 @@ rule threshold: resources: mem_mb=500000 if config["dask_scheduler"] == "distributed" else 250000, runtime=180, + params: + threshold=lambda wildcards: int(wildcards.threshold), + zarrnii_kwargs={"orientation": config["orientation"]}, script: "../scripts/threshold.py" +rule compute_spim_histogram: + """Compute an intensity histogram from a raw SPIM channel at a given level. + +Reads the OME-Zarr at the specified pyramid *level* using Dask and accumulates +a histogram over a fixed intensity range with a fixed number of bins. The bin +edges and counts are stored as a compressed NumPy archive (``.npz``) for later +aggregation and batch-wide threshold optimisation. +""" + input: + spim=inputs["spim"].path, + output: + histogram=bids( + root=root, + datatype="seg", + stain="{stain}", + level="{level}", + suffix="histogram.npz", + **inputs["spim"].wildcards, + ), + threads: 32 + resources: + mem_mb=32000, + runtime=30, + params: + level="{level}", + hist_bins=config.get("spim_hist_bins", 500), + hist_range=config.get("spim_hist_range", [0, 65535]), + zarrnii_kwargs={"orientation": config["orientation"]}, + script: + "../scripts/compute_spim_histogram.py" + + +rule aggregate_histograms: + """Aggregate per-subject histograms into a single batch histogram. + +Collects histogram ``.npz`` files from all subjects (for a given stain and +pyramid level) and produces a single combined histogram by summing counts. +All input histograms must have been computed with the same ``hist_range`` and +``hist_bins`` parameters. +""" + input: + histograms=inputs["spim"].expand( + bids( + root=root, + datatype="seg", + stain="{stain}", + level="{level}", + suffix="histogram.npz", + **inputs["spim"].wildcards, + ), + allow_missing=True, + ), + output: + histogram=bids( + root=root, + datatype="group", + stain="{stain}", + level="{level}", + suffix="histogram.npz", + ), + threads: 1 + resources: + mem_mb=4000, + runtime=10, + script: + "../scripts/aggregate_histograms.py" + + +rule multiotsu_from_batch_histogram: + """Compute Multi-Otsu thresholds from the batch-aggregated histogram. + +Reads the aggregated histogram, optionally clips to a percentile-based +intensity range, and applies multi-level Otsu thresholding. This yields a +single set of thresholds that is consistent across all subjects in the batch +(i.e. acquisitions with the same imaging parameters). + +Outputs a TSV of threshold index → value and a histogram PNG with the +threshold positions marked. +""" + input: + histogram=bids( + root=root, + datatype="group", + stain="{stain}", + level="{level}", + suffix="histogram.npz", + ), + output: + thresholds_png=bids( + root=root, + datatype="group", + stain="{stain}", + level="{level}", + desc="otsu+k{k,[0-9]+}i{i,[0-9]+}", + suffix="thresholds.png", + ), + thresholds_tsv=bids( + root=root, + datatype="group", + stain="{stain}", + level="{level}", + desc="otsu+k{k,[0-9]+}i{i,[0-9]+}", + suffix="thresholds.tsv", + ), + threads: 1 + resources: + mem_mb=4000, + runtime=10, + params: + 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), + script: + "../scripts/multiotsu_from_batch_histogram.py" + + rule clean_segmentation: """Clean segmentation mask by removing edge artifacts and small objects. - - Performs connected component analysis to identify and exclude objects that - extend too close to the image boundaries (likely artifacts). Creates both - a cleaned mask and an exclusion mask showing what was removed. - """ + +Performs connected component analysis to identify and exclude objects that +extend too close to the image boundaries (likely artifacts). Creates both +a cleaned mask and an exclusion mask showing what was removed. +""" input: mask=bids( root=root, @@ -206,10 +323,6 @@ rule clean_segmentation: suffix="mask.ozx", **inputs["spim"].wildcards, ), - params: - max_extent=0.15, - proc_level=2, #level at which to calculate conncomp - zarrnii_kwargs={"orientation": config["orientation"]}, output: exclude_mask=bids( root=root, @@ -234,5 +347,9 @@ rule clean_segmentation: mem_mb=256000, disk_mb=2097152, runtime=30, + params: + max_extent=0.15, + proc_level=2, #level at which to calculate conncomp + zarrnii_kwargs={"orientation": config["orientation"]}, script: "../scripts/clean_segmentation.py" diff --git a/spimquant/workflow/scripts/aggregate_histograms.py b/spimquant/workflow/scripts/aggregate_histograms.py new file mode 100644 index 0000000..4e32237 --- /dev/null +++ b/spimquant/workflow/scripts/aggregate_histograms.py @@ -0,0 +1,59 @@ +"""Aggregate intensity histograms across subjects. + +Loads histogram ``.npz`` files (each containing ``counts`` and ``bin_edges``) +from multiple subjects and produces a single combined histogram by summing +counts. All input histograms must share the same bin edges (i.e. they must +have been computed with the same ``hist_range`` and ``hist_bins`` parameters). + +This is a Snakemake script; the ``snakemake`` object is automatically provided +when executed as part of a Snakemake workflow. +""" + +import numpy as np + + +def main(): + histogram_files = snakemake.input.histograms + + print(f"Aggregating {len(histogram_files)} histogram(s)") + + total_counts = None + ref_bin_edges = None + + for hist_file in histogram_files: + data = np.load(hist_file) + counts = data["counts"].astype(np.float64) + bin_edges = data["bin_edges"] + + if ref_bin_edges is None: + ref_bin_edges = bin_edges + total_counts = counts.copy() + else: + if counts.shape == total_counts.shape and np.allclose( + bin_edges, ref_bin_edges + ): + total_counts += counts + else: + # Rebin by treating bin centres as weighted samples + bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2.0 + rebinned, _ = np.histogram( + bin_centers, + bins=ref_bin_edges, + weights=counts, + ) + total_counts += rebinned + + print(f" {hist_file}: total voxels = {counts.sum():.0f}") + + print(f"Aggregated total voxels = {total_counts.sum():.0f}") + + np.savez_compressed( + snakemake.output.histogram, + counts=total_counts, + bin_edges=ref_bin_edges, + ) + print(f"Saved aggregated histogram to {snakemake.output.histogram}") + + +if __name__ == "__main__": + main() diff --git a/spimquant/workflow/scripts/compute_spim_histogram.py b/spimquant/workflow/scripts/compute_spim_histogram.py new file mode 100644 index 0000000..b6b9953 --- /dev/null +++ b/spimquant/workflow/scripts/compute_spim_histogram.py @@ -0,0 +1,51 @@ +"""Compute intensity histogram from a raw SPIM dataset channel. + +Reads the raw OME-Zarr at a given pyramid level using Dask and computes an +intensity histogram over a fixed range. The resulting bin edges and counts +are saved as a compressed NumPy archive (.npz) for downstream aggregation +and threshold computation. + +This is a Snakemake script; the ``snakemake`` object is automatically provided +when executed as part of a Snakemake workflow. +""" + +import numpy as np +from dask_setup import get_dask_client +from zarrnii import ZarrNii + + +def main(): + level = snakemake.params.level + stain = snakemake.wildcards.stain + hist_bins = snakemake.params.hist_bins + hist_range = snakemake.params.hist_range + + print(f"Computing histogram for stain={stain} at level={level}") + print(f" hist_bins={hist_bins}, hist_range={hist_range}") + + with get_dask_client(snakemake.config["dask_scheduler"], snakemake.threads): + znimg = ZarrNii.from_ome_zarr( + snakemake.input.spim, + level=level, + channel_labels=[stain], + **snakemake.params.zarrnii_kwargs, + ) + + hist_counts, bin_edges = znimg.compute_histogram( + bins=hist_bins, + range=hist_range, + ) + + hist_counts = np.asarray(hist_counts, dtype=np.float64) + bin_edges = np.asarray(bin_edges, dtype=np.float64) + + np.savez_compressed( + snakemake.output.histogram, + counts=hist_counts, + bin_edges=bin_edges, + ) + print(f"Saved histogram to {snakemake.output.histogram}") + + +if __name__ == "__main__": + main() diff --git a/spimquant/workflow/scripts/multiotsu_from_batch_histogram.py b/spimquant/workflow/scripts/multiotsu_from_batch_histogram.py new file mode 100644 index 0000000..0f642bc --- /dev/null +++ b/spimquant/workflow/scripts/multiotsu_from_batch_histogram.py @@ -0,0 +1,76 @@ +"""Compute Multi-Otsu thresholds from a batch-aggregated histogram. + +Reads the aggregated histogram ``.npz`` file (containing ``counts`` and +``bin_edges``), optionally filters to a percentile-based intensity range to +exclude background and saturation artefacts, and then applies multi-level Otsu +thresholding. + +Outputs: +- A PNG visualisation of the histogram with Otsu threshold positions marked. +- A TSV file listing every threshold index and the corresponding value. + +This is a Snakemake script; the ``snakemake`` object is automatically provided +when executed as part of a Snakemake workflow. +""" + +import matplotlib +import numpy as np + +matplotlib.use("agg") + +from zarrnii.analysis import compute_otsu_thresholds + + +def main(): + data = np.load(snakemake.input.histogram) + hist_counts = data["counts"].astype(np.float64) + bin_edges = data["bin_edges"] + + bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2.0 + + pct_lo, pct_hi = snakemake.params.hist_percentile_range + + # Derive a percentile-based intensity range from the aggregated histogram + total = hist_counts.sum() + if total > 0: + cumsum_norm = np.cumsum(hist_counts) / total + idx_lo = int(np.searchsorted(cumsum_norm, pct_lo / 100.0)) + idx_hi = int(np.searchsorted(cumsum_norm, pct_hi / 100.0)) + idx_lo = max(0, min(idx_lo, len(bin_centers) - 1)) + idx_hi = max(idx_lo + 1, min(idx_hi, len(bin_centers) - 1)) + else: + idx_lo, idx_hi = 0, len(bin_centers) - 1 + + range_lo = float(bin_centers[idx_lo]) + range_hi = float(bin_centers[idx_hi]) + print( + f"Percentile range [{pct_lo}%–{pct_hi}%]: " f"[{range_lo:.3f}, {range_hi:.3f}]" + ) + + # Slice the histogram to the percentile-based range + filtered_counts = hist_counts[idx_lo : idx_hi + 1] + filtered_edges = bin_edges[idx_lo : idx_hi + 2] + + # Compute Multi-Otsu thresholds + thresholds, fig = compute_otsu_thresholds( + filtered_counts, + classes=snakemake.params.otsu_k, + bin_edges=filtered_edges, + return_figure=True, + ) + print(f"Batch Otsu thresholds: {[f'{t:.3f}' for t in thresholds]}") + + fig.savefig(snakemake.output.thresholds_png) + + # Save thresholds to TSV + with open(snakemake.output.thresholds_tsv, "w") as fh: + fh.write("threshold_index\tthreshold_value\n") + for idx, val in enumerate(thresholds): + fh.write(f"{idx}\t{val:.6f}\n") + + print(f"Saved thresholds to {snakemake.output.thresholds_tsv}") + print(f"Saved histogram PNG to {snakemake.output.thresholds_png}") + + +if __name__ == "__main__": + main() diff --git a/spimquant/workflow/scripts/qc_batch_otsu_threshold_sweep.py b/spimquant/workflow/scripts/qc_batch_otsu_threshold_sweep.py new file mode 100644 index 0000000..fa71bfa --- /dev/null +++ b/spimquant/workflow/scripts/qc_batch_otsu_threshold_sweep.py @@ -0,0 +1,274 @@ +"""Batch Otsu threshold sweep QC report. + +For a given stain and multiotsu segmentation method, uses the batch-wide +Multi-Otsu thresholds (computed from the aggregated histogram across all +subjects) to show one mid-volume 2D crop per subject at each threshold +value. This enables visual assessment of whether the same threshold works +across an entire acquisition batch. + +The aggregated histogram PNG is embedded at the top of the report, followed +by a sweep of threshold values, each showing one patch per subject side by +side. + +This is a Snakemake script; the ``snakemake`` object is automatically provided +when executed as part of a Snakemake workflow. +""" + +import base64 +import re +from io import BytesIO +from pathlib import Path + +import matplotlib + +matplotlib.use("agg") +import matplotlib.pyplot as plt +import numpy as np +from zarrnii import ZarrNii + + +def _fig_to_base64(fig): + """Convert a matplotlib figure to a base64-encoded PNG data URI.""" + buf = BytesIO() + fig.savefig(buf, format="png", dpi=100, bbox_inches="tight") + buf.seek(0) + b64 = base64.b64encode(buf.read()).decode("utf-8") + plt.close(fig) + return f"data:image/png;base64,{b64}" + + +def _file_to_base64_png(path): + """Read a PNG file and return a base64 data URI.""" + with open(path, "rb") as fh: + b64 = base64.b64encode(fh.read()).decode("utf-8") + return f"data:image/png;base64,{b64}" + + +def _norm(arr, lo, hi): + """Linearly normalise *arr* to [0, 1] using the given bounds.""" + if hi > lo: + return np.clip((arr.astype(np.float32) - lo) / (hi - lo), 0.0, 1.0) + return np.zeros_like(arr, dtype=np.float32) + + +def _subject_label(zarr_path): + """Extract a short subject label from a file path.""" + m = re.search(r"sub-([^_/\\]+)", zarr_path) + if m: + return f"sub-{m.group(1)}" + return Path(zarr_path).stem + + +def main(): + zarrnii_kwargs = snakemake.params.zarrnii_kwargs + n_thresholds = snakemake.params.n_thresholds + patch_size = snakemake.params.patch_size + level = snakemake.params.level + pct_lo, pct_hi = snakemake.params.hist_percentile_range + + stain = snakemake.wildcards.stain + desc = snakemake.wildcards.desc + + corrected_zarrs = snakemake.input.corrected_zarrs + n_subjects = len(corrected_zarrs) + subject_labels = [_subject_label(p) for p in corrected_zarrs] + + print(f"Batch threshold sweep QC: stain={stain}, desc={desc}") + print(f" {n_subjects} subject(s): {subject_labels}") + + # ------------------------------------------------------------------ + # Load one mid-volume crop per subject + # ------------------------------------------------------------------ + half_patch = patch_size // 2 + subject_crops = [] + subject_ranges = [] + + for zarr_path, label in zip(corrected_zarrs, subject_labels): + print(f" Loading crop from {label} …") + znimg = None + for ds_offset in [3, 2, 1, 0]: + try: + candidate = ZarrNii.from_ome_zarr( + zarr_path, level=level + ds_offset, **zarrnii_kwargs + ) + znimg = candidate + print(f" Loaded at pyramid level {level + ds_offset}") + break + except Exception as exc: + print(f" Level {level + ds_offset} not available: {exc}") + + if znimg is None: + znimg = ZarrNii.from_ome_zarr(zarr_path, **zarrnii_kwargs) + + data = znimg.data.compute() + if data.ndim == 4: + data = data[0] # drop channel dim → (Z, Y, X) + + flat = data.ravel().astype(np.float32) + range_lo = float(np.percentile(flat, pct_lo)) + range_hi = float(np.percentile(flat, pct_hi)) + subject_ranges.append((range_lo, range_hi)) + + Z, Y, X = data.shape + z_mid = Z // 2 + y_pos = Y // 2 + x_pos = X // 2 + y0 = max(0, y_pos - half_patch) + y1 = min(Y, y_pos + half_patch) + x0 = max(0, x_pos - half_patch) + x1 = min(X, x_pos + half_patch) + subject_crops.append(data[z_mid, y0:y1, x0:x1]) + + # ------------------------------------------------------------------ + # Determine global sweep range and threshold values + # ------------------------------------------------------------------ + global_lo = float(np.min([r[0] for r in subject_ranges])) + global_hi = float(np.max([r[1] for r in subject_ranges])) + thresholds = np.linspace(global_lo, global_hi, n_thresholds) + print(f" Sweep range: [{global_lo:.3f}, {global_hi:.3f}]") + + # ------------------------------------------------------------------ + # Build one figure per threshold: columns = subjects + # ------------------------------------------------------------------ + threshold_entries = [] + for thresh in thresholds: + fig, axes = plt.subplots( + 1, n_subjects, figsize=(max(3 * n_subjects, 6), 3), squeeze=False + ) + axes = axes[0] + fig.suptitle(f"Threshold = {thresh:.1f}", fontsize=10, fontweight="bold") + + for ax, crop, (rlo, rhi), label in zip( + axes, subject_crops, subject_ranges, subject_labels + ): + crop_norm = _norm(crop, rlo, rhi) + mask_crop = (crop > thresh).astype(np.float32) + ax.imshow(crop_norm, cmap="gray") + mask_ma = np.ma.masked_where(mask_crop < 0.5, mask_crop) + ax.imshow(mask_ma, cmap="Reds", alpha=0.6, vmin=0, vmax=1) + ax.set_xticks([]) + ax.set_yticks([]) + ax.set_title(label, fontsize=7) + + plt.tight_layout() + threshold_entries.append( + {"threshold": float(thresh), "b64": _fig_to_base64(fig)} + ) + + # ------------------------------------------------------------------ + # Embed the batch Otsu histogram PNG + # ------------------------------------------------------------------ + otsu_hist_b64 = _file_to_base64_png(snakemake.input.thresholds_png) + + # ------------------------------------------------------------------ + # Generate HTML report + # ------------------------------------------------------------------ + sweep_html_parts = [] + for entry in threshold_entries: + thresh_label = f"{entry['threshold']:.1f}" + sweep_html_parts.append( + f'
' + f"

Threshold = {thresh_label}

" + f'threshold {thresh_label}' + f"
" + ) + sweep_html = "\n".join(sweep_html_parts) + + subjects_str = ", ".join(subject_labels) + + html = f""" + + + + + Batch Otsu Threshold Sweep QC – {stain} / {desc} + + + +
+

Batch Otsu Threshold Sweep QC Report

+ +
+

Stain: {stain}

+

Method: {desc}

+

Subjects ({n_subjects}): {subjects_str}

+

Intensity sweep range + ({pct_lo}th–{pct_hi}th percentile): + {global_lo:.1f} – {global_hi:.1f}

+

Number of threshold values: {n_thresholds}

+
+ +

Batch Histogram (Multi-level Otsu Thresholds from Aggregated Data)

+

The histogram below was computed from the aggregated + intensity data across all subjects in this batch. The vertical lines + show the Multi-Otsu threshold positions. These batch-wide thresholds + are used as a fixed reference for all subjects so that segmentation is + consistent across acquisitions.

+
+ Batch Otsu Histogram +
+ +

Threshold Sweep (one crop per subject)

+

Each row corresponds to one threshold value. Within each row, every + column shows a mid-volume 2D crop from one subject with the resulting + binary mask overlaid in red. Use this to identify a threshold that + reliably captures the signal of interest across all subjects + without excessive background segmentation.

+ + {sweep_html} + +
+ +""" + + with open(snakemake.output.html, "w") as fh: + fh.write(html) + + print(f"Saved batch threshold sweep QC to {snakemake.output.html}") + + +if __name__ == "__main__": + main() From f5efd2b0f860f8e99919342877eef81d7f027f99 Mon Sep 17 00:00:00 2001 From: Ali Khan Date: Mon, 20 Apr 2026 21:23:30 -0400 Subject: [PATCH 3/3] updates to run on non-corrected image --- spimquant/workflow/Snakefile | 27 +++++++++++ spimquant/workflow/rules/qc.smk | 47 ++++++++++++++++++- spimquant/workflow/rules/segmentation.smk | 9 ++-- .../scripts/qc_batch_otsu_threshold_sweep.py | 10 ++-- 4 files changed, 83 insertions(+), 10 deletions(-) diff --git a/spimquant/workflow/Snakefile b/spimquant/workflow/Snakefile index e19f13f..9a9931b 100644 --- a/spimquant/workflow/Snakefile +++ b/spimquant/workflow/Snakefile @@ -278,6 +278,33 @@ the active ``seg_method`` contains at least one ``otsu+k{}i{}`` entry. else [], +rule all_batch_otsu_hist_qc_raw: + """Target rule to generate batch-wide Otsu threshold sweep QC reports. + +Produces one HTML report per stain/method combination (across all subjects) +for every ``otsu+k{}i{}`` segmentation method. Each report shows the +aggregated histogram with Multi-Otsu threshold positions and one mid-volume +crop per subject at each threshold, enabling visual selection of a threshold +that generalises across the entire acquisition batch. Only meaningful when +the active ``seg_method`` contains at least one ``otsu+k{}i{}`` entry. +""" + input: + expand( + bids( + root=root, + datatype="group", + stain="{stain}", + desc="{desc}", + suffix="batchotsuthreshqcraw.html", + ), + 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( diff --git a/spimquant/workflow/rules/qc.smk b/spimquant/workflow/rules/qc.smk index 1d668c1..e50e5c0 100644 --- a/spimquant/workflow/rules/qc.smk +++ b/spimquant/workflow/rules/qc.smk @@ -458,6 +458,51 @@ Only applicable when ``seg_method`` uses the ``otsu+k{}i{}`` pattern. script: "../scripts/qc_otsu_threshold_sweep.py" +rule qc_batch_otsu_threshold_sweep_raw: + """Batch threshold sweep QC HTML report for multiotsu segmentation, on raw uncorrected images. + +Uses the batch-wide Multi-Otsu thresholds (computed from the aggregated +histogram across all subjects in the dataset) to sweep over a range of +threshold values and display one mid-volume 2D crop per subject at each +threshold. Subjects are shown as columns so that the user can immediately +see whether a given threshold generalises across the entire acquisition batch. + +The aggregated histogram PNG is embedded at the top of the report. + +Only applicable when ``seg_method`` uses the ``otsu+k{}i{}`` pattern. +""" + input: + zarrs=inputs["spim"].expand(), + thresholds_png=bids( + root=root, + datatype="group", + stain="{stain}", + level=str(config["segmentation_level"]), + desc="{desc}", + suffix="thresholds.png", + ), + output: + html=bids( + root=root, + datatype="group", + stain="{stain}", + desc="{desc}", + suffix="batchotsuthreshqcraw.html", + ), + threads: 4 + resources: + mem_mb=32000, + runtime=60, + params: + n_thresholds=10, + patch_size=300, + level=config["segmentation_level"], + hist_percentile_range=[float(x) for x in config["seg_hist_percentile_range"]], + zarrnii_kwargs=lambda wildcards: {"orientation": config["orientation"],"channel_labels":[wildcards.stain]}, + script: + "../scripts/qc_batch_otsu_threshold_sweep.py" + + rule qc_batch_otsu_threshold_sweep: """Batch threshold sweep QC HTML report for multiotsu segmentation. @@ -473,7 +518,7 @@ The aggregated histogram PNG is embedded at the top of the report. Only applicable when ``seg_method`` uses the ``otsu+k{}i{}`` pattern. """ input: - corrected_zarrs=inputs["spim"].expand( + zarrs=inputs["spim"].expand( bids( root=work, datatype="seg", diff --git a/spimquant/workflow/rules/segmentation.smk b/spimquant/workflow/rules/segmentation.smk index 1cfff63..6dd067f 100644 --- a/spimquant/workflow/rules/segmentation.smk +++ b/spimquant/workflow/rules/segmentation.smk @@ -209,12 +209,13 @@ aggregation and batch-wide threshold optimisation. suffix="histogram.npz", **inputs["spim"].wildcards, ), - threads: 32 + threads: 128 if config["dask_scheduler"] == "distributed" else 32 resources: - mem_mb=32000, - runtime=30, + mem_mb=500000 if config["dask_scheduler"] == "distributed" else 250000, + disk_mb=2097152, + runtime=180, params: - level="{level}", + level=lambda wildcards: int(wildcards.level), hist_bins=config.get("spim_hist_bins", 500), hist_range=config.get("spim_hist_range", [0, 65535]), zarrnii_kwargs={"orientation": config["orientation"]}, diff --git a/spimquant/workflow/scripts/qc_batch_otsu_threshold_sweep.py b/spimquant/workflow/scripts/qc_batch_otsu_threshold_sweep.py index fa71bfa..e09c70c 100644 --- a/spimquant/workflow/scripts/qc_batch_otsu_threshold_sweep.py +++ b/spimquant/workflow/scripts/qc_batch_otsu_threshold_sweep.py @@ -69,9 +69,9 @@ def main(): stain = snakemake.wildcards.stain desc = snakemake.wildcards.desc - corrected_zarrs = snakemake.input.corrected_zarrs - n_subjects = len(corrected_zarrs) - subject_labels = [_subject_label(p) for p in corrected_zarrs] + zarrs = snakemake.input.zarrs + n_subjects = len(zarrs) + subject_labels = [_subject_label(p) for p in zarrs] print(f"Batch threshold sweep QC: stain={stain}, desc={desc}") print(f" {n_subjects} subject(s): {subject_labels}") @@ -83,10 +83,10 @@ def main(): subject_crops = [] subject_ranges = [] - for zarr_path, label in zip(corrected_zarrs, subject_labels): + for zarr_path, label in zip(zarrs, subject_labels): print(f" Loading crop from {label} …") znimg = None - for ds_offset in [3, 2, 1, 0]: + for ds_offset in [5,4,3, 2, 1, 0]: try: candidate = ZarrNii.from_ome_zarr( zarr_path, level=level + ds_offset, **zarrnii_kwargs