[Insight-users] pointsettoimageregistration - setup & transform results

[Dean Inglis]

some additional info.  Attached is a composite of the image (a
rescaled gradient magnitude image of the segmented head surface)
and the original points (green) and the result of the registration (red).
I have initialized the transform by calculating the centroid of the
image using itkImageMomentsCalculator Center of Gravity:

  transform->SetIdentity();
  // center of the fixed (red) points
  TransformType::InputPointType centerFixed;
  centerFixed[0] = cx;
  centerFixed[1] = cy;
  centerFixed[2] = cz;
  transform->SetCenter(centerFixed);

  TransformType::OutputVectorType translation;
  translation[0] = ccg[0] - cx;
  translation[1] = ccg[1] - cy;
  translation[2] = ccg[2] - cz;
  transform->SetTranslation( translation );

and stopped using transform->SetParameters( parameters );
I have also initialized the VersorRigid3DTransformOptimizer optimizer 
scales:

  OptimizerType::ScalesType scales( transform->GetNumberOfParameters() );
  scales.Fill( 1.0 );

  scales[3] = 1./1000.;
  scales[4] = 1./1000.;
  scales[5] = 1./1000.;

However, as seen in the attached image, the points are nowhere near a 
gradient
edge.  Any ideas / suggestions would be greatly appreciated.

best,
Dean


>I am registering a set of 3D points taken using a touch probe
> over a subject's face to a 3D CT scan of the subject.  My pipeline
> compiles and now runs without crashing but the results are
> not acceptable: the output points are not positioned remotely
> correctly.   I have a few issues that I would could use some
> guidance on
>
> 1) does the pipeline posted below look reasonable or are there
> better components that should be considered (e.g., different metric,
> interpolator etc)
>
> 2) what parameters need to be tuned and how should they be
> initialized?
>
> 3) since the points need to be rotated and translated, I am
> using  a VersorRigid3DTransform.  Do I need to initialize it
> with the center of the fixed points?  When I transform the input
> points to get the final output points, is it sufficient to just do
>    FixedPointType inputPt  = fixedPointSet->GetPoint( i );
>    FixedPointType outputPt = transform->TransformPoint( inputPt );
> or do I need an inverse transform and do I also need to take into
> account center?
>
> Dean
>
>  typedef itk::StatisticsImageFilter< UnsignedCharImageType > 
> StatisticsImageFilterType;
>  StatisticsImageFilterType::Pointer statisticsImageFilter
>          = StatisticsImageFilterType::New ();
>  statisticsImageFilter->SetInput( castFilter->GetOutput() );
>  statisticsImageFilter->Update();
>
>  // set the values associated with the points to be a factor of the stdev 
> + mean of the cast output
>  float defaultScaleFactor = -1.0;
>  defaultValue = static_cast< FixedPointDataType >( defaultScaleFactor * 
> statisticsImageFilter->GetSigma() + statisticsImageFilter->GetMean() );
>
>  for(int i = 0 ; i < fixedPointSet->GetNumberOfPoints(); ++i )
>  {
>    fixedPointSet->SetPointData( i, defaultValue );
>  }
>
>  typedef itk::ShrinkImageFilter <UnsignedCharImageType, 
> UnsignedCharImageType> ShrinkImageFilterType;
>  ShrinkImageFilterType::Pointer shrinkFilter = 
> ShrinkImageFilterType::New();
>  shrinkFilter->SetInput(castFilter->GetOutput());
>  shrinkFilter->SetShrinkFactor(0, 2); // shrink the dimension by a factor 
> of 2 (i.e. 100 gets changed to 50)
>  shrinkFilter->SetShrinkFactor(1, 2);
>  shrinkFilter->SetShrinkFactor(2, 2);
>  shrinkFilter->Update();
>
>  typedef double CoordinateRepresentationType;
>  typedef itk::VersorRigid3DTransform< CoordinateRepresentationType   > 
> TransformType;
>  TransformType::Pointer transform = TransformType::New();
>  transform->SetIdentity();
>
>  TransformType::InputPointType centerFixed;
>  centerFixed[0] = cx;
>  centerFixed[1] = cy;
>  centerFixed[2] = cz;
>  transform->SetCenter(centerFixed);
>
>  typedef itk::MeanSquaresPointSetToImageMetric< FixedPointSetType, 
> UnsignedCharImageType > MetricType;
>  MetricType::Pointer metric = MetricType::New();
>  typedef MetricType::TransformType TransformBaseType;
>  typedef TransformBaseType::ParametersType ParametersType;
>
>    ParametersType parameters( transform->GetNumberOfParameters() );
>
>  //  Versor type
>    typedef    TransformType::VersorType      VersorType;
>    VersorType versor;
>
>    parameters[0] = versor.GetX();   // Rotation axis * sin(t/2)
>    parameters[1] = versor.GetY();
>    parameters[2] = versor.GetZ();
>    parameters[3] = 0.0;             // Translation
>    parameters[4] = 0.0;
>    parameters[5] = 0.0;
>
>    transform->SetParameters( parameters );
>
>  typedef itk:: NearestNeighborInterpolateImageFunction<
>                                    UnsignedCharImageType, 
> CoordinateRepresentationType     >    InterpolatorType;
>  InterpolatorType::Pointer interpolator = InterpolatorType::New();
>  //interpolator->DebugOn();
>
>  typedef itk::VersorRigid3DTransformOptimizer OptimizerType;
>  OptimizerType::Pointer optimizer = OptimizerType::New();
>
>  unsigned long   numberOfIterations =   50;
>  double          maximumStepLength  =  1.0;  // no step will be larger 
> than this
>  double          minimumStepLength  =  0.01; // convergence criterion
>  double          gradientTolerance  =  1e-6; // convergence criterion
>
>  OptimizerType::ScalesType scales( transform->GetNumberOfParameters() );
>  scales.Fill( 1.0 );
>
>  optimizer->SetScales( scales );
>  optimizer->SetNumberOfIterations( numberOfIterations );
>  optimizer->SetMinimumStepLength( minimumStepLength );
>  optimizer->SetMaximumStepLength( maximumStepLength );
>  optimizer->SetGradientMagnitudeTolerance( gradientTolerance );
>  optimizer->MinimizeOn();
>
>
>  typedef IterationCallback< OptimizerType >   IterationCallbackType;
>  IterationCallbackType::Pointer callback = IterationCallbackType::New();
>  callback->SetOptimizer( optimizer );
>
>  typedef itk::PointSetToImageRegistrationMethod< FixedPointSetType, 
> UnsignedCharImageType > RegistrationType;
>  RegistrationType::Pointer registration = RegistrationType::New();
>
>  registration->SetMetric(        metric        );
>  registration->SetOptimizer(     optimizer     );
>  registration->SetTransform(     transform     );
>  registration->SetFixedPointSet( fixedPointSet );
>  registration->SetMovingImage(   shrinkFilter->GetOutput()   );
>  registration->SetInterpolator(  interpolator  );
>  registration->SetInitialTransformParameters( transform->GetParameters());
>
> 
-------------- next part --------------
A non-text attachment was scrubbed...
Name: two_views.jpg
Type: image/jpeg
Size: 25990 bytes
Desc: not available
URL: <http://www.itk.org/pipermail/insight-users/attachments/20110920/408b9507/attachment.jpg>