[ITK-users] infinite Loop in Registration

[Michael_A]

Hello dear forum users,
 
i am quite a beginner in terms of ITK, but I have been handed the task of
putting together a registration code to co-register MRT and PET 3D-Images.
So I followed the documentation and tried to put a simple code together that
would at least be able to run and output a co-registered file (doesnt have
to be any good yet, just an output).
 
So I used the Quaternion-Transform with the respetive optimizer and a Mutual
Information Metric, Linear Interpolator and an observer (which somehow
doesnt answer anything when executing?) and put it together.
I added cout stations simply to check where loop is.
 
here is my code:
 
 
#include <iostream>
 
#include "itkImageRegistrationMethod.h"
#include "itkQuaternionRigidTransform.h"
#include "itkMutualInformationImageToImageMetric.h"
#include "itkQuaternionRigidTransformGradientDescentOptimizer.h"
#include "itkResampleImageFilter.h"
 
 
#include "itkNormalizeImageFilter.h"
#include "itkDiscreteGaussianImageFilter.h"
 
#include "itkCastImageFilter.h"
#include "itkSubtractImageFilter.h"
#include "itkRescaleIntensityImageFilter.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
 
#include "itkCommand.h"
 
#include "CommandIterationUpdate.h"
 
int main(int argc, char **argv) 
{
	const unsigned int Dimension = 3;
 
	typedef unsigned short PixelType;
  	typedef itk::Image									< PixelType, Dimension > 				FixedImageType;
  	typedef itk::Image									< PixelType, Dimension > 				MovingImageType;
  	typedef itk::Image									< PixelType, Dimension > 				ImageType;
  	typedef float InternalPixelType;
  	typedef itk::Image									< InternalPixelType, Dimension > 	
InternalImageType;
  	typedef itk::QuaternionRigidTransform						< double >					TransformType;
  	typedef itk::QuaternionRigidTransformGradientDescentOptimizer							
OptimizerType;
  	typedef itk::LinearInterpolateImageFunction			<InternalImageType,double > 		
InterpolatorType;
  	typedef itk::ImageRegistrationMethod			
<InternalImageType,InternalImageType > 	RegistrationType;
  	typedef itk::MutualInformationImageToImageMetric
<InternalImageType,InternalImageType > 	MetricType;
  	typedef itk::ImageFileReader						< FixedImageType > 					
FixedImageReaderType;
  	typedef itk::ImageFileReader						< MovingImageType > 				
MovingImageReaderType;
  	typedef itk::ResampleImageFilter					<MovingImageType,FixedImageType > 	
ResampleFilterType;
 
  	MetricType::Pointer metric = MetricType::New();
  	TransformType::Pointer transform = TransformType::New();
  	OptimizerType::Pointer optimizer = OptimizerType::New();
  	InterpolatorType::Pointer interpolator = InterpolatorType::New();
  	RegistrationType::Pointer registration = RegistrationType::New();
	
  	registration->SetMetric( metric );
  	registration->SetOptimizer( optimizer );
  	registration->SetTransform( transform );
  	registration->SetInterpolator( interpolator );
	
  	metric->SetFixedImageStandardDeviation( 0.4 );
  	metric->SetMovingImageStandardDeviation( 0.4 );
	
  	FixedImageReaderType::Pointer fixedImageReader =
FixedImageReaderType::New();
  	MovingImageReaderType::Pointer movingImageReader =
MovingImageReaderType::New();
	
  	fixedImageReader->SetFileName( argv[1] );
  	movingImageReader->SetFileName( argv[2] );
	
  	fixedImageReader->Update();
  	movingImageReader->Update();
 
std::cout << 1 << std::endl;
	ResampleFilterType::Pointer resampler = ResampleFilterType::New();
  	resampler->SetInput( movingImageReader->GetOutput() );
  
  	typedef itk::NormalizeImageFilter					<FixedImageType,InternalImageType> 	
FixedNormalizeFilterType;
  	typedef itk::NormalizeImageFilter					<MovingImageType,InternalImageType> 
MovingNormalizeFilterType;
	
  	FixedNormalizeFilterType::Pointer fixedNormalizer =
FixedNormalizeFilterType::New();
  	MovingNormalizeFilterType::Pointer movingNormalizer =
MovingNormalizeFilterType::New();
  
  	typedef itk::DiscreteGaussianImageFilter		
<InternalImageType,InternalImageType> 	GaussianFilterType;
  	GaussianFilterType::Pointer fixedSmoother = GaussianFilterType::New();
  	GaussianFilterType::Pointer movingSmoother = GaussianFilterType::New();
	
  	fixedSmoother->SetVariance( 2.0 );
  	movingSmoother->SetVariance( 2.0 );
	
  std::cout << 2 << std::endl;
    	
  	fixedNormalizer->SetInput( fixedImageReader->GetOutput() );
  	movingNormalizer->SetInput( movingImageReader->GetOutput() );
	
	fixedNormalizer->Update();
	movingNormalizer->Update();
 
  	fixedSmoother->SetInput( fixedNormalizer->GetOutput() );
  	movingSmoother->SetInput( movingNormalizer->GetOutput() );
	
	fixedSmoother->Update();
	movingSmoother->Update();
 
  	registration->SetFixedImage( fixedSmoother->GetOutput() );
  	registration->SetMovingImage( movingSmoother->GetOutput() );
	
	std::cout << 3 << std::endl;
	
	ImageType::RegionType fixedImageRegion =
fixedSmoother->GetOutput()->GetBufferedRegion();
  	registration->SetFixedImageRegion(	fixedImageRegion   );
 
  	const unsigned int numberOfPixels = fixedImageRegion.GetNumberOfPixels();
  	const unsigned int numberOfSamples = static_cast< unsigned int >(
numberOfPixels * 0.5 );
  	metric->SetNumberOfSpatialSamples( numberOfSamples );
  
  	optimizer->SetLearningRate( 10 );
  	optimizer->SetNumberOfIterations( 200 );
  	optimizer->MaximizeOn();
	
	typedef RegistrationType::ParametersType ParametersType;
	ParametersType initialParameters( transform->GetNumberOfParameters() );
	for(int i = 0; i < transform->GetNumberOfParameters(); i++)
	{
	  initialParameters[i] = 0.0; 
	}
 
	registration->SetInitialTransformParameters( initialParameters );
	//transform->SetIdentity();
	
	CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
	optimizer->AddObserver(itk::IterationEvent(), observer);
	
	std::cout << 4 << std::endl;
	
  	try
  	{
  		registration->Update();
  		std::cout << "Optimizer stop condition: " <<
registration->GetOptimizer()->GetStopConditionDescription() << std::endl;
  	}
  	catch( itk::ExceptionObject & err )
  	{
  		std::cout << "ExceptionObject caught !" << std::endl;
  		std::cout << err << std::endl;
  		return EXIT_FAILURE;
  	}
  	std::cout << 5 << std::endl;
  	transform->SetParameters(registration->GetLastTransformParameters());
 
  	// Seite 299 ff der Dokumentation
  	resampler->SetTransform( registration->GetOutput()->Get() );
  	FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();
  	resampler->SetSize( fixedImage->GetLargestPossibleRegion().GetSize() );
  	resampler->SetOutputOrigin( fixedImage->GetOrigin() );
  	resampler->SetOutputSpacing( fixedImage->GetSpacing() );
  	resampler->SetOutputDirection( fixedImage->GetDirection() );
  	resampler->SetDefaultPixelValue( 100 );
	resampler->Update();
	
std::cout << 6 << std::endl;
 
  	typedef unsigned char OutputPixelType;
  	typedef itk::Image					<OutputPixelType, Dimension > 			OutputImageType;
  	typedef itk::CastImageFilter				<FixedImageType,OutputImageType> 	
CastFilterType;
  	typedef itk::ImageFileWriter				<OutputImageType> 				WriterType;
 
  	WriterType::Pointer writer = WriterType::New();
  	CastFilterType::Pointer caster = CastFilterType::New();
 
  	caster->SetInput( resampler->GetOutput() );
	caster->Update();
  	writer->SetInput( caster->GetOutput() );
	//writer->SetFileName( (const char&) argv[1] + (const char*) "-" + (const
char&) argv[2] + "-coregistered.dcm" );
	writer->SetFileName( "111.dcm" );
	writer->Update();
	std::cout <<7  << std::endl;
	resampler->SetDefaultPixelValue( 1 );
	resampler->Update();
 
