[Insight-users] Python ITK Registration Example with Scaling

[Helen Ramsden]

Hi there,

I'm attempting to convert ITK's example ImageRegistration7 into python but
I'm having a problem defining the cost function and working out how to slot
in the functions/bits of code below:

#include "itkSubtractImageFilter.h"
#include "itkRescaleIntensityImageFilter.h"
#include "itkIdentityTransform.h"

int main( int argc, char *argv[] )
{
  if( argc < 4 )
    {
    std::cerr << "Missing Parameters " << std::endl;
    std::cerr << "Usage: " << argv[0];
    std::cerr << " fixedImageFile  movingImageFile ";
    std::cerr << " outputImagefile  [differenceBeforeRegistration] ";
    std::cerr << " [differenceAfterRegistration] ";
    std::cerr << " [steplength] ";
    std::cerr << " [initialScaling] [initialAngle] ";
    std::cerr << std::endl;
    return EXIT_FAILURE;
    }

The error message I get is:

/usr/lib/pymodules/python2.6/itkImageRegistrationMethod.pyc in
StartRegistration(*args)
   1393     def SetTransform(*args): return
_itkImageRegistrationMethod.itkImageRegistrationMethodIUS2IUS2_Pointer_SetTransform(*args)
   1394     def StartOptimization(*args): return
_itkImageRegistrationMethod.itkImageRegistrationMethodIUS2IUS2_Pointer_StartOptimization(*args)
-> 1395     def StartRegistration(*args): return
_itkImageRegistrationMethod.itkImageRegistrationMethodIUS2IUS2_Pointer_StartRegistration(*args)
   1396     def AbortGenerateDataOff(*args): return
_itkImageRegistrationMethod.itkImageRegistrationMethodIUS2IUS2_Pointer_AbortGenerateDataOff(*args)
   1397     def AbortGenerateDataOn(*args): return
_itkImageRegistrationMethod.itkImageRegistrationMethodIUS2IUS2_Pointer_AbortGenerateDataOn(*args)

RuntimeError:
/tmp/buildd/insighttoolkit-3.18.0/Code/Numerics/itkRegularStepGradientDescentBaseOptimizer.cxx:197:
itk::ERROR: RegularStepGradientDescentOptimizer(0xad25f0e8): The size of
Scales is 6, but the NumberOfParameters for the CostFunction is 4.

I'm not familiar with C++ so it's been quite a bit of trial and error so
far.

Has anyone else tried to do this or write a python script which does image
registration with scaling?

Here's the python code so far:

Thanks in advance,
Helen

# For scaling

import itk

def updateProgress():
    """
    Updates progress of the registration process.
    """
    currentParameter = self.transform.GetParameters()
    print "M: %f S: %f A: %f C: (%f,%f) T: (%f,%f)"
%(self.optimizer.GetValue(),currentParameter.GetElement(0),currentParameter.GetElement(1),currentParameter.GetElement(2),currentParameter.GetElement(3),currentParameter.GetElement(4),currentParameter.GetElement(5))
    RegistrationFilters.RegistrationFilter.updateProgress(self)

image = itk.Image[itk.US, 2].New()

ffix =  'trial1.jpg'
fmove = 'trial2.jpg'
outf =  'outtrial.png'
useMoments = 'CENTEROFMASS'

maxStepLength = 0.1
minStepLength = 0.0001
maxIterations = 500
backgroundValue = 100

fixedImageReader = itk.ImageFileReader[image].New()
movingImageReader = itk.ImageFileReader[image].New()
fixedImageReader.SetFileName( ffix)
movingImageReader.SetFileName(fmove)
fixedImageReader.Update()
movingImageReader.Update()

fixedImage = fixedImageReader.GetOutput()
movingImage = movingImageReader.GetOutput()


# Create registration framework's components
registration = itk.ImageRegistrationMethod[image,image].New()
optimizer = itk.RegularStepGradientDescentOptimizer.New()
metric = itk.MeanSquaresImageToImageMetric[image,image].New()
interpolator = itk.LinearInterpolateImageFunction[ itk.Image[ itk.US, 2],
itk.D ].New()
transform = itk.Similarity2DTransform.D.New() # or Centered?

# Initialize registration framework's components
registration.SetOptimizer(optimizer.GetPointer())
registration.SetTransform(transform.GetPointer())
registration.SetInterpolator(interpolator.GetPointer())
registration.SetMetric(metric.GetPointer())
registration.SetFixedImage(fixedImage)
registration.SetMovingImage(movingImage)
registration.SetFixedImageRegion(fixedImage.GetBufferedRegion())

optimizer.SetMaximumStepLength(maxStepLength)
optimizer.SetMinimumStepLength(minStepLength)
optimizer.SetNumberOfIterations(maxIterations)

################## Set optimizer scales#######################
optimizerScales = itk.Array.D()
optimizerScales.SetSize(6)

translationScale = 1.0/100.0
scaleScale = 10
rotationScale = 1

optimizerScales.SetElement(0,scaleScale)
optimizerScales.SetElement(1,rotationScale)
for i in range (2,5):
    optimizerScales.SetElement(i,translationScale)

optimizer.SetScales(optimizerScales)
#optimizer.SetRelaxationFactor( 0.8 )

######################################## Initialise
################################


#scale = transform.GetParameters()
#transform.SetAngle(0)
#transform.SetScale(1)

# Initialize transform
initializer = itk.CenteredTransformInitializer[transform, fixedImage,
movingImage].New()
initializer.SetTransform(transform)
initializer.SetFixedImage(fixedImage)
initializer.SetMovingImage(movingImage)

if useMoments == 'CENTEROFMASS':
    initializer.MomentsOn()
else:
    initializer.GeometryOn()
initializer.InitializeTransform()

######################################## Check what's happening - must have
same number of parameters as scales? ################################

registration.SetInitialTransformParameters(transform.GetParameters())

#
# Iteration Observer
#
def iterationUpdate():
    currentParameter = transform.GetParameters()
    print "%f %f %f %f %f %f" % ( currentParameter.GetElement(0),
                                currentParameter.GetElement(1),
                                currentParameter.GetElement(2),
                                currentParameter.GetElement(3),
                                currentParameter.GetElement(4),
                                currentParameter.GetElement(5) )

iterationCommand = itk.PyCommand.New()
iterationCommand.SetCommandCallable( iterationUpdate )
optimizer.AddObserver( itk.IterationEvent(), iterationCommand.GetPointer() )


######################################## Go ################################

print "Starting registration"

#
# Start the registration process
#

registration.StartRegistration()

#
# Get the final parameters of the transformation
#
finalParameters = registration.GetLastTransformParameters()

print "Final Registration Parameters "
print "Scale  = %f" % finalParameters.GetElement(0)
print "Angle in radians  = %f" % finalParameters.GetElement(1)
print "Rotation Center X = %f" % finalParameters.GetElement(2)
print "Rotation Center Y = %f" % finalParameters.GetElement(3)
print "Translation in  X = %f" % finalParameters.GetElement(4)
print "Translation in  Y = %f" % finalParameters.GetElement(5)

################# Resample, cast, subtract, rescale, identity
###############

# Now, we use the final transform for resampling the moving image.
resampler = itk.ResampleImageFilter[image,image].New()
resampleInterpolator = itk.NearestNeighborInterpolateImageFunction[
itk.Image[ itk.US, 2], itk.D ].New()
transform.SetParameters(finalParameters)
resampler.SetTransform(transform.GetPointer())
resampler.SetInput(movingImage)
resampler.SetSize(movingImage.GetLargestPossibleRegion().GetSize())
resampler.SetOutputSpacing(movingImage.GetSpacing())
resampler.SetOutputDirection(movingImage.GetDirection())
resampler.SetOutputOrigin(movingImage.GetOrigin())
resampler.SetDefaultPixelValue(backgroundValue)
resampler.SetInterpolator(resampleInterpolator.GetPointer())
data = resampler.GetOutput()
data.Update()

# Cast for output
#
outputCast = itk.RescaleIntensityImageFilter[image,image].New()
outputCast.SetInput( resampler.GetOutput() )
outputCast.SetOutputMinimum( 0 )
outputCast.SetOutputMaximum( 65535 )

# Subtract filter?

# Rescale intensity?

# Identity transform?

writer = itk.ImageFileWriter[image].New()

writer.SetFileName( outf)
writer.SetInput( outputCast.GetOutput() )
writer.Update()
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20110608/d12c2785/attachment.htm>
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: not available
URL: <http://www.itk.org/pipermail/insight-users/attachments/20110608/d12c2785/attachment.asc>