Porting https://itk.org/Wiki/ITK/Examples/Metrics/MeanSquaresImageToImageMetric from C++ to Python

Beginner Questions
1

Hi all, this is my first question here, I hope I understood the rules.

I have been using itk (C++ version) for many years now, so I am quite familiar with it. Recently I have been exploring the possibility to use it also in Python. So far with not a lot of luck. Here is an specific question but any link to itk or simple itk resources only will be most appreciated (so far I have found quite a few tutorials but after reading those I feel I do not yet understand the way itk works in Python).

The example:

I want to port the code from this example:

https://itk.org/Wiki/ITK/Examples/Metrics/MeanSquaresImageToImageMetric

to Python. As you will see, I am using both itk and simpleitk (I am not sure that is a good idea actually).

So far I have tried this:

import SimpleITK as sitk
import itk

def computeMetrics(nameFixed,nameMoving):

    pixelType = sitk.sitkFloat32

    fixedImage = sitk.ReadImage(nameFixed, pixelType)

    movingImage = sitk.ReadImage(nameMoving, pixelType)

    metric=itk.MeanSquaresImageToImageMetric()
    interpolator=itk.LinearInterpolateImageFunction()
    transform=itk.TranslationTransform()

    interpolator.SetInputImage( fixedImage )

    metric.SetFixedImage( fixedImage )
    metric.SetMovingImage( movingImage )
    metric.SetFixedImageRegion( fixedImage.GetLargestPossibleRegion() )
    metric.SetTransform( transform )
    metric.SetInterpolator( interpolator )

…

However, whenever I try to set the input images for the metric (both fixed and moving) I get the following error:

Traceback (most recent call last):
File “./computeImageMetric.py”, line 67, in
computeMetrics(sys.argv[1],sys.argv[2])
File “./computeImageMetric.py”, line 26, in computeMetrics
interpolator.SetInputImage( fixedImage )
TypeError: in method ‘itkImageFunctionISS2DD_SetInputImage’, argument 2 of type ‘itkImageSS2 const *’
I do not understand the error and am also having a hard time finding documentation (in python) for the classes involved. Any help will be most appreciated.

2

The short answer is: you cannot mix ITK and SimpleITK like that. metric expects and ITKImage, and SimpleITK image is similar but different class. I guess you will have to pick either ITK or SimpleITK.

@blowekamp @zivy Correct me if I am wrong.

3

Thanks a lot for the answer.

I did worry a little bit about that and tried loading the images with itk (and not simple itk), for the following code:

def computeMetrics(nameFixed,nameMoving):

    pixelType = sitk.sitkFloat32

    #fixedImage = sitk.ReadImage(nameFixed, pixelType)
    fixedImage =itk.ImageFileReader.New(FileName=nameFixed).GetOutput()

    movingImage =itk.ImageFileReader.New(FileName=nameMoving).GetOutput()
    #movingImage = sitk.ReadImage(nameMoving, pixelType)

    metric=itk.MeanSquaresImageToImageMetric()
    interpolator=itk.LinearInterpolateImageFunction()
    transform=itk.TranslationTransform()

    interpolator.SetInputImage( fixedImage )

    metric.SetFixedImage( fixedImage )
    metric.SetMovingImage( movingImage )
    metric.SetFixedImageRegion( fixedImage.GetLargestPossibleRegion() )
    metric.SetTransform( transform )
    metric.SetInterpolator( interpolator )

The error seems to persist, though

Traceback (most recent call last):
File “./computeImageMetric.py”, line 68, in
computeMetrics(sys.argv[1],sys.argv[2])
File “./computeImageMetric.py”, line 27, in computeMetrics
interpolator.SetInputImage( fixedImage )
TypeError: in method ‘itkImageFunctionISS2DD_SetInputImage’, argument 2 of type ‘itkImageSS2 const *’

4

Hello,

