1 /*=========================================================================
2 Program: vv http://www.creatis.insa-lyon.fr/rio/vv
5 - University of LYON http://www.universite-lyon.fr/
6 - Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr
7 - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr
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.
13 It is distributed under dual licence
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"
23 #include "itkContinuousIndex.h"
24 #include "itkImageRegionIterator.h"
25 #include "itkImageRegionConstIterator.h"
26 #include "itkImageRegionConstIteratorWithIndex.h"
27 #include "itkIdentityTransform.h"
32 // Constructor with default arguments
33 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
34 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
35 ::ShapedBLUTSpatioTemporalDeformableTransform():Superclass(OutputDimension,0)
39 //JV default spline order
40 for ( i=0;i<InputDimension; i++)
43 //JV default LUTSamplingfactor
44 for ( i=0;i<InputDimension; i++)
45 m_LUTSamplingFactors[i]=20;
47 for (i=0;i<InputDimension; i++)
48 m_SupportSize[i] = m_SplineOrders[i]+1;
50 //Instantiate a Vector Bspline Interpolator
51 m_VectorInterpolator =VectorInterpolatorType::New();
52 m_VectorInterpolator->SetLUTSamplingFactors(m_LUTSamplingFactors);
53 m_VectorInterpolator->SetSplineOrders(m_SplineOrders);
55 // Set Bulk transform to NULL (no bulk performed)
56 m_BulkTransform = NULL;
64 // Default grid size is zero
68 //JV region containing the parameters
69 m_GridRegion.SetSize( m_NullSize);
70 m_GridRegion.SetIndex( m_NullIndex );
72 //JV use second region over the images
73 m_PaddedGridRegion.SetSize(m_NullSize);
74 m_PaddedGridRegion.SetIndex(m_NullIndex);
76 //JV Maintain two origins
77 m_GridOrigin.Fill( 0.0 ); // default origin is all zeros
78 m_PaddedGridOrigin.Fill( 0.0 ); // default origin is all zeros
80 m_GridSpacing.Fill( 1.0 ); // default spacing is all ones
81 m_GridDirection.SetIdentity(); // default spacing is all ones
83 m_InternalParametersBuffer = ParametersType(0);
84 // Make sure the parameters pointer is not NULL after construction.
85 m_InputParametersPointer = &m_InternalParametersBuffer;
87 // Initialize coeffient images
88 m_WrappedImage = CoefficientImageType::New();
89 m_WrappedImage->SetRegions( m_GridRegion );
90 m_WrappedImage->SetOrigin( m_GridOrigin.GetDataPointer() );
91 m_WrappedImage->SetSpacing( m_GridSpacing.GetDataPointer() );
92 m_WrappedImage->SetDirection( m_GridDirection );
94 // JV Initialize the padded version
95 m_PaddedCoefficientImage = CoefficientImageType::New();
96 m_PaddedCoefficientImage->SetRegions( m_PaddedGridRegion );
97 m_PaddedCoefficientImage->SetOrigin( m_GridOrigin.GetDataPointer() );
98 m_PaddedCoefficientImage->SetSpacing( m_GridSpacing.GetDataPointer() );
99 m_PaddedCoefficientImage->SetDirection( m_GridDirection );
101 m_CoefficientImage = NULL;
103 // Variables for computing interpolation
104 for (i=0; i <InputDimension;i++)
106 m_Offset[i] = m_SplineOrders[i] / 2;
107 if ( m_SplineOrders[i] % 2 )
109 m_SplineOrderOdd[i] = true;
113 m_SplineOrderOdd[i] = false;
117 //JV padded should be used when checking
118 m_ValidRegion = m_PaddedGridRegion;
120 // Initialize jacobian images
121 for (i=0; i <OutputDimension;i++)
123 m_JacobianImage[i] = JacobianImageType::New();
124 m_JacobianImage[i]->SetRegions( m_GridRegion );
125 m_JacobianImage[i]->SetOrigin( m_GridOrigin.GetDataPointer() );
126 m_JacobianImage[i]->SetSpacing( m_GridSpacing.GetDataPointer() );
127 m_JacobianImage[i]->SetDirection( m_GridDirection );
130 /** Fixed Parameters store the following information:
135 //JV we add the splineOrders, LUTsamplingfactor, m_Mask and m_BulkTransform
141 The size of these is equal to the NInputDimensions
142 *********************************************************/
143 this->m_FixedParameters.SetSize ( NInputDimensions * (NInputDimensions + 5)+3 );
144 this->m_FixedParameters.Fill ( 0.0 );
145 for ( i=0; i<NInputDimensions; i++)
147 this->m_FixedParameters[2*NInputDimensions+i] = m_GridSpacing[i];
149 for (unsigned int di=0; di<NInputDimensions; di++)
151 for (unsigned int dj=0; dj<NInputDimensions; dj++)
153 this->m_FixedParameters[3*NInputDimensions+(di*NInputDimensions+dj)] = m_GridDirection[di][dj];
157 //JV add splineOrders
158 for ( i=0; i<NInputDimensions; i++)
160 this->m_FixedParameters[ ( (3+NInputDimensions)*NInputDimensions)+i] = (this->GetSplineOrders())[i];
163 //JV add LUTsamplingFactors
164 for ( i=0; i<NInputDimensions; i++)
166 this->m_FixedParameters[( (4+NInputDimensions)*NInputDimensions)+i ] = this->GetLUTSamplingFactors()[i];
169 // JV add the mask pointer
170 this->m_FixedParameters[( (5+NInputDimensions)*NInputDimensions)]=(double)((size_t)m_Mask);
172 // JV add the bulkTransform pointer
173 this->m_FixedParameters[( (5+NInputDimensions)*NInputDimensions) +1]=(double)((size_t)m_BulkTransform);
175 // JV add the Transform shape
176 this->m_FixedParameters[( (5+NInputDimensions)*NInputDimensions) +2]=(double)(m_TransformShape);
179 // Calculate the PointToIndex matrices
181 for( unsigned int i=0; i<OutputDimension; i++)
183 scale[i][i] = m_GridSpacing[i];
186 m_IndexToPoint = m_GridDirection * scale;
187 m_PointToIndex = m_IndexToPoint.GetInverse();
189 m_LastJacobianIndex = m_ValidRegion.GetIndex();
192 //Weights to be used for the BC
194 std::vector<double> weights(5);
203 m_Weights[0]=weights;
211 m_Weights[1]=weights;
219 m_Weights[2]=weights;
227 m_Weights[3]=weights;
229 // Update the WeightRatios
230 m_WeightRatio.resize(4);
231 for (unsigned int i=0; i<4; i++)
233 m_WeightRatio[i].resize(4);
234 for (unsigned int j=0; j<4; j++)
235 m_WeightRatio[i][j]=m_Weights[m_SplineOrders[InputDimension-1]][i]/ m_Weights[m_SplineOrders[InputDimension-1]][j];
238 //JV initialize some variables for jacobian calculation
239 m_SupportRegion.SetSize(m_SupportSize);
240 m_SupportIndex.Fill(0);
241 m_SupportRegion.SetIndex(m_SupportIndex);
242 for ( i = 0; i < SpaceDimension ; i++ )
243 m_ZeroVector[i]=itk::NumericTraits<JacobianValueType>::Zero;
246 m_FirstRegion.SetSize(m_NullSize);
247 m_FirstRegion.SetIndex(m_NullIndex);
248 m_SecondRegion.SetSize(m_NullSize);
249 m_SecondRegion.SetIndex(m_NullIndex);
250 m_ThirdRegion.SetSize(m_NullSize);
251 m_ThirdRegion.SetIndex(m_NullIndex);
253 m_BCValues.resize(0);
254 m_BCRegions.resize(0);
256 m_BC2Values.resize(0);
257 m_BC2Regions.resize(0);
259 m_BC3Values.resize(0);
260 m_BC3Regions.resize(0);
264 for ( i = 0; i < SpaceDimension ; i++ )
266 m_FirstIterator[i]= IteratorType( m_JacobianImage[i], m_FirstRegion);
267 m_SecondIterator[i]= IteratorType( m_JacobianImage[i], m_SecondRegion);
268 m_ThirdIterator[i]= IteratorType( m_JacobianImage[i], m_ThirdRegion);
269 m_BCIterators[i].resize(0);
270 m_BC2Iterators[i].resize(0);
271 m_BC3Iterators[i].resize(0);
280 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
281 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
282 ::~ShapedBLUTSpatioTemporalDeformableTransform()
288 // JV set Spline Order
289 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
291 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
292 ::SetSplineOrder(const unsigned int & splineOrder)
294 SizeType splineOrders;
295 for (unsigned int i=0;i<InputDimension;i++)splineOrders[i]=splineOrder;
297 this->SetSplineOrders(splineOrders);
301 // JV set Spline Orders
302 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
304 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
305 ::SetSplineOrders(const SizeType & splineOrders)
307 if(m_SplineOrders!=splineOrders)
309 m_SplineOrders=splineOrders;
311 //update the interpolation function
312 m_VectorInterpolator->SetSplineOrders(m_SplineOrders);
314 //update the varaibles for computing interpolation
315 for (unsigned int i=0; i <InputDimension;i++)
317 m_SupportSize[i] = m_SplineOrders[i]+1;
318 m_Offset[i] = m_SplineOrders[i] / 2;
320 if ( m_SplineOrders[i] % 2 )
322 m_SplineOrderOdd[i] = true;
326 m_SplineOrderOdd[i] = false;
330 //SupportSize is updated!, update regions
331 //JV initialize some variables for jacobian calculation
332 m_SupportRegion.SetSize(m_SupportSize);
333 m_SupportIndex.Fill(0);
334 m_SupportRegion.SetIndex(m_SupportIndex);
336 // Update the WeightRatios
337 for (unsigned int i=0; i<4; i++)
338 for (unsigned int j=0; j<4; j++)
339 m_WeightRatio[i][j]=m_Weights[m_SplineOrders[InputDimension-1]][i]/ m_Weights[m_SplineOrders[InputDimension-1]][j];
346 // JV set sampling factor
347 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
349 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
350 ::SetLUTSamplingFactor( const int & samplingFactor)
352 SizeType samplingFactors;
353 for (unsigned int i=0; i<NInputDimensions; i++)
354 samplingFactors[i]=samplingFactor;
356 //update the interpolation function
357 this->SetLUTSamplingFactors(samplingFactors);
361 // JV set sampling factors
362 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
364 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
365 ::SetLUTSamplingFactors( const SizeType & samplingFactors)
367 if(m_LUTSamplingFactors!=samplingFactors)
369 for (unsigned int i=0; i<NInputDimensions; i++)
370 m_LUTSamplingFactors[i]=samplingFactors[i];
372 //update the interpolation function
373 m_VectorInterpolator->SetLUTSamplingFactors(m_LUTSamplingFactors);
380 // Get the number of parameters
381 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
383 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
384 ::GetNumberOfParameters(void) const
387 // The number of parameters equal SpaceDimension * number of
388 // of pixels in the grid region.
389 return ( static_cast<unsigned int>( SpaceDimension ) *
390 static_cast<unsigned int>( m_GridRegion.GetNumberOfPixels() ) );
395 // Get the padded number of parameters
396 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
398 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
399 ::GetPaddedNumberOfParameters(void) const
402 // The number of parameters equal SpaceDimension * number of
403 // of pixels in the grid region.
404 return ( static_cast<unsigned int>( SpaceDimension ) *
405 static_cast<unsigned int>( m_PaddedGridRegion.GetNumberOfPixels() ) );
411 // Get the number of parameters per dimension
412 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
414 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
415 ::GetNumberOfParametersPerDimension(void) const
417 // The number of parameters per dimension equal number of
418 // of pixels in the grid region.
419 return ( static_cast<unsigned int>( m_GridRegion.GetNumberOfPixels() ) );
424 // Set the grid region
425 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
427 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
428 ::SetGridRegion( const RegionType & region )
430 if ( m_GridRegion != region )
432 m_GridRegion = region;
433 m_PaddedGridRegion=region;
435 // JV set the padded region
436 typename CoefficientImageType::RegionType::SizeType paddedSize= region.GetSize();
438 //JV size dependes on shape
439 switch (m_TransformShape)
443 paddedSize[InputDimension-1]+=4;
447 paddedSize[InputDimension-1]+=4;
451 paddedSize[InputDimension-1]+=3;
455 paddedSize[InputDimension-1]+=2;
459 paddedSize[InputDimension-1]+=3;
462 paddedSize[InputDimension-1]=+1;
464 m_PaddedGridRegion.SetSize(paddedSize);
466 // Set regions for each coefficient and jacobian image
467 m_WrappedImage->SetRegions( m_GridRegion );
468 m_PaddedCoefficientImage->SetRegions( m_PaddedGridRegion );
469 m_PaddedCoefficientImage->Allocate();
470 for (unsigned int j=0; j <OutputDimension;j++)
472 m_JacobianImage[j]->SetRegions( m_GridRegion );
475 // JV used the padded version for the valid region
476 // Set the valid region
477 // If the grid spans the interval [start, last].
478 // The valid interval for evaluation is [start+offset, last-offset]
479 // when spline order is even.
480 // The valid interval for evaluation is [start+offset, last-offset)
481 // when spline order is odd.
482 // Where offset = vcl_floor(spline / 2 ).
483 // Note that the last pixel is not included in the valid region
484 // with odd spline orders.
485 typename RegionType::SizeType size = m_PaddedGridRegion.GetSize();
486 typename RegionType::IndexType index = m_PaddedGridRegion.GetIndex();
487 for ( unsigned int j = 0; j < NInputDimensions; j++ )
490 static_cast< typename RegionType::IndexValueType >( m_Offset[j] );
492 static_cast< typename RegionType::SizeValueType> ( 2 * m_Offset[j] );
493 m_ValidRegionLast[j] = index[j] +
494 static_cast< typename RegionType::IndexValueType >( size[j] ) - 1;
496 m_ValidRegion.SetSize( size );
497 m_ValidRegion.SetIndex( index );
499 // If we are using the default parameters, update their size and set to identity.
500 // Input parameters point to internal buffer => using default parameters.
501 if (m_InputParametersPointer == &m_InternalParametersBuffer)
503 // Check if we need to resize the default parameter buffer.
504 if ( m_InternalParametersBuffer.GetSize() != this->GetNumberOfParameters() )
506 m_InternalParametersBuffer.SetSize( this->GetNumberOfParameters() );
507 // Fill with zeros for identity.
508 m_InternalParametersBuffer.Fill( 0 );
517 // Set the grid spacing
518 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
520 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
521 ::SetGridSpacing( const SpacingType & spacing )
523 if ( m_GridSpacing != spacing )
525 m_GridSpacing = spacing;
527 // Set spacing for each coefficient and jacobian image
528 m_WrappedImage->SetSpacing( m_GridSpacing.GetDataPointer() );
529 m_PaddedCoefficientImage->SetSpacing( m_GridSpacing.GetDataPointer() );
530 for (unsigned int j=0; j <OutputDimension; j++) m_JacobianImage[j]->SetSpacing( m_GridSpacing.GetDataPointer() );
534 for( unsigned int i=0; i<OutputDimension; i++)
536 scale[i][i] = m_GridSpacing[i];
539 m_IndexToPoint = m_GridDirection * scale;
540 m_PointToIndex = m_IndexToPoint.GetInverse();
547 // Set the grid direction
548 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
550 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
551 ::SetGridDirection( const DirectionType & direction )
553 if ( m_GridDirection != direction )
555 m_GridDirection = direction;
557 // Set direction for each coefficient and jacobian image
558 m_WrappedImage->SetDirection( m_GridDirection );
559 m_PaddedCoefficientImage->SetDirection( m_GridDirection );
560 for (unsigned int j=0; j <OutputDimension; j++) m_JacobianImage[j]->SetDirection( m_GridDirection );
564 for( unsigned int i=0; i<OutputDimension; i++)
566 scale[i][i] = m_GridSpacing[i];
569 m_IndexToPoint = m_GridDirection * scale;
570 m_PointToIndex = m_IndexToPoint.GetInverse();
578 // Set the grid origin
579 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
581 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
582 ::SetGridOrigin( const OriginType& origin )
584 if( m_GridOrigin!=origin)
586 m_GridOrigin = origin;
588 // JV The origin depends on the shape
589 switch (m_TransformShape)
593 m_PaddedGridOrigin=origin;
594 m_PaddedGridOrigin[InputDimension-1]=origin[InputDimension-1]-2* m_GridSpacing[InputDimension-1];
598 m_PaddedGridOrigin=origin;
599 m_PaddedGridOrigin[InputDimension-1]=origin[InputDimension-1]-1* m_GridSpacing[InputDimension-1];
603 m_PaddedGridOrigin=origin;
604 m_PaddedGridOrigin[InputDimension-1]=origin[InputDimension-1]-1* m_GridSpacing[InputDimension-1];
608 m_PaddedGridOrigin=origin;
612 m_PaddedGridOrigin=origin;
613 m_PaddedGridOrigin[InputDimension-1]=origin[InputDimension-1]-1* m_GridSpacing[InputDimension-1];
616 m_PaddedGridOrigin=origin;
619 // Set origin for each coefficient and jacobianimage
620 m_WrappedImage->SetOrigin( m_GridOrigin.GetDataPointer() );
621 m_PaddedCoefficientImage->SetOrigin( m_PaddedGridOrigin.GetDataPointer() );
622 for (unsigned int j=0; j <OutputDimension; j++) m_JacobianImage[j]->SetOrigin( m_GridOrigin.GetDataPointer() );
629 // Set the parameters
630 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
632 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
635 if( m_InputParametersPointer )
637 ParametersType * parameters =
638 const_cast<ParametersType *>( m_InputParametersPointer );
639 parameters->Fill( 0.0 );
644 itkExceptionMacro( << "Input parameters for the spline haven't been set ! "
645 << "Set them using the SetParameters or SetCoefficientImage method first." );
650 // Set the parameters
651 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
653 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
654 ::SetParameters( const ParametersType & parameters )
657 // Check if the number of parameters match the
658 // Expected number of parameters
659 if ( parameters.Size() != this->GetNumberOfParameters() )
661 itkExceptionMacro(<<"Mismatched between parameters size "
663 << " and region size "
664 << m_GridRegion.GetNumberOfPixels() );
667 // Clean up buffered parameters
668 m_InternalParametersBuffer = ParametersType( 0 );
670 // Keep a reference to the input parameters
671 m_InputParametersPointer = ¶meters;
673 // Wrap flat array as images of coefficients
674 this->WrapAsImages();
676 //JV Set padded input to vector interpolator
677 m_VectorInterpolator->SetInputImage(this->GetPaddedCoefficientImage());
679 // Modified is always called since we just have a pointer to the
680 // parameters and cannot know if the parameters have changed.
685 // Set the Fixed Parameters
686 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
688 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
689 ::SetFixedParameters( const ParametersType & parameters )
692 // JV number should be exact, no defaults for spacing
693 if ( parameters.Size() != NInputDimensions * (5 + NInputDimensions)+3 )
695 itkExceptionMacro(<< "Mismatched between parameters size "
697 << " and number of fixed parameters "
698 << NInputDimensions * (5 + NInputDimensions)+3 );
700 /*********************************************************
701 Fixed Parameters store the following information:
706 // JV we add the splineOrders, LUTsamplingfactor, mask pointer and bulktransform pointer
714 The size of these is equal to the NInputDimensions
715 *********************************************************/
717 /** Set the Grid Parameters */
719 for (unsigned int i=0; i<NInputDimensions; i++)
721 gridSize[i] = static_cast<int> (parameters[i]);
723 RegionType bsplineRegion;
724 bsplineRegion.SetSize( gridSize );
726 /** Set the Origin Parameters */
728 for (unsigned int i=0; i<NInputDimensions; i++)
730 origin[i] = parameters[NInputDimensions+i];
733 /** Set the Spacing Parameters */
735 for (unsigned int i=0; i<NInputDimensions; i++)
737 spacing[i] = parameters[2*NInputDimensions+i];
740 /** Set the Direction Parameters */
741 DirectionType direction;
742 for (unsigned int di=0; di<NInputDimensions; di++)
744 for (unsigned int dj=0; dj<NInputDimensions; dj++)
746 direction[di][dj] = parameters[3*NInputDimensions+(di*NInputDimensions+dj)];
750 //JV add the spline orders
751 SizeType splineOrders;
752 for (unsigned int i=0; i<NInputDimensions; i++)
754 splineOrders[i]=parameters[(3+NInputDimensions)*NInputDimensions+i];
757 //JV add the LUT sampling factor
758 SizeType samplingFactors;
759 for (unsigned int i=0; i<NInputDimensions; i++)
761 samplingFactors[i]=parameters[(4+NInputDimensions)*NInputDimensions+i];
764 //JV add the MaskPointer
765 m_Mask=((MaskPointer)( (size_t)parameters[(5+NInputDimensions)*NInputDimensions]));
767 //JV add the MaskPointer
768 m_BulkTransform=((BulkTransformPointer)( (size_t)parameters[(5+NInputDimensions)*NInputDimensions+1]));
770 //JV add the TransformShape
771 m_TransformShape=((unsigned int)parameters[(5+NInputDimensions)*NInputDimensions+2]);
774 this->SetSplineOrders( splineOrders );
775 this->SetGridSpacing( spacing );
776 this->SetGridDirection( direction );
777 this->SetGridOrigin( origin );
778 this->SetGridRegion( bsplineRegion );
779 this->SetLUTSamplingFactors( samplingFactors );
784 // Wrap flat parameters as images
785 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
787 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
790 //JV Wrap parameter array in vectorial image, changed parameter order: A1x A1y A1z, A2x ....
791 PixelType * dataPointer =reinterpret_cast<PixelType *>( const_cast<double *>(m_InputParametersPointer->data_block() )) ;
792 unsigned int numberOfPixels = m_GridRegion.GetNumberOfPixels();
794 m_WrappedImage->GetPixelContainer()->SetImportPointer( dataPointer,numberOfPixels);//InputDimension
795 m_CoefficientImage = m_WrappedImage;
797 //=====================================
798 //JV Create padded structure adding BC
799 //=====================================
800 PadCoefficientImage();
802 //=====================================
803 //JV Wrap jacobian into OutputDimension X Vectorial images
804 //=====================================
805 this->m_Jacobian.set_size( OutputDimension, this->GetNumberOfParameters() );
807 // Use memset to set the memory
808 // JV four rows of three comps of parameters
809 JacobianPixelType * jacobianDataPointer = reinterpret_cast<JacobianPixelType *>(this->m_Jacobian.data_block());
810 memset(jacobianDataPointer, 0, OutputDimension*numberOfPixels*sizeof(JacobianPixelType));
812 for (unsigned int j=0; j<OutputDimension; j++)
814 m_JacobianImage[j]->GetPixelContainer()->
815 SetImportPointer( jacobianDataPointer, numberOfPixels );
816 jacobianDataPointer += numberOfPixels;
819 // Reset the J parameters
820 m_LastJacobianIndex = m_ValidRegion.GetIndex();
821 m_FirstRegion.SetSize(m_NullSize);
822 m_SecondRegion.SetSize(m_NullSize);
823 for ( unsigned int j = 0; j < SpaceDimension ; j++ )
825 m_FirstIterator[j]= IteratorType( m_JacobianImage[j], m_FirstRegion);
826 m_SecondIterator[j]= IteratorType( m_JacobianImage[j], m_SecondRegion);
829 m_BCValues.resize(0);
830 m_BCRegions.resize(0);
832 m_BC2Values.resize(0);
833 m_BC2Regions.resize(0);
835 m_BC3Values.resize(0);
836 m_BC3Regions.resize(0);
842 // Set the parameters by value
843 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
845 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
846 ::SetParametersByValue( const ParametersType & parameters )
849 // Check if the number of parameters match the
850 // Expected number of parameters
851 if ( parameters.Size() != this->GetNumberOfParameters() )
853 itkExceptionMacro(<<"Mismatched between parameters size "
855 << " and region size "
856 << m_GridRegion.GetNumberOfPixels() );
860 m_InternalParametersBuffer = parameters;
861 m_InputParametersPointer = &m_InternalParametersBuffer;
863 // Wrap flat array as images of coefficients
864 this->WrapAsImages();
866 //JV Set padded input to vector interpolator
867 m_VectorInterpolator->SetInputImage(this->GetPaddedCoefficientImage());
869 // Modified is always called since we just have a pointer to the
870 // Parameters and cannot know if the parameters have changed.
875 // Get the parameters
876 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
878 typename ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
880 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
881 ::GetParameters( void ) const
883 /** NOTE: For efficiency, this class does not keep a copy of the parameters -
884 * it just keeps pointer to input parameters.
886 if (NULL == m_InputParametersPointer)
888 itkExceptionMacro( <<"Cannot GetParameters() because m_InputParametersPointer is NULL. Perhaps SetCoefficientImage() has been called causing the NULL pointer." );
891 return (*m_InputParametersPointer);
895 // Get the parameters
896 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
898 typename ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
900 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
901 ::GetFixedParameters( void ) const
903 RegionType resRegion = this->GetGridRegion( );
905 for (unsigned int i=0; i<NInputDimensions; i++)
907 this->m_FixedParameters[i] = (resRegion.GetSize())[i];
909 for (unsigned int i=0; i<NInputDimensions; i++)
911 this->m_FixedParameters[NInputDimensions+i] = (this->GetGridOrigin())[i];
913 for (unsigned int i=0; i<NInputDimensions; i++)
915 this->m_FixedParameters[2*NInputDimensions+i] = (this->GetGridSpacing())[i];
917 for (unsigned int di=0; di<NInputDimensions; di++)
919 for (unsigned int dj=0; dj<NInputDimensions; dj++)
921 this->m_FixedParameters[3*NInputDimensions+(di*NInputDimensions+dj)] = (this->GetGridDirection())[di][dj];
925 //JV add splineOrders
926 for (unsigned int i=0; i<NInputDimensions; i++)
928 this->m_FixedParameters[(3+NInputDimensions)*NInputDimensions+i] = (this->GetSplineOrders())[i];
931 //JV add LUTsamplingFactor
932 for (unsigned int i=0; i<NInputDimensions; i++)
934 this->m_FixedParameters[(4+NInputDimensions)*NInputDimensions+i] = (this->GetLUTSamplingFactors())[i];
938 this->m_FixedParameters[(5+NInputDimensions)*NInputDimensions]=(double)((size_t) m_Mask);
940 //JV add the bulktransform pointer
941 this->m_FixedParameters[(5+NInputDimensions)*NInputDimensions+1]=(double)((size_t) m_BulkTransform);
943 //JV add the transform shape
944 this->m_FixedParameters[(5+NInputDimensions)*NInputDimensions+2]=(double)( m_TransformShape);
946 return (this->m_FixedParameters);
950 // Set the B-Spline coefficients using input images
951 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
953 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
954 ::SetCoefficientImage( CoefficientImagePointer image )
956 this->SetGridSpacing( image->GetSpacing() );
957 this->SetGridOrigin( image->GetOrigin() );
958 this->SetGridDirection( image->GetDirection() );
959 this->SetGridRegion( image->GetBufferedRegion() );
960 m_CoefficientImage = image;
963 // m_WrappedImage=m_CoefficientImage;
965 // Update the interpolator
966 this->PadCoefficientImage();
967 m_VectorInterpolator->SetInputImage(this->GetPaddedCoefficientImage());
969 // Clean up buffered parameters
970 m_InternalParametersBuffer = ParametersType( 0 );
971 m_InputParametersPointer = NULL;
976 // // Set the B-Spline coefficients using input images
977 // template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
979 // ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
980 // ::SetPaddedCoefficientImage( CoefficientImagePointer image )
982 // //JV modify the region
983 // typename CoefficientImageType::RegionType region=image->GetBufferedRegion();
984 // typename CoefficientImageType::RegionType::SizeType size=region.GetSize();
985 // size[InputDimension-1]-=2;
986 // region.SetSize(size);
989 // this->SetGridRegion( region );
990 // this->SetGridSpacing( image->GetSpacing() );
991 // this->SetGridDirection( image->GetDirection() );
992 // this->SetGridOrigin( image->GetOrigin() );
993 // m_PaddedCoefficientImage = image;
994 // this->ExtractCoefficientImage();
995 // m_VectorInterpolator->SetInputImage(this->GetPaddedCoefficientImage());
997 // // Clean up buffered parameters
998 // m_InternalParametersBuffer = ParametersType( 0 );
999 // m_InputParametersPointer = NULL;
1003 // Set the B-Spline coefficients using input images
1004 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
1005 typename ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>::CoefficientImageType::Pointer
1006 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
1007 ::ExtractTemporalRow(const typename CoefficientImageType::Pointer& coefficientImage, unsigned int temporalIndex)
1010 typename CoefficientImageType::RegionType sourceRegion=coefficientImage->GetLargestPossibleRegion();
1011 sourceRegion.SetSize(InputDimension-1, 1);
1012 sourceRegion.SetIndex(InputDimension-1, temporalIndex);
1015 typedef clitk::ExtractImageFilter<CoefficientImageType, CoefficientImageType> ExtractImageFilterType;
1016 typename ExtractImageFilterType::Pointer extract=ExtractImageFilterType::New();
1017 extract->SetInput(coefficientImage);
1018 extract->SetExtractionRegion(sourceRegion);
1020 typename CoefficientImageType::Pointer row= extract->GetOutput();
1022 // Set index to zero
1023 sourceRegion.SetIndex(InputDimension-1, 0);
1024 row->SetRegions(sourceRegion);
1029 // Set the B-Spline coefficients using input images
1030 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
1032 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
1033 ::PadCoefficientImage(void)
1036 // Define paste, extract and combine filters
1037 typedef itk::PasteImageFilter<CoefficientImageType, CoefficientImageType, CoefficientImageType> PasteImageFilterType;
1038 typedef clitk::ExtractImageFilter<CoefficientImageType, CoefficientImageType> ExtractImageFilterType;
1039 typedef clitk::LinearCombinationImageFilter<CoefficientImageType, CoefficientImageType> LinearCombinationFilterType;
1040 typedef itk::MultiplyByConstantImageFilter<CoefficientImageType, double, CoefficientImageType> MultiplicationFilterType;
1043 typename CoefficientImageType::RegionType sourceRegion=m_PaddedCoefficientImage->GetLargestPossibleRegion();
1044 typename CoefficientImageType::RegionType destinationRegion=m_PaddedCoefficientImage->GetLargestPossibleRegion();
1045 typename CoefficientImageType::RegionType::SizeType sourceSize=sourceRegion.GetSize();
1046 typename CoefficientImageType::RegionType::SizeType destinationSize=destinationRegion.GetSize();
1047 typename CoefficientImageType::IndexType sourceIndex=sourceRegion.GetIndex();
1048 typename CoefficientImageType::IndexType destinationIndex=destinationRegion.GetIndex();
1050 // JV Padding depends on the shape
1051 switch (m_TransformShape)
1056 2: rabbit 4 CP 3 DOF
1057 3: rabbit 5 CP 4 DOF
1058 4: sputnik 4 CP 4 DOF
1059 5: sputnik 5 CP 5 DOF
1060 6: diamond 6 CP 5 DOF
1061 7: diamond 7 CP 6 DOF
1066 // ----------------------------------------------------------------------
1067 // The egg with 4 internal CP (starting from inhale)
1068 // Periodic, constrained to zero at the reference
1069 // at position 3 and
1070 // Coeff row BC1 BC2 0 BC3 1 2 BC4
1071 // PaddedCoeff R: 0 1 2 3 4 5 6
1074 // BC3= -weights[2]/weights[0] ( R2+R4 )
1076 // ---------------------------------------------------------------------
1078 //---------------------------------
1079 // 1. First two temporal rows are identical: paste 1-2 to 0-1
1080 typename PasteImageFilterType::Pointer paster0=PasteImageFilterType::New();
1081 paster0->SetDestinationImage(m_PaddedCoefficientImage);
1082 paster0->SetSourceImage(m_CoefficientImage);
1083 sourceRegion.SetIndex(NInputDimensions-1,1);
1084 sourceRegion.SetSize(NInputDimensions-1,2);
1085 destinationIndex[InputDimension-1]=0;
1086 paster0->SetDestinationIndex(destinationIndex);
1087 paster0->SetSourceRegion(sourceRegion);
1089 //---------------------------------
1090 // 1. Next temporal row = paste 0 to 2
1091 typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1092 typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1093 paster1->SetDestinationImage(paster0->GetOutput());
1094 paster1->SetSourceImage(row0);
1095 destinationIndex[InputDimension-1]=2;
1096 paster1->SetDestinationIndex(destinationIndex);
1097 paster1->SetSourceRegion(row0->GetLargestPossibleRegion());
1099 //---------------------------------
1100 // 2. Middle row at index 3=BC
1101 // Extract row at index 1, 2 (of coeff image)
1102 typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1105 typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1106 combine1->SetInput(0,row0);
1107 combine1->SetInput(1,row1);
1108 combine1->SetA(-m_WeightRatio[2][0]);
1109 combine1->SetB(-m_WeightRatio[2][0]);
1111 typename CoefficientImageType::Pointer bc3Row=combine1->GetOutput();
1113 // Paste middleRow at index 3 (padded coeff)
1114 typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1115 paster2->SetDestinationImage(paster1->GetOutput());
1116 paster2->SetSourceImage(bc3Row);
1117 destinationIndex[InputDimension-1]=3;
1118 paster2->SetDestinationIndex(destinationIndex);
1119 paster2->SetSourceRegion(bc3Row->GetLargestPossibleRegion());
1121 //---------------------------------
1122 // 3. Next two temporal rows identical: paste 1,2 to 4,5
1123 typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1124 paster3->SetDestinationImage(paster2->GetOutput());
1125 paster3->SetSourceImage(m_CoefficientImage);
1126 sourceRegion.SetIndex(InputDimension-1,1);
1127 sourceRegion.SetSize(InputDimension-1,2);
1128 destinationIndex[InputDimension-1]=4;
1129 paster3->SetDestinationIndex(destinationIndex);
1130 paster3->SetSourceRegion(sourceRegion);
1132 // ---------------------------------
1133 // 4. Row at index 6=BC4= R2
1134 // Paste BC3 row at index 5
1135 typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
1136 paster4->SetDestinationImage(paster3->GetOutput());
1137 paster4->SetSourceImage(row0);
1138 destinationIndex[InputDimension-1]=6;
1139 paster4->SetDestinationIndex(destinationIndex);
1140 paster4->SetSourceRegion(row0->GetLargestPossibleRegion());
1142 // Update the chain!
1144 m_PaddedCoefficientImage= paster4->GetOutput();
1151 // ----------------------------------------------------------------------
1152 // The egg with 5 internal CP (starting from inhale)
1153 // Periodic, constrained to zero at the reference
1154 // at position 3 and
1155 // Coeff row BC1 BC2 0 BC3 1 2 3 BC4
1156 // PaddedCoeff R: 0 1 2 3 4 5 6 7
1159 // BC3= -weights[3]/weights[1] ( R2+R5 ) - R4
1161 // ---------------------------------------------------------------------
1162 //---------------------------------
1163 // 1. First two temporal rows are identical: paste 2-3 to 0-1
1164 typename PasteImageFilterType::Pointer paster0=PasteImageFilterType::New();
1165 paster0->SetDestinationImage(m_PaddedCoefficientImage);
1166 paster0->SetSourceImage(m_CoefficientImage);
1167 sourceRegion.SetIndex(NInputDimensions-1,2);
1168 sourceRegion.SetSize(NInputDimensions-1,2);
1169 destinationIndex[InputDimension-1]=0;
1170 paster0->SetDestinationIndex(destinationIndex);
1171 paster0->SetSourceRegion(sourceRegion);
1173 //---------------------------------
1174 // 1. Next temporal row = paste 0 to 2
1175 typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1176 typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1177 paster1->SetDestinationImage(paster0->GetOutput());
1178 paster1->SetSourceImage(row0);
1179 destinationIndex[InputDimension-1]=2;
1180 paster1->SetDestinationIndex(destinationIndex);
1181 paster1->SetSourceRegion(row0->GetLargestPossibleRegion());
1183 //---------------------------------
1184 // 2. Middle row at index 3=BC
1185 // Extract row at index 1, 2 (of coeff image)
1186 typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1187 typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
1190 typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1191 combine1->SetInput(0,row0);
1192 combine1->SetInput(1,row2);
1195 // combine1->Update();
1196 typename LinearCombinationFilterType::Pointer combine2=LinearCombinationFilterType::New();
1197 combine2->SetInput(0,row1);
1198 combine2->SetInput(1,combine1->GetOutput());
1199 combine2->SetA(-1.);
1200 combine2->SetB(-m_WeightRatio[3][1]);
1202 typename CoefficientImageType::Pointer bc3Row=combine2->GetOutput();
1204 // Paste middleRow at index 3 (padded coeff)
1205 typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1206 paster2->SetDestinationImage(paster1->GetOutput());
1207 paster2->SetSourceImage(bc3Row);
1208 destinationIndex[InputDimension-1]=3;
1209 paster2->SetDestinationIndex(destinationIndex);
1210 paster2->SetSourceRegion(bc3Row->GetLargestPossibleRegion());
1212 //---------------------------------
1213 // 3. Next three temporal rows identical: paste 1,2,3 to 4,5,6
1214 typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1215 paster3->SetDestinationImage(paster2->GetOutput());
1216 paster3->SetSourceImage(m_CoefficientImage);
1217 sourceRegion.SetIndex(InputDimension-1,1);
1218 sourceRegion.SetSize(InputDimension-1,3);
1219 destinationIndex[InputDimension-1]=4;
1220 paster3->SetDestinationIndex(destinationIndex);
1221 paster3->SetSourceRegion(sourceRegion);
1223 // ---------------------------------
1224 // 4. Row at index 7=BC4= R2
1225 // Paste BC3 row at index 5
1226 typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
1227 paster4->SetDestinationImage(paster3->GetOutput());
1228 paster4->SetSourceImage(row0);
1229 destinationIndex[InputDimension-1]=7;
1230 paster4->SetDestinationIndex(destinationIndex);
1231 paster4->SetSourceRegion(row0->GetLargestPossibleRegion());
1233 // Update the chain!
1235 m_PaddedCoefficientImage= paster4->GetOutput();
1241 // // ----------------------------------------------------------------------
1242 // // The egg with 5 internal CP:
1243 // // Periodic, C2 smooth everywhere and constrained to zero at the reference
1244 // // Coeff row R5 BC1 0 1 2 3 BC2 R2
1245 // // PaddedCoeff R: 0 1 2 3 4 5 6 7
1246 // // BC1= -weights[2]/weights[0] ( R2+R5)
1248 // // ---------------------------------------------------------------------
1250 // // Extract rows with index 0 and 3
1251 // typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1252 // typename CoefficientImageType::Pointer row3=this->ExtractTemporalRow(m_CoefficientImage, 3);
1254 // // Paste the first row
1255 // typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1256 // destinationIndex.Fill(0);
1257 // paster1->SetDestinationIndex(destinationIndex);
1258 // paster1->SetSourceRegion(row3->GetLargestPossibleRegion());
1259 // paster1->SetSourceImage(row3);
1260 // paster1->SetDestinationImage(m_PaddedCoefficientImage);
1262 // // Linearly Combine rows for BC1 and BC2
1263 // typedef clitk::LinearCombinationImageFilter<CoefficientImageType, CoefficientImageType> LinearCombinationFilterType;
1264 // typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1265 // combine1->SetFirstInput(row0);
1266 // combine1->SetSecondInput(row3);
1267 // combine1->SetA(-m_WeightRatio[2][0]);
1268 // combine1->SetB(-m_WeightRatio[2][0]);
1269 // combine1->Update();
1270 // typename CoefficientImageType::Pointer bcRow=combine1->GetOutput();
1272 // // Paste the second row
1273 // typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1274 // destinationIndex[InputDimension-1]=1;
1275 // paster2->SetDestinationIndex(destinationIndex);
1276 // paster2->SetSourceRegion(bcRow->GetLargestPossibleRegion());
1277 // paster2->SetSourceImage(bcRow);
1278 // paster2->SetDestinationImage(paster1->GetOutput());
1280 // // Paste the coefficientImage
1281 // typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1282 // destinationIndex[InputDimension-1]=2;
1283 // paster3->SetDestinationIndex(destinationIndex);
1284 // paster3->SetSourceRegion(m_CoefficientImage->GetLargestPossibleRegion());
1285 // paster3->SetSourceImage(m_CoefficientImage);
1286 // paster3->SetDestinationImage(paster2->GetOutput());
1288 // // Paste the last two rows
1289 // typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
1290 // destinationIndex[InputDimension-1]=6;
1291 // paster4->SetDestinationIndex(destinationIndex);
1292 // paster4->SetSourceRegion(bcRow->GetLargestPossibleRegion());
1293 // paster4->SetSourceImage(bcRow);
1294 // paster4->SetDestinationImage(paster3->GetOutput());
1296 // typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
1297 // destinationIndex[InputDimension-1]=7;
1298 // paster5->SetDestinationIndex(destinationIndex);
1299 // paster5->SetSourceRegion(row0->GetLargestPossibleRegion());
1300 // paster5->SetSourceImage(row0);
1301 // paster5->SetDestinationImage(paster4->GetOutput());
1303 // // Update the chain!
1304 // paster5->Update();
1305 // m_PaddedCoefficientImage= paster5->GetOutput();
1312 // ----------------------------------------------------------------------
1313 // The rabbit with 4 internal CP
1314 // Periodic, constrained to zero at the reference
1315 // at position 3 and the extremes fixed with anit-symm bc
1316 // Coeff row BC1 0 1 BC2 2 BC3 BC4
1317 // PaddedCoeff R: 0 1 2 3 4 5 6
1319 // BC2= -weights[2]/weights[0] ( R2+R4 )
1322 // ---------------------------------------------------------------------
1324 // ---------------------------------
1325 // 0. First Row =BC1
1326 typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1327 typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1328 typename LinearCombinationFilterType::Pointer combine0=LinearCombinationFilterType::New();
1329 combine0->SetInput(0,row0);
1330 combine0->SetInput(1,row1);
1332 combine0->SetB(-1.);
1334 typename CoefficientImageType::Pointer bc1Row=combine0->GetOutput();
1336 typename PasteImageFilterType::Pointer paster0=PasteImageFilterType::New();
1337 paster0->SetDestinationImage(m_PaddedCoefficientImage);
1338 paster0->SetSourceImage(bc1Row);
1339 destinationIndex[InputDimension-1]=0;
1340 paster0->SetDestinationIndex(destinationIndex);
1341 paster0->SetSourceRegion(bc1Row->GetLargestPossibleRegion());
1343 //---------------------------------
1344 // 1. Next two temporal rows are identical: paste 0-1 to 1-2
1345 typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1346 paster1->SetDestinationImage(paster0->GetOutput());
1347 paster1->SetSourceImage(m_CoefficientImage);
1348 sourceRegion.SetIndex(NInputDimensions-1,0);
1349 destinationIndex[InputDimension-1]=1;
1350 sourceRegion.SetSize(NInputDimensions-1,2);
1351 paster1->SetDestinationIndex(destinationIndex);
1352 paster1->SetSourceRegion(sourceRegion);
1354 //---------------------------------
1355 // 2. Middle row at index 3=BC
1356 // Extract row at index 1, 2 (of coeff image)
1357 typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
1360 typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1361 combine1->SetInput(0,row1);
1362 combine1->SetInput(1,row2);
1363 combine1->SetA(-m_WeightRatio[2][0]);
1364 combine1->SetB(-m_WeightRatio[2][0]);
1366 typename CoefficientImageType::Pointer bc2Row=combine1->GetOutput();
1368 // Paste middleRow at index 3 (padded coeff)
1369 typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1370 paster2->SetDestinationImage(paster1->GetOutput());
1371 paster2->SetSourceImage(bc2Row);
1372 destinationIndex[InputDimension-1]=3;
1373 paster2->SetDestinationIndex(destinationIndex);
1374 paster2->SetSourceRegion(bc2Row->GetLargestPossibleRegion());
1376 //---------------------------------
1377 // 3. Next temporal row is identical: paste 2 to 4
1378 typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1379 paster3->SetDestinationImage(paster2->GetOutput());
1380 paster3->SetSourceImage(row2);
1381 destinationIndex[InputDimension-1]=4;
1382 paster3->SetDestinationIndex(destinationIndex);
1383 paster3->SetSourceRegion(row2->GetLargestPossibleRegion());
1385 // ---------------------------------
1386 // 4. Row at index 5=BC (paddedcoeff image) R1
1387 // Paste BC3 row at index 5
1388 typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
1389 paster4->SetDestinationImage(paster3->GetOutput());
1390 paster4->SetSourceImage(row0);
1391 destinationIndex[InputDimension-1]=5;
1392 paster4->SetDestinationIndex(destinationIndex);
1393 paster4->SetSourceRegion(row0->GetLargestPossibleRegion());
1395 // ---------------------------------
1396 // 5. Paste BC4 row at index 6
1397 typename LinearCombinationFilterType::Pointer combine3=LinearCombinationFilterType::New();
1398 combine3->SetInput(0,row0);
1399 combine3->SetInput(1,row2);
1401 combine3->SetB(-1.);
1403 typename CoefficientImageType::Pointer bc4Row=combine3->GetOutput();
1404 typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
1405 paster5->SetDestinationImage(paster4->GetOutput());
1406 paster5->SetSourceImage(bc4Row);
1407 destinationIndex[InputDimension-1]=6;
1408 paster5->SetDestinationIndex(destinationIndex);
1409 paster5->SetSourceRegion(bc4Row->GetLargestPossibleRegion());
1411 // Update the chain!
1413 m_PaddedCoefficientImage= paster5->GetOutput();
1420 // ----------------------------------------------------------------------
1421 // The rabbit with 5 internal CP
1422 // Periodic, constrained to zero at the reference at position 3.5
1423 // and the extremes fixed with anti-symmetrical BC
1424 // Coeff row BC1 0 1 BC2 2 3 BC3 BC4
1425 // PaddedCoeff R: 0 1 2 3 4 5 6 7
1427 // BC2= -weights[3]/weights[1] ( R2+R5 ) - R4
1430 // ---------------------------------------------------------------------
1432 // ---------------------------------
1433 // 0. First Row =BC1
1434 typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1435 typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1436 typename LinearCombinationFilterType::Pointer combine0=LinearCombinationFilterType::New();
1437 combine0->SetInput(0,row0);
1438 combine0->SetInput(1,row1);
1440 combine0->SetB(-1.);
1442 typename CoefficientImageType::Pointer bc1Row=combine0->GetOutput();
1444 typename PasteImageFilterType::Pointer paster0=PasteImageFilterType::New();
1445 paster0->SetDestinationImage(m_PaddedCoefficientImage);
1446 paster0->SetSourceImage(bc1Row);
1447 destinationIndex[InputDimension-1]=0;
1448 paster0->SetDestinationIndex(destinationIndex);
1449 paster0->SetSourceRegion(bc1Row->GetLargestPossibleRegion());
1451 //---------------------------------
1452 // 1. Next two temporal rows are identical: paste 0-1 to 1-2
1453 typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1454 paster1->SetDestinationImage(paster0->GetOutput());
1455 paster1->SetSourceImage(m_CoefficientImage);
1456 sourceRegion.SetIndex(InputDimension-1,0);
1457 destinationIndex[InputDimension-1]=1;
1458 sourceRegion.SetSize(NInputDimensions-1,2);
1459 paster1->SetDestinationIndex(destinationIndex);
1460 paster1->SetSourceRegion(sourceRegion);
1462 //---------------------------------
1463 // 2. Middle row at index 3=BC
1464 // Extract row at index 1, 3 (of coeff image)
1465 //typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1466 typename CoefficientImageType::Pointer row3=this->ExtractTemporalRow(m_CoefficientImage, 3);
1469 typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1470 combine1->SetInput(0,row1);
1471 combine1->SetInput(1,row3);
1476 // Extract row at index 2 (coeff image)
1477 typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
1480 typename LinearCombinationFilterType::Pointer combine2=LinearCombinationFilterType::New();
1481 combine2->SetInput(0,combine1->GetOutput());
1482 combine2->SetInput(1,row2);
1483 combine2->SetA(-m_WeightRatio[3][1]);
1484 combine2->SetB(-1.);
1486 typename CoefficientImageType::Pointer bc2Row=combine2->GetOutput();
1488 // Paste middleRow at index 3 (padded coeff)
1489 typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1490 paster2->SetDestinationImage(paster1->GetOutput());
1491 paster2->SetSourceImage(bc2Row);
1492 destinationIndex[InputDimension-1]=3;
1493 paster2->SetDestinationIndex(destinationIndex);
1494 paster2->SetSourceRegion(bc2Row->GetLargestPossibleRegion());
1496 //---------------------------------
1497 // 3. Next two temporal rows are identical: paste 2,3 to 4,5
1498 typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1499 paster3->SetDestinationImage(paster2->GetOutput());
1500 paster3->SetSourceImage(m_CoefficientImage);
1501 sourceRegion.SetIndex(InputDimension-1,2);
1502 destinationIndex[InputDimension-1]=4;
1503 sourceRegion.SetSize(NInputDimensions-1,2);
1504 paster3->SetDestinationIndex(destinationIndex);
1505 paster3->SetSourceRegion(sourceRegion);
1507 // ---------------------------------
1508 // 4. Row at index 6=BC (paddedcoeff image)R1
1509 // Paste BC3 row at index 6
1510 typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
1511 paster4->SetDestinationImage(paster3->GetOutput());
1512 paster4->SetSourceImage(row0);
1513 destinationIndex[InputDimension-1]=6;
1514 paster4->SetDestinationIndex(destinationIndex);
1515 paster4->SetSourceRegion(row0->GetLargestPossibleRegion());
1517 // ---------------------------------
1518 // 5. Paste BC4 row at index 7
1519 typename LinearCombinationFilterType::Pointer combine3=LinearCombinationFilterType::New();
1520 combine3->SetInput(0,row0);
1521 combine3->SetInput(1,row3);
1523 combine3->SetB(-1.);
1525 typename CoefficientImageType::Pointer bc4Row=combine3->GetOutput();
1526 typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
1527 paster5->SetDestinationImage(paster4->GetOutput());
1528 paster5->SetSourceImage(bc4Row);
1529 destinationIndex[InputDimension-1]=7;
1530 paster5->SetDestinationIndex(destinationIndex);
1531 paster5->SetSourceRegion(bc4Row->GetLargestPossibleRegion());
1533 // Update the chain!
1535 m_PaddedCoefficientImage= paster5->GetOutput();
1542 // ----------------------------------------------------------------------
1543 // The sputnik with 4 internal CP
1544 // Periodic, constrained to zero at the reference
1545 // at position 3 and one indepent extremes copied
1546 // Coeff row BC1 0 1 BC2 2 BC2 3
1547 // PaddedCoeff R: 0 1 2 3 4 5 6
1549 // BC2= -weights[2]/weights[0] ( R2+R4 )
1550 // BC3= weights[2]/weights[0] ( R2-R4) + R1
1551 // ---------------------------------------------------------------------
1553 //---------------------------------
1554 // 1. First Row is equal to last row: paste 3 row to 0
1555 typename CoefficientImageType::Pointer row3=this->ExtractTemporalRow(m_CoefficientImage, 3);
1556 typename PasteImageFilterType::Pointer paster0=PasteImageFilterType::New();
1557 paster0->SetDestinationImage(m_PaddedCoefficientImage);
1558 paster0->SetSourceImage(row3);
1559 destinationIndex[InputDimension-1]=0;
1560 paster0->SetDestinationIndex(destinationIndex);
1561 paster0->SetSourceRegion(row3->GetLargestPossibleRegion());
1563 //---------------------------------
1564 // 1. Next two temporal rows are identical: paste 0-1 to 1-2
1565 typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1566 paster1->SetDestinationImage(paster0->GetOutput());
1567 paster1->SetSourceImage(m_CoefficientImage);
1568 sourceRegion.SetIndex(NInputDimensions-1,0);
1569 destinationIndex[InputDimension-1]=1;
1570 sourceRegion.SetSize(NInputDimensions-1,2);
1571 paster1->SetDestinationIndex(destinationIndex);
1572 paster1->SetSourceRegion(sourceRegion);
1574 //---------------------------------
1575 // 2. Middle row at index 3=BC
1576 // Extract row at index 1, 2 (of coeff image)
1577 typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1578 typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
1581 typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1582 combine1->SetInput(0,row1);
1583 combine1->SetInput(1,row2);
1584 combine1->SetA(-m_WeightRatio[2][0]);
1585 combine1->SetB(-m_WeightRatio[2][0]);
1587 typename CoefficientImageType::Pointer bc1Row=combine1->GetOutput();
1589 // Paste middleRow at index 3 (padded coeff)
1590 typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1591 paster2->SetDestinationImage(paster1->GetOutput());
1592 paster2->SetSourceImage(bc1Row);
1593 destinationIndex[InputDimension-1]=3;
1594 paster2->SetDestinationIndex(destinationIndex);
1595 paster2->SetSourceRegion(bc1Row->GetLargestPossibleRegion());
1597 //---------------------------------
1598 // 3. Next temporal row is identical: paste 2 to 4
1599 typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1600 paster3->SetDestinationImage(paster2->GetOutput());
1601 paster3->SetSourceImage(row2);
1602 destinationIndex[InputDimension-1]=4;
1603 paster3->SetDestinationIndex(destinationIndex);
1604 paster3->SetSourceRegion(row2->GetLargestPossibleRegion());
1606 //---------------------------------
1607 // 4. Final row at index 5=BC3 (paddedcoeff image)R1 R2 R4 corresponds to index in coeff image 0 1 2
1608 // Extract row at index 1, 2 (of coeff image)already done
1609 typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1612 typename LinearCombinationFilterType::Pointer combine3=LinearCombinationFilterType::New();
1613 combine3->SetInput(0,row1);
1614 combine3->SetInput(1,row2);
1616 combine3->SetB(-1.);
1619 typename LinearCombinationFilterType::Pointer combine4=LinearCombinationFilterType::New();
1620 combine4->SetInput(0,combine3->GetOutput());
1621 combine4->SetInput(1,row0);
1622 combine4->SetA(m_WeightRatio[2][0]);
1625 typename CoefficientImageType::Pointer bc2Row=combine4->GetOutput();
1627 // Paste BC row at index 5
1628 typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
1629 paster4->SetDestinationImage(paster3->GetOutput());
1630 paster4->SetSourceImage(bc2Row);
1631 destinationIndex[InputDimension-1]=5;
1632 paster4->SetDestinationIndex(destinationIndex);
1633 paster4->SetSourceRegion(bc2Row->GetLargestPossibleRegion());
1635 // Paste row 3 at index 6
1636 typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
1637 paster5->SetDestinationImage(paster4->GetOutput());
1638 paster5->SetSourceImage(row3);
1639 destinationIndex[InputDimension-1]=6;
1640 paster5->SetDestinationIndex(destinationIndex);
1641 paster5->SetSourceRegion(row3->GetLargestPossibleRegion());
1643 // Update the chain!
1645 m_PaddedCoefficientImage= paster5->GetOutput();
1653 // ----------------------------------------------------------------------
1654 // The sputnik with 5 internal CP
1655 // Periodic, constrained to zero at the reference
1656 // at position 3 and one indepent extreme
1657 // Coeff row BC1 0 1 BC2 2 3 BC3 4
1658 // PaddedCoeff R: 0 1 2 3 4 5 6 7
1659 // BC1= R2 + R5 - R7
1660 // BC2= -weights[3]/weights[1] ( R2+R5 ) - R4
1661 // BC3= R1 + 0.5 R2 - 0.5 R7
1662 // ----------------------------------------------------------------------
1663 //---------------------------------
1665 typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1666 typename CoefficientImageType::Pointer row3=this->ExtractTemporalRow(m_CoefficientImage, 3);
1667 typename CoefficientImageType::Pointer row4=this->ExtractTemporalRow(m_CoefficientImage, 4);
1670 typename LinearCombinationFilterType::Pointer combine0=LinearCombinationFilterType::New();
1671 combine0->SetInput(0,row1);
1672 combine0->SetInput(1,row3);
1675 //combine0->Update();
1676 typename LinearCombinationFilterType::Pointer combine0bis=LinearCombinationFilterType::New();
1677 combine0bis->SetInput(0,combine0->GetOutput());
1678 combine0bis->SetInput(1,row4);
1679 combine0bis->SetA(1.);
1680 combine0bis->SetB(-1.);
1681 combine0bis->Update();
1682 typename CoefficientImageType::Pointer bc1Row=combine0bis->GetOutput();
1684 typename PasteImageFilterType::Pointer paster0=PasteImageFilterType::New();
1685 paster0->SetDestinationImage(m_PaddedCoefficientImage);
1686 paster0->SetSourceImage(bc1Row);
1687 destinationIndex[InputDimension-1]=0;
1688 paster0->SetDestinationIndex(destinationIndex);
1689 paster0->SetSourceRegion(bc1Row->GetLargestPossibleRegion());
1691 //---------------------------------
1692 // 1. Next two temporal rows are identical: paste 0-1 to 1-2
1693 typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1694 paster1->SetDestinationImage(paster0->GetOutput());
1695 paster1->SetSourceImage(m_CoefficientImage);
1696 sourceRegion.SetIndex(NInputDimensions-1,0);
1697 destinationIndex[InputDimension-1]=1;
1698 sourceRegion.SetSize(NInputDimensions-1,2);
1699 paster1->SetDestinationIndex(destinationIndex);
1700 paster1->SetSourceRegion(sourceRegion);
1702 //---------------------------------
1703 // 2. Middle row at index 3=BC
1704 // Extract row at index 1, 2 (of coeff image)
1705 //typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1706 typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
1709 typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1710 combine1->SetInput(0,row1);
1711 combine1->SetInput(1,row3);
1716 typename LinearCombinationFilterType::Pointer combine2=LinearCombinationFilterType::New();
1717 combine2->SetInput(0,combine1->GetOutput());
1718 combine2->SetInput(1,row2);
1719 combine2->SetA(-m_WeightRatio[3][1]);
1720 combine2->SetB(-1.);
1722 typename CoefficientImageType::Pointer bc2Row=combine2->GetOutput();
1724 // Paste middleRow at index 3 (padded coeff)
1725 typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1726 paster2->SetDestinationImage(paster1->GetOutput());
1727 paster2->SetSourceImage(bc2Row);
1728 destinationIndex[InputDimension-1]=3;
1729 paster2->SetDestinationIndex(destinationIndex);
1730 paster2->SetSourceRegion(bc2Row->GetLargestPossibleRegion());
1732 //---------------------------------
1733 // 3. Next two temporal rows are identical: paste 2,3 to 4,5
1734 typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1735 paster3->SetDestinationImage(paster2->GetOutput());
1736 paster3->SetSourceImage(m_CoefficientImage);
1737 sourceRegion.SetIndex(InputDimension-1,2);
1738 destinationIndex[InputDimension-1]=4;
1739 sourceRegion.SetSize(NInputDimensions-1,2);
1740 paster3->SetDestinationIndex(destinationIndex);
1741 paster3->SetSourceRegion(sourceRegion);
1743 //---------------------------------
1744 // 4. Final row at index 6=BC3
1745 typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1748 typename LinearCombinationFilterType::Pointer combine3=LinearCombinationFilterType::New();
1749 combine3->SetInput(0,row1);
1750 combine3->SetInput(1,row4);
1752 combine3->SetB(-1.);
1753 typename LinearCombinationFilterType::Pointer combine4=LinearCombinationFilterType::New();
1754 combine4->SetInput(0,row0);
1755 combine4->SetInput(1,combine3->GetOutput());
1759 typename CoefficientImageType::Pointer bc3Row=combine4->GetOutput();
1761 // Paste BC row at index 6
1762 typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
1763 paster4->SetDestinationImage(paster3->GetOutput());
1764 paster4->SetSourceImage(bc3Row);
1765 destinationIndex[InputDimension-1]=6;
1766 paster4->SetDestinationIndex(destinationIndex);
1767 paster4->SetSourceRegion(bc3Row->GetLargestPossibleRegion());
1769 // Paste row 4 at index 7
1770 typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
1771 paster5->SetDestinationImage(paster4->GetOutput());
1772 paster5->SetSourceImage(row4);
1773 destinationIndex[InputDimension-1]=7;
1774 paster5->SetDestinationIndex(destinationIndex);
1775 paster5->SetSourceRegion(row4->GetLargestPossibleRegion());
1777 // Update the chain!
1779 m_PaddedCoefficientImage= paster5->GetOutput();
1786 // ----------------------------------------------------------------------
1787 // The diamond with 4 internal CP:
1788 // Periodic, constrained to zero at the reference at position 3
1789 // Coeff row 0 1 2 BC1 3 BC2 4
1790 // PaddedCoeff R:0 1 2 3 4 5 6
1791 // BC1= -weights[2]/weights[0] ( R2+R4 )
1792 // BC2= weights[2]/weights[0] ( R0+R2-R4-R6 ) + R1
1793 // ---------------------------------------------------------------------
1795 //---------------------------------
1796 // 1. First Three temporal rows are identical: paste 0-2 to 0-2
1797 typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1798 paster1->SetDestinationImage(m_PaddedCoefficientImage);
1799 paster1->SetSourceImage(m_CoefficientImage);
1800 sourceIndex.Fill(0);
1801 destinationIndex.Fill(0);
1802 sourceSize[NInputDimensions-1]=3;
1803 sourceRegion.SetSize(sourceSize);
1804 sourceRegion.SetIndex(sourceIndex);
1805 paster1->SetDestinationIndex(destinationIndex);
1806 paster1->SetSourceRegion(sourceRegion);
1808 //---------------------------------
1809 // 2. Middle row at index 3=BC
1810 // Extract row at index 0, 4 (of coeff image)
1811 typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
1812 typename CoefficientImageType::Pointer row3=this->ExtractTemporalRow(m_CoefficientImage, 3);
1815 typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1816 combine1->SetInput(0,row2);
1817 combine1->SetInput(1,row3);
1818 combine1->SetA(-m_WeightRatio[2][0]);
1819 combine1->SetB(-m_WeightRatio[2][0]);
1821 typename CoefficientImageType::Pointer bc1Row=combine1->GetOutput();
1823 // Paste middleRow at index 3 (padded coeff)
1824 typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1825 paster2->SetDestinationImage(paster1->GetOutput());
1826 paster2->SetSourceImage(bc1Row);
1827 destinationIndex[InputDimension-1]=3;
1828 paster2->SetDestinationIndex(destinationIndex);
1829 paster2->SetSourceRegion(bc1Row->GetLargestPossibleRegion());
1831 //---------------------------------
1832 // 3. Next row identical: paste 3 to 4
1833 typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1834 paster3->SetDestinationImage(paster2->GetOutput());
1835 paster3->SetSourceImage(row3);
1836 destinationIndex[InputDimension-1]=4;
1837 paster3->SetDestinationIndex(destinationIndex);
1838 paster3->SetSourceRegion(row3->GetLargestPossibleRegion());
1840 //---------------------------------
1841 // 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
1842 // Extract row at index 0, 2 (of coeff image)already done
1843 typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1844 typename CoefficientImageType::Pointer row4=this->ExtractTemporalRow(m_CoefficientImage, 4);
1845 typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1848 typename LinearCombinationFilterType::Pointer combine3=LinearCombinationFilterType::New();
1849 combine3->SetInput(0,row0);
1850 combine3->SetInput(1,row2);
1855 typename LinearCombinationFilterType::Pointer combine4=LinearCombinationFilterType::New();
1856 combine4->SetInput(0,row3);
1857 combine4->SetInput(1,row4);
1862 typename LinearCombinationFilterType::Pointer combine5=LinearCombinationFilterType::New();
1863 combine5->SetInput(0,combine3->GetOutput());
1864 combine5->SetInput(1,combine4->GetOutput());
1866 combine5->SetB(-1.);
1869 typename LinearCombinationFilterType::Pointer combine6=LinearCombinationFilterType::New();
1870 combine6->SetInput(0,combine5->GetOutput());
1871 combine6->SetInput(1,row1);
1872 combine6->SetA(m_WeightRatio[2][0]);
1875 typename CoefficientImageType::Pointer bc2Row=combine6->GetOutput();
1877 // Paste BC row at index 5
1878 typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
1879 paster4->SetDestinationImage(paster3->GetOutput());
1880 paster4->SetSourceImage(bc2Row);
1881 destinationIndex[InputDimension-1]=5;
1882 paster4->SetDestinationIndex(destinationIndex);
1883 paster4->SetSourceRegion(bc2Row->GetLargestPossibleRegion());
1885 // Paste last row at index 6
1886 typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
1887 paster5->SetDestinationImage(paster4->GetOutput());
1888 paster5->SetSourceImage(row4);
1889 destinationIndex[InputDimension-1]=6;
1890 paster5->SetDestinationIndex(destinationIndex);
1891 paster5->SetSourceRegion(row4->GetLargestPossibleRegion());
1893 // Update the chain!
1895 m_PaddedCoefficientImage= paster5->GetOutput();
1901 // ----------------------------------------------------------------------
1902 // The diamond with 5 internal CP:
1903 // periodic, constrained to zero at the reference at position 3.5
1904 // Coeff row 0 1 2 BC1 3 4 BC2 5
1905 // PaddedCoeff R:0 1 2 3 4 5 6 7
1906 // BC1= -weights[3]/weights[1] ( R2+R5 ) - R4
1907 // BC2= weights[2]/weights[0] ( R0+R2-R5-R7 ) + R1
1908 // ---------------------------------------------------------------------
1910 //---------------------------------
1911 // 1. First Three temporal rows are identical: paste 0-2 to 0-2
1912 typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1913 paster1->SetDestinationImage(m_PaddedCoefficientImage);
1914 paster1->SetSourceImage(m_CoefficientImage);
1915 sourceIndex.Fill(0);
1916 destinationIndex.Fill(0);
1917 sourceSize[NInputDimensions-1]=3;
1918 sourceRegion.SetSize(sourceSize);
1919 sourceRegion.SetIndex(sourceIndex);
1920 paster1->SetDestinationIndex(destinationIndex);
1921 paster1->SetSourceRegion(sourceRegion);
1923 //---------------------------------
1924 // 2. Middle row at index 3=BC
1925 // Extract row at index 0, 4 (of coeff image)
1926 typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
1927 typename CoefficientImageType::Pointer row4=this->ExtractTemporalRow(m_CoefficientImage, 4);
1930 typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1931 combine1->SetInput(0,row2);
1932 combine1->SetInput(1,row4);
1937 // Extract row at index 3 (coeff image)
1938 typename CoefficientImageType::Pointer row3=this->ExtractTemporalRow(m_CoefficientImage, 3);
1941 typename LinearCombinationFilterType::Pointer combine2=LinearCombinationFilterType::New();
1942 combine2->SetInput(0,combine1->GetOutput());
1943 combine2->SetInput(1,row3);
1944 combine2->SetA(-m_WeightRatio[3][1] );
1945 combine2->SetB(-1.);
1947 typename CoefficientImageType::Pointer bc1Row=combine2->GetOutput();
1949 // Paste middleRow at index 3 (padded coeff)
1950 typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1951 paster2->SetDestinationImage(paster1->GetOutput());
1952 paster2->SetSourceImage(bc1Row);
1953 destinationIndex[InputDimension-1]=3;
1954 paster2->SetDestinationIndex(destinationIndex);
1955 paster2->SetSourceRegion(bc1Row->GetLargestPossibleRegion());
1957 //---------------------------------
1958 // 3. Next two temporal rows are identical: paste 3,4 to 4,5
1959 typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1960 paster3->SetDestinationImage(paster2->GetOutput());
1961 paster3->SetSourceImage(m_CoefficientImage);
1962 sourceIndex[InputDimension-1]=3;
1963 destinationIndex[InputDimension-1]=4;
1964 sourceSize[NInputDimensions-1]=2;
1965 sourceRegion.SetSize(sourceSize);
1966 sourceRegion.SetIndex(sourceIndex);
1967 paster3->SetDestinationIndex(destinationIndex);
1968 paster3->SetSourceRegion(sourceRegion);
1970 //---------------------------------
1971 // 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
1972 // Extract row at index 0, 2 (of coeff image)already done
1973 typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1974 typename CoefficientImageType::Pointer row5=this->ExtractTemporalRow(m_CoefficientImage, 5);
1975 typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1978 typename LinearCombinationFilterType::Pointer combine3=LinearCombinationFilterType::New();
1979 combine3->SetInput(0,row0);
1980 combine3->SetInput(1,row2);
1985 typename LinearCombinationFilterType::Pointer combine4=LinearCombinationFilterType::New();
1986 combine4->SetInput(0,row4);
1987 combine4->SetInput(1,row5);
1992 typename LinearCombinationFilterType::Pointer combine5=LinearCombinationFilterType::New();
1993 combine5->SetInput(0,combine3->GetOutput());
1994 combine5->SetInput(1,combine4->GetOutput());
1996 combine5->SetB(-1.);
1999 typename LinearCombinationFilterType::Pointer combine6=LinearCombinationFilterType::New();
2000 combine6->SetInput(0,combine5->GetOutput());
2001 combine6->SetInput(1,row1);
2002 combine6->SetA(m_WeightRatio[2][0]);
2005 typename CoefficientImageType::Pointer bc2Row=combine6->GetOutput();
2007 // Paste BC row at index 6
2008 typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
2009 paster4->SetDestinationImage(paster3->GetOutput());
2010 paster4->SetSourceImage(bc2Row);
2011 destinationIndex[InputDimension-1]=6;
2012 paster4->SetDestinationIndex(destinationIndex);
2013 paster4->SetSourceRegion(bc2Row->GetLargestPossibleRegion());
2015 // Paste last row at index 7
2016 typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
2017 paster5->SetDestinationImage(paster4->GetOutput());
2018 paster5->SetSourceImage(row5);
2019 destinationIndex[InputDimension-1]=7;
2020 paster5->SetDestinationIndex(destinationIndex);
2021 paster5->SetSourceRegion(row5->GetLargestPossibleRegion());
2023 // Update the chain!
2025 m_PaddedCoefficientImage= paster5->GetOutput();
2032 // ----------------------------------------------------------------------
2033 // The sputnik with 5 internal CP T''(0)=T''(10)
2034 // Periodic, constrained to zero at the reference
2035 // at position 3 and one indepent extreme
2036 // Coeff row BC1 0 1 BC2 2 3 BC3 4
2037 // PaddedCoeff R: 0 1 2 3 4 5 6 7
2039 // BC2= -weights[3]/weights[1] ( R2+R5 ) - R4
2041 // ---------------------------------------------------------------------
2043 //---------------------------------
2045 typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
2046 typename CoefficientImageType::Pointer row3=this->ExtractTemporalRow(m_CoefficientImage, 3);
2047 typename CoefficientImageType::Pointer row4=this->ExtractTemporalRow(m_CoefficientImage, 4);
2050 typename LinearCombinationFilterType::Pointer combine0=LinearCombinationFilterType::New();
2051 combine0->SetInput(0,row3);
2052 combine0->SetInput(1,row4);
2055 typename LinearCombinationFilterType::Pointer combine0bis=LinearCombinationFilterType::New();
2056 combine0bis->SetInput(0,combine0->GetOutput());
2057 combine0bis->SetInput(1,row1);
2058 combine0bis->SetA(1.);
2059 combine0bis->SetB(-1.);
2060 combine0bis->Update();
2061 typename CoefficientImageType::Pointer bc1Row=combine0bis->GetOutput();
2063 typename PasteImageFilterType::Pointer paster0=PasteImageFilterType::New();
2064 paster0->SetDestinationImage(m_PaddedCoefficientImage);
2065 paster0->SetSourceImage(bc1Row);
2066 destinationIndex[InputDimension-1]=0;
2067 paster0->SetDestinationIndex(destinationIndex);
2068 paster0->SetSourceRegion(bc1Row->GetLargestPossibleRegion());
2070 //---------------------------------
2071 // 1. Next two temporal rows are identical: paste 0-1 to 1-2
2072 typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
2073 paster1->SetDestinationImage(paster0->GetOutput());
2074 paster1->SetSourceImage(m_CoefficientImage);
2075 sourceRegion.SetIndex(NInputDimensions-1,0);
2076 destinationIndex[InputDimension-1]=1;
2077 sourceRegion.SetSize(NInputDimensions-1,2);
2078 paster1->SetDestinationIndex(destinationIndex);
2079 paster1->SetSourceRegion(sourceRegion);
2081 //---------------------------------
2082 // 2. Middle row at index 3=BC
2083 // Extract row at index 1, 2 (of coeff image)
2084 // typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
2085 typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
2088 typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
2089 combine1->SetInput(0,row1);
2090 combine1->SetInput(1,row3);
2095 typename LinearCombinationFilterType::Pointer combine2=LinearCombinationFilterType::New();
2096 combine2->SetInput(0,combine1->GetOutput());
2097 combine2->SetInput(1,row2);
2098 combine2->SetA(-m_WeightRatio[3][1]);
2099 combine2->SetB(-1.);
2101 typename CoefficientImageType::Pointer bc2Row=combine2->GetOutput();
2103 // Paste middleRow at index 3 (padded coeff)
2104 typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
2105 paster2->SetDestinationImage(paster1->GetOutput());
2106 paster2->SetSourceImage(bc2Row);
2107 destinationIndex[InputDimension-1]=3;
2108 paster2->SetDestinationIndex(destinationIndex);
2109 paster2->SetSourceRegion(bc2Row->GetLargestPossibleRegion());
2111 //---------------------------------
2112 // 3. Next two temporal rows are identical: paste 2,3 to 4,5
2113 typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
2114 paster3->SetDestinationImage(paster2->GetOutput());
2115 paster3->SetSourceImage(m_CoefficientImage);
2116 sourceRegion.SetIndex(InputDimension-1,2);
2117 destinationIndex[InputDimension-1]=4;
2118 sourceRegion.SetSize(NInputDimensions-1,2);
2119 paster3->SetDestinationIndex(destinationIndex);
2120 paster3->SetSourceRegion(sourceRegion);
2122 //---------------------------------
2123 // 4. Final row at index 6=BC3
2124 typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
2127 // typename LinearCombinationFilterType::Pointer combine3=LinearCombinationFilterType::New();
2128 // combine3->SetInput(0,row0);
2129 // combine3->SetInput(1,row1);
2130 // combine3->SetA(1.);
2131 // combine3->SetB(0.5);
2132 // combine3->Update();
2133 // typename CoefficientImageType::Pointer bc3Row=combine3->GetOutput();
2135 // Paste BC row at index 6
2136 typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
2137 paster4->SetDestinationImage(paster3->GetOutput());
2138 paster4->SetSourceImage(row0);
2139 destinationIndex[InputDimension-1]=6;
2140 paster4->SetDestinationIndex(destinationIndex);
2141 paster4->SetSourceRegion(row0->GetLargestPossibleRegion());
2143 // Paste row 4 at index 7
2144 typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
2145 paster5->SetDestinationImage(paster4->GetOutput());
2146 paster5->SetSourceImage(row4);
2147 destinationIndex[InputDimension-1]=7;
2148 paster5->SetDestinationIndex(destinationIndex);
2149 paster5->SetSourceRegion(row4->GetLargestPossibleRegion());
2151 // Update the chain!
2153 m_PaddedCoefficientImage= paster5->GetOutput();
2159 DD ("Shape not available");
2165 // // Extract coefficients from padded version
2166 // template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
2168 // ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2169 // ::ExtractCoefficientImage( )
2171 // ////DD("extract coeff image");
2172 // typename ExtractImageFilterType::Pointer extractImageFilter = ExtractImageFilterType::New();
2173 // typename CoefficientImageType::RegionType extractionRegion=m_PaddedCoefficientImage->GetLargestPossibleRegion();
2174 // typename CoefficientImageType::RegionType::SizeType extractionSize=extractionRegion.GetSize();
2175 // typename CoefficientImageType::RegionType::IndexType extractionIndex = extractionRegion.GetIndex();
2176 // extractionSize[InputDimension-1]-=4;
2177 // extractionIndex[InputDimension-1]=2;
2178 // extractionRegion.SetSize(extractionSize);
2179 // extractionRegion.SetIndex(extractionIndex);
2180 // extractImageFilter->SetInput(m_PaddedCoefficientImage);
2181 // extractImageFilter->SetExtractionRegion(extractionRegion);
2182 // extractImageFilter->Update();
2183 // m_CoefficientImage=extractImageFilter->GetOutput();
2188 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
2190 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2191 ::PrintSelf(std::ostream &os, itk::Indent indent) const
2194 this->Superclass::PrintSelf(os, indent);
2196 os << indent << "GridRegion: " << m_GridRegion << std::endl;
2197 os << indent << "GridOrigin: " << m_GridOrigin << std::endl;
2198 os << indent << "GridSpacing: " << m_GridSpacing << std::endl;
2199 os << indent << "GridDirection: " << m_GridDirection << std::endl;
2200 os << indent << "IndexToPoint: " << m_IndexToPoint << std::endl;
2201 os << indent << "PointToIndex: " << m_PointToIndex << std::endl;
2203 os << indent << "CoefficientImage: [ ";
2204 os << m_CoefficientImage.GetPointer() << " ]" << std::endl;
2206 os << indent << "WrappedImage: [ ";
2207 os << m_WrappedImage.GetPointer() << " ]" << std::endl;
2209 os << indent << "InputParametersPointer: "
2210 << m_InputParametersPointer << std::endl;
2211 os << indent << "ValidRegion: " << m_ValidRegion << std::endl;
2212 os << indent << "LastJacobianIndex: " << m_LastJacobianIndex << std::endl;
2213 os << indent << "BulkTransform: ";
2214 os << m_BulkTransform << std::endl;
2216 if ( m_BulkTransform )
2218 os << indent << "BulkTransformType: "
2219 << m_BulkTransform->GetNameOfClass() << std::endl;
2221 os << indent << "VectorBSplineInterpolator: ";
2222 os << m_VectorInterpolator.GetPointer() << std::endl;
2223 os << indent << "Mask: ";
2224 os << m_Mask<< std::endl;
2229 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
2231 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2232 ::InsideValidRegion( const ContinuousIndexType& index ) const
2236 if ( !m_ValidRegion.IsInside( index ) )
2241 // JV verify all dimensions
2244 typedef typename ContinuousIndexType::ValueType ValueType;
2245 for( unsigned int j = 0; j < NInputDimensions; j++ )
2247 if (m_SplineOrderOdd[j])
2249 if ( index[j] >= static_cast<ValueType>( m_ValidRegionLast[j] ) )
2261 // Transform a point
2262 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
2263 typename ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2265 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2266 ::TransformPoint(const InputPointType &inputPoint) const
2269 InputPointType transformedPoint;
2270 OutputPointType outputPoint;
2273 if ( m_BulkTransform )
2275 transformedPoint = m_BulkTransform->TransformPoint( inputPoint );
2279 transformedPoint = inputPoint;
2282 // Deformable transform
2283 if ( m_PaddedCoefficientImage )
2285 // Check if inside mask
2286 if(m_Mask && !(m_Mask->IsInside(inputPoint) ) )
2288 // Outside: no (deformable) displacement
2289 return transformedPoint;
2292 // Check if inside valid region
2294 ContinuousIndexType index;
2295 this->TransformPointToContinuousIndex( inputPoint, index );
2296 inside = this->InsideValidRegion( index );
2299 // Outside: no (deformable) displacement
2300 DD("outside valid region");
2301 return transformedPoint;
2304 // Call the vector interpolator
2305 itk::Vector<TCoordRep,SpaceDimension> displacement=m_VectorInterpolator->EvaluateAtContinuousIndex(index);
2307 // JV add for the spatial dimensions
2308 outputPoint=transformedPoint;
2309 for (unsigned int i=0; i<NInputDimensions-1; i++)
2310 outputPoint[i] += displacement[i];
2316 itkWarningMacro( << "B-spline coefficients have not been set" );
2317 outputPoint = transformedPoint;
2325 //JV Deformably transform a point
2326 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
2327 typename ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2329 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2330 ::DeformablyTransformPoint(const InputPointType &inputPoint) const
2332 OutputPointType outputPoint;
2333 if ( m_PaddedCoefficientImage )
2336 // Check if inside mask
2337 if(m_Mask && !(m_Mask->IsInside(inputPoint) ) )
2339 // Outside: no (deformable) displacement
2345 ContinuousIndexType index;
2346 this->TransformPointToContinuousIndex( inputPoint, index );
2347 inside = this->InsideValidRegion( index );
2351 //outside: no (deformable) displacement
2352 outputPoint = inputPoint;
2356 // Call the vector interpolator
2357 itk::Vector<TCoordRep,SpaceDimension> displacement=m_VectorInterpolator->EvaluateAtContinuousIndex(index);
2359 // JV add for the spatial dimensions
2360 outputPoint=inputPoint;
2361 for (unsigned int i=0; i<NInputDimensions-1; i++)
2362 outputPoint[i] += displacement[i];
2365 // No coefficients available
2368 itkWarningMacro( << "B-spline coefficients have not been set" );
2369 outputPoint = inputPoint;
2376 // JV weights are identical as for transformpoint, could be done simultaneously in metric!!!!
2377 // Compute the Jacobian in one position
2378 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
2380 typename ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2382 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2383 ::GetJacobian( const InputPointType & point ) const
2386 //========================================================
2387 // Zero all components of jacobian
2388 //========================================================
2389 // JV not thread safe (m_LastJacobianIndex), instantiate N transforms
2390 // NOTE: for efficiency, we only need to zero out the coefficients
2391 // that got fill last time this method was called.
2392 unsigned int j=0,b=0;
2394 //Set the previously-set to zero
2395 for ( j = 0; j < SpaceDimension; j++ )
2397 m_FirstIterator[j].GoToBegin();
2398 while ( !m_FirstIterator[j].IsAtEnd() )
2400 m_FirstIterator[j].Set( m_ZeroVector );
2401 ++(m_FirstIterator[j]);
2405 //Set the previously-set to zero
2406 for ( j = 0; j < SpaceDimension; j++ )
2408 m_SecondIterator[j].GoToBegin();
2409 while ( ! (m_SecondIterator[j]).IsAtEnd() )
2411 m_SecondIterator[j].Set( m_ZeroVector );
2412 ++(m_SecondIterator[j]);
2416 //Set the previously-set to zero
2418 for ( j = 0; j < SpaceDimension; j++ )
2420 m_ThirdIterator[j].GoToBegin();
2421 while ( ! (m_ThirdIterator[j]).IsAtEnd() )
2423 m_ThirdIterator[j].Set( m_ZeroVector );
2424 ++(m_ThirdIterator[j]);
2428 //Set the previously-set to zero
2430 for (b=0; b<m_BCSize;b++)
2431 if ( !( m_FirstRegion.IsInside(m_BCRegions[b]) | m_SecondRegion.IsInside(m_BCRegions[b]) ) )
2432 for ( j = 0; j < SpaceDimension; j++ )
2434 m_BCIterators[j][b].GoToBegin();
2435 while ( ! (m_BCIterators[j][b]).IsAtEnd() )
2437 m_BCIterators[j][b].Set( m_ZeroVector );
2438 ++(m_BCIterators[j][b]);
2442 //Set the previously-set to zero
2444 for (b=0; b<m_BC2Size;b++)
2445 if ( !( m_FirstRegion.IsInside(m_BC2Regions[b]) | m_SecondRegion.IsInside(m_BC2Regions[b]) ) )
2446 for ( j = 0; j < SpaceDimension; j++ )
2448 m_BC2Iterators[j][b].GoToBegin();
2449 while ( ! (m_BC2Iterators[j][b]).IsAtEnd() )
2451 m_BC2Iterators[j][b].Set( m_ZeroVector );
2452 ++(m_BC2Iterators[j][b]);
2456 //Set the previously-set to zero
2458 for (b=0; b<m_BC3Size;b++)
2459 if ( !( m_FirstRegion.IsInside(m_BC3Regions[b]) | m_SecondRegion.IsInside(m_BC3Regions[b]) ) )
2460 for ( j = 0; j < SpaceDimension; j++ )
2462 m_BC3Iterators[j][b].GoToBegin();
2463 while ( ! (m_BC3Iterators[j][b]).IsAtEnd() )
2465 m_BC3Iterators[j][b].Set( m_ZeroVector );
2466 ++(m_BC3Iterators[j][b]);
2471 //========================================================
2472 // For each dimension, copy the weight to the support region
2473 //========================================================
2475 // Check if inside mask
2476 if(m_Mask && !(m_Mask->IsInside(point) ) )
2478 // Outside: no (deformable) displacement
2479 return this->m_Jacobian;
2483 this->TransformPointToContinuousIndex( point, m_Index );
2485 // NOTE: if the support region does not lie totally within the grid
2486 // we assume zero displacement and return the input point
2487 if ( !this->InsideValidRegion( m_Index ) )
2489 return this->m_Jacobian;
2492 // Compute interpolation weights
2493 const WeightsDataType *weights=NULL;
2494 m_VectorInterpolator->EvaluateWeightsAtContinuousIndex( m_Index, &weights, m_LastJacobianIndex);
2497 m_SupportRegion.SetIndex( m_LastJacobianIndex );
2498 WrapRegion(m_SupportRegion, m_FirstRegion, m_SecondRegion, m_ThirdRegion, m_BCRegions, m_BCValues, m_BC2Regions, m_BC2Values, m_BC3Regions, m_BC3Values, m_InitialOffset);
2499 m_ThirdSize=m_ThirdRegion.GetSize()[InputDimension-1];
2500 m_BCSize=m_BCRegions.size();
2501 m_BC2Size=m_BC2Regions.size();
2502 m_BC3Size=m_BC3Regions.size();
2504 // Reset the iterators
2505 for ( j = 0; j < SpaceDimension ; j++ )
2507 m_FirstIterator[j] = IteratorType( m_JacobianImage[j], m_FirstRegion);
2508 m_SecondIterator[j] = IteratorType( m_JacobianImage[j], m_SecondRegion);
2509 if(m_ThirdSize) m_ThirdIterator[j] = IteratorType( m_JacobianImage[j], m_ThirdRegion);
2511 m_BCIterators[j].resize(m_BCSize);
2512 for (b=0; b<m_BCSize;b++)
2513 m_BCIterators[j][b]= IteratorType( m_JacobianImage[j], m_BCRegions[b]);
2514 m_BC2Iterators[j].resize(m_BC2Size);
2515 for (b=0; b<m_BC2Size;b++)
2516 m_BC2Iterators[j][b]= IteratorType( m_JacobianImage[j], m_BC2Regions[b]);
2517 m_BC3Iterators[j].resize(m_BC3Size);
2518 for (b=0; b<m_BC3Size;b++)
2519 m_BC3Iterators[j][b]= IteratorType( m_JacobianImage[j], m_BC3Regions[b]);
2522 // Skip if on a fixed condition
2525 if (m_BCSize) weights+=m_InitialOffset*m_BCRegions[0].GetNumberOfPixels();
2526 else std::cerr<<"InitialOffset without BCSize: Error!!!!!!!"<<std::endl;
2529 //copy weight to jacobian image
2530 for ( j = 0; j < SpaceDimension; j++ )
2532 // For each dimension, copy the weight to the support region
2533 while ( ! (m_FirstIterator[j]).IsAtEnd() )
2535 m_ZeroVector[j]=*weights;
2536 (m_FirstIterator[j]).Set( m_ZeroVector);
2537 ++(m_FirstIterator[j]);
2541 m_ZeroVector[j]=itk::NumericTraits<JacobianValueType>::Zero;
2542 if (j != SpaceDimension-1) weights-=m_FirstRegion.GetNumberOfPixels();
2545 // Skip BC1 and go to the second region
2546 if (m_BCSize) weights+=m_BCRegions[0].GetNumberOfPixels();
2548 // For each dimension, copy the weight to the support region
2549 //copy weight to jacobian image
2550 for ( j = 0; j < SpaceDimension; j++ )
2552 while ( ! (m_SecondIterator[j]).IsAtEnd() )
2554 m_ZeroVector[j]=*weights;
2555 (m_SecondIterator[j]).Set( m_ZeroVector);
2556 ++(m_SecondIterator[j]);
2559 weights-=m_SecondRegion.GetNumberOfPixels();
2560 m_ZeroVector[j]=itk::NumericTraits<JacobianValueType>::Zero;
2566 // Put pointer in correct position
2567 weights-=m_BCRegions[0].GetNumberOfPixels();
2569 for ( j = 0; j < SpaceDimension; j++ )
2571 for ( b=0; b < m_BCSize; b++ )
2573 while ( ! (m_BCIterators[j][b]).IsAtEnd() )
2575 //copy weight to jacobian image
2576 m_ZeroVector[j]=(*weights) * m_BCValues[b];
2577 (m_BCIterators[j][b]).Value()+= m_ZeroVector;
2578 ++(m_BCIterators[j][b]);
2581 weights-=m_BCRegions[b].GetNumberOfPixels();
2583 m_ZeroVector[j]=itk::NumericTraits<JacobianValueType>::Zero;
2585 weights+=m_BCRegions[m_BCSize-1].GetNumberOfPixels();
2588 // Add the BC2 to the weights
2591 // Move further in the weights pointer
2592 weights+=m_SecondRegion.GetNumberOfPixels();
2594 for ( j = 0; j < SpaceDimension; j++ )
2596 for ( b=0; b < m_BC2Size; b++ )
2598 while ( ! (m_BC2Iterators[j][b]).IsAtEnd() )
2600 //copy weight to jacobian image
2601 m_ZeroVector[j]=(*weights) * m_BC2Values[b];
2602 (m_BC2Iterators[j][b]).Value()+= m_ZeroVector;
2603 ++(m_BC2Iterators[j][b]);
2606 weights-=m_BC2Regions[b].GetNumberOfPixels();
2608 m_ZeroVector[j]=itk::NumericTraits<JacobianValueType>::Zero;
2610 // Move further in the weights pointer
2611 weights+=m_BC2Regions[m_BC2Size-1].GetNumberOfPixels();
2617 // For each dimension, copy the weight to the support region
2618 //copy weight to jacobian image
2619 for ( j = 0; j < SpaceDimension; j++ )
2621 while ( ! (m_ThirdIterator[j]).IsAtEnd() )
2623 m_ZeroVector[j]=*weights;
2624 (m_ThirdIterator[j]).Value()+= m_ZeroVector;
2625 ++(m_ThirdIterator[j]);
2629 // Move further in the weights pointer?
2630 if (j != SpaceDimension-1) weights-=m_ThirdRegion.GetNumberOfPixels();
2631 m_ZeroVector[j]=itk::NumericTraits<JacobianValueType>::Zero;
2635 // Add the BC3 to the weights
2638 for ( j = 0; j < SpaceDimension; j++ )
2640 for ( b=0; b < m_BC3Size; b++ )
2642 while ( ! (m_BC3Iterators[j][b]).IsAtEnd() )
2644 //copy weight to jacobian image
2645 m_ZeroVector[j]=(*weights) * m_BC3Values[b];
2646 (m_BC3Iterators[j][b]).Value()+= m_ZeroVector;
2647 ++(m_BC3Iterators[j][b]);
2650 weights-=m_BC3Regions[b].GetNumberOfPixels();
2652 m_ZeroVector[j]=itk::NumericTraits<JacobianValueType>::Zero;
2656 // Return the result
2657 return this->m_Jacobian;
2661 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
2663 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2664 ::WrapRegion( const RegionType & m_SupportRegion,
2665 RegionType & m_FirstRegion,
2666 RegionType & m_SecondRegion,
2667 RegionType & m_ThirdRegion,
2668 std::vector<RegionType>& m_BCRegions,std::vector<double>& m_BCValues,
2669 std::vector<RegionType>& m_BC2Regions,std::vector<double>& m_BC2Values,
2670 std::vector<RegionType>& m_BC3Regions,std::vector<double>& m_BC3Values,
2671 unsigned int& m_InitialOffset ) const
2674 // Synchronize regions
2676 m_FirstRegion=m_SupportRegion;
2677 m_BCRegion=m_SupportRegion;
2678 m_BCRegion.SetSize(InputDimension-1,1);
2679 m_SecondRegion=m_SupportRegion;
2680 m_ThirdRegion=m_SupportRegion;
2681 m_ThirdRegion.SetSize(InputDimension-1,0);
2682 m_BC3Regions.resize(0);
2685 // BC depends on shape
2686 switch(m_TransformShape)
2691 2: rabbit 4 CP 3 DOF
2692 3: rabbit 5 CP 4 DOF
2693 4: sputnik 4 CP 4 DOF
2694 5: sputnik 5 CP 5 DOF
2695 6: diamond 6 CP 5 DOF
2696 7: diamond 7 CP 6 DOF
2701 // ----------------------------------------------------------------------
2702 // The egg with 4 internal CP (starting from inhale)
2703 // Periodic, constrained to zero at the reference
2704 // at position 3 and
2705 // Coeff row BC1 BC2 0 BC3 1 2 BC4
2706 // PaddedCoeff R: 0 1 2 3 4 5 6
2709 // BC3= -weights[2]/weights[0] ( R2+R4 )
2711 // ---------------------------------------------------------------------
2712 switch(m_SupportRegion.GetIndex(InputDimension-1))
2717 m_FirstRegion.SetSize(InputDimension-1,2);
2718 m_FirstRegion.SetIndex(InputDimension-1,1);
2721 m_BCRegions.resize(0);
2724 m_SecondRegion.SetSize(InputDimension-1,1);
2725 m_SecondRegion.SetIndex(InputDimension-1,0);
2728 m_BC2Regions.resize(2);
2729 m_BC2Values.resize(2);
2730 m_BCRegion.SetIndex(InputDimension-1,0);
2731 m_BC2Regions[0]=m_BCRegion;
2732 m_BC2Values[0]=-m_WeightRatio[2][0];
2733 m_BCRegion.SetIndex(InputDimension-1,1);
2734 m_BC2Regions[1]=m_BCRegion;
2735 m_BC2Values[1]= -m_WeightRatio[2][0];
2742 m_FirstRegion.SetSize(InputDimension-1,1);
2743 m_FirstRegion.SetIndex(InputDimension-1,2);
2746 m_BCRegions.resize(0);
2749 m_SecondRegion.SetSize(InputDimension-1,1);
2750 m_SecondRegion.SetIndex(InputDimension-1,0);
2753 m_BC2Regions.resize(2);
2754 m_BC2Values.resize(2);
2755 m_BCRegion.SetIndex(InputDimension-1,0);
2756 m_BC2Regions[0]=m_BCRegion;
2757 m_BC2Values[0]=-m_WeightRatio[2][0];
2758 m_BCRegion.SetIndex(InputDimension-1,1);
2759 m_BC2Regions[1]=m_BCRegion;
2760 m_BC2Values[1]= -m_WeightRatio[2][0];
2763 m_FirstRegion.SetSize(InputDimension-1,1);
2764 m_FirstRegion.SetIndex(InputDimension-1,1);
2771 m_FirstRegion.SetSize(InputDimension-1,1);
2772 m_FirstRegion.SetIndex(InputDimension-1,0);
2775 m_BCRegions.resize(2);
2776 m_BCValues.resize(2);
2777 m_BCRegion.SetIndex(InputDimension-1,0);
2778 m_BCRegions[0]=m_BCRegion;
2779 m_BCValues[0]=-m_WeightRatio[2][0];
2780 m_BCRegion.SetIndex(InputDimension-1,1);
2781 m_BCRegions[1]=m_BCRegion;
2782 m_BCValues[1]= -m_WeightRatio[2][0];
2785 m_SecondRegion.SetSize(InputDimension-1,2);
2786 m_SecondRegion.SetIndex(InputDimension-1,1);
2789 m_BC2Regions.resize(0);
2796 m_FirstRegion.SetSize(InputDimension-1,0);
2799 m_BCRegions.resize(2);
2800 m_BCValues.resize(2);
2801 m_BCRegion.SetIndex(InputDimension-1,0);
2802 m_BCRegions[0]=m_BCRegion;
2803 m_BCValues[0]=-m_WeightRatio[2][0];
2804 m_BCRegion.SetIndex(InputDimension-1,1);
2805 m_BCRegions[1]=m_BCRegion;
2806 m_BCValues[1]= -m_WeightRatio[2][0];
2809 m_SecondRegion.SetSize(InputDimension-1,2);
2810 m_SecondRegion.SetIndex(InputDimension-1,1);
2813 m_BC2Regions.resize(1);
2814 m_BC2Values.resize(1);
2815 m_BCRegion.SetIndex(InputDimension-1,0);
2816 m_BC2Regions[0]=m_BCRegion;
2824 DD("supportindex > 3 ???");
2825 DD(m_SupportRegion.GetIndex(InputDimension-1));
2827 } // end switch index
2834 // ----------------------------------------------------------------------
2835 // The egg with 5 internal CP (starting from inhale)
2836 // Periodic, constrained to zero at the reference
2837 // at position 3 and
2838 // Coeff row BC1 BC2 0 BC3 1 2 3 BC4
2839 // PaddedCoeff R: 0 1 2 3 4 5 6 7
2842 // BC3= -weights[3]/weights[1] ( R2+R5 ) - R4
2844 // ---------------------------------------------------------------------
2845 switch(m_SupportRegion.GetIndex(InputDimension-1))
2850 m_FirstRegion.SetSize(InputDimension-1,2);
2851 m_FirstRegion.SetIndex(InputDimension-1,2);
2854 m_BCRegions.resize(0);
2857 m_SecondRegion.SetSize(InputDimension-1,1);
2858 m_SecondRegion.SetIndex(InputDimension-1,0);
2861 m_BC2Regions.resize(3);
2862 m_BC2Values.resize(3);
2863 m_BCRegion.SetIndex(InputDimension-1,0);
2864 m_BC2Regions[0]=m_BCRegion;
2865 m_BC2Values[0]=-m_WeightRatio[3][1];
2866 m_BCRegion.SetIndex(InputDimension-1,1);
2867 m_BC2Regions[1]=m_BCRegion;
2869 m_BCRegion.SetIndex(InputDimension-1,2);
2870 m_BC2Regions[2]=m_BCRegion;
2871 m_BC2Values[2]=-m_WeightRatio[3][1];
2878 m_FirstRegion.SetSize(InputDimension-1,1);
2879 m_FirstRegion.SetIndex(InputDimension-1,3);
2882 m_BCRegions.resize(0);
2885 m_SecondRegion.SetSize(InputDimension-1,1);
2886 m_SecondRegion.SetIndex(InputDimension-1,0);
2889 m_BC2Regions.resize(3);
2890 m_BC2Values.resize(3);
2891 m_BCRegion.SetIndex(InputDimension-1,0);
2892 m_BC2Regions[0]=m_BCRegion;
2893 m_BC2Values[0]=-m_WeightRatio[3][1];
2894 m_BCRegion.SetIndex(InputDimension-1,1);
2895 m_BC2Regions[1]=m_BCRegion;
2897 m_BCRegion.SetIndex(InputDimension-1,2);
2898 m_BC2Regions[2]=m_BCRegion;
2899 m_BC2Values[2]=-m_WeightRatio[3][1];
2902 m_FirstRegion.SetSize(InputDimension-1,1);
2903 m_FirstRegion.SetIndex(InputDimension-1,1);
2910 m_FirstRegion.SetSize(InputDimension-1,1);
2911 m_FirstRegion.SetIndex(InputDimension-1,0);
2914 m_BCRegions.resize(3);
2915 m_BCValues.resize(3);
2916 m_BCRegion.SetIndex(InputDimension-1,0);
2917 m_BCRegions[0]=m_BCRegion;
2918 m_BCValues[0]=-m_WeightRatio[3][1];
2919 m_BCRegion.SetIndex(InputDimension-1,1);
2920 m_BCRegions[1]=m_BCRegion;
2922 m_BCRegion.SetIndex(InputDimension-1,2);
2923 m_BCRegions[2]=m_BCRegion;
2924 m_BCValues[2]=-m_WeightRatio[3][1];
2927 m_SecondRegion.SetSize(InputDimension-1,2);
2928 m_SecondRegion.SetIndex(InputDimension-1,1);
2931 m_BC2Regions.resize(0);
2938 m_FirstRegion.SetSize(InputDimension-1,0);
2941 m_BCRegions.resize(3);
2942 m_BCValues.resize(3);
2943 m_BCRegion.SetIndex(InputDimension-1,0);
2944 m_BCRegions[0]=m_BCRegion;
2945 m_BCValues[0]=-m_WeightRatio[3][1];
2946 m_BCRegion.SetIndex(InputDimension-1,1);
2947 m_BCRegions[1]=m_BCRegion;
2949 m_BCRegion.SetIndex(InputDimension-1,2);
2950 m_BCRegions[2]=m_BCRegion;
2951 m_BCValues[2]=-m_WeightRatio[3][1];
2954 m_SecondRegion.SetSize(InputDimension-1,3);
2955 m_SecondRegion.SetIndex(InputDimension-1,1);
2958 m_BC2Regions.resize(0);
2965 m_FirstRegion.SetSize(InputDimension-1,3);
2966 m_FirstRegion.SetIndex(InputDimension-1,1);
2969 m_BCRegions.resize(0);
2972 m_SecondRegion.SetSize(InputDimension-1,1);
2973 m_SecondRegion.SetIndex(InputDimension-1,0);
2976 m_BC2Regions.resize(0);
2983 DD("supportindex > 3 ???");
2984 DD(m_SupportRegion.GetIndex(InputDimension-1));
2986 } // end switch index
2990 // // ----------------------------------------------------------------------
2991 // // The egg with 5 internal CP:
2992 // // Periodic, C2 smooth everywhere and constrained to zero at the reference
2993 // // Coeff row R5 BC1 0 1 2 3 BC2 R2
2994 // // PaddedCoeff R: 0 1 2 3 4 5 6 7
2995 // // BC1= -weights[2]/weights[0] ( R2+R5)
2997 // // ---------------------------------------------------------------------
2998 // switch(m_SupportRegion.GetIndex(InputDimension-1))
3003 // m_FirstRegion.SetSize(InputDimension-1,1);
3004 // m_FirstRegion.SetIndex(InputDimension-1,3);
3007 // m_BCRegions.resize(2);
3008 // m_BCValues.resize(2);
3009 // m_BCRegion.SetIndex(InputDimension-1,0);
3010 // m_BCRegions[0]=m_BCRegion;
3011 // m_BCValues[0]=-m_WeightRatio[2][0];
3012 // m_BCRegion.SetIndex(InputDimension-1,3);
3013 // m_BCRegions[1]=m_BCRegion;
3014 // m_BCValues[1]=-m_WeightRatio[2][0];
3017 // m_SecondRegion.SetSize(InputDimension-1,2);
3018 // m_SecondRegion.SetIndex(InputDimension-1,0);
3021 // m_BC2Regions.resize(0);
3028 // m_FirstRegion.SetSize(InputDimension-1,0);
3031 // m_BCRegions.resize(2);
3032 // m_BCValues.resize(2);
3033 // m_BCRegion.SetIndex(InputDimension-1,0);
3034 // m_BCRegions[0]=m_BCRegion;
3035 // m_BCValues[0]=-m_WeightRatio[2][0];
3036 // m_BCRegion.SetIndex(InputDimension-1,3);
3037 // m_BCRegions[1]=m_BCRegion;
3038 // m_BCValues[1]=-m_WeightRatio[2][0];
3041 // m_SecondRegion.SetSize(InputDimension-1,3);
3042 // m_SecondRegion.SetIndex(InputDimension-1,0);
3045 // m_BC2Regions.resize(0);
3052 // m_FirstRegion.SetSize(InputDimension-1,4);
3053 // m_FirstRegion.SetIndex(InputDimension-1,0);
3056 // m_BCRegions.resize(0);
3059 // m_SecondRegion.SetSize(InputDimension-1,0);
3062 // m_BC2Regions.resize(0);
3069 // m_FirstRegion.SetSize(InputDimension-1,3);
3070 // m_FirstRegion.SetIndex(InputDimension-1,1);
3073 // m_BCRegions.resize(2);
3074 // m_BCValues.resize(2);
3075 // m_BCRegion.SetIndex(InputDimension-1,0);
3076 // m_BCRegions[0]=m_BCRegion;
3077 // m_BCValues[0]=-m_WeightRatio[2][0];
3078 // m_BCRegion.SetIndex(InputDimension-1,3);
3079 // m_BCRegions[1]=m_BCRegion;
3080 // m_BCValues[1]=-m_WeightRatio[2][0];
3083 // m_SecondRegion.SetSize(InputDimension-1,0);
3086 // m_BC2Regions.resize(0);
3094 // m_FirstRegion.SetSize(InputDimension-1,2);
3095 // m_FirstRegion.SetIndex(InputDimension-1,2);
3098 // m_BCRegions.resize(2);
3099 // m_BCValues.resize(2);
3100 // m_BCRegion.SetIndex(InputDimension-1,0);
3101 // m_BCRegions[0]=m_BCRegion;
3102 // m_BCValues[0]=-m_WeightRatio[2][0];
3103 // m_BCRegion.SetIndex(InputDimension-1,3);
3104 // m_BCRegions[1]=m_BCRegion;
3105 // m_BCValues[1]=-m_WeightRatio[2][0];
3108 // m_SecondRegion.SetSize(InputDimension-1,1);
3109 // m_SecondRegion.SetIndex(InputDimension-1,0);
3112 // m_BC2Regions.resize(0);
3119 // DD("supportindex > 4 ???");
3120 // DD(m_SupportRegion.GetIndex(InputDimension-1));
3121 // DD(m_TransformShape);
3123 // }// end swith index
3126 // } // end case 1 shape
3130 // ----------------------------------------------------------------------
3131 // The rabbit with 4 internal CP
3132 // Periodic, constrained to zero at the reference
3133 // at position 3 and the extremes fixed with anit-symm bc
3134 // Coeff row BC1 0 1 BC2 2 BC3 BC4
3135 // PaddedCoeff R: 0 1 2 3 4 5 6
3137 // BC2= -weights[2]/weights[0] ( R2+R4 )
3140 // ---------------------------------------------------------------------
3141 switch(m_SupportRegion.GetIndex(InputDimension-1))
3146 m_FirstRegion.SetSize(InputDimension-1,0);
3149 m_BCRegions.resize(2);
3150 m_BCValues.resize(2);
3151 m_BCRegion.SetIndex(InputDimension-1,0);
3152 m_BCRegions[0]=m_BCRegion;
3154 m_BCRegion.SetIndex(InputDimension-1,1);
3155 m_BCRegions[1]=m_BCRegion;
3159 m_SecondRegion.SetSize(InputDimension-1,2);
3160 m_SecondRegion.SetIndex(InputDimension-1,0);
3163 m_BC2Regions.resize(2);
3164 m_BC2Values.resize(2);
3165 m_BCRegion.SetIndex(InputDimension-1,1);
3166 m_BC2Regions[0]=m_BCRegion;
3167 m_BC2Values[0]=-m_WeightRatio[2][0];
3168 m_BCRegion.SetIndex(InputDimension-1,2);
3169 m_BC2Regions[1]=m_BCRegion;
3170 m_BC2Values[1]= -m_WeightRatio[2][0];
3177 m_FirstRegion.SetSize(InputDimension-1,2);
3178 m_FirstRegion.SetIndex(InputDimension-1,0);
3181 m_BCRegions.resize(2);
3182 m_BCValues.resize(2);
3183 m_BCRegion.SetIndex(InputDimension-1,1);
3184 m_BCRegions[0]=m_BCRegion;
3185 m_BCValues[0]=-m_WeightRatio[2][0];
3186 m_BCRegion.SetIndex(InputDimension-1,2);
3187 m_BCRegions[1]=m_BCRegion;
3188 m_BCValues[1]= -m_WeightRatio[2][0];
3191 m_SecondRegion.SetSize(InputDimension-1,1);
3192 m_SecondRegion.SetIndex(InputDimension-1,2);
3195 m_BC2Regions.resize(0);
3202 m_FirstRegion.SetSize(InputDimension-1,1);
3203 m_FirstRegion.SetIndex(InputDimension-1,1);
3206 m_BCRegions.resize(2);
3207 m_BCValues.resize(2);
3208 m_BCRegion.SetIndex(InputDimension-1,1);
3209 m_BCRegions[0]=m_BCRegion;
3210 m_BCValues[0]=-m_WeightRatio[2][0];
3211 m_BCRegion.SetIndex(InputDimension-1,2);
3212 m_BCRegions[1]=m_BCRegion;
3213 m_BCValues[1]= -m_WeightRatio[2][0];
3216 m_SecondRegion.SetSize(InputDimension-1,1);
3217 m_SecondRegion.SetIndex(InputDimension-1,2);
3220 m_BC2Regions.resize(1);
3221 m_BC2Values.resize(1);
3222 m_BCRegion.SetIndex(InputDimension-1,0);
3223 m_BC2Regions[0]=m_BCRegion;
3231 m_FirstRegion.SetSize(InputDimension-1,0);
3234 m_BCRegions.resize(2);
3235 m_BCValues.resize(2);
3236 m_BCRegion.SetIndex(InputDimension-1,1);
3237 m_BCRegions[0]=m_BCRegion;
3238 m_BCValues[0]=-m_WeightRatio[2][0];
3239 m_BCRegion.SetIndex(InputDimension-1,2);
3240 m_BCRegions[1]=m_BCRegion;
3241 m_BCValues[1]= -m_WeightRatio[2][0];
3244 m_SecondRegion.SetSize(InputDimension-1,1);
3245 m_SecondRegion.SetIndex(InputDimension-1,2);
3248 m_BC2Regions.resize(1);
3249 m_BC2Values.resize(1);
3250 m_BCRegion.SetIndex(InputDimension-1,0);
3251 m_BC2Regions[0]=m_BCRegion;
3255 m_BC3Regions.resize(2);
3256 m_BC3Values.resize(2);
3257 m_BCRegion.SetIndex(InputDimension-1,0);
3258 m_BC3Regions[0]=m_BCRegion;
3260 m_BCRegion.SetIndex(InputDimension-1,2);
3261 m_BC3Regions[1]=m_BCRegion;
3269 DD("supportindex > 3 ???");
3270 DD(m_SupportRegion.GetIndex(InputDimension-1));
3272 } // end switch index
3275 } // end case 2 shape
3278 // ----------------------------------------------------------------------
3279 // The rabbit with 5 internal CP
3280 // Periodic, constrained to zero at the reference at position 3.5
3281 // and the extremes fixed with anti-symmetrical BC
3282 // Coeff row BC1 0 1 BC2 2 3 BC3 BC4
3283 // PaddedCoeff R: 0 1 2 3 4 5 6 7
3285 // BC2= -weights[3]/weights[1] ( R2+R5 ) - R4
3288 // ---------------------------------------------------------------------
3289 switch(m_SupportRegion.GetIndex(InputDimension-1))
3294 m_FirstRegion.SetSize(InputDimension-1,0);
3297 m_BCRegions.resize(2);
3298 m_BCValues.resize(2);
3299 m_BCRegion.SetIndex(InputDimension-1,0);
3300 m_BCRegions[0]=m_BCRegion;
3302 m_BCRegion.SetIndex(InputDimension-1,1);
3303 m_BCRegions[1]=m_BCRegion;
3307 m_SecondRegion.SetSize(InputDimension-1,2);
3308 m_SecondRegion.SetIndex(InputDimension-1,0);
3311 m_BC2Regions.resize(3);
3312 m_BC2Values.resize(3);
3313 m_BCRegion.SetIndex(InputDimension-1,1);
3314 m_BC2Regions[0]=m_BCRegion;
3315 m_BC2Values[0]=-m_WeightRatio[3][1];
3316 m_BCRegion.SetIndex(InputDimension-1,2);
3317 m_BC2Regions[1]=m_BCRegion;
3319 m_BCRegion.SetIndex(InputDimension-1,3);
3320 m_BC2Regions[2]=m_BCRegion;
3321 m_BC2Values[2]=-m_WeightRatio[3][1];
3328 m_FirstRegion.SetSize(InputDimension-1,2);
3329 m_FirstRegion.SetIndex(InputDimension-1,0);
3332 m_BCRegions.resize(3);
3333 m_BCValues.resize(3);
3334 m_BCRegion.SetIndex(InputDimension-1,1);
3335 m_BCRegions[0]=m_BCRegion;
3336 m_BCValues[0]=-m_WeightRatio[3][1];
3337 m_BCRegion.SetIndex(InputDimension-1,2);
3338 m_BCRegions[1]=m_BCRegion;
3340 m_BCRegion.SetIndex(InputDimension-1,3);
3341 m_BCRegions[2]=m_BCRegion;
3342 m_BCValues[2]=-m_WeightRatio[3][1];
3345 m_SecondRegion.SetSize(InputDimension-1,1);
3346 m_SecondRegion.SetIndex(InputDimension-1,2);
3349 m_BC2Regions.resize(0);
3356 m_FirstRegion.SetSize(InputDimension-1,1);
3357 m_FirstRegion.SetIndex(InputDimension-1,1);
3360 m_BCRegions.resize(3);
3361 m_BCValues.resize(3);
3362 m_BCRegion.SetIndex(InputDimension-1,1);
3363 m_BCRegions[0]=m_BCRegion;
3364 m_BCValues[0]=-m_WeightRatio[3][1];
3365 m_BCRegion.SetIndex(InputDimension-1,2);
3366 m_BCRegions[1]=m_BCRegion;
3368 m_BCRegion.SetIndex(InputDimension-1,3);
3369 m_BCRegions[2]=m_BCRegion;
3370 m_BCValues[2]=-m_WeightRatio[3][1];
3373 m_SecondRegion.SetSize(InputDimension-1,2);
3374 m_SecondRegion.SetIndex(InputDimension-1,2);
3377 m_BC2Regions.resize(0);
3384 m_FirstRegion.SetSize(InputDimension-1,0);
3387 m_BCRegions.resize(3);
3388 m_BCValues.resize(3);
3389 m_BCRegion.SetIndex(InputDimension-1,1);
3390 m_BCRegions[0]=m_BCRegion;
3391 m_BCValues[0]=-m_WeightRatio[3][1];
3392 m_BCRegion.SetIndex(InputDimension-1,2);
3393 m_BCRegions[1]=m_BCRegion;
3395 m_BCRegion.SetIndex(InputDimension-1,3);
3396 m_BCRegions[2]=m_BCRegion;
3397 m_BCValues[2]=-m_WeightRatio[3][1];
3400 m_SecondRegion.SetSize(InputDimension-1,2);
3401 m_SecondRegion.SetIndex(InputDimension-1,2);
3404 m_BC2Regions.resize(1);
3405 m_BC2Values.resize(1);
3406 m_BCRegion.SetIndex(InputDimension-1,0);
3407 m_BC2Regions[0]=m_BCRegion;
3416 m_FirstRegion.SetSize(InputDimension-1,2);
3417 m_FirstRegion.SetIndex(InputDimension-1,2);
3420 m_BCRegions.resize(1);
3421 m_BCValues.resize(1);
3422 m_BCRegion.SetIndex(InputDimension-1,0);
3423 m_BCRegions[0]=m_BCRegion;
3427 m_SecondRegion.SetSize(InputDimension-1,0);
3430 m_BC2Regions.resize(2);
3431 m_BC2Values.resize(2);
3432 m_BCRegion.SetIndex(InputDimension-1,0);
3433 m_BC2Regions[0]=m_BCRegion;
3435 m_BCRegion.SetIndex(InputDimension-1,3);
3436 m_BC2Regions[1]=m_BCRegion;
3444 DD("supportindex > 4 ???");
3445 DD(m_SupportRegion.GetIndex(InputDimension-1));
3447 } // end switch index
3450 } // end case 3 shape
3454 // ----------------------------------------------------------------------
3455 // The sputnik with 4 internal CP
3456 // Periodic, constrained to zero at the reference
3457 // at position 3 and one indepent extremes copied
3458 // Coeff row BC1 0 1 BC2 2 BC2 3
3459 // PaddedCoeff R: 0 1 2 3 4 5 6
3461 // BC2= -weights[2]/weights[0] ( R2+R4 )
3462 // BC3= weights[2]/weights[0] ( R2-R4) + R1
3463 // ---------------------------------------------------------------------
3464 switch(m_SupportRegion.GetIndex(InputDimension-1))
3469 m_FirstRegion.SetSize(InputDimension-1,0);
3472 m_BCRegions.resize(1);
3473 m_BCValues.resize(1);
3474 m_BCRegion.SetIndex(InputDimension-1,3);
3475 m_BCRegions[0]=m_BCRegion;
3479 m_SecondRegion.SetSize(InputDimension-1,2);
3480 m_SecondRegion.SetIndex(InputDimension-1,0);
3483 m_BC2Regions.resize(2);
3484 m_BC2Values.resize(2);
3485 m_BCRegion.SetIndex(InputDimension-1,1);
3486 m_BC2Regions[0]=m_BCRegion;
3487 m_BC2Values[0]=-m_WeightRatio[2][0];
3488 m_BCRegion.SetIndex(InputDimension-1,2);
3489 m_BC2Regions[1]=m_BCRegion;
3490 m_BC2Values[1]= -m_WeightRatio[2][0];
3497 m_FirstRegion.SetSize(InputDimension-1,2);
3498 m_FirstRegion.SetIndex(InputDimension-1,0);
3501 m_BCRegions.resize(2);
3502 m_BCValues.resize(2);
3503 m_BCRegion.SetIndex(InputDimension-1,1);
3504 m_BCRegions[0]=m_BCRegion;
3505 m_BCValues[0]=-m_WeightRatio[2][0];
3506 m_BCRegion.SetIndex(InputDimension-1,2);
3507 m_BCRegions[1]=m_BCRegion;
3508 m_BCValues[1]= -m_WeightRatio[2][0];
3511 m_SecondRegion.SetSize(InputDimension-1,1);
3512 m_SecondRegion.SetIndex(InputDimension-1,2);
3515 m_BC2Regions.resize(0);
3522 m_FirstRegion.SetSize(InputDimension-1,1);
3523 m_FirstRegion.SetIndex(InputDimension-1,1);
3526 m_BCRegions.resize(2);
3527 m_BCValues.resize(2);
3528 m_BCRegion.SetIndex(InputDimension-1,1);
3529 m_BCRegions[0]=m_BCRegion;
3530 m_BCValues[0]=-m_WeightRatio[2][0];
3531 m_BCRegion.SetIndex(InputDimension-1,2);
3532 m_BCRegions[1]=m_BCRegion;
3533 m_BCValues[1]= -m_WeightRatio[2][0];
3536 m_SecondRegion.SetSize(InputDimension-1,1);
3537 m_SecondRegion.SetIndex(InputDimension-1,2);
3540 m_BC2Regions.resize(3);
3541 m_BC2Values.resize(3);
3542 m_BCRegion.SetIndex(InputDimension-1,0);
3543 m_BC2Regions[0]=m_BCRegion;
3545 m_BCRegion.SetIndex(InputDimension-1,1);
3546 m_BC2Regions[1]=m_BCRegion;
3547 m_BC2Values[1]=m_WeightRatio[2][0];
3548 m_BCRegion.SetIndex(InputDimension-1,2);
3549 m_BC2Regions[2]=m_BCRegion;
3550 m_BC2Values[2]=-m_WeightRatio[2][0];
3557 m_FirstRegion.SetSize(InputDimension-1,0);
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];
3570 m_SecondRegion.SetSize(InputDimension-1,1);
3571 m_SecondRegion.SetIndex(InputDimension-1,2);
3574 m_BC2Regions.resize(3);
3575 m_BC2Values.resize(3);
3576 m_BCRegion.SetIndex(InputDimension-1,0);
3577 m_BC2Regions[0]=m_BCRegion;
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];
3587 m_ThirdRegion.SetSize(InputDimension-1,1);
3588 m_ThirdRegion.SetIndex(InputDimension-1,3);
3595 DD("supportindex > 3 ???");
3596 DD(m_SupportRegion.GetIndex(InputDimension-1));
3598 } // end switch index
3601 } // end case 4 shape
3605 // ----------------------------------------------------------------------
3606 // The sputnik with 5 internal CP
3607 // Periodic, constrained to zero at the reference
3608 // at position 3 and one indepent extreme
3609 // Coeff row BC1 0 1 BC2 2 3 BC3 4
3610 // PaddedCoeff R: 0 1 2 3 4 5 6 7
3611 // BC1= R2 + R5 - R7
3612 // BC2= -weights[3]/weights[1] ( R2+R5 ) - R4
3613 // BC3= R1 + 0.5 R2 - 0.5 R7
3614 // ----------------------------------------------------------------------
3615 switch(m_SupportRegion.GetIndex(InputDimension-1))
3620 m_FirstRegion.SetSize(InputDimension-1,0);
3623 m_BCRegions.resize(3);
3624 m_BCValues.resize(3);
3625 m_BCRegion.SetIndex(InputDimension-1,1);
3626 m_BCRegions[0]=m_BCRegion;
3628 m_BCRegion.SetIndex(InputDimension-1,3);
3629 m_BCRegions[1]=m_BCRegion;
3631 m_BCRegion.SetIndex(InputDimension-1,4);
3632 m_BCRegions[2]=m_BCRegion;
3636 m_SecondRegion.SetSize(InputDimension-1,2);
3637 m_SecondRegion.SetIndex(InputDimension-1,0);
3640 m_BC2Regions.resize(3);
3641 m_BC2Values.resize(3);
3642 m_BCRegion.SetIndex(InputDimension-1,1);
3643 m_BC2Regions[0]=m_BCRegion;
3644 m_BC2Values[0]=-m_WeightRatio[3][1];
3645 m_BCRegion.SetIndex(InputDimension-1,2);
3646 m_BC2Regions[1]=m_BCRegion;
3648 m_BCRegion.SetIndex(InputDimension-1,3);
3649 m_BC2Regions[2]=m_BCRegion;
3650 m_BC2Values[2]=-m_WeightRatio[3][1];
3657 m_FirstRegion.SetSize(InputDimension-1,2);
3658 m_FirstRegion.SetIndex(InputDimension-1,0);
3661 m_BCRegions.resize(3);
3662 m_BCValues.resize(3);
3663 m_BCRegion.SetIndex(InputDimension-1,1);
3664 m_BCRegions[0]=m_BCRegion;
3665 m_BCValues[0]=-m_WeightRatio[3][1];
3666 m_BCRegion.SetIndex(InputDimension-1,2);
3667 m_BCRegions[1]=m_BCRegion;
3669 m_BCRegion.SetIndex(InputDimension-1,3);
3670 m_BCRegions[2]=m_BCRegion;
3671 m_BCValues[2]=-m_WeightRatio[3][1];
3674 m_SecondRegion.SetSize(InputDimension-1,1);
3675 m_SecondRegion.SetIndex(InputDimension-1,2);
3678 m_BC2Regions.resize(0);
3685 m_FirstRegion.SetSize(InputDimension-1,1);
3686 m_FirstRegion.SetIndex(InputDimension-1,1);
3689 m_BCRegions.resize(3);
3690 m_BCValues.resize(3);
3691 m_BCRegion.SetIndex(InputDimension-1,1);
3692 m_BCRegions[0]=m_BCRegion;
3693 m_BCValues[0]=-m_WeightRatio[3][1];
3694 m_BCRegion.SetIndex(InputDimension-1,2);
3695 m_BCRegions[1]=m_BCRegion;
3697 m_BCRegion.SetIndex(InputDimension-1,3);
3698 m_BCRegions[2]=m_BCRegion;
3699 m_BCValues[2]=-m_WeightRatio[3][1];
3702 m_SecondRegion.SetSize(InputDimension-1,2);
3703 m_SecondRegion.SetIndex(InputDimension-1,2);
3706 m_BC2Regions.resize(0);
3713 m_FirstRegion.SetSize(InputDimension-1,0);
3716 m_BCRegions.resize(3);
3717 m_BCValues.resize(3);
3718 m_BCRegion.SetIndex(InputDimension-1,1);
3719 m_BCRegions[0]=m_BCRegion;
3720 m_BCValues[0]=-m_WeightRatio[3][1];
3721 m_BCRegion.SetIndex(InputDimension-1,2);
3722 m_BCRegions[1]=m_BCRegion;
3724 m_BCRegion.SetIndex(InputDimension-1,3);
3725 m_BCRegions[2]=m_BCRegion;
3726 m_BCValues[2]=-m_WeightRatio[3][1];
3729 m_SecondRegion.SetSize(InputDimension-1,2);
3730 m_SecondRegion.SetIndex(InputDimension-1,2);
3733 m_BC2Regions.resize(3);
3734 m_BC2Values.resize(3);
3735 m_BCRegion.SetIndex(InputDimension-1,0);
3736 m_BC2Regions[0]=m_BCRegion;
3738 m_BCRegion.SetIndex(InputDimension-1,1);
3739 m_BC2Regions[1]=m_BCRegion;
3741 m_BCRegion.SetIndex(InputDimension-1,4);
3742 m_BC2Regions[2]=m_BCRegion;
3743 m_BC2Values[2]=-0.5;
3750 m_FirstRegion.SetSize(InputDimension-1,2);
3751 m_FirstRegion.SetIndex(InputDimension-1,2);
3754 m_BCRegions.resize(3);
3755 m_BCValues.resize(3);
3756 m_BCRegion.SetIndex(InputDimension-1,0);
3757 m_BCRegions[0]=m_BCRegion;
3759 m_BCRegion.SetIndex(InputDimension-1,1);
3760 m_BCRegions[1]=m_BCRegion;
3762 m_BCRegion.SetIndex(InputDimension-1,4);
3763 m_BCRegions[2]=m_BCRegion;
3767 m_SecondRegion.SetSize(InputDimension-1,1);
3768 m_SecondRegion.SetIndex(InputDimension-1,4);
3771 m_BC2Regions.resize(0);
3777 DD("supportindex > 4 ???");
3778 DD(m_SupportRegion.GetIndex(InputDimension-1));
3780 } // end switch index
3783 } // end case 5 shape
3790 // ----------------------------------------------------------------------
3791 // The diamond with 4 internal CP:
3792 // Periodic, constrained to zero at the reference at position 3
3793 // Coeff row 0 1 2 BC1 3 BC2 4
3794 // PaddedCoeff R:0 1 2 3 4 5 6
3795 // BC1= -weights[2]/weights[0] ( R2+R4 )
3796 // BC2= weights[2]/weights[0] ( R0+R2-R4-R6 ) + R1
3797 // ---------------------------------------------------------------------
3798 switch(m_SupportRegion.GetIndex(InputDimension-1))
3803 m_FirstRegion.SetSize(InputDimension-1,3);
3804 m_FirstRegion.SetIndex(InputDimension-1,0);
3807 m_BCRegions.resize(2);
3808 m_BCValues.resize(2);
3809 m_BCRegion.SetIndex(InputDimension-1,2);
3810 m_BCRegions[0]=m_BCRegion;
3811 m_BCValues[0]=-m_WeightRatio[2][0];
3812 m_BCRegion.SetIndex(InputDimension-1,3);
3813 m_BCRegions[1]=m_BCRegion;
3814 m_BCValues[1]=-m_WeightRatio[2][0];
3817 m_SecondRegion.SetSize(InputDimension-1,0);
3820 m_BC2Regions.resize(0);
3827 m_FirstRegion.SetSize(InputDimension-1,2);
3828 m_FirstRegion.SetIndex(InputDimension-1,1);
3831 m_BCRegions.resize(2);
3832 m_BCValues.resize(2);
3833 m_BCRegion.SetIndex(InputDimension-1,2);
3834 m_BCRegions[0]=m_BCRegion;
3835 m_BCValues[0]=-m_WeightRatio[2][0];
3836 m_BCRegion.SetIndex(InputDimension-1,3);
3837 m_BCRegions[1]=m_BCRegion;
3838 m_BCValues[1]=-m_WeightRatio[2][0];
3841 m_SecondRegion.SetSize(InputDimension-1,1);
3842 m_SecondRegion.SetIndex(InputDimension-1,3);
3845 m_BC2Regions.resize(0);
3852 m_FirstRegion.SetSize(InputDimension-1,1);
3853 m_FirstRegion.SetIndex(InputDimension-1,2);
3856 m_BCRegions.resize(2);
3857 m_BCValues.resize(2);
3858 m_BCRegion.SetIndex(InputDimension-1,2);
3859 m_BCRegions[0]=m_BCRegion;
3860 m_BCValues[0]=-m_WeightRatio[2][0];
3861 m_BCRegion.SetIndex(InputDimension-1,3);
3862 m_BCRegions[1]=m_BCRegion;
3863 m_BCValues[1]=-m_WeightRatio[2][0];
3866 m_SecondRegion.SetSize(InputDimension-1,1);
3867 m_SecondRegion.SetIndex(InputDimension-1,3);
3870 m_BC2Regions.resize(5);
3871 m_BC2Values.resize(5);
3872 m_BCRegion.SetIndex(InputDimension-1,0);
3873 m_BC2Regions[0]=m_BCRegion;
3874 m_BC2Values[0]=m_WeightRatio[2][0];
3875 m_BCRegion.SetIndex(InputDimension-1,1);
3876 m_BC2Regions[1]=m_BCRegion;
3878 m_BCRegion.SetIndex(InputDimension-1,2);
3879 m_BC2Regions[2]=m_BCRegion;
3880 m_BC2Values[2]=m_WeightRatio[2][0];
3881 m_BCRegion.SetIndex(InputDimension-1,3);
3882 m_BC2Regions[3]=m_BCRegion;
3883 m_BC2Values[3]=-m_WeightRatio[2][0];
3884 m_BCRegion.SetIndex(InputDimension-1,4);
3885 m_BC2Regions[4]=m_BCRegion;
3886 m_BC2Values[4]=-m_WeightRatio[2][0];
3893 m_FirstRegion.SetSize(InputDimension-1,0);
3896 m_BCRegions.resize(2);
3897 m_BCValues.resize(2);
3898 m_BCRegion.SetIndex(InputDimension-1,2);
3899 m_BCRegions[0]=m_BCRegion;
3900 m_BCValues[0]=-m_WeightRatio[2][0];
3901 m_BCRegion.SetIndex(InputDimension-1,3);
3902 m_BCRegions[1]=m_BCRegion;
3903 m_BCValues[1]=-m_WeightRatio[2][0];
3906 m_SecondRegion.SetSize(InputDimension-1,1);
3907 m_SecondRegion.SetIndex(InputDimension-1,3);
3910 m_BC2Regions.resize(5);
3911 m_BC2Values.resize(5);
3912 m_BCRegion.SetIndex(InputDimension-1,0);
3913 m_BC2Regions[0]=m_BCRegion;
3914 m_BC2Values[0]=m_WeightRatio[2][0];
3915 m_BCRegion.SetIndex(InputDimension-1,1);
3916 m_BC2Regions[1]=m_BCRegion;
3918 m_BCRegion.SetIndex(InputDimension-1,2);
3919 m_BC2Regions[2]=m_BCRegion;
3920 m_BC2Values[2]=m_WeightRatio[2][0];
3921 m_BCRegion.SetIndex(InputDimension-1,3);
3922 m_BC2Regions[3]=m_BCRegion;
3923 m_BC2Values[3]=-m_WeightRatio[2][0];
3924 m_BCRegion.SetIndex(InputDimension-1,4);
3925 m_BC2Regions[4]=m_BCRegion;
3926 m_BC2Values[4]=-m_WeightRatio[2][0];
3929 m_ThirdRegion.SetSize(InputDimension-1,1);
3930 m_ThirdRegion.SetIndex(InputDimension-1,4);
3937 DD("supportindex > 3 ???");
3938 DD(m_SupportRegion.GetIndex(InputDimension-1));
3940 } // end switch index
3943 } // end case 7 shape
3947 // ----------------------------------------------------------------------
3948 // The diamond with 5 internal CP:
3949 // periodic, constrained to zero at the reference at position 3.5
3950 // Coeff row 0 1 2 BC1 3 4 BC2 5
3951 // PaddedCoeff R:0 1 2 3 4 5 6 7
3952 // BC1= -weights[3]/weights[1] ( R2+R5 ) - R4
3953 // BC2= weights[2]/weights[0] ( R0+R2-R5-R7 ) + R1
3954 // ---------------------------------------------------------------------
3955 switch(m_SupportRegion.GetIndex(InputDimension-1))
3960 m_FirstRegion.SetSize(InputDimension-1,3);
3961 m_FirstRegion.SetIndex(InputDimension-1,0);
3964 m_BCRegions.resize(3);
3965 m_BCValues.resize(3);
3966 m_BCRegion.SetIndex(InputDimension-1,2);
3967 m_BCRegions[0]=m_BCRegion;
3968 m_BCValues[0]=-m_WeightRatio[3][1];
3969 m_BCRegion.SetIndex(InputDimension-1,3);
3970 m_BCRegions[1]=m_BCRegion;
3972 m_BCRegion.SetIndex(InputDimension-1,4);
3973 m_BCRegions[2]=m_BCRegion;
3974 m_BCValues[2]=-m_WeightRatio[3][1];
3977 m_SecondRegion.SetSize(InputDimension-1,0);
3980 m_BC2Regions.resize(0);
3987 m_FirstRegion.SetSize(InputDimension-1,2);
3988 m_FirstRegion.SetIndex(InputDimension-1,1);
3991 m_BCRegions.resize(3);
3992 m_BCValues.resize(3);
3993 m_BCRegion.SetIndex(InputDimension-1,2);
3994 m_BCRegions[0]=m_BCRegion;
3995 m_BCValues[0]=-m_WeightRatio[3][1];
3996 m_BCRegion.SetIndex(InputDimension-1,3);
3997 m_BCRegions[1]=m_BCRegion;
3999 m_BCRegion.SetIndex(InputDimension-1,4);
4000 m_BCRegions[2]=m_BCRegion;
4001 m_BCValues[2]=-m_WeightRatio[3][1];
4004 m_SecondRegion.SetSize(InputDimension-1,1);
4005 m_SecondRegion.SetIndex(InputDimension-1,3);
4008 m_BC2Regions.resize(0);
4015 m_FirstRegion.SetSize(InputDimension-1,1);
4016 m_FirstRegion.SetIndex(InputDimension-1,2);
4019 m_BCRegions.resize(3);
4020 m_BCValues.resize(3);
4021 m_BCRegion.SetIndex(InputDimension-1,2);
4022 m_BCRegions[0]=m_BCRegion;
4023 m_BCValues[0]=-m_WeightRatio[3][1];
4024 m_BCRegion.SetIndex(InputDimension-1,3);
4025 m_BCRegions[1]=m_BCRegion;
4027 m_BCRegion.SetIndex(InputDimension-1,4);
4028 m_BCRegions[2]=m_BCRegion;
4029 m_BCValues[2]=-m_WeightRatio[3][1];
4032 m_SecondRegion.SetSize(InputDimension-1,2);
4033 m_SecondRegion.SetIndex(InputDimension-1,3);
4036 m_BC2Regions.resize(0);
4043 m_FirstRegion.SetSize(InputDimension-1,0);
4044 m_FirstRegion.SetIndex(InputDimension-1,0);
4047 m_BCRegions.resize(3);
4048 m_BCValues.resize(3);
4049 m_BCRegion.SetIndex(InputDimension-1,2);
4050 m_BCRegions[0]=m_BCRegion;
4051 m_BCValues[0]=-m_WeightRatio[3][1];
4052 m_BCRegion.SetIndex(InputDimension-1,3);
4053 m_BCRegions[1]=m_BCRegion;
4055 m_BCRegion.SetIndex(InputDimension-1,4);
4056 m_BCRegions[2]=m_BCRegion;
4057 m_BCValues[2]=-m_WeightRatio[3][1];
4060 m_SecondRegion.SetSize(InputDimension-1,2);
4061 m_SecondRegion.SetIndex(InputDimension-1,3);
4064 m_BC2Regions.resize(5);
4065 m_BC2Values.resize(5);
4066 m_BCRegion.SetIndex(InputDimension-1,0);
4067 m_BC2Regions[0]=m_BCRegion;
4068 m_BC2Values[0]=m_WeightRatio[2][0];
4069 m_BCRegion.SetIndex(InputDimension-1,1);
4070 m_BC2Regions[1]=m_BCRegion;
4072 m_BCRegion.SetIndex(InputDimension-1,2);
4073 m_BC2Regions[2]=m_BCRegion;
4074 m_BC2Values[2]=m_WeightRatio[2][0];
4075 m_BCRegion.SetIndex(InputDimension-1,4);
4076 m_BC2Regions[3]=m_BCRegion;
4077 m_BC2Values[3]=-m_WeightRatio[2][0];
4078 m_BCRegion.SetIndex(InputDimension-1,5);
4079 m_BC2Regions[4]=m_BCRegion;
4080 m_BC2Values[4]=-m_WeightRatio[2][0];
4088 m_FirstRegion.SetSize(InputDimension-1,2);
4089 m_FirstRegion.SetIndex(InputDimension-1,3);
4092 m_BCRegions.resize(5);
4093 m_BCValues.resize(5);
4094 m_BCRegion.SetIndex(InputDimension-1,0);
4095 m_BCRegions[0]=m_BCRegion;
4096 m_BCValues[0]=m_WeightRatio[2][0];
4097 m_BCRegion.SetIndex(InputDimension-1,1);
4098 m_BCRegions[1]=m_BCRegion;
4100 m_BCRegion.SetIndex(InputDimension-1,2);
4101 m_BCRegions[2]=m_BCRegion;
4102 m_BCValues[2]=m_WeightRatio[2][0];
4103 m_BCRegion.SetIndex(InputDimension-1,4);
4104 m_BCRegions[3]=m_BCRegion;
4105 m_BCValues[3]=-m_WeightRatio[2][0];
4106 m_BCRegion.SetIndex(InputDimension-1,5);
4107 m_BCRegions[4]=m_BCRegion;
4108 m_BCValues[4]=-m_WeightRatio[2][0];
4111 m_SecondRegion.SetSize(InputDimension-1,1);
4112 m_SecondRegion.SetIndex(InputDimension-1,5);
4115 m_BC2Regions.resize(0);
4122 DD("supportindex > 4 ???");
4123 DD(m_SupportRegion.GetIndex(InputDimension-1));
4125 } // end switch index
4128 } // end case 7 shape
4132 // ----------------------------------------------------------------------
4133 // The sputnik with 5 internal CP
4134 // Periodic, constrained to zero at the reference
4135 // at position 3 and one indepent extreme
4136 // Coeff row BC1 0 1 BC2 2 3 BC3 4
4137 // PaddedCoeff R: 0 1 2 3 4 5 6 7
4139 // BC2= -weights[3]/weights[1] ( R2+R5 ) - R4
4141 // ---------------------------------------------------------------------
4142 switch(m_SupportRegion.GetIndex(InputDimension-1))
4147 m_FirstRegion.SetSize(InputDimension-1,0);
4150 m_BCRegions.resize(3);
4151 m_BCValues.resize(3);
4152 m_BCRegion.SetIndex(InputDimension-1,1);
4153 m_BCRegions[0]=m_BCRegion;
4155 m_BCRegion.SetIndex(InputDimension-1,3);
4156 m_BCRegions[1]=m_BCRegion;
4158 m_BCRegion.SetIndex(InputDimension-1,4);
4159 m_BCRegions[2]=m_BCRegion;
4163 m_SecondRegion.SetSize(InputDimension-1,2);
4164 m_SecondRegion.SetIndex(InputDimension-1,0);
4167 m_BC2Regions.resize(3);
4168 m_BC2Values.resize(3);
4169 m_BCRegion.SetIndex(InputDimension-1,1);
4170 m_BC2Regions[0]=m_BCRegion;
4171 m_BC2Values[0]=-m_WeightRatio[3][1];
4172 m_BCRegion.SetIndex(InputDimension-1,2);
4173 m_BC2Regions[1]=m_BCRegion;
4175 m_BCRegion.SetIndex(InputDimension-1,3);
4176 m_BC2Regions[2]=m_BCRegion;
4177 m_BC2Values[2]=-m_WeightRatio[3][1];
4184 m_FirstRegion.SetSize(InputDimension-1,2);
4185 m_FirstRegion.SetIndex(InputDimension-1,0);
4188 m_BCRegions.resize(3);
4189 m_BCValues.resize(3);
4190 m_BCRegion.SetIndex(InputDimension-1,1);
4191 m_BCRegions[0]=m_BCRegion;
4192 m_BCValues[0]=-m_WeightRatio[3][1];
4193 m_BCRegion.SetIndex(InputDimension-1,2);
4194 m_BCRegions[1]=m_BCRegion;
4196 m_BCRegion.SetIndex(InputDimension-1,3);
4197 m_BCRegions[2]=m_BCRegion;
4198 m_BCValues[2]=-m_WeightRatio[3][1];
4201 m_SecondRegion.SetSize(InputDimension-1,1);
4202 m_SecondRegion.SetIndex(InputDimension-1,2);
4205 m_BC2Regions.resize(0);
4212 m_FirstRegion.SetSize(InputDimension-1,1);
4213 m_FirstRegion.SetIndex(InputDimension-1,1);
4216 m_BCRegions.resize(3);
4217 m_BCValues.resize(3);
4218 m_BCRegion.SetIndex(InputDimension-1,1);
4219 m_BCRegions[0]=m_BCRegion;
4220 m_BCValues[0]=-m_WeightRatio[3][1];
4221 m_BCRegion.SetIndex(InputDimension-1,2);
4222 m_BCRegions[1]=m_BCRegion;
4224 m_BCRegion.SetIndex(InputDimension-1,3);
4225 m_BCRegions[2]=m_BCRegion;
4226 m_BCValues[2]=-m_WeightRatio[3][1];
4229 m_SecondRegion.SetSize(InputDimension-1,2);
4230 m_SecondRegion.SetIndex(InputDimension-1,2);
4233 m_BC2Regions.resize(0);
4240 m_FirstRegion.SetSize(InputDimension-1,0);
4243 m_BCRegions.resize(3);
4244 m_BCValues.resize(3);
4245 m_BCRegion.SetIndex(InputDimension-1,1);
4246 m_BCRegions[0]=m_BCRegion;
4247 m_BCValues[0]=-m_WeightRatio[3][1];
4248 m_BCRegion.SetIndex(InputDimension-1,2);
4249 m_BCRegions[1]=m_BCRegion;
4251 m_BCRegion.SetIndex(InputDimension-1,3);
4252 m_BCRegions[2]=m_BCRegion;
4253 m_BCValues[2]=-m_WeightRatio[3][1];
4256 m_SecondRegion.SetSize(InputDimension-1,1);
4257 m_SecondRegion.SetIndex(InputDimension-1,2);
4260 m_BC2Regions.resize(1);
4261 m_BC2Values.resize(1);
4262 m_BCRegion.SetIndex(InputDimension-1,0);
4263 m_BC2Regions[0]=m_BCRegion;
4271 m_FirstRegion.SetSize(InputDimension-1,2);
4272 m_FirstRegion.SetIndex(InputDimension-1,2);
4275 m_BCRegions.resize(1);
4276 m_BCValues.resize(1);
4277 m_BCRegion.SetIndex(InputDimension-1,0);
4278 m_BCRegions[0]=m_BCRegion;
4282 m_SecondRegion.SetSize(InputDimension-1,1);
4283 m_SecondRegion.SetIndex(InputDimension-1,4);
4286 m_BC2Regions.resize(0);
4292 DD("supportindex > 4 ???");
4293 DD(m_SupportRegion.GetIndex(InputDimension-1));
4295 } // end switch index
4298 } // end case 9 shape
4303 DD ("Other shapes currently not implemented");
4306 } // end switch shape
4307 } // end wrap region
4311 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
4313 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions,NOutputDimensions>
4314 ::TransformPointToContinuousIndex( const InputPointType & point, ContinuousIndexType & index ) const
4318 itk::Vector<double, OutputDimension> tvector;
4320 for ( j = 0; j < OutputDimension; j++ )
4322 //JV find the index in the PADDED version
4323 tvector[j] = point[j] - this->m_PaddedGridOrigin[j];
4326 itk::Vector<double, OutputDimension> cvector;
4328 cvector = m_PointToIndex * tvector;
4330 for ( j = 0; j < OutputDimension; j++ )
4332 index[j] = static_cast< typename ContinuousIndexType::CoordRepType >( cvector[j] );