[Insight-users] About registration and the SoftwareGuide
Luis Ibanez
luis.ibanez at kitware.com
Sat Mar 17 10:28:23 EST 2007
Hi Xavier,
Thanks for your suggestion.
We will insert the comments in the next
edition of the ITK software guide.
Regards,
Luis
--------------------
BATY Xavier wrote:
> Hi all,
>
> As I'm working on a registration process, I have the following suggestion.
> The 2 messages (cited below) about registration parameters tunning and
> the Too many samples maps outside... execption sent by Luis recently
> should be added to registration chapter of the software guide.
>
> Those messages contains answers to many questions for people who start
> working on registration.
>
> XB
>
>> ------------------------------
>>
>> Message: 3
>> Date: Thu, 15 Mar 2007 21:26:39 -0400
>> From: "Luis Ibanez" <luis.ibanez at kitware.com>
>> Subject: Re: [Insight-users] Too many samples map outside moving image
>> buffer
>> To: Goo <gtshowtime at gmail.com>
>> Cc: insight-users at itk.org
>> Message-ID:
>> <f7abd23c0703151826j6b7c23cfwb6ffe0c9583039dd at mail.gmail.com>
>> Content-Type: text/plain; charset=ISO-8859-1; format=flowed
>>
>> Hi Goo,
>>
>> This is a common error message in image registration.
>>
>> It means that at the current iteration of the optimization,
>> the two images as so off-registration that their spatial
>> overlap is not large enough for bringing them back into
>> registration.
>>
>> The common causes of this problem are:
>>
>> 1) Poor initialization: You must initialize the transform
>> properly. Please read the ITK Software Guide
>> http://www.itk.org/ItkSoftwareGuide.pdf for a description
>> of the use of the CenteredTransformInitializer class.
>>
>>
>> 2) Optimzer steps too large. If you optimizer takes steps
>> that are too large, it risks to become unstable and to
>> send the images too far appart.
>>
>> You may want to start the optimizer with a maximum
>> step lenght of 1.0, and only increase it once you have
>> managed to fine tune all other registration parameters.
>>
>> Increasing the step length makes your program faster,
>> but it also makes it more unstable.
>>
>>
>> 3) Poor set up o the transform parameters scaling.
>>
>> This is extremely critical in registration. You must make
>> sure that you balance the relative difference of scale between
>> the rotation parameters and the translation parameters.
>>
>> In typical medical datasets such as CT and MR, translations
>> are measured in millimeters, and therefore are in the range
>> of -100:100, while rotations are measured in radians, and
>> therefore they tend to be in the range of -1:1.
>>
>>
>> A rotation of 3 radians is catastrophic, while a translation
>> of 3 millimeters is rather inoffensive. That difference in scale
>> is the one that must be accounted for.
>>
>>
>>
>>
>> Regards,
>>
>>
>> Luis
>>
>>
>>
>>
>> --------------------------------------------------------------------
>> On 3/13/07, Goo <gtshowtime at gmail.com> wrote:
>>
>>
>>> Hi All :
>>>
>>> I met an error message when using Mattes MI and Affine Transform.
>>> This message is notifying me "Too many samples map outside moving image
>>> buffer"
>>>
>>> I guess this maybe caused by TransformPoint in
>>> itkMattesMutualInformationImageToImageMetric.txx
>>>
>>> this->TransformPoint( nFixedImageSamples, parameters, mappedPoint,
>>> sampleOk, movingImageValue );
>>>
>>> but I don't know how to solve this problem.
>>>
>>> My test data are two different modality images with the size 512*512.
>>>
>>> The Metric setting :
>>> unsigned int numberOfBins = 50;
>>> unsigned int numberOfSamples = 10000;
>>>
>>> The Transform setting :
>>> initialParameters[0] = 0.0; // Initial offset in mm along X
>>> initialParameters[1] = 0.0; // Initial offset in mm along Y
>>> registration->SetInitialTransformParameters( initialParameters );
>>>
>>> The optimizer setting :
>>> optimizer->MinimizeOn();
>>> optimizer->SetMaximumStepLength( 5.00 );
>>> optimizer->SetMinimumStepLength( 0.001 );
>>> optimizer->SetNumberOfIterations( 200 );
>>> optimizer->SetRelaxationFactor( 0.8 );
>>>
>>> Can anybody help me to solve this problem?
>>>
>>> Regards.
>>>
>>
>>
>>
>> ------------------------------
>>
>> Message: 4
>> Date: Thu, 15 Mar 2007 21:05:18 -0400
>> From: "Luis Ibanez" <luis.ibanez at kitware.com>
>> Subject: [Insight-users] Re: How to choose user-defined parameter for
>> combining Mattes MI metric and Affine transform?
>> To: Goo <gtshowtime at gmail.com>
>> Cc: insight-users at itk.org
>> Message-ID:
>> <f7abd23c0703151805p62fd964dr2db34e0a09bdf394 at mail.gmail.com>
>> Content-Type: text/plain; charset=ISO-8859-1; format=flowed
>>
>> Hi Goo
>>
>> Here is some advice on how to fine tune the parameters
>> of the registration process.
>>
>>
>> 1) Set Maximum step length to 0.1 and do not change it
>> until all other parameters are stable.
>>
>> 2) Set Minimum step length to 0.001 and do not change it.
>>
>> You could interpret these two parameters as if their
>> units were radians. So, 0.1 radian = 5.7 degrees.
>>
>>
>> 3) Number of histogram bins:
>>
>> First plot the histogram of your image using the
>> example program in
>>
>> Insight/Examples/Statistics/ImageHistogram2.cxx
>>
>> In that program use first a large number of bins
>> (for example 2000) and identify the different
>> populations of intensity level and to what anatomical
>> structures they correspond.
>>
>> Once you identify the anatomical structures in the
>> histogram, then rerun that same program with less
>> and less number of bins, until you reach the minimun
>> number of bins for which all the tissues that are important
>> for your application, are still distinctly differentiated in the
>> histogram. At that point, take that number of bins and
>> us it for your Mutual Information metric.
>>
>>
>> 4) Number of Samples:
>>
>> The trade-off with the number of samples is the following:
>>
>> a) computation time of registration is linearly proportional
>> to the number of samples
>> b) the samples must be enough to significantly populate
>> the joint histogram.
>> c) Once the histogram is populated, there is not much
>> use in adding more samples.
>>
>> Therefore do the following:
>>
>> Plot the joint histogram of both images, using the number
>> of bins that you selected in item (3). You can do this by
>> modifying the code of the example:
>>
>> Insight/Examples/Statistics/
>> ImageMutualInformation1.cxx
>>
>> you have to change the code to print out the values
>> of the bins. Then use a plotting program such as gnuplot,
>> or Matlab, or even Excel and look at the distribution.
>>
>> The number of samples to take must be enough
>> for producing the same "appearance" of the joint histogram.
>>
>> As an arbitrary rule of thumb you may want to start using
>> a high number of samples (80% ~ 100%). And do not
>> change it until you have mastered the other parameters
>> of the registration. Once you get your registration to converge
>> you can revisit the number of samples and reduce it in order
>> to make the registration run faster. You can simply reduce it
>> until you find that the registration becomes unstable. That's
>> your critical bound for the minimum number of samples.
>> Take that number and multiply it by the magic number 1.5,
>> to send it back to a stable region, or if your application is
>> really critical, then use an even higher magic number x2.0.
>>
>> This is just engineering: you figure out what is the minimal
>> size of a piece of steel that will support a bridge, and then
>> you enlarge it to keep it away from the critical value.
>>
>>
>>
>> 5) The MOST critical values of the registration process are the
>> scaling parameters that define the proportions between
>> the parameters of the transform. In your case, for an Affine
>> Transform in 2D, you have 6 parameters. The first four are
>> the ones of the Matrix, and the last two are the translation.
>>
>> The rotation matrix value must be in the ranges of radians
>> which is typically [ -1 to 1 ], while the translation values are
>> in the ranges of millimeters (your image size units).
>>
>> You want to start by setting the scaling of the matrix
>> parameters to 1.0, and the scaling of the Translation
>> parameters to the holy esoteric values:
>>
>> 1.0 / ( 10.0 * pixelspacing[0] * imagesize[0] )
>> 1.0 / ( 10.0 * pixelspacing[1] * imagesize[1] )
>>
>> This is telling the optimizer that you consider that rotating
>> the image by 57 degrees is as "significant" as translating
>> the image by half its physical extent.
>>
>> Note that esoteric value has included the arbitrary number
>> 10.0 in the denominator, for no other reason that we have
>> been lucky when using that factor. This of course is just a
>> supersticion, so you should feel free to experiment with
>> different values of this number.
>>
>> Just keep in mind that what the optimizer will do is to
>> "jump" in a paramteric space of 6 dimension, and that the
>> component of the jump on every dimension will be proporitional
>> to 1/scaling factor * OptimizerStepLenght. Since you put
>> the optimizer Step Length to 0.1, then the optimizer will start
>> by exploring the rotations at jumps of about 5degrees, which
>> is a conservative rotation for most medical applications.
>>
>> If you have reasons to think that your rotations are larger or
>> smaller, then you should modify the scaling factor of the matrix
>> parameters accordingly.
>>
>> In the same way, if you thinkl that 1/10 of the image size is too
>> large as the first step for exploring the translations, then you
>> should modify the scaling of translation parameters accordingly.
>>
>>
>>
>> In order to drive all these you need to analyze the feedback that
>> the observer is providing you. For example, plot the metric values,
>> and plot the translation coordinates so that you can get a feeling
>> of how the registration is behaving.
>>
>>
>> Note also that image registration is not a science. it is a pure
>> engineerig practice, and therefore, there are no correct answers,
>> nor "truths" to be found. It is all about how much quality your want,
>> and how must computation time, and development time are you
>> willing to pay for that quality. The "satisfying" answer for your
>> specific application must be found by exploring the trade-offs
>> between the different parameters that regulate the image
>> registration process.
>>
>> If you are proficient in VTK you may want to consider attaching
>> some visualization to the Event observer, so that you can have
>> a visual feedback on the progress of the registration. This is a
>> lot more productive than trying to interpret the values printed
>> out on the console by the observer.
>>
>>
>> Regards,
>>
>>
>> Luis
>>
>>
>>
>> =================================
>> On 3/14/07, Goo <gtshowtime at gmail.com> wrote:
>>
>>
>>> > Hi, All :
>>> >
>>> > I want to combine Affine transform and Mattes MI metric into my
>>> registration
>>> > process.
>>> > The process is fine to work but the result of registration is
>>> incorrect.
>>> > I think this is caused by some coefficients setting that is not
>>> suitable for
>>> > my application,
>>> > but how to determine a applicable coefficients that is very
>>> suffering for
>>> > me.
>>> > So, I wish any experienced friends can give me some indication that
>>> > including how to set:
>>> > OptimizerScales (specially translation Scale) <---How to set ???
>>> > MaximumStepLength <---How to set ???
>>> > MinimumStepLength <---How to set ???
>>> > NumberOfHistogramBins (now I set 20~50 in typically)
>>> > NumberOfSpatialSamples (now I set 20% of image pixels )
>>> >
>>> >
>>> > My code is posting as below:
>>> > Regards.
>>> >
>>> > class CommandIterationUpdate : public itk::Command
>>> > {
>>> > public:
>>> > typedef CommandIterationUpdate Self;
>>> > typedef itk::Command Superclass;
>>> > typedef itk::SmartPointer<Self> Pointer;
>>> > itkNewMacro( Self );
>>> > protected:
>>> > CommandIterationUpdate() {};
>>> > public:
>>> > typedef itk::RegularStepGradientDescentOptimizer OptimizerType;
>>> > typedef const OptimizerType * OptimizerPointer;
>>> >
>>> > void Execute(itk::Object *caller, const itk::EventObject & event)
>>> > {
>>> > Execute( (const itk::Object *)caller, event);
>>> > }
>>> >
>>> > void Execute(const itk::Object * object, const itk::EventObject &
>>> event)
>>> > {
>>> > OptimizerPointer optimizer =
>>> > dynamic_cast< OptimizerPointer >( object );
>>> > if( ! itk::IterationEvent().CheckEvent( &event ) )
>>> > {
>>> > return;
>>> > }
>>> > std::cout << optimizer->GetCurrentIteration() << " ";
>>> > std::cout << optimizer->GetValue() << " ";
>>> > std::cout << optimizer->GetCurrentPosition();
>>> >
>>> > // Print the angle for the trace plot
>>> > vnl_matrix<double> p(2, 2);
>>> > p[0][0] = (double) optimizer->GetCurrentPosition()[0];
>>> > p[0][1] = (double) optimizer->GetCurrentPosition()[1];
>>> > p[1][0] = (double) optimizer->GetCurrentPosition()[2];
>>> > p[1][1] = (double) optimizer->GetCurrentPosition()[3];
>>> > vnl_svd<double> svd(p);
>>> > vnl_matrix<double> r(2, 2);
>>> > r = svd.U() * vnl_transpose(svd.V());
>>> > double angle = asin(r[1][0]);
>>> > std::cout << " AffineAngle: " << angle * 45.0 / atan(1.0) <<
>>> > std::endl;
>>> > }
>>> > };
>>> >
>>> > int main( int argc, char *argv[] )
>>> > {
>>> > if( argc < 4 )
>>> > {
>>> > std::cerr << "Missing Parameters " << std::endl;
>>> > std::cerr << "Usage: " << argv[0];
>>> > std::cerr << " fixedImageFile movingImageFile
>>> outputImagefile" <<
>>> > std::endl;
>>> > return 1;
>>> > }
>>> >
>>> > const unsigned int Dimension = 2;
>>> > typedef float PixelType;
>>> >
>>> > typedef itk::Image< PixelType, Dimension > FixedImageType;
>>> > typedef itk::Image< PixelType, Dimension > MovingImageType;
>>> >
>>> > typedef itk::AffineTransform<
>>> > double,
>>> > Dimension > TransformType;
>>> >
>>> > typedef itk::RegularStepGradientDescentOptimizer
>>> OptimizerType;
>>> >
>>> > typedef itk::MattesMutualInformationImageToImageMetric<
>>> > FixedImageType,
>>> > MovingImageType >
>>> MetricType;
>>> >
>>> > typedef itk:: LinearInterpolateImageFunction<
>>> > MovingImageType,
>>> > double >
>>> InterpolatorType;
>>> > typedef itk::ImageRegistrationMethod<
>>> > FixedImageType,
>>> > MovingImageType >
>>> RegistrationType;
>>> >
>>> > MetricType::Pointer metric = MetricType::New();
>>> > OptimizerType::Pointer optimizer = OptimizerType::New();
>>> > InterpolatorType::Pointer interpolator = InterpolatorType::New();
>>> > RegistrationType::Pointer registration = RegistrationType::New();
>>> >
>>> > registration->SetMetric( metric );
>>> > registration->SetOptimizer( optimizer );
>>> > registration->SetInterpolator( interpolator );
>>> >
>>> > TransformType::Pointer transform = TransformType::New();
>>> > registration->SetTransform( transform );
>>> >
>>> > typedef itk::ImageFileReader< FixedImageType >
>>> FixedImageReaderType;
>>> > typedef itk::ImageFileReader< MovingImageType >
>>> MovingImageReaderType;
>>> > FixedImageReaderType::Pointer fixedImageReader =
>>> > FixedImageReaderType::New();
>>> > MovingImageReaderType::Pointer movingImageReader =
>>> > MovingImageReaderType::New();
>>> > fixedImageReader->SetFileName( argv[1] );
>>> > movingImageReader->SetFileName( argv[2] );
>>> >
>>> >
>>> > registration->SetFixedImage( fixedImageReader->GetOutput() );
>>> > registration->SetMovingImage( movingImageReader->GetOutput() );
>>> > fixedImageReader->Update();
>>> >
>>> > registration->SetFixedImageRegion(
>>> > fixedImageReader->GetOutput()->GetBufferedRegion() );
>>> >
>>> > typedef itk::CenteredTransformInitializer<
>>> > TransformType,
>>> > FixedImageType,
>>> > MovingImageType >
>>> > TransformInitializerType;
>>> > TransformInitializerType::Pointer initializer =
>>> > TransformInitializerType::New();
>>> > initializer->SetTransform( transform );
>>> > initializer->SetFixedImage( fixedImageReader->GetOutput() );
>>> > initializer->SetMovingImage( movingImageReader->GetOutput() );
>>> > initializer->MomentsOn();
>>> > initializer->InitializeTransform();
>>> >
>>> > registration->SetInitialTransformParameters(
>>> > transform->GetParameters() );
>>> >
>>> > typedef OptimizerType::ScalesType OptimizerScalesType;
>>> > OptimizerScalesType optimizerScales(
>>> transform->GetNumberOfParameters() );
>>> >
>>> >
>>> //-----------------------------------------------------------------------------
>>>
>>> > // How to determine these ?
>>> > metric->SetNumberOfHistogramBins( 50 );
>>> > metric->SetNumberOfSpatialSamples( 10000 );
>>> > metric->ReinitializeSeed( 76926294 );
>>> >
>>> > double translationScale = 1.0 / 1000.0;
>>> > optimizerScales[0] = 1.0;
>>> > optimizerScales[1] = 1.0;
>>> > optimizerScales[2] = 1.0;
>>> > optimizerScales[3] = 1.0;
>>> > optimizerScales[4] = translationScale;
>>> > optimizerScales[5] = translationScale;
>>> > optimizer->SetScales( optimizerScales );
>>> >
>>> > double steplength = 1; //0.1
>>> > unsigned int maxNumberOfIterations = 200;
>>> > optimizer->SetMaximumStepLength( steplength );
>>> > optimizer->SetMinimumStepLength( 0.0001 ); //0.0001
>>> > optimizer->SetNumberOfIterations( maxNumberOfIterations );
>>> >
>>> //-----------------------------------------------------------------------------
>>>
>>> >
>>> >
>>> > optimizer->MinimizeOn();
>>> >
>>> >
>>> > CommandIterationUpdate::Pointer observer =
>>> CommandIterationUpdate::New();
>>> > optimizer->AddObserver( itk::IterationEvent(), observer );
>>> >
>>> > try
>>> > {
>>> > registration->StartRegistration();
>>> > }
>>> > catch( itk::ExceptionObject & err )
>>> > {
>>> > std::cerr << "ExceptionObject caught !" << std::endl;
>>> > std::cerr << err << std::endl;
>>> > return -1;
>>> > }
>>> >
>>> >
>>> >
>>> > //-----------------------------------------------------------
>>> > OptimizerType::ParametersType finalParameters =
>>> > registration->GetLastTransformParameters();
>>> >
>>> > const double finalRotationCenterX = transform->GetCenter()[0];
>>> > const double finalRotationCenterY = transform->GetCenter()[1];
>>> > const double finalTranslationX = finalParameters[4];
>>> > const double finalTranslationY = finalParameters[5];
>>> >
>>> > const unsigned int numberOfIterations =
>>> optimizer->GetCurrentIteration();
>>> >
>>> > const double bestValue = optimizer->GetValue();
>>> > std::cout << "Result = " << std::endl;
>>> > std::cout << " Center X = " << finalRotationCenterX <<
>>> std::endl;
>>> > std::cout << " Center Y = " << finalRotationCenterY <<
>>> std::endl;
>>> > std::cout << " Translation X = " << finalTranslationX << std::endl;
>>> > std::cout << " Translation Y = " << finalTranslationY << std::endl;
>>> > std::cout << " Iterations = " << numberOfIterations << std::endl;
>>> > std::cout << " Metric value = " << bestValue << std::endl;
>>> >
>>>
>>
>> ------------------------------
>
>
> _______________________________________________
> Insight-users mailing list
> Insight-users at itk.org
> http://www.itk.org/mailman/listinfo/insight-users
>
More information about the Insight-users
mailing list