cell_matcher

Match cells between two segmentation outputs from the same sample.

cell_matcher reads two cell segmentations, converts both coordinate sets to microns, optionally transforms the moving segmentation into the reference coordinate system, and performs greedy one-to-one centroid matching within a maximum distance.

The package is designed for:

• comparing two cell segmentations of the same Xenium sample (for example the original Xenium segmentation vs. proseg or Baysor),
• producing matched pairs plus unmatched cells from each input.

The two inputs are referred to as segmentation A (the reference) and segmentation B (the moving segmentation).

Installation

From this repository:

pip install git+https://github.com/cellgeni/cell_matcher.git

All readers and the plotting and image tools are installed with the base package; there are no runtime extras to add.

The package requires Python 3.12 or newer.

Supported Inputs

cell_matcher supports these input formats:

• Parquet vertex tables, including Xenium cell_boundaries.parquet layout.
• Parquet geometry tables with WKB, WKT, or Shapely-like geometry values.
• GeoJSON FeatureCollection files with Polygon or MultiPolygon geometries.
• Xenium cells.zarr.zip files. The reader uses cell_summary centroids in microns and converts integer cell_id values to the standard Xenium string IDs.
• Generic Zarr label images, with one integer label per cell.
• StarDist polygon coordinate arrays saved as .npz or .npy files. These are read from the StarDist coord array returned in the predict_instances details dictionary.

For Parquet vertex tables, the reader detects common column names such as cell_id, vertex_x, and vertex_y. Column names can also be supplied explicitly through InputConfig or CLI flags.

If a Parquet input path is a Xenium bundle directory, cell_matcher searches for cell_boundaries.parquet or nucleus_boundaries.parquet inside it.

StarDist coordinate arrays use image row/column order. cell_matcher converts them to x/y internally. They are treated as pixel coordinates by default, so provide pixel_size_um unless the arrays were already converted to microns.

Coordinate Units

Matching always runs in microns.

• Use units="micron" when coordinates are already physical coordinates.
• Use units="pixel" plus pixel_size_um=(x_um, y_um) for pixel coordinates.
• Use units="auto" for format-specific defaults: most tabular geometry inputs default to microns, while source formats that are inherently pixel-based, such as StarDist arrays and label images, require a pixel size.

No range-based pixel-versus-micron guessing is used.

Python API

from cell_matcher import InputConfig, MatchConfig, match_segmentations

result = match_segmentations(
    reference=InputConfig(
        path="segmentation_a/cell_boundaries.parquet",
        label="segmentation_a",
        units="micron",
    ),
    moving=InputConfig(
        path="segmentation_b/cell_boundaries.parquet",
        label="segmentation_b",
        units="micron",
    ),
    match=MatchConfig(max_distance_um=10.0),
    out_matches="matches.csv",
    out_unmatched_reference="unmatched_segmentation_a.txt",
    out_unmatched_moving="unmatched_segmentation_b.txt",
    verbose=True,
)

print(result.matches)
print(result.unmatched_a)
print(result.unmatched_b)
print(result.summary)

Two segmentations of the same Xenium sample already share a coordinate frame, so no transform is needed. If the inputs are in different coordinate frames, supply an affine transform that maps segmentation B into segmentation A:

result = match_segmentations(
    reference=InputConfig(path="segmentation_a/cell_boundaries.parquet", label="segmentation_a"),
    moving=InputConfig(path="segmentation_b/cell_boundaries.parquet", label="segmentation_b"),
    match=MatchConfig(max_distance_um=8.0, transform_path="registration_affine.csv"),
)

Pass verbose=True to match_segmentations (or --verbose on the command line) to print each stage (reading, validating, transforming, matching) to stderr as it runs, followed by a summary at the end.

Assignment Modes

MatchConfig.mode selects how pairs are chosen:

