Examples/RegistrationITKv4/IterativeClosestPoint2.cxx

/*=========================================================================
 *
 *  Copyright Insight Software Consortium
 *
 *  Licensed under the Apache License, Version 2.0 (the "License");
 *  you may not use this file except in compliance with the License.
 *  You may obtain a copy of the License at
 *
 *         http://www.apache.org/licenses/LICENSE-2.0.txt
 *
 *  Unless required by applicable law or agreed to in writing, software
 *  distributed under the License is distributed on an "AS IS" BASIS,
 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *  See the License for the specific language governing permissions and
 *  limitations under the License.
 *
 *=========================================================================*/

// Software Guide : BeginLatex
//
// This example illustrates how to perform Iterative Closest Point (ICP)
// registration in ITK using sets of 3D points.
//
// Software Guide : EndLatex

// Software Guide : BeginLatex
//
// The first step is to include the relevant headers.
//
// Software Guide : EndLatex

// Software Guide : BeginCodeSnippet
#include "itkEuler3DTransform.h"
#include "itkEuclideanDistancePointMetric.h"
#include "itkLevenbergMarquardtOptimizer.h"
#include "itkPointSetToPointSetRegistrationMethod.h"
#include <iostream>
#include <fstream>
// Software Guide : EndCodeSnippet

class CommandIterationUpdate : public itk::Command
{
public:
  using Self = CommandIterationUpdate;
  using Superclass = itk::Command;
  using Pointer = itk::SmartPointer<Self>;
  itkNewMacro( Self );

protected:
  CommandIterationUpdate() = default;

public:
  using OptimizerType = itk::LevenbergMarquardtOptimizer;
  using OptimizerPointer = const OptimizerType *;

  void Execute(itk::Object *caller, const itk::EventObject & event) override
    {
    Execute( (const itk::Object *)caller, event);
    }

  void Execute(const itk::Object * object, const itk::EventObject & event) override
    {
    auto optimizer = dynamic_cast< OptimizerPointer >( object );
    if( optimizer == nullptr )
      {
      itkExceptionMacro( "Could not cast optimizer." );
      }

    if( ! itk::IterationEvent().CheckEvent( &event ) )
      {
      return;
      }

    std::cout << "Value = " << optimizer->GetCachedValue() << std::endl;
    std::cout << "Position = "  << optimizer->GetCachedCurrentPosition();
    std::cout << std::endl << std::endl;
    }
};


