X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=registration%2FclitkBLUTDIRGenericFilter.cxx;h=ca10c2d04d4fc4adc52e035d160b66843808a663;hb=2b4ec898fdd6dc8ea5780c0e30c7124272a4202e;hp=c8c42f0f422f039bd30814bac37bcad99a44cd86;hpb=df51e4850f784f8f2fbaead7bb82bd6897fb9fab;p=clitk.git diff --git a/registration/clitkBLUTDIRGenericFilter.cxx b/registration/clitkBLUTDIRGenericFilter.cxx old mode 100755 new mode 100644 index c8c42f0..ca10c2d --- a/registration/clitkBLUTDIRGenericFilter.cxx +++ b/registration/clitkBLUTDIRGenericFilter.cxx @@ -3,7 +3,7 @@ Program: vv http://www.creatis.insa-lyon.fr/rio/vv Authors belong to: - University of LYON http://www.universite-lyon.fr/ -- Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr +- Léon Bérard cancer center http://www.centreleonberard.fr - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr This software is distributed WITHOUT ANY WARRANTY; without even @@ -14,7 +14,7 @@ It is distributed under dual licence - BSD See included LICENSE.txt file - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html -======================================================================-====*/ +===========================================================================**/ #ifndef clitkBLUTDIRGenericFilter_cxx #define clitkBLUTDIRGenericFilter_cxx @@ -28,6 +28,7 @@ It is distributed under dual licence ===================================================*/ #include "clitkBLUTDIRGenericFilter.h" +#include "clitkBLUTDIRCommandIterationUpdateDVF.h" namespace clitk { @@ -78,7 +79,7 @@ namespace clitk { InitializeImageType<2>(); InitializeImageType<3>(); - m_Verbose=true; + m_Verbose=false; } //=========================================================================// @@ -93,6 +94,8 @@ namespace clitk } if (m_ArgsInfo.output_given) SetOutputFilename(m_ArgsInfo.output_arg); + + if (m_ArgsInfo.verbose_given) m_Verbose=true; } //=========================================================================// @@ -273,6 +276,8 @@ namespace clitk template void BLUTDIRGenericFilter::UpdateWithInputImageType() { + if (m_Verbose) std::cout << "BLUTDIRGenericFilter::UpdateWithInputImageType()" << std::endl; + //============================================================================= //Input //============================================================================= @@ -317,7 +322,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(); @@ -382,6 +386,9 @@ namespace clitk // Crop the fixedImage to the bounding box to facilitate multi-resolution typedef itk::ExtractImageFilter ExtractImageFilterType; typename ExtractImageFilterType::Pointer extractImageFilter=ExtractImageFilterType::New(); +#if ITK_VERSION_MAJOR == 4 + extractImageFilter->SetDirectionCollapseToSubmatrix(); +#endif extractImageFilter->SetInput(fixedImage); extractImageFilter->SetExtractionRegion(transformRegion); extractImageFilter->Update(); @@ -579,6 +586,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 @@ -603,7 +629,7 @@ namespace clitk typename MetricType::Pointer metric=genericMetric->GetMetricPointer(); if (movingMask) metric->SetMovingImageMask(movingMask); -#ifdef ITK_USE_OPTIMIZED_REGISTRATION_METHODS +#if defined(ITK_USE_OPTIMIZED_REGISTRATION_METHODS) || ITK_VERSION_MAJOR >= 4 if (threadsGiven) { metric->SetNumberOfThreads( threads ); if (m_Verbose) std::cout<< "Using " << threads << " threads." << std::endl; @@ -620,7 +646,7 @@ namespace clitk GenericOptimizerType::Pointer genericOptimizer = GenericOptimizerType::New(); genericOptimizer->SetArgsInfo(m_ArgsInfo); genericOptimizer->SetMaximize(genericMetric->GetMaximize()); - genericOptimizer->SetNumberOfParameters(transform->GetNumberOfParameters()); + genericOptimizer->SetNumberOfParameters(regTransform->GetNumberOfParameters()); typedef itk::SingleValuedNonLinearOptimizer OptimizerType; OptimizerType::Pointer optimizer = genericOptimizer->GetOptimizerPointer(); @@ -633,7 +659,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; @@ -643,7 +669,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 "<SetMaximize(genericMetric->GetMaximize()); command->SetMetricRegion(metricRegion); registration->AddObserver( itk::IterationEvent(), command ); + + if (m_ArgsInfo.coeff_given) + { + if(m_ArgsInfo.itkbspline_flag) { + itkExceptionMacro("--coeff and --itkbpline are incompatible"); + } + + std::cout << std::endl << "Output coefficient images every " << m_ArgsInfo.coeffEveryN_arg << " iterations." << std::endl; + typedef CommandIterationUpdateDVF DVFCommandType; + typename DVFCommandType::Pointer observerdvf = DVFCommandType::New(); + observerdvf->SetFixedImage(fixedImage); + observerdvf->SetTransform(transform); + observerdvf->SetArgsInfo(m_ArgsInfo); + optimizer->AddObserver( itk::IterationEvent(), observerdvf ); + } } @@ -691,7 +732,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: " @@ -731,21 +772,28 @@ namespace clitk // Compute the DVF (only deformable transform) //======================================================= typedef itk::Vector< float, SpaceDimension > DisplacementType; - typedef itk::Image< DisplacementType, InputImageType::ImageDimension > DeformationFieldType; - typedef itk::TransformToDeformationFieldSource ConvertorType; + typedef itk::Image< DisplacementType, InputImageType::ImageDimension > DisplacementFieldType; +#if ITK_VERSION_MAJOR >= 4 + typedef itk::TransformToDisplacementFieldSource ConvertorType; +#else + typedef itk::TransformToDeformationFieldSource ConvertorType; +#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); filter->SetOutputParametersFromImage(fixedImage); filter->Update(); - typename DeformationFieldType::Pointer field = filter->GetOutput(); + typename DisplacementFieldType::Pointer field = filter->GetOutput(); //======================================================= // Write the DVF //======================================================= - typedef itk::ImageFileWriter< DeformationFieldType > FieldWriterType; + typedef itk::ImageFileWriter< DisplacementFieldType > FieldWriterType; typename FieldWriterType::Pointer fieldWriter = FieldWriterType::New(); fieldWriter->SetFileName( m_ArgsInfo.vf_arg ); fieldWriter->SetInput( field ); @@ -766,8 +814,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();