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.