Examples/RegistrationITKv4/LandmarkWarping2.cxx

/*=========================================================================
 *
 *  Copyright NumFOCUS
 *
 *  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 deform an image using a KernelBase spline
//  and two sets of landmarks.
//
//  \index{WarpImageFilter}
//  \index{LandmarkDisplacementFieldSource}
//
//  In addition to standard headers included in previous examples, this
//  example requires the following includes:
//
//  Software Guide : EndLatex
 
#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkDisplacementFieldTransform.h"
#include "itkResampleImageFilter.h"
// Software Guide : BeginCodeSnippet
#include "itkVector.h"
#include "itkLandmarkDisplacementFieldSource.h"
#include <fstream>
// Software Guide : EndCodeSnippet
 
int
main(int argc, char * argv[])
{
 
  if (argc < 3)
  {
    std::cerr << "Missing Parameters " << std::endl;
    std::cerr << "Usage: " << argv[0];
    std::cerr << " landmarksFile fixedImage ";
    std::cerr << "movingImage deformedMovingImage" << std::endl;
    return EXIT_FAILURE;
  }
 
  constexpr unsigned int Dimension = 2;
  using VectorComponentType = float;
 
  using VectorType = itk::Vector<VectorComponentType, Dimension>;
 
  using DisplacementFieldType = itk::Image<VectorType, Dimension>;
 
  using PixelType = unsigned char;
  using FixedImageType = itk::Image<PixelType, Dimension>;
  using MovingImageType = itk::Image<PixelType, Dimension>;
 
  using FixedReaderType = itk::ImageFileReader<FixedImageType>;
  using MovingReaderType = itk::ImageFileReader<MovingImageType>;
 
  using MovingWriterType = itk::ImageFileWriter<MovingImageType>;
 
 
  FixedReaderType::Pointer fixedReader = FixedReaderType::New();
  fixedReader->SetFileName(argv[2]);
 
  try
  {
    fixedReader->Update();
  }
  catch (const itk::ExceptionObject & excp)
  {
    std::cerr << "Exception thrown " << std::endl;
    std::cerr << excp << std::endl;
    return EXIT_FAILURE;
  }
 
 
  MovingReaderType::Pointer movingReader = MovingReaderType::New();
  MovingWriterType::Pointer movingWriter = MovingWriterType::New();
 
  movingReader->SetFileName(argv[3]);
  movingWriter->SetFileName(argv[4]);
 
 
  FixedImageType::ConstPointer fixedImage = fixedReader->GetOutput();
 
  //  Software Guide : BeginLatex
  //
  //  After reading in the fixed and moving images, the \code{deformer} object
  //  is instantiated from the \code{itk::LandmarkDisplacementFieldSource}
  //  class, and parameters of the image space and orientation are set.
  //
  //  Software Guide : EndLatex
 
  // Software Guide : BeginCodeSnippet
 
  using DisplacementSourceType =
    itk::LandmarkDisplacementFieldSource<DisplacementFieldType>;
 
  DisplacementSourceType::Pointer deformer = DisplacementSourceType::New();
 
  deformer->SetOutputSpacing(fixedImage->GetSpacing());
  deformer->SetOutputOrigin(fixedImage->GetOrigin());
  deformer->SetOutputRegion(fixedImage->GetLargestPossibleRegion());
  deformer->SetOutputDirection(fixedImage->GetDirection());
 
  // Software Guide : EndCodeSnippet
 
  //  Software Guide : BeginLatex
  //
  //  Source and target landmarks are then created, and the points themselves
  //  are read in from a file stream.
  //
  //  Software Guide : EndLatex
 
  // Software Guide : BeginCodeSnippet
 
  using LandmarkContainerType = DisplacementSourceType::LandmarkContainer;
  using LandmarkPointType = DisplacementSourceType::LandmarkPointType;
 
  LandmarkContainerType::Pointer sourceLandmarks =
    LandmarkContainerType::New();
  LandmarkContainerType::Pointer targetLandmarks =
    LandmarkContainerType::New();
 
  LandmarkPointType sourcePoint;
  LandmarkPointType targetPoint;
 
  std::ifstream pointsFile;
  pointsFile.open(argv[1]);
 
  unsigned int pointId = 0;
 
  pointsFile >> sourcePoint;
  pointsFile >> targetPoint;
 
  while (!pointsFile.fail())
  {
    sourceLandmarks->InsertElement(pointId, sourcePoint);
    targetLandmarks->InsertElement(pointId, targetPoint);
    ++pointId;
 
    pointsFile >> sourcePoint;
    pointsFile >> targetPoint;
  }
 
  pointsFile.close();
 
  //  Software Guide : EndCodeSnippet
 
  //  Software Guide : BeginLatex
  //
  //  The source and target landmark objects are then assigned to
  //  \code{deformer}.
  //
  //  Software Guide : EndLatex
 
  //  Software Guide : BeginCodeSnippet
 
  deformer->SetSourceLandmarks(sourceLandmarks);
  deformer->SetTargetLandmarks(targetLandmarks);
 
  //  Software Guide : EndCodeSnippet
 
  //  Software Guide : BeginLatex
  //
  //  After calling \code{UpdateLargestPossibleRegion()} on the
  //  \code{deformer}, the displacement field may be obtained via the
  //  \code{GetOutput()} method.
  //
  //  Software Guide : EndLatex
 
  try
  {
    deformer->UpdateLargestPossibleRegion();
  }
  catch (const itk::ExceptionObject & excp)
  {
    std::cerr << "Exception thrown " << std::endl;
    std::cerr << excp << std::endl;
    return EXIT_FAILURE;
  }
 
  DisplacementFieldType::Pointer displacementField = deformer->GetOutput();
 
  using InterpolatorPrecisionType = double;
  using TransformPrecisionType = float;
  using FilterType = itk::ResampleImageFilter<MovingImageType,
                                              MovingImageType,
                                              InterpolatorPrecisionType,
                                              TransformPrecisionType>;
  FilterType::Pointer warper = FilterType::New();
 
  using InterpolatorType =
    itk::LinearInterpolateImageFunction<MovingImageType,
                                        InterpolatorPrecisionType>;
 
  InterpolatorType::Pointer interpolator = InterpolatorType::New();
 
  warper->SetInterpolator(interpolator);
 
  using DisplacementFieldTransformType =
    itk::DisplacementFieldTransform<TransformPrecisionType, Dimension>;
  auto displacementTransform = DisplacementFieldTransformType::New();
  displacementTransform->SetDisplacementField(displacementField);
  warper->SetTransform(displacementTransform);
 
  warper->SetOutputSpacing(displacementField->GetSpacing());
  warper->SetOutputOrigin(displacementField->GetOrigin());
 
  warper->SetOutputDirection(displacementField->GetDirection());
  warper->SetSize(displacementField->GetLargestPossibleRegion().GetSize());
 
  warper->SetInput(movingReader->GetOutput());
 
  movingWriter->SetInput(warper->GetOutput());
 
  try
  {
    movingWriter->Update();
  }
  catch (const itk::ExceptionObject & excp)
  {
    std::cerr << "Exception thrown " << std::endl;
    std::cerr << excp << std::endl;
    return EXIT_FAILURE;
  }
 
  return EXIT_SUCCESS;
}