[Insight-users] Can not reproduce metric value

[Bill Lorensen]

The difference could be that the internal values have higher precision
than the printed values.

You should be getting an exception with the try/catch.


On Sun, Oct 4, 2009 at 1:51 PM, Quy Pham Sy <phamsyquybk at gmail.com> wrote:
> hi,  Bill
> Thanks very much for you help. I tried adding your code. Actually your code
> outputs the metric value to screen, am i right?
> To make all my question clearer
> here is the summary, there're issues:
>
> 1 - When I used the output (final transform parameters) from
> ImageRegistration7 example (in ITK source code)
> to reproduce the metric  value directly from Metric Object (same as metric
> object type used in ImageRegistration7 example)
> I can not reproduce the same metric value as the example. Why?
>
> Output from ImageRegistration7
> -------------------------------------
>  Scale         = 1.08438
>  Angle (radians) 0.154531
>  Angle (degrees) 8.85397
>  Center X      = 170.327
>  Center Y      = 148.556
>  Translation X = -26.6822
>  Translation Y = -93.4428
>  Iterations    = 500
>  Metric value  = 4551.86
> ---------------------------------------
> Two input image: fixed.bmp; moving.bmp
> Set params code for transformation of metric
> ---------------------------------------------------
> transform->SetScale(1.08438);
> transform->SetAngle(0.154531);
> TransformType::TranslationType translation;
> translation[0] = -26.6822;
> translation[1] = -93.6822;
> transform->SetTranslation(translation);
> MetricType::TransformParametersType transformParams =
> transform->GetParameters();
>
> double metricval = metric->GetValue(transformParams);
> /*result: metricval = 4570.285 still != 4551.86*/
> ----------------------------------------------------
>
> 2. In theory,  the metric value from registration (determined by: transform
> type, metric type, interpolate type) of two images (fixed image, moving
> image) is supposed to be "MINIMUM METRIC VALUE".
> But when i try a set of params (as following) i get a smaller value, event
> down to 0 (still use two input image as above) the value which is impossible
> since two images are not identical !!!
> ----------------------------------------------
> transform->SetScale(0.92632);
> transform->SetAngle(0.15095);
> TransformType::TranslationType translation;
> translation[0] = 28.29828;
> translation[1] = 94.07266;
> transform->SetTranslation(translation);
> MetricType::TransformParametersType transformParams =
> transform->GetParameters();
>
> double metricval = metric->GetValue(transformParams);
> /*result: metricval == 0 !!!! */
> -------------------------------------------------
> Really appreciate any suggest and help!
> //========================================================================
> //              HERE IS THE SOURCE CODE TO REPRODUCE THESE PROBLEMS
> //========================================================================
> int main( int argc, char * argv[] )
> {
> const     unsigned int   Dimension = 2;
> typedef   unsigned char  PixelType;
> typedef itk::Image< PixelType, Dimension >   ImageType;
> typedef itk::ImageFileReader< ImageType >  ReaderType;
> ReaderType::Pointer fixedReader  = ReaderType::New();
> ReaderType::Pointer movingReader = ReaderType::New();
>   fixedReader->SetFileName(  "fixed.bmp" );
>   movingReader->SetFileName( "moving.bmp" );
>
> try
> {
> fixedReader->Update();
> movingReader->Update();
> }
> catch( itk::ExceptionObject & excep )
> {
> std::cerr << "Exception catched !" << std::endl;
> std::cerr << excep << std::endl;
> }
>
> typedef itk::MeanSquaresImageToImageMetric< ImageType, ImageType >
>  MetricType;
> MetricType::Pointer metric = MetricType::New();
> typedef itk::CenteredSimilarity2DTransform < double >  TransformType;
> TransformType::Pointer transform = TransformType::New();
> typedef itk::LinearInterpolateImageFunction< ImageType, double >
>  InterpolatorType;
> InterpolatorType::Pointer interpolator = InterpolatorType::New();
> ImageType::Pointer fixedImage  = fixedReader->GetOutput();
> ImageType::Pointer movingImage = movingReader->GetOutput();
> typedef itk::CenteredTransformInitializer< TransformType, ImageType,
> ImageType >  TransformInitializerType;
> TransformInitializerType::Pointer initializer =
> TransformInitializerType::New();
> initializer->SetTransform(   transform );
> initializer->SetFixedImage(  fixedImage );
> initializer->SetMovingImage( movingImage );
> initializer->MomentsOn();
> initializer->InitializeTransform();
> metric->SetTransform( transform );
> metric->SetInterpolator( interpolator );
> metric->SetFixedImage(  fixedImage  );
> metric->SetMovingImage( movingImage );
> metric->SetFixedImageRegion(  fixedImage->GetBufferedRegion()  );
> try
> {
> metric->Initialize();
> }
> catch( itk::ExceptionObject & excep )
> {
> std::cerr << "Exception catched !" << std::endl;
> std::cerr << excep << std::endl;
> return EXIT_FAILURE;
> }
> //************************************************************************/
>         //  CHANGE PARAMETERS HERE TO REPRODUCE PROBLEM
> transform->SetScale(0.92632);
> transform->SetAngle(0.15095);
> TransformType::TranslationType translation;
> translation[0] = 28.29828;
> translation[1] = 94.07266;
>
> //*************************************************************************/
>
> transform->SetTranslation(translation);
> MetricType::TransformParametersType transformParams =
> transform->GetParameters();
> double metricval = 0.0;
> try
> {
> metricval = metric->GetValue(transformParams);
> std::cout << "metricval: " << metricval << std::endl;
> }
> catch( itk::ExceptionObject & excep )
> {
> std::cerr << "Exception catched !" << std::endl;
> std::cerr << excep << std::endl;
> std::cout << "metricval: " << metricval << std::endl;
> return EXIT_FAILURE;
> }
> return EXIT_SUCCESS;
> }
> =======================================================================================
>