[Insight-users] PointSetToImageRegistrationMethod

[Luis Ibanez]

Hi Alon,

Thanks for the detailed description of your problem.


    1)  Why are you using RGB as pixel type ?
        Specially if you are using a CT scan as input...

    2)  You are using the NormalizedCorrelation metric.
        (I'm surprised to see that it compiles for an
        RGB pixel type...)


Please let us know...

In the meantime:


     a) Yes, what you are doing is a reasonable approach.
        However, please note that the border of an object
        is exactly the location where intensity values are
        the most variable.

        Instead of registering against the CT image, you
        may want to register against the gradient magnitude
        of the image.


     b) Extracting a surface from the CT and doing
        PointSetToPointSet registration is also a reasonable
        option. (at the price of requiring a registration).

     c) The final value of -0.5 for the normalized correlation
        seems to be an accident. The ideal value of a solution
        would be -1.0.


     d) Could you post to the list the sequence of values
        of the metric vs the iterations of the optimizer ?

     e) Not having seen the output of (d), from the parameters
        of your optimizer, the initial step of "20" seems to
        be very large. Specially considering that the scales
        that you are using for the angles are = 1.0.

        Therefore, when you set the value of the step, you
        should consider that the units are radians, and that
        20 radians is extremely large.

        You should start with something like 0.1 radias.



   Please let us know about the questions above.


     Regards,


        Luis


--------------------
Alon Mozes wrote:
> Hello,
>  
> I have been trying to register a pre-defined shape to a CT scan of a 
> patient with that shape also in the scan.  My approach has been to 
> create an .STL model of that shape, read it in as a series of points, 
> and use that as the fixedPointSet and the CT as the movingImage in the 
> PointSetToImageRegistrationMethod.  I am also using the 
> VersorRigid3DTransform and VersorRigid3DTransformOptimizer as part of 
> the registration (based on examples I'd seen).  The value passed in to 
> the method is a fixed value for each point representing the intensity of 
> the shape in the scan.  The params pointer contains the initial 
> transform (a 3D offset) based on the user picking a point in 3D.  I 
> display the shape point set overlaid on the contoured CT: before any 
> transformation, with the initial transformation, and with the final 
> transformation.  The initial transform puts the point set in 
> approximately the right place, but the resulting transform does not line 
> the point set up well with the contoured shape from the CT scan.  The 
> registration often causes it to wander a bit, converging after 100 or so 
> iterations.  I'm wondering:
>  
> 1)  Is this the best approach for registering what I need?  Maybe 
> intensity based comparisons are not ideal.  Perhaps I should try 
> registering the point set to the contoured CT using some distance metric 
> instead of intensity (like a more traditional surface match?)
>  
> 2) The optimizer->GetValue result once the optimizer is finished is 
> always about -0.5.  My understanding is that this is equivalent to an 
> RMS error, but then this value makes no sense.  What does this value 
> really mean?
>  
> 3) The CT has negative intensity values.  Does this affect the 
> registration? Are all positive intensities assumed?
>  
> Any help is greatly appreciated.
> Thanks,
> Alon
>  
> Here's my registration method:
>  
>  From the .h:
>  
> typedef itk::PointSetToImageRegistrationMethod<PointSetType, 
> ITKImageType > RegistrationType;
> RegistrationType::Pointer m_Registration;
>  
> The .cpp method:
>  
> void RegisterReference::Register(double params[6], int value)
> {
> const unsigned int Dimension = 3;
> typedef itk::RGBPixel< float > PixelType;
> //typedef itk::PointSet< PixelType, 3 > PointSetType;
> PointSetType::Pointer pointSet = PointSetType::New();
> vtkSTLReader * reader = vtkSTLReader::New();
> reader->SetFileName ("fiducial_fine.STL");
> reader->Update();
> vtkPolyData *inputMesh = reader->GetOutput();
> unsigned int pointId = 0;
> vtkPoints *points = inputMesh->GetPoints();
> for (int i=0; i<points->GetNumberOfPoints(); i++) {
> PointSetType::PointType p;
> p[0] = points->GetPoint(i)[0];
> p[1] = points->GetPoint(i)[1];
> p[2] = points->GetPoint(i)[2];
> pointSet->SetPoint( pointId, p );
> pointSet->SetPointData( pointId, value);
> ++pointId;
> }
> m_Registration = RegistrationType::New();
> 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;
> MetricType::Pointer metric = MetricType::New();
> OptimizerType::Pointer optimizer = OptimizerType::New();
> InterpolatorType::Pointer interpolator = InterpolatorType::New();
> m_Registration->SetMetric( metric );
> m_Registration->SetOptimizer( optimizer );
> m_Registration->SetInterpolator( interpolator );
> TransformType::Pointer transform = TransformType::New();
> transform->SetIdentity();
> m_Registration->SetInitialTransformParameters( transform->GetParameters() );
> RegistrationType::ParametersType 
> initialParameters(transform->GetNumberOfParameters() );
> initialParameters.Fill( 0.0 );
> initialParameters[3] = params[3];
> initialParameters[4] = params[4];
> initialParameters[5] = params[5];
> metric->SetTransform( transform );
> metric->SetInterpolator( interpolator );
> metric->SetFixedPointSet( pointSet );
> metric->SetMovingImage( m_MovingImage );
> m_Registration->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 / 1000.0;
> optimizerScales[0] = 1.0;
> optimizerScales[1] = 1.0;
> optimizerScales[2] = 1.0;
> optimizerScales[3] = translationScale;
> optimizerScales[4] = translationScale;
> optimizerScales[5] = translationScale;
> optimizer->SetMaximumStepLength( 20.0 );
> optimizer->SetMinimumStepLength( 0.001 );
> optimizer->SetNumberOfIterations( 500 );
> optimizer->SetRelaxationFactor( 0.90 );
> optimizer->SetGradientMagnitudeTolerance( 0.05 );
> optimizer->MinimizeOn();
> m_Registration->SetFixedPointSet( pointSet );
> m_Registration->SetMovingImage( m_MovingImage );
> ParametersType parameters( transform->GetNumberOfParameters() );
> // initialize the offset/vector part
> for( unsigned int k = 0; k < Dimension; k++ ) {
> parameters[k]= 0.0f;
> }
> //Prepare the observer to listen for registration iterations
> RegistrationObserver::Pointer observer = RegistrationObserver::New();
> optimizer->AddObserver( itk::IterationEvent(), observer );
> optimizer->AddObserver( itk::StartEvent(), observer );
> try
> {
> m_Registration->SetInitialTransformParameters( initialParameters );
> m_Registration->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;
> char* s = ostr.str();
> AfxMessageBox(s);
> }
> OptimizerType::ParametersType finalParameters 
> =m_Registration->GetLastTransformParameters();
> params[0] = finalParameters[0];
> params[1] = finalParameters[1];
> params[2] = finalParameters[2];
> params[3] = finalParameters[3];
> params[4] = finalParameters[4];
> params[5] = finalParameters[5];
> const unsigned int nmbrOfIterations = optimizer->GetCurrentIteration();
> const double bestValue = optimizer->GetValue();
> }
>  
> 
> 
> 
> ------------------------------------------------------------------------
> 
> _______________________________________________
> Insight-users mailing list
> Insight-users at itk.org
> http://www.itk.org/mailman/listinfo/insight-users