From 8362c2e3e92e0364649318366d93a7118f9629aa Mon Sep 17 00:00:00 2001
From: Leonardo Schwarz <leonardo.schwarz@fgcz.ethz.ch>
Date: Fri, 5 Jul 2024 11:51:23 +0200
Subject: [PATCH] plot_intensity_threshold.py

---
 .../pipeline_config/artifacts_mapping.py      |  2 +
 .../workflow/qc/plot_intensity_threshold.py   | 78 +++++++++++++++++++
 .../workflow/rules/rules_qc.smk               | 13 ++++
 3 files changed, 93 insertions(+)
 create mode 100644 src/depiction_targeted_preproc/workflow/qc/plot_intensity_threshold.py

diff --git a/src/depiction_targeted_preproc/pipeline_config/artifacts_mapping.py b/src/depiction_targeted_preproc/pipeline_config/artifacts_mapping.py
index 0fe2521..5421ef2 100644
--- a/src/depiction_targeted_preproc/pipeline_config/artifacts_mapping.py
+++ b/src/depiction_targeted_preproc/pipeline_config/artifacts_mapping.py
@@ -29,6 +29,8 @@
         #        "cluster_default_stats_kmeans.csv",
         "cluster_default_kmeans.png",
         "cluster_default_hdbscan.png",
+        "qc/plot_intensity_threshold_all.pdf",
+        "qc/plot_intensity_threshold_fg.pdf",
         # "exp_plot_map_comparison.pdf",
         # "qc/plot_marker_presence_mini.pdf",
     ],
diff --git a/src/depiction_targeted_preproc/workflow/qc/plot_intensity_threshold.py b/src/depiction_targeted_preproc/workflow/qc/plot_intensity_threshold.py
new file mode 100644
index 0000000..e8d93e6
--- /dev/null
+++ b/src/depiction_targeted_preproc/workflow/qc/plot_intensity_threshold.py
@@ -0,0 +1,78 @@
+from __future__ import annotations
+
+from pathlib import Path
+
+import cyclopts
+import matplotlib.pyplot as plt
+import numpy as np
+import polars as pl
+from numpy.typing import NDArray
+from tqdm import tqdm
+
+from depiction.image.multi_channel_image import MultiChannelImage
+
+app = cyclopts.App()
+
+
+@app.default
+def main(image_hdf5: Path, output_all_pixels_pdf: Path, output_foreground_pixels_pdf: Path) -> None:
+    image = MultiChannelImage.read_hdf5(image_hdf5)
+    data_flat = image.data_flat.values.ravel()
+
+    # compute thresholds to evaluate
+    v_min = 0
+    v_max = np.percentile(data_flat, 99)
+    thresholds = np.logspace(np.log10(v_min + 1), np.log10(v_max + 1), 500) - 1
+    plot_threshold_all_pixels(image=image, thresholds=thresholds, output_pdf=output_all_pixels_pdf)
+    plot_threshold_foreground_only(image=image, thresholds=thresholds, output_pdf=output_foreground_pixels_pdf)
+
+
+def plot_threshold_all_pixels(image: MultiChannelImage, thresholds: NDArray[float], output_pdf: Path) -> None:
+    collect = []
+    data_flat = image.data_flat
+    for threshold in tqdm(thresholds):
+        counts = (data_flat > threshold).sum("c")
+        collect.append(
+            {"threshold": threshold, "mean": counts.mean(), "p25": counts.quantile(0.25), "p75": counts.quantile(0.75)}
+        )
+    df = pl.DataFrame(collect)
+    plt.figure()
+    plt.title(f"Detected Targets by Intensity Threshold (N={np.prod(data_flat.shape):,})")
+    plt.fill_between(df["threshold"], df["p25"], df["p75"], color="gray", alpha=0.5, label="p25-p75")
+    plt.plot(df["threshold"], df["mean"], label="mean")
+    plt.xlabel("Threshold")
+    plt.ylabel("Detected Targets (agg. over pixels)")
+    plt.grid()
+    plt.legend()
+    plt.savefig(output_pdf, bbox_inches="tight")
+
+
+def plot_threshold_foreground_only(image: MultiChannelImage, thresholds: NDArray[float], output_pdf: Path) -> None:
+    background_targets = 5
+    threshold = 25
+    data = image.data_spatial
+    bg = (data > threshold).sum("c") <= background_targets
+    fg = ~bg
+    fg.plot.imshow(yincrease=False, cmap="gray")
+    # TODO if this actually gets used it will have to be supplemented with information about the foreground mask
+    collect = []
+    data_fg_flat = data.where(fg).stack(i=("y", "x")).dropna("i", how="all")
+    for threshold in tqdm(thresholds):
+        counts = (data_fg_flat > threshold).sum("c")
+        collect.append(
+            {"threshold": threshold, "mean": counts.mean(), "p25": counts.quantile(0.25), "p75": counts.quantile(0.75)}
+        )
+    df = pl.DataFrame(collect)
+    plt.figure()
+    plt.title(f"Detected Targets by Intensity Threshold (N={np.prod(data_fg_flat.shape):,})")
+    plt.fill_between(df["threshold"], df["p25"], df["p75"], color="gray", alpha=0.5, label="p25-p75")
+    plt.plot(df["threshold"], df["mean"], label="mean")
+    plt.xlabel("Threshold")
+    plt.ylabel("Detected Targets (agg. over pixels)")
+    plt.grid()
+    plt.legend()
+    plt.savefig(output_pdf, bbox_inches="tight")
+
+
+if __name__ == "__main__":
+    app()
diff --git a/src/depiction_targeted_preproc/workflow/rules/rules_qc.smk b/src/depiction_targeted_preproc/workflow/rules/rules_qc.smk
index 69e8519..39b93e5 100644
--- a/src/depiction_targeted_preproc/workflow/rules/rules_qc.smk
+++ b/src/depiction_targeted_preproc/workflow/rules/rules_qc.smk
@@ -162,3 +162,16 @@ rule qc_plot_scan_direction:
         "python -m depiction_targeted_preproc.workflow.qc.plot_scan_direction"
         " --input-imzml-path {input.imzml[0]}"
         " --output-pdf {output.pdf}"
+
+
+rule qc_plot_intensity_threshold:
+    input:
+        image="{sample}/images_default.hdf5",
+    output:
+        pdf_all="{sample}/qc/plot_intensity_threshold_all.pdf",
+        pdf_fg="{sample}/qc/plot_intensity_threshold_fg.pdf",
+    shell:
+        "python -m depiction_targeted_preproc.workflow.qc.plot_intensity_threshold"
+        " --image-hdf5 {input.image}"
+        " --output-all-pixels-pdf {output.pdf_all}"
+        " --output-foreground-pixels-pdf {output.pdf_fg}"