[Insight-users] PointSetToImageRegistrationMethod
Alon Mozes
amozes77 at yahoo.com
Sun Jan 4 22:48:48 EST 2009
Hi Luis,
Thanks so much for your reply! I think your advice is spot on. Specifically...
1) You're right. RGB is not actually used beyond the typedef. Unfortunately, it's messy code that has some leftovers from all the variations I've tried. I neglected to delete that line.
2) Normalized correlation metric seems to work with the code I'm using.
a) I verified in VolView that the edge is, as you mentioned, quite variable. I have now created a solid point cloud (instead of just the surface point model) to try to mitigate that. Along with other fixes you suggested (see e), it seems to work reasonably well (at least a 90% solution) but I may need to try something else if I can't get that last 10%. I looked at the gradient magnitude also in VolView and the edge looks a little wide which I'm guessing could lead to the same uncertainty in the registration. I'm guessing the perfect match can only be achieved if the edge in the scan is sharp regardless of what method I use.
b) I do extract a surface using contours (just for visualization) but it doesn't look especially accurate. I'd like to avoid a more involved segmentation if I don't need an accurate model.
c) I can get values of -0.7 to -0.9 now, depending on the initial guess. Since the edge falls off in value in the scan, I'm guessing a perfect -1.0 is too much to ask for.
d) Since it's working, I haven't had the chance to dump the value at each iteration.
e) This was it!! My initial rotation step value of 20 was too high (I must've played around with that value assuming it was degrees and never revisited it). Changing it to 0.1 radians did the trick. The points line up reasonably well with the scan. Now I just have to get it to be rock solid every time. I'm thinking improving the initial guess will help a lot. Also, if my final result value is worse than -0.8, I should bump the initial guess and try to register again.
I'm open to any other ideas to perfecting the solution. Thanks again for your help!
Alon
--- On Wed, 12/31/08, Luis Ibanez <luis.ibanez at kitware.com> wrote:
From: Luis Ibanez <luis.ibanez at kitware.com>
Subject: Re: [Insight-users] PointSetToImageRegistrationMethod
To: "Alon Mozes" <amozes77 at yahoo.com>
Cc: insight-users at itk.org
Date: Wednesday, December 31, 2008, 11:31 AM
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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20090104/f20502d9/attachment.htm>
More information about the Insight-users
mailing list