Image Processing¶

This is a minimal example of how to use the ngio library for applying some basic image processing techniques.

For this example we will apply gaussian blur to an image.

Step 1: Setup¶

We will first create a simple function to apply gaussian blur to an image. This function will take an image and a sigma value as input and return the blurred image.

In [1]:

import numpy as np
import skimage


def gaussian_blur(image: np.ndarray, sigma: float) -> np.ndarray:
    """Apply gaussian blur to an image."""
    original_type = image.dtype
    image = skimage.filters.gaussian(
        image, sigma=sigma, channel_axis=0, preserve_range=True
    )
    # Convert the image back to the original type
    image = image.astype(original_type)
    return image

import numpy as np import skimage def gaussian_blur(image: np.ndarray, sigma: float) -> np.ndarray: """Apply gaussian blur to an image.""" original_type = image.dtype image = skimage.filters.gaussian( image, sigma=sigma, channel_axis=0, preserve_range=True ) # Convert the image back to the original type image = image.astype(original_type) return image

Step 2: Open the OmeZarr container¶

In [2]:

from pathlib import Path

from ngio import open_ome_zarr_container
from ngio.utils import download_ome_zarr_dataset

# Download the dataset
download_dir = Path(".").absolute().parent.parent / "data"

hcs_path = download_ome_zarr_dataset("CardiomyocyteTiny", download_dir=download_dir)
image_path = hcs_path / "B" / "03" / "0"

# Open the ome-zarr container
ome_zarr = open_ome_zarr_container(image_path)

from pathlib import Path from ngio import open_ome_zarr_container from ngio.utils import download_ome_zarr_dataset # Download the dataset download_dir = Path(".").absolute().parent.parent / "data" hcs_path = download_ome_zarr_dataset("CardiomyocyteTiny", download_dir=download_dir) image_path = hcs_path / "B" / "03" / "0" # Open the ome-zarr container ome_zarr = open_ome_zarr_container(image_path)

Unzipping contents of '/home/runner/work/ngio/ngio/data/20200812-CardiomyocyteDifferentiation14-Cycle1-tiny.zarr.zip' to '/home/runner/work/ngio/ngio/data/tmp'

Step 3: Create a new empty omeZarr container¶

ngio provide a simple way to "derive" a new container from an existing one. This is useful when you want to apply some processing to an image and save the results in a new container that preserves the original metadata and dimensions (unless explicitly changed when deriving).

In [3]:

# First we will need the image object
image = ome_zarr.get_image()

# Second we need to derive a new ome-zarr image where we will store
# the processed image

blurred_omezarr_path = image_path.parent / "0_blurred"
blurred_omezarr = ome_zarr.derive_image(
    store=blurred_omezarr_path, name="Blurred Image", overwrite=True
)
blurred_image = blurred_omezarr.get_image()

# First we will need the image object image = ome_zarr.get_image() # Second we need to derive a new ome-zarr image where we will store # the processed image blurred_omezarr_path = image_path.parent / "0_blurred" blurred_omezarr = ome_zarr.derive_image( store=blurred_omezarr_path, name="Blurred Image", overwrite=True ) blurred_image = blurred_omezarr.get_image()

Step 4: Apply the gaussian blur and consolidate the processed image¶

In [4]:

# We can use the axes order to specify how we query the image data.
# Here we will reorder the axes to be ["c", "z", "y", "x"].
# So that it will be compatible with the gaussian blur function
# which expects the channel axis to be the first one.
image_data = image.get_as_numpy(axes_order=["c", "z", "y", "x"])
# Apply gaussian blur to the image
sigma = 5.0
blurred_image_data = gaussian_blur(image_data, sigma=sigma)

# Set the processed image data back to the ome-zarr image
blurred_image.set_array(patch=blurred_image_data, axes_order=["c", "z", "y", "x"])

# The `set_array` method only saved the blurred image to the container at a specific
# resolution level. So all other resolution levels are still empty.
# To propagate the changes to all resolution levels,
# we can use the `consolidate` method.
blurred_image.consolidate()

# We can use the axes order to specify how we query the image data. # Here we will reorder the axes to be ["c", "z", "y", "x"]. # So that it will be compatible with the gaussian blur function # which expects the channel axis to be the first one. image_data = image.get_as_numpy(axes_order=["c", "z", "y", "x"]) # Apply gaussian blur to the image sigma = 5.0 blurred_image_data = gaussian_blur(image_data, sigma=sigma) # Set the processed image data back to the ome-zarr image blurred_image.set_array(patch=blurred_image_data, axes_order=["c", "z", "y", "x"]) # The `set_array` method only saved the blurred image to the container at a specific # resolution level. So all other resolution levels are still empty. # To propagate the changes to all resolution levels, # we can use the `consolidate` method. blurred_image.consolidate()

Plot the results¶

Finally, we can visualize the original and blurred images using matplotlib.

In [5]:

import matplotlib.pyplot as plt
import numpy as np
from matplotlib.colors import ListedColormap

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

fig, axs = plt.subplots(2, 1, figsize=(8, 4))
axs[0].set_title("Original image")
axs[0].imshow(image.get_as_numpy(c=0, z=1, axes_order=["y", "x"]), cmap="gray")
axs[1].set_title("Blurred image")
axs[1].imshow(blurred_image.get_as_numpy(c=0, z=1, axes_order=["y", "x"]), cmap="gray")
for ax in axs:
    ax.axis("off")
plt.tight_layout()
plt.show()

import matplotlib.pyplot as plt import numpy as np from matplotlib.colors import ListedColormap rand_cmap = np.random.rand(1000, 3) rand_cmap[0] = 0 rand_cmap = ListedColormap(rand_cmap) fig, axs = plt.subplots(2, 1, figsize=(8, 4)) axs[0].set_title("Original image") axs[0].imshow(image.get_as_numpy(c=0, z=1, axes_order=["y", "x"]), cmap="gray") axs[1].set_title("Blurred image") axs[1].imshow(blurred_image.get_as_numpy(c=0, z=1, axes_order=["y", "x"]), cmap="gray") for ax in axs: ax.axis("off") plt.tight_layout() plt.show()

Step 5: Out of memory processing¶

Sometimes we want to apply some simple processing to larger than memory images. In this case, we can use the dask library to process the image in chunks. In ngio we can simply query the data as a dask array and apply the desired processing function to it.

In [6]:

from dask import array as da


def dask_gaussian_blur(image: da.Array, sigma: float) -> da.Array:
    """Apply gaussian blur to a dask array."""
    # This will itroduce some edge artifacts at chunk boundaries
    # In a real application, consider using map_overlap to mitigate this
    # With appropriate depth based on sigma
    return da.map_blocks(gaussian_blur, image, dtype=image.dtype, sigma=sigma)


image_dask = image.get_as_dask(axes_order=["c", "z", "y", "x"])
blurred_image_dask = dask_gaussian_blur(image_dask, sigma=sigma)
blurred_image_dask

from dask import array as da def dask_gaussian_blur(image: da.Array, sigma: float) -> da.Array: """Apply gaussian blur to a dask array.""" # This will itroduce some edge artifacts at chunk boundaries # In a real application, consider using map_overlap to mitigate this # With appropriate depth based on sigma return da.map_blocks(gaussian_blur, image, dtype=image.dtype, sigma=sigma) image_dask = image.get_as_dask(axes_order=["c", "z", "y", "x"]) blurred_image_dask = dask_gaussian_blur(image_dask, sigma=sigma) blurred_image_dask

Out[6]:

Array Chunk Bytes 42.19 MiB 10.55 MiB Shape (1, 2, 2160, 5120) (1, 1, 2160, 2560) Dask graph 4 chunks in 3 graph layers Data type uint16 numpy.ndarray

Step 6. Image Processing Iterators¶

ngio provides an alternative way to process large images using iterators. This API is not meant to replace dask but to provide a simple way to iterate over arbitrary regions, moreover it provides a simple way to implement default broadcasting behaviors.

In [7]:

from ngio.experimental.iterators import ImageProcessingIterator

iterator = ImageProcessingIterator(
    input_image=image,
    output_image=blurred_image,
    axes_order=["c", "z", "y", "x"],
)

# After initializing the iterator, the iterator will have created
# will iterate over the entire image.
print(f"Iterator after initialization: {iterator}")

# Iterate over an arbitrary region of interest table
# We can use the product method that performs a cartesian product
# between the iterator and the table.
table = ome_zarr.get_roi_table("FOV_ROI_table")
iterator = iterator.product(table)
print(f"Iterator after product with table: {iterator}")

# We can explicitly set a broadcasting behavior
# For example we can iterate over all zyx planes, and broadcast all the other
# spatial dimensions
iterator = iterator.by_zyx()

# Finally (if needed) we can check if the regions are not-overlapping
iterator.require_no_regions_overlap()
# We can also check if the regions lay on non-overlapping chunks
iterator.require_no_chunks_overlap()

# Now we can map the gaussian blur function to the iterator
iterator.map_as_numpy(lambda x: gaussian_blur(x, sigma=sigma))

# No need to consolidate, the iterator takes care of that
# after all the regions have been processed

from ngio.experimental.iterators import ImageProcessingIterator iterator = ImageProcessingIterator( input_image=image, output_image=blurred_image, axes_order=["c", "z", "y", "x"], ) # After initializing the iterator, the iterator will have created # will iterate over the entire image. print(f"Iterator after initialization: {iterator}") # Iterate over an arbitrary region of interest table # We can use the product method that performs a cartesian product # between the iterator and the table. table = ome_zarr.get_roi_table("FOV_ROI_table") iterator = iterator.product(table) print(f"Iterator after product with table: {iterator}") # We can explicitly set a broadcasting behavior # For example we can iterate over all zyx planes, and broadcast all the other # spatial dimensions iterator = iterator.by_zyx() # Finally (if needed) we can check if the regions are not-overlapping iterator.require_no_regions_overlap() # We can also check if the regions lay on non-overlapping chunks iterator.require_no_chunks_overlap() # Now we can map the gaussian blur function to the iterator iterator.map_as_numpy(lambda x: gaussian_blur(x, sigma=sigma)) # No need to consolidate, the iterator takes care of that # after all the regions have been processed

Iterator after initialization: ImageProcessingIterator(regions=1)
Iterator after product with table: ImageProcessingIterator(regions=2)