[Insight-users] About PointSetToImageRegistrationMethod

[xiaojun chen]

Hi,all,
    I need to register a CAD model to CT image of a patient with this
model. First of all, I
got an initial transform matrix with fiducial registration method, which
already puts the CAD model in approximately the right place, and then I
obtained the gradient magnitude image of the CT data, and tried using
PointSetToImageRegistrationMethod to register the
point set to it.
    But when registration is performed, only rotation parameters are
changed during registration, and the fractional translation around zero is
observed. I tried changing translation scale maximum and minimum step
length but it is not translating at all and only rotation is effected by
change of these parameters, and the registration result is totally wrong.
    Could anyone please help me take a look at the source codes and give me
some suggestions?Thanks!

P.S. the source codes are as follows:

void registerCADmodelToImage()
{
       typedef itk::Image<double, 3> ITKImageType;
       typedef
itk::PointSetToImageRegistrationMethod<PointSetType,ITKImageType>
RegistrationType;
       typedef itk::PointSet< double, 3 > PointSetType;
       typedef itk::VersorRigid3DTransform< double > TransformType;
       typedef itk::VersorRigid3DTransformOptimizer OptimizerType;
       typedef
itk::NormalizedCorrelationPointSetToImageMetric<PointSetType,ITKImageType>
MetricType;
       typedef itk:: LinearInterpolateImageFunction<ITKImageType,double >
InterpolatorType;
       typedef MetricType::TransformType TransformBaseType;
       typedef TransformBaseType::ParametersType ParametersType;

       RegistrationType::Pointer RegistrationMethod=
RegistrationType::New();
       PointSetType::Pointer pointSet= PointSetType::New();
       MetricType::Pointer metric = MetricType::New();
       OptimizerType::Pointer optimizer = OptimizerType::New();
       InterpolatorType::Pointer interpolator = InterpolatorType::New();
       TransformType::Pointer transform = TransformType::New();
       transform->SetIdentity();

//--------------------------------------------------------------------------------------------------
// Set the initial registration matrix for the CAD model
//--------------------------------------------------------------------------------------------------
       itk::Matrix<double, 3, 3> transformMatrix,transformMatrixIni;
       itk::Vector<double, 3>    vec,vecIni;
       static const int D = 3;
       int i, j;

       transformMatrix[0][0]=0.1675;
transformMatrix[0][1]=-0.985872;transformMatrix[0][2]=0;
       transformMatrix[1][0]=-0.985872;
transformMatrix[1][1]=-0.1675;transformMatrix[1][2]=0;
       transformMatrix[2][0]=0;
transformMatrix[2][1]=0;transformMatrix[2][2]=-1;
       vec[0]=8.9;vec[1]=41.47;vec[2]=-103.88;

       vtkSmartPointer<vtkMatrix4x4> vtkmatInitial =
vtkSmartPointer<vtkMatrix4x4>::New();
       vtkmatInitial->Identity();

       for (i=0; i < D; i++)
               {
               for (j=0; j < D; j++)
                       {
                       (*vtkmatInitial)[i][j] = transformMatrix[i][j];
                       }
               (*vtkmatInitial)[i][D] = vec[i];
               }

//--------------------------------------------------------------------------------------------------
// Set the point data set for the CAD model
//--------------------------------------------------------------------------------------------------
       unsigned int pointId=0;
       setPointData(35,34,vtkmatInitial,&pointId,pointSet);
       setPointData(25,36.679,vtkmatInitial,&pointId,pointSet);

//----------------------------------------------------------------------------
// Load the moving image
//----------------------------------------------------------------------------
 typedef itk::ImageFileReader<ITKImageType> MovingImageReaderType;

 MovingImageReaderType::Pointer movingReader =
MovingImageReaderType::New();
 movingReader->SetFileName( "1.mhd" );

 try
   {
   movingReader->Update();
   }
 catch( itk::ExceptionObject & e )
   {
   std::cout << e.GetDescription() << std::endl;
   return ;
   }

 typedef itk::GradientMagnitudeImageFilter<
                 ITKImageType, ITKImageType >  filterType;

  // Create and setup a gradient filter
 filterType::Pointer gradientFilter = filterType::New();
 gradientFilter->SetInput( movingReader->GetOutput() );
 gradientFilter->Update();


 RegistrationMethod->SetInitialTransformParameters(transform->GetParameters()
);
       RegistrationMethod->SetTransform( transform );

       // Scale the translation components of the Transform in the Optimizer
       typedef OptimizerType::ScalesType OptimizerScalesType;
       OptimizerScalesType optimizerScales(
transform->GetNumberOfParameters());

       const double translationScale = 1.0 ;
       optimizerScales[0] = 0.1;
       optimizerScales[1] = 0.1;
       optimizerScales[2] = 0.1;
       optimizerScales[3] = translationScale;
       optimizerScales[4] = translationScale;
       optimizerScales[5] = translationScale;
       optimizer->SetMaximumStepLength( 0.1 );
       optimizer->SetMinimumStepLength( 0.01 );
       optimizer->SetNumberOfIterations( 500 );
       optimizer->SetRelaxationFactor( 0.90 );
       optimizer->SetGradientMagnitudeTolerance( 0.005 );
       optimizer->MinimizeOn();
       RegistrationMethod->SetFixedPointSet( pointSet );
       RegistrationMethod->SetMetric( metric );
       RegistrationMethod->SetOptimizer( optimizer );
       RegistrationMethod->SetInterpolator( interpolator );
       RegistrationMethod->SetMovingImage( gradientFilter->GetOutput() );

       try
               {
               RegistrationMethod->Update();
               }

       catch( itk::ExceptionObject & e )
               {
               std::cout << "testing!!!" << std::endl;
               std::cerr << e << std::endl;
               ostrstream ostr;
               cout.rdbuf(ostr.rdbuf());
               std::cout << e << std::endl;
               }

       OptimizerType::ParametersType
finalParameters=RegistrationMethod->GetLastTransformParameters();

       double versorX              = finalParameters[0];
       double versorY              = finalParameters[1];
       double versorZ              = finalParameters[2];
       double finalTranslationX    = finalParameters[3];
       double finalTranslationY    = finalParameters[4];
       double finalTranslationZ    = finalParameters[5];

       const unsigned int numberOfIterations =
optimizer->GetCurrentIteration();
       const double bestValue = optimizer->GetValue();

 // Print out results
 //
 std::cout << std::endl << std::endl;
 std::cout << "Result = " << std::endl;
 std::cout << " versor X      = " << versorX  << std::endl;
 std::cout << " versor Y      = " << versorY  << std::endl;
 std::cout << " versor Z      = " << versorZ  << std::endl;
 std::cout << " Translation X = " << finalTranslationX  << std::endl;
 std::cout << " Translation Y = " << finalTranslationY  << std::endl;
 std::cout << " Translation Z = " << finalTranslationZ  << std::endl;
 std::cout << " Iterations    = " << numberOfIterations << std::endl;
 std::cout << " Metric value  = " << bestValue          << std::endl;
}
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20120430/a0c8c998/attachment.htm>