registration/clitkShapedBLUTSpatioTemporalDeformableTransform.txx

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