register 3d images (X,Y,Z) and register from various objectives (x20 /x63) ?

[Nal44]

Hi,

I have a set of IF multiplex that I want to register but they are X.Y.Z (3D volume) ,
Using the code you sent me before (using the ch00 then match to all other channels with the same name) can I use the code or modify it to work with 3D images ?

Thanks a lot,
Antho

PS: also can I use valis to register images acquired from different objectives (x20 and x63) ? so scaling up/down before registering ?

Thanks

[cdgatenbee]

Hi @Nal44,
Could you tell me a bit more about how the data is structured? Are you wanting to align several 3D volumes, each in a different file (i.e. each file has several Z planes)? Is a set of images, where each is a single Z-plane? Is everything in 1 file, but needs to be aligned a bit better? Regardless, I think that there should be some way to align 3D IF images.

Best,

• Chandler

[Nal44]

Hi Chandler,

The first set is composed of several files 3D volumes (X,Y,Z), such as 5 individuals tif files (5 rounds per ROI) with 4 channels each (DAPI, beads, marker 1, marker 2).
So each file as: 4 channels Z,X,Z (12 planes approx.), so need to align all rounds to let's say to the DAPI channel (00), then apply the transform to other channels per round.

The cells shifted slightly between the rounds since it was acquired using multiple ROIS...

For the other set, I have a 2d image (low mag ) and a higher one (x63 oil) on similar ROI.

The idea is to register the 2d image (lower magnification ) onto the 3D image (higher magnification), then save as a OME.TIFF with 5 channels aligned (3d as 3 channels + 2d (dapi+marker) =5 channels total).

Hence, it needs to be scaled then registered , so I can use a segmentation algorithm the get the masks on the 3D image and the registered 2d image and cross all markers for quantification .

Let me know if this is clear or need more explanation,
Thanks a lot,
Antho

[cdgatenbee]

Hi @Nal44,
While it may be a little tricky, I think this should be possible. In the first case, if the individual planes within each volume do not need to be aligned, then it should be possible to align the DAPI channel in last plane of volume 1 to the DAPI channel in first plane volume 2, then last DAPI in volume 2 to first DAPI in volume 3, etc... If this sounds like what you had in mind, let me know, and I can try to put together some code. Even better would be if you could share a few images (at least 3, and they could be lower resolution), which would be very helpful in developing the code, since I don't have any images/volumes like this.

For the second set of images. I think what you would want to do is make the 3D image your reference image (Valis(reference_img_f="path/to/3D_volume")), so that the 2D image will be aligned to it. Since the 2D image is the one being warped, it should get scaled up to the same size as the 3D one. During the registration, both images will be rescaled to have similar sizes, so any differences in resolution in the full scale images shouldn't be an issue. So, for this dataset, maybe try that approach and let me know how it goes.

Best,
-Chandler

[Nal44]

Hi Chandler,

Thanks a lot,
I don't think it will since the first and last planes are out of focus in generally when we do Z stack acquisitions, hence the middle plane might work better (focus plane).

The way I did was to take the first Serie, segment the nucleus in 3d, then use this mask to manually align the other series to the first (translations/rotations) and applied the the same operations to the corresponding channels for each round (not very effective).
The result is ok, but not perfect as expected.

These images are from collaborators, hence I will check first with them if I can share them (these are cells it should be ok I believe).

For the second point I will try next week , I got delayed with too much work this week,
Thanks a lot ;D.
antho

[Nal44]

Hi Chandler,

Files are uploading now.
there are 5 files from one F02.
The task is to align r2,r3r,4,r5 to R1 (all channels ) , files are X,Y,Z 4 channels.

Please let me know once received.

thanks a lot for the help /support !
Antho

[cdgatenbee]

Hi @Nal44,
Just wanted to let you know that I was able to download the images (thanks for sharing!). Also, thank you for letting me know that using the middle plane may be best, which should actually be easier than trying to use the first or last plane. There are a few small things I need to take care of, but I'll get to this asap, and will let you know how it goes.

Best,
-Chandler

[Nal44]

Hi Chandler,

No probs, thank you for helping out ; D

Please let me know,
antho

[cdgatenbee]

Hi @Nal44,
Thanks again for sharing the images, I was able to use them to put together the following script that will align and merge your 3D volumes. The "trick" was to make a new SlideReader class that by default reads the center plane of the volume, which can be done by subclassing the BioFormatsSlideReader and tweaking the slide2vips and slide2image methods:

from valis import registration, warp_tools, slide_io, valtils

class Read3D(slide_io.BioFormatsSlideReader):
    """Reads middle plane of 3D volume

    """

    def __init__(self, src_f, series=None, *args, **kwargs):
        """
        Parameters
        -----------
        src_f : str
            Path to slide

        series : int
            The series to be read. If `series` is None, the the `series`
            will be set to the series associated with the largest image.

        """
        slide_io.init_jvm()
        super().__init__(src_f=src_f, *args, **kwargs)


    def slide2vips(self, level, series=None, xywh=None, tile_wh=None, z=None, t=0, *args, **kwargs):
        """Convert slide to pyvips.Image
        If `z` is None, the center plane of the volume will be returned
        """

        if z is None:
            z = self.metadata.n_z // 2 # Get middle plane

        vips_slide = super().slide2vips(level, series=series, xywh=xywh, tile_wh=tile_wh, z=z, t=0, *args, **kwargs)

        return vips_slide

    def slide2image(self, level, series=None, xywh=None, tile_wh=None, z=None, t=0, *args, **kwargs):
        """Convert slide to numpy array
        If `z` is None, the center plane of the volume will be returned
        """

        if z is None:
            z = self.metadata.n_z // 2 # Get middle plane

        img = super().slide2image(level, series=series, xywh=xywh, tile_wh=tile_wh, z=z, t=0, *args, **kwargs)

        return img

