itksimilarity2dtransformtest.cxx

DTMK软件开发包,此为开源软件,是一款很好的医学图像开发资源.

/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkSimilarity2DTransformTest.cxx,v $
  Language:  C++
  Date:      $Date: 2008-06-29 13:25:41 $
  Version:   $Revision: 1.17 $

  Copyright (c) Insight Software Consortium. All rights reserved.
  See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even 
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
     PURPOSE.  See the above copyright notices for more information.

=========================================================================*/
#if defined(_MSC_VER)
#pragma warning ( disable : 4786 )
#endif

#include <iostream>

#include "itkCenteredSimilarity2DTransform.h"
#include "vnl/vnl_vector_fixed.h"
#include "itkVector.h"

namespace 
{
bool CheckEqual( 
 itk::Point<double,2> p1,
 itk::Point<double,2> p2 )
{
  const double epsilon = 1e-10;
  for( unsigned int i = 0; i < 2; i++ )
    {
    if( fabs( p1[i] - p2[i] ) > epsilon )
      {
      std::cout << p1 << " != " << p2 << ":[ FAILED ]" << std::endl;
      return false;
      }
    }
  std::cout << p1 << " == " << p2 << ":[ PASSED ]" << std::endl;
  return true;
}
}


int itkSimilarity2DTransformTest(int ,char *[] )
{

  std::cout << "==================================" << std::endl;
  std::cout << "Testing Similarity 2D Transform" << std::endl << std::endl;

  const double epsilon = 1e-10;
  const unsigned int N = 2;
  bool Ok = true;

  typedef itk::CenteredSimilarity2DTransform<double>  SimilarityTransformType;
  SimilarityTransformType::Pointer transform = SimilarityTransformType::New();

  // Test the identity transform
  std::cout << "Testing Identity:";
  transform->SetIdentity();
 
  SimilarityTransformType::InputPointType::ValueType pInit[2] = {10,10};
  SimilarityTransformType::InputPointType p = pInit;
  SimilarityTransformType::OutputPointType r;


  r = transform->TransformPoint( p );
  for(unsigned int i=0; i<N; i++)
    {
    if( fabs( p[i]- r[i] ) > epsilon )
      {
      Ok = false;
      break;    
      }
    }
  if( !Ok )
    { 
    std::cerr << "Error with Identity transform" << std::endl;
    return EXIT_FAILURE;
    }
  else
    {
    std::cout << " [ PASSED ] " << std::endl;
    }

  // Test SetAngle/GetAngle
  SimilarityTransformType::Pointer transform1 = SimilarityTransformType::New();
  SimilarityTransformType::Pointer transform2 = SimilarityTransformType::New();
  transform1->SetIdentity();
  double angle1 = .125;
  transform1->SetAngle(angle1);
  transform2->SetMatrix(transform1->GetMatrix());
  std::cout << "Testing SetAngle(" << angle1 << ")/GetAngle():";
  const double epsilon2 = 1e-5;
  if (fabs(transform2->GetAngle() - angle1) > epsilon2)
    {
    std::cerr << "Error with SetAngle/GetAngle:" << std::endl;
    std::cerr << "transform1->SetAngle: " << angle1 << std::endl;
    std::cerr << "transform2->GetAngle: " << transform2->GetAngle() << std::endl;
    return EXIT_FAILURE;
    }
  else
    {
    std::cout << " [ PASSED ] " << std::endl;
    }

  transform1->SetIdentity();
  transform1->SetAngle(-angle1);
  transform2->SetMatrix(transform1->GetMatrix());
  std::cout << "Testing SetAngle(" << -angle1 << ")/GetAngle():";
  if (fabs(transform2->GetAngle() - (-angle1)) > epsilon2)
    {
    std::cerr << "Error with SetAngle/GetAngle:" << std::endl;
    std::cerr << "transform1->SetAngle: " << -angle1 << std::endl;
    std::cerr << "transform2->GetAngle: " << transform2->GetAngle() << std::endl;
    return EXIT_FAILURE;
    }
  else
    {
    std::cout << " [ PASSED ] " << std::endl;
    }
  // Test the Set/Get Parameters
  std::cout << "Testing Set/GetParameters():" << std::endl;
  SimilarityTransformType::ParametersType params(6);

  for(unsigned int i=0;i<6;i++)
    {
    params[i]=i+1;
    }

  std::cout << "Input Parameters = " << params << std::endl;

  transform->SetParameters(params);
  SimilarityTransformType::ParametersType outputParams(6);
  
  outputParams = transform->GetParameters();
  
  std::cout << "Output Parameters = " << outputParams << std::endl;
  
  for(unsigned int i=0; i<4; i++) // do not test for the offset
    {
    if( fabs( outputParams[i]-params[i] ) > epsilon )
      {
      Ok = false;
      break;    
      }
    }
  if( !Ok )
    { 
    std::cerr << "Error with Set/GetParameters:" << std::endl;
    std::cerr << "Input:" << params << std::endl;
    std::cerr << "Output:" << outputParams << std::endl;
    return EXIT_FAILURE;
    }
  else
    {
    std::cout << " [ PASSED ] " << std::endl;
    }
  
  // 15 degrees in radians
  transform->SetIdentity();
  const double angle = 15.0 * atan( 1.0f ) / 45.0; 
  const double sinth = sin( angle );
  const double costh = cos( angle );


  std::cout << "Testing Rotation:";
  transform->SetAngle(angle);

  // Rotate an itk::Point
  SimilarityTransformType::InputPointType q;
  p = pInit;
  q[0] =  p[0] * costh - p[1] * sinth;
  q[1] =  p[0] * sinth + p[1] * costh;
  r = transform->TransformPoint( p );
  for(unsigned int i=0; i<N; i++)
    {
    if( fabs( q[i]- r[i] ) > epsilon )
      {
      Ok = false;
      break;    
      }
    }
  if( !Ok )
    { 
    std::cerr << "Error rotating point   : " << p << std::endl;
    std::cerr << "Result should be       : " << q << std::endl;
    std::cerr << "Reported Result is     : " << r << std::endl;
    return EXIT_FAILURE;
    }
  else
    {
    std::cout << " [ PASSED ] " << std::endl;
    }
  
  std::cout << "Testing Translation:";

  transform->SetAngle(0);
  
  SimilarityTransformType::OffsetType::ValueType ioffsetInit[2] = {1,4};
  SimilarityTransformType::OffsetType ioffset = ioffsetInit;

  transform->SetOffset( ioffset );

  q = p + ioffset;
      
  r = transform->TransformPoint( p );
  for(unsigned int i=0; i<N; i++)
    {
    if( fabs( q[i]- r[i] ) > epsilon )
      {
      Ok = false;
      break;    
      }
    }
  if( !Ok )
    { 
    std::cerr << "Error translating point: " << p << std::endl;
    std::cerr << "Result should be       : " << q << std::endl;
    std::cerr << "Reported Result is     : " << r << std::endl;
    return EXIT_FAILURE;
    }
  else
    {
    std::cout << " [ PASSED ] " << std::endl;
    }

  // Testing the Jacobian
  std::cout << "Testing Jacobian:";
  SimilarityTransformType::JacobianType jacobian0 =  transform->GetJacobian(p);

  if( 
    (jacobian0[0][0] != 10) || (jacobian0[0][1] != -10) || (jacobian0[0][2] != 0) || (jacobian0[0][3] != 0) 
    || (jacobian0[0][4] != 1) || (jacobian0[0][5] != 0) || (jacobian0[1][0] != 10) || (jacobian0[1][1] != 10)
    || (jacobian0[1][2] !=0 ) || (jacobian0[1][3] != 0) || (jacobian0[1][4] !=0) || (jacobian0[1][5] != 1)
    ) 
    {
    std::cerr << "Error with Jacobian: " << jacobian0 << std::endl;
    return EXIT_FAILURE;
    }
  else
    {
    std::cout << " [ PASSED ] " << std::endl;
    }

  {
  // Test instantiation, inverse computation, back transform etc.
  typedef itk::Similarity2DTransform<double>  TransformType;
  TransformType::Pointer t1 = TransformType::New();

  // Set parameters
  TransformType::ParametersType parameters( t1->GetNumberOfParameters() );

  parameters[0] = 2.0;
  parameters[1] = -21.0 / 180.0 * vnl_math::pi;
  parameters[2] = 12.0;
  parameters[3] = -8.9;

  TransformType::InputPointType center;
  center[0] = 67.8;
  center[1] = -0.2;

  t1->SetParameters( parameters );
  t1->SetCenter( center );

  TransformType::InputPointType p1;
  p1[0] = 96.8;
  p1[1] = -3.2;

  TransformType::InputPointType p2;
  p2 = t1->TransformPoint( p1 );

  // Test inverse
  TransformType::Pointer t2;
  t1->CloneInverseTo( t2 );

  TransformType::InputPointType p3;
  p3 = t2->TransformPoint( p2 );

  std::cout << "Test CloneInverseTo(): ";
  if( !CheckEqual( p1, p3 ) )
    {
    return EXIT_FAILURE;
    }

  TransformType::Pointer t2dash = TransformType::New();
  t1->GetInverse( t2dash );
  TransformType::InputPointType p3dash;
  p3dash = t2dash->TransformPoint( p2 );

  std::cout << "Test GetInverse(): ";
  if( !CheckEqual( p1, p3dash ) )
    {
    return EXIT_FAILURE;
    }


  // Test clone
  TransformType::Pointer t3;
  t1->CloneTo( t3 );


  TransformType::InputPointType p4;
  p4 = t3->TransformPoint( p1 );