<html>
<head>
<style>
.hmmessage P
{
margin:0px;
padding:0px
}
body.hmmessage
{
font-size: 10pt;
font-family:Verdana
}
</style>
</head>
<body class='hmmessage'>
Hi everyone,<br><br>I have spent days trying to resolve this problem and hope someone out there has dealt with this before. I want to rigidly register an image (moving component) to a mesh (fixed component). The image in question is a vertebra I segmented using ITK-SNAP and and exported as MHD file (so 3d image). The mesh is VTK Polydata mesh exported from ITK-SNAP from another segmented vertebra.<br><br>Before registration I decimate the mesh (vtkDecimatePro), write it back to file and read it as an ITK Mesh. I take that mesh an set all PointData to zero. I do this so I can register a signed danielsson distance map to the mesh. The purpose of registering a distance map to a zeroed mesh is to reduce computational overhead. However, when I get the transformation it appears that it's not using the distance map/zeroed mesh properly. I should be seeing a translation of (10,45,35) but instead get a translation that approaches 0. The same goes for rotation and scaling, the "better" the optimizer settings are, the worse the transformation is (eg extreme shrinking, rotation that results in perpendicular alignment). It always goes till the max iteration setting. I verified all this using Paraview and even hard coding the correct translation (the correct rotation/scaling already near identity) it screws things up.<br><br>I have looked at the ICP examples and the PointSettoImage test examples and obviously none of them deal with this problem. The closest thing I found was in ICP example 3, where distance map is used for the fixed pointset (EuclideanMetric->SetDistanceMap). The problem is, my overall project depends on having different vertebre images registering to the same fixed mesh.<br><br>Any help would be GREATLY appreciated, I've exhausted all other resources.<br><br>PS cleaned up relevant source code included below<br><br>//-----------------------------------------------------------<br>// setPointsToZero: reset ITK Mesh PointData to zero<br>//-----------------------------------------------------------<br>void setPointsToZero (MeshType::Pointer &templateMesh){<br><br> PixelTypeDouble zero = (0.0,0.0,0.0);<br><br> // Get the number of points in the mesh<br> int numPoints = templateMesh->GetNumberOfPoints();<br> if(numPoints == 0)<br> {<br> templateMesh->Print(std::cerr << "Mesh has no points." << std::endl);<br> } <br> <br> // Iterate over all the points in the itk mesh seting values to 0 <br> MeshType::PointsContainer::Pointer points = templateMesh->GetPoints();<br> for(MeshType::PointsContainer::Iterator i = points->Begin(); i != points->End(); i++){<br> <br> // Get the point index from the point container iterator<br> templateMesh->SetPointData(i->Index(), zero);<br><br> }<br> <br>}<br><br>//--------------------------------------------------------------<br>// registerMeshToImg: 3d rigid register distance map to ITK Mesh<br>//--------------------------------------------------------------<br>void registerMeshToImg(MeshType::Pointer &meshTemplate, ImageTypeDouble::Pointer &imgDistanceMap)<br>{<br><br> //-----------------------------------------------------------<br> // Set up the Metric<br> //-----------------------------------------------------------<br> typedef itk::MeanSquaresPointSetToImageMetric< <br> MeshType, <br> ImageTypeDouble > <br> MetricType;<br><br> typedef MetricType::TransformType TransformBaseType;<br> typedef TransformBaseType::ParametersType ParametersType;<br> typedef TransformBaseType::JacobianType JacobianType;<br><br> MetricType::Pointer metric = MetricType::New();<br><br> //-----------------------------------------------------------<br> // Set up a Transform<br> //-----------------------------------------------------------<br> typedef itk::Similarity3DTransform< PixelTypeDouble > TransformType;<br><br> TransformType::Pointer transform = TransformType::New();<br><br> //------------------------------------------------------------<br> // Set up an Interpolator<br> //------------------------------------------------------------<br> typedef itk::LinearInterpolateImageFunction< <br> ImageTypeDouble,<br> PixelTypeDouble > InterpolatorType;<br><br> InterpolatorType::Pointer interpolator = InterpolatorType::New();<br><br> // Optimizer Type<br> typedef itk::RegularStepGradientDescentOptimizer OptimizerType;<br><br> OptimizerType::Pointer optimizer = OptimizerType::New();<br><br> // Registration Method<br> typedef itk::PointSetToImageRegistrationMethod< <br> MeshType, <br> ImageTypeDouble > RegistrationType;<br><br><br> RegistrationType::Pointer registration = RegistrationType::New();<br><br> // Scale the components of the Transform in the Optimizer<br> OptimizerType::ScalesType scales( transform->GetNumberOfParameters() );<br> scales.Fill( 1.0 );<br><br> unsigned long numberOfIterations = 500;<br> double maximumStepLenght = 0.001; // no step will be larger than this<br> double minimumStepLenght = 0.0001; // convergence criterion<br> double gradientTolerance = 1e-6; // convergence criterion<br> <br> optimizer->SetScales( scales );<br> optimizer->SetNumberOfIterations( numberOfIterations );<br> optimizer->SetMinimumStepLength( minimumStepLenght );<br> optimizer->SetMaximumStepLength( maximumStepLenght );<br> optimizer->SetGradientMagnitudeTolerance( gradientTolerance );<br> optimizer->MinimizeOn();<br> <br> // Start from an Identity transform<br> transform->SetIdentity();<br> transform->SetCenter( /* Center of Gravity of Mesh */ );<br> <br> registration->SetInitialTransformParameters( transform->GetParameters() );<br> <br> //------------------------------------------------------<br> // Connect all the components required for Registration<br> //------------------------------------------------------<br> registration->SetMetric( metric );<br> registration->SetOptimizer( optimizer );<br> registration->SetTransform( transform );<br> registration->SetFixedPointSet( meshTemplate );<br> registration->SetMovingImage( imgDistanceMap );<br> registration->SetInterpolator( interpolator );<br> <br> //------------------------------------------------------------<br> // Set up transform parameters<br> //------------------------------------------------------------<br> ParametersType parameters( transform->GetNumberOfParameters() );<br><br> try <br> {<br> registration->StartRegistration();<br> }<br> catch( itk::ExceptionObject & e )<br> {<br> std::cout << e << std::endl;<br> }<br> <br>}<br><br><br /><hr />Windows Live™ Hotmail®:…more than just e-mail. <a href='http://windowslive.com/online/hotmail?ocid=TXT_TAGLM_WL_HM_more_042009' target='_new'>Check it out.</a></body>
</html>