updates to figures incl customization of montaging

right now only seg=coarse uses custom setting..

Author: akhanf

spimquant/workflow/rules/qc.smk

@@ -46,7 +46,7 @@ saturation/clip fraction (percentage of voxels at the maximum bin).
             suffix="histogram.png",
             **inputs["spim"].wildcards,
         ),
-    threads: 8
+    threads: 4
     resources:
         mem_mb=16000,
         runtime=30,
@@ -88,7 +88,7 @@ Aspect ratio is corrected using voxel spacings from ``ZarrNii.get_zooms()``.
             suffix="segslices.png",
             **inputs["spim"].wildcards,
         ),
-    threads: 8
+    threads: 4
     resources:
         mem_mb=16000,
         runtime=30,
@@ -127,7 +127,7 @@ for isotropic display and physically correct aspect ratio.
             suffix="vesselslices.png",
             **inputs["spim"].wildcards,
         ),
-    threads: 8
+    threads: 4
     resources:
         mem_mb=16000,
         runtime=30,
@@ -185,13 +185,14 @@ quality within individual brain regions.
             suffix="roimontage.png",
             **inputs["spim"].wildcards,
         ),
-    threads: 32
+    threads: 4
     resources:
         mem_mb=32000,
         runtime=15,
     params:
-        max_rois=25,
-        n_cols=5,
+        max_rois=lambda wildcards: 25 if wildcards.seg == "coarse" else 100,
+        n_cols=lambda wildcards: 5 if wildcards.seg == "coarse" else 10,
+        patch_size=lambda wildcards: 2000 if wildcards.seg == "coarse" else 500,
         level=config["segmentation_level"],
     script:
         "../scripts/qc_segmentation_roi_zoom.py"
@@ -241,13 +242,14 @@ ZarrNiiAtlas for atlas-based ROI cropping.
             suffix="vesselroimontage.png",
             **inputs["spim"].wildcards,
         ),
-    threads: 32
+    threads: 4
     resources:
         mem_mb=32000,
         runtime=15,
     params:
-        max_rois=25,
-        n_cols=5,
+        max_rois=lambda wildcards: 25 if wildcards.seg == "coarse" else 100,
+        n_cols=lambda wildcards: 5 if wildcards.seg == "coarse" else 10,
+        patch_size=lambda wildcards: 2000 if wildcards.seg == "coarse" else 500,
         level=config["segmentation_level"],
     script:
         "../scripts/qc_segmentation_roi_zoom.py"

spimquant/workflow/scripts/qc_intensity_histogram.py

@@ -24,7 +24,7 @@ def main():
     hist_bins = snakemake.params.hist_bins
     hist_range = snakemake.params.hist_range
 
