- /*=========================================================================
+ /*=========================================================================
Program: vv http://www.creatis.insa-lyon.fr/rio/vv
- Authors belong to:
+ Authors belong to:
- University of LYON http://www.universite-lyon.fr/
- Léon Bérard cancer center http://www.centreleonberard.fr
- CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr
namespace clitk
{
-//==============================================================================
-//Creating an observer class that allows us to change parameters at subsequent levels
-//==============================================================================
-template <typename TRegistration>
-class RegistrationInterfaceCommand : public itk::Command
-{
-public:
- typedef RegistrationInterfaceCommand Self;
-
-
- typedef itk::Command Superclass;
- typedef itk::SmartPointer<Self> Pointer;
- itkNewMacro( Self );
-protected:
- RegistrationInterfaceCommand() {};
-public:
- typedef TRegistration RegistrationType;
- typedef RegistrationType * RegistrationPointer;
-
- // Two arguments are passed to the Execute() method: the first
- // is the pointer to the object which invoked the event and the
- // second is the event that was invoked.
- void Execute(itk::Object * object, const itk::EventObject & event) {
- if( !(itk::IterationEvent().CheckEvent( &event )) ) {
- return;
+ //==============================================================================
+ //Creating an observer class that allows us to change parameters at subsequent levels
+ //==============================================================================
+ template <typename TRegistration>
+ class RegistrationInterfaceCommand : public itk::Command
+ {
+ public:
+ typedef RegistrationInterfaceCommand Self;
+
+
+ typedef itk::Command Superclass;
+ typedef itk::SmartPointer<Self> Pointer;
+ itkNewMacro( Self );
+ protected:
+ RegistrationInterfaceCommand() {};
+ public:
+ typedef TRegistration RegistrationType;
+ typedef RegistrationType * RegistrationPointer;
+
+ // Two arguments are passed to the Execute() method: the first
+ // is the pointer to the object which invoked the event and the
+ // second is the event that was invoked.
+ void Execute(itk::Object * object, const itk::EventObject & event)
+ {
+ if( !(itk::IterationEvent().CheckEvent( &event )) )
+ {
+ return;
+ }
+ RegistrationPointer registration = dynamic_cast<RegistrationPointer>( object );
+ unsigned int numberOfLevels=registration->GetNumberOfLevels();
+ unsigned int currentLevel=registration->GetCurrentLevel()+1;
+ std::cout<<std::endl;
+ std::cout<<"========================================"<<std::endl;
+ std::cout<<"Starting resolution level "<<currentLevel<<" of "<<numberOfLevels<<"..."<<std::endl;
+ std::cout<<"========================================"<<std::endl;
+ std::cout<<std::endl;
}
- RegistrationPointer registration = dynamic_cast<RegistrationPointer>( object );
- unsigned int numberOfLevels=registration->GetNumberOfLevels();
- unsigned int currentLevel=registration->GetCurrentLevel()+1;
- std::cout<<std::endl;
- std::cout<<"========================================"<<std::endl;
- std::cout<<"Starting resolution level "<<currentLevel<<" of "<<numberOfLevels<<"..."<<std::endl;
- std::cout<<"========================================"<<std::endl;
- std::cout<<std::endl;
- }
-
- void Execute(const itk::Object * , const itk::EventObject & ) {
- return;
- }
-
-};
-
-
-//==============================================================================
-// Creating an observer class that allows output at each iteration
-//==============================================================================
-class CommandIterationUpdate : public itk::Command
-{
-public:
- typedef CommandIterationUpdate Self;
- typedef itk::Command Superclass;
- typedef itk::SmartPointer<Self> Pointer;
- itkNewMacro( Self );
-protected:
- CommandIterationUpdate() {};
-public:
- typedef clitk::GenericOptimizer<args_info_clitkBSplineDeformableRegistration> OptimizerType;
- typedef const OptimizerType * OptimizerPointer;
-
- // We set the generic optimizer
- void SetOptimizer(OptimizerPointer o) {
- m_Optimizer=o;
- }
-
- // Execute
- void Execute(itk::Object *caller, const itk::EventObject & event) {
- Execute( (const itk::Object *)caller, event);
- }
-
- void Execute(const itk::Object * object, const itk::EventObject & event) {
- if( !(itk::IterationEvent().CheckEvent( &event )) ) {
- return;
+
+ void Execute(const itk::Object * , const itk::EventObject & )
+ { return; }
+
+ };
+
+
+ //==============================================================================
+ // Creating an observer class that allows output at each iteration
+ //==============================================================================
+ class CommandIterationUpdate : public itk::Command
+ {
+ public:
+ typedef CommandIterationUpdate Self;
+ typedef itk::Command Superclass;
+ typedef itk::SmartPointer<Self> Pointer;
+ itkNewMacro( Self );
+ protected:
+ CommandIterationUpdate() {};
+ public:
+ typedef clitk::GenericOptimizer<args_info_clitkBSplineDeformableRegistration> OptimizerType;
+ typedef const OptimizerType * OptimizerPointer;
+
+ // We set the generic optimizer
+ void SetOptimizer(OptimizerPointer o){m_Optimizer=o;}
+
+ // Execute
+ void Execute(itk::Object *caller, const itk::EventObject & event)
+ {
+ Execute( (const itk::Object *)caller, event);
}
-
- m_Optimizer->OutputIterationInfo();
- }
-
- OptimizerPointer m_Optimizer;
-};
-
-
-//==============================================================================
-// Update with the number of dimensions
-//==============================================================================
-template<unsigned int Dimension>
-void BSplineDeformableRegistrationGenericFilter::UpdateWithDim(std::string PixelType)
-{
-
- if (m_Verbose) std::cout << "Images were detected to be "<< Dimension << "D and " << PixelType << "..." << std::endl;
-
- if(PixelType == "short") {
- if (m_Verbose) std::cout << "Launching warp in "<< Dimension <<"D and signed short..." << std::endl;
- UpdateWithDimAndPixelType<Dimension, signed short>();
- }
- // else if(PixelType == "unsigned_short"){
- // if (m_Verbose) std::cout << "Launching warp in "<< Dimension <<"D and unsigned_short..." << std::endl;
- // UpdateWithDimAndPixelType<Dimension, unsigned short>();
- // }
-
- // else if (PixelType == "unsigned_char"){
- // if (m_Verbose) std::cout << "Launching warp in "<< Dimension <<"D and unsigned_char..." << std::endl;
- // UpdateWithDimAndPixelType<Dimension, unsigned char>();
- // }
-
- // else if (PixelType == "char"){
- // if (m_Verbose) std::cout << "Launching warp in "<< Dimension <<"D and signed_char..." << std::endl;
- // UpdateWithDimAndPixelType<Dimension, signed char>();
- // }
- else {
- if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and float..." << std::endl;
- UpdateWithDimAndPixelType<Dimension, float>();
- }
-}
-
-
-
-//==============================================================================
-// Update with the number of dimensions and pixeltype
-//==============================================================================
-template<unsigned int ImageDimension, class PixelType>
-void BSplineDeformableRegistrationGenericFilter::UpdateWithDimAndPixelType()
-{
-
- //=======================================================
- // Run-time
- //=======================================================
- bool threadsGiven=m_ArgsInfo.threads_given;
- int threads=m_ArgsInfo.threads_arg;
-
- typedef itk::Image< PixelType, ImageDimension > FixedImageType;
- typedef itk::Image< PixelType, ImageDimension > MovingImageType;
- const unsigned int SpaceDimension = ImageDimension;
- typedef double TCoordRep;
-
-
- //=======================================================
- //Input
- //=======================================================
- typedef itk::ImageFileReader< FixedImageType > FixedImageReaderType;
- typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;
-
- typename FixedImageReaderType::Pointer fixedImageReader = FixedImageReaderType::New();
- typename MovingImageReaderType::Pointer movingImageReader = MovingImageReaderType::New();
-
- fixedImageReader->SetFileName( m_ArgsInfo.reference_arg );
- movingImageReader->SetFileName( m_ArgsInfo.target_arg );
- if (m_Verbose) std::cout<<"Reading images..."<<std::endl;
- fixedImageReader->Update();
- movingImageReader->Update();
-
- typename FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();
- typename MovingImageType::Pointer movingImage =movingImageReader->GetOutput();
- typename FixedImageType::RegionType fixedImageRegion = fixedImage->GetLargestPossibleRegion();
-
- // The metric region: where should the metric be CALCULATED (depends on mask)
- 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();
-
- // The transform region: where should the transform be DEFINED (depends on mask)
- typename FixedImageType::RegionType transformRegion = fixedImage->GetLargestPossibleRegion();
- typename FixedImageType::RegionType::SizeType transformRegionSize=transformRegion.GetSize();
- typename FixedImageType::RegionType::IndexType transformRegionIndex=transformRegion.GetIndex();
- typename FixedImageType::PointType transformRegionOrigin=fixedImage->GetOrigin();
-
-
- //=======================================================
- // If given, we connect a mask to the fixed image
- //======================================================
- typedef itk::ImageMaskSpatialObject< ImageDimension > MaskType;
- typename MaskType::Pointer spatialObjectMask=NULL;
-
- if (m_ArgsInfo.mask_given) {
- typedef itk::Image< unsigned char, ImageDimension > ImageMaskType;
- typedef itk::ImageFileReader< ImageMaskType > MaskReaderType;
- typename MaskReaderType::Pointer maskReader = MaskReaderType::New();
- maskReader->SetFileName(m_ArgsInfo.mask_arg);
-
- try {
- maskReader->Update();
- } catch( itk::ExceptionObject & err ) {
- std::cerr << "ExceptionObject caught while reading mask !" << std::endl;
- std::cerr << err << std::endl;
- return;
- }
- if (m_Verbose)std::cout <<"Reference image mask was read..." <<std::endl;
-
-
- // Set the image to the spatialObject
- spatialObjectMask = MaskType::New();
- spatialObjectMask->SetImage( maskReader->GetOutput() );
-
- // Find the bounding box of the "inside" label
- typedef itk::LabelStatisticsImageFilter<ImageMaskType, ImageMaskType> StatisticsImageFilterType;
- typename StatisticsImageFilterType::Pointer statisticsImageFilter=StatisticsImageFilterType::New();
- statisticsImageFilter->SetInput(maskReader->GetOutput());
- statisticsImageFilter->SetLabelInput(maskReader->GetOutput());
- statisticsImageFilter->Update();
- typename StatisticsImageFilterType::BoundingBoxType boundingBox = statisticsImageFilter->GetBoundingBox(1);
-
- // Limit the transform region to the mask
- for (unsigned int i=0; i<ImageDimension; i++) {
- transformRegionIndex[i]=boundingBox[2*i];
- transformRegionSize[i]=boundingBox[2*i+1]-boundingBox[2*i]+1;
+
+ void Execute(const itk::Object * object, const itk::EventObject & event)
+ {
+ if( !(itk::IterationEvent().CheckEvent( &event )) )
+ {
+ return;
+ }
+
+ m_Optimizer->OutputIterationInfo();
}
- transformRegion.SetSize(transformRegionSize);
- transformRegion.SetIndex(transformRegionIndex);
- fixedImage->TransformIndexToPhysicalPoint(transformRegion.GetIndex(), transformRegionOrigin);
-
- // Limit the metric region to the mask
- metricRegion=transformRegion;
- fixedImage->TransformIndexToPhysicalPoint(metricRegion.GetIndex(), metricRegionOrigin);
- }
-
-
-
- //=======================================================
- // Regions
- //=====================================================
- if (m_Verbose) {
- // Fixed image region
- std::cout<<"The fixed image has its origin at "<<fixedImage->GetOrigin()<<std::endl
- <<"The fixed image region starts at index "<<fixedImageRegion.GetIndex()<<std::endl
- <<"The fixed image region has size "<< fixedImageRegion.GetSize()<<std::endl;
-
- // Transform region
- std::cout<<"The transform has its origin at "<<transformRegionOrigin<<std::endl
- <<"The transform region will start at index "<<transformRegion.GetIndex()<<std::endl
- <<"The transform region has size "<< transformRegion.GetSize()<<std::endl;
-
- // Metric region
- std::cout<<"The metric region has its origin at "<<metricRegionOrigin<<std::endl
- <<"The metric region will start at index "<<metricRegion.GetIndex()<<std::endl
- <<"The metric region has size "<< metricRegion.GetSize()<<std::endl;
-
- }
-
-
- //=======================================================
- //Pyramids
- //=======================================================
- typedef itk::RecursiveMultiResolutionPyramidImageFilter< FixedImageType, FixedImageType> FixedImagePyramidType;
- typedef itk::RecursiveMultiResolutionPyramidImageFilter< MovingImageType, MovingImageType> MovingImagePyramidType;
- typename FixedImagePyramidType::Pointer fixedImagePyramid = FixedImagePyramidType::New();
- typename MovingImagePyramidType::Pointer movingImagePyramid = MovingImagePyramidType::New();
-
-
- // //=======================================================
- // // Rigid Transform
- // //=======================================================
- // typedef itk::Euler3DTransform <double> RigidTransformType;
- // RigidTransformType::Pointer rigidTransform=RigidTransformType::New();
-
- // if (m_ArgsInfo.rigid_given)
- // {
- // itk::Matrix<double,4,4> rigidTransformMatrix=clitk::ReadMatrix3D(m_ArgsInfo.rigid_arg);
-
- // //Set the rotation
- // itk::Matrix<double,3,3> finalRotation = clitk::GetRotationalPartMatrix3D(rigidTransformMatrix);
- // rigidTransform->SetMatrix(finalRotation);
-
- // //Set the translation
- // itk::Vector<double,3> finalTranslation = clitk::GetTranslationPartMatrix3D(rigidTransformMatrix);
- // rigidTransform->SetTranslation(finalTranslation);
-
- // }
-
-
- //=======================================================
- // BSpline Transform
- //=======================================================
- typename FixedImageType::RegionType::SizeType splineOrders ;
-
- //Default is cubic splines
- splineOrders.Fill(3);
- if (m_ArgsInfo.order_given)
- for(unsigned int i=0; i<ImageDimension; i++)
- splineOrders[i]=m_ArgsInfo.order_arg[i];
-
- // BLUT or ITK FFD
- typedef itk::Transform<TCoordRep, ImageDimension, SpaceDimension> TransformType;
- typename TransformType::Pointer transform;
- typedef itk::BSplineDeformableTransform<TCoordRep,SpaceDimension, 3> BSplineTransformType;
- typedef BSplineTransformType* BSplineTransformPointer;
- typedef clitk::BSplineDeformableTransform<TCoordRep,ImageDimension, SpaceDimension > BLUTTransformType;
- typedef BLUTTransformType* BLUTTransformPointer;
-
- // JV parameter array is passed by reference, create outside context so it exists afterwards!!!!!
- typedef typename TransformType::ParametersType ParametersType;
- ParametersType parameters;
-
-
- // CLITK BLUT transform
- if(m_ArgsInfo.wlut_flag) {
- typename BLUTTransformType::Pointer bsplineTransform = BLUTTransformType::New();
- if (m_Verbose) std::cout<<"Setting the spline orders to "<<splineOrders<<"..."<<std::endl;
- bsplineTransform->SetSplineOrders(splineOrders);
-
- //-------------------------------------------------------------------------
- // Define the region: Either the spacing or the number of CP should be given
- //-------------------------------------------------------------------------
-
- // Region
- typedef typename BSplineTransformType::RegionType RegionType;
- RegionType bsplineRegion;
- typename RegionType::SizeType gridSizeOnImage;
- typename RegionType::SizeType gridBorderSize;
- typename RegionType::SizeType totalGridSize;
-
- // Spacing
- typedef typename BSplineTransformType::SpacingType SpacingType;
- SpacingType fixedImageSpacing, chosenSpacing, adaptedSpacing;
- fixedImageSpacing = fixedImage->GetSpacing();
-
- // Only spacing given: adjust if necessary
- if (m_ArgsInfo.spacing_given && !m_ArgsInfo.control_given) {
- for(unsigned int r=0; r<ImageDimension; r++) {
- chosenSpacing[r]= m_ArgsInfo.spacing_arg[r];
- gridSizeOnImage[r] = ceil( (double) transformRegion.GetSize()[r] / ( itk::Math::Round<double>(chosenSpacing[r]/fixedImageSpacing[r]) ) );
- adaptedSpacing[r]= ( itk::Math::Round<double>(chosenSpacing[r]/fixedImageSpacing[r]) *fixedImageSpacing[r] ) ;
- }
- if (m_Verbose) std::cout<<"The chosen control point spacing "<<chosenSpacing<<"..."<<std::endl;
- if (m_Verbose) std::cout<<"The control points spacing was adapted to "<<adaptedSpacing<<"..."<<std::endl;
- if (m_Verbose) std::cout<<"The number of (internal) control points is "<<gridSizeOnImage<<"..."<<std::endl;
- }
+ OptimizerPointer m_Optimizer;
+ };
- // Only number of CP given: adjust if necessary
- else if (m_ArgsInfo.control_given && !m_ArgsInfo.spacing_given) {
- for(unsigned int r=0; r<ImageDimension; r++) {
- gridSizeOnImage[r]= m_ArgsInfo.control_arg[r];
- chosenSpacing[r]=fixedImageSpacing[r]*( (double)(transformRegion.GetSize()[r]) /
- (double)(gridSizeOnImage[r]) );
- adaptedSpacing[r]= fixedImageSpacing[r]* ceil( (double)(transformRegion.GetSize()[r] - 1) /
- (double)(gridSizeOnImage[r] - 1) );
- }
- if (m_Verbose) std::cout<<"The chosen control point spacing "<<chosenSpacing<<"..."<<std::endl;
- if (m_Verbose) std::cout<<"The control points spacing was adapted to "<<adaptedSpacing<<"..."<<std::endl;
- if (m_Verbose) std::cout<<"The number of (internal) control points is "<<gridSizeOnImage<<"..."<<std::endl;
- }
- // Spacing and number of CP given: no adjustment adjust, just warnings
- else if (m_ArgsInfo.control_given && m_ArgsInfo.spacing_given) {
- for(unsigned int r=0; r<ImageDimension; r++) {
- adaptedSpacing[r]= m_ArgsInfo.spacing_arg[r];
- gridSizeOnImage[r] = m_ArgsInfo.control_arg[r];
- if (gridSizeOnImage[r]*adaptedSpacing[r]< transformRegion.GetSize()[r]*fixedImageSpacing[r]) {
- std::cout<<"WARNING: Specified control point region ("<<gridSizeOnImage[r]*adaptedSpacing[r]
- <<"mm) does not cover the transform region ("<< transformRegion.GetSize()[r]*fixedImageSpacing[r]
- <<"mm) for dimension "<<r<<"!" <<std::endl
- <<"Specify only --spacing or --control for automatic adjustment..."<<std::endl;
- }
- if ( fmod(adaptedSpacing[r], fixedImageSpacing[r]) ) {
- std::cout<<"WARNING: Specified control point spacing for dimension "<<r
- <<" does not allow exact representation of BLUT FFD!"<<std::endl
- <<"Spacing ratio is non-integer: "<<adaptedSpacing[r]/ fixedImageSpacing[r]<<std::endl
- <<"Specify only --spacing or --control for automatic adjustment..."<<std::endl;
- }
- }
- if (m_Verbose) std::cout<<"The control points spacing was set to "<<adaptedSpacing<<"..."<<std::endl;
- if (m_Verbose) std::cout<<"The number of (internal) control points spacing is "<<gridSizeOnImage<<"..."<<std::endl;
- }
+ //==============================================================================
+ // Update with the number of dimensions
+ //==============================================================================
+ template<unsigned int Dimension>
+ void BSplineDeformableRegistrationGenericFilter::UpdateWithDim(std::string PixelType)
+ {
- //JV border size should depend on spline order
- for(unsigned int r=0; r<ImageDimension; r++) gridBorderSize[r]=splineOrders[r]; // Border for spline order = 3 ( 1 lower, 2 upper )
- totalGridSize = gridSizeOnImage + gridBorderSize;
- bsplineRegion.SetSize( totalGridSize );
- if (m_Verbose) std::cout<<"The total control point grid size was set to "<<totalGridSize<<"..."<<std::endl;
-
- // Direction
- typename FixedImageType::DirectionType gridDirection = fixedImage->GetDirection();
- SpacingType gridOriginOffset = gridDirection * adaptedSpacing;
-
- // Origin: 1 CP border for spatial dimensions
- typedef typename BSplineTransformType::OriginType OriginType;
- OriginType gridOrigin = transformRegionOrigin - gridOriginOffset;
- if (m_Verbose) std::cout<<"The control point grid origin was set to "<<gridOrigin<<"..."<<std::endl;
-
- // Set
- bsplineTransform->SetGridSpacing( adaptedSpacing );
- bsplineTransform->SetGridOrigin( gridOrigin );
- bsplineTransform->SetGridRegion( bsplineRegion );
- bsplineTransform->SetGridDirection( gridDirection );
-
- //Bulk transform
- //if (m_Verbose) std::cout<<"Setting rigid transform..."<<std::endl;
- //bsplineTransform->SetBulkTransform( rigidTransform );
-
- //Vector BSpline interpolator
- //bsplineTransform->SetOutputSpacing(fixedImage->GetSpacing());
- typename RegionType::SizeType samplingFactors;
- for (unsigned int i=0; i< ImageDimension; i++) {
- if (m_Verbose) std::cout<<"For dimension "<<i<<", the ideal sampling factor (if integer) is a multitude of "
- << (double)adaptedSpacing[i]/ (double) fixedImageSpacing[i]<<"..."<<std::endl;
- if (m_ArgsInfo.samplingFactor_given) samplingFactors[i]=m_ArgsInfo.samplingFactor_arg[i];
- else samplingFactors[i]=(int) ((double)adaptedSpacing[i]/ (double) movingImage->GetSpacing()[i]);
- if (m_Verbose) std::cout<<"Setting sampling factor "<<i<<" to "<<samplingFactors[i]<<"..."<<std::endl;
+ if (m_Verbose) std::cout << "Images were detected to be "<< Dimension << "D and " << PixelType << "..." << std::endl;
+
+ if(PixelType == "short"){
+ if (m_Verbose) std::cout << "Launching warp in "<< Dimension <<"D and signed short..." << std::endl;
+ UpdateWithDimAndPixelType<Dimension, signed short>();
}
- bsplineTransform->SetLUTSamplingFactors(samplingFactors);
-
- //initial parameters
- if (m_ArgsInfo.init_given) {
- typedef itk::ImageFileReader<typename BLUTTransformType::CoefficientImageType> CoefficientReaderType;
- typename CoefficientReaderType::Pointer coeffReader=CoefficientReaderType::New();
- coeffReader->SetFileName(m_ArgsInfo.init_arg[0]);
- coeffReader->Update();
- bsplineTransform->SetCoefficientImage(coeffReader->GetOutput());
- } else {
- //typedef typename TransformType::ParametersType ParametersType;
- const unsigned int numberOfParameters = bsplineTransform->GetNumberOfParameters();
- parameters=ParametersType( numberOfParameters );
- parameters.Fill( 0.0 );
- bsplineTransform->SetParameters( parameters );
+ // else if(PixelType == "unsigned_short"){
+ // if (m_Verbose) std::cout << "Launching warp in "<< Dimension <<"D and unsigned_short..." << std::endl;
+ // UpdateWithDimAndPixelType<Dimension, unsigned short>();
+ // }
+
+ // else if (PixelType == "unsigned_char"){
+ // if (m_Verbose) std::cout << "Launching warp in "<< Dimension <<"D and unsigned_char..." << std::endl;
+ // UpdateWithDimAndPixelType<Dimension, unsigned char>();
+ // }
+
+ // else if (PixelType == "char"){
+ // if (m_Verbose) std::cout << "Launching warp in "<< Dimension <<"D and signed_char..." << std::endl;
+ // UpdateWithDimAndPixelType<Dimension, signed char>();
+ // }
+ else {
+ if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and float..." << std::endl;
+ UpdateWithDimAndPixelType<Dimension, float>();
}
-
- // Mask
- if (spatialObjectMask) bsplineTransform->SetMask( spatialObjectMask );
-
- // Pass
- transform=bsplineTransform;
}
- //ITK BSpline transform
- else {
- typename BSplineTransformType::Pointer bsplineTransform = BSplineTransformType::New();
-
- // Define the region
- typedef typename BSplineTransformType::RegionType RegionType;
- RegionType bsplineRegion;
- typename RegionType::SizeType gridSizeOnImage;
- typename RegionType::SizeType gridBorderSize;
- typename RegionType::SizeType totalGridSize;
-
- //Set the number of control points
- for(unsigned int r=0; r<ImageDimension; r++) gridSizeOnImage[r]=m_ArgsInfo.control_arg[r];
- if (m_Verbose) std::cout<<"Setting the number of internal control points "<<gridSizeOnImage<<"..."<<std::endl;
- gridBorderSize.Fill( 3 ); // Border for spline order = 3 ( 1 lower, 2 upper )
- totalGridSize = gridSizeOnImage + gridBorderSize;
- bsplineRegion.SetSize( totalGridSize );
-
- // Spacing
- typedef typename BSplineTransformType::SpacingType SpacingType;
- SpacingType spacing = fixedImage->GetSpacing();
- typename FixedImageType::SizeType fixedImageSize = fixedImageRegion.GetSize();
- if (m_ArgsInfo.spacing_given) {
-
- for(unsigned int r=0; r<ImageDimension; r++) {
- spacing[r] =m_ArgsInfo.spacing_arg[r];
+
+
+ //==============================================================================
+ // Update with the number of dimensions and pixeltype
+ //==============================================================================
+ template<unsigned int ImageDimension, class PixelType>
+ void BSplineDeformableRegistrationGenericFilter::UpdateWithDimAndPixelType()
+ {
+
+ //=======================================================
+ // Run-time
+ //=======================================================
+ bool threadsGiven=m_ArgsInfo.threads_given;
+ int threads=m_ArgsInfo.threads_arg;
+
+ typedef itk::Image< PixelType, ImageDimension > FixedImageType;
+ typedef itk::Image< PixelType, ImageDimension > MovingImageType;
+ const unsigned int SpaceDimension = ImageDimension;
+ typedef double TCoordRep;
+
+
+ //=======================================================
+ //Input
+ //=======================================================
+ typedef itk::ImageFileReader< FixedImageType > FixedImageReaderType;
+ typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;
+
+ typename FixedImageReaderType::Pointer fixedImageReader = FixedImageReaderType::New();
+ typename MovingImageReaderType::Pointer movingImageReader = MovingImageReaderType::New();
+
+ fixedImageReader->SetFileName( m_ArgsInfo.reference_arg );
+ movingImageReader->SetFileName( m_ArgsInfo.target_arg );
+ if (m_Verbose) std::cout<<"Reading images..."<<std::endl;
+ fixedImageReader->Update();
+ movingImageReader->Update();
+
+ typename FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();
+ typename MovingImageType::Pointer movingImage =movingImageReader->GetOutput();
+ typename FixedImageType::RegionType fixedImageRegion = fixedImage->GetLargestPossibleRegion();
+
+ // The metric region: where should the metric be CALCULATED (depends on mask)
+ 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();
+
+ // The transform region: where should the transform be DEFINED (depends on mask)
+ typename FixedImageType::RegionType transformRegion = fixedImage->GetLargestPossibleRegion();
+ typename FixedImageType::RegionType::SizeType transformRegionSize=transformRegion.GetSize();
+ typename FixedImageType::RegionType::IndexType transformRegionIndex=transformRegion.GetIndex();
+ typename FixedImageType::PointType transformRegionOrigin=fixedImage->GetOrigin();
+
+
+ //=======================================================
+ // If given, we connect a mask to the fixed image
+ //======================================================
+ typedef itk::ImageMaskSpatialObject< ImageDimension > MaskType;
+ typename MaskType::Pointer spatialObjectMask=NULL;
+
+ if (m_ArgsInfo.mask_given)
+ {
+ typedef itk::Image< unsigned char, ImageDimension > ImageMaskType;
+ typedef itk::ImageFileReader< ImageMaskType > MaskReaderType;
+ typename MaskReaderType::Pointer maskReader = MaskReaderType::New();
+ maskReader->SetFileName(m_ArgsInfo.mask_arg);
+
+ try
+ {
+ maskReader->Update();
+ }
+ catch( itk::ExceptionObject & err )
+ {
+ std::cerr << "ExceptionObject caught while reading mask !" << std::endl;
+ std::cerr << err << std::endl;
+ return;
+ }
+ if (m_Verbose)std::cout <<"Reference image mask was read..." <<std::endl;
+
+
+ // Set the image to the spatialObject
+ spatialObjectMask = MaskType::New();
+ spatialObjectMask->SetImage( maskReader->GetOutput() );
+
+ // Find the bounding box of the "inside" label
+ typedef itk::LabelStatisticsImageFilter<ImageMaskType, ImageMaskType> StatisticsImageFilterType;
+ typename StatisticsImageFilterType::Pointer statisticsImageFilter=StatisticsImageFilterType::New();
+ statisticsImageFilter->SetInput(maskReader->GetOutput());
+ statisticsImageFilter->SetLabelInput(maskReader->GetOutput());
+ statisticsImageFilter->Update();
+ typename StatisticsImageFilterType::BoundingBoxType boundingBox = statisticsImageFilter->GetBoundingBox(1);
+
+ // Limit the transform region to the mask
+ for (unsigned int i=0; i<ImageDimension; i++)
+ {
+ transformRegionIndex[i]=boundingBox[2*i];
+ transformRegionSize[i]=boundingBox[2*i+1]-boundingBox[2*i]+1;
+ }
+ transformRegion.SetSize(transformRegionSize);
+ transformRegion.SetIndex(transformRegionIndex);
+ fixedImage->TransformIndexToPhysicalPoint(transformRegion.GetIndex(), transformRegionOrigin);
+
+ // Limit the metric region to the mask
+ metricRegion=transformRegion;
+ fixedImage->TransformIndexToPhysicalPoint(metricRegion.GetIndex(), metricRegionOrigin);
+
+ }
+
+
+
+ //=======================================================
+ // Regions
+ //=====================================================
+ if (m_Verbose)
+ {
+ // Fixed image region
+ std::cout<<"The fixed image has its origin at "<<fixedImage->GetOrigin()<<std::endl
+ <<"The fixed image region starts at index "<<fixedImageRegion.GetIndex()<<std::endl
+ <<"The fixed image region has size "<< fixedImageRegion.GetSize()<<std::endl;
+
+ // Transform region
+ std::cout<<"The transform has its origin at "<<transformRegionOrigin<<std::endl
+ <<"The transform region will start at index "<<transformRegion.GetIndex()<<std::endl
+ <<"The transform region has size "<< transformRegion.GetSize()<<std::endl;
+
+ // Metric region
+ std::cout<<"The metric region has its origin at "<<metricRegionOrigin<<std::endl
+ <<"The metric region will start at index "<<metricRegion.GetIndex()<<std::endl
+ <<"The metric region has size "<< metricRegion.GetSize()<<std::endl;
+
+ }
+
+
+ //=======================================================
+ //Pyramids
+ //=======================================================
+ typedef itk::RecursiveMultiResolutionPyramidImageFilter< FixedImageType, FixedImageType> FixedImagePyramidType;
+ typedef itk::RecursiveMultiResolutionPyramidImageFilter< MovingImageType, MovingImageType> MovingImagePyramidType;
+ typename FixedImagePyramidType::Pointer fixedImagePyramid = FixedImagePyramidType::New();
+ typename MovingImagePyramidType::Pointer movingImagePyramid = MovingImagePyramidType::New();
+
+
+ // //=======================================================
+ // // Rigid Transform
+ // //=======================================================
+ // typedef itk::Euler3DTransform <double> RigidTransformType;
+ // RigidTransformType::Pointer rigidTransform=RigidTransformType::New();
+
+ // if (m_ArgsInfo.rigid_given)
+ // {
+ // itk::Matrix<double,4,4> rigidTransformMatrix=clitk::ReadMatrix3D(m_ArgsInfo.rigid_arg);
+
+ // //Set the rotation
+ // itk::Matrix<double,3,3> finalRotation = clitk::GetRotationalPartMatrix3D(rigidTransformMatrix);
+ // rigidTransform->SetMatrix(finalRotation);
+
+ // //Set the translation
+ // itk::Vector<double,3> finalTranslation = clitk::GetTranslationPartMatrix3D(rigidTransformMatrix);
+ // rigidTransform->SetTranslation(finalTranslation);
+
+ // }
+
+
+ //=======================================================
+ // BSpline Transform
+ //=======================================================
+ typename FixedImageType::RegionType::SizeType splineOrders ;
+
+ //Default is cubic splines
+ splineOrders.Fill(3);
+ if (m_ArgsInfo.order_given)
+ for(unsigned int i=0; i<ImageDimension;i++)
+ splineOrders[i]=m_ArgsInfo.order_arg[i];
+
+ // BLUT or ITK FFD
+ typedef itk::Transform<TCoordRep, ImageDimension, SpaceDimension> TransformType;
+ typename TransformType::Pointer transform;
+ typedef itk::BSplineDeformableTransform<TCoordRep,SpaceDimension, 3> BSplineTransformType;
+ typedef BSplineTransformType* BSplineTransformPointer;
+ typedef clitk::BSplineDeformableTransform<TCoordRep,ImageDimension, SpaceDimension > BLUTTransformType;
+ typedef BLUTTransformType* BLUTTransformPointer;
+
+ // JV parameter array is passed by reference, create outside context so it exists afterwards!!!!!
+ typedef typename TransformType::ParametersType ParametersType;
+ ParametersType parameters;
+
+
+ // CLITK BLUT transform
+ if(m_ArgsInfo.wlut_flag)
+ {
+ typename BLUTTransformType::Pointer bsplineTransform = BLUTTransformType::New();
+ if (m_Verbose) std::cout<<"Setting the spline orders to "<<splineOrders<<"..."<<std::endl;
+ bsplineTransform->SetSplineOrders(splineOrders);
+
+ //-------------------------------------------------------------------------
+ // Define the region: Either the spacing or the number of CP should be given
+ //-------------------------------------------------------------------------
+
+ // Region
+ typedef typename BSplineTransformType::RegionType RegionType;
+ RegionType bsplineRegion;
+ typename RegionType::SizeType gridSizeOnImage;
+ typename RegionType::SizeType gridBorderSize;
+ typename RegionType::SizeType totalGridSize;
+
+ // Spacing
+ typedef typename BSplineTransformType::SpacingType SpacingType;
+ SpacingType fixedImageSpacing, chosenSpacing, adaptedSpacing;
+ fixedImageSpacing = fixedImage->GetSpacing();
+
+ // Only spacing given: adjust if necessary
+ if (m_ArgsInfo.spacing_given && !m_ArgsInfo.control_given)
+ {
+ for(unsigned int r=0; r<ImageDimension; r++)
+ {
+ chosenSpacing[r]= m_ArgsInfo.spacing_arg[r];
- gridSizeOnImage[r] = ceil( (double) transformRegion.GetSize()[r] / ( round(chosenSpacing[r]/fixedImageSpacing[r]) ) );
- adaptedSpacing[r]= ( round(chosenSpacing[r]/fixedImageSpacing[r]) *fixedImageSpacing[r] ) ;
++ gridSizeOnImage[r] = ceil( (double) transformRegion.GetSize()[r] / ( itk::Math::Round<double>(chosenSpacing[r]/fixedImageSpacing[r]) ) );
++ adaptedSpacing[r]= ( itk::Math::Round<double>(chosenSpacing[r]/fixedImageSpacing[r]) *fixedImageSpacing[r] ) ;
+ }
+ if (m_Verbose) std::cout<<"The chosen control point spacing "<<chosenSpacing<<"..."<<std::endl;
+ if (m_Verbose) std::cout<<"The control points spacing was adapted to "<<adaptedSpacing<<"..."<<std::endl;
+ if (m_Verbose) std::cout<<"The number of (internal) control points is "<<gridSizeOnImage<<"..."<<std::endl;
+ }
+
+ // Only number of CP given: adjust if necessary
+ else if (m_ArgsInfo.control_given && !m_ArgsInfo.spacing_given)
+ {
+ for(unsigned int r=0; r<ImageDimension; r++)
+ {
+ gridSizeOnImage[r]= m_ArgsInfo.control_arg[r];
+ chosenSpacing[r]=fixedImageSpacing[r]*( (double)(transformRegion.GetSize()[r]) /
+ (double)(gridSizeOnImage[r]) );
+ adaptedSpacing[r]= fixedImageSpacing[r]* ceil( (double)(transformRegion.GetSize()[r] - 1) /
+ (double)(gridSizeOnImage[r] - 1) );
+ }
+ if (m_Verbose) std::cout<<"The chosen control point spacing "<<chosenSpacing<<"..."<<std::endl;
+ if (m_Verbose) std::cout<<"The control points spacing was adapted to "<<adaptedSpacing<<"..."<<std::endl;
+ if (m_Verbose) std::cout<<"The number of (internal) control points is "<<gridSizeOnImage<<"..."<<std::endl;
+ }
+
+ // Spacing and number of CP given: no adjustment adjust, just warnings
+ else if (m_ArgsInfo.control_given && m_ArgsInfo.spacing_given)
+ {
+ for(unsigned int r=0; r<ImageDimension; r++)
+ {
+ adaptedSpacing[r]= m_ArgsInfo.spacing_arg[r];
+ gridSizeOnImage[r] = m_ArgsInfo.control_arg[r];
+ if (gridSizeOnImage[r]*adaptedSpacing[r]< transformRegion.GetSize()[r]*fixedImageSpacing[r])
+ {
+ std::cout<<"WARNING: Specified control point region ("<<gridSizeOnImage[r]*adaptedSpacing[r]
+ <<"mm) does not cover the transform region ("<< transformRegion.GetSize()[r]*fixedImageSpacing[r]
+ <<"mm) for dimension "<<r<<"!" <<std::endl
+ <<"Specify only --spacing or --control for automatic adjustment..."<<std::endl;
+ }
+ if ( fmod(adaptedSpacing[r], fixedImageSpacing[r]) )
+ {
+ std::cout<<"WARNING: Specified control point spacing for dimension "<<r
+ <<" does not allow exact representation of BLUT FFD!"<<std::endl
+ <<"Spacing ratio is non-integer: "<<adaptedSpacing[r]/ fixedImageSpacing[r]<<std::endl
+ <<"Specify only --spacing or --control for automatic adjustment..."<<std::endl;
+ }
+ }
+ if (m_Verbose) std::cout<<"The control points spacing was set to "<<adaptedSpacing<<"..."<<std::endl;
+ if (m_Verbose) std::cout<<"The number of (internal) control points spacing is "<<gridSizeOnImage<<"..."<<std::endl;
+ }
+
+ //JV border size should depend on spline order
+ for(unsigned int r=0; r<ImageDimension; r++) gridBorderSize[r]=splineOrders[r]; // Border for spline order = 3 ( 1 lower, 2 upper )
+ totalGridSize = gridSizeOnImage + gridBorderSize;
+ bsplineRegion.SetSize( totalGridSize );
+ if (m_Verbose) std::cout<<"The total control point grid size was set to "<<totalGridSize<<"..."<<std::endl;
+
+ // Direction
+ typename FixedImageType::DirectionType gridDirection = fixedImage->GetDirection();
+ SpacingType gridOriginOffset = gridDirection * adaptedSpacing;
+
+ // Origin: 1 CP border for spatial dimensions
+ typedef typename BSplineTransformType::OriginType OriginType;
+ OriginType gridOrigin = transformRegionOrigin - gridOriginOffset;
+ if (m_Verbose) std::cout<<"The control point grid origin was set to "<<gridOrigin<<"..."<<std::endl;
+
+ // Set
+ bsplineTransform->SetGridSpacing( adaptedSpacing );
+ bsplineTransform->SetGridOrigin( gridOrigin );
+ bsplineTransform->SetGridRegion( bsplineRegion );
+ bsplineTransform->SetGridDirection( gridDirection );
+
+ //Bulk transform
+ //if (m_Verbose) std::cout<<"Setting rigid transform..."<<std::endl;
+ //bsplineTransform->SetBulkTransform( rigidTransform );
+
+ //Vector BSpline interpolator
+ //bsplineTransform->SetOutputSpacing(fixedImage->GetSpacing());
+ typename RegionType::SizeType samplingFactors;
+ for (unsigned int i=0; i< ImageDimension; i++)
+ {
+ if (m_Verbose) std::cout<<"For dimension "<<i<<", the ideal sampling factor (if integer) is a multitude of "
+ << (double)adaptedSpacing[i]/ (double) fixedImageSpacing[i]<<"..."<<std::endl;
+ if (m_ArgsInfo.samplingFactor_given) samplingFactors[i]=m_ArgsInfo.samplingFactor_arg[i];
+ else samplingFactors[i]=(int) ((double)adaptedSpacing[i]/ (double) movingImage->GetSpacing()[i]);
+ if (m_Verbose) std::cout<<"Setting sampling factor "<<i<<" to "<<samplingFactors[i]<<"..."<<std::endl;
+ }
+ bsplineTransform->SetLUTSamplingFactors(samplingFactors);
+
+ //initial parameters
+ if (m_ArgsInfo.init_given)
+ {
+ typedef itk::ImageFileReader<typename BLUTTransformType::CoefficientImageType> CoefficientReaderType;
+ typename CoefficientReaderType::Pointer coeffReader=CoefficientReaderType::New();
+ coeffReader->SetFileName(m_ArgsInfo.init_arg[0]);
+ coeffReader->Update();
+ bsplineTransform->SetCoefficientImage(coeffReader->GetOutput());
+ }
+ else
+ {
+ //typedef typename TransformType::ParametersType ParametersType;
+ const unsigned int numberOfParameters = bsplineTransform->GetNumberOfParameters();
+ parameters=ParametersType( numberOfParameters );
+ parameters.Fill( 0.0 );
+ bsplineTransform->SetParameters( parameters );
+ }
+
+ // Mask
+ if (spatialObjectMask) bsplineTransform->SetMask( spatialObjectMask );
+
+ // Pass
+ transform=bsplineTransform;
}
- } else {
- for(unsigned int r=0; r<ImageDimension; r++) {
- spacing[r] *= static_cast<double>(fixedImageSize[r] - 1) /
- static_cast<double>(gridSizeOnImage[r] - 1);
- }
- }
- if (m_Verbose) std::cout<<"The control points spacing was set to "<<spacing<<"..."<<std::endl;
-
- // Direction
- typename FixedImageType::DirectionType gridDirection = fixedImage->GetDirection();
- SpacingType gridOriginOffset = gridDirection * spacing;
-
- // Origin
- typedef typename BSplineTransformType::OriginType OriginType;
- OriginType origin = fixedImage->GetOrigin();
- OriginType gridOrigin = origin - gridOriginOffset;
-
- // Set
- bsplineTransform->SetGridSpacing( spacing );
- bsplineTransform->SetGridOrigin( gridOrigin );
- bsplineTransform->SetGridRegion( bsplineRegion );
- bsplineTransform->SetGridDirection( gridDirection );
-
- // Bulk transform
- // if (m_Verbose) std::cout<<"Setting rigid transform..."<<std::endl;
- // bsplineTransform->SetBulkTransform( rigidTransform );
-
- // Initial parameters
- if (m_ArgsInfo.init_given) {
- typedef itk::ImageFileReader<typename BSplineTransformType::ImageType> CoefficientReaderType;
- typename BSplineTransformType::ImageType::Pointer coeffImages[SpaceDimension];
- for(unsigned int i=0; i<SpaceDimension; i++) {
- typename CoefficientReaderType::Pointer coeffReader=CoefficientReaderType::New();
- coeffReader->SetFileName(m_ArgsInfo.init_arg[i]);
- coeffReader->Update();
- coeffImages[i]=coeffReader->GetOutput();
- }
- bsplineTransform->SetCoefficientImage(coeffImages);
- } else {
- const unsigned int numberOfParameters = bsplineTransform->GetNumberOfParameters();
- parameters=ParametersType( numberOfParameters );
- parameters.Fill( 0.0 );
- bsplineTransform->SetParameters( parameters );
- }
-
- // Pass
- transform=bsplineTransform;
- }
+ //ITK BSpline transform
+ else
+ {
+ typename BSplineTransformType::Pointer bsplineTransform = BSplineTransformType::New();
+
+ // Define the region
+ typedef typename BSplineTransformType::RegionType RegionType;
+ RegionType bsplineRegion;
+ typename RegionType::SizeType gridSizeOnImage;
+ typename RegionType::SizeType gridBorderSize;
+ typename RegionType::SizeType totalGridSize;
+
+#if 0
+ //Set the number of control points
+ for(unsigned int r=0; r<ImageDimension; r++) gridSizeOnImage[r]=m_ArgsInfo.control_arg[r];
+ if (m_Verbose) std::cout<<"Setting the number of internal control points "<<gridSizeOnImage<<"..."<<std::endl;
+ gridBorderSize.Fill( 3 ); // Border for spline order = 3 ( 1 lower, 2 upper )
+ totalGridSize = gridSizeOnImage + gridBorderSize;
+ bsplineRegion.SetSize( totalGridSize );
+
+ // Spacing
+ typedef typename BSplineTransformType::SpacingType SpacingType;
+ SpacingType spacing = fixedImage->GetSpacing();
+ typename FixedImageType::SizeType fixedImageSize = fixedImageRegion.GetSize();
+ if (m_ArgsInfo.spacing_given)
+ {
+
+ for(unsigned int r=0; r<ImageDimension; r++)
+ {
+ spacing[r] =m_ArgsInfo.spacing_arg[r];
+ }
+ }
+ else
+ {
+ for(unsigned int r=0; r<ImageDimension; r++)
+ {
+ spacing[r] *= static_cast<double>(fixedImageSize[r] - 1) /
+ static_cast<double>(gridSizeOnImage[r] - 1);
+ }
+ }
+ if (m_Verbose) std::cout<<"The control points spacing was set to "<<spacing<<"..."<<std::endl;
+
+ // Direction
+ typename FixedImageType::DirectionType gridDirection = fixedImage->GetDirection();
+ SpacingType gridOriginOffset = gridDirection * spacing;
+
+ // Origin
+ typedef typename BSplineTransformType::OriginType OriginType;
+ OriginType origin = fixedImage->GetOrigin();
+ OriginType gridOrigin = origin - gridOriginOffset;
+
+ // Set
+ bsplineTransform->SetGridSpacing( spacing );
+ bsplineTransform->SetGridOrigin( gridOrigin );
+ bsplineTransform->SetGridRegion( bsplineRegion );
+ bsplineTransform->SetGridDirection( gridDirection );
+
+ // Bulk transform
+ // if (m_Verbose) std::cout<<"Setting rigid transform..."<<std::endl;
+ // bsplineTransform->SetBulkTransform( rigidTransform );
+
+ // Initial parameters
+ if (m_ArgsInfo.init_given)
+ {
+ typedef itk::ImageFileReader<typename BSplineTransformType::ImageType> CoefficientReaderType;
+ typename BSplineTransformType::ImageType::Pointer coeffImages[SpaceDimension];
+ for(unsigned int i=0; i<SpaceDimension; i++)
+ {
+ typename CoefficientReaderType::Pointer coeffReader=CoefficientReaderType::New();
+ coeffReader->SetFileName(m_ArgsInfo.init_arg[i]);
+ coeffReader->Update();
+ coeffImages[i]=coeffReader->GetOutput();
+ }
+ bsplineTransform->SetCoefficientImage(coeffImages);
+ }
+ else
+ {
+ const unsigned int numberOfParameters = bsplineTransform->GetNumberOfParameters();
+ parameters=ParametersType( numberOfParameters );
+ parameters.Fill( 0.0 );
+ bsplineTransform->SetParameters( parameters );
+ }
+
+ // Pass
+ transform=bsplineTransform;
+#else
+ // Spacing
+ typedef typename BSplineTransformType::SpacingType SpacingType;
+ SpacingType fixedImageSpacing, chosenSpacing, adaptedSpacing;
+ fixedImageSpacing = fixedImage->GetSpacing();
+
+ // Only spacing given: adjust if necessary
+ if (m_ArgsInfo.spacing_given && !m_ArgsInfo.control_given)
+ {
+ for(unsigned int r=0; r<ImageDimension; r++)
+ {
+ chosenSpacing[r]= m_ArgsInfo.spacing_arg[r];
+ gridSizeOnImage[r] = ceil( (double) transformRegion.GetSize()[r] / ( round(chosenSpacing[r]/fixedImageSpacing[r]) ) );
+ adaptedSpacing[r]= ( round(chosenSpacing[r]/fixedImageSpacing[r]) *fixedImageSpacing[r] ) ;
+ }
+ if (m_Verbose) std::cout<<"The chosen control point spacing "<<chosenSpacing<<"..."<<std::endl;
+ if (m_Verbose) std::cout<<"The control points spacing was adapted to "<<adaptedSpacing<<"..."<<std::endl;
+ if (m_Verbose) std::cout<<"The number of (internal) control points is "<<gridSizeOnImage<<"..."<<std::endl;
+ }
+
+ // Only number of CP given: adjust if necessary
+ else if (m_ArgsInfo.control_given && !m_ArgsInfo.spacing_given)
+ {
+ for(unsigned int r=0; r<ImageDimension; r++)
+ {
+ gridSizeOnImage[r]= m_ArgsInfo.control_arg[r];
+ chosenSpacing[r]=fixedImageSpacing[r]*( (double)(transformRegion.GetSize()[r]) /
+ (double)(gridSizeOnImage[r]) );
+ adaptedSpacing[r]= fixedImageSpacing[r]* ceil( (double)(transformRegion.GetSize()[r] - 1) /
+ (double)(gridSizeOnImage[r] - 1) );
+ }
+ if (m_Verbose) std::cout<<"The chosen control point spacing "<<chosenSpacing<<"..."<<std::endl;
+ if (m_Verbose) std::cout<<"The control points spacing was adapted to "<<adaptedSpacing<<"..."<<std::endl;
+ if (m_Verbose) std::cout<<"The number of (internal) control points is "<<gridSizeOnImage<<"..."<<std::endl;
+ }
+
+ // Spacing and number of CP given: no adjustment adjust, just warnings
+ else if (m_ArgsInfo.control_given && m_ArgsInfo.spacing_given)
+ {
+ for(unsigned int r=0; r<ImageDimension; r++)
+ {
+ adaptedSpacing[r]= m_ArgsInfo.spacing_arg[r];
+ gridSizeOnImage[r] = m_ArgsInfo.control_arg[r];
+ if (gridSizeOnImage[r]*adaptedSpacing[r]< transformRegion.GetSize()[r]*fixedImageSpacing[r])
+ {
+ std::cout<<"WARNING: Specified control point region ("<<gridSizeOnImage[r]*adaptedSpacing[r]
+ <<"mm) does not cover the transform region ("<< transformRegion.GetSize()[r]*fixedImageSpacing[r]
+ <<"mm) for dimension "<<r<<"!" <<std::endl
+ <<"Specify only --spacing or --control for automatic adjustment..."<<std::endl;
+ }
+ if ( fmod(adaptedSpacing[r], fixedImageSpacing[r]) )
+ {
+ std::cout<<"WARNING: Specified control point spacing for dimension "<<r
+ <<" does not allow exact representation of BLUT FFD!"<<std::endl
+ <<"Spacing ratio is non-integer: "<<adaptedSpacing[r]/ fixedImageSpacing[r]<<std::endl
+ <<"Specify only --spacing or --control for automatic adjustment..."<<std::endl;
+ }
+ }
+ if (m_Verbose) std::cout<<"The control points spacing was set to "<<adaptedSpacing<<"..."<<std::endl;
+ if (m_Verbose) std::cout<<"The number of (internal) control points spacing is "<<gridSizeOnImage<<"..."<<std::endl;
+ }
+
+ //JV border size should depend on spline order
+ for(unsigned int r=0; r<ImageDimension; r++) gridBorderSize[r]=splineOrders[r]; // Border for spline order = 3 ( 1 lower, 2 upper )
+ totalGridSize = gridSizeOnImage + gridBorderSize;
+ bsplineRegion.SetSize( totalGridSize );
+ if (m_Verbose) std::cout<<"The total control point grid size was set to "<<totalGridSize<<"..."<<std::endl;
+
+ // Direction
+ typename FixedImageType::DirectionType gridDirection = fixedImage->GetDirection();
+ SpacingType gridOriginOffset = gridDirection * adaptedSpacing;
+
+ // Origin: 1 CP border for spatial dimensions
+ typedef typename BSplineTransformType::OriginType OriginType;
+ OriginType gridOrigin = transformRegionOrigin - gridOriginOffset;
+ if (m_Verbose) std::cout<<"The control point grid origin was set to "<<gridOrigin<<"..."<<std::endl;
+
+ // Set
+ bsplineTransform->SetGridSpacing( adaptedSpacing );
+ bsplineTransform->SetGridOrigin( gridOrigin );
+ bsplineTransform->SetGridRegion( bsplineRegion );
+ bsplineTransform->SetGridDirection( gridDirection );
+
+ //Bulk transform
+ //if (m_Verbose) std::cout<<"Setting rigid transform..."<<std::endl;
+ //bsplineTransform->SetBulkTransform( rigidTransform );
+
+ //Vector BSpline interpolator
+ //bsplineTransform->SetOutputSpacing(fixedImage->GetSpacing());
+
+ // Initial parameters
+ if (m_ArgsInfo.init_given)
+ {
+ typedef itk::ImageFileReader<typename BSplineTransformType::ImageType> CoefficientReaderType;
+ typename BSplineTransformType::ImageType::Pointer coeffImages[SpaceDimension];
+ for(unsigned int i=0; i<SpaceDimension; i++)
+ {
+ typename CoefficientReaderType::Pointer coeffReader=CoefficientReaderType::New();
+ coeffReader->SetFileName(m_ArgsInfo.init_arg[i]);
+ coeffReader->Update();
+ coeffImages[i]=coeffReader->GetOutput();
+ }
+ bsplineTransform->SetCoefficientImage(coeffImages);
+ }
+ else
+ {
+ const unsigned int numberOfParameters = bsplineTransform->GetNumberOfParameters();
+ parameters=ParametersType( numberOfParameters );
+ parameters.Fill( 0.0 );
+ bsplineTransform->SetParameters( parameters );
+ }
+
+ // Pass
+ transform=bsplineTransform;
+#endif
- //=======================================================
- // Interpolator
- //=======================================================
- typedef clitk::GenericInterpolator<args_info_clitkBSplineDeformableRegistration, FixedImageType, TCoordRep > GenericInterpolatorType;
- typename GenericInterpolatorType::Pointer genericInterpolator=GenericInterpolatorType::New();
- genericInterpolator->SetArgsInfo(m_ArgsInfo);
- typedef itk::InterpolateImageFunction< FixedImageType, TCoordRep > InterpolatorType;
- typename InterpolatorType::Pointer interpolator=genericInterpolator->GetInterpolatorPointer();
-
-
- //=======================================================
- // Metric
- //=======================================================
- typedef clitk::GenericMetric<args_info_clitkBSplineDeformableRegistration, FixedImageType,MovingImageType > GenericMetricType;
- typename GenericMetricType::Pointer genericMetric=GenericMetricType::New();
- genericMetric->SetArgsInfo(m_ArgsInfo);
- genericMetric->SetFixedImage(fixedImage);
- genericMetric->SetFixedImageRegion(metricRegion);
- typedef itk::ImageToImageMetric< FixedImageType, MovingImageType > MetricType;
- typename MetricType::Pointer metric=genericMetric->GetMetricPointer();
- if (spatialObjectMask) metric->SetFixedImageMask( spatialObjectMask );
+ }
+
+
+ //=======================================================
+ // Interpolator
+ //=======================================================
+ typedef clitk::GenericInterpolator<args_info_clitkBSplineDeformableRegistration, FixedImageType, TCoordRep > GenericInterpolatorType;
+ typename GenericInterpolatorType::Pointer genericInterpolator=GenericInterpolatorType::New();
+ genericInterpolator->SetArgsInfo(m_ArgsInfo);
+ typedef itk::InterpolateImageFunction< FixedImageType, TCoordRep > InterpolatorType;
+ typename InterpolatorType::Pointer interpolator=genericInterpolator->GetInterpolatorPointer();
+
+
+ //=======================================================
+ // Metric
+ //=======================================================
+ typedef clitk::GenericMetric<args_info_clitkBSplineDeformableRegistration, FixedImageType,MovingImageType > GenericMetricType;
+ typename GenericMetricType::Pointer genericMetric=GenericMetricType::New();
+ genericMetric->SetArgsInfo(m_ArgsInfo);
+ genericMetric->SetFixedImage(fixedImage);
+ if (spatialObjectMask) genericMetric->SetFixedImageMask( spatialObjectMask );
+ genericMetric->SetFixedImageRegion(metricRegion);
+ typedef itk::ImageToImageMetric< FixedImageType, MovingImageType > MetricType;
+ typename MetricType::Pointer metric=genericMetric->GetMetricPointer();
#ifdef ITK_USE_OPTIMIZED_REGISTRATION_METHODS
- if (threadsGiven) metric->SetNumberOfThreads( threads );
+ if (threadsGiven) metric->SetNumberOfThreads( threads );
#else
- if (m_Verbose) std::cout<<"Not setting the number of threads (not compiled with USE_OPTIMIZED_REGISTRATION_METHODS)..."<<std::endl;
+ if (m_Verbose) std::cout<<"Not setting the number of threads (not compiled with USE_OPTIMIZED_REGISTRATION_METHODS)..."<<std::endl;
#endif
git://git.creatis.insa-lyon.fr / clitk.git / commitdiff