Interpolation problem for sparse image data? #351

Activity

1. LucaMarconato commented on May 7, 2026

   

   LucaMarconato

   on May 7, 2026

   MemberAuthor

   More actions

   @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")
2. LucaMarconato commented on May 7, 2026

   

   LucaMarconato

   on May 7, 2026

   MemberAuthor

   More actions

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