Merge pull request #4 from lambda-science/dev

Dev

myoquant/commands/run_sdh.py

@@ -164,21 +164,23 @@ def sdh_analysis(
         progress.add_task(description="Reading all inputs...", total=None)
         image_ndarray_sdh = imread(image_path)
 
-    if mask_path is not None:
-        mask_ndarray = imread(mask_path)
-        if np.unique(mask_ndarray).shape[0] != 2:
-            console.print(
-                "The mask image should be a binary image with only 2 values (0 and 1).",
-                style="red",
-            )
-            raise ValueError
-        if len(image_ndarray_sdh.shape) > 2:
-            mask_ndarray = np.repeat(
-                mask_ndarray.reshape(mask_ndarray.shape[0], mask_ndarray.shape[1], 1),
-                image_ndarray_sdh.shape[2],
-                axis=2,
-            )
-        image_ndarray_sdh = image_ndarray_sdh * mask_ndarray
+        if mask_path is not None:
+            mask_ndarray = imread(mask_path)
+            if np.unique(mask_ndarray).shape[0] != 2:
+                console.print(
+                    "The mask image should be a binary image with only 2 values (0 and 1).",
+                    style="red",
+                )
+                raise ValueError
+            if len(image_ndarray_sdh.shape) > 2:
+                mask_ndarray = np.repeat(
+                    mask_ndarray.reshape(
+                        mask_ndarray.shape[0], mask_ndarray.shape[1], 1
+                    ),
+                    image_ndarray_sdh.shape[2],
+                    axis=2,
+                )
+            image_ndarray_sdh = image_ndarray_sdh * mask_ndarray
         if cellpose_path is not None:
             mask_cellpose = imread(cellpose_path)
 

myoquant/src/ATP_analysis.py

@@ -1,9 +1,9 @@
 import pandas as pd
-from skimage.measure import regionprops_table
 from scipy.stats import gaussian_kde
 from sklearn.mixture import GaussianMixture
-
+from .common_func import extract_single_image, df_from_cellpose_mask
 import numpy as np
+import matplotlib.pyplot as plt
 
 labels_predict = {1: "fiber type 1", 2: "fiber type 2"}
 np.random.seed(42)
@@ -12,13 +12,8 @@
 def get_all_intensity(image_array, df_cellpose):
     all_cell_median_intensity = []
     for index in range(len(df_cellpose)):
-        single_cell_img = image_array[
-            df_cellpose.iloc[index, 5] : df_cellpose.iloc[index, 7],
-            df_cellpose.iloc[index, 6] : df_cellpose.iloc[index, 8],
-        ].copy()
+        single_cell_img = extract_single_image(image_array, df_cellpose, index)
 
-        single_cell_mask = df_cellpose.iloc[index, 9].copy()
-        single_cell_img[~single_cell_mask] = 0
         # Calculate median pixel intensity of the cell but ignore 0 values
         single_cell_median_intensity = np.median(single_cell_img[single_cell_img > 0])
         all_cell_median_intensity.append(single_cell_median_intensity)
@@ -44,6 +39,25 @@ def estimate_threshold(intensity_list):
     return threshold
 
 
+def plot_density(all_cell_median_intensity, intensity_threshold):
+    if intensity_threshold == 0:
+        intensity_threshold = estimate_threshold(all_cell_median_intensity)
+    fig, ax = plt.subplots(figsize=(10, 5))
+    density = gaussian_kde(all_cell_median_intensity)
+    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)
+    ax.plot(xs, density_xs_values, label="Estimated Density")
+    ax.axvline(x=intensity_threshold, 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)
     if intensity_threshold is None:
@@ -76,19 +90,7 @@ def paint_full_image(image_atp, df_cellpose, class_predicted_all):
 
 
 def run_atp_analysis(image_array, mask_cellpose, intensity_threshold=None):
-    props_cellpose = regionprops_table(
-        mask_cellpose,
-        properties=[
-            "label",
-            "area",
-            "centroid",
-            "eccentricity",
-            "bbox",
-            "image",
-            "perimeter",
-        ],
-    )
-    df_cellpose = pd.DataFrame(props_cellpose)
+    df_cellpose = df_from_cellpose_mask(mask_cellpose)
     class_predicted_all, intensity_all = predict_all_cells(
         image_array, df_cellpose, intensity_threshold
     )

myoquant/src/HE_analysis.py

@@ -1,39 +1,27 @@
 import numpy as np
 import pandas as pd
 from skimage.draw import line
