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://www.centreleonberard.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;
1042 typename CoefficientImageType::RegionType sourceRegion=m_PaddedCoefficientImage->GetLargestPossibleRegion();
1043 typename CoefficientImageType::RegionType destinationRegion=m_PaddedCoefficientImage->GetLargestPossibleRegion();
1044 typename CoefficientImageType::RegionType::SizeType sourceSize=sourceRegion.GetSize();
1045 typename CoefficientImageType::RegionType::SizeType destinationSize=destinationRegion.GetSize();
1046 typename CoefficientImageType::IndexType sourceIndex=sourceRegion.GetIndex();
1047 typename CoefficientImageType::IndexType destinationIndex=destinationRegion.GetIndex();
1049 // JV Padding depends on the shape
1050 switch (m_TransformShape)
1055 2: rabbit 4 CP 3 DOF
1056 3: rabbit 5 CP 4 DOF
1057 4: sputnik 4 CP 4 DOF
1058 5: sputnik 5 CP 5 DOF
1059 6: diamond 6 CP 5 DOF
1060 7: diamond 7 CP 6 DOF
1065 // ----------------------------------------------------------------------
1066 // The egg with 4 internal CP (starting from inhale)
1067 // Periodic, constrained to zero at the reference
1068 // at position 3 and
1069 // Coeff row BC1 BC2 0 BC3 1 2 BC4
1070 // PaddedCoeff R: 0 1 2 3 4 5 6
1073 // BC3= -weights[2]/weights[0] ( R2+R4 )
1075 // ---------------------------------------------------------------------
1077 //---------------------------------
1078 // 1. First two temporal rows are identical: paste 1-2 to 0-1
1079 typename PasteImageFilterType::Pointer paster0=PasteImageFilterType::New();
1080 paster0->SetDestinationImage(m_PaddedCoefficientImage);
1081 paster0->SetSourceImage(m_CoefficientImage);
1082 sourceRegion.SetIndex(NInputDimensions-1,1);
1083 sourceRegion.SetSize(NInputDimensions-1,2);
1084 destinationIndex[InputDimension-1]=0;
1085 paster0->SetDestinationIndex(destinationIndex);
1086 paster0->SetSourceRegion(sourceRegion);
1088 //---------------------------------
1089 // 1. Next temporal row = paste 0 to 2
1090 typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1091 typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1092 paster1->SetDestinationImage(paster0->GetOutput());
1093 paster1->SetSourceImage(row0);
1094 destinationIndex[InputDimension-1]=2;
1095 paster1->SetDestinationIndex(destinationIndex);
1096 paster1->SetSourceRegion(row0->GetLargestPossibleRegion());
1098 //---------------------------------
1099 // 2. Middle row at index 3=BC
1100 // Extract row at index 1, 2 (of coeff image)
1101 typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1104 typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1105 combine1->SetInput(0,row0);
1106 combine1->SetInput(1,row1);
1107 combine1->SetA(-m_WeightRatio[2][0]);
1108 combine1->SetB(-m_WeightRatio[2][0]);
1110 typename CoefficientImageType::Pointer bc3Row=combine1->GetOutput();
1112 // Paste middleRow at index 3 (padded coeff)
1113 typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1114 paster2->SetDestinationImage(paster1->GetOutput());
1115 paster2->SetSourceImage(bc3Row);
1116 destinationIndex[InputDimension-1]=3;
1117 paster2->SetDestinationIndex(destinationIndex);
1118 paster2->SetSourceRegion(bc3Row->GetLargestPossibleRegion());
1120 //---------------------------------
1121 // 3. Next two temporal rows identical: paste 1,2 to 4,5
1122 typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1123 paster3->SetDestinationImage(paster2->GetOutput());
1124 paster3->SetSourceImage(m_CoefficientImage);
1125 sourceRegion.SetIndex(InputDimension-1,1);
1126 sourceRegion.SetSize(InputDimension-1,2);
1127 destinationIndex[InputDimension-1]=4;
1128 paster3->SetDestinationIndex(destinationIndex);
1129 paster3->SetSourceRegion(sourceRegion);
1131 // ---------------------------------
1132 // 4. Row at index 6=BC4= R2
1133 // Paste BC3 row at index 5
1134 typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
1135 paster4->SetDestinationImage(paster3->GetOutput());
1136 paster4->SetSourceImage(row0);
1137 destinationIndex[InputDimension-1]=6;
1138 paster4->SetDestinationIndex(destinationIndex);
1139 paster4->SetSourceRegion(row0->GetLargestPossibleRegion());
1141 // Update the chain!
1143 m_PaddedCoefficientImage= paster4->GetOutput();
1150 // ----------------------------------------------------------------------
1151 // The egg with 5 internal CP (starting from inhale)
1152 // Periodic, constrained to zero at the reference
1153 // at position 3 and
1154 // Coeff row BC1 BC2 0 BC3 1 2 3 BC4
1155 // PaddedCoeff R: 0 1 2 3 4 5 6 7
1158 // BC3= -weights[3]/weights[1] ( R2+R5 ) - R4
1160 // ---------------------------------------------------------------------
1161 //---------------------------------
1162 // 1. First two temporal rows are identical: paste 2-3 to 0-1
1163 typename PasteImageFilterType::Pointer paster0=PasteImageFilterType::New();
1164 paster0->SetDestinationImage(m_PaddedCoefficientImage);
1165 paster0->SetSourceImage(m_CoefficientImage);
1166 sourceRegion.SetIndex(NInputDimensions-1,2);
1167 sourceRegion.SetSize(NInputDimensions-1,2);
1168 destinationIndex[InputDimension-1]=0;
1169 paster0->SetDestinationIndex(destinationIndex);
1170 paster0->SetSourceRegion(sourceRegion);
1172 //---------------------------------
1173 // 1. Next temporal row = paste 0 to 2
1174 typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1175 typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1176 paster1->SetDestinationImage(paster0->GetOutput());
1177 paster1->SetSourceImage(row0);
1178 destinationIndex[InputDimension-1]=2;
1179 paster1->SetDestinationIndex(destinationIndex);
1180 paster1->SetSourceRegion(row0->GetLargestPossibleRegion());
1182 //---------------------------------
1183 // 2. Middle row at index 3=BC
1184 // Extract row at index 1, 2 (of coeff image)
1185 typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1186 typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
1189 typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1190 combine1->SetInput(0,row0);
1191 combine1->SetInput(1,row2);
1194 // combine1->Update();
1195 typename LinearCombinationFilterType::Pointer combine2=LinearCombinationFilterType::New();
1196 combine2->SetInput(0,row1);
1197 combine2->SetInput(1,combine1->GetOutput());
1198 combine2->SetA(-1.);
1199 combine2->SetB(-m_WeightRatio[3][1]);
1201 typename CoefficientImageType::Pointer bc3Row=combine2->GetOutput();
1203 // Paste middleRow at index 3 (padded coeff)
1204 typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1205 paster2->SetDestinationImage(paster1->GetOutput());
1206 paster2->SetSourceImage(bc3Row);
1207 destinationIndex[InputDimension-1]=3;
1208 paster2->SetDestinationIndex(destinationIndex);
1209 paster2->SetSourceRegion(bc3Row->GetLargestPossibleRegion());
1211 //---------------------------------
1212 // 3. Next three temporal rows identical: paste 1,2,3 to 4,5,6
1213 typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1214 paster3->SetDestinationImage(paster2->GetOutput());
1215 paster3->SetSourceImage(m_CoefficientImage);
1216 sourceRegion.SetIndex(InputDimension-1,1);
1217 sourceRegion.SetSize(InputDimension-1,3);
1218 destinationIndex[InputDimension-1]=4;
1219 paster3->SetDestinationIndex(destinationIndex);
1220 paster3->SetSourceRegion(sourceRegion);
1222 // ---------------------------------
1223 // 4. Row at index 7=BC4= R2
1224 // Paste BC3 row at index 5
1225 typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
1226 paster4->SetDestinationImage(paster3->GetOutput());
1227 paster4->SetSourceImage(row0);
1228 destinationIndex[InputDimension-1]=7;
1229 paster4->SetDestinationIndex(destinationIndex);
1230 paster4->SetSourceRegion(row0->GetLargestPossibleRegion());
1232 // Update the chain!
1234 m_PaddedCoefficientImage= paster4->GetOutput();
1240 // // ----------------------------------------------------------------------
1241 // // The egg with 5 internal CP:
1242 // // Periodic, C2 smooth everywhere and constrained to zero at the reference
1243 // // Coeff row R5 BC1 0 1 2 3 BC2 R2
1244 // // PaddedCoeff R: 0 1 2 3 4 5 6 7
1245 // // BC1= -weights[2]/weights[0] ( R2+R5)
1247 // // ---------------------------------------------------------------------
1249 // // Extract rows with index 0 and 3
1250 // typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1251 // typename CoefficientImageType::Pointer row3=this->ExtractTemporalRow(m_CoefficientImage, 3);
1253 // // Paste the first row
1254 // typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1255 // destinationIndex.Fill(0);
1256 // paster1->SetDestinationIndex(destinationIndex);
1257 // paster1->SetSourceRegion(row3->GetLargestPossibleRegion());
1258 // paster1->SetSourceImage(row3);
1259 // paster1->SetDestinationImage(m_PaddedCoefficientImage);
1261 // // Linearly Combine rows for BC1 and BC2
1262 // typedef clitk::LinearCombinationImageFilter<CoefficientImageType, CoefficientImageType> LinearCombinationFilterType;
1263 // typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1264 // combine1->SetFirstInput(row0);
1265 // combine1->SetSecondInput(row3);
1266 // combine1->SetA(-m_WeightRatio[2][0]);
1267 // combine1->SetB(-m_WeightRatio[2][0]);
1268 // combine1->Update();
1269 // typename CoefficientImageType::Pointer bcRow=combine1->GetOutput();
1271 // // Paste the second row
1272 // typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1273 // destinationIndex[InputDimension-1]=1;
1274 // paster2->SetDestinationIndex(destinationIndex);
1275 // paster2->SetSourceRegion(bcRow->GetLargestPossibleRegion());
1276 // paster2->SetSourceImage(bcRow);
1277 // paster2->SetDestinationImage(paster1->GetOutput());
1279 // // Paste the coefficientImage
1280 // typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1281 // destinationIndex[InputDimension-1]=2;
1282 // paster3->SetDestinationIndex(destinationIndex);
1283 // paster3->SetSourceRegion(m_CoefficientImage->GetLargestPossibleRegion());
1284 // paster3->SetSourceImage(m_CoefficientImage);
1285 // paster3->SetDestinationImage(paster2->GetOutput());
1287 // // Paste the last two rows
1288 // typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
1289 // destinationIndex[InputDimension-1]=6;
1290 // paster4->SetDestinationIndex(destinationIndex);
1291 // paster4->SetSourceRegion(bcRow->GetLargestPossibleRegion());
1292 // paster4->SetSourceImage(bcRow);
1293 // paster4->SetDestinationImage(paster3->GetOutput());
1295 // typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
1296 // destinationIndex[InputDimension-1]=7;
1297 // paster5->SetDestinationIndex(destinationIndex);
1298 // paster5->SetSourceRegion(row0->GetLargestPossibleRegion());
1299 // paster5->SetSourceImage(row0);
1300 // paster5->SetDestinationImage(paster4->GetOutput());
1302 // // Update the chain!
1303 // paster5->Update();
1304 // m_PaddedCoefficientImage= paster5->GetOutput();
1311 // ----------------------------------------------------------------------
1312 // The rabbit with 4 internal CP
1313 // Periodic, constrained to zero at the reference
1314 // at position 3 and the extremes fixed with anit-symm bc
1315 // Coeff row BC1 0 1 BC2 2 BC3 BC4
1316 // PaddedCoeff R: 0 1 2 3 4 5 6
1318 // BC2= -weights[2]/weights[0] ( R2+R4 )
1321 // ---------------------------------------------------------------------
1323 // ---------------------------------
1324 // 0. First Row =BC1
1325 typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1326 typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1327 typename LinearCombinationFilterType::Pointer combine0=LinearCombinationFilterType::New();
1328 combine0->SetInput(0,row0);
1329 combine0->SetInput(1,row1);
1331 combine0->SetB(-1.);
1333 typename CoefficientImageType::Pointer bc1Row=combine0->GetOutput();
1335 typename PasteImageFilterType::Pointer paster0=PasteImageFilterType::New();
1336 paster0->SetDestinationImage(m_PaddedCoefficientImage);
1337 paster0->SetSourceImage(bc1Row);
1338 destinationIndex[InputDimension-1]=0;
1339 paster0->SetDestinationIndex(destinationIndex);
1340 paster0->SetSourceRegion(bc1Row->GetLargestPossibleRegion());
1342 //---------------------------------
1343 // 1. Next two temporal rows are identical: paste 0-1 to 1-2
1344 typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1345 paster1->SetDestinationImage(paster0->GetOutput());
1346 paster1->SetSourceImage(m_CoefficientImage);
1347 sourceRegion.SetIndex(NInputDimensions-1,0);
1348 destinationIndex[InputDimension-1]=1;
1349 sourceRegion.SetSize(NInputDimensions-1,2);
1350 paster1->SetDestinationIndex(destinationIndex);
1351 paster1->SetSourceRegion(sourceRegion);
1353 //---------------------------------
1354 // 2. Middle row at index 3=BC
1355 // Extract row at index 1, 2 (of coeff image)
1356 typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
1359 typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1360 combine1->SetInput(0,row1);
1361 combine1->SetInput(1,row2);
1362 combine1->SetA(-m_WeightRatio[2][0]);
1363 combine1->SetB(-m_WeightRatio[2][0]);
1365 typename CoefficientImageType::Pointer bc2Row=combine1->GetOutput();
1367 // Paste middleRow at index 3 (padded coeff)
1368 typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1369 paster2->SetDestinationImage(paster1->GetOutput());
1370 paster2->SetSourceImage(bc2Row);
1371 destinationIndex[InputDimension-1]=3;
1372 paster2->SetDestinationIndex(destinationIndex);
1373 paster2->SetSourceRegion(bc2Row->GetLargestPossibleRegion());
1375 //---------------------------------
1376 // 3. Next temporal row is identical: paste 2 to 4
1377 typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1378 paster3->SetDestinationImage(paster2->GetOutput());
1379 paster3->SetSourceImage(row2);
1380 destinationIndex[InputDimension-1]=4;
1381 paster3->SetDestinationIndex(destinationIndex);
1382 paster3->SetSourceRegion(row2->GetLargestPossibleRegion());
1384 // ---------------------------------
1385 // 4. Row at index 5=BC (paddedcoeff image) R1
1386 // Paste BC3 row at index 5
1387 typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
1388 paster4->SetDestinationImage(paster3->GetOutput());
1389 paster4->SetSourceImage(row0);
1390 destinationIndex[InputDimension-1]=5;
1391 paster4->SetDestinationIndex(destinationIndex);
1392 paster4->SetSourceRegion(row0->GetLargestPossibleRegion());
1394 // ---------------------------------
1395 // 5. Paste BC4 row at index 6
1396 typename LinearCombinationFilterType::Pointer combine3=LinearCombinationFilterType::New();
1397 combine3->SetInput(0,row0);
1398 combine3->SetInput(1,row2);
1400 combine3->SetB(-1.);
1402 typename CoefficientImageType::Pointer bc4Row=combine3->GetOutput();
1403 typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
1404 paster5->SetDestinationImage(paster4->GetOutput());
1405 paster5->SetSourceImage(bc4Row);
1406 destinationIndex[InputDimension-1]=6;
1407 paster5->SetDestinationIndex(destinationIndex);
1408 paster5->SetSourceRegion(bc4Row->GetLargestPossibleRegion());
1410 // Update the chain!
1412 m_PaddedCoefficientImage= paster5->GetOutput();
1419 // ----------------------------------------------------------------------
1420 // The rabbit with 5 internal CP
1421 // Periodic, constrained to zero at the reference at position 3.5
1422 // and the extremes fixed with anti-symmetrical BC
1423 // Coeff row BC1 0 1 BC2 2 3 BC3 BC4
1424 // PaddedCoeff R: 0 1 2 3 4 5 6 7
1426 // BC2= -weights[3]/weights[1] ( R2+R5 ) - R4
1429 // ---------------------------------------------------------------------
1431 // ---------------------------------
1432 // 0. First Row =BC1
1433 typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1434 typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1435 typename LinearCombinationFilterType::Pointer combine0=LinearCombinationFilterType::New();
1436 combine0->SetInput(0,row0);
1437 combine0->SetInput(1,row1);
1439 combine0->SetB(-1.);
1441 typename CoefficientImageType::Pointer bc1Row=combine0->GetOutput();
1443 typename PasteImageFilterType::Pointer paster0=PasteImageFilterType::New();
1444 paster0->SetDestinationImage(m_PaddedCoefficientImage);
1445 paster0->SetSourceImage(bc1Row);
1446 destinationIndex[InputDimension-1]=0;
1447 paster0->SetDestinationIndex(destinationIndex);
1448 paster0->SetSourceRegion(bc1Row->GetLargestPossibleRegion());
1450 //---------------------------------
1451 // 1. Next two temporal rows are identical: paste 0-1 to 1-2
1452 typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1453 paster1->SetDestinationImage(paster0->GetOutput());
1454 paster1->SetSourceImage(m_CoefficientImage);
1455 sourceRegion.SetIndex(InputDimension-1,0);
1456 destinationIndex[InputDimension-1]=1;
1457 sourceRegion.SetSize(NInputDimensions-1,2);
1458 paster1->SetDestinationIndex(destinationIndex);
1459 paster1->SetSourceRegion(sourceRegion);
1461 //---------------------------------
1462 // 2. Middle row at index 3=BC
1463 // Extract row at index 1, 3 (of coeff image)
1464 //typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1465 typename CoefficientImageType::Pointer row3=this->ExtractTemporalRow(m_CoefficientImage, 3);
1468 typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1469 combine1->SetInput(0,row1);
1470 combine1->SetInput(1,row3);
1475 // Extract row at index 2 (coeff image)
1476 typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
1479 typename LinearCombinationFilterType::Pointer combine2=LinearCombinationFilterType::New();
1480 combine2->SetInput(0,combine1->GetOutput());
1481 combine2->SetInput(1,row2);
1482 combine2->SetA(-m_WeightRatio[3][1]);
1483 combine2->SetB(-1.);
1485 typename CoefficientImageType::Pointer bc2Row=combine2->GetOutput();
1487 // Paste middleRow at index 3 (padded coeff)
1488 typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1489 paster2->SetDestinationImage(paster1->GetOutput());
1490 paster2->SetSourceImage(bc2Row);
1491 destinationIndex[InputDimension-1]=3;
1492 paster2->SetDestinationIndex(destinationIndex);
1493 paster2->SetSourceRegion(bc2Row->GetLargestPossibleRegion());
1495 //---------------------------------
1496 // 3. Next two temporal rows are identical: paste 2,3 to 4,5
1497 typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1498 paster3->SetDestinationImage(paster2->GetOutput());
1499 paster3->SetSourceImage(m_CoefficientImage);
1500 sourceRegion.SetIndex(InputDimension-1,2);
1501 destinationIndex[InputDimension-1]=4;
1502 sourceRegion.SetSize(NInputDimensions-1,2);
1503 paster3->SetDestinationIndex(destinationIndex);
1504 paster3->SetSourceRegion(sourceRegion);
1506 // ---------------------------------
1507 // 4. Row at index 6=BC (paddedcoeff image)R1
1508 // Paste BC3 row at index 6
1509 typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
1510 paster4->SetDestinationImage(paster3->GetOutput());
1511 paster4->SetSourceImage(row0);
1512 destinationIndex[InputDimension-1]=6;
1513 paster4->SetDestinationIndex(destinationIndex);
1514 paster4->SetSourceRegion(row0->GetLargestPossibleRegion());
1516 // ---------------------------------
1517 // 5. Paste BC4 row at index 7
1518 typename LinearCombinationFilterType::Pointer combine3=LinearCombinationFilterType::New();
1519 combine3->SetInput(0,row0);
1520 combine3->SetInput(1,row3);
1522 combine3->SetB(-1.);
1524 typename CoefficientImageType::Pointer bc4Row=combine3->GetOutput();
1525 typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
1526 paster5->SetDestinationImage(paster4->GetOutput());
1527 paster5->SetSourceImage(bc4Row);
1528 destinationIndex[InputDimension-1]=7;
1529 paster5->SetDestinationIndex(destinationIndex);
1530 paster5->SetSourceRegion(bc4Row->GetLargestPossibleRegion());
1532 // Update the chain!
1534 m_PaddedCoefficientImage= paster5->GetOutput();
1541 // ----------------------------------------------------------------------
1542 // The sputnik with 4 internal CP
1543 // Periodic, constrained to zero at the reference
1544 // at position 3 and one indepent extremes copied
1545 // Coeff row BC1 0 1 BC2 2 BC2 3
1546 // PaddedCoeff R: 0 1 2 3 4 5 6
1548 // BC2= -weights[2]/weights[0] ( R2+R4 )
1549 // BC3= weights[2]/weights[0] ( R2-R4) + R1
1550 // ---------------------------------------------------------------------
1552 //---------------------------------
1553 // 1. First Row is equal to last row: paste 3 row to 0
1554 typename CoefficientImageType::Pointer row3=this->ExtractTemporalRow(m_CoefficientImage, 3);
1555 typename PasteImageFilterType::Pointer paster0=PasteImageFilterType::New();
1556 paster0->SetDestinationImage(m_PaddedCoefficientImage);
1557 paster0->SetSourceImage(row3);
1558 destinationIndex[InputDimension-1]=0;
1559 paster0->SetDestinationIndex(destinationIndex);
1560 paster0->SetSourceRegion(row3->GetLargestPossibleRegion());
1562 //---------------------------------
1563 // 1. Next two temporal rows are identical: paste 0-1 to 1-2
1564 typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1565 paster1->SetDestinationImage(paster0->GetOutput());
1566 paster1->SetSourceImage(m_CoefficientImage);
1567 sourceRegion.SetIndex(NInputDimensions-1,0);
1568 destinationIndex[InputDimension-1]=1;
1569 sourceRegion.SetSize(NInputDimensions-1,2);
1570 paster1->SetDestinationIndex(destinationIndex);
1571 paster1->SetSourceRegion(sourceRegion);
1573 //---------------------------------
1574 // 2. Middle row at index 3=BC
1575 // Extract row at index 1, 2 (of coeff image)
1576 typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1577 typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
1580 typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1581 combine1->SetInput(0,row1);
1582 combine1->SetInput(1,row2);
1583 combine1->SetA(-m_WeightRatio[2][0]);
1584 combine1->SetB(-m_WeightRatio[2][0]);
1586 typename CoefficientImageType::Pointer bc1Row=combine1->GetOutput();
1588 // Paste middleRow at index 3 (padded coeff)
1589 typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1590 paster2->SetDestinationImage(paster1->GetOutput());
1591 paster2->SetSourceImage(bc1Row);
1592 destinationIndex[InputDimension-1]=3;
1593 paster2->SetDestinationIndex(destinationIndex);
1594 paster2->SetSourceRegion(bc1Row->GetLargestPossibleRegion());
1596 //---------------------------------
1597 // 3. Next temporal row is identical: paste 2 to 4
1598 typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1599 paster3->SetDestinationImage(paster2->GetOutput());
1600 paster3->SetSourceImage(row2);
1601 destinationIndex[InputDimension-1]=4;
1602 paster3->SetDestinationIndex(destinationIndex);
1603 paster3->SetSourceRegion(row2->GetLargestPossibleRegion());
1605 //---------------------------------
1606 // 4. Final row at index 5=BC3 (paddedcoeff image)R1 R2 R4 corresponds to index in coeff image 0 1 2
1607 // Extract row at index 1, 2 (of coeff image)already done
1608 typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1611 typename LinearCombinationFilterType::Pointer combine3=LinearCombinationFilterType::New();
1612 combine3->SetInput(0,row1);
1613 combine3->SetInput(1,row2);
1615 combine3->SetB(-1.);
1618 typename LinearCombinationFilterType::Pointer combine4=LinearCombinationFilterType::New();
1619 combine4->SetInput(0,combine3->GetOutput());
1620 combine4->SetInput(1,row0);
1621 combine4->SetA(m_WeightRatio[2][0]);
1624 typename CoefficientImageType::Pointer bc2Row=combine4->GetOutput();
1626 // Paste BC row at index 5
1627 typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
1628 paster4->SetDestinationImage(paster3->GetOutput());
1629 paster4->SetSourceImage(bc2Row);
1630 destinationIndex[InputDimension-1]=5;
1631 paster4->SetDestinationIndex(destinationIndex);
1632 paster4->SetSourceRegion(bc2Row->GetLargestPossibleRegion());
1634 // Paste row 3 at index 6
1635 typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
1636 paster5->SetDestinationImage(paster4->GetOutput());
1637 paster5->SetSourceImage(row3);
1638 destinationIndex[InputDimension-1]=6;
1639 paster5->SetDestinationIndex(destinationIndex);
1640 paster5->SetSourceRegion(row3->GetLargestPossibleRegion());
1642 // Update the chain!
1644 m_PaddedCoefficientImage= paster5->GetOutput();
1652 // ----------------------------------------------------------------------
1653 // The sputnik with 5 internal CP
1654 // Periodic, constrained to zero at the reference
1655 // at position 3 and one indepent extreme
1656 // Coeff row BC1 0 1 BC2 2 3 BC3 4
1657 // PaddedCoeff R: 0 1 2 3 4 5 6 7
1658 // BC1= R2 + R5 - R7
1659 // BC2= -weights[3]/weights[1] ( R2+R5 ) - R4
1660 // BC3= R1 + 0.5 R2 - 0.5 R7
1661 // ----------------------------------------------------------------------
1662 //---------------------------------
1664 typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1665 typename CoefficientImageType::Pointer row3=this->ExtractTemporalRow(m_CoefficientImage, 3);
1666 typename CoefficientImageType::Pointer row4=this->ExtractTemporalRow(m_CoefficientImage, 4);
1669 typename LinearCombinationFilterType::Pointer combine0=LinearCombinationFilterType::New();
1670 combine0->SetInput(0,row1);
1671 combine0->SetInput(1,row3);
1674 //combine0->Update();
1675 typename LinearCombinationFilterType::Pointer combine0bis=LinearCombinationFilterType::New();
1676 combine0bis->SetInput(0,combine0->GetOutput());
1677 combine0bis->SetInput(1,row4);
1678 combine0bis->SetA(1.);
1679 combine0bis->SetB(-1.);
1680 combine0bis->Update();
1681 typename CoefficientImageType::Pointer bc1Row=combine0bis->GetOutput();
1683 typename PasteImageFilterType::Pointer paster0=PasteImageFilterType::New();
1684 paster0->SetDestinationImage(m_PaddedCoefficientImage);
1685 paster0->SetSourceImage(bc1Row);
1686 destinationIndex[InputDimension-1]=0;
1687 paster0->SetDestinationIndex(destinationIndex);
1688 paster0->SetSourceRegion(bc1Row->GetLargestPossibleRegion());
1690 //---------------------------------
1691 // 1. Next two temporal rows are identical: paste 0-1 to 1-2
1692 typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1693 paster1->SetDestinationImage(paster0->GetOutput());
1694 paster1->SetSourceImage(m_CoefficientImage);
1695 sourceRegion.SetIndex(NInputDimensions-1,0);
1696 destinationIndex[InputDimension-1]=1;
1697 sourceRegion.SetSize(NInputDimensions-1,2);
1698 paster1->SetDestinationIndex(destinationIndex);
1699 paster1->SetSourceRegion(sourceRegion);
1701 //---------------------------------
1702 // 2. Middle row at index 3=BC
1703 // Extract row at index 1, 2 (of coeff image)
1704 //typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1705 typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
1708 typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1709 combine1->SetInput(0,row1);
1710 combine1->SetInput(1,row3);
1715 typename LinearCombinationFilterType::Pointer combine2=LinearCombinationFilterType::New();
1716 combine2->SetInput(0,combine1->GetOutput());
1717 combine2->SetInput(1,row2);
1718 combine2->SetA(-m_WeightRatio[3][1]);
1719 combine2->SetB(-1.);
1721 typename CoefficientImageType::Pointer bc2Row=combine2->GetOutput();
1723 // Paste middleRow at index 3 (padded coeff)
1724 typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1725 paster2->SetDestinationImage(paster1->GetOutput());
1726 paster2->SetSourceImage(bc2Row);
1727 destinationIndex[InputDimension-1]=3;
1728 paster2->SetDestinationIndex(destinationIndex);
1729 paster2->SetSourceRegion(bc2Row->GetLargestPossibleRegion());
1731 //---------------------------------
1732 // 3. Next two temporal rows are identical: paste 2,3 to 4,5
1733 typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1734 paster3->SetDestinationImage(paster2->GetOutput());
1735 paster3->SetSourceImage(m_CoefficientImage);
1736 sourceRegion.SetIndex(InputDimension-1,2);
1737 destinationIndex[InputDimension-1]=4;
1738 sourceRegion.SetSize(NInputDimensions-1,2);
1739 paster3->SetDestinationIndex(destinationIndex);
1740 paster3->SetSourceRegion(sourceRegion);
1742 //---------------------------------
1743 // 4. Final row at index 6=BC3
1744 typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1747 typename LinearCombinationFilterType::Pointer combine3=LinearCombinationFilterType::New();
1748 combine3->SetInput(0,row1);
1749 combine3->SetInput(1,row4);
1751 combine3->SetB(-1.);
1752 typename LinearCombinationFilterType::Pointer combine4=LinearCombinationFilterType::New();
1753 combine4->SetInput(0,row0);
1754 combine4->SetInput(1,combine3->GetOutput());
1758 typename CoefficientImageType::Pointer bc3Row=combine4->GetOutput();
1760 // Paste BC row at index 6
1761 typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
1762 paster4->SetDestinationImage(paster3->GetOutput());
1763 paster4->SetSourceImage(bc3Row);
1764 destinationIndex[InputDimension-1]=6;
1765 paster4->SetDestinationIndex(destinationIndex);
1766 paster4->SetSourceRegion(bc3Row->GetLargestPossibleRegion());
1768 // Paste row 4 at index 7
1769 typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
1770 paster5->SetDestinationImage(paster4->GetOutput());
1771 paster5->SetSourceImage(row4);
1772 destinationIndex[InputDimension-1]=7;
1773 paster5->SetDestinationIndex(destinationIndex);
1774 paster5->SetSourceRegion(row4->GetLargestPossibleRegion());
1776 // Update the chain!
1778 m_PaddedCoefficientImage= paster5->GetOutput();
1785 // ----------------------------------------------------------------------
1786 // The diamond with 4 internal CP:
1787 // Periodic, constrained to zero at the reference at position 3
1788 // Coeff row 0 1 2 BC1 3 BC2 4
1789 // PaddedCoeff R:0 1 2 3 4 5 6
1790 // BC1= -weights[2]/weights[0] ( R2+R4 )
1791 // BC2= weights[2]/weights[0] ( R0+R2-R4-R6 ) + R1
1792 // ---------------------------------------------------------------------
1794 //---------------------------------
1795 // 1. First Three temporal rows are identical: paste 0-2 to 0-2
1796 typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1797 paster1->SetDestinationImage(m_PaddedCoefficientImage);
1798 paster1->SetSourceImage(m_CoefficientImage);
1799 sourceIndex.Fill(0);
1800 destinationIndex.Fill(0);
1801 sourceSize[NInputDimensions-1]=3;
1802 sourceRegion.SetSize(sourceSize);
1803 sourceRegion.SetIndex(sourceIndex);
1804 paster1->SetDestinationIndex(destinationIndex);
1805 paster1->SetSourceRegion(sourceRegion);
1807 //---------------------------------
1808 // 2. Middle row at index 3=BC
1809 // Extract row at index 0, 4 (of coeff image)
1810 typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
1811 typename CoefficientImageType::Pointer row3=this->ExtractTemporalRow(m_CoefficientImage, 3);
1814 typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1815 combine1->SetInput(0,row2);
1816 combine1->SetInput(1,row3);
1817 combine1->SetA(-m_WeightRatio[2][0]);
1818 combine1->SetB(-m_WeightRatio[2][0]);
1820 typename CoefficientImageType::Pointer bc1Row=combine1->GetOutput();
1822 // Paste middleRow at index 3 (padded coeff)
1823 typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1824 paster2->SetDestinationImage(paster1->GetOutput());
1825 paster2->SetSourceImage(bc1Row);
1826 destinationIndex[InputDimension-1]=3;
1827 paster2->SetDestinationIndex(destinationIndex);
1828 paster2->SetSourceRegion(bc1Row->GetLargestPossibleRegion());
1830 //---------------------------------
1831 // 3. Next row identical: paste 3 to 4
1832 typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1833 paster3->SetDestinationImage(paster2->GetOutput());
1834 paster3->SetSourceImage(row3);
1835 destinationIndex[InputDimension-1]=4;
1836 paster3->SetDestinationIndex(destinationIndex);
1837 paster3->SetSourceRegion(row3->GetLargestPossibleRegion());
1839 //---------------------------------
1840 // 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
1841 // Extract row at index 0, 2 (of coeff image)already done
1842 typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1843 typename CoefficientImageType::Pointer row4=this->ExtractTemporalRow(m_CoefficientImage, 4);
1844 typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1847 typename LinearCombinationFilterType::Pointer combine3=LinearCombinationFilterType::New();
1848 combine3->SetInput(0,row0);
1849 combine3->SetInput(1,row2);
1854 typename LinearCombinationFilterType::Pointer combine4=LinearCombinationFilterType::New();
1855 combine4->SetInput(0,row3);
1856 combine4->SetInput(1,row4);
1861 typename LinearCombinationFilterType::Pointer combine5=LinearCombinationFilterType::New();
1862 combine5->SetInput(0,combine3->GetOutput());
1863 combine5->SetInput(1,combine4->GetOutput());
1865 combine5->SetB(-1.);
1868 typename LinearCombinationFilterType::Pointer combine6=LinearCombinationFilterType::New();
1869 combine6->SetInput(0,combine5->GetOutput());
1870 combine6->SetInput(1,row1);
1871 combine6->SetA(m_WeightRatio[2][0]);
1874 typename CoefficientImageType::Pointer bc2Row=combine6->GetOutput();
1876 // Paste BC row at index 5
1877 typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
1878 paster4->SetDestinationImage(paster3->GetOutput());
1879 paster4->SetSourceImage(bc2Row);
1880 destinationIndex[InputDimension-1]=5;
1881 paster4->SetDestinationIndex(destinationIndex);
1882 paster4->SetSourceRegion(bc2Row->GetLargestPossibleRegion());
1884 // Paste last row at index 6
1885 typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
1886 paster5->SetDestinationImage(paster4->GetOutput());
1887 paster5->SetSourceImage(row4);
1888 destinationIndex[InputDimension-1]=6;
1889 paster5->SetDestinationIndex(destinationIndex);
1890 paster5->SetSourceRegion(row4->GetLargestPossibleRegion());
1892 // Update the chain!
1894 m_PaddedCoefficientImage= paster5->GetOutput();
1900 // ----------------------------------------------------------------------
1901 // The diamond with 5 internal CP:
1902 // periodic, constrained to zero at the reference at position 3.5
1903 // Coeff row 0 1 2 BC1 3 4 BC2 5
1904 // PaddedCoeff R:0 1 2 3 4 5 6 7
1905 // BC1= -weights[3]/weights[1] ( R2+R5 ) - R4
1906 // BC2= weights[2]/weights[0] ( R0+R2-R5-R7 ) + R1
1907 // ---------------------------------------------------------------------
1909 //---------------------------------
1910 // 1. First Three temporal rows are identical: paste 0-2 to 0-2
1911 typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1912 paster1->SetDestinationImage(m_PaddedCoefficientImage);
1913 paster1->SetSourceImage(m_CoefficientImage);
1914 sourceIndex.Fill(0);
1915 destinationIndex.Fill(0);
1916 sourceSize[NInputDimensions-1]=3;
1917 sourceRegion.SetSize(sourceSize);
1918 sourceRegion.SetIndex(sourceIndex);
1919 paster1->SetDestinationIndex(destinationIndex);
1920 paster1->SetSourceRegion(sourceRegion);
1922 //---------------------------------
1923 // 2. Middle row at index 3=BC
1924 // Extract row at index 0, 4 (of coeff image)
1925 typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
1926 typename CoefficientImageType::Pointer row4=this->ExtractTemporalRow(m_CoefficientImage, 4);
1929 typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1930 combine1->SetInput(0,row2);
1931 combine1->SetInput(1,row4);
1936 // Extract row at index 3 (coeff image)
1937 typename CoefficientImageType::Pointer row3=this->ExtractTemporalRow(m_CoefficientImage, 3);
1940 typename LinearCombinationFilterType::Pointer combine2=LinearCombinationFilterType::New();
1941 combine2->SetInput(0,combine1->GetOutput());
1942 combine2->SetInput(1,row3);
1943 combine2->SetA(-m_WeightRatio[3][1] );
1944 combine2->SetB(-1.);
1946 typename CoefficientImageType::Pointer bc1Row=combine2->GetOutput();
1948 // Paste middleRow at index 3 (padded coeff)
1949 typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1950 paster2->SetDestinationImage(paster1->GetOutput());
1951 paster2->SetSourceImage(bc1Row);
1952 destinationIndex[InputDimension-1]=3;
1953 paster2->SetDestinationIndex(destinationIndex);
1954 paster2->SetSourceRegion(bc1Row->GetLargestPossibleRegion());
1956 //---------------------------------
1957 // 3. Next two temporal rows are identical: paste 3,4 to 4,5
1958 typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1959 paster3->SetDestinationImage(paster2->GetOutput());
1960 paster3->SetSourceImage(m_CoefficientImage);
1961 sourceIndex[InputDimension-1]=3;
1962 destinationIndex[InputDimension-1]=4;
1963 sourceSize[NInputDimensions-1]=2;
1964 sourceRegion.SetSize(sourceSize);
1965 sourceRegion.SetIndex(sourceIndex);
1966 paster3->SetDestinationIndex(destinationIndex);
1967 paster3->SetSourceRegion(sourceRegion);
1969 //---------------------------------
1970 // 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
1971 // Extract row at index 0, 2 (of coeff image)already done
1972 typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1973 typename CoefficientImageType::Pointer row5=this->ExtractTemporalRow(m_CoefficientImage, 5);
1974 typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1977 typename LinearCombinationFilterType::Pointer combine3=LinearCombinationFilterType::New();
1978 combine3->SetInput(0,row0);
1979 combine3->SetInput(1,row2);
1984 typename LinearCombinationFilterType::Pointer combine4=LinearCombinationFilterType::New();
1985 combine4->SetInput(0,row4);
1986 combine4->SetInput(1,row5);
1991 typename LinearCombinationFilterType::Pointer combine5=LinearCombinationFilterType::New();
1992 combine5->SetInput(0,combine3->GetOutput());
1993 combine5->SetInput(1,combine4->GetOutput());
1995 combine5->SetB(-1.);
1998 typename LinearCombinationFilterType::Pointer combine6=LinearCombinationFilterType::New();
1999 combine6->SetInput(0,combine5->GetOutput());
2000 combine6->SetInput(1,row1);
2001 combine6->SetA(m_WeightRatio[2][0]);
2004 typename CoefficientImageType::Pointer bc2Row=combine6->GetOutput();
2006 // Paste BC row at index 6
2007 typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
2008 paster4->SetDestinationImage(paster3->GetOutput());
2009 paster4->SetSourceImage(bc2Row);
2010 destinationIndex[InputDimension-1]=6;
2011 paster4->SetDestinationIndex(destinationIndex);
2012 paster4->SetSourceRegion(bc2Row->GetLargestPossibleRegion());
2014 // Paste last row at index 7
2015 typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
2016 paster5->SetDestinationImage(paster4->GetOutput());
2017 paster5->SetSourceImage(row5);
2018 destinationIndex[InputDimension-1]=7;
2019 paster5->SetDestinationIndex(destinationIndex);
2020 paster5->SetSourceRegion(row5->GetLargestPossibleRegion());
2022 // Update the chain!
2024 m_PaddedCoefficientImage= paster5->GetOutput();
2031 // ----------------------------------------------------------------------
2032 // The sputnik with 5 internal CP T''(0)=T''(10)
2033 // Periodic, constrained to zero at the reference
2034 // at position 3 and one indepent extreme
2035 // Coeff row BC1 0 1 BC2 2 3 BC3 4
2036 // PaddedCoeff R: 0 1 2 3 4 5 6 7
2038 // BC2= -weights[3]/weights[1] ( R2+R5 ) - R4
2040 // ---------------------------------------------------------------------
2042 //---------------------------------
2044 typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
2045 typename CoefficientImageType::Pointer row3=this->ExtractTemporalRow(m_CoefficientImage, 3);
2046 typename CoefficientImageType::Pointer row4=this->ExtractTemporalRow(m_CoefficientImage, 4);
2049 typename LinearCombinationFilterType::Pointer combine0=LinearCombinationFilterType::New();
2050 combine0->SetInput(0,row3);
2051 combine0->SetInput(1,row4);
2054 typename LinearCombinationFilterType::Pointer combine0bis=LinearCombinationFilterType::New();
2055 combine0bis->SetInput(0,combine0->GetOutput());
2056 combine0bis->SetInput(1,row1);
2057 combine0bis->SetA(1.);
2058 combine0bis->SetB(-1.);
2059 combine0bis->Update();
2060 typename CoefficientImageType::Pointer bc1Row=combine0bis->GetOutput();
2062 typename PasteImageFilterType::Pointer paster0=PasteImageFilterType::New();
2063 paster0->SetDestinationImage(m_PaddedCoefficientImage);
2064 paster0->SetSourceImage(bc1Row);
2065 destinationIndex[InputDimension-1]=0;
2066 paster0->SetDestinationIndex(destinationIndex);
2067 paster0->SetSourceRegion(bc1Row->GetLargestPossibleRegion());
2069 //---------------------------------
2070 // 1. Next two temporal rows are identical: paste 0-1 to 1-2
2071 typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
2072 paster1->SetDestinationImage(paster0->GetOutput());
2073 paster1->SetSourceImage(m_CoefficientImage);
2074 sourceRegion.SetIndex(NInputDimensions-1,0);
2075 destinationIndex[InputDimension-1]=1;
2076 sourceRegion.SetSize(NInputDimensions-1,2);
2077 paster1->SetDestinationIndex(destinationIndex);
2078 paster1->SetSourceRegion(sourceRegion);
2080 //---------------------------------
2081 // 2. Middle row at index 3=BC
2082 // Extract row at index 1, 2 (of coeff image)
2083 // typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
2084 typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
2087 typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
2088 combine1->SetInput(0,row1);
2089 combine1->SetInput(1,row3);
2094 typename LinearCombinationFilterType::Pointer combine2=LinearCombinationFilterType::New();
2095 combine2->SetInput(0,combine1->GetOutput());
2096 combine2->SetInput(1,row2);
2097 combine2->SetA(-m_WeightRatio[3][1]);
2098 combine2->SetB(-1.);
2100 typename CoefficientImageType::Pointer bc2Row=combine2->GetOutput();
2102 // Paste middleRow at index 3 (padded coeff)
2103 typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
2104 paster2->SetDestinationImage(paster1->GetOutput());
2105 paster2->SetSourceImage(bc2Row);
2106 destinationIndex[InputDimension-1]=3;
2107 paster2->SetDestinationIndex(destinationIndex);
2108 paster2->SetSourceRegion(bc2Row->GetLargestPossibleRegion());
2110 //---------------------------------
2111 // 3. Next two temporal rows are identical: paste 2,3 to 4,5
2112 typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
2113 paster3->SetDestinationImage(paster2->GetOutput());
2114 paster3->SetSourceImage(m_CoefficientImage);
2115 sourceRegion.SetIndex(InputDimension-1,2);
2116 destinationIndex[InputDimension-1]=4;
2117 sourceRegion.SetSize(NInputDimensions-1,2);
2118 paster3->SetDestinationIndex(destinationIndex);
2119 paster3->SetSourceRegion(sourceRegion);
2121 //---------------------------------
2122 // 4. Final row at index 6=BC3
2123 typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
2126 // typename LinearCombinationFilterType::Pointer combine3=LinearCombinationFilterType::New();
2127 // combine3->SetInput(0,row0);
2128 // combine3->SetInput(1,row1);
2129 // combine3->SetA(1.);
2130 // combine3->SetB(0.5);
2131 // combine3->Update();
2132 // typename CoefficientImageType::Pointer bc3Row=combine3->GetOutput();
2134 // Paste BC row at index 6
2135 typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
2136 paster4->SetDestinationImage(paster3->GetOutput());
2137 paster4->SetSourceImage(row0);
2138 destinationIndex[InputDimension-1]=6;
2139 paster4->SetDestinationIndex(destinationIndex);
2140 paster4->SetSourceRegion(row0->GetLargestPossibleRegion());
2142 // Paste row 4 at index 7
2143 typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
2144 paster5->SetDestinationImage(paster4->GetOutput());
2145 paster5->SetSourceImage(row4);
2146 destinationIndex[InputDimension-1]=7;
2147 paster5->SetDestinationIndex(destinationIndex);
2148 paster5->SetSourceRegion(row4->GetLargestPossibleRegion());
2150 // Update the chain!
2152 m_PaddedCoefficientImage= paster5->GetOutput();
2158 DD ("Shape not available");
2164 // // Extract coefficients from padded version
2165 // template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
2167 // ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2168 // ::ExtractCoefficientImage( )
2170 // ////DD("extract coeff image");
2171 // typename ExtractImageFilterType::Pointer extractImageFilter = ExtractImageFilterType::New();
2172 // typename CoefficientImageType::RegionType extractionRegion=m_PaddedCoefficientImage->GetLargestPossibleRegion();
2173 // typename CoefficientImageType::RegionType::SizeType extractionSize=extractionRegion.GetSize();
2174 // typename CoefficientImageType::RegionType::IndexType extractionIndex = extractionRegion.GetIndex();
2175 // extractionSize[InputDimension-1]-=4;
2176 // extractionIndex[InputDimension-1]=2;
2177 // extractionRegion.SetSize(extractionSize);
2178 // extractionRegion.SetIndex(extractionIndex);
2179 // extractImageFilter->SetInput(m_PaddedCoefficientImage);
2180 // extractImageFilter->SetExtractionRegion(extractionRegion);
2181 // extractImageFilter->Update();
2182 // m_CoefficientImage=extractImageFilter->GetOutput();
2187 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
2189 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2190 ::PrintSelf(std::ostream &os, itk::Indent indent) const
2193 this->Superclass::PrintSelf(os, indent);
2195 os << indent << "GridRegion: " << m_GridRegion << std::endl;
2196 os << indent << "GridOrigin: " << m_GridOrigin << std::endl;
2197 os << indent << "GridSpacing: " << m_GridSpacing << std::endl;
2198 os << indent << "GridDirection: " << m_GridDirection << std::endl;
2199 os << indent << "IndexToPoint: " << m_IndexToPoint << std::endl;
2200 os << indent << "PointToIndex: " << m_PointToIndex << std::endl;
2202 os << indent << "CoefficientImage: [ ";
2203 os << m_CoefficientImage.GetPointer() << " ]" << std::endl;
2205 os << indent << "WrappedImage: [ ";
2206 os << m_WrappedImage.GetPointer() << " ]" << std::endl;
2208 os << indent << "InputParametersPointer: "
2209 << m_InputParametersPointer << std::endl;
2210 os << indent << "ValidRegion: " << m_ValidRegion << std::endl;
2211 os << indent << "LastJacobianIndex: " << m_LastJacobianIndex << std::endl;
2212 os << indent << "BulkTransform: ";
2213 os << m_BulkTransform << std::endl;
2215 if ( m_BulkTransform )
2217 os << indent << "BulkTransformType: "
2218 << m_BulkTransform->GetNameOfClass() << std::endl;
2220 os << indent << "VectorBSplineInterpolator: ";
2221 os << m_VectorInterpolator.GetPointer() << std::endl;
2222 os << indent << "Mask: ";
2223 os << m_Mask<< std::endl;
2228 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
2230 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2231 ::InsideValidRegion( const ContinuousIndexType& index ) const
2235 if ( !m_ValidRegion.IsInside( index ) )
2240 // JV verify all dimensions
2243 typedef typename ContinuousIndexType::ValueType ValueType;
2244 for( unsigned int j = 0; j < NInputDimensions; j++ )
2246 if (m_SplineOrderOdd[j])
2248 if ( index[j] >= static_cast<ValueType>( m_ValidRegionLast[j] ) )
2260 // Transform a point
2261 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
2262 typename ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2264 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2265 ::TransformPoint(const InputPointType &inputPoint) const
2268 InputPointType transformedPoint;
2269 OutputPointType outputPoint;
2272 if ( m_BulkTransform )
2274 transformedPoint = m_BulkTransform->TransformPoint( inputPoint );
2278 transformedPoint = inputPoint;
2281 // Deformable transform
2282 if ( m_PaddedCoefficientImage )
2284 // Check if inside mask
2285 if(m_Mask && !(m_Mask->IsInside(inputPoint) ) )
2287 // Outside: no (deformable) displacement
2288 return transformedPoint;
2291 // Check if inside valid region
2293 ContinuousIndexType index;
2294 this->TransformPointToContinuousIndex( inputPoint, index );
2295 inside = this->InsideValidRegion( index );
2298 // Outside: no (deformable) displacement
2299 DD("outside valid region");
2300 return transformedPoint;
2303 // Call the vector interpolator
2304 itk::Vector<TCoordRep,SpaceDimension> displacement=m_VectorInterpolator->EvaluateAtContinuousIndex(index);
2306 // JV add for the spatial dimensions
2307 outputPoint=transformedPoint;
2308 for (unsigned int i=0; i<NInputDimensions-1; i++)
2309 outputPoint[i] += displacement[i];
2315 itkWarningMacro( << "B-spline coefficients have not been set" );
2316 outputPoint = transformedPoint;
2324 //JV Deformably transform a point
2325 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
2326 typename ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2328 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2329 ::DeformablyTransformPoint(const InputPointType &inputPoint) const
2331 OutputPointType outputPoint;
2332 if ( m_PaddedCoefficientImage )
2335 // Check if inside mask
2336 if(m_Mask && !(m_Mask->IsInside(inputPoint) ) )
2338 // Outside: no (deformable) displacement
2344 ContinuousIndexType index;
2345 this->TransformPointToContinuousIndex( inputPoint, index );
2346 inside = this->InsideValidRegion( index );
2350 //outside: no (deformable) displacement
2351 outputPoint = inputPoint;
2355 // Call the vector interpolator
2356 itk::Vector<TCoordRep,SpaceDimension> displacement=m_VectorInterpolator->EvaluateAtContinuousIndex(index);
2358 // JV add for the spatial dimensions
2359 outputPoint=inputPoint;
2360 for (unsigned int i=0; i<NInputDimensions-1; i++)
2361 outputPoint[i] += displacement[i];
2364 // No coefficients available
2367 itkWarningMacro( << "B-spline coefficients have not been set" );
2368 outputPoint = inputPoint;
2375 // JV weights are identical as for transformpoint, could be done simultaneously in metric!!!!
2376 // Compute the Jacobian in one position
2377 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
2379 typename ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2381 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2382 ::GetJacobian( const InputPointType & point ) const
2385 //========================================================
2386 // Zero all components of jacobian
2387 //========================================================
2388 // JV not thread safe (m_LastJacobianIndex), instantiate N transforms
2389 // NOTE: for efficiency, we only need to zero out the coefficients
2390 // that got fill last time this method was called.
2391 unsigned int j=0,b=0;
2393 //Set the previously-set to zero
2394 for ( j = 0; j < SpaceDimension; j++ )
2396 m_FirstIterator[j].GoToBegin();
2397 while ( !m_FirstIterator[j].IsAtEnd() )
2399 m_FirstIterator[j].Set( m_ZeroVector );
2400 ++(m_FirstIterator[j]);
2404 //Set the previously-set to zero
2405 for ( j = 0; j < SpaceDimension; j++ )
2407 m_SecondIterator[j].GoToBegin();
2408 while ( ! (m_SecondIterator[j]).IsAtEnd() )
2410 m_SecondIterator[j].Set( m_ZeroVector );
2411 ++(m_SecondIterator[j]);
2415 //Set the previously-set to zero
2417 for ( j = 0; j < SpaceDimension; j++ )
2419 m_ThirdIterator[j].GoToBegin();
2420 while ( ! (m_ThirdIterator[j]).IsAtEnd() )
2422 m_ThirdIterator[j].Set( m_ZeroVector );
2423 ++(m_ThirdIterator[j]);
2427 //Set the previously-set to zero
2429 for (b=0; b<m_BCSize;b++)
2430 if ( !( m_FirstRegion.IsInside(m_BCRegions[b]) | m_SecondRegion.IsInside(m_BCRegions[b]) ) )
2431 for ( j = 0; j < SpaceDimension; j++ )
2433 m_BCIterators[j][b].GoToBegin();
2434 while ( ! (m_BCIterators[j][b]).IsAtEnd() )
2436 m_BCIterators[j][b].Set( m_ZeroVector );
2437 ++(m_BCIterators[j][b]);
2441 //Set the previously-set to zero
2443 for (b=0; b<m_BC2Size;b++)
2444 if ( !( m_FirstRegion.IsInside(m_BC2Regions[b]) | m_SecondRegion.IsInside(m_BC2Regions[b]) ) )
2445 for ( j = 0; j < SpaceDimension; j++ )
2447 m_BC2Iterators[j][b].GoToBegin();
2448 while ( ! (m_BC2Iterators[j][b]).IsAtEnd() )
2450 m_BC2Iterators[j][b].Set( m_ZeroVector );
2451 ++(m_BC2Iterators[j][b]);
2455 //Set the previously-set to zero
2457 for (b=0; b<m_BC3Size;b++)
2458 if ( !( m_FirstRegion.IsInside(m_BC3Regions[b]) | m_SecondRegion.IsInside(m_BC3Regions[b]) ) )
2459 for ( j = 0; j < SpaceDimension; j++ )
2461 m_BC3Iterators[j][b].GoToBegin();
2462 while ( ! (m_BC3Iterators[j][b]).IsAtEnd() )
2464 m_BC3Iterators[j][b].Set( m_ZeroVector );
2465 ++(m_BC3Iterators[j][b]);
2470 //========================================================
2471 // For each dimension, copy the weight to the support region
2472 //========================================================
2474 // Check if inside mask
2475 if(m_Mask && !(m_Mask->IsInside(point) ) )
2477 // Outside: no (deformable) displacement
2478 return this->m_Jacobian;
2482 this->TransformPointToContinuousIndex( point, m_Index );
2484 // NOTE: if the support region does not lie totally within the grid
2485 // we assume zero displacement and return the input point
2486 if ( !this->InsideValidRegion( m_Index ) )
2488 return this->m_Jacobian;
2491 // Compute interpolation weights
2492 const WeightsDataType *weights=NULL;
2493 m_VectorInterpolator->EvaluateWeightsAtContinuousIndex( m_Index, &weights, m_LastJacobianIndex);
2496 m_SupportRegion.SetIndex( m_LastJacobianIndex );
2497 WrapRegion(m_SupportRegion, m_FirstRegion, m_SecondRegion, m_ThirdRegion, m_BCRegions, m_BCValues, m_BC2Regions, m_BC2Values, m_BC3Regions, m_BC3Values, m_InitialOffset);
2498 m_ThirdSize=m_ThirdRegion.GetSize()[InputDimension-1];
2499 m_BCSize=m_BCRegions.size();
2500 m_BC2Size=m_BC2Regions.size();
2501 m_BC3Size=m_BC3Regions.size();
2503 // Reset the iterators
2504 for ( j = 0; j < SpaceDimension ; j++ )
2506 m_FirstIterator[j] = IteratorType( m_JacobianImage[j], m_FirstRegion);
2507 m_SecondIterator[j] = IteratorType( m_JacobianImage[j], m_SecondRegion);
2508 if(m_ThirdSize) m_ThirdIterator[j] = IteratorType( m_JacobianImage[j], m_ThirdRegion);
2510 m_BCIterators[j].resize(m_BCSize);
2511 for (b=0; b<m_BCSize;b++)
2512 m_BCIterators[j][b]= IteratorType( m_JacobianImage[j], m_BCRegions[b]);
2513 m_BC2Iterators[j].resize(m_BC2Size);
2514 for (b=0; b<m_BC2Size;b++)
2515 m_BC2Iterators[j][b]= IteratorType( m_JacobianImage[j], m_BC2Regions[b]);
2516 m_BC3Iterators[j].resize(m_BC3Size);
2517 for (b=0; b<m_BC3Size;b++)
2518 m_BC3Iterators[j][b]= IteratorType( m_JacobianImage[j], m_BC3Regions[b]);
2521 // Skip if on a fixed condition
2524 if (m_BCSize) weights+=m_InitialOffset*m_BCRegions[0].GetNumberOfPixels();
2525 else std::cerr<<"InitialOffset without BCSize: Error!!!!!!!"<<std::endl;
2528 //copy weight to jacobian image
2529 for ( j = 0; j < SpaceDimension; j++ )
2531 // For each dimension, copy the weight to the support region
2532 while ( ! (m_FirstIterator[j]).IsAtEnd() )
2534 m_ZeroVector[j]=*weights;
2535 (m_FirstIterator[j]).Set( m_ZeroVector);
2536 ++(m_FirstIterator[j]);
2540 m_ZeroVector[j]=itk::NumericTraits<JacobianValueType>::Zero;
2541 if (j != SpaceDimension-1) weights-=m_FirstRegion.GetNumberOfPixels();
2544 // Skip BC1 and go to the second region
2545 if (m_BCSize) weights+=m_BCRegions[0].GetNumberOfPixels();
2547 // For each dimension, copy the weight to the support region
2548 //copy weight to jacobian image
2549 for ( j = 0; j < SpaceDimension; j++ )
2551 while ( ! (m_SecondIterator[j]).IsAtEnd() )
2553 m_ZeroVector[j]=*weights;
2554 (m_SecondIterator[j]).Set( m_ZeroVector);
2555 ++(m_SecondIterator[j]);
2558 weights-=m_SecondRegion.GetNumberOfPixels();
2559 m_ZeroVector[j]=itk::NumericTraits<JacobianValueType>::Zero;
2565 // Put pointer in correct position
2566 weights-=m_BCRegions[0].GetNumberOfPixels();
2568 for ( j = 0; j < SpaceDimension; j++ )
2570 for ( b=0; b < m_BCSize; b++ )
2572 while ( ! (m_BCIterators[j][b]).IsAtEnd() )
2574 //copy weight to jacobian image
2575 m_ZeroVector[j]=(*weights) * m_BCValues[b];
2576 (m_BCIterators[j][b]).Value()+= m_ZeroVector;
2577 ++(m_BCIterators[j][b]);
2580 weights-=m_BCRegions[b].GetNumberOfPixels();
2582 m_ZeroVector[j]=itk::NumericTraits<JacobianValueType>::Zero;
2584 weights+=m_BCRegions[m_BCSize-1].GetNumberOfPixels();
2587 // Add the BC2 to the weights
2590 // Move further in the weights pointer
2591 weights+=m_SecondRegion.GetNumberOfPixels();
2593 for ( j = 0; j < SpaceDimension; j++ )
2595 for ( b=0; b < m_BC2Size; b++ )
2597 while ( ! (m_BC2Iterators[j][b]).IsAtEnd() )
2599 //copy weight to jacobian image
2600 m_ZeroVector[j]=(*weights) * m_BC2Values[b];
2601 (m_BC2Iterators[j][b]).Value()+= m_ZeroVector;
2602 ++(m_BC2Iterators[j][b]);
2605 weights-=m_BC2Regions[b].GetNumberOfPixels();
2607 m_ZeroVector[j]=itk::NumericTraits<JacobianValueType>::Zero;
2609 // Move further in the weights pointer
2610 weights+=m_BC2Regions[m_BC2Size-1].GetNumberOfPixels();
2616 // For each dimension, copy the weight to the support region
2617 //copy weight to jacobian image
2618 for ( j = 0; j < SpaceDimension; j++ )
2620 while ( ! (m_ThirdIterator[j]).IsAtEnd() )
2622 m_ZeroVector[j]=*weights;
2623 (m_ThirdIterator[j]).Value()+= m_ZeroVector;
2624 ++(m_ThirdIterator[j]);
2628 // Move further in the weights pointer?
2629 if (j != SpaceDimension-1) weights-=m_ThirdRegion.GetNumberOfPixels();
2630 m_ZeroVector[j]=itk::NumericTraits<JacobianValueType>::Zero;
2634 // Add the BC3 to the weights
2637 for ( j = 0; j < SpaceDimension; j++ )
2639 for ( b=0; b < m_BC3Size; b++ )
2641 while ( ! (m_BC3Iterators[j][b]).IsAtEnd() )
2643 //copy weight to jacobian image
2644 m_ZeroVector[j]=(*weights) * m_BC3Values[b];
2645 (m_BC3Iterators[j][b]).Value()+= m_ZeroVector;
2646 ++(m_BC3Iterators[j][b]);
2649 weights-=m_BC3Regions[b].GetNumberOfPixels();
2651 m_ZeroVector[j]=itk::NumericTraits<JacobianValueType>::Zero;
2655 // Return the result
2656 return this->m_Jacobian;
2660 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
2662 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2663 ::WrapRegion( const RegionType & m_SupportRegion,
2664 RegionType & m_FirstRegion,
2665 RegionType & m_SecondRegion,
2666 RegionType & m_ThirdRegion,
2667 std::vector<RegionType>& m_BCRegions,std::vector<double>& m_BCValues,
2668 std::vector<RegionType>& m_BC2Regions,std::vector<double>& m_BC2Values,
2669 std::vector<RegionType>& m_BC3Regions,std::vector<double>& m_BC3Values,
2670 unsigned int& m_InitialOffset ) const
2673 // Synchronize regions
2675 m_FirstRegion=m_SupportRegion;
2676 m_BCRegion=m_SupportRegion;
2677 m_BCRegion.SetSize(InputDimension-1,1);
2678 m_SecondRegion=m_SupportRegion;
2679 m_ThirdRegion=m_SupportRegion;
2680 m_ThirdRegion.SetSize(InputDimension-1,0);
2681 m_BC3Regions.resize(0);
2684 // BC depends on shape
2685 switch(m_TransformShape)
2690 2: rabbit 4 CP 3 DOF
2691 3: rabbit 5 CP 4 DOF
2692 4: sputnik 4 CP 4 DOF
2693 5: sputnik 5 CP 5 DOF
2694 6: diamond 6 CP 5 DOF
2695 7: diamond 7 CP 6 DOF
2700 // ----------------------------------------------------------------------
2701 // The egg with 4 internal CP (starting from inhale)
2702 // Periodic, constrained to zero at the reference
2703 // at position 3 and
2704 // Coeff row BC1 BC2 0 BC3 1 2 BC4
2705 // PaddedCoeff R: 0 1 2 3 4 5 6
2708 // BC3= -weights[2]/weights[0] ( R2+R4 )
2710 // ---------------------------------------------------------------------
2711 switch(m_SupportRegion.GetIndex(InputDimension-1))
2716 m_FirstRegion.SetSize(InputDimension-1,2);
2717 m_FirstRegion.SetIndex(InputDimension-1,1);
2720 m_BCRegions.resize(0);
2723 m_SecondRegion.SetSize(InputDimension-1,1);
2724 m_SecondRegion.SetIndex(InputDimension-1,0);
2727 m_BC2Regions.resize(2);
2728 m_BC2Values.resize(2);
2729 m_BCRegion.SetIndex(InputDimension-1,0);
2730 m_BC2Regions[0]=m_BCRegion;
2731 m_BC2Values[0]=-m_WeightRatio[2][0];
2732 m_BCRegion.SetIndex(InputDimension-1,1);
2733 m_BC2Regions[1]=m_BCRegion;
2734 m_BC2Values[1]= -m_WeightRatio[2][0];
2741 m_FirstRegion.SetSize(InputDimension-1,1);
2742 m_FirstRegion.SetIndex(InputDimension-1,2);
2745 m_BCRegions.resize(0);
2748 m_SecondRegion.SetSize(InputDimension-1,1);
2749 m_SecondRegion.SetIndex(InputDimension-1,0);
2752 m_BC2Regions.resize(2);
2753 m_BC2Values.resize(2);
2754 m_BCRegion.SetIndex(InputDimension-1,0);
2755 m_BC2Regions[0]=m_BCRegion;
2756 m_BC2Values[0]=-m_WeightRatio[2][0];
2757 m_BCRegion.SetIndex(InputDimension-1,1);
2758 m_BC2Regions[1]=m_BCRegion;
2759 m_BC2Values[1]= -m_WeightRatio[2][0];
2762 m_FirstRegion.SetSize(InputDimension-1,1);
2763 m_FirstRegion.SetIndex(InputDimension-1,1);
2770 m_FirstRegion.SetSize(InputDimension-1,1);
2771 m_FirstRegion.SetIndex(InputDimension-1,0);
2774 m_BCRegions.resize(2);
2775 m_BCValues.resize(2);
2776 m_BCRegion.SetIndex(InputDimension-1,0);
2777 m_BCRegions[0]=m_BCRegion;
2778 m_BCValues[0]=-m_WeightRatio[2][0];
2779 m_BCRegion.SetIndex(InputDimension-1,1);
2780 m_BCRegions[1]=m_BCRegion;
2781 m_BCValues[1]= -m_WeightRatio[2][0];
2784 m_SecondRegion.SetSize(InputDimension-1,2);
2785 m_SecondRegion.SetIndex(InputDimension-1,1);
2788 m_BC2Regions.resize(0);
2795 m_FirstRegion.SetSize(InputDimension-1,0);
2798 m_BCRegions.resize(2);
2799 m_BCValues.resize(2);
2800 m_BCRegion.SetIndex(InputDimension-1,0);
2801 m_BCRegions[0]=m_BCRegion;
2802 m_BCValues[0]=-m_WeightRatio[2][0];
2803 m_BCRegion.SetIndex(InputDimension-1,1);
2804 m_BCRegions[1]=m_BCRegion;
2805 m_BCValues[1]= -m_WeightRatio[2][0];
2808 m_SecondRegion.SetSize(InputDimension-1,2);
2809 m_SecondRegion.SetIndex(InputDimension-1,1);
2812 m_BC2Regions.resize(1);
2813 m_BC2Values.resize(1);
2814 m_BCRegion.SetIndex(InputDimension-1,0);
2815 m_BC2Regions[0]=m_BCRegion;
2823 DD("supportindex > 3 ???");
2824 DD(m_SupportRegion.GetIndex(InputDimension-1));
2826 } // end switch index
2833 // ----------------------------------------------------------------------
2834 // The egg with 5 internal CP (starting from inhale)
2835 // Periodic, constrained to zero at the reference
2836 // at position 3 and
2837 // Coeff row BC1 BC2 0 BC3 1 2 3 BC4
2838 // PaddedCoeff R: 0 1 2 3 4 5 6 7
2841 // BC3= -weights[3]/weights[1] ( R2+R5 ) - R4
2843 // ---------------------------------------------------------------------
2844 switch(m_SupportRegion.GetIndex(InputDimension-1))
2849 m_FirstRegion.SetSize(InputDimension-1,2);
2850 m_FirstRegion.SetIndex(InputDimension-1,2);
2853 m_BCRegions.resize(0);
2856 m_SecondRegion.SetSize(InputDimension-1,1);
2857 m_SecondRegion.SetIndex(InputDimension-1,0);
2860 m_BC2Regions.resize(3);
2861 m_BC2Values.resize(3);
2862 m_BCRegion.SetIndex(InputDimension-1,0);
2863 m_BC2Regions[0]=m_BCRegion;
2864 m_BC2Values[0]=-m_WeightRatio[3][1];
2865 m_BCRegion.SetIndex(InputDimension-1,1);
2866 m_BC2Regions[1]=m_BCRegion;
2868 m_BCRegion.SetIndex(InputDimension-1,2);
2869 m_BC2Regions[2]=m_BCRegion;
2870 m_BC2Values[2]=-m_WeightRatio[3][1];
2877 m_FirstRegion.SetSize(InputDimension-1,1);
2878 m_FirstRegion.SetIndex(InputDimension-1,3);
2881 m_BCRegions.resize(0);
2884 m_SecondRegion.SetSize(InputDimension-1,1);
2885 m_SecondRegion.SetIndex(InputDimension-1,0);
2888 m_BC2Regions.resize(3);
2889 m_BC2Values.resize(3);
2890 m_BCRegion.SetIndex(InputDimension-1,0);
2891 m_BC2Regions[0]=m_BCRegion;
2892 m_BC2Values[0]=-m_WeightRatio[3][1];
2893 m_BCRegion.SetIndex(InputDimension-1,1);
2894 m_BC2Regions[1]=m_BCRegion;
2896 m_BCRegion.SetIndex(InputDimension-1,2);
2897 m_BC2Regions[2]=m_BCRegion;
2898 m_BC2Values[2]=-m_WeightRatio[3][1];
2901 m_FirstRegion.SetSize(InputDimension-1,1);
2902 m_FirstRegion.SetIndex(InputDimension-1,1);
2909 m_FirstRegion.SetSize(InputDimension-1,1);
2910 m_FirstRegion.SetIndex(InputDimension-1,0);
2913 m_BCRegions.resize(3);
2914 m_BCValues.resize(3);
2915 m_BCRegion.SetIndex(InputDimension-1,0);
2916 m_BCRegions[0]=m_BCRegion;
2917 m_BCValues[0]=-m_WeightRatio[3][1];
2918 m_BCRegion.SetIndex(InputDimension-1,1);
2919 m_BCRegions[1]=m_BCRegion;
2921 m_BCRegion.SetIndex(InputDimension-1,2);
2922 m_BCRegions[2]=m_BCRegion;
2923 m_BCValues[2]=-m_WeightRatio[3][1];
2926 m_SecondRegion.SetSize(InputDimension-1,2);
2927 m_SecondRegion.SetIndex(InputDimension-1,1);
2930 m_BC2Regions.resize(0);
2937 m_FirstRegion.SetSize(InputDimension-1,0);
2940 m_BCRegions.resize(3);
2941 m_BCValues.resize(3);
2942 m_BCRegion.SetIndex(InputDimension-1,0);
2943 m_BCRegions[0]=m_BCRegion;
2944 m_BCValues[0]=-m_WeightRatio[3][1];
2945 m_BCRegion.SetIndex(InputDimension-1,1);
2946 m_BCRegions[1]=m_BCRegion;
2948 m_BCRegion.SetIndex(InputDimension-1,2);
2949 m_BCRegions[2]=m_BCRegion;
2950 m_BCValues[2]=-m_WeightRatio[3][1];
2953 m_SecondRegion.SetSize(InputDimension-1,3);
2954 m_SecondRegion.SetIndex(InputDimension-1,1);
2957 m_BC2Regions.resize(0);
2964 m_FirstRegion.SetSize(InputDimension-1,3);
2965 m_FirstRegion.SetIndex(InputDimension-1,1);
2968 m_BCRegions.resize(0);
2971 m_SecondRegion.SetSize(InputDimension-1,1);
2972 m_SecondRegion.SetIndex(InputDimension-1,0);
2975 m_BC2Regions.resize(0);
2982 DD("supportindex > 3 ???");
2983 DD(m_SupportRegion.GetIndex(InputDimension-1));
2985 } // end switch index
2989 // // ----------------------------------------------------------------------
2990 // // The egg with 5 internal CP:
2991 // // Periodic, C2 smooth everywhere and constrained to zero at the reference
2992 // // Coeff row R5 BC1 0 1 2 3 BC2 R2
2993 // // PaddedCoeff R: 0 1 2 3 4 5 6 7
2994 // // BC1= -weights[2]/weights[0] ( R2+R5)
2996 // // ---------------------------------------------------------------------
2997 // switch(m_SupportRegion.GetIndex(InputDimension-1))
3002 // m_FirstRegion.SetSize(InputDimension-1,1);
3003 // m_FirstRegion.SetIndex(InputDimension-1,3);
3006 // m_BCRegions.resize(2);
3007 // m_BCValues.resize(2);
3008 // m_BCRegion.SetIndex(InputDimension-1,0);
3009 // m_BCRegions[0]=m_BCRegion;
3010 // m_BCValues[0]=-m_WeightRatio[2][0];
3011 // m_BCRegion.SetIndex(InputDimension-1,3);
3012 // m_BCRegions[1]=m_BCRegion;
3013 // m_BCValues[1]=-m_WeightRatio[2][0];
3016 // m_SecondRegion.SetSize(InputDimension-1,2);
3017 // m_SecondRegion.SetIndex(InputDimension-1,0);
3020 // m_BC2Regions.resize(0);
3027 // m_FirstRegion.SetSize(InputDimension-1,0);
3030 // m_BCRegions.resize(2);
3031 // m_BCValues.resize(2);
3032 // m_BCRegion.SetIndex(InputDimension-1,0);
3033 // m_BCRegions[0]=m_BCRegion;
3034 // m_BCValues[0]=-m_WeightRatio[2][0];
3035 // m_BCRegion.SetIndex(InputDimension-1,3);
3036 // m_BCRegions[1]=m_BCRegion;
3037 // m_BCValues[1]=-m_WeightRatio[2][0];
3040 // m_SecondRegion.SetSize(InputDimension-1,3);
3041 // m_SecondRegion.SetIndex(InputDimension-1,0);
3044 // m_BC2Regions.resize(0);
3051 // m_FirstRegion.SetSize(InputDimension-1,4);
3052 // m_FirstRegion.SetIndex(InputDimension-1,0);
3055 // m_BCRegions.resize(0);
3058 // m_SecondRegion.SetSize(InputDimension-1,0);
3061 // m_BC2Regions.resize(0);
3068 // m_FirstRegion.SetSize(InputDimension-1,3);
3069 // m_FirstRegion.SetIndex(InputDimension-1,1);
3072 // m_BCRegions.resize(2);
3073 // m_BCValues.resize(2);
3074 // m_BCRegion.SetIndex(InputDimension-1,0);
3075 // m_BCRegions[0]=m_BCRegion;
3076 // m_BCValues[0]=-m_WeightRatio[2][0];
3077 // m_BCRegion.SetIndex(InputDimension-1,3);
3078 // m_BCRegions[1]=m_BCRegion;
3079 // m_BCValues[1]=-m_WeightRatio[2][0];
3082 // m_SecondRegion.SetSize(InputDimension-1,0);
3085 // m_BC2Regions.resize(0);
3093 // m_FirstRegion.SetSize(InputDimension-1,2);
3094 // m_FirstRegion.SetIndex(InputDimension-1,2);
3097 // m_BCRegions.resize(2);
3098 // m_BCValues.resize(2);
3099 // m_BCRegion.SetIndex(InputDimension-1,0);
3100 // m_BCRegions[0]=m_BCRegion;
3101 // m_BCValues[0]=-m_WeightRatio[2][0];
3102 // m_BCRegion.SetIndex(InputDimension-1,3);
3103 // m_BCRegions[1]=m_BCRegion;
3104 // m_BCValues[1]=-m_WeightRatio[2][0];
3107 // m_SecondRegion.SetSize(InputDimension-1,1);
3108 // m_SecondRegion.SetIndex(InputDimension-1,0);
3111 // m_BC2Regions.resize(0);
3118 // DD("supportindex > 4 ???");
3119 // DD(m_SupportRegion.GetIndex(InputDimension-1));
3120 // DD(m_TransformShape);
3122 // }// end swith index
3125 // } // end case 1 shape
3129 // ----------------------------------------------------------------------
3130 // The rabbit with 4 internal CP
3131 // Periodic, constrained to zero at the reference
3132 // at position 3 and the extremes fixed with anit-symm bc
3133 // Coeff row BC1 0 1 BC2 2 BC3 BC4
3134 // PaddedCoeff R: 0 1 2 3 4 5 6
3136 // BC2= -weights[2]/weights[0] ( R2+R4 )
3139 // ---------------------------------------------------------------------
3140 switch(m_SupportRegion.GetIndex(InputDimension-1))
3145 m_FirstRegion.SetSize(InputDimension-1,0);
3148 m_BCRegions.resize(2);
3149 m_BCValues.resize(2);
3150 m_BCRegion.SetIndex(InputDimension-1,0);
3151 m_BCRegions[0]=m_BCRegion;
3153 m_BCRegion.SetIndex(InputDimension-1,1);
3154 m_BCRegions[1]=m_BCRegion;
3158 m_SecondRegion.SetSize(InputDimension-1,2);
3159 m_SecondRegion.SetIndex(InputDimension-1,0);
3162 m_BC2Regions.resize(2);
3163 m_BC2Values.resize(2);
3164 m_BCRegion.SetIndex(InputDimension-1,1);
3165 m_BC2Regions[0]=m_BCRegion;
3166 m_BC2Values[0]=-m_WeightRatio[2][0];
3167 m_BCRegion.SetIndex(InputDimension-1,2);
3168 m_BC2Regions[1]=m_BCRegion;
3169 m_BC2Values[1]= -m_WeightRatio[2][0];
3176 m_FirstRegion.SetSize(InputDimension-1,2);
3177 m_FirstRegion.SetIndex(InputDimension-1,0);
3180 m_BCRegions.resize(2);
3181 m_BCValues.resize(2);
3182 m_BCRegion.SetIndex(InputDimension-1,1);
3183 m_BCRegions[0]=m_BCRegion;
3184 m_BCValues[0]=-m_WeightRatio[2][0];
3185 m_BCRegion.SetIndex(InputDimension-1,2);
3186 m_BCRegions[1]=m_BCRegion;
3187 m_BCValues[1]= -m_WeightRatio[2][0];
3190 m_SecondRegion.SetSize(InputDimension-1,1);
3191 m_SecondRegion.SetIndex(InputDimension-1,2);
3194 m_BC2Regions.resize(0);
3201 m_FirstRegion.SetSize(InputDimension-1,1);
3202 m_FirstRegion.SetIndex(InputDimension-1,1);
3205 m_BCRegions.resize(2);
3206 m_BCValues.resize(2);
3207 m_BCRegion.SetIndex(InputDimension-1,1);
3208 m_BCRegions[0]=m_BCRegion;
3209 m_BCValues[0]=-m_WeightRatio[2][0];
3210 m_BCRegion.SetIndex(InputDimension-1,2);
3211 m_BCRegions[1]=m_BCRegion;
3212 m_BCValues[1]= -m_WeightRatio[2][0];
3215 m_SecondRegion.SetSize(InputDimension-1,1);
3216 m_SecondRegion.SetIndex(InputDimension-1,2);
3219 m_BC2Regions.resize(1);
3220 m_BC2Values.resize(1);
3221 m_BCRegion.SetIndex(InputDimension-1,0);
3222 m_BC2Regions[0]=m_BCRegion;
3230 m_FirstRegion.SetSize(InputDimension-1,0);
3233 m_BCRegions.resize(2);
3234 m_BCValues.resize(2);
3235 m_BCRegion.SetIndex(InputDimension-1,1);
3236 m_BCRegions[0]=m_BCRegion;
3237 m_BCValues[0]=-m_WeightRatio[2][0];
3238 m_BCRegion.SetIndex(InputDimension-1,2);
3239 m_BCRegions[1]=m_BCRegion;
3240 m_BCValues[1]= -m_WeightRatio[2][0];
3243 m_SecondRegion.SetSize(InputDimension-1,1);
3244 m_SecondRegion.SetIndex(InputDimension-1,2);
3247 m_BC2Regions.resize(1);
3248 m_BC2Values.resize(1);
3249 m_BCRegion.SetIndex(InputDimension-1,0);
3250 m_BC2Regions[0]=m_BCRegion;
3254 m_BC3Regions.resize(2);
3255 m_BC3Values.resize(2);
3256 m_BCRegion.SetIndex(InputDimension-1,0);
3257 m_BC3Regions[0]=m_BCRegion;
3259 m_BCRegion.SetIndex(InputDimension-1,2);
3260 m_BC3Regions[1]=m_BCRegion;
3268 DD("supportindex > 3 ???");
3269 DD(m_SupportRegion.GetIndex(InputDimension-1));
3271 } // end switch index
3274 } // end case 2 shape
3277 // ----------------------------------------------------------------------
3278 // The rabbit with 5 internal CP
3279 // Periodic, constrained to zero at the reference at position 3.5
3280 // and the extremes fixed with anti-symmetrical BC
3281 // Coeff row BC1 0 1 BC2 2 3 BC3 BC4
3282 // PaddedCoeff R: 0 1 2 3 4 5 6 7
3284 // BC2= -weights[3]/weights[1] ( R2+R5 ) - R4
3287 // ---------------------------------------------------------------------
3288 switch(m_SupportRegion.GetIndex(InputDimension-1))
3293 m_FirstRegion.SetSize(InputDimension-1,0);
3296 m_BCRegions.resize(2);
3297 m_BCValues.resize(2);
3298 m_BCRegion.SetIndex(InputDimension-1,0);
3299 m_BCRegions[0]=m_BCRegion;
3301 m_BCRegion.SetIndex(InputDimension-1,1);
3302 m_BCRegions[1]=m_BCRegion;
3306 m_SecondRegion.SetSize(InputDimension-1,2);
3307 m_SecondRegion.SetIndex(InputDimension-1,0);
3310 m_BC2Regions.resize(3);
3311 m_BC2Values.resize(3);
3312 m_BCRegion.SetIndex(InputDimension-1,1);
3313 m_BC2Regions[0]=m_BCRegion;
3314 m_BC2Values[0]=-m_WeightRatio[3][1];
3315 m_BCRegion.SetIndex(InputDimension-1,2);
3316 m_BC2Regions[1]=m_BCRegion;
3318 m_BCRegion.SetIndex(InputDimension-1,3);
3319 m_BC2Regions[2]=m_BCRegion;
3320 m_BC2Values[2]=-m_WeightRatio[3][1];
3327 m_FirstRegion.SetSize(InputDimension-1,2);
3328 m_FirstRegion.SetIndex(InputDimension-1,0);
3331 m_BCRegions.resize(3);
3332 m_BCValues.resize(3);
3333 m_BCRegion.SetIndex(InputDimension-1,1);
3334 m_BCRegions[0]=m_BCRegion;
3335 m_BCValues[0]=-m_WeightRatio[3][1];
3336 m_BCRegion.SetIndex(InputDimension-1,2);
3337 m_BCRegions[1]=m_BCRegion;
3339 m_BCRegion.SetIndex(InputDimension-1,3);
3340 m_BCRegions[2]=m_BCRegion;
3341 m_BCValues[2]=-m_WeightRatio[3][1];
3344 m_SecondRegion.SetSize(InputDimension-1,1);
3345 m_SecondRegion.SetIndex(InputDimension-1,2);
3348 m_BC2Regions.resize(0);
3355 m_FirstRegion.SetSize(InputDimension-1,1);
3356 m_FirstRegion.SetIndex(InputDimension-1,1);
3359 m_BCRegions.resize(3);
3360 m_BCValues.resize(3);
3361 m_BCRegion.SetIndex(InputDimension-1,1);
3362 m_BCRegions[0]=m_BCRegion;
3363 m_BCValues[0]=-m_WeightRatio[3][1];
3364 m_BCRegion.SetIndex(InputDimension-1,2);
3365 m_BCRegions[1]=m_BCRegion;
3367 m_BCRegion.SetIndex(InputDimension-1,3);
3368 m_BCRegions[2]=m_BCRegion;
3369 m_BCValues[2]=-m_WeightRatio[3][1];
3372 m_SecondRegion.SetSize(InputDimension-1,2);
3373 m_SecondRegion.SetIndex(InputDimension-1,2);
3376 m_BC2Regions.resize(0);
3383 m_FirstRegion.SetSize(InputDimension-1,0);
3386 m_BCRegions.resize(3);
3387 m_BCValues.resize(3);
3388 m_BCRegion.SetIndex(InputDimension-1,1);
3389 m_BCRegions[0]=m_BCRegion;
3390 m_BCValues[0]=-m_WeightRatio[3][1];
3391 m_BCRegion.SetIndex(InputDimension-1,2);
3392 m_BCRegions[1]=m_BCRegion;
3394 m_BCRegion.SetIndex(InputDimension-1,3);
3395 m_BCRegions[2]=m_BCRegion;
3396 m_BCValues[2]=-m_WeightRatio[3][1];
3399 m_SecondRegion.SetSize(InputDimension-1,2);
3400 m_SecondRegion.SetIndex(InputDimension-1,2);
3403 m_BC2Regions.resize(1);
3404 m_BC2Values.resize(1);
3405 m_BCRegion.SetIndex(InputDimension-1,0);
3406 m_BC2Regions[0]=m_BCRegion;
3415 m_FirstRegion.SetSize(InputDimension-1,2);
3416 m_FirstRegion.SetIndex(InputDimension-1,2);
3419 m_BCRegions.resize(1);
3420 m_BCValues.resize(1);
3421 m_BCRegion.SetIndex(InputDimension-1,0);
3422 m_BCRegions[0]=m_BCRegion;
3426 m_SecondRegion.SetSize(InputDimension-1,0);
3429 m_BC2Regions.resize(2);
3430 m_BC2Values.resize(2);
3431 m_BCRegion.SetIndex(InputDimension-1,0);
3432 m_BC2Regions[0]=m_BCRegion;
3434 m_BCRegion.SetIndex(InputDimension-1,3);
3435 m_BC2Regions[1]=m_BCRegion;
3443 DD("supportindex > 4 ???");
3444 DD(m_SupportRegion.GetIndex(InputDimension-1));
3446 } // end switch index
3449 } // end case 3 shape
3453 // ----------------------------------------------------------------------
3454 // The sputnik with 4 internal CP
3455 // Periodic, constrained to zero at the reference
3456 // at position 3 and one indepent extremes copied
3457 // Coeff row BC1 0 1 BC2 2 BC2 3
3458 // PaddedCoeff R: 0 1 2 3 4 5 6
3460 // BC2= -weights[2]/weights[0] ( R2+R4 )
3461 // BC3= weights[2]/weights[0] ( R2-R4) + R1
3462 // ---------------------------------------------------------------------
3463 switch(m_SupportRegion.GetIndex(InputDimension-1))
3468 m_FirstRegion.SetSize(InputDimension-1,0);
3471 m_BCRegions.resize(1);
3472 m_BCValues.resize(1);
3473 m_BCRegion.SetIndex(InputDimension-1,3);
3474 m_BCRegions[0]=m_BCRegion;
3478 m_SecondRegion.SetSize(InputDimension-1,2);
3479 m_SecondRegion.SetIndex(InputDimension-1,0);
3482 m_BC2Regions.resize(2);
3483 m_BC2Values.resize(2);
3484 m_BCRegion.SetIndex(InputDimension-1,1);
3485 m_BC2Regions[0]=m_BCRegion;
3486 m_BC2Values[0]=-m_WeightRatio[2][0];
3487 m_BCRegion.SetIndex(InputDimension-1,2);
3488 m_BC2Regions[1]=m_BCRegion;
3489 m_BC2Values[1]= -m_WeightRatio[2][0];
3496 m_FirstRegion.SetSize(InputDimension-1,2);
3497 m_FirstRegion.SetIndex(InputDimension-1,0);
3500 m_BCRegions.resize(2);
3501 m_BCValues.resize(2);
3502 m_BCRegion.SetIndex(InputDimension-1,1);
3503 m_BCRegions[0]=m_BCRegion;
3504 m_BCValues[0]=-m_WeightRatio[2][0];
3505 m_BCRegion.SetIndex(InputDimension-1,2);
3506 m_BCRegions[1]=m_BCRegion;
3507 m_BCValues[1]= -m_WeightRatio[2][0];
3510 m_SecondRegion.SetSize(InputDimension-1,1);
3511 m_SecondRegion.SetIndex(InputDimension-1,2);
3514 m_BC2Regions.resize(0);
3521 m_FirstRegion.SetSize(InputDimension-1,1);
3522 m_FirstRegion.SetIndex(InputDimension-1,1);
3525 m_BCRegions.resize(2);
3526 m_BCValues.resize(2);
3527 m_BCRegion.SetIndex(InputDimension-1,1);
3528 m_BCRegions[0]=m_BCRegion;
3529 m_BCValues[0]=-m_WeightRatio[2][0];
3530 m_BCRegion.SetIndex(InputDimension-1,2);
3531 m_BCRegions[1]=m_BCRegion;
3532 m_BCValues[1]= -m_WeightRatio[2][0];
3535 m_SecondRegion.SetSize(InputDimension-1,1);
3536 m_SecondRegion.SetIndex(InputDimension-1,2);
3539 m_BC2Regions.resize(3);
3540 m_BC2Values.resize(3);
3541 m_BCRegion.SetIndex(InputDimension-1,0);
3542 m_BC2Regions[0]=m_BCRegion;
3544 m_BCRegion.SetIndex(InputDimension-1,1);
3545 m_BC2Regions[1]=m_BCRegion;
3546 m_BC2Values[1]=m_WeightRatio[2][0];
3547 m_BCRegion.SetIndex(InputDimension-1,2);
3548 m_BC2Regions[2]=m_BCRegion;
3549 m_BC2Values[2]=-m_WeightRatio[2][0];
3556 m_FirstRegion.SetSize(InputDimension-1,0);
3559 m_BCRegions.resize(2);
3560 m_BCValues.resize(2);
3561 m_BCRegion.SetIndex(InputDimension-1,1);
3562 m_BCRegions[0]=m_BCRegion;
3563 m_BCValues[0]=-m_WeightRatio[2][0];
3564 m_BCRegion.SetIndex(InputDimension-1,2);
3565 m_BCRegions[1]=m_BCRegion;
3566 m_BCValues[1]= -m_WeightRatio[2][0];
3569 m_SecondRegion.SetSize(InputDimension-1,1);
3570 m_SecondRegion.SetIndex(InputDimension-1,2);
3573 m_BC2Regions.resize(3);
3574 m_BC2Values.resize(3);
3575 m_BCRegion.SetIndex(InputDimension-1,0);
3576 m_BC2Regions[0]=m_BCRegion;
3578 m_BCRegion.SetIndex(InputDimension-1,1);
3579 m_BC2Regions[1]=m_BCRegion;
3580 m_BC2Values[1]=m_WeightRatio[2][0];
3581 m_BCRegion.SetIndex(InputDimension-1,2);
3582 m_BC2Regions[2]=m_BCRegion;
3583 m_BC2Values[2]=-m_WeightRatio[2][0];
3586 m_ThirdRegion.SetSize(InputDimension-1,1);
3587 m_ThirdRegion.SetIndex(InputDimension-1,3);
3594 DD("supportindex > 3 ???");
3595 DD(m_SupportRegion.GetIndex(InputDimension-1));
3597 } // end switch index
3600 } // end case 4 shape
3604 // ----------------------------------------------------------------------
3605 // The sputnik with 5 internal CP
3606 // Periodic, constrained to zero at the reference
3607 // at position 3 and one indepent extreme
3608 // Coeff row BC1 0 1 BC2 2 3 BC3 4
3609 // PaddedCoeff R: 0 1 2 3 4 5 6 7
3610 // BC1= R2 + R5 - R7
3611 // BC2= -weights[3]/weights[1] ( R2+R5 ) - R4
3612 // BC3= R1 + 0.5 R2 - 0.5 R7
3613 // ----------------------------------------------------------------------
3614 switch(m_SupportRegion.GetIndex(InputDimension-1))
3619 m_FirstRegion.SetSize(InputDimension-1,0);
3622 m_BCRegions.resize(3);
3623 m_BCValues.resize(3);
3624 m_BCRegion.SetIndex(InputDimension-1,1);
3625 m_BCRegions[0]=m_BCRegion;
3627 m_BCRegion.SetIndex(InputDimension-1,3);
3628 m_BCRegions[1]=m_BCRegion;
3630 m_BCRegion.SetIndex(InputDimension-1,4);
3631 m_BCRegions[2]=m_BCRegion;
3635 m_SecondRegion.SetSize(InputDimension-1,2);
3636 m_SecondRegion.SetIndex(InputDimension-1,0);
3639 m_BC2Regions.resize(3);
3640 m_BC2Values.resize(3);
3641 m_BCRegion.SetIndex(InputDimension-1,1);
3642 m_BC2Regions[0]=m_BCRegion;
3643 m_BC2Values[0]=-m_WeightRatio[3][1];
3644 m_BCRegion.SetIndex(InputDimension-1,2);
3645 m_BC2Regions[1]=m_BCRegion;
3647 m_BCRegion.SetIndex(InputDimension-1,3);
3648 m_BC2Regions[2]=m_BCRegion;
3649 m_BC2Values[2]=-m_WeightRatio[3][1];
3656 m_FirstRegion.SetSize(InputDimension-1,2);
3657 m_FirstRegion.SetIndex(InputDimension-1,0);
3660 m_BCRegions.resize(3);
3661 m_BCValues.resize(3);
3662 m_BCRegion.SetIndex(InputDimension-1,1);
3663 m_BCRegions[0]=m_BCRegion;
3664 m_BCValues[0]=-m_WeightRatio[3][1];
3665 m_BCRegion.SetIndex(InputDimension-1,2);
3666 m_BCRegions[1]=m_BCRegion;
3668 m_BCRegion.SetIndex(InputDimension-1,3);
3669 m_BCRegions[2]=m_BCRegion;
3670 m_BCValues[2]=-m_WeightRatio[3][1];
3673 m_SecondRegion.SetSize(InputDimension-1,1);
3674 m_SecondRegion.SetIndex(InputDimension-1,2);
3677 m_BC2Regions.resize(0);
3684 m_FirstRegion.SetSize(InputDimension-1,1);
3685 m_FirstRegion.SetIndex(InputDimension-1,1);
3688 m_BCRegions.resize(3);
3689 m_BCValues.resize(3);
3690 m_BCRegion.SetIndex(InputDimension-1,1);
3691 m_BCRegions[0]=m_BCRegion;
3692 m_BCValues[0]=-m_WeightRatio[3][1];
3693 m_BCRegion.SetIndex(InputDimension-1,2);
3694 m_BCRegions[1]=m_BCRegion;
3696 m_BCRegion.SetIndex(InputDimension-1,3);
3697 m_BCRegions[2]=m_BCRegion;
3698 m_BCValues[2]=-m_WeightRatio[3][1];
3701 m_SecondRegion.SetSize(InputDimension-1,2);
3702 m_SecondRegion.SetIndex(InputDimension-1,2);
3705 m_BC2Regions.resize(0);
3712 m_FirstRegion.SetSize(InputDimension-1,0);
3715 m_BCRegions.resize(3);
3716 m_BCValues.resize(3);
3717 m_BCRegion.SetIndex(InputDimension-1,1);
3718 m_BCRegions[0]=m_BCRegion;
3719 m_BCValues[0]=-m_WeightRatio[3][1];
3720 m_BCRegion.SetIndex(InputDimension-1,2);
3721 m_BCRegions[1]=m_BCRegion;
3723 m_BCRegion.SetIndex(InputDimension-1,3);
3724 m_BCRegions[2]=m_BCRegion;
3725 m_BCValues[2]=-m_WeightRatio[3][1];
3728 m_SecondRegion.SetSize(InputDimension-1,2);
3729 m_SecondRegion.SetIndex(InputDimension-1,2);
3732 m_BC2Regions.resize(3);
3733 m_BC2Values.resize(3);
3734 m_BCRegion.SetIndex(InputDimension-1,0);
3735 m_BC2Regions[0]=m_BCRegion;
3737 m_BCRegion.SetIndex(InputDimension-1,1);
3738 m_BC2Regions[1]=m_BCRegion;
3740 m_BCRegion.SetIndex(InputDimension-1,4);
3741 m_BC2Regions[2]=m_BCRegion;
3742 m_BC2Values[2]=-0.5;
3749 m_FirstRegion.SetSize(InputDimension-1,2);
3750 m_FirstRegion.SetIndex(InputDimension-1,2);
3753 m_BCRegions.resize(3);
3754 m_BCValues.resize(3);
3755 m_BCRegion.SetIndex(InputDimension-1,0);
3756 m_BCRegions[0]=m_BCRegion;
3758 m_BCRegion.SetIndex(InputDimension-1,1);
3759 m_BCRegions[1]=m_BCRegion;
3761 m_BCRegion.SetIndex(InputDimension-1,4);
3762 m_BCRegions[2]=m_BCRegion;
3766 m_SecondRegion.SetSize(InputDimension-1,1);
3767 m_SecondRegion.SetIndex(InputDimension-1,4);
3770 m_BC2Regions.resize(0);
3776 DD("supportindex > 4 ???");
3777 DD(m_SupportRegion.GetIndex(InputDimension-1));
3779 } // end switch index
3782 } // end case 5 shape
3789 // ----------------------------------------------------------------------
3790 // The diamond with 4 internal CP:
3791 // Periodic, constrained to zero at the reference at position 3
3792 // Coeff row 0 1 2 BC1 3 BC2 4
3793 // PaddedCoeff R:0 1 2 3 4 5 6
3794 // BC1= -weights[2]/weights[0] ( R2+R4 )
3795 // BC2= weights[2]/weights[0] ( R0+R2-R4-R6 ) + R1
3796 // ---------------------------------------------------------------------
3797 switch(m_SupportRegion.GetIndex(InputDimension-1))
3802 m_FirstRegion.SetSize(InputDimension-1,3);
3803 m_FirstRegion.SetIndex(InputDimension-1,0);
3806 m_BCRegions.resize(2);
3807 m_BCValues.resize(2);
3808 m_BCRegion.SetIndex(InputDimension-1,2);
3809 m_BCRegions[0]=m_BCRegion;
3810 m_BCValues[0]=-m_WeightRatio[2][0];
3811 m_BCRegion.SetIndex(InputDimension-1,3);
3812 m_BCRegions[1]=m_BCRegion;
3813 m_BCValues[1]=-m_WeightRatio[2][0];
3816 m_SecondRegion.SetSize(InputDimension-1,0);
3819 m_BC2Regions.resize(0);
3826 m_FirstRegion.SetSize(InputDimension-1,2);
3827 m_FirstRegion.SetIndex(InputDimension-1,1);
3830 m_BCRegions.resize(2);
3831 m_BCValues.resize(2);
3832 m_BCRegion.SetIndex(InputDimension-1,2);
3833 m_BCRegions[0]=m_BCRegion;
3834 m_BCValues[0]=-m_WeightRatio[2][0];
3835 m_BCRegion.SetIndex(InputDimension-1,3);
3836 m_BCRegions[1]=m_BCRegion;
3837 m_BCValues[1]=-m_WeightRatio[2][0];
3840 m_SecondRegion.SetSize(InputDimension-1,1);
3841 m_SecondRegion.SetIndex(InputDimension-1,3);
3844 m_BC2Regions.resize(0);
3851 m_FirstRegion.SetSize(InputDimension-1,1);
3852 m_FirstRegion.SetIndex(InputDimension-1,2);
3855 m_BCRegions.resize(2);
3856 m_BCValues.resize(2);
3857 m_BCRegion.SetIndex(InputDimension-1,2);
3858 m_BCRegions[0]=m_BCRegion;
3859 m_BCValues[0]=-m_WeightRatio[2][0];
3860 m_BCRegion.SetIndex(InputDimension-1,3);
3861 m_BCRegions[1]=m_BCRegion;
3862 m_BCValues[1]=-m_WeightRatio[2][0];
3865 m_SecondRegion.SetSize(InputDimension-1,1);
3866 m_SecondRegion.SetIndex(InputDimension-1,3);
3869 m_BC2Regions.resize(5);
3870 m_BC2Values.resize(5);
3871 m_BCRegion.SetIndex(InputDimension-1,0);
3872 m_BC2Regions[0]=m_BCRegion;
3873 m_BC2Values[0]=m_WeightRatio[2][0];
3874 m_BCRegion.SetIndex(InputDimension-1,1);
3875 m_BC2Regions[1]=m_BCRegion;
3877 m_BCRegion.SetIndex(InputDimension-1,2);
3878 m_BC2Regions[2]=m_BCRegion;
3879 m_BC2Values[2]=m_WeightRatio[2][0];
3880 m_BCRegion.SetIndex(InputDimension-1,3);
3881 m_BC2Regions[3]=m_BCRegion;
3882 m_BC2Values[3]=-m_WeightRatio[2][0];
3883 m_BCRegion.SetIndex(InputDimension-1,4);
3884 m_BC2Regions[4]=m_BCRegion;
3885 m_BC2Values[4]=-m_WeightRatio[2][0];
3892 m_FirstRegion.SetSize(InputDimension-1,0);
3895 m_BCRegions.resize(2);
3896 m_BCValues.resize(2);
3897 m_BCRegion.SetIndex(InputDimension-1,2);
3898 m_BCRegions[0]=m_BCRegion;
3899 m_BCValues[0]=-m_WeightRatio[2][0];
3900 m_BCRegion.SetIndex(InputDimension-1,3);
3901 m_BCRegions[1]=m_BCRegion;
3902 m_BCValues[1]=-m_WeightRatio[2][0];
3905 m_SecondRegion.SetSize(InputDimension-1,1);
3906 m_SecondRegion.SetIndex(InputDimension-1,3);
3909 m_BC2Regions.resize(5);
3910 m_BC2Values.resize(5);
3911 m_BCRegion.SetIndex(InputDimension-1,0);
3912 m_BC2Regions[0]=m_BCRegion;
3913 m_BC2Values[0]=m_WeightRatio[2][0];
3914 m_BCRegion.SetIndex(InputDimension-1,1);
3915 m_BC2Regions[1]=m_BCRegion;
3917 m_BCRegion.SetIndex(InputDimension-1,2);
3918 m_BC2Regions[2]=m_BCRegion;
3919 m_BC2Values[2]=m_WeightRatio[2][0];
3920 m_BCRegion.SetIndex(InputDimension-1,3);
3921 m_BC2Regions[3]=m_BCRegion;
3922 m_BC2Values[3]=-m_WeightRatio[2][0];
3923 m_BCRegion.SetIndex(InputDimension-1,4);
3924 m_BC2Regions[4]=m_BCRegion;
3925 m_BC2Values[4]=-m_WeightRatio[2][0];
3928 m_ThirdRegion.SetSize(InputDimension-1,1);
3929 m_ThirdRegion.SetIndex(InputDimension-1,4);
3936 DD("supportindex > 3 ???");
3937 DD(m_SupportRegion.GetIndex(InputDimension-1));
3939 } // end switch index
3942 } // end case 7 shape
3946 // ----------------------------------------------------------------------
3947 // The diamond with 5 internal CP:
3948 // periodic, constrained to zero at the reference at position 3.5
3949 // Coeff row 0 1 2 BC1 3 4 BC2 5
3950 // PaddedCoeff R:0 1 2 3 4 5 6 7
3951 // BC1= -weights[3]/weights[1] ( R2+R5 ) - R4
3952 // BC2= weights[2]/weights[0] ( R0+R2-R5-R7 ) + R1
3953 // ---------------------------------------------------------------------
3954 switch(m_SupportRegion.GetIndex(InputDimension-1))
3959 m_FirstRegion.SetSize(InputDimension-1,3);
3960 m_FirstRegion.SetIndex(InputDimension-1,0);
3963 m_BCRegions.resize(3);
3964 m_BCValues.resize(3);
3965 m_BCRegion.SetIndex(InputDimension-1,2);
3966 m_BCRegions[0]=m_BCRegion;
3967 m_BCValues[0]=-m_WeightRatio[3][1];
3968 m_BCRegion.SetIndex(InputDimension-1,3);
3969 m_BCRegions[1]=m_BCRegion;
3971 m_BCRegion.SetIndex(InputDimension-1,4);
3972 m_BCRegions[2]=m_BCRegion;
3973 m_BCValues[2]=-m_WeightRatio[3][1];
3976 m_SecondRegion.SetSize(InputDimension-1,0);
3979 m_BC2Regions.resize(0);
3986 m_FirstRegion.SetSize(InputDimension-1,2);
3987 m_FirstRegion.SetIndex(InputDimension-1,1);
3990 m_BCRegions.resize(3);
3991 m_BCValues.resize(3);
3992 m_BCRegion.SetIndex(InputDimension-1,2);
3993 m_BCRegions[0]=m_BCRegion;
3994 m_BCValues[0]=-m_WeightRatio[3][1];
3995 m_BCRegion.SetIndex(InputDimension-1,3);
3996 m_BCRegions[1]=m_BCRegion;
3998 m_BCRegion.SetIndex(InputDimension-1,4);
3999 m_BCRegions[2]=m_BCRegion;
4000 m_BCValues[2]=-m_WeightRatio[3][1];
4003 m_SecondRegion.SetSize(InputDimension-1,1);
4004 m_SecondRegion.SetIndex(InputDimension-1,3);
4007 m_BC2Regions.resize(0);
4014 m_FirstRegion.SetSize(InputDimension-1,1);
4015 m_FirstRegion.SetIndex(InputDimension-1,2);
4018 m_BCRegions.resize(3);
4019 m_BCValues.resize(3);
4020 m_BCRegion.SetIndex(InputDimension-1,2);
4021 m_BCRegions[0]=m_BCRegion;
4022 m_BCValues[0]=-m_WeightRatio[3][1];
4023 m_BCRegion.SetIndex(InputDimension-1,3);
4024 m_BCRegions[1]=m_BCRegion;
4026 m_BCRegion.SetIndex(InputDimension-1,4);
4027 m_BCRegions[2]=m_BCRegion;
4028 m_BCValues[2]=-m_WeightRatio[3][1];
4031 m_SecondRegion.SetSize(InputDimension-1,2);
4032 m_SecondRegion.SetIndex(InputDimension-1,3);
4035 m_BC2Regions.resize(0);
4042 m_FirstRegion.SetSize(InputDimension-1,0);
4043 m_FirstRegion.SetIndex(InputDimension-1,0);
4046 m_BCRegions.resize(3);
4047 m_BCValues.resize(3);
4048 m_BCRegion.SetIndex(InputDimension-1,2);
4049 m_BCRegions[0]=m_BCRegion;
4050 m_BCValues[0]=-m_WeightRatio[3][1];
4051 m_BCRegion.SetIndex(InputDimension-1,3);
4052 m_BCRegions[1]=m_BCRegion;
4054 m_BCRegion.SetIndex(InputDimension-1,4);
4055 m_BCRegions[2]=m_BCRegion;
4056 m_BCValues[2]=-m_WeightRatio[3][1];
4059 m_SecondRegion.SetSize(InputDimension-1,2);
4060 m_SecondRegion.SetIndex(InputDimension-1,3);
4063 m_BC2Regions.resize(5);
4064 m_BC2Values.resize(5);
4065 m_BCRegion.SetIndex(InputDimension-1,0);
4066 m_BC2Regions[0]=m_BCRegion;
4067 m_BC2Values[0]=m_WeightRatio[2][0];
4068 m_BCRegion.SetIndex(InputDimension-1,1);
4069 m_BC2Regions[1]=m_BCRegion;
4071 m_BCRegion.SetIndex(InputDimension-1,2);
4072 m_BC2Regions[2]=m_BCRegion;
4073 m_BC2Values[2]=m_WeightRatio[2][0];
4074 m_BCRegion.SetIndex(InputDimension-1,4);
4075 m_BC2Regions[3]=m_BCRegion;
4076 m_BC2Values[3]=-m_WeightRatio[2][0];
4077 m_BCRegion.SetIndex(InputDimension-1,5);
4078 m_BC2Regions[4]=m_BCRegion;
4079 m_BC2Values[4]=-m_WeightRatio[2][0];
4087 m_FirstRegion.SetSize(InputDimension-1,2);
4088 m_FirstRegion.SetIndex(InputDimension-1,3);
4091 m_BCRegions.resize(5);
4092 m_BCValues.resize(5);
4093 m_BCRegion.SetIndex(InputDimension-1,0);
4094 m_BCRegions[0]=m_BCRegion;
4095 m_BCValues[0]=m_WeightRatio[2][0];
4096 m_BCRegion.SetIndex(InputDimension-1,1);
4097 m_BCRegions[1]=m_BCRegion;
4099 m_BCRegion.SetIndex(InputDimension-1,2);
4100 m_BCRegions[2]=m_BCRegion;
4101 m_BCValues[2]=m_WeightRatio[2][0];
4102 m_BCRegion.SetIndex(InputDimension-1,4);
4103 m_BCRegions[3]=m_BCRegion;
4104 m_BCValues[3]=-m_WeightRatio[2][0];
4105 m_BCRegion.SetIndex(InputDimension-1,5);
4106 m_BCRegions[4]=m_BCRegion;
4107 m_BCValues[4]=-m_WeightRatio[2][0];
4110 m_SecondRegion.SetSize(InputDimension-1,1);
4111 m_SecondRegion.SetIndex(InputDimension-1,5);
4114 m_BC2Regions.resize(0);
4121 DD("supportindex > 4 ???");
4122 DD(m_SupportRegion.GetIndex(InputDimension-1));
4124 } // end switch index
4127 } // end case 7 shape
4131 // ----------------------------------------------------------------------
4132 // The sputnik with 5 internal CP
4133 // Periodic, constrained to zero at the reference
4134 // at position 3 and one indepent extreme
4135 // Coeff row BC1 0 1 BC2 2 3 BC3 4
4136 // PaddedCoeff R: 0 1 2 3 4 5 6 7
4138 // BC2= -weights[3]/weights[1] ( R2+R5 ) - R4
4140 // ---------------------------------------------------------------------
4141 switch(m_SupportRegion.GetIndex(InputDimension-1))
4146 m_FirstRegion.SetSize(InputDimension-1,0);
4149 m_BCRegions.resize(3);
4150 m_BCValues.resize(3);
4151 m_BCRegion.SetIndex(InputDimension-1,1);
4152 m_BCRegions[0]=m_BCRegion;
4154 m_BCRegion.SetIndex(InputDimension-1,3);
4155 m_BCRegions[1]=m_BCRegion;
4157 m_BCRegion.SetIndex(InputDimension-1,4);
4158 m_BCRegions[2]=m_BCRegion;
4162 m_SecondRegion.SetSize(InputDimension-1,2);
4163 m_SecondRegion.SetIndex(InputDimension-1,0);
4166 m_BC2Regions.resize(3);
4167 m_BC2Values.resize(3);
4168 m_BCRegion.SetIndex(InputDimension-1,1);
4169 m_BC2Regions[0]=m_BCRegion;
4170 m_BC2Values[0]=-m_WeightRatio[3][1];
4171 m_BCRegion.SetIndex(InputDimension-1,2);
4172 m_BC2Regions[1]=m_BCRegion;
4174 m_BCRegion.SetIndex(InputDimension-1,3);
4175 m_BC2Regions[2]=m_BCRegion;
4176 m_BC2Values[2]=-m_WeightRatio[3][1];
4183 m_FirstRegion.SetSize(InputDimension-1,2);
4184 m_FirstRegion.SetIndex(InputDimension-1,0);
4187 m_BCRegions.resize(3);
4188 m_BCValues.resize(3);
4189 m_BCRegion.SetIndex(InputDimension-1,1);
4190 m_BCRegions[0]=m_BCRegion;
4191 m_BCValues[0]=-m_WeightRatio[3][1];
4192 m_BCRegion.SetIndex(InputDimension-1,2);
4193 m_BCRegions[1]=m_BCRegion;
4195 m_BCRegion.SetIndex(InputDimension-1,3);
4196 m_BCRegions[2]=m_BCRegion;
4197 m_BCValues[2]=-m_WeightRatio[3][1];
4200 m_SecondRegion.SetSize(InputDimension-1,1);
4201 m_SecondRegion.SetIndex(InputDimension-1,2);
4204 m_BC2Regions.resize(0);
4211 m_FirstRegion.SetSize(InputDimension-1,1);
4212 m_FirstRegion.SetIndex(InputDimension-1,1);
4215 m_BCRegions.resize(3);
4216 m_BCValues.resize(3);
4217 m_BCRegion.SetIndex(InputDimension-1,1);
4218 m_BCRegions[0]=m_BCRegion;
4219 m_BCValues[0]=-m_WeightRatio[3][1];
4220 m_BCRegion.SetIndex(InputDimension-1,2);
4221 m_BCRegions[1]=m_BCRegion;
4223 m_BCRegion.SetIndex(InputDimension-1,3);
4224 m_BCRegions[2]=m_BCRegion;
4225 m_BCValues[2]=-m_WeightRatio[3][1];
4228 m_SecondRegion.SetSize(InputDimension-1,2);
4229 m_SecondRegion.SetIndex(InputDimension-1,2);
4232 m_BC2Regions.resize(0);
4239 m_FirstRegion.SetSize(InputDimension-1,0);
4242 m_BCRegions.resize(3);
4243 m_BCValues.resize(3);
4244 m_BCRegion.SetIndex(InputDimension-1,1);
4245 m_BCRegions[0]=m_BCRegion;
4246 m_BCValues[0]=-m_WeightRatio[3][1];
4247 m_BCRegion.SetIndex(InputDimension-1,2);
4248 m_BCRegions[1]=m_BCRegion;
4250 m_BCRegion.SetIndex(InputDimension-1,3);
4251 m_BCRegions[2]=m_BCRegion;
4252 m_BCValues[2]=-m_WeightRatio[3][1];
4255 m_SecondRegion.SetSize(InputDimension-1,1);
4256 m_SecondRegion.SetIndex(InputDimension-1,2);
4259 m_BC2Regions.resize(1);
4260 m_BC2Values.resize(1);
4261 m_BCRegion.SetIndex(InputDimension-1,0);
4262 m_BC2Regions[0]=m_BCRegion;
4270 m_FirstRegion.SetSize(InputDimension-1,2);
4271 m_FirstRegion.SetIndex(InputDimension-1,2);
4274 m_BCRegions.resize(1);
4275 m_BCValues.resize(1);
4276 m_BCRegion.SetIndex(InputDimension-1,0);
4277 m_BCRegions[0]=m_BCRegion;
4281 m_SecondRegion.SetSize(InputDimension-1,1);
4282 m_SecondRegion.SetIndex(InputDimension-1,4);
4285 m_BC2Regions.resize(0);
4291 DD("supportindex > 4 ???");
4292 DD(m_SupportRegion.GetIndex(InputDimension-1));
4294 } // end switch index
4297 } // end case 9 shape
4302 DD ("Other shapes currently not implemented");
4305 } // end switch shape
4306 } // end wrap region
4310 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
4312 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions,NOutputDimensions>
4313 ::TransformPointToContinuousIndex( const InputPointType & point, ContinuousIndexType & index ) const
4317 itk::Vector<double, OutputDimension> tvector;
4319 for ( j = 0; j < OutputDimension; j++ )
4321 //JV find the index in the PADDED version
4322 tvector[j] = point[j] - this->m_PaddedGridOrigin[j];
4325 itk::Vector<double, OutputDimension> cvector;
4327 cvector = m_PointToIndex * tvector;
4329 for ( j = 0; j < OutputDimension; j++ )
4331 index[j] = static_cast< typename ContinuousIndexType::CoordRepType >( cvector[j] );