Offset when rasterizing data

When rasterizing images or labels (basically resampling), there is an offset (see picture) that I think is due to the fact that somewhere in the code for rasterization the coordinates of the pixels are considered to be associated with the top left corner, and in other parts of the code for rasterization the coordinates refer to the center of the pixels.

This half pixel difference is not negligible since when we align together images of different resolutions, one pixel for one image can be many pixels for another.

The screenshot is from this example:

spatialdata/examples/dev-examples/spatial_query_and_rasterization.py

Line 15 in eb99aeb

def _visualize_crop_affine_labels_2d() -> None:

Another important point, the rasterization functions could be used for computing aggregation functions. Until this offset problem is addressed, the aggregation would be imprecise.

The function transform_to_data_extent() (which is coming soon in a new PR) is affected by this bug. I had to create a test with permissive tolerance, which reads as assert np.allclose(a, first_a, rtol=0.005).

When this bug is addressed we should reduce the tolerance.

The bug is not just of rasterize, but multiple factors are involved. First, napari considers the coordinates as the center of pixels. Here is for instance the effect of rotating an image by -45 degrees, and showing also the (0, 0) coordinate before and after rotation. Note that the rotation doesn't preserve the location of the (0, 0)-pixel top-left corner.

I need to check, but I suspect that ndinterp() used in transform() and rasterize() may not operate as in napari, leading to some misalignment.

I confirm that one of the reasons for the bug is the fact that (0, 0) is not the top-left corner but the center of the top-left pixel, as shown in the image above. I am fixing this in #453. Nevertheless, there is also another problem, apparently due to ndinterp.affine_transform(): scipy/scipy#20038.

@timtreis

After digging deeper into this I found out the causes of the bug, some of which can be fixed (as I did partially did in #453), but one of them requires a design choice to be made. I will explain this here.

The bugs are the following (they all lead to an error to up to 1 pixel, which is noticeable when the pixel sizes are large):

ndinterp.affine_transform() adds a black border padding. Tracked here BUG: ndinterp.affine_transform() creates an image with black border, in each direction scipy/scipy#20038. This is probably due to the interpolation algorithm, even if theoretically, since we don't use interpolation, this could be avoided. A workaround (computationally expensive) could be an optional flag that induces a padding of the image before calling ndinterp.affine_transform(), and followed by a appropriate cropping. Anyway, even if there is a noticeable border padding added, the alignment of the non-black pixels is not compromised.
The bounding box query function adds the wrong translation vector, this is fixed in Bug ndinterp.affine_transform() non-pixel perfect #453 (not merged).
The functions rasterize() and _transform_raster(), which both call ndinterp.affine_transform(), need to adjust the transformation matrix to account to how a raster image is positioned with respect to the coordinates of the cartesian axes. Practically, once the point 4. below will be addressed, we have to adjust the implementation. Also, the offset that is added to the rasterized image (in rasterize()) and to the transformed image (in _set_transformation_for_transformed_elements()), needs to take into account for this).

The design choice to be made is due to the following:
4. We have to choose where the center of the pixel (0, 0) is, with respect to the cartesian point (0, 0). Napari chooses to (a) have the center of the pixel (0, 0) in the point (0, 0), I propose instead to (b) have it in (0.5, 0.5).

@timtreis @Sonja-Stockhaus I kindly ask for feedback, how does matplotlib operates? With the approach (a) or (b)?

Some example implications of the case (a) are the following:

