Feature extraction¶
Measure per-label features and store them as a table.
Measure regionprops features from a segmented image with ngio and skimage, and write
them back as a feature table in the OME-Zarr container. By the end the container holds a
table with one row per label, ready to be read back or aggregated across a plate.
Step 1: write the measurement function¶
Start with the function that does the measuring — here a thin wrapper around
skimage.measure.regionprops_table, taking one image patch, one label patch, and
the region's Roi. It can return a DataFrame or a plain dict of columns; either
way the rows must carry the object id in a label column.
import numpy as np
import pandas as pd
from skimage import measure
from ngio import Roi
def extract_features(image: np.ndarray, label: np.ndarray, roi: Roi) -> pd.DataFrame:
"""Basic feature extraction using skimage.measure.regionprops_table."""
label = label.squeeze(-1) # Remove the channel axis if present
roi_feat_table = measure.regionprops_table(
label_image=label,
intensity_image=image,
properties=[
"label",
"area",
"mean_intensity",
"max_intensity",
"min_intensity",
],
)
return pd.DataFrame(roi_feat_table)
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("CardiomyocyteTinyMip", 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: set up the inputs¶
from ngio.transforms import ZoomTransform
# Take the image to measure
image = ome_zarr.get_image()
# Get the nuclei label
nuclei = ome_zarr.get_label("nuclei")
# Here the image is stored at a higher resolution than the nuclei label
print(f"Image dimensions: {image.dimensions}, pixel size: {image.pixel_size}")
print(f"Nuclei dimensions: {nuclei.dimensions}, pixel size: {nuclei.pixel_size}")
# So resample the label up to the image resolution with a transform
zoom_transform = ZoomTransform(
input_image=nuclei,
target_image=image,
order="nearest", # Nearest-neighbour interpolation, so label ids stay intact
)
Step 4: use the FeatureExtractorIterator to create a feature table¶
measure runs the measurement over every region and joins the results
into a FeatureTable referencing the input label — one call, one table. The
per-region measurements schedule exactly like reduce, so a
mapper=ThreadedMapper("auto") parallelizes them; the join still happens once,
at the end. The iterator writes nothing: storing the table is your explicit
add_table call, where the name, backend and overwrite policy belong.
from ngio.iterators import FeatureExtractorIterator
iterator = FeatureExtractorIterator(
input_image=image,
input_label=nuclei,
label_transforms=[zoom_transform],
axes_order=["y", "x", "c"],
)
# Measure every region and join the per-region results into ONE FeatureTable.
# Pass `mapper=ThreadedMapper("auto")` to fan the measurements out in parallel;
# the join always happens once, at the end. Nothing is written yet.
feat_table = iterator.measure(extract_features)
assert feat_table is not None # a serial run always returns the table
# Storing the table is a separate, explicit step.
ome_zarr.add_table("nuclei_regionprops", feat_table, overwrite=True)
For flows the default join does not fit — a different table type, filtering, or
aggregation — either declare a custom join with with_join(...), or drop down to the loop that
measure replaces:
from ngio.tables import FeatureTable
feat_frames = []
for image_data, label_data, roi in iterator.iter_as_numpy():
print(f"Processing ROI: {roi}")
feat_frames.append(extract_features(image_data, label_data, roi))
# Concatenate the per-region frames into one table
manual_table = FeatureTable(table_data=pd.concat(feat_frames), reference_label="nuclei")
ome_zarr.add_table("nuclei_regionprops_manual", manual_table, overwrite=True)
Measuring with a halo¶
Tiling a large image (by_grid, by_blocks) cuts objects at the tile
edges — a border nucleus would be measured on half its pixels. with_halo
fixes that by reading a margin of context around each tile: the function
receives the grown patches (and the grown roi), so a border object is seen
whole by at least one tile. The price is that every tile that sees it
measures it, so the same label appears more than once. Deduplicating is
your declared join's job, and every normalized row carries two provenance
columns for exactly that: roi_index (the region's global index) and
roi_name. The default join keeps the duplicates as-is; a declared join
picks one row per object — here the one from the tile that saw the largest
piece. (roi_index, roi_name, and _ngio_index are reserved: a
measurement function returning one of them is refused.)
from ngio.tables import Table
def keep_most_complete(results: list[pd.DataFrame]) -> Table:
"""One row per object: the measurement from the tile that saw most of it."""
joined = pd.concat([frame for frame in results if len(frame)])
joined = (
joined.sort_values("area", ascending=False)
.drop_duplicates("label")
.drop(columns=["roi_index", "roi_name"])
.set_index("label")
.sort_index()
)
return FeatureTable(table_data=joined, reference_label="nuclei")
# Four tiles, each reading 32 px of context past its edges: a border nucleus
# is measured whole by every tile that sees it, so its label shows up more
# than once — each row stamped with the `roi_index`/`roi_name` it came from.
tiled = iterator.by_blocks(num_y=2, num_x=2).with_halo(y=32, x=32)
tiled_table = tiled.with_join(keep_most_complete).measure(extract_features)
assert tiled_table is not None
ome_zarr.add_table("nuclei_regionprops_tiled", tiled_table, overwrite=True)
Sanity check: read the table back¶
| label | area | mean_intensity-0 | max_intensity-0 | min_intensity-0 | roi_name | roi_index |
|---|---|---|---|---|---|---|
| 1 | 1360.00 | 184.58 | 268.00 | 125.00 | roi_0 | 0 |
| 2 | 2464.00 | 273.25 | 461.00 | 132.00 | roi_0 | 0 |
| 3 | 1968.00 | 277.29 | 429.00 | 143.00 | roi_0 | 0 |
| 4 | 5120.00 | 279.04 | 413.00 | 118.00 | roi_0 | 0 |
| 5 | 288.00 | 243.32 | 341.00 | 147.00 | roi_0 | 0 |
Plot the features¶
The table is made for exactly this — one column against another (the area converted to µm²), every dot one nucleus:
df = ome_zarr.get_table("nuclei_regionprops").dataframe
area_um2 = df["area"] * image.pixel_size.x**2
fig, ax = plt.subplots(figsize=(6.4, 4.2))
# rasterized=True bakes the ~1500 dots into a small embedded image (set_dpi keeps
# them crisp); the axes and labels stay vector, so the figure embeds light.
ax.scatter(
area_um2,
df["mean_intensity-0"],
s=10,
alpha=0.45,
color="#22a699",
linewidths=0,
rasterized=True,
)
ax.set_xlabel("nucleus area (µm²)")
ax.set_ylabel("mean DAPI intensity")
ax.set_title("one dot per nucleus")
ax.spines[["top", "right"]].set_visible(False)
fig.set_dpi(200)
print(
figure_html(
fig,
alt="Scatter of nucleus area against mean DAPI intensity: one dense "
"cloud around 130 square micrometers, and a tail of small, bright "
"nuclei in the upper left.",
)
)
Most nuclei sit in one cloud; the small, bright ones in the upper left — condensed chromatin, dividing or dying — are the kind of subpopulation you measure to find.
Next steps¶
- Iterators — halos, joins, and the read-only iterators in full.
- Distributed processing — the same measurement split across jobs.
- HCS exploration — aggregate feature tables across a plate.
- Table specifications — how feature tables are stored.