I am trying to read in two different sets of DICOM image series in two different folders and perform a registration. The series seems to be read in correctly, everything goes smoothly until registration->Update() is being called. It crushes directly and likely an abort function is called inside Update(). The program was working fine with 2D images. Here is the code
#include "itkImageRegistrationMethod.h"
#include "itkMeanSquaresImageToImageMetric.h"
#include "itkTimeProbesCollectorBase.h"
#include "itkMemoryProbesCollectorBase.h"
#include "itkBSplineTransform.h"
#include "itkLBFGSBOptimizer.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkResampleImageFilter.h"
#include "itkCastImageFilter.h"
#include "itkSquaredDifferenceImageFilter.h"
#include "itkGDCMImageIO.h"
#include "itkGDCMSeriesFileNames.h"
#include "itkNumericSeriesFileNames.h"
#include "gdcmUIDGenerator.h"
#include "itkImageSeriesReader.h"
#include "itkImageSeriesWriter.h"
#define SRI
#include "itkCommand.h"
class CommandIterationUpdate : public itk::Command
{
public:
typedef CommandIterationUpdate Self;
typedef itk::Command Superclass;
typedef itk::SmartPointer<Self> Pointer;
itkNewMacro( Self );
protected:
CommandIterationUpdate() {};
public:
typedef itk::LBFGSBOptimizer OptimizerType;
typedef const OptimizerType * OptimizerPointer;
void Execute(itk::Object *caller, const itk::EventObject & event)
{
Execute( (const itk::Object *)caller, event);
}
void Execute(const itk::Object * object, const itk::EventObject & event)
{
OptimizerPointer optimizer =
dynamic_cast< OptimizerPointer >( object );
if( !(itk::IterationEvent().CheckEvent( &event )) )
{
return;
}
std::cout << optimizer->GetCurrentIteration() << " ";
std::cout << optimizer->GetValue() << " ";
std::cout << optimizer->GetInfinityNormOfProjectedGradient() << std::endl;
}
};
int main( int argc, char *argv[] )
{
if( argc < 4 )
{
std::cerr << "Missing Parameters " << std::endl;
std::cerr << "Usage: " << argv[0];
std::cerr << " fixedImageFile movingImageFile outputImagefile ";
std::cerr << " [differenceOutputfile] [differenceBeforeRegistration] ";
std::cerr << " [deformationField] ";
return EXIT_FAILURE;
}
const unsigned int ImageDimension = 3;
typedef float PixelType;
typedef itk::Image< PixelType, ImageDimension > FixedImageType;
typedef itk::Image< PixelType, ImageDimension > MovingImageType;
const unsigned int SpaceDimension = ImageDimension;
const unsigned int SplineOrder = 3;
typedef double CoordinateRepType;
typedef itk::BSplineTransform<
CoordinateRepType,
SpaceDimension,
SplineOrder > TransformType;
typedef itk::LBFGSBOptimizer OptimizerType;
typedef itk::MeanSquaresImageToImageMetric<
FixedImageType,
MovingImageType > MetricType;
typedef itk:: LinearInterpolateImageFunction<
MovingImageType,
double > InterpolatorType;
typedef itk::ImageRegistrationMethod<
FixedImageType,
MovingImageType > 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 );
#ifdef SRI
typedef itk::ImageSeriesReader< FixedImageType > FixedImageReaderType;
typedef itk::ImageSeriesReader< MovingImageType > MovingImageReaderType;
typedef itk::GDCMImageIO ImageIOType;
typedef itk::GDCMSeriesFileNames InputNamesGeneratorType;
typedef itk::NumericSeriesFileNames OutputNameGeneratorType;
typedef itk::ImageSeriesWriter< MovingImageType, MovingImageType > SeriesWriterType;
ImageIOType::Pointer gdcmIO = ImageIOType::New();
InputNamesGeneratorType::Pointer fixedImageNames = InputNamesGeneratorType::New();
InputNamesGeneratorType::Pointer movingImageNames = InputNamesGeneratorType::New();
fixedImageNames->SetInputDirectory( argv[1] );
movingImageNames->SetInputDirectory( argv[2] );
const FixedImageReaderType::FileNamesContainer & fixedNames = fixedImageNames- >GetInputFileNames();
const MovingImageReaderType::FileNamesContainer & movingNames = movingImageNames->GetInputFileNames();
#else
typedef itk::ImageFileReader< FixedImageType > FixedImageReaderType;
typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;
#endif
FixedImageReaderType::Pointer fixedImageReader = FixedImageReaderType::New();
MovingImageReaderType::Pointer movingImageReader = MovingImageReaderType::New();
#ifdef SRI
fixedImageReader->SetImageIO( gdcmIO );
movingImageReader->SetImageIO( gdcmIO );
fixedImageReader->SetFileNames( fixedNames );
movingImageReader->SetFileNames( movingNames );
#else
fixedImageReader->SetFileName( argv[1] );
movingImageReader->SetFileName( argv[2] );
#endif
FixedImageType::ConstPointer fixedImage = fixedImageReader->GetOutput();
registration->SetFixedImage( fixedImage );
registration->SetMovingImage( movingImageReader->GetOutput() );
fixedImageReader->Update();
FixedImageType::RegionType fixedRegion = fixedImage->GetBufferedRegion();
registration->SetFixedImageRegion( fixedRegion );
typedef TransformType::RegionType RegionType;
unsigned int numberOfGridNodes = 8;
TransformType::PhysicalDimensionsType fixedPhysicalDimensions;
TransformType::MeshSizeType meshSize;
TransformType::OriginType fixedOrigin;
for( unsigned int i=0; i< SpaceDimension; i++ )
{
fixedOrigin = fixedImage->GetOrigin()[i];
fixedPhysicalDimensions[i] = fixedImage->GetSpacing()[i] *
static_cast<double>(
fixedImage->GetLargestPossibleRegion().GetSize()[i] - 1 );
}
meshSize.Fill( numberOfGridNodes - SplineOrder );
transform->SetTransformDomainOrigin( fixedOrigin );
transform->SetTransformDomainPhysicalDimensions(
fixedPhysicalDimensions );
transform->SetTransformDomainMeshSize( meshSize );
transform->SetTransformDomainDirection( fixedImage->GetDirection() );
typedef TransformType::ParametersType ParametersType;
const unsigned int numberOfParameters =
transform->GetNumberOfParameters();
ParametersType parameters( numberOfParameters );
parameters.Fill( 0.0 );
transform->SetParameters( parameters );
registration->SetInitialTransformParameters( transform->GetParameters() );
std::cout << "Intial Parameters = " << std::endl;
std::cout << transform->GetParameters() << std::endl;
OptimizerType::BoundSelectionType boundSelect( transform->GetNumberOfParameters() );
OptimizerType::BoundValueType upperBound( transform->GetNumberOfParameters() );
OptimizerType::BoundValueType lowerBound( transform->GetNumberOfParameters() );
boundSelect.Fill( 0 );
upperBound.Fill( 0.0 );
lowerBound.Fill( 0.0 );
optimizer->SetBoundSelection( boundSelect );
optimizer->SetUpperBound( upperBound );
optimizer->SetLowerBound( lowerBound );
optimizer->SetCostFunctionConvergenceFactor( 1e+12 );
optimizer->SetProjectedGradientTolerance( 1.0 );
optimizer->SetMaximumNumberOfIterations( 500 );
optimizer->SetMaximumNumberOfEvaluations( 500 );
optimizer->SetMaximumNumberOfCorrections( 5 );
CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
optimizer->AddObserver( itk::IterationEvent(), observer );
itk::TimeProbesCollectorBase chronometer;
itk::MemoryProbesCollectorBase memorymeter;
std::cout << std::endl << "Starting Registration" << std::endl;
try
{
memorymeter.Start( "Registration" );
chronometer.Start( "Registration" );
registration->Update();
chronometer.Stop( "Registration" );
memorymeter.Stop( "Registration" );
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;
// Report the time taken by the registration
chronometer.Report( std::cout );
memorymeter.Report( std::cout );
transform->SetParameters( finalParameters );
// resample and output
I have been struggling on this for weeks and still couldn't figure what is the problem. I tried to look up in the user guide and examples but none of them is reading DICOM image series.
It would also be helpful if anyone could provide me with an example of ITK registration on image series.
Thanks in advance.
