Processing¶

In this notebook we will implement a couple of mock image analysis workflows using ngio.

Maximum intensity projection¶

In this workflow we will read a volumetric image and create a maximum intensity projection (MIP) along the z-axis.

step 1: Create a ngff image¶

For this example we will use the following publicly available image

In [1]:

import matplotlib.pyplot as plt

from ngio.core import NgffImage

ngff_image = NgffImage("../../data/20200812-CardiomyocyteDifferentiation14-Cycle1.zarr/B/03/0")

import matplotlib.pyplot as plt from ngio.core import NgffImage ngff_image = NgffImage("../../data/20200812-CardiomyocyteDifferentiation14-Cycle1.zarr/B/03/0")

step 2: Create a new ngff image to store the MIP¶

In [2]:

mip_ngff = ngff_image.derive_new_image("../../data/20200812-CardiomyocyteDifferentiation14-Cycle1.zarr/B/03/0_mip",
                                       name="MIP",
                                       shape=(1, 1, 2160, 5120))

mip_ngff = ngff_image.derive_new_image("../../data/20200812-CardiomyocyteDifferentiation14-Cycle1.zarr/B/03/0_mip", name="MIP", shape=(1, 1, 2160, 5120))

step 3: Run the workflow¶

For each roi in the image, create a MIP and store it in the new image

In [3]:

# Get the source imag
source_image = ngff_image.get_image()
print("Source image loaded with shape:", source_image.shape)

# Get the MIP image
mip_image = mip_ngff.get_image()
print("MIP image loaded with shape:", mip_image.shape)

# Get a ROI table
roi_table = ngff_image.table.get_table("FOV_ROI_table")
print("ROI table loaded with", len(roi_table.rois), "ROIs")

# For each ROI in the table
# - get the data from the source image
# - calculate the MIP
# - set the data in the MIP image
for roi in roi_table.rois:
    print(f" - Processing ROI {roi.infos.get('field_index')}")
    patch = source_image.get_array_from_roi(roi)
    mip_patch = patch.max(axis=1, keepdims=True)
    mip_image.set_array_from_roi(patch=mip_patch, roi=roi)
    
print("MIP image saved")

plt.figure(figsize=(5, 5))
plt.title("Mip")
plt.imshow(mip_image.on_disk_array[0, 0, :, :], cmap="gray")
plt.axis('off')
plt.tight_layout()
plt.show()

# Get the source imag source_image = ngff_image.get_image() print("Source image loaded with shape:", source_image.shape) # Get the MIP image mip_image = mip_ngff.get_image() print("MIP image loaded with shape:", mip_image.shape) # Get a ROI table roi_table = ngff_image.table.get_table("FOV_ROI_table") print("ROI table loaded with", len(roi_table.rois), "ROIs") # For each ROI in the table # - get the data from the source image # - calculate the MIP # - set the data in the MIP image for roi in roi_table.rois: print(f" - Processing ROI {roi.infos.get('field_index')}") patch = source_image.get_array_from_roi(roi) mip_patch = patch.max(axis=1, keepdims=True) mip_image.set_array_from_roi(patch=mip_patch, roi=roi) print("MIP image saved") plt.figure(figsize=(5, 5)) plt.title("Mip") plt.imshow(mip_image.on_disk_array[0, 0, :, :], cmap="gray") plt.axis('off') plt.tight_layout() plt.show()

Source image loaded with shape: (1, 2, 2160, 5120)
MIP image loaded with shape: (1, 1, 2160, 5120)
ROI table loaded with 2 ROIs
 - Processing ROI None
 - Processing ROI None
MIP image saved

step 4: Consolidate the results (!!! Important)¶

In this we wrote the mip image to a single level of the image pyramid. To truly consolidate the results we would need to write the mip to all levels of the image pyramid. We can do this by calling the .consolidate() method on the image.

In [4]:

# Get the MIP image at a lower resolution
mip_image_2 = mip_ngff.get_image(path="2")

image_before_consolidation = mip_image_2.get_array(c=0, z=0)

# Consolidate the pyramid
mip_image.consolidate()

image_after_consolidation = mip_image_2.get_array(c=0, z=0)

fig, axs = plt.subplots(2, 1, figsize=(10, 5))
axs[0].set_title("Before consolidation")
axs[0].imshow(image_before_consolidation, cmap="gray")
axs[1].set_title("After consolidation")
axs[1].imshow(image_after_consolidation, cmap="gray")
for ax in axs:
    ax.axis('off')
plt.tight_layout()
plt.show()

# Get the MIP image at a lower resolution mip_image_2 = mip_ngff.get_image(path="2") image_before_consolidation = mip_image_2.get_array(c=0, z=0) # Consolidate the pyramid mip_image.consolidate() image_after_consolidation = mip_image_2.get_array(c=0, z=0) fig, axs = plt.subplots(2, 1, figsize=(10, 5)) axs[0].set_title("Before consolidation") axs[0].imshow(image_before_consolidation, cmap="gray") axs[1].set_title("After consolidation") axs[1].imshow(image_after_consolidation, cmap="gray") for ax in axs: ax.axis('off') plt.tight_layout() plt.show()

step 5: Create a new ROI table¶

As a final step we will create a new ROI table that contains the MIPs as ROIs. Where we correct the z bounds of the ROIs to reflect the MIP.

In [5]:

mip_roi_table = mip_ngff.table.new("FOV_ROI_table", overwrite=True)

roi_list = []
for roi in roi_table.rois:
    print(f" - Processing ROI {roi.infos.get('field_index')}")
    roi.z_length = 1 # In the MIP image, the z dimension is 1
    roi_list.append(roi)