You may find this SimpleITK Example useful: https://simpleitk.readthedocs.io/en/master/Examples/ImageRegistrationMethodExhaustive/Documentation.html

It utilizes an exhaustive optimizer to evaluate the metric on a grid pattern. Additionally, you may be interested in the sitk::ImageRegistrationMethod::MetricEvaluate method.

1 Like
5

Hello @nicill,

I wouldn’t recommend mixing and matching ITK and SimpleITK, but it should be easy to do via their numpy array import export facilities:

import SimpleITK as sitk
import itk

in_file_name = 'my_image.jpg'

def itk2sitk(img):
    return sitk.GetImageFromArray(itk.GetArrayViewFromImage(img), img.GetNumberOfComponentsPerPixel()>1)

def sitk2itk(img):
    return itk.GetImageFromArray(sitk.GetArrayViewFromImage(img), img.GetNumberOfComponentsPerPixel()>1)


# start with ITK, do some SimpleITK stuff, go back to ITK
img_itk = itk.imread(in_file_name)
img_sitk = sitk.Median(itk2sitk(img_itk))
itk.imwrite(sitk2itk(img_sitk), 'itk2sitk2itk.jpg')


# start with SimpleITK, do some ITK stuff, go back to SimpleITK
img_sitk = sitk.ReadImage(in_file_name)
median_filter = itk.MedianImageFilter.New(sitk2itk(img_sitk))
median_filter.Update()
img_itk = median_filter.GetOutput()
sitk.WriteImage(itk2sitk(img_itk), 'sitk2itk2sitk.jpg')

This code works for grayscale images, but I believe I uncovered a bug in ITK Python with respect to vector images. So, if nothing else, you have helped improve the toolkits.

2 Likes
6

As a side note, porting the example to Python and contributing to the ITK and/or SimpleITK Examples would be a great contribution:
https://itk.org/ITKExamples/Documentation/Contribute/index.html
https://simpleitk.readthedocs.io/en/master/Examples/index.html

Both toolkits are built in a large part thanks to the great contributions from the community. And this would definitely be such one :100:.

2 Likes
7

Thanks for your answer, it helped a lot.

I am not sure that is what you meant but your post made me think of computing the metric by setting up a registration scenario and then getting the metric (with the sitk::ImageRegistrationMethod::MetricEvaluate method that you mention by calling it before the registration is executed.

It looks like this:

def computeMetric(fixed,moving,mode):

# First, set up "phony" registration
R = sitk.ImageRegistrationMethod()
R.SetOptimizerAsGradientDescentLineSearch(learningRate=1.0,numberOfIterations=1,convergenceMinimumValue=1e-5,convergenceWindowSize=5)
R.SetInitialTransform(sitk.Transform(2,sitk.sitkIdentity)) # Transformation deliberately not using any initializer
R.SetInterpolator(sitk.sitkLinear)

#second, choose metric
if(mode==0): R.SetMetricAsJointHistogramMutualInformation()
elif(mode==1): R.SetMetricAsMattesMutualInformation(numberOfHistogramBins = 50)
elif(mode==2): R.SetMetricAsMeanSquares()
elif(mode==3): R.SetMetricAsCorrelation()
else:
    print("Error in computeMetric, unrecognized metric")
    sys.exit ( 1 )

#third, get the metric value
print(R.MetricEvaluate(fixed, moving))

I have been testing a few things (metric with same image, when the image order changes, after some registration methods) and I am getting some weird values from the mutual information metrics (modes 0 and 1) but least squares and correlation seem to work fine.

1 Like
8

Thanks for the tip about working with itk and simpleItk, at this moment I can do all I need to do with only simpleItk but I will definetely keep it mind.

9

Very nice!

That is how the “MetricEvaluate” method is designed to be used. It is not so much a “phony” registration as to evaluate the metric the interpolator and transforms are needed. Only the optimizer is not used since the registration method is just being evaluated at one set of transform parameters.