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