Xenium–CODEX alignment (interactive Napari landmark registration)¶

Purpose¶

This notebook performs interactive spatial registration between:

• Xenium (RNA; using Xenium morphology image + nucleus boundaries), and
• CODEX (protein; multiplex IF image + per-cell table)

It creates a combined SpatialData object, lets you place corresponding landmarks in Napari, computes an affine transform, and writes aligned outputs to .zarr.
A second section optionally aligns an externally generated binary mask (e.g., from HALO) into the same aligned coordinate system.

Inputs¶

CODEX inputs (protein)¶

• A multi-channel CODEX image (e.g., *.tif)
• A per-cell measurement table (CSV) with bounding box columns XMin/XMax/YMin/YMax (used to compute centroids)
  • Optional columns used as metadata: cell/nucleus area fields, etc.
  • Marker intensity columns (features) are detected automatically (columns containing "Cell Intensity" unless list_of_marker_names is provided)

Xenium inputs (RNA / morphology)¶

• A Xenium SpatialData .zarr produced by Xenium processing / upstream pipeline
  • Expected elements used here:
    • images["morphology_focus"] (reference image)
    • shapes["nucleus_boundaries"] (reference shapes)
    • tables["table"] (cell/nucleus table)

Optional mask input¶

• A 2D label mask image (e.g., HALO-exported TIFF) where one or more integer pixel values represent foreground (e.g., collagen-positive).
  • The notebook converts this into a SpatialData labels/image element and registers it to the aligned coordinate frame.

Outputs¶

Paths are currently hard-coded to lab/HPC locations; for external users, set an OUTDIR and write there.

Typical outputs written by this notebook:

• slide_1_correction.zarr
  CODEX-only SpatialData created from the CODEX image + table.
• sdata_align_slide_1.zarr
  Combined Xenium+CODEX SpatialData after landmark-based registration.
• slide_1_left_masked.zarr
  Aligned SpatialData containing the registered binary mask element (name: "mask_fg").

Dependencies¶

Core Python packages:

• spatialdata, spatialdata-plot, anndata, scanpy
• napari-spatialdata (interactive UI; required for placing landmarks)
• tifffile, numpy, pandas, geopandas, shapely, xarray, scipy, matplotlib

Recommended run order¶

1. Build CODEX SpatialData from image + table (SptaialData_object)
2. Write CODEX SpatialData to .zarr (optional but recommended)
3. Create merged Xenium+CODEX SpatialData (Allign_interactive)
4. Launch Napari (Interactive(...)) and create landmark shapes:
   • "landmark_rna" on Xenium (reference)
   • "landmark_protein" on CODEX (moving)
5. Compute affine transform with align_elements_using_landmarks(...)
6. Write aligned .zarr
7. (Optional) Load a binary mask, add it as "mask_fg", and repeat landmark registration to align mask to the aligned data.

Notes / common pitfalls¶

• This workflow is manual/interactive: results depend on where you click landmarks in Napari.
• Coordinate systems must be consistent. The notebook sometimes renames the protein "global" coordinate system (see protein_coor).

In [ ]:

import tifffile
import os
from spatialdata.models import Image2DModel
import pandas as pd
import geopandas as gpd
from shapely.geometry import Point
import numpy as np
from spatialdata.models import ShapesModel
import anndata
import scanpy as sc
from spatialdata import SpatialData
from spatialdata import polygon_query
from spatialdata import bounding_box_query
from napari_spatialdata import Interactive
from spatialdata.models import TableModel
from spatialdata_plot.pl.utils import set_zero_in_cmap_to_transparent
from spatialdata import SpatialData
import spatialdata as sd
from scipy import ndimage as ndi
import seaborn as sns
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
import matplotlib as mpl
from spatialdata.transformations import (
    BaseTransformation,
    Identity,
    Sequence,
    align_elements_using_landmarks,
    get_transformation,
    set_transformation,
    Affine
)

import tifffile import os from spatialdata.models import Image2DModel import pandas as pd import geopandas as gpd from shapely.geometry import Point import numpy as np from spatialdata.models import ShapesModel import anndata import scanpy as sc from spatialdata import SpatialData from spatialdata import polygon_query from spatialdata import bounding_box_query from napari_spatialdata import Interactive from spatialdata.models import TableModel from spatialdata_plot.pl.utils import set_zero_in_cmap_to_transparent from spatialdata import SpatialData import spatialdata as sd from scipy import ndimage as ndi import seaborn as sns import matplotlib.pyplot as plt import matplotlib.colors as mcolors import matplotlib as mpl from spatialdata.transformations import ( BaseTransformation, Identity, Sequence, align_elements_using_landmarks, get_transformation, set_transformation, Affine )

