Segment Actin Channel - batch process a folder of images

from skimage.io import imread
import pyclesperanto_prototype as cle  # version 0.24.3
import napari_segment_blobs_and_things_with_membranes as nsbatwm  # version 0.3.8
from skimage.segmentation import watershed
import napari_simpleitk_image_processing as nsitk  # version 0.4.5
from skimage.filters import median, sobel
from skimage.morphology import disk
import matplotlib.pyplot as plt
import numpy as np
import tifffile as tif
from tqdm import tqdm
import napari
import os

Helper functions

def seeded_watershed(image, seeds, mask, outline_sigma=0):

    image = np.asarray(image)

    if outline_sigma != 0:
        outline_blurred = gaussian(image, sigma=outline_sigma)
    else:
        outline_blurred = image

    wt =  watershed(outline_blurred, markers=seeds, mask=mask)
    
    return cle.smooth_labels(wt, None, 1.0)

# create napari viewer instance
if 'viewer' not in globals():
    viewer = napari.Viewer()

# a path where nuclei segmentation images were saved
path = "/Users/laura/projects/Bio-image_analysis_school_ScadsAI/prepared_dataset"

# take a look what is in that directory
prepared_images = os.listdir(path)
prepared_images

['aphidicolin_timelapse.tif',
 '.DS_Store',
 'nocodazole_timelapse.tif',
 'AZ-I_timelapse.tif',
 'taxol_timelapse.tif',
 'labels_nuclei',
 'latrunculin B_timelapse.tif',
 'epothilone B_timelapse.tif',
 'monastrol_timelapse.tif',
 'AZ-H_timelapse.tif',
 'DMSO_timelapse.tif',
 'colchicine_timelapse.tif',
 'doxorubicin_timelapse.tif',
 'cytochalasin B_timelapse.tif',
 'labels_tubulin',
 'labels_actin',
 'AZ-A_timelapse.tif']

Process one image as an example

filename = "aphidicolin_timelapse"

img = tif.imread(os.path.join(path, f'{filename}.tif'))

# viewer.add_image(img) # and then in the viewer right click on the layer - split RGB or:
viewer.add_image(
    img, 
    name=["tubulin", "actin", "nuclei"],
    colormap=["magenta", "green", "blue"],
    channel_axis=3)

[<Image layer 'tubulin' at 0x166ad5730>,
 <Image layer 'actin' at 0x166ad5340>,
 <Image layer 'nuclei' at 0x166b67f70>]

# load nuclei segmentation
labels = tif.imread(os.path.join(path, "labels_nuclei", 'aphidicolin_timelapse_labels_dapi.tif'))

# add nuclei segmentation to the viewer
viewer.add_labels(labels, name="nuclei_segmented")

<Labels layer 'nuclei_segmented' at 0x166ad5c40>

# we get actin channel as a numpy array from the viewer
image = viewer.layers['actin'].data

# initialize a new list where we will put our segmented actin images
binarized_imgs = []
labels_actin = []

# iterate through "time" dimension
for t in tqdm(range(image.shape[0])):
    
    # get the binary image
    binary_image = nsitk.threshold_huang(image[t])
    binarized_imgs.append(binary_image)
    
    # noise removal without removing so much of the edges as gaussian blur
    noise_removed_img = cle.median_sphere(image[t], None, 1.0, 1.0, 0.0)
    
    # detect edges
    sobel_edges = sobel(noise_removed_img)
    normalized_sobel = (sobel_edges - sobel_edges.min()) / (sobel_edges.max() - sobel_edges.min())
    
    # perform segmentation with watershed using nuclei as seeds and binary image as a mask
    labels_ws = watershed(normalized_sobel, labels[t], mask=binary_image)
    labels_actin.append(labels_ws)

100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 8/8 [00:02<00:00,  3.75it/s]

timelapse = np.stack(binarized_imgs, axis=0)
viewer.add_labels(timelapse, name="binary_mask")

<Labels layer 'binary_mask' at 0x16a1c62e0>

# note: some layers visibility and/or contrast limits were adjusted interactively in the viewer
screenshot = viewer.screenshot(canvas_only=False)

plt.figure(figsize=(15, 10))
plt.imshow(screenshot)
plt.axis('off')
plt.show()

labels_timelapse = np.stack(labels_actin, axis=0)
viewer.add_labels(labels_timelapse, name="actin_segmented")

<Labels layer 'actin_segmented' at 0x16b675d00>

screenshot = viewer.screenshot(canvas_only=False)

plt.figure(figsize=(15, 10))
plt.imshow(screenshot)
plt.axis('off')
plt.show()

save_path = os.path.join(path, "labels_actin")

if not os.path.exists(save_path):
    os.makedirs(save_path)
    
    tif.imwrite(os.path.join(save_path, f'{filename}_labels_actin.tif'), labels_timelapse)

Process a batch of images from a folder

save_path = os.path.join(path, "labels_actin")
nuclei_labels_path = os.path.join(path, "labels_nuclei")

# create save directory if it does not exist yet
if not os.path.exists(save_path):
    os.makedirs(save_path)

for filename in tqdm(os.listdir(path)):
    
    if not filename.endswith("tif"):
        continue
        
    # print(f"Processing image {filename}")
    rgb_image = tif.imread(os.path.join(path, filename))
    
    # take the green (actin) channel only
    img = rgb_image[:, :, :, 1]
    
    # load nuclei labels
    labels = tif.imread(os.path.join(nuclei_labels_path, f'{filename.split(".tif")[0]}_labels_dapi.tif'))
    
    # initialize an empty list where we will put obtained label images
    labels_actin = []
    binarized_imgs = []
    
    # iterate though the "time" dimension of the image
    for t in range(img.shape[0]):
        
        # get the binary image
        binary_image = nsitk.threshold_huang(img[t])
        binarized_imgs.append(binary_image)
        
        # noise removal without removing so much of the edges as gaussian blur
        noise_removed_img = cle.median_sphere(img[t], None, 1.0, 1.0, 0.0)

        # detect edges
        sobel_edges = sobel(noise_removed_img)
        normalized_sobel = (sobel_edges - sobel_edges.min()) / (sobel_edges.max() - sobel_edges.min())

        # perform segmentation with watershed using nuclei as seeds and binary image as a mask
        labels_ws = watershed(normalized_sobel, labels[t], mask=binary_image)
        labels_actin.append(labels_ws)
        
    # timelapse = np.stack(binarized_imgs, axis=0)
    labels_timelapse = np.stack(labels_actin, axis=0)
    
    # save the labels image
    tif.imwrite(os.path.join(save_path, filename.split(".tif")[0] + "_labels_actin.tif"), labels_timelapse)

100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 17/17 [00:20<00:00,  1.18s/it]