ITK  4.13.0
Insight Segmentation and Registration Toolkit
SphinxExamples/src/Registration/Metricsv4/PerformRegistrationOnVectorImages/Code.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.
*
*=========================================================================*/
#include "itkCommand.h"
#include <iomanip>
int main(int argc, char *argv[])
{
if( argc < 4 )
{
std::cerr << "Missing Parameters " << std::endl;
std::cerr << "Usage: " << argv[0];
std::cerr << " fixedImageFile movingImageFile ";
std::cerr << " resampledMovingImageFile ";
std::cerr << " [numberOfAffineIterations numberOfDisplacementIterations] ";
std::cerr << std::endl;
return EXIT_FAILURE;
}
const char * fixedImageFile = argv[1];
const char * movingImageFile = argv[2];
const char * resampledMovingImageFile = argv[3];
unsigned int numberOfAffineIterations = 2;
unsigned int numberOfDisplacementIterations = 2;
if( argc >= 5 )
{
numberOfAffineIterations = atoi( argv[4] );
}
if( argc >= 6 )
{
numberOfDisplacementIterations = atoi( argv[5] );
}
const unsigned int Dimension = 2;
// The input images have red, blue, and green pixel components.
const unsigned int NumberOfPixelComponents = 3;
typedef float PixelComponentType;
typedef double ParametersValueType;
ReaderType::Pointer fixedImageReader = ReaderType::New();
fixedImageReader->SetFileName( fixedImageFile );
ReaderType::Pointer movingImageReader = ReaderType::New();
movingImageReader->SetFileName( movingImageFile );
// Get the input images
try
{
fixedImageReader->Update();
movingImageReader->Update();
}
catch( itk::ExceptionObject & error )
{
std::cerr << "Error while reading images: " << error << std::endl;
return EXIT_FAILURE;
}
InputImageType::Pointer fixedImage = fixedImageReader->GetOutput();
InputImageType::Pointer movingImage = movingImageReader->GetOutput();
AffineTransformType;
AffineTransformType::Pointer affineTransform = AffineTransformType::New();
ParametersValueType, Dimension > DisplacementTransformType;
DisplacementTransformType::Pointer displacementTransform =
DisplacementTransformType::New();
typedef DisplacementTransformType::DisplacementFieldType
DisplacementFieldType;
DisplacementFieldType::Pointer displacementField =
DisplacementFieldType::New();
// Set the displacement field to be the same as the fixed image region, which will
// act by default as the virtual domain in this example.
displacementField->SetRegions( fixedImage->GetLargestPossibleRegion() );
// make sure the displacement field has the same spatial information as the image
displacementField->CopyInformation( fixedImage );
std::cout << "fixedImage->GetLargestPossibleRegion(): "
<< fixedImage->GetLargestPossibleRegion() << std::endl;
displacementField->Allocate();
// Fill it with 0's
DisplacementTransformType::OutputVectorType zeroVector;
zeroVector.Fill( 0.0 );
displacementField->FillBuffer( zeroVector );
// Assign to transform
displacementTransform->SetDisplacementField( displacementField );
displacementTransform->SetGaussianSmoothingVarianceForTheUpdateField( 5.0 );
displacementTransform->SetGaussianSmoothingVarianceForTheTotalField( 6.0 );
// Identity transform for fixed image
IdentityTransformType;
IdentityTransformType::Pointer identityTransform =
IdentityTransformType::New();
// The metric
VirtualImageType;
InputImageType, InputImageType, VirtualImageType, InputPixelType::Length >
MetricTraitsType;
InputImageType, InputImageType, VirtualImageType, double, MetricTraitsType >
MetricType;
typedef MetricType::FixedSampledPointSetType
PointSetType;
MetricType::Pointer metric = MetricType::New();
PointSetType::Pointer pointSet = PointSetType::New();
fixedImage->GetLargestPossibleRegion() );
itk::SizeValueType index = 0;
itk::SizeValueType count = 0;
for( fixedImageIt.GoToBegin(); !fixedImageIt.IsAtEnd(); ++fixedImageIt )
{
// take every N^th point
if ( count % 2 == 0 )
{
PointType point;
fixedImage->TransformIndexToPhysicalPoint( fixedImageIt.GetIndex(), point );
pointSet->SetPoint( index, point );
++index;
}
++count;
}
metric->SetFixedSampledPointSet( pointSet );
metric->SetUseFixedSampledPointSet( true );
// Assign images and transforms.
// By not setting a virtual domain image or virtual domain settings,
// the metric will use the fixed image for the virtual domain.
metric->SetFixedImage( fixedImage );
metric->SetMovingImage( movingImage );
metric->SetFixedTransform( identityTransform );
metric->SetMovingTransform( affineTransform );
const bool gaussian = false;
metric->SetUseMovingImageGradientFilter( gaussian );
metric->SetUseFixedImageGradientFilter( gaussian );
metric->Initialize();
RegistrationParameterScalesFromShiftType;
RegistrationParameterScalesFromShiftType::Pointer shiftScaleEstimator =
RegistrationParameterScalesFromShiftType::New();
shiftScaleEstimator->SetMetric( metric );
//
// Affine registration
//
std::cout << "First do an affine registration..." << std::endl;
typedef itk::GradientDescentOptimizerv4 OptimizerType;
OptimizerType::Pointer optimizer = OptimizerType::New();
optimizer->SetMetric( metric );
optimizer->SetNumberOfIterations( numberOfAffineIterations );
optimizer->SetScalesEstimator( shiftScaleEstimator );
optimizer->StartOptimization();
std::cout << "After optimization affine parameters are: "
<< affineTransform->GetParameters() << std::endl;
//
// Deformable registration
//
std::cout << "Now, add the displacement field to the composite transform..."
<< std::endl;
CompositeType::Pointer compositeTransform = CompositeType::New();
compositeTransform->AddTransform( affineTransform );
compositeTransform->AddTransform( displacementTransform );
// Optimize only the displacement field, but still using the
// previously computed affine transformation.
compositeTransform->SetOnlyMostRecentTransformToOptimizeOn();
metric->SetMovingTransform( compositeTransform );
metric->SetUseFixedSampledPointSet( false );
metric->Initialize();
// Optimizer
optimizer->SetScalesEstimator( shiftScaleEstimator );
optimizer->SetMetric( metric );
optimizer->SetNumberOfIterations( numberOfDisplacementIterations );
try
{
if( numberOfDisplacementIterations > 0 )
{
optimizer->StartOptimization();
}
else
{
std::cout << "*** Skipping displacement field optimization.\n";
}
}
catch( itk::ExceptionObject & error )
{
std::cerr << "Error: " << error << std::endl;
return EXIT_FAILURE;
}
// Warp the image with the displacement field
ResampleFilterType;
ResampleFilterType::Pointer resampler = ResampleFilterType::New();
resampler->SetTransform( compositeTransform );
resampler->SetInput( movingImageReader->GetOutput() );
resampler->SetReferenceImage( fixedImage );
resampler->SetUseReferenceImage( true );
// Write the warped image into a file
OutputPixelType;
OutputImageType;
CastFilterType;
CastFilterType::Pointer caster = CastFilterType::New();
caster->SetInput( resampler->GetOutput() );
WriterType::Pointer writer = WriterType::New();
writer->SetFileName( resampledMovingImageFile );
writer->SetInput( caster->GetOutput() );
try
{
writer->UpdateLargestPossibleRegion();
}
catch( itk::ExceptionObject & error )
{
std::cerr << "Error: " << error << std::endl;
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}