[ITK Community] [Insight-users] registration using Mutual Information

Luis Ibanez luis.ibanez at kitware.com
Sat Dec 14 10:46:18 EST 2013


Massinissa,

The fact that you are compiling for Debug, explains why the code is taking
a long time to run.

Please give it a try at compiling for Release, that can easily make a
difference of 10X or 20X
in the execution time. Most of the C++ template optimizations that ITK
relies on, only happen
when you compile for Release.

Once you run it in Release mode, please share the output of the Command
observer class.
This should help illuminate the progression of the transform and the metric
values.

    Thanks

        Luis




On Fri, Dec 13, 2013 at 2:59 PM, Massinissa Bandou <
Massinissa.Bandou at usherbrooke.ca> wrote:

> Thank you for replying Mr Luis!!!!
>
> I'm using win7 x64 and I'm compiling on debug mode.
>
> The *source* image is from a series of dicom
> format: float
> image dimension: 120 x 120 x 137
> voxel dimension: 0.5 x 0.5 x 0.596756
>
> The *target* image is a raw file.
> format: short
> image dimension: 256 x 256 x 256
> voxel dimension: 0.308 x 0.308 x 0.308
>
> I decided to register 2 same images and use the same code as
> *Registration/ImageRegistration8.cxx*
> but I got the following error:
>
> <http://itk-users.7.n7.nabble.com/file/n32985/error.png>
>
> and my code is:
>
> int main(int argc, char *argv[])
> {
>
> vtkSmartPointer<vtkImageReader> imagetarget =
> vtkSmartPointer<vtkImageReader>::New();
>
>
> imagetarget->SetFileName("C:/Users/Massi/Desktop/Engineering/IRM_CT_FIDUCIALS/CT_(130802_m_FDGF18_004ct)/image/CT_volume.raw");
>         imagetarget->SetFileDimensionality(3);
>         imagetarget->SetDataExtent(0,511,0,511,0,511);
>         imagetarget->SetDataSpacing(0.154,0.154,0.154);
>         imagetarget->SetDataByteOrderToLittleEndian();
>         imagetarget->SetDataScalarTypeToShort();
>         imagetarget->Update();
>
>         vtkSmartPointer<vtkImageReslice> resliceMyTargetImage =
> vtkSmartPointer<vtkImageReslice>::New();
>         resliceMyTargetImage
> ->SetInputConnection(imagetarget->GetOutputPort());
>         resliceMyTargetImage->SetOutputExtent(0,255, 0,255, 0,255);
>         double p[3];
>         for(unsigned int i=0;i<3;i++){
>                 p[i] = (512 * 0.154) / 256;
>         }
>         resliceMyTargetImage->SetOutputSpacing(p[0],p[1],p[2]);
>         resliceMyTargetImage->SetOutputDimensionality(3);
>         resliceMyTargetImage->SetInterpolationModeToCubic();
>         resliceMyTargetImage->AutoCropOutputOn();
>         resliceMyTargetImage->SetResliceAxesOrigin(0,0,0);
>         resliceMyTargetImage->Update();
>
>         vtkSmartPointer<vtkDICOMImageReader> imagesource =
> vtkSmartPointer<vtkDICOMImageReader>::New();
>
>
> imagesource->SetDirectoryName("C:/Users/Massi/Downloads/M016_Labpet8_20130425/T_6");
>         imagesource->Update();
>
>
>         typedef itk::Image<float, 3> FixedImageType;
>         typedef itk::Image<short, 3> MovingImageType;
>
>         /*transform vtk image to itk image*/
>         typedef itk::VTKImageToImageFilter<FixedImageType>
> FixedImageReaderType;
>         typedef itk::VTKImageToImageFilter<MovingImageType>
> MovingImageReaderType;
>
>     FixedImageReaderType::Pointer fixedImageReader =
> FixedImageReaderType::New();
>         MovingImageReaderType::Pointer movingImageReader =
> MovingImageReaderType::New();
>         fixedImageReader->SetInput(imagesource->GetOutput());  // source
>         movingImageReader->SetInput(resliceMyTargetImage->GetOutput());
> //target
>
>         typedef itk::VersorRigid3DTransform< double > TransformType;
>         typedef itk::VersorRigid3DTransformOptimizer
> 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 );
>
>         fixedImageReader->Update();
>
>
> registration->SetFixedImageRegion(fixedImageReader->GetOutput()->GetBufferedRegion());
>
>         typedef
>
> itk::CenteredTransformInitializer<TransformType,FixedImageType,MovingImageType>
> TransformInitializerType;
>         TransformInitializerType::Pointer initializer =
> TransformInitializerType::New();
>
>         initializer->SetTransform(transform);
>         initializer->SetFixedImage(fixedImageReader->GetOutput());
>         initializer->SetMovingImage(movingImageReader->GetOutput());
>         initializer->MomentsOn();
>         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 / 1000.0;
>         optimizerScales[0] = 1.0;
>         optimizerScales[1] = 1.0;
>         optimizerScales[2] = 1.0;
>         optimizerScales[3] = translationScale;
>         optimizerScales[4] = translationScale;
>         optimizerScales[5] = translationScale;
>         optimizer->SetScales( optimizerScales );
>         optimizer->SetMaximumStepLength( 0.2000  );
>         optimizer->SetMinimumStepLength( 0.0001 );
>         optimizer->SetNumberOfIterations( 200 );
>
>         CommandIterationUpdate::Pointer observer =
> CommandIterationUpdate::New();
>         optimizer->AddObserver( itk::IterationEvent(), observer );
>
>         try{
>                 registration->Update();
>                 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;
>     }
>
>         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;
>
>         transform->SetParameters( finalParameters );
>         TransformType::MatrixType matrix = transform->GetMatrix();
>         TransformType::OffsetType offset = transform->GetOffset();
>         std::cout << "Matrix = " << std::endl << matrix << std::endl;
>         std::cout << "Offset = " << std::endl << offset << std::endl;
>
>         typedef itk::ResampleImageFilter< MovingImageType,FixedImageType >
> ResampleFilterType;
>         TransformType::Pointer finalTransform = TransformType::New();
>         finalTransform->SetCenter( transform->GetCenter() );
>         finalTransform->SetParameters( finalParameters );
>         finalTransform->SetFixedParameters(
> transform->GetFixedParameters() );
>         ResampleFilterType::Pointer resampler = ResampleFilterType::New();
>         resampler->SetTransform( finalTransform );
>         resampler->SetInput( movingImageReader->GetOutput() );
>         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 );
>
>         typedef unsigned char
> OutputPixelType;
>         typedef itk::Image< OutputPixelType, 3 >
>  OutputImageType;
>         typedef itk::CastImageFilter< FixedImageType, OutputImageType >
> CastFilterType;
>         typedef itk::ImageFileWriter< OutputImageType >
> WriterType;
>         //WriterType::Pointer      writer =  WriterType::New();
>         CastFilterType::Pointer  caster =  CastFilterType::New();
>         //writer->SetFileName( "result.raw" );
>         caster->SetInput( resampler->GetOutput() );
>         //writer->SetInput( caster->GetOutput()   );
>         //writer->Update();
>
>         typedef
> itk::SubtractImageFilter<FixedImageType,FixedImageType,FixedImageType >
> DifferenceFilterType;
>         DifferenceFilterType::Pointer difference =
> DifferenceFilterType::New();
>         typedef
> itk::RescaleIntensityImageFilter<FixedImageType,OutputImageType >
> RescalerType;
>         RescalerType::Pointer intensityRescaler = RescalerType::New();
>         intensityRescaler->SetInput( difference->GetOutput() );
>         intensityRescaler->SetOutputMinimum(   0 );
>         intensityRescaler->SetOutputMaximum( 255 );
>         difference->SetInput1( fixedImageReader->GetOutput() );
>         difference->SetInput2( resampler->GetOutput() );
>         resampler->SetDefaultPixelValue( 1 );
>
>
>         /*connect to VTK*/
>         typedef itk::ImageToVTKImageFilter<OutputImageType> ConnectorType;
>         ConnectorType::Pointer connector = ConnectorType::New();
>         connector->SetInput(caster->GetOutput());
>         connector->Update();
>
>         vtkSmartPointer<vtkGPUVolumeRayCastMapper> volumeMapper =
> vtkSmartPointer<vtkGPUVolumeRayCastMapper>::New();
>         volumeMapper->SetInput(connector->GetOutput());
>         volumeMapper->AutoAdjustSampleDistancesOff();
>         volumeMapper->SetBlendModeToComposite();
>
>         /********************set color & opacity*******************/
>     vtkSmartPointer<vtkVolumeProperty> volumeProperty =
> vtkSmartPointer<vtkVolumeProperty>::New();
>     vtkSmartPointer<vtkPiecewiseFunction> compositeOpacity =
> vtkSmartPointer<vtkPiecewiseFunction>::New();
>         vtkSmartPointer<vtkColorTransferFunction> color =
> vtkSmartPointer<vtkColorTransferFunction>::New();
>
>
>
> compositeOpacity->AddPoint(connector->GetOutput()->GetScalarRange()[0],0.0);
>
>
> compositeOpacity->AddPoint(connector->GetOutput()->GetScalarRange()[1]*3/255,0.0);
>
>
> compositeOpacity->AddPoint(connector->GetOutput()->GetScalarRange()[1]*15/255,0.1);
>
>
> compositeOpacity->AddPoint(connector->GetOutput()->GetScalarRange()[1]*70/255,0.2);
>
>
> compositeOpacity->AddPoint(connector->GetOutput()->GetScalarRange()[1]*75/255,0.3);
>
>
>
> color->AddRGBPoint(connector->GetOutput()->GetScalarRange()[1]*0/255,0.0,0.0,0.0);
>
>
> color->AddRGBPoint(connector->GetOutput()->GetScalarRange()[1]*3/255,1.0,0.3,0.3);//
>
>
> color->AddRGBPoint(connector->GetOutput()->GetScalarRange()[1]*20/255,0.0,0.0,1.0);
>
>
> color->AddRGBPoint(connector->GetOutput()->GetScalarRange()[1]*50/255,0.9,0.37,0.27);
>
>
> color->AddRGBPoint(connector->GetOutput()->GetScalarRange()[1]*80.0/255,0.15,0.15,0.5);
>
>
> color->AddRGBPoint(connector->GetOutput()->GetScalarRange()[1]*120.0/255,1.0,1.0,1.0);
>
>                         volumeProperty->SetAmbient(0.5);
>                         volumeProperty->SetDiffuse(0.8);
>                         volumeProperty->SetSpecular(8.0);
>                         volumeProperty->SetInterpolationTypeToLinear();
>     volumeProperty->SetColor(color);
>         volumeProperty->SetScalarOpacity(compositeOpacity);
>     volumeProperty->ShadeOn();
>
>     vtkSmartPointer<vtkVolume> volume = vtkSmartPointer<vtkVolume>::New();
>     volume->SetMapper(volumeMapper);
>     volume->SetProperty(volumeProperty);
>
>         vtkSmartPointer<vtkRenderWindow> renderWindow =
> vtkSmartPointer<vtkRenderWindow>::New();
>         vtkSmartPointer<vtkRenderer> renderer =
> vtkSmartPointer<vtkRenderer>::New();
>         renderWindow->AddRenderer(renderer);
>         renderer->AddVolume(volume);
>         renderer->SetBackground(0.9,0.9,0.9);
>         renderer->GradientBackgroundOn();
>         vtkSmartPointer<vtkRenderWindowInteractor> renderWindowInteractor =
>     vtkSmartPointer<vtkRenderWindowInteractor>::New();
>   renderWindowInteractor->SetRenderWindow(renderWindow);
>
>         /*Set camera position*/
>         renderer->GetActiveCamera()->Azimuth(30);
>         renderer->GetActiveCamera()->Elevation(30);
>         renderer->ResetCamera();
>           renderWindow->Render();
>   renderWindowInteractor->Start();
>
>         return EXIT_SUCCESS;
> }
>
>
>
>
> --
> View this message in context:
> http://itk-users.7.n7.nabble.com/registration-using-Mutual-Information-tp32972p32985.html
> Sent from the ITK - Users mailing list archive at Nabble.com.
> _____________________________________
> Powered by www.kitware.com
>
> Visit other Kitware open-source projects at
> http://www.kitware.com/opensource/opensource.html
>
> Kitware offers ITK Training Courses, for more information visit:
> http://www.kitware.com/products/protraining.php
>
> 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
> _______________________________________________
> Community mailing list
> Community at itk.org
> http://public.kitware.com/cgi-bin/mailman/listinfo/community
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/community/attachments/20131214/6fbccd60/attachment.html>
-------------- next part --------------
_____________________________________
Powered by www.kitware.com

Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html

Kitware offers ITK Training Courses, for more information visit:
http://www.kitware.com/products/protraining.php

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 Community mailing list