/*=========================================================================
  Program:   vv                     http://www.creatis.insa-lyon.fr/rio/vv

  Authors belong to: 
  - University of LYON              http://www.universite-lyon.fr/
  - Léon Bérard cancer center       http://www.centreleonberard.fr
  - CREATIS CNRS laboratory         http://www.creatis.insa-lyon.fr

  This software is distributed WITHOUT ANY WARRANTY; without even
  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  PURPOSE.  See the copyright notices for more information.

  It is distributed under dual licence

  - BSD        See included LICENSE.txt file
  - CeCILL-B   http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
===========================================================================**/
#ifndef __clitkShapedBLUTSpatioTemporalDeformableTransform_txx
#define __clitkShapedBLUTSpatioTemporalDeformableTransform_txx
#include "clitkShapedBLUTSpatioTemporalDeformableTransform.h"

//itk
#include "itkContinuousIndex.h"
#include "itkImageRegionIterator.h"
#include "itkImageRegionConstIterator.h"
#include "itkImageRegionConstIteratorWithIndex.h"
#include "itkIdentityTransform.h"

namespace clitk
{

  // Constructor with default arguments
  template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
  ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>::ShapedBLUTSpatioTemporalDeformableTransform():Superclass(0)
  {
    unsigned int i;
    
    //JV default spline order
    for ( i=0;i<InputDimension; i++)
      m_SplineOrders[i]=3;

    //JV default LUTSamplingfactor
    for ( i=0;i<InputDimension; i++)
      m_LUTSamplingFactors[i]=20;

    for (i=0;i<InputDimension; i++)
      m_SupportSize[i] = m_SplineOrders[i]+1;

    //Instantiate a Vector Bspline Interpolator
    m_VectorInterpolator =VectorInterpolatorType::New();
    m_VectorInterpolator->SetLUTSamplingFactors(m_LUTSamplingFactors);
    m_VectorInterpolator->SetSplineOrders(m_SplineOrders);

    // Set Bulk transform to NULL (no bulk performed)
    m_BulkTransform = NULL;

    // Mask
    m_Mask=NULL;

    // Shape
    m_TransformShape=0;

    // Default grid size is zero
    m_NullSize.Fill(0);
    m_NullIndex.Fill(0);
  
    //JV region containing the parameters
    m_GridRegion.SetSize( m_NullSize);
    m_GridRegion.SetIndex( m_NullIndex );

    //JV use second region over the images
    m_PaddedGridRegion.SetSize(m_NullSize);
    m_PaddedGridRegion.SetIndex(m_NullIndex);
    
    //JV Maintain two origins
    m_GridOrigin.Fill( 0.0 );  // default origin is all zeros
    m_PaddedGridOrigin.Fill( 0.0 );  // default origin is all zeros
    
    m_GridSpacing.Fill( 1.0 ); // default spacing is all ones
    m_GridDirection.SetIdentity(); // default spacing is all ones

    m_InternalParametersBuffer = ParametersType(0);
    // Make sure the parameters pointer is not NULL after construction.
    m_InputParametersPointer = &m_InternalParametersBuffer;

    // Initialize coeffient images
    m_WrappedImage = CoefficientImageType::New();
    m_WrappedImage->SetRegions( m_GridRegion );
    m_WrappedImage->SetOrigin( m_GridOrigin.GetDataPointer() );
    m_WrappedImage->SetSpacing( m_GridSpacing.GetDataPointer() );
    m_WrappedImage->SetDirection( m_GridDirection );

    // JV Initialize the padded version
    m_PaddedCoefficientImage = CoefficientImageType::New();
    m_PaddedCoefficientImage->SetRegions( m_PaddedGridRegion );
    m_PaddedCoefficientImage->SetOrigin( m_GridOrigin.GetDataPointer() );
    m_PaddedCoefficientImage->SetSpacing( m_GridSpacing.GetDataPointer() );
    m_PaddedCoefficientImage->SetDirection( m_GridDirection );

    m_CoefficientImage = NULL;
  
    // Variables for computing interpolation
    for (i=0; i <InputDimension;i++)
      {
	m_Offset[i] = m_SplineOrders[i] / 2;
	if ( m_SplineOrders[i] % 2 ) 
	  {
	    m_SplineOrderOdd[i] = true;
	  }
	else
	  {
	    m_SplineOrderOdd[i] = false;
	  }
      }

    //JV padded should be used when checking
    m_ValidRegion = m_PaddedGridRegion;
	
    // Initialize jacobian images 
    for (i=0; i <OutputDimension;i++)
      {
	m_JacobianImage[i] = JacobianImageType::New();
	m_JacobianImage[i]->SetRegions( m_GridRegion );
	m_JacobianImage[i]->SetOrigin( m_GridOrigin.GetDataPointer() );
	m_JacobianImage[i]->SetSpacing( m_GridSpacing.GetDataPointer() );
	m_JacobianImage[i]->SetDirection( m_GridDirection );
      }
    
    /** Fixed Parameters store the following information:
     *     Grid Size
     *     Grid Origin
     *     Grid Spacing
     *     Grid Direction */
    //JV we add the splineOrders, LUTsamplingfactor, m_Mask  and m_BulkTransform
    /*
      Spline orders
      Sampling factors
      m_Mask
      m_BulkTransform
      The size of these is equal to the  NInputDimensions
    *********************************************************/
    this->m_FixedParameters.SetSize ( NInputDimensions * (NInputDimensions + 5)+3 );
    this->m_FixedParameters.Fill ( 0.0 );
    for ( i=0; i<NInputDimensions; i++)
      {
	this->m_FixedParameters[2*NInputDimensions+i] = m_GridSpacing[i];
      }
    for (unsigned int di=0; di<NInputDimensions; di++)
      {
	for (unsigned int dj=0; dj<NInputDimensions; dj++)
	  {
	    this->m_FixedParameters[3*NInputDimensions+(di*NInputDimensions+dj)] = m_GridDirection[di][dj];
	  }
      }

    //JV add splineOrders
    for ( i=0; i<NInputDimensions; i++)
      {
	this->m_FixedParameters[ ( (3+NInputDimensions)*NInputDimensions)+i] = (this->GetSplineOrders())[i];
      }
    
    //JV add LUTsamplingFactors
    for ( i=0; i<NInputDimensions; i++)
      {
	this->m_FixedParameters[( (4+NInputDimensions)*NInputDimensions)+i ] = this->GetLUTSamplingFactors()[i];
      }

    // JV add the mask pointer
    this->m_FixedParameters[( (5+NInputDimensions)*NInputDimensions)]=(double)((size_t)m_Mask);

    // JV add the bulkTransform pointer
    this->m_FixedParameters[( (5+NInputDimensions)*NInputDimensions) +1]=(double)((size_t)m_BulkTransform);

    // JV add the Transform shape
    this->m_FixedParameters[( (5+NInputDimensions)*NInputDimensions) +2]=(double)(m_TransformShape);


    // Calculate the PointToIndex matrices
    DirectionType scale;
    for( unsigned int i=0; i<OutputDimension; i++)
      {
	scale[i][i] = m_GridSpacing[i];
      }

    m_IndexToPoint = m_GridDirection * scale;
    m_PointToIndex = m_IndexToPoint.GetInverse();
 
    m_LastJacobianIndex = m_ValidRegion.GetIndex();


    //Weights to be used for the BC
    m_Weights.resize(4);
    std::vector<double> weights(5);
    for (i=0; i<4; i++) 
      weights[0]=0;
    // NN
    weights[0]=1;
    weights[1]=1./2;
    weights[2]=0;
    weights[3]=0;
    weights[4]=0;
    m_Weights[0]=weights;

    // Linear
    weights[0]=1;
    weights[1]=1./2;
    weights[2]=0;
    weights[3]=0;
    weights[4]=0;
    m_Weights[1]=weights;

    // quadratic
    weights[0]=3./4.;
    weights[1]=1./2.;
    weights[2]=1./8.;
    weights[3]=0;
    weights[4]=0;
    m_Weights[2]=weights;

    // cubic
    weights[0]=2./3.;
    weights[1]=23./48. ;
    weights[2]=1./6.;
    weights[3]=1./48.;
    weights[4]=0;
    m_Weights[3]=weights;

    // Update the WeightRatios
    m_WeightRatio.resize(4);
    for (unsigned int i=0; i<4; i++)
      {
	m_WeightRatio[i].resize(4);
	for (unsigned int j=0; j<4; j++)
	  m_WeightRatio[i][j]=m_Weights[m_SplineOrders[InputDimension-1]][i]/ m_Weights[m_SplineOrders[InputDimension-1]][j];
      }

    //JV initialize some variables for jacobian calculation
    m_SupportRegion.SetSize(m_SupportSize);
    m_SupportIndex.Fill(0);
    m_SupportRegion.SetIndex(m_SupportIndex);
    for (  i = 0; i < SpaceDimension ; i++ ) 
      m_ZeroVector[i]=itk::NumericTraits<JacobianValueType>::Zero;
    
    m_InitialOffset=0;
    m_FirstRegion.SetSize(m_NullSize);
    m_FirstRegion.SetIndex(m_NullIndex);
    m_SecondRegion.SetSize(m_NullSize);
    m_SecondRegion.SetIndex(m_NullIndex);
    m_ThirdRegion.SetSize(m_NullSize);
    m_ThirdRegion.SetIndex(m_NullIndex);
        
    m_BCValues.resize(0);
    m_BCRegions.resize(0);
    m_BCSize=0;
    m_BC2Values.resize(0);
    m_BC2Regions.resize(0);
    m_BC2Size=0;
    m_BC3Values.resize(0);
    m_BC3Regions.resize(0);
    m_BC3Size=0;


    for (  i = 0; i < SpaceDimension ; i++ )
      {
  	m_FirstIterator[i]= IteratorType( m_JacobianImage[i], m_FirstRegion);
    	m_SecondIterator[i]= IteratorType( m_JacobianImage[i], m_SecondRegion);
    	m_ThirdIterator[i]= IteratorType( m_JacobianImage[i], m_ThirdRegion);
      	m_BCIterators[i].resize(0);
	m_BC2Iterators[i].resize(0);
	m_BC3Iterators[i].resize(0);
      }

    this->Modified();

  }
    

  // Destructor
  template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
  ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
  ::~ShapedBLUTSpatioTemporalDeformableTransform()
  {

  }


  // JV set Spline Order
  template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
  void
  ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
  ::SetSplineOrder(const unsigned int & splineOrder)
  {
    SizeType splineOrders;
    for (unsigned int i=0;i<InputDimension;i++)splineOrders[i]=splineOrder;

    this->SetSplineOrders(splineOrders);
  }
   

  // JV set Spline Orders
  template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
  void
  ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
  ::SetSplineOrders(const SizeType & splineOrders)
  {
    if(m_SplineOrders!=splineOrders)
      {
	m_SplineOrders=splineOrders;
	
	//update the interpolation function
	m_VectorInterpolator->SetSplineOrders(m_SplineOrders);
	
	//update the varaibles for computing interpolation
	for (unsigned int i=0; i <InputDimension;i++)
	  {
	    m_SupportSize[i] = m_SplineOrders[i]+1;
	    m_Offset[i] = m_SplineOrders[i] / 2;
	    
	    if ( m_SplineOrders[i] % 2 ) 
	      {
		m_SplineOrderOdd[i] = true;
	      }
	    else
	      {
		m_SplineOrderOdd[i] = false;
	      }
	  }

	//SupportSize is updated!, update regions
	//JV initialize some variables for jacobian calculation
	m_SupportRegion.SetSize(m_SupportSize);
	m_SupportIndex.Fill(0);
	m_SupportRegion.SetIndex(m_SupportIndex);
	
	// Update the WeightRatios
	for (unsigned int i=0; i<4; i++)
	  for (unsigned int j=0; j<4; j++)
	    m_WeightRatio[i][j]=m_Weights[m_SplineOrders[InputDimension-1]][i]/ m_Weights[m_SplineOrders[InputDimension-1]][j];

	this->Modified();
      }
  }
   

  // JV set sampling factor
  template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
  void
  ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
  ::SetLUTSamplingFactor( const int & samplingFactor)
  {
    SizeType samplingFactors;
    for (unsigned int i=0; i<NInputDimensions; i++)
      samplingFactors[i]=samplingFactor;
    
    //update the interpolation function
    this->SetLUTSamplingFactors(samplingFactors);
  }


  // JV set sampling factors
  template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
  void
  ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
  ::SetLUTSamplingFactors( const SizeType & samplingFactors)
  {
    if(m_LUTSamplingFactors!=samplingFactors)
      {
	for (unsigned int i=0; i<NInputDimensions; i++)
	  m_LUTSamplingFactors[i]=samplingFactors[i];
	
	//update the interpolation function
	m_VectorInterpolator->SetLUTSamplingFactors(m_LUTSamplingFactors);
	
	this->Modified();
      }
  }


  // Get the number of parameters
  template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
  typename ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>::NumberOfParametersType
  ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
  ::GetNumberOfParameters(void) const
  {

    // The number of parameters equal SpaceDimension * number of
    // of pixels in the grid region.
    return ( static_cast<unsigned int>( SpaceDimension ) *
	     static_cast<unsigned int>( m_GridRegion.GetNumberOfPixels() ) );

  }


  // Get the padded number of parameters
  template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
  unsigned int
  ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
  ::GetPaddedNumberOfParameters(void) const
  {

    // The number of parameters equal SpaceDimension * number of
    // of pixels in the grid region.
    return ( static_cast<unsigned int>( SpaceDimension ) *
	     static_cast<unsigned int>( m_PaddedGridRegion.GetNumberOfPixels() ) );

  }



  // Get the number of parameters per dimension
  template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
  unsigned int
  ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
  ::GetNumberOfParametersPerDimension(void) const
  {
    // The number of parameters per dimension equal number of
    // of pixels in the grid region.
    return ( static_cast<unsigned int>( m_GridRegion.GetNumberOfPixels() ) );

  }


  // Set the grid region
  template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
  void
  ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
  ::SetGridRegion( const RegionType & region )
  {
    if ( m_GridRegion != region )
      {
	m_GridRegion = region;
	m_PaddedGridRegion=region;
	
	// JV set the padded region
	typename CoefficientImageType::RegionType::SizeType paddedSize= region.GetSize();
	
	//JV   size dependes on shape
	switch (m_TransformShape)
	  {
	  case 0:
	  case 1:
	    paddedSize[InputDimension-1]+=4;
	    break;
	  case 2:
	  case 3:
	    paddedSize[InputDimension-1]+=4;
	    break;
	  case 4:
	  case 5:
	    paddedSize[InputDimension-1]+=3;
	    break;
	  case 6:
	  case 7:
	    paddedSize[InputDimension-1]+=2;
	    break;
	  case 8:
	  case 9:
	    paddedSize[InputDimension-1]+=3;
	    break;
	  default:
	    paddedSize[InputDimension-1]=+1;
	  }
	m_PaddedGridRegion.SetSize(paddedSize);
	
	// Set regions for each coefficient and jacobian image
	m_WrappedImage->SetRegions( m_GridRegion );
	m_PaddedCoefficientImage->SetRegions( m_PaddedGridRegion );
	m_PaddedCoefficientImage->Allocate();
	for (unsigned int j=0; j <OutputDimension;j++)
	  {
	    m_JacobianImage[j]->SetRegions( m_GridRegion );
	  }
	
	// JV used the padded version for the valid region
	// Set the valid region
	// If the grid spans the interval [start, last].
	// The valid interval for evaluation is [start+offset, last-offset]
	// when spline order is even.
	// The valid interval for evaluation is [start+offset, last-offset)
	// when spline order is odd.
	// Where offset = std::floor(spline / 2 ).
	// Note that the last pixel is not included in the valid region
	// with odd spline orders.
	typename RegionType::SizeType size = m_PaddedGridRegion.GetSize();
	typename RegionType::IndexType index = m_PaddedGridRegion.GetIndex();
	for ( unsigned int j = 0; j < NInputDimensions; j++ )
	  {
	    index[j] += 
	      static_cast< typename RegionType::IndexValueType >( m_Offset[j] );
	    size[j] -= 
	      static_cast< typename RegionType::SizeValueType> ( 2 * m_Offset[j] );
	    m_ValidRegionLast[j] = index[j] +
	      static_cast< typename RegionType::IndexValueType >( size[j] ) - 1;
	  }
	m_ValidRegion.SetSize( size );
	m_ValidRegion.SetIndex( index );
	
	// If we are using the default parameters, update their size and set to identity.
	// Input parameters point to internal buffer => using default parameters.
	if (m_InputParametersPointer == &m_InternalParametersBuffer)
	  {
	    // Check if we need to resize the default parameter buffer.
	    if ( m_InternalParametersBuffer.GetSize() != this->GetNumberOfParameters() )
	      {
		m_InternalParametersBuffer.SetSize( this->GetNumberOfParameters() );
		// Fill with zeros for identity.
		m_InternalParametersBuffer.Fill( 0 );
	      }
	  }

	this->Modified();
      }
  }


  // Set the grid spacing
  template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
  void
  ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
  ::SetGridSpacing( const SpacingType & spacing )
  {
    if ( m_GridSpacing != spacing )
      {
	m_GridSpacing = spacing;

	// Set spacing for each coefficient and jacobian image
	m_WrappedImage->SetSpacing( m_GridSpacing.GetDataPointer() );
	m_PaddedCoefficientImage->SetSpacing( m_GridSpacing.GetDataPointer() );
	for (unsigned int j=0; j <OutputDimension; j++) m_JacobianImage[j]->SetSpacing( m_GridSpacing.GetDataPointer() );
	
	// Set scale
	DirectionType scale;
	for( unsigned int i=0; i<OutputDimension; i++)
	  {
	    scale[i][i] = m_GridSpacing[i];
	  }

	m_IndexToPoint = m_GridDirection * scale;
	m_PointToIndex = m_IndexToPoint.GetInverse();

	this->Modified();
      }

  }

  // Set the grid direction
  template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
  void
  ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
  ::SetGridDirection( const DirectionType & direction )
  {
    if ( m_GridDirection != direction )
      {
	m_GridDirection = direction;

	// Set direction for each coefficient and jacobian image
	m_WrappedImage->SetDirection( m_GridDirection );
	m_PaddedCoefficientImage->SetDirection( m_GridDirection );
	for (unsigned int j=0; j <OutputDimension; j++) m_JacobianImage[j]->SetDirection( m_GridDirection );
	
	// Set scale
	DirectionType scale;
	for( unsigned int i=0; i<OutputDimension; i++)
	  {
	    scale[i][i] = m_GridSpacing[i];
	  }

	m_IndexToPoint = m_GridDirection * scale;
	m_PointToIndex = m_IndexToPoint.GetInverse();

	this->Modified();
      }

  }


  // Set the grid origin
  template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
  void
  ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
  ::SetGridOrigin( const OriginType& origin )
  {
    if( m_GridOrigin!=origin)
      {
	m_GridOrigin = origin;
	
	// JV The origin depends on the shape 
	switch (m_TransformShape)
	  {
	  case 0:
	  case 1:
	    m_PaddedGridOrigin=origin;
	    m_PaddedGridOrigin[InputDimension-1]=origin[InputDimension-1]-2* m_GridSpacing[InputDimension-1];
	    break;
	  case 2:
	  case 3:
	    m_PaddedGridOrigin=origin;
	    m_PaddedGridOrigin[InputDimension-1]=origin[InputDimension-1]-1* m_GridSpacing[InputDimension-1];
	    break;
	  case 4:
	  case 5:
	    m_PaddedGridOrigin=origin;
	    m_PaddedGridOrigin[InputDimension-1]=origin[InputDimension-1]-1* m_GridSpacing[InputDimension-1];
	    break;
	  case 6:
	  case 7:
	    m_PaddedGridOrigin=origin;
	    break;
	  case 8:
	  case 9:
	    m_PaddedGridOrigin=origin;
	    m_PaddedGridOrigin[InputDimension-1]=origin[InputDimension-1]-1* m_GridSpacing[InputDimension-1];
	    break;
	  default: 
	    m_PaddedGridOrigin=origin;
	  }
    
	// Set origin for each coefficient and jacobianimage
	m_WrappedImage->SetOrigin( m_GridOrigin.GetDataPointer() );
	m_PaddedCoefficientImage->SetOrigin( m_PaddedGridOrigin.GetDataPointer() );
	for (unsigned int j=0; j <OutputDimension; j++) m_JacobianImage[j]->SetOrigin( m_GridOrigin.GetDataPointer() );
	
	this->Modified();
      }
  }


  // Set the parameters
  template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
  void
  ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
  ::SetIdentity()
  {
    if( m_InputParametersPointer )
      {
	ParametersType * parameters =
	  const_cast<ParametersType *>( m_InputParametersPointer );
	parameters->Fill( 0.0 );
	this->Modified();
      }
    else 
      {
	itkExceptionMacro( << "Input parameters for the spline haven't been set ! "
			   << "Set them using the SetParameters or SetCoefficientImage method first." );
      }
  }


  // Set the parameters
  template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
  void
  ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
  ::SetParameters( const ParametersType & parameters )
  {

    // Check if the number of parameters match the
    // Expected number of parameters
    if ( parameters.Size() != this->GetNumberOfParameters() )
      {
	itkExceptionMacro(<<"Mismatched between parameters size "
			  << parameters.size() 
			  << " and region size "
			  << m_GridRegion.GetNumberOfPixels() );
      }

    // Clean up buffered parameters
    m_InternalParametersBuffer = ParametersType( 0 );

    // Keep a reference to the input parameters
    m_InputParametersPointer = &parameters;

    // Wrap flat array as images of coefficients
    this->WrapAsImages();

    //JV Set padded input to vector interpolator
    m_VectorInterpolator->SetInputImage(this->GetPaddedCoefficientImage());

    // Modified is always called since we just have a pointer to the
    // parameters and cannot know if the parameters have changed.
    this->Modified();
  }


  // Set the Fixed Parameters
  template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
  void
  ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
  ::SetFixedParameters( const ParametersType & parameters )
  {
 
    // JV number should be exact, no defaults for spacing
    if ( parameters.Size() != NInputDimensions * (5 + NInputDimensions)+3 )
      {
	itkExceptionMacro(<< "Mismatched between parameters size "
			  << parameters.size() 
			  << " and number of fixed parameters "
			  << NInputDimensions * (5 + NInputDimensions)+3 );
      }
    /********************************************************* 
    Fixed Parameters store the following information:
        Grid Size
        Grid Origin
        Grid Spacing
        Grid Direction */
    // JV we add the splineOrders, LUTsamplingfactor, mask pointer and bulktransform pointer
    /*
      Spline orders
      Sampling factors
      m_Mask
      m_BulkTransform
      m_TransformShape

      The size of these is equal to the  NInputDimensions
    *********************************************************/
  
    /** Set the Grid Parameters */
    SizeType   gridSize;
    for (unsigned int i=0; i<NInputDimensions; i++)
      {
	gridSize[i] = static_cast<int> (parameters[i]);
      }
    RegionType bsplineRegion;
    bsplineRegion.SetSize( gridSize );
  
    /** Set the Origin Parameters */
    OriginType origin;
    for (unsigned int i=0; i<NInputDimensions; i++)
      {
	origin[i] = parameters[NInputDimensions+i];
      }
  
    /** Set the Spacing Parameters */
    SpacingType spacing;
    for (unsigned int i=0; i<NInputDimensions; i++)
      {
	spacing[i] = parameters[2*NInputDimensions+i];
      }

    /** Set the Direction Parameters */
    DirectionType direction;
    for (unsigned int di=0; di<NInputDimensions; di++)
      {
	for (unsigned int dj=0; dj<NInputDimensions; dj++)
	  {
	    direction[di][dj] = parameters[3*NInputDimensions+(di*NInputDimensions+dj)];
	  }
      }

    //JV add the spline orders
    SizeType splineOrders;
    for (unsigned int i=0; i<NInputDimensions; i++)
      {
	splineOrders[i]=parameters[(3+NInputDimensions)*NInputDimensions+i];
      }

    //JV add the LUT sampling factor
    SizeType  samplingFactors;
    for (unsigned int i=0; i<NInputDimensions; i++)
      {
	samplingFactors[i]=parameters[(4+NInputDimensions)*NInputDimensions+i];
      }

    //JV add the MaskPointer
    m_Mask=((MaskPointer)( (size_t)parameters[(5+NInputDimensions)*NInputDimensions]));

    //JV add the MaskPointer
    m_BulkTransform=((BulkTransformPointer)( (size_t)parameters[(5+NInputDimensions)*NInputDimensions+1]));

    //JV add the TransformShape
    m_TransformShape=((unsigned int)parameters[(5+NInputDimensions)*NInputDimensions+2]);

    // Set the members 
    this->SetSplineOrders( splineOrders );
    this->SetGridSpacing( spacing );
    this->SetGridDirection( direction );
    this->SetGridOrigin( origin );
    this->SetGridRegion( bsplineRegion );
    this->SetLUTSamplingFactors( samplingFactors );
      
  }


  // Wrap flat parameters as images
  template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
  void
  ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
  ::WrapAsImages()
  {
    //JV Wrap parameter array in vectorial image, changed parameter order: A1x A1y A1z, A2x ....
    PixelType * dataPointer =reinterpret_cast<PixelType *>( const_cast<double *>(m_InputParametersPointer->data_block() )) ;
    unsigned int numberOfPixels = m_GridRegion.GetNumberOfPixels();
    
    m_WrappedImage->GetPixelContainer()->SetImportPointer( dataPointer,numberOfPixels);//InputDimension
    m_CoefficientImage = m_WrappedImage;

    //=====================================
    //JV Create padded structure adding BC
    //=====================================
    PadCoefficientImage();
        
    //=====================================
    //JV Wrap jacobian into OutputDimension X Vectorial images
    //=====================================
    this->m_SharedDataBSplineJacobian.set_size( OutputDimension, this->GetNumberOfParameters() );

    // Use memset to set the memory
    // JV four rows of three comps of parameters
    JacobianPixelType * jacobianDataPointer = reinterpret_cast<JacobianPixelType *>(this->m_SharedDataBSplineJacobian.data_block());
    memset(jacobianDataPointer, 0,  OutputDimension*numberOfPixels*sizeof(JacobianPixelType));

    for (unsigned int j=0; j<OutputDimension; j++)
      {
	m_JacobianImage[j]->GetPixelContainer()->
	  SetImportPointer( jacobianDataPointer, numberOfPixels );
	jacobianDataPointer += numberOfPixels;
      }

    // Reset the J parameters 
    m_LastJacobianIndex = m_ValidRegion.GetIndex();
    m_FirstRegion.SetSize(m_NullSize);
    m_SecondRegion.SetSize(m_NullSize);
    for ( unsigned int j = 0; j < SpaceDimension ; j++ )
      {
	m_FirstIterator[j]= IteratorType( m_JacobianImage[j], m_FirstRegion);
	m_SecondIterator[j]= IteratorType( m_JacobianImage[j], m_SecondRegion);
      }

    m_BCValues.resize(0);
    m_BCRegions.resize(0);
    m_BCSize=0;
    m_BC2Values.resize(0);
    m_BC2Regions.resize(0);
    m_BC2Size=0;
    m_BC3Values.resize(0);
    m_BC3Regions.resize(0);
    m_BC3Size=0;

  }


  // Set the parameters by value
  template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
  void
  ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
  ::SetParametersByValue( const ParametersType & parameters )
  {

    // Check if the number of parameters match the
    // Expected number of parameters
    if ( parameters.Size() != this->GetNumberOfParameters() )
      {
	itkExceptionMacro(<<"Mismatched between parameters size "
			  << parameters.size() 
			  << " and region size "
			  << m_GridRegion.GetNumberOfPixels() );
      }

    // Copy it
    m_InternalParametersBuffer = parameters;
    m_InputParametersPointer = &m_InternalParametersBuffer;

    // Wrap flat array as images of coefficients
    this->WrapAsImages();

    //JV Set padded input to vector interpolator
    m_VectorInterpolator->SetInputImage(this->GetPaddedCoefficientImage());

    // Modified is always called since we just have a pointer to the
    // Parameters and cannot know if the parameters have changed.
    this->Modified();

  }

  // Get the parameters
  template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
  const 
  typename ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
  ::ParametersType &
  ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
  ::GetParameters( void ) const
  {
    /** NOTE: For efficiency, this class does not keep a copy of the parameters - 
     * it just keeps pointer to input parameters. 
     */
    if (NULL == m_InputParametersPointer)
      {
	itkExceptionMacro( <<"Cannot GetParameters() because m_InputParametersPointer is NULL. Perhaps SetCoefficientImage() has been called causing the NULL pointer." );
      }

    return (*m_InputParametersPointer);
  }


  // Get the parameters
  template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
  const 
  typename ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
  ::ParametersType &
  ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
  ::GetFixedParameters( void ) const
  {
    RegionType resRegion = this->GetGridRegion(  );
  
    for (unsigned int i=0; i<NInputDimensions; i++)
      {
	this->m_FixedParameters[i] = (resRegion.GetSize())[i];
      }
    for (unsigned int i=0; i<NInputDimensions; i++)
      {
	this->m_FixedParameters[NInputDimensions+i] = (this->GetGridOrigin())[i];
      } 
    for (unsigned int i=0; i<NInputDimensions; i++)
      {
	this->m_FixedParameters[2*NInputDimensions+i] =  (this->GetGridSpacing())[i];
      }
    for (unsigned int di=0; di<NInputDimensions; di++)
      {
	for (unsigned int dj=0; dj<NInputDimensions; dj++)
	  {
	    this->m_FixedParameters[3*NInputDimensions+(di*NInputDimensions+dj)] = (this->GetGridDirection())[di][dj];
	  }
      }

    //JV add splineOrders
    for (unsigned int i=0; i<NInputDimensions; i++)
      {
	this->m_FixedParameters[(3+NInputDimensions)*NInputDimensions+i] = (this->GetSplineOrders())[i];
      }
    
    //JV add LUTsamplingFactor
    for (unsigned int i=0; i<NInputDimensions; i++)
      {
	this->m_FixedParameters[(4+NInputDimensions)*NInputDimensions+i] = (this->GetLUTSamplingFactors())[i];
      }
    
    //JV add the mask
    this->m_FixedParameters[(5+NInputDimensions)*NInputDimensions]=(double)((size_t) m_Mask);

    //JV add the bulktransform pointer
    this->m_FixedParameters[(5+NInputDimensions)*NInputDimensions+1]=(double)((size_t) m_BulkTransform);

    //JV add the transform shape
    this->m_FixedParameters[(5+NInputDimensions)*NInputDimensions+2]=(double)( m_TransformShape);
    
    return (this->m_FixedParameters);
  }

  
  // Set the B-Spline coefficients using input images
  template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
  void 
  ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
  ::SetCoefficientImage( CoefficientImagePointer image )
  {
    this->SetGridSpacing( image->GetSpacing() );
    this->SetGridOrigin( image->GetOrigin() );
    this->SetGridDirection( image->GetDirection() );
    this->SetGridRegion( image->GetBufferedRegion() );
    m_CoefficientImage = image;
    
    //JV
    //   m_WrappedImage=m_CoefficientImage;

    // Update the interpolator
    this->PadCoefficientImage();
    m_VectorInterpolator->SetInputImage(this->GetPaddedCoefficientImage());

    // Clean up buffered parameters
    m_InternalParametersBuffer = ParametersType( 0 );
    m_InputParametersPointer  = NULL;
      
  }  


  //   // Set the B-Spline coefficients using input images
  //   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
  //   void 
  //   ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
  //   ::SetPaddedCoefficientImage( CoefficientImagePointer image )
  //   {  
  //     //JV modify the region
  //     typename CoefficientImageType::RegionType region=image->GetBufferedRegion();
  //     typename CoefficientImageType::RegionType::SizeType size=region.GetSize();
  //     size[InputDimension-1]-=2;
  //     region.SetSize(size);
    
  //     //Set properties
  //     this->SetGridRegion( region );
  //     this->SetGridSpacing( image->GetSpacing() );
  //     this->SetGridDirection( image->GetDirection() );
  //     this->SetGridOrigin( image->GetOrigin() );
  //     m_PaddedCoefficientImage = image;
  //     this->ExtractCoefficientImage();
  //     m_VectorInterpolator->SetInputImage(this->GetPaddedCoefficientImage());
    
  //     // Clean up buffered parameters
  //     m_InternalParametersBuffer = ParametersType( 0 );
  //     m_InputParametersPointer  = NULL;
      
  //   } 

  // Set the B-Spline coefficients using input images
  template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
  typename   ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>::CoefficientImageType::Pointer
  ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
  ::ExtractTemporalRow(const typename CoefficientImageType::Pointer& coefficientImage, unsigned int temporalIndex)
  {
    // Region
    typename   CoefficientImageType::RegionType sourceRegion=coefficientImage->GetLargestPossibleRegion();
    sourceRegion.SetSize(InputDimension-1, 1);
    sourceRegion.SetIndex(InputDimension-1, temporalIndex);

    // Extract
    typedef clitk::ExtractImageFilter<CoefficientImageType, CoefficientImageType> ExtractImageFilterType;
    typename ExtractImageFilterType::Pointer extract=ExtractImageFilterType::New();
    extract->SetInput(coefficientImage);
    extract->SetExtractionRegion(sourceRegion);
    extract->Update();
    typename CoefficientImageType::Pointer row= extract->GetOutput();

    // Set index to zero
    sourceRegion.SetIndex(InputDimension-1, 0);
    row->SetRegions(sourceRegion);
    return row;

  }