get_extent() would have to be modified to return a bounding box that has origin in (-0.5, -0.5) and applying an operation to an image would lead to a new different origin (for instance a 10x scale would lead to (-5, -5). I find this unintuitive;
If an image has a 10x scale operation, calling transform() would return the scaled image with a translation vector attached, to account for the offset between the old bounding box origin and the new bounding box origin. More precisely, this translation vector is the sum of 3 vectors: the offset between the (0, 0) point and the image origin pre-transformation; the offset between the (0, 0) point and the image origin post-transformation; the offset between the transformed (0, 0) point and the transformed origin (the last two numbers don't coincide: the former depends only on the resolution of the transformed image and in particular takes into account for the "corner" padding in case of rotations, the second only on the transformation itself). In other words, if before the transformation we have 1 pixel = 2x3 µm, now we may have 1 pixel = 20x30µm (which is taken into account by the latter translation) but the number of pixels in the new image to represent this could be something different, always when the pixel is rotated (which is taken into account by the former translation). The implementation around this is very tedious.

Some implications of the case (b) are the following:

get_extent() keeps returning (0, 0)
If you apply transform to an image, the resulted transformation doesn't require the addition of any translation vector. The implementation becomes much much easier.
If you make a napari layer passing affine=my_transformation.to_affine_matrix(...), you get the wrong result; you need to use napari-spatialdata instead, which will add a tiny translation to account for the different interpretation of pixel centers.

The approach (a) is motivated by http://alvyray.com/Memos/CG/Microsoft/6_pixel.pdf. After reading this (just done) I am more inclined for this approach. We should carefully what this implies for the result of the transform() and rasterize() APIs, and how matplotlib operates before continuing.

Using extent and origin I can pretty much place the pixel anywhere I want but the default also seems to be the center of the pixel, like Napari. I personally also find this more intuitive.

import matplotlib.pyplot as plt
import numpy as np

matrix = np.array([[2]])

plt.imshow(matrix, cmap='gray')
plt.grid(color='red', linestyle='-', linewidth=2)
plt.xticks(np.arange(-1, 2))
plt.yticks(np.arange(-1, 2))

plt.show()

The tests for get_centroids() (get_centroids() #465) needs to be adjusted when this issue is solved.

The new "shapes labels interchangeability" notebook which shows rasterize() is also affected by this bug.

the atol in the test in test_rasterize_bins() can be set to a smaller value after this bug is fixed.

@t-a-m-i

When the issue is addressed, we should double check against this example from @t-a-m-i that all is working correctly: https://scverse.zulipchat.com/#narrow/channel/315824-spatial/topic/SpatialData/near/479568220

Hi, did any fix happen to get implemented? I have a Xenium dataset where the cell polygons all seem increasingly offset in lower resolution images. But it's a bit unclear to me if this is just a natural consequence of aliasing or some kind of offset error I should manually fix (mentioned in the thread above).

Hi, unfortunately a fix is not available yet because we will change the transformation system soon, and the problem should not appear anynore.

This whole message is AI-generated. The issue was automatically discovered and reported by an AI agent (Claude) during an autonomous bug hunt on the spatialdata code base. It has not been verified or triaged by a human yet; the needs: triage label is set so that a maintainer can confirm it. The reproduction script below was executed by the agent in an isolated environment (see Environment) and its output is pasted verbatim.

I re-attempted to reproduce this issue on the current main branch (commit ccf1ea048d054b6624214bf618008a9f9ae223e0) as part of a triage pass over open issues, and confirmed it reproduces.

Reproduction script (PEP 723 inline metadata; run with uv run repro.py, no other setup needed):

# /// script
# requires-python = ">=3.12"
# dependencies = [
#     "spatialdata @ git+https://github.com/scverse/spatialdata.git@ccf1ea048d054b6624214bf618008a9f9ae223e0",
# ]
# ///
"""Issue #165: half-pixel offset when rasterizing raster data at a different resolution.
Same probe as scverse/spatialdata#1220 (findings/issues/19): a single bright pixel is rasterized at
coarser resolutions and the centre of the bright output pixel is mapped back to global coordinates."""
import warnings

import numpy as np
import pandas as pd

from spatialdata import get_extent, rasterize, transform
from spatialdata.models import Image2DModel, PointsModel

warnings.simplefilter("ignore")
arr = np.zeros((1, 100, 100), dtype=np.float32)
arr[0, 40, 60] = 1.0  # bright pixel covering x in [60, 61], y in [40, 41]
img = Image2DModel.parse(arr)
bug = False
for units_per_pixel in [1.0, 2.0, 4.0]:
    out = rasterize(img, axes=("x", "y"), min_coordinate=[0, 0], max_coordinate=[100, 100],
                    target_coordinate_system="global", target_unit_to_pixels=1 / units_per_pixel)
    o = np.asarray(out.data.compute())[0]
    ys, xs = np.nonzero(o)
    centre_px = pd.DataFrame({"x": xs.astype(float) + 0.5, "y": ys.astype(float) + 0.5})
    centre_units = transform(PointsModel.parse(centre_px, transformations={"global": out.attrs["transform"]["global"]}),
                             to_coordinate_system="global").compute()
    expected_x = (xs[0] + 0.5) * units_per_pixel
    shift = centre_units.x.iloc[0] - expected_x
    ext = get_extent(out)["x"]
    print(f"output pixel size s={units_per_pixel}: bright output pixel x={xs[0]}, centre maps to x={centre_units.x.iloc[0]:.2f} "
          f"(expected {expected_x:.2f}), shift={shift:+.2f}; extent x=({ext[0]:.1f}, {ext[1]:.1f}) (expected (0, 100))")
    bug |= abs(shift) > 1e-9
print("VERDICT:", "REPRODUCED (offset still present; duplicate of #1220)" if bug else "NOT REPRODUCED")

Environment: uv run --no-project --python 3.13, which resolves spatialdata from git+https://github.com/scverse/spatialdata.git@ccf1ea048d054b6624214bf618008a9f9ae223e0 (extra dependencies, if any, are pinned in the script's inline metadata above).

Output (verbatim, only uv package-resolution/install log lines removed):

output pixel size s=1.0: bright output pixel x=60, centre maps to x=60.50 (expected 60.50), shift=+0.00; extent x=(0.0, 100.0) (expected (0, 100))
output pixel size s=2.0: bright output pixel x=30, centre maps to x=60.50 (expected 61.00), shift=-0.50; extent x=(-0.5, 99.5) (expected (0, 100))
output pixel size s=4.0: bright output pixel x=15, centre maps to x=60.50 (expected 62.00), shift=-1.50; extent x=(-1.5, 98.5) (expected (0, 100))
VERDICT: REPRODUCED (offset still present; duplicate of #1220)
exit=0

Agent's severity assessment: medium (see the priority: * label added to this issue).

Automatically generated; discovered by an AI agent (Claude) and not yet reviewed by a human.

Navigation Menu

Repository navigation

Offset when rasterizing data #165

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Offset when rasterizing data #165

Description

Activity

LucaMarconato commented on Mar 1, 2023

LucaMarconato commented on Jan 30, 2024

LucaMarconato commented on Feb 6, 2024

LucaMarconato commented on Feb 6, 2024

LucaMarconato commented on Feb 12, 2024

LucaMarconato commented on Feb 12, 2024

timtreis commented on Feb 12, 2024

LucaMarconato commented on Feb 22, 2024

LucaMarconato commented on May 28, 2024

LucaMarconato commented on Jun 14, 2024

LucaMarconato commented on Dec 13, 2024

kushalnimkar commented on Jan 8, 2026

LucaMarconato commented on Jan 9, 2026

LucaMarconato commented on Sep 13, 2026

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions