[Insight-users] On MultiResolutionImageRegistrationMethod
Lydia Ng
lng@insightful.com
Mon, 3 Feb 2003 14:00:44 -0800
Sateesh,
In what way was the MultiResImageRegistration2.cxx not working? =
Crashing, throwing exception or producing the expected results?
In your earlier email, you had a test case which had a 2 degree =
rotation. Did the example work for that case?
Is the code in your email all of your program or just a snippet?
I believe that the MultiResImageRegistration2.cxx example also has a =
Command (RegistrationInterfaceCommand) that adjusts the some of the =
optimization parameters at each resolution level. Tuning of the =
optimization parameters is very important to the success of the =
registration and it seems to be missing in your code.
- Lydia
-----Original Message-----
From: Sri Valli [mailto:valli_gummadi@yahoo.co.in]=20
Sent: Saturday, February 01, 2003 5:42 AM
To: insight-users@public.kitware.com
Subject: [Insight-users] On MultiResolutionImageRegistrationMethod
Hi Lydia,
=A0=A0=A0=A0 Thanks for ur suggestions. Well i could correct the problem =
which was coming before. Now again i face problem on rotation. This =
example MultiResImageRegistration2.cxx is not working for registering =
with a rotated image. I have pasted my code at the end of this mail. =
Please go through the code and send me the code which have to be added =
to it or to be modified. I am using a data set of (256 x 256 x 1). The =
pixel values are from 0 to 120. I have rotated the image in z-direction =
with angle 10 degrees. I gave this rotated image as moving image and =
normal image as fixed image to the application with code below. Please =
go through it and send me the missing code.
=A0=A0=A0 Thanks in advance.
-Regards,
=A0=A0 Sateesh.
=A0
My code......Follows.....
=A0 if( argc < 3 )
=A0=A0=A0 {
=A0=A0=A0 std::cerr << "Missing Parameters " << std::endl;
=A0=A0=A0 std::cerr << "Usage: " << argv[0];
=A0=A0=A0 std::cerr << " fixedImageFile=A0 movingImageFile ";
=A0=A0=A0 std::cerr << " outputImagefile"=A0=A0=A0=A0 << std::endl;
=A0=A0=A0 return 1;
=A0=A0=A0 }
=A0=20
=A0 const=A0=A0=A0 unsigned int=A0=A0=A0 Dimension =3D 2;
=A0 typedef=A0 unsigned short=A0 PixelType;
=A0=20
=A0 typedef itk::Image< PixelType, Dimension >=A0 FixedImageType;
=A0 typedef itk::Image< PixelType, Dimension >=A0 MovingImageType;
=A0 typedef=A0=A0 float=A0=A0=A0=A0 InternalPixelType;
=A0 typedef itk::AffineTransform< double, Dimension > TransformType;
=A0 // Software Guide : EndCodeSnippet
=A0
=A0 typedef itk::RegularStepGradientDescentOptimizer=A0=A0=A0=A0=A0=A0 =
OptimizerType;
=A0 typedef itk::LinearInterpolateImageFunction<=20
=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=
=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0 InternalImageType,
=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=
=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0 =
double=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0 > InterpolatorType;
=A0 typedef itk::MattesMutualInformationImageToImageMetric<=20
=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=
=A0=A0=A0=A0=A0=A0=A0=A0& nbsp;=A0=A0=A0=A0=A0=A0=A0 InternalImageType,=20
=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=
=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0 InternalImageType =
>=A0=A0=A0 MetricType;
=A0 typedef OptimizerType::ScalesType=A0=A0=A0=A0=A0=A0 =
OptimizerScalesType;
=A0 typedef itk::MultiResolutionImageRegistrationMethod<=20
=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=
=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0 InternalImageType,=20
=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=
=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0 InternalImageType=A0=A0=A0 > =
RegistrationType;
=A0 typedef itk::RecursiveMultiResolutionPyramidImageFilter<
=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=
=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0 InternalImageType,
=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=
=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0 InternalImageType=A0 >=A0=A0=A0 =
FixedImagePyramidType;
=A0
=A0 typedef itk::RecursiveMultiResolutionPyramidImageFilter<
=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=
=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0 InternalImageType,
=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=
=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0 InternalImageType=A0 >=A0=A0 =
MovingImagePyramidType;
=A0 OptimizerType::Pointer=A0=A0=A0=A0=A0 optimizer=A0=A0=A0=A0 =3D =
OptimizerType::New();
=A0 InterpolatorType::Pointer=A0=A0 interpolator=A0 =3D =
InterpolatorType::New();
=A0 RegistrationType::Pointer=A0=A0 registration=A0 =3D =
RegistrationType::New();
=A0 MetricType::Pointer=A0=A0=A0=A0=A0=A0=A0=A0 =
metric=A0=A0=A0=A0=A0=A0=A0 =3D MetricType::New();
=A0 registration->SetOptimizer(=A0=A0=A0=A0 optimizer=A0=A0=A0=A0 );
=A0 registration->SetInterpolator(=A0 interpolator=A0 );
=A0 registration->SetMetric( metric=A0 );
=A0 std::cout << "\n Set :: Optimizer, interpolator, metric ";
=A0=20
=A0 TransformType::Pointer=A0=A0 transform=A0 =3D TransformType::New();
=A0 registration->SetTransform( transform );
=A0=A0 // Software Guide : EndCodeSnippet
=A0 std::cout << "\n Set :: Transform ";
=A0
=A0
=A0 FixedImagePyramidType::Pointer fixedImagePyramid =3D=20
=A0=A0=A0=A0=A0 FixedImagePyramidType::New();
=A0 MovingImagePyramidType::Pointer movingImagePyramid =3D
=A0=A0=A0=A0=A0 MovingImagePyramidType::New();
=A0 typedef itk::ImageFileReader< FixedImageType=A0 > =
FixedImageReaderType;
=A0 typedef itk::ImageFileReader< MovingImageType > =
MovingImageReaderType;
=A0 FixedImageReaderType::Pointer=A0 fixedImageReader=A0 =3D =
FixedImageReaderType::New();
=A0 MovingImageReaderType::Pointer movingImageReader =3D =
MovingImageReaderType::New();
=A0 fixedImageReader->SetFileName(=A0 argv[1] );
=A0 movingImageReader->SetFileName( argv[2] );
=A0 int nx =3D 256, ny =3D 256, nz =3D 5;
=A0 typedef=A0 itk::RawImageIO<PixelType,Dimension>=A0=A0 RawReaderType;
=A0 RawReaderType::Pointer=A0 rawReader=A0 =3D RawReaderType::New();
=A0 rawReader->SetDimensions( 0, nx );
=A0 rawReader->SetDimensions( 1, ny );
=A0// rawReader->SetDimensions( 2, nz );
=A0 fixedImageReader->SetImageIO( rawReader );
=A0 try
=A0 {
=A0=A0=A0 fixedImageReader->Update();
=A0 }
=A0 catch( itk::ExceptionObject & e )
=A0 {
=A0=A0=A0 std::cerr << "Exception caught during Raw file reading " << =
std::endl;
=A0=A0=A0 std::cerr << e << std::endl;
=A0=A0=A0 return -1;
=A0 }
=A0 movingImageReader->SetImageIO( rawReader );
=A0 try
=A0 {
=A0=A0 movingImageReader->Update();
=A0 }
=A0 catch( itk::ExceptionObject & e )
=A0 {
=A0=A0=A0 std::cerr << "Exception caught during Raw file reading " << =
std::endl;
=A0=A0=A0 std::cerr << e << std::endl;
=A0=A0=A0 return -1;
=A0 }
=A0 std::cout << "\n First file : " << argv[1] << std::endl;
=A0 std::cout << "Second file : " << argv[2] << std::endl;
=A0 typedef itk::CastImageFilter<=20
=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0 =
FixedImageType, InternalImageType > FixedCastFilterType;
=A0=20
=A0 typedef itk::CastImageFilter<=20
=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0 =
MovingImageType, InternalImageType > MovingCastFilterType;
=A0 FixedCastFilterType::Pointer fixedCaster=A0=A0 =3D =
FixedCastFilterType::New();
=A0 MovingCastFilterType::Pointer movingCaster =3D =
MovingCastFilterType::New();
=A0 fixedCaster->SetInput(=A0 fixedImageReader->GetOutput() );
=A0 std::cout << "\n\n Set Fixed Caster";
=A0 movingCaster->SetInput( movingImageReader->GetOutput() );
=A0=20
=A0 std::cout << "\n\n Set Moving Caster";
=A0 registration->SetFixedImage(=A0=A0=A0 =
fixedCaster->GetOutput()=A0=A0=A0 );
=A0 std::cout << "\n\n SetFixedImage";
=A0 registration->SetMovingImage(=A0=A0 movingCaster->GetOutput()=A0=A0 =
);
=A0 std::cout << "\n\n SetMovingImage";
=A0 fixedCaster->Update();
=A0 std::cout << "\n\n Update";
=A0 =
registration->SetFixedImageRegion(fixedCaster->GetOutput()->GetBufferedRe=
gion() );
=A0 std::cout << "\n\n SetFixedImageRegion";
=A0=A0=20
=A0transform->SetIdentity();
=A0registration->SetInitialTransformParameters( =
transform->GetParameters() );
=A0OptimizerScalesType optimizerScales( =
transform->GetNumberOfParameters() );
=A0 std::cout << "\nGet No of Params : " << =
transform->GetNumberOfParameters() << std::endl;
=A0 optimizerScales[0] =3D 1.0;// / 100.0; // scale for M11
=A0 optimizerScales[1] =3D 1.0;// / 100.0; // scale for M12
=A0 optimizerScales[2] =3D 1.0;// / 100.0; // scale for M21
=A0 optimizerScales[3] =3D 1.0;// / 100.0; // scale for M22
=A0 typedef itk::Image< InternalPixelType, Dimension > =
InternalImageType;
=A0optimizerScales[4] =3D 1.0 / 1000000.0; // scale for translation on Y
=A0 optimizerScales[5] =3D 1.0 / 1000000.0; // scale for translation on =
Z
=A0optimizer->SetScales( optimizerScales );
=A0metric->SetNumberOfHistogramBins( 50 );
=A0 metric->SetNumberOfSpatialSamples( 1000 );
=A0optimizer->SetNumberOfIterations(=A0=A0=A0 200=A0=A0 );
=A0CommandIterationUpdate::Pointer observer =3D =
CommandIterationUpdate::New();
=A0 optimizer->AddObserver( itk::IterationEvent(), observer );
=A0typedef RegistrationInterfaceCommand<RegistrationType> CommandType;
=A0 CommandType::Pointer command =3D CommandType::New();
=A0 registration->AddObserver( itk::IterationEvent(), command );
=A0
=A0 registration->SetNumberOfLevels( 3 );
=A0 std::cout << "\n Before Registration ";
=A0 try=20
=A0=A0=A0 {=20
=A0=A0=A0 registration->StartRegistration();=20
=A0=A0=A0 }=20
=A0 catch( itk::ExceptionObject & err )=20
=A0=A0=A0 {=20
=A0=A0=A0 std::cout << "ExceptionObject caught !" << std::endl;=20
=A0=A0=A0 std::cout << err << std::endl;=20
=A0=A0=A0 return -1;
=A0=A0=A0 }=20
=A0 std::cout << "Hello";
=A0 typedef RegistrationType::ParametersType ParametersType;
=A0 ParametersType finalParameters =3D =
registration->GetLastTransformParameters();
=A0=20
=A0 double TranslationAlongX =3D finalParameters[4];
=A0 double TranslationAlongY =3D finalParameters[5];
=A0=20
=A0 unsigned int numberOfIterations =3D =
optimizer->GetCurrentIteration();
=A0=20
=A0 double bestValue =3D optimizer->GetValue();
=A0=20
=A0 //
=A0 // Print out results
=A0 //
=A0 std::cout << "Result =3D " << std::endl;
=A0 std::cout << " Translation X =3D " << TranslationAlongX=A0 << =
std::endl;
=A0 std::cout << " Translation Y =3D " << TranslationAlongY=A0 << =
std::endl;
=A0 std::cout << " Iterations=A0=A0=A0 =3D " << numberOfIterations << =
std::endl;
=A0 std::cout << " Metric value=A0 =3D " << =
bestValue=A0=A0=A0=A0=A0=A0=A0=A0=A0 << std::endl;
typedef itk::ResampleImageFilter<=20
=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=
=A0=A0 MovingImageType,=20
=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=
=A0=A0 FixedIm ageType >=A0=A0=A0 ResampleFilterType;
=A0 TransformType::Pointer finalTransform =3D TransformType::New();
=A0 finalTransform->SetParameters( finalParameters );
=A0 finalTransform->PrintSelf(std::cout,NULL);
=A0 ResampleFilterType::Pointer resample =3D ResampleFilterType::New();
=A0 resample->SetTransform( finalTransform );
=A0 resample->SetInput( movingImageReader->GetOutput() );
=A0 FixedImageType::Pointer fixedImage =3D =
fixedImageReader->GetOutput();
=A0 resample->SetSize(=A0=A0=A0 =
fixedImage->GetLargestPossibleRegion().GetSize() );
=A0 resample->SetOutputOrigin(=A0 fixedImage->GetOrigin() );
=A0 resample->SetOutputSpacing( fixedImage->GetSpacing() );
=A0 resample->SetDefaultPixelValue( 100 );
=A0 typedef=A0 short=A0 OutputPixelType;
=A0 typedef itk::Image< OutputPixelType, Dimension > OutputImageType;
=A0=20
=A0 typedef itk::CastImageFilter<=20
=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0 =
FixedImageType,
=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0 =
OutputImageType > CastFilterType;
=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=20
=A0 typedef itk::ImageFileWriter< OutputImageType >=A0 WriterType;
=A0 WriterType::Pointer=A0=A0=A0=A0=A0 writer =3D=A0 WriterType::New();
=A0 CastFilterType::Pointer=A0 caster =3D=A0 CastFilterType::New();
=A0 writer->SetFileName( argv[3] );
=A0 std::cout << "\n OutPut : " << argv[3] << std::endl;
=A0=20
=A0 caster->SetInput( resample->GetOutput() );
=A0 writer->SetInput( caster->GetOutput()=A0=A0 );
=A0 writer->SetImageIO(rawReader);
=A0 std::cout << "\n Writer : SetInput, SetImageIO";
=A0 try
=A0{
=A0=A0writer->Write();
=A0=A0std::cout << "\n Write()";
=A0}
=A0 catch( itk::ExceptionObject & e )
=A0=A0=A0 {
=A0=A0=A0 std::cerr << "Exception caught during Raw file writing " << =
std::endl;
=A0=A0=A0 std::cerr << e << std::endl;
=A0=A0=A0 return -1;
=A0=A0=A0 }
=A0 std::cout << "File succesfully writen ! " << std::endl;