• greedy (default): accept the shortest available pair first. Fast, but not globally optimal.
• hungarian: minimum total distance assignment. Globally optimal for the distance cost. Best when accuracy matters more than speed.
• mutual_nearest: keep a pair only when each cell is the other's nearest cell. Stricter, fewer false positives, more unmatched cells.
• overlap: match polygons by intersection-over-union (IoU). Use only when both inputs share one coordinate frame and carry polygon geometries. Set min_iou to drop weak overlaps. A transform is not supported in this mode; register the inputs first.

result = match_segmentations(
    reference=InputConfig(path="segmentation_a/cell_boundaries.parquet", label="segmentation_a"),
    moving=InputConfig(path="segmentation_b/cell_boundaries.parquet", label="segmentation_b"),
    match=MatchConfig(mode="hungarian", max_distance_um=8.0),
)

Confidence Metrics

With MatchConfig.compute_confidence=True (default), centroid matches (modes greedy, hungarian, mutual_nearest) include per-match quality columns:

• distance_score: 1 - distance / max_distance_um, clipped to 0-1.
• second_distance_um: distance from the reference cell to its second-nearest moving cell.
• margin_um: second_distance_um - distance_um. A small margin means an ambiguous match.
• ambiguous_candidate_count: number of moving cells within max_distance_um of the reference cell.
• iou_score: polygon intersection over union (Jaccard) for the matched pair.
• iom_score: polygon intersection over the smaller area (Szymkiewicz–Simpson).

iou_score and iom_score need polygons for both cells. When compute_confidence=True, the readers load boundaries automatically (parquet vertex/geometry tables, GeoJSON, Xenium cells.zarr.zip polygons, or traced label-image contours), and the moving polygons are carried through the transform so the scores stay aligned. The scores are NaN only when the input has no boundaries to load (for example centroid-only cells.parquet).

Overlap matches (mode overlap) include iou, second_best_iou, iou_margin, and overlap_candidate_count.

Validation

Each input is checked before matching. MatchConfig.on_duplicate_ids controls the response to duplicate ids (warn, raise, or ignore). The result summary carries a report per input under reference_validation and moving_validation with counts for duplicate ids, non-finite centroids, empty geometries, and the resolved units.

Command Line

cell_matcher \
  --reference segmentation_a/cell_boundaries.parquet \
  --reference-label segmentation_a \
  --moving segmentation_b/cell_boundaries.parquet \
  --moving-label segmentation_b \
  --max-distance-um 10 \
  --out-matches matches.csv \
  --out-unmatched-reference unmatched_segmentation_a.txt \
  --out-unmatched-moving unmatched_segmentation_b.txt \
  --verbose

Add --transform registration_affine.csv only when the two segmentations are in different coordinate frames.

Useful CLI options:

• --reference-format and --moving-format: auto, parquet, geojson, zarr, or stardist.
• --reference-units and --moving-units: auto, micron, or pixel.
• --reference-pixel-size and --moving-pixel-size: pixel size as X_UM Y_UM.
• --invert-transform: invert the supplied affine before applying it.
• --mode: greedy, hungarian, mutual_nearest, or overlap.
• --min-iou: minimum IoU for a candidate pair in overlap mode.
• --no-confidence: skip per-match confidence columns.
• --on-duplicate-ids: warn, raise, or ignore.
• --reference-id-column, --reference-x-column, --reference-y-column, and matching --moving-* flags for explicit column names.

Transform Format

The transform maps moving centroids into the reference coordinate system. It is applied after unit conversion to microns.

Supported transform files:

• CSV containing a numeric 3x3 matrix,
• .npy file containing a 3x3 array,
• JSON file containing a 3x3 nested list.

When no transform is supplied, the identity matrix is used.

Creating a Transform

cell_matcher expects a 3x3 matrix that maps moving coordinates to reference coordinates after both inputs have been converted to microns:

import numpy as np

moving_h = np.array([x_moving_um, y_moving_um, 1.0])
x_reference_um, y_reference_um, _ = transform @ moving_h

If a registration tool gives you a pixel-space matrix, convert it to microns before passing it to cell_matcher:

import numpy as np


