adding erosion setting

Author: Corentin

.gitignore

@@ -169,4 +169,5 @@ debug_data/
 !nuclei.tif
 !cytoplasm.tif
 !binary_mask_sdh.tif
-data/*
+data/*
+sample_atp_intensity_plot.png

.vscode/launch.json

@@ -26,7 +26,7 @@
             "request": "launch",
             "module": "myoquant",
             "justMyCode": true,
-            "args": ["atp-analysis", "sample_img/sample_atp.jpg", "--cellpose-path", "sample_img/sample_atp_cellpose_mask.tiff"],
+            "args": ["atp-analysis", "sample_img/sample_atp.jpg", "--cellpose-path", "sample_img/sample_atp_cellpose_mask.tiff", "--intensity-method", "mean", "--n-classes", "2", "--erosion"],
         }
     ]
 }

myoquant/commands/run_atp.py

@@ -59,6 +59,23 @@ def atp_analysis(
         None,
         help="Approximative single cell diameter in pixel for CellPose detection. If not specified, Cellpose will try to deduce it.",
     ),
+    channel: int = typer.Option(
+        None,
+        help="Image channel to use for the analysis. If not specified, the analysis will be performed on all three channels.",
+    ),
+    n_classes: int = typer.Option(
+        2,
+        max=10,
+        help="The number of classes of cell to detect. If not specified this is defaulted to two classes.",
+    ),
+    intensity_method: str = typer.Option(
+        "median",
+        help="The method to use to compute the intensity of the cell. Can be either 'median' or 'mean'.",
+    ),
+    erosion: bool = typer.Option(
+        False,
+        help="Perform an erosion on the cells images to remove signal in the cell membrane (usefull for fluo)",
+    ),
     export_map: bool = typer.Option(
         True,
         help="Export the original image with cells painted by classification label.",
@@ -135,6 +152,8 @@ def atp_analysis(
     ) as progress:
         progress.add_task(description="Reading all inputs...", total=None)
         image_ndarray = imread(image_path)
+        if channel is not None:
+            image_ndarray = image_ndarray[:, :, channel]
 
         if mask_path is not None:
             mask_ndarray = imread(mask_path)
@@ -200,8 +219,13 @@ def atp_analysis(
         transient=False,
     ) as progress:
         progress.add_task(description="Detecting fiber types...", total=None)
-        result_df, full_label_map, df_cellpose_details = run_atp_analysis(
-            image_ndarray, mask_cellpose, intensity_threshold
+        result_df, full_label_map, df_cellpose_details, fig = run_atp_analysis(
+            image_ndarray,
+            mask_cellpose,
+            intensity_threshold,
+            n_classes,
+            intensity_method,
+            erosion,
         )
     if export_map:
         with Progress(
@@ -214,7 +238,12 @@ def atp_analysis(
                 description="Blending label and original image together...", total=None
             )
             labelRGB_map = label2rgb(image_ndarray, full_label_map)
-            overlay_img = blend_image_with_label(image_ndarray, labelRGB_map)
+            if channel is not None:
+                overlay_img = blend_image_with_label(
+                    image_ndarray, labelRGB_map, fluo=True
+                )
+            else:
+                overlay_img = blend_image_with_label(image_ndarray, labelRGB_map)
             overlay_filename = image_path.stem + "_label_blend.tiff"
             overlay_img.save(output_path / overlay_filename)
 
@@ -239,6 +268,11 @@ def atp_analysis(
         f"💾 OUTPUT: Summary Table saved as {output_path/csv_name}",
         style="green",
     )
+    plot_name = image_path.stem + "_intensity_plot.png"
+    fig.savefig(output_path / plot_name)
+    console.print(
+        f"💾 OUTPUT: Intensity Plot saved as {output_path/plot_name}", style="green"
+    )
     if export_map:
         console.print(
             f"💾 OUTPUT: Overlay image saved as {output_path/overlay_filename}",

myoquant/src/ATP_analysis.py

@@ -3,45 +3,65 @@
 from sklearn.mixture import GaussianMixture
 from .common_func import extract_single_image, df_from_cellpose_mask
 import numpy as np
+import matplotlib
+
+matplotlib.use("agg")
+
 import matplotlib.pyplot as plt
 
 labels_predict = {1: "fiber type 1", 2: "fiber type 2"}
 np.random.seed(42)
 
 
-def get_all_intensity(image_array, df_cellpose):
+def get_all_intensity(
+    image_array, df_cellpose, intensity_method="median", erosion=False
+):
     all_cell_median_intensity = []
     for index in range(len(df_cellpose)):
-        single_cell_img = extract_single_image(image_array, df_cellpose, index)
+        single_cell_img = extract_single_image(image_array, df_cellpose, index, erosion)
 
         # Calculate median pixel intensity of the cell but ignore 0 values
-        single_cell_median_intensity = np.median(single_cell_img[single_cell_img > 0])
+        if intensity_method == "median":
+            single_cell_median_intensity = np.median(
+                single_cell_img[single_cell_img > 0]
+            )
+        elif intensity_method == "mean":
+            single_cell_median_intensity = np.mean(single_cell_img[single_cell_img > 0])
         all_cell_median_intensity.append(single_cell_median_intensity)
     return all_cell_median_intensity
 
 
-def estimate_threshold(intensity_list):
+def estimate_threshold(intensity_list, n_classes=2):
     density = gaussian_kde(intensity_list)
     density.covariance_factor = lambda: 0.25
     density._compute_covariance()
 
     # Create a vector of 256 values going from 0 to 256:
     xs = np.linspace(0, 255, 256)
     density_xs_values = density(xs)
-    gmm = GaussianMixture(n_components=2).fit(np.array(intensity_list).reshape(-1, 1))
+    gmm = GaussianMixture(n_components=n_classes).fit(
+        np.array(intensity_list).reshape(-1, 1)
+    )
 
     # Find the x values of the two peaks
     peaks_x = np.sort(gmm.means_.flatten())
     # Find the minimum point between the two peaks
-    min_index = np.argmin(density_xs_values[(xs > peaks_x[0]) & (xs < peaks_x[1])])
-    threshold = peaks_x[0] + xs[min_index]
 
-    return threshold
+    threshold_list = []
+    length = len(peaks_x)
+    for index, peaks in enumerate(peaks_x):
+        if index == length - 1:
+            break
+        min_index = np.argmin(
+            density_xs_values[(xs > peaks) & (xs < peaks_x[index + 1])]
+        )
+        threshold_list.append(peaks + xs[min_index])
+    return threshold_list
 
 
-def plot_density(all_cell_median_intensity, intensity_threshold):
+def plot_density(all_cell_median_intensity, intensity_threshold, n_classes=2):
     if intensity_threshold == 0:
-        intensity_threshold = estimate_threshold(all_cell_median_intensity)
+        intensity_threshold = estimate_threshold(all_cell_median_intensity, n_classes)
     fig, ax = plt.subplots(figsize=(10, 5))
     density = gaussian_kde(all_cell_median_intensity)
     density.covariance_factor = lambda: 0.25
@@ -51,22 +71,44 @@ def plot_density(all_cell_median_intensity, intensity_threshold):
     xs = np.linspace(0, 255, 256)
     density_xs_values = density(xs)
     ax.plot(xs, density_xs_values, label="Estimated Density")
-    ax.axvline(x=intensity_threshold, color="red", label="Threshold")
+    for values in intensity_threshold:
+        ax.axvline(x=values, color="red", label="Threshold")
     ax.set_xlabel("Pixel Intensity")
     ax.set_ylabel("Density")
     ax.legend()
     return fig
 
 
-def predict_all_cells(histo_img, cellpose_df, intensity_threshold):
-    all_cell_median_intensity = get_all_intensity(histo_img, cellpose_df)
+def merge_peaks_too_close(peak_list):
+    pass
+
+
+def classify_cells_intensity(all_cell_median_intensity, intensity_threshold):
+    muscle_fiber_type_all = []
+    for intensity in all_cell_median_intensity:
+        class_cell = np.searchsorted(intensity_threshold, intensity, side="right")
+        muscle_fiber_type_all.append(class_cell)
+    return muscle_fiber_type_all
+
+
+def predict_all_cells(
+    histo_img,
+    cellpose_df,
+    intensity_threshold,
+    n_classes=2,
+    intensity_method="median",
+    erosion=False,
+):
+    all_cell_median_intensity = get_all_intensity(
+        histo_img, cellpose_df, intensity_method, erosion
+    )
     if intensity_threshold is None:
-        intensity_threshold = estimate_threshold(all_cell_median_intensity)
+        intensity_threshold = estimate_threshold(all_cell_median_intensity, n_classes)
 
-    muscle_fiber_type_all = [
-        1 if x > intensity_threshold else 2 for x in all_cell_median_intensity
-    ]
-    return muscle_fiber_type_all, all_cell_median_intensity
+    muscle_fiber_type_all = classify_cells_intensity(
+        all_cell_median_intensity, intensity_threshold
+    )
+    return muscle_fiber_type_all, all_cell_median_intensity, intensity_threshold
 
 
 def paint_full_image(image_atp, df_cellpose, class_predicted_all):
@@ -76,24 +118,46 @@ def paint_full_image(image_atp, df_cellpose, class_predicted_all):
     # for index in track(range(len(df_cellpose)), description="Painting cells"):
     for index in range(len(df_cellpose)):
         single_cell_mask = df_cellpose.iloc[index, 9].copy()
-        if class_predicted_all[index] == 1:
-            image_atp_paint[
-                df_cellpose.iloc[index, 5] : df_cellpose.iloc[index, 7],
-                df_cellpose.iloc[index, 6] : df_cellpose.iloc[index, 8],
-            ][single_cell_mask] = 1
-        elif class_predicted_all[index] == 2:
-            image_atp_paint[
-                df_cellpose.iloc[index, 5] : df_cellpose.iloc[index, 7],
-                df_cellpose.iloc[index, 6] : df_cellpose.iloc[index, 8],
-            ][single_cell_mask] = 2
+        image_atp_paint[
+            df_cellpose.iloc[index, 5] : df_cellpose.iloc[index, 7],
+            df_cellpose.iloc[index, 6] : df_cellpose.iloc[index, 8],
+        ][single_cell_mask] = (
+            class_predicted_all[index] + 1
+        )
     return image_atp_paint
 
 
-def run_atp_analysis(image_array, mask_cellpose, intensity_threshold=None):
+def label_list_from_threhsold(threshold_list):
+    label_list = []
+    length = len(threshold_list)
+    for index, threshold in enumerate(threshold_list):
+        if index == 0:
+            label_list.append(f"<{threshold}")
+        if index == length - 1:
+            label_list.append(f">{threshold}")
+        else:
+            label_list.append(f">{threshold} & <{threshold_list[index+1]}")
+    return label_list
+
+
+def run_atp_analysis(
+    image_array,
+    mask_cellpose,
+    intensity_threshold=None,
+    n_classes=2,
+    intensity_method="median",
+    erosion=False,
+):
     df_cellpose = df_from_cellpose_mask(mask_cellpose)
-    class_predicted_all, intensity_all = predict_all_cells(
-        image_array, df_cellpose, intensity_threshold
+    class_predicted_all, intensity_all, intensity_threshold = predict_all_cells(
+        image_array,
+        df_cellpose,
+        intensity_threshold,
+        n_classes,
+        intensity_method,
+        erosion,
     )
+    fig = plot_density(intensity_all, intensity_threshold, n_classes)
     df_cellpose["muscle_cell_type"] = class_predicted_all
     df_cellpose["cell_intensity"] = intensity_all
     count_per_label = np.unique(class_predicted_all, return_counts=True)
@@ -102,15 +166,16 @@ def run_atp_analysis(image_array, mask_cellpose, intensity_threshold=None):
     headers = ["Feature", "Raw Count", "Proportion (%)"]
     data = []
     data.append(["Muscle Fibers", len(class_predicted_all), 100])
+    label_list = label_list_from_threhsold(intensity_threshold)
     for index, elem in enumerate(count_per_label[0]):
         data.append(
             [
-                labels_predict[int(elem)],
+                label_list[int(elem)],
                 count_per_label[1][int(index)],
                 100 * count_per_label[1][int(index)] / len(class_predicted_all),
             ]
         )
     result_df = pd.DataFrame(columns=headers, data=data)
     # Paint The Full Image
     full_label_map = paint_full_image(image_array, df_cellpose, class_predicted_all)
-    return result_df, full_label_map, df_cellpose
+    return result_df, full_label_map, df_cellpose, fig

myoquant/src/common_func.py

@@ -5,6 +5,8 @@
 
 os.environ["TF_CPP_MIN_LOG_LEVEL"] = "3"
 import sys
+import math
+
 import tensorflow as tf
 import torch
 from cellpose.models import Cellpose
@@ -14,11 +16,14 @@
 import numpy as np
 from PIL import Image
 from skimage.measure import regionprops_table
+from skimage.morphology import binary_erosion
 import pandas as pd
 
 # from .gradcam import make_gradcam_heatmap, save_and_display_gradcam
 from .random_brightness import RandomBrightness
 
+import imageio
+
 tf.random.set_seed(42)
 np.random.seed(42)
 
@@ -168,12 +173,21 @@ def df_from_stardist_mask(mask, intensity_image=None):
     return df_stardist
 
 
-def extract_single_image(raw_image, df_props, index):
+def extract_single_image(raw_image, df_props, index, erosion=False):
     single_entity_img = raw_image[
         df_props.iloc[index, 5] : df_props.iloc[index, 7],
         df_props.iloc[index, 6] : df_props.iloc[index, 8],
     ].copy()
+    surface_area = df_props.iloc[index, 1]
+    cell_radius = math.sqrt(surface_area / math.pi)
     single_entity_mask = df_props.iloc[index, 9]
+    erosion_size = int(cell_radius / 5)  # 20% of the cell
+    if erosion:
+        for i in range(erosion_size):
+            single_entity_mask = binary_erosion(
+                single_entity_mask, out=single_entity_mask
+            )
+
     single_entity_img[~single_entity_mask] = 0
     return single_entity_img
 
@@ -192,6 +206,14 @@ def label2rgb(img_ndarray, label_map):
         0: [255, 255, 255],
         1: [15, 157, 88],
         2: [219, 68, 55],
+        3: [100, 128, 170],
+        4: [231, 204, 143],
+        5: [202, 137, 115],
+        6: [178, 143, 172],
+        7: [144, 191, 207],
+        8: [148, 187, 187],
+        9: [78, 86, 105],
+        10: [245, 90, 87],
     }
     img_rgb = np.zeros((img_ndarray.shape[0], img_ndarray.shape[1], 3), dtype=np.uint8)