No subject
Mon Dec 15 09:51:36 EST 2008
=A0
typedef itk::PointSetToImageRegistrationMethod<PointSetType, ITKImageType >=
RegistrationType;
RegistrationType::Pointer m_Registration;
=A0
The .cpp method:
=A0
void RegisterReference::Register(double params[6], int value)
{
const unsigned int Dimension =3D 3;
typedef itk::RGBPixel< float > PixelType;
//typedef itk::PointSet< PixelType, 3 > PointSetType;
PointSetType::Pointer pointSet =3D PointSetType::New();
vtkSTLReader * reader =3D vtkSTLReader::New();
reader->SetFileName ("fiducial_fine.STL");
reader->Update();
vtkPolyData *inputMesh =3D reader->GetOutput();
unsigned int pointId =3D 0;
vtkPoints *points =3D inputMesh->GetPoints();
for (int i=3D0; i<points->GetNumberOfPoints(); i++) {
PointSetType::PointType p;
p[0] =3D points->GetPoint(i)[0];
p[1] =3D points->GetPoint(i)[1];
p[2] =3D points->GetPoint(i)[2];
pointSet->SetPoint( pointId, p );
pointSet->SetPointData( pointId, value);
++pointId;
}
m_Registration =3D RegistrationType::New();
typedef itk::VersorRigid3DTransform< double > TransformType;
typedef itk::VersorRigid3DTransformOptimizer OptimizerType;
typedef itk::NormalizedCorrelationPointSetToImageMetric<=20
PointSetType,=20
ITKImageType > MetricType;=20
typedef itk:: LinearInterpolateImageFunction<=20
ITKImageType,
double > InterpolatorType;
typedef MetricType::TransformType TransformBaseType;
typedef TransformBaseType::ParametersType ParametersType;
MetricType::Pointer metric =3D MetricType::New();
OptimizerType::Pointer optimizer =3D OptimizerType::New();
InterpolatorType::Pointer interpolator =3D InterpolatorType::New();
m_Registration->SetMetric( metric );
m_Registration->SetOptimizer( optimizer );
m_Registration->SetInterpolator( interpolator );
TransformType::Pointer transform =3D TransformType::New();
transform->SetIdentity();
m_Registration->SetInitialTransformParameters( transform->GetParameters() )=
;
RegistrationType::ParametersType initialParameters(transform->GetNumberOfPa=
rameters() );
initialParameters.Fill( 0.0 );
initialParameters[3] =3D params[3];=20
initialParameters[4] =3D params[4];=20
initialParameters[5] =3D params[5];
metric->SetTransform( transform );=20
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 =3D 1.0 / 1000.0;
optimizerScales[0] =3D 1.0;
optimizerScales[1] =3D 1.0;
optimizerScales[2] =3D 1.0;
optimizerScales[3] =3D translationScale;
optimizerScales[4] =3D translationScale;
optimizerScales[5] =3D 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 =3D 0; k < Dimension; k++ ) {
parameters[k]=3D 0.0f;
}
//Prepare the observer to listen for registration iterations
RegistrationObserver::Pointer observer =3D RegistrationObserver::New();
optimizer->AddObserver( itk::IterationEvent(), observer );
optimizer->AddObserver( itk::StartEvent(), observer );
try=20
{
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());=20
std::cout << e << std::endl;
char* s =3D ostr.str();=20
AfxMessageBox(s);
}
OptimizerType::ParametersType finalParameters =3Dm_Registration->GetLastTra=
nsformParameters();
params[0] =3D finalParameters[0];
params[1] =3D finalParameters[1];
params[2] =3D finalParameters[2];
params[3] =3D finalParameters[3];
params[4] =3D finalParameters[4];
params[5] =3D finalParameters[5];
const unsigned int nmbrOfIterations =3D optimizer->GetCurrentIteration();
const double bestValue =3D optimizer->GetValue();
}
=A0=0A=0A=0A
--0-112149324-1230443879=:27989
Content-Type: text/html; charset=us-ascii
<table cellspacing="0" cellpadding="0" border="0" ><tr><td valign="top" style="font: inherit;"><DIV>Hello,</DIV>
<DIV> </DIV>
<DIV>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:</DIV>
<DIV> </DIV>
<DIV>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?)</DIV>
<DIV> </DIV>
<DIV>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?</DIV>
<DIV> </DIV>
<DIV>3) The CT has negative intensity values. Does this affect the registration? Are all positive intensities assumed?</DIV>
<DIV> </DIV>
<DIV>Any help is greatly appreciated.</DIV>
<DIV>Thanks,</DIV>
<DIV>Alon</DIV>
<DIV> </DIV>
<DIV>Here's my registration method:</DIV>
<DIV> </DIV>
<DIV>From the .h:</DIV>
<DIV> </DIV><FONT size=2>
<DIV></FONT><FONT color=#0000ff size=2>typedef</FONT><FONT size=2> itk::PointSetToImageRegistrationMethod<PointSetType, ITKImageType > RegistrationType;</DIV></FONT><FONT size=2>
<DIV>RegistrationType::Pointer m_Registration;</DIV></FONT>
<DIV> </DIV>
<DIV>The .cpp method:</DIV>
<DIV> </DIV><FONT color=#0000ff size=2>
<DIV>void</FONT><FONT size=2> RegisterReference::Register(</FONT><FONT color=#0000ff size=2>double</FONT><FONT size=2> params[6], </FONT><FONT color=#0000ff size=2>int</FONT><FONT size=2> value)</DIV>
<DIV>{</DIV>
<DIV></FONT><FONT color=#0000ff size=2>const</FONT><FONT size=2> </FONT><FONT color=#0000ff size=2>unsigned</FONT><FONT size=2> </FONT><FONT color=#0000ff size=2>int</FONT><FONT size=2> Dimension = 3;</DIV>
<DIV></FONT><FONT color=#0000ff size=2>typedef</FONT><FONT size=2> itk::RGBPixel< </FONT><FONT color=#0000ff size=2>float</FONT><FONT size=2> > PixelType;</DIV>
<DIV></FONT><FONT color=#008000 size=2>//typedef itk::PointSet< PixelType, 3 > PointSetType;</DIV></FONT><FONT size=2>
<DIV>PointSetType::Pointer pointSet = PointSetType::New();</DIV>
<DIV>vtkSTLReader * reader = vtkSTLReader::New();</DIV>
<DIV>reader->SetFileName ("fiducial_fine.STL");</DIV>
<DIV>reader->Update();</DIV>
<DIV>vtkPolyData *inputMesh = reader->GetOutput();</DIV>
<DIV></FONT><FONT color=#0000ff size=2>unsigned</FONT><FONT size=2> </FONT><FONT color=#0000ff size=2>int</FONT><FONT size=2> pointId = 0;</DIV>
<DIV>vtkPoints *points = inputMesh->GetPoints();</DIV>
<DIV></FONT><FONT color=#0000ff size=2>for</FONT><FONT size=2> (</FONT><FONT color=#0000ff size=2>int</FONT><FONT size=2> i=0; i<points->GetNumberOfPoints(); i++) {</DIV>
<DIV>PointSetType::PointType p;</DIV>
<DIV>p[0] = points->GetPoint(i)[0];</DIV>
<DIV>p[1] = points->GetPoint(i)[1];</DIV>
<DIV>p[2] = points->GetPoint(i)[2];</DIV>
<DIV>pointSet->SetPoint( pointId, p );</DIV>
<DIV>pointSet->SetPointData( pointId, value);</DIV>
<DIV>++pointId;</DIV>
<DIV>}</DIV>
<DIV>m_Registration = RegistrationType::New();</DIV>
<DIV></FONT><FONT color=#0000ff size=2>typedef</FONT><FONT size=2> itk::VersorRigid3DTransform< </FONT><FONT color=#0000ff size=2>double</FONT><FONT size=2> > TransformType;</DIV>
<DIV></FONT><FONT color=#0000ff size=2>typedef</FONT><FONT size=2> itk::VersorRigid3DTransformOptimizer OptimizerType;</DIV>
<DIV></FONT><FONT color=#0000ff size=2>typedef</FONT><FONT size=2> itk::NormalizedCorrelationPointSetToImageMetric< </DIV>
<DIV>PointSetType, </DIV>
<DIV>ITKImageType > MetricType; </DIV>
<DIV></DIV>
<DIV></FONT><FONT color=#0000ff size=2>typedef</FONT><FONT size=2> itk:: LinearInterpolateImageFunction< </DIV>
<DIV>ITKImageType,</DIV>
<DIV></FONT><FONT color=#0000ff size=2>double</FONT><FONT size=2> > InterpolatorType;</DIV>
<DIV></FONT><FONT color=#0000ff size=2>typedef</FONT><FONT size=2> MetricType::TransformType TransformBaseType;</DIV>
<DIV></FONT><FONT color=#0000ff size=2>typedef</FONT><FONT size=2> TransformBaseType::ParametersType ParametersType;</DIV>
<DIV>MetricType::Pointer metric = MetricType::New();</DIV>
<DIV>OptimizerType::Pointer optimizer = OptimizerType::New();</DIV>
<DIV>InterpolatorType::Pointer interpolator = InterpolatorType::New();</DIV>
<DIV></DIV>
<DIV>m_Registration->SetMetric( metric );</DIV>
<DIV>m_Registration->SetOptimizer( optimizer );</DIV>
<DIV>m_Registration->SetInterpolator( interpolator );</DIV>
<DIV>TransformType::Pointer transform = TransformType::New();</DIV>
<DIV>transform->SetIdentity();</DIV>
<DIV>m_Registration->SetInitialTransformParameters( transform->GetParameters() );</DIV>
<DIV>RegistrationType::ParametersType initialParameters(transform->GetNumberOfParameters() );</DIV>
<DIV>initialParameters.Fill( 0.0 );</DIV>
<DIV>initialParameters[3] = params[3]; </DIV>
<DIV>initialParameters[4] = params[4]; </DIV>
<DIV>initialParameters[5] = params[5];</DIV>
<DIV>metric->SetTransform( transform ); </DIV>
<DIV>metric->SetInterpolator( interpolator );</DIV>
<DIV>metric->SetFixedPointSet( pointSet );</DIV>
<DIV>metric->SetMovingImage( m_MovingImage );</DIV>
<DIV>m_Registration->SetTransform( transform );</DIV>
<DIV></FONT><FONT color=#008000 size=2>// Scale the translation components of the Transform in the Optimizer</DIV></FONT><FONT size=2>
<DIV></FONT><FONT color=#0000ff size=2>typedef</FONT><FONT size=2> OptimizerType::ScalesType OptimizerScalesType;</DIV>
<DIV>OptimizerScalesType optimizerScales( transform->GetNumberOfParameters() );</DIV>
<DIV></FONT><FONT color=#0000ff size=2>const</FONT><FONT size=2> </FONT><FONT color=#0000ff size=2>double</FONT><FONT size=2> translationScale = 1.0 / 1000.0;</DIV>
<DIV>optimizerScales[0] = 1.0;</DIV>
<DIV>optimizerScales[1] = 1.0;</DIV>
<DIV>optimizerScales[2] = 1.0;</DIV>
<DIV>optimizerScales[3] = translationScale;</DIV>
<DIV>optimizerScales[4] = translationScale;</DIV>
<DIV>optimizerScales[5] = translationScale;</DIV>
<DIV>optimizer->SetMaximumStepLength( 20.0 );</DIV>
<DIV>optimizer->SetMinimumStepLength( 0.001 );</DIV>
<DIV>optimizer->SetNumberOfIterations( 500 );</DIV>
<DIV>optimizer->SetRelaxationFactor( 0.90 );</DIV>
<DIV>optimizer->SetGradientMagnitudeTolerance( 0.05 );</DIV>
<DIV>optimizer->MinimizeOn();</DIV>
<DIV>m_Registration->SetFixedPointSet( pointSet );</DIV>
<DIV>m_Registration->SetMovingImage( m_MovingImage );</DIV>
<DIV>ParametersType parameters( transform->GetNumberOfParameters() );</DIV>
<DIV></FONT><FONT color=#008000 size=2>// initialize the offset/vector part</DIV></FONT><FONT size=2>
<DIV></FONT><FONT color=#0000ff size=2>for</FONT><FONT size=2>( </FONT><FONT color=#0000ff size=2>unsigned</FONT><FONT size=2> </FONT><FONT color=#0000ff size=2>int</FONT><FONT size=2> k = 0; k < Dimension; k++ ) {</DIV>
<DIV>parameters[k]= 0.0f;</DIV>
<DIV>}</DIV>
<DIV></FONT><FONT color=#008000 size=2>//Prepare the observer to listen for registration iterations</DIV></FONT><FONT size=2>
<DIV>RegistrationObserver::Pointer observer = RegistrationObserver::New();</DIV>
<DIV>optimizer->AddObserver( itk::IterationEvent(), observer );</DIV>
<DIV>optimizer->AddObserver( itk::StartEvent(), observer );</DIV>
<DIV></FONT><FONT color=#0000ff size=2>try</FONT><FONT size=2> </DIV>
<DIV>{</DIV>
<DIV>m_Registration->SetInitialTransformParameters( initialParameters );</DIV>
<DIV>m_Registration->Update();</DIV>
<DIV>}</DIV>
<DIV></FONT><FONT color=#0000ff size=2>catch</FONT><FONT size=2>( itk::ExceptionObject & e )</DIV>
<DIV>{</DIV>
<DIV>std::cout << "testing!!!" << std::endl;</DIV>
<DIV>std::cerr << e << std::endl;</DIV>
<DIV>ostrstream ostr;</DIV>
<DIV>cout.rdbuf(ostr.rdbuf()); </DIV>
<DIV>std::cout << e << std::endl;</DIV>
<DIV></FONT><FONT color=#0000ff size=2>char</FONT><FONT size=2>* s = ostr.str(); </DIV>
<DIV>AfxMessageBox(s);</DIV>
<DIV>}</DIV>
<DIV>OptimizerType::ParametersType finalParameters =m_Registration->GetLastTransformParameters();</DIV>
<DIV>params[0] = finalParameters[0];</DIV>
<DIV>params[1] = finalParameters[1];</DIV>
<DIV>params[2] = finalParameters[2];</DIV>
<DIV>params[3] = finalParameters[3];</DIV>
<DIV>params[4] = finalParameters[4];</DIV>
<DIV>params[5] = finalParameters[5];</DIV>
<DIV></FONT><FONT color=#0000ff size=2>const</FONT><FONT size=2> </FONT><FONT color=#0000ff size=2>unsigned</FONT><FONT size=2> </FONT><FONT color=#0000ff size=2>int</FONT><FONT size=2> nmbrOfIterations = optimizer->GetCurrentIteration();</DIV>
<DIV></FONT><FONT color=#0000ff size=2>const</FONT><FONT size=2> </FONT><FONT color=#0000ff size=2>double</FONT><FONT size=2> bestValue = optimizer->GetValue();</DIV>
<DIV>}</DIV></FONT>
<DIV> </DIV></td></tr></table><br>
--0-112149324-1230443879=:27989--
More information about the Insight-users
mailing list