[ITK-users] Fwd: Metric Evaluation

[Matt McCormick]

Hi Suneeza,

As show here:

  http://www.itk.org/Doxygen/html/Examples_2RegistrationITKv4_2MeanSquaresImageMetric1_8cxx-example.html

Does calling

  metric->SetVirtualDomainFromImage( fixedImageReader->GetOutput() );

improve the result?

HTH,
Matt

On Thu, Jan 7, 2016 at 4:52 AM, suneeza hanif <suneezahanif02 at gmail.com> wrote:
> Hello Matt,
>
> Here is the full version of code. Can you please guide me where am I doing
> something wrong because I am unable to get it. Even If I stop after one
> iteration the values are very different.
>
> I would be highly thankful to you.
>
> const    unsigned int    Dimension = 2;
> typedef  float          PixelType;
>
> typedef itk::Image< PixelType, Dimension >  FixedImageType;
> typedef itk::Image< PixelType, Dimension >  MovingImageType;
>
> typedef itk::TranslationTransform< double, Dimension > TransformType;
>
> typedef itk::RegularStepGradientDescentOptimizerv4<double> OptimizerType;
>
> typedef itk::MeanSquaresImageToImageMetricv4< FixedImageType,MovingImageType
>> MetricType;
>
> typedef itk::ImageRegistrationMethodv4< FixedImageType, MovingImageType,
> TransformType > RegistrationType;
>
> MetricType::Pointer         metric        = MetricType::New();
> OptimizerType::Pointer      optimizer     = OptimizerType::New();
> RegistrationType::Pointer   registration  = RegistrationType::New();
>
> registration->SetMetric(        metric        );
> registration->SetOptimizer(     optimizer     );
>
> typedef itk::LinearInterpolateImageFunction< FixedImageType, double >
> FixedLinearInterpolatorType;
>
> typedef itk::LinearInterpolateImageFunction< MovingImageType, double >
> MovingLinearInterpolatorType;
>
> FixedLinearInterpolatorType::Pointer fixedInterpolator =
> FixedLinearInterpolatorType::New();
> MovingLinearInterpolatorType::Pointer movingInterpolator =
> MovingLinearInterpolatorType::New();
>
> metric->SetFixedInterpolator(  fixedInterpolator  );
> metric->SetMovingInterpolator(  movingInterpolator  );
>
> typedef itk::ImageFileReader< FixedImageType  >   FixedImageReaderType;
> typedef itk::ImageFileReader< MovingImageType >   MovingImageReaderType;
>
> FixedImageReaderType::Pointer   fixedImageReader     =
> FixedImageReaderType::New();
> MovingImageReaderType::Pointer  movingImageReader    =
> MovingImageReaderType::New();
>
> fixedImageReader->SetFileName(  input1 );
> movingImageReader->SetFileName( input2 );
>
> registration->SetFixedImage( fixedImageReader->GetOutput() );
> registration->SetMovingImage( movingImageReader->GetOutput() );
> typedef float                    OutputPixelType;
>
> typedef itk::Image< PixelType, Dimension > OutputImageType;
> TransformType::Pointer transform = TransformType::New();
>
> transform->SetIdentity();
> registration->SetInitialTransform( transform );
>
> optimizer->SetLearningRate( 4 );
> optimizer->SetMinimumStepLength( 0.001 );
> optimizer->SetRelaxationFactor( 0.5 );
>
> std::cout<<"Transform before registration: "<<transform<<endl;
> RegistrationType::MetricSamplingStrategyType  samplingStrategy  =
> RegistrationType::RANDOM;
> double samplingPercentage = 0.20;
> registration->SetMetricSamplingStrategy( samplingStrategy );
> registration->SetMetricSamplingPercentage( samplingPercentage );
>
> optimizer->SetNumberOfIterations( 1 );
>
> CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
> optimizer->AddObserver( itk::IterationEvent(), observer );
>
> const unsigned int numberOfLevels = 1;
> registration->SetNumberOfLevels ( numberOfLevels );
> try
> {
> std::cout<<"Evaluating registration parameters..."<<std::endl;
> registration->Update();
> std::cout << "Optimizer stop condition: ";
> std::cout<< registration->GetOptimizer()->GetStopConditionDescription();
> std::cout << std::endl;
> }
> catch( itk::ExceptionObject & err )
> {
> std::cerr << "ExceptionObject caught !" << std::endl;
> std::cerr << err << std::endl;
> system("pause");
> return EXIT_FAILURE;
> }
>
> TransformType::ParametersType finalParameters =
> registration->GetTransform()->GetParameters();
> const double TranslationAlongX = finalParameters[0];
> const double TranslationAlongY = finalParameters[1];
>
> const unsigned int numberOfIterations = optimizer->GetCurrentIteration();
>
> const double bestValue = optimizer->GetValue();
>
> std::cout << "Result = " << std::endl;
> std::cout << " Translation X = " << TranslationAlongX  << std::endl;
> std::cout << " Translation Y = " << TranslationAlongY  << std::endl;
> std::cout << " Iterations    = " << numberOfIterations << std::endl;
> std::cout << " Metric value  = " << bestValue          << std::endl;
>
>
>
> typedef itk::ResampleImageFilter< MovingImageType, FixedImageType >
> ResampleFilterType;
>
> ResampleFilterType::Pointer resampler = ResampleFilterType::New();
> resampler->SetInput( movingImageReader->GetOutput() );
> resampler->SetTransform( registration->GetTransform() );
> FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();
>
> resampler->SetSize( fixedImage->GetLargestPossibleRegion().GetSize() );
> resampler->SetOutputOrigin(  fixedImage->GetOrigin() );
> resampler->SetOutputSpacing( fixedImage->GetSpacing() );
> resampler->SetOutputDirection( fixedImage->GetDirection() );
> resampler->SetDefaultPixelValue( 0 );
>
>
> ////////////////////////////////////////////////////////////////////////////////////////////
>
> std::cout<<"Transform after registration: "<<transform<<endl;
>
> MetricType::Pointer metric2 = MetricType::New();
>
> TransformType::Pointer transformCopy = transform->Clone();
> std::cout << "Transform copy: " << transformCopy << std::endl;
>
> metric2->SetFixedImage( fixedImageReader->GetOutput() );
> metric2->SetMovingImage(movingImageReader->GetOutput() );
> metric2->SetFixedInterpolator(  fixedInterpolator  );
> metric2->SetMovingInterpolator(  movingInterpolator  );
> metric2->SetTransform(transformCopy);
>
> try
> {
> metric2->Initialize();
> }
> catch( itk::ExceptionObject & err )
> {
> std::cerr << "ExceptionObject caught !" << std::endl;
> std::cerr << err << std::endl;
> }
>
> cout<<endl<<"Metric Value (Mean Square Image to Image):
> "<<metric2->GetValue()<<endl;
>
> ////////////////////////////////////////////////////////////////////////////////////////////
>
> Best Regards,
> Suneeza
>
> On Wed, Jan 6, 2016 at 11:53 PM, suneeza hanif <suneezahanif02 at gmail.com>
> wrote:
>>
>> Hello Matt,
>>
>> Thanks a lot for your help.
>>
>> I am still unable to get same metric even after using this.
>>
>> I am using Centered2DTransform and during registration I have not
>> initialized it so it will take (0,0) as center of rotation. Then I defined
>> FixedParameters again with (0,0) but metric value is still quite different.
>>
>> Best Regards,
>> Suneeza
>>
>