From: Vivien Delmon <vivien.delmon@creatis.insa-lyon.fr>
Date: Tue, 30 Aug 2011 13:24:27 +0000 (+0200)
Subject: Merge commit '488f24aa984ae24adc9458bf5fbf3b2351415742'
X-Git-Tag: v1.3.0~238^2~6
X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=commitdiff_plain;h=6194949c0beb1589904e22381b9aba1bbface172;p=clitk.git

Merge commit '488f24aa984ae24adc9458bf5fbf3b2351415742'

Conflicts:
	registration/clitkBSplineDeformableRegistrationGenericFilter.txx
---

6194949c0beb1589904e22381b9aba1bbface172
diff --cc registration/clitkBSplineDeformableRegistrationGenericFilter.txx
index c49a873,c1747a6..a0f0d7a
--- a/registration/clitkBSplineDeformableRegistrationGenericFilter.txx
+++ b/registration/clitkBSplineDeformableRegistrationGenericFilter.txx
@@@ -1,7 -1,7 +1,7 @@@
- /*=========================================================================
+ /*=========================================================================
    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
@@@ -23,680 -23,532 +23,680 @@@
  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