def pixel_affine_to_micron_affine(
    matrix_px,
    moving_pixel_size_um,
    reference_pixel_size_um,
):
    moving_pixel_size_um = np.asarray(moving_pixel_size_um, dtype=float)
    reference_pixel_size_um = np.asarray(reference_pixel_size_um, dtype=float)

    moving_um_to_px = np.diag([1 / moving_pixel_size_um[0], 1 / moving_pixel_size_um[1], 1])
    reference_px_to_um = np.diag([reference_pixel_size_um[0], reference_pixel_size_um[1], 1])

    return reference_px_to_um @ matrix_px @ moving_um_to_px

Palom

Palom can estimate a global affine from thumbnails and then scale that affine to the selected pyramid level. Use aligner.affine_matrix when you want the global moving-to-reference affine. Do not pass block_affine_matrices_da to cell_matcher; that is a local, piecewise transform and this package currently accepts one global affine.

import numpy as np
import palom


def level_pixel_size_um(reader, level):
    pixel_size = np.asarray(reader.pixel_size, dtype=float)

    if pixel_size.ndim == 0:
        pixel_size = np.repeat(pixel_size, 2)

    return tuple(pixel_size * reader.level_downsamples[level])


LEVEL = 1
THUMBNAIL_LEVEL = 3

reference_reader = palom.reader.OmePyramidReader("reference.ome.tif")
moving_reader = palom.reader.OmePyramidReader("moving.ome.tif")

aligner = palom.align.Aligner(
    ref_img=reference_reader.read_level_channels(LEVEL, 0),
    moving_img=moving_reader.read_level_channels(LEVEL, 0),
    ref_thumbnail=reference_reader.read_level_channels(THUMBNAIL_LEVEL, 0).compute(),
    moving_thumbnail=moving_reader.read_level_channels(THUMBNAIL_LEVEL, 0).compute(),
    ref_thumbnail_down_factor=(
        reference_reader.level_downsamples[THUMBNAIL_LEVEL] / reference_reader.level_downsamples[LEVEL]
    ),
    moving_thumbnail_down_factor=(
        moving_reader.level_downsamples[THUMBNAIL_LEVEL] / moving_reader.level_downsamples[LEVEL]
    ),
)

aligner.coarse_register_affine(n_keypoints=4000)
matrix_px = aligner.affine_matrix

matrix_um = pixel_affine_to_micron_affine(
    matrix_px,
    moving_pixel_size_um=level_pixel_size_um(moving_reader, LEVEL),
    reference_pixel_size_um=level_pixel_size_um(reference_reader, LEVEL),
)
np.savetxt("palom_affine_um.csv", matrix_um, delimiter=",")

Then use it as the moving-to-reference transform:

from cell_matcher import InputConfig, MatchConfig, match_segmentations

result = match_segmentations(
    reference=InputConfig(path="reference/cells.zarr.zip", label="reference"),
    moving=InputConfig(path="moving/cells.zarr.zip", label="moving"),
    match=MatchConfig(
        max_distance_um=10,
        transform_path="palom_affine_um.csv",
    ),
)

scikit-image registration

For translation-only registration, use skimage.registration.phase_cross_correlation. It returns shifts in array axis order (row, column), so the affine translation is (x=column, y=row).

import numpy as np
from skimage.registration import phase_cross_correlation

shift_rc, error, phase = phase_cross_correlation(
    reference_image,
    moving_image,
    upsample_factor=20,
    normalization=None,
)

shift_row, shift_col = shift_rc
matrix_px = np.array(
    [
        [1, 0, shift_col],
        [0, 1, shift_row],
        [0, 0, 1],
    ],
    dtype=float,
)

matrix_um = pixel_affine_to_micron_affine(
    matrix_px,
    moving_pixel_size_um=(0.5, 0.5),
    reference_pixel_size_um=(0.5, 0.5),
)
np.savetxt("skimage_translation_um.csv", matrix_um, delimiter=",")

