<div class="gmail_quote"><div class="gmail_quote">Hi,all,<br>
I need to register a CAD model to CT image of a patient with this model. First of all, I<br>
got an initial transform matrix with fiducial registration method, which<br>
already puts the CAD model in approximately the right place, and then I<br>
obtained the gradient magnitude image of the
CT data, and tried using PointSetToImageRegistrationMethod to register the<br>
point set to it.<br>
But when registration is performed, only rotation parameters are<br>
changed during registration, and the fractional translation around zero is<br>
observed. I tried changing translation scale maximum and minimum step<br>
length but it is not translating at all and only rotation is effected by<br>
change of these parameters, and the registration result is totally wrong.<br>
Could anyone please help me take a look at the source codes and give me<br>
some suggestions?Thanks!<br>
<br>
P.S. the source codes are as follows:<br>
<br>
void registerCADmodelToImage()<br>
{<br>
typedef itk::Image<double, 3> ITKImageType;<br>
typedef itk::PointSetToImageRegistrationMethod<PointSetType,ITKImageType><br>
RegistrationType;<br>
typedef itk::PointSet< double, 3 > PointSetType;<br>
typedef itk::VersorRigid3DTransform< double > TransformType;<br>
typedef itk::VersorRigid3DTransformOptimizer OptimizerType;<br>
typedef<br>
itk::NormalizedCorrelationPointSetToImageMetric<PointSetType,ITKImageType><br>
MetricType;<br>
typedef itk:: LinearInterpolateImageFunction<ITKImageType,double ><br>
InterpolatorType;<br>
typedef MetricType::TransformType TransformBaseType;<br>
typedef TransformBaseType::ParametersType ParametersType;<br>
<br>
RegistrationType::Pointer RegistrationMethod= RegistrationType::New();<br>
PointSetType::Pointer pointSet= PointSetType::New();<br>
MetricType::Pointer metric = MetricType::New();<br>
OptimizerType::Pointer optimizer = OptimizerType::New();<br>
InterpolatorType::Pointer interpolator = InterpolatorType::New();<br>
TransformType::Pointer transform = TransformType::New();<br>
transform->SetIdentity();<br>
<br>
//--------------------------------------------------------------------------------------------------<br>
// Set the initial registration matrix for the CAD model<br>
//--------------------------------------------------------------------------------------------------<br>
itk::Matrix<double, 3, 3> transformMatrix,transformMatrixIni;<br>
itk::Vector<double, 3> vec,vecIni;<br>
static const int D = 3;<br>
int i, j;<br>
<br>
transformMatrix[0][0]=0.1675;<br>
transformMatrix[0][1]=-0.985872;transformMatrix[0][2]=0;<br>
transformMatrix[1][0]=-0.985872;<br>
transformMatrix[1][1]=-0.1675;transformMatrix[1][2]=0;<br>
transformMatrix[2][0]=0; transformMatrix[2][1]=0;transformMatrix[2][2]=-1;<br>
vec[0]=8.9;vec[1]=41.47;vec[2]=-103.88;<br>
<br>
vtkSmartPointer<vtkMatrix4x4> vtkmatInitial =<br>
vtkSmartPointer<vtkMatrix4x4>::New();<br>
vtkmatInitial->Identity();<br>
<br>
for (i=0; i < D; i++)<br>
{<br>
for (j=0; j < D; j++)<br>
{<br>
(*vtkmatInitial)[i][j] = transformMatrix[i][j];<br>
}<br>
(*vtkmatInitial)[i][D] = vec[i];<br>
}<br>
<br>
//--------------------------------------------------------------------------------------------------<br>
// Set the point data set for the CAD model<br>
//--------------------------------------------------------------------------------------------------<br>
unsigned int pointId=0;<br>
setPointData(35,34,vtkmatInitial,&pointId,pointSet);<br>
setPointData(25,36.679,vtkmatInitial,&pointId,pointSet);<br>
<br><div><div class="h5">
//----------------------------------------------------------------------------<br>
// Load the moving image<br>
//----------------------------------------------------------------------------<br>
typedef itk::ImageFileReader<ITKImageType> MovingImageReaderType;<br>
<br>
MovingImageReaderType::Pointer movingReader =<br>
MovingImageReaderType::New();<br>
movingReader->SetFileName( "1.mhd" );<br>
<br>
try<br>
{<br>
movingReader->Update();<br>
}<br>
catch( itk::ExceptionObject & e )<br>
{<br>
std::cout << e.GetDescription() << std::endl;<br>
return ;<br>
}<br>
<br>
typedef itk::GradientMagnitudeImageFilter<<br>
ITKImageType, ITKImageType > filterType;<br>
<br>
// Create and setup a gradient filter<br>
filterType::Pointer gradientFilter = filterType::New();<br>
gradientFilter->SetInput( movingReader->GetOutput() );<br>
gradientFilter->Update();<br>
<br>
RegistrationMethod->SetInitialTransformParameters(transform->GetParameters()<br>
);<br>
RegistrationMethod->SetTransform( transform );<br>
<br>
// Scale the translation components of the Transform in the Optimizer<br>
typedef OptimizerType::ScalesType OptimizerScalesType;<br>
OptimizerScalesType optimizerScales( transform->GetNumberOfParameters());<br>
<br>
const double translationScale = 1.0 ;<br>
optimizerScales[0] = 0.1;<br>
optimizerScales[1] = 0.1;<br>
optimizerScales[2] = 0.1;<br>
optimizerScales[3] = translationScale;<br>
optimizerScales[4] = translationScale;<br>
optimizerScales[5] = translationScale;<br>
optimizer->SetMaximumStepLength( 0.1 );<br>
optimizer->SetMinimumStepLength( 0.01 );<br>
optimizer->SetNumberOfIterations( 500 );<br>
optimizer->SetRelaxationFactor( 0.90 );<br>
optimizer->SetGradientMagnitudeTolerance( 0.005 );<br>
optimizer->MinimizeOn();<br>
RegistrationMethod->SetFixedPointSet( pointSet );<br>
RegistrationMethod->SetMetric( metric );<br>
RegistrationMethod->SetOptimizer( optimizer );<br>
RegistrationMethod->SetInterpolator( interpolator );<br>
RegistrationMethod->SetMovingImage( gradientFilter->GetOutput() );<br>
<br>
try<br>
{<br>
RegistrationMethod->Update();<br>
}<br>
<br>
catch( itk::ExceptionObject & e )<br>
{<br>
std::cout << "testing!!!" << std::endl;<br>
std::cerr << e << std::endl;<br>
ostrstream ostr;<br>
cout.rdbuf(ostr.rdbuf());<br>
std::cout << e << std::endl;<br>
}<br>
<br>
OptimizerType::ParametersType<br>
finalParameters=RegistrationMethod->GetLastTransformParameters();<br>
<br>
double versorX = finalParameters[0];<br>
double versorY = finalParameters[1];<br>
double versorZ = finalParameters[2];<br>
double finalTranslationX = finalParameters[3];<br>
double finalTranslationY = finalParameters[4];<br>
double finalTranslationZ = finalParameters[5];<br>
<br>
const unsigned int numberOfIterations = optimizer->GetCurrentIteration();<br>
const double bestValue = optimizer->GetValue();<br>
<br>
// Print out results<br>
//<br>
std::cout << std::endl << std::endl;<br>
std::cout << "Result = " << std::endl;<br>
std::cout << " versor X = " << versorX << std::endl;<br>
std::cout << " versor Y = " << versorY << std::endl;<br>
std::cout << " versor Z = " << versorZ << std::endl;<br>
std::cout << " Translation X = " << finalTranslationX << std::endl;<br>
std::cout << " Translation Y = " << finalTranslationY << std::endl;<br>
std::cout << " Translation Z = " << finalTranslationZ << std::endl;<br>
std::cout << " Iterations = " << numberOfIterations << std::endl;<br>
std::cout << " Metric value = " << bestValue << std::endl;<br>
}<br>
<br>
<br>
<br>
</div></div></div><br>
</div><br>