X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=registration%2FclitkBLUTDIRGenericFilter.cxx;h=9f7b84f456f36496f1be156400a0c43441dd2983;hb=refs%2Fheads%2FextentSimon;hp=0822bd9f3f0c434325d4995ed711d8e3e6133ba3;hpb=6194949c0beb1589904e22381b9aba1bbface172;p=clitk.git diff --git a/registration/clitkBLUTDIRGenericFilter.cxx b/registration/clitkBLUTDIRGenericFilter.cxx index 0822bd9..9f7b84f 100644 --- a/registration/clitkBLUTDIRGenericFilter.cxx +++ b/registration/clitkBLUTDIRGenericFilter.cxx @@ -29,6 +29,16 @@ It is distributed under dual licence #include "clitkBLUTDIRGenericFilter.h" #include "clitkBLUTDIRCommandIterationUpdateDVF.h" +#include "itkCenteredTransformInitializer.h" +#if ITK_VERSION_MAJOR >= 4 +# if ITK_VERSION_MINOR < 6 +# include "itkTransformToDisplacementFieldSource.h" +# else +# include "itkTransformToDisplacementFieldFilter.h" +# endif +#else +# include "itkTransformToDeformationFieldSource.h" +#endif namespace clitk { @@ -322,7 +332,6 @@ namespace clitk // The metric region with respect to the extracted transform region: // where should the metric be CALCULATED (depends on transform) typename FixedImageType::RegionType metricRegion = fixedImage->GetLargestPossibleRegion(); - typename FixedImageType::RegionType::SizeType metricRegionSize=metricRegion.GetSize(); typename FixedImageType::RegionType::IndexType metricRegionIndex=metricRegion.GetIndex(); typename FixedImageType::PointType metricRegionOrigin=fixedImage->GetOrigin(); @@ -387,6 +396,7 @@ namespace clitk // Crop the fixedImage to the bounding box to facilitate multi-resolution typedef itk::ExtractImageFilter ExtractImageFilterType; typename ExtractImageFilterType::Pointer extractImageFilter=ExtractImageFilterType::New(); + extractImageFilter->SetDirectionCollapseToSubmatrix(); extractImageFilter->SetInput(fixedImage); extractImageFilter->SetExtractionRegion(transformRegion); extractImageFilter->Update(); @@ -493,6 +503,20 @@ namespace clitk itk::Vector finalTranslation = clitk::GetTranslationPartMatrix3D(rigidTransformMatrix); rigidTransform->SetTranslation(finalTranslation); } + else if (m_ArgsInfo.centre_flag) + { + if(m_Verbose) std::cout<<"No itinial matrix given and \"centre\" flag switched on. Centering all images..."< TransformInitializerType; + typename TransformInitializerType::Pointer initializer = TransformInitializerType::New(); + initializer->SetTransform( rigidTransform ); + initializer->SetFixedImage( fixedImage ); + initializer->SetMovingImage( movingImage ); + initializer->GeometryOn(); + initializer->InitializeTransform(); + } //======================================================= @@ -584,6 +608,25 @@ namespace clitk transform->SetParameters( parameters ); if (m_ArgsInfo.initCoeff_given) initializer->SetInitialParameters(m_ArgsInfo.initCoeff_arg, parameters); + //------------------------------------------------------------------------- + // DEBUG: use an itk BSpline instead of multilabel BLUTs + //------------------------------------------------------------------------- + typedef itk::Transform< TCoordRep, 3, 3 > RegistrationTransformType; + RegistrationTransformType::Pointer regTransform(transform); + typedef itk::BSplineDeformableTransform SingleBSplineTransformType; + typename SingleBSplineTransformType::Pointer sTransform; + if(m_ArgsInfo.itkbspline_flag) { + if( transform->GetTransforms().size()>1) + itkExceptionMacro(<< "invalid --itkbspline option if there is more than 1 label") + sTransform = SingleBSplineTransformType::New(); + sTransform->SetBulkTransform( transform->GetTransforms()[0]->GetBulkTransform() ); + sTransform->SetGridSpacing( transform->GetTransforms()[0]->GetGridSpacing() ); + sTransform->SetGridOrigin( transform->GetTransforms()[0]->GetGridOrigin() ); + sTransform->SetGridRegion( transform->GetTransforms()[0]->GetGridRegion() ); + sTransform->SetParameters( transform->GetTransforms()[0]->GetParameters() ); + regTransform = sTransform; + transform = NULL; // free memory + } //======================================================= // Interpolator @@ -607,16 +650,10 @@ namespace clitk typedef itk::ImageToImageMetric< FixedImageType, MovingImageType > MetricType; typename MetricType::Pointer metric=genericMetric->GetMetricPointer(); if (movingMask) metric->SetMovingImageMask(movingMask); - -#ifdef ITK_USE_OPTIMIZED_REGISTRATION_METHODS if (threadsGiven) { metric->SetNumberOfThreads( threads ); if (m_Verbose) std::cout<< "Using " << threads << " threads." << std::endl; } -#else - if (m_Verbose) std::cout<<"Not setting the number of threads (not compiled with USE_OPTIMIZED_REGISTRATION_METHODS)..."<SetArgsInfo(m_ArgsInfo); genericOptimizer->SetMaximize(genericMetric->GetMaximize()); - genericOptimizer->SetNumberOfParameters(transform->GetNumberOfParameters()); + genericOptimizer->SetNumberOfParameters(regTransform->GetNumberOfParameters()); typedef itk::SingleValuedNonLinearOptimizer OptimizerType; OptimizerType::Pointer optimizer = genericOptimizer->GetOptimizerPointer(); @@ -638,7 +675,7 @@ namespace clitk registration->SetMetric( metric ); registration->SetOptimizer( optimizer ); registration->SetInterpolator( interpolator ); - registration->SetTransform (transform); + registration->SetTransform (regTransform ); if(threadsGiven) { registration->SetNumberOfThreads(threads); if (m_Verbose) std::cout<< "Using " << threads << " threads." << std::endl; @@ -648,7 +685,7 @@ namespace clitk registration->SetFixedImageRegion( metricRegion ); registration->SetFixedImagePyramid( fixedImagePyramid ); registration->SetMovingImagePyramid( movingImagePyramid ); - registration->SetInitialTransformParameters( transform->GetParameters() ); + registration->SetInitialTransformParameters( regTransform->GetParameters() ); registration->SetNumberOfLevels( m_ArgsInfo.levels_arg ); if (m_Verbose) std::cout<<"Setting the number of resolution levels to "< DVFCommandType; typename DVFCommandType::Pointer observerdvf = DVFCommandType::New(); @@ -693,7 +734,7 @@ namespace clitk try { - registration->StartRegistration(); + registration->Update(); } catch( itk::ExceptionObject & err ) { @@ -707,7 +748,7 @@ namespace clitk // Get the result //======================================================= OptimizerType::ParametersType finalParameters = registration->GetLastTransformParameters(); - transform->SetParameters( finalParameters ); + regTransform->SetParameters( finalParameters ); if (m_Verbose) { std::cout<<"Stop condition description: " @@ -749,15 +790,24 @@ namespace clitk typedef itk::Vector< float, SpaceDimension > DisplacementType; typedef itk::Image< DisplacementType, InputImageType::ImageDimension > DisplacementFieldType; #if ITK_VERSION_MAJOR >= 4 +# if ITK_VERSION_MINOR < 6 typedef itk::TransformToDisplacementFieldSource ConvertorType; -#else - typedef itk::TransformToDeformationFieldSource ConvertorType; +# else + typedef itk::TransformToDisplacementFieldFilter ConvertorType; +# endif #endif typename ConvertorType::Pointer filter= ConvertorType::New(); filter->SetNumberOfThreads(1); - transform->SetBulkTransform(NULL); - filter->SetTransform(transform); + if(m_ArgsInfo.itkbspline_flag) + sTransform->SetBulkTransform(NULL); + else + transform->SetBulkTransform(NULL); + filter->SetTransform(regTransform); +#if ITK_VERSION_MAJOR > 4 || (ITK_VERSION_MAJOR == 4 && ITK_VERSION_MINOR >= 6) + filter->SetReferenceImage(fixedImage); +#else filter->SetOutputParametersFromImage(fixedImage); +#endif filter->Update(); typename DisplacementFieldType::Pointer field = filter->GetOutput(); @@ -786,8 +836,13 @@ namespace clitk //======================================================= typedef itk::ResampleImageFilter< MovingImageType, FixedImageType > ResampleFilterType; typename ResampleFilterType::Pointer resampler = ResampleFilterType::New(); - if (rigidTransform) transform->SetBulkTransform(rigidTransform); - resampler->SetTransform( transform ); + if (rigidTransform) { + if(m_ArgsInfo.itkbspline_flag) + sTransform->SetBulkTransform(rigidTransform); + else + transform->SetBulkTransform(rigidTransform); + } + resampler->SetTransform( regTransform ); resampler->SetInput( movingImage); resampler->SetOutputParametersFromImage(fixedImage); resampler->Update();