[Insight-users] Mesh to Distance Map registration problem

[Luis Ibanez]

Hi "SWAP 2000 - 2001"

Thanks for detailed description of the problem.

Your approach seems to be quite reasonable.

The first question I have, after looking at the code,
is:

     How are you initializing the translation
     of the transform ?

It seems that you start with an identity,...

Can you post to the list the output generated by the
Observer that should be connected to the optimizer ?

This is the output that lists the Metric values at
every iteration along with the transform parameters.

--

In general, my suggestion will be to first debug
the code by using a trivial example:  E.g. a small
image of a Sphere in a 100x100x100 pixels image,
and a mesh with a sphere of the matching size and
position.


Once you get that  example to work, you could move
into a more realistic scenario, such as the one with
the vertebral image.


     Regards,


         Luis



------------------------
SWAP 2000 - 2001 wrote:
> Hi everyone,
> 
> 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.
> 
> 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.
> 
> 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.
> 
> Any help would be GREATLY appreciated, I've exhausted all other resources.
> 
> PS cleaned up relevant source code included below
> 
> //-----------------------------------------------------------
> // setPointsToZero: reset ITK Mesh PointData to zero
> //-----------------------------------------------------------
> void setPointsToZero (MeshType::Pointer &templateMesh){
> 
>     PixelTypeDouble zero = (0.0,0.0,0.0);
> 
>     // Get the number of points in the mesh
>    int numPoints = templateMesh->GetNumberOfPoints();
>    if(numPoints == 0)
>     {
>         templateMesh->Print(std::cerr << "Mesh has no points." << 
> std::endl);
>     }   
>   
>     // Iterate over all the points in the itk mesh seting values to 0   
>     MeshType::PointsContainer::Pointer points = templateMesh->GetPoints();
>     for(MeshType::PointsContainer::Iterator i = points->Begin(); i != 
> points->End(); i++){
>        
>         // Get the point index from the point container iterator
>         templateMesh->SetPointData(i->Index(), zero);
> 
>     }
>    
> }
> 
> //--------------------------------------------------------------
> // registerMeshToImg: 3d rigid register distance map to ITK Mesh
> //--------------------------------------------------------------
> void registerMeshToImg(MeshType::Pointer &meshTemplate, 
> ImageTypeDouble::Pointer &imgDistanceMap)
> {
> 
>   //-----------------------------------------------------------
>   // Set up  the Metric
>   //-----------------------------------------------------------
>   typedef itk::MeanSquaresPointSetToImageMetric< 
>                                        MeshType,
>                                        ImageTypeDouble >  
>                                                     MetricType;
> 
>   typedef MetricType::TransformType                 TransformBaseType;
>   typedef TransformBaseType::ParametersType         ParametersType;
>   typedef TransformBaseType::JacobianType           JacobianType;
> 
>   MetricType::Pointer  metric = MetricType::New();
> 
>   //-----------------------------------------------------------
>   // Set up a Transform
>   //-----------------------------------------------------------
>   typedef itk::Similarity3DTransform< PixelTypeDouble > TransformType;
> 
>   TransformType::Pointer transform = TransformType::New();
> 
>   //------------------------------------------------------------
>   // Set up an Interpolator
>   //------------------------------------------------------------
>   typedef itk::LinearInterpolateImageFunction<
>     ImageTypeDouble,
>     PixelTypeDouble > InterpolatorType;
> 
>   InterpolatorType::Pointer interpolator = InterpolatorType::New();
> 
>   // Optimizer Type
>   typedef itk::RegularStepGradientDescentOptimizer       OptimizerType;
> 
>   OptimizerType::Pointer      optimizer     = OptimizerType::New();
> 
>   // Registration Method
>   typedef itk::PointSetToImageRegistrationMethod<
>     MeshType,
>     ImageTypeDouble >    RegistrationType;
> 
> 
>   RegistrationType::Pointer   registration  = RegistrationType::New();
> 
>   // Scale the components of the Transform in the Optimizer
>   OptimizerType::ScalesType scales( transform->GetNumberOfParameters() );
>   scales.Fill( 1.0 );
> 
>   unsigned long   numberOfIterations =   500;
>   double          maximumStepLenght  =  0.001;  // no step will be 
> larger than this
>   double          minimumStepLenght  =  0.0001; // convergence criterion
>   double         gradientTolerance  =  1e-6; // convergence criterion
>  
>   optimizer->SetScales( scales );
>   optimizer->SetNumberOfIterations( numberOfIterations );
>   optimizer->SetMinimumStepLength( minimumStepLenght );
>   optimizer->SetMaximumStepLength( maximumStepLenght );
>   optimizer->SetGradientMagnitudeTolerance( gradientTolerance );
>   optimizer->MinimizeOn();
>  
>   // Start from an Identity transform
>   transform->SetIdentity();
>   transform->SetCenter( /* Center of Gravity of Mesh */ );
>  
>   registration->SetInitialTransformParameters( transform->GetParameters() );
>  
>   //------------------------------------------------------
>   // Connect all the components required for Registration
>   //------------------------------------------------------
>   registration->SetMetric(        metric        );
>   registration->SetOptimizer(     optimizer     );
>   registration->SetTransform(     transform     );
>   registration->SetFixedPointSet( meshTemplate );
>   registration->SetMovingImage(   imgDistanceMap   );
>   registration->SetInterpolator(  interpolator  );
>  
>   //------------------------------------------------------------
>   // Set up transform parameters
>   //------------------------------------------------------------
>   ParametersType parameters( transform->GetNumberOfParameters() );
> 
>   try
>     {
>     registration->StartRegistration();
>     }
>   catch( itk::ExceptionObject & e )
>     {
>     std::cout << e << std::endl;
>     }
>    
> }
> _____________________________________
> Powered by www.kitware.com
> 
> Visit other Kitware open-source projects at
> http://www.kitware.com/opensource/opensource.html
> 
> Please keep messages on-topic and check the ITK FAQ at: http://www.itk.org/Wiki/ITK_FAQ
> 
> Follow this link to subscribe/unsubscribe:
> http://www.itk.org/mailman/listinfo/insight-users