One can now align each volume to R1 by passing the above Read3D class to the reader_cls argument of Valis.register:

img_list = list(pathlib.Path(src_dir).rglob("*.tiff"))

ref_img = [str(x) for x in img_list if str(x).endswith("R1.tiff")][0]
registrar = registration.Valis(src_dir, dst_dir, reference_img_f=ref_img, align_to_reference=True)
rigid_registrar, non_rigid_registrar, error_df = registrar.register(reader_cls=Read3D)

Running the above produce the following overlap image, which I think looks pretty good:

The Valis.warp_and_merge_slides method wasn't really designed to merge 3D images, but you can use the following method to warp and merge each plane, which can then be stacked to build the volume:

def warp_and_merge_volumes(registrar, level=0, non_rigid=True, crop=True, channel_names=None, interp_method="bicubic"):

    ref_slide = registrar.get_ref_slide()
    ref_reader = ref_slide.reader

    # Aligned to reference, but did not assume images were ordered.
    # Want to merge channels in the same order as the rounds
    src_f_list = registrar.get_sorted_img_f_list()
    valtils.sort_nicely(src_f_list)

    if level != 0:
        if not np.issubdtype(type(level), np.integer):
            msg = "Need slide level to be an integer indicating pyramid level"
            valtils.print_warning(msg)
        aligned_slide_shape = ref_slide.val_obj.get_aligned_slide_shape(level)
    else:
        aligned_slide_shape = ref_slide.aligned_slide_shape_rc

    do_crop = False
    if isinstance(crop, bool) or isinstance(crop, str):
        do_crop = True
        crop_method = ref_slide.get_crop_method(crop)
        if crop_method is not False:
            if crop_method == registration.CROP_REF:
                scaled_aligned_shape_rc = ref_slide.slide_dimensions_wh[level][::-1]

            elif crop_method == registration.CROP_OVERLAP:
                scaled_aligned_shape_rc = aligned_slide_shape

    # Merge each plane, and then stack to create volume
    n_planes = ref_reader.metadata.n_z
    all_merged_planes = [None] * n_planes
    for i in range(n_planes):
        merged_plane = None
        for src_f in src_f_list:
            slide_obj = registrar.get_slide(src_f)

            # Get plane and warp using registration parameters
            plane = slide_obj.reader.slide2vips(level, z=i)
            if non_rigid:
                bk_dxdy = slide_obj.bk_dxdy
            else:
                bk_dxdy = None

            if do_crop:
                crop_bbox_xywh, _ = slide_obj.get_crop_xywh(crop=crop_method, out_shape_rc=scaled_aligned_shape_rc)
            else:
                crop_bbox_xywh = None

            warped_plane = warp_tools.warp_img(img=plane,
                                               M=slide_obj.M,
                                               bk_dxdy=bk_dxdy,
                                               transformation_src_shape_rc=slide_obj.processed_img_shape_rc,
                                               transformation_dst_shape_rc=slide_obj.reg_img_shape_rc,
                                               out_shape_rc=aligned_slide_shape,
                                               bbox_xywh=crop_bbox_xywh,
                                               interp_method=interp_method)
            if merged_plane is None:
                merged_plane = warped_plane
            else:
                merged_plane = merged_plane.bandjoin(warped_plane)

        all_merged_planes[i] = merged_plane

    merged_volume = all_merged_planes[0].bandjoin(all_merged_planes[1:])

    # Create ome-xml for merged volume
    size_x = merged_plane.width
    size_y = merged_plane.height
    size_z = n_planes
    size_c = merged_plane.bands
    size_t = 1
    bf_dtype = slide_io.vips2bf_dtype(merged_volume.format)
    shape_xyzct = [size_x, size_y, size_z, size_c, size_t]
    pixel_physical_size_xyu = ref_slide.reader.metadata.pixel_physical_size_xyu
    channel_names = [registrar.get_slide(src_f).reader.metadata.channel_names for src_f in src_f_list]
    channel_names = [j for i in channel_names for j in i]


    ome_xml_obj = slide_io.create_ome_xml(shape_xyzct=shape_xyzct,
                                      bf_dtype=bf_dtype,
                                      is_rgb=False,
                                      pixel_physical_size_xyu=pixel_physical_size_xyu,
                                      channel_names=channel_names
                                      )

    ome_xml = ome_xml_obj.to_xml()

    return merged_volume, ome_xml

Running the above creates a volume with 16 planes, and 20 channels per plane. It also produces the ome-xml metadata which you can use to save the image using slide_io.save_ome_tiff:

merged_volume, ome_xml = warp_and_merge_volumes(registrar)

out_f = os.path.join(dst_dir, f"{registrar.name}_merged.ome.tiff")
slide_io.save_ome_tiff(merged_volume, out_f, ome_xml=ome_xml)

It did take quite a bit of time to save this image (1.5h), but there may be ways to speed this up (different compression method, lower Q, etc...?). Anyhow, once saved, the registered volume looks like this in QuPath:

Please try out the above script and let me know how it goes.

Best,
-Chandler

[Nal44]

Hi Chandler,

Thanks a lot !!
This looks great :) , and what I need.
I will give it a go this week, when I get some free time (busy week) .

Great support and very helpful :D
Thanks
antho