mip_roi_table.set_rois(roi_list, overwrite=True)
mip_roi_table.write()

mip_roi_table.table

mip_roi_table = mip_ngff.table.new("FOV_ROI_table", overwrite=True) roi_list = [] for roi in roi_table.rois: print(f" - Processing ROI {roi.infos.get('field_index')}") roi.z_length = 1 # In the MIP image, the z dimension is 1 roi_list.append(roi) mip_roi_table.set_rois(roi_list, overwrite=True) mip_roi_table.write() mip_roi_table.table

 - Processing ROI None
 - Processing ROI None

Out[5]:

	x_micrometer	y_micrometer	z_micrometer	len_x_micrometer	len_y_micrometer	len_z_micrometer	x_micrometer_original	y_micrometer_original
FieldIndex
FOV_1	0.0	0.0	0.0	416.0	351.0	1	-1448.300049	-1517.699951
FOV_2	416.0	0.0	0.0	416.0	351.0	1	-1032.300049	-1517.699951

Image segmentation¶

Now we can use the MIP image to segment the image using a simple thresholding algorithm.

In [6]:

# Setup a simple segmentation function

import numpy as np
from matplotlib.colors import ListedColormap
from skimage.filters import threshold_otsu
from skimage.measure import label

rand_cmap = np.random.rand(1000, 3)
rand_cmap[0] = 0
rand_cmap = ListedColormap(rand_cmap)


def otsu_threshold_segmentation(image: np.ndarray, max_label:int) -> np.ndarray:
    """Simple segmentation using Otsu thresholding."""
    threshold = threshold_otsu(image)
    binary = image > threshold
    label_image = label(binary)
    label_image += max_label
    label_image = np.where(binary, label_image, 0)
    return label_image

# Setup a simple segmentation function import numpy as np from matplotlib.colors import ListedColormap from skimage.filters import threshold_otsu from skimage.measure import label rand_cmap = np.random.rand(1000, 3) rand_cmap[0] = 0 rand_cmap = ListedColormap(rand_cmap) def otsu_threshold_segmentation(image: np.ndarray, max_label:int) -> np.ndarray: """Simple segmentation using Otsu thresholding.""" threshold = threshold_otsu(image) binary = image > threshold label_image = label(binary) label_image += max_label label_image = np.where(binary, label_image, 0) return label_image

step 1: Derive a new label image from the MIP image¶

In [7]:

nuclei_image = mip_ngff.label.derive(name="nuclei", overwrite=True)

nuclei_image = mip_ngff.label.derive(name="nuclei", overwrite=True)

step 2: Run the workflow¶

In [8]:

# Get the source imag
source_image = mip_ngff.get_image()
print("Source image loaded with shape:", source_image.shape)

# Get a ROI table
roi_table = mip_ngff.table.get_table("FOV_ROI_table")
print("ROI table loaded with", len(roi_table.rois), "ROIs")

# Find the DAPI channel
dapi_idx = source_image.get_channel_idx(label="DAPI")

# For each ROI in the table
# - get the data from the source image
# - calculate the Segmentation
# - set the data in segmentation image
max_label = 0
for roi in roi_table.rois:
    print(f" - Processing ROI {roi.infos.get('field_index')}")
    patch = source_image.get_array_from_roi(roi, c=dapi_idx)
    segmentation = otsu_threshold_segmentation(patch, max_label)

    # Add the max label of the previous segmentation to avoid overlapping labels
    max_label = segmentation.max()

    nuclei_image.set_array_from_roi(patch=segmentation, roi=roi)

# Consolidate the segmentation image
nuclei_image.consolidate()

print("Segmentation image saved")
fig, axs = plt.subplots(2, 1, figsize=(10, 5))
axs[0].set_title("MIP")
axs[0].imshow(source_image.on_disk_array[0, 0], cmap="gray")
axs[1].set_title("Nuclei segmentation")
axs[1].imshow(nuclei_image.on_disk_array[0], cmap=rand_cmap, interpolation='nearest')
for ax in axs:
    ax.axis('off')
plt.tight_layout()
plt.show()

# Get the source imag source_image = mip_ngff.get_image() print("Source image loaded with shape:", source_image.shape) # Get a ROI table roi_table = mip_ngff.table.get_table("FOV_ROI_table") print("ROI table loaded with", len(roi_table.rois), "ROIs") # Find the DAPI channel dapi_idx = source_image.get_channel_idx(label="DAPI") # For each ROI in the table # - get the data from the source image # - calculate the Segmentation # - set the data in segmentation image max_label = 0 for roi in roi_table.rois: print(f" - Processing ROI {roi.infos.get('field_index')}") patch = source_image.get_array_from_roi(roi, c=dapi_idx) segmentation = otsu_threshold_segmentation(patch, max_label) # Add the max label of the previous segmentation to avoid overlapping labels max_label = segmentation.max() nuclei_image.set_array_from_roi(patch=segmentation, roi=roi) # Consolidate the segmentation image nuclei_image.consolidate() print("Segmentation image saved") fig, axs = plt.subplots(2, 1, figsize=(10, 5)) axs[0].set_title("MIP") axs[0].imshow(source_image.on_disk_array[0, 0], cmap="gray") axs[1].set_title("Nuclei segmentation") axs[1].imshow(nuclei_image.on_disk_array[0], cmap=rand_cmap, interpolation='nearest') for ax in axs: ax.axis('off') plt.tight_layout() plt.show()

Source image loaded with shape: (1, 1, 2160, 5120)
ROI table loaded with 2 ROIs
 - Processing ROI None
 - Processing ROI None

Segmentation image saved

In [ ]: