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