[Insight-users] Re: Non-rigid Registration

[Luis Ibanez]

Hi Andy,


A) In your registration process that returns the
    transform with a value equal to your initial
    guess.... for how many iterations did the optimizer
    actually run  ?


B) In your multiresolution example, I found strange that
    the spacing and origin of the wrapper filter is using
    values from the matched->Output() image, you probably
    your use the values from the fixed image.

    Also your heuristic for setting the number of iterations
    is a bit arbitrary...  you should verify if at every
    level of the pyramid, the registration stops because the
    minimum error has been achieved or because all the alloted
    iterations have been used.  If it happens to stop because
    of full use of iterations, this means that the level didn't
    really converge,... and it is very unlikely that the next
    resolution level could do better.



Regards,



    Luis



-----------------
Andy Eow wrote:

> Hi Luis,
> 
> Thanks for your advice ... i think AffineTransform is more applicable to me.
> Anyhow, I tried implementing registration using AffineTransform but it
> always returns my initial guess as the solution after registration
> completes. I've attached my code below, could you help me take a look if I'm
> setting everything up correctly?
> 
> ImageType::Pointer registerAffine( ImageType::Pointer fixedImg,
> ImageType::Pointer movingImg )
> {
>     typedef itk::AffineTransform<double, dimension> TransformType;
>     typedef itk::RegularStepGradientDescentBaseOptimizer OptimizerType;
>     typedef itk::LinearInterpolateImageFunction <ImageType, double>
> InterpolatorType;
>     typedef itk::ImageRegistrationMethod <ImageType, ImageType>
> RegistrationType;
>     typedef itk::MutualInformationImageToImageMetric<ImageType, ImageType>
> MetricType;
>     MetricType::Pointer metric = MetricType::New();
>     TransformType::Pointer transform = TransformType::New();
>     OptimizerType::Pointer optimizer = OptimizerType::New();
>     InterpolatorType::Pointer interpolator = InterpolatorType::New();
>     RegistrationType::Pointer registration  = RegistrationType::New();
> 
>     RegistrationType::ParametersType
> guess(transform->GetNumberOfParameters());
>     guess[0] = 1.0;
>     guess[1] = 0.0;
>     guess[2] = 0.0;
>     guess[3] = 1.0;
>     guess[4] = 0.0;
>     guess[5] = 0.0;
>     registration->SetInitialTransformParameters (guess);
> 
>     metric->SetMovingImageStandardDeviation(1.0);
>     metric->SetFixedImageStandardDeviation(1.0);
>     metric->SetFixedImageRegion(fixedImg->GetBufferedRegion());
>     metric->SetNumberOfSpatialSamples(1000);
> 
>     typedef OptimizerType::ScalesType ScaleType;
>     ScaleType scales(transform->GetNumberOfParameters());
>     scales.Fill(1.0);
>     optimizer->SetScales(scales);
>     optimizer->SetNumberOfIterations(300);
> 
>     registration->SetMetric(metric);
>     registration->SetOptimizer(optimizer);
>     registration->SetTransform(transform);
>     registration->SetInterpolator(interpolator);
>     registration->SetFixedImage(movingImg);
>     registration->SetMovingImage(fixedImg);
>     registration->StartRegistration();
>     RegistrationType::ParametersType solution =
> registration->GetLastTransformParameters();
>     cout << solution << endl;
> 
>     // function not complete
>     return NULL;
> }
> 
> 
> The other question is more pertaining to my old code where i was trying to
> find the deformation field using a multiresolution approach. Code is as
> follows but the registered output (after warping) is a highly shrunk version
> of the moving image at the bottom right corner of the registered image.
> Could you also help me take a look at it to see if I've setup everything
> correctly?
> 
> ImageType::Pointer registerImage( ImageType::Pointer fixedImg,
> ImageType::Pointer movingImg )
> {
>     // Cast input images to internal type
>     typedef float InternalPixelType;
>     typedef itk::Image<InternalPixelType, dimension> InternalImageType;
>     typedef itk::CastImageFilter<ImageType, InternalImageType>
> ImageCasterType;
>     ImageCasterType::Pointer fixedCaster = ImageCasterType::New();
>     fixedCaster->SetInput(fixedImg);
>     fixedCaster->Update();
>     ImageCasterType::Pointer movingCaster = ImageCasterType::New();
>     movingCaster->SetInput(movingImg);
>     movingCaster->Update();
> 
>     // Correct for intensity differences
>     typedef itk::HistogramMatchingImageFilter<InternalImageType,
> InternalImageType> MatchingFilterType;
>     MatchingFilterType::Pointer matchFilter = MatchingFilterType::New();
>     matchFilter->SetInput(movingCaster->GetOutput());
>     matchFilter->SetReferenceImage(fixedCaster->GetOutput());
>     matchFilter->SetNumberOfHistogramLevels(32);
>     matchFilter->SetNumberOfMatchPoints(100);
>     matchFilter->ThresholdAtMeanIntensityOn();
>     matchFilter->Update();
> 
>     // Register
>     typedef float VectorPixelType;
>     typedef itk::Vector<VectorPixelType, dimension> VectorType;
>     typedef itk::Image<VectorType, dimension> DeformationFieldType;
>     typedef itk::MultiResolutionPDEDeformableRegistration<InternalImageType,
> InternalImageType, DeformationFieldType>
> 
> RegistrationFilterType;
>     RegistrationFilterType::Pointer registration =
> RegistrationFilterType::New();
>     registration->SetFixedImage(fixedCaster->GetOutput());
>     registration->SetMovingImage(matchFilter->GetOutput());
>     unsigned int numLevel = 5;
>     unsigned int numIterations[10];
>     numIterations[0] = 128;
>     for (unsigned int iLevel = 1; iLevel < numLevel; iLevel++)
>         numIterations[iLevel] = numIterations[iLevel-1] / 2;
>     registration->SetNumberOfLevels(numLevel);
>     registration->SetNumberOfIterations(numIterations);
>     registration->Print(cout);
>     registration->Update();
> 
>     // Apply result (Warpping)
>     typedef itk::WarpImageFilter<ImageType, ImageType, DeformationFieldType>
> WarpFilterType;
>     WarpFilterType::Pointer warper = WarpFilterType::New();
>     typedef WarpFilterType::CoordRepType CoordRepType;
>     typedef itk::NearestNeighborInterpolateImageFunction<ImageType,
> CoordRepType> InterpolatorType;
>     InterpolatorType::Pointer interpolator = InterpolatorType::New();
>     warper->SetInput(movingImg);
>     warper->SetInterpolator(interpolator);
>     warper->SetOutputSpacing(matchFilter->GetOutput()->GetSpacing());
>     warper->SetOutputOrigin(matchFilter->GetOutput()->GetOrigin());
>     warper->SetDeformationField(registration->GetOutput());
>     warper->Update();
>     return warper->GetOutput();
> }
> 
> Thank you very much Luis!!!
> 
> Cheers,
> Andy
> 
> 
> 
>>----- Original Message ----- 
>>From: "Luis Ibanez" <luis.ibanez at kitware.com>
>>To: "Andy Eow" <vtkitk at hotmail.com>
>>Cc: <Insight-users at itk.org>
>>Sent: Tuesday, March 09, 2004 4:55 PM
>>Subject: Re: Non-rigid Registration
>>
>>
>>
>>>Hi Andy,
>>>
>>>ITK doesn't have implemented a straightforward
>>>method for extracting the Affine part of a
>>>deformation field.
>>>
>>>The deformation field is much richer in information
>>>than a simple Affine transform.
>>>
>>>You could compute the components of the Affine
>>>transform by using the KernelTransform.
>>>http://www.itk.org/Insight/Doxygen/html/classitk_1_1KernelTransform.html
>>>
>>>Simply take a number of samples from the deformation
>>>field and generate with them source and targe
>>>landmarks. Feed those in any of the subclasses
>>>of the KernelTransform and compute its matrix.
>>>
>>>This transform has the nice property of computing
>>>the Affine part of the desformation.
>>>
>>>
>>>   For details you can look at the IEEE TMI paper by
>>>   Davis, Khotanzad, Flamig, and Harms,
>>>   Vol. 16, No. 3 June 1997.
>>>
>>>
>>>Now, note that it is an overkill to use a deformable
>>>registration method if at the end the only thing that
>>>you want to use is the Affine part. You will be computing
>>>hundreds/thousands of parameters, just to reduce them later
>>>to only 12. It will be much simpler to use a basic registration
>>>algorithm with an AffineTransform...
>>>
>>>
>>>--------------
>>>
>>>
>>>About using MutualInformation for non-rigid
>>>registration, You can do this right now in ITK.
>>>
>>>One option is to use BSpline transform
>>>
>>
> http://www.itk.org/Insight/Doxygen/html/classitk_1_1BSplineDeformableTransform.html
> 
>>>along with any of the MutualInformation metric classes
>>>available in ITK (Viola Wells, Mattes or Histogram).
>>>
>>>
>>>and the other is to use the FEMRegistration filter
>>>and in the configuration file select Mutual Information
>>>as metric.
>>>
>>>
>>>
>>>
>>>   Regards,
>>>
>>>
>>>
>>>        Luis
>>>
>>>
>>>
>>>-----------------
>>>Andy Eow wrote:
>>>
>>>
>>>>Hi Luis,
>>>>
>>>>I'm currently performing a non-rigid registration using
>>>>DemonsRegistrationFilter to obtain a DeformationField. However, I'd
> 
> like
> 
>>>>to recover the affine rotation matrix representation of the
> 
> deformation,
> 
>>>>modify this matrix and reapply the affine matrix to the target image.
> 
> Is
> 
>>>>this possible within itk? Are there any examples that demonstrate
> 
> this?
> 
>>>>Also, does itk currently have a non-rigid registration filter based on
> 
> a
> 
>>>>mutual information metric?
>>>>
>>>>Thank you very much.
>>>>
>>>>Cheers,
>>>>Andy
>>>
>>>
>>>
>>>
>