ITK/Examples/Registration/ImageRegistrationMethodBSpline: Difference between revisions

From KitwarePublic
< ITK‎ | Examples
Jump to navigationJump to search
(Deprecated content that is moved to sphinx)
 
(9 intermediate revisions by one other user not shown)
Line 1: Line 1:
==ImageRegistrationMethodBSpline.cxx==
{{warning|1=The media wiki content on this page is no longer maintainedThe examples presented on the https://itk.org/Wiki/* pages likely require ITK version 4.13 or earlier releasesIn many cases, the examples on this page no longer conform to the best practices for modern ITK versions.}}
<source lang="cpp">
#include "itkImageRegistrationMethod.h"
#include "itkMeanSquaresImageToImageMetric.h"
#include "itkTimeProbesCollectorBase.h"
#include "itkSpatialObjectToImageFilter.h"
#include "itkEllipseSpatialObject.h"
 
#include "itkBSplineDeformableTransform.h"
#include "itkLBFGSOptimizer.h"
#include "itkImageFileWriter.h"
#include "itkResampleImageFilter.h"
#include "itkCastImageFilter.h"
#include "itkSquaredDifferenceImageFilter.h"
 
const    unsigned int    ImageDimension = 2;
typedef float          PixelType;
 
typedef itk::Image< PixelType, ImageDimension >  ImageType;
 
static void CreateEllipseImage(ImageType::Pointer image);
static void CreateCircleImage(ImageType::Pointer image);
 
int main( int argc, char *argv[] )
{
 
  const unsigned int SpaceDimension = ImageDimension;
  const unsigned int SplineOrder = 3;
  typedef double CoordinateRepType;
 
  typedef itk::BSplineDeformableTransform<
                            CoordinateRepType,
                            SpaceDimension,
                            SplineOrder >    TransformType;
 
  typedef itk::LBFGSOptimizer      OptimizerType;
 
 
  typedef itk::MeanSquaresImageToImageMetric<
                                    ImageType,
                                    ImageType >    MetricType;
 
  typedef itk:: LinearInterpolateImageFunction<
                                    ImageType,
                                    double          >    InterpolatorType;
 
  typedef itk::ImageRegistrationMethod<
                                    ImageType,
                                    ImageType >    RegistrationType;
 
  MetricType::Pointer        metric        = MetricType::New();
  OptimizerType::Pointer      optimizer    = OptimizerType::New();
  InterpolatorType::Pointer  interpolator  = InterpolatorType::New();
  RegistrationType::Pointer  registration  = RegistrationType::New();
 
 
  registration->SetMetric(        metric        );
  registration->SetOptimizer(    optimizer    );
  registration->SetInterpolator(  interpolator  );
 
  TransformType::Pointer  transform = TransformType::New();
  registration->SetTransform( transform );
 
  // Create the synthetic images
  ImageType::Pointer  fixedImage  = ImageType::New();
  CreateCircleImage(fixedImage);
 
  ImageType::Pointer movingImage = ImageType::New();
  CreateEllipseImage(movingImage);
 
  // Write the images
  typedef itk::ImageFileWriter< ImageType >  WriterType;
/*
  WriterType::Pointer      fixedWriter =  WriterType::New();
  fixedWriter->SetFileName("fixed.png");
  fixedWriter->SetInput( fixedImage);
  fixedWriter->Update();
 
  WriterType::Pointer      movingWriter =  WriterType::New();
  movingWriter->SetFileName("moving.png");
  movingWriter->SetInput( movingImage);
  movingWriter->Update();
*/
 
  // Setup the registration
  registration->SetFixedImage( fixedImage  );
  registration->SetMovingImage(  movingImage);
 
  ImageType::RegionType fixedRegion = fixedImage->GetBufferedRegion();
 
  registration->SetFixedImageRegion( fixedRegion );
 
  //  Here we define the parameters of the BSplineDeformableTransform grid. We
  //  arbitrarily decide to use a grid with $5 \times 5$ nodes within the image.
   //  The reader should note that the BSpline computation requires a
  //  finite support region ( 1 grid node at the lower borders and 2
  //  grid nodes at upper borders). Therefore in this example, we set
  //  the grid size to be $8 \times 8$ and place the grid origin such that
  //  grid node (1,1) coincides with the first pixel in the fixed image.
 
  typedef TransformType::RegionType RegionType;
  RegionType bsplineRegion;
  RegionType::SizeType  gridSizeOnImage;
  RegionType::SizeType  gridBorderSize;
  RegionType::SizeType  totalGridSize;
 
  gridSizeOnImage.Fill( 5 );
  gridBorderSize.Fill( 3 );    // Border for spline order = 3 ( 1 lower, 2 upper )
  totalGridSize = gridSizeOnImage + gridBorderSize;
 
  bsplineRegion.SetSize( totalGridSize );
 
  typedef TransformType::SpacingType SpacingType;
  SpacingType spacing = fixedImage->GetSpacing();
 
  typedef TransformType::OriginType OriginType;
  OriginType origin = fixedImage->GetOrigin();
 
  ImageType::SizeType fixedImageSize = fixedRegion.GetSize();
 
  for(unsigned int r=0; r<ImageDimension; r++)
    {
    spacing[r] *= static_cast<double>(fixedImageSize[r] - 1)  /
                  static_cast<double>(gridSizeOnImage[r] - 1);
    }
 
  ImageType::DirectionType gridDirection = fixedImage->GetDirection();
  SpacingType gridOriginOffset = gridDirection * spacing;
 
  OriginType gridOrigin = origin - gridOriginOffset;
 
  transform->SetGridSpacing( spacing );
  transform->SetGridOrigin( gridOrigin );
  transform->SetGridRegion( bsplineRegion );
  transform->SetGridDirection( gridDirection );
 
  typedef TransformType::ParametersType    ParametersType;
 
  const unsigned int numberOfParameters =
              transform->GetNumberOfParameters();
 
  ParametersType parameters( numberOfParameters );
 
  parameters.Fill( 0.0 );
 
  transform->SetParameters( parameters );
 
  //  We now pass the parameters of the current transform as the initial
  //  parameters to be used when the registration process starts.
 
  registration->SetInitialTransformParameters( transform->GetParameters() );
 
  std::cout << "Intial Parameters = " << std::endl;
  std::cout << transform->GetParameters() << std::endl;
 
  //  Next we set the parameters of the LBFGS Optimizer.
 
  optimizer->SetGradientConvergenceTolerance( 0.05 );
  optimizer->SetLineSearchAccuracy( 0.9 );
  optimizer->SetDefaultStepLength( 1.5 );
  optimizer->TraceOn();
  optimizer->SetMaximumNumberOfFunctionEvaluations( 1000 );
 
  std::cout << std::endl << "Starting Registration" << std::endl;
 
  try
    {
    registration->StartRegistration();
    std::cout << "Optimizer stop condition = "
              << registration->GetOptimizer()->GetStopConditionDescription()
              << std::endl;
    }
  catch( itk::ExceptionObject & err )
    {
    std::cerr << "ExceptionObject caught !" << std::endl;
    std::cerr << err << std::endl;
    return EXIT_FAILURE;
    }
 
  OptimizerType::ParametersType finalParameters =
                    registration->GetLastTransformParameters();
 
  std::cout << "Last Transform Parameters" << std::endl;
  std::cout << finalParameters << std::endl;
 
  transform->SetParameters( finalParameters );
 
  typedef itk::ResampleImageFilter<
                            ImageType,
                            ImageType >    ResampleFilterType;
 
  ResampleFilterType::Pointer resample = ResampleFilterType::New();
 
  resample->SetTransform( transform );
  resample->SetInput( movingImage );
 
  resample->SetSize(    fixedImage->GetLargestPossibleRegion().GetSize() );
  resample->SetOutputOrigin(  fixedImage->GetOrigin() );
  resample->SetOutputSpacing( fixedImage->GetSpacing() );
  resample->SetOutputDirection( fixedImage->GetDirection() );
  resample->SetDefaultPixelValue( 100 );
 
  typedef  unsigned char  OutputPixelType;
 
  typedef itk::Image< OutputPixelType, ImageDimension > OutputImageType;
 
  typedef itk::CastImageFilter<
                        ImageType,
                        OutputImageType > CastFilterType;
 
  typedef itk::ImageFileWriter< OutputImageType >  OutputWriterType;
                       
  OutputWriterType::Pointer      writer =  OutputWriterType::New();
  CastFilterType::Pointer  caster =  CastFilterType::New();
 
 
  writer->SetFileName("output.png");
 
  caster->SetInput( resample->GetOutput() );
  writer->SetInput( caster->GetOutput()  );
 
  try
    {
    writer->Update();
    }
  catch( itk::ExceptionObject & err )
    {
    std::cerr << "ExceptionObject caught !" << std::endl;
    std::cerr << err << std::endl;
    return EXIT_FAILURE;
    }
 
  return EXIT_SUCCESS;
}
 
void CreateEllipseImage(ImageType::Pointer image)
{
  typedef itk::EllipseSpatialObject< ImageDimension >  EllipseType;
 
  typedef itk::SpatialObjectToImageFilter<
    EllipseType, ImageType >  SpatialObjectToImageFilterType;
 
  SpatialObjectToImageFilterType::Pointer imageFilter =
    SpatialObjectToImageFilterType::New();
 
  ImageType::SizeType size;
  size[ 0 ] =  100;
  size[ 1 ] =  100;
 
  imageFilter->SetSize( size );
 
  ImageType::SpacingType spacing;
  spacing.Fill(1);
  imageFilter->SetSpacing(spacing);
 
  EllipseType::Pointer ellipse    = EllipseType::New();
  EllipseType::ArrayType radiusArray;
  radiusArray[0] = 10;
  radiusArray[1] = 20;
  ellipse->SetRadius(radiusArray);
 
  typedef EllipseType::TransformType                TransformType;
  TransformType::Pointer transform = TransformType::New();
  transform->SetIdentity();
 
  TransformType::OutputVectorType  translation;
  TransformType::CenterType        center;
 
  translation[ 0 ] =  65;
  translation[ 1 ] =  45;
  transform->Translate( translation, false );
 
  ellipse->SetObjectToParentTransform( transform );
 
  imageFilter->SetInput(ellipse);
 
  ellipse->SetDefaultInsideValue(255);
  ellipse->SetDefaultOutsideValue(0);
  imageFilter->SetUseObjectValue( true );
  imageFilter->SetOutsideValue( 0 );
 
  imageFilter->Update();
 
  image->Graft(imageFilter->GetOutput());
 
}
 
void CreateCircleImage(ImageType::Pointer image)
{
typedef itk::EllipseSpatialObject< ImageDimension >  EllipseType;
 
  typedef itk::SpatialObjectToImageFilter<
    EllipseType, ImageType >  SpatialObjectToImageFilterType;
 
  SpatialObjectToImageFilterType::Pointer imageFilter =
    SpatialObjectToImageFilterType::New();
 
  ImageType::SizeType size;
  size[ 0 ] =  100;
  size[ 1 ] =  100;
 
  imageFilter->SetSize( size );
 
  ImageType::SpacingType spacing;
  spacing.Fill(1);
  imageFilter->SetSpacing(spacing);
 
  EllipseType::Pointer ellipse    = EllipseType::New();
  EllipseType::ArrayType radiusArray;
  radiusArray[0] = 10;
  radiusArray[1] = 10;
  ellipse->SetRadius(radiusArray);
 
  typedef EllipseType::TransformType                TransformType;
  TransformType::Pointer transform = TransformType::New();
  transform->SetIdentity();
 
  TransformType::OutputVectorType  translation;
  TransformType::CenterType        center;
 
  translation[ 0 ] =  50;
  translation[ 1 ] =  50;
  transform->Translate( translation, false );
 
  ellipse->SetObjectToParentTransform( transform );
 
  imageFilter->SetInput(ellipse);
 
  ellipse->SetDefaultInsideValue(255);
  ellipse->SetDefaultOutsideValue(0);
  imageFilter->SetUseObjectValue( true );
  imageFilter->SetOutsideValue( 0 );
 
  imageFilter->Update();
 
  image->Graft(imageFilter->GetOutput());
}
</source>
 
==CMakeLists.txt==
<source lang="cmake">
cmake_minimum_required(VERSION 2.6)
 
PROJECT(BSpline)
 
FIND_PACKAGE(ITK REQUIRED)
INCLUDE(${ITK_USE_FILE})
 
ADD_EXECUTABLE(BSpline BSpline.cxx)
TARGET_LINK_LIBRARIES(BSpline ITKIO ITKNumerics)
 
 
</source>

Latest revision as of 14:57, 7 June 2019

Warning: The media wiki content on this page is no longer maintained. The examples presented on the https://itk.org/Wiki/* pages likely require ITK version 4.13 or earlier releases. In many cases, the examples on this page no longer conform to the best practices for modern ITK versions.