move polar map scripts to a new model and add more

mu_map/polar_map/eval_perfusion.py

+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+
+from mu_map.polar_map.prepare import headers
+
+data = pd.read_csv("data/polar_maps/perfusion.csv")
+
+baseline = data[data[headers.ac] & (data[headers.type] == "symbia")]
+
+correction_none = data[~data[headers.ac]]
+correction_syn = data[data[headers.ac] & (data[headers.type] == "synthetic")]
+correction_ct = data[data[headers.ac] & (data[headers.type] == "ct")]
+
+keys_segments = [f"segment_{i}" for i in range(1, 18)]
+
+_baseline = baseline[keys_segments].values
+_correction_none = correction_none[keys_segments].values
+_correction_syn = correction_syn[keys_segments].values
+_correction_ct = correction_ct[keys_segments].values
+
+
+def absolute_percent_error(prediction: np.ndarray, target: np.ndarray) -> float:
+    mean_p = prediction.mean(axis=0)
+    mean_t = target.mean(axis=0)
+    diff = np.absolute(mean_p - mean_t)
+    return 100.0 * diff / mean_t
+
+
+def apc(prediction: np.ndarray, target: np.ndarray) -> float:
+    return absolute_percent_error(prediction, target)
+
+
+combinations = [
+    ("NoAC to AC", _correction_none, _baseline),
+    ("SYN to AC", _correction_syn, _baseline),
+    ("CT to AC", _correction_ct, _baseline),
+    ("NoAC to CT", _correction_none, _correction_ct),
+    ("SYN to CT", _correction_syn, _correction_ct),
+]
+
+for label, prediction, target in combinations:
+    apcs = absolute_percent_error(prediction, target)
+    _mean = apcs.mean()
+    _std = apcs.std()
+    print(f"APE - {label:>12}: {_mean:.3f}±{_std:.3f}")
+
+
+n_plots = 5
+plot_width = 4
+fig, axs = plt.subplots(1, n_plots, figsize=(n_plots * plot_width + 1, plot_width + 1))
+fig.suptitle("Correlation of Segment Perfusion")
+
+
+def jitter(array: np.ndarray, scale=1):
+    return array + (np.random.random(array.shape) * scale - 0.5 * scale)
+
+
+axs[0].scatter(
+    jitter(_baseline.flatten()),
+    jitter(_correction_none.flatten()),
+    c="#67a9cf",
+    s=20,
+    alpha=0.6,
+    edgecolor="black",
+)
+axs[1].scatter(
+    jitter(_baseline.flatten()),
+    jitter(_correction_syn.flatten()),
+    c="#67a9cf",
+    s=20,
+    alpha=0.6,
+    edgecolor="black",
+)
+axs[2].scatter(
+    jitter(_baseline.flatten()),
+    jitter(_correction_ct.flatten()),
+    c="#67a9cf",
+    s=20,
+    alpha=0.6,
+    edgecolor="black",
+)
+
+axs[3].scatter(
+    jitter(_correction_ct.flatten()),
+    jitter(_correction_none.flatten()),
+    c="#67a9cf",
+    s=20,
+    alpha=0.6,
+    edgecolor="black",
+)
+axs[4].scatter(
+    jitter(_correction_ct.flatten()),
+    jitter(_correction_syn.flatten()),
+    c="#67a9cf",
+    s=20,
+    alpha=0.6,
+    edgecolor="black",
+)
+
+axs[0].set_xlabel("Perfusion AC")
+axs[1].set_xlabel("Perfusion AC")
+axs[2].set_xlabel("Perfusion AC")
+axs[3].set_xlabel("Perfusion AC CT")
+axs[4].set_xlabel("Perfusion AC CT")
+
+axs[0].set_ylabel("Perfusion NoAC")
+axs[1].set_ylabel("Perfusion AC SYN")
+axs[2].set_ylabel("Perfusion AC CT")
+axs[3].set_ylabel("Perfusion NoAC")
+axs[4].set_ylabel("Perfusion AC SYN")
+
+axs[0].set_title("  NoAC vs. AC")
+axs[1].set_title("AC SYN vs. AC")
+axs[2].set_title(" AC CT vs. AC")
+axs[3].set_title("  NoAC vs. AC CT")
+axs[4].set_title("AC SYN vs. AC CT")
+
+for ax in axs:
+    ax.set_xlim((30, 100))
+    ax.set_ylim((30, 100))
+    ax.plot((30, 100), (30, 100), color="#ef8a62")
+
+plt.tight_layout()
+plt.show()

mu_map/data/get_perfusion_scores.py → mu_map/polar_map/get_perfusion.py

@@ -46,7 +46,7 @@ if __name__ == "__main__":
 
     import pandas as pd
 
-    from mu_map.data.prepare_polar_maps import headers
+    from mu_map.polar_map.prepare import headers
 
     parser = argparse.ArgumentParser(
         formatter_class=argparse.ArgumentDefaultsHelpFormatter,
@@ -85,6 +85,8 @@ if __name__ == "__main__":
     args.perfusion_csv = os.path.join(args.polar_map_dir, args.perfusion_csv)
 
     data = pd.read_csv(args.csv)
+    data = data[data[headers.segments]]
+
     rows = []
     for i in range(len(data)):
         row = data.iloc[i]

mu_map/data/prepare_polar_maps.py → mu_map/polar_map/prepare.py

@@ -31,9 +31,11 @@ def is_scatter_corrected(dcm: dcm_type) -> bool:
 def is_attenuation_corrected(dcm: dcm_type) -> bool:
     return not ("NoAC" in dcm.SeriesDescription)
 
+
 def shows_segments(dcm: dcm_type) -> bool:
     return "PolarMapSeg" in dcm.SeriesDescription
 
+
 def get_type(dcm: dcm_type) -> str:
     description = dcm.SeriesDescription.lower()
     if "syn" in description:
@@ -183,6 +185,17 @@ if __name__ == "__main__":
         polar_map = img[top:bottom, left:right]
         polar_map = cv.cvtColor(polar_map, cv.COLOR_RGB2BGR)
 
+        # Note: the images created by Syngo.Via show polar maps of rest studies
+        # further down per default. This little bit of code selects this area
+        # further down if correcting this was forgotten
+        if polar_map.mean() < 10:
+            offset = 649
+            top = top + offset
+            bottom = bottom + offset
+
+            polar_map = img[top:bottom, left:right]
+            polar_map = cv.cvtColor(polar_map, cv.COLOR_RGB2BGR)
+
         _ac = "ac" if row[headers.ac] else "nac"
         _sc = "sc" if row[headers.sc] else "nsc"
         _seg = "segments" if row[headers.segments] else "gray"

mu_map/polar_map/visualize.py

+import argparse
+import os
+
+import cv2 as cv
+import matplotlib as mlp
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+
+from mu_map.polar_map.prepare import headers
+
+
+def get_circular_mark(shape: np.ndarray, channels:int=1):
+    mask = np.full((*shape, channels), 0, np.uint8)
+    cx, cy = np.array(mask.shape[:2]) // 2
+    mask = cv.circle(
+        mask,
+        center=(cx-1, cy),
+        radius=cx - 2,
+        color=(255,) * channels,
+        thickness=cv.FILLED,
+    )
+    mask = mask == 255
+    return mask[:, :, 0] if channels == 1 else mask
+
+
+parser = argparse.ArgumentParser(description="Visualize polar maps of different reconstructions", formatter_class=argparse.ArgumentDefaultsHelpFormatter)
+parser.add_argument("--polar_map_dir", type=str, default="data/polar_maps", help="directory containing the polar map images")
+parser.add_argument("--images_dir", type=str, default="images", help="sub-directory under <polar_map_dir> containing the actual image files")
+parser.add_argument("--csv", type=str, default="polar_maps.csv", help="file under <polar_map_dir> containing meta information about the polar maps")
+parser.add_argument("--baseline", choices=["symbia", "stir"], default="symbia", help="select the polar map treated as the baseline")
+parser.add_argument("--id", type=int, help="select a specific study to show by its id")
+parser.add_argument("--relative_difference", action="store_true", help="show the difference as a percentage difference to the baseline")
+args = parser.parse_args()
+
+args.images_dir = os.path.join(args.polar_map_dir, args.images_dir)
+args.csv = os.path.join(args.polar_map_dir, args.csv)
+
+meta = pd.read_csv(args.csv)
+ids = meta[headers.id].unique()
+
+if args.id:
+    assert args.id in ids, f"Id {args.id} is not available. Chose one of {ids}."
+    ids = [args.id]
+
+
+for _id in ids:
+    print(f"Show id {_id:03d}")
+    _meta = meta[(meta[headers.id] == _id) & ~(meta[headers.segments])]
+
+    file_recon_ac = _meta[(_meta[headers.type] == "symbia") & _meta[headers.ac]][headers.file].values[0]
+    file_recon_nac = _meta[~_meta[headers.ac]][headers.file].values[0]
+    file_recon_syn = _meta[_meta[headers.type] == "synthetic"][headers.file].values[0]
+    file_recon_ct = _meta[_meta[headers.type] == "ct"][headers.file].values[0]
+
+    recon_ac = cv.imread(os.path.join(args.images_dir, file_recon_ac), cv.IMREAD_GRAYSCALE)
+    recon_nac = cv.imread(os.path.join(args.images_dir, file_recon_nac), cv.IMREAD_GRAYSCALE)
+    recon_syn = cv.imread(os.path.join(args.images_dir, file_recon_syn), cv.IMREAD_GRAYSCALE)
+    recon_ct = cv.imread(os.path.join(args.images_dir, file_recon_ct), cv.IMREAD_GRAYSCALE)
+
+    baseline = recon_ac.copy() if args.baseline == "symbia" else recon_ct.copy()
+    recons = [recon_nac, recon_syn, recon_ct]
+    labels = ["No AC", "AC SYN", "AC CT"]
+
+
+    diffs = map(lambda recon: (recon.astype(float) - baseline) * 100.0 / 255.0, recons)
+
+    if args.relative_difference:
+        divider = np.where(baseline > 0, baseline, 1)
+        diffs = map(lambda diff: 100 * diff / divider, diffs)
+
+    diffs = list(diffs)
+
+    diff_min = min(map(lambda diff: diff.min(), diffs))
+    diff_max = max(map(lambda diff: diff.max(), diffs))
+    diff_min = -max(abs(diff_min), diff_max)
+    diff_max = max(abs(diff_min), diff_max)
+
+    diffs = map(lambda diff: (diff - diff_min) / (diff_max - diff_min), diffs)
+    diffs = list(diffs)
+
+    fig, axs = plt.subplots(2, 4, figsize=(16, 8))
+    for ax in axs.flatten():
+        ax.set_axis_off()
+
+    mask = np.full((*baseline.shape, 4), 0, np.uint8)
+    cx, cy = np.array(mask.shape[:2]) // 2
+    mask = cv.circle(
+        mask,
+        center=(cx - 1, cy + 1),
+        radius=cx - 2,
+        color=(255, 255, 255, 255),
+        thickness=cv.FILLED,
+    )
+    mask = mask == 255
+    black = np.zeros(mask.shape, np.uint8)
+
+
+    def colormap_and_mask(img, cm):
+        _img = cm(img)
+        return np.where(mask, _img, black)
+
+    cm_plasma = mlp.colormaps["plasma"]
+    cm_divergent = mlp.colormaps["RdYlBu"].reversed()
+    axs[0, 0].imshow(colormap_and_mask(baseline, cm_plasma))
+    axs[0, 0].set_title("AC" if args.baseline == "symbia" else "AC - CT")
+    for i, (recon, diff, label) in enumerate(zip(recons, diffs, labels), start=1):
+        axs[0, i].imshow(colormap_and_mask(recon, cm_plasma))
+        axs[0, i].set_title(label)
+        axs[1, i].imshow(colormap_and_mask(diff, cm_divergent))
+
+    fig.colorbar(
+        mlp.cm.ScalarMappable(
+            norm=mlp.colors.Normalize(vmin=diff_min, vmax=diff_max), cmap=cm_divergent
+        ),
+        ax=axs[1, 1:4],
+    )
+    fig.colorbar(
+        mlp.cm.ScalarMappable(
+            norm=mlp.colors.Normalize(vmin=0, vmax=100), cmap=cm_plasma
+        ),
+        ax=axs[0, 1:4]
+    )
+    plt.show()