Histology image registration

[rockdeme]

Hey!

I'm interested in registering some histology images and I was wondering if Elastix could be used for this purpose. I tried to register some demo images, which were cropped from the same tissue section and I introduced a 30° rotation to one of them. As the images are almost identical, in theory it should not be a difficult task but I can't get any meaningful results.

Demo images can be downloaded from here: https://drive.google.com/file/d/1iSplSMc6WaANrDq7UilynyVLgrvr4Jpt/view?usp=sharing

import itk
import numpy as np
from PIL import Image
import matplotlib.pyplot as plt

# open and convert images to grayscale
moving_image = Image.open("demo_1.png").convert('L')
moving_image = np.array(moving_image).astype(np.float32)
moving_image /= 255.0
moving_image = itk.image_view_from_array(moving_image)

fixed_image = Image.open("demo_2.png").convert('L')
fixed_image = np.array(fixed_image, np.float32)
fixed_image /= 255.0
fixed_image = itk.image_view_from_array(fixed_image)

# initialize Elastix with affine transform
parameter_object = itk.ParameterObject.New()
default_affine_parameter_map = parameter_object.GetDefaultParameterMap('affine')
default_affine_parameter_map['FinalBSplineInterpolationOrder'] = ['0']
parameter_object.AddParameterMap(default_affine_parameter_map)

result_image_affine, result_transform_parameters = itk.elastix_registration_method(fixed_image,
                                                                                   moving_image,
                                                                                   parameter_object=parameter_object,
                                                                                   log_to_console=True)

# plot the results
im = np.zeros((fixed_image.shape[0], fixed_image.shape[1], 3))
im[:, :, 0] = fixed_image
im[:, :, 1] = result_image_affine

fig, axs = plt.subplots(1, 4, sharey=True, figsize=[15, 5])
axs[0].imshow(result_image_affine)
axs[0].set_title('Result', fontsize=30)
axs[1].imshow(fixed_image)
axs[1].set_title('Fixed', fontsize=30)
axs[2].imshow(moving_image)
axs[2].set_title('Moving', fontsize=30)
axs[3].imshow(im)
axs[3].set_title('Overlap', fontsize=30)
plt.show()

I don't see any improvements in the error metric when looking at the output.

These are the results :

[thewtex]

Hi,

First, please optimize the rigid and affine transform before trying a bspline, which will be unnecessary in this case regardless.

[rockdeme]

Thanks for the quick reply @thewtex!

I've been trying to optimize a simple Euler transform for starters, but there is just no improvement no matter what parameter / optimizer / metric I choose. Honestly I'm a bit lost on how to proceed from here. Also tried chaining with a subsequent affine transform and got the same results.

parameter_object = itk.ParameterObject.New()
parameter_map_rigid = parameter_object.GetDefaultParameterMap('rigid')

parameter_map_rigid['FixedInternalImagePixelType'] = ['float']
parameter_map_rigid['FixedImageDimension'] = ['2']
parameter_map_rigid['MovingInternalImagePixelType'] = ['float']
parameter_map_rigid['MovingImageDimension'] = ['2']

parameter_map_rigid['Registration'] = ['MultiResolutionRegistration']
parameter_map_rigid['Interpolator'] = ['BSplineInterpolator']
parameter_map_rigid['ResampleInterpolator'] = ['FinalBSplineInterpolator']
parameter_map_rigid['Resampler'] = ['DefaultResampler']

parameter_map_rigid['FixedImagePyramid'] = ['FixedRecursiveImagePyramid']
parameter_map_rigid['MovingImagePyramid'] = ['MovingRecursiveImagePyramid']

parameter_map_rigid['Optimizer'] = ['AdaptiveStochasticGradientDescent']
parameter_map_rigid['Transform'] = ['EulerTransform']
parameter_map_rigid['Metric'] = ['AdvancedMattesMutualInformation']

parameter_map_rigid['AutomaticScalesEstimation'] = ['true']
parameter_map_rigid['NumberOfResolutions'] = ['2']

parameter_map_rigid['MaximumNumberOfIterations'] = ['1000']
parameter_map_rigid['NumberOfSpatialSamples'] = ['2048']
parameter_map_rigid['NewSamplesEveryIteration'] = ['true']
parameter_map_rigid['ImageSampler'] = ['Random']

parameter_map_rigid['BSplineInterpolationOrder'] = ['1']
parameter_map_rigid['FinalBSplineInterpolationOrder'] = ['3']

parameter_object.AddParameterMap(parameter_map_rigid)

result_image_affine, result_transform_parameters = itk.elastix_registration_method(fixed_image,
                                                                                   moving_image,
                                                                                   parameter_object=parameter_object,
                                                                                   log_to_console=True)

Log snippet

1:ItNr  2:Metric    3a:Time 3b:StepSize 4:||Gradient||  Time[ms]
0   -0.076791   0.000000    7.235377    0.021713    267.9
1   -0.071761   0.000000    7.235377    0.019584    1.1
2   -0.075115   0.000000    7.235377    0.035574    1.1
3   -0.073115   0.000000    7.235377    0.006367    1.0
4   -0.080415   0.936323    6.926545    0.031977    0.9
5   -0.073495   0.194758    7.168891    0.011891    1.1
6   -0.081848   1.168299    6.854063    0.032334    1.0
7   -0.071832   2.167861    6.558349    0.012248    0.9
8   -0.087272   1.405984    6.781354    0.032551    0.9
9   -0.081111   0.682835    7.007521    0.020713    1.0
10  -0.073901   0.000000    7.235377    0.035399    0.9
11  -0.080938   0.000000    7.235377    0.024698    1.0
12  -0.082372   0.000000    7.235377    0.017510    1.0
13  -0.079852   0.998593    6.906938    0.009833    1.3
...
990 -0.076338   84.885894   1.434968    0.014674    0.9
991 -0.075349   84.610046   1.438717    0.022607    0.9
992 -0.072099   84.195247   1.444390    0.018278    0.9
993 -0.074220   83.624382   1.452271    0.014025    0.9
994 -0.079642   83.141341   1.459007    0.009778    0.9
995 -0.073838   82.469599   1.468479    0.033366    0.9
996 -0.073229   83.224529   1.457842    0.011144    0.9
997 -0.081512   84.221872   1.444024    0.028603    0.9
998 -0.088392   83.461767   1.454531    0.008984    0.9
999 -0.076497   84.422828   1.441271    0.031268    0.9
Time spent in resolution 1 (ITK initialization and iterating): 1.25
Stopping condition: Maximum number of iterations has been reached.

Visual output is the same as the one in the original post.

[thewtex]

Hi, yes, first it is essential to get a Euler transform optimized.

parameter_map_rigid['NumberOfResolutions'] = ['2']

Given the high frequency content of the images, this should likely be increased.

You could also try a different similarity metric.

[dzenanz]

Perhaps fiddle with parameter scales to get the optimization more rotation-prone?

[mstaring]

All of those, and also increase (MaximumStepLength 1.0) the default is 1, try 2 or 3.

I often look at the column stepsize and the column ||gradient|| and aim for the product to be ~1 in the beginning.