Interpolation problem for sparse image data?

[LucaMarconato]

While improving the Visium HD notebook after #303 has been addressed I have noticed that there must be a problem with interpolations of images when the target canvas is small, because I can't replicate the plot shown in figure with the commented code (when I tried it shows everything black; if I put a large figsize e.g. (30, 30), I start seeing some white pixels).

To reproduce please run the Visium HD notebook.

[LucaMarconato]

@timtreis I asked Claude to turn this into a self-contained script (runnable with uv run my_script.py. It can save you some credits 😁

# /// script
# requires-python = ">=3.10"
# dependencies = [
#   "spatialdata>=0.7.0",
#   "spatialdata-plot>=0.3.0",
#   "matplotlib",
#   "numpy",
#   "sparse<0.16",
# ]
# ///
"""
Reproducer for https://github.com/scverse/spatialdata-plot/issues/351
"Interpolation problem for sparse image data?"

Root cause
----------
spatialdata-plot loses most non-zero pixels when rendering large, sparse images
onto a small canvas.  Two distinct code paths both exhibit the problem:

  scale=None (default)
    _rasterize_if_necessary() pre-downsamples with spatialdata.rasterize(order=0)
    (nearest-neighbour).  The output grid samples one input pixel per cell, so
    an n× downsampling retains only ~1/n² of the non-zero pixels.

  scale="full"
    Skips pre-rasterisation and passes the full array to ax.imshow().
    Matplotlib's "antialiased" interpolation averages each output pixel over a
    wide kernel; a lone pixel with count=1 is diluted to ≈1/k² where k is the
    downsampling factor.  With norm=Normalize(vmin=0, vmax=0.01) and a binary
    black/white colormap the diluted value falls in the black half ([0, 0.5))
    and nothing is visible.

Synthetic data
--------------
Mimics rasterized_002um from the Visium HD notebook:
  - 4000×4000 pixels (≈ 2 µm/px resolution over an 8 mm tissue section)
  - Very sparse: 200 non-zero pixels out of 16 M total (0.00125 %)
  - Non-zero values are integer-like counts 1–4, as produced by rasterize_bins
"""

import numpy as np
import spatialdata as sd
import spatialdata_plot  # noqa: F401
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap, Normalize
from spatialdata.models import Image2DModel

# ── Synthetic sparse rasterized-bin image ────────────────────────────────────
HEIGHT, WIDTH = 4000, 4000
N_NONZERO = 200
CHANNEL = "AA986860"  # gene name from the Visium HD notebook

rng = np.random.default_rng(42)
dense = np.zeros((1, HEIGHT, WIDTH), dtype=np.float32)
ys = rng.integers(0, HEIGHT, N_NONZERO)
xs = rng.integers(0, WIDTH, N_NONZERO)
dense[0, ys, xs] = rng.integers(1, 5, N_NONZERO).astype(np.float32)

image = Image2DModel.parse(dense, dims=("c", "y", "x"), c_coords=[CHANNEL])
sdata = sd.SpatialData(images={"rasterized_002um": image})

# colormap and norm from notebook Cell 34 (binary: black=0, white>0)
cmap_bw = ListedColormap(["#000000", "#ffffff"])
norm = Normalize(vmin=0, vmax=0.01)

# ── Quantitative diagnostic ───────────────────────────────────────────────────
# Simulate _rasterize_if_necessary to measure how many non-zero pixels
# survive the nearest-neighbour pre-downsampling step.
from spatialdata._core.operations.rasterize import rasterize as sd_rasterize
from spatialdata import get_extent

extent = get_extent(image)
nz_source = int(np.count_nonzero(dense))

def survival_after_rasterize(img_da, fig_w_inches: float, fig_h_inches: float, dpi: float = 100) -> int:
    scale = min(dpi * fig_h_inches / img_da.shape[1],
                dpi * fig_w_inches / img_da.shape[2])
    rast = sd_rasterize(
        img_da,
        ("y", "x"),
        [extent["y"][0], extent["x"][0]],
        [extent["y"][1], extent["x"][1]],
        "global",
        target_unit_to_pixels=scale,
    )
    return int(np.count_nonzero(rast.values))

nz_small  = survival_after_rasterize(image, 5,  5)   # figsize=(5,5)
nz_medium = survival_after_rasterize(image, 15, 15)  # figsize=(15,15)
nz_large  = survival_after_rasterize(image, 30, 30)  # figsize=(30,30)

print("── Pixel survival after spatialdata.rasterize(order=0) ─────────────")
print(f"  Source image:              {nz_source:>4} non-zero pixels  ({HEIGHT}×{WIDTH})")
print(f"  scale=None, figsize=(5,5): {nz_small:>4} survive  "
      f"({100*nz_small/nz_source:.1f} %)   8× downsampling")
print(f"  scale=None, figsize=(15,15):{nz_medium:>3} survive  "
      f"({100*nz_medium/nz_source:.1f} %)   2.7× downsampling")
print(f"  scale=None, figsize=(30,30):{nz_large:>3} survive  "
      f"({100*nz_large/nz_source:.1f} %)   1.3× downsampling")
print()

# ── Figure layout: 2×2 grid ──────────────────────────────────────────────────
fig, axes = plt.subplots(2, 2, figsize=(12, 10))
fig.suptitle(
    "Issue #351 — sparse image rendering\n"
    f"Source: {N_NONZERO} non-zero pixels in a {HEIGHT}×{WIDTH} image",
    fontsize=13,
)

# ── Panel 1: Bug A — scale=None (default), small canvas ─────────────────────
# figsize=(5,5) @ 100 dpi → 500 px canvas → 8× downsampling with order=0
# Only ~3 of 200 pixels survive; they appear as sporadic white specks.
inner_fig_a, inner_ax_a = plt.subplots(figsize=(5, 5), dpi=100)
sdata.pl.render_images(
    "rasterized_002um", channel=[CHANNEL], cmap=cmap_bw, norm=norm
).pl.show(ax=inner_ax_a)
inner_ax_a.set_title(
    f"Bug A: scale=None, figsize=(5,5)\n"
    f"{nz_source} source px → {nz_small} survive → nearly all black"
)
inner_fig_a.savefig("bug_A_scale_none.png", dpi=100, bbox_inches="tight")
plt.close(inner_fig_a)

# ── Panel 2: Bug B — scale="full", small canvas ──────────────────────────────
# Full 4000×4000 array passed to ax.imshow(); matplotlib's antialiased
# downsampling dilutes the sparse values below the binary-cmap threshold.
inner_fig_b, inner_ax_b = plt.subplots(figsize=(5, 5), dpi=100)
sdata.pl.render_images(
    "rasterized_002um", channel=[CHANNEL], scale="full", cmap=cmap_bw, norm=norm
).pl.show(ax=inner_ax_b)
inner_ax_b.set_title(
    "Bug B: scale='full', figsize=(5,5)\n"
    "Antialiased avg dilutes count=1 to ≈0.004 < vmax=0.01 → all black"
)
inner_fig_b.savefig("bug_B_scale_full.png", dpi=100, bbox_inches="tight")
plt.close(inner_fig_b)

# ── Panel 3: Reference — plt.spy() shows ALL non-zero pixels ────────────────
# The notebook uses plt.spy() as a workaround (Cell 34).
# spy() plots each non-zero matrix entry as a fixed-size marker,
# circumventing interpolation entirely.
inner_fig_c, inner_ax_c = plt.subplots(figsize=(5, 5), dpi=100)
inner_ax_c.spy(
    sdata["rasterized_002um"].sel(c=CHANNEL).values,
    markersize=1.5,
    color="white",
    markerfacecolor="white",
)
inner_ax_c.set_facecolor("black")
inner_ax_c.set_title(
    f"Reference: plt.spy()\nAll {nz_source} non-zero pixels visible (notebook workaround)"
)
inner_fig_c.savefig("reference_spy.png", dpi=100, bbox_inches="tight")
plt.close(inner_fig_c)

# ── Panel 4: Reference — direct imshow, interpolation='none' ─────────────────
# The other notebook workaround (Cell 34): bypasses spatialdata-plot's
# rasterisation; nearest-neighbour display keeps non-zero values intact.
inner_fig_d, inner_ax_d = plt.subplots(figsize=(5, 5), dpi=100)
inner_ax_d.imshow(
    sdata["rasterized_002um"].sel(c=CHANNEL).values,
    cmap=cmap_bw,
    norm=norm,
    interpolation="none",
)
inner_ax_d.set_title(
    "Reference: direct imshow, interpolation='none'\n"
    "Same aliasing as Bug A — notebook Cell 34 notes this also 'suffers'"
)
inner_fig_d.savefig("reference_imshow_none.png", dpi=100, bbox_inches="tight")
plt.close(inner_fig_d)

print("Saved: bug_A_scale_none.png, bug_B_scale_full.png,")
print("       reference_spy.png, reference_imshow_none.png")

[LucaMarconato]

(Not an important issue to fix, but I stumbled upon this today and I didn't like that it was not reproducible).