Skip to content

Instantly share code, notes, and snippets.

@thewtex
Forked from anonymous/Registration.hpp
Last active August 29, 2015 13:56
Show Gist options
  • Save thewtex/8984576 to your computer and use it in GitHub Desktop.
Save thewtex/8984576 to your computer and use it in GitHub Desktop.
cmake_minimum_required( VERSION 2.8 )
project( Registration )
find_package( ITK REQUIRED )
include( ${ITK_USE_FILE} )
include_directories( ${CMAKE_CURRENT_SOURCE_DIR} )
add_executable( ${PROJECT_NAME} main.cpp )
target_link_libraries( ${PROJECT_NAME} ${ITK_LIBRARIES} )
#include "Registration.hpp"
#include "itkImageFileReader.h"
int main( int argc, char * argv [] )
{
if( argc < 2 )
{
std::cerr << "Usage: " << argv[0] << " inputImage" << std::endl;
}
const char * inputImageFileName = argv[1];
typedef float PixelType;
const unsigned int Dimension = 2;
typedef itk::Image< PixelType, Dimension > ImageType;
typedef itk::ImageFileReader< ImageType > ReaderType;
ReaderType::Pointer reader = ReaderType::New();
reader->SetFileName( inputImageFileName );
try
{
reader->Update();
}
catch( itk::ExceptionObject & err )
{
std::cerr << "Error during reading: " << err << std::endl;
return EXIT_FAILURE;
}
ImageType::Pointer inputImage = reader->GetOutput();
/*
ImageType is a 3D double image
input image is one such image previously read from a netCDF file
*/
typedef itk::RegistrationFilter< ImageType > RegistrationFilterType;
RegistrationFilterType::Pointer registration = RegistrationFilterType::New();
// this is just a test: output transform should be identity
registration->SetFixedImage( inputImage );
registration->SetMovingImage( inputImage );
registration->ObserveRegistration( true );
try
{
registration->Update();
}
catch( itk::ExceptionObject & err )
{
std::cerr << "Error: " << err << std::endl;
}
return 0;
}
#ifndef _Registration_hpp
#define _Registration_hpp
#include <iostream>
//#include "Types.hpp"
#include <itkImageRegistrationMethod.h>
#include <itkTranslationTransform.h>
#include <itkMeanSquaresImageToImageMetric.h>
#include <itkLinearInterpolateImageFunction.h>
#include <itkRegularStepGradientDescentOptimizer.h>
#include <itkResampleImageFilter.h>
#include <itkCommand.h>
namespace itk
{
class CommandIteration : public itk::Command {
public:
typedef CommandIteration Self;
typedef itk::Command SuperClass;
typedef itk::SmartPointer< Self > Pointer;
itkNewMacro( Self );
protected:
CommandIteration()
{
}
public:
typedef itk::RegularStepGradientDescentOptimizer OptimizerType;
typedef const OptimizerType * OptimizerPointer;
void Execute( itk::Object * caller, const itk::EventObject & event )
{
this-> Execute( (const itk::Object *) caller, event );
}
void Execute( const itk::Object * caller, const itk::EventObject & event )
{ OptimizerPointer optimizer = dynamic_cast< OptimizerPointer >( caller );
if( typeid( event ) == typeid( itk::IterationEvent ) )
{
std::cout << optimizer->GetCurrentIteration() << " : ";
std::cout << optimizer->GetValue() << " : ";
std::cout << optimizer->GetCurrentPosition() << std::endl;
}
}
};
template< class T>
class RegistrationFilter : public ImageToImageFilter< T, T >
{
private:
const static unsigned int Dimension = T::ImageDimension;
typedef T ImageType;
typedef itk::TranslationTransform< double, T::ImageDimension > TransformType;
typedef itk::RegularStepGradientDescentOptimizer OptimizerType;
typedef itk::LinearInterpolateImageFunction< ImageType , double > InterpolatorType;
typedef itk::MeanSquaresImageToImageMetric< ImageType , ImageType > MetricType;
typedef itk::ImageRegistrationMethod< ImageType , ImageType > RegistrationType;
typedef typename RegistrationType::ParametersType ParametersType;
typedef itk::ResampleImageFilter< ImageType , ImageType > ResamplerType;
public:
/** Standard class typedefs. */
typedef RegistrationFilter Self;
typedef ImageToImageFilter< T, T > Superclass;
typedef SmartPointer< Self > Pointer;
typedef SmartPointer< const Self > ConstPointer;
/** Method for creation through the object factory. */
itkNewMacro(Self);
/** Run-time type information (and related methods). */
itkTypeMacro(RegistrationFilter, ImageToImageFilter);
/** Parameter personalization */
void SetFixedImage(typename T::Pointer fi)
{
_fixedImage = fi;
}
void SetMovingImage(typename T::Pointer mi)
{
_movingImage = mi;
}
void ObserveRegistration(bool obr)
{
_observeRegistration = obr;
}
protected:
RegistrationFilter()
{
_transform = TransformType::New();
_transform->SetIdentity();
_optimizer = OptimizerType::New();
_optimizer->SetMaximumStepLength( 4.00 );
_optimizer->SetMinimumStepLength( 0.01 );
_optimizer->SetNumberOfIterations( 100 );
_optimizer->MaximizeOff();
_observeRegistration = true;
_observer = CommandIteration::New();
_interpolator = InterpolatorType::New();
_metric = MetricType::New();
_registrator = RegistrationType::New();
_registrator->SetTransform( _transform );
_registrator->SetOptimizer( _optimizer );
_registrator->SetInterpolator( _interpolator );
_registrator->SetMetric( _metric );
_registrator->SetInitialTransformParameters( _transform->GetParameters() );
_resampler = ResamplerType::New();
_resampler->SetTransform ( _transform );
}
// Does the real work.
virtual void GenerateData()
{
std::cout<<" --- BEGIN REGISTRATION --- "<<std::endl;
std::cout<<"Setting parameters"<<std::endl;
_registrator->SetFixedImage( _fixedImage );
_registrator->SetMovingImage( _movingImage );
_registrator->SetFixedImageRegion( _fixedImage->GetLargestPossibleRegion() );
_resampler->SetInput( _movingImage );
_resampler->SetOutputOrigin( _fixedImage->GetOrigin() );
_resampler->SetOutputSpacing( _fixedImage->GetSpacing() );
_resampler->SetSize( _fixedImage->GetLargestPossibleRegion().GetSize() );
if ( _observeRegistration )
_optimizer->AddObserver( itk::IterationEvent(), _observer);
try
{
std::cout<<"Optimizing transform"<<std::endl;
_registrator->Update();
}
catch( itk::ExceptionObject & excp )
{
std::cerr << "Error in registration" << std::endl;
std::cerr << excp << std::endl;
}
_transform->SetParameters( _registrator->GetLastTransformParameters() );
std::cout<<"Resampling moving image from optimized transform"<<std::endl;
_resampler->Update();
std::cout<<" --- END REGISTRATION --- "<<std::endl;
this->GraftOutput( _fixedImage );
}
RegistrationFilter(const Self &); //purposely not implemented
void operator=(const Self &); //purposely not implemented
//Component filters
typename TransformType::Pointer _transform;
typename OptimizerType::Pointer _optimizer;
typename InterpolatorType::Pointer _interpolator;
typename MetricType::Pointer _metric;
typename RegistrationType::Pointer _registrator;
typename ResamplerType::Pointer _resampler;
typename CommandIteration::Pointer _observer;
typename T::Pointer _fixedImage, _movingImage;
bool _observeRegistration;
};
} //namespace ITK
#endif /* defined(____Registration__) */
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment