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     typedef itk::MultiplyByConstantImageFilter<CoefficientImageType, double, CoefficientImageType> MultiplicationFilterType;
1041
1042     // Regions
1043     typename   CoefficientImageType::RegionType sourceRegion=m_PaddedCoefficientImage->GetLargestPossibleRegion();
1044     typename   CoefficientImageType::RegionType destinationRegion=m_PaddedCoefficientImage->GetLargestPossibleRegion();
1045     typename   CoefficientImageType::RegionType::SizeType sourceSize=sourceRegion.GetSize();
1046     typename   CoefficientImageType::RegionType::SizeType destinationSize=destinationRegion.GetSize();
1047     typename   CoefficientImageType::IndexType sourceIndex=sourceRegion.GetIndex();
1048     typename   CoefficientImageType::IndexType destinationIndex=destinationRegion.GetIndex();
1049
1050     // JV Padding depends on the shape
1051     switch (m_TransformShape)
1052       {
1053         /*  The shapes are
1054             0: egg     4 CP  3 DOF
1055             1: egg     5 CP  4 DOF
1056             2: rabbit  4 CP  3 DOF
1057             3: rabbit  5 CP  4 DOF
1058             4: sputnik 4 CP  4 DOF
1059             5: sputnik 5 CP  5 DOF
1060             6: diamond 6 CP  5 DOF
1061             7: diamond 7 CP  6 DOF
1062         */
1063
1064       case 0:
1065         {
1066           // ----------------------------------------------------------------------
1067           // The egg with 4 internal CP (starting from inhale)
1068           // Periodic, constrained to zero at the reference
1069           // at position 3 and
1070           // Coeff      row  BC1 BC2  0  BC3  1   2  BC4
1071           // PaddedCoeff  R:  0   1   2   3   4   5   6
1072           // BC1= R4
1073           // BC2= R5
1074           // BC3= -weights[2]/weights[0] ( R2+R4 )
1075           // BC4= R2
1076           // ---------------------------------------------------------------------
1077
1078           //---------------------------------
1079           // 1. First two temporal rows are identical: paste 1-2 to 0-1
1080           typename PasteImageFilterType::Pointer paster0=PasteImageFilterType::New();
1081           paster0->SetDestinationImage(m_PaddedCoefficientImage);
1082           paster0->SetSourceImage(m_CoefficientImage);
1083           sourceRegion.SetIndex(NInputDimensions-1,1);
1084           sourceRegion.SetSize(NInputDimensions-1,2);
1085           destinationIndex[InputDimension-1]=0;
1086           paster0->SetDestinationIndex(destinationIndex);
1087           paster0->SetSourceRegion(sourceRegion);
1088
1089           //---------------------------------
1090           // 1. Next temporal row = paste 0 to 2
1091           typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1092           typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1093           paster1->SetDestinationImage(paster0->GetOutput());
1094           paster1->SetSourceImage(row0);
1095           destinationIndex[InputDimension-1]=2;
1096           paster1->SetDestinationIndex(destinationIndex);
1097           paster1->SetSourceRegion(row0->GetLargestPossibleRegion());
1098
1099           //---------------------------------
1100           // 2. Middle row at index 3=BC
1101           // Extract row at index 1, 2 (of coeff image)
1102           typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1103
1104           // Combine
1105           typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1106           combine1->SetInput(0,row0);
1107           combine1->SetInput(1,row1);
1108           combine1->SetA(-m_WeightRatio[2][0]);
1109           combine1->SetB(-m_WeightRatio[2][0]);
1110           combine1->Update();
1111           typename CoefficientImageType::Pointer bc3Row=combine1->GetOutput();
1112
1113           // Paste middleRow at index 3 (padded coeff)
1114           typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1115           paster2->SetDestinationImage(paster1->GetOutput());
1116           paster2->SetSourceImage(bc3Row);
1117           destinationIndex[InputDimension-1]=3;
1118           paster2->SetDestinationIndex(destinationIndex);
1119           paster2->SetSourceRegion(bc3Row->GetLargestPossibleRegion());
1120
1121           //---------------------------------
1122           // 3. Next two temporal rows identical: paste 1,2 to 4,5
1123           typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1124           paster3->SetDestinationImage(paster2->GetOutput());
1125           paster3->SetSourceImage(m_CoefficientImage);
1126           sourceRegion.SetIndex(InputDimension-1,1);
1127           sourceRegion.SetSize(InputDimension-1,2);
1128           destinationIndex[InputDimension-1]=4;
1129           paster3->SetDestinationIndex(destinationIndex);
1130           paster3->SetSourceRegion(sourceRegion);
1131
1132           // ---------------------------------
1133           // 4. Row at index 6=BC4= R2
1134           // Paste BC3 row at index 5
1135           typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
1136           paster4->SetDestinationImage(paster3->GetOutput());
1137           paster4->SetSourceImage(row0);
1138           destinationIndex[InputDimension-1]=6;
1139           paster4->SetDestinationIndex(destinationIndex);
1140           paster4->SetSourceRegion(row0->GetLargestPossibleRegion());
1141
1142           // Update the chain!
1143           paster4->Update();
1144           m_PaddedCoefficientImage= paster4->GetOutput();
1145
1146           break;
1147         }
1148
1149       case 1:
1150         {
1151           // ----------------------------------------------------------------------
1152           // The egg with 5 internal CP (starting from inhale)
1153           // Periodic, constrained to zero at the reference
1154           // at position 3 and
1155           // Coeff      row  BC1 BC2  0  BC3  1   2   3  BC4
1156           // PaddedCoeff  R:  0   1   2   3   4   5   6   7
1157           // BC1= R5
1158           // BC2= R6
1159           // BC3= -weights[3]/weights[1] ( R2+R5 ) - R4
1160           // BC4= R2
1161           // ---------------------------------------------------------------------
1162           //---------------------------------
1163           // 1. First two temporal rows are identical: paste 2-3 to 0-1
1164           typename PasteImageFilterType::Pointer paster0=PasteImageFilterType::New();
1165           paster0->SetDestinationImage(m_PaddedCoefficientImage);
1166           paster0->SetSourceImage(m_CoefficientImage);
1167           sourceRegion.SetIndex(NInputDimensions-1,2);
1168           sourceRegion.SetSize(NInputDimensions-1,2);
1169           destinationIndex[InputDimension-1]=0;
1170           paster0->SetDestinationIndex(destinationIndex);
1171           paster0->SetSourceRegion(sourceRegion);
1172
1173           //---------------------------------
1174           // 1. Next temporal row = paste 0 to 2
1175           typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1176           typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1177           paster1->SetDestinationImage(paster0->GetOutput());
1178           paster1->SetSourceImage(row0);
1179           destinationIndex[InputDimension-1]=2;
1180           paster1->SetDestinationIndex(destinationIndex);
1181           paster1->SetSourceRegion(row0->GetLargestPossibleRegion());
1182
1183           //---------------------------------
1184           // 2. Middle row at index 3=BC
1185           // Extract row at index 1, 2 (of coeff image)
1186           typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1187           typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
1188
1189           // Combine
1190           typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1191           combine1->SetInput(0,row0);
1192           combine1->SetInput(1,row2);
1193           combine1->SetA(1);
1194           combine1->SetB(1);
1195           // combine1->Update();
1196           typename LinearCombinationFilterType::Pointer combine2=LinearCombinationFilterType::New();
1197           combine2->SetInput(0,row1);
1198           combine2->SetInput(1,combine1->GetOutput());
1199           combine2->SetA(-1.);
1200           combine2->SetB(-m_WeightRatio[3][1]);
1201           combine2->Update();
1202           typename CoefficientImageType::Pointer bc3Row=combine2->GetOutput();
1203
1204           // Paste middleRow at index 3 (padded coeff)
1205           typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1206           paster2->SetDestinationImage(paster1->GetOutput());
1207           paster2->SetSourceImage(bc3Row);
1208           destinationIndex[InputDimension-1]=3;
1209           paster2->SetDestinationIndex(destinationIndex);
1210           paster2->SetSourceRegion(bc3Row->GetLargestPossibleRegion());
1211
1212           //---------------------------------
1213           // 3. Next three temporal rows identical: paste 1,2,3 to 4,5,6
1214           typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1215           paster3->SetDestinationImage(paster2->GetOutput());
1216           paster3->SetSourceImage(m_CoefficientImage);
1217           sourceRegion.SetIndex(InputDimension-1,1);
1218           sourceRegion.SetSize(InputDimension-1,3);
1219           destinationIndex[InputDimension-1]=4;
1220           paster3->SetDestinationIndex(destinationIndex);
1221           paster3->SetSourceRegion(sourceRegion);
1222
1223           // ---------------------------------
1224           // 4. Row at index 7=BC4= R2
1225           // Paste BC3 row at index 5
1226           typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
1227           paster4->SetDestinationImage(paster3->GetOutput());
1228           paster4->SetSourceImage(row0);
1229           destinationIndex[InputDimension-1]=7;
1230           paster4->SetDestinationIndex(destinationIndex);
1231           paster4->SetSourceRegion(row0->GetLargestPossibleRegion());
1232
1233           // Update the chain!
1234           paster4->Update();
1235           m_PaddedCoefficientImage= paster4->GetOutput();
1236
1237           break;
1238         }
1239
1240
1241 //        // ----------------------------------------------------------------------
1242 //        // The egg with 5 internal CP:
1243 //        // Periodic, C2 smooth everywhere and constrained to zero at the reference
1244 //        // Coeff       row R5 BC1  0   1   2   3  BC2  R2
1245 //        // PaddedCoeff R:  0   1   2   3   4   5   6   7
1246 //        // BC1= -weights[2]/weights[0] ( R2+R5)
1247 //        // BC2= BC1
1248 //        // ---------------------------------------------------------------------
1249
1250 //        // Extract rows with index 0 and 3
1251 //        typename    CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1252 //        typename    CoefficientImageType::Pointer row3=this->ExtractTemporalRow(m_CoefficientImage, 3);
1253
1254 //        // Paste the first row
1255 //        typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1256 //        destinationIndex.Fill(0);
1257 //        paster1->SetDestinationIndex(destinationIndex);
1258 //        paster1->SetSourceRegion(row3->GetLargestPossibleRegion());
1259 //        paster1->SetSourceImage(row3);
1260 //        paster1->SetDestinationImage(m_PaddedCoefficientImage);
1261
1262 //        // Linearly Combine rows for BC1 and BC2
1263 //        typedef clitk::LinearCombinationImageFilter<CoefficientImageType, CoefficientImageType> LinearCombinationFilterType;
1264 //        typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1265 //        combine1->SetFirstInput(row0);
1266 //        combine1->SetSecondInput(row3);
1267 //        combine1->SetA(-m_WeightRatio[2][0]);
1268 //        combine1->SetB(-m_WeightRatio[2][0]);
1269 //        combine1->Update();
1270 //        typename CoefficientImageType::Pointer bcRow=combine1->GetOutput();
1271
1272 //        // Paste the second row
1273 //        typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1274 //        destinationIndex[InputDimension-1]=1;
1275 //        paster2->SetDestinationIndex(destinationIndex);
1276 //        paster2->SetSourceRegion(bcRow->GetLargestPossibleRegion());
1277 //        paster2->SetSourceImage(bcRow);
1278 //        paster2->SetDestinationImage(paster1->GetOutput());
1279
1280 //        // Paste the coefficientImage
1281 //        typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1282 //        destinationIndex[InputDimension-1]=2;
1283 //        paster3->SetDestinationIndex(destinationIndex);
1284 //        paster3->SetSourceRegion(m_CoefficientImage->GetLargestPossibleRegion());
1285 //        paster3->SetSourceImage(m_CoefficientImage);
1286 //        paster3->SetDestinationImage(paster2->GetOutput());
1287
1288 //        // Paste the last two rows
1289 //        typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
1290 //        destinationIndex[InputDimension-1]=6;
1291 //        paster4->SetDestinationIndex(destinationIndex);
1292 //        paster4->SetSourceRegion(bcRow->GetLargestPossibleRegion());
1293 //        paster4->SetSourceImage(bcRow);
1294 //        paster4->SetDestinationImage(paster3->GetOutput());
1295
1296 //        typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
1297 //        destinationIndex[InputDimension-1]=7;
1298 //        paster5->SetDestinationIndex(destinationIndex);
1299 //        paster5->SetSourceRegion(row0->GetLargestPossibleRegion());
1300 //        paster5->SetSourceImage(row0);
1301 //        paster5->SetDestinationImage(paster4->GetOutput());
1302
1303 //        // Update the chain!
1304 //        paster5->Update();
1305 //        m_PaddedCoefficientImage= paster5->GetOutput();
1306
1307 //        break;
1308 //      }
1309
1310       case 2:
1311         {
1312           // ----------------------------------------------------------------------
1313           // The rabbit with 4 internal CP
1314           // Periodic, constrained to zero at the reference
1315           // at position 3 and the extremes fixed with anit-symm bc
1316           // Coeff      row  BC1  0   1  BC2  2  BC3 BC4
1317           // PaddedCoeff  R:  0   1   2   3   4   5   6
1318           // BC1= 2*R1-R2
1319           // BC2= -weights[2]/weights[0] ( R2+R4 )
1320           // BC3=  R1
1321           // BC4= 2*R1-R4
1322           // ---------------------------------------------------------------------
1323
1324           // ---------------------------------
1325           // 0. First Row =BC1
1326           typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1327           typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1328           typename LinearCombinationFilterType::Pointer combine0=LinearCombinationFilterType::New();
1329           combine0->SetInput(0,row0);
1330           combine0->SetInput(1,row1);
1331           combine0->SetA(2.);
1332           combine0->SetB(-1.);
1333           combine0->Update();
1334           typename CoefficientImageType::Pointer bc1Row=combine0->GetOutput();
1335
1336           typename PasteImageFilterType::Pointer paster0=PasteImageFilterType::New();
1337           paster0->SetDestinationImage(m_PaddedCoefficientImage);
1338           paster0->SetSourceImage(bc1Row);
1339           destinationIndex[InputDimension-1]=0;
1340           paster0->SetDestinationIndex(destinationIndex);
1341           paster0->SetSourceRegion(bc1Row->GetLargestPossibleRegion());
1342
1343           //---------------------------------
1344           // 1. Next two temporal rows are identical: paste 0-1 to 1-2
1345           typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1346           paster1->SetDestinationImage(paster0->GetOutput());
1347           paster1->SetSourceImage(m_CoefficientImage);
1348           sourceRegion.SetIndex(NInputDimensions-1,0);
1349           destinationIndex[InputDimension-1]=1;
1350           sourceRegion.SetSize(NInputDimensions-1,2);
1351           paster1->SetDestinationIndex(destinationIndex);
1352           paster1->SetSourceRegion(sourceRegion);
1353
1354           //---------------------------------
1355           // 2. Middle row at index 3=BC
1356           // Extract row at index 1, 2 (of coeff image)
1357           typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
1358
1359           // Combine
1360           typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1361           combine1->SetInput(0,row1);
1362           combine1->SetInput(1,row2);
1363           combine1->SetA(-m_WeightRatio[2][0]);
1364           combine1->SetB(-m_WeightRatio[2][0]);
1365           combine1->Update();
1366           typename CoefficientImageType::Pointer bc2Row=combine1->GetOutput();
1367
1368           // Paste middleRow at index 3 (padded coeff)
1369           typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1370           paster2->SetDestinationImage(paster1->GetOutput());
1371           paster2->SetSourceImage(bc2Row);
1372           destinationIndex[InputDimension-1]=3;
1373           paster2->SetDestinationIndex(destinationIndex);
1374           paster2->SetSourceRegion(bc2Row->GetLargestPossibleRegion());
1375
1376           //---------------------------------
1377           // 3. Next temporal row is identical: paste 2 to 4
1378           typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1379           paster3->SetDestinationImage(paster2->GetOutput());
1380           paster3->SetSourceImage(row2);
1381           destinationIndex[InputDimension-1]=4;
1382           paster3->SetDestinationIndex(destinationIndex);
1383           paster3->SetSourceRegion(row2->GetLargestPossibleRegion());
1384
1385           // ---------------------------------
1386           // 4. Row at index 5=BC (paddedcoeff image) R1
1387           // Paste BC3 row at index 5
1388           typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
1389           paster4->SetDestinationImage(paster3->GetOutput());
1390           paster4->SetSourceImage(row0);
1391           destinationIndex[InputDimension-1]=5;
1392           paster4->SetDestinationIndex(destinationIndex);
1393           paster4->SetSourceRegion(row0->GetLargestPossibleRegion());
1394
1395           // ---------------------------------
1396           // 5. Paste BC4 row at index 6
1397           typename LinearCombinationFilterType::Pointer combine3=LinearCombinationFilterType::New();
1398           combine3->SetInput(0,row0);
1399           combine3->SetInput(1,row2);
1400           combine3->SetA(2.);
1401           combine3->SetB(-1.);
1402           combine3->Update();
1403           typename CoefficientImageType::Pointer bc4Row=combine3->GetOutput();
1404           typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
1405           paster5->SetDestinationImage(paster4->GetOutput());
1406           paster5->SetSourceImage(bc4Row);
1407           destinationIndex[InputDimension-1]=6;
1408           paster5->SetDestinationIndex(destinationIndex);
1409           paster5->SetSourceRegion(bc4Row->GetLargestPossibleRegion());
1410
1411           // Update the chain!
1412           paster5->Update();
1413           m_PaddedCoefficientImage= paster5->GetOutput();
1414
1415           break;
1416         }
1417
1418       case 3:
1419         {
1420           // ----------------------------------------------------------------------
1421           // The rabbit with 5 internal CP
1422           // Periodic, constrained to zero at the reference at position 3.5
1423           // and the extremes fixed with anti-symmetrical BC
1424           // Coeff      row  BC1  0   1  BC2  2   3  BC3 BC4
1425           // PaddedCoeff  R:  0   1   2   3   4   5   6   7
1426           // BC1= 2*R1-R2
1427           // BC2= -weights[3]/weights[1] ( R2+R5 ) - R4
1428           // BC3= R1
1429           // BC4= 2*R1-R5
1430           // ---------------------------------------------------------------------
1431
1432           // ---------------------------------
1433           // 0. First Row =BC1
1434           typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1435           typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1436           typename LinearCombinationFilterType::Pointer combine0=LinearCombinationFilterType::New();
1437           combine0->SetInput(0,row0);
1438           combine0->SetInput(1,row1);
1439           combine0->SetA(2.);
1440           combine0->SetB(-1.);
1441           combine0->Update();
1442           typename CoefficientImageType::Pointer bc1Row=combine0->GetOutput();
1443
1444           typename PasteImageFilterType::Pointer paster0=PasteImageFilterType::New();
1445           paster0->SetDestinationImage(m_PaddedCoefficientImage);
1446           paster0->SetSourceImage(bc1Row);
1447           destinationIndex[InputDimension-1]=0;
1448           paster0->SetDestinationIndex(destinationIndex);
1449           paster0->SetSourceRegion(bc1Row->GetLargestPossibleRegion());
1450
1451           //---------------------------------
1452           // 1. Next two temporal rows are identical: paste 0-1 to 1-2
1453           typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1454           paster1->SetDestinationImage(paster0->GetOutput());
1455           paster1->SetSourceImage(m_CoefficientImage);
1456           sourceRegion.SetIndex(InputDimension-1,0);
1457           destinationIndex[InputDimension-1]=1;
1458           sourceRegion.SetSize(NInputDimensions-1,2);
1459           paster1->SetDestinationIndex(destinationIndex);
1460           paster1->SetSourceRegion(sourceRegion);
1461
1462           //---------------------------------
1463           // 2. Middle row at index 3=BC
1464           // Extract row at index 1, 3 (of coeff image)
1465           //typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1466           typename CoefficientImageType::Pointer row3=this->ExtractTemporalRow(m_CoefficientImage, 3);
1467
1468           // Combine
1469           typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1470           combine1->SetInput(0,row1);
1471           combine1->SetInput(1,row3);
1472           combine1->SetA(1.);
1473           combine1->SetB(1.);
1474           combine1->Update();
1475
1476           // Extract row at index 2 (coeff image)
1477           typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
1478
1479           // Combine
1480           typename LinearCombinationFilterType::Pointer combine2=LinearCombinationFilterType::New();
1481           combine2->SetInput(0,combine1->GetOutput());
1482           combine2->SetInput(1,row2);
1483           combine2->SetA(-m_WeightRatio[3][1]);
1484           combine2->SetB(-1.);
1485           combine2->Update();
1486           typename CoefficientImageType::Pointer bc2Row=combine2->GetOutput();
1487
1488           // Paste middleRow at index 3 (padded coeff)
1489           typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1490           paster2->SetDestinationImage(paster1->GetOutput());
1491           paster2->SetSourceImage(bc2Row);
1492           destinationIndex[InputDimension-1]=3;
1493           paster2->SetDestinationIndex(destinationIndex);
1494           paster2->SetSourceRegion(bc2Row->GetLargestPossibleRegion());
1495
1496           //---------------------------------
1497           // 3. Next two temporal rows are identical: paste 2,3 to 4,5
1498           typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1499           paster3->SetDestinationImage(paster2->GetOutput());
1500           paster3->SetSourceImage(m_CoefficientImage);
1501           sourceRegion.SetIndex(InputDimension-1,2);
1502           destinationIndex[InputDimension-1]=4;
1503           sourceRegion.SetSize(NInputDimensions-1,2);
1504           paster3->SetDestinationIndex(destinationIndex);
1505           paster3->SetSourceRegion(sourceRegion);
1506
1507           // ---------------------------------
1508           // 4. Row at index 6=BC (paddedcoeff image)R1
1509           // Paste BC3 row at index 6
1510           typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
1511           paster4->SetDestinationImage(paster3->GetOutput());
1512           paster4->SetSourceImage(row0);
1513           destinationIndex[InputDimension-1]=6;
1514           paster4->SetDestinationIndex(destinationIndex);
1515           paster4->SetSourceRegion(row0->GetLargestPossibleRegion());
1516
1517           // ---------------------------------
1518           // 5. Paste BC4 row at index 7
1519           typename LinearCombinationFilterType::Pointer combine3=LinearCombinationFilterType::New();
1520           combine3->SetInput(0,row0);
1521           combine3->SetInput(1,row3);
1522           combine3->SetA(2.);
1523           combine3->SetB(-1.);
1524           combine3->Update();
1525           typename CoefficientImageType::Pointer bc4Row=combine3->GetOutput();
1526           typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
1527           paster5->SetDestinationImage(paster4->GetOutput());
1528           paster5->SetSourceImage(bc4Row);
1529           destinationIndex[InputDimension-1]=7;
1530           paster5->SetDestinationIndex(destinationIndex);
1531           paster5->SetSourceRegion(bc4Row->GetLargestPossibleRegion());
1532
1533           // Update the chain!
1534           paster5->Update();
1535           m_PaddedCoefficientImage= paster5->GetOutput();
1536
1537           break;
1538         }
1539
1540       case 4:
1541         {
1542           // ----------------------------------------------------------------------
1543           // The sputnik with 4 internal CP
1544           // Periodic, constrained to zero at the reference
1545           // at position 3 and one indepent extremes copied
1546           // Coeff     row BC1  0   1  BC2  2  BC2  3
1547           // PaddedCoeff R: 0   1   2   3   4   5   6
1548           // BC1= R6
1549           // BC2= -weights[2]/weights[0] ( R2+R4 )
1550           // BC3=  weights[2]/weights[0] ( R2-R4) + R1
1551           // ---------------------------------------------------------------------
1552
1553           //---------------------------------
1554           // 1. First Row is equal to last row: paste 3 row to 0
1555           typename CoefficientImageType::Pointer row3=this->ExtractTemporalRow(m_CoefficientImage, 3);
1556           typename PasteImageFilterType::Pointer paster0=PasteImageFilterType::New();
1557           paster0->SetDestinationImage(m_PaddedCoefficientImage);
1558           paster0->SetSourceImage(row3);
1559           destinationIndex[InputDimension-1]=0;
1560           paster0->SetDestinationIndex(destinationIndex);
1561           paster0->SetSourceRegion(row3->GetLargestPossibleRegion());
1562
1563           //---------------------------------
1564           // 1. Next two temporal rows are identical: paste 0-1 to 1-2
1565           typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1566           paster1->SetDestinationImage(paster0->GetOutput());
1567           paster1->SetSourceImage(m_CoefficientImage);
1568           sourceRegion.SetIndex(NInputDimensions-1,0);
1569           destinationIndex[InputDimension-1]=1;
1570           sourceRegion.SetSize(NInputDimensions-1,2);
1571           paster1->SetDestinationIndex(destinationIndex);
1572           paster1->SetSourceRegion(sourceRegion);
1573
1574           //---------------------------------
1575           // 2. Middle row at index 3=BC
1576           // Extract row at index 1, 2 (of coeff image)
1577           typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1578           typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
1579
1580           // Combine
1581           typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1582           combine1->SetInput(0,row1);
1583           combine1->SetInput(1,row2);
1584           combine1->SetA(-m_WeightRatio[2][0]);
1585           combine1->SetB(-m_WeightRatio[2][0]);
1586           combine1->Update();
1587           typename CoefficientImageType::Pointer bc1Row=combine1->GetOutput();
1588
1589           // Paste middleRow at index 3 (padded coeff)
1590           typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1591           paster2->SetDestinationImage(paster1->GetOutput());
1592           paster2->SetSourceImage(bc1Row);
1593           destinationIndex[InputDimension-1]=3;
1594           paster2->SetDestinationIndex(destinationIndex);
1595           paster2->SetSourceRegion(bc1Row->GetLargestPossibleRegion());
1596
1597           //---------------------------------
1598           // 3. Next temporal row is identical: paste 2 to 4
1599           typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1600           paster3->SetDestinationImage(paster2->GetOutput());
1601           paster3->SetSourceImage(row2);
1602           destinationIndex[InputDimension-1]=4;
1603           paster3->SetDestinationIndex(destinationIndex);
1604           paster3->SetSourceRegion(row2->GetLargestPossibleRegion());
1605
1606           //---------------------------------
1607           // 4. Final row at index 5=BC3 (paddedcoeff image)R1 R2 R4 corresponds to index in coeff image 0 1 2
1608           // Extract row at index 1, 2 (of coeff image)already done
1609           typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1610
1611           // Combine
1612           typename LinearCombinationFilterType::Pointer combine3=LinearCombinationFilterType::New();
1613           combine3->SetInput(0,row1);
1614           combine3->SetInput(1,row2);
1615           combine3->SetA(1.);
1616           combine3->SetB(-1.);
1617
1618           // Combine
1619           typename LinearCombinationFilterType::Pointer combine4=LinearCombinationFilterType::New();
1620           combine4->SetInput(0,combine3->GetOutput());
1621           combine4->SetInput(1,row0);
1622           combine4->SetA(m_WeightRatio[2][0]);
1623           combine4->SetB(1.);
1624           combine4->Update();
1625           typename CoefficientImageType::Pointer bc2Row=combine4->GetOutput();
1626
1627           // Paste BC row at index 5
1628           typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
1629           paster4->SetDestinationImage(paster3->GetOutput());
1630           paster4->SetSourceImage(bc2Row);
1631           destinationIndex[InputDimension-1]=5;
1632           paster4->SetDestinationIndex(destinationIndex);
1633           paster4->SetSourceRegion(bc2Row->GetLargestPossibleRegion());
1634
1635           // Paste row 3  at index 6
1636           typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
1637           paster5->SetDestinationImage(paster4->GetOutput());
1638           paster5->SetSourceImage(row3);
1639           destinationIndex[InputDimension-1]=6;
1640           paster5->SetDestinationIndex(destinationIndex);
1641           paster5->SetSourceRegion(row3->GetLargestPossibleRegion());
1642
1643           // Update the chain!
1644           paster5->Update();
1645           m_PaddedCoefficientImage= paster5->GetOutput();
1646
1647           break;
1648         }
1649
1650       case 5:
1651         {
1652
1653           // ----------------------------------------------------------------------
1654           // The sputnik with 5 internal CP
1655           // Periodic, constrained to zero at the reference
1656           // at position 3 and one indepent extreme
1657           // Coeff     row BC1  0   1  BC2  2   3  BC3  4
1658           // PaddedCoeff R: 0   1   2   3   4   5   6   7
1659           // BC1= R2 + R5 - R7
1660           // BC2= -weights[3]/weights[1] ( R2+R5 ) - R4
1661           // BC3= R1 + 0.5 R2 - 0.5 R7
1662           // ----------------------------------------------------------------------
1663           //---------------------------------
1664           // 1. First Row =BC
1665           typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1666           typename CoefficientImageType::Pointer row3=this->ExtractTemporalRow(m_CoefficientImage, 3);
1667           typename CoefficientImageType::Pointer row4=this->ExtractTemporalRow(m_CoefficientImage, 4);
1668
1669           // Combine
1670           typename LinearCombinationFilterType::Pointer combine0=LinearCombinationFilterType::New();
1671           combine0->SetInput(0,row1);
1672           combine0->SetInput(1,row3);
1673           combine0->SetA(1.);
1674           combine0->SetB(1.);
1675           //combine0->Update();
1676           typename LinearCombinationFilterType::Pointer combine0bis=LinearCombinationFilterType::New();
1677           combine0bis->SetInput(0,combine0->GetOutput());
1678           combine0bis->SetInput(1,row4);
1679           combine0bis->SetA(1.);
1680           combine0bis->SetB(-1.);
1681           combine0bis->Update();
1682           typename CoefficientImageType::Pointer bc1Row=combine0bis->GetOutput();
1683
1684           typename PasteImageFilterType::Pointer paster0=PasteImageFilterType::New();
1685           paster0->SetDestinationImage(m_PaddedCoefficientImage);
1686           paster0->SetSourceImage(bc1Row);
1687           destinationIndex[InputDimension-1]=0;
1688           paster0->SetDestinationIndex(destinationIndex);
1689           paster0->SetSourceRegion(bc1Row->GetLargestPossibleRegion());
1690
1691           //---------------------------------
1692           // 1. Next two temporal rows are identical: paste 0-1 to 1-2
1693           typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1694           paster1->SetDestinationImage(paster0->GetOutput());
1695           paster1->SetSourceImage(m_CoefficientImage);
1696           sourceRegion.SetIndex(NInputDimensions-1,0);
1697           destinationIndex[InputDimension-1]=1;
1698           sourceRegion.SetSize(NInputDimensions-1,2);
1699           paster1->SetDestinationIndex(destinationIndex);
1700           paster1->SetSourceRegion(sourceRegion);
1701
1702           //---------------------------------
1703           // 2. Middle row at index 3=BC
1704           // Extract row at index 1, 2 (of coeff image)
1705           //typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1706           typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
1707
1708           // Combine
1709           typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1710           combine1->SetInput(0,row1);
1711           combine1->SetInput(1,row3);
1712           combine1->SetA(1.);
1713           combine1->SetB(1.);
1714           combine1->Update();
1715
1716           typename LinearCombinationFilterType::Pointer combine2=LinearCombinationFilterType::New();
1717           combine2->SetInput(0,combine1->GetOutput());
1718           combine2->SetInput(1,row2);
1719           combine2->SetA(-m_WeightRatio[3][1]);
1720           combine2->SetB(-1.);
1721           combine2->Update();
1722           typename CoefficientImageType::Pointer bc2Row=combine2->GetOutput();
1723
1724           // Paste middleRow at index 3 (padded coeff)
1725           typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1726           paster2->SetDestinationImage(paster1->GetOutput());
1727           paster2->SetSourceImage(bc2Row);
1728           destinationIndex[InputDimension-1]=3;
1729           paster2->SetDestinationIndex(destinationIndex);
1730           paster2->SetSourceRegion(bc2Row->GetLargestPossibleRegion());
1731
1732           //---------------------------------
1733           // 3. Next two temporal rows are identical: paste 2,3 to 4,5
1734           typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1735           paster3->SetDestinationImage(paster2->GetOutput());
1736           paster3->SetSourceImage(m_CoefficientImage);
1737           sourceRegion.SetIndex(InputDimension-1,2);
1738           destinationIndex[InputDimension-1]=4;
1739           sourceRegion.SetSize(NInputDimensions-1,2);
1740           paster3->SetDestinationIndex(destinationIndex);
1741           paster3->SetSourceRegion(sourceRegion);
1742
1743           //---------------------------------
1744           // 4. Final row at index 6=BC3
1745           typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1746
1747           // Combine
1748           typename LinearCombinationFilterType::Pointer combine3=LinearCombinationFilterType::New();
1749           combine3->SetInput(0,row1);
1750           combine3->SetInput(1,row4);
1751           combine3->SetA(1.);
1752           combine3->SetB(-1.);
1753           typename LinearCombinationFilterType::Pointer combine4=LinearCombinationFilterType::New();
1754           combine4->SetInput(0,row0);
1755           combine4->SetInput(1,combine3->GetOutput());
1756           combine4->SetA(1.);
1757           combine4->SetB(.5);
1758           combine4->Update();
1759           typename CoefficientImageType::Pointer bc3Row=combine4->GetOutput();
1760
1761           // Paste BC row at index 6
1762           typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
1763           paster4->SetDestinationImage(paster3->GetOutput());
1764           paster4->SetSourceImage(bc3Row);
1765           destinationIndex[InputDimension-1]=6;
1766           paster4->SetDestinationIndex(destinationIndex);
1767           paster4->SetSourceRegion(bc3Row->GetLargestPossibleRegion());
1768
1769           // Paste row 4  at index 7
1770           typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
1771           paster5->SetDestinationImage(paster4->GetOutput());
1772           paster5->SetSourceImage(row4);
1773           destinationIndex[InputDimension-1]=7;
1774           paster5->SetDestinationIndex(destinationIndex);
1775           paster5->SetSourceRegion(row4->GetLargestPossibleRegion());
1776
1777           // Update the chain!
1778           paster5->Update();
1779           m_PaddedCoefficientImage= paster5->GetOutput();
1780
1781           break;
1782         }
1783
1784       case 6:
1785         {
1786           // ----------------------------------------------------------------------
1787           // The diamond with 4 internal CP:
1788           // Periodic, constrained to zero at the reference at position 3
1789           // Coeff      row 0   1   2  BC1  3  BC2  4
1790           // PaddedCoeff  R:0   1   2   3   4   5   6
1791           // BC1= -weights[2]/weights[0] ( R2+R4 )
1792           // BC2= weights[2]/weights[0]  ( R0+R2-R4-R6 ) + R1
1793           // ---------------------------------------------------------------------
1794
1795           //---------------------------------
1796           // 1. First Three temporal rows are identical: paste 0-2 to 0-2
1797           typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1798           paster1->SetDestinationImage(m_PaddedCoefficientImage);
1799           paster1->SetSourceImage(m_CoefficientImage);
1800           sourceIndex.Fill(0);
1801           destinationIndex.Fill(0);
1802           sourceSize[NInputDimensions-1]=3;
1803           sourceRegion.SetSize(sourceSize);
1804           sourceRegion.SetIndex(sourceIndex);
1805           paster1->SetDestinationIndex(destinationIndex);
1806           paster1->SetSourceRegion(sourceRegion);
1807
1808           //---------------------------------
1809           // 2. Middle row at index 3=BC
1810           // Extract row at index 0, 4 (of coeff image)
1811           typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
1812           typename CoefficientImageType::Pointer row3=this->ExtractTemporalRow(m_CoefficientImage, 3);
1813
1814           // Combine
1815           typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1816           combine1->SetInput(0,row2);
1817           combine1->SetInput(1,row3);
1818           combine1->SetA(-m_WeightRatio[2][0]);
1819           combine1->SetB(-m_WeightRatio[2][0]);
1820           combine1->Update();
1821           typename CoefficientImageType::Pointer bc1Row=combine1->GetOutput();
1822
1823           // Paste middleRow at index 3 (padded coeff)
1824           typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1825           paster2->SetDestinationImage(paster1->GetOutput());
1826           paster2->SetSourceImage(bc1Row);
1827           destinationIndex[InputDimension-1]=3;
1828           paster2->SetDestinationIndex(destinationIndex);
1829           paster2->SetSourceRegion(bc1Row->GetLargestPossibleRegion());
1830
1831           //---------------------------------
1832           // 3. Next  row identical: paste 3 to 4
1833           typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1834           paster3->SetDestinationImage(paster2->GetOutput());
1835           paster3->SetSourceImage(row3);
1836           destinationIndex[InputDimension-1]=4;
1837           paster3->SetDestinationIndex(destinationIndex);
1838           paster3->SetSourceRegion(row3->GetLargestPossibleRegion());
1839
1840           //---------------------------------
1841           // 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
1842           // Extract row at index 0, 2 (of coeff image)already done
1843           typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1844           typename CoefficientImageType::Pointer row4=this->ExtractTemporalRow(m_CoefficientImage, 4);
1845           typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1846
1847           // Combine
1848           typename LinearCombinationFilterType::Pointer combine3=LinearCombinationFilterType::New();
1849           combine3->SetInput(0,row0);
1850           combine3->SetInput(1,row2);
1851           combine3->SetA(1.);
1852           combine3->SetB(1.);
1853
1854           // Combine
1855           typename LinearCombinationFilterType::Pointer combine4=LinearCombinationFilterType::New();
1856           combine4->SetInput(0,row3);
1857           combine4->SetInput(1,row4);
1858           combine4->SetA(1.);
1859           combine4->SetB(1.);
1860
1861           // Combine
1862           typename LinearCombinationFilterType::Pointer combine5=LinearCombinationFilterType::New();
1863           combine5->SetInput(0,combine3->GetOutput());
1864           combine5->SetInput(1,combine4->GetOutput());
1865           combine5->SetA(1.);
1866           combine5->SetB(-1.);
1867
1868           // Combine
1869           typename LinearCombinationFilterType::Pointer combine6=LinearCombinationFilterType::New();
1870           combine6->SetInput(0,combine5->GetOutput());
1871           combine6->SetInput(1,row1);
1872           combine6->SetA(m_WeightRatio[2][0]);
1873           combine6->SetB(1.);
1874           combine6->Update();
1875           typename CoefficientImageType::Pointer bc2Row=combine6->GetOutput();
1876
1877           // Paste BC row at index 5
1878           typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
1879           paster4->SetDestinationImage(paster3->GetOutput());
1880           paster4->SetSourceImage(bc2Row);
1881           destinationIndex[InputDimension-1]=5;
1882           paster4->SetDestinationIndex(destinationIndex);
1883           paster4->SetSourceRegion(bc2Row->GetLargestPossibleRegion());
1884
1885           // Paste last row at index 6
1886           typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
1887           paster5->SetDestinationImage(paster4->GetOutput());
1888           paster5->SetSourceImage(row4);
1889           destinationIndex[InputDimension-1]=6;
1890           paster5->SetDestinationIndex(destinationIndex);
1891           paster5->SetSourceRegion(row4->GetLargestPossibleRegion());
1892
1893           // Update the chain!
1894           paster5->Update();
1895           m_PaddedCoefficientImage= paster5->GetOutput();
1896
1897           break;
1898         }
1899       case 7:
1900         {
1901           // ----------------------------------------------------------------------
1902           // The diamond with 5 internal CP:
1903           // periodic, constrained to zero at the reference at position 3.5
1904           // Coeff      row 0   1   2  BC1  3   4  BC2  5
1905           // PaddedCoeff  R:0   1   2   3   4   5   6   7
1906           // BC1= -weights[3]/weights[1] ( R2+R5 ) - R4
1907           // BC2= weights[2]/weights[0] ( R0+R2-R5-R7 ) + R1
1908           // ---------------------------------------------------------------------
1909
1910           //---------------------------------
1911           // 1. First Three temporal rows are identical: paste 0-2 to 0-2
1912           typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1913           paster1->SetDestinationImage(m_PaddedCoefficientImage);
1914           paster1->SetSourceImage(m_CoefficientImage);
1915           sourceIndex.Fill(0);
1916           destinationIndex.Fill(0);
1917           sourceSize[NInputDimensions-1]=3;
1918           sourceRegion.SetSize(sourceSize);
1919           sourceRegion.SetIndex(sourceIndex);
1920           paster1->SetDestinationIndex(destinationIndex);
1921           paster1->SetSourceRegion(sourceRegion);
1922
1923           //---------------------------------
1924           // 2. Middle row at index 3=BC
1925           // Extract row at index 0, 4 (of coeff image)
1926           typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
1927           typename CoefficientImageType::Pointer row4=this->ExtractTemporalRow(m_CoefficientImage, 4);
1928
1929           // Combine
1930           typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1931           combine1->SetInput(0,row2);
1932           combine1->SetInput(1,row4);
1933           combine1->SetA(1.);
1934           combine1->SetB(1.);
1935           combine1->Update();
1936
1937           // Extract row at index 3 (coeff image)
1938           typename CoefficientImageType::Pointer row3=this->ExtractTemporalRow(m_CoefficientImage, 3);
1939
1940           // Combine
1941           typename LinearCombinationFilterType::Pointer combine2=LinearCombinationFilterType::New();
1942           combine2->SetInput(0,combine1->GetOutput());
1943           combine2->SetInput(1,row3);
1944           combine2->SetA(-m_WeightRatio[3][1] );
1945           combine2->SetB(-1.);
1946           combine2->Update();
1947           typename CoefficientImageType::Pointer bc1Row=combine2->GetOutput();
1948
1949           // Paste middleRow at index 3 (padded coeff)
1950           typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1951           paster2->SetDestinationImage(paster1->GetOutput());
1952           paster2->SetSourceImage(bc1Row);
1953           destinationIndex[InputDimension-1]=3;
1954           paster2->SetDestinationIndex(destinationIndex);
1955           paster2->SetSourceRegion(bc1Row->GetLargestPossibleRegion());
1956
1957           //---------------------------------
1958           // 3. Next two temporal rows are identical: paste 3,4 to 4,5
1959           typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1960           paster3->SetDestinationImage(paster2->GetOutput());
1961           paster3->SetSourceImage(m_CoefficientImage);
1962           sourceIndex[InputDimension-1]=3;
1963           destinationIndex[InputDimension-1]=4;
1964           sourceSize[NInputDimensions-1]=2;
1965           sourceRegion.SetSize(sourceSize);
1966           sourceRegion.SetIndex(sourceIndex);
1967           paster3->SetDestinationIndex(destinationIndex);
1968           paster3->SetSourceRegion(sourceRegion);
1969
1970           //---------------------------------
1971           // 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
1972           // Extract row at index 0, 2 (of coeff image)already done
1973           typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1974           typename CoefficientImageType::Pointer row5=this->ExtractTemporalRow(m_CoefficientImage, 5);
1975           typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1976
1977           // Combine
1978           typename LinearCombinationFilterType::Pointer combine3=LinearCombinationFilterType::New();
1979           combine3->SetInput(0,row0);
1980           combine3->SetInput(1,row2);
1981           combine3->SetA(1.);
1982           combine3->SetB(1.);
1983
1984           // Combine
1985           typename LinearCombinationFilterType::Pointer combine4=LinearCombinationFilterType::New();
1986           combine4->SetInput(0,row4);
1987           combine4->SetInput(1,row5);
1988           combine4->SetA(1.);
1989           combine4->SetB(1.);
1990
1991           // Combine
1992           typename LinearCombinationFilterType::Pointer combine5=LinearCombinationFilterType::New();
1993           combine5->SetInput(0,combine3->GetOutput());
1994           combine5->SetInput(1,combine4->GetOutput());
1995           combine5->SetA(1.);
1996           combine5->SetB(-1.);
1997
1998           // Combine
1999           typename LinearCombinationFilterType::Pointer combine6=LinearCombinationFilterType::New();
2000           combine6->SetInput(0,combine5->GetOutput());
2001           combine6->SetInput(1,row1);
2002           combine6->SetA(m_WeightRatio[2][0]);
2003           combine6->SetB(1.);
2004           combine6->Update();
2005           typename CoefficientImageType::Pointer bc2Row=combine6->GetOutput();
2006
2007           // Paste BC row at index 6
2008           typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
2009           paster4->SetDestinationImage(paster3->GetOutput());
2010           paster4->SetSourceImage(bc2Row);
2011           destinationIndex[InputDimension-1]=6;
2012           paster4->SetDestinationIndex(destinationIndex);
2013           paster4->SetSourceRegion(bc2Row->GetLargestPossibleRegion());
2014
2015           // Paste last row at index 7
2016           typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
2017           paster5->SetDestinationImage(paster4->GetOutput());
2018           paster5->SetSourceImage(row5);
2019           destinationIndex[InputDimension-1]=7;
2020           paster5->SetDestinationIndex(destinationIndex);
2021           paster5->SetSourceRegion(row5->GetLargestPossibleRegion());
2022
2023           // Update the chain!
2024           paster5->Update();
2025           m_PaddedCoefficientImage= paster5->GetOutput();
2026
2027           break;
2028         }
2029
2030       case 9:
2031         {
2032           // ----------------------------------------------------------------------
2033           // The sputnik with 5 internal CP T''(0)=T''(10)
2034           // Periodic, constrained to zero at the reference
2035           // at position 3 and one indepent extreme
2036           // Coeff     row BC1  0   1  BC2  2   3  BC3  4
2037           // PaddedCoeff R: 0   1   2   3   4   5   6   7
2038           // BC1= -R2+R5+R7
2039           // BC2= -weights[3]/weights[1] ( R2+R5 ) - R4
2040           // BC3= R1
2041           // ---------------------------------------------------------------------
2042
2043           //---------------------------------
2044           // 1. First Row =BC
2045           typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
2046           typename CoefficientImageType::Pointer row3=this->ExtractTemporalRow(m_CoefficientImage, 3);
2047           typename CoefficientImageType::Pointer row4=this->ExtractTemporalRow(m_CoefficientImage, 4);
2048
2049           // Combine
2050           typename LinearCombinationFilterType::Pointer combine0=LinearCombinationFilterType::New();
2051           combine0->SetInput(0,row3);
2052           combine0->SetInput(1,row4);
2053           combine0->SetA(1.);
2054           combine0->SetB(1.);
2055           typename LinearCombinationFilterType::Pointer combine0bis=LinearCombinationFilterType::New();
2056           combine0bis->SetInput(0,combine0->GetOutput());
2057           combine0bis->SetInput(1,row1);
2058           combine0bis->SetA(1.);
2059           combine0bis->SetB(-1.);
2060           combine0bis->Update();
2061           typename CoefficientImageType::Pointer bc1Row=combine0bis->GetOutput();
2062
2063           typename PasteImageFilterType::Pointer paster0=PasteImageFilterType::New();
2064           paster0->SetDestinationImage(m_PaddedCoefficientImage);
2065           paster0->SetSourceImage(bc1Row);
2066           destinationIndex[InputDimension-1]=0;
2067           paster0->SetDestinationIndex(destinationIndex);
2068           paster0->SetSourceRegion(bc1Row->GetLargestPossibleRegion());
2069
2070           //---------------------------------
2071           // 1. Next two temporal rows are identical: paste 0-1 to 1-2
2072           typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
2073           paster1->SetDestinationImage(paster0->GetOutput());
2074           paster1->SetSourceImage(m_CoefficientImage);
2075           sourceRegion.SetIndex(NInputDimensions-1,0);
2076           destinationIndex[InputDimension-1]=1;
2077           sourceRegion.SetSize(NInputDimensions-1,2);
2078           paster1->SetDestinationIndex(destinationIndex);
2079           paster1->SetSourceRegion(sourceRegion);
2080
2081           //---------------------------------
2082           // 2. Middle row at index 3=BC
2083           // Extract row at index 1, 2 (of coeff image)
2084           // typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
2085           typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
2086
2087           // Combine
2088           typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
2089           combine1->SetInput(0,row1);
2090           combine1->SetInput(1,row3);
2091           combine1->SetA(1.);
2092           combine1->SetB(1.);
2093           combine1->Update();
2094
2095           typename LinearCombinationFilterType::Pointer combine2=LinearCombinationFilterType::New();
2096           combine2->SetInput(0,combine1->GetOutput());
2097           combine2->SetInput(1,row2);
2098           combine2->SetA(-m_WeightRatio[3][1]);
2099           combine2->SetB(-1.);
2100           combine2->Update();
2101           typename CoefficientImageType::Pointer bc2Row=combine2->GetOutput();
2102
2103           // Paste middleRow at index 3 (padded coeff)
2104           typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
2105           paster2->SetDestinationImage(paster1->GetOutput());
2106           paster2->SetSourceImage(bc2Row);
2107           destinationIndex[InputDimension-1]=3;
2108           paster2->SetDestinationIndex(destinationIndex);
2109           paster2->SetSourceRegion(bc2Row->GetLargestPossibleRegion());
2110
2111           //---------------------------------
2112           // 3. Next two temporal rows are identical: paste 2,3 to 4,5
2113           typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
2114           paster3->SetDestinationImage(paster2->GetOutput());
2115           paster3->SetSourceImage(m_CoefficientImage);
2116           sourceRegion.SetIndex(InputDimension-1,2);
2117           destinationIndex[InputDimension-1]=4;
2118           sourceRegion.SetSize(NInputDimensions-1,2);
2119           paster3->SetDestinationIndex(destinationIndex);
2120           paster3->SetSourceRegion(sourceRegion);
2121
2122           //---------------------------------
2123           // 4. Final row at index 6=BC3
2124           typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
2125
2126           //      // Combine
2127           //      typename LinearCombinationFilterType::Pointer combine3=LinearCombinationFilterType::New();
2128           //      combine3->SetInput(0,row0);
2129           //      combine3->SetInput(1,row1);
2130           //      combine3->SetA(1.);
2131           //      combine3->SetB(0.5);
2132           //      combine3->Update();
2133           //      typename CoefficientImageType::Pointer bc3Row=combine3->GetOutput();
2134
2135           // Paste BC row at index 6
2136           typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
2137           paster4->SetDestinationImage(paster3->GetOutput());
2138           paster4->SetSourceImage(row0);
2139           destinationIndex[InputDimension-1]=6;
2140           paster4->SetDestinationIndex(destinationIndex);
2141           paster4->SetSourceRegion(row0->GetLargestPossibleRegion());
2142
2143           // Paste row 4  at index 7
2144           typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
2145           paster5->SetDestinationImage(paster4->GetOutput());
2146           paster5->SetSourceImage(row4);
2147           destinationIndex[InputDimension-1]=7;
2148           paster5->SetDestinationIndex(destinationIndex);
2149           paster5->SetSourceRegion(row4->GetLargestPossibleRegion());
2150
2151           // Update the chain!
2152           paster5->Update();
2153           m_PaddedCoefficientImage= paster5->GetOutput();
2154
2155           break;
2156         }
2157
2158       default:
2159         DD ("Shape not available");
2160       }
2161
2162   }
2163
2164
2165   // // Extract coefficients from padded version
2166   // template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
2167   // void
2168   // ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2169   // ::ExtractCoefficientImage( )
2170   // {
2171   //   ////DD("extract coeff image");
2172   //   typename ExtractImageFilterType::Pointer extractImageFilter = ExtractImageFilterType::New();
2173   //   typename CoefficientImageType::RegionType extractionRegion=m_PaddedCoefficientImage->GetLargestPossibleRegion();
2174   //   typename CoefficientImageType::RegionType::SizeType extractionSize=extractionRegion.GetSize();
2175   //   typename CoefficientImageType::RegionType::IndexType extractionIndex = extractionRegion.GetIndex();
2176   //   extractionSize[InputDimension-1]-=4;
2177   //   extractionIndex[InputDimension-1]=2;
2178   //   extractionRegion.SetSize(extractionSize);
2179   //   extractionRegion.SetIndex(extractionIndex);
2180   //   extractImageFilter->SetInput(m_PaddedCoefficientImage);
2181   //   extractImageFilter->SetExtractionRegion(extractionRegion);
2182   //   extractImageFilter->Update();
2183   //   m_CoefficientImage=extractImageFilter->GetOutput();
2184   // }
2185
2186
2187   // Print self
2188   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
2189   void
2190   ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2191   ::PrintSelf(std::ostream &os, itk::Indent indent) const
2192   {
2193
2194     this->Superclass::PrintSelf(os, indent);
2195
2196     os << indent << "GridRegion: " << m_GridRegion << std::endl;
2197     os << indent << "GridOrigin: " << m_GridOrigin << std::endl;
2198     os << indent << "GridSpacing: " << m_GridSpacing << std::endl;
2199     os << indent << "GridDirection: " << m_GridDirection << std::endl;
2200     os << indent << "IndexToPoint: " << m_IndexToPoint << std::endl;
2201     os << indent << "PointToIndex: " << m_PointToIndex << std::endl;
2202
2203     os << indent << "CoefficientImage: [ ";
2204     os << m_CoefficientImage.GetPointer() << " ]" << std::endl;
2205
2206     os << indent << "WrappedImage: [ ";
2207     os << m_WrappedImage.GetPointer() << " ]" << std::endl;
2208
2209     os << indent << "InputParametersPointer: "
2210        << m_InputParametersPointer << std::endl;
2211     os << indent << "ValidRegion: " << m_ValidRegion << std::endl;
2212     os << indent << "LastJacobianIndex: " << m_LastJacobianIndex << std::endl;
2213     os << indent << "BulkTransform: ";
2214     os << m_BulkTransform << std::endl;
2215
2216     if ( m_BulkTransform )
2217       {
2218         os << indent << "BulkTransformType: "
2219            << m_BulkTransform->GetNameOfClass() << std::endl;
2220       }
2221     os << indent << "VectorBSplineInterpolator: ";
2222     os << m_VectorInterpolator.GetPointer() << std::endl;
2223     os << indent << "Mask: ";
2224     os << m_Mask<< std::endl;
2225   }
2226
2227
2228   // Verify
2229   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
2230   bool
2231   ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2232   ::InsideValidRegion( const ContinuousIndexType& index ) const
2233   {
2234     bool inside = true;
2235
2236     if ( !m_ValidRegion.IsInside( index ) )
2237       {
2238         inside = false;
2239       }
2240
2241     // JV verify all dimensions
2242     if ( inside)
2243       {
2244         typedef typename ContinuousIndexType::ValueType ValueType;
2245         for( unsigned int j = 0; j < NInputDimensions; j++ )
2246           {
2247             if (m_SplineOrderOdd[j])
2248               {
2249                 if ( index[j] >= static_cast<ValueType>( m_ValidRegionLast[j] ) )
2250                   {
2251                     inside = false;
2252                     break;
2253                   }
2254               }
2255           }
2256       }
2257     return inside;
2258   }
2259
2260
2261   // Transform a point
2262   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
2263   typename ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2264   ::OutputPointType
2265   ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2266   ::TransformPoint(const InputPointType &inputPoint) const
2267   {
2268
2269     InputPointType transformedPoint;
2270     OutputPointType outputPoint;
2271
2272     // BulkTransform
2273     if ( m_BulkTransform )
2274       {
2275         transformedPoint = m_BulkTransform->TransformPoint( inputPoint );
2276       }
2277     else
2278       {
2279         transformedPoint = inputPoint;
2280       }
2281
2282     // Deformable transform
2283     if ( m_PaddedCoefficientImage )
2284       {
2285         // Check if inside mask
2286         if(m_Mask &&  !(m_Mask->IsInside(inputPoint) ) )
2287           {
2288             // Outside: no (deformable) displacement
2289             return transformedPoint;
2290           }
2291
2292         // Check if inside valid region
2293         bool inside = true;
2294         ContinuousIndexType index;
2295         this->TransformPointToContinuousIndex( inputPoint, index );
2296         inside = this->InsideValidRegion( index );
2297         if ( !inside )
2298           {
2299             // Outside: no (deformable) displacement
2300             DD("outside valid region");
2301             return transformedPoint;
2302           }
2303
2304         // Call the vector interpolator
2305         itk::Vector<TCoordRep,SpaceDimension> displacement=m_VectorInterpolator->EvaluateAtContinuousIndex(index);
2306
2307         // JV add for the spatial dimensions
2308         outputPoint=transformedPoint;
2309         for (unsigned int i=0; i<NInputDimensions-1; i++)
2310           outputPoint[i] += displacement[i];
2311
2312       }
2313
2314     else
2315       {
2316         itkWarningMacro( << "B-spline coefficients have not been set" );
2317         outputPoint = transformedPoint;
2318       }
2319
2320     return outputPoint;
2321   }
2322
2323
2324
2325   //JV Deformably transform a point
2326   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
2327   typename ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2328   ::OutputPointType
2329   ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2330   ::DeformablyTransformPoint(const InputPointType &inputPoint) const
2331   {
2332     OutputPointType outputPoint;
2333     if ( m_PaddedCoefficientImage )
2334       {
2335
2336         // Check if inside mask
2337         if(m_Mask && !(m_Mask->IsInside(inputPoint) ) )
2338           {
2339             // Outside: no (deformable) displacement
2340             return inputPoint;
2341           }
2342
2343         // Check if inside
2344         bool inside = true;
2345         ContinuousIndexType index;
2346         this->TransformPointToContinuousIndex( inputPoint, index );
2347         inside = this->InsideValidRegion( index );
2348
2349         if ( !inside )
2350           {
2351             //outside: no (deformable) displacement
2352             outputPoint = inputPoint;
2353             return outputPoint;
2354           }
2355
2356         // Call the vector interpolator
2357         itk::Vector<TCoordRep,SpaceDimension> displacement=m_VectorInterpolator->EvaluateAtContinuousIndex(index);
2358
2359         // JV add for the spatial dimensions
2360         outputPoint=inputPoint;
2361         for (unsigned int i=0; i<NInputDimensions-1; i++)
2362           outputPoint[i] += displacement[i];
2363       }
2364
2365     // No coefficients available
2366     else
2367       {
2368         itkWarningMacro( << "B-spline coefficients have not been set" );
2369         outputPoint = inputPoint;
2370       }
2371
2372     return outputPoint;
2373   }
2374
2375
2376   // JV weights are identical as for transformpoint, could be done simultaneously in metric!!!!
2377   // Compute the Jacobian in one position
2378   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
2379   const
2380   typename ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2381   ::JacobianType &
2382   ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2383   ::GetJacobian( const InputPointType & point ) const
2384   {
2385
2386     //========================================================
2387     // Zero all components of jacobian
2388     //========================================================
2389     // JV not thread safe (m_LastJacobianIndex), instantiate N transforms
2390     // NOTE: for efficiency, we only need to zero out the coefficients
2391     // that got fill last time this method was called.
2392     unsigned int j=0,b=0;
2393
2394     //Set the previously-set to zero
2395     for ( j = 0; j < SpaceDimension; j++ )
2396       {
2397         m_FirstIterator[j].GoToBegin();
2398         while ( !m_FirstIterator[j].IsAtEnd() )
2399           {
2400             m_FirstIterator[j].Set( m_ZeroVector );
2401             ++(m_FirstIterator[j]);
2402           }
2403       }
2404
2405     //Set the previously-set to zero
2406     for ( j = 0; j < SpaceDimension; j++ )
2407       {
2408         m_SecondIterator[j].GoToBegin();
2409         while ( ! (m_SecondIterator[j]).IsAtEnd() )
2410           {
2411             m_SecondIterator[j].Set( m_ZeroVector );
2412             ++(m_SecondIterator[j]);
2413           }
2414       }
2415
2416     //Set the previously-set to zero
2417     if (m_ThirdSize)
2418       for ( j = 0; j < SpaceDimension; j++ )
2419         {
2420           m_ThirdIterator[j].GoToBegin();
2421           while ( ! (m_ThirdIterator[j]).IsAtEnd() )
2422             {
2423               m_ThirdIterator[j].Set( m_ZeroVector );
2424               ++(m_ThirdIterator[j]);
2425             }
2426         }
2427
2428     //Set the previously-set to zero
2429     if (m_BCSize)
2430       for (b=0; b<m_BCSize;b++)
2431         if ( !( m_FirstRegion.IsInside(m_BCRegions[b]) | m_SecondRegion.IsInside(m_BCRegions[b]) )  )
2432           for ( j = 0; j < SpaceDimension; j++ )
2433             {
2434               m_BCIterators[j][b].GoToBegin();
2435               while ( ! (m_BCIterators[j][b]).IsAtEnd() )
2436                 {
2437                   m_BCIterators[j][b].Set( m_ZeroVector );
2438                   ++(m_BCIterators[j][b]);
2439                 }
2440             }
2441
2442     //Set the previously-set to zero
2443     if (m_BC2Size)
2444       for (b=0; b<m_BC2Size;b++)
2445         if ( !( m_FirstRegion.IsInside(m_BC2Regions[b]) | m_SecondRegion.IsInside(m_BC2Regions[b]) )  )
2446           for ( j = 0; j < SpaceDimension; j++ )
2447             {
2448               m_BC2Iterators[j][b].GoToBegin();
2449               while ( ! (m_BC2Iterators[j][b]).IsAtEnd() )
2450                 {
2451                   m_BC2Iterators[j][b].Set( m_ZeroVector );
2452                   ++(m_BC2Iterators[j][b]);
2453                 }
2454             }
2455
2456     //Set the previously-set to zero
2457     if (m_BC3Size)
2458       for (b=0; b<m_BC3Size;b++)
2459         if ( !( m_FirstRegion.IsInside(m_BC3Regions[b]) | m_SecondRegion.IsInside(m_BC3Regions[b]) )  )
2460           for ( j = 0; j < SpaceDimension; j++ )
2461             {
2462               m_BC3Iterators[j][b].GoToBegin();
2463               while ( ! (m_BC3Iterators[j][b]).IsAtEnd() )
2464                 {
2465                   m_BC3Iterators[j][b].Set( m_ZeroVector );
2466                   ++(m_BC3Iterators[j][b]);
2467                 }
2468             }
2469
2470
2471     //========================================================
2472     // For each dimension, copy the weight to the support region
2473     //========================================================
2474
2475     // Check if inside mask
2476     if(m_Mask &&  !(m_Mask->IsInside(point) ) )
2477       {
2478         // Outside: no (deformable) displacement
2479         return this->m_Jacobian;
2480       }
2481
2482     // Get index
2483     this->TransformPointToContinuousIndex( point, m_Index );
2484
2485     // NOTE: if the support region does not lie totally within the grid
2486     // we assume zero displacement and return the input point
2487     if ( !this->InsideValidRegion( m_Index ) )
2488       {
2489         return this->m_Jacobian;
2490       }
2491
2492     // Compute interpolation weights
2493     const WeightsDataType *weights=NULL;
2494     m_VectorInterpolator->EvaluateWeightsAtContinuousIndex( m_Index, &weights, m_LastJacobianIndex);
2495
2496     // Get support
2497     m_SupportRegion.SetIndex( m_LastJacobianIndex );
2498     WrapRegion(m_SupportRegion, m_FirstRegion, m_SecondRegion, m_ThirdRegion, m_BCRegions, m_BCValues, m_BC2Regions, m_BC2Values, m_BC3Regions, m_BC3Values, m_InitialOffset);
2499     m_ThirdSize=m_ThirdRegion.GetSize()[InputDimension-1];
2500     m_BCSize=m_BCRegions.size();
2501     m_BC2Size=m_BC2Regions.size();
2502     m_BC3Size=m_BC3Regions.size();
2503
2504     // Reset the iterators
2505     for ( j = 0; j < SpaceDimension ; j++ )
2506       {
2507         m_FirstIterator[j] = IteratorType( m_JacobianImage[j], m_FirstRegion);
2508         m_SecondIterator[j] = IteratorType( m_JacobianImage[j], m_SecondRegion);
2509         if(m_ThirdSize) m_ThirdIterator[j] = IteratorType( m_JacobianImage[j], m_ThirdRegion);
2510
2511         m_BCIterators[j].resize(m_BCSize);
2512         for (b=0; b<m_BCSize;b++)
2513           m_BCIterators[j][b]= IteratorType( m_JacobianImage[j], m_BCRegions[b]);
2514         m_BC2Iterators[j].resize(m_BC2Size);
2515         for (b=0; b<m_BC2Size;b++)
2516           m_BC2Iterators[j][b]= IteratorType( m_JacobianImage[j], m_BC2Regions[b]);
2517         m_BC3Iterators[j].resize(m_BC3Size);
2518         for (b=0; b<m_BC3Size;b++)
2519           m_BC3Iterators[j][b]= IteratorType( m_JacobianImage[j], m_BC3Regions[b]);
2520       }
2521
2522     // Skip if on a fixed condition
2523     if(m_InitialOffset)
2524       {
2525         if (m_BCSize)  weights+=m_InitialOffset*m_BCRegions[0].GetNumberOfPixels();
2526         else std::cerr<<"InitialOffset without BCSize: Error!!!!!!!"<<std::endl;
2527       }
2528
2529     //copy weight to jacobian image
2530     for ( j = 0; j < SpaceDimension; j++ )
2531       {
2532         // For each dimension, copy the weight to the support region
2533         while ( ! (m_FirstIterator[j]).IsAtEnd() )
2534           {
2535             m_ZeroVector[j]=*weights;
2536             (m_FirstIterator[j]).Set( m_ZeroVector);
2537             ++(m_FirstIterator[j]);
2538             weights++;
2539           }
2540
2541         m_ZeroVector[j]=itk::NumericTraits<JacobianValueType>::Zero;
2542         if (j != SpaceDimension-1)  weights-=m_FirstRegion.GetNumberOfPixels();
2543       }
2544
2545     // Skip BC1 and go to the second region
2546     if (m_BCSize) weights+=m_BCRegions[0].GetNumberOfPixels();
2547
2548     // For each dimension, copy the weight to the support region
2549     //copy weight to jacobian image
2550     for ( j = 0; j < SpaceDimension; j++ )
2551       {
2552         while ( ! (m_SecondIterator[j]).IsAtEnd() )
2553           {
2554             m_ZeroVector[j]=*weights;
2555             (m_SecondIterator[j]).Set( m_ZeroVector);
2556             ++(m_SecondIterator[j]);
2557             weights++;
2558           }
2559         weights-=m_SecondRegion.GetNumberOfPixels();
2560         m_ZeroVector[j]=itk::NumericTraits<JacobianValueType>::Zero;
2561       }
2562
2563     // Now do BC1:
2564     if (m_BCSize)
2565       {
2566         //  Put pointer in correct position
2567         weights-=m_BCRegions[0].GetNumberOfPixels();
2568
2569         for ( j = 0; j < SpaceDimension; j++ )
2570           {
2571             for ( b=0; b < m_BCSize; b++ )
2572               {
2573                 while ( ! (m_BCIterators[j][b]).IsAtEnd() )
2574                   {
2575                     //copy weight to jacobian image
2576                     m_ZeroVector[j]=(*weights) * m_BCValues[b];
2577                     (m_BCIterators[j][b]).Value()+= m_ZeroVector;
2578                     ++(m_BCIterators[j][b]);
2579                     weights++;
2580                   }
2581                 weights-=m_BCRegions[b].GetNumberOfPixels();
2582               }
2583             m_ZeroVector[j]=itk::NumericTraits<JacobianValueType>::Zero;
2584           }
2585         weights+=m_BCRegions[m_BCSize-1].GetNumberOfPixels();
2586       }
2587
2588     // Add the BC2 to the weights
2589     if (m_BC2Size)
2590       {
2591         // Move further in the weights pointer
2592         weights+=m_SecondRegion.GetNumberOfPixels();
2593
2594         for ( j = 0; j < SpaceDimension; j++ )
2595           {
2596             for ( b=0; b < m_BC2Size; b++ )
2597               {
2598                 while ( ! (m_BC2Iterators[j][b]).IsAtEnd() )
2599                   {
2600                     //copy weight to jacobian image
2601                     m_ZeroVector[j]=(*weights) * m_BC2Values[b];
2602                     (m_BC2Iterators[j][b]).Value()+= m_ZeroVector;
2603                     ++(m_BC2Iterators[j][b]);
2604                     weights++;
2605                   }
2606                 weights-=m_BC2Regions[b].GetNumberOfPixels();
2607               }
2608             m_ZeroVector[j]=itk::NumericTraits<JacobianValueType>::Zero;
2609           }
2610         // Move further in the weights pointer
2611         weights+=m_BC2Regions[m_BC2Size-1].GetNumberOfPixels();
2612       }
2613
2614     // Third Region
2615     if(m_ThirdSize)
2616       {
2617         // For each dimension, copy the weight to the support region
2618         //copy weight to jacobian image
2619         for ( j = 0; j < SpaceDimension; j++ )
2620           {
2621             while ( ! (m_ThirdIterator[j]).IsAtEnd() )
2622               {
2623                 m_ZeroVector[j]=*weights;
2624                 (m_ThirdIterator[j]).Value()+= m_ZeroVector;
2625                 ++(m_ThirdIterator[j]);
2626                 weights++;
2627               }
2628
2629             // Move further in the weights pointer?
2630             if (j != SpaceDimension-1) weights-=m_ThirdRegion.GetNumberOfPixels();
2631             m_ZeroVector[j]=itk::NumericTraits<JacobianValueType>::Zero;
2632           }
2633       }
2634
2635     // Add the BC3 to the weights
2636     if (m_BC3Size)
2637       {
2638         for ( j = 0; j < SpaceDimension; j++ )
2639           {
2640             for ( b=0; b < m_BC3Size; b++ )
2641               {
2642                 while ( ! (m_BC3Iterators[j][b]).IsAtEnd() )
2643                   {
2644                     //copy weight to jacobian image
2645                     m_ZeroVector[j]=(*weights) * m_BC3Values[b];
2646                     (m_BC3Iterators[j][b]).Value()+= m_ZeroVector;
2647                     ++(m_BC3Iterators[j][b]);
2648                     weights++;
2649                   }
2650                 weights-=m_BC3Regions[b].GetNumberOfPixels();
2651               }
2652             m_ZeroVector[j]=itk::NumericTraits<JacobianValueType>::Zero;
2653           }
2654       }
2655
2656     // Return the result
2657     return this->m_Jacobian;
2658   }
2659
2660
2661   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
2662   inline void
2663   ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2664   ::WrapRegion( const RegionType & m_SupportRegion,
2665                 RegionType & m_FirstRegion,
2666                 RegionType & m_SecondRegion,
2667                 RegionType & m_ThirdRegion,
2668                 std::vector<RegionType>& m_BCRegions,std::vector<double>& m_BCValues,
2669                 std::vector<RegionType>& m_BC2Regions,std::vector<double>& m_BC2Values,
2670                 std::vector<RegionType>& m_BC3Regions,std::vector<double>& m_BC3Values,
2671                 unsigned int& m_InitialOffset ) const
2672   {
2673
2674     // Synchronize regions
2675     m_InitialOffset=0;
2676     m_FirstRegion=m_SupportRegion;
2677     m_BCRegion=m_SupportRegion;
2678     m_BCRegion.SetSize(InputDimension-1,1);
2679     m_SecondRegion=m_SupportRegion;
2680     m_ThirdRegion=m_SupportRegion;
2681     m_ThirdRegion.SetSize(InputDimension-1,0);
2682     m_BC3Regions.resize(0);
2683
2684
2685     // BC depends on shape
2686     switch(m_TransformShape)
2687       {
2688         /*  The shapes are
2689             0: egg     4 CP  3 DOF
2690             1: egg     5 CP  4 DOF
2691             2: rabbit  4 CP  3 DOF
2692             3: rabbit  5 CP  4 DOF
2693             4: sputnik 4 CP  4 DOF
2694             5: sputnik 5 CP  5 DOF
2695             6: diamond 6 CP  5 DOF
2696             7: diamond 7 CP  6 DOF
2697         */
2698
2699       case 0:
2700         {
2701           // ----------------------------------------------------------------------
2702           // The egg with 4 internal CP (starting from inhale)
2703           // Periodic, constrained to zero at the reference
2704           // at position 3 and
2705           // Coeff      row  BC1 BC2  0  BC3  1   2  BC4
2706           // PaddedCoeff  R:  0   1   2   3   4   5   6
2707           // BC1= R4
2708           // BC2= R5
2709           // BC3= -weights[2]/weights[0] ( R2+R4 )
2710           // BC4= R2
2711           // ---------------------------------------------------------------------
2712           switch(m_SupportRegion.GetIndex(InputDimension-1))
2713             {
2714             case 0:
2715               {
2716                 // Lower end
2717                 m_FirstRegion.SetSize(InputDimension-1,2);
2718                 m_FirstRegion.SetIndex(InputDimension-1,1);
2719
2720                 // BC
2721                 m_BCRegions.resize(0);
2722
2723                 // Second part
2724                 m_SecondRegion.SetSize(InputDimension-1,1);
2725                 m_SecondRegion.SetIndex(InputDimension-1,0);
2726
2727                 // BC
2728                 m_BC2Regions.resize(2);
2729                 m_BC2Values.resize(2);
2730                 m_BCRegion.SetIndex(InputDimension-1,0);
2731                 m_BC2Regions[0]=m_BCRegion;
2732                 m_BC2Values[0]=-m_WeightRatio[2][0];
2733                 m_BCRegion.SetIndex(InputDimension-1,1);
2734                 m_BC2Regions[1]=m_BCRegion;
2735                 m_BC2Values[1]= -m_WeightRatio[2][0];
2736
2737                 break;
2738               }
2739             case 1:
2740               {
2741                 // Lower end
2742                 m_FirstRegion.SetSize(InputDimension-1,1);
2743                 m_FirstRegion.SetIndex(InputDimension-1,2);
2744
2745                 // BC
2746                 m_BCRegions.resize(0);
2747
2748                 // Second part
2749                 m_SecondRegion.SetSize(InputDimension-1,1);
2750                 m_SecondRegion.SetIndex(InputDimension-1,0);
2751
2752                 // BC
2753                 m_BC2Regions.resize(2);
2754                 m_BC2Values.resize(2);
2755                 m_BCRegion.SetIndex(InputDimension-1,0);
2756                 m_BC2Regions[0]=m_BCRegion;
2757                 m_BC2Values[0]=-m_WeightRatio[2][0];
2758                 m_BCRegion.SetIndex(InputDimension-1,1);
2759                 m_BC2Regions[1]=m_BCRegion;
2760                 m_BC2Values[1]= -m_WeightRatio[2][0];
2761
2762                 // Third Part
2763                 m_FirstRegion.SetSize(InputDimension-1,1);
2764                 m_FirstRegion.SetIndex(InputDimension-1,1);
2765
2766                 break;
2767               }
2768             case 2:
2769               {
2770                 // Lower end
2771                 m_FirstRegion.SetSize(InputDimension-1,1);
2772                 m_FirstRegion.SetIndex(InputDimension-1,0);
2773
2774                 // BC
2775                 m_BCRegions.resize(2);
2776                 m_BCValues.resize(2);
2777                 m_BCRegion.SetIndex(InputDimension-1,0);
2778                 m_BCRegions[0]=m_BCRegion;
2779                 m_BCValues[0]=-m_WeightRatio[2][0];
2780                 m_BCRegion.SetIndex(InputDimension-1,1);
2781                 m_BCRegions[1]=m_BCRegion;
2782                 m_BCValues[1]= -m_WeightRatio[2][0];
2783
2784                 // Second part
2785                 m_SecondRegion.SetSize(InputDimension-1,2);
2786                 m_SecondRegion.SetIndex(InputDimension-1,1);
2787
2788                 // BC2
2789                 m_BC2Regions.resize(0);
2790
2791                 break;
2792               }
2793             case 3:
2794               {
2795                 // Lower end
2796                 m_FirstRegion.SetSize(InputDimension-1,0);
2797
2798                 // BC
2799                 m_BCRegions.resize(2);
2800                 m_BCValues.resize(2);
2801                 m_BCRegion.SetIndex(InputDimension-1,0);
2802                 m_BCRegions[0]=m_BCRegion;
2803                 m_BCValues[0]=-m_WeightRatio[2][0];
2804                 m_BCRegion.SetIndex(InputDimension-1,1);
2805                 m_BCRegions[1]=m_BCRegion;
2806                 m_BCValues[1]= -m_WeightRatio[2][0];
2807
2808                 // Second part
2809                 m_SecondRegion.SetSize(InputDimension-1,2);
2810                 m_SecondRegion.SetIndex(InputDimension-1,1);
2811
2812                 // BC2
2813                 m_BC2Regions.resize(1);
2814                 m_BC2Values.resize(1);
2815                 m_BCRegion.SetIndex(InputDimension-1,0);
2816                 m_BC2Regions[0]=m_BCRegion;
2817                 m_BC2Values[0]=1;
2818
2819                 break;
2820               }
2821
2822             default:
2823               {
2824                 DD("supportindex > 3 ???");
2825                 DD(m_SupportRegion.GetIndex(InputDimension-1));
2826               }
2827             } // end switch index
2828
2829           break;
2830         }
2831
2832       case 1:
2833         {
2834           // ----------------------------------------------------------------------
2835           // The egg with 5 internal CP (starting from inhale)
2836           // Periodic, constrained to zero at the reference
2837           // at position 3 and
2838           // Coeff      row  BC1 BC2  0  BC3  1   2   3  BC4
2839           // PaddedCoeff  R:  0   1   2   3   4   5   6   7
2840           // BC1= R5
2841           // BC2= R6
2842           // BC3= -weights[3]/weights[1] ( R2+R5 ) - R4
2843           // BC4= R2
2844           // ---------------------------------------------------------------------
2845           switch(m_SupportRegion.GetIndex(InputDimension-1))
2846             {
2847             case 0:
2848               {
2849                 // Lower end
2850                 m_FirstRegion.SetSize(InputDimension-1,2);
2851                 m_FirstRegion.SetIndex(InputDimension-1,2);
2852
2853                 // BC
2854                 m_BCRegions.resize(0);
2855
2856                 // Second part
2857                 m_SecondRegion.SetSize(InputDimension-1,1);
2858                 m_SecondRegion.SetIndex(InputDimension-1,0);
2859
2860                 // BC
2861                 m_BC2Regions.resize(3);
2862                 m_BC2Values.resize(3);
2863                 m_BCRegion.SetIndex(InputDimension-1,0);
2864                 m_BC2Regions[0]=m_BCRegion;
2865                 m_BC2Values[0]=-m_WeightRatio[3][1];
2866                 m_BCRegion.SetIndex(InputDimension-1,1);
2867                 m_BC2Regions[1]=m_BCRegion;
2868                 m_BC2Values[1]= -1;
2869                 m_BCRegion.SetIndex(InputDimension-1,2);
2870                 m_BC2Regions[2]=m_BCRegion;
2871                 m_BC2Values[2]=-m_WeightRatio[3][1];
2872
2873                 break;
2874               }
2875             case 1:
2876               {
2877                 // Lower end
2878                 m_FirstRegion.SetSize(InputDimension-1,1);
2879                 m_FirstRegion.SetIndex(InputDimension-1,3);
2880
2881                 // BC
2882                 m_BCRegions.resize(0);
2883
2884                 // Second part
2885                 m_SecondRegion.SetSize(InputDimension-1,1);
2886                 m_SecondRegion.SetIndex(InputDimension-1,0);
2887
2888                 // BC
2889                 m_BC2Regions.resize(3);
2890                 m_BC2Values.resize(3);
2891                 m_BCRegion.SetIndex(InputDimension-1,0);
2892                 m_BC2Regions[0]=m_BCRegion;
2893                 m_BC2Values[0]=-m_WeightRatio[3][1];
2894                 m_BCRegion.SetIndex(InputDimension-1,1);
2895                 m_BC2Regions[1]=m_BCRegion;
2896                 m_BC2Values[1]= -1;
2897                 m_BCRegion.SetIndex(InputDimension-1,2);
2898                 m_BC2Regions[2]=m_BCRegion;
2899                 m_BC2Values[2]=-m_WeightRatio[3][1];
2900
2901                 // Third Part
2902                 m_FirstRegion.SetSize(InputDimension-1,1);
2903                 m_FirstRegion.SetIndex(InputDimension-1,1);
2904
2905                 break;
2906               }
2907             case 2:
2908               {
2909                 // Lower end
2910                 m_FirstRegion.SetSize(InputDimension-1,1);
2911                 m_FirstRegion.SetIndex(InputDimension-1,0);
2912
2913                 // BC
2914                 m_BCRegions.resize(3);
2915                 m_BCValues.resize(3);
2916                 m_BCRegion.SetIndex(InputDimension-1,0);
2917                 m_BCRegions[0]=m_BCRegion;
2918                 m_BCValues[0]=-m_WeightRatio[3][1];
2919                 m_BCRegion.SetIndex(InputDimension-1,1);
2920                 m_BCRegions[1]=m_BCRegion;
2921                 m_BCValues[1]= -1;
2922                 m_BCRegion.SetIndex(InputDimension-1,2);
2923                 m_BCRegions[2]=m_BCRegion;
2924                 m_BCValues[2]=-m_WeightRatio[3][1];
2925
2926                 // Second part
2927                 m_SecondRegion.SetSize(InputDimension-1,2);
2928                 m_SecondRegion.SetIndex(InputDimension-1,1);
2929
2930                 // BC2
2931                 m_BC2Regions.resize(0);
2932
2933                 break;
2934               }
2935             case 3:
2936               {
2937                 // Lower end
2938                 m_FirstRegion.SetSize(InputDimension-1,0);
2939
2940                 // BC
2941                 m_BCRegions.resize(3);
2942                 m_BCValues.resize(3);
2943                 m_BCRegion.SetIndex(InputDimension-1,0);
2944                 m_BCRegions[0]=m_BCRegion;
2945                 m_BCValues[0]=-m_WeightRatio[3][1];
2946                 m_BCRegion.SetIndex(InputDimension-1,1);
2947                 m_BCRegions[1]=m_BCRegion;
2948                 m_BCValues[1]= -1;
2949                 m_BCRegion.SetIndex(InputDimension-1,2);
2950                 m_BCRegions[2]=m_BCRegion;
2951                 m_BCValues[2]=-m_WeightRatio[3][1];
2952
2953                 // Second part
2954                 m_SecondRegion.SetSize(InputDimension-1,3);
2955                 m_SecondRegion.SetIndex(InputDimension-1,1);
2956
2957                 // BC2
2958                 m_BC2Regions.resize(0);
2959
2960                 break;
2961               }
2962             case 4:
2963               {
2964                 // Lower end
2965                 m_FirstRegion.SetSize(InputDimension-1,3);
2966                 m_FirstRegion.SetIndex(InputDimension-1,1);
2967
2968                 // BC
2969                 m_BCRegions.resize(0);
2970
2971                 // Second part
2972                 m_SecondRegion.SetSize(InputDimension-1,1);
2973                 m_SecondRegion.SetIndex(InputDimension-1,0);
2974
2975                 // BC2
2976                 m_BC2Regions.resize(0);
2977
2978                 break;
2979               }
2980
2981             default:
2982               {
2983                 DD("supportindex > 3 ???");
2984                 DD(m_SupportRegion.GetIndex(InputDimension-1));
2985               }
2986             } // end switch index
2987
2988           break;
2989         }
2990 //        // ----------------------------------------------------------------------
2991 //        // The egg with 5 internal CP:
2992 //        // Periodic, C2 smooth everywhere and constrained to zero at the reference
2993 //        // Coeff       row R5 BC1  0   1   2   3  BC2  R2
2994 //        // PaddedCoeff R:  0   1   2   3   4   5   6   7
2995 //        // BC1= -weights[2]/weights[0] ( R2+R5)
2996 //        // BC2= BC1
2997 //        // ---------------------------------------------------------------------
2998 //        switch(m_SupportRegion.GetIndex(InputDimension-1))
2999 //          {
3000 //          case 0:
3001 //            {
3002 //              // Lower end
3003 //              m_FirstRegion.SetSize(InputDimension-1,1);
3004 //              m_FirstRegion.SetIndex(InputDimension-1,3);
3005
3006 //              // BC
3007 //              m_BCRegions.resize(2);
3008 //              m_BCValues.resize(2);
3009 //              m_BCRegion.SetIndex(InputDimension-1,0);
3010 //              m_BCRegions[0]=m_BCRegion;
3011 //              m_BCValues[0]=-m_WeightRatio[2][0];
3012 //              m_BCRegion.SetIndex(InputDimension-1,3);
3013 //              m_BCRegions[1]=m_BCRegion;
3014 //              m_BCValues[1]=-m_WeightRatio[2][0];
3015
3016 //              // Second part
3017 //              m_SecondRegion.SetSize(InputDimension-1,2);
3018 //              m_SecondRegion.SetIndex(InputDimension-1,0);
3019
3020 //              // BC2
3021 //              m_BC2Regions.resize(0);
3022
3023 //              break;
3024 //            }
3025 //          case 1:
3026 //            {
3027 //              // Lower end
3028 //              m_FirstRegion.SetSize(InputDimension-1,0);
3029
3030 //              // BC
3031 //              m_BCRegions.resize(2);
3032 //              m_BCValues.resize(2);
3033 //              m_BCRegion.SetIndex(InputDimension-1,0);
3034 //              m_BCRegions[0]=m_BCRegion;
3035 //              m_BCValues[0]=-m_WeightRatio[2][0];
3036 //              m_BCRegion.SetIndex(InputDimension-1,3);
3037 //              m_BCRegions[1]=m_BCRegion;
3038 //              m_BCValues[1]=-m_WeightRatio[2][0];
3039
3040 //              // Second part
3041 //              m_SecondRegion.SetSize(InputDimension-1,3);
3042 //              m_SecondRegion.SetIndex(InputDimension-1,0);
3043
3044 //              // BC2
3045 //              m_BC2Regions.resize(0);
3046
3047 //              break;
3048 //            }
3049 //          case 2:
3050 //            {
3051 //              // Lower end
3052 //              m_FirstRegion.SetSize(InputDimension-1,4);
3053 //              m_FirstRegion.SetIndex(InputDimension-1,0);
3054
3055 //              // BC
3056 //              m_BCRegions.resize(0);
3057
3058 //              // Second part
3059 //              m_SecondRegion.SetSize(InputDimension-1,0);
3060
3061 //              // BC2
3062 //              m_BC2Regions.resize(0);
3063
3064 //              break;
3065 //            }
3066 //          case 3:
3067 //            {
3068 //              // Lower end
3069 //              m_FirstRegion.SetSize(InputDimension-1,3);
3070 //              m_FirstRegion.SetIndex(InputDimension-1,1);
3071
3072 //              // BC
3073 //              m_BCRegions.resize(2);
3074 //              m_BCValues.resize(2);
3075 //              m_BCRegion.SetIndex(InputDimension-1,0);
3076 //              m_BCRegions[0]=m_BCRegion;
3077 //              m_BCValues[0]=-m_WeightRatio[2][0];
3078 //              m_BCRegion.SetIndex(InputDimension-1,3);
3079 //              m_BCRegions[1]=m_BCRegion;
3080 //              m_BCValues[1]=-m_WeightRatio[2][0];
3081
3082 //              // Second part
3083 //              m_SecondRegion.SetSize(InputDimension-1,0);
3084
3085 //              // BC2
3086 //              m_BC2Regions.resize(0);
3087
3088 //              break;
3089 //            }
3090
3091 //          case 4:
3092 //            {
3093 //              // Lower end
3094 //              m_FirstRegion.SetSize(InputDimension-1,2);
3095 //              m_FirstRegion.SetIndex(InputDimension-1,2);
3096
3097 //              // BC
3098 //              m_BCRegions.resize(2);
3099 //              m_BCValues.resize(2);
3100 //              m_BCRegion.SetIndex(InputDimension-1,0);
3101 //              m_BCRegions[0]=m_BCRegion;
3102 //              m_BCValues[0]=-m_WeightRatio[2][0];
3103 //              m_BCRegion.SetIndex(InputDimension-1,3);
3104 //              m_BCRegions[1]=m_BCRegion;
3105 //              m_BCValues[1]=-m_WeightRatio[2][0];
3106
3107 //              // Second part
3108 //              m_SecondRegion.SetSize(InputDimension-1,1);
3109 //              m_SecondRegion.SetIndex(InputDimension-1,0);
3110
3111 //              // BC2
3112 //              m_BC2Regions.resize(0);
3113
3114 //              break;
3115 //            }
3116
3117 //          default:
3118 //            {
3119 //              DD("supportindex > 4 ???");
3120 //              DD(m_SupportRegion.GetIndex(InputDimension-1));
3121 //              DD(m_TransformShape);
3122 //            }
3123 //          }// end swith index
3124 //        break;
3125
3126 //      } // end case 1 shape
3127
3128       case 2:
3129         {
3130           // ----------------------------------------------------------------------
3131           // The rabbit with 4 internal CP
3132           // Periodic, constrained to zero at the reference
3133           // at position 3 and the extremes fixed with anit-symm bc
3134           // Coeff      row  BC1  0   1  BC2  2  BC3 BC4
3135           // PaddedCoeff  R:  0   1   2   3   4   5   6
3136           // BC1= 2*R1-R0
3137           // BC2= -weights[2]/weights[0] ( R2+R4 )
3138           // BC3=  R1
3139           // BC4= 2*R1-R4
3140           // ---------------------------------------------------------------------
3141           switch(m_SupportRegion.GetIndex(InputDimension-1))
3142             {
3143             case 0:
3144               {
3145                 // Lower end
3146                 m_FirstRegion.SetSize(InputDimension-1,0);
3147
3148                 // BC
3149                 m_BCRegions.resize(2);
3150                 m_BCValues.resize(2);
3151                 m_BCRegion.SetIndex(InputDimension-1,0);
3152                 m_BCRegions[0]=m_BCRegion;
3153                 m_BCValues[0]=2;
3154                 m_BCRegion.SetIndex(InputDimension-1,1);
3155                 m_BCRegions[1]=m_BCRegion;
3156                 m_BCValues[1]= -1;
3157
3158                 // Second part
3159                 m_SecondRegion.SetSize(InputDimension-1,2);
3160                 m_SecondRegion.SetIndex(InputDimension-1,0);
3161
3162                 // BC
3163                 m_BC2Regions.resize(2);
3164                 m_BC2Values.resize(2);
3165                 m_BCRegion.SetIndex(InputDimension-1,1);
3166                 m_BC2Regions[0]=m_BCRegion;
3167                 m_BC2Values[0]=-m_WeightRatio[2][0];
3168                 m_BCRegion.SetIndex(InputDimension-1,2);
3169                 m_BC2Regions[1]=m_BCRegion;
3170                 m_BC2Values[1]= -m_WeightRatio[2][0];
3171
3172                 break;
3173               }
3174             case 1:
3175               {
3176                 // Lower end
3177                 m_FirstRegion.SetSize(InputDimension-1,2);
3178                 m_FirstRegion.SetIndex(InputDimension-1,0);
3179
3180                 // BC
3181                 m_BCRegions.resize(2);
3182                 m_BCValues.resize(2);
3183                 m_BCRegion.SetIndex(InputDimension-1,1);
3184                 m_BCRegions[0]=m_BCRegion;
3185                 m_BCValues[0]=-m_WeightRatio[2][0];
3186                 m_BCRegion.SetIndex(InputDimension-1,2);
3187                 m_BCRegions[1]=m_BCRegion;
3188                 m_BCValues[1]= -m_WeightRatio[2][0];
3189
3190                 // Second part
3191                 m_SecondRegion.SetSize(InputDimension-1,1);
3192                 m_SecondRegion.SetIndex(InputDimension-1,2);
3193
3194                 // BC2
3195                 m_BC2Regions.resize(0);
3196
3197                 break;
3198               }
3199             case 2:
3200               {
3201                 // Lower end
3202                 m_FirstRegion.SetSize(InputDimension-1,1);
3203                 m_FirstRegion.SetIndex(InputDimension-1,1);
3204
3205                 // BC
3206                 m_BCRegions.resize(2);
3207                 m_BCValues.resize(2);
3208                 m_BCRegion.SetIndex(InputDimension-1,1);
3209                 m_BCRegions[0]=m_BCRegion;
3210                 m_BCValues[0]=-m_WeightRatio[2][0];
3211                 m_BCRegion.SetIndex(InputDimension-1,2);
3212                 m_BCRegions[1]=m_BCRegion;
3213                 m_BCValues[1]= -m_WeightRatio[2][0];
3214
3215                 // Second part
3216                 m_SecondRegion.SetSize(InputDimension-1,1);
3217                 m_SecondRegion.SetIndex(InputDimension-1,2);
3218
3219                 // BC2
3220                 m_BC2Regions.resize(1);
3221                 m_BC2Values.resize(1);
3222                 m_BCRegion.SetIndex(InputDimension-1,0);
3223                 m_BC2Regions[0]=m_BCRegion;
3224                 m_BC2Values[0]=1;
3225
3226                 break;
3227               }
3228             case 3:
3229               {
3230                 // Lower end
3231                 m_FirstRegion.SetSize(InputDimension-1,0);
3232
3233                 // BC
3234                 m_BCRegions.resize(2);
3235                 m_BCValues.resize(2);
3236                 m_BCRegion.SetIndex(InputDimension-1,1);
3237                 m_BCRegions[0]=m_BCRegion;
3238                 m_BCValues[0]=-m_WeightRatio[2][0];
3239                 m_BCRegion.SetIndex(InputDimension-1,2);
3240                 m_BCRegions[1]=m_BCRegion;
3241                 m_BCValues[1]= -m_WeightRatio[2][0];
3242
3243                 // Second part
3244                 m_SecondRegion.SetSize(InputDimension-1,1);
3245                 m_SecondRegion.SetIndex(InputDimension-1,2);
3246
3247                 // BC2
3248                 m_BC2Regions.resize(1);
3249                 m_BC2Values.resize(1);
3250                 m_BCRegion.SetIndex(InputDimension-1,0);
3251                 m_BC2Regions[0]=m_BCRegion;
3252                 m_BC2Values[0]=1;
3253
3254                 // BC3
3255                 m_BC3Regions.resize(2);
3256                 m_BC3Values.resize(2);
3257                 m_BCRegion.SetIndex(InputDimension-1,0);
3258                 m_BC3Regions[0]=m_BCRegion;
3259                 m_BC3Values[0]=2;
3260                 m_BCRegion.SetIndex(InputDimension-1,2);
3261                 m_BC3Regions[1]=m_BCRegion;
3262                 m_BC3Values[1]= -1;
3263
3264                 break;
3265               }
3266
3267             default:
3268               {
3269                 DD("supportindex > 3 ???");
3270                 DD(m_SupportRegion.GetIndex(InputDimension-1));
3271               }
3272             } // end switch index
3273
3274           break;
3275         } // end case 2 shape
3276       case 3:
3277         {
3278           // ----------------------------------------------------------------------
3279           // The rabbit with 5 internal CP
3280           // Periodic, constrained to zero at the reference at position 3.5
3281           // and the extremes fixed with anti-symmetrical BC
3282           // Coeff      row  BC1  0   1  BC2  2   3  BC3 BC4
3283           // PaddedCoeff  R:  0   1   2   3   4   5   6   7
3284           // BC1= 2*R1-R2
3285           // BC2= -weights[3]/weights[1] ( R2+R5 ) - R4
3286           // BC3= R1
3287           // BC4= 2*R1-R5
3288           // ---------------------------------------------------------------------
3289           switch(m_SupportRegion.GetIndex(InputDimension-1))
3290             {
3291             case 0:
3292               {
3293                 // Lower end
3294                 m_FirstRegion.SetSize(InputDimension-1,0);
3295
3296                 // BC
3297                 m_BCRegions.resize(2);
3298                 m_BCValues.resize(2);
3299                 m_BCRegion.SetIndex(InputDimension-1,0);
3300                 m_BCRegions[0]=m_BCRegion;
3301                 m_BCValues[0]=2.;
3302                 m_BCRegion.SetIndex(InputDimension-1,1);
3303                 m_BCRegions[1]=m_BCRegion;
3304                 m_BCValues[1]= -1.;
3305
3306                 // Second part
3307                 m_SecondRegion.SetSize(InputDimension-1,2);
3308                 m_SecondRegion.SetIndex(InputDimension-1,0);
3309
3310                 // BC2
3311                 m_BC2Regions.resize(3);
3312                 m_BC2Values.resize(3);
3313                 m_BCRegion.SetIndex(InputDimension-1,1);
3314                 m_BC2Regions[0]=m_BCRegion;
3315                 m_BC2Values[0]=-m_WeightRatio[3][1];
3316                 m_BCRegion.SetIndex(InputDimension-1,2);
3317                 m_BC2Regions[1]=m_BCRegion;
3318                 m_BC2Values[1]= -1;
3319                 m_BCRegion.SetIndex(InputDimension-1,3);
3320                 m_BC2Regions[2]=m_BCRegion;
3321                 m_BC2Values[2]=-m_WeightRatio[3][1];
3322
3323                 break;
3324               }
3325             case 1:
3326               {
3327                 // Lower end
3328                 m_FirstRegion.SetSize(InputDimension-1,2);
3329                 m_FirstRegion.SetIndex(InputDimension-1,0);
3330
3331                 // BC
3332                 m_BCRegions.resize(3);
3333                 m_BCValues.resize(3);
3334                 m_BCRegion.SetIndex(InputDimension-1,1);
3335                 m_BCRegions[0]=m_BCRegion;
3336                 m_BCValues[0]=-m_WeightRatio[3][1];
3337                 m_BCRegion.SetIndex(InputDimension-1,2);
3338                 m_BCRegions[1]=m_BCRegion;
3339                 m_BCValues[1]= -1;
3340                 m_BCRegion.SetIndex(InputDimension-1,3);
3341                 m_BCRegions[2]=m_BCRegion;
3342                 m_BCValues[2]=-m_WeightRatio[3][1];
3343
3344                 // Second part
3345                 m_SecondRegion.SetSize(InputDimension-1,1);
3346                 m_SecondRegion.SetIndex(InputDimension-1,2);
3347
3348                 // BC2
3349                 m_BC2Regions.resize(0);
3350
3351                 break;
3352               }
3353             case 2:
3354               {
3355                 // Lower end
3356                 m_FirstRegion.SetSize(InputDimension-1,1);
3357                 m_FirstRegion.SetIndex(InputDimension-1,1);
3358
3359                 // BC
3360                 m_BCRegions.resize(3);
3361                 m_BCValues.resize(3);
3362                 m_BCRegion.SetIndex(InputDimension-1,1);
3363                 m_BCRegions[0]=m_BCRegion;
3364                 m_BCValues[0]=-m_WeightRatio[3][1];
3365                 m_BCRegion.SetIndex(InputDimension-1,2);
3366                 m_BCRegions[1]=m_BCRegion;
3367                 m_BCValues[1]= -1;
3368                 m_BCRegion.SetIndex(InputDimension-1,3);
3369                 m_BCRegions[2]=m_BCRegion;
3370                 m_BCValues[2]=-m_WeightRatio[3][1];
3371
3372                 // Second part
3373                 m_SecondRegion.SetSize(InputDimension-1,2);
3374                 m_SecondRegion.SetIndex(InputDimension-1,2);
3375
3376                 // BC2
3377                 m_BC2Regions.resize(0);
3378
3379                 break;
3380               }
3381             case 3:
3382               {
3383                 // Lower end
3384                 m_FirstRegion.SetSize(InputDimension-1,0);
3385
3386                 // BC
3387                 m_BCRegions.resize(3);
3388                 m_BCValues.resize(3);
3389                 m_BCRegion.SetIndex(InputDimension-1,1);
3390                 m_BCRegions[0]=m_BCRegion;
3391                 m_BCValues[0]=-m_WeightRatio[3][1];
3392                 m_BCRegion.SetIndex(InputDimension-1,2);
3393                 m_BCRegions[1]=m_BCRegion;
3394                 m_BCValues[1]= -1;
3395                 m_BCRegion.SetIndex(InputDimension-1,3);
3396                 m_BCRegions[2]=m_BCRegion;
3397                 m_BCValues[2]=-m_WeightRatio[3][1];
3398
3399                 // Second part
3400                 m_SecondRegion.SetSize(InputDimension-1,2);
3401                 m_SecondRegion.SetIndex(InputDimension-1,2);
3402
3403                 // BC2
3404                 m_BC2Regions.resize(1);
3405                 m_BC2Values.resize(1);
3406                 m_BCRegion.SetIndex(InputDimension-1,0);
3407                 m_BC2Regions[0]=m_BCRegion;
3408                 m_BC2Values[0]=1.;
3409
3410                 break;
3411               }
3412
3413             case 4:
3414               {
3415                 // Lower end
3416                 m_FirstRegion.SetSize(InputDimension-1,2);
3417                 m_FirstRegion.SetIndex(InputDimension-1,2);
3418
3419                 // BC2
3420                 m_BCRegions.resize(1);
3421                 m_BCValues.resize(1);
3422                 m_BCRegion.SetIndex(InputDimension-1,0);
3423                 m_BCRegions[0]=m_BCRegion;
3424                 m_BCValues[0]=1.;
3425
3426                 // Second part
3427                 m_SecondRegion.SetSize(InputDimension-1,0);
3428
3429                 // BC2
3430                 m_BC2Regions.resize(2);
3431                 m_BC2Values.resize(2);
3432                 m_BCRegion.SetIndex(InputDimension-1,0);
3433                 m_BC2Regions[0]=m_BCRegion;
3434                 m_BC2Values[0]=2;
3435                 m_BCRegion.SetIndex(InputDimension-1,3);
3436                 m_BC2Regions[1]=m_BCRegion;
3437                 m_BC2Values[1]=-1;
3438
3439                 break;
3440
3441               }
3442             default:
3443               {
3444                 DD("supportindex > 4 ???");
3445                 DD(m_SupportRegion.GetIndex(InputDimension-1));
3446               }
3447             } // end switch index
3448
3449           break;
3450         } // end case 3 shape
3451
3452       case 4:
3453         {
3454           // ----------------------------------------------------------------------
3455           // The sputnik with 4 internal CP
3456           // Periodic, constrained to zero at the reference
3457           // at position 3 and one indepent extremes copied
3458           // Coeff     row BC1  0   1  BC2  2  BC2  3
3459           // PaddedCoeff R: 0   1   2   3   4   5   6
3460           // BC1= R6
3461           // BC2= -weights[2]/weights[0] ( R2+R4 )
3462           // BC3=  weights[2]/weights[0] ( R2-R4) + R1
3463           // ---------------------------------------------------------------------
3464           switch(m_SupportRegion.GetIndex(InputDimension-1))
3465             {
3466             case 0:
3467               {
3468                 // Lower end
3469                 m_FirstRegion.SetSize(InputDimension-1,0);
3470
3471                 // BC
3472                 m_BCRegions.resize(1);
3473                 m_BCValues.resize(1);
3474                 m_BCRegion.SetIndex(InputDimension-1,3);
3475                 m_BCRegions[0]=m_BCRegion;
3476                 m_BCValues[0]=1;
3477
3478                 // Second part
3479                 m_SecondRegion.SetSize(InputDimension-1,2);
3480                 m_SecondRegion.SetIndex(InputDimension-1,0);
3481
3482                 // BC2
3483                 m_BC2Regions.resize(2);
3484                 m_BC2Values.resize(2);
3485                 m_BCRegion.SetIndex(InputDimension-1,1);
3486                 m_BC2Regions[0]=m_BCRegion;
3487                 m_BC2Values[0]=-m_WeightRatio[2][0];
3488                 m_BCRegion.SetIndex(InputDimension-1,2);
3489                 m_BC2Regions[1]=m_BCRegion;
3490                 m_BC2Values[1]= -m_WeightRatio[2][0];
3491
3492                 break;
3493               }
3494             case 1:
3495               {
3496                 // Lower end
3497                 m_FirstRegion.SetSize(InputDimension-1,2);
3498                 m_FirstRegion.SetIndex(InputDimension-1,0);
3499
3500                 // BC
3501                 m_BCRegions.resize(2);
3502                 m_BCValues.resize(2);
3503                 m_BCRegion.SetIndex(InputDimension-1,1);
3504                 m_BCRegions[0]=m_BCRegion;
3505                 m_BCValues[0]=-m_WeightRatio[2][0];
3506                 m_BCRegion.SetIndex(InputDimension-1,2);
3507                 m_BCRegions[1]=m_BCRegion;
3508                 m_BCValues[1]= -m_WeightRatio[2][0];
3509
3510                 // Second part
3511                 m_SecondRegion.SetSize(InputDimension-1,1);
3512                 m_SecondRegion.SetIndex(InputDimension-1,2);
3513
3514                 // BC2
3515                 m_BC2Regions.resize(0);
3516
3517                 break;
3518               }
3519             case 2:
3520               {
3521                 // Lower end
3522                 m_FirstRegion.SetSize(InputDimension-1,1);
3523                 m_FirstRegion.SetIndex(InputDimension-1,1);
3524
3525                 // BC
3526                 m_BCRegions.resize(2);
3527                 m_BCValues.resize(2);
3528                 m_BCRegion.SetIndex(InputDimension-1,1);
3529                 m_BCRegions[0]=m_BCRegion;
3530                 m_BCValues[0]=-m_WeightRatio[2][0];
3531                 m_BCRegion.SetIndex(InputDimension-1,2);
3532                 m_BCRegions[1]=m_BCRegion;
3533                 m_BCValues[1]= -m_WeightRatio[2][0];
3534
3535                 // Second part
3536                 m_SecondRegion.SetSize(InputDimension-1,1);
3537                 m_SecondRegion.SetIndex(InputDimension-1,2);
3538
3539                 // BC2
3540                 m_BC2Regions.resize(3);
3541                 m_BC2Values.resize(3);
3542                 m_BCRegion.SetIndex(InputDimension-1,0);
3543                 m_BC2Regions[0]=m_BCRegion;
3544                 m_BC2Values[0]=1;
3545                 m_BCRegion.SetIndex(InputDimension-1,1);
3546                 m_BC2Regions[1]=m_BCRegion;
3547                 m_BC2Values[1]=m_WeightRatio[2][0];
3548                 m_BCRegion.SetIndex(InputDimension-1,2);
3549                 m_BC2Regions[2]=m_BCRegion;
3550                 m_BC2Values[2]=-m_WeightRatio[2][0];
3551
3552                 break;
3553               }
3554             case 3:
3555               {
3556                 // Lower end
3557                 m_FirstRegion.SetSize(InputDimension-1,0);
3558
3559                 // BC
3560                 m_BCRegions.resize(2);
3561                 m_BCValues.resize(2);
3562                 m_BCRegion.SetIndex(InputDimension-1,1);
3563                 m_BCRegions[0]=m_BCRegion;
3564                 m_BCValues[0]=-m_WeightRatio[2][0];
3565                 m_BCRegion.SetIndex(InputDimension-1,2);
3566                 m_BCRegions[1]=m_BCRegion;
3567                 m_BCValues[1]= -m_WeightRatio[2][0];
3568
3569                 // Second part
3570                 m_SecondRegion.SetSize(InputDimension-1,1);
3571                 m_SecondRegion.SetIndex(InputDimension-1,2);
3572
3573                 // BC2
3574                 m_BC2Regions.resize(3);
3575                 m_BC2Values.resize(3);
3576                 m_BCRegion.SetIndex(InputDimension-1,0);
3577                 m_BC2Regions[0]=m_BCRegion;
3578                 m_BC2Values[0]=1;
3579                 m_BCRegion.SetIndex(InputDimension-1,1);
3580                 m_BC2Regions[1]=m_BCRegion;
3581                 m_BC2Values[1]=m_WeightRatio[2][0];
3582                 m_BCRegion.SetIndex(InputDimension-1,2);
3583                 m_BC2Regions[2]=m_BCRegion;
3584                 m_BC2Values[2]=-m_WeightRatio[2][0];
3585
3586                 // Third part
3587                 m_ThirdRegion.SetSize(InputDimension-1,1);
3588                 m_ThirdRegion.SetIndex(InputDimension-1,3);
3589
3590                 break;
3591               }
3592
3593             default:
3594               {
3595                 DD("supportindex > 3 ???");
3596                 DD(m_SupportRegion.GetIndex(InputDimension-1));
3597               }
3598             } // end switch index
3599
3600           break;
3601         } // end case 4 shape
3602
3603       case 5:
3604         {
3605           // ----------------------------------------------------------------------
3606           // The sputnik with 5 internal CP
3607           // Periodic, constrained to zero at the reference
3608           // at position 3 and one indepent extreme
3609           // Coeff     row BC1  0   1  BC2  2   3  BC3  4
3610           // PaddedCoeff R: 0   1   2   3   4   5   6   7
3611           // BC1= R2 + R5 - R7
3612           // BC2= -weights[3]/weights[1] ( R2+R5 ) - R4
3613           // BC3= R1 + 0.5 R2 - 0.5 R7
3614           // ----------------------------------------------------------------------
3615           switch(m_SupportRegion.GetIndex(InputDimension-1))
3616             {
3617             case 0:
3618               {
3619                 // Lower end
3620                 m_FirstRegion.SetSize(InputDimension-1,0);
3621
3622                 // BC
3623                 m_BCRegions.resize(3);
3624                 m_BCValues.resize(3);
3625                 m_BCRegion.SetIndex(InputDimension-1,1);
3626                 m_BCRegions[0]=m_BCRegion;
3627                 m_BCValues[0]=1.;
3628                 m_BCRegion.SetIndex(InputDimension-1,3);
3629                 m_BCRegions[1]=m_BCRegion;
3630                 m_BCValues[1]=1.;
3631                 m_BCRegion.SetIndex(InputDimension-1,4);
3632                 m_BCRegions[2]=m_BCRegion;
3633                 m_BCValues[2]=-1.;
3634
3635                 // Second part
3636                 m_SecondRegion.SetSize(InputDimension-1,2);
3637                 m_SecondRegion.SetIndex(InputDimension-1,0);
3638
3639                 // BC2
3640                 m_BC2Regions.resize(3);
3641                 m_BC2Values.resize(3);
3642                 m_BCRegion.SetIndex(InputDimension-1,1);
3643                 m_BC2Regions[0]=m_BCRegion;
3644                 m_BC2Values[0]=-m_WeightRatio[3][1];
3645                 m_BCRegion.SetIndex(InputDimension-1,2);
3646                 m_BC2Regions[1]=m_BCRegion;
3647                 m_BC2Values[1]= -1;
3648                 m_BCRegion.SetIndex(InputDimension-1,3);
3649                 m_BC2Regions[2]=m_BCRegion;
3650                 m_BC2Values[2]=-m_WeightRatio[3][1];
3651
3652                 break;
3653               }
3654             case 1:
3655               {
3656                 // Lower end
3657                 m_FirstRegion.SetSize(InputDimension-1,2);
3658                 m_FirstRegion.SetIndex(InputDimension-1,0);
3659
3660                 // BC
3661                 m_BCRegions.resize(3);
3662                 m_BCValues.resize(3);
3663                 m_BCRegion.SetIndex(InputDimension-1,1);
3664                 m_BCRegions[0]=m_BCRegion;
3665                 m_BCValues[0]=-m_WeightRatio[3][1];
3666                 m_BCRegion.SetIndex(InputDimension-1,2);
3667                 m_BCRegions[1]=m_BCRegion;
3668                 m_BCValues[1]= -1;
3669                 m_BCRegion.SetIndex(InputDimension-1,3);
3670                 m_BCRegions[2]=m_BCRegion;
3671                 m_BCValues[2]=-m_WeightRatio[3][1];
3672
3673                 // Second part
3674                 m_SecondRegion.SetSize(InputDimension-1,1);
3675                 m_SecondRegion.SetIndex(InputDimension-1,2);
3676
3677                 // BC2
3678                 m_BC2Regions.resize(0);
3679
3680                 break;
3681               }
3682             case 2:
3683               {
3684                 // Lower end
3685                 m_FirstRegion.SetSize(InputDimension-1,1);
3686                 m_FirstRegion.SetIndex(InputDimension-1,1);
3687
3688                 // BC
3689                 m_BCRegions.resize(3);
3690                 m_BCValues.resize(3);
3691                 m_BCRegion.SetIndex(InputDimension-1,1);
3692                 m_BCRegions[0]=m_BCRegion;
3693                 m_BCValues[0]=-m_WeightRatio[3][1];
3694                 m_BCRegion.SetIndex(InputDimension-1,2);
3695                 m_BCRegions[1]=m_BCRegion;
3696                 m_BCValues[1]= -1;
3697                 m_BCRegion.SetIndex(InputDimension-1,3);
3698                 m_BCRegions[2]=m_BCRegion;
3699                 m_BCValues[2]=-m_WeightRatio[3][1];
3700
3701                 // Second part
3702                 m_SecondRegion.SetSize(InputDimension-1,2);
3703                 m_SecondRegion.SetIndex(InputDimension-1,2);
3704
3705                 // BC2
3706                 m_BC2Regions.resize(0);
3707
3708                 break;
3709               }
3710             case 3:
3711               {
3712                 // Lower end
3713                 m_FirstRegion.SetSize(InputDimension-1,0);
3714
3715                 // BC
3716                 m_BCRegions.resize(3);
3717                 m_BCValues.resize(3);
3718                 m_BCRegion.SetIndex(InputDimension-1,1);
3719                 m_BCRegions[0]=m_BCRegion;
3720                 m_BCValues[0]=-m_WeightRatio[3][1];
3721                 m_BCRegion.SetIndex(InputDimension-1,2);
3722                 m_BCRegions[1]=m_BCRegion;
3723                 m_BCValues[1]= -1;
3724                 m_BCRegion.SetIndex(InputDimension-1,3);
3725                 m_BCRegions[2]=m_BCRegion;
3726                 m_BCValues[2]=-m_WeightRatio[3][1];
3727
3728                 // Second part
3729                 m_SecondRegion.SetSize(InputDimension-1,2);
3730                 m_SecondRegion.SetIndex(InputDimension-1,2);
3731
3732                 // BC2
3733                 m_BC2Regions.resize(3);
3734                 m_BC2Values.resize(3);
3735                 m_BCRegion.SetIndex(InputDimension-1,0);
3736                 m_BC2Regions[0]=m_BCRegion;
3737                 m_BC2Values[0]=1;
3738                 m_BCRegion.SetIndex(InputDimension-1,1);
3739                 m_BC2Regions[1]=m_BCRegion;
3740                 m_BC2Values[1]=0.5;
3741                 m_BCRegion.SetIndex(InputDimension-1,4);
3742                 m_BC2Regions[2]=m_BCRegion;
3743                 m_BC2Values[2]=-0.5;
3744
3745                 break;
3746               }
3747             case 4:
3748               {
3749                 // Lower end
3750                 m_FirstRegion.SetSize(InputDimension-1,2);
3751                 m_FirstRegion.SetIndex(InputDimension-1,2);
3752
3753                 // BC
3754                 m_BCRegions.resize(3);
3755                 m_BCValues.resize(3);
3756                 m_BCRegion.SetIndex(InputDimension-1,0);
3757                 m_BCRegions[0]=m_BCRegion;
3758                 m_BCValues[0]=1;
3759                 m_BCRegion.SetIndex(InputDimension-1,1);
3760                 m_BCRegions[1]=m_BCRegion;
3761                 m_BCValues[1]=0.5;
3762                 m_BCRegion.SetIndex(InputDimension-1,4);
3763                 m_BCRegions[2]=m_BCRegion;
3764                 m_BCValues[2]=-0.5;
3765
3766                 // Second part
3767                 m_SecondRegion.SetSize(InputDimension-1,1);
3768                 m_SecondRegion.SetIndex(InputDimension-1,4);
3769
3770                 // BC2
3771                 m_BC2Regions.resize(0);
3772
3773                 break;
3774               }
3775             default:
3776               {
3777                 DD("supportindex > 4 ???");
3778                 DD(m_SupportRegion.GetIndex(InputDimension-1));
3779               }
3780             } // end switch index
3781
3782           break;
3783         } // end case 5 shape
3784
3785
3786
3787
3788       case 6:
3789         {
3790           // ----------------------------------------------------------------------
3791           // The diamond with 4 internal CP:
3792           // Periodic, constrained to zero at the reference at position 3
3793           // Coeff      row 0   1   2  BC1  3  BC2  4
3794           // PaddedCoeff  R:0   1   2   3   4   5   6
3795           // BC1= -weights[2]/weights[0] ( R2+R4 )
3796           // BC2= weights[2]/weights[0]  ( R0+R2-R4-R6 ) + R1
3797           // ---------------------------------------------------------------------
3798           switch(m_SupportRegion.GetIndex(InputDimension-1))
3799             {
3800             case 0:
3801               {
3802                 // Lower end
3803                 m_FirstRegion.SetSize(InputDimension-1,3);
3804                 m_FirstRegion.SetIndex(InputDimension-1,0);
3805
3806                 // BC
3807                 m_BCRegions.resize(2);
3808                 m_BCValues.resize(2);
3809                 m_BCRegion.SetIndex(InputDimension-1,2);
3810                 m_BCRegions[0]=m_BCRegion;
3811                 m_BCValues[0]=-m_WeightRatio[2][0];
3812                 m_BCRegion.SetIndex(InputDimension-1,3);
3813                 m_BCRegions[1]=m_BCRegion;
3814                 m_BCValues[1]=-m_WeightRatio[2][0];
3815
3816                 // Second part
3817                 m_SecondRegion.SetSize(InputDimension-1,0);
3818
3819                 // BC2
3820                 m_BC2Regions.resize(0);
3821
3822                 break;
3823               }
3824             case 1:
3825               {
3826                 // Lower end
3827                 m_FirstRegion.SetSize(InputDimension-1,2);
3828                 m_FirstRegion.SetIndex(InputDimension-1,1);
3829
3830                 // BC
3831                 m_BCRegions.resize(2);
3832                 m_BCValues.resize(2);
3833                 m_BCRegion.SetIndex(InputDimension-1,2);
3834                 m_BCRegions[0]=m_BCRegion;
3835                 m_BCValues[0]=-m_WeightRatio[2][0];
3836                 m_BCRegion.SetIndex(InputDimension-1,3);
3837                 m_BCRegions[1]=m_BCRegion;
3838                 m_BCValues[1]=-m_WeightRatio[2][0];
3839
3840                 // Second part
3841                 m_SecondRegion.SetSize(InputDimension-1,1);
3842                 m_SecondRegion.SetIndex(InputDimension-1,3);
3843
3844                 // BC2
3845                 m_BC2Regions.resize(0);
3846
3847                 break;
3848               }
3849             case 2:
3850               {
3851                 // Lower end
3852                 m_FirstRegion.SetSize(InputDimension-1,1);
3853                 m_FirstRegion.SetIndex(InputDimension-1,2);
3854
3855                 // BC
3856                 m_BCRegions.resize(2);
3857                 m_BCValues.resize(2);
3858                 m_BCRegion.SetIndex(InputDimension-1,2);
3859                 m_BCRegions[0]=m_BCRegion;
3860                 m_BCValues[0]=-m_WeightRatio[2][0];
3861                 m_BCRegion.SetIndex(InputDimension-1,3);
3862                 m_BCRegions[1]=m_BCRegion;
3863                 m_BCValues[1]=-m_WeightRatio[2][0];
3864
3865                 // Second part
3866                 m_SecondRegion.SetSize(InputDimension-1,1);
3867                 m_SecondRegion.SetIndex(InputDimension-1,3);
3868
3869                 // BC2
3870                 m_BC2Regions.resize(5);
3871                 m_BC2Values.resize(5);
3872                 m_BCRegion.SetIndex(InputDimension-1,0);
3873                 m_BC2Regions[0]=m_BCRegion;
3874                 m_BC2Values[0]=m_WeightRatio[2][0];
3875                 m_BCRegion.SetIndex(InputDimension-1,1);
3876                 m_BC2Regions[1]=m_BCRegion;
3877                 m_BC2Values[1]=1;
3878                 m_BCRegion.SetIndex(InputDimension-1,2);
3879                 m_BC2Regions[2]=m_BCRegion;
3880                 m_BC2Values[2]=m_WeightRatio[2][0];
3881                 m_BCRegion.SetIndex(InputDimension-1,3);
3882                 m_BC2Regions[3]=m_BCRegion;
3883                 m_BC2Values[3]=-m_WeightRatio[2][0];
3884                 m_BCRegion.SetIndex(InputDimension-1,4);
3885                 m_BC2Regions[4]=m_BCRegion;
3886                 m_BC2Values[4]=-m_WeightRatio[2][0];
3887
3888                 break;
3889               }
3890             case 3:
3891               {
3892                 // Lower end
3893                 m_FirstRegion.SetSize(InputDimension-1,0);
3894
3895                 // BC
3896                 m_BCRegions.resize(2);
3897                 m_BCValues.resize(2);
3898                 m_BCRegion.SetIndex(InputDimension-1,2);
3899                 m_BCRegions[0]=m_BCRegion;
3900                 m_BCValues[0]=-m_WeightRatio[2][0];
3901                 m_BCRegion.SetIndex(InputDimension-1,3);
3902                 m_BCRegions[1]=m_BCRegion;
3903                 m_BCValues[1]=-m_WeightRatio[2][0];
3904
3905                 // Second part
3906                 m_SecondRegion.SetSize(InputDimension-1,1);
3907                 m_SecondRegion.SetIndex(InputDimension-1,3);
3908
3909                 // BC2
3910                 m_BC2Regions.resize(5);
3911                 m_BC2Values.resize(5);
3912                 m_BCRegion.SetIndex(InputDimension-1,0);
3913                 m_BC2Regions[0]=m_BCRegion;
3914                 m_BC2Values[0]=m_WeightRatio[2][0];
3915                 m_BCRegion.SetIndex(InputDimension-1,1);
3916                 m_BC2Regions[1]=m_BCRegion;
3917                 m_BC2Values[1]=1;
3918                 m_BCRegion.SetIndex(InputDimension-1,2);
3919                 m_BC2Regions[2]=m_BCRegion;
3920                 m_BC2Values[2]=m_WeightRatio[2][0];
3921                 m_BCRegion.SetIndex(InputDimension-1,3);
3922                 m_BC2Regions[3]=m_BCRegion;
3923                 m_BC2Values[3]=-m_WeightRatio[2][0];
3924                 m_BCRegion.SetIndex(InputDimension-1,4);
3925                 m_BC2Regions[4]=m_BCRegion;
3926                 m_BC2Values[4]=-m_WeightRatio[2][0];
3927
3928                 // Third part
3929                 m_ThirdRegion.SetSize(InputDimension-1,1);
3930                 m_ThirdRegion.SetIndex(InputDimension-1,4);
3931
3932                 break;
3933               }
3934
3935             default:
3936               {
3937                 DD("supportindex > 3 ???");
3938                 DD(m_SupportRegion.GetIndex(InputDimension-1));
3939               }
3940             } // end switch index
3941
3942           break;
3943         } // end case 7 shape
3944
3945       case 7:
3946         {
3947           // ----------------------------------------------------------------------
3948           // The diamond with 5 internal CP:
3949           // periodic, constrained to zero at the reference at position 3.5
3950           // Coeff      row 0   1   2  BC1  3   4  BC2  5
3951           // PaddedCoeff  R:0   1   2   3   4   5   6   7
3952           // BC1= -weights[3]/weights[1] ( R2+R5 ) - R4
3953           // BC2= weights[2]/weights[0] ( R0+R2-R5-R7 ) + R1
3954           // ---------------------------------------------------------------------
3955           switch(m_SupportRegion.GetIndex(InputDimension-1))
3956             {
3957             case 0:
3958               {
3959                 // Lower end
3960                 m_FirstRegion.SetSize(InputDimension-1,3);
3961                 m_FirstRegion.SetIndex(InputDimension-1,0);
3962
3963                 // BC
3964                 m_BCRegions.resize(3);
3965                 m_BCValues.resize(3);
3966                 m_BCRegion.SetIndex(InputDimension-1,2);
3967                 m_BCRegions[0]=m_BCRegion;
3968                 m_BCValues[0]=-m_WeightRatio[3][1];
3969                 m_BCRegion.SetIndex(InputDimension-1,3);
3970                 m_BCRegions[1]=m_BCRegion;
3971                 m_BCValues[1]= -1;
3972                 m_BCRegion.SetIndex(InputDimension-1,4);
3973                 m_BCRegions[2]=m_BCRegion;
3974                 m_BCValues[2]=-m_WeightRatio[3][1];
3975
3976                 // Second part
3977                 m_SecondRegion.SetSize(InputDimension-1,0);
3978
3979                 // BC2
3980                 m_BC2Regions.resize(0);
3981
3982                 break;
3983               }
3984             case 1:
3985               {
3986                 // Lower end
3987                 m_FirstRegion.SetSize(InputDimension-1,2);
3988                 m_FirstRegion.SetIndex(InputDimension-1,1);
3989
3990                 // BC
3991                 m_BCRegions.resize(3);
3992                 m_BCValues.resize(3);
3993                 m_BCRegion.SetIndex(InputDimension-1,2);
3994                 m_BCRegions[0]=m_BCRegion;
3995                 m_BCValues[0]=-m_WeightRatio[3][1];
3996                 m_BCRegion.SetIndex(InputDimension-1,3);
3997                 m_BCRegions[1]=m_BCRegion;
3998                 m_BCValues[1]= -1;
3999                 m_BCRegion.SetIndex(InputDimension-1,4);
4000                 m_BCRegions[2]=m_BCRegion;
4001                 m_BCValues[2]=-m_WeightRatio[3][1];
4002
4003                 // Second part
4004                 m_SecondRegion.SetSize(InputDimension-1,1);
4005                 m_SecondRegion.SetIndex(InputDimension-1,3);
4006
4007                 // BC2
4008                 m_BC2Regions.resize(0);
4009
4010                 break;
4011               }
4012             case 2:
4013               {
4014                 // Lower end
4015                 m_FirstRegion.SetSize(InputDimension-1,1);
4016                 m_FirstRegion.SetIndex(InputDimension-1,2);
4017
4018                 // BC
4019                 m_BCRegions.resize(3);
4020                 m_BCValues.resize(3);
4021                 m_BCRegion.SetIndex(InputDimension-1,2);
4022                 m_BCRegions[0]=m_BCRegion;
4023                 m_BCValues[0]=-m_WeightRatio[3][1];
4024                 m_BCRegion.SetIndex(InputDimension-1,3);
4025                 m_BCRegions[1]=m_BCRegion;
4026                 m_BCValues[1]= -1;
4027                 m_BCRegion.SetIndex(InputDimension-1,4);
4028                 m_BCRegions[2]=m_BCRegion;
4029                 m_BCValues[2]=-m_WeightRatio[3][1];
4030
4031                 // Second part
4032                 m_SecondRegion.SetSize(InputDimension-1,2);
4033                 m_SecondRegion.SetIndex(InputDimension-1,3);
4034
4035                 // BC2
4036                 m_BC2Regions.resize(0);
4037
4038                 break;
4039               }
4040             case 3:
4041               {
4042                 // Lower end
4043                 m_FirstRegion.SetSize(InputDimension-1,0);
4044                 m_FirstRegion.SetIndex(InputDimension-1,0);
4045
4046                 // BC
4047                 m_BCRegions.resize(3);
4048                 m_BCValues.resize(3);
4049                 m_BCRegion.SetIndex(InputDimension-1,2);
4050                 m_BCRegions[0]=m_BCRegion;
4051                 m_BCValues[0]=-m_WeightRatio[3][1];
4052                 m_BCRegion.SetIndex(InputDimension-1,3);
4053                 m_BCRegions[1]=m_BCRegion;
4054                 m_BCValues[1]= -1;
4055                 m_BCRegion.SetIndex(InputDimension-1,4);
4056                 m_BCRegions[2]=m_BCRegion;
4057                 m_BCValues[2]=-m_WeightRatio[3][1];
4058
4059                 // Second part
4060                 m_SecondRegion.SetSize(InputDimension-1,2);
4061                 m_SecondRegion.SetIndex(InputDimension-1,3);
4062
4063                 // BC2
4064                 m_BC2Regions.resize(5);
4065                 m_BC2Values.resize(5);
4066                 m_BCRegion.SetIndex(InputDimension-1,0);
4067                 m_BC2Regions[0]=m_BCRegion;
4068                 m_BC2Values[0]=m_WeightRatio[2][0];
4069                 m_BCRegion.SetIndex(InputDimension-1,1);
4070                 m_BC2Regions[1]=m_BCRegion;
4071                 m_BC2Values[1]=1;
4072                 m_BCRegion.SetIndex(InputDimension-1,2);
4073                 m_BC2Regions[2]=m_BCRegion;
4074                 m_BC2Values[2]=m_WeightRatio[2][0];
4075                 m_BCRegion.SetIndex(InputDimension-1,4);
4076                 m_BC2Regions[3]=m_BCRegion;
4077                 m_BC2Values[3]=-m_WeightRatio[2][0];
4078                 m_BCRegion.SetIndex(InputDimension-1,5);
4079                 m_BC2Regions[4]=m_BCRegion;
4080                 m_BC2Values[4]=-m_WeightRatio[2][0];
4081
4082                 break;
4083               }
4084
4085             case 4:
4086               {
4087                 // Lower end
4088                 m_FirstRegion.SetSize(InputDimension-1,2);
4089                 m_FirstRegion.SetIndex(InputDimension-1,3);
4090
4091                 // BC
4092                 m_BCRegions.resize(5);
4093                 m_BCValues.resize(5);
4094                 m_BCRegion.SetIndex(InputDimension-1,0);
4095                 m_BCRegions[0]=m_BCRegion;
4096                 m_BCValues[0]=m_WeightRatio[2][0];
4097                 m_BCRegion.SetIndex(InputDimension-1,1);
4098                 m_BCRegions[1]=m_BCRegion;
4099                 m_BCValues[1]=1;
4100                 m_BCRegion.SetIndex(InputDimension-1,2);
4101                 m_BCRegions[2]=m_BCRegion;
4102                 m_BCValues[2]=m_WeightRatio[2][0];
4103                 m_BCRegion.SetIndex(InputDimension-1,4);
4104                 m_BCRegions[3]=m_BCRegion;
4105                 m_BCValues[3]=-m_WeightRatio[2][0];
4106                 m_BCRegion.SetIndex(InputDimension-1,5);
4107                 m_BCRegions[4]=m_BCRegion;
4108                 m_BCValues[4]=-m_WeightRatio[2][0];
4109
4110                 // Second part
4111                 m_SecondRegion.SetSize(InputDimension-1,1);
4112                 m_SecondRegion.SetIndex(InputDimension-1,5);
4113
4114                 // BC2
4115                 m_BC2Regions.resize(0);
4116
4117                 break;
4118               }
4119
4120             default:
4121               {
4122                 DD("supportindex > 4 ???");
4123                 DD(m_SupportRegion.GetIndex(InputDimension-1));
4124               }
4125             } // end switch index
4126
4127           break;
4128         } // end case 7 shape
4129
4130       case 9:
4131         {
4132           // ----------------------------------------------------------------------
4133           // The sputnik with 5 internal CP
4134           // Periodic, constrained to zero at the reference
4135           // at position 3 and one indepent extreme
4136           // Coeff     row BC1  0   1  BC2  2   3  BC3  4
4137           // PaddedCoeff R: 0   1   2   3   4   5   6   7
4138           // BC1= -R2+R5-R7
4139           // BC2= -weights[3]/weights[1] ( R2+R5 ) - R4
4140           // BC3= R1
4141           // ---------------------------------------------------------------------
4142           switch(m_SupportRegion.GetIndex(InputDimension-1))
4143             {
4144             case 0:
4145               {
4146                 // Lower end
4147                 m_FirstRegion.SetSize(InputDimension-1,0);
4148
4149                 // BC
4150                 m_BCRegions.resize(3);
4151                 m_BCValues.resize(3);
4152                 m_BCRegion.SetIndex(InputDimension-1,1);
4153                 m_BCRegions[0]=m_BCRegion;
4154                 m_BCValues[0]=-1;
4155                 m_BCRegion.SetIndex(InputDimension-1,3);
4156                 m_BCRegions[1]=m_BCRegion;
4157                 m_BCValues[1]=1;
4158                 m_BCRegion.SetIndex(InputDimension-1,4);
4159                 m_BCRegions[2]=m_BCRegion;
4160                 m_BCValues[2]=-1;
4161
4162                 // Second part
4163                 m_SecondRegion.SetSize(InputDimension-1,2);
4164                 m_SecondRegion.SetIndex(InputDimension-1,0);
4165
4166                 // BC2
4167                 m_BC2Regions.resize(3);
4168                 m_BC2Values.resize(3);
4169                 m_BCRegion.SetIndex(InputDimension-1,1);
4170                 m_BC2Regions[0]=m_BCRegion;
4171                 m_BC2Values[0]=-m_WeightRatio[3][1];
4172                 m_BCRegion.SetIndex(InputDimension-1,2);
4173                 m_BC2Regions[1]=m_BCRegion;
4174                 m_BC2Values[1]= -1;
4175                 m_BCRegion.SetIndex(InputDimension-1,3);
4176                 m_BC2Regions[2]=m_BCRegion;
4177                 m_BC2Values[2]=-m_WeightRatio[3][1];
4178
4179                 break;
4180               }
4181             case 1:
4182               {
4183                 // Lower end
4184                 m_FirstRegion.SetSize(InputDimension-1,2);
4185                 m_FirstRegion.SetIndex(InputDimension-1,0);
4186
4187                 // BC
4188                 m_BCRegions.resize(3);
4189                 m_BCValues.resize(3);
4190                 m_BCRegion.SetIndex(InputDimension-1,1);
4191                 m_BCRegions[0]=m_BCRegion;
4192                 m_BCValues[0]=-m_WeightRatio[3][1];
4193                 m_BCRegion.SetIndex(InputDimension-1,2);
4194                 m_BCRegions[1]=m_BCRegion;
4195                 m_BCValues[1]= -1;
4196                 m_BCRegion.SetIndex(InputDimension-1,3);
4197                 m_BCRegions[2]=m_BCRegion;
4198                 m_BCValues[2]=-m_WeightRatio[3][1];
4199
4200                 // Second part
4201                 m_SecondRegion.SetSize(InputDimension-1,1);
4202                 m_SecondRegion.SetIndex(InputDimension-1,2);
4203
4204                 // BC2
4205                 m_BC2Regions.resize(0);
4206
4207                 break;
4208               }
4209             case 2:
4210               {
4211                 // Lower end
4212                 m_FirstRegion.SetSize(InputDimension-1,1);
4213                 m_FirstRegion.SetIndex(InputDimension-1,1);
4214
4215                 // BC
4216                 m_BCRegions.resize(3);
4217                 m_BCValues.resize(3);
4218                 m_BCRegion.SetIndex(InputDimension-1,1);
4219                 m_BCRegions[0]=m_BCRegion;
4220                 m_BCValues[0]=-m_WeightRatio[3][1];
4221                 m_BCRegion.SetIndex(InputDimension-1,2);
4222                 m_BCRegions[1]=m_BCRegion;
4223                 m_BCValues[1]= -1;
4224                 m_BCRegion.SetIndex(InputDimension-1,3);
4225                 m_BCRegions[2]=m_BCRegion;
4226                 m_BCValues[2]=-m_WeightRatio[3][1];
4227
4228                 // Second part
4229                 m_SecondRegion.SetSize(InputDimension-1,2);
4230                 m_SecondRegion.SetIndex(InputDimension-1,2);
4231
4232                 // BC2
4233                 m_BC2Regions.resize(0);
4234
4235                 break;
4236               }
4237             case 3:
4238               {
4239                 // Lower end
4240                 m_FirstRegion.SetSize(InputDimension-1,0);
4241
4242                 // BC
4243                 m_BCRegions.resize(3);
4244                 m_BCValues.resize(3);
4245                 m_BCRegion.SetIndex(InputDimension-1,1);
4246                 m_BCRegions[0]=m_BCRegion;
4247                 m_BCValues[0]=-m_WeightRatio[3][1];
4248                 m_BCRegion.SetIndex(InputDimension-1,2);
4249                 m_BCRegions[1]=m_BCRegion;
4250                 m_BCValues[1]= -1;
4251                 m_BCRegion.SetIndex(InputDimension-1,3);
4252                 m_BCRegions[2]=m_BCRegion;
4253                 m_BCValues[2]=-m_WeightRatio[3][1];
4254
4255                 // Second part
4256                 m_SecondRegion.SetSize(InputDimension-1,1);
4257                 m_SecondRegion.SetIndex(InputDimension-1,2);
4258
4259                 // BC
4260                 m_BC2Regions.resize(1);
4261                 m_BC2Values.resize(1);
4262                 m_BCRegion.SetIndex(InputDimension-1,0);
4263                 m_BC2Regions[0]=m_BCRegion;
4264                 m_BC2Values[0]=1;
4265
4266                 break;
4267               }
4268             case 4:
4269               {
4270                 // Lower end
4271                 m_FirstRegion.SetSize(InputDimension-1,2);
4272                 m_FirstRegion.SetIndex(InputDimension-1,2);
4273
4274                 // BC
4275                 m_BCRegions.resize(1);
4276                 m_BCValues.resize(1);
4277                 m_BCRegion.SetIndex(InputDimension-1,0);
4278                 m_BCRegions[0]=m_BCRegion;
4279                 m_BCValues[0]=1;
4280
4281                 // Second part
4282                 m_SecondRegion.SetSize(InputDimension-1,1);
4283                 m_SecondRegion.SetIndex(InputDimension-1,4);
4284
4285                 // BC2
4286                 m_BC2Regions.resize(0);
4287
4288                 break;
4289               }
4290             default:
4291               {
4292                 DD("supportindex > 4 ???");
4293                 DD(m_SupportRegion.GetIndex(InputDimension-1));
4294               }
4295             } // end switch index
4296
4297           break;
4298         } // end case 9 shape
4299
4300
4301       default:
4302         {
4303           DD ("Other shapes currently not implemented");
4304         }
4305
4306       } // end switch shape
4307   } // end wrap region
4308
4309
4310
4311   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
4312   void
4313   ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions,NOutputDimensions>
4314   ::TransformPointToContinuousIndex( const InputPointType & point, ContinuousIndexType & index ) const
4315   {
4316     unsigned int j;
4317
4318     itk::Vector<double, OutputDimension> tvector;
4319
4320     for ( j = 0; j < OutputDimension; j++ )
4321       {
4322         //JV find the index in the PADDED version
4323         tvector[j] = point[j] - this->m_PaddedGridOrigin[j];
4324       }
4325
4326     itk::Vector<double, OutputDimension> cvector;
4327
4328     cvector = m_PointToIndex * tvector;
4329
4330     for ( j = 0; j < OutputDimension; j++ )
4331       {
4332         index[j] = static_cast< typename ContinuousIndexType::CoordRepType >( cvector[j] );
4333       }
4334   }
4335
4336
4337 } // namespace
4338
4339 #endif
4340