[Insight-users] Problem with Rigid3DTransform
Luis Ibanez
luis.ibanez at kitware.com
Thu Mar 12 09:58:28 EDT 2009
Hi Serena,
Thanks for your detailed report.
A) By "n.160" do you mean the "Iteration Number 160" ?
or are you referring to slice number 160 in the resampled image ?
B) What do you mean by
"the other have half full statistics" ??
C) It is perfectly fine to combine the VersorRigid3DTransform
with the MattesMutualInformationMetrics. There is no reason
why they couldn't work together. We used this combination
on a regular basis.
D) It is strange for the resampled image to have a corrupted
slice...
Can you describe what you mean by "corrupted" ?
It may indicate a memory management problem ....
It will be ideal if you could post this output image
in a public web site, so that we can take a look at it.
Please let us know about items (A) and (B).
Thanks
Luis
------------------------
Serena Fabbri wrote:
>
> Hi All,
>
> I am trying to register T1-MRI and T2-MRI images with
>
> itkVersorRigid3DTransform
> itkVersorRigid3DTransformOptimizer
> itkMattesMutualInformationImageToImageMetric
>
> fixed: T2
> moving: T1
>
> but i find that T1-reg has the last 4 slides corrupted.
> the n. 160 is empty and the other have half full statistics.
> I don't think that this task is difficult for MI because the 2 images
> look very similar and I did a test: i registered T1 with T1 and i found
> T1-reg corrupted.
> I think I made a mistake but i don't recognize where it is.
> maybe it is not possible to use
> MattesMutualInformationImageToImageMetric with VersorRigid3DTransform,
> is it?
>
>
>
> T1
> dim 176 256 160
> spacing 1 1 1
> origin 0 0 0
>
> T2
> dim 392 512 160
> spacing 0.5 0.5 1 origin 0 0 0
>
> they are mha format
>
> this is my code.....any idea?
>
> Thanks.
>
> Serena.
>
>
>
>
> const unsigned int Dimension = 3;
> typedef float PixelType;
>
> typedef itk::Image< PixelType, Dimension > FixedImageType;
> typedef itk::Image< PixelType, Dimension > MovingImageType;
>
>
>
> typedef itk::VersorRigid3DTransform< double > TransformType;
>
>
>
>
> typedef itk::VersorRigid3DTransformOptimizer OptimizerType;
>
> typedef itk::MattesMutualInformationImageToImageMetric<
> 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->SetOptimizer( optimizer );
> registration->SetInterpolator( interpolator );
>
>
>
>
> registration->SetMetric( metric );
> metric->SetNumberOfHistogramBins( 32 );
> metric->SetNumberOfSpatialSamples( 20000 );
>
> if( argc > 9 )
> {
>
> metric->SetUseExplicitPDFDerivatives( atoi( argv[10] ) );
> }
>
> if( argc > 10 )
> {
>
> metric->SetUseCachingOfBSplineWeights( atoi( argv[11] ) );
> }
>
>
>
> TransformType::Pointer transform = TransformType::New();
> registration->SetTransform( transform );
>
>
>
> typedef itk::ImageFileReader< FixedImageType > FixedImageReaderType;
> typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;
>
> FixedImageReaderType::Pointer fixedImageReader =
> FixedImageReaderType::New();
> MovingImageReaderType::Pointer movingImageReader =
> MovingImageReaderType::New();
>
> fixedImageReader->SetFileName( argv[1] ); fixedImageReader->Update();
> movingImageReader->SetFileName( argv[2] ); movingImageReader->Update();
>
>
>
>
>
> typedef itk::RescaleIntensityImageFilter< FixedImageType, FixedImageType
> > RescalerType;
>
> RescalerType::Pointer intensityRescaler1 = RescalerType::New();
>
> intensityRescaler1->SetInput( fixedImageReader->GetOutput() );
> intensityRescaler1->SetOutputMinimum( 0 );
> intensityRescaler1->SetOutputMaximum( 255 );
> intensityRescaler1->Update();
>
> RescalerType::Pointer intensityRescaler2 = RescalerType::New();
>
> intensityRescaler2->SetInput( movingImageReader->GetOutput() );
> intensityRescaler2->SetOutputMinimum( 0 );
> intensityRescaler2->SetOutputMaximum( 255 );
> intensityRescaler2->Update();
>
> registration->SetFixedImage( intensityRescaler1->GetOutput() );
> registration->SetMovingImage( intensityRescaler2->GetOutput() );
>
>
>
>
>
>
> registration->SetFixedImageRegion(
> intensityRescaler1->GetOutput()->GetBufferedRegion() );
>
>
>
> typedef itk::CenteredTransformInitializer< TransformType,
> FixedImageType,
> MovingImageType
> >
> TransformInitializerType;
>
> TransformInitializerType::Pointer initializer =
> TransformInitializerType::New();
>
> initializer->SetTransform( transform );
> initializer->SetFixedImage( intensityRescaler1->GetOutput() );
> initializer->SetMovingImage( intensityRescaler2->GetOutput() );
> initializer->GeometryOn();
>
> initializer->InitializeTransform();
>
>
> typedef TransformType::VersorType VersorType;
> typedef VersorType::VectorType VectorType;
>
> VersorType rotation;
> VectorType axis;
>
> axis[0] = 0.0;
> axis[1] = 0.0;
> axis[2] = 1.0;
>
> const double angle = 0;
>
> rotation.Set( axis, angle );
>
> transform->SetRotation( rotation );
>
> registration->SetInitialTransformParameters(
> transform->GetParameters() );
>
> typedef OptimizerType::ScalesType OptimizerScalesType;
> OptimizerScalesType optimizerScales(
> transform->GetNumberOfParameters() );
> const double translationScale = 1.0 / 2000.0;
>
> optimizerScales[0] = 1.0;//original
> optimizerScales[1] = 1.0;
> optimizerScales[2] = 1.0;
> optimizerScales[3] = translationScale;
> optimizerScales[4] = translationScale;
> optimizerScales[5] = translationScale;
>
> optimizer->SetScales( optimizerScales );
>
> optimizer->SetMaximumStepLength( 0.1000 );//original 0.2000
> optimizer->SetMinimumStepLength( 0.0001 );
>
> optimizer->SetNumberOfIterations( 500 );
>
>
>
> CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
> optimizer->AddObserver( itk::IterationEvent(), observer );
>
>
> try
> {
> registration->StartRegistration();
> }
> catch( itk::ExceptionObject & err )
> {
> std::cerr << "ExceptionObject caught !" << std::endl;
> std::cerr << err << std::endl;
> return EXIT_FAILURE;
> }
>
> OptimizerType::ParametersType finalParameters =
> registration->GetLastTransformParameters();
>
>
> const double versorX = finalParameters[0];
> const double versorY = finalParameters[1];
> const double versorZ = finalParameters[2];
> const double finalTranslationX = finalParameters[3];
> const double finalTranslationY = finalParameters[4];
> const double finalTranslationZ = finalParameters[5];
>
> const unsigned int numberOfIterations = optimizer->GetCurrentIteration();
>
> const double bestValue = optimizer->GetValue();
>
>
>
> std::cout << std::endl << std::endl;
> std::cout << "Result = " << std::endl;
> std::cout << " versor X = " << versorX << std::endl;
> std::cout << " versor Y = " << versorY << std::endl;
> std::cout << " versor Z = " << versorZ << std::endl;
> std::cout << " Translation X = " << finalTranslationX << std::endl;
> std::cout << " Translation Y = " << finalTranslationY << std::endl;
> std::cout << " Translation Z = " << finalTranslationZ << std::endl;
> std::cout << " Iterations = " << numberOfIterations << std::endl;
> std::cout << " Metric value = " << bestValue << std::endl;
>
>
>
>
>
> typedef itk::ResampleImageFilter<
> MovingImageType,
> FixedImageType > ResampleFilterType;
>
> TransformType::Pointer finalTransform = TransformType::New();
>
> finalTransform->SetCenter( transform->GetCenter() );
>
> finalTransform->SetParameters( finalParameters );
>
> ResampleFilterType::Pointer resampler = ResampleFilterType::New();
>
> resampler->SetTransform( finalTransform );
>
> resampler->SetInput( intensityRescaler2->GetOutput() );
>
>
> FixedImageType::Pointer fixedImage = intensityRescaler1->GetOutput();
>
> resampler->SetSize(
> fixedImage->GetLargestPossibleRegion().GetSize() );
> resampler->SetOutputOrigin( fixedImage->GetOrigin() );
> resampler->SetOutputSpacing( fixedImage->GetSpacing() );
> resampler->SetOutputDirection( fixedImage->GetDirection() );
> resampler->SetDefaultPixelValue( 100 );
>
> 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();
>
>
> writer->SetFileName( argv[3] );
>
>
>
> caster->SetInput( resampler->GetOutput() );
> writer->SetInput( caster->GetOutput() );
> writer->Update();
>
>
>
>
>
> _____________________________________
> Powered by www.kitware.com
>
> Visit other Kitware open-source projects at
> http://www.kitware.com/opensource/opensource.html
>
> Please keep messages on-topic and check the ITK FAQ at:
> http://www.itk.org/Wiki/ITK_FAQ
>
> Follow this link to subscribe/unsubscribe:
> http://www.itk.org/mailman/listinfo/insight-users
>
More information about the Insight-users
mailing list