-    with get_dask_client(snakemake.config["dask_scheduler"], snakemake.threads):
+    with get_dask_client("threads", snakemake.threads):
         znimg = ZarrNii.from_ome_zarr(
             snakemake.input.spim,
             level=level,

spimquant/workflow/scripts/qc_roi_summary.py

@@ -2,7 +2,7 @@
 
 Loads the merged segmentation-statistics TSV (all stains) together with
 the atlas label table, then produces bar-chart visualisations of the
-top brain-regions ranked by field fraction and count for each stain.
+top brain-regions ranked by field fraction and density for each stain.
 
 This is a Snakemake script that expects the ``snakemake`` object to be
 available, which is automatically provided when executed as part of a
@@ -19,7 +19,6 @@
 # Suffixes used to identify stain-prefixed metric columns in mergedsegstats TSV.
 # Columns follow the pattern "{stain}+{metric}", e.g. "Abeta+fieldfrac".
 _SUFFIX_FIELDFRAC = "+fieldfrac"
-_SUFFIX_COUNT = "+count"
 _SUFFIX_DENSITY = "+density"
 
 
@@ -62,15 +61,13 @@ def main():
 
     # Identify stain-prefixed metric columns (pattern: "{stain}+{metric}")
     ff_cols = [c for c in stats_df.columns if c.endswith(_SUFFIX_FIELDFRAC)]
-    count_cols = [c for c in stats_df.columns if c.endswith(_SUFFIX_COUNT)]
     density_cols = [c for c in stats_df.columns if c.endswith(_SUFFIX_DENSITY)]
 
-    # Determine number of rows: 1 row per metric type (ff, count, density)
+    # Determine number of rows: 1 row per metric type (ff,  density)
     # with one subplot per stain within each row
     n_ff = len(ff_cols)
-    n_count = len(count_cols)
     n_density = len(density_cols)
-    n_rows = (1 if n_ff else 0) + (1 if n_count else 0) + (1 if n_density else 0)
+    n_rows = (1 if n_ff else 0) + (1 if n_density else 0)
 
     if n_rows == 0:
         fig, ax = plt.subplots(figsize=(8, 4))
@@ -95,8 +92,6 @@ def main():
     row_specs = []
     if n_ff:
         row_specs.append(("Field Fraction (%)", ff_cols, "steelblue"))
-    if n_count:
-        row_specs.append(("Count (objects)", count_cols, "darkorange"))
     if n_density:
         row_specs.append(("Density (objects/vol)", density_cols, "forestgreen"))
 

spimquant/workflow/scripts/qc_segmentation_overview.py

@@ -113,8 +113,8 @@ def main():
         for col, sl in enumerate(slice_indices):
             idx = [slice(None)] * 3
             idx[ax_idx] = int(sl)
-            spim_sl = np.rot90(spim_norm[tuple(idx)])
-            mask_sl = np.rot90(mask_norm[tuple(idx)])
+            spim_sl = spim_norm[tuple(idx)]
+            mask_sl = mask_norm[tuple(idx)]
 
             ax = axes[row, col]
             ax.imshow(spim_sl, cmap="gray", vmin=0, vmax=1, aspect=aspect)

spimquant/workflow/scripts/qc_segmentation_roi_zoom.py

@@ -23,47 +23,23 @@
 from scipy.ndimage import zoom
 
 
-def _percentile_norm(arr, pct_low=1, pct_high=99):
-    """Percentile-normalise *arr* to the range [0, 1] using global statistics."""
+def _estimate_global_percentiles(
+    arr,
+    pct_low=1,
+    pct_high=99,
+):
+    """
+    Estimate percentile normalization bounds from an image
+    """
     lo = np.percentile(arr, pct_low)
     hi = np.percentile(arr, pct_high)
-    if hi > lo:
-        return np.clip((arr.astype(float) - lo) / (hi - lo), 0.0, 1.0)
-    return np.zeros_like(arr, dtype=float)
-
-
-def _match_shape(source, target_shape, order=1):
-    """Zoom *source* array to *target_shape* if shapes differ."""
-    if source.shape == target_shape:
-        return source
-    factors = [t / s for t, s in zip(target_shape, source.shape)]
-    return zoom(source, factors, order=order)
-
+    return float(lo), float(hi)
 
-def _select_best_z_slice(ff_crop):
-    """Return the Z-index with the most field-fraction signal.
 
-    Falls back to the central slice when no signal is present.
-    """
-    ff_per_z = ff_crop.sum(axis=(0, 1))
-    if ff_per_z.max() > 0:
-        return int(ff_per_z.argmax())
-    return ff_crop.shape[2] // 2
-
-
-def _get_bounding_box(mask, pad=5):
-    """Return index slices for the bounding box of a boolean *mask* with padding.
-
-    Returns ``None`` when the mask is empty.
-    """
-    indices = np.where(mask)
-    if not indices[0].size:
-        return None
-    shape = mask.shape
-    return tuple(
-        slice(max(0, int(idx.min()) - pad), min(sz, int(idx.max()) + pad + 1))
-        for idx, sz in zip(indices, shape)
-    )
+def _apply_fixed_percentile_norm(arr, lo, hi):
+    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 main():
@@ -73,20 +49,33 @@ def main():
     max_rois = snakemake.params.max_rois
     n_cols = snakemake.params.n_cols
 
-    spim_img = ZarrNii.from_ome_zarr(snakemake.input.spim,level=snakemake.params.level, downsample_near_isotropic=True,channel_labels=[snakemake.wildcards.stain])
-    mask_img = ZarrNii.from_ome_zarr(snakemake.input.mask,level=0)
-    
-    atlas = ZarrNiiAtlas.from_files(snakemake.input.dseg_nii,snakemake.input.label_tsv)
+    spim_img = ZarrNii.from_ome_zarr(
+        snakemake.input.spim,
+        level=snakemake.params.level,
+        downsample_near_isotropic=True,
+        channel_labels=[snakemake.wildcards.stain],
+    )
+    mask_img = ZarrNii.from_ome_zarr(snakemake.input.mask, level=0)
 
-    dseg_data = atlas.dseg.data.compute()
+    atlas = ZarrNiiAtlas.from_files(snakemake.input.dseg_nii, snakemake.input.label_tsv)
 
+    dseg_data = atlas.dseg.data.compute()
 
     # Voxel dimensions (mm) for physical aspect-ratio correction - not implemented yet
-    # but should be easy with ZarrNii image .scale 
+    # but should be easy with ZarrNii image .scale
     aspect_axial = 1
 
+    spim_img_ds = ZarrNii.from_ome_zarr(
+        snakemake.input.spim,
+        level=(int(snakemake.params.level) + 5),
+        downsample_near_isotropic=True,
+        channel_labels=[snakemake.wildcards.stain],
+    )
 
-    # Global normalisations  - not implemented yet
+    # estimate once globally, from a coarse version of the full image
+    glob_lo, glob_hi = _estimate_global_percentiles(
+        spim_img_ds.data.compute(), pct_low=1, pct_high=99
+    )
 
     # Load atlas label table
     label_df = atlas.labels_df
@@ -145,20 +134,27 @@ def main():
         label_id = int(row["index"])
         label_name = str(row.get("name", label_id))
 
-
-        #get cropped images for this label
+        # get cropped images for this label
         bbox_min, bbox_max = atlas.get_region_bounding_box(region_ids=label_id)
         center_coord = tuple((x + y) / 2 for x, y in zip(bbox_min, bbox_max))
-        spim_crop = spim_img.crop_centered(center_coord, patch_size=(2000,2000,1))
-        mask_crop = mask_img.crop_centered(center_coord, patch_size=(2000,2000,1))
-
+        spim_crop = spim_img.crop_centered(
+            center_coord,
+            patch_size=(snakemake.params.patch_size, snakemake.params.patch_size, 1),
+        )
+        mask_crop = mask_img.crop_centered(
+            center_coord,
+            patch_size=(snakemake.params.patch_size, snakemake.params.patch_size, 1),
+        )
 
-        spim_sl = np.rot90(spim_crop.data[0, :, :].squeeze().compute())
-        mask_sl = np.rot90(mask_crop.data[0, :, :].squeeze().compute())
+        spim_sl = spim_crop.data[0, :, :].squeeze().compute()
+        spim_sl = _apply_fixed_percentile_norm(spim_sl, glob_lo, glob_hi)
+        mask_sl = mask_crop.data[0, :, :].squeeze().compute()
 
         ax.imshow(spim_sl, cmap="gray")
         mask_masked = np.ma.masked_where(mask_sl < 100, mask_sl)
-        ax.imshow(mask_masked, cmap="spring", alpha=0.6, vmin=0, vmax=100, aspect=aspect_axial)
+        ax.imshow(
+            mask_masked, cmap="spring", alpha=0.6, vmin=0, vmax=100, aspect=aspect_axial
+        )
         ax.set_title(label_name, fontsize=7, pad=2)
         ax.set_xticks([])
         ax.set_yticks([])