Transformation on large image results in very large chunks

[yyi17]

I am analyzing the 10x Genomics Xenium example datasets.

https://cf.10xgenomics.com/samples/xenium/4.0.0/Human_Breast_Biomarkers_S1_Top/Human_Breast_Biomarkers_S1_Top_outs.zip
https://cf.10xgenomics.com/samples/xenium/4.0.0/Human_Breast_Biomarkers_S1_Top/Human_Breast_Biomarkers_S1_Top_he_image.ome.tif
https://cf.10xgenomics.com/samples/xenium/4.0.0/Human_Breast_Biomarkers_S1_Top/Human_Breast_Biomarkers_S1_Top_he_imagealignment.csv

sdata = spatialdata_io.xenium(
    data_dir,
    nucleus_boundaries=False,
    cells_labels=False,
    nucleus_labels=False,
    transcripts=False,
    aligned_images=True,
)
print(sdata)
print(sd.transformations.get_transformation(sdata.images["he_image"], get_all=True))

The dataset contains an additional H&E image with an Affine to the global coordinate system.

SpatialData object
├── Images
│     ├── 'he_image': DataTree[cyx] (3, 43995, 20515), (3, 21997, 10257), (3, 10998, 5128), (3, 5499, 2564), (3, 2749, 1282)
│     └── 'morphology_focus': DataTree[cyx] (4, 27420, 53994), (4, 13710, 26997), (4, 6855, 13498), (4, 3427, 6749), (4, 1713, 3374)
├── Shapes
│     └── 'cell_boundaries': GeoDataFrame shape: (209467, 1) (2D shapes)
└── Tables
      └── 'table': AnnData (209467, 280)
with coordinate systems:
    ▸ 'global', with elements:
        he_image (Images), morphology_focus (Images), cell_boundaries (Shapes)

{'global': Affine (x, y -> x, y)
     [ 9.88824121e-03  1.28956723e+00 -1.10808655e+03]
     [-1.28956723e+00  9.88824121e-03  2.67224245e+04]
     [0. 0. 1.]}

After transforming the image to the global coordinate system, the image contains only 3 large chunks.

image = sdata_full.images["he_image"]["scale0"].image
print(image)

image_trans = sd.transform(sdata.images["he_image"], to_coordinate_system="global")["scale0"].image
print(image_trans)

Before transform:

	Array	Chunk
Bytes	2.52 GiB	16.00 MiB
Shape	(3, 43995, 20515)	(1, 4096, 4096)
Dask graph	198 chunks in 5 graph layers
Data type	uint8 numpy.ndarray

After transform:

	Array	Chunk
Bytes	4.28 GiB	1.54 GiB
Shape	(3, 26890, 56937)	(3, 26890, 20515)
Dask graph	3 chunks in 9 graph layers
Data type	uint8 numpy.ndarray

for k, v in image.data[:, :1024, :1024].dask.layers.items():
    print(k, type(v), len(v))

for k, v in image_trans.data[:, :1024, :1024].dask.layers.items():
    print(k, type(v), len(v))

Before transform:

array-1f6906e306bdd2a15af122a29f1cc176 <class 'dask.highlevelgraph.MaterializedLayer'> 1
_map_read_frame-e4934b9d23dc66be228a602e126e7d51 <class 'dask.blockwise.Blockwise'> 1
getitem-2e03aac7dab88134e57e5b5893266c0c <class 'dask.highlevelgraph.MaterializedLayer'> 1
transpose-557e428061a6933b28b8db44ad0ed7d8 <class 'dask.blockwise.Blockwise'> 1
rechunk-merge-7640fe91c4d24108ec25cbaaedaeb750 <class 'dask.highlevelgraph.MaterializedLayer'> 396
getitem-9694c8fd0ff4f4b05c07751859e8b88a <class 'dask.highlevelgraph.MaterializedLayer'> 3

After transform:

affine_transform-e20839f41827529f6fcd15c7d7795e6b <class 'dask.highlevelgraph.MaterializedLayer'> 3
array-1f6906e306bdd2a15af122a29f1cc176 <class 'dask.highlevelgraph.MaterializedLayer'> 1
_map_read_frame-e4934b9d23dc66be228a602e126e7d51 <class 'dask.blockwise.Blockwise'> 1
getitem-2e03aac7dab88134e57e5b5893266c0c <class 'dask.highlevelgraph.MaterializedLayer'> 1
transpose-557e428061a6933b28b8db44ad0ed7d8 <class 'dask.blockwise.Blockwise'> 1
rechunk-merge-7640fe91c4d24108ec25cbaaedaeb750 <class 'dask.highlevelgraph.MaterializedLayer'> 396
getitem-bb69ef301749c714cbc84798ebeaeca8 <class 'dask.highlevelgraph.MaterializedLayer'> 72
getitem-f9415ec8602e80f5ee1d861ed3a40acc <class 'dask.highlevelgraph.MaterializedLayer'> 90
getitem-b95f72c34293e0492a54f7416d743a1a <class 'dask.highlevelgraph.MaterializedLayer'> 72
getitem-8f49fea7021e419e7c77ceea26ee6933 <class 'dask.highlevelgraph.MaterializedLayer'> 1

It seems that affine_transform will cause the reading of large blocks of the image.

Is any solutions of slicing the transformed image data that minimizes the reading?

Desktop:

• OS: Windows 10 22H2 19045.2604
• Version: spatialdata 0.7.3, spatialdata_io 0.7.0

[LucaMarconato]

Hey, I suggest to save the object to Zarr and re-read it. The object returned from spatialdata_io.xenium() is lazily constructed from OME-TIFF and it doesn't take advantage of the Zarr chunking.

[LucaMarconato]

Saving and reloading is the general recommendation for SpatialData objects coming from spatialdata-io. Nevertheless there is more to say here. In the case of Xenium >= 2.x.x and the morphology_focus image, we could actually take advantage of the pyramidal OME-TIFF format and construct the multiscale image directly from each scale and chunk (instead of lazily loading the large image and downscaling to a custom chunk size). This is currently not done as the spatialdata.models.Image2DModel.parse() method does not accept a DataTree object, but only a DataArray. If we extend this we can then also add this feature to spatialdata-io. Saving and reloading the object would still be preferable, but it would not be necessary for the case you reported.