int main(int argc, char * argv[] )
{

  if( argc < 3 )
    {
    std::cerr << "Arguments Missing. " << std::endl;
    std::cerr <<
      "Usage:  IterativeClosestPoint2   fixedPointsFile  movingPointsFile "
      << std::endl;
    return EXIT_FAILURE;
    }

  constexpr unsigned int Dimension = 3;

// Software Guide : BeginLatex
//
// First, define the necessary types for the moving and fixed point sets.
//
// Software Guide : EndLatex

// Software Guide : BeginCodeSnippet
  using PointSetType = itk::PointSet< float, Dimension >;

  PointSetType::Pointer fixedPointSet  = PointSetType::New();
  PointSetType::Pointer movingPointSet = PointSetType::New();

  using PointType = PointSetType::PointType;

  using PointsContainer = PointSetType::PointsContainer;

  PointsContainer::Pointer fixedPointContainer  = PointsContainer::New();
  PointsContainer::Pointer movingPointContainer = PointsContainer::New();

  PointType fixedPoint;
  PointType movingPoint;
// Software Guide : EndCodeSnippet

  // Read the file containing coordinates of fixed points.
  std::ifstream   fixedFile;
  fixedFile.open( argv[1] );
  if( fixedFile.fail() )
    {
    std::cerr << "Error opening points file with name : " << std::endl;
    std::cerr << argv[1] << std::endl;
    return EXIT_FAILURE;
    }

  unsigned int pointId = 0;
  fixedFile >> fixedPoint;
  while( !fixedFile.eof() )
    {
    fixedPointContainer->InsertElement( pointId, fixedPoint );
    fixedFile >> fixedPoint;
    pointId++;
    }
  fixedPointSet->SetPoints( fixedPointContainer );
  std::cout <<
    "Number of fixed Points = " << fixedPointSet->GetNumberOfPoints()
    << std::endl;

  // Read the file containing coordinates of moving points.
  std::ifstream   movingFile;
  movingFile.open( argv[2] );
  if( movingFile.fail() )
    {
    std::cerr << "Error opening points file with name : " << std::endl;
    std::cerr << argv[2] << std::endl;
    return EXIT_FAILURE;
    }

  pointId = 0;
  movingFile >> movingPoint;
  while( !movingFile.eof() )
    {
    movingPointContainer->InsertElement( pointId, movingPoint );
    movingFile >> movingPoint;
    pointId++;
    }
  movingPointSet->SetPoints( movingPointContainer );
  std::cout <<
    "Number of moving Points = "
    << movingPointSet->GetNumberOfPoints() << std::endl;


// Software Guide : BeginLatex
//
// After the points are read in from files, setup the metric to be used
// later by the registration.
//
// Software Guide : EndLatex

// Software Guide : BeginCodeSnippet
  using MetricType = itk::EuclideanDistancePointMetric<
                          PointSetType, PointSetType >;

  MetricType::Pointer  metric = MetricType::New();
// Software Guide : EndCodeSnippet

// Software Guide : BeginLatex
//
// Next, setup the tranform, optimizers, and registration.
//
// Software Guide : EndLatex

// Software Guide : BeginCodeSnippet
  using TransformType = itk::Euler3DTransform< double >;

  TransformType::Pointer transform = TransformType::New();


  // Optimizer Type
  using OptimizerType = itk::LevenbergMarquardtOptimizer;

  OptimizerType::Pointer      optimizer     = OptimizerType::New();
  optimizer->SetUseCostFunctionGradient(false);

  // Registration Method
  using RegistrationType = itk::PointSetToPointSetRegistrationMethod<
                            PointSetType, PointSetType >;


  RegistrationType::Pointer   registration  = RegistrationType::New();
// Software Guide : EndCodeSnippet

// Software Guide : BeginLatex
//
// Scale the translation components of the Transform in the Optimizer
//
// Software Guide : EndLatex

// Software Guide : BeginCodeSnippet
  OptimizerType::ScalesType scales( transform->GetNumberOfParameters() );
// Software Guide : EndCodeSnippet

// Software Guide : BeginLatex
//
// Next, set the scales and ranges for translations and rotations in the
// transform. Also, set the convergence criteria and number of iterations
// to be used by the optimizer.
//
// Software Guide : EndLatex

// Software Guide : BeginCodeSnippet
  constexpr double translationScale = 1000.0; // dynamic range of translations
  constexpr double rotationScale = 1.0;       // dynamic range of rotations

  scales[0] = 1.0 / rotationScale;
  scales[1] = 1.0 / rotationScale;
  scales[2] = 1.0 / rotationScale;
  scales[3] = 1.0 / translationScale;
  scales[4] = 1.0 / translationScale;
  scales[5] = 1.0 / translationScale;

  unsigned long   numberOfIterations =  2000;
  double          gradientTolerance  =  1e-4;   // convergence criterion
  double          valueTolerance     =  1e-4;   // convergence criterion
  double          epsilonFunction    =  1e-5;   // convergence criterion


  optimizer->SetScales( scales );
  optimizer->SetNumberOfIterations( numberOfIterations );
  optimizer->SetValueTolerance( valueTolerance );
  optimizer->SetGradientTolerance( gradientTolerance );
  optimizer->SetEpsilonFunction( epsilonFunction );
// Software Guide : EndCodeSnippet

// Software Guide : BeginLatex
//
// Here we start with an identity transform, although the user will usually
// be able to provide a better guess than this.
//
// Software Guide : EndLatex

// Software Guide : BeginCodeSnippet
  transform->SetIdentity();
// Software Guide : EndCodeSnippet

  registration->SetInitialTransformParameters( transform->GetParameters() );

// Software Guide : BeginLatex
//
// Connect all the components required for the registration.
//
// Software Guide : EndLatex

// Software Guide : BeginCodeSnippet
  registration->SetMetric(        metric        );
  registration->SetOptimizer(     optimizer     );
  registration->SetTransform(     transform     );
  registration->SetFixedPointSet( fixedPointSet );
  registration->SetMovingPointSet(   movingPointSet   );
// Software Guide : EndCodeSnippet
//
  // Connect an observer
  CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
  optimizer->AddObserver( itk::IterationEvent(), observer );

  try
    {
    registration->Update();
    }
  catch( itk::ExceptionObject & e )
    {
    std::cerr << e << std::endl;
    return EXIT_FAILURE;
    }

  std::cout << "Solution = " << transform->GetParameters() << std::endl;
  std::cout << "Stopping condition: " << optimizer->GetStopConditionDescription() << std::endl;

  return EXIT_SUCCESS;
}