itksimilarity3dtransformtest.cxx

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

    ParametersType parameters( np ); // Number of parameters

    VersorType versor;

    parameters[0] = versor.GetX();   // Rotation axis * sin(t/2)
    parameters[1] = versor.GetY();
    parameters[2] = versor.GetZ();
    parameters[3] = 8.0;             // Translation
    parameters[4] = 7.0;
    parameters[5] = 6.0;
    parameters[6] = 1.0;             // Scale

    transform->SetParameters( parameters );

    ParametersType parameters2 = transform->GetParameters();

    const double tolerance = 1e-8;

    for(unsigned int p=0; p<np; p++)
      {
      if( fabs( parameters[p] - parameters2[p] ) > tolerance )
        {
        std::cout << "Output parameter does not match input " << std::endl;
        return EXIT_FAILURE;
        }
      }
     std::cout << "Input/Output parameter check Passed !"  << std::endl;

     // Try the GetJacobian method
     TransformType::InputPointType  aPoint;
     aPoint[0] = 10.0;
     aPoint[1] = 20.0;
     aPoint[2] = -10.0;
     JacobianType   jacobian = transform->GetJacobian( aPoint );
     std::cout << "Jacobian: "  << std::endl;
     std::cout << jacobian << std::endl;

     // copy the read one just for getting the right matrix size
     JacobianType   TheoreticalJacobian = jacobian;
     
     TheoreticalJacobian[0][0] =    0.0;
     TheoreticalJacobian[1][0] =  206.0;
     TheoreticalJacobian[2][0] =  -84.0;

     TheoreticalJacobian[0][1] = -206.0;
     TheoreticalJacobian[1][1] =    0.0;
     TheoreticalJacobian[2][1] =   42.0;

     TheoreticalJacobian[0][2] =   84.0;
     TheoreticalJacobian[1][2] =  -42.0;
     TheoreticalJacobian[2][2] =    0.0;

     TheoreticalJacobian[0][3] = 1.0;
     TheoreticalJacobian[1][3] = 0.0;
     TheoreticalJacobian[2][3] = 0.0;

     TheoreticalJacobian[0][4] = 0.0;
     TheoreticalJacobian[1][4] = 1.0;
     TheoreticalJacobian[2][4] = 0.0;

     TheoreticalJacobian[0][5] = 0.0;
     TheoreticalJacobian[1][5] = 0.0;
     TheoreticalJacobian[2][5] = 1.0;

     TheoreticalJacobian[0][6] =  -21.0;
     TheoreticalJacobian[1][6] =  -42.0;
     TheoreticalJacobian[2][6] = -103.0;

     for(unsigned int ii=0; ii < 3; ii++)
       {
       for(unsigned int jj=0; jj < 7; jj++)
         {
         if( vnl_math_abs( TheoreticalJacobian[ii][jj] - jacobian[ii][jj] ) > 1e-5 )
           {
           std::cout << "Jacobian components differ from expected values ";
           std::cout << std::endl << std::endl;
           std::cout << "Expected Jacobian = " << std::endl;
           std::cout << TheoreticalJacobian << std::endl << std::endl;
           std::cout << "Computed Jacobian = " << std::endl;
           std::cout << jacobian << std::endl << std::endl;
           std::cout << std::endl << "Test FAILED ! " << std::endl;
           return EXIT_FAILURE;
           }
         }
       }
  }

  {
  std::cout << " Exercise the SetIdentity() method " << std::endl; 
  TransformType::Pointer  transform = TransformType::New();

  itk::Vector<double,3> axis(1);

  const double angle = (atan(1.0)/45.0)*30.0; // turn 30 degrees

  transform->SetRotation( axis, angle );

  TransformType::InputPointType  center;
  center[0] = 31;
  center[1] = 62;
  center[2] = 93;
  
  transform->SetCenter( center );

  transform->SetIdentity();

  const unsigned int np = transform->GetNumberOfParameters();

  ParametersType parameters( np ); // Number of parameters

  VersorType versor;

  parameters[0] = versor.GetX();   // Rotation axis * sin(t/2)
  parameters[1] = versor.GetY();
  parameters[2] = versor.GetZ();
  parameters[3] = 0.0;             // Translation
  parameters[4] = 0.0;
  parameters[5] = 0.0;
  parameters[6] = 1.0;


  ParametersType parameters2 = transform->GetParameters();

  const double tolerance = 1e-8;

  for(unsigned int p=0; p<np; p++)
    {
    if( fabs( parameters[p] - parameters2[p] ) > tolerance )
      {
      std::cout << "Output parameter does not match input " << std::endl;
      return EXIT_FAILURE;
      }
    }
  std::cout << "Input/Output parameter check Passed !"  << std::endl;
  }

  {
  std::cout << " Exercise the Scaling methods " << std::endl; 
  TransformType::Pointer  transform = TransformType::New();

  itk::Vector<double,3> axis(1);

  const double angle = (atan(1.0)/45.0)*30.0; // turn 30 degrees

  transform->SetRotation( axis, angle );

  TransformType::InputPointType  center;
  center[0] = 31;
  center[1] = 62;
  center[2] = 93;
  
  transform->SetCenter( center );

  TransformType::OutputVectorType translation;
  translation[0] = 17;
  translation[1] = 19;
  translation[2] = 23;

  transform->SetTranslation( translation );


  const double scale = 2.5;

  transform->SetScale( scale );

  const double rscale = transform->GetScale();

  const double tolerance = 1e-8;

  if( fabs( rscale - scale ) > tolerance )
    {
    std::cout << "Error in Set/Get Scale() " << std::endl;
    return EXIT_FAILURE;
    }

  const unsigned int np = transform->GetNumberOfParameters();

  ParametersType parameters( np ); // Number of parameters

  VersorType versor;
  versor.Set( axis, angle );

  parameters[0] = versor.GetX();   // Rotation axis * sin(t/2)
  parameters[1] = versor.GetY();
  parameters[2] = versor.GetZ();
  parameters[3] = translation[0];
  parameters[4] = translation[1];
  parameters[5] = translation[2];
  parameters[6] = scale;


  ParametersType parameters2 = transform->GetParameters();

  for(unsigned int p=0; p<np; p++)
    {
    if( fabs( parameters[p] - parameters2[p] ) > tolerance )
      {
      std::cout << "Output parameter does not match input " << std::endl;
      return EXIT_FAILURE;
      }
    }
  std::cout << "Input/Output parameter check Passed !"  << std::endl;
  }


  {
     // Testing SetMatrix()
     std::cout << "Testing SetMatrix() ... ";
     bool Ok;
     unsigned int par;

     typedef TransformType::MatrixType MatrixType;
     MatrixType matrix;

     TransformType::Pointer t = TransformType::New();
      
     // attempt to set an non-orthogonal matrix
     par = 0;
     for( unsigned int row = 0; row < 3; row++ )
        {
        for( unsigned int col = 0; col < 3; col++ )
          {
          matrix[row][col] = static_cast<double>( par + 1 );
          ++par;
          }
        }

     Ok = false;

     std::cout << "Setting non-orthogonal matrix = " << std::endl;
     std::cout << matrix << std::endl;

     try
      {
      t->SetMatrix( matrix );
      }
     catch ( itk::ExceptionObject & err )
      {
      std::cout << "Caught expected exception" << err << std::endl;
      Ok = true;
      }
     catch( ... )
      {
      std::cout << "Caught unknown exception" << std::endl;
      }

     if( !Ok )
      {
      std::cout << "Error: expected to catch an exception when attempting";
      std::cout << " to set an non-orthogonal matrix." << std::endl;
      return EXIT_FAILURE;
      }

      t = TransformType::New();

      // attempt to set an (orthogonal + scale) matrix
      matrix.GetVnlMatrix().set_identity();

      double a = 1.0 / 180.0 * vnl_math::pi;
      double s = 0.5;
      matrix[0][0] =        cos( a ) * s;
      matrix[0][1] = -1.0 * sin( a ) * s;
      matrix[1][0] =        sin( a ) * s; 
      matrix[1][1] =        cos( a ) * s;
      matrix[2][2] =                   s;

     Ok = true;
     try
      {
      t->SetMatrix( matrix );
      }
     catch ( itk::ExceptionObject & err )
      {
      std::cout << err << std::endl;
      Ok = false;
      }
     catch( ... )
      {
      std::cout << "Caught unknown exception" << std::endl;
      Ok = false;
      }

     if( !Ok )
      {
      std::cout << "Error: caught unexpected exception" << std::endl;
      return EXIT_FAILURE;
      }

   // Check the computed parameters
    TransformType::InputPointType center;
    center[0] = 15.0;
    center[1] = 16.0;
    center[2] = 17.0;

    typedef TransformType::ParametersType ParametersType;
    ParametersType e( t->GetNumberOfParameters() );
    e.Fill( 0.0 );
    e[2] = sin(0.5 * a);
    e[6] = 0.5;

    t = TransformType::New();
    t->SetCenter( center );
    t->SetParameters( e );

    TransformType::Pointer t2 = TransformType::New();
    t2->SetCenter( center );
    t2->SetMatrix( t->GetMatrix() );

    ParametersType p = t2->GetParameters();

    for( unsigned int k = 0; k < e.GetSize(); k++ )
      {
      if( fabs( e[k] - p[k] ) > epsilon )
        {
        std::cout << " [ FAILED ] " << std::endl;
        std::cout << "Expected parameters: " << e << std::endl;
        std::cout << "but got: " << p << std::endl;
        return EXIT_FAILURE; 
        }
      }

    std::cout << "[ PASSED ]" << std::endl;
    }

  std::cout << std::endl << "Test PASSED ! " << std::endl;

  return EXIT_SUCCESS;

}