-from skimage.measure import regionprops_table
-
-from .draw_line import *
+from .common_func import (
+    extract_single_image,
+    df_from_cellpose_mask,
+    df_from_stardist_mask,
+)
+from .draw_line import (
+    line_equation,
+    calculate_intersection,
+    calculate_closest_point,
+    calculate_distance,
+)
+import matplotlib.pyplot as plt
 
 
 def extract_ROIs(histo_img, index, cellpose_df, mask_stardist):
-    single_cell_img = histo_img[
-        cellpose_df.iloc[index, 5] : cellpose_df.iloc[index, 7],
-        cellpose_df.iloc[index, 6] : cellpose_df.iloc[index, 8],
-    ].copy()
-    nucleus_single_cell_img = mask_stardist[
-        cellpose_df.iloc[index, 5] : cellpose_df.iloc[index, 7],
-        cellpose_df.iloc[index, 6] : cellpose_df.iloc[index, 8],
-    ].copy()
+    single_cell_img = extract_single_image(histo_img, cellpose_df, index)
+    nucleus_single_cell_img = extract_single_image(mask_stardist, cellpose_df, index)
     single_cell_mask = cellpose_df.iloc[index, 9]
-    single_cell_img[~single_cell_mask] = 0
-    nucleus_single_cell_img[~single_cell_mask] = 0
-
-    props_nuc_single = regionprops_table(
-        nucleus_single_cell_img,
-        intensity_image=single_cell_img,
-        properties=[
-            "label",
-            "area",
-            "centroid",
-            "eccentricity",
-            "bbox",
-            "image",
-            "perimeter",
-            "feret_diameter_max",
-        ],
+    df_nuc_single = df_from_stardist_mask(
+        nucleus_single_cell_img, intensity_image=single_cell_img
     )
-    df_nuc_single = pd.DataFrame(props_nuc_single)
     return single_cell_img, nucleus_single_cell_img, single_cell_mask, df_nuc_single
 
 
@@ -45,14 +33,22 @@ def single_cell_analysis(
     y_fiber,
     cell_label,
     internalised_threshold=0.75,
+    draw_and_return=False,
 ):
+    if draw_and_return:
+        fig_size = (5, 5)
+        f, ax = plt.subplots(figsize=fig_size)
+        ax.scatter(x_fiber, y_fiber, color="white")
+
     n_nuc, n_nuc_intern, n_nuc_periph = 0, 0, 0
     new_row_lst = []
     new_col_names = df_nuc_single.columns.tolist()
     new_col_names.append("internalised")
     new_col_names.append("score_eccentricity")
     new_col_names.append("cell label")
     for _, value in df_nuc_single.iterrows():
+        if draw_and_return:
+            ax.scatter(value[3], value[2], color="red")
         n_nuc += 1
         value = value.tolist()
         # Extend line and find closest point
@@ -68,6 +64,21 @@ def single_cell_analysis(
             m, b, (single_cell_img.shape[0], single_cell_img.shape[1])
         )
         border_point = calculate_closest_point(value[3], value[2], intersections_lst)
+        if draw_and_return:
+            ax.plot(
+                (x_fiber, border_point[0]),
+                (y_fiber, border_point[1]),
+                "ro--",
+                linewidth=1,
+                markersize=1,
+            )
+            ax.plot(
+                (x_fiber, value[3]),
+                (y_fiber, value[2]),
+                "go--",
+                linewidth=1,
+                markersize=1,
+            )
         rr, cc = line(
             int(y_fiber),
             int(x_fiber),
@@ -95,6 +106,8 @@ def single_cell_analysis(
                         value.append(ratio_dist)
                         value.append(cell_label)
                     new_row_lst.append(value)
+                    if draw_and_return:
+                        ax.scatter(coords[1], coords[0], color="red", s=10)
                     break
             except IndexError:
                 coords = list(zip(rr, cc))[index3 - 1]
@@ -114,8 +127,13 @@ def single_cell_analysis(
                     value.append(ratio_dist)
                     value.append(cell_label)
                 new_row_lst.append(value)
+                if draw_and_return:
+                    ax.scatter(coords[1], coords[0], color="red", s=10)
+                    ax.axis("off")
                 break
     df_nuc_single_stats = pd.DataFrame(new_row_lst, columns=new_col_names)
+    if draw_and_return:
+        return n_nuc, n_nuc_intern, n_nuc_periph, df_nuc_single_stats, ax
     return n_nuc, n_nuc_intern, n_nuc_periph, df_nuc_single_stats
 
 
@@ -172,20 +190,7 @@ def paint_histo_img(histo_img, cellpose_df, prediction_df):
 
 
 def run_he_analysis(image_ndarray, mask_cellpose, mask_stardist, eccentricity_thresh):
-    props_cellpose = regionprops_table(
-        mask_cellpose,
-        properties=[
-            "label",
-            "area",
-            "centroid",
-            "eccentricity",
-            "bbox",
-            "image",
-            "perimeter",
-            "feret_diameter_max",
-        ],
-    )
-    df_cellpose = pd.DataFrame(props_cellpose)
+    df_cellpose = df_from_cellpose_mask(mask_cellpose)
     df_nuc_analysis, all_nuc_df_stats = predict_all_cells(
         image_ndarray, df_cellpose, mask_stardist, eccentricity_thresh
     )

myoquant/src/SDH_analysis.py

@@ -4,9 +4,7 @@
 
 import pandas as pd
 import tensorflow as tf
-from rich.progress import track
-from skimage.measure import regionprops_table
-
+from .common_func import extract_single_image, df_from_cellpose_mask
 from .gradcam import make_gradcam_heatmap, save_and_display_gradcam
 import numpy as np
 
@@ -32,14 +30,7 @@ def resize_batch_cells(histo_img, cellpose_df):
     img_array_full = np.empty((len(cellpose_df), 256, 256, 3))
     # for index in track(range(len(cellpose_df)), description="Resizing cells"):
     for index in range(len(cellpose_df)):
-        single_cell_img = histo_img[
-            cellpose_df.iloc[index, 5] : cellpose_df.iloc[index, 7],
-            cellpose_df.iloc[index, 6] : cellpose_df.iloc[index, 8],
-        ].copy()
-
-        single_cell_mask = cellpose_df.iloc[index, 9].copy()
-        single_cell_img[~single_cell_mask] = 0
-
+        single_cell_img = extract_single_image(histo_img, cellpose_df, index)
         img_array_full[index] = tf.image.resize(single_cell_img, (256, 256))
     return img_array_full
 
@@ -54,7 +45,9 @@ def predict_all_cells(histo_img, cellpose_df, _model_SDH):
         predicted_class_array[index_counter] = prediction_result.argmax()
         predicted_proba_array[index_counter] = np.amax(prediction_result)
         index_counter += 1
-    return predicted_class_array, predicted_proba_array
+    cellpose_df["class_predicted"] = predicted_class_array
+    cellpose_df["proba_predicted"] = predicted_proba_array
+    return cellpose_df
 
 
 def paint_full_image(image_sdh, df_cellpose, class_predicted_all):
@@ -78,43 +71,28 @@ def paint_full_image(image_sdh, df_cellpose, class_predicted_all):
 
 
 def run_sdh_analysis(image_array, model_SDH, mask_cellpose):
-    props_cellpose = regionprops_table(
-        mask_cellpose,
-        properties=[
-            "label",
-            "area",
-            "centroid",
-            "eccentricity",
-            "bbox",
-            "image",
-            "perimeter",
-        ],
-    )
-    df_cellpose = pd.DataFrame(props_cellpose)
-    class_predicted_all, proba_predicted_all = predict_all_cells(
-        image_array, df_cellpose, model_SDH
-    )
-    df_cellpose["class_predicted"] = class_predicted_all
-    df_cellpose["proba_predicted"] = proba_predicted_all
-    count_per_label = np.unique(class_predicted_all, return_counts=True)
-    class_and_proba_df = pd.DataFrame(
-        list(zip(class_predicted_all, proba_predicted_all)),
-        columns=["class", "proba"],
-    )
+    df_cellpose = df_from_cellpose_mask(mask_cellpose)
+
+    df_cellpose = predict_all_cells(image_array, df_cellpose, model_SDH)
+    count_per_label = np.unique(df_cellpose["class_predicted"], return_counts=True)
 
     # Result table dict
     headers = ["Feature", "Raw Count", "Proportion (%)"]
     data = []
-    data.append(["Muscle Fibers", len(class_predicted_all), 100])
+    data.append(["Muscle Fibers", len(df_cellpose["class_predicted"]), 100])
     for elem in count_per_label[0]:
         data.append(
             [
                 labels_predict[int(elem)],
                 count_per_label[1][int(elem)],
-                100 * count_per_label[1][int(elem)] / len(class_predicted_all),
+                100
+                * count_per_label[1][int(elem)]
+                / len(df_cellpose["class_predicted"]),
             ]
         )
     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)
+    full_label_map = paint_full_image(
+        image_array, df_cellpose, df_cellpose["class_predicted"]
+    )
     return result_df, full_label_map, df_cellpose