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