For an affine estimated from matched landmarks, use points in x/y order and estimate a transform from moving points to reference points:

import numpy as np
from skimage.transform import AffineTransform

moving_points_xy = np.array(
    [
        [100, 200],
        [250, 210],
        [120, 360],
        [300, 390],
    ],
    dtype=float,
)
reference_points_xy = np.array(
    [
        [115, 188],
        [265, 196],
        [139, 348],
        [320, 375],
    ],
    dtype=float,
)

tform = AffineTransform.from_estimate(moving_points_xy, reference_points_xy)

if not tform:
    raise RuntimeError("Affine estimation failed.")

matrix_um = pixel_affine_to_micron_affine(
    tform.params,
    moving_pixel_size_um=(0.5, 0.5),
    reference_pixel_size_um=(0.5, 0.5),
)
np.savetxt("skimage_affine_um.csv", matrix_um, delimiter=",")

Outputs

The returned MatchResult contains:

• matches: a DataFrame with cell_id_<reference_label>, cell_id_<moving_label>, and distance_um.
• unmatched_a: reference cell IDs with no match.
• unmatched_b: moving cell IDs with no match.
• summary: counts and matched-distance statistics.
• reference_config, moving_config, match_config: the configs used.
• reference_cells, moving_cells: the cell sets used for matching (moving is after transform), so plots need only the result.

Optional output files can be written for matches and unmatched IDs.

Saving Results

Save the full result (matches, cell sets, and configs) to one file with pickle:

result.save("result.pkl")  # a ".pkl" suffix is added if omitted
loaded = MatchResult.load("result.pkl")

Visualization

The plotting helpers live in cell_matcher.plot. Each takes only the MatchResult, because the result carries the cell sets and the distance cutoff:

from cell_matcher import InputConfig, MatchConfig, match_segmentations
from cell_matcher.plot import centroid_assignments, match_overview

reference = InputConfig(path="segmentation_a/cells.zarr.zip", label="segmentation_a")
moving = InputConfig(path="segmentation_b/cells.zarr.zip", label="segmentation_b")
match = MatchConfig(max_distance_um=10)

result = match_segmentations(reference=reference, moving=moving, match=match)

match_overview(result, save_path="match_overview.png")
centroid_assignments(result, save_path="assignment_panels.png")

centroid_assignments draws three panels that share axis limits:

1. reference: reference cells are green when matched, grey when not.
2. moving: moving cells are green when matched, grey when not.
3. overlap: matched cells are green, reference-only and moving-only cells get distinct colors, so the direction of each mismatch is visible.

Points are small and semi-transparent so dots do not fully hide each other.

cell_closeup zooms in on one cell:

from cell_matcher.plot import cell_closeup

cell_closeup(result, "cell-123", save_path="cell_closeup.png")

It draws the target cell boundary and centroid at full resolution, the nearby cells from both segmentations, and the matched cell, with a connector line. A text box lists the match metrics for the pair (distance_um, match_score, second_distance_um, margin_um, ambiguous_candidate_count, or the overlap equivalents). source selects reference or moving; it is inferred from the cell id when omitted.

Boundaries are drawn only when the cell sets carry polygons. In centroid mode the moving boundaries are dropped during the transform, so by default only the match centroid is drawn. Pass show_segmentation=True to re-read the inputs from result.reference and result.moving, fetch the full-resolution boundaries, and apply the same unit conversion and transform, so the match boundary lines up with the reference:

cell_closeup(result, "cell-123", show_segmentation=True, save_path="cell_closeup.png")

show_segmentation fetches boundaries from parquet and zarr inputs. For parquet it uses the vertex or geometry columns; for a zarr.zip it reads the Xenium polygon vertices, and for a generic label-image zarr it traces each label's contour into a polygon (pixel coordinates are scaled to microns). If the source holds only centroids (for example Xenium cells.parquet) there are no polygons to fetch, so cell_closeup warns and draws centroids only. Point it at boundary files instead:

