[Insight-users] Problem with Rigid3DTransform
Luis Ibanez
luis.ibanez at kitware.com
Thu Mar 12 17:07:47 EDT 2009
Hi Serena,
Thanks a lot for uploading your resampled image.
The good news is that:
The image is perfectly fine. :-)
What you are observing is the normal process of
filling up the borders of the image with the
default value "100" that you set in the resample
filter with:
resampler->SetDefaultPixelValue( 100 );
You will find this explained in the ITK Software Guide
http://www.itk.org/ItkSoftwareGuide.pdf
In the Resample Image sections, as well as in the
Image Registration Chapter.
What seems to have happened is that your fixed
image is a bit larger than the moving image, and
therefore when you map the moving image into the
coordinate system of the fixed image, there is room
left along the sides. Slices 157 and 158 looks cut
long a diagonal, simply because your final transform
rotated the image a bit and therefore the planes
of the moving image are not parallel to the ones
of the fixed image.
Disregard, my comments about a potential memory
problem. The slices in your image are not really
corrupted, they are simply filled up with the
"default" intensity value used by the resample
image filter.
Regards,
Luis
-----------------------
Serena Fabbri wrote:
>
> Hi Luis, Hi All,
>
> Thanks for your reply.
>
> A) n.160 is the slice number 160
> B)slice n. 160 is empty
> slice n.159, n.158 have a corner empty.
>
> This is the output :
>
>
>
> Result =
> versor X = -0.00406936
> versor Y = -0.0015104
> versor Z = 0.00378959
> Translation X = -9.60539
> Translation Y = -14.3073
> Translation Z = 1.52914
> Iterations = 203
> Metric value = -0.204644
> Matrix = 0.999967 -0.00756676 -0.00305159
> 0.00759135 0.999938 0.00812714
> 0.0029899 -0.00815003 0.999962
>
> Offset = [-8.39288, -15.6875, 2.28104]
>
>
> What do you mean with "memory problem"?
> what do I have to check?
> (I am using ITK 3.10.0)
>
>
>
> thanks again.
>
> Serena.
>
>
> On Thu, 12 Mar 2009, Luis Ibanez wrote:
>
>>
>> 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