Repository navigation
to_polygons() on labels is not pixel-exact: half-pixel shift, cut corners, filled holes, chunk-dependent output #1301
Description
Description
I expected to_polygons(labels) to give polygons that match the label pixels exactly, so that assigning transcripts by point-in-polygon gives the same result as the pixel lookup labels[floor(y), floor(x)]. That doesn't hold. On a synthetic dataset (60 cells, 200k random points), ~3.1% of points were assigned differently.
The root cause seems to be _region_props_to_polygon, which uses skimage.measure.find_contours(mask, 0.5). That traces outlines through points between pixel centres (pixel centre = integer coordinate). But spatialdata itself (e.g. rasterize) and spatialdata-plot use the convention that pixel [i, j] covers [j, j+1) x [i, i+1). Consequences:
- Half-pixel shift + cut corners: a single pixel at
[1, 1]becomesPOLYGON ((1 1.5, 0.5 1, 1 0.5, 1.5 1, 1 1.5))(area 0.5) instead ofbox(1, 1, 2, 2)(area 1). - Holes are filled: the hole contour is added as a second shell of a
MultiPolygoninstead of an interior ring. The result is invalid (Nested shells) and covers the hole. - Nested labels overlap: a label lying in another label's hole is covered by both polygons.
- Diagonally connected pixels of one label become separate diamonds (area 1 instead of 2).
- Output depends on dask chunking: a 2x2 label in one chunk gives an octagon (area 3.5). Split over 2x2 chunks, it gives 4 diamonds touching only at their tips (area 2). The true area is 4 in both cases.
Figure
Labels are filled, to_polygons() output is red, and test points (4x4 per pixel) are coloured by where pixel lookup and point-in-polygon disagree. Blue = inside a labelled pixel but outside its polygon; orange = inside a polygon of another label/background. Green dashed lines = dask chunk edges.
Minimal example (no plotting)
import numpy as np, dask.array as da, shapely
from spatialdata import to_polygons
from spatialdata.models import Labels2DModel
def polys(lab, chunks=None):
arr = da.from_array(lab, chunks=chunks or lab.shape)
return to_polygons(Labels2DModel.parse(arr, dims=("y", "x"))).geometry
# 1. single pixel
lab = np.zeros((3, 3), np.int32); lab[1, 1] = 1
print(polys(lab).iloc[0]) # diamond, area 0.5; expected POLYGON ((1 1, 2 1, 2 2, 1 2, 1 1))
# 2. hole
lab = np.zeros((5, 5), np.int32); lab[1:4, 1:4] = 1; lab[2, 2] = 0
g = polys(lab).iloc[0]
print(g.is_valid, shapely.is_valid_reason(g)) # False, Nested shells
# 3. nested label
lab[2, 2] = 2
gs = polys(lab)
print(gs[1].intersection(gs[2]).area) # > 0: polygons overlap
# 4. diagonal pixels
lab = np.zeros((4, 4), np.int32); lab[1, 1] = lab[2, 2] = 1
print(polys(lab).iloc[0].area) # 1.0, expected 2
# 5. chunking
lab = np.zeros((4, 4), np.int32); lab[1:3, 1:3] = 1
print(polys(lab).iloc[0].area, polys(lab, chunks=(2, 2)).iloc[0].area) # 3.5 vs 2.0, expected 4 and 4Expected behaviour
Pixel-edge polygons with interior rings, the same convention as rasterize/spatialdata-plot, valid geometry, and output that does not depend on chunking. For comparison, rasterio.features.shapes(arr, mask=arr > 0, connectivity=4) followed by a union_all per label gives exact results (0 mismatches on the same 200k points).
Script to reproduce the figure (spatialdata-plot)
"""Minimal examples: spatialdata.to_polygons(labels) is not consistent with the labels' pixels.
Labels are filled, to_polygons() outlines are red. Test points on a regular grid show
pixel lookup labels[floor(y), floor(x)] vs. point-in-polygon:
grey = same result
blue = inside a labelled pixel but outside that label's polygon
orange = inside a polygon of a different label (background or another cell)
"""
import warnings
import dask.array as da
import geopandas as gpd
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import shapely
import spatialdata_plot # noqa: F401
from spatialdata import SpatialData, to_polygons
from spatialdata.models import Labels2DModel, PointsModel
warnings.filterwarnings("ignore")
matplotlib.use("Agg")
CATS = ["agree", "pixel only", "wrong polygon"]
COLORS = {"agree": "lightgrey", "pixel only": "tab:blue", "wrong polygon": "tab:orange"}
def make_example(lab, chunks=None):
labels = Labels2DModel.parse(da.from_array(lab, chunks=chunks or lab.shape), dims=("y", "x"))
shapes = to_polygons(labels)
# Test points: a 4x4 grid inside every pixel
h, w = lab.shape
xs, ys = np.meshgrid(np.arange(0.125, w, 0.25), np.arange(0.125, h, 0.25))
xs, ys = xs.ravel(), ys.ravel()
truth = lab[ys.astype(int), xs.astype(int)]
cats = []
for x, y, t in zip(xs, ys, truth):
hits = {int(i) for i, g in shapes.geometry.items() if g.covers(shapely.Point(x, y))}
if t != 0 and t not in hits:
cats.append("pixel only")
elif hits - {t}:
cats.append("wrong polygon")
else:
cats.append("agree")
df = pd.DataFrame({"x": xs, "y": ys, "check": pd.Categorical(cats, categories=CATS)})
points = PointsModel.parse(df)
sdata = SpatialData(labels={"labels": labels}, shapes={"polygons": shapes}, points={"points": points})
return sdata, shapes
def plot(sdata, shapes, title, ax, chunks=None):
(
sdata.pl.render_labels("labels", fill_alpha=0.6)
.pl.render_shapes("polygons", fill_alpha=0, outline_alpha=1, outline_color="red", outline_width=2)
.pl.render_points("points", color="check", groups=[c for c in CATS if c != "agree"] ,
palette=[COLORS[c] for c in CATS if c != "agree"], size=12)
.pl.show(ax=ax, title=title)
)
from matplotlib.lines import Line2D
ax.legend(handles=[Line2D([], [], color="red", lw=2, label="to_polygons()")]
+ [Line2D([], [], ls="", marker="o", color=COLORS[c], label=c) for c in CATS[1:]],
loc="lower right", fontsize=7, framealpha=0.9)
h, w = sdata["labels"].shape
ax.set_xticks(range(w + 1)); ax.set_yticks(range(h + 1)); ax.grid(True, color="black", lw=0.4, alpha=0.4)
if chunks:
for y in range(0, h + 1, chunks[0]): ax.axhline(y, color="green", ls="--", lw=1.5)
for x in range(0, w + 1, chunks[1]): ax.axvline(x, color="green", ls="--", lw=1.5)
true_area = {int(i): int((np.asarray(sdata["labels"]) == i).sum()) for i in shapes.index}
info = ", ".join(f"{i}: area {g.area:g} (true {true_area[int(i)]}), valid={g.is_valid}" for i, g in shapes.geometry.items())
ax.set_xlabel(info, fontsize=7)
examples = []
# 1. A single pixel becomes a diamond, shifted by half a pixel
lab = np.zeros((3, 3), np.int32); lab[1, 1] = 1
examples.append(("1. single pixel: half-pixel shift + cut corners", lab, None))
# 2. A hole in a label is filled in, and the MultiPolygon is invalid
lab = np.zeros((5, 5), np.int32); lab[1:4, 1:4] = 1; lab[2, 2] = 0
examples.append(("2. hole is filled (invalid MultiPolygon)", lab, None))
# 3. A label inside another label's hole: polygons overlap
lab = np.zeros((5, 5), np.int32); lab[1:4, 1:4] = 1; lab[2, 2] = 2
examples.append(("3. nested label: overlapping polygons", lab, None))
# 4. Diagonally touching pixels of one label
lab = np.zeros((4, 4), np.int32); lab[1, 1] = 1; lab[2, 2] = 1
examples.append(("4. diagonal pixels of one label", lab, None))
# 5. The result depends on dask chunking (chunk edges = green dashed lines)
lab = np.zeros((4, 4), np.int32); lab[1:3, 1:3] = 1
examples.append(("5a. 2x2 label, one chunk", lab, None))
examples.append(("5b. same label, 2x2 chunks", lab, (2, 2)))
fig, axes = plt.subplots(2, 3, figsize=(15, 10))
for ax, (title, lab, chunks) in zip(axes.ravel(), examples):
sdata, shapes = make_example(lab, chunks)
print(title, "->", [g.wkt for g in shapes.geometry])
plot(sdata, shapes, title, ax, chunks)
single_fig, single_ax = plt.subplots(figsize=(5, 5))
plot(sdata, shapes, title, single_ax, chunks)
single_fig.savefig(f"figs/example_{title.split('.')[0]}.png", dpi=120, bbox_inches="tight")
plt.close(single_fig)
fig.tight_layout()
fig.savefig("figs/all_examples.png", dpi=120, bbox_inches="tight")Versions
spatialdata 0.8.0, spatialdata-plot 0.4.3, scikit-image 0.26.0, shapely 2.2.0, geopandas 1.2.0, dask 2026.8.0, Python 3.13.18
Activity