segmentation_lite.py

import os

import datetime
import PIL
import numpy as np
import pandas as pd

from skimage.io import imsave
from skimage.filters import threshold_multiotsu
from skimage.measure import regionprops, label
from skimage.segmentation import mark_boundaries, expand_labels, watershed


PIL.Image.MAX_IMAGE_PIXELS = 10000000000
data_folder = ""
nuclei_marker = "nuclei"
nuclei_expansion = 20
bin_percentage = 100.0


def cell_segmentation(nuclei_img):
    sd_labels = label(nuclei_img)

    if nuclei_expansion >= 0:
        sd_labels = expand_labels(sd_labels, distance=nuclei_expansion)

    np.save(data_folder + '/analysis/segmentation_mask.npy', sd_labels)
    print(">>>   Final segmentation result numpy binary data saved =", datetime.datetime.now().strftime("%H:%M:%S"), flush=True)

    boundaries = np.uint8(mark_boundaries(np.zeros((len(nuclei_img), len(nuclei_img[0]), 3)), sd_labels, color=(0, 1, 0)) * 255)
    rgba = np.dstack([boundaries, np.full(boundaries.shape[:2], 255, dtype=np.uint8)])
    black_mask = (boundaries[:, :, 0] == 0) & (boundaries[:, :, 1] == 0) & (boundaries[:, :, 2] == 0)
    rgba[black_mask, 3] = 0
    imsave(data_folder + '/analysis/segmentation_mask_show.png', rgba)

    return sd_labels


# function to calculate the marker intensities for each cell
def marker_calculation(marker, marker_img, cellLabels, data_table, marker_distribution):
    # applying segmentation mask over the marker image
    marker_img_otsu = (marker_img - np.amin(marker_img)) / (np.amax(marker_img) - np.amin(marker_img))
    marker_img_otsu = np.uint8(marker_img_otsu * 255)
    c_otsu_level = 3
    c_otsu_level_min = 1
    c_otsu_level_max = 2
    if 'weak' in marker_distribution:
        c_otsu_level = 5
        c_otsu_level_min = 2
        c_otsu_level_max = 4
    if 'localized' in marker_distribution:
        c_otsu_level_min = c_otsu_level - 2
        c_otsu_level_max = c_otsu_level - 1
    c_otsu = threshold_multiotsu(marker_img_otsu, c_otsu_level)
    temp_otsu = [c_otsu[c_otsu_level_min - 1], c_otsu[c_otsu_level_max - 1]]
    cell_binarized_threshold = max(0.0, temp_otsu[0] - 1) * bin_percentage / 100.0
    marker_properties = regionprops(cellLabels, marker_img_otsu)
    # obtaining mean intensity per cell
    for cell in marker_properties:
        data_table[cell.label][marker] = cell.intensity_mean
        cell_image = cell.image_intensity
        cell_image = cell_image[(cell_image != 0) & (~np.isnan(cell_image))]
        data_table[cell.label][marker + '_local_90'] = np.quantile(cell_image, 0.9) if len(cell_image) > 0 else 0
        #if (marker == "protein"): print("sample",cell_binarized_threshold,data_table[cell.label][marker + '_local_90'],cell.intensity_mean)
        data_table[cell.label][marker + '_bin_low'] = 1 if data_table[cell.label][marker + '_local_90'] >= cell_binarized_threshold else 0
        data_table[cell.label][marker + '_bin_high'] = 1 if cell.intensity_mean >= cell_binarized_threshold else 0

    print(">>>   Marker " + marker + " calculated =", datetime.datetime.now().strftime("%H:%M:%S"), flush=True)


def pipex_lite(basefolder, nuclei_blobs, nuclei_img, protein_img, protein_type="unknown", protein_distribution="unknown", full_run=True):
    global data_folder
    data_folder = basefolder

    os.makedirs(os.path.join(data_folder, 'analysis'), exist_ok=True)

    cell_labels = cell_segmentation(nuclei_blobs)

    # creating data table with segmented cell information
    cell_properties = regionprops(cell_labels)
    data_table = {}
    for cell in cell_properties:
        data_cell = {}
        data_cell['cell_id'] = cell.label
        data_cell['size'] = cell.area
        data_cell['x'] = int(cell.centroid[1])
        data_cell['y'] = int(cell.centroid[0])
        data_table[cell.label] = data_cell

    marker_calculation('nuclei', nuclei_img, cell_labels, data_table, "unknown")
    marker_calculation('protein', protein_img, cell_labels, data_table, protein_distribution)

    # dumpming data_table in cell_data.csv file
    df = pd.DataFrame.from_dict(data_table, orient='index')
    binarized_marker_columns = []
    measure_markers = ['nuclei', 'protein']
    for marker in measure_markers:
        binarized_marker_columns.append(marker + "_local_90")
        binarized_marker_columns.append(marker + "_bin_low")
        binarized_marker_columns.append(marker + "_bin_high")
    measure_markers.extend(binarized_marker_columns)
    measure_markers.insert(0, 'y')
    measure_markers.insert(0, 'x')
    measure_markers.insert(0, 'size')
    measure_markers.insert(0, 'cell_id')
    df = df.reindex(measure_markers, axis=1)
    df.to_csv(data_folder + '/analysis/cell_data.csv', index=False)

    print(">>>   FCS file saved =", datetime.datetime.now().strftime("%H:%M:%S"), flush=True)

    if full_run:
        for index, row in df[df['protein_bin_low'] == 0].iterrows():
            cell_labels[cell_labels == row['cell_id']] = 0

        rgba = np.zeros((len(cell_labels), len(cell_labels[0]), 4), dtype=np.uint8)
        pos_value_mask = cell_labels > 0
        rgba[pos_value_mask] = [255, 128, 0, 128]
        imsave(data_folder + '/analysis/segmentation_mask_binarized_low_show.png', np.uint8(rgba))

        for index, row in df[df['protein_bin_high'] == 0].iterrows():
            cell_labels[cell_labels == row['cell_id']] = 0

        rgba = np.zeros((len(cell_labels), len(cell_labels[0]), 4), dtype=np.uint8)
        pos_value_mask = cell_labels > 0
        rgba[pos_value_mask] = [255, 0, 0, 128]
        imsave(data_folder + '/analysis/segmentation_mask_binarized_high_show.png', np.uint8(rgba))

    print(">>>   Binarized masks saved =", datetime.datetime.now().strftime("%H:%M:%S"), flush=True)

    print(">>>   End time segmentation =", datetime.datetime.now().strftime("%H:%M:%S"), flush=True)