[Insight-users] [ITK Community] 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://www.itk.org/pipermail/insight-users/attachments/20131214/6fbccd60/attachment.htm>
More information about the Insight-users
mailing list