registration/clitkShapedBLUTSpatioTemporalDeformableTransform.txx

   1 /*=========================================================================
   2   Program:   vv                     http://www.creatis.insa-lyon.fr/rio/vv
   3
   4   Authors belong to:
   5   - University of LYON              http://www.universite-lyon.fr/
   6   - Léon Bérard cancer center       http://www.centreleonberard.fr
   7   - CREATIS CNRS laboratory         http://www.creatis.insa-lyon.fr
   8
   9   This software is distributed WITHOUT ANY WARRANTY; without even
  10   the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  11   PURPOSE.  See the copyright notices for more information.
  12
  13   It is distributed under dual licence
  14
  15   - BSD        See included LICENSE.txt file
  16   - CeCILL-B   http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
  17 ===========================================================================**/
  18 #ifndef __clitkShapedBLUTSpatioTemporalDeformableTransform_txx
  19 #define __clitkShapedBLUTSpatioTemporalDeformableTransform_txx
  20 #include "clitkShapedBLUTSpatioTemporalDeformableTransform.h"
  21
  22 //itk
  23 #include "itkContinuousIndex.h"
  24 #include "itkImageRegionIterator.h"
  25 #include "itkImageRegionConstIterator.h"
  26 #include "itkImageRegionConstIteratorWithIndex.h"
  27 #include "itkIdentityTransform.h"
  28
  29 namespace clitk
  30 {
  31
  32   // Constructor with default arguments
  33   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
  34   ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
  35   ::ShapedBLUTSpatioTemporalDeformableTransform():Superclass(OutputDimension,0)
  36   {
  37     unsigned int i;
  38
  39     //JV default spline order
  40     for ( i=0;i<InputDimension; i++)
  41       m_SplineOrders[i]=3;
  42
  43     //JV default LUTSamplingfactor
  44     for ( i=0;i<InputDimension; i++)
  45       m_LUTSamplingFactors[i]=20;
  46
  47     for (i=0;i<InputDimension; i++)
  48       m_SupportSize[i] = m_SplineOrders[i]+1;
  49
  50     //Instantiate a Vector Bspline Interpolator
  51     m_VectorInterpolator =VectorInterpolatorType::New();
  52     m_VectorInterpolator->SetLUTSamplingFactors(m_LUTSamplingFactors);
  53     m_VectorInterpolator->SetSplineOrders(m_SplineOrders);
  54
  55     // Set Bulk transform to NULL (no bulk performed)
  56     m_BulkTransform = NULL;
  57
  58     // Mask
  59     m_Mask=NULL;
  60
  61     // Shape
  62     m_TransformShape=0;
  63
  64     // Default grid size is zero
  65     m_NullSize.Fill(0);
  66     m_NullIndex.Fill(0);
  67
  68     //JV region containing the parameters
  69     m_GridRegion.SetSize( m_NullSize);
  70     m_GridRegion.SetIndex( m_NullIndex );
  71
  72     //JV use second region over the images
  73     m_PaddedGridRegion.SetSize(m_NullSize);
  74     m_PaddedGridRegion.SetIndex(m_NullIndex);
  75
  76     //JV Maintain two origins
  77     m_GridOrigin.Fill( 0.0 );  // default origin is all zeros
  78     m_PaddedGridOrigin.Fill( 0.0 );  // default origin is all zeros
  79
  80     m_GridSpacing.Fill( 1.0 ); // default spacing is all ones
  81     m_GridDirection.SetIdentity(); // default spacing is all ones
  82
  83     m_InternalParametersBuffer = ParametersType(0);
  84     // Make sure the parameters pointer is not NULL after construction.
  85     m_InputParametersPointer = &m_InternalParametersBuffer;
  86
  87     // Initialize coeffient images
  88     m_WrappedImage = CoefficientImageType::New();
  89     m_WrappedImage->SetRegions( m_GridRegion );
  90     m_WrappedImage->SetOrigin( m_GridOrigin.GetDataPointer() );
  91     m_WrappedImage->SetSpacing( m_GridSpacing.GetDataPointer() );
  92     m_WrappedImage->SetDirection( m_GridDirection );
  93
  94     // JV Initialize the padded version
  95     m_PaddedCoefficientImage = CoefficientImageType::New();
  96     m_PaddedCoefficientImage->SetRegions( m_PaddedGridRegion );
  97     m_PaddedCoefficientImage->SetOrigin( m_GridOrigin.GetDataPointer() );
  98     m_PaddedCoefficientImage->SetSpacing( m_GridSpacing.GetDataPointer() );
  99     m_PaddedCoefficientImage->SetDirection( m_GridDirection );
 100
 101     m_CoefficientImage = NULL;
 102
 103     // Variables for computing interpolation
 104     for (i=0; i <InputDimension;i++)
 105       {
 106         m_Offset[i] = m_SplineOrders[i] / 2;
 107         if ( m_SplineOrders[i] % 2 )
 108           {
 109             m_SplineOrderOdd[i] = true;
 110           }
 111         else
 112           {
 113             m_SplineOrderOdd[i] = false;
 114           }
 115       }
 116
 117     //JV padded should be used when checking
 118     m_ValidRegion = m_PaddedGridRegion;
 119
 120     // Initialize jacobian images
 121     for (i=0; i <OutputDimension;i++)
 122       {
 123         m_JacobianImage[i] = JacobianImageType::New();
 124         m_JacobianImage[i]->SetRegions( m_GridRegion );
 125         m_JacobianImage[i]->SetOrigin( m_GridOrigin.GetDataPointer() );
 126         m_JacobianImage[i]->SetSpacing( m_GridSpacing.GetDataPointer() );
 127         m_JacobianImage[i]->SetDirection( m_GridDirection );
 128       }
 129
 130     /** Fixed Parameters store the following information:
 131      *     Grid Size
 132      *     Grid Origin
 133      *     Grid Spacing
 134      *     Grid Direction */
 135     //JV we add the splineOrders, LUTsamplingfactor, m_Mask  and m_BulkTransform
 136     /*
 137       Spline orders
 138       Sampling factors
 139       m_Mask
 140       m_BulkTransform
 141       The size of these is equal to the  NInputDimensions
 142     *********************************************************/
 143     this->m_FixedParameters.SetSize ( NInputDimensions * (NInputDimensions + 5)+3 );
 144     this->m_FixedParameters.Fill ( 0.0 );
 145     for ( i=0; i<NInputDimensions; i++)
 146       {
 147         this->m_FixedParameters[2*NInputDimensions+i] = m_GridSpacing[i];
 148       }
 149     for (unsigned int di=0; di<NInputDimensions; di++)
 150       {
 151         for (unsigned int dj=0; dj<NInputDimensions; dj++)
 152           {
 153             this->m_FixedParameters[3*NInputDimensions+(di*NInputDimensions+dj)] = m_GridDirection[di][dj];
 154           }
 155       }
 156
 157     //JV add splineOrders
 158     for ( i=0; i<NInputDimensions; i++)
 159       {
 160         this->m_FixedParameters[ ( (3+NInputDimensions)*NInputDimensions)+i] = (this->GetSplineOrders())[i];
 161       }
 162
 163     //JV add LUTsamplingFactors
 164     for ( i=0; i<NInputDimensions; i++)
 165       {
 166         this->m_FixedParameters[( (4+NInputDimensions)*NInputDimensions)+i ] = this->GetLUTSamplingFactors()[i];
 167       }
 168
 169     // JV add the mask pointer
 170     this->m_FixedParameters[( (5+NInputDimensions)*NInputDimensions)]=(double)((size_t)m_Mask);
 171
 172     // JV add the bulkTransform pointer
 173     this->m_FixedParameters[( (5+NInputDimensions)*NInputDimensions) +1]=(double)((size_t)m_BulkTransform);
 174
 175     // JV add the Transform shape
 176     this->m_FixedParameters[( (5+NInputDimensions)*NInputDimensions) +2]=(double)(m_TransformShape);
 177
 178
 179     // Calculate the PointToIndex matrices
 180     DirectionType scale;
 181     for( unsigned int i=0; i<OutputDimension; i++)
 182       {
 183         scale[i][i] = m_GridSpacing[i];
 184       }
 185
 186     m_IndexToPoint = m_GridDirection * scale;
 187     m_PointToIndex = m_IndexToPoint.GetInverse();
 188
 189     m_LastJacobianIndex = m_ValidRegion.GetIndex();
 190
 191
 192     //Weights to be used for the BC
 193     m_Weights.resize(4);
 194     std::vector<double> weights(5);
 195     for (i=0; i<4; i++)
 196       weights[0]=0;
 197     // NN
 198     weights[0]=1;
 199     weights[1]=1./2;
 200     weights[2]=0;
 201     weights[3]=0;
 202     weights[4]=0;
 203     m_Weights[0]=weights;
 204
 205     // Linear
 206     weights[0]=1;
 207     weights[1]=1./2;
 208     weights[2]=0;
 209     weights[3]=0;
 210     weights[4]=0;
 211     m_Weights[1]=weights;
 212
 213     // quadratic
 214     weights[0]=3./4.;
 215     weights[1]=1./2.;
 216     weights[2]=1./8.;
 217     weights[3]=0;
 218     weights[4]=0;
 219     m_Weights[2]=weights;
 220
 221     // cubic
 222     weights[0]=2./3.;
 223     weights[1]=23./48. ;
 224     weights[2]=1./6.;
 225     weights[3]=1./48.;
 226     weights[4]=0;
 227     m_Weights[3]=weights;
 228
 229     // Update the WeightRatios
 230     m_WeightRatio.resize(4);
 231     for (unsigned int i=0; i<4; i++)
 232       {
 233         m_WeightRatio[i].resize(4);
 234         for (unsigned int j=0; j<4; j++)
 235           m_WeightRatio[i][j]=m_Weights[m_SplineOrders[InputDimension-1]][i]/ m_Weights[m_SplineOrders[InputDimension-1]][j];
 236       }
 237
 238     //JV initialize some variables for jacobian calculation
 239     m_SupportRegion.SetSize(m_SupportSize);
 240     m_SupportIndex.Fill(0);
 241     m_SupportRegion.SetIndex(m_SupportIndex);
 242     for (  i = 0; i < SpaceDimension ; i++ )
 243       m_ZeroVector[i]=itk::NumericTraits<JacobianValueType>::Zero;
 244
 245     m_InitialOffset=0;
 246     m_FirstRegion.SetSize(m_NullSize);
 247     m_FirstRegion.SetIndex(m_NullIndex);
 248     m_SecondRegion.SetSize(m_NullSize);
 249     m_SecondRegion.SetIndex(m_NullIndex);
 250     m_ThirdRegion.SetSize(m_NullSize);
 251     m_ThirdRegion.SetIndex(m_NullIndex);
 252
 253     m_BCValues.resize(0);
 254     m_BCRegions.resize(0);
 255     m_BCSize=0;
 256     m_BC2Values.resize(0);
 257     m_BC2Regions.resize(0);
 258     m_BC2Size=0;
 259     m_BC3Values.resize(0);
 260     m_BC3Regions.resize(0);
 261     m_BC3Size=0;
 262
 263
 264     for (  i = 0; i < SpaceDimension ; i++ )
 265       {
 266         m_FirstIterator[i]= IteratorType( m_JacobianImage[i], m_FirstRegion);
 267         m_SecondIterator[i]= IteratorType( m_JacobianImage[i], m_SecondRegion);
 268         m_ThirdIterator[i]= IteratorType( m_JacobianImage[i], m_ThirdRegion);
 269         m_BCIterators[i].resize(0);
 270         m_BC2Iterators[i].resize(0);
 271         m_BC3Iterators[i].resize(0);
 272       }
 273
 274     this->Modified();
 275
 276   }
 277
 278
 279   // Destructor
 280   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
 281   ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
 282   ::~ShapedBLUTSpatioTemporalDeformableTransform()
 283   {
 284
 285   }
 286
 287
 288   // JV set Spline Order
 289   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
 290   void
 291   ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
 292   ::SetSplineOrder(const unsigned int & splineOrder)
 293   {
 294     SizeType splineOrders;
 295     for (unsigned int i=0;i<InputDimension;i++)splineOrders[i]=splineOrder;
 296
 297     this->SetSplineOrders(splineOrders);
 298   }
 299
 300
 301   // JV set Spline Orders
 302   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
 303   void
 304   ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
 305   ::SetSplineOrders(const SizeType & splineOrders)
 306   {
 307     if(m_SplineOrders!=splineOrders)
 308       {
 309         m_SplineOrders=splineOrders;
 310
 311         //update the interpolation function
 312         m_VectorInterpolator->SetSplineOrders(m_SplineOrders);
 313
 314         //update the varaibles for computing interpolation
 315         for (unsigned int i=0; i <InputDimension;i++)
 316           {
 317             m_SupportSize[i] = m_SplineOrders[i]+1;
 318             m_Offset[i] = m_SplineOrders[i] / 2;
 319
 320             if ( m_SplineOrders[i] % 2 )
 321               {
 322                 m_SplineOrderOdd[i] = true;
 323               }
 324             else
 325               {
 326                 m_SplineOrderOdd[i] = false;
 327               }
 328           }
 329
 330         //SupportSize is updated!, update regions
 331         //JV initialize some variables for jacobian calculation
 332         m_SupportRegion.SetSize(m_SupportSize);
 333         m_SupportIndex.Fill(0);
 334         m_SupportRegion.SetIndex(m_SupportIndex);
 335
 336         // Update the WeightRatios
 337         for (unsigned int i=0; i<4; i++)
 338           for (unsigned int j=0; j<4; j++)
 339             m_WeightRatio[i][j]=m_Weights[m_SplineOrders[InputDimension-1]][i]/ m_Weights[m_SplineOrders[InputDimension-1]][j];
 340
 341         this->Modified();
 342       }
 343   }
 344
 345
 346   // JV set sampling factor
 347   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
 348   void
 349   ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
 350   ::SetLUTSamplingFactor( const int & samplingFactor)
 351   {
 352     SizeType samplingFactors;
 353     for (unsigned int i=0; i<NInputDimensions; i++)
 354       samplingFactors[i]=samplingFactor;
 355
 356     //update the interpolation function
 357     this->SetLUTSamplingFactors(samplingFactors);
 358   }
 359
 360
 361   // JV set sampling factors
 362   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
 363   void
 364   ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
 365   ::SetLUTSamplingFactors( const SizeType & samplingFactors)
 366   {
 367     if(m_LUTSamplingFactors!=samplingFactors)
 368       {
 369         for (unsigned int i=0; i<NInputDimensions; i++)
 370           m_LUTSamplingFactors[i]=samplingFactors[i];
 371
 372         //update the interpolation function
 373         m_VectorInterpolator->SetLUTSamplingFactors(m_LUTSamplingFactors);
 374
 375         this->Modified();
 376       }
 377   }
 378
 379
 380   // Get the number of parameters
 381   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
 382   unsigned int
 383   ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
 384   ::GetNumberOfParameters(void) const
 385   {
 386
 387     // The number of parameters equal SpaceDimension * number of
 388     // of pixels in the grid region.
 389     return ( static_cast<unsigned int>( SpaceDimension ) *
 390              static_cast<unsigned int>( m_GridRegion.GetNumberOfPixels() ) );
 391
 392   }
 393
 394
 395   // Get the padded number of parameters
 396   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
 397   unsigned int
 398   ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
 399   ::GetPaddedNumberOfParameters(void) const
 400   {
 401
 402     // The number of parameters equal SpaceDimension * number of
 403     // of pixels in the grid region.
 404     return ( static_cast<unsigned int>( SpaceDimension ) *
 405              static_cast<unsigned int>( m_PaddedGridRegion.GetNumberOfPixels() ) );
 406
 407   }
 408
 409
 410
 411   // Get the number of parameters per dimension
 412   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
 413   unsigned int
 414   ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
 415   ::GetNumberOfParametersPerDimension(void) const
 416   {
 417     // The number of parameters per dimension equal number of
 418     // of pixels in the grid region.
 419     return ( static_cast<unsigned int>( m_GridRegion.GetNumberOfPixels() ) );
 420
 421   }
 422
 423
 424   // Set the grid region
 425   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
 426   void
 427   ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
 428   ::SetGridRegion( const RegionType & region )
 429   {
 430     if ( m_GridRegion != region )
 431       {
 432         m_GridRegion = region;
 433         m_PaddedGridRegion=region;
 434
 435         // JV set the padded region
 436         typename CoefficientImageType::RegionType::SizeType paddedSize= region.GetSize();
 437
 438         //JV   size dependes on shape
 439         switch (m_TransformShape)
 440           {
 441           case 0:
 442           case 1:
 443             paddedSize[InputDimension-1]+=4;
 444             break;
 445           case 2:
 446           case 3:
 447             paddedSize[InputDimension-1]+=4;
 448             break;
 449           case 4:
 450           case 5:
 451             paddedSize[InputDimension-1]+=3;
 452             break;
 453           case 6:
 454           case 7:
 455             paddedSize[InputDimension-1]+=2;
 456             break;
 457           case 8:
 458           case 9:
 459             paddedSize[InputDimension-1]+=3;
 460             break;
 461           default:
 462             paddedSize[InputDimension-1]=+1;
 463           }
 464         m_PaddedGridRegion.SetSize(paddedSize);
 465
 466         // Set regions for each coefficient and jacobian image
 467         m_WrappedImage->SetRegions( m_GridRegion );
 468         m_PaddedCoefficientImage->SetRegions( m_PaddedGridRegion );
 469         m_PaddedCoefficientImage->Allocate();
 470         for (unsigned int j=0; j <OutputDimension;j++)
 471           {
 472             m_JacobianImage[j]->SetRegions( m_GridRegion );
 473           }
 474
 475         // JV used the padded version for the valid region
 476         // Set the valid region
 477         // If the grid spans the interval [start, last].
 478         // The valid interval for evaluation is [start+offset, last-offset]
 479         // when spline order is even.
 480         // The valid interval for evaluation is [start+offset, last-offset)
 481         // when spline order is odd.
 482         // Where offset = vcl_floor(spline / 2 ).
 483         // Note that the last pixel is not included in the valid region
 484         // with odd spline orders.
 485         typename RegionType::SizeType size = m_PaddedGridRegion.GetSize();
 486         typename RegionType::IndexType index = m_PaddedGridRegion.GetIndex();
 487         for ( unsigned int j = 0; j < NInputDimensions; j++ )
 488           {
 489             index[j] +=
 490               static_cast< typename RegionType::IndexValueType >( m_Offset[j] );
 491             size[j] -=
 492               static_cast< typename RegionType::SizeValueType> ( 2 * m_Offset[j] );
 493             m_ValidRegionLast[j] = index[j] +
 494               static_cast< typename RegionType::IndexValueType >( size[j] ) - 1;
 495           }
 496         m_ValidRegion.SetSize( size );
 497         m_ValidRegion.SetIndex( index );
 498
 499         // If we are using the default parameters, update their size and set to identity.
 500         // Input parameters point to internal buffer => using default parameters.
 501         if (m_InputParametersPointer == &m_InternalParametersBuffer)
 502           {
 503             // Check if we need to resize the default parameter buffer.
 504             if ( m_InternalParametersBuffer.GetSize() != this->GetNumberOfParameters() )
 505               {
 506                 m_InternalParametersBuffer.SetSize( this->GetNumberOfParameters() );
 507                 // Fill with zeros for identity.
 508                 m_InternalParametersBuffer.Fill( 0 );
 509               }
 510           }
 511
 512         this->Modified();
 513       }
 514   }
 515
 516
 517   // Set the grid spacing
 518   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
 519   void
 520   ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
 521   ::SetGridSpacing( const SpacingType & spacing )
 522   {
 523     if ( m_GridSpacing != spacing )
 524       {
 525         m_GridSpacing = spacing;
 526
 527         // Set spacing for each coefficient and jacobian image
 528         m_WrappedImage->SetSpacing( m_GridSpacing.GetDataPointer() );
 529         m_PaddedCoefficientImage->SetSpacing( m_GridSpacing.GetDataPointer() );
 530         for (unsigned int j=0; j <OutputDimension; j++) m_JacobianImage[j]->SetSpacing( m_GridSpacing.GetDataPointer() );
 531
 532         // Set scale
 533         DirectionType scale;
 534         for( unsigned int i=0; i<OutputDimension; i++)
 535           {
 536             scale[i][i] = m_GridSpacing[i];
 537           }
 538
 539         m_IndexToPoint = m_GridDirection * scale;
 540         m_PointToIndex = m_IndexToPoint.GetInverse();
 541
 542         this->Modified();
 543       }
 544
 545   }
 546
 547   // Set the grid direction
 548   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
 549   void
 550   ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
 551   ::SetGridDirection( const DirectionType & direction )
 552   {
 553     if ( m_GridDirection != direction )
 554       {
 555         m_GridDirection = direction;
 556
 557         // Set direction for each coefficient and jacobian image
 558         m_WrappedImage->SetDirection( m_GridDirection );
 559         m_PaddedCoefficientImage->SetDirection( m_GridDirection );
 560         for (unsigned int j=0; j <OutputDimension; j++) m_JacobianImage[j]->SetDirection( m_GridDirection );
 561
 562         // Set scale
 563         DirectionType scale;
 564         for( unsigned int i=0; i<OutputDimension; i++)
 565           {
 566             scale[i][i] = m_GridSpacing[i];
 567           }
 568
 569         m_IndexToPoint = m_GridDirection * scale;
 570         m_PointToIndex = m_IndexToPoint.GetInverse();
 571
 572         this->Modified();
 573       }
 574
 575   }
 576
 577
 578   // Set the grid origin
 579   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
 580   void
 581   ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
 582   ::SetGridOrigin( const OriginType& origin )
 583   {
 584     if( m_GridOrigin!=origin)
 585       {
 586         m_GridOrigin = origin;
 587
 588         // JV The origin depends on the shape
 589         switch (m_TransformShape)
 590           {
 591           case 0:
 592           case 1:
 593             m_PaddedGridOrigin=origin;
 594             m_PaddedGridOrigin[InputDimension-1]=origin[InputDimension-1]-2* m_GridSpacing[InputDimension-1];
 595             break;
 596           case 2:
 597           case 3:
 598             m_PaddedGridOrigin=origin;
 599             m_PaddedGridOrigin[InputDimension-1]=origin[InputDimension-1]-1* m_GridSpacing[InputDimension-1];
 600             break;
 601           case 4:
 602           case 5:
 603             m_PaddedGridOrigin=origin;
 604             m_PaddedGridOrigin[InputDimension-1]=origin[InputDimension-1]-1* m_GridSpacing[InputDimension-1];
 605             break;
 606           case 6:
 607           case 7:
 608             m_PaddedGridOrigin=origin;
 609             break;
 610           case 8:
 611           case 9:
 612             m_PaddedGridOrigin=origin;
 613             m_PaddedGridOrigin[InputDimension-1]=origin[InputDimension-1]-1* m_GridSpacing[InputDimension-1];
 614             break;
 615           default:
 616             m_PaddedGridOrigin=origin;
 617           }
 618
 619         // Set origin for each coefficient and jacobianimage
 620         m_WrappedImage->SetOrigin( m_GridOrigin.GetDataPointer() );
 621         m_PaddedCoefficientImage->SetOrigin( m_PaddedGridOrigin.GetDataPointer() );
 622         for (unsigned int j=0; j <OutputDimension; j++) m_JacobianImage[j]->SetOrigin( m_GridOrigin.GetDataPointer() );
 623
 624         this->Modified();
 625       }
 626   }
 627
 628
 629   // Set the parameters
 630   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
 631   void
 632   ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
 633   ::SetIdentity()
 634   {
 635     if( m_InputParametersPointer )
 636       {
 637         ParametersType * parameters =
 638           const_cast<ParametersType *>( m_InputParametersPointer );
 639         parameters->Fill( 0.0 );
 640         this->Modified();
 641       }
 642     else
 643       {
 644         itkExceptionMacro( << "Input parameters for the spline haven't been set ! "
 645                            << "Set them using the SetParameters or SetCoefficientImage method first." );
 646       }
 647   }
 648
 649
 650   // Set the parameters
 651   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
 652   void
 653   ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
 654   ::SetParameters( const ParametersType & parameters )
 655   {
 656
 657     // Check if the number of parameters match the
 658     // Expected number of parameters
 659     if ( parameters.Size() != this->GetNumberOfParameters() )
 660       {
 661         itkExceptionMacro(<<"Mismatched between parameters size "
 662                           << parameters.size()
 663                           << " and region size "
 664                           << m_GridRegion.GetNumberOfPixels() );
 665       }
 666
 667     // Clean up buffered parameters
 668     m_InternalParametersBuffer = ParametersType( 0 );
 669
 670     // Keep a reference to the input parameters
 671     m_InputParametersPointer = &parameters;
 672
 673     // Wrap flat array as images of coefficients
 674     this->WrapAsImages();
 675
 676     //JV Set padded input to vector interpolator
 677     m_VectorInterpolator->SetInputImage(this->GetPaddedCoefficientImage());
 678
 679     // Modified is always called since we just have a pointer to the
 680     // parameters and cannot know if the parameters have changed.
 681     this->Modified();
 682   }
 683
 684
 685   // Set the Fixed Parameters
 686   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
 687   void
 688   ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
 689   ::SetFixedParameters( const ParametersType & parameters )
 690   {
 691
 692     // JV number should be exact, no defaults for spacing
 693     if ( parameters.Size() != NInputDimensions * (5 + NInputDimensions)+3 )
 694       {
 695         itkExceptionMacro(<< "Mismatched between parameters size "
 696                           << parameters.size()
 697                           << " and number of fixed parameters "
 698                           << NInputDimensions * (5 + NInputDimensions)+3 );
 699       }
 700     /*********************************************************
 701     Fixed Parameters store the following information:
 702         Grid Size
 703         Grid Origin
 704         Grid Spacing
 705         Grid Direction */
 706     // JV we add the splineOrders, LUTsamplingfactor, mask pointer and bulktransform pointer
 707     /*
 708       Spline orders
 709       Sampling factors
 710       m_Mask
 711       m_BulkTransform
 712       m_TransformShape
 713
 714       The size of these is equal to the  NInputDimensions
 715     *********************************************************/
 716
 717     /** Set the Grid Parameters */
 718     SizeType   gridSize;
 719     for (unsigned int i=0; i<NInputDimensions; i++)
 720       {
 721         gridSize[i] = static_cast<int> (parameters[i]);
 722       }
 723     RegionType bsplineRegion;
 724     bsplineRegion.SetSize( gridSize );
 725
 726     /** Set the Origin Parameters */
 727     OriginType origin;
 728     for (unsigned int i=0; i<NInputDimensions; i++)
 729       {
 730         origin[i] = parameters[NInputDimensions+i];
 731       }
 732
 733     /** Set the Spacing Parameters */
 734     SpacingType spacing;
 735     for (unsigned int i=0; i<NInputDimensions; i++)
 736       {
 737         spacing[i] = parameters[2*NInputDimensions+i];
 738       }
 739
 740     /** Set the Direction Parameters */
 741     DirectionType direction;
 742     for (unsigned int di=0; di<NInputDimensions; di++)
 743       {
 744         for (unsigned int dj=0; dj<NInputDimensions; dj++)
 745           {
 746             direction[di][dj] = parameters[3*NInputDimensions+(di*NInputDimensions+dj)];
 747           }
 748       }
 749
 750     //JV add the spline orders
 751     SizeType splineOrders;
 752     for (unsigned int i=0; i<NInputDimensions; i++)
 753       {
 754         splineOrders[i]=parameters[(3+NInputDimensions)*NInputDimensions+i];
 755       }
 756
 757     //JV add the LUT sampling factor
 758     SizeType  samplingFactors;
 759     for (unsigned int i=0; i<NInputDimensions; i++)
 760       {
 761         samplingFactors[i]=parameters[(4+NInputDimensions)*NInputDimensions+i];
 762       }
 763
 764     //JV add the MaskPointer
 765     m_Mask=((MaskPointer)( (size_t)parameters[(5+NInputDimensions)*NInputDimensions]));
 766
 767     //JV add the MaskPointer
 768     m_BulkTransform=((BulkTransformPointer)( (size_t)parameters[(5+NInputDimensions)*NInputDimensions+1]));
 769
 770     //JV add the TransformShape
 771     m_TransformShape=((unsigned int)parameters[(5+NInputDimensions)*NInputDimensions+2]);
 772
 773     // Set the members
 774     this->SetSplineOrders( splineOrders );
 775     this->SetGridSpacing( spacing );
 776     this->SetGridDirection( direction );
 777     this->SetGridOrigin( origin );
 778     this->SetGridRegion( bsplineRegion );
 779     this->SetLUTSamplingFactors( samplingFactors );
 780
 781   }
 782
 783
 784   // Wrap flat parameters as images
 785   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
 786   void
 787   ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
 788   ::WrapAsImages()
 789   {
 790     //JV Wrap parameter array in vectorial image, changed parameter order: A1x A1y A1z, A2x ....
 791     PixelType * dataPointer =reinterpret_cast<PixelType *>( const_cast<double *>(m_InputParametersPointer->data_block() )) ;
 792     unsigned int numberOfPixels = m_GridRegion.GetNumberOfPixels();
 793
 794     m_WrappedImage->GetPixelContainer()->SetImportPointer( dataPointer,numberOfPixels);//InputDimension
 795     m_CoefficientImage = m_WrappedImage;
 796
 797     //=====================================
 798     //JV Create padded structure adding BC
 799     //=====================================
 800     PadCoefficientImage();
 801
 802     //=====================================
 803     //JV Wrap jacobian into OutputDimension X Vectorial images
 804     //=====================================
 805     this->m_Jacobian.set_size( OutputDimension, this->GetNumberOfParameters() );
 806
 807     // Use memset to set the memory
 808     // JV four rows of three comps of parameters
 809     JacobianPixelType * jacobianDataPointer = reinterpret_cast<JacobianPixelType *>(this->m_Jacobian.data_block());
 810     memset(jacobianDataPointer, 0,  OutputDimension*numberOfPixels*sizeof(JacobianPixelType));
 811
 812     for (unsigned int j=0; j<OutputDimension; j++)
 813       {
 814         m_JacobianImage[j]->GetPixelContainer()->
 815           SetImportPointer( jacobianDataPointer, numberOfPixels );
 816         jacobianDataPointer += numberOfPixels;
 817       }
 818
 819     // Reset the J parameters
 820     m_LastJacobianIndex = m_ValidRegion.GetIndex();
 821     m_FirstRegion.SetSize(m_NullSize);
 822     m_SecondRegion.SetSize(m_NullSize);
 823     for ( unsigned int j = 0; j < SpaceDimension ; j++ )
 824       {
 825         m_FirstIterator[j]= IteratorType( m_JacobianImage[j], m_FirstRegion);
 826         m_SecondIterator[j]= IteratorType( m_JacobianImage[j], m_SecondRegion);
 827       }
 828
 829     m_BCValues.resize(0);
 830     m_BCRegions.resize(0);
 831     m_BCSize=0;
 832     m_BC2Values.resize(0);
 833     m_BC2Regions.resize(0);
 834     m_BC2Size=0;
 835     m_BC3Values.resize(0);
 836     m_BC3Regions.resize(0);
 837     m_BC3Size=0;
 838
 839   }
 840
 841
 842   // Set the parameters by value
 843   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
 844   void
 845   ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
 846   ::SetParametersByValue( const ParametersType & parameters )
 847   {
 848
 849     // Check if the number of parameters match the
 850     // Expected number of parameters
 851     if ( parameters.Size() != this->GetNumberOfParameters() )
 852       {
 853         itkExceptionMacro(<<"Mismatched between parameters size "
 854                           << parameters.size()
 855                           << " and region size "
 856                           << m_GridRegion.GetNumberOfPixels() );
 857       }
 858
 859     // Copy it
 860     m_InternalParametersBuffer = parameters;
 861     m_InputParametersPointer = &m_InternalParametersBuffer;
 862
 863     // Wrap flat array as images of coefficients
 864     this->WrapAsImages();
 865
 866     //JV Set padded input to vector interpolator
 867     m_VectorInterpolator->SetInputImage(this->GetPaddedCoefficientImage());
 868
 869     // Modified is always called since we just have a pointer to the
 870     // Parameters and cannot know if the parameters have changed.
 871     this->Modified();
 872
 873   }
 874
 875   // Get the parameters
 876   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
 877   const
 878   typename ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
 879   ::ParametersType &
 880   ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
 881   ::GetParameters( void ) const
 882   {
 883     /** NOTE: For efficiency, this class does not keep a copy of the parameters -
 884      * it just keeps pointer to input parameters.
 885      */
 886     if (NULL == m_InputParametersPointer)
 887       {
 888         itkExceptionMacro( <<"Cannot GetParameters() because m_InputParametersPointer is NULL. Perhaps SetCoefficientImage() has been called causing the NULL pointer." );
 889       }
 890
 891     return (*m_InputParametersPointer);
 892   }
 893
 894
 895   // Get the parameters
 896   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
 897   const
 898   typename ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
 899   ::ParametersType &
 900   ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
 901   ::GetFixedParameters( void ) const
 902   {
 903     RegionType resRegion = this->GetGridRegion(  );
 904
 905     for (unsigned int i=0; i<NInputDimensions; i++)
 906       {
 907         this->m_FixedParameters[i] = (resRegion.GetSize())[i];
 908       }
 909     for (unsigned int i=0; i<NInputDimensions; i++)
 910       {
 911         this->m_FixedParameters[NInputDimensions+i] = (this->GetGridOrigin())[i];
 912       }
 913     for (unsigned int i=0; i<NInputDimensions; i++)
 914       {
 915         this->m_FixedParameters[2*NInputDimensions+i] =  (this->GetGridSpacing())[i];
 916       }
 917     for (unsigned int di=0; di<NInputDimensions; di++)
 918       {
 919         for (unsigned int dj=0; dj<NInputDimensions; dj++)
 920           {
 921             this->m_FixedParameters[3*NInputDimensions+(di*NInputDimensions+dj)] = (this->GetGridDirection())[di][dj];
 922           }
 923       }
 924
 925     //JV add splineOrders
 926     for (unsigned int i=0; i<NInputDimensions; i++)
 927       {
 928         this->m_FixedParameters[(3+NInputDimensions)*NInputDimensions+i] = (this->GetSplineOrders())[i];
 929       }
 930
 931     //JV add LUTsamplingFactor
 932     for (unsigned int i=0; i<NInputDimensions; i++)
 933       {
 934         this->m_FixedParameters[(4+NInputDimensions)*NInputDimensions+i] = (this->GetLUTSamplingFactors())[i];
 935       }
 936
 937     //JV add the mask
 938     this->m_FixedParameters[(5+NInputDimensions)*NInputDimensions]=(double)((size_t) m_Mask);
 939
 940     //JV add the bulktransform pointer
 941     this->m_FixedParameters[(5+NInputDimensions)*NInputDimensions+1]=(double)((size_t) m_BulkTransform);
 942
 943     //JV add the transform shape
 944     this->m_FixedParameters[(5+NInputDimensions)*NInputDimensions+2]=(double)( m_TransformShape);
 945
 946     return (this->m_FixedParameters);
 947   }
 948
 949
 950   // Set the B-Spline coefficients using input images
 951   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
 952   void
 953   ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
 954   ::SetCoefficientImage( CoefficientImagePointer image )
 955   {
 956     this->SetGridSpacing( image->GetSpacing() );
 957     this->SetGridOrigin( image->GetOrigin() );
 958     this->SetGridDirection( image->GetDirection() );
 959     this->SetGridRegion( image->GetBufferedRegion() );
 960     m_CoefficientImage = image;
 961
 962     //JV
 963     //   m_WrappedImage=m_CoefficientImage;
 964
 965     // Update the interpolator
 966     this->PadCoefficientImage();
 967     m_VectorInterpolator->SetInputImage(this->GetPaddedCoefficientImage());
 968
 969     // Clean up buffered parameters
 970     m_InternalParametersBuffer = ParametersType( 0 );
 971     m_InputParametersPointer  = NULL;
 972
 973   }
 974
 975
 976   //   // Set the B-Spline coefficients using input images
 977   //   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
 978   //   void
 979   //   ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
 980   //   ::SetPaddedCoefficientImage( CoefficientImagePointer image )
 981   //   {
 982   //     //JV modify the region
 983   //     typename CoefficientImageType::RegionType region=image->GetBufferedRegion();
 984   //     typename CoefficientImageType::RegionType::SizeType size=region.GetSize();
 985   //     size[InputDimension-1]-=2;
 986   //     region.SetSize(size);
 987
 988   //     //Set properties
 989   //     this->SetGridRegion( region );
 990   //     this->SetGridSpacing( image->GetSpacing() );
 991   //     this->SetGridDirection( image->GetDirection() );
 992   //     this->SetGridOrigin( image->GetOrigin() );
 993   //     m_PaddedCoefficientImage = image;
 994   //     this->ExtractCoefficientImage();
 995   //     m_VectorInterpolator->SetInputImage(this->GetPaddedCoefficientImage());
 996
 997   //     // Clean up buffered parameters
 998   //     m_InternalParametersBuffer = ParametersType( 0 );
 999   //     m_InputParametersPointer  = NULL;
1000
1001   //   }
1002
1003   // Set the B-Spline coefficients using input images
1004   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
1005   typename   ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>::CoefficientImageType::Pointer
1006   ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
1007   ::ExtractTemporalRow(const typename CoefficientImageType::Pointer& coefficientImage, unsigned int temporalIndex)
1008   {
1009     // Region
1010     typename   CoefficientImageType::RegionType sourceRegion=coefficientImage->GetLargestPossibleRegion();
1011     sourceRegion.SetSize(InputDimension-1, 1);
1012     sourceRegion.SetIndex(InputDimension-1, temporalIndex);
1013
1014     // Extract
1015     typedef clitk::ExtractImageFilter<CoefficientImageType, CoefficientImageType> ExtractImageFilterType;
1016     typename ExtractImageFilterType::Pointer extract=ExtractImageFilterType::New();
1017     extract->SetInput(coefficientImage);
1018     extract->SetExtractionRegion(sourceRegion);
1019     extract->Update();
1020     typename CoefficientImageType::Pointer row= extract->GetOutput();
1021
1022     // Set index to zero
1023     sourceRegion.SetIndex(InputDimension-1, 0);
1024     row->SetRegions(sourceRegion);
1025     return row;
1026
1027   }
1028
1029   // Set the B-Spline coefficients using input images
1030   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
1031   void
1032   ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
1033   ::PadCoefficientImage(void)
1034   {
1035
1036     // Define paste, extract and combine filters
1037     typedef itk::PasteImageFilter<CoefficientImageType, CoefficientImageType, CoefficientImageType> PasteImageFilterType;
1038     typedef clitk::ExtractImageFilter<CoefficientImageType, CoefficientImageType> ExtractImageFilterType;
1039     typedef clitk::LinearCombinationImageFilter<CoefficientImageType, CoefficientImageType> LinearCombinationFilterType;
1040
1041     // Regions
1042     typename   CoefficientImageType::RegionType sourceRegion=m_PaddedCoefficientImage->GetLargestPossibleRegion();
1043     typename   CoefficientImageType::RegionType destinationRegion=m_PaddedCoefficientImage->GetLargestPossibleRegion();
1044     typename   CoefficientImageType::RegionType::SizeType sourceSize=sourceRegion.GetSize();
1045     typename   CoefficientImageType::RegionType::SizeType destinationSize=destinationRegion.GetSize();
1046     typename   CoefficientImageType::IndexType sourceIndex=sourceRegion.GetIndex();
1047     typename   CoefficientImageType::IndexType destinationIndex=destinationRegion.GetIndex();
1048
1049     // JV Padding depends on the shape
1050     switch (m_TransformShape)
1051       {
1052         /*  The shapes are
1053             0: egg     4 CP  3 DOF
1054             1: egg     5 CP  4 DOF
1055             2: rabbit  4 CP  3 DOF
1056             3: rabbit  5 CP  4 DOF
1057             4: sputnik 4 CP  4 DOF
1058             5: sputnik 5 CP  5 DOF
1059             6: diamond 6 CP  5 DOF
1060             7: diamond 7 CP  6 DOF
1061         */
1062
1063       case 0:
1064         {
1065           // ----------------------------------------------------------------------
1066           // The egg with 4 internal CP (starting from inhale)
1067           // Periodic, constrained to zero at the reference
1068           // at position 3 and
1069           // Coeff      row  BC1 BC2  0  BC3  1   2  BC4
1070           // PaddedCoeff  R:  0   1   2   3   4   5   6
1071           // BC1= R4
1072           // BC2= R5
1073           // BC3= -weights[2]/weights[0] ( R2+R4 )
1074           // BC4= R2
1075           // ---------------------------------------------------------------------
1076
1077           //---------------------------------
1078           // 1. First two temporal rows are identical: paste 1-2 to 0-1
1079           typename PasteImageFilterType::Pointer paster0=PasteImageFilterType::New();
1080           paster0->SetDestinationImage(m_PaddedCoefficientImage);
1081           paster0->SetSourceImage(m_CoefficientImage);
1082           sourceRegion.SetIndex(NInputDimensions-1,1);
1083           sourceRegion.SetSize(NInputDimensions-1,2);
1084           destinationIndex[InputDimension-1]=0;
1085           paster0->SetDestinationIndex(destinationIndex);
1086           paster0->SetSourceRegion(sourceRegion);
1087
1088           //---------------------------------
1089           // 1. Next temporal row = paste 0 to 2
1090           typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1091           typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1092           paster1->SetDestinationImage(paster0->GetOutput());
1093           paster1->SetSourceImage(row0);
1094           destinationIndex[InputDimension-1]=2;
1095           paster1->SetDestinationIndex(destinationIndex);
1096           paster1->SetSourceRegion(row0->GetLargestPossibleRegion());
1097
1098           //---------------------------------
1099           // 2. Middle row at index 3=BC
1100           // Extract row at index 1, 2 (of coeff image)
1101           typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1102
1103           // Combine
1104           typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1105           combine1->SetInput(0,row0);
1106           combine1->SetInput(1,row1);
1107           combine1->SetA(-m_WeightRatio[2][0]);
1108           combine1->SetB(-m_WeightRatio[2][0]);
1109           combine1->Update();
1110           typename CoefficientImageType::Pointer bc3Row=combine1->GetOutput();
1111
1112           // Paste middleRow at index 3 (padded coeff)
1113           typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1114           paster2->SetDestinationImage(paster1->GetOutput());
1115           paster2->SetSourceImage(bc3Row);
1116           destinationIndex[InputDimension-1]=3;
1117           paster2->SetDestinationIndex(destinationIndex);
1118           paster2->SetSourceRegion(bc3Row->GetLargestPossibleRegion());
1119
1120           //---------------------------------
1121           // 3. Next two temporal rows identical: paste 1,2 to 4,5
1122           typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1123           paster3->SetDestinationImage(paster2->GetOutput());
1124           paster3->SetSourceImage(m_CoefficientImage);
1125           sourceRegion.SetIndex(InputDimension-1,1);
1126           sourceRegion.SetSize(InputDimension-1,2);
1127           destinationIndex[InputDimension-1]=4;
1128           paster3->SetDestinationIndex(destinationIndex);
1129           paster3->SetSourceRegion(sourceRegion);
1130
1131           // ---------------------------------
1132           // 4. Row at index 6=BC4= R2
1133           // Paste BC3 row at index 5
1134           typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
1135           paster4->SetDestinationImage(paster3->GetOutput());
1136           paster4->SetSourceImage(row0);
1137           destinationIndex[InputDimension-1]=6;
1138           paster4->SetDestinationIndex(destinationIndex);
1139           paster4->SetSourceRegion(row0->GetLargestPossibleRegion());
1140
1141           // Update the chain!
1142           paster4->Update();
1143           m_PaddedCoefficientImage= paster4->GetOutput();
1144
1145           break;
1146         }
1147
1148       case 1:
1149         {
1150           // ----------------------------------------------------------------------
1151           // The egg with 5 internal CP (starting from inhale)
1152           // Periodic, constrained to zero at the reference
1153           // at position 3 and
1154           // Coeff      row  BC1 BC2  0  BC3  1   2   3  BC4
1155           // PaddedCoeff  R:  0   1   2   3   4   5   6   7
1156           // BC1= R5
1157           // BC2= R6
1158           // BC3= -weights[3]/weights[1] ( R2+R5 ) - R4
1159           // BC4= R2
1160           // ---------------------------------------------------------------------
1161           //---------------------------------
1162           // 1. First two temporal rows are identical: paste 2-3 to 0-1
1163           typename PasteImageFilterType::Pointer paster0=PasteImageFilterType::New();
1164           paster0->SetDestinationImage(m_PaddedCoefficientImage);
1165           paster0->SetSourceImage(m_CoefficientImage);
1166           sourceRegion.SetIndex(NInputDimensions-1,2);
1167           sourceRegion.SetSize(NInputDimensions-1,2);
1168           destinationIndex[InputDimension-1]=0;
1169           paster0->SetDestinationIndex(destinationIndex);
1170           paster0->SetSourceRegion(sourceRegion);
1171
1172           //---------------------------------
1173           // 1. Next temporal row = paste 0 to 2
1174           typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1175           typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1176           paster1->SetDestinationImage(paster0->GetOutput());
1177           paster1->SetSourceImage(row0);
1178           destinationIndex[InputDimension-1]=2;
1179           paster1->SetDestinationIndex(destinationIndex);
1180           paster1->SetSourceRegion(row0->GetLargestPossibleRegion());
1181
1182           //---------------------------------
1183           // 2. Middle row at index 3=BC
1184           // Extract row at index 1, 2 (of coeff image)
1185           typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1186           typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
1187
1188           // Combine
1189           typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1190           combine1->SetInput(0,row0);
1191           combine1->SetInput(1,row2);
1192           combine1->SetA(1);
1193           combine1->SetB(1);
1194           // combine1->Update();
1195           typename LinearCombinationFilterType::Pointer combine2=LinearCombinationFilterType::New();
1196           combine2->SetInput(0,row1);
1197           combine2->SetInput(1,combine1->GetOutput());
1198           combine2->SetA(-1.);
1199           combine2->SetB(-m_WeightRatio[3][1]);
1200           combine2->Update();
1201           typename CoefficientImageType::Pointer bc3Row=combine2->GetOutput();
1202
1203           // Paste middleRow at index 3 (padded coeff)
1204           typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1205           paster2->SetDestinationImage(paster1->GetOutput());
1206           paster2->SetSourceImage(bc3Row);
1207           destinationIndex[InputDimension-1]=3;
1208           paster2->SetDestinationIndex(destinationIndex);
1209           paster2->SetSourceRegion(bc3Row->GetLargestPossibleRegion());
1210
1211           //---------------------------------
1212           // 3. Next three temporal rows identical: paste 1,2,3 to 4,5,6
1213           typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1214           paster3->SetDestinationImage(paster2->GetOutput());
1215           paster3->SetSourceImage(m_CoefficientImage);
1216           sourceRegion.SetIndex(InputDimension-1,1);
1217           sourceRegion.SetSize(InputDimension-1,3);
1218           destinationIndex[InputDimension-1]=4;
1219           paster3->SetDestinationIndex(destinationIndex);
1220           paster3->SetSourceRegion(sourceRegion);
1221
1222           // ---------------------------------
1223           // 4. Row at index 7=BC4= R2
1224           // Paste BC3 row at index 5
1225           typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
1226           paster4->SetDestinationImage(paster3->GetOutput());
1227           paster4->SetSourceImage(row0);
1228           destinationIndex[InputDimension-1]=7;
1229           paster4->SetDestinationIndex(destinationIndex);
1230           paster4->SetSourceRegion(row0->GetLargestPossibleRegion());
1231
1232           // Update the chain!
1233           paster4->Update();
1234           m_PaddedCoefficientImage= paster4->GetOutput();
1235
1236           break;
1237         }
1238
1239
1240 //        // ----------------------------------------------------------------------
1241 //        // The egg with 5 internal CP:
1242 //        // Periodic, C2 smooth everywhere and constrained to zero at the reference
1243 //        // Coeff       row R5 BC1  0   1   2   3  BC2  R2
1244 //        // PaddedCoeff R:  0   1   2   3   4   5   6   7
1245 //        // BC1= -weights[2]/weights[0] ( R2+R5)
1246 //        // BC2= BC1
1247 //        // ---------------------------------------------------------------------
1248
1249 //        // Extract rows with index 0 and 3
1250 //        typename    CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1251 //        typename    CoefficientImageType::Pointer row3=this->ExtractTemporalRow(m_CoefficientImage, 3);
1252
1253 //        // Paste the first row
1254 //        typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1255 //        destinationIndex.Fill(0);
1256 //        paster1->SetDestinationIndex(destinationIndex);
1257 //        paster1->SetSourceRegion(row3->GetLargestPossibleRegion());
1258 //        paster1->SetSourceImage(row3);
1259 //        paster1->SetDestinationImage(m_PaddedCoefficientImage);
1260
1261 //        // Linearly Combine rows for BC1 and BC2
1262 //        typedef clitk::LinearCombinationImageFilter<CoefficientImageType, CoefficientImageType> LinearCombinationFilterType;
1263 //        typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1264 //        combine1->SetFirstInput(row0);
1265 //        combine1->SetSecondInput(row3);
1266 //        combine1->SetA(-m_WeightRatio[2][0]);
1267 //        combine1->SetB(-m_WeightRatio[2][0]);
1268 //        combine1->Update();
1269 //        typename CoefficientImageType::Pointer bcRow=combine1->GetOutput();
1270
1271 //        // Paste the second row
1272 //        typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1273 //        destinationIndex[InputDimension-1]=1;
1274 //        paster2->SetDestinationIndex(destinationIndex);
1275 //        paster2->SetSourceRegion(bcRow->GetLargestPossibleRegion());
1276 //        paster2->SetSourceImage(bcRow);
1277 //        paster2->SetDestinationImage(paster1->GetOutput());
1278
1279 //        // Paste the coefficientImage
1280 //        typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1281 //        destinationIndex[InputDimension-1]=2;
1282 //        paster3->SetDestinationIndex(destinationIndex);
1283 //        paster3->SetSourceRegion(m_CoefficientImage->GetLargestPossibleRegion());
1284 //        paster3->SetSourceImage(m_CoefficientImage);
1285 //        paster3->SetDestinationImage(paster2->GetOutput());
1286
1287 //        // Paste the last two rows
1288 //        typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
1289 //        destinationIndex[InputDimension-1]=6;
1290 //        paster4->SetDestinationIndex(destinationIndex);
1291 //        paster4->SetSourceRegion(bcRow->GetLargestPossibleRegion());
1292 //        paster4->SetSourceImage(bcRow);
1293 //        paster4->SetDestinationImage(paster3->GetOutput());
1294
1295 //        typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
1296 //        destinationIndex[InputDimension-1]=7;
1297 //        paster5->SetDestinationIndex(destinationIndex);
1298 //        paster5->SetSourceRegion(row0->GetLargestPossibleRegion());
1299 //        paster5->SetSourceImage(row0);
1300 //        paster5->SetDestinationImage(paster4->GetOutput());
1301
1302 //        // Update the chain!
1303 //        paster5->Update();
1304 //        m_PaddedCoefficientImage= paster5->GetOutput();
1305
1306 //        break;
1307 //      }
1308
1309       case 2:
1310         {
1311           // ----------------------------------------------------------------------
1312           // The rabbit with 4 internal CP
1313           // Periodic, constrained to zero at the reference
1314           // at position 3 and the extremes fixed with anit-symm bc
1315           // Coeff      row  BC1  0   1  BC2  2  BC3 BC4
1316           // PaddedCoeff  R:  0   1   2   3   4   5   6
1317           // BC1= 2*R1-R2
1318           // BC2= -weights[2]/weights[0] ( R2+R4 )
1319           // BC3=  R1
1320           // BC4= 2*R1-R4
1321           // ---------------------------------------------------------------------
1322
1323           // ---------------------------------
1324           // 0. First Row =BC1
1325           typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1326           typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1327           typename LinearCombinationFilterType::Pointer combine0=LinearCombinationFilterType::New();
1328           combine0->SetInput(0,row0);
1329           combine0->SetInput(1,row1);
1330           combine0->SetA(2.);
1331           combine0->SetB(-1.);
1332           combine0->Update();
1333           typename CoefficientImageType::Pointer bc1Row=combine0->GetOutput();
1334
1335           typename PasteImageFilterType::Pointer paster0=PasteImageFilterType::New();
1336           paster0->SetDestinationImage(m_PaddedCoefficientImage);
1337           paster0->SetSourceImage(bc1Row);
1338           destinationIndex[InputDimension-1]=0;
1339           paster0->SetDestinationIndex(destinationIndex);
1340           paster0->SetSourceRegion(bc1Row->GetLargestPossibleRegion());
1341
1342           //---------------------------------
1343           // 1. Next two temporal rows are identical: paste 0-1 to 1-2
1344           typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1345           paster1->SetDestinationImage(paster0->GetOutput());
1346           paster1->SetSourceImage(m_CoefficientImage);
1347           sourceRegion.SetIndex(NInputDimensions-1,0);
1348           destinationIndex[InputDimension-1]=1;
1349           sourceRegion.SetSize(NInputDimensions-1,2);
1350           paster1->SetDestinationIndex(destinationIndex);
1351           paster1->SetSourceRegion(sourceRegion);
1352
1353           //---------------------------------
1354           // 2. Middle row at index 3=BC
1355           // Extract row at index 1, 2 (of coeff image)
1356           typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
1357
1358           // Combine
1359           typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1360           combine1->SetInput(0,row1);
1361           combine1->SetInput(1,row2);
1362           combine1->SetA(-m_WeightRatio[2][0]);
1363           combine1->SetB(-m_WeightRatio[2][0]);
1364           combine1->Update();
1365           typename CoefficientImageType::Pointer bc2Row=combine1->GetOutput();
1366
1367           // Paste middleRow at index 3 (padded coeff)
1368           typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1369           paster2->SetDestinationImage(paster1->GetOutput());
1370           paster2->SetSourceImage(bc2Row);
1371           destinationIndex[InputDimension-1]=3;
1372           paster2->SetDestinationIndex(destinationIndex);
1373           paster2->SetSourceRegion(bc2Row->GetLargestPossibleRegion());
1374
1375           //---------------------------------
1376           // 3. Next temporal row is identical: paste 2 to 4
1377           typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1378           paster3->SetDestinationImage(paster2->GetOutput());
1379           paster3->SetSourceImage(row2);
1380           destinationIndex[InputDimension-1]=4;
1381           paster3->SetDestinationIndex(destinationIndex);
1382           paster3->SetSourceRegion(row2->GetLargestPossibleRegion());
1383
1384           // ---------------------------------
1385           // 4. Row at index 5=BC (paddedcoeff image) R1
1386           // Paste BC3 row at index 5
1387           typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
1388           paster4->SetDestinationImage(paster3->GetOutput());
1389           paster4->SetSourceImage(row0);
1390           destinationIndex[InputDimension-1]=5;
1391           paster4->SetDestinationIndex(destinationIndex);
1392           paster4->SetSourceRegion(row0->GetLargestPossibleRegion());
1393
1394           // ---------------------------------
1395           // 5. Paste BC4 row at index 6
1396           typename LinearCombinationFilterType::Pointer combine3=LinearCombinationFilterType::New();
1397           combine3->SetInput(0,row0);
1398           combine3->SetInput(1,row2);
1399           combine3->SetA(2.);
1400           combine3->SetB(-1.);
1401           combine3->Update();
1402           typename CoefficientImageType::Pointer bc4Row=combine3->GetOutput();
1403           typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
1404           paster5->SetDestinationImage(paster4->GetOutput());
1405           paster5->SetSourceImage(bc4Row);
1406           destinationIndex[InputDimension-1]=6;
1407           paster5->SetDestinationIndex(destinationIndex);
1408           paster5->SetSourceRegion(bc4Row->GetLargestPossibleRegion());
1409
1410           // Update the chain!
1411           paster5->Update();
1412           m_PaddedCoefficientImage= paster5->GetOutput();
1413
1414           break;
1415         }
1416
1417       case 3:
1418         {
1419           // ----------------------------------------------------------------------
1420           // The rabbit with 5 internal CP
1421           // Periodic, constrained to zero at the reference at position 3.5
1422           // and the extremes fixed with anti-symmetrical BC
1423           // Coeff      row  BC1  0   1  BC2  2   3  BC3 BC4
1424           // PaddedCoeff  R:  0   1   2   3   4   5   6   7
1425           // BC1= 2*R1-R2
1426           // BC2= -weights[3]/weights[1] ( R2+R5 ) - R4
1427           // BC3= R1
1428           // BC4= 2*R1-R5
1429           // ---------------------------------------------------------------------
1430
1431           // ---------------------------------
1432           // 0. First Row =BC1
1433           typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1434           typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1435           typename LinearCombinationFilterType::Pointer combine0=LinearCombinationFilterType::New();
1436           combine0->SetInput(0,row0);
1437           combine0->SetInput(1,row1);
1438           combine0->SetA(2.);
1439           combine0->SetB(-1.);
1440           combine0->Update();
1441           typename CoefficientImageType::Pointer bc1Row=combine0->GetOutput();
1442
1443           typename PasteImageFilterType::Pointer paster0=PasteImageFilterType::New();
1444           paster0->SetDestinationImage(m_PaddedCoefficientImage);
1445           paster0->SetSourceImage(bc1Row);
1446           destinationIndex[InputDimension-1]=0;
1447           paster0->SetDestinationIndex(destinationIndex);
1448           paster0->SetSourceRegion(bc1Row->GetLargestPossibleRegion());
1449
1450           //---------------------------------
1451           // 1. Next two temporal rows are identical: paste 0-1 to 1-2
1452           typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1453           paster1->SetDestinationImage(paster0->GetOutput());
1454           paster1->SetSourceImage(m_CoefficientImage);
1455           sourceRegion.SetIndex(InputDimension-1,0);
1456           destinationIndex[InputDimension-1]=1;
1457           sourceRegion.SetSize(NInputDimensions-1,2);
1458           paster1->SetDestinationIndex(destinationIndex);
1459           paster1->SetSourceRegion(sourceRegion);
1460
1461           //---------------------------------
1462           // 2. Middle row at index 3=BC
1463           // Extract row at index 1, 3 (of coeff image)
1464           //typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1465           typename CoefficientImageType::Pointer row3=this->ExtractTemporalRow(m_CoefficientImage, 3);
1466
1467           // Combine
1468           typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1469           combine1->SetInput(0,row1);
1470           combine1->SetInput(1,row3);
1471           combine1->SetA(1.);
1472           combine1->SetB(1.);
1473           combine1->Update();
1474
1475           // Extract row at index 2 (coeff image)
1476           typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
1477
1478           // Combine
1479           typename LinearCombinationFilterType::Pointer combine2=LinearCombinationFilterType::New();
1480           combine2->SetInput(0,combine1->GetOutput());
1481           combine2->SetInput(1,row2);
1482           combine2->SetA(-m_WeightRatio[3][1]);
1483           combine2->SetB(-1.);
1484           combine2->Update();
1485           typename CoefficientImageType::Pointer bc2Row=combine2->GetOutput();
1486
1487           // Paste middleRow at index 3 (padded coeff)
1488           typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1489           paster2->SetDestinationImage(paster1->GetOutput());
1490           paster2->SetSourceImage(bc2Row);
1491           destinationIndex[InputDimension-1]=3;
1492           paster2->SetDestinationIndex(destinationIndex);
1493           paster2->SetSourceRegion(bc2Row->GetLargestPossibleRegion());
1494
1495           //---------------------------------
1496           // 3. Next two temporal rows are identical: paste 2,3 to 4,5
1497           typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1498           paster3->SetDestinationImage(paster2->GetOutput());
1499           paster3->SetSourceImage(m_CoefficientImage);
1500           sourceRegion.SetIndex(InputDimension-1,2);
1501           destinationIndex[InputDimension-1]=4;
1502           sourceRegion.SetSize(NInputDimensions-1,2);
1503           paster3->SetDestinationIndex(destinationIndex);
1504           paster3->SetSourceRegion(sourceRegion);
1505
1506           // ---------------------------------
1507           // 4. Row at index 6=BC (paddedcoeff image)R1
1508           // Paste BC3 row at index 6
1509           typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
1510           paster4->SetDestinationImage(paster3->GetOutput());
1511           paster4->SetSourceImage(row0);
1512           destinationIndex[InputDimension-1]=6;
1513           paster4->SetDestinationIndex(destinationIndex);
1514           paster4->SetSourceRegion(row0->GetLargestPossibleRegion());
1515
1516           // ---------------------------------
1517           // 5. Paste BC4 row at index 7
1518           typename LinearCombinationFilterType::Pointer combine3=LinearCombinationFilterType::New();
1519           combine3->SetInput(0,row0);
1520           combine3->SetInput(1,row3);
1521           combine3->SetA(2.);
1522           combine3->SetB(-1.);
1523           combine3->Update();
1524           typename CoefficientImageType::Pointer bc4Row=combine3->GetOutput();
1525           typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
1526           paster5->SetDestinationImage(paster4->GetOutput());
1527           paster5->SetSourceImage(bc4Row);
1528           destinationIndex[InputDimension-1]=7;
1529           paster5->SetDestinationIndex(destinationIndex);
1530           paster5->SetSourceRegion(bc4Row->GetLargestPossibleRegion());
1531
1532           // Update the chain!
1533           paster5->Update();
1534           m_PaddedCoefficientImage= paster5->GetOutput();
1535
1536           break;
1537         }
1538
1539       case 4:
1540         {
1541           // ----------------------------------------------------------------------
1542           // The sputnik with 4 internal CP
1543           // Periodic, constrained to zero at the reference
1544           // at position 3 and one indepent extremes copied
1545           // Coeff     row BC1  0   1  BC2  2  BC2  3
1546           // PaddedCoeff R: 0   1   2   3   4   5   6
1547           // BC1= R6
1548           // BC2= -weights[2]/weights[0] ( R2+R4 )
1549           // BC3=  weights[2]/weights[0] ( R2-R4) + R1
1550           // ---------------------------------------------------------------------
1551
1552           //---------------------------------
1553           // 1. First Row is equal to last row: paste 3 row to 0
1554           typename CoefficientImageType::Pointer row3=this->ExtractTemporalRow(m_CoefficientImage, 3);
1555           typename PasteImageFilterType::Pointer paster0=PasteImageFilterType::New();
1556           paster0->SetDestinationImage(m_PaddedCoefficientImage);
1557           paster0->SetSourceImage(row3);
1558           destinationIndex[InputDimension-1]=0;
1559           paster0->SetDestinationIndex(destinationIndex);
1560           paster0->SetSourceRegion(row3->GetLargestPossibleRegion());
1561
1562           //---------------------------------
1563           // 1. Next two temporal rows are identical: paste 0-1 to 1-2
1564           typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1565           paster1->SetDestinationImage(paster0->GetOutput());
1566           paster1->SetSourceImage(m_CoefficientImage);
1567           sourceRegion.SetIndex(NInputDimensions-1,0);
1568           destinationIndex[InputDimension-1]=1;
1569           sourceRegion.SetSize(NInputDimensions-1,2);
1570           paster1->SetDestinationIndex(destinationIndex);
1571           paster1->SetSourceRegion(sourceRegion);
1572
1573           //---------------------------------
1574           // 2. Middle row at index 3=BC
1575           // Extract row at index 1, 2 (of coeff image)
1576           typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1577           typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
1578
1579           // Combine
1580           typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1581           combine1->SetInput(0,row1);
1582           combine1->SetInput(1,row2);
1583           combine1->SetA(-m_WeightRatio[2][0]);
1584           combine1->SetB(-m_WeightRatio[2][0]);
1585           combine1->Update();
1586           typename CoefficientImageType::Pointer bc1Row=combine1->GetOutput();
1587
1588           // Paste middleRow at index 3 (padded coeff)
1589           typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1590           paster2->SetDestinationImage(paster1->GetOutput());
1591           paster2->SetSourceImage(bc1Row);
1592           destinationIndex[InputDimension-1]=3;
1593           paster2->SetDestinationIndex(destinationIndex);
1594           paster2->SetSourceRegion(bc1Row->GetLargestPossibleRegion());
1595
1596           //---------------------------------
1597           // 3. Next temporal row is identical: paste 2 to 4
1598           typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1599           paster3->SetDestinationImage(paster2->GetOutput());
1600           paster3->SetSourceImage(row2);
1601           destinationIndex[InputDimension-1]=4;
1602           paster3->SetDestinationIndex(destinationIndex);
1603           paster3->SetSourceRegion(row2->GetLargestPossibleRegion());
1604
1605           //---------------------------------
1606           // 4. Final row at index 5=BC3 (paddedcoeff image)R1 R2 R4 corresponds to index in coeff image 0 1 2
1607           // Extract row at index 1, 2 (of coeff image)already done
1608           typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1609
1610           // Combine
1611           typename LinearCombinationFilterType::Pointer combine3=LinearCombinationFilterType::New();
1612           combine3->SetInput(0,row1);
1613           combine3->SetInput(1,row2);
1614           combine3->SetA(1.);
1615           combine3->SetB(-1.);
1616
1617           // Combine
1618           typename LinearCombinationFilterType::Pointer combine4=LinearCombinationFilterType::New();
1619           combine4->SetInput(0,combine3->GetOutput());
1620           combine4->SetInput(1,row0);
1621           combine4->SetA(m_WeightRatio[2][0]);
1622           combine4->SetB(1.);
1623           combine4->Update();
1624           typename CoefficientImageType::Pointer bc2Row=combine4->GetOutput();
1625
1626           // Paste BC row at index 5
1627           typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
1628           paster4->SetDestinationImage(paster3->GetOutput());
1629           paster4->SetSourceImage(bc2Row);
1630           destinationIndex[InputDimension-1]=5;
1631           paster4->SetDestinationIndex(destinationIndex);
1632           paster4->SetSourceRegion(bc2Row->GetLargestPossibleRegion());
1633
1634           // Paste row 3  at index 6
1635           typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
1636           paster5->SetDestinationImage(paster4->GetOutput());
1637           paster5->SetSourceImage(row3);
1638           destinationIndex[InputDimension-1]=6;
1639           paster5->SetDestinationIndex(destinationIndex);
1640           paster5->SetSourceRegion(row3->GetLargestPossibleRegion());
1641
1642           // Update the chain!
1643           paster5->Update();
1644           m_PaddedCoefficientImage= paster5->GetOutput();
1645
1646           break;
1647         }
1648
1649       case 5:
1650         {
1651
1652           // ----------------------------------------------------------------------
1653           // The sputnik with 5 internal CP
1654           // Periodic, constrained to zero at the reference
1655           // at position 3 and one indepent extreme
1656           // Coeff     row BC1  0   1  BC2  2   3  BC3  4
1657           // PaddedCoeff R: 0   1   2   3   4   5   6   7
1658           // BC1= R2 + R5 - R7
1659           // BC2= -weights[3]/weights[1] ( R2+R5 ) - R4
1660           // BC3= R1 + 0.5 R2 - 0.5 R7
1661           // ----------------------------------------------------------------------
1662           //---------------------------------
1663           // 1. First Row =BC
1664           typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1665           typename CoefficientImageType::Pointer row3=this->ExtractTemporalRow(m_CoefficientImage, 3);
1666           typename CoefficientImageType::Pointer row4=this->ExtractTemporalRow(m_CoefficientImage, 4);
1667
1668           // Combine
1669           typename LinearCombinationFilterType::Pointer combine0=LinearCombinationFilterType::New();
1670           combine0->SetInput(0,row1);
1671           combine0->SetInput(1,row3);
1672           combine0->SetA(1.);
1673           combine0->SetB(1.);
1674           //combine0->Update();
1675           typename LinearCombinationFilterType::Pointer combine0bis=LinearCombinationFilterType::New();
1676           combine0bis->SetInput(0,combine0->GetOutput());
1677           combine0bis->SetInput(1,row4);
1678           combine0bis->SetA(1.);
1679           combine0bis->SetB(-1.);
1680           combine0bis->Update();
1681           typename CoefficientImageType::Pointer bc1Row=combine0bis->GetOutput();
1682
1683           typename PasteImageFilterType::Pointer paster0=PasteImageFilterType::New();
1684           paster0->SetDestinationImage(m_PaddedCoefficientImage);
1685           paster0->SetSourceImage(bc1Row);
1686           destinationIndex[InputDimension-1]=0;
1687           paster0->SetDestinationIndex(destinationIndex);
1688           paster0->SetSourceRegion(bc1Row->GetLargestPossibleRegion());
1689
1690           //---------------------------------
1691           // 1. Next two temporal rows are identical: paste 0-1 to 1-2
1692           typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1693           paster1->SetDestinationImage(paster0->GetOutput());
1694           paster1->SetSourceImage(m_CoefficientImage);
1695           sourceRegion.SetIndex(NInputDimensions-1,0);
1696           destinationIndex[InputDimension-1]=1;
1697           sourceRegion.SetSize(NInputDimensions-1,2);
1698           paster1->SetDestinationIndex(destinationIndex);
1699           paster1->SetSourceRegion(sourceRegion);
1700
1701           //---------------------------------
1702           // 2. Middle row at index 3=BC
1703           // Extract row at index 1, 2 (of coeff image)
1704           //typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1705           typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
1706
1707           // Combine
1708           typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1709           combine1->SetInput(0,row1);
1710           combine1->SetInput(1,row3);
1711           combine1->SetA(1.);
1712           combine1->SetB(1.);
1713           combine1->Update();
1714
1715           typename LinearCombinationFilterType::Pointer combine2=LinearCombinationFilterType::New();
1716           combine2->SetInput(0,combine1->GetOutput());
1717           combine2->SetInput(1,row2);
1718           combine2->SetA(-m_WeightRatio[3][1]);
1719           combine2->SetB(-1.);
1720           combine2->Update();
1721           typename CoefficientImageType::Pointer bc2Row=combine2->GetOutput();
1722
1723           // Paste middleRow at index 3 (padded coeff)
1724           typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1725           paster2->SetDestinationImage(paster1->GetOutput());
1726           paster2->SetSourceImage(bc2Row);
1727           destinationIndex[InputDimension-1]=3;
1728           paster2->SetDestinationIndex(destinationIndex);
1729           paster2->SetSourceRegion(bc2Row->GetLargestPossibleRegion());
1730
1731           //---------------------------------
1732           // 3. Next two temporal rows are identical: paste 2,3 to 4,5
1733           typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1734           paster3->SetDestinationImage(paster2->GetOutput());
1735           paster3->SetSourceImage(m_CoefficientImage);
1736           sourceRegion.SetIndex(InputDimension-1,2);
1737           destinationIndex[InputDimension-1]=4;
1738           sourceRegion.SetSize(NInputDimensions-1,2);
1739           paster3->SetDestinationIndex(destinationIndex);
1740           paster3->SetSourceRegion(sourceRegion);
1741
1742           //---------------------------------
1743           // 4. Final row at index 6=BC3
1744           typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1745
1746           // Combine
1747           typename LinearCombinationFilterType::Pointer combine3=LinearCombinationFilterType::New();
1748           combine3->SetInput(0,row1);
1749           combine3->SetInput(1,row4);
1750           combine3->SetA(1.);
1751           combine3->SetB(-1.);
1752           typename LinearCombinationFilterType::Pointer combine4=LinearCombinationFilterType::New();
1753           combine4->SetInput(0,row0);
1754           combine4->SetInput(1,combine3->GetOutput());
1755           combine4->SetA(1.);
1756           combine4->SetB(.5);
1757           combine4->Update();
1758           typename CoefficientImageType::Pointer bc3Row=combine4->GetOutput();
1759
1760           // Paste BC row at index 6
1761           typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
1762           paster4->SetDestinationImage(paster3->GetOutput());
1763           paster4->SetSourceImage(bc3Row);
1764           destinationIndex[InputDimension-1]=6;
1765           paster4->SetDestinationIndex(destinationIndex);
1766           paster4->SetSourceRegion(bc3Row->GetLargestPossibleRegion());
1767
1768           // Paste row 4  at index 7
1769           typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
1770           paster5->SetDestinationImage(paster4->GetOutput());
1771           paster5->SetSourceImage(row4);
1772           destinationIndex[InputDimension-1]=7;
1773           paster5->SetDestinationIndex(destinationIndex);
1774           paster5->SetSourceRegion(row4->GetLargestPossibleRegion());
1775
1776           // Update the chain!
1777           paster5->Update();
1778           m_PaddedCoefficientImage= paster5->GetOutput();
1779
1780           break;
1781         }
1782
1783       case 6:
1784         {
1785           // ----------------------------------------------------------------------
1786           // The diamond with 4 internal CP:
1787           // Periodic, constrained to zero at the reference at position 3
1788           // Coeff      row 0   1   2  BC1  3  BC2  4
1789           // PaddedCoeff  R:0   1   2   3   4   5   6
1790           // BC1= -weights[2]/weights[0] ( R2+R4 )
1791           // BC2= weights[2]/weights[0]  ( R0+R2-R4-R6 ) + R1
1792           // ---------------------------------------------------------------------
1793
1794           //---------------------------------
1795           // 1. First Three temporal rows are identical: paste 0-2 to 0-2
1796           typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1797           paster1->SetDestinationImage(m_PaddedCoefficientImage);
1798           paster1->SetSourceImage(m_CoefficientImage);
1799           sourceIndex.Fill(0);
1800           destinationIndex.Fill(0);
1801           sourceSize[NInputDimensions-1]=3;
1802           sourceRegion.SetSize(sourceSize);
1803           sourceRegion.SetIndex(sourceIndex);
1804           paster1->SetDestinationIndex(destinationIndex);
1805           paster1->SetSourceRegion(sourceRegion);
1806
1807           //---------------------------------
1808           // 2. Middle row at index 3=BC
1809           // Extract row at index 0, 4 (of coeff image)
1810           typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
1811           typename CoefficientImageType::Pointer row3=this->ExtractTemporalRow(m_CoefficientImage, 3);
1812
1813           // Combine
1814           typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1815           combine1->SetInput(0,row2);
1816           combine1->SetInput(1,row3);
1817           combine1->SetA(-m_WeightRatio[2][0]);
1818           combine1->SetB(-m_WeightRatio[2][0]);
1819           combine1->Update();
1820           typename CoefficientImageType::Pointer bc1Row=combine1->GetOutput();
1821
1822           // Paste middleRow at index 3 (padded coeff)
1823           typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1824           paster2->SetDestinationImage(paster1->GetOutput());
1825           paster2->SetSourceImage(bc1Row);
1826           destinationIndex[InputDimension-1]=3;
1827           paster2->SetDestinationIndex(destinationIndex);
1828           paster2->SetSourceRegion(bc1Row->GetLargestPossibleRegion());
1829
1830           //---------------------------------
1831           // 3. Next  row identical: paste 3 to 4
1832           typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1833           paster3->SetDestinationImage(paster2->GetOutput());
1834           paster3->SetSourceImage(row3);
1835           destinationIndex[InputDimension-1]=4;
1836           paster3->SetDestinationIndex(destinationIndex);
1837           paster3->SetSourceRegion(row3->GetLargestPossibleRegion());
1838
1839           //---------------------------------
1840           // 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
1841           // Extract row at index 0, 2 (of coeff image)already done
1842           typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1843           typename CoefficientImageType::Pointer row4=this->ExtractTemporalRow(m_CoefficientImage, 4);
1844           typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1845
1846           // Combine
1847           typename LinearCombinationFilterType::Pointer combine3=LinearCombinationFilterType::New();
1848           combine3->SetInput(0,row0);
1849           combine3->SetInput(1,row2);
1850           combine3->SetA(1.);
1851           combine3->SetB(1.);
1852
1853           // Combine
1854           typename LinearCombinationFilterType::Pointer combine4=LinearCombinationFilterType::New();
1855           combine4->SetInput(0,row3);
1856           combine4->SetInput(1,row4);
1857           combine4->SetA(1.);
1858           combine4->SetB(1.);
1859
1860           // Combine
1861           typename LinearCombinationFilterType::Pointer combine5=LinearCombinationFilterType::New();
1862           combine5->SetInput(0,combine3->GetOutput());
1863           combine5->SetInput(1,combine4->GetOutput());
1864           combine5->SetA(1.);
1865           combine5->SetB(-1.);
1866
1867           // Combine
1868           typename LinearCombinationFilterType::Pointer combine6=LinearCombinationFilterType::New();
1869           combine6->SetInput(0,combine5->GetOutput());
1870           combine6->SetInput(1,row1);
1871           combine6->SetA(m_WeightRatio[2][0]);
1872           combine6->SetB(1.);
1873           combine6->Update();
1874           typename CoefficientImageType::Pointer bc2Row=combine6->GetOutput();
1875
1876           // Paste BC row at index 5
1877           typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
1878           paster4->SetDestinationImage(paster3->GetOutput());
1879           paster4->SetSourceImage(bc2Row);
1880           destinationIndex[InputDimension-1]=5;
1881           paster4->SetDestinationIndex(destinationIndex);
1882           paster4->SetSourceRegion(bc2Row->GetLargestPossibleRegion());
1883
1884           // Paste last row at index 6
1885           typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
1886           paster5->SetDestinationImage(paster4->GetOutput());
1887           paster5->SetSourceImage(row4);
1888           destinationIndex[InputDimension-1]=6;
1889           paster5->SetDestinationIndex(destinationIndex);
1890           paster5->SetSourceRegion(row4->GetLargestPossibleRegion());
1891
1892           // Update the chain!
1893           paster5->Update();
1894           m_PaddedCoefficientImage= paster5->GetOutput();
1895
1896           break;
1897         }
1898       case 7:
1899         {
1900           // ----------------------------------------------------------------------
1901           // The diamond with 5 internal CP:
1902           // periodic, constrained to zero at the reference at position 3.5
1903           // Coeff      row 0   1   2  BC1  3   4  BC2  5
1904           // PaddedCoeff  R:0   1   2   3   4   5   6   7
1905           // BC1= -weights[3]/weights[1] ( R2+R5 ) - R4
1906           // BC2= weights[2]/weights[0] ( R0+R2-R5-R7 ) + R1
1907           // ---------------------------------------------------------------------
1908
1909           //---------------------------------
1910           // 1. First Three temporal rows are identical: paste 0-2 to 0-2
1911           typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1912           paster1->SetDestinationImage(m_PaddedCoefficientImage);
1913           paster1->SetSourceImage(m_CoefficientImage);
1914           sourceIndex.Fill(0);
1915           destinationIndex.Fill(0);
1916           sourceSize[NInputDimensions-1]=3;
1917           sourceRegion.SetSize(sourceSize);
1918           sourceRegion.SetIndex(sourceIndex);
1919           paster1->SetDestinationIndex(destinationIndex);
1920           paster1->SetSourceRegion(sourceRegion);
1921
1922           //---------------------------------
1923           // 2. Middle row at index 3=BC
1924           // Extract row at index 0, 4 (of coeff image)
1925           typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
1926           typename CoefficientImageType::Pointer row4=this->ExtractTemporalRow(m_CoefficientImage, 4);
1927
1928           // Combine
1929           typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1930           combine1->SetInput(0,row2);
1931           combine1->SetInput(1,row4);
1932           combine1->SetA(1.);
1933           combine1->SetB(1.);
1934           combine1->Update();
1935
1936           // Extract row at index 3 (coeff image)
1937           typename CoefficientImageType::Pointer row3=this->ExtractTemporalRow(m_CoefficientImage, 3);
1938
1939           // Combine
1940           typename LinearCombinationFilterType::Pointer combine2=LinearCombinationFilterType::New();
1941           combine2->SetInput(0,combine1->GetOutput());
1942           combine2->SetInput(1,row3);
1943           combine2->SetA(-m_WeightRatio[3][1] );
1944           combine2->SetB(-1.);
1945           combine2->Update();
1946           typename CoefficientImageType::Pointer bc1Row=combine2->GetOutput();
1947
1948           // Paste middleRow at index 3 (padded coeff)
1949           typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1950           paster2->SetDestinationImage(paster1->GetOutput());
1951           paster2->SetSourceImage(bc1Row);
1952           destinationIndex[InputDimension-1]=3;
1953           paster2->SetDestinationIndex(destinationIndex);
1954           paster2->SetSourceRegion(bc1Row->GetLargestPossibleRegion());
1955
1956           //---------------------------------
1957           // 3. Next two temporal rows are identical: paste 3,4 to 4,5
1958           typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1959           paster3->SetDestinationImage(paster2->GetOutput());
1960           paster3->SetSourceImage(m_CoefficientImage);
1961           sourceIndex[InputDimension-1]=3;
1962           destinationIndex[InputDimension-1]=4;
1963           sourceSize[NInputDimensions-1]=2;
1964           sourceRegion.SetSize(sourceSize);
1965           sourceRegion.SetIndex(sourceIndex);
1966           paster3->SetDestinationIndex(destinationIndex);
1967           paster3->SetSourceRegion(sourceRegion);
1968
1969           //---------------------------------
1970           // 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
1971           // Extract row at index 0, 2 (of coeff image)already done
1972           typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1973           typename CoefficientImageType::Pointer row5=this->ExtractTemporalRow(m_CoefficientImage, 5);
1974           typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1975
1976           // Combine
1977           typename LinearCombinationFilterType::Pointer combine3=LinearCombinationFilterType::New();
1978           combine3->SetInput(0,row0);
1979           combine3->SetInput(1,row2);
1980           combine3->SetA(1.);
1981           combine3->SetB(1.);
1982
1983           // Combine
1984           typename LinearCombinationFilterType::Pointer combine4=LinearCombinationFilterType::New();
1985           combine4->SetInput(0,row4);
1986           combine4->SetInput(1,row5);
1987           combine4->SetA(1.);
1988           combine4->SetB(1.);
1989
1990           // Combine
1991           typename LinearCombinationFilterType::Pointer combine5=LinearCombinationFilterType::New();
1992           combine5->SetInput(0,combine3->GetOutput());
1993           combine5->SetInput(1,combine4->GetOutput());
1994           combine5->SetA(1.);
1995           combine5->SetB(-1.);
1996
1997           // Combine
1998           typename LinearCombinationFilterType::Pointer combine6=LinearCombinationFilterType::New();
1999           combine6->SetInput(0,combine5->GetOutput());
2000           combine6->SetInput(1,row1);
2001           combine6->SetA(m_WeightRatio[2][0]);
2002           combine6->SetB(1.);
2003           combine6->Update();
2004           typename CoefficientImageType::Pointer bc2Row=combine6->GetOutput();
2005
2006           // Paste BC row at index 6
2007           typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
2008           paster4->SetDestinationImage(paster3->GetOutput());
2009           paster4->SetSourceImage(bc2Row);
2010           destinationIndex[InputDimension-1]=6;
2011           paster4->SetDestinationIndex(destinationIndex);
2012           paster4->SetSourceRegion(bc2Row->GetLargestPossibleRegion());
2013
2014           // Paste last row at index 7
2015           typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
2016           paster5->SetDestinationImage(paster4->GetOutput());
2017           paster5->SetSourceImage(row5);
2018           destinationIndex[InputDimension-1]=7;
2019           paster5->SetDestinationIndex(destinationIndex);
2020           paster5->SetSourceRegion(row5->GetLargestPossibleRegion());
2021
2022           // Update the chain!
2023           paster5->Update();
2024           m_PaddedCoefficientImage= paster5->GetOutput();
2025
2026           break;
2027         }
2028
2029       case 9:
2030         {
2031           // ----------------------------------------------------------------------
2032           // The sputnik with 5 internal CP T''(0)=T''(10)
2033           // Periodic, constrained to zero at the reference
2034           // at position 3 and one indepent extreme
2035           // Coeff     row BC1  0   1  BC2  2   3  BC3  4
2036           // PaddedCoeff R: 0   1   2   3   4   5   6   7
2037           // BC1= -R2+R5+R7
2038           // BC2= -weights[3]/weights[1] ( R2+R5 ) - R4
2039           // BC3= R1
2040           // ---------------------------------------------------------------------
2041
2042           //---------------------------------
2043           // 1. First Row =BC
2044           typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
2045           typename CoefficientImageType::Pointer row3=this->ExtractTemporalRow(m_CoefficientImage, 3);
2046           typename CoefficientImageType::Pointer row4=this->ExtractTemporalRow(m_CoefficientImage, 4);
2047
2048           // Combine
2049           typename LinearCombinationFilterType::Pointer combine0=LinearCombinationFilterType::New();
2050           combine0->SetInput(0,row3);
2051           combine0->SetInput(1,row4);
2052           combine0->SetA(1.);
2053           combine0->SetB(1.);
2054           typename LinearCombinationFilterType::Pointer combine0bis=LinearCombinationFilterType::New();
2055           combine0bis->SetInput(0,combine0->GetOutput());
2056           combine0bis->SetInput(1,row1);
2057           combine0bis->SetA(1.);
2058           combine0bis->SetB(-1.);
2059           combine0bis->Update();
2060           typename CoefficientImageType::Pointer bc1Row=combine0bis->GetOutput();
2061
2062           typename PasteImageFilterType::Pointer paster0=PasteImageFilterType::New();
2063           paster0->SetDestinationImage(m_PaddedCoefficientImage);
2064           paster0->SetSourceImage(bc1Row);
2065           destinationIndex[InputDimension-1]=0;
2066           paster0->SetDestinationIndex(destinationIndex);
2067           paster0->SetSourceRegion(bc1Row->GetLargestPossibleRegion());
2068
2069           //---------------------------------
2070           // 1. Next two temporal rows are identical: paste 0-1 to 1-2
2071           typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
2072           paster1->SetDestinationImage(paster0->GetOutput());
2073           paster1->SetSourceImage(m_CoefficientImage);
2074           sourceRegion.SetIndex(NInputDimensions-1,0);
2075           destinationIndex[InputDimension-1]=1;
2076           sourceRegion.SetSize(NInputDimensions-1,2);
2077           paster1->SetDestinationIndex(destinationIndex);
2078           paster1->SetSourceRegion(sourceRegion);
2079
2080           //---------------------------------
2081           // 2. Middle row at index 3=BC
2082           // Extract row at index 1, 2 (of coeff image)
2083           // typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
2084           typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
2085
2086           // Combine
2087           typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
2088           combine1->SetInput(0,row1);
2089           combine1->SetInput(1,row3);
2090           combine1->SetA(1.);
2091           combine1->SetB(1.);
2092           combine1->Update();
2093
2094           typename LinearCombinationFilterType::Pointer combine2=LinearCombinationFilterType::New();
2095           combine2->SetInput(0,combine1->GetOutput());
2096           combine2->SetInput(1,row2);
2097           combine2->SetA(-m_WeightRatio[3][1]);
2098           combine2->SetB(-1.);
2099           combine2->Update();
2100           typename CoefficientImageType::Pointer bc2Row=combine2->GetOutput();
2101
2102           // Paste middleRow at index 3 (padded coeff)
2103           typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
2104           paster2->SetDestinationImage(paster1->GetOutput());
2105           paster2->SetSourceImage(bc2Row);
2106           destinationIndex[InputDimension-1]=3;
2107           paster2->SetDestinationIndex(destinationIndex);
2108           paster2->SetSourceRegion(bc2Row->GetLargestPossibleRegion());
2109
2110           //---------------------------------
2111           // 3. Next two temporal rows are identical: paste 2,3 to 4,5
2112           typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
2113           paster3->SetDestinationImage(paster2->GetOutput());
2114           paster3->SetSourceImage(m_CoefficientImage);
2115           sourceRegion.SetIndex(InputDimension-1,2);
2116           destinationIndex[InputDimension-1]=4;
2117           sourceRegion.SetSize(NInputDimensions-1,2);
2118           paster3->SetDestinationIndex(destinationIndex);
2119           paster3->SetSourceRegion(sourceRegion);
2120
2121           //---------------------------------
2122           // 4. Final row at index 6=BC3
2123           typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
2124
2125           //      // Combine
2126           //      typename LinearCombinationFilterType::Pointer combine3=LinearCombinationFilterType::New();
2127           //      combine3->SetInput(0,row0);
2128           //      combine3->SetInput(1,row1);
2129           //      combine3->SetA(1.);
2130           //      combine3->SetB(0.5);
2131           //      combine3->Update();
2132           //      typename CoefficientImageType::Pointer bc3Row=combine3->GetOutput();
2133
2134           // Paste BC row at index 6
2135           typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
2136           paster4->SetDestinationImage(paster3->GetOutput());
2137           paster4->SetSourceImage(row0);
2138           destinationIndex[InputDimension-1]=6;
2139           paster4->SetDestinationIndex(destinationIndex);
2140           paster4->SetSourceRegion(row0->GetLargestPossibleRegion());
2141
2142           // Paste row 4  at index 7
2143           typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
2144           paster5->SetDestinationImage(paster4->GetOutput());
2145           paster5->SetSourceImage(row4);
2146           destinationIndex[InputDimension-1]=7;
2147           paster5->SetDestinationIndex(destinationIndex);
2148           paster5->SetSourceRegion(row4->GetLargestPossibleRegion());
2149
2150           // Update the chain!
2151           paster5->Update();
2152           m_PaddedCoefficientImage= paster5->GetOutput();
2153
2154           break;
2155         }
2156
2157       default:
2158         DD ("Shape not available");
2159       }
2160
2161   }
2162
2163
2164   // // Extract coefficients from padded version
2165   // template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
2166   // void
2167   // ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2168   // ::ExtractCoefficientImage( )
2169   // {
2170   //   ////DD("extract coeff image");
2171   //   typename ExtractImageFilterType::Pointer extractImageFilter = ExtractImageFilterType::New();
2172   //   typename CoefficientImageType::RegionType extractionRegion=m_PaddedCoefficientImage->GetLargestPossibleRegion();
2173   //   typename CoefficientImageType::RegionType::SizeType extractionSize=extractionRegion.GetSize();
2174   //   typename CoefficientImageType::RegionType::IndexType extractionIndex = extractionRegion.GetIndex();
2175   //   extractionSize[InputDimension-1]-=4;
2176   //   extractionIndex[InputDimension-1]=2;
2177   //   extractionRegion.SetSize(extractionSize);
2178   //   extractionRegion.SetIndex(extractionIndex);
2179   //   extractImageFilter->SetInput(m_PaddedCoefficientImage);
2180   //   extractImageFilter->SetExtractionRegion(extractionRegion);
2181   //   extractImageFilter->Update();
2182   //   m_CoefficientImage=extractImageFilter->GetOutput();
2183   // }
2184
2185
2186   // Print self
2187   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
2188   void
2189   ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2190   ::PrintSelf(std::ostream &os, itk::Indent indent) const
2191   {
2192
2193     this->Superclass::PrintSelf(os, indent);
2194
2195     os << indent << "GridRegion: " << m_GridRegion << std::endl;
2196     os << indent << "GridOrigin: " << m_GridOrigin << std::endl;
2197     os << indent << "GridSpacing: " << m_GridSpacing << std::endl;
2198     os << indent << "GridDirection: " << m_GridDirection << std::endl;
2199     os << indent << "IndexToPoint: " << m_IndexToPoint << std::endl;
2200     os << indent << "PointToIndex: " << m_PointToIndex << std::endl;
2201
2202     os << indent << "CoefficientImage: [ ";
2203     os << m_CoefficientImage.GetPointer() << " ]" << std::endl;
2204
2205     os << indent << "WrappedImage: [ ";
2206     os << m_WrappedImage.GetPointer() << " ]" << std::endl;
2207
2208     os << indent << "InputParametersPointer: "
2209        << m_InputParametersPointer << std::endl;
2210     os << indent << "ValidRegion: " << m_ValidRegion << std::endl;
2211     os << indent << "LastJacobianIndex: " << m_LastJacobianIndex << std::endl;
2212     os << indent << "BulkTransform: ";
2213     os << m_BulkTransform << std::endl;
2214
2215     if ( m_BulkTransform )
2216       {
2217         os << indent << "BulkTransformType: "
2218            << m_BulkTransform->GetNameOfClass() << std::endl;
2219       }
2220     os << indent << "VectorBSplineInterpolator: ";
2221     os << m_VectorInterpolator.GetPointer() << std::endl;
2222     os << indent << "Mask: ";
2223     os << m_Mask<< std::endl;
2224   }
2225
2226
2227   // Verify
2228   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
2229   bool
2230   ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2231   ::InsideValidRegion( const ContinuousIndexType& index ) const
2232   {
2233     bool inside = true;
2234
2235     if ( !m_ValidRegion.IsInside( index ) )
2236       {
2237         inside = false;
2238       }
2239
2240     // JV verify all dimensions
2241     if ( inside)
2242       {
2243         typedef typename ContinuousIndexType::ValueType ValueType;
2244         for( unsigned int j = 0; j < NInputDimensions; j++ )
2245           {
2246             if (m_SplineOrderOdd[j])
2247               {
2248                 if ( index[j] >= static_cast<ValueType>( m_ValidRegionLast[j] ) )
2249                   {
2250                     inside = false;
2251                     break;
2252                   }
2253               }
2254           }
2255       }
2256     return inside;
2257   }
2258
2259
2260   // Transform a point
2261   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
2262   typename ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2263   ::OutputPointType
2264   ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2265   ::TransformPoint(const InputPointType &inputPoint) const
2266   {
2267
2268     InputPointType transformedPoint;
2269     OutputPointType outputPoint;
2270
2271     // BulkTransform
2272     if ( m_BulkTransform )
2273       {
2274         transformedPoint = m_BulkTransform->TransformPoint( inputPoint );
2275       }
2276     else
2277       {
2278         transformedPoint = inputPoint;
2279       }
2280
2281     // Deformable transform
2282     if ( m_PaddedCoefficientImage )
2283       {
2284         // Check if inside mask
2285         if(m_Mask &&  !(m_Mask->IsInside(inputPoint) ) )
2286           {
2287             // Outside: no (deformable) displacement
2288             return transformedPoint;
2289           }
2290
2291         // Check if inside valid region
2292         bool inside = true;
2293         ContinuousIndexType index;
2294         this->TransformPointToContinuousIndex( inputPoint, index );
2295         inside = this->InsideValidRegion( index );
2296         if ( !inside )
2297           {
2298             // Outside: no (deformable) displacement
2299             DD("outside valid region");
2300             return transformedPoint;
2301           }
2302
2303         // Call the vector interpolator
2304         itk::Vector<TCoordRep,SpaceDimension> displacement=m_VectorInterpolator->EvaluateAtContinuousIndex(index);
2305
2306         // JV add for the spatial dimensions
2307         outputPoint=transformedPoint;
2308         for (unsigned int i=0; i<NInputDimensions-1; i++)
2309           outputPoint[i] += displacement[i];
2310
2311       }
2312
2313     else
2314       {
2315         itkWarningMacro( << "B-spline coefficients have not been set" );
2316         outputPoint = transformedPoint;
2317       }
2318
2319     return outputPoint;
2320   }
2321
2322
2323
2324   //JV Deformably transform a point
2325   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
2326   typename ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2327   ::OutputPointType
2328   ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2329   ::DeformablyTransformPoint(const InputPointType &inputPoint) const
2330   {
2331     OutputPointType outputPoint;
2332     if ( m_PaddedCoefficientImage )
2333       {
2334
2335         // Check if inside mask
2336         if(m_Mask && !(m_Mask->IsInside(inputPoint) ) )
2337           {
2338             // Outside: no (deformable) displacement
2339             return inputPoint;
2340           }
2341
2342         // Check if inside
2343         bool inside = true;
2344         ContinuousIndexType index;
2345         this->TransformPointToContinuousIndex( inputPoint, index );
2346         inside = this->InsideValidRegion( index );
2347
2348         if ( !inside )
2349           {
2350             //outside: no (deformable) displacement
2351             outputPoint = inputPoint;
2352             return outputPoint;
2353           }
2354
2355         // Call the vector interpolator
2356         itk::Vector<TCoordRep,SpaceDimension> displacement=m_VectorInterpolator->EvaluateAtContinuousIndex(index);
2357
2358         // JV add for the spatial dimensions
2359         outputPoint=inputPoint;
2360         for (unsigned int i=0; i<NInputDimensions-1; i++)
2361           outputPoint[i] += displacement[i];
2362       }
2363
2364     // No coefficients available
2365     else
2366       {
2367         itkWarningMacro( << "B-spline coefficients have not been set" );
2368         outputPoint = inputPoint;
2369       }
2370
2371     return outputPoint;
2372   }
2373
2374
2375   // JV weights are identical as for transformpoint, could be done simultaneously in metric!!!!
2376   // Compute the Jacobian in one position
2377   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
2378   const
2379   typename ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2380   ::JacobianType &
2381   ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2382   ::GetJacobian( const InputPointType & point ) const
2383   {
2384
2385     //========================================================
2386     // Zero all components of jacobian
2387     //========================================================
2388     // JV not thread safe (m_LastJacobianIndex), instantiate N transforms
2389     // NOTE: for efficiency, we only need to zero out the coefficients
2390     // that got fill last time this method was called.
2391     unsigned int j=0,b=0;
2392
2393     //Set the previously-set to zero
2394     for ( j = 0; j < SpaceDimension; j++ )
2395       {
2396         m_FirstIterator[j].GoToBegin();
2397         while ( !m_FirstIterator[j].IsAtEnd() )
2398           {
2399             m_FirstIterator[j].Set( m_ZeroVector );
2400             ++(m_FirstIterator[j]);
2401           }
2402       }
2403
2404     //Set the previously-set to zero
2405     for ( j = 0; j < SpaceDimension; j++ )
2406       {
2407         m_SecondIterator[j].GoToBegin();
2408         while ( ! (m_SecondIterator[j]).IsAtEnd() )
2409           {
2410             m_SecondIterator[j].Set( m_ZeroVector );
2411             ++(m_SecondIterator[j]);
2412           }
2413       }
2414
2415     //Set the previously-set to zero
2416     if (m_ThirdSize)
2417       for ( j = 0; j < SpaceDimension; j++ )
2418         {
2419           m_ThirdIterator[j].GoToBegin();
2420           while ( ! (m_ThirdIterator[j]).IsAtEnd() )
2421             {
2422               m_ThirdIterator[j].Set( m_ZeroVector );
2423               ++(m_ThirdIterator[j]);
2424             }
2425         }
2426
2427     //Set the previously-set to zero
2428     if (m_BCSize)
2429       for (b=0; b<m_BCSize;b++)
2430         if ( !( m_FirstRegion.IsInside(m_BCRegions[b]) | m_SecondRegion.IsInside(m_BCRegions[b]) )  )
2431           for ( j = 0; j < SpaceDimension; j++ )
2432             {
2433               m_BCIterators[j][b].GoToBegin();
2434               while ( ! (m_BCIterators[j][b]).IsAtEnd() )
2435                 {
2436                   m_BCIterators[j][b].Set( m_ZeroVector );
2437                   ++(m_BCIterators[j][b]);
2438                 }
2439             }
2440
2441     //Set the previously-set to zero
2442     if (m_BC2Size)
2443       for (b=0; b<m_BC2Size;b++)
2444         if ( !( m_FirstRegion.IsInside(m_BC2Regions[b]) | m_SecondRegion.IsInside(m_BC2Regions[b]) )  )
2445           for ( j = 0; j < SpaceDimension; j++ )
2446             {
2447               m_BC2Iterators[j][b].GoToBegin();
2448               while ( ! (m_BC2Iterators[j][b]).IsAtEnd() )
2449                 {
2450                   m_BC2Iterators[j][b].Set( m_ZeroVector );
2451                   ++(m_BC2Iterators[j][b]);
2452                 }
2453             }
2454
2455     //Set the previously-set to zero
2456     if (m_BC3Size)
2457       for (b=0; b<m_BC3Size;b++)
2458         if ( !( m_FirstRegion.IsInside(m_BC3Regions[b]) | m_SecondRegion.IsInside(m_BC3Regions[b]) )  )
2459           for ( j = 0; j < SpaceDimension; j++ )
2460             {
2461               m_BC3Iterators[j][b].GoToBegin();
2462               while ( ! (m_BC3Iterators[j][b]).IsAtEnd() )
2463                 {
2464                   m_BC3Iterators[j][b].Set( m_ZeroVector );
2465                   ++(m_BC3Iterators[j][b]);
2466                 }
2467             }
2468
2469
2470     //========================================================
2471     // For each dimension, copy the weight to the support region
2472     //========================================================
2473
2474     // Check if inside mask
2475     if(m_Mask &&  !(m_Mask->IsInside(point) ) )
2476       {
2477         // Outside: no (deformable) displacement
2478         return this->m_Jacobian;
2479       }
2480
2481     // Get index
2482     this->TransformPointToContinuousIndex( point, m_Index );
2483
2484     // NOTE: if the support region does not lie totally within the grid
2485     // we assume zero displacement and return the input point
2486     if ( !this->InsideValidRegion( m_Index ) )
2487       {
2488         return this->m_Jacobian;
2489       }
2490
2491     // Compute interpolation weights
2492     const WeightsDataType *weights=NULL;
2493     m_VectorInterpolator->EvaluateWeightsAtContinuousIndex( m_Index, &weights, m_LastJacobianIndex);
2494
2495     // Get support
2496     m_SupportRegion.SetIndex( m_LastJacobianIndex );
2497     WrapRegion(m_SupportRegion, m_FirstRegion, m_SecondRegion, m_ThirdRegion, m_BCRegions, m_BCValues, m_BC2Regions, m_BC2Values, m_BC3Regions, m_BC3Values, m_InitialOffset);
2498     m_ThirdSize=m_ThirdRegion.GetSize()[InputDimension-1];
2499     m_BCSize=m_BCRegions.size();
2500     m_BC2Size=m_BC2Regions.size();
2501     m_BC3Size=m_BC3Regions.size();
2502
2503     // Reset the iterators
2504     for ( j = 0; j < SpaceDimension ; j++ )
2505       {
2506         m_FirstIterator[j] = IteratorType( m_JacobianImage[j], m_FirstRegion);
2507         m_SecondIterator[j] = IteratorType( m_JacobianImage[j], m_SecondRegion);
2508         if(m_ThirdSize) m_ThirdIterator[j] = IteratorType( m_JacobianImage[j], m_ThirdRegion);
2509
2510         m_BCIterators[j].resize(m_BCSize);
2511         for (b=0; b<m_BCSize;b++)
2512           m_BCIterators[j][b]= IteratorType( m_JacobianImage[j], m_BCRegions[b]);
2513         m_BC2Iterators[j].resize(m_BC2Size);
2514         for (b=0; b<m_BC2Size;b++)
2515           m_BC2Iterators[j][b]= IteratorType( m_JacobianImage[j], m_BC2Regions[b]);
2516         m_BC3Iterators[j].resize(m_BC3Size);
2517         for (b=0; b<m_BC3Size;b++)
2518           m_BC3Iterators[j][b]= IteratorType( m_JacobianImage[j], m_BC3Regions[b]);
2519       }
2520
2521     // Skip if on a fixed condition
2522     if(m_InitialOffset)
2523       {
2524         if (m_BCSize)  weights+=m_InitialOffset*m_BCRegions[0].GetNumberOfPixels();
2525         else std::cerr<<"InitialOffset without BCSize: Error!!!!!!!"<<std::endl;
2526       }
2527
2528     //copy weight to jacobian image
2529     for ( j = 0; j < SpaceDimension; j++ )
2530       {
2531         // For each dimension, copy the weight to the support region
2532         while ( ! (m_FirstIterator[j]).IsAtEnd() )
2533           {
2534             m_ZeroVector[j]=*weights;
2535             (m_FirstIterator[j]).Set( m_ZeroVector);
2536             ++(m_FirstIterator[j]);
2537             weights++;
2538           }
2539
2540         m_ZeroVector[j]=itk::NumericTraits<JacobianValueType>::Zero;
2541         if (j != SpaceDimension-1)  weights-=m_FirstRegion.GetNumberOfPixels();
2542       }
2543
2544     // Skip BC1 and go to the second region
2545     if (m_BCSize) weights+=m_BCRegions[0].GetNumberOfPixels();
2546
2547     // For each dimension, copy the weight to the support region
2548     //copy weight to jacobian image
2549     for ( j = 0; j < SpaceDimension; j++ )
2550       {
2551         while ( ! (m_SecondIterator[j]).IsAtEnd() )
2552           {
2553             m_ZeroVector[j]=*weights;
2554             (m_SecondIterator[j]).Set( m_ZeroVector);
2555             ++(m_SecondIterator[j]);
2556             weights++;
2557           }
2558         weights-=m_SecondRegion.GetNumberOfPixels();
2559         m_ZeroVector[j]=itk::NumericTraits<JacobianValueType>::Zero;
2560       }
2561
2562     // Now do BC1:
2563     if (m_BCSize)
2564       {
2565         //  Put pointer in correct position
2566         weights-=m_BCRegions[0].GetNumberOfPixels();
2567
2568         for ( j = 0; j < SpaceDimension; j++ )
2569           {
2570             for ( b=0; b < m_BCSize; b++ )
2571               {
2572                 while ( ! (m_BCIterators[j][b]).IsAtEnd() )
2573                   {
2574                     //copy weight to jacobian image
2575                     m_ZeroVector[j]=(*weights) * m_BCValues[b];
2576                     (m_BCIterators[j][b]).Value()+= m_ZeroVector;
2577                     ++(m_BCIterators[j][b]);
2578                     weights++;
2579                   }
2580                 weights-=m_BCRegions[b].GetNumberOfPixels();
2581               }
2582             m_ZeroVector[j]=itk::NumericTraits<JacobianValueType>::Zero;
2583           }
2584         weights+=m_BCRegions[m_BCSize-1].GetNumberOfPixels();
2585       }
2586
2587     // Add the BC2 to the weights
2588     if (m_BC2Size)
2589       {
2590         // Move further in the weights pointer
2591         weights+=m_SecondRegion.GetNumberOfPixels();
2592
2593         for ( j = 0; j < SpaceDimension; j++ )
2594           {
2595             for ( b=0; b < m_BC2Size; b++ )
2596               {
2597                 while ( ! (m_BC2Iterators[j][b]).IsAtEnd() )
2598                   {
2599                     //copy weight to jacobian image
2600                     m_ZeroVector[j]=(*weights) * m_BC2Values[b];
2601                     (m_BC2Iterators[j][b]).Value()+= m_ZeroVector;
2602                     ++(m_BC2Iterators[j][b]);
2603                     weights++;
2604                   }
2605                 weights-=m_BC2Regions[b].GetNumberOfPixels();
2606               }
2607             m_ZeroVector[j]=itk::NumericTraits<JacobianValueType>::Zero;
2608           }
2609         // Move further in the weights pointer
2610         weights+=m_BC2Regions[m_BC2Size-1].GetNumberOfPixels();
2611       }
2612
2613     // Third Region
2614     if(m_ThirdSize)
2615       {
2616         // For each dimension, copy the weight to the support region
2617         //copy weight to jacobian image
2618         for ( j = 0; j < SpaceDimension; j++ )
2619           {
2620             while ( ! (m_ThirdIterator[j]).IsAtEnd() )
2621               {
2622                 m_ZeroVector[j]=*weights;
2623                 (m_ThirdIterator[j]).Value()+= m_ZeroVector;
2624                 ++(m_ThirdIterator[j]);
2625                 weights++;
2626               }
2627
2628             // Move further in the weights pointer?
2629             if (j != SpaceDimension-1) weights-=m_ThirdRegion.GetNumberOfPixels();
2630             m_ZeroVector[j]=itk::NumericTraits<JacobianValueType>::Zero;
2631           }
2632       }
2633
2634     // Add the BC3 to the weights
2635     if (m_BC3Size)
2636       {
2637         for ( j = 0; j < SpaceDimension; j++ )
2638           {
2639             for ( b=0; b < m_BC3Size; b++ )
2640               {
2641                 while ( ! (m_BC3Iterators[j][b]).IsAtEnd() )
2642                   {
2643                     //copy weight to jacobian image
2644                     m_ZeroVector[j]=(*weights) * m_BC3Values[b];
2645                     (m_BC3Iterators[j][b]).Value()+= m_ZeroVector;
2646                     ++(m_BC3Iterators[j][b]);
2647                     weights++;
2648                   }
2649                 weights-=m_BC3Regions[b].GetNumberOfPixels();
2650               }
2651             m_ZeroVector[j]=itk::NumericTraits<JacobianValueType>::Zero;
2652           }
2653       }
2654
2655     // Return the result
2656     return this->m_Jacobian;
2657   }
2658
2659
2660   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
2661   inline void
2662   ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2663   ::WrapRegion( const RegionType & m_SupportRegion,
2664                 RegionType & m_FirstRegion,
2665                 RegionType & m_SecondRegion,
2666                 RegionType & m_ThirdRegion,
2667                 std::vector<RegionType>& m_BCRegions,std::vector<double>& m_BCValues,
2668                 std::vector<RegionType>& m_BC2Regions,std::vector<double>& m_BC2Values,
2669                 std::vector<RegionType>& m_BC3Regions,std::vector<double>& m_BC3Values,
2670                 unsigned int& m_InitialOffset ) const
2671   {
2672
2673     // Synchronize regions
2674     m_InitialOffset=0;
2675     m_FirstRegion=m_SupportRegion;
2676     m_BCRegion=m_SupportRegion;
2677     m_BCRegion.SetSize(InputDimension-1,1);
2678     m_SecondRegion=m_SupportRegion;
2679     m_ThirdRegion=m_SupportRegion;
2680     m_ThirdRegion.SetSize(InputDimension-1,0);
2681     m_BC3Regions.resize(0);
2682
2683
2684     // BC depends on shape
2685     switch(m_TransformShape)
2686       {
2687         /*  The shapes are
2688             0: egg     4 CP  3 DOF
2689             1: egg     5 CP  4 DOF
2690             2: rabbit  4 CP  3 DOF
2691             3: rabbit  5 CP  4 DOF
2692             4: sputnik 4 CP  4 DOF
2693             5: sputnik 5 CP  5 DOF
2694             6: diamond 6 CP  5 DOF
2695             7: diamond 7 CP  6 DOF
2696         */
2697
2698       case 0:
2699         {
2700           // ----------------------------------------------------------------------
2701           // The egg with 4 internal CP (starting from inhale)
2702           // Periodic, constrained to zero at the reference
2703           // at position 3 and
2704           // Coeff      row  BC1 BC2  0  BC3  1   2  BC4
2705           // PaddedCoeff  R:  0   1   2   3   4   5   6
2706           // BC1= R4
2707           // BC2= R5
2708           // BC3= -weights[2]/weights[0] ( R2+R4 )
2709           // BC4= R2
2710           // ---------------------------------------------------------------------
2711           switch(m_SupportRegion.GetIndex(InputDimension-1))
2712             {
2713             case 0:
2714               {
2715                 // Lower end
2716                 m_FirstRegion.SetSize(InputDimension-1,2);
2717                 m_FirstRegion.SetIndex(InputDimension-1,1);
2718
2719                 // BC
2720                 m_BCRegions.resize(0);
2721
2722                 // Second part
2723                 m_SecondRegion.SetSize(InputDimension-1,1);
2724                 m_SecondRegion.SetIndex(InputDimension-1,0);
2725
2726                 // BC
2727                 m_BC2Regions.resize(2);
2728                 m_BC2Values.resize(2);
2729                 m_BCRegion.SetIndex(InputDimension-1,0);
2730                 m_BC2Regions[0]=m_BCRegion;
2731                 m_BC2Values[0]=-m_WeightRatio[2][0];
2732                 m_BCRegion.SetIndex(InputDimension-1,1);
2733                 m_BC2Regions[1]=m_BCRegion;
2734                 m_BC2Values[1]= -m_WeightRatio[2][0];
2735
2736                 break;
2737               }
2738             case 1:
2739               {
2740                 // Lower end
2741                 m_FirstRegion.SetSize(InputDimension-1,1);
2742                 m_FirstRegion.SetIndex(InputDimension-1,2);
2743
2744                 // BC
2745                 m_BCRegions.resize(0);
2746
2747                 // Second part
2748                 m_SecondRegion.SetSize(InputDimension-1,1);
2749                 m_SecondRegion.SetIndex(InputDimension-1,0);
2750
2751                 // BC
2752                 m_BC2Regions.resize(2);
2753                 m_BC2Values.resize(2);
2754                 m_BCRegion.SetIndex(InputDimension-1,0);
2755                 m_BC2Regions[0]=m_BCRegion;
2756                 m_BC2Values[0]=-m_WeightRatio[2][0];
2757                 m_BCRegion.SetIndex(InputDimension-1,1);
2758                 m_BC2Regions[1]=m_BCRegion;
2759                 m_BC2Values[1]= -m_WeightRatio[2][0];
2760
2761                 // Third Part
2762                 m_FirstRegion.SetSize(InputDimension-1,1);
2763                 m_FirstRegion.SetIndex(InputDimension-1,1);
2764
2765                 break;
2766               }
2767             case 2:
2768               {
2769                 // Lower end
2770                 m_FirstRegion.SetSize(InputDimension-1,1);
2771                 m_FirstRegion.SetIndex(InputDimension-1,0);
2772
2773                 // BC
2774                 m_BCRegions.resize(2);
2775                 m_BCValues.resize(2);
2776                 m_BCRegion.SetIndex(InputDimension-1,0);
2777                 m_BCRegions[0]=m_BCRegion;
2778                 m_BCValues[0]=-m_WeightRatio[2][0];
2779                 m_BCRegion.SetIndex(InputDimension-1,1);
2780                 m_BCRegions[1]=m_BCRegion;
2781                 m_BCValues[1]= -m_WeightRatio[2][0];
2782
2783                 // Second part
2784                 m_SecondRegion.SetSize(InputDimension-1,2);
2785                 m_SecondRegion.SetIndex(InputDimension-1,1);
2786
2787                 // BC2
2788                 m_BC2Regions.resize(0);
2789
2790                 break;
2791               }
2792             case 3:
2793               {
2794                 // Lower end
2795                 m_FirstRegion.SetSize(InputDimension-1,0);
2796
2797                 // BC
2798                 m_BCRegions.resize(2);
2799                 m_BCValues.resize(2);
2800                 m_BCRegion.SetIndex(InputDimension-1,0);
2801                 m_BCRegions[0]=m_BCRegion;
2802                 m_BCValues[0]=-m_WeightRatio[2][0];
2803                 m_BCRegion.SetIndex(InputDimension-1,1);
2804                 m_BCRegions[1]=m_BCRegion;
2805                 m_BCValues[1]= -m_WeightRatio[2][0];
2806
2807                 // Second part
2808                 m_SecondRegion.SetSize(InputDimension-1,2);
2809                 m_SecondRegion.SetIndex(InputDimension-1,1);
2810
2811                 // BC2
2812                 m_BC2Regions.resize(1);
2813                 m_BC2Values.resize(1);
2814                 m_BCRegion.SetIndex(InputDimension-1,0);
2815                 m_BC2Regions[0]=m_BCRegion;
2816                 m_BC2Values[0]=1;
2817
2818                 break;
2819               }
2820
2821             default:
2822               {
2823                 DD("supportindex > 3 ???");
2824                 DD(m_SupportRegion.GetIndex(InputDimension-1));
2825               }
2826             } // end switch index
2827
2828           break;
2829         }
2830
2831       case 1:
2832         {
2833           // ----------------------------------------------------------------------
2834           // The egg with 5 internal CP (starting from inhale)
2835           // Periodic, constrained to zero at the reference
2836           // at position 3 and
2837           // Coeff      row  BC1 BC2  0  BC3  1   2   3  BC4
2838           // PaddedCoeff  R:  0   1   2   3   4   5   6   7
2839           // BC1= R5
2840           // BC2= R6
2841           // BC3= -weights[3]/weights[1] ( R2+R5 ) - R4
2842           // BC4= R2
2843           // ---------------------------------------------------------------------
2844           switch(m_SupportRegion.GetIndex(InputDimension-1))
2845             {
2846             case 0:
2847               {
2848                 // Lower end
2849                 m_FirstRegion.SetSize(InputDimension-1,2);
2850                 m_FirstRegion.SetIndex(InputDimension-1,2);
2851
2852                 // BC
2853                 m_BCRegions.resize(0);
2854
2855                 // Second part
2856                 m_SecondRegion.SetSize(InputDimension-1,1);
2857                 m_SecondRegion.SetIndex(InputDimension-1,0);
2858
2859                 // BC
2860                 m_BC2Regions.resize(3);
2861                 m_BC2Values.resize(3);
2862                 m_BCRegion.SetIndex(InputDimension-1,0);
2863                 m_BC2Regions[0]=m_BCRegion;
2864                 m_BC2Values[0]=-m_WeightRatio[3][1];
2865                 m_BCRegion.SetIndex(InputDimension-1,1);
2866                 m_BC2Regions[1]=m_BCRegion;
2867                 m_BC2Values[1]= -1;
2868                 m_BCRegion.SetIndex(InputDimension-1,2);
2869                 m_BC2Regions[2]=m_BCRegion;
2870                 m_BC2Values[2]=-m_WeightRatio[3][1];
2871
2872                 break;
2873               }
2874             case 1:
2875               {
2876                 // Lower end
2877                 m_FirstRegion.SetSize(InputDimension-1,1);
2878                 m_FirstRegion.SetIndex(InputDimension-1,3);
2879
2880                 // BC
2881                 m_BCRegions.resize(0);
2882
2883                 // Second part
2884                 m_SecondRegion.SetSize(InputDimension-1,1);
2885                 m_SecondRegion.SetIndex(InputDimension-1,0);
2886
2887                 // BC
2888                 m_BC2Regions.resize(3);
2889                 m_BC2Values.resize(3);
2890                 m_BCRegion.SetIndex(InputDimension-1,0);
2891                 m_BC2Regions[0]=m_BCRegion;
2892                 m_BC2Values[0]=-m_WeightRatio[3][1];
2893                 m_BCRegion.SetIndex(InputDimension-1,1);
2894                 m_BC2Regions[1]=m_BCRegion;
2895                 m_BC2Values[1]= -1;
2896                 m_BCRegion.SetIndex(InputDimension-1,2);
2897                 m_BC2Regions[2]=m_BCRegion;
2898                 m_BC2Values[2]=-m_WeightRatio[3][1];
2899
2900                 // Third Part
2901                 m_FirstRegion.SetSize(InputDimension-1,1);
2902                 m_FirstRegion.SetIndex(InputDimension-1,1);
2903
2904                 break;
2905               }
2906             case 2:
2907               {
2908                 // Lower end
2909                 m_FirstRegion.SetSize(InputDimension-1,1);
2910                 m_FirstRegion.SetIndex(InputDimension-1,0);
2911
2912                 // BC
2913                 m_BCRegions.resize(3);
2914                 m_BCValues.resize(3);
2915                 m_BCRegion.SetIndex(InputDimension-1,0);
2916                 m_BCRegions[0]=m_BCRegion;
2917                 m_BCValues[0]=-m_WeightRatio[3][1];
2918                 m_BCRegion.SetIndex(InputDimension-1,1);
2919                 m_BCRegions[1]=m_BCRegion;
2920                 m_BCValues[1]= -1;
2921                 m_BCRegion.SetIndex(InputDimension-1,2);
2922                 m_BCRegions[2]=m_BCRegion;
2923                 m_BCValues[2]=-m_WeightRatio[3][1];
2924
2925                 // Second part
2926                 m_SecondRegion.SetSize(InputDimension-1,2);
2927                 m_SecondRegion.SetIndex(InputDimension-1,1);
2928
2929                 // BC2
2930                 m_BC2Regions.resize(0);
2931
2932                 break;
2933               }
2934             case 3:
2935               {
2936                 // Lower end
2937                 m_FirstRegion.SetSize(InputDimension-1,0);
2938
2939                 // BC
2940                 m_BCRegions.resize(3);
2941                 m_BCValues.resize(3);
2942                 m_BCRegion.SetIndex(InputDimension-1,0);
2943                 m_BCRegions[0]=m_BCRegion;
2944                 m_BCValues[0]=-m_WeightRatio[3][1];
2945                 m_BCRegion.SetIndex(InputDimension-1,1);
2946                 m_BCRegions[1]=m_BCRegion;
2947                 m_BCValues[1]= -1;
2948                 m_BCRegion.SetIndex(InputDimension-1,2);
2949                 m_BCRegions[2]=m_BCRegion;
2950                 m_BCValues[2]=-m_WeightRatio[3][1];
2951
2952                 // Second part
2953                 m_SecondRegion.SetSize(InputDimension-1,3);
2954                 m_SecondRegion.SetIndex(InputDimension-1,1);
2955
2956                 // BC2
2957                 m_BC2Regions.resize(0);
2958
2959                 break;
2960               }
2961             case 4:
2962               {
2963                 // Lower end
2964                 m_FirstRegion.SetSize(InputDimension-1,3);
2965                 m_FirstRegion.SetIndex(InputDimension-1,1);
2966
2967                 // BC
2968                 m_BCRegions.resize(0);
2969
2970                 // Second part
2971                 m_SecondRegion.SetSize(InputDimension-1,1);
2972                 m_SecondRegion.SetIndex(InputDimension-1,0);
2973
2974                 // BC2
2975                 m_BC2Regions.resize(0);
2976
2977                 break;
2978               }
2979
2980             default:
2981               {
2982                 DD("supportindex > 3 ???");
2983                 DD(m_SupportRegion.GetIndex(InputDimension-1));
2984               }
2985             } // end switch index
2986
2987           break;
2988         }
2989 //        // ----------------------------------------------------------------------
2990 //        // The egg with 5 internal CP:
2991 //        // Periodic, C2 smooth everywhere and constrained to zero at the reference
2992 //        // Coeff       row R5 BC1  0   1   2   3  BC2  R2
2993 //        // PaddedCoeff R:  0   1   2   3   4   5   6   7
2994 //        // BC1= -weights[2]/weights[0] ( R2+R5)
2995 //        // BC2= BC1
2996 //        // ---------------------------------------------------------------------
2997 //        switch(m_SupportRegion.GetIndex(InputDimension-1))
2998 //          {
2999 //          case 0:
3000 //            {
3001 //              // Lower end
3002 //              m_FirstRegion.SetSize(InputDimension-1,1);
3003 //              m_FirstRegion.SetIndex(InputDimension-1,3);
3004
3005 //              // BC
3006 //              m_BCRegions.resize(2);
3007 //              m_BCValues.resize(2);
3008 //              m_BCRegion.SetIndex(InputDimension-1,0);
3009 //              m_BCRegions[0]=m_BCRegion;
3010 //              m_BCValues[0]=-m_WeightRatio[2][0];
3011 //              m_BCRegion.SetIndex(InputDimension-1,3);
3012 //              m_BCRegions[1]=m_BCRegion;
3013 //              m_BCValues[1]=-m_WeightRatio[2][0];
3014
3015 //              // Second part
3016 //              m_SecondRegion.SetSize(InputDimension-1,2);
3017 //              m_SecondRegion.SetIndex(InputDimension-1,0);
3018
3019 //              // BC2
3020 //              m_BC2Regions.resize(0);
3021
3022 //              break;
3023 //            }
3024 //          case 1:
3025 //            {
3026 //              // Lower end
3027 //              m_FirstRegion.SetSize(InputDimension-1,0);
3028
3029 //              // BC
3030 //              m_BCRegions.resize(2);
3031 //              m_BCValues.resize(2);
3032 //              m_BCRegion.SetIndex(InputDimension-1,0);
3033 //              m_BCRegions[0]=m_BCRegion;
3034 //              m_BCValues[0]=-m_WeightRatio[2][0];
3035 //              m_BCRegion.SetIndex(InputDimension-1,3);
3036 //              m_BCRegions[1]=m_BCRegion;
3037 //              m_BCValues[1]=-m_WeightRatio[2][0];
3038
3039 //              // Second part
3040 //              m_SecondRegion.SetSize(InputDimension-1,3);
3041 //              m_SecondRegion.SetIndex(InputDimension-1,0);
3042
3043 //              // BC2
3044 //              m_BC2Regions.resize(0);
3045
3046 //              break;
3047 //            }
3048 //          case 2:
3049 //            {
3050 //              // Lower end
3051 //              m_FirstRegion.SetSize(InputDimension-1,4);
3052 //              m_FirstRegion.SetIndex(InputDimension-1,0);
3053
3054 //              // BC
3055 //              m_BCRegions.resize(0);
3056
3057 //              // Second part
3058 //              m_SecondRegion.SetSize(InputDimension-1,0);
3059
3060 //              // BC2
3061 //              m_BC2Regions.resize(0);
3062
3063 //              break;
3064 //            }
3065 //          case 3:
3066 //            {
3067 //              // Lower end
3068 //              m_FirstRegion.SetSize(InputDimension-1,3);
3069 //              m_FirstRegion.SetIndex(InputDimension-1,1);
3070
3071 //              // BC
3072 //              m_BCRegions.resize(2);
3073 //              m_BCValues.resize(2);
3074 //              m_BCRegion.SetIndex(InputDimension-1,0);
3075 //              m_BCRegions[0]=m_BCRegion;
3076 //              m_BCValues[0]=-m_WeightRatio[2][0];
3077 //              m_BCRegion.SetIndex(InputDimension-1,3);
3078 //              m_BCRegions[1]=m_BCRegion;
3079 //              m_BCValues[1]=-m_WeightRatio[2][0];
3080
3081 //              // Second part
3082 //              m_SecondRegion.SetSize(InputDimension-1,0);
3083
3084 //              // BC2
3085 //              m_BC2Regions.resize(0);
3086
3087 //              break;
3088 //            }
3089
3090 //          case 4:
3091 //            {
3092 //              // Lower end
3093 //              m_FirstRegion.SetSize(InputDimension-1,2);
3094 //              m_FirstRegion.SetIndex(InputDimension-1,2);
3095
3096 //              // BC
3097 //              m_BCRegions.resize(2);
3098 //              m_BCValues.resize(2);
3099 //              m_BCRegion.SetIndex(InputDimension-1,0);
3100 //              m_BCRegions[0]=m_BCRegion;
3101 //              m_BCValues[0]=-m_WeightRatio[2][0];
3102 //              m_BCRegion.SetIndex(InputDimension-1,3);
3103 //              m_BCRegions[1]=m_BCRegion;
3104 //              m_BCValues[1]=-m_WeightRatio[2][0];
3105
3106 //              // Second part
3107 //              m_SecondRegion.SetSize(InputDimension-1,1);
3108 //              m_SecondRegion.SetIndex(InputDimension-1,0);
3109
3110 //              // BC2
3111 //              m_BC2Regions.resize(0);
3112
3113 //              break;
3114 //            }
3115
3116 //          default:
3117 //            {
3118 //              DD("supportindex > 4 ???");
3119 //              DD(m_SupportRegion.GetIndex(InputDimension-1));
3120 //              DD(m_TransformShape);
3121 //            }
3122 //          }// end swith index
3123 //        break;
3124
3125 //      } // end case 1 shape
3126
3127       case 2:
3128         {
3129           // ----------------------------------------------------------------------
3130           // The rabbit with 4 internal CP
3131           // Periodic, constrained to zero at the reference
3132           // at position 3 and the extremes fixed with anit-symm bc
3133           // Coeff      row  BC1  0   1  BC2  2  BC3 BC4
3134           // PaddedCoeff  R:  0   1   2   3   4   5   6
3135           // BC1= 2*R1-R0
3136           // BC2= -weights[2]/weights[0] ( R2+R4 )
3137           // BC3=  R1
3138           // BC4= 2*R1-R4
3139           // ---------------------------------------------------------------------
3140           switch(m_SupportRegion.GetIndex(InputDimension-1))
3141             {
3142             case 0:
3143               {
3144                 // Lower end
3145                 m_FirstRegion.SetSize(InputDimension-1,0);
3146
3147                 // BC
3148                 m_BCRegions.resize(2);
3149                 m_BCValues.resize(2);
3150                 m_BCRegion.SetIndex(InputDimension-1,0);
3151                 m_BCRegions[0]=m_BCRegion;
3152                 m_BCValues[0]=2;
3153                 m_BCRegion.SetIndex(InputDimension-1,1);
3154                 m_BCRegions[1]=m_BCRegion;
3155                 m_BCValues[1]= -1;
3156
3157                 // Second part
3158                 m_SecondRegion.SetSize(InputDimension-1,2);
3159                 m_SecondRegion.SetIndex(InputDimension-1,0);
3160
3161                 // BC
3162                 m_BC2Regions.resize(2);
3163                 m_BC2Values.resize(2);
3164                 m_BCRegion.SetIndex(InputDimension-1,1);
3165                 m_BC2Regions[0]=m_BCRegion;
3166                 m_BC2Values[0]=-m_WeightRatio[2][0];
3167                 m_BCRegion.SetIndex(InputDimension-1,2);
3168                 m_BC2Regions[1]=m_BCRegion;
3169                 m_BC2Values[1]= -m_WeightRatio[2][0];
3170
3171                 break;
3172               }
3173             case 1:
3174               {
3175                 // Lower end
3176                 m_FirstRegion.SetSize(InputDimension-1,2);
3177                 m_FirstRegion.SetIndex(InputDimension-1,0);
3178
3179                 // BC
3180                 m_BCRegions.resize(2);
3181                 m_BCValues.resize(2);
3182                 m_BCRegion.SetIndex(InputDimension-1,1);
3183                 m_BCRegions[0]=m_BCRegion;
3184                 m_BCValues[0]=-m_WeightRatio[2][0];
3185                 m_BCRegion.SetIndex(InputDimension-1,2);
3186                 m_BCRegions[1]=m_BCRegion;
3187                 m_BCValues[1]= -m_WeightRatio[2][0];
3188
3189                 // Second part
3190                 m_SecondRegion.SetSize(InputDimension-1,1);
3191                 m_SecondRegion.SetIndex(InputDimension-1,2);
3192
3193                 // BC2
3194                 m_BC2Regions.resize(0);
3195
3196                 break;
3197               }
3198             case 2:
3199               {
3200                 // Lower end
3201                 m_FirstRegion.SetSize(InputDimension-1,1);
3202                 m_FirstRegion.SetIndex(InputDimension-1,1);
3203
3204                 // BC
3205                 m_BCRegions.resize(2);
3206                 m_BCValues.resize(2);
3207                 m_BCRegion.SetIndex(InputDimension-1,1);
3208                 m_BCRegions[0]=m_BCRegion;
3209                 m_BCValues[0]=-m_WeightRatio[2][0];
3210                 m_BCRegion.SetIndex(InputDimension-1,2);
3211                 m_BCRegions[1]=m_BCRegion;
3212                 m_BCValues[1]= -m_WeightRatio[2][0];
3213
3214                 // Second part
3215                 m_SecondRegion.SetSize(InputDimension-1,1);
3216                 m_SecondRegion.SetIndex(InputDimension-1,2);
3217
3218                 // BC2
3219                 m_BC2Regions.resize(1);
3220                 m_BC2Values.resize(1);
3221                 m_BCRegion.SetIndex(InputDimension-1,0);
3222                 m_BC2Regions[0]=m_BCRegion;
3223                 m_BC2Values[0]=1;
3224
3225                 break;
3226               }
3227             case 3:
3228               {
3229                 // Lower end
3230                 m_FirstRegion.SetSize(InputDimension-1,0);
3231
3232                 // BC
3233                 m_BCRegions.resize(2);
3234                 m_BCValues.resize(2);
3235                 m_BCRegion.SetIndex(InputDimension-1,1);
3236                 m_BCRegions[0]=m_BCRegion;
3237                 m_BCValues[0]=-m_WeightRatio[2][0];
3238                 m_BCRegion.SetIndex(InputDimension-1,2);
3239                 m_BCRegions[1]=m_BCRegion;
3240                 m_BCValues[1]= -m_WeightRatio[2][0];
3241
3242                 // Second part
3243                 m_SecondRegion.SetSize(InputDimension-1,1);
3244                 m_SecondRegion.SetIndex(InputDimension-1,2);
3245
3246                 // BC2
3247                 m_BC2Regions.resize(1);
3248                 m_BC2Values.resize(1);
3249                 m_BCRegion.SetIndex(InputDimension-1,0);
3250                 m_BC2Regions[0]=m_BCRegion;
3251                 m_BC2Values[0]=1;
3252
3253                 // BC3
3254                 m_BC3Regions.resize(2);
3255                 m_BC3Values.resize(2);
3256                 m_BCRegion.SetIndex(InputDimension-1,0);
3257                 m_BC3Regions[0]=m_BCRegion;
3258                 m_BC3Values[0]=2;
3259                 m_BCRegion.SetIndex(InputDimension-1,2);
3260                 m_BC3Regions[1]=m_BCRegion;
3261                 m_BC3Values[1]= -1;
3262
3263                 break;
3264               }
3265
3266             default:
3267               {
3268                 DD("supportindex > 3 ???");
3269                 DD(m_SupportRegion.GetIndex(InputDimension-1));
3270               }
3271             } // end switch index
3272
3273           break;
3274         } // end case 2 shape
3275       case 3:
3276         {
3277           // ----------------------------------------------------------------------
3278           // The rabbit with 5 internal CP
3279           // Periodic, constrained to zero at the reference at position 3.5
3280           // and the extremes fixed with anti-symmetrical BC
3281           // Coeff      row  BC1  0   1  BC2  2   3  BC3 BC4
3282           // PaddedCoeff  R:  0   1   2   3   4   5   6   7
3283           // BC1= 2*R1-R2
3284           // BC2= -weights[3]/weights[1] ( R2+R5 ) - R4
3285           // BC3= R1
3286           // BC4= 2*R1-R5
3287           // ---------------------------------------------------------------------
3288           switch(m_SupportRegion.GetIndex(InputDimension-1))
3289             {
3290             case 0:
3291               {
3292                 // Lower end
3293                 m_FirstRegion.SetSize(InputDimension-1,0);
3294
3295                 // BC
3296                 m_BCRegions.resize(2);
3297                 m_BCValues.resize(2);
3298                 m_BCRegion.SetIndex(InputDimension-1,0);
3299                 m_BCRegions[0]=m_BCRegion;
3300                 m_BCValues[0]=2.;
3301                 m_BCRegion.SetIndex(InputDimension-1,1);
3302                 m_BCRegions[1]=m_BCRegion;
3303                 m_BCValues[1]= -1.;
3304
3305                 // Second part
3306                 m_SecondRegion.SetSize(InputDimension-1,2);
3307                 m_SecondRegion.SetIndex(InputDimension-1,0);
3308
3309                 // BC2
3310                 m_BC2Regions.resize(3);
3311                 m_BC2Values.resize(3);
3312                 m_BCRegion.SetIndex(InputDimension-1,1);
3313                 m_BC2Regions[0]=m_BCRegion;
3314                 m_BC2Values[0]=-m_WeightRatio[3][1];
3315                 m_BCRegion.SetIndex(InputDimension-1,2);
3316                 m_BC2Regions[1]=m_BCRegion;
3317                 m_BC2Values[1]= -1;
3318                 m_BCRegion.SetIndex(InputDimension-1,3);
3319                 m_BC2Regions[2]=m_BCRegion;
3320                 m_BC2Values[2]=-m_WeightRatio[3][1];
3321
3322                 break;
3323               }
3324             case 1:
3325               {
3326                 // Lower end
3327                 m_FirstRegion.SetSize(InputDimension-1,2);
3328                 m_FirstRegion.SetIndex(InputDimension-1,0);
3329
3330                 // BC
3331                 m_BCRegions.resize(3);
3332                 m_BCValues.resize(3);
3333                 m_BCRegion.SetIndex(InputDimension-1,1);
3334                 m_BCRegions[0]=m_BCRegion;
3335                 m_BCValues[0]=-m_WeightRatio[3][1];
3336                 m_BCRegion.SetIndex(InputDimension-1,2);
3337                 m_BCRegions[1]=m_BCRegion;
3338                 m_BCValues[1]= -1;
3339                 m_BCRegion.SetIndex(InputDimension-1,3);
3340                 m_BCRegions[2]=m_BCRegion;
3341                 m_BCValues[2]=-m_WeightRatio[3][1];
3342
3343                 // Second part
3344                 m_SecondRegion.SetSize(InputDimension-1,1);
3345                 m_SecondRegion.SetIndex(InputDimension-1,2);
3346
3347                 // BC2
3348                 m_BC2Regions.resize(0);
3349
3350                 break;
3351               }
3352             case 2:
3353               {
3354                 // Lower end
3355                 m_FirstRegion.SetSize(InputDimension-1,1);
3356                 m_FirstRegion.SetIndex(InputDimension-1,1);
3357
3358                 // BC
3359                 m_BCRegions.resize(3);
3360                 m_BCValues.resize(3);
3361                 m_BCRegion.SetIndex(InputDimension-1,1);
3362                 m_BCRegions[0]=m_BCRegion;
3363                 m_BCValues[0]=-m_WeightRatio[3][1];
3364                 m_BCRegion.SetIndex(InputDimension-1,2);
3365                 m_BCRegions[1]=m_BCRegion;
3366                 m_BCValues[1]= -1;
3367                 m_BCRegion.SetIndex(InputDimension-1,3);
3368                 m_BCRegions[2]=m_BCRegion;
3369                 m_BCValues[2]=-m_WeightRatio[3][1];
3370
3371                 // Second part
3372                 m_SecondRegion.SetSize(InputDimension-1,2);
3373                 m_SecondRegion.SetIndex(InputDimension-1,2);
3374
3375                 // BC2
3376                 m_BC2Regions.resize(0);
3377
3378                 break;
3379               }
3380             case 3:
3381               {
3382                 // Lower end
3383                 m_FirstRegion.SetSize(InputDimension-1,0);
3384
3385                 // BC
3386                 m_BCRegions.resize(3);
3387                 m_BCValues.resize(3);
3388                 m_BCRegion.SetIndex(InputDimension-1,1);
3389                 m_BCRegions[0]=m_BCRegion;
3390                 m_BCValues[0]=-m_WeightRatio[3][1];
3391                 m_BCRegion.SetIndex(InputDimension-1,2);
3392                 m_BCRegions[1]=m_BCRegion;
3393                 m_BCValues[1]= -1;
3394                 m_BCRegion.SetIndex(InputDimension-1,3);
3395                 m_BCRegions[2]=m_BCRegion;
3396                 m_BCValues[2]=-m_WeightRatio[3][1];
3397
3398                 // Second part
3399                 m_SecondRegion.SetSize(InputDimension-1,2);
3400                 m_SecondRegion.SetIndex(InputDimension-1,2);
3401
3402                 // BC2
3403                 m_BC2Regions.resize(1);
3404                 m_BC2Values.resize(1);
3405                 m_BCRegion.SetIndex(InputDimension-1,0);
3406                 m_BC2Regions[0]=m_BCRegion;
3407                 m_BC2Values[0]=1.;
3408
3409                 break;
3410               }
3411
3412             case 4:
3413               {
3414                 // Lower end
3415                 m_FirstRegion.SetSize(InputDimension-1,2);
3416                 m_FirstRegion.SetIndex(InputDimension-1,2);
3417
3418                 // BC2
3419                 m_BCRegions.resize(1);
3420                 m_BCValues.resize(1);
3421                 m_BCRegion.SetIndex(InputDimension-1,0);
3422                 m_BCRegions[0]=m_BCRegion;
3423                 m_BCValues[0]=1.;
3424
3425                 // Second part
3426                 m_SecondRegion.SetSize(InputDimension-1,0);
3427
3428                 // BC2
3429                 m_BC2Regions.resize(2);
3430                 m_BC2Values.resize(2);
3431                 m_BCRegion.SetIndex(InputDimension-1,0);
3432                 m_BC2Regions[0]=m_BCRegion;
3433                 m_BC2Values[0]=2;
3434                 m_BCRegion.SetIndex(InputDimension-1,3);
3435                 m_BC2Regions[1]=m_BCRegion;
3436                 m_BC2Values[1]=-1;
3437
3438                 break;
3439
3440               }
3441             default:
3442               {
3443                 DD("supportindex > 4 ???");
3444                 DD(m_SupportRegion.GetIndex(InputDimension-1));
3445               }
3446             } // end switch index
3447
3448           break;
3449         } // end case 3 shape
3450
3451       case 4:
3452         {
3453           // ----------------------------------------------------------------------
3454           // The sputnik with 4 internal CP
3455           // Periodic, constrained to zero at the reference
3456           // at position 3 and one indepent extremes copied
3457           // Coeff     row BC1  0   1  BC2  2  BC2  3
3458           // PaddedCoeff R: 0   1   2   3   4   5   6
3459           // BC1= R6
3460           // BC2= -weights[2]/weights[0] ( R2+R4 )
3461           // BC3=  weights[2]/weights[0] ( R2-R4) + R1
3462           // ---------------------------------------------------------------------
3463           switch(m_SupportRegion.GetIndex(InputDimension-1))
3464             {
3465             case 0:
3466               {
3467                 // Lower end
3468                 m_FirstRegion.SetSize(InputDimension-1,0);
3469
3470                 // BC
3471                 m_BCRegions.resize(1);
3472                 m_BCValues.resize(1);
3473                 m_BCRegion.SetIndex(InputDimension-1,3);
3474                 m_BCRegions[0]=m_BCRegion;
3475                 m_BCValues[0]=1;
3476
3477                 // Second part
3478                 m_SecondRegion.SetSize(InputDimension-1,2);
3479                 m_SecondRegion.SetIndex(InputDimension-1,0);
3480
3481                 // BC2
3482                 m_BC2Regions.resize(2);
3483                 m_BC2Values.resize(2);
3484                 m_BCRegion.SetIndex(InputDimension-1,1);
3485                 m_BC2Regions[0]=m_BCRegion;
3486                 m_BC2Values[0]=-m_WeightRatio[2][0];
3487                 m_BCRegion.SetIndex(InputDimension-1,2);
3488                 m_BC2Regions[1]=m_BCRegion;
3489                 m_BC2Values[1]= -m_WeightRatio[2][0];
3490
3491                 break;
3492               }
3493             case 1:
3494               {
3495                 // Lower end
3496                 m_FirstRegion.SetSize(InputDimension-1,2);
3497                 m_FirstRegion.SetIndex(InputDimension-1,0);
3498
3499                 // BC
3500                 m_BCRegions.resize(2);
3501                 m_BCValues.resize(2);
3502                 m_BCRegion.SetIndex(InputDimension-1,1);
3503                 m_BCRegions[0]=m_BCRegion;
3504                 m_BCValues[0]=-m_WeightRatio[2][0];
3505                 m_BCRegion.SetIndex(InputDimension-1,2);
3506                 m_BCRegions[1]=m_BCRegion;
3507                 m_BCValues[1]= -m_WeightRatio[2][0];
3508
3509                 // Second part
3510                 m_SecondRegion.SetSize(InputDimension-1,1);
3511                 m_SecondRegion.SetIndex(InputDimension-1,2);
3512
3513                 // BC2
3514                 m_BC2Regions.resize(0);
3515
3516                 break;
3517               }
3518             case 2:
3519               {
3520                 // Lower end
3521                 m_FirstRegion.SetSize(InputDimension-1,1);
3522                 m_FirstRegion.SetIndex(InputDimension-1,1);
3523
3524                 // BC
3525                 m_BCRegions.resize(2);
3526                 m_BCValues.resize(2);
3527                 m_BCRegion.SetIndex(InputDimension-1,1);
3528                 m_BCRegions[0]=m_BCRegion;
3529                 m_BCValues[0]=-m_WeightRatio[2][0];
3530                 m_BCRegion.SetIndex(InputDimension-1,2);
3531                 m_BCRegions[1]=m_BCRegion;
3532                 m_BCValues[1]= -m_WeightRatio[2][0];
3533
3534                 // Second part
3535                 m_SecondRegion.SetSize(InputDimension-1,1);
3536                 m_SecondRegion.SetIndex(InputDimension-1,2);
3537
3538                 // BC2
3539                 m_BC2Regions.resize(3);
3540                 m_BC2Values.resize(3);
3541                 m_BCRegion.SetIndex(InputDimension-1,0);
3542                 m_BC2Regions[0]=m_BCRegion;
3543                 m_BC2Values[0]=1;
3544                 m_BCRegion.SetIndex(InputDimension-1,1);
3545                 m_BC2Regions[1]=m_BCRegion;
3546                 m_BC2Values[1]=m_WeightRatio[2][0];
3547                 m_BCRegion.SetIndex(InputDimension-1,2);
3548                 m_BC2Regions[2]=m_BCRegion;
3549                 m_BC2Values[2]=-m_WeightRatio[2][0];
3550
3551                 break;
3552               }
3553             case 3:
3554               {
3555                 // Lower end
3556                 m_FirstRegion.SetSize(InputDimension-1,0);
3557
3558                 // BC
3559                 m_BCRegions.resize(2);
3560                 m_BCValues.resize(2);
3561                 m_BCRegion.SetIndex(InputDimension-1,1);
3562                 m_BCRegions[0]=m_BCRegion;
3563                 m_BCValues[0]=-m_WeightRatio[2][0];
3564                 m_BCRegion.SetIndex(InputDimension-1,2);
3565                 m_BCRegions[1]=m_BCRegion;
3566                 m_BCValues[1]= -m_WeightRatio[2][0];
3567
3568                 // Second part
3569                 m_SecondRegion.SetSize(InputDimension-1,1);
3570                 m_SecondRegion.SetIndex(InputDimension-1,2);
3571
3572                 // BC2
3573                 m_BC2Regions.resize(3);
3574                 m_BC2Values.resize(3);
3575                 m_BCRegion.SetIndex(InputDimension-1,0);
3576                 m_BC2Regions[0]=m_BCRegion;
3577                 m_BC2Values[0]=1;
3578                 m_BCRegion.SetIndex(InputDimension-1,1);
3579                 m_BC2Regions[1]=m_BCRegion;
3580                 m_BC2Values[1]=m_WeightRatio[2][0];
3581                 m_BCRegion.SetIndex(InputDimension-1,2);
3582                 m_BC2Regions[2]=m_BCRegion;
3583                 m_BC2Values[2]=-m_WeightRatio[2][0];
3584
3585                 // Third part
3586                 m_ThirdRegion.SetSize(InputDimension-1,1);
3587                 m_ThirdRegion.SetIndex(InputDimension-1,3);
3588
3589                 break;
3590               }
3591
3592             default:
3593               {
3594                 DD("supportindex > 3 ???");
3595                 DD(m_SupportRegion.GetIndex(InputDimension-1));
3596               }
3597             } // end switch index
3598
3599           break;
3600         } // end case 4 shape
3601
3602       case 5:
3603         {
3604           // ----------------------------------------------------------------------
3605           // The sputnik with 5 internal CP
3606           // Periodic, constrained to zero at the reference
3607           // at position 3 and one indepent extreme
3608           // Coeff     row BC1  0   1  BC2  2   3  BC3  4
3609           // PaddedCoeff R: 0   1   2   3   4   5   6   7
3610           // BC1= R2 + R5 - R7
3611           // BC2= -weights[3]/weights[1] ( R2+R5 ) - R4
3612           // BC3= R1 + 0.5 R2 - 0.5 R7
3613           // ----------------------------------------------------------------------
3614           switch(m_SupportRegion.GetIndex(InputDimension-1))
3615             {
3616             case 0:
3617               {
3618                 // Lower end
3619                 m_FirstRegion.SetSize(InputDimension-1,0);
3620
3621                 // BC
3622                 m_BCRegions.resize(3);
3623                 m_BCValues.resize(3);
3624                 m_BCRegion.SetIndex(InputDimension-1,1);
3625                 m_BCRegions[0]=m_BCRegion;
3626                 m_BCValues[0]=1.;
3627                 m_BCRegion.SetIndex(InputDimension-1,3);
3628                 m_BCRegions[1]=m_BCRegion;
3629                 m_BCValues[1]=1.;
3630                 m_BCRegion.SetIndex(InputDimension-1,4);
3631                 m_BCRegions[2]=m_BCRegion;
3632                 m_BCValues[2]=-1.;
3633
3634                 // Second part
3635                 m_SecondRegion.SetSize(InputDimension-1,2);
3636                 m_SecondRegion.SetIndex(InputDimension-1,0);
3637
3638                 // BC2
3639                 m_BC2Regions.resize(3);
3640                 m_BC2Values.resize(3);
3641                 m_BCRegion.SetIndex(InputDimension-1,1);
3642                 m_BC2Regions[0]=m_BCRegion;
3643                 m_BC2Values[0]=-m_WeightRatio[3][1];
3644                 m_BCRegion.SetIndex(InputDimension-1,2);
3645                 m_BC2Regions[1]=m_BCRegion;
3646                 m_BC2Values[1]= -1;
3647                 m_BCRegion.SetIndex(InputDimension-1,3);
3648                 m_BC2Regions[2]=m_BCRegion;
3649                 m_BC2Values[2]=-m_WeightRatio[3][1];
3650
3651                 break;
3652               }
3653             case 1:
3654               {
3655                 // Lower end
3656                 m_FirstRegion.SetSize(InputDimension-1,2);
3657                 m_FirstRegion.SetIndex(InputDimension-1,0);
3658
3659                 // BC
3660                 m_BCRegions.resize(3);
3661                 m_BCValues.resize(3);
3662                 m_BCRegion.SetIndex(InputDimension-1,1);
3663                 m_BCRegions[0]=m_BCRegion;
3664                 m_BCValues[0]=-m_WeightRatio[3][1];
3665                 m_BCRegion.SetIndex(InputDimension-1,2);
3666                 m_BCRegions[1]=m_BCRegion;
3667                 m_BCValues[1]= -1;
3668                 m_BCRegion.SetIndex(InputDimension-1,3);
3669                 m_BCRegions[2]=m_BCRegion;
3670                 m_BCValues[2]=-m_WeightRatio[3][1];
3671
3672                 // Second part
3673                 m_SecondRegion.SetSize(InputDimension-1,1);
3674                 m_SecondRegion.SetIndex(InputDimension-1,2);
3675
3676                 // BC2
3677                 m_BC2Regions.resize(0);
3678
3679                 break;
3680               }
3681             case 2:
3682               {
3683                 // Lower end
3684                 m_FirstRegion.SetSize(InputDimension-1,1);
3685                 m_FirstRegion.SetIndex(InputDimension-1,1);
3686
3687                 // BC
3688                 m_BCRegions.resize(3);
3689                 m_BCValues.resize(3);
3690                 m_BCRegion.SetIndex(InputDimension-1,1);
3691                 m_BCRegions[0]=m_BCRegion;
3692                 m_BCValues[0]=-m_WeightRatio[3][1];
3693                 m_BCRegion.SetIndex(InputDimension-1,2);
3694                 m_BCRegions[1]=m_BCRegion;
3695                 m_BCValues[1]= -1;
3696                 m_BCRegion.SetIndex(InputDimension-1,3);
3697                 m_BCRegions[2]=m_BCRegion;
3698                 m_BCValues[2]=-m_WeightRatio[3][1];
3699
3700                 // Second part
3701                 m_SecondRegion.SetSize(InputDimension-1,2);
3702                 m_SecondRegion.SetIndex(InputDimension-1,2);
3703
3704                 // BC2
3705                 m_BC2Regions.resize(0);
3706
3707                 break;
3708               }
3709             case 3:
3710               {
3711                 // Lower end
3712                 m_FirstRegion.SetSize(InputDimension-1,0);
3713
3714                 // BC
3715                 m_BCRegions.resize(3);
3716                 m_BCValues.resize(3);
3717                 m_BCRegion.SetIndex(InputDimension-1,1);
3718                 m_BCRegions[0]=m_BCRegion;
3719                 m_BCValues[0]=-m_WeightRatio[3][1];
3720                 m_BCRegion.SetIndex(InputDimension-1,2);
3721                 m_BCRegions[1]=m_BCRegion;
3722                 m_BCValues[1]= -1;
3723                 m_BCRegion.SetIndex(InputDimension-1,3);
3724                 m_BCRegions[2]=m_BCRegion;
3725                 m_BCValues[2]=-m_WeightRatio[3][1];
3726
3727                 // Second part
3728                 m_SecondRegion.SetSize(InputDimension-1,2);
3729                 m_SecondRegion.SetIndex(InputDimension-1,2);
3730
3731                 // BC2
3732                 m_BC2Regions.resize(3);
3733                 m_BC2Values.resize(3);
3734                 m_BCRegion.SetIndex(InputDimension-1,0);
3735                 m_BC2Regions[0]=m_BCRegion;
3736                 m_BC2Values[0]=1;
3737                 m_BCRegion.SetIndex(InputDimension-1,1);
3738                 m_BC2Regions[1]=m_BCRegion;
3739                 m_BC2Values[1]=0.5;
3740                 m_BCRegion.SetIndex(InputDimension-1,4);
3741                 m_BC2Regions[2]=m_BCRegion;
3742                 m_BC2Values[2]=-0.5;
3743
3744                 break;
3745               }
3746             case 4:
3747               {
3748                 // Lower end
3749                 m_FirstRegion.SetSize(InputDimension-1,2);
3750                 m_FirstRegion.SetIndex(InputDimension-1,2);
3751
3752                 // BC
3753                 m_BCRegions.resize(3);
3754                 m_BCValues.resize(3);
3755                 m_BCRegion.SetIndex(InputDimension-1,0);
3756                 m_BCRegions[0]=m_BCRegion;
3757                 m_BCValues[0]=1;
3758                 m_BCRegion.SetIndex(InputDimension-1,1);
3759                 m_BCRegions[1]=m_BCRegion;
3760                 m_BCValues[1]=0.5;
3761                 m_BCRegion.SetIndex(InputDimension-1,4);
3762                 m_BCRegions[2]=m_BCRegion;
3763                 m_BCValues[2]=-0.5;
3764
3765                 // Second part
3766                 m_SecondRegion.SetSize(InputDimension-1,1);
3767                 m_SecondRegion.SetIndex(InputDimension-1,4);
3768
3769                 // BC2
3770                 m_BC2Regions.resize(0);
3771
3772                 break;
3773               }
3774             default:
3775               {
3776                 DD("supportindex > 4 ???");
3777                 DD(m_SupportRegion.GetIndex(InputDimension-1));
3778               }
3779             } // end switch index
3780
3781           break;
3782         } // end case 5 shape
3783
3784
3785
3786
3787       case 6:
3788         {
3789           // ----------------------------------------------------------------------
3790           // The diamond with 4 internal CP:
3791           // Periodic, constrained to zero at the reference at position 3
3792           // Coeff      row 0   1   2  BC1  3  BC2  4
3793           // PaddedCoeff  R:0   1   2   3   4   5   6
3794           // BC1= -weights[2]/weights[0] ( R2+R4 )
3795           // BC2= weights[2]/weights[0]  ( R0+R2-R4-R6 ) + R1
3796           // ---------------------------------------------------------------------
3797           switch(m_SupportRegion.GetIndex(InputDimension-1))
3798             {
3799             case 0:
3800               {
3801                 // Lower end
3802                 m_FirstRegion.SetSize(InputDimension-1,3);
3803                 m_FirstRegion.SetIndex(InputDimension-1,0);
3804
3805                 // BC
3806                 m_BCRegions.resize(2);
3807                 m_BCValues.resize(2);
3808                 m_BCRegion.SetIndex(InputDimension-1,2);
3809                 m_BCRegions[0]=m_BCRegion;
3810                 m_BCValues[0]=-m_WeightRatio[2][0];
3811                 m_BCRegion.SetIndex(InputDimension-1,3);
3812                 m_BCRegions[1]=m_BCRegion;
3813                 m_BCValues[1]=-m_WeightRatio[2][0];
3814
3815                 // Second part
3816                 m_SecondRegion.SetSize(InputDimension-1,0);
3817
3818                 // BC2
3819                 m_BC2Regions.resize(0);
3820
3821                 break;
3822               }
3823             case 1:
3824               {
3825                 // Lower end
3826                 m_FirstRegion.SetSize(InputDimension-1,2);
3827                 m_FirstRegion.SetIndex(InputDimension-1,1);
3828
3829                 // BC
3830                 m_BCRegions.resize(2);
3831                 m_BCValues.resize(2);
3832                 m_BCRegion.SetIndex(InputDimension-1,2);
3833                 m_BCRegions[0]=m_BCRegion;
3834                 m_BCValues[0]=-m_WeightRatio[2][0];
3835                 m_BCRegion.SetIndex(InputDimension-1,3);
3836                 m_BCRegions[1]=m_BCRegion;
3837                 m_BCValues[1]=-m_WeightRatio[2][0];
3838
3839                 // Second part
3840                 m_SecondRegion.SetSize(InputDimension-1,1);
3841                 m_SecondRegion.SetIndex(InputDimension-1,3);
3842
3843                 // BC2
3844                 m_BC2Regions.resize(0);
3845
3846                 break;
3847               }
3848             case 2:
3849               {
3850                 // Lower end
3851                 m_FirstRegion.SetSize(InputDimension-1,1);
3852                 m_FirstRegion.SetIndex(InputDimension-1,2);
3853
3854                 // BC
3855                 m_BCRegions.resize(2);
3856                 m_BCValues.resize(2);
3857                 m_BCRegion.SetIndex(InputDimension-1,2);
3858                 m_BCRegions[0]=m_BCRegion;
3859                 m_BCValues[0]=-m_WeightRatio[2][0];
3860                 m_BCRegion.SetIndex(InputDimension-1,3);
3861                 m_BCRegions[1]=m_BCRegion;
3862                 m_BCValues[1]=-m_WeightRatio[2][0];
3863
3864                 // Second part
3865                 m_SecondRegion.SetSize(InputDimension-1,1);
3866                 m_SecondRegion.SetIndex(InputDimension-1,3);
3867
3868                 // BC2
3869                 m_BC2Regions.resize(5);
3870                 m_BC2Values.resize(5);
3871                 m_BCRegion.SetIndex(InputDimension-1,0);
3872                 m_BC2Regions[0]=m_BCRegion;
3873                 m_BC2Values[0]=m_WeightRatio[2][0];
3874                 m_BCRegion.SetIndex(InputDimension-1,1);
3875                 m_BC2Regions[1]=m_BCRegion;
3876                 m_BC2Values[1]=1;
3877                 m_BCRegion.SetIndex(InputDimension-1,2);
3878                 m_BC2Regions[2]=m_BCRegion;
3879                 m_BC2Values[2]=m_WeightRatio[2][0];
3880                 m_BCRegion.SetIndex(InputDimension-1,3);
3881                 m_BC2Regions[3]=m_BCRegion;
3882                 m_BC2Values[3]=-m_WeightRatio[2][0];
3883                 m_BCRegion.SetIndex(InputDimension-1,4);
3884                 m_BC2Regions[4]=m_BCRegion;
3885                 m_BC2Values[4]=-m_WeightRatio[2][0];
3886
3887                 break;
3888               }
3889             case 3:
3890               {
3891                 // Lower end
3892                 m_FirstRegion.SetSize(InputDimension-1,0);
3893
3894                 // BC
3895                 m_BCRegions.resize(2);
3896                 m_BCValues.resize(2);
3897                 m_BCRegion.SetIndex(InputDimension-1,2);
3898                 m_BCRegions[0]=m_BCRegion;
3899                 m_BCValues[0]=-m_WeightRatio[2][0];
3900                 m_BCRegion.SetIndex(InputDimension-1,3);
3901                 m_BCRegions[1]=m_BCRegion;
3902                 m_BCValues[1]=-m_WeightRatio[2][0];
3903
3904                 // Second part
3905                 m_SecondRegion.SetSize(InputDimension-1,1);
3906                 m_SecondRegion.SetIndex(InputDimension-1,3);
3907
3908                 // BC2
3909                 m_BC2Regions.resize(5);
3910                 m_BC2Values.resize(5);
3911                 m_BCRegion.SetIndex(InputDimension-1,0);
3912                 m_BC2Regions[0]=m_BCRegion;
3913                 m_BC2Values[0]=m_WeightRatio[2][0];
3914                 m_BCRegion.SetIndex(InputDimension-1,1);
3915                 m_BC2Regions[1]=m_BCRegion;
3916                 m_BC2Values[1]=1;
3917                 m_BCRegion.SetIndex(InputDimension-1,2);
3918                 m_BC2Regions[2]=m_BCRegion;
3919                 m_BC2Values[2]=m_WeightRatio[2][0];
3920                 m_BCRegion.SetIndex(InputDimension-1,3);
3921                 m_BC2Regions[3]=m_BCRegion;
3922                 m_BC2Values[3]=-m_WeightRatio[2][0];
3923                 m_BCRegion.SetIndex(InputDimension-1,4);
3924                 m_BC2Regions[4]=m_BCRegion;
3925                 m_BC2Values[4]=-m_WeightRatio[2][0];
3926
3927                 // Third part
3928                 m_ThirdRegion.SetSize(InputDimension-1,1);
3929                 m_ThirdRegion.SetIndex(InputDimension-1,4);
3930
3931                 break;
3932               }
3933
3934             default:
3935               {
3936                 DD("supportindex > 3 ???");
3937                 DD(m_SupportRegion.GetIndex(InputDimension-1));
3938               }
3939             } // end switch index
3940
3941           break;
3942         } // end case 7 shape
3943
3944       case 7:
3945         {
3946           // ----------------------------------------------------------------------
3947           // The diamond with 5 internal CP:
3948           // periodic, constrained to zero at the reference at position 3.5
3949           // Coeff      row 0   1   2  BC1  3   4  BC2  5
3950           // PaddedCoeff  R:0   1   2   3   4   5   6   7
3951           // BC1= -weights[3]/weights[1] ( R2+R5 ) - R4
3952           // BC2= weights[2]/weights[0] ( R0+R2-R5-R7 ) + R1
3953           // ---------------------------------------------------------------------
3954           switch(m_SupportRegion.GetIndex(InputDimension-1))
3955             {
3956             case 0:
3957               {
3958                 // Lower end
3959                 m_FirstRegion.SetSize(InputDimension-1,3);
3960                 m_FirstRegion.SetIndex(InputDimension-1,0);
3961
3962                 // BC
3963                 m_BCRegions.resize(3);
3964                 m_BCValues.resize(3);
3965                 m_BCRegion.SetIndex(InputDimension-1,2);
3966                 m_BCRegions[0]=m_BCRegion;
3967                 m_BCValues[0]=-m_WeightRatio[3][1];
3968                 m_BCRegion.SetIndex(InputDimension-1,3);
3969                 m_BCRegions[1]=m_BCRegion;
3970                 m_BCValues[1]= -1;
3971                 m_BCRegion.SetIndex(InputDimension-1,4);
3972                 m_BCRegions[2]=m_BCRegion;
3973                 m_BCValues[2]=-m_WeightRatio[3][1];
3974
3975                 // Second part
3976                 m_SecondRegion.SetSize(InputDimension-1,0);
3977
3978                 // BC2
3979                 m_BC2Regions.resize(0);
3980
3981                 break;
3982               }
3983             case 1:
3984               {
3985                 // Lower end
3986                 m_FirstRegion.SetSize(InputDimension-1,2);
3987                 m_FirstRegion.SetIndex(InputDimension-1,1);
3988
3989                 // BC
3990                 m_BCRegions.resize(3);
3991                 m_BCValues.resize(3);
3992                 m_BCRegion.SetIndex(InputDimension-1,2);
3993                 m_BCRegions[0]=m_BCRegion;
3994                 m_BCValues[0]=-m_WeightRatio[3][1];
3995                 m_BCRegion.SetIndex(InputDimension-1,3);
3996                 m_BCRegions[1]=m_BCRegion;
3997                 m_BCValues[1]= -1;
3998                 m_BCRegion.SetIndex(InputDimension-1,4);
3999                 m_BCRegions[2]=m_BCRegion;
4000                 m_BCValues[2]=-m_WeightRatio[3][1];
4001
4002                 // Second part
4003                 m_SecondRegion.SetSize(InputDimension-1,1);
4004                 m_SecondRegion.SetIndex(InputDimension-1,3);
4005
4006                 // BC2
4007                 m_BC2Regions.resize(0);
4008
4009                 break;
4010               }
4011             case 2:
4012               {
4013                 // Lower end
4014                 m_FirstRegion.SetSize(InputDimension-1,1);
4015                 m_FirstRegion.SetIndex(InputDimension-1,2);
4016
4017                 // BC
4018                 m_BCRegions.resize(3);
4019                 m_BCValues.resize(3);
4020                 m_BCRegion.SetIndex(InputDimension-1,2);
4021                 m_BCRegions[0]=m_BCRegion;
4022                 m_BCValues[0]=-m_WeightRatio[3][1];
4023                 m_BCRegion.SetIndex(InputDimension-1,3);
4024                 m_BCRegions[1]=m_BCRegion;
4025                 m_BCValues[1]= -1;
4026                 m_BCRegion.SetIndex(InputDimension-1,4);
4027                 m_BCRegions[2]=m_BCRegion;
4028                 m_BCValues[2]=-m_WeightRatio[3][1];
4029
4030                 // Second part
4031                 m_SecondRegion.SetSize(InputDimension-1,2);
4032                 m_SecondRegion.SetIndex(InputDimension-1,3);
4033
4034                 // BC2
4035                 m_BC2Regions.resize(0);
4036
4037                 break;
4038               }
4039             case 3:
4040               {
4041                 // Lower end
4042                 m_FirstRegion.SetSize(InputDimension-1,0);
4043                 m_FirstRegion.SetIndex(InputDimension-1,0);
4044
4045                 // BC
4046                 m_BCRegions.resize(3);
4047                 m_BCValues.resize(3);
4048                 m_BCRegion.SetIndex(InputDimension-1,2);
4049                 m_BCRegions[0]=m_BCRegion;
4050                 m_BCValues[0]=-m_WeightRatio[3][1];
4051                 m_BCRegion.SetIndex(InputDimension-1,3);
4052                 m_BCRegions[1]=m_BCRegion;
4053                 m_BCValues[1]= -1;
4054                 m_BCRegion.SetIndex(InputDimension-1,4);
4055                 m_BCRegions[2]=m_BCRegion;
4056                 m_BCValues[2]=-m_WeightRatio[3][1];
4057
4058                 // Second part
4059                 m_SecondRegion.SetSize(InputDimension-1,2);
4060                 m_SecondRegion.SetIndex(InputDimension-1,3);
4061
4062                 // BC2
4063                 m_BC2Regions.resize(5);
4064                 m_BC2Values.resize(5);
4065                 m_BCRegion.SetIndex(InputDimension-1,0);
4066                 m_BC2Regions[0]=m_BCRegion;
4067                 m_BC2Values[0]=m_WeightRatio[2][0];
4068                 m_BCRegion.SetIndex(InputDimension-1,1);
4069                 m_BC2Regions[1]=m_BCRegion;
4070                 m_BC2Values[1]=1;
4071                 m_BCRegion.SetIndex(InputDimension-1,2);
4072                 m_BC2Regions[2]=m_BCRegion;
4073                 m_BC2Values[2]=m_WeightRatio[2][0];
4074                 m_BCRegion.SetIndex(InputDimension-1,4);
4075                 m_BC2Regions[3]=m_BCRegion;
4076                 m_BC2Values[3]=-m_WeightRatio[2][0];
4077                 m_BCRegion.SetIndex(InputDimension-1,5);
4078                 m_BC2Regions[4]=m_BCRegion;
4079                 m_BC2Values[4]=-m_WeightRatio[2][0];
4080
4081                 break;
4082               }
4083
4084             case 4:
4085               {
4086                 // Lower end
4087                 m_FirstRegion.SetSize(InputDimension-1,2);
4088                 m_FirstRegion.SetIndex(InputDimension-1,3);
4089
4090                 // BC
4091                 m_BCRegions.resize(5);
4092                 m_BCValues.resize(5);
4093                 m_BCRegion.SetIndex(InputDimension-1,0);
4094                 m_BCRegions[0]=m_BCRegion;
4095                 m_BCValues[0]=m_WeightRatio[2][0];
4096                 m_BCRegion.SetIndex(InputDimension-1,1);
4097                 m_BCRegions[1]=m_BCRegion;
4098                 m_BCValues[1]=1;
4099                 m_BCRegion.SetIndex(InputDimension-1,2);
4100                 m_BCRegions[2]=m_BCRegion;
4101                 m_BCValues[2]=m_WeightRatio[2][0];
4102                 m_BCRegion.SetIndex(InputDimension-1,4);
4103                 m_BCRegions[3]=m_BCRegion;
4104                 m_BCValues[3]=-m_WeightRatio[2][0];
4105                 m_BCRegion.SetIndex(InputDimension-1,5);
4106                 m_BCRegions[4]=m_BCRegion;
4107                 m_BCValues[4]=-m_WeightRatio[2][0];
4108
4109                 // Second part
4110                 m_SecondRegion.SetSize(InputDimension-1,1);
4111                 m_SecondRegion.SetIndex(InputDimension-1,5);
4112
4113                 // BC2
4114                 m_BC2Regions.resize(0);
4115
4116                 break;
4117               }
4118
4119             default:
4120               {
4121                 DD("supportindex > 4 ???");
4122                 DD(m_SupportRegion.GetIndex(InputDimension-1));
4123               }
4124             } // end switch index
4125
4126           break;
4127         } // end case 7 shape
4128
4129       case 9:
4130         {
4131           // ----------------------------------------------------------------------
4132           // The sputnik with 5 internal CP
4133           // Periodic, constrained to zero at the reference
4134           // at position 3 and one indepent extreme
4135           // Coeff     row BC1  0   1  BC2  2   3  BC3  4
4136           // PaddedCoeff R: 0   1   2   3   4   5   6   7
4137           // BC1= -R2+R5-R7
4138           // BC2= -weights[3]/weights[1] ( R2+R5 ) - R4
4139           // BC3= R1
4140           // ---------------------------------------------------------------------
4141           switch(m_SupportRegion.GetIndex(InputDimension-1))
4142             {
4143             case 0:
4144               {
4145                 // Lower end
4146                 m_FirstRegion.SetSize(InputDimension-1,0);
4147
4148                 // BC
4149                 m_BCRegions.resize(3);
4150                 m_BCValues.resize(3);
4151                 m_BCRegion.SetIndex(InputDimension-1,1);
4152                 m_BCRegions[0]=m_BCRegion;
4153                 m_BCValues[0]=-1;
4154                 m_BCRegion.SetIndex(InputDimension-1,3);
4155                 m_BCRegions[1]=m_BCRegion;
4156                 m_BCValues[1]=1;
4157                 m_BCRegion.SetIndex(InputDimension-1,4);
4158                 m_BCRegions[2]=m_BCRegion;
4159                 m_BCValues[2]=-1;
4160
4161                 // Second part
4162                 m_SecondRegion.SetSize(InputDimension-1,2);
4163                 m_SecondRegion.SetIndex(InputDimension-1,0);
4164
4165                 // BC2
4166                 m_BC2Regions.resize(3);
4167                 m_BC2Values.resize(3);
4168                 m_BCRegion.SetIndex(InputDimension-1,1);
4169                 m_BC2Regions[0]=m_BCRegion;
4170                 m_BC2Values[0]=-m_WeightRatio[3][1];
4171                 m_BCRegion.SetIndex(InputDimension-1,2);
4172                 m_BC2Regions[1]=m_BCRegion;
4173                 m_BC2Values[1]= -1;
4174                 m_BCRegion.SetIndex(InputDimension-1,3);
4175                 m_BC2Regions[2]=m_BCRegion;
4176                 m_BC2Values[2]=-m_WeightRatio[3][1];
4177
4178                 break;
4179               }
4180             case 1:
4181               {
4182                 // Lower end
4183                 m_FirstRegion.SetSize(InputDimension-1,2);
4184                 m_FirstRegion.SetIndex(InputDimension-1,0);
4185
4186                 // BC
4187                 m_BCRegions.resize(3);
4188                 m_BCValues.resize(3);
4189                 m_BCRegion.SetIndex(InputDimension-1,1);
4190                 m_BCRegions[0]=m_BCRegion;
4191                 m_BCValues[0]=-m_WeightRatio[3][1];
4192                 m_BCRegion.SetIndex(InputDimension-1,2);
4193                 m_BCRegions[1]=m_BCRegion;
4194                 m_BCValues[1]= -1;
4195                 m_BCRegion.SetIndex(InputDimension-1,3);
4196                 m_BCRegions[2]=m_BCRegion;
4197                 m_BCValues[2]=-m_WeightRatio[3][1];
4198
4199                 // Second part
4200                 m_SecondRegion.SetSize(InputDimension-1,1);
4201                 m_SecondRegion.SetIndex(InputDimension-1,2);
4202
4203                 // BC2
4204                 m_BC2Regions.resize(0);
4205
4206                 break;
4207               }
4208             case 2:
4209               {
4210                 // Lower end
4211                 m_FirstRegion.SetSize(InputDimension-1,1);
4212                 m_FirstRegion.SetIndex(InputDimension-1,1);
4213
4214                 // BC
4215                 m_BCRegions.resize(3);
4216                 m_BCValues.resize(3);
4217                 m_BCRegion.SetIndex(InputDimension-1,1);
4218                 m_BCRegions[0]=m_BCRegion;
4219                 m_BCValues[0]=-m_WeightRatio[3][1];
4220                 m_BCRegion.SetIndex(InputDimension-1,2);
4221                 m_BCRegions[1]=m_BCRegion;
4222                 m_BCValues[1]= -1;
4223                 m_BCRegion.SetIndex(InputDimension-1,3);
4224                 m_BCRegions[2]=m_BCRegion;
4225                 m_BCValues[2]=-m_WeightRatio[3][1];
4226
4227                 // Second part
4228                 m_SecondRegion.SetSize(InputDimension-1,2);
4229                 m_SecondRegion.SetIndex(InputDimension-1,2);
4230
4231                 // BC2
4232                 m_BC2Regions.resize(0);
4233
4234                 break;
4235               }
4236             case 3:
4237               {
4238                 // Lower end
4239                 m_FirstRegion.SetSize(InputDimension-1,0);
4240
4241                 // BC
4242                 m_BCRegions.resize(3);
4243                 m_BCValues.resize(3);
4244                 m_BCRegion.SetIndex(InputDimension-1,1);
4245                 m_BCRegions[0]=m_BCRegion;
4246                 m_BCValues[0]=-m_WeightRatio[3][1];
4247                 m_BCRegion.SetIndex(InputDimension-1,2);
4248                 m_BCRegions[1]=m_BCRegion;
4249                 m_BCValues[1]= -1;
4250                 m_BCRegion.SetIndex(InputDimension-1,3);
4251                 m_BCRegions[2]=m_BCRegion;
4252                 m_BCValues[2]=-m_WeightRatio[3][1];
4253
4254                 // Second part
4255                 m_SecondRegion.SetSize(InputDimension-1,1);
4256                 m_SecondRegion.SetIndex(InputDimension-1,2);
4257
4258                 // BC
4259                 m_BC2Regions.resize(1);
4260                 m_BC2Values.resize(1);
4261                 m_BCRegion.SetIndex(InputDimension-1,0);
4262                 m_BC2Regions[0]=m_BCRegion;
4263                 m_BC2Values[0]=1;
4264
4265                 break;
4266               }
4267             case 4:
4268               {
4269                 // Lower end
4270                 m_FirstRegion.SetSize(InputDimension-1,2);
4271                 m_FirstRegion.SetIndex(InputDimension-1,2);
4272
4273                 // BC
4274                 m_BCRegions.resize(1);
4275                 m_BCValues.resize(1);
4276                 m_BCRegion.SetIndex(InputDimension-1,0);
4277                 m_BCRegions[0]=m_BCRegion;
4278                 m_BCValues[0]=1;
4279
4280                 // Second part
4281                 m_SecondRegion.SetSize(InputDimension-1,1);
4282                 m_SecondRegion.SetIndex(InputDimension-1,4);
4283
4284                 // BC2
4285                 m_BC2Regions.resize(0);
4286
4287                 break;
4288               }
4289             default:
4290               {
4291                 DD("supportindex > 4 ???");
4292                 DD(m_SupportRegion.GetIndex(InputDimension-1));
4293               }
4294             } // end switch index
4295
4296           break;
4297         } // end case 9 shape
4298
4299
4300       default:
4301         {
4302           DD ("Other shapes currently not implemented");
4303         }
4304
4305       } // end switch shape
4306   } // end wrap region
4307
4308
4309
4310   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
4311   void
4312   ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions,NOutputDimensions>
4313   ::TransformPointToContinuousIndex( const InputPointType & point, ContinuousIndexType & index ) const
4314   {
4315     unsigned int j;
4316
4317     itk::Vector<double, OutputDimension> tvector;
4318
4319     for ( j = 0; j < OutputDimension; j++ )
4320       {
4321         //JV find the index in the PADDED version
4322         tvector[j] = point[j] - this->m_PaddedGridOrigin[j];
4323       }
4324
4325     itk::Vector<double, OutputDimension> cvector;
4326
4327     cvector = m_PointToIndex * tvector;
4328
4329     for ( j = 0; j < OutputDimension; j++ )
4330       {
4331         index[j] = static_cast< typename ContinuousIndexType::CoordRepType >( cvector[j] );
4332       }
4333   }
4334
4335
4336 } // namespace
4337
4338 #endif
4339