How to do 2D image registration?

[risan-raja]

Hi, I am trying to register two histology sections and I am stumped on to figure out, how to do 2D image registration. Could you help me out here? @rohitrango

[rohitrango]

can you provide me with some more details, i.e. what are you trying to register, what data format are you using, have you tried adopting the basic tutorial to your data, what kind of registration (affine, rigid, deformable) do you want?

if you provide a minimal reproducible example, we can go from there

[risan-raja]

Hi @rohitrango , Thanks for getting back to me. Yes I have followed your basic tutorial. I am trying to register two histological WSIs. I am trying to register them in a downsampled manner as they are more than 100k x 100k in size. Will this help? I am primarily looking forword to use non-rigid deformation ie I want to extract the DVF from the syn registration method. I have initially aligned them using SIFT, the size of the downsampled images are in the size 8192px . I have divided the image into 4x4 grid and now i am trying to register against the patch corresponding to the fixed/moving image.

[rohitrango]

have you tried registering them with ANTs first? Or do you have any pipeline for any baseline? I'm unfamiliar with literature on deformable registration for WSI, can you share some papers or public datasets?

If you followed the basic tutorial, then FireANTs should work out of the box. The image loader and registration algorithms will detect it is a 2D image and perform intensity-based registration in 2D. Let me know if you face any error.

[risan-raja]

Hi, Yes I have experimented with ANTsPy first, however it is quite slow and doesn't make use of GPU. Whenever I try to load a single image, the inverse function throws an error in your code. Could I mail you ?

[risan-raja]

This is an example of what I am trying to align

patches_moving: list[NDArray]
patches_fixed: list[NDArray]
sitk_patches_moving = [sitk.GetImageFromArray(p) for p in patches_moving]
sitk_patches_fixed = [sitk.GetImageFromArray(p) for p in patches_fixed]
fire_ants_moving = [Image(x) for x in sitk_patches_moving]
fire_ants_fixed = [Image(x) for x in sitk_patches_fixed]
b1 = BatchedImages([fire_ants_moving[9]])
b2 = BatchedImages([fire_ants_fixed[9]])
# specify some values
scales = [4, 2, 1]  # scales at which to perform registration
iterations = [100, 50, 20]
optim = "Adam"
lr = 3e-3
loss_type = "cc"
loss_device = "cuda:1"
# create affine registration object
affine = AffineRegistration(
    scales,
    iterations,
    b1,
    b2,
    optimizer=optim,
    optimizer_lr=lr,
    cc_kernel_size=5,
    loss_type=loss_type,
    loss_device=loss_device,
)

Now when I try to register them using the function below

transformed_images = affine.optimize(save_transformed=True)

I get the following error.

{
    "name": "NotImplementedError",
    "message": "Got 4D input, but trilinear mode needs 5D input",
    "stack": "---------------------------------------------------------------------------
NotImplementedError                       Traceback (most recent call last)
Cell In[16], line 2
      1 start = time()
----> 2 transformed_images = affine.optimize(save_transformed=True)
      3 torch.cuda.empty_cache()
      4 end = time()

File /apps/FireANTS/fireants/registration/affine.py:82, in AffineRegistration.optimize(self, save_transformed)
     80     sigmas = 0.5 * torch.tensor([sz/szdown for sz, szdown in zip(fixed_size, size_down)], device=fixed_arrays.device)
     81     gaussians = [gaussian_1d(s, truncated=2) for s in sigmas]
---> 82     fixed_image_down = downsample(fixed_arrays, size=size_down, mode=self.fixed_images.interpolate_mode, gaussians=gaussians)
     83     moving_image_blur = separable_filtering(moving_arrays, gaussians)
     84 else:

File /apps/FireANTS/fireants/utils/imageutils.py:109, in downsample(image, size, mode, sigma, gaussians)
    107 # otherwise gaussians is given, just downsample
    108 image_smooth = separable_filtering(image, gaussians)
--> 109 image_down = F.interpolate(image_smooth, size=size, mode=mode, align_corners=True)
    110 return image_down

File /usr/local/lib/python3.10/dist-packages/torch/nn/functional.py:4082, in interpolate(input, size, scale_factor, mode, align_corners, recompute_scale_factor, antialias)
   4080     raise NotImplementedError(\"Got 4D input, but linear mode needs 3D input\")
   4081 if input.dim() == 4 and mode == \"trilinear\":
-> 4082     raise NotImplementedError(\"Got 4D input, but trilinear mode needs 5D input\")
   4083 if input.dim() == 5 and mode == \"linear\":
   4084     raise NotImplementedError(\"Got 5D input, but linear mode needs 3D input\")

NotImplementedError: Got 4D input, but trilinear mode needs 5D input"
}

I believe it's because torch.interpolate doesn't support 2D arrays and I have to give it 3D arrays. But when I try to expand the dimension by adding a dummy axis, the inverse of the image doesn't compute and I get an error. Or should I make a 2d Patch $(H, W)$ to $(1,1, H, W)$

Now when I try this:

patches_moving: list[NDArray]
patches_fixed: list[NDArray]
sitk_patches_moving = [sitk.GetImageFromArray(np.expand_dims(p, 0)) for p in patches_moving]
fire_ants_moving = [Image(x) for x in sitk_patches_moving]

The following error appears:

{
    "name": "_LinAlgError",
    "message": "linalg.inv: The diagonal element 3 is zero, the inversion could not be completed because the input matrix is singular.",
    "stack": "---------------------------------------------------------------------------
_LinAlgError                              Traceback (most recent call last)
Cell In[22], line 3
      1 sitk_patches_moving = [sitk.GetImageFromArray(np.expand_dims(p, 0)) for p in patches_moving]
      2 sitk_patches_fixed = [sitk.GetImageFromArray(np.expand_dims(p, 0)) for p in patches_fixed]
----> 3 fire_ants_moving = [Image(x) for x in sitk_patches_moving]
      4 fire_ants_fixed = [Image(x) for x in sitk_patches_fixed]

Cell In[22], line 3, in <listcomp>(.0)
      1 sitk_patches_moving = [sitk.GetImageFromArray(np.expand_dims(p, 0)) for p in patches_moving]
      2 sitk_patches_fixed = [sitk.GetImageFromArray(np.expand_dims(p, 0)) for p in patches_fixed]
----> 3 fire_ants_moving = [Image(x) for x in sitk_patches_moving]
      4 fire_ants_fixed = [Image(x) for x in sitk_patches_fixed]

File /apps/FireANTS/fireants/io/image.py:59, in Image.__init__(self, itk_image, device, is_segmentation, max_seg_label, background_seg_label, seg_preprocessor, spacing, direction, origin, center)
     57 # save the mapping from physical to torch and vice versa
     58 self.torch2phy = torch.from_numpy(np.matmul(px2phy, torch2px)).to(device).float().unsqueeze(0)
---> 59 self.phy2torch = torch.inverse(self.torch2phy[0]).float().unsqueeze(0)
     60 # also save intermediates just in case (as numpy arrays)
     61 self._torch2px = torch2px

_LinAlgError: linalg.inv: The diagonal element 3 is zero, the inversion could not be completed because the input matrix is singular."
}

[rohitrango]

This is a bug in the code, where interpolate_mode is always set to trilinear.

In fireants/io/image.py, Line 90 can you please change

        self.interpolate_mode = 'bilinear' if self.images[0] == 2 else 'trilinear'

to

        self.interpolate_mode = 'bilinear' if len(self.images[0].shape) == 4 else 'trilinear'

this should work, if it does you can submit a PR for the bug fix. Meanwhile, I should set up some unittests, I never really tested the code for 2D images.