std::cout << 8 << std::endl;
 
  	typedef itk::SubtractImageFilter		
<FixedImageType,FixedImageType,FixedImageType > DifferenceFilterType;
  	DifferenceFilterType::Pointer difference = DifferenceFilterType::New();
  	difference->SetInput1( fixedImageReader->GetOutput() );
  	difference->SetInput2( resampler->GetOutput() );
	difference->Update();
	
	typedef itk::RescaleIntensityImageFilter		<FixedImageType,OutputImageType>	
RescalerType;
	RescalerType::Pointer intensityRescaler = RescalerType::New();
	intensityRescaler->SetInput(   difference->GetOutput()  );
	intensityRescaler->SetOutputMinimum( 0 );
	intensityRescaler->SetOutputMaximum( 255 );
	intensityRescaler->Update();
	
	std::cout << 9 << std::endl;
 
	WriterType::Pointer writer2 = WriterType::New();
	writer2->SetInput( intensityRescaler->GetOutput() );
	//writer2->SetFileName((const char&) argv[1] + "-" + (const char&) argv[2]
+ "-Difference.dcm");
	writer2->SetFileName("123.dcm");
	writer2->Update();
	std::cout << 10 << std::endl;
	return EXIT_SUCCESS;
}
 
 
the observer class has been put into a header file and included at the
beginning. It does the exact same as in the documentation, yet failes to
cout anything :/
 
the Test run with image data of a MRT and a PET loops at 4, meaning the
terminal says:
 
1
2
3
4
 
and then it simpy calculates to no end...
 
Does someone know why? What did I miss to implement or has something been
done crucially wrong?
 
best regards and thanks in advance,
Michael



--
View this message in context: http://itk-users.7.n7.nabble.com/infinite-Loop-in-Registration-tp33864.html
Sent from the ITK - Users mailing list archive at Nabble.com.