Image segmentation¶

Segment an image one field of view at a time.

Segment an OME-Zarr image with ngio and skimage, one field of view at a time, and write the result back as a label. The second half repeats the segmentation inside a mask, so it only runs where you want it to.

Step 1: set up¶

Start with a function that segments an image, using skimage to do the work.

# Setup a simple segmentation function
import numpy as np
import skimage


def otsu_threshold_segmentation(image: np.ndarray) -> np.ndarray:
    """Simple segmentation using Otsu thresholding.

    Note there is no bookkeeping here: the function numbers its objects from
    1 like any segmenter, and keeping ids unique across regions is the
    iterator's job, not the function's.
    """
    threshold = skimage.filters.threshold_otsu(image)
    binary = image > threshold
    return skimage.measure.label(binary).astype(np.uint32)

Step 2: open the OME-Zarr container¶

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("./data").absolute()
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)

Step 3: segment the image¶

Rather than segmenting the image all at once, map the function over its FOVs. Two problems come free with tiling — every FOV numbers its objects from 1, and an object crossing a FOV boundary comes out as two objects — and with_stitch() solves both: each FOV's ids land in a block of their own during the map, the halo overlap is used to merge split objects afterwards, and the surviving ids are renumbered to a dense 1..N.

from ngio.iterators import SegmentationIterator

# Take the image to read from, and the FOV table naming the regions to walk
image = ome_zarr.get_image()
roi_table = ome_zarr.get_roi_table("FOV_ROI_table")

# Derive an empty label image to write the segmentation into
label = ome_zarr.derive_label("new_label", overwrite=True)

# Setup the segmentation iterator
seg_iterator = SegmentationIterator(
    input_image=image,
    output_label=label,
    channel_selection="DAPI",
    axes_order=["z", "y", "x"],
    consolidation_mode="auto",
)
seg_iterator = seg_iterator.with_stitch().product(roi_table)

# Split any remaining time axis, so each step yields one whole ZYX volume
seg_iterator = seg_iterator.by_zyx()

# Each FOV reads a halo past its edge; `with_stitch()` uses that overlap to give
# every FOV its own id block and to merge objects split by a FOV boundary,
# then renumbers everything to a dense 1..N.
seg_iterator = seg_iterator.with_halo(x=16, y=16)

seg_iterator.map(otsu_threshold_segmentation)

# No need to consolidate, the iterator does it automatically after the map

Plot the segmentation¶

rand_cmap = random_label_cmap()
original = image.get_as_numpy(c=0, z=1, axes_order=["y", "x"])

# The data does not fill its uint16 range, so window it on percentiles.
vmin, vmax = np.percentile(original, (1, 99.8))

fig, axs = plt.subplots(2, 1, figsize=(8, 6))
axs[0].set_title("Original image")
axs[0].imshow(original, cmap="gray", vmin=vmin, vmax=vmax)
axs[1].set_title("Final segmentation")
axs[1].imshow(
    label.get_as_numpy(z=1, axes_order=["y", "x"]),
    cmap=rand_cmap,
    interpolation="nearest",
)
for ax in axs:
    ax.axis("off")
fig.tight_layout()
print(figure_html(fig))

Step 4: masked image segmentation¶

Now use a mask to restrict the segmentation to certain areas of the image. Here you create the mask by hand for illustration, but in a real pipeline it would usually come from another segmentation.

# Create a basic mask for illustration purposes
mask = ome_zarr.derive_label("mask", overwrite=True)
mask_data = mask.get_as_numpy(axes_order=["z", "y", "x"])
mask_data[:, 200:-200, 500:2000] = 1
mask_data[:, 200:-200, 3000:-500] = 2
mask_data[:, 600:-600, 1200:-1000] = 0
mask_data[:, 700:-700, 1600:-1500] = 3
mask.set_array(mask_data, axes_order=["z", "y", "x"])
mask.consolidate(mode="auto")

fig, axs = plt.subplots(2, 1, figsize=(8, 6))
axs[0].set_title("Original image")
axs[0].imshow(original, cmap="gray", vmin=vmin, vmax=vmax)
axs[1].set_title("Mask")
axs[1].imshow(
    mask.get_as_numpy(z=1, axes_order=["y", "x"]),
    cmap=rand_cmap,
    interpolation="nearest",
)
for ax in axs:
    ax.axis("off")
fig.tight_layout()
print(figure_html(fig))

Note that the next step rebinds image to the masked image, so the plot below shows the masked image rather than the original one. Masks are not a tile grid, so there is no stitch here; ids stay unique across masks with UniqueLabelsTransform, where each ROI's own label picks the id block.

from ngio.iterators import MaskedSegmentationIterator

# Take a masked image, which carries its masking ROI table with it
image = ome_zarr.get_masked_image(masking_label_name="mask")

# Derive an empty label image to write the segmentation into
label = ome_zarr.derive_label("masked_new_label", overwrite=True)

# Setup the masked segmentation iterator
from ngio.transforms import UniqueLabelsTransform

# Each mask's ids land in a block of their own — the ROI's label picks the
# block, so the transform needs no per-region state and the map can even run
# under a parallel mapper.
seg_iterator = MaskedSegmentationIterator(
    input_image=image,
    output_label=label,
    channel_selection="DAPI",
    axes_order=["z", "y", "x"],
    consolidation_mode="auto",
    output_transforms=[UniqueLabelsTransform(10_000)],
)

# Split any remaining time axis, so each step yields one whole ZYX volume
seg_iterator = seg_iterator.by_zyx()

seg_iterator.map(otsu_threshold_segmentation)

# No need to consolidate, the iterator does it automatically after the map

fig, axs = plt.subplots(2, 1, figsize=(8, 6))
axs[0].set_title("Original image")
axs[0].imshow(
    image.get_as_numpy(c=0, z=1, axes_order=["y", "x"]),
    cmap="gray",
    vmin=vmin,
    vmax=vmax,
)
axs[1].set_title("Final segmentation")
axs[1].imshow(
    label.get_as_numpy(z=1, axes_order=["y", "x"]),
    cmap=rand_cmap,
    interpolation="nearest",
)
for ax in axs:
    ax.axis("off")
fig.tight_layout()
print(figure_html(fig))

Next steps¶

Feature extraction — measure the objects you just segmented.
Stitching — the same mechanism on a plain tile grid, tuned.
Iterators — halos, stitching and parallel mappers in full.
Masked images and labels — read data object-by-object.

Beyond the tutorials¶

The ngio workshop has hands-on marimo notebooks covering containers, images, labels and tables, and the processing iterators. Run them locally with uv, in the browser via molab, or read them as static pages.