multiview-stitcher / multiview-stitcher

A toolbox for registering / fusing / stitching large multi-view / multi-positioning image datasets in 2-3D.
https://multiview-stitcher.github.io/multiview-stitcher/
BSD 3-Clause "New" or "Revised" License
47 stars 6 forks source link

Best way to get global coordinates of fused image? #21

Closed dpshepherd closed 5 months ago

dpshepherd commented 5 months ago

Hi @m-albert,

What is the best way to get the global coordinates of the fused image? The use case is that we segment the fused image using Cellpose. Using the resulting label image from Cellpose, which doesn't know about global coordinates, we then need to iterate over the segmented cells, grab the coordinates of that label, and check if any of our point localizations (processed separately) fall within the label.

Normally I would do this by shifting everything using a reference tile and then calculating from there, but I realized I'm not entirely sure if there is a tile held in a constant position as reference during the global optimization steps.

Thanks!!

m-albert commented 5 months ago

Normally I would do this by shifting everything using a reference tile and then calculating from there, but I realized I'm not entirely sure if there is a tile held in a constant position as reference during the global optimization steps.

@dpshepherd The easiest will be to use the transformations stored in a given msims to transform its pixel coordinates into physical coordinates in the registered space (transform_key='registered'). You can do this as follows independently of whether your pixel locations correspond to input tiles or a fused image:

# extract transformations from a msim
affine = msi_utils.get_transform_from_msim(msim, transform_key='registered').data.squeeze()
origin = si_utils.get_origin_from_sim(msi_utils.get_sim_from_msim(msim), asarray=True)
spacing = si_utils.get_spacing_from_sim(msi_utils.get_sim_from_msim(msim), asarray=True)

# example pixel coordinate in msim
pt_px = np.array([100, 150])

# what's the physical coordinate of this pixel location corresponding to transform_key='registered'?

# variant 1: use skimage to do the transformation
import skimage

T_px_to_phys_image = skimage.transform.AffineTransform(
    scale=spacing, translation=origin)
T_phys_image_to_phys_registered = skimage.transform.AffineTransform(
    matrix=affine)
T_px_to_phys_registered = T_phys_image_to_phys_registered + T_px_to_phys_image

pt_phys_registered_skimage = T_px_to_phys_registered(pt_px)

# variant 2: pure numpy version of the above code snippet:
pt_phys_image = pt_px * spacing + origin
pt_phys_registered_np = (np.array(affine) @ np.array(list(pt_phys_image) + [1]))[:-1]

print(pt_phys_registered_skimage[0], pt_phys_registered_np)

Here's a example script recapitulating what I understood from the description of your workflow.

Small registration / fusion example:

import numpy as np
import skimage
from scipy import ndimage
import tifffile
import dask.array as da

from multiview_stitcher import spatial_image_utils as si_utils
from multiview_stitcher import msi_utils, registration, fusion, transformation, param_utils

im = skimage.data.human_mitosis()

tile1 = im[:300]
tile2 = im[260:]

isotropic_spacing = 10

msims = [
    msi_utils.get_msim_from_sim(
        si_utils.get_sim_from_array(
        tile, dims=['y', 'x'],
        translation={'y': y_offset, 'x': 0},
        scale={'y': isotropic_spacing, 'x': isotropic_spacing},
        transform_key='metadata'
        ), scale_factors=[]
    ) for tile, y_offset in zip([tile1, tile2], [0, 250 * isotropic_spacing])
]

params = registration.register(
    msims,
    reg_channel_index=0,
    transform_key='metadata', new_transform_key='registered',
    plot_summary=True,
    )

fused = fusion.fuse(
    [msi_utils.get_sim_from_msim(msim) for msim in msims],
    transform_key='registered').sel(t=0, c=0)

from matplotlib import pyplot as plt
plt.figure()
plt.imshow(fused)

image

image

Extract some pixel coordinates independently from the tiles and the fused image.

point_locs_pixel_space = [
    np.array([p.centroid for p in skimage.measure.regionprops(ndimage.label(im > 50)[0])])
    for im in [tile1, tile2, fused]
]

Transform the pixel locations into physical coordinates in "registered" space:

point_locs_pixel_space = [
    np.array([p.centroid for p in skimage.measure.regionprops(ndimage.label(im > 50)[0])])
    for im in [tile1, tile2, fused]
]

point_locs_phys_registered_space = []
for itile, (pts, msim) in enumerate(
    zip(
        point_locs_pixel_space,
        msims + [msi_utils.get_msim_from_sim(fused, scale_factors=[])]
        )
    ):

    # get origin, spacing and affine into registered coordinate system
    origin = si_utils.get_origin_from_sim(msi_utils.get_sim_from_msim(msim), asarray=True)
    spacing = si_utils.get_spacing_from_sim(msi_utils.get_sim_from_msim(msim), asarray=True)
    affine = msi_utils.get_transform_from_msim(msim, transform_key='registered').data.squeeze()

    T_px_to_phys_image = skimage.transform.AffineTransform(
        scale=spacing, translation=origin)
    T_phys_image_to_phys_registered = skimage.transform.AffineTransform(
        matrix=affine)

    point_locs_phys_registered_space.append(
        T_phys_image_to_phys_registered(T_px_to_phys_image(pts))
    )

Plot all points in the registered coordinate system:

from matplotlib import pyplot as plt

plt.figure()
for ipt, pts in enumerate(point_locs_phys_registered_space[:]):
    plt.scatter(
        pts[:, 1], pts[:, 0],
        label=['tile1', 'tile2', 'fused'][ipt],
        marker=['x', 'o', '+'][ipt])
plt.legend()
plt.title('Detections from tiles and fused images in registered physical space')

image

dpshepherd commented 5 months ago

Thank you!! We'll work to implement and circle back around with confirmation or more questions.

dpshepherd commented 5 months ago

This all worked perfectly, using the origin from the msim was the key (was not properly accounting for that before). The context here is placing 3D localizations (for MERFISH) from many individual tiles into their proper global coordinates from the full-scale data, while also placing cell segmentations from a fused, downsampled image onto the same global coordinate system.

Here all of the localizations shifted into place using the transforms for 40ish tiles spanning a 5 mm x 5 mm x 50 um acquisition.

Thanks!

image

m-albert commented 5 months ago

Great to hear it worked!

Thanks a lot for the feedback, explaining your use case and even posting an image of the result. Very appreciated!