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