Hi there, <br><br>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:<br><br>
<div style="margin-left: 40px;">#include "itkSubtractImageFilter.h"<br>#include "itkRescaleIntensityImageFilter.h"<br>#include "itkIdentityTransform.h"<br><br>int main( int argc, char *argv[] )<br>
{<br> if( argc < 4 )<br> {<br> std::cerr << "Missing Parameters " << std::endl;<br> std::cerr << "Usage: " << argv[0];<br> std::cerr << " fixedImageFile movingImageFile ";<br>
std::cerr << " outputImagefile [differenceBeforeRegistration] ";<br> std::cerr << " [differenceAfterRegistration] ";<br> std::cerr << " [steplength] ";<br> std::cerr << " [initialScaling] [initialAngle] ";<br>
std::cerr << std::endl;<br> return EXIT_FAILURE;<br> }<br></div> <br>The error message I get is:<br><br><div style="margin-left: 40px;">/usr/lib/pymodules/python2.6/itkImageRegistrationMethod.pyc in StartRegistration(*args)<br>
1393 def SetTransform(*args): return _itkImageRegistrationMethod.itkImageRegistrationMethodIUS2IUS2_Pointer_SetTransform(*args)<br> 1394 def StartOptimization(*args): return _itkImageRegistrationMethod.itkImageRegistrationMethodIUS2IUS2_Pointer_StartOptimization(*args)<br>
-> 1395 def StartRegistration(*args): return _itkImageRegistrationMethod.itkImageRegistrationMethodIUS2IUS2_Pointer_StartRegistration(*args)<br> 1396 def AbortGenerateDataOff(*args): return _itkImageRegistrationMethod.itkImageRegistrationMethodIUS2IUS2_Pointer_AbortGenerateDataOff(*args)<br>
1397 def AbortGenerateDataOn(*args): return _itkImageRegistrationMethod.itkImageRegistrationMethodIUS2IUS2_Pointer_AbortGenerateDataOn(*args)<br><br>RuntimeError: /tmp/buildd/insighttoolkit-3.18.0/Code/Numerics/itkRegularStepGradientDescentBaseOptimizer.cxx:197:<br>
itk::ERROR: RegularStepGradientDescentOptimizer(0xad25f0e8): The size of Scales is 6, but the NumberOfParameters for the CostFunction is 4.<br></div><br>I'm not familiar with C++ so it's been quite a bit of trial and error so far. <br>
<br>Has anyone else tried to do this or write a python script which does image registration with scaling?<br><br>Here's the python code so far:<br><br>Thanks in advance, <br>Helen<br><br># For scaling<br><br>import itk<br>
<br>def updateProgress():<br> """<br> Updates progress of the registration process.<br> """<br> currentParameter = self.transform.GetParameters()<br> 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))<br>
RegistrationFilters.RegistrationFilter.updateProgress(self)<br><br>image = itk.Image[itk.US, 2].New()<br><br>ffix = 'trial1.jpg'<br>fmove = 'trial2.jpg'<br>outf = 'outtrial.png' <br>useMoments = 'CENTEROFMASS'<br>
<br>maxStepLength = 0.1<br>minStepLength = 0.0001<br>maxIterations = 500<br>backgroundValue = 100<br><br>fixedImageReader = itk.ImageFileReader[image].New()<br>movingImageReader = itk.ImageFileReader[image].New()<br>fixedImageReader.SetFileName( ffix)<br>
movingImageReader.SetFileName(fmove)<br>fixedImageReader.Update()<br>movingImageReader.Update()<br><br>fixedImage = fixedImageReader.GetOutput() <br>movingImage = movingImageReader.GetOutput()<br><br><br># Create registration framework's components<br>
registration = itk.ImageRegistrationMethod[image,image].New()<br>optimizer = itk.RegularStepGradientDescentOptimizer.New()<br>metric = itk.MeanSquaresImageToImageMetric[image,image].New()<br>interpolator = itk.LinearInterpolateImageFunction[ itk.Image[ itk.US, 2], itk.D ].New()<br>
transform = itk.Similarity2DTransform.D.New() # or Centered?<br><br># Initialize registration framework's components<br>registration.SetOptimizer(optimizer.GetPointer())<br>registration.SetTransform(transform.GetPointer())<br>
registration.SetInterpolator(interpolator.GetPointer())<br>registration.SetMetric(metric.GetPointer())<br>registration.SetFixedImage(fixedImage)<br>registration.SetMovingImage(movingImage)<br>registration.SetFixedImageRegion(fixedImage.GetBufferedRegion())<br>
<br>optimizer.SetMaximumStepLength(maxStepLength)<br>optimizer.SetMinimumStepLength(minStepLength)<br>optimizer.SetNumberOfIterations(maxIterations)<br><br>################## Set optimizer scales#######################<br>
optimizerScales = itk.Array.D()<br>optimizerScales.SetSize(6)<br><br>translationScale = 1.0/100.0<br>scaleScale = 10<br>rotationScale = 1<br><br>optimizerScales.SetElement(0,scaleScale)<br>optimizerScales.SetElement(1,rotationScale)<br>
for i in range (2,5):<br> optimizerScales.SetElement(i,translationScale)<br><br>optimizer.SetScales(optimizerScales)<br>#optimizer.SetRelaxationFactor( 0.8 )<br><br>######################################## Initialise ################################<br>
<br><br>#scale = transform.GetParameters()<br>#transform.SetAngle(0)<br>#transform.SetScale(1)<br><br># Initialize transform<br>initializer = itk.CenteredTransformInitializer[transform, fixedImage, movingImage].New()<br>initializer.SetTransform(transform)<br>
initializer.SetFixedImage(fixedImage)<br>initializer.SetMovingImage(movingImage)<br><br>if useMoments == 'CENTEROFMASS':<br> initializer.MomentsOn()<br>else:<br> initializer.GeometryOn()<br>initializer.InitializeTransform()<br>
<br>######################################## Check what's happening - must have same number of parameters as scales? ################################<br><br>registration.SetInitialTransformParameters(transform.GetParameters())<br>
<br>#<br># Iteration Observer<br>#<br>def iterationUpdate():<br> currentParameter = transform.GetParameters()<br> print "%f %f %f %f %f %f" % ( currentParameter.GetElement(0),<br> currentParameter.GetElement(1),<br>
currentParameter.GetElement(2),<br> currentParameter.GetElement(3),<br> currentParameter.GetElement(4),<br> currentParameter.GetElement(5) )<br>
<br>iterationCommand = itk.PyCommand.New()<br>iterationCommand.SetCommandCallable( iterationUpdate )<br>optimizer.AddObserver( itk.IterationEvent(), iterationCommand.GetPointer() )<br><br><br>######################################## Go ################################<br>
<br>print "Starting registration"<br><br>#<br># Start the registration process<br>#<br><br>registration.StartRegistration()<br><br>#<br># Get the final parameters of the transformation<br>#<br>finalParameters = registration.GetLastTransformParameters()<br>
<br>print "Final Registration Parameters "<br>print "Scale = %f" % finalParameters.GetElement(0)<br>print "Angle in radians = %f" % finalParameters.GetElement(1)<br>print "Rotation Center X = %f" % finalParameters.GetElement(2)<br>
print "Rotation Center Y = %f" % finalParameters.GetElement(3)<br>print "Translation in X = %f" % finalParameters.GetElement(4)<br>print "Translation in Y = %f" % finalParameters.GetElement(5)<br>
<br>################# Resample, cast, subtract, rescale, identity ###############<br><br># Now, we use the final transform for resampling the moving image.<br>resampler = itk.ResampleImageFilter[image,image].New()<br>resampleInterpolator = itk.NearestNeighborInterpolateImageFunction[ itk.Image[ itk.US, 2], itk.D ].New()<br>
transform.SetParameters(finalParameters)<br>resampler.SetTransform(transform.GetPointer())<br>resampler.SetInput(movingImage)<br>resampler.SetSize(movingImage.GetLargestPossibleRegion().GetSize())<br>resampler.SetOutputSpacing(movingImage.GetSpacing())<br>
resampler.SetOutputDirection(movingImage.GetDirection())<br>resampler.SetOutputOrigin(movingImage.GetOrigin())<br>resampler.SetDefaultPixelValue(backgroundValue)<br>resampler.SetInterpolator(resampleInterpolator.GetPointer())<br>
data = resampler.GetOutput()<br>data.Update()<br><br># Cast for output<br>#<br>outputCast = itk.RescaleIntensityImageFilter[image,image].New()<br>outputCast.SetInput( resampler.GetOutput() )<br>outputCast.SetOutputMinimum( 0 )<br>
outputCast.SetOutputMaximum( 65535 )<br><br># Subtract filter?<br><br># Rescale intensity?<br><br># Identity transform?<br><br>writer = itk.ImageFileWriter[image].New()<br><br>writer.SetFileName( outf)<br>writer.SetInput( outputCast.GetOutput() )<br>
writer.Update()<br><br><br><br><br><br><br><br><br>