/vf/users/kanferg/conda_v1/envs/spatialdata_env/lib/python3.12/site-packages/dask/dataframe/__init__.py:31: FutureWarning: The legacy Dask DataFrame implementation is deprecated and will be removed in a future version. Set the configuration option `dataframe.query-planning` to `True` or None to enable the new Dask Dataframe implementation and silence this warning.
  warnings.warn(
/vf/users/kanferg/conda_v1/envs/spatialdata_env/lib/python3.12/site-packages/xarray_schema/__init__.py:1: UserWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html. The pkg_resources package is slated for removal as early as 2025-11-30. Refrain from using this package or pin to Setuptools<81.
  from pkg_resources import DistributionNotFound, get_distribution

Step 1 — Create a CODEX SpatialData object¶

This section reads the CODEX image and the associated per-cell table, converts the table to an AnnData (marker intensities as features), and packages image + shapes + table into a SpatialData object.

In [ ]:

class CodexProcessing:
    def __init__(self,path_slide,slide_name,tableFile_name,keep_list = None,list_of_marker_names = None, radius = 2, table_export = None):
        self.path_slide = path_slide
        self.slide_name = slide_name
        self.tableFile_name = tableFile_name
        self.keep_list = keep_list #obs
        self.list_of_marker_names = list_of_marker_names #vars
        self.radius = radius
        self.table_export = table_export

class GenerComp(CodexProcessing):
    def __init__(self,*args, **kwargs):
        super().__init__(*args, **kwargs)
        
    def read_codex_image(self):
        img_array = tifffile.imread(os.path.join(self.path_slide,self.slide_name))
        print(f"image dimentions: {np.shape(img_array)}")
        return img_array
    
    def read_codex_table(self):
        if self.table_export:
            df =self.table_export
        else:
            df = pd.read_csv(os.path.join(self.path_slide,self.tableFile_name))
        df['centroid_x'] = (df['XMin'] + df['XMax']) / 2
        df['centroid_y'] = (df['YMin'] + df['YMax']) / 2
        return df
    
    def table_to_adata(self):
        self.df.rename(columns={
                'Cell Area (µm²)': 'Area',
                'Cytoplasm Area (µm²)': 'Cyto_Area',
                'Nucleus Area (µm²)': 'Nucleus_Area',
                'Nucleus Perimeter (µm)': 'Nucleus_Perimeter',
                'Classifier Label':'Classifier_Label',
                'DAPI Positive Classification':'DAPI'
            }, inplace=True)
        if self.keep_list:
            keep_list = self.keep_list
        else:
            keep_list =['Area','Cyto_Area','Nucleus_Area','Nucleus_Perimeter']
        if not self.list_of_marker_names:
            list_of_marker_names = [col for col in df.columns.to_numpy() if 'Cell Intensity' in col]
        else:
            list_of_marker_names = self.list_of_marker_names
        marker_cols = [col for col in self.df.columns if col in list_of_marker_names]
        X = self.df[marker_cols].to_numpy()         # shape (n_cells, n_markers)
        obs = self.df.loc[:,keep_list].copy()     # non-marker metadata as obs
        adata = anndata.AnnData(X=X, obs=obs)
        adata.var['markers'] = marker_cols     # naming the features (optional)
        return adata 

class SptaialData_object(GenerComp):
    def __init__(self,*args, **kwargs):
        super().__init__(*args, **kwargs)
        self.img_array = self.read_codex_image()
        self.df = self.read_codex_table()
        self.adata = self.table_to_adata()
        self.image_element = self.spatial_image()
        self.shapes_element = self.spatial_shape()
        self.table_element = self.spatial_table()
        self.sdata = self.generate()
    
    def spatial_image(self):
        image_element = Image2DModel.parse(self.img_array, scale_factors=(2,2,2))
        return image_element
        
    def spatial_shape(self):
        geometry = [Point(x, y) for x, y in zip(self.df['centroid_x'], self.df['centroid_y'])]
        gdf = gpd.GeoDataFrame(self.df[['centroid_x','centroid_y']], geometry=geometry)
        gdf['radius'] = self.radius
        shapes_element = ShapesModel.parse(gdf)
        return shapes_element
    
    def spatial_table(self):
        adata_sdata = TableModel.parse(self.adata)
        adata_sdata.uns["spatialdata_attrs"] = {"region": "cells",          # name of the shapes element
                                                "region_key": "region",     # obs column that indicates region name
                                                "instance_key": "cell_id"}
        adata_sdata.obs["region"] = "cells"
        adata_sdata.obs["cell_id"] = self.shapes_element.index.astype(str)
        adata_sdata.obs["region"] = adata_sdata.obs["region"].astype("category")
        adata_sdata.obs["cell_id"] = adata_sdata.obs["cell_id"].astype(str)
        return adata_sdata
    
    def generate(self):
        sdata = sd.SpatialData(
                images={"codex_image": self.image_element},
                shapes={"cells": self.shapes_element},
                tables={"cells_table": self.table_element})
        return sdata                

class CodexProcessing: def __init__(self,path_slide,slide_name,tableFile_name,keep_list = None,list_of_marker_names = None, radius = 2, table_export = None): self.path_slide = path_slide self.slide_name = slide_name self.tableFile_name = tableFile_name self.keep_list = keep_list #obs self.list_of_marker_names = list_of_marker_names #vars self.radius = radius self.table_export = table_export class GenerComp(CodexProcessing): def __init__(self,*args, **kwargs): super().__init__(*args, **kwargs) def read_codex_image(self): img_array = tifffile.imread(os.path.join(self.path_slide,self.slide_name)) print(f"image dimentions: {np.shape(img_array)}") return img_array def read_codex_table(self): if self.table_export: df =self.table_export else: df = pd.read_csv(os.path.join(self.path_slide,self.tableFile_name)) df['centroid_x'] = (df['XMin'] + df['XMax']) / 2 df['centroid_y'] = (df['YMin'] + df['YMax']) / 2 return df def table_to_adata(self): self.df.rename(columns={ 'Cell Area (µm²)': 'Area', 'Cytoplasm Area (µm²)': 'Cyto_Area', 'Nucleus Area (µm²)': 'Nucleus_Area', 'Nucleus Perimeter (µm)': 'Nucleus_Perimeter', 'Classifier Label':'Classifier_Label', 'DAPI Positive Classification':'DAPI' }, inplace=True) if self.keep_list: keep_list = self.keep_list else: keep_list =['Area','Cyto_Area','Nucleus_Area','Nucleus_Perimeter'] if not self.list_of_marker_names: list_of_marker_names = [col for col in df.columns.to_numpy() if 'Cell Intensity' in col] else: list_of_marker_names = self.list_of_marker_names marker_cols = [col for col in self.df.columns if col in list_of_marker_names] X = self.df[marker_cols].to_numpy() # shape (n_cells, n_markers) obs = self.df.loc[:,keep_list].copy() # non-marker metadata as obs adata = anndata.AnnData(X=X, obs=obs) adata.var['markers'] = marker_cols # naming the features (optional) return adata class SptaialData_object(GenerComp): def __init__(self,*args, **kwargs): super().__init__(*args, **kwargs) self.img_array = self.read_codex_image() self.df = self.read_codex_table() self.adata = self.table_to_adata() self.image_element = self.spatial_image() self.shapes_element = self.spatial_shape() self.table_element = self.spatial_table() self.sdata = self.generate() def spatial_image(self): image_element = Image2DModel.parse(self.img_array, scale_factors=(2,2,2)) return image_element def spatial_shape(self): geometry = [Point(x, y) for x, y in zip(self.df['centroid_x'], self.df['centroid_y'])] gdf = gpd.GeoDataFrame(self.df[['centroid_x','centroid_y']], geometry=geometry) gdf['radius'] = self.radius shapes_element = ShapesModel.parse(gdf) return shapes_element def spatial_table(self): adata_sdata = TableModel.parse(self.adata) adata_sdata.uns["spatialdata_attrs"] = {"region": "cells", # name of the shapes element "region_key": "region", # obs column that indicates region name "instance_key": "cell_id"} adata_sdata.obs["region"] = "cells" adata_sdata.obs["cell_id"] = self.shapes_element.index.astype(str) adata_sdata.obs["region"] = adata_sdata.obs["region"].astype("category") adata_sdata.obs["cell_id"] = adata_sdata.obs["cell_id"].astype(str) return adata_sdata def generate(self): sdata = sd.SpatialData( images={"codex_image": self.image_element}, shapes={"cells": self.shapes_element}, tables={"cells_table": self.table_element}) return sdata

In [ ]:

# Large input files for this notebook are hosted on Zenodo; download and unpack the Zenodo bundle locally and set the base data path below to your local folder.
DATA_ROOT = "data_external/nguyenlab-il15act-dataset-v1"
path = os.path.join(DATA_ROOT, "04_codex_slide1", "Xenium1_repeat_02232025_40844.csv")
df = pd.read_csv(
    path,
    engine="python",              # more forgiving than C engine
    on_bad_lines="warn",          # or "skip" to suppress warnings
    sep=",",
)
df = df.loc[df['DAPI Positive Classification']==1]

# Large input files for this notebook are hosted on Zenodo; download and unpack the Zenodo bundle locally and set the base data path below to your local folder. DATA_ROOT = "data_external/nguyenlab-il15act-dataset-v1" path = os.path.join(DATA_ROOT, "04_codex_slide1", "Xenium1_repeat_02232025_40844.csv") df = pd.read_csv( path, engine="python", # more forgiving than C engine on_bad_lines="warn", # or "skip" to suppress warnings sep=",", ) df = df.loc[df['DAPI Positive Classification']==1]

In [ ]:

path_slide1 = os.path.join(DATA_ROOT, "04_codex_slide1")
file_name_slide1 = 'xenium1_02232025.tif'
tableFile_name = 'Xenium1_repeat_02232025_40844_melanoma.er.qptiff_object_Data.csv'
keep_list = ['Area','Cyto_Area','Nucleus_Area','Nucleus_Perimeter','Classifier_Label','DAPI']
list_of_marker_names = ['Col1A1 Cell Intensity']
table_export = df

path_slide1 = os.path.join(DATA_ROOT, "04_codex_slide1") file_name_slide1 = 'xenium1_02232025.tif' tableFile_name = 'Xenium1_repeat_02232025_40844_melanoma.er.qptiff_object_Data.csv' keep_list = ['Area','Cyto_Area','Nucleus_Area','Nucleus_Perimeter','Classifier_Label','DAPI'] list_of_marker_names = ['Col1A1 Cell Intensity'] table_export = df

In [ ]:

sdata_slide1_protein = SptaialData_object(path_slide = path_slide1,slide_name = file_name_slide1,tableFile_name = tableFile_name,list_of_marker_names = list_of_marker_names,keep_list = keep_list,table_export = df) 
sdata_protein = sdata_slide1_protein.generate()
sdata_protein['cells_table'].obs['Classifier_Label'] = sdata_protein['cells_table'].obs['Classifier_Label'].astype(str)

sdata_slide1_protein = SptaialData_object(path_slide = path_slide1,slide_name = file_name_slide1,tableFile_name = tableFile_name,list_of_marker_names = list_of_marker_names,keep_list = keep_list,table_export = df) sdata_protein = sdata_slide1_protein.generate() sdata_protein['cells_table'].obs['Classifier_Label'] = sdata_protein['cells_table'].obs['Classifier_Label'].astype(str)

In [ ]:

path_sdata = os.path.join(DATA_ROOT, "05_spatialdata_zarr")
path_zarr = os.path.join(path_sdata,'slide_1_correction.zarr')
sdata_protein.write(path_zarr,overwrite=True)

path_sdata = os.path.join(DATA_ROOT, "05_spatialdata_zarr") path_zarr = os.path.join(path_sdata,'slide_1_correction.zarr') sdata_protein.write(path_zarr,overwrite=True)

Step 2 — Merge Xenium + CODEX and place landmarks in Napari¶

We build a combined SpatialData object with:

• Xenium morphology image + nucleus boundaries (reference)
• CODEX image + per-cell shapes/table (moving)

Then we open Napari (Interactive) to create corresponding landmark shapes:

• landmark_rna (reference coords on Xenium)
• landmark_protein (moving coords on CODEX)

These landmarks are used to compute an affine registration.

In [ ]:

class AllignmentInit:
    def __init__(self,path_sdata_protein,sdata_name_protein,path_sdata_rna,sdata_name_rna):
        self.path_sdata_protein = path_sdata_protein # end with /
        self.sdata_name_protein = sdata_name_protein
        self.path_sdata_rna = path_sdata_rna
        self.sdata_name_rna = sdata_name_rna
        self.sdata_protein,self.sdata_rna = self.read_data()    
    
    def read_data(self):
        return SpatialData.read(self.path_sdata_protein + self.sdata_name_protein),SpatialData.read(self.path_sdata_rna + self.sdata_name_rna)

        
class Allign_interactive(AllignmentInit):
    def __init__(self,protein_coor = None, *args, **kwargs):
        self.protein_coor = protein_coor #e.g. "codex_coor"
        super().__init__(*args, **kwargs)
        self.sdata_align = self.alignObject()
        
    def alignObject(self):
        if self.protein_coor:
            self.sdata_protein.rename_coordinate_systems({"global": self.protein_coor})
        sdata_align = sd.SpatialData(images={
                                            "xenium": self.sdata_rna.images["morphology_focus"],
                                            "codex": self.sdata_protein.images["codex_image"],
                                                },
                                                shapes = {"xenium_shape":self.sdata_rna.shapes["nucleus_boundaries"],
                                                          "codex_shape":self.sdata_protein.shapes["cells"],
                                                            },
                                                tables = {"xenium_table":self.sdata_rna.tables["table"],
                                                          "codex_table":self.sdata_protein.tables["cells_table"],})
        return sdata_align

class AllignmentInit: def __init__(self,path_sdata_protein,sdata_name_protein,path_sdata_rna,sdata_name_rna): self.path_sdata_protein = path_sdata_protein # end with / self.sdata_name_protein = sdata_name_protein self.path_sdata_rna = path_sdata_rna self.sdata_name_rna = sdata_name_rna self.sdata_protein,self.sdata_rna = self.read_data() def read_data(self): return SpatialData.read(self.path_sdata_protein + self.sdata_name_protein),SpatialData.read(self.path_sdata_rna + self.sdata_name_rna) class Allign_interactive(AllignmentInit): def __init__(self,protein_coor = None, *args, **kwargs): self.protein_coor = protein_coor #e.g. "codex_coor" super().__init__(*args, **kwargs) self.sdata_align = self.alignObject() def alignObject(self): if self.protein_coor: self.sdata_protein.rename_coordinate_systems({"global": self.protein_coor}) sdata_align = sd.SpatialData(images={ "xenium": self.sdata_rna.images["morphology_focus"], "codex": self.sdata_protein.images["codex_image"], }, shapes = {"xenium_shape":self.sdata_rna.shapes["nucleus_boundaries"], "codex_shape":self.sdata_protein.shapes["cells"], }, tables = {"xenium_table":self.sdata_rna.tables["table"], "codex_table":self.sdata_protein.tables["cells_table"],}) return sdata_align

In [ ]:

path_sdata_protein = os.path.join(DATA_ROOT, "05_spatialdata_zarr") + "/"
sdata_name_protein = 'slide_1_correction.zarr'
# path_sdata_rna = os.path.join(DATA_ROOT, "05_spatialdata_zarr") + "/"
# sdata_name_rna = 'output-XETG00202__0040847_Right__SCAF04317_Right_R1__20241121__163918.zarr'
path_sdata_rna = os.path.join(DATA_ROOT, "05_spatialdata_zarr") + "/"
sdata_name_rna = 'output-XETG00202__0040844_Left__SCAF04316_Left_R1__20241121__163918.zarr'
protein_coor = "codex_coor" # note sure if this is correct
#test_sdata = AllignmentInit(path_sdata_protein = path_sdata_protein, sdata_name_protein = sdata_name_protein, path_sdata_rna = path_sdata_rna, sdata_name_rna = sdata_name_rna)
sdata_align = Allign_interactive(path_sdata_protein = path_sdata_protein, sdata_name_protein = sdata_name_protein, path_sdata_rna = path_sdata_rna, sdata_name_rna = sdata_name_rna, protein_coor = protein_coor)

path_sdata_protein = os.path.join(DATA_ROOT, "05_spatialdata_zarr") + "/" sdata_name_protein = 'slide_1_correction.zarr' # path_sdata_rna = os.path.join(DATA_ROOT, "05_spatialdata_zarr") + "/" # sdata_name_rna = 'output-XETG00202__0040847_Right__SCAF04317_Right_R1__20241121__163918.zarr' path_sdata_rna = os.path.join(DATA_ROOT, "05_spatialdata_zarr") + "/" sdata_name_rna = 'output-XETG00202__0040844_Left__SCAF04316_Left_R1__20241121__163918.zarr' protein_coor = "codex_coor" # note sure if this is correct #test_sdata = AllignmentInit(path_sdata_protein = path_sdata_protein, sdata_name_protein = sdata_name_protein, path_sdata_rna = path_sdata_rna, sdata_name_rna = sdata_name_rna) sdata_align = Allign_interactive(path_sdata_protein = path_sdata_protein, sdata_name_protein = sdata_name_protein, path_sdata_rna = path_sdata_rna, sdata_name_rna = sdata_name_rna, protein_coor = protein_coor)

In [ ]:

sdata_align = sdata_align.sdata_align

sdata_align = sdata_align.sdata_align

In [ ]:

Interactive(sdata_align)

Interactive(sdata_align)

Step 3 — Compute affine transform from landmarks¶

Uses align_elements_using_landmarks(...) to compute an affine transformation that maps CODEX → Xenium coordinates and stores the result in a new coordinate system (e.g., "aligned").

In [ ]:

affine = align_elements_using_landmarks(
    references_coords=sdata_align["landmark_rna"],
    moving_coords=sdata_align["landmark_protein"],
    reference_element=sdata_align["xenium"],
    moving_element=sdata_align["codex"],
    reference_coordinate_system="global",
    moving_coordinate_system="codex_coor",
    new_coordinate_system="aligned")

affine = align_elements_using_landmarks( references_coords=sdata_align["landmark_rna"], moving_coords=sdata_align["landmark_protein"], reference_element=sdata_align["xenium"], moving_element=sdata_align["codex"], reference_coordinate_system="global", moving_coordinate_system="codex_coor", new_coordinate_system="aligned")

In [ ]:

path_sdata = '/data/HiTIF/data/spatialomics/melanoma/data/spatialdata'
path_zarr = os.path.join(path_sdata,'sdata_align_slide_1.zarr')
sdata_align.write(path_zarr,overwrite=True)

path_sdata = '/data/HiTIF/data/spatialomics/melanoma/data/spatialdata' path_zarr = os.path.join(path_sdata,'sdata_align_slide_1.zarr') sdata_align.write(path_zarr,overwrite=True)

Optional — Register an external binary mask (e.g., HALO collagen mask)¶

Loads a label mask TIFF, converts it into a SpatialData element (mask_fg), and registers it into the aligned Xenium/CODEX coordinate system using a second set of landmarks.

In [ ]:

def image_mask_to_labels_sdata(arr, foreground_val):
    # 1) make 0/1 binary mask
    bin_mask = np.isin(arr, [foreground_vals]).astype(np.uint8)
    
    # 2) labels are 2D (y, x) — NO channel axis
    da_lbl = xr.DataArray(bin_mask, dims=("y", "x"))
    
    # 3) parse into SpatialData labels element
    labels_elem = Labels2DModel.parse(da_lbl)
    sdata_labels = SpatialData(labels={"mask_label": labels_elem})
    
    # optional rename coord system
    sdata_labels.rename_coordinate_systems({"global": "mask_coo"})
    return sdata_labels,da_lbl

def arr_to_spatialdata_binary_labels(arr, foreground_vals, name="mask_binary"):
    """
    arr: 2D numpy array (label mask)
    foreground_vals: int or list of ints (values in arr to keep as foreground)
    name: name of the labels layer in SpatialData

    Returns: SpatialData with labels[name] being a uint8 {0,1} mask.
    """
    fg = np.atleast_1d(foreground_vals)

    # 1) Make binary mask: foreground -> 1, else -> 0
    bin_mask = np.isin(arr, np.atleast_1d(foreground_vals)).astype(np.uint8)

    # 2) Wrap as xarray (SpatialData models expect labeled dims)
    da = xr.DataArray(
                        bin_mask[None, :, :],                 # shape (1, y, x)
                        dims=("c", "y", "x"),
                        coords={"c": ["mask"]}                # channel name
                    )

    img = Image2DModel.parse(da)

    # 4) build SpatialData with an image element
    sdata = SpatialData(images={"mask_fg": img})
    sdata.rename_coordinate_systems({"global": "mask_coo"})
    return sdata

def load_tiff_level(src, series=0, level=-1):
    """
    Load a pyramid level from a TIFF using tifffile.
    level=0 is largest, level=-1 is smallest.
    """
    with tf.TiffFile(src) as tif:
        s = tif.series[series]
        if hasattr(s, "levels") and len(s.levels) > 0:
            lvl = s.levels[level]
            arr = lvl.asarray()  # loads into RAM; pick a small level unless you know it's manageable
        else:
            # no pyramid levels: fall back to the series itself
            arr = s.asarray()
    return arr

def top_labels_from_thumb(arr, thumb_step=32, topk=4):
    """
    Compute top-k most frequent label values on a downsampled thumbnail.
    Returns (thumb, values_sorted, counts_sorted).
    """
    thumb = arr[::thumb_step, ::thumb_step]
    vals, counts = np.unique(thumb, return_counts=True)
    order = np.argsort(counts)[::-1]
    vals = vals[order]
    counts = counts[order]
    return thumb, vals[:topk], counts[:topk]

def show_label_preview(thumb, title="Label preview (thumb)"):
    plt.figure(figsize=(7, 7))
    plt.imshow(thumb, interpolation="nearest")  # default colormap is fine for labels
    plt.title(title)
    plt.axis("off")
    plt.show()

def show_binary_masks(thumb, bg_val, label_vals):
    """
    Display background mask + masks for each label in label_vals (in the thumb space).
    """
    # Background
    plt.figure(figsize=(6, 6))
    plt.imshow((thumb == bg_val), cmap="gray", interpolation="nearest")
    plt.title(f"Background mask (value == {bg_val})")
    plt.axis("off")
    plt.show()

    # Foreground labels
    for v in label_vals:
        plt.figure(figsize=(6, 6))
        plt.imshow((thumb == v), cmap="gray", interpolation="nearest")
        plt.title(f"Mask (value == {v})")
        plt.axis("off")
        plt.show()

def show_bg_plus_label(thumb, bg_val, v):
    """
    Show a 2-class image: background vs one selected label.
    """
    two_class = np.where(thumb == v, 1, 0)  # label=v -> 1, everything else -> 0
    plt.figure(figsize=(6, 6))
    plt.imshow(two_class, cmap="gray", interpolation="nearest")
    plt.title(f"2-class: value {v} vs rest (bg assumed {bg_val})")
    plt.axis("off")
    plt.show()

def image_mask_to_labels_sdata(arr, foreground_val): # 1) make 0/1 binary mask bin_mask = np.isin(arr, [foreground_vals]).astype(np.uint8) # 2) labels are 2D (y, x) — NO channel axis da_lbl = xr.DataArray(bin_mask, dims=("y", "x")) # 3) parse into SpatialData labels element labels_elem = Labels2DModel.parse(da_lbl) sdata_labels = SpatialData(labels={"mask_label": labels_elem}) # optional rename coord system sdata_labels.rename_coordinate_systems({"global": "mask_coo"}) return sdata_labels,da_lbl def arr_to_spatialdata_binary_labels(arr, foreground_vals, name="mask_binary"): """ arr: 2D numpy array (label mask) foreground_vals: int or list of ints (values in arr to keep as foreground) name: name of the labels layer in SpatialData Returns: SpatialData with labels[name] being a uint8 {0,1} mask. """ fg = np.atleast_1d(foreground_vals) # 1) Make binary mask: foreground -> 1, else -> 0 bin_mask = np.isin(arr, np.atleast_1d(foreground_vals)).astype(np.uint8) # 2) Wrap as xarray (SpatialData models expect labeled dims) da = xr.DataArray( bin_mask[None, :, :], # shape (1, y, x) dims=("c", "y", "x"), coords={"c": ["mask"]} # channel name ) img = Image2DModel.parse(da) # 4) build SpatialData with an image element sdata = SpatialData(images={"mask_fg": img}) sdata.rename_coordinate_systems({"global": "mask_coo"}) return sdata def load_tiff_level(src, series=0, level=-1): """ Load a pyramid level from a TIFF using tifffile. level=0 is largest, level=-1 is smallest. """ with tf.TiffFile(src) as tif: s = tif.series[series] if hasattr(s, "levels") and len(s.levels) > 0: lvl = s.levels[level] arr = lvl.asarray() # loads into RAM; pick a small level unless you know it's manageable else: # no pyramid levels: fall back to the series itself arr = s.asarray() return arr def top_labels_from_thumb(arr, thumb_step=32, topk=4): """ Compute top-k most frequent label values on a downsampled thumbnail. Returns (thumb, values_sorted, counts_sorted). """ thumb = arr[::thumb_step, ::thumb_step] vals, counts = np.unique(thumb, return_counts=True) order = np.argsort(counts)[::-1] vals = vals[order] counts = counts[order] return thumb, vals[:topk], counts[:topk] def show_label_preview(thumb, title="Label preview (thumb)"): plt.figure(figsize=(7, 7)) plt.imshow(thumb, interpolation="nearest") # default colormap is fine for labels plt.title(title) plt.axis("off") plt.show() def show_binary_masks(thumb, bg_val, label_vals): """ Display background mask + masks for each label in label_vals (in the thumb space). """ # Background plt.figure(figsize=(6, 6)) plt.imshow((thumb == bg_val), cmap="gray", interpolation="nearest") plt.title(f"Background mask (value == {bg_val})") plt.axis("off") plt.show() # Foreground labels for v in label_vals: plt.figure(figsize=(6, 6)) plt.imshow((thumb == v), cmap="gray", interpolation="nearest") plt.title(f"Mask (value == {v})") plt.axis("off") plt.show() def show_bg_plus_label(thumb, bg_val, v): """ Show a 2-class image: background vs one selected label. """ two_class = np.where(thumb == v, 1, 0) # label=v -> 1, everything else -> 0 plt.figure(figsize=(6, 6)) plt.imshow(two_class, cmap="gray", interpolation="nearest") plt.title(f"2-class: value {v} vs rest (bg assumed {bg_val})") plt.axis("off") plt.show()

In [ ]:

path_codex = os.path.join(DATA_ROOT, "04_codex_slide1")
file_xen1 = f'{path_codex}/Tina_xenium1_large.tif'
arr_1 = load_tiff_level(file_xen1, series=0, level=-1)
print("Loaded level array:", arr_1.shape, arr_1.dtype)

path_codex = os.path.join(DATA_ROOT, "04_codex_slide1") file_xen1 = f'{path_codex}/Tina_xenium1_large.tif' arr_1 = load_tiff_level(file_xen1, series=0, level=-1) print("Loaded level array:", arr_1.shape, arr_1.dtype)

In [ ]:

sdata_mask = arr_to_spatialdata_binary_labels(arr = arr_1, foreground_vals = 1)

bin_mask = np.isin(arr_1, [1]).astype(np.uint8)
    
# 2) labels are 2D (y, x) — NO channel axis
da_lbl = xr.DataArray(bin_mask, dims=("y", "x"))

# 3) parse into SpatialData labels element
labels_elem = Labels2DModel.parse(da_lbl)

sdata_mask = arr_to_spatialdata_binary_labels(arr = arr_1, foreground_vals = 1) bin_mask = np.isin(arr_1, [1]).astype(np.uint8) # 2) labels are 2D (y, x) — NO channel axis da_lbl = xr.DataArray(bin_mask, dims=("y", "x")) # 3) parse into SpatialData labels element labels_elem = Labels2DModel.parse(da_lbl)

preprocessed sdata (previously aligned) file will be used for registering the mask too

In [ ]:

path_sdata = os.path.join(DATA_ROOT, "05_spatialdata_zarr")
path_zarr = os.path.join(path_sdata,'slide_1_left_maked.zarr')
sdata_1 = sd.read_zarr(path_zarr)

path_sdata = os.path.join(DATA_ROOT, "05_spatialdata_zarr") path_zarr = os.path.join(path_sdata,'slide_1_left_maked.zarr') sdata_1 = sd.read_zarr(path_zarr)

In [ ]:

sdata_1.labels["mask_fg"] = labels_elem

sdata_1.labels["mask_fg"] = labels_elem

In [ ]:

sdata_1

sdata_1

In [ ]:

t = get_transformation(sdata_1.images["codex_image_L"], to_coordinate_system="slide1_left")
set_transformation(sdata_1.labels["mask_fg"], transformation=t, to_coordinate_system="global")

t = get_transformation(sdata_1.images["codex_image_L"], to_coordinate_system="slide1_left") set_transformation(sdata_1.labels["mask_fg"], transformation=t, to_coordinate_system="global")

In [ ]:

sdata_1.write_element("mask_fg",overwrite=True)

sdata_1.write_element("mask_fg",overwrite=True)

Using spatialdata Interactive napari we can create new shapes, e.g. "codex_ref" and "mask_fg_ref" for register mask with xenium and codex.

In [ ]:

Interactive(sdata_1)

Interactive(sdata_1)

In [ ]:

affine = align_elements_using_landmarks(
    references_coords=sdata_1["codex_ref"],
    moving_coords=sdata_1["mask_fg_ref"],
    reference_element=sdata_1["codex_image_L"],
    moving_element=sdata_1["mask_fg"],
    reference_coordinate_system="slide1_left",
    moving_coordinate_system="global", #### not sure ??????
    new_coordinate_system="slide1_left")

affine = align_elements_using_landmarks( references_coords=sdata_1["codex_ref"], moving_coords=sdata_1["mask_fg_ref"], reference_element=sdata_1["codex_image_L"], moving_element=sdata_1["mask_fg"], reference_coordinate_system="slide1_left", moving_coordinate_system="global", #### not sure ?????? new_coordinate_system="slide1_left")

In [ ]:

path_sdata = os.path.join(DATA_ROOT, "05_spatialdata_zarr")
path_zarr = os.path.join(path_sdata,'slide_1_left_masked.zarr')
sdata_1.write(path_zarr,overwrite=True)

path_sdata = os.path.join(DATA_ROOT, "05_spatialdata_zarr") path_zarr = os.path.join(path_sdata,'slide_1_left_masked.zarr') sdata_1.write(path_zarr,overwrite=True)