[Insight-users] Mesh to Distance Map registration problem

[SWAP 2000 - 2001]

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;
    }
    
}


_________________________________________________________________
Windows Live™ Hotmail®:…more than just e-mail.
http://windowslive.com/online/hotmail?ocid=TXT_TAGLM_WL_HM_more_042009
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20090426/0c1b5fcb/attachment.htm>