cell_closeup(
    result,
    "cell-123",
    reference_segmentation="segmentation_a/cell_boundaries.parquet",
    moving_segmentation="segmentation_b/cell_boundaries.parquet",
    save_path="cell_closeup.png",
)

Pass reference_cells, moving_cells, or max_distance_um as keyword overrides when you want to plot custom cell sets or a different cutoff.

cell_closeup can also draw a background image under the cells with image. It is treated as Xenium pixel data (origin top-left, y down) and mapped to microns with image_pixel_size_um. Leave image_pixel_size_um=None (default) to read the pixel size from the image metadata (OME-TIFF PhysicalSizeX/Y, or the TIFF resolution tags), falling back to the Xenium default of 0.2125 um/pixel (a notice is printed whenever this fallback is used); pass a scalar or (x, y) to override. image may be a numpy array or any image path; for a TIFF/OME-TIFF it is opened lazily through tifffile's zarr store and only the crop covering the zoom window is read into memory, so the full image is never loaded. tifffile assembles a multi-file OME-TIFF such as the Xenium morphology_focus set into the full multi-channel stack. image_channel selects the stain plane (0 = DAPI). When an image is given (and no ax is supplied), cell_closeup draws two panels side by side: outlines only on the left (as without an image), and the image in the background with bright, full-alpha boundary lines drawn over a black halo so they contrast against the image. It then returns (fig, (ax_left, ax_right)).

For a false-color composite, pass an RGB array, or a callable image(row0, row1, col0, col1) that returns the crop for those pixel bounds. The callable composites only the zoom-window crop at full resolution, which suits blending the morphology_focus channels into RGB:

CH_COLORS = np.array([[0, 0.4, 1.0], [0, 1.0, 0.3], [1.0, 0.5, 0], [1.0, 0.2, 0.6]], np.float32)
channel_files = [
    "ch0000_dapi.ome.tif",
    "ch0001_atp1a1_cd45_e-cadherin.ome.tif",
    "ch0002_18s.ome.tif",
    "ch0003_alphasma_vimentin.ome.tif",
]

def composite(row0, row1, col0, col1):
    rgb = np.zeros((row1 - row0, col1 - col0, 3), np.float32)
    for path, color in zip(channel_files, CH_COLORS):
        with tifffile.TiffFile(f"morphology_focus/{path}", is_ome=False) as tif:
            plane = tif.series[0].levels[0].asarray()[row0:row1, col0:col1].astype(np.float32)
        lo, hi = np.percentile(plane, (1, 99))
        norm = np.clip((plane - lo) / (hi - lo + 1e-6), 0, 1)
        rgb += norm[..., None] * color
    return np.clip(rgb, 0, 1)

cell_closeup(result, "cell-123", image=composite, image_pixel_size_um=0.2125,
             show_segmentation=True)

image_pixel_size_um is the full-resolution pixel size, so the callback reads full-resolution channel crops (levels[0]). Percentile normalization here is per-crop; precompute global lo/hi in the closure if you want fixed contrast.

cell_closeup(
    result,
    "cell-123",
    image="morphology_focus",  # folder, an .ome.tif, or a numpy array
    image_channel=0,  # 0 = DAPI
    show_segmentation=True,
    save_path="cell_closeup.png",
)

Reading TIFF images uses tifffile and zarr, both core dependencies, so lazy region reads work out of the box.

The overview plot shows assignment counts, the matched-distance distribution, and distance-ranked pairs. When the matches carry overlap scores it also adds an IoU panel and an IoM panel (iou_score/iom_score for centroid modes, or the iou column for overlap mode), each a 0-1 histogram with the median marked.

Development

Run the test suite with:

uv run --extra dev pytest -q

Format and lint the code with Ruff:

uv run --extra dev ruff format .
uv run --extra dev ruff check .

The tests cover centroid matching, Parquet readers, GeoJSON reading, Zarr label images, Xenium cells.zarr.zip files, StarDist coordinate arrays, affine transforms, unit conversion, visualization helpers, end-to-end matching, and the CLI.