  // Set the B-Spline coefficients using input images
  template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
  void 
  ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
  ::PadCoefficientImage(void)
  {
  
    // Define paste, extract and combine filters 
    typedef itk::PasteImageFilter<CoefficientImageType, CoefficientImageType, CoefficientImageType> PasteImageFilterType;
    typedef clitk::ExtractImageFilter<CoefficientImageType, CoefficientImageType> ExtractImageFilterType;
    typedef clitk::LinearCombinationImageFilter<CoefficientImageType, CoefficientImageType> LinearCombinationFilterType;

    // Regions
    typename   CoefficientImageType::RegionType sourceRegion=m_PaddedCoefficientImage->GetLargestPossibleRegion();
    typename   CoefficientImageType::RegionType destinationRegion=m_PaddedCoefficientImage->GetLargestPossibleRegion();
    typename   CoefficientImageType::RegionType::SizeType sourceSize=sourceRegion.GetSize();
    typename   CoefficientImageType::IndexType sourceIndex=sourceRegion.GetIndex();
    typename   CoefficientImageType::IndexType destinationIndex=destinationRegion.GetIndex();
 
    // JV Padding depends on the shape 
    switch (m_TransformShape)
      {
	/*  The shapes are
	    0: egg     4 CP  3 DOF
	    1: egg     5 CP  4 DOF 
	    2: rabbit  4 CP  3 DOF 
	    3: rabbit  5 CP  4 DOF
	    4: sputnik 4 CP  4 DOF
	    5: sputnik 5 CP  5 DOF
	    6: diamond 6 CP  5 DOF
	    7: diamond 7 CP  6 DOF
	*/
	
      case 0:
	{
	  // ----------------------------------------------------------------------
	  // The egg with 4 internal CP (starting from inhale)
	  // Periodic, constrained to zero at the reference 
	  // at position 3 and
	  // Coeff      row  BC1 BC2  0  BC3  1   2  BC4
	  // PaddedCoeff  R:  0   1   2   3   4   5   6 
	  // BC1= R4
	  // BC2= R5
	  // BC3= -weights[2]/weights[0] ( R2+R4 ) 
	  // BC4= R2
	  // ---------------------------------------------------------------------
	
	  //---------------------------------
	  // 1. First two temporal rows are identical: paste 1-2 to 0-1
	  typename PasteImageFilterType::Pointer paster0=PasteImageFilterType::New();
	  paster0->SetDestinationImage(m_PaddedCoefficientImage);
	  paster0->SetSourceImage(m_CoefficientImage);
	  sourceRegion.SetIndex(NInputDimensions-1,1);
	  sourceRegion.SetSize(NInputDimensions-1,2);
	  destinationIndex[InputDimension-1]=0;
	  paster0->SetDestinationIndex(destinationIndex);
	  paster0->SetSourceRegion(sourceRegion);
	  
	  //---------------------------------
	  // 1. Next temporal row = paste 0 to 2 
	  typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
	  typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
	  paster1->SetDestinationImage(paster0->GetOutput());
	  paster1->SetSourceImage(row0);
	  destinationIndex[InputDimension-1]=2;
	  paster1->SetDestinationIndex(destinationIndex);
	  paster1->SetSourceRegion(row0->GetLargestPossibleRegion());
	  
	  //---------------------------------
	  // 2. Middle row at index 3=BC
	  // Extract row at index 1, 2 (of coeff image)
	  typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
	  
	  // Combine
	  typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
	  combine1->SetInput(0,row0);
	  combine1->SetInput(1,row1);
	  combine1->SetA(-m_WeightRatio[2][0]);
	  combine1->SetB(-m_WeightRatio[2][0]);
	  combine1->Update();
	  typename CoefficientImageType::Pointer bc3Row=combine1->GetOutput();

	  // Paste middleRow at index 3 (padded coeff)
	  typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
	  paster2->SetDestinationImage(paster1->GetOutput());
	  paster2->SetSourceImage(bc3Row);
	  destinationIndex[InputDimension-1]=3;
	  paster2->SetDestinationIndex(destinationIndex);
	  paster2->SetSourceRegion(bc3Row->GetLargestPossibleRegion());
	
	  //---------------------------------
	  // 3. Next two temporal rows identical: paste 1,2 to 4,5
	  typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
	  paster3->SetDestinationImage(paster2->GetOutput());
	  paster3->SetSourceImage(m_CoefficientImage);
	  sourceRegion.SetIndex(InputDimension-1,1);
	  sourceRegion.SetSize(InputDimension-1,2);
	  destinationIndex[InputDimension-1]=4;
	  paster3->SetDestinationIndex(destinationIndex);
	  paster3->SetSourceRegion(sourceRegion);
	
	  // ---------------------------------
	  // 4. Row at index 6=BC4= R2
	  // Paste BC3 row at index 5
	  typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
	  paster4->SetDestinationImage(paster3->GetOutput());
	  paster4->SetSourceImage(row0);
	  destinationIndex[InputDimension-1]=6;
	  paster4->SetDestinationIndex(destinationIndex);
	  paster4->SetSourceRegion(row0->GetLargestPossibleRegion());
	
	  // Update the chain!
	  paster4->Update();
	  m_PaddedCoefficientImage= paster4->GetOutput();

	  break;
	}
	
      case 1:
	{
	  // ----------------------------------------------------------------------
	  // The egg with 5 internal CP (starting from inhale)
	  // Periodic, constrained to zero at the reference 
	  // at position 3 and
	  // Coeff      row  BC1 BC2  0  BC3  1   2   3  BC4
	  // PaddedCoeff  R:  0   1   2   3   4   5   6   7
	  // BC1= R5
	  // BC2= R6
	  // BC3= -weights[3]/weights[1] ( R2+R5 ) - R4 
	  // BC4= R2
	  // ---------------------------------------------------------------------
	  //---------------------------------
	  // 1. First two temporal rows are identical: paste 2-3 to 0-1
	  typename PasteImageFilterType::Pointer paster0=PasteImageFilterType::New();
	  paster0->SetDestinationImage(m_PaddedCoefficientImage);
	  paster0->SetSourceImage(m_CoefficientImage);
	  sourceRegion.SetIndex(NInputDimensions-1,2);
	  sourceRegion.SetSize(NInputDimensions-1,2);
	  destinationIndex[InputDimension-1]=0;
	  paster0->SetDestinationIndex(destinationIndex);
	  paster0->SetSourceRegion(sourceRegion);
	  
	  //---------------------------------
	  // 1. Next temporal row = paste 0 to 2 
	  typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
	  typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
	  paster1->SetDestinationImage(paster0->GetOutput());
	  paster1->SetSourceImage(row0);
	  destinationIndex[InputDimension-1]=2;
	  paster1->SetDestinationIndex(destinationIndex);
	  paster1->SetSourceRegion(row0->GetLargestPossibleRegion());
	  
	  //---------------------------------
	  // 2. Middle row at index 3=BC
	  // Extract row at index 1, 2 (of coeff image)
	  typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
	  typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
	  
	  // Combine
	  typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
	  combine1->SetInput(0,row0);
	  combine1->SetInput(1,row2);
	  combine1->SetA(1);
	  combine1->SetB(1);
	  // combine1->Update();
	  typename LinearCombinationFilterType::Pointer combine2=LinearCombinationFilterType::New();
	  combine2->SetInput(0,row1);
	  combine2->SetInput(1,combine1->GetOutput());
	  combine2->SetA(-1.);
	  combine2->SetB(-m_WeightRatio[3][1]);
	  combine2->Update();
	  typename CoefficientImageType::Pointer bc3Row=combine2->GetOutput();

	  // Paste middleRow at index 3 (padded coeff)
	  typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
	  paster2->SetDestinationImage(paster1->GetOutput());
	  paster2->SetSourceImage(bc3Row);
	  destinationIndex[InputDimension-1]=3;
	  paster2->SetDestinationIndex(destinationIndex);
	  paster2->SetSourceRegion(bc3Row->GetLargestPossibleRegion());
	
	  //---------------------------------
	  // 3. Next three temporal rows identical: paste 1,2,3 to 4,5,6
	  typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
	  paster3->SetDestinationImage(paster2->GetOutput());
	  paster3->SetSourceImage(m_CoefficientImage);
	  sourceRegion.SetIndex(InputDimension-1,1);
	  sourceRegion.SetSize(InputDimension-1,3);
	  destinationIndex[InputDimension-1]=4;
	  paster3->SetDestinationIndex(destinationIndex);
	  paster3->SetSourceRegion(sourceRegion);
	
	  // ---------------------------------
	  // 4. Row at index 7=BC4= R2
	  // Paste BC3 row at index 5
	  typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
	  paster4->SetDestinationImage(paster3->GetOutput());
	  paster4->SetSourceImage(row0);
	  destinationIndex[InputDimension-1]=7;
	  paster4->SetDestinationIndex(destinationIndex);
	  paster4->SetSourceRegion(row0->GetLargestPossibleRegion());
	
	  // Update the chain!
	  paster4->Update();
	  m_PaddedCoefficientImage= paster4->GetOutput();

	  break;
	}


// 	  // ----------------------------------------------------------------------
// 	  // The egg with 5 internal CP: 
// 	  // Periodic, C2 smooth everywhere and constrained to zero at the reference
// 	  // Coeff       row R5 BC1  0   1   2   3  BC2  R2  
// 	  // PaddedCoeff R:  0   1   2   3   4   5   6   7    
// 	  // BC1= -weights[2]/weights[0] ( R2+R5)
// 	  // BC2= BC1  
// 	  // ---------------------------------------------------------------------

// 	  // Extract rows with index 0 and 3
// 	  typename    CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
// 	  typename    CoefficientImageType::Pointer row3=this->ExtractTemporalRow(m_CoefficientImage, 3);
	  
// 	  // Paste the first row
// 	  typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
// 	  destinationIndex.Fill(0);
// 	  paster1->SetDestinationIndex(destinationIndex);
// 	  paster1->SetSourceRegion(row3->GetLargestPossibleRegion());
// 	  paster1->SetSourceImage(row3);
// 	  paster1->SetDestinationImage(m_PaddedCoefficientImage);

// 	  // Linearly Combine rows for BC1 and BC2
// 	  typedef clitk::LinearCombinationImageFilter<CoefficientImageType, CoefficientImageType> LinearCombinationFilterType;
// 	  typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
// 	  combine1->SetFirstInput(row0);
// 	  combine1->SetSecondInput(row3);
// 	  combine1->SetA(-m_WeightRatio[2][0]);
// 	  combine1->SetB(-m_WeightRatio[2][0]);
// 	  combine1->Update();
// 	  typename CoefficientImageType::Pointer bcRow=combine1->GetOutput();

// 	  // Paste the second row
// 	  typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
// 	  destinationIndex[InputDimension-1]=1;
// 	  paster2->SetDestinationIndex(destinationIndex);
// 	  paster2->SetSourceRegion(bcRow->GetLargestPossibleRegion());
// 	  paster2->SetSourceImage(bcRow);
// 	  paster2->SetDestinationImage(paster1->GetOutput());

// 	  // Paste the coefficientImage
// 	  typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
// 	  destinationIndex[InputDimension-1]=2;
// 	  paster3->SetDestinationIndex(destinationIndex);
// 	  paster3->SetSourceRegion(m_CoefficientImage->GetLargestPossibleRegion());
// 	  paster3->SetSourceImage(m_CoefficientImage);
// 	  paster3->SetDestinationImage(paster2->GetOutput());

// 	  // Paste the last two rows
// 	  typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
// 	  destinationIndex[InputDimension-1]=6;
// 	  paster4->SetDestinationIndex(destinationIndex);
// 	  paster4->SetSourceRegion(bcRow->GetLargestPossibleRegion());
// 	  paster4->SetSourceImage(bcRow);
// 	  paster4->SetDestinationImage(paster3->GetOutput());

// 	  typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
// 	  destinationIndex[InputDimension-1]=7;
// 	  paster5->SetDestinationIndex(destinationIndex);
// 	  paster5->SetSourceRegion(row0->GetLargestPossibleRegion());
// 	  paster5->SetSourceImage(row0);
// 	  paster5->SetDestinationImage(paster4->GetOutput());
	  
// 	  // Update the chain!
// 	  paster5->Update();
// 	  m_PaddedCoefficientImage= paster5->GetOutput();
	  
// 	  break;
// 	}

      case 2:
	{
	  // ----------------------------------------------------------------------
	  // The rabbit with 4 internal CP
	  // Periodic, constrained to zero at the reference 
	  // at position 3 and the extremes fixed with anit-symm bc
	  // Coeff      row  BC1  0   1  BC2  2  BC3 BC4 
	  // PaddedCoeff  R:  0   1   2   3   4   5   6 
	  // BC1= 2*R1-R2
	  // BC2= -weights[2]/weights[0] ( R2+R4 )
	  // BC3=  R1
	  // BC4= 2*R1-R4
	  // ---------------------------------------------------------------------

	  // ---------------------------------
	  // 0. First Row =BC1
	  typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
	  typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
	  typename LinearCombinationFilterType::Pointer combine0=LinearCombinationFilterType::New();
	  combine0->SetInput(0,row0);
	  combine0->SetInput(1,row1);
	  combine0->SetA(2.);
	  combine0->SetB(-1.);
	  combine0->Update();
	  typename CoefficientImageType::Pointer bc1Row=combine0->GetOutput();

	  typename PasteImageFilterType::Pointer paster0=PasteImageFilterType::New();
	  paster0->SetDestinationImage(m_PaddedCoefficientImage);
	  paster0->SetSourceImage(bc1Row);
	  destinationIndex[InputDimension-1]=0;
	  paster0->SetDestinationIndex(destinationIndex);
	  paster0->SetSourceRegion(bc1Row->GetLargestPossibleRegion());

	  //---------------------------------
	  // 1. Next two temporal rows are identical: paste 0-1 to 1-2
	  typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
	  paster1->SetDestinationImage(paster0->GetOutput());
	  paster1->SetSourceImage(m_CoefficientImage);
	  sourceRegion.SetIndex(NInputDimensions-1,0);
	  destinationIndex[InputDimension-1]=1;
	  sourceRegion.SetSize(NInputDimensions-1,2);
	  paster1->SetDestinationIndex(destinationIndex);
	  paster1->SetSourceRegion(sourceRegion);

	  //---------------------------------
	  // 2. Middle row at index 3=BC
	  // Extract row at index 1, 2 (of coeff image)
	  typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
  
	  // Combine
	  typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
	  combine1->SetInput(0,row1);
	  combine1->SetInput(1,row2);
	  combine1->SetA(-m_WeightRatio[2][0]);
	  combine1->SetB(-m_WeightRatio[2][0]);
	  combine1->Update();
	  typename CoefficientImageType::Pointer bc2Row=combine1->GetOutput();

	  // Paste middleRow at index 3 (padded coeff)
	  typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
	  paster2->SetDestinationImage(paster1->GetOutput());
	  paster2->SetSourceImage(bc2Row);
	  destinationIndex[InputDimension-1]=3;
	  paster2->SetDestinationIndex(destinationIndex);
	  paster2->SetSourceRegion(bc2Row->GetLargestPossibleRegion());
	  
	  //---------------------------------
	  // 3. Next temporal row is identical: paste 2 to 4
	  typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
	  paster3->SetDestinationImage(paster2->GetOutput());
	  paster3->SetSourceImage(row2);
	  destinationIndex[InputDimension-1]=4;
	  paster3->SetDestinationIndex(destinationIndex);
	  paster3->SetSourceRegion(row2->GetLargestPossibleRegion());

	  // ---------------------------------
	  // 4. Row at index 5=BC (paddedcoeff image) R1
	  // Paste BC3 row at index 5
	  typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
	  paster4->SetDestinationImage(paster3->GetOutput());
	  paster4->SetSourceImage(row0);
	  destinationIndex[InputDimension-1]=5;
	  paster4->SetDestinationIndex(destinationIndex);
	  paster4->SetSourceRegion(row0->GetLargestPossibleRegion());
    
	  // ---------------------------------
	  // 5. Paste BC4 row at index 6
	  typename LinearCombinationFilterType::Pointer combine3=LinearCombinationFilterType::New();
	  combine3->SetInput(0,row0);
	  combine3->SetInput(1,row2);
	  combine3->SetA(2.);
	  combine3->SetB(-1.);
	  combine3->Update();
	  typename CoefficientImageType::Pointer bc4Row=combine3->GetOutput();
	  typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
	  paster5->SetDestinationImage(paster4->GetOutput());
	  paster5->SetSourceImage(bc4Row);
	  destinationIndex[InputDimension-1]=6;
	  paster5->SetDestinationIndex(destinationIndex);
	  paster5->SetSourceRegion(bc4Row->GetLargestPossibleRegion());
	  
	  // Update the chain!
	  paster5->Update();
	  m_PaddedCoefficientImage= paster5->GetOutput();

	  break;
	}
     
      case 3: 
	{
	  // ----------------------------------------------------------------------
	  // The rabbit with 5 internal CP
	  // Periodic, constrained to zero at the reference at position 3.5 
	  // and the extremes fixed with anti-symmetrical BC
	  // Coeff      row  BC1  0   1  BC2  2   3  BC3 BC4
	  // PaddedCoeff  R:  0   1   2   3   4   5   6   7    
	  // BC1= 2*R1-R2
	  // BC2= -weights[3]/weights[1] ( R2+R5 ) - R4
	  // BC3= R1
	  // BC4= 2*R1-R5
	  // ---------------------------------------------------------------------

	  // ---------------------------------
	  // 0. First Row =BC1
	  typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
	  typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
	  typename LinearCombinationFilterType::Pointer combine0=LinearCombinationFilterType::New();
	  combine0->SetInput(0,row0);
	  combine0->SetInput(1,row1);
	  combine0->SetA(2.);
	  combine0->SetB(-1.);
	  combine0->Update();
	  typename CoefficientImageType::Pointer bc1Row=combine0->GetOutput();

	  typename PasteImageFilterType::Pointer paster0=PasteImageFilterType::New();
	  paster0->SetDestinationImage(m_PaddedCoefficientImage);
	  paster0->SetSourceImage(bc1Row);
	  destinationIndex[InputDimension-1]=0;
	  paster0->SetDestinationIndex(destinationIndex);
	  paster0->SetSourceRegion(bc1Row->GetLargestPossibleRegion());

	  //---------------------------------
	  // 1. Next two temporal rows are identical: paste 0-1 to 1-2
	  typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
	  paster1->SetDestinationImage(paster0->GetOutput());
	  paster1->SetSourceImage(m_CoefficientImage);
	  sourceRegion.SetIndex(InputDimension-1,0);
	  destinationIndex[InputDimension-1]=1;
	  sourceRegion.SetSize(NInputDimensions-1,2);
	  paster1->SetDestinationIndex(destinationIndex);
	  paster1->SetSourceRegion(sourceRegion);

	  //---------------------------------
	  // 2. Middle row at index 3=BC
	  // Extract row at index 1, 3 (of coeff image)
	  //typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
	  typename CoefficientImageType::Pointer row3=this->ExtractTemporalRow(m_CoefficientImage, 3);
  
	  // Combine
	  typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
	  combine1->SetInput(0,row1);
	  combine1->SetInput(1,row3);
	  combine1->SetA(1.);
	  combine1->SetB(1.);
	  combine1->Update();
	  
	  // Extract row at index 2 (coeff image)
	  typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);

	  // Combine
	  typename LinearCombinationFilterType::Pointer combine2=LinearCombinationFilterType::New();
	  combine2->SetInput(0,combine1->GetOutput());
	  combine2->SetInput(1,row2);
	  combine2->SetA(-m_WeightRatio[3][1]);
	  combine2->SetB(-1.);
	  combine2->Update();
	  typename CoefficientImageType::Pointer bc2Row=combine2->GetOutput();

	  // Paste middleRow at index 3 (padded coeff)
	  typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
	  paster2->SetDestinationImage(paster1->GetOutput());
	  paster2->SetSourceImage(bc2Row);
	  destinationIndex[InputDimension-1]=3;
	  paster2->SetDestinationIndex(destinationIndex);
	  paster2->SetSourceRegion(bc2Row->GetLargestPossibleRegion());
	  
	  //---------------------------------
	  // 3. Next two temporal rows are identical: paste 2,3 to 4,5
	  typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
	  paster3->SetDestinationImage(paster2->GetOutput());
	  paster3->SetSourceImage(m_CoefficientImage);
	  sourceRegion.SetIndex(InputDimension-1,2);
	  destinationIndex[InputDimension-1]=4;
	  sourceRegion.SetSize(NInputDimensions-1,2);
	  paster3->SetDestinationIndex(destinationIndex);
	  paster3->SetSourceRegion(sourceRegion);

	  // ---------------------------------
	  // 4. Row at index 6=BC (paddedcoeff image)R1
	  // Paste BC3 row at index 6
	  typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
	  paster4->SetDestinationImage(paster3->GetOutput());
	  paster4->SetSourceImage(row0);
	  destinationIndex[InputDimension-1]=6;
	  paster4->SetDestinationIndex(destinationIndex);
	  paster4->SetSourceRegion(row0->GetLargestPossibleRegion());

	  // ---------------------------------
	  // 5. Paste BC4 row at index 7
	  typename LinearCombinationFilterType::Pointer combine3=LinearCombinationFilterType::New();
	  combine3->SetInput(0,row0);
	  combine3->SetInput(1,row3);
	  combine3->SetA(2.);
	  combine3->SetB(-1.);
	  combine3->Update();
	  typename CoefficientImageType::Pointer bc4Row=combine3->GetOutput();
	  typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
	  paster5->SetDestinationImage(paster4->GetOutput());
	  paster5->SetSourceImage(bc4Row);
	  destinationIndex[InputDimension-1]=7;
	  paster5->SetDestinationIndex(destinationIndex);
	  paster5->SetSourceRegion(bc4Row->GetLargestPossibleRegion());
	  
	  // Update the chain!
	  paster5->Update();
	  m_PaddedCoefficientImage= paster5->GetOutput();

	  break;
	}
	
      case 4: 
	{
	  // ----------------------------------------------------------------------
	  // The sputnik with 4 internal CP
	  // Periodic, constrained to zero at the reference 
	  // at position 3 and one indepent extremes copied
	  // Coeff     row BC1  0   1  BC2  2  BC2  3
	  // PaddedCoeff R: 0   1   2   3   4   5   6      
	  // BC1= R6
	  // BC2= -weights[2]/weights[0] ( R2+R4 )
	  // BC3=  weights[2]/weights[0] ( R2-R4) + R1
	  // ---------------------------------------------------------------------

	  //---------------------------------
	  // 1. First Row is equal to last row: paste 3 row to 0
	  typename CoefficientImageType::Pointer row3=this->ExtractTemporalRow(m_CoefficientImage, 3);
	  typename PasteImageFilterType::Pointer paster0=PasteImageFilterType::New();
	  paster0->SetDestinationImage(m_PaddedCoefficientImage);
	  paster0->SetSourceImage(row3);
	  destinationIndex[InputDimension-1]=0;
	  paster0->SetDestinationIndex(destinationIndex);
	  paster0->SetSourceRegion(row3->GetLargestPossibleRegion());

	  //---------------------------------
	  // 1. Next two temporal rows are identical: paste 0-1 to 1-2
	  typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
	  paster1->SetDestinationImage(paster0->GetOutput());
	  paster1->SetSourceImage(m_CoefficientImage);
	  sourceRegion.SetIndex(NInputDimensions-1,0);
	  destinationIndex[InputDimension-1]=1;
	  sourceRegion.SetSize(NInputDimensions-1,2);
	  paster1->SetDestinationIndex(destinationIndex);
	  paster1->SetSourceRegion(sourceRegion);

	  //---------------------------------
	  // 2. Middle row at index 3=BC
	  // Extract row at index 1, 2 (of coeff image)
	  typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
	  typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
  
	  // Combine
	  typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
	  combine1->SetInput(0,row1);
	  combine1->SetInput(1,row2);
	  combine1->SetA(-m_WeightRatio[2][0]);
	  combine1->SetB(-m_WeightRatio[2][0]);
	  combine1->Update();
	  typename CoefficientImageType::Pointer bc1Row=combine1->GetOutput();

	  // Paste middleRow at index 3 (padded coeff)
	  typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
	  paster2->SetDestinationImage(paster1->GetOutput());
	  paster2->SetSourceImage(bc1Row);
	  destinationIndex[InputDimension-1]=3;
	  paster2->SetDestinationIndex(destinationIndex);
	  paster2->SetSourceRegion(bc1Row->GetLargestPossibleRegion());
	  
	  //---------------------------------
	  // 3. Next temporal row is identical: paste 2 to 4
	  typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
	  paster3->SetDestinationImage(paster2->GetOutput());
	  paster3->SetSourceImage(row2);
	  destinationIndex[InputDimension-1]=4;
	  paster3->SetDestinationIndex(destinationIndex);
	  paster3->SetSourceRegion(row2->GetLargestPossibleRegion());

	  //---------------------------------
	  // 4. Final row at index 5=BC3 (paddedcoeff image)R1 R2 R4 corresponds to index in coeff image 0 1 2
 	  // Extract row at index 1, 2 (of coeff image)already done
	  typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
    
	  // Combine
	  typename LinearCombinationFilterType::Pointer combine3=LinearCombinationFilterType::New();
	  combine3->SetInput(0,row1);
	  combine3->SetInput(1,row2);
	  combine3->SetA(1.);
	  combine3->SetB(-1.);

	  // Combine
	  typename LinearCombinationFilterType::Pointer combine4=LinearCombinationFilterType::New();
	  combine4->SetInput(0,combine3->GetOutput());
	  combine4->SetInput(1,row0);
	  combine4->SetA(m_WeightRatio[2][0]);
	  combine4->SetB(1.);
	  combine4->Update();
	  typename CoefficientImageType::Pointer bc2Row=combine4->GetOutput();

	  // Paste BC row at index 5
	  typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
	  paster4->SetDestinationImage(paster3->GetOutput());
	  paster4->SetSourceImage(bc2Row);
	  destinationIndex[InputDimension-1]=5;
	  paster4->SetDestinationIndex(destinationIndex);
	  paster4->SetSourceRegion(bc2Row->GetLargestPossibleRegion());
    
	  // Paste row 3  at index 6
	  typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
	  paster5->SetDestinationImage(paster4->GetOutput());
	  paster5->SetSourceImage(row3);
	  destinationIndex[InputDimension-1]=6;
	  paster5->SetDestinationIndex(destinationIndex);
	  paster5->SetSourceRegion(row3->GetLargestPossibleRegion());
	  
	  // Update the chain!
	  paster5->Update();
	  m_PaddedCoefficientImage= paster5->GetOutput();

	  break;
	}

      case 5: 
	{

	  // ----------------------------------------------------------------------
	  // The sputnik with 5 internal CP
	  // Periodic, constrained to zero at the reference 
	  // at position 3 and one indepent extreme
	  // Coeff     row BC1  0   1  BC2  2   3  BC3  4
	  // PaddedCoeff R: 0   1   2   3   4   5   6   7   
	  // BC1= R2 + R5 - R7
	  // BC2= -weights[3]/weights[1] ( R2+R5 ) - R4
	  // BC3= R1 + 0.5 R2 - 0.5 R7
	  // ----------------------------------------------------------------------
	  //---------------------------------
	  // 1. First Row =BC 
	  typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
	  typename CoefficientImageType::Pointer row3=this->ExtractTemporalRow(m_CoefficientImage, 3);
	  typename CoefficientImageType::Pointer row4=this->ExtractTemporalRow(m_CoefficientImage, 4);

	  // Combine
	  typename LinearCombinationFilterType::Pointer combine0=LinearCombinationFilterType::New();
	  combine0->SetInput(0,row1);
	  combine0->SetInput(1,row3);
	  combine0->SetA(1.);
	  combine0->SetB(1.);
	  //combine0->Update();
	  typename LinearCombinationFilterType::Pointer combine0bis=LinearCombinationFilterType::New();
	  combine0bis->SetInput(0,combine0->GetOutput());
	  combine0bis->SetInput(1,row4);
	  combine0bis->SetA(1.);
	  combine0bis->SetB(-1.);
	  combine0bis->Update();
	  typename CoefficientImageType::Pointer bc1Row=combine0bis->GetOutput();

	  typename PasteImageFilterType::Pointer paster0=PasteImageFilterType::New();
	  paster0->SetDestinationImage(m_PaddedCoefficientImage);
	  paster0->SetSourceImage(bc1Row);
	  destinationIndex[InputDimension-1]=0;
	  paster0->SetDestinationIndex(destinationIndex);
	  paster0->SetSourceRegion(bc1Row->GetLargestPossibleRegion());

	  //---------------------------------
	  // 1. Next two temporal rows are identical: paste 0-1 to 1-2
	  typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
	  paster1->SetDestinationImage(paster0->GetOutput());
	  paster1->SetSourceImage(m_CoefficientImage);
	  sourceRegion.SetIndex(NInputDimensions-1,0);
	  destinationIndex[InputDimension-1]=1;
	  sourceRegion.SetSize(NInputDimensions-1,2);
	  paster1->SetDestinationIndex(destinationIndex);
	  paster1->SetSourceRegion(sourceRegion);

	  //---------------------------------
	  // 2. Middle row at index 3=BC
	  // Extract row at index 1, 2 (of coeff image)
	  //typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
	  typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
  
	  // Combine
	  typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
	  combine1->SetInput(0,row1);
	  combine1->SetInput(1,row3);
	  combine1->SetA(1.);
	  combine1->SetB(1.);
	  combine1->Update();

	  typename LinearCombinationFilterType::Pointer combine2=LinearCombinationFilterType::New();
	  combine2->SetInput(0,combine1->GetOutput());
	  combine2->SetInput(1,row2);
	  combine2->SetA(-m_WeightRatio[3][1]);
	  combine2->SetB(-1.);
	  combine2->Update();
	  typename CoefficientImageType::Pointer bc2Row=combine2->GetOutput();

	  // Paste middleRow at index 3 (padded coeff)
	  typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
	  paster2->SetDestinationImage(paster1->GetOutput());
	  paster2->SetSourceImage(bc2Row);
	  destinationIndex[InputDimension-1]=3;
	  paster2->SetDestinationIndex(destinationIndex);
	  paster2->SetSourceRegion(bc2Row->GetLargestPossibleRegion());
	  
	  //---------------------------------
	  // 3. Next two temporal rows are identical: paste 2,3 to 4,5
	  typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
	  paster3->SetDestinationImage(paster2->GetOutput());
	  paster3->SetSourceImage(m_CoefficientImage);
	  sourceRegion.SetIndex(InputDimension-1,2);
	  destinationIndex[InputDimension-1]=4;
	  sourceRegion.SetSize(NInputDimensions-1,2);
	  paster3->SetDestinationIndex(destinationIndex);
	  paster3->SetSourceRegion(sourceRegion);

	  //---------------------------------
	  // 4. Final row at index 6=BC3
	  typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
    
	  // Combine
	  typename LinearCombinationFilterType::Pointer combine3=LinearCombinationFilterType::New();
	  combine3->SetInput(0,row1);
	  combine3->SetInput(1,row4);
	  combine3->SetA(1.);
	  combine3->SetB(-1.);
	  typename LinearCombinationFilterType::Pointer combine4=LinearCombinationFilterType::New();
	  combine4->SetInput(0,row0);
	  combine4->SetInput(1,combine3->GetOutput());
	  combine4->SetA(1.);
	  combine4->SetB(.5);
	  combine4->Update();
	  typename CoefficientImageType::Pointer bc3Row=combine4->GetOutput();

	  // Paste BC row at index 6
	  typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
	  paster4->SetDestinationImage(paster3->GetOutput());
	  paster4->SetSourceImage(bc3Row);
	  destinationIndex[InputDimension-1]=6;
	  paster4->SetDestinationIndex(destinationIndex);
	  paster4->SetSourceRegion(bc3Row->GetLargestPossibleRegion());
    
	  // Paste row 4  at index 7
	  typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
	  paster5->SetDestinationImage(paster4->GetOutput());
	  paster5->SetSourceImage(row4);
	  destinationIndex[InputDimension-1]=7;
	  paster5->SetDestinationIndex(destinationIndex);
	  paster5->SetSourceRegion(row4->GetLargestPossibleRegion());
	  
	  // Update the chain!
	  paster5->Update();
	  m_PaddedCoefficientImage= paster5->GetOutput();

	  break;
	}

      case 6: 
	{
	  // ----------------------------------------------------------------------
	  // The diamond with 4 internal CP:
	  // Periodic, constrained to zero at the reference at position 3
	  // Coeff      row 0   1   2  BC1  3  BC2  4
	  // PaddedCoeff  R:0   1   2   3   4   5   6    
	  // BC1= -weights[2]/weights[0] ( R2+R4 )
	  // BC2= weights[2]/weights[0]  ( R0+R2-R4-R6 ) + R1
	  // ---------------------------------------------------------------------

	  //---------------------------------
	  // 1. First Three temporal rows are identical: paste 0-2 to 0-2
	  typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
	  paster1->SetDestinationImage(m_PaddedCoefficientImage);
	  paster1->SetSourceImage(m_CoefficientImage);
	  sourceIndex.Fill(0);
	  destinationIndex.Fill(0);
	  sourceSize[NInputDimensions-1]=3;
	  sourceRegion.SetSize(sourceSize);
	  sourceRegion.SetIndex(sourceIndex);
	  paster1->SetDestinationIndex(destinationIndex);
	  paster1->SetSourceRegion(sourceRegion);

	  //---------------------------------
	  // 2. Middle row at index 3=BC
	  // Extract row at index 0, 4 (of coeff image)
	  typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
	  typename CoefficientImageType::Pointer row3=this->ExtractTemporalRow(m_CoefficientImage, 3);
  
	  // Combine
	  typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
	  combine1->SetInput(0,row2);
	  combine1->SetInput(1,row3);
	  combine1->SetA(-m_WeightRatio[2][0]);
	  combine1->SetB(-m_WeightRatio[2][0]);
	  combine1->Update();
	  typename CoefficientImageType::Pointer bc1Row=combine1->GetOutput();

	  // Paste middleRow at index 3 (padded coeff)
	  typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
	  paster2->SetDestinationImage(paster1->GetOutput());
	  paster2->SetSourceImage(bc1Row);
	  destinationIndex[InputDimension-1]=3;
	  paster2->SetDestinationIndex(destinationIndex);
	  paster2->SetSourceRegion(bc1Row->GetLargestPossibleRegion());
	  
	  //---------------------------------
	  // 3. Next  row identical: paste 3 to 4
	  typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
	  paster3->SetDestinationImage(paster2->GetOutput());
	  paster3->SetSourceImage(row3);
	  destinationIndex[InputDimension-1]=4;
	  paster3->SetDestinationIndex(destinationIndex);
	  paster3->SetSourceRegion(row3->GetLargestPossibleRegion());

	  //---------------------------------
	  // 4. Final row at index 6=BC (paddedcoeff image)R0 R2 R5 R7 R1 corresponds to index in coeff image 0 2 4 5 1
 	  // Extract row at index 0, 2 (of coeff image)already done
	  typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
	  typename CoefficientImageType::Pointer row4=this->ExtractTemporalRow(m_CoefficientImage, 4);
	  typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
    
	  // Combine
	  typename LinearCombinationFilterType::Pointer combine3=LinearCombinationFilterType::New();
	  combine3->SetInput(0,row0);
	  combine3->SetInput(1,row2);
	  combine3->SetA(1.);
	  combine3->SetB(1.);

	  // Combine
	  typename LinearCombinationFilterType::Pointer combine4=LinearCombinationFilterType::New();
	  combine4->SetInput(0,row3);
	  combine4->SetInput(1,row4);
	  combine4->SetA(1.);
	  combine4->SetB(1.);
	  
	  // Combine
	  typename LinearCombinationFilterType::Pointer combine5=LinearCombinationFilterType::New();
	  combine5->SetInput(0,combine3->GetOutput());
	  combine5->SetInput(1,combine4->GetOutput());
	  combine5->SetA(1.);
	  combine5->SetB(-1.);

	  // Combine
	  typename LinearCombinationFilterType::Pointer combine6=LinearCombinationFilterType::New();
	  combine6->SetInput(0,combine5->GetOutput());
	  combine6->SetInput(1,row1);
	  combine6->SetA(m_WeightRatio[2][0]);
	  combine6->SetB(1.);
	  combine6->Update();
	  typename CoefficientImageType::Pointer bc2Row=combine6->GetOutput();

	  // Paste BC row at index 5
	  typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
	  paster4->SetDestinationImage(paster3->GetOutput());
	  paster4->SetSourceImage(bc2Row);
	  destinationIndex[InputDimension-1]=5;
	  paster4->SetDestinationIndex(destinationIndex);
	  paster4->SetSourceRegion(bc2Row->GetLargestPossibleRegion());
    
	  // Paste last row at index 6
	  typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
	  paster5->SetDestinationImage(paster4->GetOutput());
	  paster5->SetSourceImage(row4);
	  destinationIndex[InputDimension-1]=6;
	  paster5->SetDestinationIndex(destinationIndex);
	  paster5->SetSourceRegion(row4->GetLargestPossibleRegion());
	  
	  // Update the chain!
	  paster5->Update();
	  m_PaddedCoefficientImage= paster5->GetOutput();

	  break;
	}
      case 7: 
	{
	  // ----------------------------------------------------------------------
	  // The diamond with 5 internal CP: 
	  // periodic, constrained to zero at the reference at position 3.5
	  // Coeff      row 0   1   2  BC1  3   4  BC2  5
	  // PaddedCoeff  R:0   1   2   3   4   5   6   7    
	  // BC1= -weights[3]/weights[1] ( R2+R5 ) - R4
	  // BC2= weights[2]/weights[0] ( R0+R2-R5-R7 ) + R1
	  // ---------------------------------------------------------------------

	  //---------------------------------
	  // 1. First Three temporal rows are identical: paste 0-2 to 0-2
	  typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
	  paster1->SetDestinationImage(m_PaddedCoefficientImage);
	  paster1->SetSourceImage(m_CoefficientImage);
	  sourceIndex.Fill(0);
	  destinationIndex.Fill(0);
	  sourceSize[NInputDimensions-1]=3;
	  sourceRegion.SetSize(sourceSize);
	  sourceRegion.SetIndex(sourceIndex);
	  paster1->SetDestinationIndex(destinationIndex);
	  paster1->SetSourceRegion(sourceRegion);

	  //---------------------------------
	  // 2. Middle row at index 3=BC
	  // Extract row at index 0, 4 (of coeff image)
	  typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
	  typename CoefficientImageType::Pointer row4=this->ExtractTemporalRow(m_CoefficientImage, 4);
  
	  // Combine
	  typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
	  combine1->SetInput(0,row2);
	  combine1->SetInput(1,row4);
	  combine1->SetA(1.);
	  combine1->SetB(1.);
	  combine1->Update();
	  
	  // Extract row at index 3 (coeff image)
	  typename CoefficientImageType::Pointer row3=this->ExtractTemporalRow(m_CoefficientImage, 3);

	  // Combine
	  typename LinearCombinationFilterType::Pointer combine2=LinearCombinationFilterType::New();
	  combine2->SetInput(0,combine1->GetOutput());
	  combine2->SetInput(1,row3);
	  combine2->SetA(-m_WeightRatio[3][1] );
	  combine2->SetB(-1.);
	  combine2->Update();
	  typename CoefficientImageType::Pointer bc1Row=combine2->GetOutput();

	  // Paste middleRow at index 3 (padded coeff)
	  typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
	  paster2->SetDestinationImage(paster1->GetOutput());
	  paster2->SetSourceImage(bc1Row);
	  destinationIndex[InputDimension-1]=3;
	  paster2->SetDestinationIndex(destinationIndex);
	  paster2->SetSourceRegion(bc1Row->GetLargestPossibleRegion());
	  
	  //---------------------------------
	  // 3. Next two temporal rows are identical: paste 3,4 to 4,5
	  typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
	  paster3->SetDestinationImage(paster2->GetOutput());
	  paster3->SetSourceImage(m_CoefficientImage);
	  sourceIndex[InputDimension-1]=3;
	  destinationIndex[InputDimension-1]=4;
	  sourceSize[NInputDimensions-1]=2;
	  sourceRegion.SetSize(sourceSize);
	  sourceRegion.SetIndex(sourceIndex);
	  paster3->SetDestinationIndex(destinationIndex);
	  paster3->SetSourceRegion(sourceRegion);

	  //---------------------------------
	  // 4. Final row at index 6=BC (paddedcoeff image)R0 R2 R5 R7 R1 corresponds to index in coeff image 0 2 4 5 1
 	  // Extract row at index 0, 2 (of coeff image)already done
	  typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
	  typename CoefficientImageType::Pointer row5=this->ExtractTemporalRow(m_CoefficientImage, 5);
	  typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
    
	  // Combine
	  typename LinearCombinationFilterType::Pointer combine3=LinearCombinationFilterType::New();
	  combine3->SetInput(0,row0);
	  combine3->SetInput(1,row2);
	  combine3->SetA(1.);
	  combine3->SetB(1.);

	  // Combine
	  typename LinearCombinationFilterType::Pointer combine4=LinearCombinationFilterType::New();
	  combine4->SetInput(0,row4);
	  combine4->SetInput(1,row5);
	  combine4->SetA(1.);
	  combine4->SetB(1.);
	  
	  // Combine
	  typename LinearCombinationFilterType::Pointer combine5=LinearCombinationFilterType::New();
	  combine5->SetInput(0,combine3->GetOutput());
	  combine5->SetInput(1,combine4->GetOutput());
	  combine5->SetA(1.);
	  combine5->SetB(-1.);

	  // Combine
	  typename LinearCombinationFilterType::Pointer combine6=LinearCombinationFilterType::New();
	  combine6->SetInput(0,combine5->GetOutput());
	  combine6->SetInput(1,row1);
	  combine6->SetA(m_WeightRatio[2][0]);
	  combine6->SetB(1.);
	  combine6->Update();
	  typename CoefficientImageType::Pointer bc2Row=combine6->GetOutput();

	  // Paste BC row at index 6
	  typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
	  paster4->SetDestinationImage(paster3->GetOutput());
	  paster4->SetSourceImage(bc2Row);
	  destinationIndex[InputDimension-1]=6;
	  paster4->SetDestinationIndex(destinationIndex);
	  paster4->SetSourceRegion(bc2Row->GetLargestPossibleRegion());
    
	  // Paste last row at index 7
	  typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
	  paster5->SetDestinationImage(paster4->GetOutput());
	  paster5->SetSourceImage(row5);
	  destinationIndex[InputDimension-1]=7;
	  paster5->SetDestinationIndex(destinationIndex);
	  paster5->SetSourceRegion(row5->GetLargestPossibleRegion());
	  
	  // Update the chain!
	  paster5->Update();
	  m_PaddedCoefficientImage= paster5->GetOutput();

	  break;
	}

      case 9:
	{
	  // ----------------------------------------------------------------------
	  // The sputnik with 5 internal CP T''(0)=T''(10)
	  // Periodic, constrained to zero at the reference 
	  // at position 3 and one indepent extreme
	  // Coeff     row BC1  0   1  BC2  2   3  BC3  4
	  // PaddedCoeff R: 0   1   2   3   4   5   6   7   
	  // BC1= -R2+R5+R7
	  // BC2= -weights[3]/weights[1] ( R2+R5 ) - R4
	  // BC3= R1
	  // ---------------------------------------------------------------------
	 
	  //---------------------------------
	  // 1. First Row =BC 
	  typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
	  typename CoefficientImageType::Pointer row3=this->ExtractTemporalRow(m_CoefficientImage, 3);
	  typename CoefficientImageType::Pointer row4=this->ExtractTemporalRow(m_CoefficientImage, 4);

	  // Combine
	  typename LinearCombinationFilterType::Pointer combine0=LinearCombinationFilterType::New();
	  combine0->SetInput(0,row3);
	  combine0->SetInput(1,row4);
	  combine0->SetA(1.);
	  combine0->SetB(1.);
	  typename LinearCombinationFilterType::Pointer combine0bis=LinearCombinationFilterType::New();
	  combine0bis->SetInput(0,combine0->GetOutput());
	  combine0bis->SetInput(1,row1);
	  combine0bis->SetA(1.);
	  combine0bis->SetB(-1.);
	  combine0bis->Update();
	  typename CoefficientImageType::Pointer bc1Row=combine0bis->GetOutput();

	  typename PasteImageFilterType::Pointer paster0=PasteImageFilterType::New();
	  paster0->SetDestinationImage(m_PaddedCoefficientImage);
	  paster0->SetSourceImage(bc1Row);
	  destinationIndex[InputDimension-1]=0;
	  paster0->SetDestinationIndex(destinationIndex);
	  paster0->SetSourceRegion(bc1Row->GetLargestPossibleRegion());

	  //---------------------------------
	  // 1. Next two temporal rows are identical: paste 0-1 to 1-2
	  typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
	  paster1->SetDestinationImage(paster0->GetOutput());
	  paster1->SetSourceImage(m_CoefficientImage);
	  sourceRegion.SetIndex(NInputDimensions-1,0);
	  destinationIndex[InputDimension-1]=1;
	  sourceRegion.SetSize(NInputDimensions-1,2);
	  paster1->SetDestinationIndex(destinationIndex);
	  paster1->SetSourceRegion(sourceRegion);

	  //---------------------------------
	  // 2. Middle row at index 3=BC
	  // Extract row at index 1, 2 (of coeff image)
	  // typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
	  typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
  
	  // Combine
	  typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
	  combine1->SetInput(0,row1);
	  combine1->SetInput(1,row3);
	  combine1->SetA(1.);
	  combine1->SetB(1.);
	  combine1->Update();

	  typename LinearCombinationFilterType::Pointer combine2=LinearCombinationFilterType::New();
	  combine2->SetInput(0,combine1->GetOutput());
	  combine2->SetInput(1,row2);
	  combine2->SetA(-m_WeightRatio[3][1]);
	  combine2->SetB(-1.);
	  combine2->Update();
	  typename CoefficientImageType::Pointer bc2Row=combine2->GetOutput();

	  // Paste middleRow at index 3 (padded coeff)
	  typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
	  paster2->SetDestinationImage(paster1->GetOutput());
	  paster2->SetSourceImage(bc2Row);
	  destinationIndex[InputDimension-1]=3;
	  paster2->SetDestinationIndex(destinationIndex);
	  paster2->SetSourceRegion(bc2Row->GetLargestPossibleRegion());
	  
	  //---------------------------------
	  // 3. Next two temporal rows are identical: paste 2,3 to 4,5
	  typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
	  paster3->SetDestinationImage(paster2->GetOutput());
	  paster3->SetSourceImage(m_CoefficientImage);
	  sourceRegion.SetIndex(InputDimension-1,2);
	  destinationIndex[InputDimension-1]=4;
	  sourceRegion.SetSize(NInputDimensions-1,2);
	  paster3->SetDestinationIndex(destinationIndex);
	  paster3->SetSourceRegion(sourceRegion);

	  //---------------------------------
	  // 4. Final row at index 6=BC3
	  typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
    
	  // 	  // Combine
	  // 	  typename LinearCombinationFilterType::Pointer combine3=LinearCombinationFilterType::New();
	  // 	  combine3->SetInput(0,row0);
	  // 	  combine3->SetInput(1,row1);
	  // 	  combine3->SetA(1.);
	  // 	  combine3->SetB(0.5);
	  // 	  combine3->Update();
	  // 	  typename CoefficientImageType::Pointer bc3Row=combine3->GetOutput();

	  // Paste BC row at index 6
	  typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
	  paster4->SetDestinationImage(paster3->GetOutput());
	  paster4->SetSourceImage(row0);
	  destinationIndex[InputDimension-1]=6;
	  paster4->SetDestinationIndex(destinationIndex);
	  paster4->SetSourceRegion(row0->GetLargestPossibleRegion());
    
	  // Paste row 4  at index 7
	  typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
	  paster5->SetDestinationImage(paster4->GetOutput());
	  paster5->SetSourceImage(row4);
	  destinationIndex[InputDimension-1]=7;
	  paster5->SetDestinationIndex(destinationIndex);
	  paster5->SetSourceRegion(row4->GetLargestPossibleRegion());
	  
	  // Update the chain!
	  paster5->Update();
	  m_PaddedCoefficientImage= paster5->GetOutput();

	  break;
	}

      default:
	DD ("Shape not available");
      }
    
  }
  
  
  // // Extract coefficients from padded version
  // template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
  // void 
  // ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
  // ::ExtractCoefficientImage( )
  // {
  //   ////DD("extract coeff image");
  //   typename ExtractImageFilterType::Pointer extractImageFilter = ExtractImageFilterType::New();
  //   typename CoefficientImageType::RegionType extractionRegion=m_PaddedCoefficientImage->GetLargestPossibleRegion();
  //   typename CoefficientImageType::RegionType::SizeType extractionSize=extractionRegion.GetSize();
  //   typename CoefficientImageType::RegionType::IndexType extractionIndex = extractionRegion.GetIndex();
  //   extractionSize[InputDimension-1]-=4;    
  //   extractionIndex[InputDimension-1]=2;
  //   extractionRegion.SetSize(extractionSize);
  //   extractionRegion.SetIndex(extractionIndex);
  //   extractImageFilter->SetInput(m_PaddedCoefficientImage);
  //   extractImageFilter->SetExtractionRegion(extractionRegion);
  //   extractImageFilter->Update();
  //   m_CoefficientImage=extractImageFilter->GetOutput();
  // } 


  // Print self
  template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
  void
  ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
  ::PrintSelf(std::ostream &os, itk::Indent indent) const
  {

    this->Superclass::PrintSelf(os, indent);

    os << indent << "GridRegion: " << m_GridRegion << std::endl;
    os << indent << "GridOrigin: " << m_GridOrigin << std::endl;
    os << indent << "GridSpacing: " << m_GridSpacing << std::endl;
    os << indent << "GridDirection: " << m_GridDirection << std::endl;
    os << indent << "IndexToPoint: " << m_IndexToPoint << std::endl;
    os << indent << "PointToIndex: " << m_PointToIndex << std::endl;

    os << indent << "CoefficientImage: [ ";
    os << m_CoefficientImage.GetPointer() << " ]" << std::endl;

    os << indent << "WrappedImage: [ ";
    os << m_WrappedImage.GetPointer() << " ]" << std::endl;
 
    os << indent << "InputParametersPointer: " 
       << m_InputParametersPointer << std::endl;
    os << indent << "ValidRegion: " << m_ValidRegion << std::endl;
    os << indent << "LastJacobianIndex: " << m_LastJacobianIndex << std::endl;
    os << indent << "BulkTransform: ";
    os << m_BulkTransform << std::endl;

    if ( m_BulkTransform )
      {
	os << indent << "BulkTransformType: " 
	   << m_BulkTransform->GetNameOfClass() << std::endl;
      }
    os << indent << "VectorBSplineInterpolator: ";
    os << m_VectorInterpolator.GetPointer() << std::endl;
    os << indent << "Mask: ";
    os << m_Mask<< std::endl;
  }


  // Verify 
  template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
  bool 
  ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
  ::InsideValidRegion( const ContinuousIndexType& index ) const
  {
    bool inside = true;

    if ( !m_ValidRegion.IsInside( index ) )
      {
	inside = false;
      }

    // JV verify all dimensions 
    if ( inside)
      {
	typedef typename ContinuousIndexType::ValueType ValueType;
	for( unsigned int j = 0; j < NInputDimensions; j++ )
	  {
	    if (m_SplineOrderOdd[j])
	      {
		if ( index[j] >= static_cast<ValueType>( m_ValidRegionLast[j] ) )
		  { 
		    inside = false;
		    break;
		  }
	      }
	  }
      }
    return inside;
  }


  // Transform a point
  template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
  typename ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
  ::OutputPointType
  ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
  ::TransformPoint(const InputPointType &inputPoint) const 
  {
  
    InputPointType transformedPoint;
    OutputPointType outputPoint;

    // BulkTransform
    if ( m_BulkTransform )
      {
	transformedPoint = m_BulkTransform->TransformPoint( inputPoint );
      }
    else
      {
	transformedPoint = inputPoint;
      }
    
    // Deformable transform
    if ( m_PaddedCoefficientImage )
      {
	// Check if inside mask
	if(m_Mask &&  !(m_Mask->IsInside(inputPoint) ) )
	  {
	    // Outside: no (deformable) displacement
	    return transformedPoint;
	  }		
	
	// Check if inside valid region
	bool inside = true;
	ContinuousIndexType index;
	this->TransformPointToContinuousIndex( inputPoint, index );
	inside = this->InsideValidRegion( index );
	if ( !inside )
	  {
	    // Outside: no (deformable) displacement
	    DD("outside valid region");
	    return transformedPoint;
	  }
		
	// Call the vector interpolator
	itk::Vector<TCoordRep,SpaceDimension> displacement=m_VectorInterpolator->EvaluateAtContinuousIndex(index);
	
	// JV add for the spatial dimensions
	outputPoint=transformedPoint;
	for (unsigned int i=0; i<NInputDimensions-1; i++)
	  outputPoint[i] += displacement[i];

      }

    else
      {
	itkWarningMacro( << "B-spline coefficients have not been set" );
	outputPoint = transformedPoint;
      }
    
    return outputPoint;
  }



  //JV Deformably transform a point
  template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
  typename ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
  ::OutputPointType
  ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
  ::DeformablyTransformPoint(const InputPointType &inputPoint) const 
  {
    OutputPointType outputPoint;
    if ( m_PaddedCoefficientImage )
      {

	// Check if inside mask
	if(m_Mask && !(m_Mask->IsInside(inputPoint) ) )
	  {
	    // Outside: no (deformable) displacement
	    return inputPoint;
	  }	
	
	// Check if inside
	bool inside = true;
	ContinuousIndexType index;
	this->TransformPointToContinuousIndex( inputPoint, index );
	inside = this->InsideValidRegion( index );

	if ( !inside )
	  {
	    //outside: no (deformable) displacement
	    outputPoint = inputPoint;
	    return outputPoint;
	  }

	// Call the vector interpolator
	itk::Vector<TCoordRep,SpaceDimension> displacement=m_VectorInterpolator->EvaluateAtContinuousIndex(index);
	
	// JV add for the spatial dimensions
	outputPoint=inputPoint;
	for (unsigned int i=0; i<NInputDimensions-1; i++)
	  outputPoint[i] += displacement[i];
      }

    // No coefficients available
    else
      {
	itkWarningMacro( << "B-spline coefficients have not been set" );
	outputPoint = inputPoint;
      }
   
    return outputPoint;
  }
  

  // JV weights are identical as for transformpoint, could be done simultaneously in metric!!!!
  // Compute the Jacobian in one position 
  template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
  void
  ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
  ::ComputeJacobianWithRespectToParameters( const InputPointType & point, JacobianType & jacobian) const
  {
  
    //========================================================
    // Zero all components of jacobian
    //========================================================
    // JV not thread safe (m_LastJacobianIndex), instantiate N transforms
    // NOTE: for efficiency, we only need to zero out the coefficients
    // that got fill last time this method was called.
    unsigned int j=0,b=0;
      
    //Set the previously-set to zero 
    for ( j = 0; j < SpaceDimension; j++ )
      {
	m_FirstIterator[j].GoToBegin();
	while ( !m_FirstIterator[j].IsAtEnd() )
	  {
	    m_FirstIterator[j].Set( m_ZeroVector );
	    ++(m_FirstIterator[j]);
	  }
      }
    
    //Set the previously-set to zero 
    for ( j = 0; j < SpaceDimension; j++ )
      {
	m_SecondIterator[j].GoToBegin();
	while ( ! (m_SecondIterator[j]).IsAtEnd() )
	  {
	    m_SecondIterator[j].Set( m_ZeroVector );
	    ++(m_SecondIterator[j]);
	  }
      }

    //Set the previously-set to zero 
    if (m_ThirdSize)
      for ( j = 0; j < SpaceDimension; j++ )
	{
	  m_ThirdIterator[j].GoToBegin();
	  while ( ! (m_ThirdIterator[j]).IsAtEnd() )
	    {
	      m_ThirdIterator[j].Set( m_ZeroVector );
	      ++(m_ThirdIterator[j]);
	    }
	}

    //Set the previously-set to zero 
    if (m_BCSize)
      for (b=0; b<m_BCSize;b++)
	if ( !( m_FirstRegion.IsInside(m_BCRegions[b]) | m_SecondRegion.IsInside(m_BCRegions[b]) )  )
	  for ( j = 0; j < SpaceDimension; j++ ) 
	    {
	      m_BCIterators[j][b].GoToBegin();
	      while ( ! (m_BCIterators[j][b]).IsAtEnd() )
		{
		  m_BCIterators[j][b].Set( m_ZeroVector );
		  ++(m_BCIterators[j][b]);
		}
	    }
	
    //Set the previously-set to zero 
    if (m_BC2Size)
      for (b=0; b<m_BC2Size;b++)
	if ( !( m_FirstRegion.IsInside(m_BC2Regions[b]) | m_SecondRegion.IsInside(m_BC2Regions[b]) )  )
	  for ( j = 0; j < SpaceDimension; j++ ) 
	    {
	      m_BC2Iterators[j][b].GoToBegin();
	      while ( ! (m_BC2Iterators[j][b]).IsAtEnd() )
		{
		  m_BC2Iterators[j][b].Set( m_ZeroVector );
		  ++(m_BC2Iterators[j][b]);
		}
	    }

    //Set the previously-set to zero 
    if (m_BC3Size)
      for (b=0; b<m_BC3Size;b++)
	if ( !( m_FirstRegion.IsInside(m_BC3Regions[b]) | m_SecondRegion.IsInside(m_BC3Regions[b]) )  )
	  for ( j = 0; j < SpaceDimension; j++ ) 
	    {
	      m_BC3Iterators[j][b].GoToBegin();
	      while ( ! (m_BC3Iterators[j][b]).IsAtEnd() )
		{
		  m_BC3Iterators[j][b].Set( m_ZeroVector );
		  ++(m_BC3Iterators[j][b]);
		}
	    }


    //========================================================
    // For each dimension, copy the weight to the support region
    //========================================================

    // Check if inside mask
    if(m_Mask &&  !(m_Mask->IsInside(point) ) )
      {
	// Outside: no (deformable) displacement
        jacobian = m_SharedDataBSplineJacobian;
        return;
      }	

    // Get index   
    this->TransformPointToContinuousIndex( point, m_Index );

    // NOTE: if the support region does not lie totally within the grid
    // we assume zero displacement and return the input point
    if ( !this->InsideValidRegion( m_Index ) )
      {
        jacobian = m_SharedDataBSplineJacobian;
        return;
      }

    // Compute interpolation weights
    const WeightsDataType *weights=NULL;
    m_VectorInterpolator->EvaluateWeightsAtContinuousIndex( m_Index, &weights, m_LastJacobianIndex);

    // Get support
    m_SupportRegion.SetIndex( m_LastJacobianIndex );
    WrapRegion(m_SupportRegion, m_FirstRegion, m_SecondRegion, m_ThirdRegion, m_BCRegions, m_BCValues, m_BC2Regions, m_BC2Values, m_BC3Regions, m_BC3Values, m_InitialOffset);
    m_ThirdSize=m_ThirdRegion.GetSize()[InputDimension-1]; 
    m_BCSize=m_BCRegions.size();
    m_BC2Size=m_BC2Regions.size();
    m_BC3Size=m_BC3Regions.size();

    // Reset the iterators
    for ( j = 0; j < SpaceDimension ; j++ ) 
      {
	m_FirstIterator[j] = IteratorType( m_JacobianImage[j], m_FirstRegion);
	m_SecondIterator[j] = IteratorType( m_JacobianImage[j], m_SecondRegion);
	if(m_ThirdSize)	m_ThirdIterator[j] = IteratorType( m_JacobianImage[j], m_ThirdRegion);

	m_BCIterators[j].resize(m_BCSize);
	for (b=0; b<m_BCSize;b++)
	  m_BCIterators[j][b]= IteratorType( m_JacobianImage[j], m_BCRegions[b]);
	m_BC2Iterators[j].resize(m_BC2Size);
	for (b=0; b<m_BC2Size;b++)
	  m_BC2Iterators[j][b]= IteratorType( m_JacobianImage[j], m_BC2Regions[b]);
	m_BC3Iterators[j].resize(m_BC3Size);
	for (b=0; b<m_BC3Size;b++)
	  m_BC3Iterators[j][b]= IteratorType( m_JacobianImage[j], m_BC3Regions[b]);
      }

    // Skip if on a fixed condition
    if(m_InitialOffset)
      {
	if (m_BCSize)  weights+=m_InitialOffset*m_BCRegions[0].GetNumberOfPixels();
	else std::cerr<<"InitialOffset without BCSize: Error!!!!!!!"<<std::endl;
      }

    //copy weight to jacobian image
    for ( j = 0; j < SpaceDimension; j++ )
      {
	// For each dimension, copy the weight to the support region
	while ( ! (m_FirstIterator[j]).IsAtEnd() )
	  {
	    m_ZeroVector[j]=*weights;
	    (m_FirstIterator[j]).Set( m_ZeroVector);
	    ++(m_FirstIterator[j]);
	    weights++;
	  }
	
	m_ZeroVector[j]=itk::NumericTraits<JacobianValueType>::Zero;
	if (j != SpaceDimension-1)  weights-=m_FirstRegion.GetNumberOfPixels();
      }
            
    // Skip BC1 and go to the second region
    if (m_BCSize) weights+=m_BCRegions[0].GetNumberOfPixels();
      
    // For each dimension, copy the weight to the support region
    //copy weight to jacobian image
    for ( j = 0; j < SpaceDimension; j++ )
      {
	while ( ! (m_SecondIterator[j]).IsAtEnd() )
	  {
	    m_ZeroVector[j]=*weights;
	    (m_SecondIterator[j]).Set( m_ZeroVector);
	    ++(m_SecondIterator[j]);
	    weights++;
	  }
	weights-=m_SecondRegion.GetNumberOfPixels();
	m_ZeroVector[j]=itk::NumericTraits<JacobianValueType>::Zero;
      }

    // Now do BC1:
    if (m_BCSize)
      {	
	//  Put pointer in correct position
	weights-=m_BCRegions[0].GetNumberOfPixels();
	
	for ( j = 0; j < SpaceDimension; j++ )
	  {
	    for ( b=0; b < m_BCSize; b++ )
	      {
		while ( ! (m_BCIterators[j][b]).IsAtEnd() )
		  {
		    //copy weight to jacobian image
		    m_ZeroVector[j]=(*weights) * m_BCValues[b];
		    (m_BCIterators[j][b]).Value()+= m_ZeroVector;
		    ++(m_BCIterators[j][b]);
		    weights++;
		  }
		weights-=m_BCRegions[b].GetNumberOfPixels();
	      }
	    m_ZeroVector[j]=itk::NumericTraits<JacobianValueType>::Zero;
	  }
	weights+=m_BCRegions[m_BCSize-1].GetNumberOfPixels();
      }
        
    // Add the BC2 to the weights
    if (m_BC2Size)
      {
	// Move further in the weights pointer
	weights+=m_SecondRegion.GetNumberOfPixels();

	for ( j = 0; j < SpaceDimension; j++ )
	  {
	    for ( b=0; b < m_BC2Size; b++ )
	      {
		while ( ! (m_BC2Iterators[j][b]).IsAtEnd() )
		  {
		    //copy weight to jacobian image
		    m_ZeroVector[j]=(*weights) * m_BC2Values[b];
		    (m_BC2Iterators[j][b]).Value()+= m_ZeroVector;
		    ++(m_BC2Iterators[j][b]);
		    weights++;
		  }
		weights-=m_BC2Regions[b].GetNumberOfPixels();
	      }
	    m_ZeroVector[j]=itk::NumericTraits<JacobianValueType>::Zero;
	  }
	// Move further in the weights pointer
	weights+=m_BC2Regions[m_BC2Size-1].GetNumberOfPixels();
      }
  
    // Third Region
    if(m_ThirdSize)
      {
	// For each dimension, copy the weight to the support region
	//copy weight to jacobian image
	for ( j = 0; j < SpaceDimension; j++ )
	  {
	    while ( ! (m_ThirdIterator[j]).IsAtEnd() )
	      {
		m_ZeroVector[j]=*weights;
		(m_ThirdIterator[j]).Value()+= m_ZeroVector;
		++(m_ThirdIterator[j]);
		weights++;
	      }

	    // Move further in the weights pointer?
	    if (j != SpaceDimension-1) weights-=m_ThirdRegion.GetNumberOfPixels();
	    m_ZeroVector[j]=itk::NumericTraits<JacobianValueType>::Zero;
	  }
      }

    // Add the BC3 to the weights
    if (m_BC3Size)
      {
	for ( j = 0; j < SpaceDimension; j++ )
	  {
	    for ( b=0; b < m_BC3Size; b++ )
	      {
		while ( ! (m_BC3Iterators[j][b]).IsAtEnd() )
		  {
		    //copy weight to jacobian image
		    m_ZeroVector[j]=(*weights) * m_BC3Values[b];
		    (m_BC3Iterators[j][b]).Value()+= m_ZeroVector;
		    ++(m_BC3Iterators[j][b]);
		    weights++;
		  }
		weights-=m_BC3Regions[b].GetNumberOfPixels();
	      }
	    m_ZeroVector[j]=itk::NumericTraits<JacobianValueType>::Zero;
	  }
      }

    // Return the result
    jacobian = m_SharedDataBSplineJacobian;
  }

 
  template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
  inline void 
  ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
  ::WrapRegion( const RegionType & m_SupportRegion, 
		RegionType & m_FirstRegion, 
		RegionType & m_SecondRegion,
		RegionType & m_ThirdRegion,  
		std::vector<RegionType>& m_BCRegions,std::vector<double>& m_BCValues, 
		std::vector<RegionType>& m_BC2Regions,std::vector<double>& m_BC2Values,
		std::vector<RegionType>& m_BC3Regions,std::vector<double>& m_BC3Values,
		unsigned int& m_InitialOffset ) const
  {

    // Synchronize regions
    m_InitialOffset=0;
    m_FirstRegion=m_SupportRegion;
    m_BCRegion=m_SupportRegion;
    m_BCRegion.SetSize(InputDimension-1,1);
    m_SecondRegion=m_SupportRegion;
    m_ThirdRegion=m_SupportRegion;
    m_ThirdRegion.SetSize(InputDimension-1,0);
    m_BC3Regions.resize(0);


    // BC depends on shape
    switch(m_TransformShape)
      {
	/*  The shapes are
	    0: egg     4 CP  3 DOF
	    1: egg     5 CP  4 DOF 
	    2: rabbit  4 CP  3 DOF 
	    3: rabbit  5 CP  4 DOF
	    4: sputnik 4 CP  4 DOF
	    5: sputnik 5 CP  5 DOF
	    6: diamond 6 CP  5 DOF
	    7: diamond 7 CP  6 DOF
	*/
  
      case 0:
	{	  
	  // ----------------------------------------------------------------------
	  // The egg with 4 internal CP (starting from inhale)
	  // Periodic, constrained to zero at the reference 
	  // at position 3 and
	  // Coeff      row  BC1 BC2  0  BC3  1   2  BC4
	  // PaddedCoeff  R:  0   1   2   3   4   5   6 
	  // BC1= R4
	  // BC2= R5
	  // BC3= -weights[2]/weights[0] ( R2+R4 ) 
	  // BC4= R2
	  // ---------------------------------------------------------------------
	  switch(m_SupportRegion.GetIndex(InputDimension-1)) 
	    {
	    case 0:
	      {
		// Lower end 
		m_FirstRegion.SetSize(InputDimension-1,2);
		m_FirstRegion.SetIndex(InputDimension-1,1);
		
		// BC
		m_BCRegions.resize(0);
		     
		// Second part
		m_SecondRegion.SetSize(InputDimension-1,1);
		m_SecondRegion.SetIndex(InputDimension-1,0);
 
		// BC
		m_BC2Regions.resize(2);
		m_BC2Values.resize(2);
		m_BCRegion.SetIndex(InputDimension-1,0);
		m_BC2Regions[0]=m_BCRegion;
		m_BC2Values[0]=-m_WeightRatio[2][0];
		m_BCRegion.SetIndex(InputDimension-1,1);
		m_BC2Regions[1]=m_BCRegion;
		m_BC2Values[1]= -m_WeightRatio[2][0];
		    		      
		break;
	      }
	    case 1:
	      {
		// Lower end 
		m_FirstRegion.SetSize(InputDimension-1,1);
		m_FirstRegion.SetIndex(InputDimension-1,2);
		      
		// BC
		m_BCRegions.resize(0);
		      
		// Second part
		m_SecondRegion.SetSize(InputDimension-1,1);
		m_SecondRegion.SetIndex(InputDimension-1,0);
		      
		// BC
		m_BC2Regions.resize(2);
		m_BC2Values.resize(2);
		m_BCRegion.SetIndex(InputDimension-1,0);
		m_BC2Regions[0]=m_BCRegion;
		m_BC2Values[0]=-m_WeightRatio[2][0];
		m_BCRegion.SetIndex(InputDimension-1,1);
		m_BC2Regions[1]=m_BCRegion;
		m_BC2Values[1]= -m_WeightRatio[2][0];

		// Third Part
		m_FirstRegion.SetSize(InputDimension-1,1);
		m_FirstRegion.SetIndex(InputDimension-1,1);
		      
		break;
	      }
	    case 2:
	      {
		// Lower end 
		m_FirstRegion.SetSize(InputDimension-1,1);
		m_FirstRegion.SetIndex(InputDimension-1,0);
		      
		// BC
		m_BCRegions.resize(2);
		m_BCValues.resize(2);
		m_BCRegion.SetIndex(InputDimension-1,0);
		m_BCRegions[0]=m_BCRegion;
		m_BCValues[0]=-m_WeightRatio[2][0];
		m_BCRegion.SetIndex(InputDimension-1,1);
		m_BCRegions[1]=m_BCRegion;
		m_BCValues[1]= -m_WeightRatio[2][0];

		// Second part
		m_SecondRegion.SetSize(InputDimension-1,2);
		m_SecondRegion.SetIndex(InputDimension-1,1);
		      
		// BC2
		m_BC2Regions.resize(0);

		break;
	      }
	    case 3:
	      {
		// Lower end 
		m_FirstRegion.SetSize(InputDimension-1,0);
		      
		// BC
		m_BCRegions.resize(2);
		m_BCValues.resize(2);
		m_BCRegion.SetIndex(InputDimension-1,0);
		m_BCRegions[0]=m_BCRegion;
		m_BCValues[0]=-m_WeightRatio[2][0];
		m_BCRegion.SetIndex(InputDimension-1,1);
		m_BCRegions[1]=m_BCRegion;
		m_BCValues[1]= -m_WeightRatio[2][0];

		// Second part
		m_SecondRegion.SetSize(InputDimension-1,2);
		m_SecondRegion.SetIndex(InputDimension-1,1);
		      
		// BC2
		m_BC2Regions.resize(1);
		m_BC2Values.resize(1);
		m_BCRegion.SetIndex(InputDimension-1,0);
		m_BC2Regions[0]=m_BCRegion;
		m_BC2Values[0]=1;

		break;
	      }
		    
	    default:
	      {
		DD("supportindex > 3 ???");
		DD(m_SupportRegion.GetIndex(InputDimension-1));
	      }
	    } // end switch index
		
	  break;
	}

      case 1:
	{
	  // ----------------------------------------------------------------------
	  // The egg with 5 internal CP (starting from inhale)
	  // Periodic, constrained to zero at the reference 
	  // at position 3 and
	  // Coeff      row  BC1 BC2  0  BC3  1   2   3  BC4
	  // PaddedCoeff  R:  0   1   2   3   4   5   6   7
	  // BC1= R5
	  // BC2= R6
	  // BC3= -weights[3]/weights[1] ( R2+R5 ) - R4 
	  // BC4= R2
	  // ---------------------------------------------------------------------
	  switch(m_SupportRegion.GetIndex(InputDimension-1)) 
	    {
	    case 0:
	      {
		// Lower end 
		m_FirstRegion.SetSize(InputDimension-1,2);
		m_FirstRegion.SetIndex(InputDimension-1,2);
		
		// BC
		m_BCRegions.resize(0);
		     
		// Second part
		m_SecondRegion.SetSize(InputDimension-1,1);
		m_SecondRegion.SetIndex(InputDimension-1,0);
 
		// BC
		m_BC2Regions.resize(3);
		m_BC2Values.resize(3);
		m_BCRegion.SetIndex(InputDimension-1,0);
		m_BC2Regions[0]=m_BCRegion;
		m_BC2Values[0]=-m_WeightRatio[3][1];
		m_BCRegion.SetIndex(InputDimension-1,1);
		m_BC2Regions[1]=m_BCRegion;
		m_BC2Values[1]= -1;
		m_BCRegion.SetIndex(InputDimension-1,2);
		m_BC2Regions[2]=m_BCRegion;
		m_BC2Values[2]=-m_WeightRatio[3][1];
		    		      
		break;
	      }
	    case 1:
	      {
		// Lower end 
		m_FirstRegion.SetSize(InputDimension-1,1);
		m_FirstRegion.SetIndex(InputDimension-1,3);
		      
		// BC
		m_BCRegions.resize(0);
		      
		// Second part
		m_SecondRegion.SetSize(InputDimension-1,1);
		m_SecondRegion.SetIndex(InputDimension-1,0);
		      
		// BC
		m_BC2Regions.resize(3);
		m_BC2Values.resize(3);
		m_BCRegion.SetIndex(InputDimension-1,0);
		m_BC2Regions[0]=m_BCRegion;
		m_BC2Values[0]=-m_WeightRatio[3][1];
		m_BCRegion.SetIndex(InputDimension-1,1);
		m_BC2Regions[1]=m_BCRegion;
		m_BC2Values[1]= -1;
		m_BCRegion.SetIndex(InputDimension-1,2);
		m_BC2Regions[2]=m_BCRegion;
		m_BC2Values[2]=-m_WeightRatio[3][1];

		// Third Part
		m_FirstRegion.SetSize(InputDimension-1,1);
		m_FirstRegion.SetIndex(InputDimension-1,1);
		      
		break;
	      }
	    case 2:
	      {
		// Lower end 
		m_FirstRegion.SetSize(InputDimension-1,1);
		m_FirstRegion.SetIndex(InputDimension-1,0);
		      
		// BC
		m_BCRegions.resize(3);
		m_BCValues.resize(3);
		m_BCRegion.SetIndex(InputDimension-1,0);
		m_BCRegions[0]=m_BCRegion;
		m_BCValues[0]=-m_WeightRatio[3][1];
		m_BCRegion.SetIndex(InputDimension-1,1);
		m_BCRegions[1]=m_BCRegion;
		m_BCValues[1]= -1;
		m_BCRegion.SetIndex(InputDimension-1,2);
		m_BCRegions[2]=m_BCRegion;
		m_BCValues[2]=-m_WeightRatio[3][1];

		// Second part
		m_SecondRegion.SetSize(InputDimension-1,2);
		m_SecondRegion.SetIndex(InputDimension-1,1);
		      
		// BC2
		m_BC2Regions.resize(0);

		break;
	      }
	    case 3:
	      {
		// Lower end 
		m_FirstRegion.SetSize(InputDimension-1,0);
		      
		// BC
		m_BCRegions.resize(3);
		m_BCValues.resize(3);
		m_BCRegion.SetIndex(InputDimension-1,0);
		m_BCRegions[0]=m_BCRegion;
		m_BCValues[0]=-m_WeightRatio[3][1];
		m_BCRegion.SetIndex(InputDimension-1,1);
		m_BCRegions[1]=m_BCRegion;
		m_BCValues[1]= -1;
		m_BCRegion.SetIndex(InputDimension-1,2);
		m_BCRegions[2]=m_BCRegion;
		m_BCValues[2]=-m_WeightRatio[3][1];

		// Second part
		m_SecondRegion.SetSize(InputDimension-1,3);
		m_SecondRegion.SetIndex(InputDimension-1,1);
		      
		// BC2
		m_BC2Regions.resize(0);
	
		break;
	      }
	    case 4:
	      {
		// Lower end 
		m_FirstRegion.SetSize(InputDimension-1,3);
		m_FirstRegion.SetIndex(InputDimension-1,1);
		      
		// BC
		m_BCRegions.resize(0);
	
		// Second part
		m_SecondRegion.SetSize(InputDimension-1,1);
		m_SecondRegion.SetIndex(InputDimension-1,0);
		      
		// BC2
		m_BC2Regions.resize(0);
	
		break;
	      }
		    
	    default:
	      {
		DD("supportindex > 3 ???");
		DD(m_SupportRegion.GetIndex(InputDimension-1));
	      }
	    } // end switch index
		
	  break;
	}
// 	  // ----------------------------------------------------------------------
// 	  // The egg with 5 internal CP: 
// 	  // Periodic, C2 smooth everywhere and constrained to zero at the reference
// 	  // Coeff       row R5 BC1  0   1   2   3  BC2  R2  
// 	  // PaddedCoeff R:  0   1   2   3   4   5   6   7    
// 	  // BC1= -weights[2]/weights[0] ( R2+R5)
// 	  // BC2= BC1  
// 	  // ---------------------------------------------------------------------
// 	  switch(m_SupportRegion.GetIndex(InputDimension-1))
// 	    {
// 	    case 0:
// 	      {
// 		// Lower end 
// 		m_FirstRegion.SetSize(InputDimension-1,1);
// 		m_FirstRegion.SetIndex(InputDimension-1,3);
			
// 		// BC
// 		m_BCRegions.resize(2);
// 		m_BCValues.resize(2);
// 		m_BCRegion.SetIndex(InputDimension-1,0);
// 		m_BCRegions[0]=m_BCRegion;
// 		m_BCValues[0]=-m_WeightRatio[2][0];
// 		m_BCRegion.SetIndex(InputDimension-1,3);
// 		m_BCRegions[1]=m_BCRegion;
// 		m_BCValues[1]=-m_WeightRatio[2][0];
	      
// 		// Second part
// 		m_SecondRegion.SetSize(InputDimension-1,2);
// 		m_SecondRegion.SetIndex(InputDimension-1,0);

// 		// BC2
// 		m_BC2Regions.resize(0);

// 		break;
// 	      }
// 	    case 1:
// 	      {
// 		// Lower end 
// 		m_FirstRegion.SetSize(InputDimension-1,0);
		      
// 		// BC
// 		m_BCRegions.resize(2);
// 		m_BCValues.resize(2);
// 		m_BCRegion.SetIndex(InputDimension-1,0);
// 		m_BCRegions[0]=m_BCRegion;
// 		m_BCValues[0]=-m_WeightRatio[2][0];
// 		m_BCRegion.SetIndex(InputDimension-1,3);
// 		m_BCRegions[1]=m_BCRegion;
// 		m_BCValues[1]=-m_WeightRatio[2][0];
	      
// 		// Second part
// 		m_SecondRegion.SetSize(InputDimension-1,3);
// 		m_SecondRegion.SetIndex(InputDimension-1,0);

// 		// BC2
// 		m_BC2Regions.resize(0);

// 		break;
// 	      }
// 	    case 2:
// 	      {
// 		// Lower end 
// 		m_FirstRegion.SetSize(InputDimension-1,4);
// 		m_FirstRegion.SetIndex(InputDimension-1,0);
	      
// 		// BC
// 		m_BCRegions.resize(0);
	    	      
// 		// Second part
// 		m_SecondRegion.SetSize(InputDimension-1,0);
	
// 		// BC2
// 		m_BC2Regions.resize(0);
	      
// 		break;
// 	      }
// 	    case 3:
// 	      {
// 		// Lower end 
// 		m_FirstRegion.SetSize(InputDimension-1,3);
// 		m_FirstRegion.SetIndex(InputDimension-1,1);
		
// 		// BC
// 		m_BCRegions.resize(2);
// 		m_BCValues.resize(2);
// 		m_BCRegion.SetIndex(InputDimension-1,0);
// 		m_BCRegions[0]=m_BCRegion;
// 		m_BCValues[0]=-m_WeightRatio[2][0];
// 		m_BCRegion.SetIndex(InputDimension-1,3);
// 		m_BCRegions[1]=m_BCRegion;
// 		m_BCValues[1]=-m_WeightRatio[2][0];
		
// 		// Second part
// 		m_SecondRegion.SetSize(InputDimension-1,0);

// 		// BC2
// 		m_BC2Regions.resize(0);
		 
// 		break;
// 	      }

// 	    case 4:
// 	      {
// 		// Lower end 
// 		m_FirstRegion.SetSize(InputDimension-1,2);
// 		m_FirstRegion.SetIndex(InputDimension-1,2);
		
// 		// BC
// 		m_BCRegions.resize(2);
// 		m_BCValues.resize(2);
// 		m_BCRegion.SetIndex(InputDimension-1,0);
// 		m_BCRegions[0]=m_BCRegion;
// 		m_BCValues[0]=-m_WeightRatio[2][0];
// 		m_BCRegion.SetIndex(InputDimension-1,3);
// 		m_BCRegions[1]=m_BCRegion;
// 		m_BCValues[1]=-m_WeightRatio[2][0];

// 		// Second part
// 		m_SecondRegion.SetSize(InputDimension-1,1);
// 		m_SecondRegion.SetIndex(InputDimension-1,0);

// 		// BC2
// 		m_BC2Regions.resize(0);

// 		break;
// 	      }

// 	    default:
// 	      {
// 		DD("supportindex > 4 ???");
// 		DD(m_SupportRegion.GetIndex(InputDimension-1));
// 		DD(m_TransformShape);
// 	      }
// 	    }// end swith index
// 	  break;
	  
// 	} // end case 1 shape

      case 2:
	{
	  // ----------------------------------------------------------------------
	  // The rabbit with 4 internal CP
	  // Periodic, constrained to zero at the reference 
	  // at position 3 and the extremes fixed with anit-symm bc
	  // Coeff      row  BC1  0   1  BC2  2  BC3 BC4 
	  // PaddedCoeff  R:  0   1   2   3   4   5   6 
	  // BC1= 2*R1-R0
	  // BC2= -weights[2]/weights[0] ( R2+R4 )
	  // BC3=  R1
	  // BC4= 2*R1-R4
	  // ---------------------------------------------------------------------
	  switch(m_SupportRegion.GetIndex(InputDimension-1)) 
	    {
	    case 0:
	      {
		// Lower end 
		m_FirstRegion.SetSize(InputDimension-1,0);

		// BC
		m_BCRegions.resize(2);
		m_BCValues.resize(2);
		m_BCRegion.SetIndex(InputDimension-1,0);
		m_BCRegions[0]=m_BCRegion;
		m_BCValues[0]=2;
		m_BCRegion.SetIndex(InputDimension-1,1);
		m_BCRegions[1]=m_BCRegion;
		m_BCValues[1]= -1;
		     
		// Second part
		m_SecondRegion.SetSize(InputDimension-1,2);
		m_SecondRegion.SetIndex(InputDimension-1,0);
 
		// BC
		m_BC2Regions.resize(2);
		m_BC2Values.resize(2);
		m_BCRegion.SetIndex(InputDimension-1,1);
		m_BC2Regions[0]=m_BCRegion;
		m_BC2Values[0]=-m_WeightRatio[2][0];
		m_BCRegion.SetIndex(InputDimension-1,2);
		m_BC2Regions[1]=m_BCRegion;
		m_BC2Values[1]= -m_WeightRatio[2][0];
		    		      
		break;
	      }
	    case 1:
	      {
		// Lower end 
		m_FirstRegion.SetSize(InputDimension-1,2);
		m_FirstRegion.SetIndex(InputDimension-1,0);
		      
		// BC
		m_BCRegions.resize(2);
		m_BCValues.resize(2);
		m_BCRegion.SetIndex(InputDimension-1,1);
		m_BCRegions[0]=m_BCRegion;
		m_BCValues[0]=-m_WeightRatio[2][0];
		m_BCRegion.SetIndex(InputDimension-1,2);
		m_BCRegions[1]=m_BCRegion;
		m_BCValues[1]= -m_WeightRatio[2][0];
		      
		// Second part
		m_SecondRegion.SetSize(InputDimension-1,1);
		m_SecondRegion.SetIndex(InputDimension-1,2);
		      
		// BC2
		m_BC2Regions.resize(0);
		      
		break;
	      }
	    case 2:
	      {
		// Lower end 
		m_FirstRegion.SetSize(InputDimension-1,1);
		m_FirstRegion.SetIndex(InputDimension-1,1);
		      
		// BC
		m_BCRegions.resize(2);
		m_BCValues.resize(2);
		m_BCRegion.SetIndex(InputDimension-1,1);
		m_BCRegions[0]=m_BCRegion;
		m_BCValues[0]=-m_WeightRatio[2][0];
		m_BCRegion.SetIndex(InputDimension-1,2);
		m_BCRegions[1]=m_BCRegion;
		m_BCValues[1]= -m_WeightRatio[2][0];

		// Second part
		m_SecondRegion.SetSize(InputDimension-1,1);
		m_SecondRegion.SetIndex(InputDimension-1,2);
		      
		// BC2
		m_BC2Regions.resize(1);
		m_BC2Values.resize(1);
		m_BCRegion.SetIndex(InputDimension-1,0);
		m_BC2Regions[0]=m_BCRegion;
		m_BC2Values[0]=1;
			   		      
		break;
	      }
	    case 3:
	      {
		// Lower end 
		m_FirstRegion.SetSize(InputDimension-1,0);
		      
		// BC
		m_BCRegions.resize(2);
		m_BCValues.resize(2);
		m_BCRegion.SetIndex(InputDimension-1,1);
		m_BCRegions[0]=m_BCRegion;
		m_BCValues[0]=-m_WeightRatio[2][0];
		m_BCRegion.SetIndex(InputDimension-1,2);
		m_BCRegions[1]=m_BCRegion;
		m_BCValues[1]= -m_WeightRatio[2][0];

		// Second part
		m_SecondRegion.SetSize(InputDimension-1,1);
		m_SecondRegion.SetIndex(InputDimension-1,2);
		      
		// BC2
		m_BC2Regions.resize(1);
		m_BC2Values.resize(1);
		m_BCRegion.SetIndex(InputDimension-1,0);
		m_BC2Regions[0]=m_BCRegion;
		m_BC2Values[0]=1;

		// BC3
		m_BC3Regions.resize(2);
		m_BC3Values.resize(2);
		m_BCRegion.SetIndex(InputDimension-1,0);
		m_BC3Regions[0]=m_BCRegion;
		m_BC3Values[0]=2;
		m_BCRegion.SetIndex(InputDimension-1,2);
		m_BC3Regions[1]=m_BCRegion;
		m_BC3Values[1]= -1;

		break;
	      }
		    
	    default:
	      {
		DD("supportindex > 3 ???");
		DD(m_SupportRegion.GetIndex(InputDimension-1));
	      }
	    } // end switch index
		
	  break;
	} // end case 2 shape
      case 3: 
	{
	  // ----------------------------------------------------------------------
	  // The rabbit with 5 internal CP
	  // Periodic, constrained to zero at the reference at position 3.5 
	  // and the extremes fixed with anti-symmetrical BC
	  // Coeff      row  BC1  0   1  BC2  2   3  BC3 BC4
	  // PaddedCoeff  R:  0   1   2   3   4   5   6   7    
	  // BC1= 2*R1-R2
	  // BC2= -weights[3]/weights[1] ( R2+R5 ) - R4
	  // BC3= R1
	  // BC4= 2*R1-R5
	  // ---------------------------------------------------------------------
	  switch(m_SupportRegion.GetIndex(InputDimension-1)) 
	    {
	    case 0:
	      {
		// Lower end 
		m_FirstRegion.SetSize(InputDimension-1,0);
			    
		// BC
		m_BCRegions.resize(2);
		m_BCValues.resize(2);
		m_BCRegion.SetIndex(InputDimension-1,0);
		m_BCRegions[0]=m_BCRegion;
		m_BCValues[0]=2.;
		m_BCRegion.SetIndex(InputDimension-1,1);
		m_BCRegions[1]=m_BCRegion;
		m_BCValues[1]= -1.;

		// Second part
		m_SecondRegion.SetSize(InputDimension-1,2);
		m_SecondRegion.SetIndex(InputDimension-1,0);
	
		// BC2
		m_BC2Regions.resize(3);
		m_BC2Values.resize(3);
		m_BCRegion.SetIndex(InputDimension-1,1);
		m_BC2Regions[0]=m_BCRegion;
		m_BC2Values[0]=-m_WeightRatio[3][1];
		m_BCRegion.SetIndex(InputDimension-1,2);
		m_BC2Regions[1]=m_BCRegion;
		m_BC2Values[1]= -1;
		m_BCRegion.SetIndex(InputDimension-1,3);
		m_BC2Regions[2]=m_BCRegion;
		m_BC2Values[2]=-m_WeightRatio[3][1];	  
		      
		break;
	      }
	    case 1:
	      {
		// Lower end 
		m_FirstRegion.SetSize(InputDimension-1,2);
		m_FirstRegion.SetIndex(InputDimension-1,0);
		      
		// BC
		m_BCRegions.resize(3);
		m_BCValues.resize(3);
		m_BCRegion.SetIndex(InputDimension-1,1);
		m_BCRegions[0]=m_BCRegion;
		m_BCValues[0]=-m_WeightRatio[3][1];
		m_BCRegion.SetIndex(InputDimension-1,2);
		m_BCRegions[1]=m_BCRegion;
		m_BCValues[1]= -1;
		m_BCRegion.SetIndex(InputDimension-1,3);
		m_BCRegions[2]=m_BCRegion;
		m_BCValues[2]=-m_WeightRatio[3][1];
		      
		// Second part
		m_SecondRegion.SetSize(InputDimension-1,1);
		m_SecondRegion.SetIndex(InputDimension-1,2);
		      
		// BC2
		m_BC2Regions.resize(0);
		      
		break;
	      }
	    case 2:
	      {
		// Lower end 
		m_FirstRegion.SetSize(InputDimension-1,1);
		m_FirstRegion.SetIndex(InputDimension-1,1);
		      
		// BC
		m_BCRegions.resize(3);
		m_BCValues.resize(3);
		m_BCRegion.SetIndex(InputDimension-1,1);
		m_BCRegions[0]=m_BCRegion;
		m_BCValues[0]=-m_WeightRatio[3][1];
		m_BCRegion.SetIndex(InputDimension-1,2);
		m_BCRegions[1]=m_BCRegion;
		m_BCValues[1]= -1;
		m_BCRegion.SetIndex(InputDimension-1,3);
		m_BCRegions[2]=m_BCRegion;
		m_BCValues[2]=-m_WeightRatio[3][1];	  
		      
		// Second part
		m_SecondRegion.SetSize(InputDimension-1,2);
		m_SecondRegion.SetIndex(InputDimension-1,2);
		      
		// BC2
		m_BC2Regions.resize(0);
		      
		break;
	      }
	    case 3:
	      {
		// Lower end 
		m_FirstRegion.SetSize(InputDimension-1,0);
		      
		// BC
		m_BCRegions.resize(3);
		m_BCValues.resize(3);
		m_BCRegion.SetIndex(InputDimension-1,1);
		m_BCRegions[0]=m_BCRegion;
		m_BCValues[0]=-m_WeightRatio[3][1];
		m_BCRegion.SetIndex(InputDimension-1,2);
		m_BCRegions[1]=m_BCRegion;
		m_BCValues[1]= -1;
		m_BCRegion.SetIndex(InputDimension-1,3);
		m_BCRegions[2]=m_BCRegion;
		m_BCValues[2]=-m_WeightRatio[3][1];	  	  
		      
		// Second part
		m_SecondRegion.SetSize(InputDimension-1,2);
		m_SecondRegion.SetIndex(InputDimension-1,2);
		      
		// BC2
		m_BC2Regions.resize(1);
		m_BC2Values.resize(1);
		m_BCRegion.SetIndex(InputDimension-1,0);
		m_BC2Regions[0]=m_BCRegion;
		m_BC2Values[0]=1.;
				      
		break;
	      }
		    
	    case 4:
	      {
		// Lower end 
		m_FirstRegion.SetSize(InputDimension-1,2);
		m_FirstRegion.SetIndex(InputDimension-1,2);
		      
		// BC2
		m_BCRegions.resize(1);
		m_BCValues.resize(1);
		m_BCRegion.SetIndex(InputDimension-1,0);
		m_BCRegions[0]=m_BCRegion;
		m_BCValues[0]=1.;
			      
		// Second part
		m_SecondRegion.SetSize(InputDimension-1,0);
				      
		// BC2
		m_BC2Regions.resize(2);
		m_BC2Values.resize(2);
		m_BCRegion.SetIndex(InputDimension-1,0);
		m_BC2Regions[0]=m_BCRegion;
		m_BC2Values[0]=2;
		m_BCRegion.SetIndex(InputDimension-1,3);
		m_BC2Regions[1]=m_BCRegion;
		m_BC2Values[1]=-1;
		      
		break;
		      
	      }
	    default:
	      {
		DD("supportindex > 4 ???");
		DD(m_SupportRegion.GetIndex(InputDimension-1));
	      }
	    } // end switch index
	  
	  break;
	} // end case 3 shape

      case 4: 
	{
	  // ----------------------------------------------------------------------
	  // The sputnik with 4 internal CP
	  // Periodic, constrained to zero at the reference 
	  // at position 3 and one indepent extremes copied
	  // Coeff     row BC1  0   1  BC2  2  BC2  3
	  // PaddedCoeff R: 0   1   2   3   4   5   6      
	  // BC1= R6
	  // BC2= -weights[2]/weights[0] ( R2+R4 )
	  // BC3=  weights[2]/weights[0] ( R2-R4) + R1
	  // ---------------------------------------------------------------------
	  switch(m_SupportRegion.GetIndex(InputDimension-1)) 
	    {
	    case 0:
	      {
		// Lower end 
		m_FirstRegion.SetSize(InputDimension-1,0);

		// BC
		m_BCRegions.resize(1);
		m_BCValues.resize(1);
		m_BCRegion.SetIndex(InputDimension-1,3);
		m_BCRegions[0]=m_BCRegion;
		m_BCValues[0]=1;

		// Second part
		m_SecondRegion.SetSize(InputDimension-1,2);
		m_SecondRegion.SetIndex(InputDimension-1,0);		
		
		// BC2
		m_BC2Regions.resize(2);
		m_BC2Values.resize(2);
		m_BCRegion.SetIndex(InputDimension-1,1);
		m_BC2Regions[0]=m_BCRegion;
		m_BC2Values[0]=-m_WeightRatio[2][0];
		m_BCRegion.SetIndex(InputDimension-1,2);
		m_BC2Regions[1]=m_BCRegion;
		m_BC2Values[1]= -m_WeightRatio[2][0];
		      
		break;
	      }
	    case 1:
	      {
		// Lower end 
		m_FirstRegion.SetSize(InputDimension-1,2);
		m_FirstRegion.SetIndex(InputDimension-1,0);
		      
		// BC
		m_BCRegions.resize(2);
		m_BCValues.resize(2);
		m_BCRegion.SetIndex(InputDimension-1,1);
		m_BCRegions[0]=m_BCRegion;
		m_BCValues[0]=-m_WeightRatio[2][0];
		m_BCRegion.SetIndex(InputDimension-1,2);
		m_BCRegions[1]=m_BCRegion;
		m_BCValues[1]= -m_WeightRatio[2][0];
		      
		// Second part
		m_SecondRegion.SetSize(InputDimension-1,1);
		m_SecondRegion.SetIndex(InputDimension-1,2);
		      
		// BC2
		m_BC2Regions.resize(0);
		      
		break;
	      }
	    case 2:
	      {
		// Lower end 
		m_FirstRegion.SetSize(InputDimension-1,1);
		m_FirstRegion.SetIndex(InputDimension-1,1);
		      
		// BC
		m_BCRegions.resize(2);
		m_BCValues.resize(2);
		m_BCRegion.SetIndex(InputDimension-1,1);
		m_BCRegions[0]=m_BCRegion;
		m_BCValues[0]=-m_WeightRatio[2][0];
		m_BCRegion.SetIndex(InputDimension-1,2);
		m_BCRegions[1]=m_BCRegion;
		m_BCValues[1]= -m_WeightRatio[2][0];

		// Second part
		m_SecondRegion.SetSize(InputDimension-1,1);
		m_SecondRegion.SetIndex(InputDimension-1,2);
		      
		// BC2
		m_BC2Regions.resize(3);
		m_BC2Values.resize(3);
		m_BCRegion.SetIndex(InputDimension-1,0);
		m_BC2Regions[0]=m_BCRegion;
		m_BC2Values[0]=1;
		m_BCRegion.SetIndex(InputDimension-1,1);
		m_BC2Regions[1]=m_BCRegion;
		m_BC2Values[1]=m_WeightRatio[2][0];
		m_BCRegion.SetIndex(InputDimension-1,2);
		m_BC2Regions[2]=m_BCRegion;
		m_BC2Values[2]=-m_WeightRatio[2][0];	  
		   		      
		break;
	      }
	    case 3:
	      {
		// Lower end 
		m_FirstRegion.SetSize(InputDimension-1,0);
		      
		// BC
		m_BCRegions.resize(2);
		m_BCValues.resize(2);
		m_BCRegion.SetIndex(InputDimension-1,1);
		m_BCRegions[0]=m_BCRegion;
		m_BCValues[0]=-m_WeightRatio[2][0];
		m_BCRegion.SetIndex(InputDimension-1,2);
		m_BCRegions[1]=m_BCRegion;
		m_BCValues[1]= -m_WeightRatio[2][0];

		// Second part
		m_SecondRegion.SetSize(InputDimension-1,1);
		m_SecondRegion.SetIndex(InputDimension-1,2);
		      
		// BC2
		m_BC2Regions.resize(3);
		m_BC2Values.resize(3);
		m_BCRegion.SetIndex(InputDimension-1,0);
		m_BC2Regions[0]=m_BCRegion;
		m_BC2Values[0]=1;
		m_BCRegion.SetIndex(InputDimension-1,1);
		m_BC2Regions[1]=m_BCRegion;
		m_BC2Values[1]=m_WeightRatio[2][0];
		m_BCRegion.SetIndex(InputDimension-1,2);
		m_BC2Regions[2]=m_BCRegion;
		m_BC2Values[2]=-m_WeightRatio[2][0];

		// Third part
		m_ThirdRegion.SetSize(InputDimension-1,1);
		m_ThirdRegion.SetIndex(InputDimension-1,3);

		break;
	      }
		    
	    default:
	      {
		DD("supportindex > 3 ???");
		DD(m_SupportRegion.GetIndex(InputDimension-1));
	      }
	    } // end switch index
		
	  break;
	} // end case 4 shape
       
      case 5: 
	{
	  // ----------------------------------------------------------------------
	  // The sputnik with 5 internal CP
	  // Periodic, constrained to zero at the reference 
	  // at position 3 and one indepent extreme
	  // Coeff     row BC1  0   1  BC2  2   3  BC3  4
	  // PaddedCoeff R: 0   1   2   3   4   5   6   7   
	  // BC1= R2 + R5 - R7
	  // BC2= -weights[3]/weights[1] ( R2+R5 ) - R4
	  // BC3= R1 + 0.5 R2 - 0.5 R7
	  // ----------------------------------------------------------------------
	  switch(m_SupportRegion.GetIndex(InputDimension-1)) 
	    {
	    case 0:
	      {
		// Lower end 
		m_FirstRegion.SetSize(InputDimension-1,0);

		// BC
		m_BCRegions.resize(3);
		m_BCValues.resize(3);
		m_BCRegion.SetIndex(InputDimension-1,1);
		m_BCRegions[0]=m_BCRegion;
		m_BCValues[0]=1.;
		m_BCRegion.SetIndex(InputDimension-1,3);
		m_BCRegions[1]=m_BCRegion;
		m_BCValues[1]=1.;
		m_BCRegion.SetIndex(InputDimension-1,4);
		m_BCRegions[2]=m_BCRegion;
		m_BCValues[2]=-1.;

		// Second part
		m_SecondRegion.SetSize(InputDimension-1,2);
		m_SecondRegion.SetIndex(InputDimension-1,0);		
		
		// BC2
		m_BC2Regions.resize(3);
		m_BC2Values.resize(3);
		m_BCRegion.SetIndex(InputDimension-1,1);
		m_BC2Regions[0]=m_BCRegion;
		m_BC2Values[0]=-m_WeightRatio[3][1];
		m_BCRegion.SetIndex(InputDimension-1,2);
		m_BC2Regions[1]=m_BCRegion;
		m_BC2Values[1]= -1;
		m_BCRegion.SetIndex(InputDimension-1,3);
		m_BC2Regions[2]=m_BCRegion;
		m_BC2Values[2]=-m_WeightRatio[3][1];
		      
		break;
	      }
	    case 1:
	      {
		// Lower end 
		m_FirstRegion.SetSize(InputDimension-1,2);
		m_FirstRegion.SetIndex(InputDimension-1,0);
		      
		// BC
		m_BCRegions.resize(3);
		m_BCValues.resize(3);
		m_BCRegion.SetIndex(InputDimension-1,1);
		m_BCRegions[0]=m_BCRegion;
		m_BCValues[0]=-m_WeightRatio[3][1];
		m_BCRegion.SetIndex(InputDimension-1,2);
		m_BCRegions[1]=m_BCRegion;
		m_BCValues[1]= -1;
		m_BCRegion.SetIndex(InputDimension-1,3);
		m_BCRegions[2]=m_BCRegion;
		m_BCValues[2]=-m_WeightRatio[3][1];
		      
		// Second part
		m_SecondRegion.SetSize(InputDimension-1,1);
		m_SecondRegion.SetIndex(InputDimension-1,2);
		      
		// BC2
		m_BC2Regions.resize(0);
		      
		break;
	      }
	    case 2:
	      {
		// Lower end 
		m_FirstRegion.SetSize(InputDimension-1,1);
		m_FirstRegion.SetIndex(InputDimension-1,1);
		      
		// BC
		m_BCRegions.resize(3);
		m_BCValues.resize(3);
		m_BCRegion.SetIndex(InputDimension-1,1);
		m_BCRegions[0]=m_BCRegion;
		m_BCValues[0]=-m_WeightRatio[3][1];
		m_BCRegion.SetIndex(InputDimension-1,2);
		m_BCRegions[1]=m_BCRegion;
		m_BCValues[1]= -1;
		m_BCRegion.SetIndex(InputDimension-1,3);
		m_BCRegions[2]=m_BCRegion;
		m_BCValues[2]=-m_WeightRatio[3][1];

		// Second part
		m_SecondRegion.SetSize(InputDimension-1,2);
		m_SecondRegion.SetIndex(InputDimension-1,2);
		      
		// BC2
		m_BC2Regions.resize(0);
			   		      
		break;
	      }
	    case 3:
	      {
		// Lower end 
		m_FirstRegion.SetSize(InputDimension-1,0);
		      
		// BC
		m_BCRegions.resize(3);
		m_BCValues.resize(3);
		m_BCRegion.SetIndex(InputDimension-1,1);
		m_BCRegions[0]=m_BCRegion;
		m_BCValues[0]=-m_WeightRatio[3][1];
		m_BCRegion.SetIndex(InputDimension-1,2);
		m_BCRegions[1]=m_BCRegion;
		m_BCValues[1]= -1;
		m_BCRegion.SetIndex(InputDimension-1,3);
		m_BCRegions[2]=m_BCRegion;
		m_BCValues[2]=-m_WeightRatio[3][1];

		// Second part
		m_SecondRegion.SetSize(InputDimension-1,2);
		m_SecondRegion.SetIndex(InputDimension-1,2);
		      
		// BC2
		m_BC2Regions.resize(3);
		m_BC2Values.resize(3);
		m_BCRegion.SetIndex(InputDimension-1,0);
		m_BC2Regions[0]=m_BCRegion;
		m_BC2Values[0]=1;
		m_BCRegion.SetIndex(InputDimension-1,1);
		m_BC2Regions[1]=m_BCRegion;
		m_BC2Values[1]=0.5;
		m_BCRegion.SetIndex(InputDimension-1,4);
		m_BC2Regions[2]=m_BCRegion;
		m_BC2Values[2]=-0.5;
	
		break;
	      }
	    case 4:
	      {
		// Lower end 
		m_FirstRegion.SetSize(InputDimension-1,2);
		m_FirstRegion.SetIndex(InputDimension-1,2);		
      
		// BC
		m_BCRegions.resize(3);
		m_BCValues.resize(3);
		m_BCRegion.SetIndex(InputDimension-1,0);
		m_BCRegions[0]=m_BCRegion;
		m_BCValues[0]=1;
		m_BCRegion.SetIndex(InputDimension-1,1);
		m_BCRegions[1]=m_BCRegion;
		m_BCValues[1]=0.5;
		m_BCRegion.SetIndex(InputDimension-1,4);
		m_BCRegions[2]=m_BCRegion;
		m_BCValues[2]=-0.5;
	
		// Second part
		m_SecondRegion.SetSize(InputDimension-1,1);
		m_SecondRegion.SetIndex(InputDimension-1,4);
		      
		// BC2
		m_BC2Regions.resize(0);

		break;
	      }
	    default:
	      {
		DD("supportindex > 4 ???");
		DD(m_SupportRegion.GetIndex(InputDimension-1));
	      }
	    } // end switch index
		
	  break;
	} // end case 5 shape



     
      case 6: 
	{
	  // ----------------------------------------------------------------------
	  // The diamond with 4 internal CP:
	  // Periodic, constrained to zero at the reference at position 3
	  // Coeff      row 0   1   2  BC1  3  BC2  4
	  // PaddedCoeff  R:0   1   2   3   4   5   6    
	  // BC1= -weights[2]/weights[0] ( R2+R4 )
	  // BC2= weights[2]/weights[0]  ( R0+R2-R4-R6 ) + R1
	  // ---------------------------------------------------------------------
	  switch(m_SupportRegion.GetIndex(InputDimension-1)) 
	    {
	    case 0:
	      {
		// Lower end 
		m_FirstRegion.SetSize(InputDimension-1,3);
		m_FirstRegion.SetIndex(InputDimension-1,0);
		
		// BC
		m_BCRegions.resize(2);
		m_BCValues.resize(2);
		m_BCRegion.SetIndex(InputDimension-1,2);
		m_BCRegions[0]=m_BCRegion;
		m_BCValues[0]=-m_WeightRatio[2][0];
		m_BCRegion.SetIndex(InputDimension-1,3);
		m_BCRegions[1]=m_BCRegion;
		m_BCValues[1]=-m_WeightRatio[2][0];
		      
		// Second part
		m_SecondRegion.SetSize(InputDimension-1,0);

		// BC2
		m_BC2Regions.resize(0);
		
		break;
	      }
	    case 1:
	      {
		// Lower end 
		m_FirstRegion.SetSize(InputDimension-1,2);
		m_FirstRegion.SetIndex(InputDimension-1,1);
		
		// BC
		m_BCRegions.resize(2);
		m_BCValues.resize(2);
		m_BCRegion.SetIndex(InputDimension-1,2);
		m_BCRegions[0]=m_BCRegion;
		m_BCValues[0]=-m_WeightRatio[2][0];
		m_BCRegion.SetIndex(InputDimension-1,3);
		m_BCRegions[1]=m_BCRegion;
		m_BCValues[1]=-m_WeightRatio[2][0];  
	      
		// Second part
		m_SecondRegion.SetSize(InputDimension-1,1);
		m_SecondRegion.SetIndex(InputDimension-1,3);

		// BC2
		m_BC2Regions.resize(0);

		break;
	      }
	    case 2:
	      {
		// Lower end 
		m_FirstRegion.SetSize(InputDimension-1,1);
		m_FirstRegion.SetIndex(InputDimension-1,2);
		
		// BC
		m_BCRegions.resize(2);
		m_BCValues.resize(2);
		m_BCRegion.SetIndex(InputDimension-1,2);
		m_BCRegions[0]=m_BCRegion;
		m_BCValues[0]=-m_WeightRatio[2][0];
		m_BCRegion.SetIndex(InputDimension-1,3);
		m_BCRegions[1]=m_BCRegion;
		m_BCValues[1]=-m_WeightRatio[2][0];	  	  
	      
		// Second part
		m_SecondRegion.SetSize(InputDimension-1,1);
		m_SecondRegion.SetIndex(InputDimension-1,3);

		// BC2
		m_BC2Regions.resize(5);
		m_BC2Values.resize(5);
		m_BCRegion.SetIndex(InputDimension-1,0);
		m_BC2Regions[0]=m_BCRegion;
		m_BC2Values[0]=m_WeightRatio[2][0];
		m_BCRegion.SetIndex(InputDimension-1,1);
		m_BC2Regions[1]=m_BCRegion;
		m_BC2Values[1]=1;
		m_BCRegion.SetIndex(InputDimension-1,2);
		m_BC2Regions[2]=m_BCRegion;
		m_BC2Values[2]=m_WeightRatio[2][0];	  
		m_BCRegion.SetIndex(InputDimension-1,3);
		m_BC2Regions[3]=m_BCRegion;
		m_BC2Values[3]=-m_WeightRatio[2][0];	
		m_BCRegion.SetIndex(InputDimension-1,4);
		m_BC2Regions[4]=m_BCRegion;
		m_BC2Values[4]=-m_WeightRatio[2][0];
		
		break;
	      }
	    case 3:
	      {
		// Lower end 
		m_FirstRegion.SetSize(InputDimension-1,0);
		
		// BC
		m_BCRegions.resize(2);
		m_BCValues.resize(2);
		m_BCRegion.SetIndex(InputDimension-1,2);
		m_BCRegions[0]=m_BCRegion;
		m_BCValues[0]=-m_WeightRatio[2][0];
		m_BCRegion.SetIndex(InputDimension-1,3);
		m_BCRegions[1]=m_BCRegion;
		m_BCValues[1]=-m_WeightRatio[2][0];	  	  
	      
		// Second part
		m_SecondRegion.SetSize(InputDimension-1,1);
		m_SecondRegion.SetIndex(InputDimension-1,3);

		// BC2
		m_BC2Regions.resize(5);
		m_BC2Values.resize(5);
		m_BCRegion.SetIndex(InputDimension-1,0);
		m_BC2Regions[0]=m_BCRegion;
		m_BC2Values[0]=m_WeightRatio[2][0];
		m_BCRegion.SetIndex(InputDimension-1,1);
		m_BC2Regions[1]=m_BCRegion;
		m_BC2Values[1]=1;
		m_BCRegion.SetIndex(InputDimension-1,2);
		m_BC2Regions[2]=m_BCRegion;
		m_BC2Values[2]=m_WeightRatio[2][0];	  
		m_BCRegion.SetIndex(InputDimension-1,3);
		m_BC2Regions[3]=m_BCRegion;
		m_BC2Values[3]=-m_WeightRatio[2][0];	
		m_BCRegion.SetIndex(InputDimension-1,4);
		m_BC2Regions[4]=m_BCRegion;
		m_BC2Values[4]=-m_WeightRatio[2][0];
	
		// Third part
		m_ThirdRegion.SetSize(InputDimension-1,1);
		m_ThirdRegion.SetIndex(InputDimension-1,4);
		
		break;
	      }
	      
	    default:
	      {
		DD("supportindex > 3 ???");
		DD(m_SupportRegion.GetIndex(InputDimension-1));
	      }
	    } // end switch index
		
	  break;
	} // end case 7 shape

      case 7: 
	{
	  // ----------------------------------------------------------------------
	  // The diamond with 5 internal CP: 
	  // periodic, constrained to zero at the reference at position 3.5
	  // Coeff      row 0   1   2  BC1  3   4  BC2  5
	  // PaddedCoeff  R:0   1   2   3   4   5   6   7    
	  // BC1= -weights[3]/weights[1] ( R2+R5 ) - R4
	  // BC2= weights[2]/weights[0] ( R0+R2-R5-R7 ) + R1
	  // ---------------------------------------------------------------------
	  switch(m_SupportRegion.GetIndex(InputDimension-1)) 
	    {
	    case 0:
	      {
		// Lower end 
		m_FirstRegion.SetSize(InputDimension-1,3);
		m_FirstRegion.SetIndex(InputDimension-1,0);
		
		// BC
		m_BCRegions.resize(3);
		m_BCValues.resize(3);
		m_BCRegion.SetIndex(InputDimension-1,2);
		m_BCRegions[0]=m_BCRegion;
		m_BCValues[0]=-m_WeightRatio[3][1];
		m_BCRegion.SetIndex(InputDimension-1,3);
		m_BCRegions[1]=m_BCRegion;
		m_BCValues[1]= -1;
		m_BCRegion.SetIndex(InputDimension-1,4);
		m_BCRegions[2]=m_BCRegion;
		m_BCValues[2]=-m_WeightRatio[3][1];	  
	      
		// Second part
		m_SecondRegion.SetSize(InputDimension-1,0);

		// BC2
		m_BC2Regions.resize(0);
		
		break;
	      }
	    case 1:
	      {
		// Lower end 
		m_FirstRegion.SetSize(InputDimension-1,2);
		m_FirstRegion.SetIndex(InputDimension-1,1);
		
		// BC
		m_BCRegions.resize(3);
		m_BCValues.resize(3);
		m_BCRegion.SetIndex(InputDimension-1,2);
		m_BCRegions[0]=m_BCRegion;
		m_BCValues[0]=-m_WeightRatio[3][1];
		m_BCRegion.SetIndex(InputDimension-1,3);
		m_BCRegions[1]=m_BCRegion;
		m_BCValues[1]= -1;
		m_BCRegion.SetIndex(InputDimension-1,4);
		m_BCRegions[2]=m_BCRegion;
		m_BCValues[2]=-m_WeightRatio[3][1];	  
	      
		// Second part
		m_SecondRegion.SetSize(InputDimension-1,1);
		m_SecondRegion.SetIndex(InputDimension-1,3);

		// BC2
		m_BC2Regions.resize(0);

		break;
	      }
	    case 2:
	      {
		// Lower end 
		m_FirstRegion.SetSize(InputDimension-1,1);
		m_FirstRegion.SetIndex(InputDimension-1,2);
		
		// BC
		m_BCRegions.resize(3);
		m_BCValues.resize(3);
		m_BCRegion.SetIndex(InputDimension-1,2);
		m_BCRegions[0]=m_BCRegion;
		m_BCValues[0]=-m_WeightRatio[3][1];
		m_BCRegion.SetIndex(InputDimension-1,3);
		m_BCRegions[1]=m_BCRegion;
		m_BCValues[1]= -1;
		m_BCRegion.SetIndex(InputDimension-1,4);
		m_BCRegions[2]=m_BCRegion;
		m_BCValues[2]=-m_WeightRatio[3][1];	  	  
	      
		// Second part
		m_SecondRegion.SetSize(InputDimension-1,2);
		m_SecondRegion.SetIndex(InputDimension-1,3);

		// BC2
		m_BC2Regions.resize(0);
		
		break;
	      }
	    case 3:
	      {
		// Lower end 
		m_FirstRegion.SetSize(InputDimension-1,0);
		m_FirstRegion.SetIndex(InputDimension-1,0);
		
		// BC
		m_BCRegions.resize(3);
		m_BCValues.resize(3);
		m_BCRegion.SetIndex(InputDimension-1,2);
		m_BCRegions[0]=m_BCRegion;
		m_BCValues[0]=-m_WeightRatio[3][1];
		m_BCRegion.SetIndex(InputDimension-1,3);
		m_BCRegions[1]=m_BCRegion;
		m_BCValues[1]= -1;
		m_BCRegion.SetIndex(InputDimension-1,4);
		m_BCRegions[2]=m_BCRegion;
		m_BCValues[2]=-m_WeightRatio[3][1];	  	  
	      
		// Second part
		m_SecondRegion.SetSize(InputDimension-1,2);
		m_SecondRegion.SetIndex(InputDimension-1,3);

		// BC2
		m_BC2Regions.resize(5);
		m_BC2Values.resize(5);
		m_BCRegion.SetIndex(InputDimension-1,0);
		m_BC2Regions[0]=m_BCRegion;
		m_BC2Values[0]=m_WeightRatio[2][0];
		m_BCRegion.SetIndex(InputDimension-1,1);
		m_BC2Regions[1]=m_BCRegion;
		m_BC2Values[1]=1;
		m_BCRegion.SetIndex(InputDimension-1,2);
		m_BC2Regions[2]=m_BCRegion;
		m_BC2Values[2]=m_WeightRatio[2][0];	  
		m_BCRegion.SetIndex(InputDimension-1,4);
		m_BC2Regions[3]=m_BCRegion;
		m_BC2Values[3]=-m_WeightRatio[2][0];	
		m_BCRegion.SetIndex(InputDimension-1,5);
		m_BC2Regions[4]=m_BCRegion;
		m_BC2Values[4]=-m_WeightRatio[2][0];	
		
		break;
	      }
	      
	    case 4:
	      {
		// Lower end 
		m_FirstRegion.SetSize(InputDimension-1,2);
		m_FirstRegion.SetIndex(InputDimension-1,3);
		
		// BC
		m_BCRegions.resize(5);
		m_BCValues.resize(5);
		m_BCRegion.SetIndex(InputDimension-1,0);
		m_BCRegions[0]=m_BCRegion;
		m_BCValues[0]=m_WeightRatio[2][0];
		m_BCRegion.SetIndex(InputDimension-1,1);
		m_BCRegions[1]=m_BCRegion;
		m_BCValues[1]=1;
		m_BCRegion.SetIndex(InputDimension-1,2);
		m_BCRegions[2]=m_BCRegion;
		m_BCValues[2]=m_WeightRatio[2][0];	  
		m_BCRegion.SetIndex(InputDimension-1,4);
		m_BCRegions[3]=m_BCRegion;
		m_BCValues[3]=-m_WeightRatio[2][0];	
		m_BCRegion.SetIndex(InputDimension-1,5);
		m_BCRegions[4]=m_BCRegion;
		m_BCValues[4]=-m_WeightRatio[2][0];	
	      
		// Second part
		m_SecondRegion.SetSize(InputDimension-1,1);
		m_SecondRegion.SetIndex(InputDimension-1,5);

		// BC2
		m_BC2Regions.resize(0);
			
		break;
	      }

	    default:
	      {
		DD("supportindex > 4 ???");
		DD(m_SupportRegion.GetIndex(InputDimension-1));
	      }
	    } // end switch index
		
	  break;
	} // end case 7 shape

      case 9:
	{
	  // ----------------------------------------------------------------------
	  // The sputnik with 5 internal CP
	  // Periodic, constrained to zero at the reference 
	  // at position 3 and one indepent extreme
	  // Coeff     row BC1  0   1  BC2  2   3  BC3  4
	  // PaddedCoeff R: 0   1   2   3   4   5   6   7   
	  // BC1= -R2+R5-R7
	  // BC2= -weights[3]/weights[1] ( R2+R5 ) - R4
	  // BC3= R1
	  // ---------------------------------------------------------------------
	  switch(m_SupportRegion.GetIndex(InputDimension-1)) 
	    {
	    case 0:
	      {
		// Lower end 
		m_FirstRegion.SetSize(InputDimension-1,0);

		// BC
		m_BCRegions.resize(3);
		m_BCValues.resize(3);
		m_BCRegion.SetIndex(InputDimension-1,1);
		m_BCRegions[0]=m_BCRegion;
		m_BCValues[0]=-1;
		m_BCRegion.SetIndex(InputDimension-1,3);
		m_BCRegions[1]=m_BCRegion;
		m_BCValues[1]=1;
		m_BCRegion.SetIndex(InputDimension-1,4);
		m_BCRegions[2]=m_BCRegion;
		m_BCValues[2]=-1;

		// Second part
		m_SecondRegion.SetSize(InputDimension-1,2);
		m_SecondRegion.SetIndex(InputDimension-1,0);		
		
		// BC2
		m_BC2Regions.resize(3);
		m_BC2Values.resize(3);
		m_BCRegion.SetIndex(InputDimension-1,1);
		m_BC2Regions[0]=m_BCRegion;
		m_BC2Values[0]=-m_WeightRatio[3][1];
		m_BCRegion.SetIndex(InputDimension-1,2);
		m_BC2Regions[1]=m_BCRegion;
		m_BC2Values[1]= -1;
		m_BCRegion.SetIndex(InputDimension-1,3);
		m_BC2Regions[2]=m_BCRegion;
		m_BC2Values[2]=-m_WeightRatio[3][1];
		      
		break;
	      }
	    case 1:
	      {
		// Lower end 
		m_FirstRegion.SetSize(InputDimension-1,2);
		m_FirstRegion.SetIndex(InputDimension-1,0);
		      
		// BC
		m_BCRegions.resize(3);
		m_BCValues.resize(3);
		m_BCRegion.SetIndex(InputDimension-1,1);
		m_BCRegions[0]=m_BCRegion;
		m_BCValues[0]=-m_WeightRatio[3][1];
		m_BCRegion.SetIndex(InputDimension-1,2);
		m_BCRegions[1]=m_BCRegion;
		m_BCValues[1]= -1;
		m_BCRegion.SetIndex(InputDimension-1,3);
		m_BCRegions[2]=m_BCRegion;
		m_BCValues[2]=-m_WeightRatio[3][1];
		      
		// Second part
		m_SecondRegion.SetSize(InputDimension-1,1);
		m_SecondRegion.SetIndex(InputDimension-1,2);
		      
		// BC2
		m_BC2Regions.resize(0);
		      
		break;
	      }
	    case 2:
	      {
		// Lower end 
		m_FirstRegion.SetSize(InputDimension-1,1);
		m_FirstRegion.SetIndex(InputDimension-1,1);
		      
		// BC
		m_BCRegions.resize(3);
		m_BCValues.resize(3);
		m_BCRegion.SetIndex(InputDimension-1,1);
		m_BCRegions[0]=m_BCRegion;
		m_BCValues[0]=-m_WeightRatio[3][1];
		m_BCRegion.SetIndex(InputDimension-1,2);
		m_BCRegions[1]=m_BCRegion;
		m_BCValues[1]= -1;
		m_BCRegion.SetIndex(InputDimension-1,3);
		m_BCRegions[2]=m_BCRegion;
		m_BCValues[2]=-m_WeightRatio[3][1];

		// Second part
		m_SecondRegion.SetSize(InputDimension-1,2);
		m_SecondRegion.SetIndex(InputDimension-1,2);
		      
		// BC2
		m_BC2Regions.resize(0);
			   		      
		break;
	      }
	    case 3:
	      {
		// Lower end 
		m_FirstRegion.SetSize(InputDimension-1,0);
		      
		// BC
		m_BCRegions.resize(3);
		m_BCValues.resize(3);
		m_BCRegion.SetIndex(InputDimension-1,1);
		m_BCRegions[0]=m_BCRegion;
		m_BCValues[0]=-m_WeightRatio[3][1];
		m_BCRegion.SetIndex(InputDimension-1,2);
		m_BCRegions[1]=m_BCRegion;
		m_BCValues[1]= -1;
		m_BCRegion.SetIndex(InputDimension-1,3);
		m_BCRegions[2]=m_BCRegion;
		m_BCValues[2]=-m_WeightRatio[3][1];

		// Second part
		m_SecondRegion.SetSize(InputDimension-1,1);
		m_SecondRegion.SetIndex(InputDimension-1,2);
		      
		// BC
		m_BC2Regions.resize(1);
		m_BC2Values.resize(1);
		m_BCRegion.SetIndex(InputDimension-1,0);
		m_BC2Regions[0]=m_BCRegion;
		m_BC2Values[0]=1;
	
		break;
	      }
	    case 4:
	      {
		// Lower end 
		m_FirstRegion.SetSize(InputDimension-1,2);
		m_FirstRegion.SetIndex(InputDimension-1,2);		
      
		// BC
		m_BCRegions.resize(1);
		m_BCValues.resize(1);
		m_BCRegion.SetIndex(InputDimension-1,0);
		m_BCRegions[0]=m_BCRegion;
		m_BCValues[0]=1;
		
		// Second part
		m_SecondRegion.SetSize(InputDimension-1,1);
		m_SecondRegion.SetIndex(InputDimension-1,4);
		      
		// BC2
		m_BC2Regions.resize(0);

		break;
	      }
	    default:
	      {
		DD("supportindex > 4 ???");
		DD(m_SupportRegion.GetIndex(InputDimension-1));
	      }
	    } // end switch index
		
	  break;
	} // end case 9 shape
       

      default:
	{
	  DD ("Other shapes currently not implemented");
	}
	      
      } // end switch shape
  } // end wrap region
	
	

  template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
  void
  ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions,NOutputDimensions>
  ::TransformPointToContinuousIndex( const InputPointType & point, ContinuousIndexType & index ) const
  {
    unsigned int j;
    
    itk::Vector<double, OutputDimension> tvector;
    
    for ( j = 0; j < OutputDimension; j++ )
      {
	//JV find the index in the PADDED version
	tvector[j] = point[j] - this->m_PaddedGridOrigin[j];
      }
    
    itk::Vector<double, OutputDimension> cvector;
    
    cvector = m_PointToIndex * tvector;
    
    for ( j = 0; j < OutputDimension; j++ )
      {
	index[j] = static_cast< typename ContinuousIndexType::CoordRepType >( cvector[j] );
      }
  }
  
  
} // namespace

#endif