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 #if ITK_VERSION_MAJOR >= 4
36 ::ShapedBLUTSpatioTemporalDeformableTransform():Superclass(0)
38 ::ShapedBLUTSpatioTemporalDeformableTransform():Superclass(OutputDimension,0)
43 //JV default spline order
44 for ( i=0;i<InputDimension; i++)
47 //JV default LUTSamplingfactor
48 for ( i=0;i<InputDimension; i++)
49 m_LUTSamplingFactors[i]=20;
51 for (i=0;i<InputDimension; i++)
52 m_SupportSize[i] = m_SplineOrders[i]+1;
54 //Instantiate a Vector Bspline Interpolator
55 m_VectorInterpolator =VectorInterpolatorType::New();
56 m_VectorInterpolator->SetLUTSamplingFactors(m_LUTSamplingFactors);
57 m_VectorInterpolator->SetSplineOrders(m_SplineOrders);
59 // Set Bulk transform to NULL (no bulk performed)
60 m_BulkTransform = NULL;
68 // Default grid size is zero
72 //JV region containing the parameters
73 m_GridRegion.SetSize( m_NullSize);
74 m_GridRegion.SetIndex( m_NullIndex );
76 //JV use second region over the images
77 m_PaddedGridRegion.SetSize(m_NullSize);
78 m_PaddedGridRegion.SetIndex(m_NullIndex);
80 //JV Maintain two origins
81 m_GridOrigin.Fill( 0.0 ); // default origin is all zeros
82 m_PaddedGridOrigin.Fill( 0.0 ); // default origin is all zeros
84 m_GridSpacing.Fill( 1.0 ); // default spacing is all ones
85 m_GridDirection.SetIdentity(); // default spacing is all ones
87 m_InternalParametersBuffer = ParametersType(0);
88 // Make sure the parameters pointer is not NULL after construction.
89 m_InputParametersPointer = &m_InternalParametersBuffer;
91 // Initialize coeffient images
92 m_WrappedImage = CoefficientImageType::New();
93 m_WrappedImage->SetRegions( m_GridRegion );
94 m_WrappedImage->SetOrigin( m_GridOrigin.GetDataPointer() );
95 m_WrappedImage->SetSpacing( m_GridSpacing.GetDataPointer() );
96 m_WrappedImage->SetDirection( m_GridDirection );
98 // JV Initialize the padded version
99 m_PaddedCoefficientImage = CoefficientImageType::New();
100 m_PaddedCoefficientImage->SetRegions( m_PaddedGridRegion );
101 m_PaddedCoefficientImage->SetOrigin( m_GridOrigin.GetDataPointer() );
102 m_PaddedCoefficientImage->SetSpacing( m_GridSpacing.GetDataPointer() );
103 m_PaddedCoefficientImage->SetDirection( m_GridDirection );
105 m_CoefficientImage = NULL;
107 // Variables for computing interpolation
108 for (i=0; i <InputDimension;i++)
110 m_Offset[i] = m_SplineOrders[i] / 2;
111 if ( m_SplineOrders[i] % 2 )
113 m_SplineOrderOdd[i] = true;
117 m_SplineOrderOdd[i] = false;
121 //JV padded should be used when checking
122 m_ValidRegion = m_PaddedGridRegion;
124 // Initialize jacobian images
125 for (i=0; i <OutputDimension;i++)
127 m_JacobianImage[i] = JacobianImageType::New();
128 m_JacobianImage[i]->SetRegions( m_GridRegion );
129 m_JacobianImage[i]->SetOrigin( m_GridOrigin.GetDataPointer() );
130 m_JacobianImage[i]->SetSpacing( m_GridSpacing.GetDataPointer() );
131 m_JacobianImage[i]->SetDirection( m_GridDirection );
134 /** Fixed Parameters store the following information:
139 //JV we add the splineOrders, LUTsamplingfactor, m_Mask and m_BulkTransform
145 The size of these is equal to the NInputDimensions
146 *********************************************************/
147 this->m_FixedParameters.SetSize ( NInputDimensions * (NInputDimensions + 5)+3 );
148 this->m_FixedParameters.Fill ( 0.0 );
149 for ( i=0; i<NInputDimensions; i++)
151 this->m_FixedParameters[2*NInputDimensions+i] = m_GridSpacing[i];
153 for (unsigned int di=0; di<NInputDimensions; di++)
155 for (unsigned int dj=0; dj<NInputDimensions; dj++)
157 this->m_FixedParameters[3*NInputDimensions+(di*NInputDimensions+dj)] = m_GridDirection[di][dj];
161 //JV add splineOrders
162 for ( i=0; i<NInputDimensions; i++)
164 this->m_FixedParameters[ ( (3+NInputDimensions)*NInputDimensions)+i] = (this->GetSplineOrders())[i];
167 //JV add LUTsamplingFactors
168 for ( i=0; i<NInputDimensions; i++)
170 this->m_FixedParameters[( (4+NInputDimensions)*NInputDimensions)+i ] = this->GetLUTSamplingFactors()[i];
173 // JV add the mask pointer
174 this->m_FixedParameters[( (5+NInputDimensions)*NInputDimensions)]=(double)((size_t)m_Mask);
176 // JV add the bulkTransform pointer
177 this->m_FixedParameters[( (5+NInputDimensions)*NInputDimensions) +1]=(double)((size_t)m_BulkTransform);
179 // JV add the Transform shape
180 this->m_FixedParameters[( (5+NInputDimensions)*NInputDimensions) +2]=(double)(m_TransformShape);
183 // Calculate the PointToIndex matrices
185 for( unsigned int i=0; i<OutputDimension; i++)
187 scale[i][i] = m_GridSpacing[i];
190 m_IndexToPoint = m_GridDirection * scale;
191 m_PointToIndex = m_IndexToPoint.GetInverse();
193 m_LastJacobianIndex = m_ValidRegion.GetIndex();
196 //Weights to be used for the BC
198 std::vector<double> weights(5);
207 m_Weights[0]=weights;
215 m_Weights[1]=weights;
223 m_Weights[2]=weights;
231 m_Weights[3]=weights;
233 // Update the WeightRatios
234 m_WeightRatio.resize(4);
235 for (unsigned int i=0; i<4; i++)
237 m_WeightRatio[i].resize(4);
238 for (unsigned int j=0; j<4; j++)
239 m_WeightRatio[i][j]=m_Weights[m_SplineOrders[InputDimension-1]][i]/ m_Weights[m_SplineOrders[InputDimension-1]][j];
242 //JV initialize some variables for jacobian calculation
243 m_SupportRegion.SetSize(m_SupportSize);
244 m_SupportIndex.Fill(0);
245 m_SupportRegion.SetIndex(m_SupportIndex);
246 for ( i = 0; i < SpaceDimension ; i++ )
247 m_ZeroVector[i]=itk::NumericTraits<JacobianValueType>::Zero;
250 m_FirstRegion.SetSize(m_NullSize);
251 m_FirstRegion.SetIndex(m_NullIndex);
252 m_SecondRegion.SetSize(m_NullSize);
253 m_SecondRegion.SetIndex(m_NullIndex);
254 m_ThirdRegion.SetSize(m_NullSize);
255 m_ThirdRegion.SetIndex(m_NullIndex);
257 m_BCValues.resize(0);
258 m_BCRegions.resize(0);
260 m_BC2Values.resize(0);
261 m_BC2Regions.resize(0);
263 m_BC3Values.resize(0);
264 m_BC3Regions.resize(0);
268 for ( i = 0; i < SpaceDimension ; i++ )
270 m_FirstIterator[i]= IteratorType( m_JacobianImage[i], m_FirstRegion);
271 m_SecondIterator[i]= IteratorType( m_JacobianImage[i], m_SecondRegion);
272 m_ThirdIterator[i]= IteratorType( m_JacobianImage[i], m_ThirdRegion);
273 m_BCIterators[i].resize(0);
274 m_BC2Iterators[i].resize(0);
275 m_BC3Iterators[i].resize(0);
284 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
285 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
286 ::~ShapedBLUTSpatioTemporalDeformableTransform()
292 // JV set Spline Order
293 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
295 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
296 ::SetSplineOrder(const unsigned int & splineOrder)
298 SizeType splineOrders;
299 for (unsigned int i=0;i<InputDimension;i++)splineOrders[i]=splineOrder;
301 this->SetSplineOrders(splineOrders);
305 // JV set Spline Orders
306 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
308 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
309 ::SetSplineOrders(const SizeType & splineOrders)
311 if(m_SplineOrders!=splineOrders)
313 m_SplineOrders=splineOrders;
315 //update the interpolation function
316 m_VectorInterpolator->SetSplineOrders(m_SplineOrders);
318 //update the varaibles for computing interpolation
319 for (unsigned int i=0; i <InputDimension;i++)
321 m_SupportSize[i] = m_SplineOrders[i]+1;
322 m_Offset[i] = m_SplineOrders[i] / 2;
324 if ( m_SplineOrders[i] % 2 )
326 m_SplineOrderOdd[i] = true;
330 m_SplineOrderOdd[i] = false;
334 //SupportSize is updated!, update regions
335 //JV initialize some variables for jacobian calculation
336 m_SupportRegion.SetSize(m_SupportSize);
337 m_SupportIndex.Fill(0);
338 m_SupportRegion.SetIndex(m_SupportIndex);
340 // Update the WeightRatios
341 for (unsigned int i=0; i<4; i++)
342 for (unsigned int j=0; j<4; j++)
343 m_WeightRatio[i][j]=m_Weights[m_SplineOrders[InputDimension-1]][i]/ m_Weights[m_SplineOrders[InputDimension-1]][j];
350 // JV set sampling factor
351 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
353 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
354 ::SetLUTSamplingFactor( const int & samplingFactor)
356 SizeType samplingFactors;
357 for (unsigned int i=0; i<NInputDimensions; i++)
358 samplingFactors[i]=samplingFactor;
360 //update the interpolation function
361 this->SetLUTSamplingFactors(samplingFactors);
365 // JV set sampling factors
366 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
368 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
369 ::SetLUTSamplingFactors( const SizeType & samplingFactors)
371 if(m_LUTSamplingFactors!=samplingFactors)
373 for (unsigned int i=0; i<NInputDimensions; i++)
374 m_LUTSamplingFactors[i]=samplingFactors[i];
376 //update the interpolation function
377 m_VectorInterpolator->SetLUTSamplingFactors(m_LUTSamplingFactors);
384 // Get the number of parameters
385 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
386 #if ITK_VERSION_MAJOR >= 4
387 typename ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>::NumberOfParametersType
391 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
392 ::GetNumberOfParameters(void) const
395 // The number of parameters equal SpaceDimension * number of
396 // of pixels in the grid region.
397 return ( static_cast<unsigned int>( SpaceDimension ) *
398 static_cast<unsigned int>( m_GridRegion.GetNumberOfPixels() ) );
403 // Get the padded number of parameters
404 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
406 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
407 ::GetPaddedNumberOfParameters(void) const
410 // The number of parameters equal SpaceDimension * number of
411 // of pixels in the grid region.
412 return ( static_cast<unsigned int>( SpaceDimension ) *
413 static_cast<unsigned int>( m_PaddedGridRegion.GetNumberOfPixels() ) );
419 // Get the number of parameters per dimension
420 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
422 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
423 ::GetNumberOfParametersPerDimension(void) const
425 // The number of parameters per dimension equal number of
426 // of pixels in the grid region.
427 return ( static_cast<unsigned int>( m_GridRegion.GetNumberOfPixels() ) );
432 // Set the grid region
433 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
435 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
436 ::SetGridRegion( const RegionType & region )
438 if ( m_GridRegion != region )
440 m_GridRegion = region;
441 m_PaddedGridRegion=region;
443 // JV set the padded region
444 typename CoefficientImageType::RegionType::SizeType paddedSize= region.GetSize();
446 //JV size dependes on shape
447 switch (m_TransformShape)
451 paddedSize[InputDimension-1]+=4;
455 paddedSize[InputDimension-1]+=4;
459 paddedSize[InputDimension-1]+=3;
463 paddedSize[InputDimension-1]+=2;
467 paddedSize[InputDimension-1]+=3;
470 paddedSize[InputDimension-1]=+1;
472 m_PaddedGridRegion.SetSize(paddedSize);
474 // Set regions for each coefficient and jacobian image
475 m_WrappedImage->SetRegions( m_GridRegion );
476 m_PaddedCoefficientImage->SetRegions( m_PaddedGridRegion );
477 m_PaddedCoefficientImage->Allocate();
478 for (unsigned int j=0; j <OutputDimension;j++)
480 m_JacobianImage[j]->SetRegions( m_GridRegion );
483 // JV used the padded version for the valid region
484 // Set the valid region
485 // If the grid spans the interval [start, last].
486 // The valid interval for evaluation is [start+offset, last-offset]
487 // when spline order is even.
488 // The valid interval for evaluation is [start+offset, last-offset)
489 // when spline order is odd.
490 // Where offset = vcl_floor(spline / 2 ).
491 // Note that the last pixel is not included in the valid region
492 // with odd spline orders.
493 typename RegionType::SizeType size = m_PaddedGridRegion.GetSize();
494 typename RegionType::IndexType index = m_PaddedGridRegion.GetIndex();
495 for ( unsigned int j = 0; j < NInputDimensions; j++ )
498 static_cast< typename RegionType::IndexValueType >( m_Offset[j] );
500 static_cast< typename RegionType::SizeValueType> ( 2 * m_Offset[j] );
501 m_ValidRegionLast[j] = index[j] +
502 static_cast< typename RegionType::IndexValueType >( size[j] ) - 1;
504 m_ValidRegion.SetSize( size );
505 m_ValidRegion.SetIndex( index );
507 // If we are using the default parameters, update their size and set to identity.
508 // Input parameters point to internal buffer => using default parameters.
509 if (m_InputParametersPointer == &m_InternalParametersBuffer)
511 // Check if we need to resize the default parameter buffer.
512 if ( m_InternalParametersBuffer.GetSize() != this->GetNumberOfParameters() )
514 m_InternalParametersBuffer.SetSize( this->GetNumberOfParameters() );
515 // Fill with zeros for identity.
516 m_InternalParametersBuffer.Fill( 0 );
525 // Set the grid spacing
526 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
528 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
529 ::SetGridSpacing( const SpacingType & spacing )
531 if ( m_GridSpacing != spacing )
533 m_GridSpacing = spacing;
535 // Set spacing for each coefficient and jacobian image
536 m_WrappedImage->SetSpacing( m_GridSpacing.GetDataPointer() );
537 m_PaddedCoefficientImage->SetSpacing( m_GridSpacing.GetDataPointer() );
538 for (unsigned int j=0; j <OutputDimension; j++) m_JacobianImage[j]->SetSpacing( m_GridSpacing.GetDataPointer() );
542 for( unsigned int i=0; i<OutputDimension; i++)
544 scale[i][i] = m_GridSpacing[i];
547 m_IndexToPoint = m_GridDirection * scale;
548 m_PointToIndex = m_IndexToPoint.GetInverse();
555 // Set the grid direction
556 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
558 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
559 ::SetGridDirection( const DirectionType & direction )
561 if ( m_GridDirection != direction )
563 m_GridDirection = direction;
565 // Set direction for each coefficient and jacobian image
566 m_WrappedImage->SetDirection( m_GridDirection );
567 m_PaddedCoefficientImage->SetDirection( m_GridDirection );
568 for (unsigned int j=0; j <OutputDimension; j++) m_JacobianImage[j]->SetDirection( m_GridDirection );
572 for( unsigned int i=0; i<OutputDimension; i++)
574 scale[i][i] = m_GridSpacing[i];
577 m_IndexToPoint = m_GridDirection * scale;
578 m_PointToIndex = m_IndexToPoint.GetInverse();
586 // Set the grid origin
587 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
589 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
590 ::SetGridOrigin( const OriginType& origin )
592 if( m_GridOrigin!=origin)
594 m_GridOrigin = origin;
596 // JV The origin depends on the shape
597 switch (m_TransformShape)
601 m_PaddedGridOrigin=origin;
602 m_PaddedGridOrigin[InputDimension-1]=origin[InputDimension-1]-2* m_GridSpacing[InputDimension-1];
606 m_PaddedGridOrigin=origin;
607 m_PaddedGridOrigin[InputDimension-1]=origin[InputDimension-1]-1* m_GridSpacing[InputDimension-1];
611 m_PaddedGridOrigin=origin;
612 m_PaddedGridOrigin[InputDimension-1]=origin[InputDimension-1]-1* m_GridSpacing[InputDimension-1];
616 m_PaddedGridOrigin=origin;
620 m_PaddedGridOrigin=origin;
621 m_PaddedGridOrigin[InputDimension-1]=origin[InputDimension-1]-1* m_GridSpacing[InputDimension-1];
624 m_PaddedGridOrigin=origin;
627 // Set origin for each coefficient and jacobianimage
628 m_WrappedImage->SetOrigin( m_GridOrigin.GetDataPointer() );
629 m_PaddedCoefficientImage->SetOrigin( m_PaddedGridOrigin.GetDataPointer() );
630 for (unsigned int j=0; j <OutputDimension; j++) m_JacobianImage[j]->SetOrigin( m_GridOrigin.GetDataPointer() );
637 // Set the parameters
638 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
640 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
643 if( m_InputParametersPointer )
645 ParametersType * parameters =
646 const_cast<ParametersType *>( m_InputParametersPointer );
647 parameters->Fill( 0.0 );
652 itkExceptionMacro( << "Input parameters for the spline haven't been set ! "
653 << "Set them using the SetParameters or SetCoefficientImage method first." );
658 // Set the parameters
659 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
661 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
662 ::SetParameters( const ParametersType & parameters )
665 // Check if the number of parameters match the
666 // Expected number of parameters
667 if ( parameters.Size() != this->GetNumberOfParameters() )
669 itkExceptionMacro(<<"Mismatched between parameters size "
671 << " and region size "
672 << m_GridRegion.GetNumberOfPixels() );
675 // Clean up buffered parameters
676 m_InternalParametersBuffer = ParametersType( 0 );
678 // Keep a reference to the input parameters
679 m_InputParametersPointer = ¶meters;
681 // Wrap flat array as images of coefficients
682 this->WrapAsImages();
684 //JV Set padded input to vector interpolator
685 m_VectorInterpolator->SetInputImage(this->GetPaddedCoefficientImage());
687 // Modified is always called since we just have a pointer to the
688 // parameters and cannot know if the parameters have changed.
693 // Set the Fixed Parameters
694 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
696 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
697 ::SetFixedParameters( const ParametersType & parameters )
700 // JV number should be exact, no defaults for spacing
701 if ( parameters.Size() != NInputDimensions * (5 + NInputDimensions)+3 )
703 itkExceptionMacro(<< "Mismatched between parameters size "
705 << " and number of fixed parameters "
706 << NInputDimensions * (5 + NInputDimensions)+3 );
708 /*********************************************************
709 Fixed Parameters store the following information:
714 // JV we add the splineOrders, LUTsamplingfactor, mask pointer and bulktransform pointer
722 The size of these is equal to the NInputDimensions
723 *********************************************************/
725 /** Set the Grid Parameters */
727 for (unsigned int i=0; i<NInputDimensions; i++)
729 gridSize[i] = static_cast<int> (parameters[i]);
731 RegionType bsplineRegion;
732 bsplineRegion.SetSize( gridSize );
734 /** Set the Origin Parameters */
736 for (unsigned int i=0; i<NInputDimensions; i++)
738 origin[i] = parameters[NInputDimensions+i];
741 /** Set the Spacing Parameters */
743 for (unsigned int i=0; i<NInputDimensions; i++)
745 spacing[i] = parameters[2*NInputDimensions+i];
748 /** Set the Direction Parameters */
749 DirectionType direction;
750 for (unsigned int di=0; di<NInputDimensions; di++)
752 for (unsigned int dj=0; dj<NInputDimensions; dj++)
754 direction[di][dj] = parameters[3*NInputDimensions+(di*NInputDimensions+dj)];
758 //JV add the spline orders
759 SizeType splineOrders;
760 for (unsigned int i=0; i<NInputDimensions; i++)
762 splineOrders[i]=parameters[(3+NInputDimensions)*NInputDimensions+i];
765 //JV add the LUT sampling factor
766 SizeType samplingFactors;
767 for (unsigned int i=0; i<NInputDimensions; i++)
769 samplingFactors[i]=parameters[(4+NInputDimensions)*NInputDimensions+i];
772 //JV add the MaskPointer
773 m_Mask=((MaskPointer)( (size_t)parameters[(5+NInputDimensions)*NInputDimensions]));
775 //JV add the MaskPointer
776 m_BulkTransform=((BulkTransformPointer)( (size_t)parameters[(5+NInputDimensions)*NInputDimensions+1]));
778 //JV add the TransformShape
779 m_TransformShape=((unsigned int)parameters[(5+NInputDimensions)*NInputDimensions+2]);
782 this->SetSplineOrders( splineOrders );
783 this->SetGridSpacing( spacing );
784 this->SetGridDirection( direction );
785 this->SetGridOrigin( origin );
786 this->SetGridRegion( bsplineRegion );
787 this->SetLUTSamplingFactors( samplingFactors );
792 // Wrap flat parameters as images
793 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
795 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
798 //JV Wrap parameter array in vectorial image, changed parameter order: A1x A1y A1z, A2x ....
799 PixelType * dataPointer =reinterpret_cast<PixelType *>( const_cast<double *>(m_InputParametersPointer->data_block() )) ;
800 unsigned int numberOfPixels = m_GridRegion.GetNumberOfPixels();
802 m_WrappedImage->GetPixelContainer()->SetImportPointer( dataPointer,numberOfPixels);//InputDimension
803 m_CoefficientImage = m_WrappedImage;
805 //=====================================
806 //JV Create padded structure adding BC
807 //=====================================
808 PadCoefficientImage();
810 //=====================================
811 //JV Wrap jacobian into OutputDimension X Vectorial images
812 //=====================================
813 #if ITK_VERSION_MAJOR >= 4
814 this->m_SharedDataBSplineJacobian.set_size( OutputDimension, this->GetNumberOfParameters() );
816 this->m_Jacobian.set_size( OutputDimension, this->GetNumberOfParameters() );
819 // Use memset to set the memory
820 // JV four rows of three comps of parameters
821 #if ITK_VERSION_MAJOR >= 4
822 JacobianPixelType * jacobianDataPointer = reinterpret_cast<JacobianPixelType *>(this->m_SharedDataBSplineJacobian.data_block());
824 JacobianPixelType * jacobianDataPointer = reinterpret_cast<JacobianPixelType *>(this->m_Jacobian.data_block());
826 memset(jacobianDataPointer, 0, OutputDimension*numberOfPixels*sizeof(JacobianPixelType));
828 for (unsigned int j=0; j<OutputDimension; j++)
830 m_JacobianImage[j]->GetPixelContainer()->
831 SetImportPointer( jacobianDataPointer, numberOfPixels );
832 jacobianDataPointer += numberOfPixels;
835 // Reset the J parameters
836 m_LastJacobianIndex = m_ValidRegion.GetIndex();
837 m_FirstRegion.SetSize(m_NullSize);
838 m_SecondRegion.SetSize(m_NullSize);
839 for ( unsigned int j = 0; j < SpaceDimension ; j++ )
841 m_FirstIterator[j]= IteratorType( m_JacobianImage[j], m_FirstRegion);
842 m_SecondIterator[j]= IteratorType( m_JacobianImage[j], m_SecondRegion);
845 m_BCValues.resize(0);
846 m_BCRegions.resize(0);
848 m_BC2Values.resize(0);
849 m_BC2Regions.resize(0);
851 m_BC3Values.resize(0);
852 m_BC3Regions.resize(0);
858 // Set the parameters by value
859 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
861 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
862 ::SetParametersByValue( const ParametersType & parameters )
865 // Check if the number of parameters match the
866 // Expected number of parameters
867 if ( parameters.Size() != this->GetNumberOfParameters() )
869 itkExceptionMacro(<<"Mismatched between parameters size "
871 << " and region size "
872 << m_GridRegion.GetNumberOfPixels() );
876 m_InternalParametersBuffer = parameters;
877 m_InputParametersPointer = &m_InternalParametersBuffer;
879 // Wrap flat array as images of coefficients
880 this->WrapAsImages();
882 //JV Set padded input to vector interpolator
883 m_VectorInterpolator->SetInputImage(this->GetPaddedCoefficientImage());
885 // Modified is always called since we just have a pointer to the
886 // Parameters and cannot know if the parameters have changed.
891 // Get the parameters
892 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
894 typename ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
896 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
897 ::GetParameters( void ) const
899 /** NOTE: For efficiency, this class does not keep a copy of the parameters -
900 * it just keeps pointer to input parameters.
902 if (NULL == m_InputParametersPointer)
904 itkExceptionMacro( <<"Cannot GetParameters() because m_InputParametersPointer is NULL. Perhaps SetCoefficientImage() has been called causing the NULL pointer." );
907 return (*m_InputParametersPointer);
911 // Get the parameters
912 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
914 typename ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
916 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
917 ::GetFixedParameters( void ) const
919 RegionType resRegion = this->GetGridRegion( );
921 for (unsigned int i=0; i<NInputDimensions; i++)
923 this->m_FixedParameters[i] = (resRegion.GetSize())[i];
925 for (unsigned int i=0; i<NInputDimensions; i++)
927 this->m_FixedParameters[NInputDimensions+i] = (this->GetGridOrigin())[i];
929 for (unsigned int i=0; i<NInputDimensions; i++)
931 this->m_FixedParameters[2*NInputDimensions+i] = (this->GetGridSpacing())[i];
933 for (unsigned int di=0; di<NInputDimensions; di++)
935 for (unsigned int dj=0; dj<NInputDimensions; dj++)
937 this->m_FixedParameters[3*NInputDimensions+(di*NInputDimensions+dj)] = (this->GetGridDirection())[di][dj];
941 //JV add splineOrders
942 for (unsigned int i=0; i<NInputDimensions; i++)
944 this->m_FixedParameters[(3+NInputDimensions)*NInputDimensions+i] = (this->GetSplineOrders())[i];
947 //JV add LUTsamplingFactor
948 for (unsigned int i=0; i<NInputDimensions; i++)
950 this->m_FixedParameters[(4+NInputDimensions)*NInputDimensions+i] = (this->GetLUTSamplingFactors())[i];
954 this->m_FixedParameters[(5+NInputDimensions)*NInputDimensions]=(double)((size_t) m_Mask);
956 //JV add the bulktransform pointer
957 this->m_FixedParameters[(5+NInputDimensions)*NInputDimensions+1]=(double)((size_t) m_BulkTransform);
959 //JV add the transform shape
960 this->m_FixedParameters[(5+NInputDimensions)*NInputDimensions+2]=(double)( m_TransformShape);
962 return (this->m_FixedParameters);
966 // Set the B-Spline coefficients using input images
967 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
969 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
970 ::SetCoefficientImage( CoefficientImagePointer image )
972 this->SetGridSpacing( image->GetSpacing() );
973 this->SetGridOrigin( image->GetOrigin() );
974 this->SetGridDirection( image->GetDirection() );
975 this->SetGridRegion( image->GetBufferedRegion() );
976 m_CoefficientImage = image;
979 // m_WrappedImage=m_CoefficientImage;
981 // Update the interpolator
982 this->PadCoefficientImage();
983 m_VectorInterpolator->SetInputImage(this->GetPaddedCoefficientImage());
985 // Clean up buffered parameters
986 m_InternalParametersBuffer = ParametersType( 0 );
987 m_InputParametersPointer = NULL;
992 // // Set the B-Spline coefficients using input images
993 // template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
995 // ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
996 // ::SetPaddedCoefficientImage( CoefficientImagePointer image )
998 // //JV modify the region
999 // typename CoefficientImageType::RegionType region=image->GetBufferedRegion();
1000 // typename CoefficientImageType::RegionType::SizeType size=region.GetSize();
1001 // size[InputDimension-1]-=2;
1002 // region.SetSize(size);
1005 // this->SetGridRegion( region );
1006 // this->SetGridSpacing( image->GetSpacing() );
1007 // this->SetGridDirection( image->GetDirection() );
1008 // this->SetGridOrigin( image->GetOrigin() );
1009 // m_PaddedCoefficientImage = image;
1010 // this->ExtractCoefficientImage();
1011 // m_VectorInterpolator->SetInputImage(this->GetPaddedCoefficientImage());
1013 // // Clean up buffered parameters
1014 // m_InternalParametersBuffer = ParametersType( 0 );
1015 // m_InputParametersPointer = NULL;
1019 // Set the B-Spline coefficients using input images
1020 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
1021 typename ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>::CoefficientImageType::Pointer
1022 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
1023 ::ExtractTemporalRow(const typename CoefficientImageType::Pointer& coefficientImage, unsigned int temporalIndex)
1026 typename CoefficientImageType::RegionType sourceRegion=coefficientImage->GetLargestPossibleRegion();
1027 sourceRegion.SetSize(InputDimension-1, 1);
1028 sourceRegion.SetIndex(InputDimension-1, temporalIndex);
1031 typedef clitk::ExtractImageFilter<CoefficientImageType, CoefficientImageType> ExtractImageFilterType;
1032 typename ExtractImageFilterType::Pointer extract=ExtractImageFilterType::New();
1033 extract->SetInput(coefficientImage);
1034 extract->SetExtractionRegion(sourceRegion);
1036 typename CoefficientImageType::Pointer row= extract->GetOutput();
1038 // Set index to zero
1039 sourceRegion.SetIndex(InputDimension-1, 0);
1040 row->SetRegions(sourceRegion);
1045 // Set the B-Spline coefficients using input images
1046 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
1048 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
1049 ::PadCoefficientImage(void)
1052 // Define paste, extract and combine filters
1053 typedef itk::PasteImageFilter<CoefficientImageType, CoefficientImageType, CoefficientImageType> PasteImageFilterType;
1054 typedef clitk::ExtractImageFilter<CoefficientImageType, CoefficientImageType> ExtractImageFilterType;
1055 typedef clitk::LinearCombinationImageFilter<CoefficientImageType, CoefficientImageType> LinearCombinationFilterType;
1058 typename CoefficientImageType::RegionType sourceRegion=m_PaddedCoefficientImage->GetLargestPossibleRegion();
1059 typename CoefficientImageType::RegionType destinationRegion=m_PaddedCoefficientImage->GetLargestPossibleRegion();
1060 typename CoefficientImageType::RegionType::SizeType sourceSize=sourceRegion.GetSize();
1061 typename CoefficientImageType::RegionType::SizeType destinationSize=destinationRegion.GetSize();
1062 typename CoefficientImageType::IndexType sourceIndex=sourceRegion.GetIndex();
1063 typename CoefficientImageType::IndexType destinationIndex=destinationRegion.GetIndex();
1065 // JV Padding depends on the shape
1066 switch (m_TransformShape)
1071 2: rabbit 4 CP 3 DOF
1072 3: rabbit 5 CP 4 DOF
1073 4: sputnik 4 CP 4 DOF
1074 5: sputnik 5 CP 5 DOF
1075 6: diamond 6 CP 5 DOF
1076 7: diamond 7 CP 6 DOF
1081 // ----------------------------------------------------------------------
1082 // The egg with 4 internal CP (starting from inhale)
1083 // Periodic, constrained to zero at the reference
1084 // at position 3 and
1085 // Coeff row BC1 BC2 0 BC3 1 2 BC4
1086 // PaddedCoeff R: 0 1 2 3 4 5 6
1089 // BC3= -weights[2]/weights[0] ( R2+R4 )
1091 // ---------------------------------------------------------------------
1093 //---------------------------------
1094 // 1. First two temporal rows are identical: paste 1-2 to 0-1
1095 typename PasteImageFilterType::Pointer paster0=PasteImageFilterType::New();
1096 paster0->SetDestinationImage(m_PaddedCoefficientImage);
1097 paster0->SetSourceImage(m_CoefficientImage);
1098 sourceRegion.SetIndex(NInputDimensions-1,1);
1099 sourceRegion.SetSize(NInputDimensions-1,2);
1100 destinationIndex[InputDimension-1]=0;
1101 paster0->SetDestinationIndex(destinationIndex);
1102 paster0->SetSourceRegion(sourceRegion);
1104 //---------------------------------
1105 // 1. Next temporal row = paste 0 to 2
1106 typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1107 typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1108 paster1->SetDestinationImage(paster0->GetOutput());
1109 paster1->SetSourceImage(row0);
1110 destinationIndex[InputDimension-1]=2;
1111 paster1->SetDestinationIndex(destinationIndex);
1112 paster1->SetSourceRegion(row0->GetLargestPossibleRegion());
1114 //---------------------------------
1115 // 2. Middle row at index 3=BC
1116 // Extract row at index 1, 2 (of coeff image)
1117 typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1120 typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1121 combine1->SetInput(0,row0);
1122 combine1->SetInput(1,row1);
1123 combine1->SetA(-m_WeightRatio[2][0]);
1124 combine1->SetB(-m_WeightRatio[2][0]);
1126 typename CoefficientImageType::Pointer bc3Row=combine1->GetOutput();
1128 // Paste middleRow at index 3 (padded coeff)
1129 typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1130 paster2->SetDestinationImage(paster1->GetOutput());
1131 paster2->SetSourceImage(bc3Row);
1132 destinationIndex[InputDimension-1]=3;
1133 paster2->SetDestinationIndex(destinationIndex);
1134 paster2->SetSourceRegion(bc3Row->GetLargestPossibleRegion());
1136 //---------------------------------
1137 // 3. Next two temporal rows identical: paste 1,2 to 4,5
1138 typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1139 paster3->SetDestinationImage(paster2->GetOutput());
1140 paster3->SetSourceImage(m_CoefficientImage);
1141 sourceRegion.SetIndex(InputDimension-1,1);
1142 sourceRegion.SetSize(InputDimension-1,2);
1143 destinationIndex[InputDimension-1]=4;
1144 paster3->SetDestinationIndex(destinationIndex);
1145 paster3->SetSourceRegion(sourceRegion);
1147 // ---------------------------------
1148 // 4. Row at index 6=BC4= R2
1149 // Paste BC3 row at index 5
1150 typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
1151 paster4->SetDestinationImage(paster3->GetOutput());
1152 paster4->SetSourceImage(row0);
1153 destinationIndex[InputDimension-1]=6;
1154 paster4->SetDestinationIndex(destinationIndex);
1155 paster4->SetSourceRegion(row0->GetLargestPossibleRegion());
1157 // Update the chain!
1159 m_PaddedCoefficientImage= paster4->GetOutput();
1166 // ----------------------------------------------------------------------
1167 // The egg with 5 internal CP (starting from inhale)
1168 // Periodic, constrained to zero at the reference
1169 // at position 3 and
1170 // Coeff row BC1 BC2 0 BC3 1 2 3 BC4
1171 // PaddedCoeff R: 0 1 2 3 4 5 6 7
1174 // BC3= -weights[3]/weights[1] ( R2+R5 ) - R4
1176 // ---------------------------------------------------------------------
1177 //---------------------------------
1178 // 1. First two temporal rows are identical: paste 2-3 to 0-1
1179 typename PasteImageFilterType::Pointer paster0=PasteImageFilterType::New();
1180 paster0->SetDestinationImage(m_PaddedCoefficientImage);
1181 paster0->SetSourceImage(m_CoefficientImage);
1182 sourceRegion.SetIndex(NInputDimensions-1,2);
1183 sourceRegion.SetSize(NInputDimensions-1,2);
1184 destinationIndex[InputDimension-1]=0;
1185 paster0->SetDestinationIndex(destinationIndex);
1186 paster0->SetSourceRegion(sourceRegion);
1188 //---------------------------------
1189 // 1. Next temporal row = paste 0 to 2
1190 typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1191 typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1192 paster1->SetDestinationImage(paster0->GetOutput());
1193 paster1->SetSourceImage(row0);
1194 destinationIndex[InputDimension-1]=2;
1195 paster1->SetDestinationIndex(destinationIndex);
1196 paster1->SetSourceRegion(row0->GetLargestPossibleRegion());
1198 //---------------------------------
1199 // 2. Middle row at index 3=BC
1200 // Extract row at index 1, 2 (of coeff image)
1201 typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1202 typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
1205 typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1206 combine1->SetInput(0,row0);
1207 combine1->SetInput(1,row2);
1210 // combine1->Update();
1211 typename LinearCombinationFilterType::Pointer combine2=LinearCombinationFilterType::New();
1212 combine2->SetInput(0,row1);
1213 combine2->SetInput(1,combine1->GetOutput());
1214 combine2->SetA(-1.);
1215 combine2->SetB(-m_WeightRatio[3][1]);
1217 typename CoefficientImageType::Pointer bc3Row=combine2->GetOutput();
1219 // Paste middleRow at index 3 (padded coeff)
1220 typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1221 paster2->SetDestinationImage(paster1->GetOutput());
1222 paster2->SetSourceImage(bc3Row);
1223 destinationIndex[InputDimension-1]=3;
1224 paster2->SetDestinationIndex(destinationIndex);
1225 paster2->SetSourceRegion(bc3Row->GetLargestPossibleRegion());
1227 //---------------------------------
1228 // 3. Next three temporal rows identical: paste 1,2,3 to 4,5,6
1229 typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1230 paster3->SetDestinationImage(paster2->GetOutput());
1231 paster3->SetSourceImage(m_CoefficientImage);
1232 sourceRegion.SetIndex(InputDimension-1,1);
1233 sourceRegion.SetSize(InputDimension-1,3);
1234 destinationIndex[InputDimension-1]=4;
1235 paster3->SetDestinationIndex(destinationIndex);
1236 paster3->SetSourceRegion(sourceRegion);
1238 // ---------------------------------
1239 // 4. Row at index 7=BC4= R2
1240 // Paste BC3 row at index 5
1241 typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
1242 paster4->SetDestinationImage(paster3->GetOutput());
1243 paster4->SetSourceImage(row0);
1244 destinationIndex[InputDimension-1]=7;
1245 paster4->SetDestinationIndex(destinationIndex);
1246 paster4->SetSourceRegion(row0->GetLargestPossibleRegion());
1248 // Update the chain!
1250 m_PaddedCoefficientImage= paster4->GetOutput();
1256 // // ----------------------------------------------------------------------
1257 // // The egg with 5 internal CP:
1258 // // Periodic, C2 smooth everywhere and constrained to zero at the reference
1259 // // Coeff row R5 BC1 0 1 2 3 BC2 R2
1260 // // PaddedCoeff R: 0 1 2 3 4 5 6 7
1261 // // BC1= -weights[2]/weights[0] ( R2+R5)
1263 // // ---------------------------------------------------------------------
1265 // // Extract rows with index 0 and 3
1266 // typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1267 // typename CoefficientImageType::Pointer row3=this->ExtractTemporalRow(m_CoefficientImage, 3);
1269 // // Paste the first row
1270 // typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1271 // destinationIndex.Fill(0);
1272 // paster1->SetDestinationIndex(destinationIndex);
1273 // paster1->SetSourceRegion(row3->GetLargestPossibleRegion());
1274 // paster1->SetSourceImage(row3);
1275 // paster1->SetDestinationImage(m_PaddedCoefficientImage);
1277 // // Linearly Combine rows for BC1 and BC2
1278 // typedef clitk::LinearCombinationImageFilter<CoefficientImageType, CoefficientImageType> LinearCombinationFilterType;
1279 // typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1280 // combine1->SetFirstInput(row0);
1281 // combine1->SetSecondInput(row3);
1282 // combine1->SetA(-m_WeightRatio[2][0]);
1283 // combine1->SetB(-m_WeightRatio[2][0]);
1284 // combine1->Update();
1285 // typename CoefficientImageType::Pointer bcRow=combine1->GetOutput();
1287 // // Paste the second row
1288 // typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1289 // destinationIndex[InputDimension-1]=1;
1290 // paster2->SetDestinationIndex(destinationIndex);
1291 // paster2->SetSourceRegion(bcRow->GetLargestPossibleRegion());
1292 // paster2->SetSourceImage(bcRow);
1293 // paster2->SetDestinationImage(paster1->GetOutput());
1295 // // Paste the coefficientImage
1296 // typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1297 // destinationIndex[InputDimension-1]=2;
1298 // paster3->SetDestinationIndex(destinationIndex);
1299 // paster3->SetSourceRegion(m_CoefficientImage->GetLargestPossibleRegion());
1300 // paster3->SetSourceImage(m_CoefficientImage);
1301 // paster3->SetDestinationImage(paster2->GetOutput());
1303 // // Paste the last two rows
1304 // typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
1305 // destinationIndex[InputDimension-1]=6;
1306 // paster4->SetDestinationIndex(destinationIndex);
1307 // paster4->SetSourceRegion(bcRow->GetLargestPossibleRegion());
1308 // paster4->SetSourceImage(bcRow);
1309 // paster4->SetDestinationImage(paster3->GetOutput());
1311 // typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
1312 // destinationIndex[InputDimension-1]=7;
1313 // paster5->SetDestinationIndex(destinationIndex);
1314 // paster5->SetSourceRegion(row0->GetLargestPossibleRegion());
1315 // paster5->SetSourceImage(row0);
1316 // paster5->SetDestinationImage(paster4->GetOutput());
1318 // // Update the chain!
1319 // paster5->Update();
1320 // m_PaddedCoefficientImage= paster5->GetOutput();
1327 // ----------------------------------------------------------------------
1328 // The rabbit with 4 internal CP
1329 // Periodic, constrained to zero at the reference
1330 // at position 3 and the extremes fixed with anit-symm bc
1331 // Coeff row BC1 0 1 BC2 2 BC3 BC4
1332 // PaddedCoeff R: 0 1 2 3 4 5 6
1334 // BC2= -weights[2]/weights[0] ( R2+R4 )
1337 // ---------------------------------------------------------------------
1339 // ---------------------------------
1340 // 0. First Row =BC1
1341 typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1342 typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1343 typename LinearCombinationFilterType::Pointer combine0=LinearCombinationFilterType::New();
1344 combine0->SetInput(0,row0);
1345 combine0->SetInput(1,row1);
1347 combine0->SetB(-1.);
1349 typename CoefficientImageType::Pointer bc1Row=combine0->GetOutput();
1351 typename PasteImageFilterType::Pointer paster0=PasteImageFilterType::New();
1352 paster0->SetDestinationImage(m_PaddedCoefficientImage);
1353 paster0->SetSourceImage(bc1Row);
1354 destinationIndex[InputDimension-1]=0;
1355 paster0->SetDestinationIndex(destinationIndex);
1356 paster0->SetSourceRegion(bc1Row->GetLargestPossibleRegion());
1358 //---------------------------------
1359 // 1. Next two temporal rows are identical: paste 0-1 to 1-2
1360 typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1361 paster1->SetDestinationImage(paster0->GetOutput());
1362 paster1->SetSourceImage(m_CoefficientImage);
1363 sourceRegion.SetIndex(NInputDimensions-1,0);
1364 destinationIndex[InputDimension-1]=1;
1365 sourceRegion.SetSize(NInputDimensions-1,2);
1366 paster1->SetDestinationIndex(destinationIndex);
1367 paster1->SetSourceRegion(sourceRegion);
1369 //---------------------------------
1370 // 2. Middle row at index 3=BC
1371 // Extract row at index 1, 2 (of coeff image)
1372 typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
1375 typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1376 combine1->SetInput(0,row1);
1377 combine1->SetInput(1,row2);
1378 combine1->SetA(-m_WeightRatio[2][0]);
1379 combine1->SetB(-m_WeightRatio[2][0]);
1381 typename CoefficientImageType::Pointer bc2Row=combine1->GetOutput();
1383 // Paste middleRow at index 3 (padded coeff)
1384 typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1385 paster2->SetDestinationImage(paster1->GetOutput());
1386 paster2->SetSourceImage(bc2Row);
1387 destinationIndex[InputDimension-1]=3;
1388 paster2->SetDestinationIndex(destinationIndex);
1389 paster2->SetSourceRegion(bc2Row->GetLargestPossibleRegion());
1391 //---------------------------------
1392 // 3. Next temporal row is identical: paste 2 to 4
1393 typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1394 paster3->SetDestinationImage(paster2->GetOutput());
1395 paster3->SetSourceImage(row2);
1396 destinationIndex[InputDimension-1]=4;
1397 paster3->SetDestinationIndex(destinationIndex);
1398 paster3->SetSourceRegion(row2->GetLargestPossibleRegion());
1400 // ---------------------------------
1401 // 4. Row at index 5=BC (paddedcoeff image) R1
1402 // Paste BC3 row at index 5
1403 typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
1404 paster4->SetDestinationImage(paster3->GetOutput());
1405 paster4->SetSourceImage(row0);
1406 destinationIndex[InputDimension-1]=5;
1407 paster4->SetDestinationIndex(destinationIndex);
1408 paster4->SetSourceRegion(row0->GetLargestPossibleRegion());
1410 // ---------------------------------
1411 // 5. Paste BC4 row at index 6
1412 typename LinearCombinationFilterType::Pointer combine3=LinearCombinationFilterType::New();
1413 combine3->SetInput(0,row0);
1414 combine3->SetInput(1,row2);
1416 combine3->SetB(-1.);
1418 typename CoefficientImageType::Pointer bc4Row=combine3->GetOutput();
1419 typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
1420 paster5->SetDestinationImage(paster4->GetOutput());
1421 paster5->SetSourceImage(bc4Row);
1422 destinationIndex[InputDimension-1]=6;
1423 paster5->SetDestinationIndex(destinationIndex);
1424 paster5->SetSourceRegion(bc4Row->GetLargestPossibleRegion());
1426 // Update the chain!
1428 m_PaddedCoefficientImage= paster5->GetOutput();
1435 // ----------------------------------------------------------------------
1436 // The rabbit with 5 internal CP
1437 // Periodic, constrained to zero at the reference at position 3.5
1438 // and the extremes fixed with anti-symmetrical BC
1439 // Coeff row BC1 0 1 BC2 2 3 BC3 BC4
1440 // PaddedCoeff R: 0 1 2 3 4 5 6 7
1442 // BC2= -weights[3]/weights[1] ( R2+R5 ) - R4
1445 // ---------------------------------------------------------------------
1447 // ---------------------------------
1448 // 0. First Row =BC1
1449 typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1450 typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1451 typename LinearCombinationFilterType::Pointer combine0=LinearCombinationFilterType::New();
1452 combine0->SetInput(0,row0);
1453 combine0->SetInput(1,row1);
1455 combine0->SetB(-1.);
1457 typename CoefficientImageType::Pointer bc1Row=combine0->GetOutput();
1459 typename PasteImageFilterType::Pointer paster0=PasteImageFilterType::New();
1460 paster0->SetDestinationImage(m_PaddedCoefficientImage);
1461 paster0->SetSourceImage(bc1Row);
1462 destinationIndex[InputDimension-1]=0;
1463 paster0->SetDestinationIndex(destinationIndex);
1464 paster0->SetSourceRegion(bc1Row->GetLargestPossibleRegion());
1466 //---------------------------------
1467 // 1. Next two temporal rows are identical: paste 0-1 to 1-2
1468 typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1469 paster1->SetDestinationImage(paster0->GetOutput());
1470 paster1->SetSourceImage(m_CoefficientImage);
1471 sourceRegion.SetIndex(InputDimension-1,0);
1472 destinationIndex[InputDimension-1]=1;
1473 sourceRegion.SetSize(NInputDimensions-1,2);
1474 paster1->SetDestinationIndex(destinationIndex);
1475 paster1->SetSourceRegion(sourceRegion);
1477 //---------------------------------
1478 // 2. Middle row at index 3=BC
1479 // Extract row at index 1, 3 (of coeff image)
1480 //typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1481 typename CoefficientImageType::Pointer row3=this->ExtractTemporalRow(m_CoefficientImage, 3);
1484 typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1485 combine1->SetInput(0,row1);
1486 combine1->SetInput(1,row3);
1491 // Extract row at index 2 (coeff image)
1492 typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
1495 typename LinearCombinationFilterType::Pointer combine2=LinearCombinationFilterType::New();
1496 combine2->SetInput(0,combine1->GetOutput());
1497 combine2->SetInput(1,row2);
1498 combine2->SetA(-m_WeightRatio[3][1]);
1499 combine2->SetB(-1.);
1501 typename CoefficientImageType::Pointer bc2Row=combine2->GetOutput();
1503 // Paste middleRow at index 3 (padded coeff)
1504 typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1505 paster2->SetDestinationImage(paster1->GetOutput());
1506 paster2->SetSourceImage(bc2Row);
1507 destinationIndex[InputDimension-1]=3;
1508 paster2->SetDestinationIndex(destinationIndex);
1509 paster2->SetSourceRegion(bc2Row->GetLargestPossibleRegion());
1511 //---------------------------------
1512 // 3. Next two temporal rows are identical: paste 2,3 to 4,5
1513 typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1514 paster3->SetDestinationImage(paster2->GetOutput());
1515 paster3->SetSourceImage(m_CoefficientImage);
1516 sourceRegion.SetIndex(InputDimension-1,2);
1517 destinationIndex[InputDimension-1]=4;
1518 sourceRegion.SetSize(NInputDimensions-1,2);
1519 paster3->SetDestinationIndex(destinationIndex);
1520 paster3->SetSourceRegion(sourceRegion);
1522 // ---------------------------------
1523 // 4. Row at index 6=BC (paddedcoeff image)R1
1524 // Paste BC3 row at index 6
1525 typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
1526 paster4->SetDestinationImage(paster3->GetOutput());
1527 paster4->SetSourceImage(row0);
1528 destinationIndex[InputDimension-1]=6;
1529 paster4->SetDestinationIndex(destinationIndex);
1530 paster4->SetSourceRegion(row0->GetLargestPossibleRegion());
1532 // ---------------------------------
1533 // 5. Paste BC4 row at index 7
1534 typename LinearCombinationFilterType::Pointer combine3=LinearCombinationFilterType::New();
1535 combine3->SetInput(0,row0);
1536 combine3->SetInput(1,row3);
1538 combine3->SetB(-1.);
1540 typename CoefficientImageType::Pointer bc4Row=combine3->GetOutput();
1541 typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
1542 paster5->SetDestinationImage(paster4->GetOutput());
1543 paster5->SetSourceImage(bc4Row);
1544 destinationIndex[InputDimension-1]=7;
1545 paster5->SetDestinationIndex(destinationIndex);
1546 paster5->SetSourceRegion(bc4Row->GetLargestPossibleRegion());
1548 // Update the chain!
1550 m_PaddedCoefficientImage= paster5->GetOutput();
1557 // ----------------------------------------------------------------------
1558 // The sputnik with 4 internal CP
1559 // Periodic, constrained to zero at the reference
1560 // at position 3 and one indepent extremes copied
1561 // Coeff row BC1 0 1 BC2 2 BC2 3
1562 // PaddedCoeff R: 0 1 2 3 4 5 6
1564 // BC2= -weights[2]/weights[0] ( R2+R4 )
1565 // BC3= weights[2]/weights[0] ( R2-R4) + R1
1566 // ---------------------------------------------------------------------
1568 //---------------------------------
1569 // 1. First Row is equal to last row: paste 3 row to 0
1570 typename CoefficientImageType::Pointer row3=this->ExtractTemporalRow(m_CoefficientImage, 3);
1571 typename PasteImageFilterType::Pointer paster0=PasteImageFilterType::New();
1572 paster0->SetDestinationImage(m_PaddedCoefficientImage);
1573 paster0->SetSourceImage(row3);
1574 destinationIndex[InputDimension-1]=0;
1575 paster0->SetDestinationIndex(destinationIndex);
1576 paster0->SetSourceRegion(row3->GetLargestPossibleRegion());
1578 //---------------------------------
1579 // 1. Next two temporal rows are identical: paste 0-1 to 1-2
1580 typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1581 paster1->SetDestinationImage(paster0->GetOutput());
1582 paster1->SetSourceImage(m_CoefficientImage);
1583 sourceRegion.SetIndex(NInputDimensions-1,0);
1584 destinationIndex[InputDimension-1]=1;
1585 sourceRegion.SetSize(NInputDimensions-1,2);
1586 paster1->SetDestinationIndex(destinationIndex);
1587 paster1->SetSourceRegion(sourceRegion);
1589 //---------------------------------
1590 // 2. Middle row at index 3=BC
1591 // Extract row at index 1, 2 (of coeff image)
1592 typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1593 typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
1596 typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1597 combine1->SetInput(0,row1);
1598 combine1->SetInput(1,row2);
1599 combine1->SetA(-m_WeightRatio[2][0]);
1600 combine1->SetB(-m_WeightRatio[2][0]);
1602 typename CoefficientImageType::Pointer bc1Row=combine1->GetOutput();
1604 // Paste middleRow at index 3 (padded coeff)
1605 typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1606 paster2->SetDestinationImage(paster1->GetOutput());
1607 paster2->SetSourceImage(bc1Row);
1608 destinationIndex[InputDimension-1]=3;
1609 paster2->SetDestinationIndex(destinationIndex);
1610 paster2->SetSourceRegion(bc1Row->GetLargestPossibleRegion());
1612 //---------------------------------
1613 // 3. Next temporal row is identical: paste 2 to 4
1614 typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1615 paster3->SetDestinationImage(paster2->GetOutput());
1616 paster3->SetSourceImage(row2);
1617 destinationIndex[InputDimension-1]=4;
1618 paster3->SetDestinationIndex(destinationIndex);
1619 paster3->SetSourceRegion(row2->GetLargestPossibleRegion());
1621 //---------------------------------
1622 // 4. Final row at index 5=BC3 (paddedcoeff image)R1 R2 R4 corresponds to index in coeff image 0 1 2
1623 // Extract row at index 1, 2 (of coeff image)already done
1624 typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1627 typename LinearCombinationFilterType::Pointer combine3=LinearCombinationFilterType::New();
1628 combine3->SetInput(0,row1);
1629 combine3->SetInput(1,row2);
1631 combine3->SetB(-1.);
1634 typename LinearCombinationFilterType::Pointer combine4=LinearCombinationFilterType::New();
1635 combine4->SetInput(0,combine3->GetOutput());
1636 combine4->SetInput(1,row0);
1637 combine4->SetA(m_WeightRatio[2][0]);
1640 typename CoefficientImageType::Pointer bc2Row=combine4->GetOutput();
1642 // Paste BC row at index 5
1643 typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
1644 paster4->SetDestinationImage(paster3->GetOutput());
1645 paster4->SetSourceImage(bc2Row);
1646 destinationIndex[InputDimension-1]=5;
1647 paster4->SetDestinationIndex(destinationIndex);
1648 paster4->SetSourceRegion(bc2Row->GetLargestPossibleRegion());
1650 // Paste row 3 at index 6
1651 typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
1652 paster5->SetDestinationImage(paster4->GetOutput());
1653 paster5->SetSourceImage(row3);
1654 destinationIndex[InputDimension-1]=6;
1655 paster5->SetDestinationIndex(destinationIndex);
1656 paster5->SetSourceRegion(row3->GetLargestPossibleRegion());
1658 // Update the chain!
1660 m_PaddedCoefficientImage= paster5->GetOutput();
1668 // ----------------------------------------------------------------------
1669 // The sputnik with 5 internal CP
1670 // Periodic, constrained to zero at the reference
1671 // at position 3 and one indepent extreme
1672 // Coeff row BC1 0 1 BC2 2 3 BC3 4
1673 // PaddedCoeff R: 0 1 2 3 4 5 6 7
1674 // BC1= R2 + R5 - R7
1675 // BC2= -weights[3]/weights[1] ( R2+R5 ) - R4
1676 // BC3= R1 + 0.5 R2 - 0.5 R7
1677 // ----------------------------------------------------------------------
1678 //---------------------------------
1680 typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1681 typename CoefficientImageType::Pointer row3=this->ExtractTemporalRow(m_CoefficientImage, 3);
1682 typename CoefficientImageType::Pointer row4=this->ExtractTemporalRow(m_CoefficientImage, 4);
1685 typename LinearCombinationFilterType::Pointer combine0=LinearCombinationFilterType::New();
1686 combine0->SetInput(0,row1);
1687 combine0->SetInput(1,row3);
1690 //combine0->Update();
1691 typename LinearCombinationFilterType::Pointer combine0bis=LinearCombinationFilterType::New();
1692 combine0bis->SetInput(0,combine0->GetOutput());
1693 combine0bis->SetInput(1,row4);
1694 combine0bis->SetA(1.);
1695 combine0bis->SetB(-1.);
1696 combine0bis->Update();
1697 typename CoefficientImageType::Pointer bc1Row=combine0bis->GetOutput();
1699 typename PasteImageFilterType::Pointer paster0=PasteImageFilterType::New();
1700 paster0->SetDestinationImage(m_PaddedCoefficientImage);
1701 paster0->SetSourceImage(bc1Row);
1702 destinationIndex[InputDimension-1]=0;
1703 paster0->SetDestinationIndex(destinationIndex);
1704 paster0->SetSourceRegion(bc1Row->GetLargestPossibleRegion());
1706 //---------------------------------
1707 // 1. Next two temporal rows are identical: paste 0-1 to 1-2
1708 typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1709 paster1->SetDestinationImage(paster0->GetOutput());
1710 paster1->SetSourceImage(m_CoefficientImage);
1711 sourceRegion.SetIndex(NInputDimensions-1,0);
1712 destinationIndex[InputDimension-1]=1;
1713 sourceRegion.SetSize(NInputDimensions-1,2);
1714 paster1->SetDestinationIndex(destinationIndex);
1715 paster1->SetSourceRegion(sourceRegion);
1717 //---------------------------------
1718 // 2. Middle row at index 3=BC
1719 // Extract row at index 1, 2 (of coeff image)
1720 //typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1721 typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
1724 typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1725 combine1->SetInput(0,row1);
1726 combine1->SetInput(1,row3);
1731 typename LinearCombinationFilterType::Pointer combine2=LinearCombinationFilterType::New();
1732 combine2->SetInput(0,combine1->GetOutput());
1733 combine2->SetInput(1,row2);
1734 combine2->SetA(-m_WeightRatio[3][1]);
1735 combine2->SetB(-1.);
1737 typename CoefficientImageType::Pointer bc2Row=combine2->GetOutput();
1739 // Paste middleRow at index 3 (padded coeff)
1740 typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1741 paster2->SetDestinationImage(paster1->GetOutput());
1742 paster2->SetSourceImage(bc2Row);
1743 destinationIndex[InputDimension-1]=3;
1744 paster2->SetDestinationIndex(destinationIndex);
1745 paster2->SetSourceRegion(bc2Row->GetLargestPossibleRegion());
1747 //---------------------------------
1748 // 3. Next two temporal rows are identical: paste 2,3 to 4,5
1749 typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1750 paster3->SetDestinationImage(paster2->GetOutput());
1751 paster3->SetSourceImage(m_CoefficientImage);
1752 sourceRegion.SetIndex(InputDimension-1,2);
1753 destinationIndex[InputDimension-1]=4;
1754 sourceRegion.SetSize(NInputDimensions-1,2);
1755 paster3->SetDestinationIndex(destinationIndex);
1756 paster3->SetSourceRegion(sourceRegion);
1758 //---------------------------------
1759 // 4. Final row at index 6=BC3
1760 typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1763 typename LinearCombinationFilterType::Pointer combine3=LinearCombinationFilterType::New();
1764 combine3->SetInput(0,row1);
1765 combine3->SetInput(1,row4);
1767 combine3->SetB(-1.);
1768 typename LinearCombinationFilterType::Pointer combine4=LinearCombinationFilterType::New();
1769 combine4->SetInput(0,row0);
1770 combine4->SetInput(1,combine3->GetOutput());
1774 typename CoefficientImageType::Pointer bc3Row=combine4->GetOutput();
1776 // Paste BC row at index 6
1777 typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
1778 paster4->SetDestinationImage(paster3->GetOutput());
1779 paster4->SetSourceImage(bc3Row);
1780 destinationIndex[InputDimension-1]=6;
1781 paster4->SetDestinationIndex(destinationIndex);
1782 paster4->SetSourceRegion(bc3Row->GetLargestPossibleRegion());
1784 // Paste row 4 at index 7
1785 typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
1786 paster5->SetDestinationImage(paster4->GetOutput());
1787 paster5->SetSourceImage(row4);
1788 destinationIndex[InputDimension-1]=7;
1789 paster5->SetDestinationIndex(destinationIndex);
1790 paster5->SetSourceRegion(row4->GetLargestPossibleRegion());
1792 // Update the chain!
1794 m_PaddedCoefficientImage= paster5->GetOutput();
1801 // ----------------------------------------------------------------------
1802 // The diamond with 4 internal CP:
1803 // Periodic, constrained to zero at the reference at position 3
1804 // Coeff row 0 1 2 BC1 3 BC2 4
1805 // PaddedCoeff R:0 1 2 3 4 5 6
1806 // BC1= -weights[2]/weights[0] ( R2+R4 )
1807 // BC2= weights[2]/weights[0] ( R0+R2-R4-R6 ) + R1
1808 // ---------------------------------------------------------------------
1810 //---------------------------------
1811 // 1. First Three temporal rows are identical: paste 0-2 to 0-2
1812 typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1813 paster1->SetDestinationImage(m_PaddedCoefficientImage);
1814 paster1->SetSourceImage(m_CoefficientImage);
1815 sourceIndex.Fill(0);
1816 destinationIndex.Fill(0);
1817 sourceSize[NInputDimensions-1]=3;
1818 sourceRegion.SetSize(sourceSize);
1819 sourceRegion.SetIndex(sourceIndex);
1820 paster1->SetDestinationIndex(destinationIndex);
1821 paster1->SetSourceRegion(sourceRegion);
1823 //---------------------------------
1824 // 2. Middle row at index 3=BC
1825 // Extract row at index 0, 4 (of coeff image)
1826 typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
1827 typename CoefficientImageType::Pointer row3=this->ExtractTemporalRow(m_CoefficientImage, 3);
1830 typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1831 combine1->SetInput(0,row2);
1832 combine1->SetInput(1,row3);
1833 combine1->SetA(-m_WeightRatio[2][0]);
1834 combine1->SetB(-m_WeightRatio[2][0]);
1836 typename CoefficientImageType::Pointer bc1Row=combine1->GetOutput();
1838 // Paste middleRow at index 3 (padded coeff)
1839 typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1840 paster2->SetDestinationImage(paster1->GetOutput());
1841 paster2->SetSourceImage(bc1Row);
1842 destinationIndex[InputDimension-1]=3;
1843 paster2->SetDestinationIndex(destinationIndex);
1844 paster2->SetSourceRegion(bc1Row->GetLargestPossibleRegion());
1846 //---------------------------------
1847 // 3. Next row identical: paste 3 to 4
1848 typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1849 paster3->SetDestinationImage(paster2->GetOutput());
1850 paster3->SetSourceImage(row3);
1851 destinationIndex[InputDimension-1]=4;
1852 paster3->SetDestinationIndex(destinationIndex);
1853 paster3->SetSourceRegion(row3->GetLargestPossibleRegion());
1855 //---------------------------------
1856 // 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
1857 // Extract row at index 0, 2 (of coeff image)already done
1858 typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1859 typename CoefficientImageType::Pointer row4=this->ExtractTemporalRow(m_CoefficientImage, 4);
1860 typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1863 typename LinearCombinationFilterType::Pointer combine3=LinearCombinationFilterType::New();
1864 combine3->SetInput(0,row0);
1865 combine3->SetInput(1,row2);
1870 typename LinearCombinationFilterType::Pointer combine4=LinearCombinationFilterType::New();
1871 combine4->SetInput(0,row3);
1872 combine4->SetInput(1,row4);
1877 typename LinearCombinationFilterType::Pointer combine5=LinearCombinationFilterType::New();
1878 combine5->SetInput(0,combine3->GetOutput());
1879 combine5->SetInput(1,combine4->GetOutput());
1881 combine5->SetB(-1.);
1884 typename LinearCombinationFilterType::Pointer combine6=LinearCombinationFilterType::New();
1885 combine6->SetInput(0,combine5->GetOutput());
1886 combine6->SetInput(1,row1);
1887 combine6->SetA(m_WeightRatio[2][0]);
1890 typename CoefficientImageType::Pointer bc2Row=combine6->GetOutput();
1892 // Paste BC row at index 5
1893 typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
1894 paster4->SetDestinationImage(paster3->GetOutput());
1895 paster4->SetSourceImage(bc2Row);
1896 destinationIndex[InputDimension-1]=5;
1897 paster4->SetDestinationIndex(destinationIndex);
1898 paster4->SetSourceRegion(bc2Row->GetLargestPossibleRegion());
1900 // Paste last row at index 6
1901 typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
1902 paster5->SetDestinationImage(paster4->GetOutput());
1903 paster5->SetSourceImage(row4);
1904 destinationIndex[InputDimension-1]=6;
1905 paster5->SetDestinationIndex(destinationIndex);
1906 paster5->SetSourceRegion(row4->GetLargestPossibleRegion());
1908 // Update the chain!
1910 m_PaddedCoefficientImage= paster5->GetOutput();
1916 // ----------------------------------------------------------------------
1917 // The diamond with 5 internal CP:
1918 // periodic, constrained to zero at the reference at position 3.5
1919 // Coeff row 0 1 2 BC1 3 4 BC2 5
1920 // PaddedCoeff R:0 1 2 3 4 5 6 7
1921 // BC1= -weights[3]/weights[1] ( R2+R5 ) - R4
1922 // BC2= weights[2]/weights[0] ( R0+R2-R5-R7 ) + R1
1923 // ---------------------------------------------------------------------
1925 //---------------------------------
1926 // 1. First Three temporal rows are identical: paste 0-2 to 0-2
1927 typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1928 paster1->SetDestinationImage(m_PaddedCoefficientImage);
1929 paster1->SetSourceImage(m_CoefficientImage);
1930 sourceIndex.Fill(0);
1931 destinationIndex.Fill(0);
1932 sourceSize[NInputDimensions-1]=3;
1933 sourceRegion.SetSize(sourceSize);
1934 sourceRegion.SetIndex(sourceIndex);
1935 paster1->SetDestinationIndex(destinationIndex);
1936 paster1->SetSourceRegion(sourceRegion);
1938 //---------------------------------
1939 // 2. Middle row at index 3=BC
1940 // Extract row at index 0, 4 (of coeff image)
1941 typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
1942 typename CoefficientImageType::Pointer row4=this->ExtractTemporalRow(m_CoefficientImage, 4);
1945 typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1946 combine1->SetInput(0,row2);
1947 combine1->SetInput(1,row4);
1952 // Extract row at index 3 (coeff image)
1953 typename CoefficientImageType::Pointer row3=this->ExtractTemporalRow(m_CoefficientImage, 3);
1956 typename LinearCombinationFilterType::Pointer combine2=LinearCombinationFilterType::New();
1957 combine2->SetInput(0,combine1->GetOutput());
1958 combine2->SetInput(1,row3);
1959 combine2->SetA(-m_WeightRatio[3][1] );
1960 combine2->SetB(-1.);
1962 typename CoefficientImageType::Pointer bc1Row=combine2->GetOutput();
1964 // Paste middleRow at index 3 (padded coeff)
1965 typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1966 paster2->SetDestinationImage(paster1->GetOutput());
1967 paster2->SetSourceImage(bc1Row);
1968 destinationIndex[InputDimension-1]=3;
1969 paster2->SetDestinationIndex(destinationIndex);
1970 paster2->SetSourceRegion(bc1Row->GetLargestPossibleRegion());
1972 //---------------------------------
1973 // 3. Next two temporal rows are identical: paste 3,4 to 4,5
1974 typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1975 paster3->SetDestinationImage(paster2->GetOutput());
1976 paster3->SetSourceImage(m_CoefficientImage);
1977 sourceIndex[InputDimension-1]=3;
1978 destinationIndex[InputDimension-1]=4;
1979 sourceSize[NInputDimensions-1]=2;
1980 sourceRegion.SetSize(sourceSize);
1981 sourceRegion.SetIndex(sourceIndex);
1982 paster3->SetDestinationIndex(destinationIndex);
1983 paster3->SetSourceRegion(sourceRegion);
1985 //---------------------------------
1986 // 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
1987 // Extract row at index 0, 2 (of coeff image)already done
1988 typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1989 typename CoefficientImageType::Pointer row5=this->ExtractTemporalRow(m_CoefficientImage, 5);
1990 typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1993 typename LinearCombinationFilterType::Pointer combine3=LinearCombinationFilterType::New();
1994 combine3->SetInput(0,row0);
1995 combine3->SetInput(1,row2);
2000 typename LinearCombinationFilterType::Pointer combine4=LinearCombinationFilterType::New();
2001 combine4->SetInput(0,row4);
2002 combine4->SetInput(1,row5);
2007 typename LinearCombinationFilterType::Pointer combine5=LinearCombinationFilterType::New();
2008 combine5->SetInput(0,combine3->GetOutput());
2009 combine5->SetInput(1,combine4->GetOutput());
2011 combine5->SetB(-1.);
2014 typename LinearCombinationFilterType::Pointer combine6=LinearCombinationFilterType::New();
2015 combine6->SetInput(0,combine5->GetOutput());
2016 combine6->SetInput(1,row1);
2017 combine6->SetA(m_WeightRatio[2][0]);
2020 typename CoefficientImageType::Pointer bc2Row=combine6->GetOutput();
2022 // Paste BC row at index 6
2023 typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
2024 paster4->SetDestinationImage(paster3->GetOutput());
2025 paster4->SetSourceImage(bc2Row);
2026 destinationIndex[InputDimension-1]=6;
2027 paster4->SetDestinationIndex(destinationIndex);
2028 paster4->SetSourceRegion(bc2Row->GetLargestPossibleRegion());
2030 // Paste last row at index 7
2031 typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
2032 paster5->SetDestinationImage(paster4->GetOutput());
2033 paster5->SetSourceImage(row5);
2034 destinationIndex[InputDimension-1]=7;
2035 paster5->SetDestinationIndex(destinationIndex);
2036 paster5->SetSourceRegion(row5->GetLargestPossibleRegion());
2038 // Update the chain!
2040 m_PaddedCoefficientImage= paster5->GetOutput();
2047 // ----------------------------------------------------------------------
2048 // The sputnik with 5 internal CP T''(0)=T''(10)
2049 // Periodic, constrained to zero at the reference
2050 // at position 3 and one indepent extreme
2051 // Coeff row BC1 0 1 BC2 2 3 BC3 4
2052 // PaddedCoeff R: 0 1 2 3 4 5 6 7
2054 // BC2= -weights[3]/weights[1] ( R2+R5 ) - R4
2056 // ---------------------------------------------------------------------
2058 //---------------------------------
2060 typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
2061 typename CoefficientImageType::Pointer row3=this->ExtractTemporalRow(m_CoefficientImage, 3);
2062 typename CoefficientImageType::Pointer row4=this->ExtractTemporalRow(m_CoefficientImage, 4);
2065 typename LinearCombinationFilterType::Pointer combine0=LinearCombinationFilterType::New();
2066 combine0->SetInput(0,row3);
2067 combine0->SetInput(1,row4);
2070 typename LinearCombinationFilterType::Pointer combine0bis=LinearCombinationFilterType::New();
2071 combine0bis->SetInput(0,combine0->GetOutput());
2072 combine0bis->SetInput(1,row1);
2073 combine0bis->SetA(1.);
2074 combine0bis->SetB(-1.);
2075 combine0bis->Update();
2076 typename CoefficientImageType::Pointer bc1Row=combine0bis->GetOutput();
2078 typename PasteImageFilterType::Pointer paster0=PasteImageFilterType::New();
2079 paster0->SetDestinationImage(m_PaddedCoefficientImage);
2080 paster0->SetSourceImage(bc1Row);
2081 destinationIndex[InputDimension-1]=0;
2082 paster0->SetDestinationIndex(destinationIndex);
2083 paster0->SetSourceRegion(bc1Row->GetLargestPossibleRegion());
2085 //---------------------------------
2086 // 1. Next two temporal rows are identical: paste 0-1 to 1-2
2087 typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
2088 paster1->SetDestinationImage(paster0->GetOutput());
2089 paster1->SetSourceImage(m_CoefficientImage);
2090 sourceRegion.SetIndex(NInputDimensions-1,0);
2091 destinationIndex[InputDimension-1]=1;
2092 sourceRegion.SetSize(NInputDimensions-1,2);
2093 paster1->SetDestinationIndex(destinationIndex);
2094 paster1->SetSourceRegion(sourceRegion);
2096 //---------------------------------
2097 // 2. Middle row at index 3=BC
2098 // Extract row at index 1, 2 (of coeff image)
2099 // typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
2100 typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
2103 typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
2104 combine1->SetInput(0,row1);
2105 combine1->SetInput(1,row3);
2110 typename LinearCombinationFilterType::Pointer combine2=LinearCombinationFilterType::New();
2111 combine2->SetInput(0,combine1->GetOutput());
2112 combine2->SetInput(1,row2);
2113 combine2->SetA(-m_WeightRatio[3][1]);
2114 combine2->SetB(-1.);
2116 typename CoefficientImageType::Pointer bc2Row=combine2->GetOutput();
2118 // Paste middleRow at index 3 (padded coeff)
2119 typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
2120 paster2->SetDestinationImage(paster1->GetOutput());
2121 paster2->SetSourceImage(bc2Row);
2122 destinationIndex[InputDimension-1]=3;
2123 paster2->SetDestinationIndex(destinationIndex);
2124 paster2->SetSourceRegion(bc2Row->GetLargestPossibleRegion());
2126 //---------------------------------
2127 // 3. Next two temporal rows are identical: paste 2,3 to 4,5
2128 typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
2129 paster3->SetDestinationImage(paster2->GetOutput());
2130 paster3->SetSourceImage(m_CoefficientImage);
2131 sourceRegion.SetIndex(InputDimension-1,2);
2132 destinationIndex[InputDimension-1]=4;
2133 sourceRegion.SetSize(NInputDimensions-1,2);
2134 paster3->SetDestinationIndex(destinationIndex);
2135 paster3->SetSourceRegion(sourceRegion);
2137 //---------------------------------
2138 // 4. Final row at index 6=BC3
2139 typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
2142 // typename LinearCombinationFilterType::Pointer combine3=LinearCombinationFilterType::New();
2143 // combine3->SetInput(0,row0);
2144 // combine3->SetInput(1,row1);
2145 // combine3->SetA(1.);
2146 // combine3->SetB(0.5);
2147 // combine3->Update();
2148 // typename CoefficientImageType::Pointer bc3Row=combine3->GetOutput();
2150 // Paste BC row at index 6
2151 typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
2152 paster4->SetDestinationImage(paster3->GetOutput());
2153 paster4->SetSourceImage(row0);
2154 destinationIndex[InputDimension-1]=6;
2155 paster4->SetDestinationIndex(destinationIndex);
2156 paster4->SetSourceRegion(row0->GetLargestPossibleRegion());
2158 // Paste row 4 at index 7
2159 typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
2160 paster5->SetDestinationImage(paster4->GetOutput());
2161 paster5->SetSourceImage(row4);
2162 destinationIndex[InputDimension-1]=7;
2163 paster5->SetDestinationIndex(destinationIndex);
2164 paster5->SetSourceRegion(row4->GetLargestPossibleRegion());
2166 // Update the chain!
2168 m_PaddedCoefficientImage= paster5->GetOutput();
2174 DD ("Shape not available");
2180 // // Extract coefficients from padded version
2181 // template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
2183 // ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2184 // ::ExtractCoefficientImage( )
2186 // ////DD("extract coeff image");
2187 // typename ExtractImageFilterType::Pointer extractImageFilter = ExtractImageFilterType::New();
2188 // typename CoefficientImageType::RegionType extractionRegion=m_PaddedCoefficientImage->GetLargestPossibleRegion();
2189 // typename CoefficientImageType::RegionType::SizeType extractionSize=extractionRegion.GetSize();
2190 // typename CoefficientImageType::RegionType::IndexType extractionIndex = extractionRegion.GetIndex();
2191 // extractionSize[InputDimension-1]-=4;
2192 // extractionIndex[InputDimension-1]=2;
2193 // extractionRegion.SetSize(extractionSize);
2194 // extractionRegion.SetIndex(extractionIndex);
2195 // extractImageFilter->SetInput(m_PaddedCoefficientImage);
2196 // extractImageFilter->SetExtractionRegion(extractionRegion);
2197 // extractImageFilter->Update();
2198 // m_CoefficientImage=extractImageFilter->GetOutput();
2203 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
2205 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2206 ::PrintSelf(std::ostream &os, itk::Indent indent) const
2209 this->Superclass::PrintSelf(os, indent);
2211 os << indent << "GridRegion: " << m_GridRegion << std::endl;
2212 os << indent << "GridOrigin: " << m_GridOrigin << std::endl;
2213 os << indent << "GridSpacing: " << m_GridSpacing << std::endl;
2214 os << indent << "GridDirection: " << m_GridDirection << std::endl;
2215 os << indent << "IndexToPoint: " << m_IndexToPoint << std::endl;
2216 os << indent << "PointToIndex: " << m_PointToIndex << std::endl;
2218 os << indent << "CoefficientImage: [ ";
2219 os << m_CoefficientImage.GetPointer() << " ]" << std::endl;
2221 os << indent << "WrappedImage: [ ";
2222 os << m_WrappedImage.GetPointer() << " ]" << std::endl;
2224 os << indent << "InputParametersPointer: "
2225 << m_InputParametersPointer << std::endl;
2226 os << indent << "ValidRegion: " << m_ValidRegion << std::endl;
2227 os << indent << "LastJacobianIndex: " << m_LastJacobianIndex << std::endl;
2228 os << indent << "BulkTransform: ";
2229 os << m_BulkTransform << std::endl;
2231 if ( m_BulkTransform )
2233 os << indent << "BulkTransformType: "
2234 << m_BulkTransform->GetNameOfClass() << std::endl;
2236 os << indent << "VectorBSplineInterpolator: ";
2237 os << m_VectorInterpolator.GetPointer() << std::endl;
2238 os << indent << "Mask: ";
2239 os << m_Mask<< std::endl;
2244 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
2246 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2247 ::InsideValidRegion( const ContinuousIndexType& index ) const
2251 if ( !m_ValidRegion.IsInside( index ) )
2256 // JV verify all dimensions
2259 typedef typename ContinuousIndexType::ValueType ValueType;
2260 for( unsigned int j = 0; j < NInputDimensions; j++ )
2262 if (m_SplineOrderOdd[j])
2264 if ( index[j] >= static_cast<ValueType>( m_ValidRegionLast[j] ) )
2276 // Transform a point
2277 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
2278 typename ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2280 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2281 ::TransformPoint(const InputPointType &inputPoint) const
2284 InputPointType transformedPoint;
2285 OutputPointType outputPoint;
2288 if ( m_BulkTransform )
2290 transformedPoint = m_BulkTransform->TransformPoint( inputPoint );
2294 transformedPoint = inputPoint;
2297 // Deformable transform
2298 if ( m_PaddedCoefficientImage )
2300 // Check if inside mask
2301 if(m_Mask && !(m_Mask->IsInside(inputPoint) ) )
2303 // Outside: no (deformable) displacement
2304 return transformedPoint;
2307 // Check if inside valid region
2309 ContinuousIndexType index;
2310 this->TransformPointToContinuousIndex( inputPoint, index );
2311 inside = this->InsideValidRegion( index );
2314 // Outside: no (deformable) displacement
2315 DD("outside valid region");
2316 return transformedPoint;
2319 // Call the vector interpolator
2320 itk::Vector<TCoordRep,SpaceDimension> displacement=m_VectorInterpolator->EvaluateAtContinuousIndex(index);
2322 // JV add for the spatial dimensions
2323 outputPoint=transformedPoint;
2324 for (unsigned int i=0; i<NInputDimensions-1; i++)
2325 outputPoint[i] += displacement[i];
2331 itkWarningMacro( << "B-spline coefficients have not been set" );
2332 outputPoint = transformedPoint;
2340 //JV Deformably transform a point
2341 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
2342 typename ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2344 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2345 ::DeformablyTransformPoint(const InputPointType &inputPoint) const
2347 OutputPointType outputPoint;
2348 if ( m_PaddedCoefficientImage )
2351 // Check if inside mask
2352 if(m_Mask && !(m_Mask->IsInside(inputPoint) ) )
2354 // Outside: no (deformable) displacement
2360 ContinuousIndexType index;
2361 this->TransformPointToContinuousIndex( inputPoint, index );
2362 inside = this->InsideValidRegion( index );
2366 //outside: no (deformable) displacement
2367 outputPoint = inputPoint;
2371 // Call the vector interpolator
2372 itk::Vector<TCoordRep,SpaceDimension> displacement=m_VectorInterpolator->EvaluateAtContinuousIndex(index);
2374 // JV add for the spatial dimensions
2375 outputPoint=inputPoint;
2376 for (unsigned int i=0; i<NInputDimensions-1; i++)
2377 outputPoint[i] += displacement[i];
2380 // No coefficients available
2383 itkWarningMacro( << "B-spline coefficients have not been set" );
2384 outputPoint = inputPoint;
2391 // JV weights are identical as for transformpoint, could be done simultaneously in metric!!!!
2392 // Compute the Jacobian in one position
2393 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
2394 #if ITK_VERSION_MAJOR >= 4
2396 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2397 ::ComputeJacobianWithRespectToParameters( const InputPointType & point, JacobianType & jacobian) const
2400 typename ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2402 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2403 ::GetJacobian( const InputPointType & point ) const
2407 //========================================================
2408 // Zero all components of jacobian
2409 //========================================================
2410 // JV not thread safe (m_LastJacobianIndex), instantiate N transforms
2411 // NOTE: for efficiency, we only need to zero out the coefficients
2412 // that got fill last time this method was called.
2413 unsigned int j=0,b=0;
2415 //Set the previously-set to zero
2416 for ( j = 0; j < SpaceDimension; j++ )
2418 m_FirstIterator[j].GoToBegin();
2419 while ( !m_FirstIterator[j].IsAtEnd() )
2421 m_FirstIterator[j].Set( m_ZeroVector );
2422 ++(m_FirstIterator[j]);
2426 //Set the previously-set to zero
2427 for ( j = 0; j < SpaceDimension; j++ )
2429 m_SecondIterator[j].GoToBegin();
2430 while ( ! (m_SecondIterator[j]).IsAtEnd() )
2432 m_SecondIterator[j].Set( m_ZeroVector );
2433 ++(m_SecondIterator[j]);
2437 //Set the previously-set to zero
2439 for ( j = 0; j < SpaceDimension; j++ )
2441 m_ThirdIterator[j].GoToBegin();
2442 while ( ! (m_ThirdIterator[j]).IsAtEnd() )
2444 m_ThirdIterator[j].Set( m_ZeroVector );
2445 ++(m_ThirdIterator[j]);
2449 //Set the previously-set to zero
2451 for (b=0; b<m_BCSize;b++)
2452 if ( !( m_FirstRegion.IsInside(m_BCRegions[b]) | m_SecondRegion.IsInside(m_BCRegions[b]) ) )
2453 for ( j = 0; j < SpaceDimension; j++ )
2455 m_BCIterators[j][b].GoToBegin();
2456 while ( ! (m_BCIterators[j][b]).IsAtEnd() )
2458 m_BCIterators[j][b].Set( m_ZeroVector );
2459 ++(m_BCIterators[j][b]);
2463 //Set the previously-set to zero
2465 for (b=0; b<m_BC2Size;b++)
2466 if ( !( m_FirstRegion.IsInside(m_BC2Regions[b]) | m_SecondRegion.IsInside(m_BC2Regions[b]) ) )
2467 for ( j = 0; j < SpaceDimension; j++ )
2469 m_BC2Iterators[j][b].GoToBegin();
2470 while ( ! (m_BC2Iterators[j][b]).IsAtEnd() )
2472 m_BC2Iterators[j][b].Set( m_ZeroVector );
2473 ++(m_BC2Iterators[j][b]);
2477 //Set the previously-set to zero
2479 for (b=0; b<m_BC3Size;b++)
2480 if ( !( m_FirstRegion.IsInside(m_BC3Regions[b]) | m_SecondRegion.IsInside(m_BC3Regions[b]) ) )
2481 for ( j = 0; j < SpaceDimension; j++ )
2483 m_BC3Iterators[j][b].GoToBegin();
2484 while ( ! (m_BC3Iterators[j][b]).IsAtEnd() )
2486 m_BC3Iterators[j][b].Set( m_ZeroVector );
2487 ++(m_BC3Iterators[j][b]);
2492 //========================================================
2493 // For each dimension, copy the weight to the support region
2494 //========================================================
2496 // Check if inside mask
2497 if(m_Mask && !(m_Mask->IsInside(point) ) )
2499 // Outside: no (deformable) displacement
2500 #if ITK_VERSION_MAJOR >= 4
2501 jacobian = m_SharedDataBSplineJacobian;
2504 return this->m_Jacobian;
2509 this->TransformPointToContinuousIndex( point, m_Index );
2511 // NOTE: if the support region does not lie totally within the grid
2512 // we assume zero displacement and return the input point
2513 if ( !this->InsideValidRegion( m_Index ) )
2515 #if ITK_VERSION_MAJOR >= 4
2516 jacobian = m_SharedDataBSplineJacobian;
2519 return this->m_Jacobian;
2523 // Compute interpolation weights
2524 const WeightsDataType *weights=NULL;
2525 m_VectorInterpolator->EvaluateWeightsAtContinuousIndex( m_Index, &weights, m_LastJacobianIndex);
2528 m_SupportRegion.SetIndex( m_LastJacobianIndex );
2529 WrapRegion(m_SupportRegion, m_FirstRegion, m_SecondRegion, m_ThirdRegion, m_BCRegions, m_BCValues, m_BC2Regions, m_BC2Values, m_BC3Regions, m_BC3Values, m_InitialOffset);
2530 m_ThirdSize=m_ThirdRegion.GetSize()[InputDimension-1];
2531 m_BCSize=m_BCRegions.size();
2532 m_BC2Size=m_BC2Regions.size();
2533 m_BC3Size=m_BC3Regions.size();
2535 // Reset the iterators
2536 for ( j = 0; j < SpaceDimension ; j++ )
2538 m_FirstIterator[j] = IteratorType( m_JacobianImage[j], m_FirstRegion);
2539 m_SecondIterator[j] = IteratorType( m_JacobianImage[j], m_SecondRegion);
2540 if(m_ThirdSize) m_ThirdIterator[j] = IteratorType( m_JacobianImage[j], m_ThirdRegion);
2542 m_BCIterators[j].resize(m_BCSize);
2543 for (b=0; b<m_BCSize;b++)
2544 m_BCIterators[j][b]= IteratorType( m_JacobianImage[j], m_BCRegions[b]);
2545 m_BC2Iterators[j].resize(m_BC2Size);
2546 for (b=0; b<m_BC2Size;b++)
2547 m_BC2Iterators[j][b]= IteratorType( m_JacobianImage[j], m_BC2Regions[b]);
2548 m_BC3Iterators[j].resize(m_BC3Size);
2549 for (b=0; b<m_BC3Size;b++)
2550 m_BC3Iterators[j][b]= IteratorType( m_JacobianImage[j], m_BC3Regions[b]);
2553 // Skip if on a fixed condition
2556 if (m_BCSize) weights+=m_InitialOffset*m_BCRegions[0].GetNumberOfPixels();
2557 else std::cerr<<"InitialOffset without BCSize: Error!!!!!!!"<<std::endl;
2560 //copy weight to jacobian image
2561 for ( j = 0; j < SpaceDimension; j++ )
2563 // For each dimension, copy the weight to the support region
2564 while ( ! (m_FirstIterator[j]).IsAtEnd() )
2566 m_ZeroVector[j]=*weights;
2567 (m_FirstIterator[j]).Set( m_ZeroVector);
2568 ++(m_FirstIterator[j]);
2572 m_ZeroVector[j]=itk::NumericTraits<JacobianValueType>::Zero;
2573 if (j != SpaceDimension-1) weights-=m_FirstRegion.GetNumberOfPixels();
2576 // Skip BC1 and go to the second region
2577 if (m_BCSize) weights+=m_BCRegions[0].GetNumberOfPixels();
2579 // For each dimension, copy the weight to the support region
2580 //copy weight to jacobian image
2581 for ( j = 0; j < SpaceDimension; j++ )
2583 while ( ! (m_SecondIterator[j]).IsAtEnd() )
2585 m_ZeroVector[j]=*weights;
2586 (m_SecondIterator[j]).Set( m_ZeroVector);
2587 ++(m_SecondIterator[j]);
2590 weights-=m_SecondRegion.GetNumberOfPixels();
2591 m_ZeroVector[j]=itk::NumericTraits<JacobianValueType>::Zero;
2597 // Put pointer in correct position
2598 weights-=m_BCRegions[0].GetNumberOfPixels();
2600 for ( j = 0; j < SpaceDimension; j++ )
2602 for ( b=0; b < m_BCSize; b++ )
2604 while ( ! (m_BCIterators[j][b]).IsAtEnd() )
2606 //copy weight to jacobian image
2607 m_ZeroVector[j]=(*weights) * m_BCValues[b];
2608 (m_BCIterators[j][b]).Value()+= m_ZeroVector;
2609 ++(m_BCIterators[j][b]);
2612 weights-=m_BCRegions[b].GetNumberOfPixels();
2614 m_ZeroVector[j]=itk::NumericTraits<JacobianValueType>::Zero;
2616 weights+=m_BCRegions[m_BCSize-1].GetNumberOfPixels();
2619 // Add the BC2 to the weights
2622 // Move further in the weights pointer
2623 weights+=m_SecondRegion.GetNumberOfPixels();
2625 for ( j = 0; j < SpaceDimension; j++ )
2627 for ( b=0; b < m_BC2Size; b++ )
2629 while ( ! (m_BC2Iterators[j][b]).IsAtEnd() )
2631 //copy weight to jacobian image
2632 m_ZeroVector[j]=(*weights) * m_BC2Values[b];
2633 (m_BC2Iterators[j][b]).Value()+= m_ZeroVector;
2634 ++(m_BC2Iterators[j][b]);
2637 weights-=m_BC2Regions[b].GetNumberOfPixels();
2639 m_ZeroVector[j]=itk::NumericTraits<JacobianValueType>::Zero;
2641 // Move further in the weights pointer
2642 weights+=m_BC2Regions[m_BC2Size-1].GetNumberOfPixels();
2648 // For each dimension, copy the weight to the support region
2649 //copy weight to jacobian image
2650 for ( j = 0; j < SpaceDimension; j++ )
2652 while ( ! (m_ThirdIterator[j]).IsAtEnd() )
2654 m_ZeroVector[j]=*weights;
2655 (m_ThirdIterator[j]).Value()+= m_ZeroVector;
2656 ++(m_ThirdIterator[j]);
2660 // Move further in the weights pointer?
2661 if (j != SpaceDimension-1) weights-=m_ThirdRegion.GetNumberOfPixels();
2662 m_ZeroVector[j]=itk::NumericTraits<JacobianValueType>::Zero;
2666 // Add the BC3 to the weights
2669 for ( j = 0; j < SpaceDimension; j++ )
2671 for ( b=0; b < m_BC3Size; b++ )
2673 while ( ! (m_BC3Iterators[j][b]).IsAtEnd() )
2675 //copy weight to jacobian image
2676 m_ZeroVector[j]=(*weights) * m_BC3Values[b];
2677 (m_BC3Iterators[j][b]).Value()+= m_ZeroVector;
2678 ++(m_BC3Iterators[j][b]);
2681 weights-=m_BC3Regions[b].GetNumberOfPixels();
2683 m_ZeroVector[j]=itk::NumericTraits<JacobianValueType>::Zero;
2687 // Return the result
2688 #if ITK_VERSION_MAJOR >= 4
2689 jacobian = m_SharedDataBSplineJacobian;
2691 return this->m_Jacobian;
2696 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
2698 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2699 ::WrapRegion( const RegionType & m_SupportRegion,
2700 RegionType & m_FirstRegion,
2701 RegionType & m_SecondRegion,
2702 RegionType & m_ThirdRegion,
2703 std::vector<RegionType>& m_BCRegions,std::vector<double>& m_BCValues,
2704 std::vector<RegionType>& m_BC2Regions,std::vector<double>& m_BC2Values,
2705 std::vector<RegionType>& m_BC3Regions,std::vector<double>& m_BC3Values,
2706 unsigned int& m_InitialOffset ) const
2709 // Synchronize regions
2711 m_FirstRegion=m_SupportRegion;
2712 m_BCRegion=m_SupportRegion;
2713 m_BCRegion.SetSize(InputDimension-1,1);
2714 m_SecondRegion=m_SupportRegion;
2715 m_ThirdRegion=m_SupportRegion;
2716 m_ThirdRegion.SetSize(InputDimension-1,0);
2717 m_BC3Regions.resize(0);
2720 // BC depends on shape
2721 switch(m_TransformShape)
2726 2: rabbit 4 CP 3 DOF
2727 3: rabbit 5 CP 4 DOF
2728 4: sputnik 4 CP 4 DOF
2729 5: sputnik 5 CP 5 DOF
2730 6: diamond 6 CP 5 DOF
2731 7: diamond 7 CP 6 DOF
2736 // ----------------------------------------------------------------------
2737 // The egg with 4 internal CP (starting from inhale)
2738 // Periodic, constrained to zero at the reference
2739 // at position 3 and
2740 // Coeff row BC1 BC2 0 BC3 1 2 BC4
2741 // PaddedCoeff R: 0 1 2 3 4 5 6
2744 // BC3= -weights[2]/weights[0] ( R2+R4 )
2746 // ---------------------------------------------------------------------
2747 switch(m_SupportRegion.GetIndex(InputDimension-1))
2752 m_FirstRegion.SetSize(InputDimension-1,2);
2753 m_FirstRegion.SetIndex(InputDimension-1,1);
2756 m_BCRegions.resize(0);
2759 m_SecondRegion.SetSize(InputDimension-1,1);
2760 m_SecondRegion.SetIndex(InputDimension-1,0);
2763 m_BC2Regions.resize(2);
2764 m_BC2Values.resize(2);
2765 m_BCRegion.SetIndex(InputDimension-1,0);
2766 m_BC2Regions[0]=m_BCRegion;
2767 m_BC2Values[0]=-m_WeightRatio[2][0];
2768 m_BCRegion.SetIndex(InputDimension-1,1);
2769 m_BC2Regions[1]=m_BCRegion;
2770 m_BC2Values[1]= -m_WeightRatio[2][0];
2777 m_FirstRegion.SetSize(InputDimension-1,1);
2778 m_FirstRegion.SetIndex(InputDimension-1,2);
2781 m_BCRegions.resize(0);
2784 m_SecondRegion.SetSize(InputDimension-1,1);
2785 m_SecondRegion.SetIndex(InputDimension-1,0);
2788 m_BC2Regions.resize(2);
2789 m_BC2Values.resize(2);
2790 m_BCRegion.SetIndex(InputDimension-1,0);
2791 m_BC2Regions[0]=m_BCRegion;
2792 m_BC2Values[0]=-m_WeightRatio[2][0];
2793 m_BCRegion.SetIndex(InputDimension-1,1);
2794 m_BC2Regions[1]=m_BCRegion;
2795 m_BC2Values[1]= -m_WeightRatio[2][0];
2798 m_FirstRegion.SetSize(InputDimension-1,1);
2799 m_FirstRegion.SetIndex(InputDimension-1,1);
2806 m_FirstRegion.SetSize(InputDimension-1,1);
2807 m_FirstRegion.SetIndex(InputDimension-1,0);
2810 m_BCRegions.resize(2);
2811 m_BCValues.resize(2);
2812 m_BCRegion.SetIndex(InputDimension-1,0);
2813 m_BCRegions[0]=m_BCRegion;
2814 m_BCValues[0]=-m_WeightRatio[2][0];
2815 m_BCRegion.SetIndex(InputDimension-1,1);
2816 m_BCRegions[1]=m_BCRegion;
2817 m_BCValues[1]= -m_WeightRatio[2][0];
2820 m_SecondRegion.SetSize(InputDimension-1,2);
2821 m_SecondRegion.SetIndex(InputDimension-1,1);
2824 m_BC2Regions.resize(0);
2831 m_FirstRegion.SetSize(InputDimension-1,0);
2834 m_BCRegions.resize(2);
2835 m_BCValues.resize(2);
2836 m_BCRegion.SetIndex(InputDimension-1,0);
2837 m_BCRegions[0]=m_BCRegion;
2838 m_BCValues[0]=-m_WeightRatio[2][0];
2839 m_BCRegion.SetIndex(InputDimension-1,1);
2840 m_BCRegions[1]=m_BCRegion;
2841 m_BCValues[1]= -m_WeightRatio[2][0];
2844 m_SecondRegion.SetSize(InputDimension-1,2);
2845 m_SecondRegion.SetIndex(InputDimension-1,1);
2848 m_BC2Regions.resize(1);
2849 m_BC2Values.resize(1);
2850 m_BCRegion.SetIndex(InputDimension-1,0);
2851 m_BC2Regions[0]=m_BCRegion;
2859 DD("supportindex > 3 ???");
2860 DD(m_SupportRegion.GetIndex(InputDimension-1));
2862 } // end switch index
2869 // ----------------------------------------------------------------------
2870 // The egg with 5 internal CP (starting from inhale)
2871 // Periodic, constrained to zero at the reference
2872 // at position 3 and
2873 // Coeff row BC1 BC2 0 BC3 1 2 3 BC4
2874 // PaddedCoeff R: 0 1 2 3 4 5 6 7
2877 // BC3= -weights[3]/weights[1] ( R2+R5 ) - R4
2879 // ---------------------------------------------------------------------
2880 switch(m_SupportRegion.GetIndex(InputDimension-1))
2885 m_FirstRegion.SetSize(InputDimension-1,2);
2886 m_FirstRegion.SetIndex(InputDimension-1,2);
2889 m_BCRegions.resize(0);
2892 m_SecondRegion.SetSize(InputDimension-1,1);
2893 m_SecondRegion.SetIndex(InputDimension-1,0);
2896 m_BC2Regions.resize(3);
2897 m_BC2Values.resize(3);
2898 m_BCRegion.SetIndex(InputDimension-1,0);
2899 m_BC2Regions[0]=m_BCRegion;
2900 m_BC2Values[0]=-m_WeightRatio[3][1];
2901 m_BCRegion.SetIndex(InputDimension-1,1);
2902 m_BC2Regions[1]=m_BCRegion;
2904 m_BCRegion.SetIndex(InputDimension-1,2);
2905 m_BC2Regions[2]=m_BCRegion;
2906 m_BC2Values[2]=-m_WeightRatio[3][1];
2913 m_FirstRegion.SetSize(InputDimension-1,1);
2914 m_FirstRegion.SetIndex(InputDimension-1,3);
2917 m_BCRegions.resize(0);
2920 m_SecondRegion.SetSize(InputDimension-1,1);
2921 m_SecondRegion.SetIndex(InputDimension-1,0);
2924 m_BC2Regions.resize(3);
2925 m_BC2Values.resize(3);
2926 m_BCRegion.SetIndex(InputDimension-1,0);
2927 m_BC2Regions[0]=m_BCRegion;
2928 m_BC2Values[0]=-m_WeightRatio[3][1];
2929 m_BCRegion.SetIndex(InputDimension-1,1);
2930 m_BC2Regions[1]=m_BCRegion;
2932 m_BCRegion.SetIndex(InputDimension-1,2);
2933 m_BC2Regions[2]=m_BCRegion;
2934 m_BC2Values[2]=-m_WeightRatio[3][1];
2937 m_FirstRegion.SetSize(InputDimension-1,1);
2938 m_FirstRegion.SetIndex(InputDimension-1,1);
2945 m_FirstRegion.SetSize(InputDimension-1,1);
2946 m_FirstRegion.SetIndex(InputDimension-1,0);
2949 m_BCRegions.resize(3);
2950 m_BCValues.resize(3);
2951 m_BCRegion.SetIndex(InputDimension-1,0);
2952 m_BCRegions[0]=m_BCRegion;
2953 m_BCValues[0]=-m_WeightRatio[3][1];
2954 m_BCRegion.SetIndex(InputDimension-1,1);
2955 m_BCRegions[1]=m_BCRegion;
2957 m_BCRegion.SetIndex(InputDimension-1,2);
2958 m_BCRegions[2]=m_BCRegion;
2959 m_BCValues[2]=-m_WeightRatio[3][1];
2962 m_SecondRegion.SetSize(InputDimension-1,2);
2963 m_SecondRegion.SetIndex(InputDimension-1,1);
2966 m_BC2Regions.resize(0);
2973 m_FirstRegion.SetSize(InputDimension-1,0);
2976 m_BCRegions.resize(3);
2977 m_BCValues.resize(3);
2978 m_BCRegion.SetIndex(InputDimension-1,0);
2979 m_BCRegions[0]=m_BCRegion;
2980 m_BCValues[0]=-m_WeightRatio[3][1];
2981 m_BCRegion.SetIndex(InputDimension-1,1);
2982 m_BCRegions[1]=m_BCRegion;
2984 m_BCRegion.SetIndex(InputDimension-1,2);
2985 m_BCRegions[2]=m_BCRegion;
2986 m_BCValues[2]=-m_WeightRatio[3][1];
2989 m_SecondRegion.SetSize(InputDimension-1,3);
2990 m_SecondRegion.SetIndex(InputDimension-1,1);
2993 m_BC2Regions.resize(0);
3000 m_FirstRegion.SetSize(InputDimension-1,3);
3001 m_FirstRegion.SetIndex(InputDimension-1,1);
3004 m_BCRegions.resize(0);
3007 m_SecondRegion.SetSize(InputDimension-1,1);
3008 m_SecondRegion.SetIndex(InputDimension-1,0);
3011 m_BC2Regions.resize(0);
3018 DD("supportindex > 3 ???");
3019 DD(m_SupportRegion.GetIndex(InputDimension-1));
3021 } // end switch index
3025 // // ----------------------------------------------------------------------
3026 // // The egg with 5 internal CP:
3027 // // Periodic, C2 smooth everywhere and constrained to zero at the reference
3028 // // Coeff row R5 BC1 0 1 2 3 BC2 R2
3029 // // PaddedCoeff R: 0 1 2 3 4 5 6 7
3030 // // BC1= -weights[2]/weights[0] ( R2+R5)
3032 // // ---------------------------------------------------------------------
3033 // switch(m_SupportRegion.GetIndex(InputDimension-1))
3038 // m_FirstRegion.SetSize(InputDimension-1,1);
3039 // m_FirstRegion.SetIndex(InputDimension-1,3);
3042 // m_BCRegions.resize(2);
3043 // m_BCValues.resize(2);
3044 // m_BCRegion.SetIndex(InputDimension-1,0);
3045 // m_BCRegions[0]=m_BCRegion;
3046 // m_BCValues[0]=-m_WeightRatio[2][0];
3047 // m_BCRegion.SetIndex(InputDimension-1,3);
3048 // m_BCRegions[1]=m_BCRegion;
3049 // m_BCValues[1]=-m_WeightRatio[2][0];
3052 // m_SecondRegion.SetSize(InputDimension-1,2);
3053 // m_SecondRegion.SetIndex(InputDimension-1,0);
3056 // m_BC2Regions.resize(0);
3063 // m_FirstRegion.SetSize(InputDimension-1,0);
3066 // m_BCRegions.resize(2);
3067 // m_BCValues.resize(2);
3068 // m_BCRegion.SetIndex(InputDimension-1,0);
3069 // m_BCRegions[0]=m_BCRegion;
3070 // m_BCValues[0]=-m_WeightRatio[2][0];
3071 // m_BCRegion.SetIndex(InputDimension-1,3);
3072 // m_BCRegions[1]=m_BCRegion;
3073 // m_BCValues[1]=-m_WeightRatio[2][0];
3076 // m_SecondRegion.SetSize(InputDimension-1,3);
3077 // m_SecondRegion.SetIndex(InputDimension-1,0);
3080 // m_BC2Regions.resize(0);
3087 // m_FirstRegion.SetSize(InputDimension-1,4);
3088 // m_FirstRegion.SetIndex(InputDimension-1,0);
3091 // m_BCRegions.resize(0);
3094 // m_SecondRegion.SetSize(InputDimension-1,0);
3097 // m_BC2Regions.resize(0);
3104 // m_FirstRegion.SetSize(InputDimension-1,3);
3105 // m_FirstRegion.SetIndex(InputDimension-1,1);
3108 // m_BCRegions.resize(2);
3109 // m_BCValues.resize(2);
3110 // m_BCRegion.SetIndex(InputDimension-1,0);
3111 // m_BCRegions[0]=m_BCRegion;
3112 // m_BCValues[0]=-m_WeightRatio[2][0];
3113 // m_BCRegion.SetIndex(InputDimension-1,3);
3114 // m_BCRegions[1]=m_BCRegion;
3115 // m_BCValues[1]=-m_WeightRatio[2][0];
3118 // m_SecondRegion.SetSize(InputDimension-1,0);
3121 // m_BC2Regions.resize(0);
3129 // m_FirstRegion.SetSize(InputDimension-1,2);
3130 // m_FirstRegion.SetIndex(InputDimension-1,2);
3133 // m_BCRegions.resize(2);
3134 // m_BCValues.resize(2);
3135 // m_BCRegion.SetIndex(InputDimension-1,0);
3136 // m_BCRegions[0]=m_BCRegion;
3137 // m_BCValues[0]=-m_WeightRatio[2][0];
3138 // m_BCRegion.SetIndex(InputDimension-1,3);
3139 // m_BCRegions[1]=m_BCRegion;
3140 // m_BCValues[1]=-m_WeightRatio[2][0];
3143 // m_SecondRegion.SetSize(InputDimension-1,1);
3144 // m_SecondRegion.SetIndex(InputDimension-1,0);
3147 // m_BC2Regions.resize(0);
3154 // DD("supportindex > 4 ???");
3155 // DD(m_SupportRegion.GetIndex(InputDimension-1));
3156 // DD(m_TransformShape);
3158 // }// end swith index
3161 // } // end case 1 shape
3165 // ----------------------------------------------------------------------
3166 // The rabbit with 4 internal CP
3167 // Periodic, constrained to zero at the reference
3168 // at position 3 and the extremes fixed with anit-symm bc
3169 // Coeff row BC1 0 1 BC2 2 BC3 BC4
3170 // PaddedCoeff R: 0 1 2 3 4 5 6
3172 // BC2= -weights[2]/weights[0] ( R2+R4 )
3175 // ---------------------------------------------------------------------
3176 switch(m_SupportRegion.GetIndex(InputDimension-1))
3181 m_FirstRegion.SetSize(InputDimension-1,0);
3184 m_BCRegions.resize(2);
3185 m_BCValues.resize(2);
3186 m_BCRegion.SetIndex(InputDimension-1,0);
3187 m_BCRegions[0]=m_BCRegion;
3189 m_BCRegion.SetIndex(InputDimension-1,1);
3190 m_BCRegions[1]=m_BCRegion;
3194 m_SecondRegion.SetSize(InputDimension-1,2);
3195 m_SecondRegion.SetIndex(InputDimension-1,0);
3198 m_BC2Regions.resize(2);
3199 m_BC2Values.resize(2);
3200 m_BCRegion.SetIndex(InputDimension-1,1);
3201 m_BC2Regions[0]=m_BCRegion;
3202 m_BC2Values[0]=-m_WeightRatio[2][0];
3203 m_BCRegion.SetIndex(InputDimension-1,2);
3204 m_BC2Regions[1]=m_BCRegion;
3205 m_BC2Values[1]= -m_WeightRatio[2][0];
3212 m_FirstRegion.SetSize(InputDimension-1,2);
3213 m_FirstRegion.SetIndex(InputDimension-1,0);
3216 m_BCRegions.resize(2);
3217 m_BCValues.resize(2);
3218 m_BCRegion.SetIndex(InputDimension-1,1);
3219 m_BCRegions[0]=m_BCRegion;
3220 m_BCValues[0]=-m_WeightRatio[2][0];
3221 m_BCRegion.SetIndex(InputDimension-1,2);
3222 m_BCRegions[1]=m_BCRegion;
3223 m_BCValues[1]= -m_WeightRatio[2][0];
3226 m_SecondRegion.SetSize(InputDimension-1,1);
3227 m_SecondRegion.SetIndex(InputDimension-1,2);
3230 m_BC2Regions.resize(0);
3237 m_FirstRegion.SetSize(InputDimension-1,1);
3238 m_FirstRegion.SetIndex(InputDimension-1,1);
3241 m_BCRegions.resize(2);
3242 m_BCValues.resize(2);
3243 m_BCRegion.SetIndex(InputDimension-1,1);
3244 m_BCRegions[0]=m_BCRegion;
3245 m_BCValues[0]=-m_WeightRatio[2][0];
3246 m_BCRegion.SetIndex(InputDimension-1,2);
3247 m_BCRegions[1]=m_BCRegion;
3248 m_BCValues[1]= -m_WeightRatio[2][0];
3251 m_SecondRegion.SetSize(InputDimension-1,1);
3252 m_SecondRegion.SetIndex(InputDimension-1,2);
3255 m_BC2Regions.resize(1);
3256 m_BC2Values.resize(1);
3257 m_BCRegion.SetIndex(InputDimension-1,0);
3258 m_BC2Regions[0]=m_BCRegion;
3266 m_FirstRegion.SetSize(InputDimension-1,0);
3269 m_BCRegions.resize(2);
3270 m_BCValues.resize(2);
3271 m_BCRegion.SetIndex(InputDimension-1,1);
3272 m_BCRegions[0]=m_BCRegion;
3273 m_BCValues[0]=-m_WeightRatio[2][0];
3274 m_BCRegion.SetIndex(InputDimension-1,2);
3275 m_BCRegions[1]=m_BCRegion;
3276 m_BCValues[1]= -m_WeightRatio[2][0];
3279 m_SecondRegion.SetSize(InputDimension-1,1);
3280 m_SecondRegion.SetIndex(InputDimension-1,2);
3283 m_BC2Regions.resize(1);
3284 m_BC2Values.resize(1);
3285 m_BCRegion.SetIndex(InputDimension-1,0);
3286 m_BC2Regions[0]=m_BCRegion;
3290 m_BC3Regions.resize(2);
3291 m_BC3Values.resize(2);
3292 m_BCRegion.SetIndex(InputDimension-1,0);
3293 m_BC3Regions[0]=m_BCRegion;
3295 m_BCRegion.SetIndex(InputDimension-1,2);
3296 m_BC3Regions[1]=m_BCRegion;
3304 DD("supportindex > 3 ???");
3305 DD(m_SupportRegion.GetIndex(InputDimension-1));
3307 } // end switch index
3310 } // end case 2 shape
3313 // ----------------------------------------------------------------------
3314 // The rabbit with 5 internal CP
3315 // Periodic, constrained to zero at the reference at position 3.5
3316 // and the extremes fixed with anti-symmetrical BC
3317 // Coeff row BC1 0 1 BC2 2 3 BC3 BC4
3318 // PaddedCoeff R: 0 1 2 3 4 5 6 7
3320 // BC2= -weights[3]/weights[1] ( R2+R5 ) - R4
3323 // ---------------------------------------------------------------------
3324 switch(m_SupportRegion.GetIndex(InputDimension-1))
3329 m_FirstRegion.SetSize(InputDimension-1,0);
3332 m_BCRegions.resize(2);
3333 m_BCValues.resize(2);
3334 m_BCRegion.SetIndex(InputDimension-1,0);
3335 m_BCRegions[0]=m_BCRegion;
3337 m_BCRegion.SetIndex(InputDimension-1,1);
3338 m_BCRegions[1]=m_BCRegion;
3342 m_SecondRegion.SetSize(InputDimension-1,2);
3343 m_SecondRegion.SetIndex(InputDimension-1,0);
3346 m_BC2Regions.resize(3);
3347 m_BC2Values.resize(3);
3348 m_BCRegion.SetIndex(InputDimension-1,1);
3349 m_BC2Regions[0]=m_BCRegion;
3350 m_BC2Values[0]=-m_WeightRatio[3][1];
3351 m_BCRegion.SetIndex(InputDimension-1,2);
3352 m_BC2Regions[1]=m_BCRegion;
3354 m_BCRegion.SetIndex(InputDimension-1,3);
3355 m_BC2Regions[2]=m_BCRegion;
3356 m_BC2Values[2]=-m_WeightRatio[3][1];
3363 m_FirstRegion.SetSize(InputDimension-1,2);
3364 m_FirstRegion.SetIndex(InputDimension-1,0);
3367 m_BCRegions.resize(3);
3368 m_BCValues.resize(3);
3369 m_BCRegion.SetIndex(InputDimension-1,1);
3370 m_BCRegions[0]=m_BCRegion;
3371 m_BCValues[0]=-m_WeightRatio[3][1];
3372 m_BCRegion.SetIndex(InputDimension-1,2);
3373 m_BCRegions[1]=m_BCRegion;
3375 m_BCRegion.SetIndex(InputDimension-1,3);
3376 m_BCRegions[2]=m_BCRegion;
3377 m_BCValues[2]=-m_WeightRatio[3][1];
3380 m_SecondRegion.SetSize(InputDimension-1,1);
3381 m_SecondRegion.SetIndex(InputDimension-1,2);
3384 m_BC2Regions.resize(0);
3391 m_FirstRegion.SetSize(InputDimension-1,1);
3392 m_FirstRegion.SetIndex(InputDimension-1,1);
3395 m_BCRegions.resize(3);
3396 m_BCValues.resize(3);
3397 m_BCRegion.SetIndex(InputDimension-1,1);
3398 m_BCRegions[0]=m_BCRegion;
3399 m_BCValues[0]=-m_WeightRatio[3][1];
3400 m_BCRegion.SetIndex(InputDimension-1,2);
3401 m_BCRegions[1]=m_BCRegion;
3403 m_BCRegion.SetIndex(InputDimension-1,3);
3404 m_BCRegions[2]=m_BCRegion;
3405 m_BCValues[2]=-m_WeightRatio[3][1];
3408 m_SecondRegion.SetSize(InputDimension-1,2);
3409 m_SecondRegion.SetIndex(InputDimension-1,2);
3412 m_BC2Regions.resize(0);
3419 m_FirstRegion.SetSize(InputDimension-1,0);
3422 m_BCRegions.resize(3);
3423 m_BCValues.resize(3);
3424 m_BCRegion.SetIndex(InputDimension-1,1);
3425 m_BCRegions[0]=m_BCRegion;
3426 m_BCValues[0]=-m_WeightRatio[3][1];
3427 m_BCRegion.SetIndex(InputDimension-1,2);
3428 m_BCRegions[1]=m_BCRegion;
3430 m_BCRegion.SetIndex(InputDimension-1,3);
3431 m_BCRegions[2]=m_BCRegion;
3432 m_BCValues[2]=-m_WeightRatio[3][1];
3435 m_SecondRegion.SetSize(InputDimension-1,2);
3436 m_SecondRegion.SetIndex(InputDimension-1,2);
3439 m_BC2Regions.resize(1);
3440 m_BC2Values.resize(1);
3441 m_BCRegion.SetIndex(InputDimension-1,0);
3442 m_BC2Regions[0]=m_BCRegion;
3451 m_FirstRegion.SetSize(InputDimension-1,2);
3452 m_FirstRegion.SetIndex(InputDimension-1,2);
3455 m_BCRegions.resize(1);
3456 m_BCValues.resize(1);
3457 m_BCRegion.SetIndex(InputDimension-1,0);
3458 m_BCRegions[0]=m_BCRegion;
3462 m_SecondRegion.SetSize(InputDimension-1,0);
3465 m_BC2Regions.resize(2);
3466 m_BC2Values.resize(2);
3467 m_BCRegion.SetIndex(InputDimension-1,0);
3468 m_BC2Regions[0]=m_BCRegion;
3470 m_BCRegion.SetIndex(InputDimension-1,3);
3471 m_BC2Regions[1]=m_BCRegion;
3479 DD("supportindex > 4 ???");
3480 DD(m_SupportRegion.GetIndex(InputDimension-1));
3482 } // end switch index
3485 } // end case 3 shape
3489 // ----------------------------------------------------------------------
3490 // The sputnik with 4 internal CP
3491 // Periodic, constrained to zero at the reference
3492 // at position 3 and one indepent extremes copied
3493 // Coeff row BC1 0 1 BC2 2 BC2 3
3494 // PaddedCoeff R: 0 1 2 3 4 5 6
3496 // BC2= -weights[2]/weights[0] ( R2+R4 )
3497 // BC3= weights[2]/weights[0] ( R2-R4) + R1
3498 // ---------------------------------------------------------------------
3499 switch(m_SupportRegion.GetIndex(InputDimension-1))
3504 m_FirstRegion.SetSize(InputDimension-1,0);
3507 m_BCRegions.resize(1);
3508 m_BCValues.resize(1);
3509 m_BCRegion.SetIndex(InputDimension-1,3);
3510 m_BCRegions[0]=m_BCRegion;
3514 m_SecondRegion.SetSize(InputDimension-1,2);
3515 m_SecondRegion.SetIndex(InputDimension-1,0);
3518 m_BC2Regions.resize(2);
3519 m_BC2Values.resize(2);
3520 m_BCRegion.SetIndex(InputDimension-1,1);
3521 m_BC2Regions[0]=m_BCRegion;
3522 m_BC2Values[0]=-m_WeightRatio[2][0];
3523 m_BCRegion.SetIndex(InputDimension-1,2);
3524 m_BC2Regions[1]=m_BCRegion;
3525 m_BC2Values[1]= -m_WeightRatio[2][0];
3532 m_FirstRegion.SetSize(InputDimension-1,2);
3533 m_FirstRegion.SetIndex(InputDimension-1,0);
3536 m_BCRegions.resize(2);
3537 m_BCValues.resize(2);
3538 m_BCRegion.SetIndex(InputDimension-1,1);
3539 m_BCRegions[0]=m_BCRegion;
3540 m_BCValues[0]=-m_WeightRatio[2][0];
3541 m_BCRegion.SetIndex(InputDimension-1,2);
3542 m_BCRegions[1]=m_BCRegion;
3543 m_BCValues[1]= -m_WeightRatio[2][0];
3546 m_SecondRegion.SetSize(InputDimension-1,1);
3547 m_SecondRegion.SetIndex(InputDimension-1,2);
3550 m_BC2Regions.resize(0);
3557 m_FirstRegion.SetSize(InputDimension-1,1);
3558 m_FirstRegion.SetIndex(InputDimension-1,1);
3561 m_BCRegions.resize(2);
3562 m_BCValues.resize(2);
3563 m_BCRegion.SetIndex(InputDimension-1,1);
3564 m_BCRegions[0]=m_BCRegion;
3565 m_BCValues[0]=-m_WeightRatio[2][0];
3566 m_BCRegion.SetIndex(InputDimension-1,2);
3567 m_BCRegions[1]=m_BCRegion;
3568 m_BCValues[1]= -m_WeightRatio[2][0];
3571 m_SecondRegion.SetSize(InputDimension-1,1);
3572 m_SecondRegion.SetIndex(InputDimension-1,2);
3575 m_BC2Regions.resize(3);
3576 m_BC2Values.resize(3);
3577 m_BCRegion.SetIndex(InputDimension-1,0);
3578 m_BC2Regions[0]=m_BCRegion;
3580 m_BCRegion.SetIndex(InputDimension-1,1);
3581 m_BC2Regions[1]=m_BCRegion;
3582 m_BC2Values[1]=m_WeightRatio[2][0];
3583 m_BCRegion.SetIndex(InputDimension-1,2);
3584 m_BC2Regions[2]=m_BCRegion;
3585 m_BC2Values[2]=-m_WeightRatio[2][0];
3592 m_FirstRegion.SetSize(InputDimension-1,0);
3595 m_BCRegions.resize(2);
3596 m_BCValues.resize(2);
3597 m_BCRegion.SetIndex(InputDimension-1,1);
3598 m_BCRegions[0]=m_BCRegion;
3599 m_BCValues[0]=-m_WeightRatio[2][0];
3600 m_BCRegion.SetIndex(InputDimension-1,2);
3601 m_BCRegions[1]=m_BCRegion;
3602 m_BCValues[1]= -m_WeightRatio[2][0];
3605 m_SecondRegion.SetSize(InputDimension-1,1);
3606 m_SecondRegion.SetIndex(InputDimension-1,2);
3609 m_BC2Regions.resize(3);
3610 m_BC2Values.resize(3);
3611 m_BCRegion.SetIndex(InputDimension-1,0);
3612 m_BC2Regions[0]=m_BCRegion;
3614 m_BCRegion.SetIndex(InputDimension-1,1);
3615 m_BC2Regions[1]=m_BCRegion;
3616 m_BC2Values[1]=m_WeightRatio[2][0];
3617 m_BCRegion.SetIndex(InputDimension-1,2);
3618 m_BC2Regions[2]=m_BCRegion;
3619 m_BC2Values[2]=-m_WeightRatio[2][0];
3622 m_ThirdRegion.SetSize(InputDimension-1,1);
3623 m_ThirdRegion.SetIndex(InputDimension-1,3);
3630 DD("supportindex > 3 ???");
3631 DD(m_SupportRegion.GetIndex(InputDimension-1));
3633 } // end switch index
3636 } // end case 4 shape
3640 // ----------------------------------------------------------------------
3641 // The sputnik with 5 internal CP
3642 // Periodic, constrained to zero at the reference
3643 // at position 3 and one indepent extreme
3644 // Coeff row BC1 0 1 BC2 2 3 BC3 4
3645 // PaddedCoeff R: 0 1 2 3 4 5 6 7
3646 // BC1= R2 + R5 - R7
3647 // BC2= -weights[3]/weights[1] ( R2+R5 ) - R4
3648 // BC3= R1 + 0.5 R2 - 0.5 R7
3649 // ----------------------------------------------------------------------
3650 switch(m_SupportRegion.GetIndex(InputDimension-1))
3655 m_FirstRegion.SetSize(InputDimension-1,0);
3658 m_BCRegions.resize(3);
3659 m_BCValues.resize(3);
3660 m_BCRegion.SetIndex(InputDimension-1,1);
3661 m_BCRegions[0]=m_BCRegion;
3663 m_BCRegion.SetIndex(InputDimension-1,3);
3664 m_BCRegions[1]=m_BCRegion;
3666 m_BCRegion.SetIndex(InputDimension-1,4);
3667 m_BCRegions[2]=m_BCRegion;
3671 m_SecondRegion.SetSize(InputDimension-1,2);
3672 m_SecondRegion.SetIndex(InputDimension-1,0);
3675 m_BC2Regions.resize(3);
3676 m_BC2Values.resize(3);
3677 m_BCRegion.SetIndex(InputDimension-1,1);
3678 m_BC2Regions[0]=m_BCRegion;
3679 m_BC2Values[0]=-m_WeightRatio[3][1];
3680 m_BCRegion.SetIndex(InputDimension-1,2);
3681 m_BC2Regions[1]=m_BCRegion;
3683 m_BCRegion.SetIndex(InputDimension-1,3);
3684 m_BC2Regions[2]=m_BCRegion;
3685 m_BC2Values[2]=-m_WeightRatio[3][1];
3692 m_FirstRegion.SetSize(InputDimension-1,2);
3693 m_FirstRegion.SetIndex(InputDimension-1,0);
3696 m_BCRegions.resize(3);
3697 m_BCValues.resize(3);
3698 m_BCRegion.SetIndex(InputDimension-1,1);
3699 m_BCRegions[0]=m_BCRegion;
3700 m_BCValues[0]=-m_WeightRatio[3][1];
3701 m_BCRegion.SetIndex(InputDimension-1,2);
3702 m_BCRegions[1]=m_BCRegion;
3704 m_BCRegion.SetIndex(InputDimension-1,3);
3705 m_BCRegions[2]=m_BCRegion;
3706 m_BCValues[2]=-m_WeightRatio[3][1];
3709 m_SecondRegion.SetSize(InputDimension-1,1);
3710 m_SecondRegion.SetIndex(InputDimension-1,2);
3713 m_BC2Regions.resize(0);
3720 m_FirstRegion.SetSize(InputDimension-1,1);
3721 m_FirstRegion.SetIndex(InputDimension-1,1);
3724 m_BCRegions.resize(3);
3725 m_BCValues.resize(3);
3726 m_BCRegion.SetIndex(InputDimension-1,1);
3727 m_BCRegions[0]=m_BCRegion;
3728 m_BCValues[0]=-m_WeightRatio[3][1];
3729 m_BCRegion.SetIndex(InputDimension-1,2);
3730 m_BCRegions[1]=m_BCRegion;
3732 m_BCRegion.SetIndex(InputDimension-1,3);
3733 m_BCRegions[2]=m_BCRegion;
3734 m_BCValues[2]=-m_WeightRatio[3][1];
3737 m_SecondRegion.SetSize(InputDimension-1,2);
3738 m_SecondRegion.SetIndex(InputDimension-1,2);
3741 m_BC2Regions.resize(0);
3748 m_FirstRegion.SetSize(InputDimension-1,0);
3751 m_BCRegions.resize(3);
3752 m_BCValues.resize(3);
3753 m_BCRegion.SetIndex(InputDimension-1,1);
3754 m_BCRegions[0]=m_BCRegion;
3755 m_BCValues[0]=-m_WeightRatio[3][1];
3756 m_BCRegion.SetIndex(InputDimension-1,2);
3757 m_BCRegions[1]=m_BCRegion;
3759 m_BCRegion.SetIndex(InputDimension-1,3);
3760 m_BCRegions[2]=m_BCRegion;
3761 m_BCValues[2]=-m_WeightRatio[3][1];
3764 m_SecondRegion.SetSize(InputDimension-1,2);
3765 m_SecondRegion.SetIndex(InputDimension-1,2);
3768 m_BC2Regions.resize(3);
3769 m_BC2Values.resize(3);
3770 m_BCRegion.SetIndex(InputDimension-1,0);
3771 m_BC2Regions[0]=m_BCRegion;
3773 m_BCRegion.SetIndex(InputDimension-1,1);
3774 m_BC2Regions[1]=m_BCRegion;
3776 m_BCRegion.SetIndex(InputDimension-1,4);
3777 m_BC2Regions[2]=m_BCRegion;
3778 m_BC2Values[2]=-0.5;
3785 m_FirstRegion.SetSize(InputDimension-1,2);
3786 m_FirstRegion.SetIndex(InputDimension-1,2);
3789 m_BCRegions.resize(3);
3790 m_BCValues.resize(3);
3791 m_BCRegion.SetIndex(InputDimension-1,0);
3792 m_BCRegions[0]=m_BCRegion;
3794 m_BCRegion.SetIndex(InputDimension-1,1);
3795 m_BCRegions[1]=m_BCRegion;
3797 m_BCRegion.SetIndex(InputDimension-1,4);
3798 m_BCRegions[2]=m_BCRegion;
3802 m_SecondRegion.SetSize(InputDimension-1,1);
3803 m_SecondRegion.SetIndex(InputDimension-1,4);
3806 m_BC2Regions.resize(0);
3812 DD("supportindex > 4 ???");
3813 DD(m_SupportRegion.GetIndex(InputDimension-1));
3815 } // end switch index
3818 } // end case 5 shape
3825 // ----------------------------------------------------------------------
3826 // The diamond with 4 internal CP:
3827 // Periodic, constrained to zero at the reference at position 3
3828 // Coeff row 0 1 2 BC1 3 BC2 4
3829 // PaddedCoeff R:0 1 2 3 4 5 6
3830 // BC1= -weights[2]/weights[0] ( R2+R4 )
3831 // BC2= weights[2]/weights[0] ( R0+R2-R4-R6 ) + R1
3832 // ---------------------------------------------------------------------
3833 switch(m_SupportRegion.GetIndex(InputDimension-1))
3838 m_FirstRegion.SetSize(InputDimension-1,3);
3839 m_FirstRegion.SetIndex(InputDimension-1,0);
3842 m_BCRegions.resize(2);
3843 m_BCValues.resize(2);
3844 m_BCRegion.SetIndex(InputDimension-1,2);
3845 m_BCRegions[0]=m_BCRegion;
3846 m_BCValues[0]=-m_WeightRatio[2][0];
3847 m_BCRegion.SetIndex(InputDimension-1,3);
3848 m_BCRegions[1]=m_BCRegion;
3849 m_BCValues[1]=-m_WeightRatio[2][0];
3852 m_SecondRegion.SetSize(InputDimension-1,0);
3855 m_BC2Regions.resize(0);
3862 m_FirstRegion.SetSize(InputDimension-1,2);
3863 m_FirstRegion.SetIndex(InputDimension-1,1);
3866 m_BCRegions.resize(2);
3867 m_BCValues.resize(2);
3868 m_BCRegion.SetIndex(InputDimension-1,2);
3869 m_BCRegions[0]=m_BCRegion;
3870 m_BCValues[0]=-m_WeightRatio[2][0];
3871 m_BCRegion.SetIndex(InputDimension-1,3);
3872 m_BCRegions[1]=m_BCRegion;
3873 m_BCValues[1]=-m_WeightRatio[2][0];
3876 m_SecondRegion.SetSize(InputDimension-1,1);
3877 m_SecondRegion.SetIndex(InputDimension-1,3);
3880 m_BC2Regions.resize(0);
3887 m_FirstRegion.SetSize(InputDimension-1,1);
3888 m_FirstRegion.SetIndex(InputDimension-1,2);
3891 m_BCRegions.resize(2);
3892 m_BCValues.resize(2);
3893 m_BCRegion.SetIndex(InputDimension-1,2);
3894 m_BCRegions[0]=m_BCRegion;
3895 m_BCValues[0]=-m_WeightRatio[2][0];
3896 m_BCRegion.SetIndex(InputDimension-1,3);
3897 m_BCRegions[1]=m_BCRegion;
3898 m_BCValues[1]=-m_WeightRatio[2][0];
3901 m_SecondRegion.SetSize(InputDimension-1,1);
3902 m_SecondRegion.SetIndex(InputDimension-1,3);
3905 m_BC2Regions.resize(5);
3906 m_BC2Values.resize(5);
3907 m_BCRegion.SetIndex(InputDimension-1,0);
3908 m_BC2Regions[0]=m_BCRegion;
3909 m_BC2Values[0]=m_WeightRatio[2][0];
3910 m_BCRegion.SetIndex(InputDimension-1,1);
3911 m_BC2Regions[1]=m_BCRegion;
3913 m_BCRegion.SetIndex(InputDimension-1,2);
3914 m_BC2Regions[2]=m_BCRegion;
3915 m_BC2Values[2]=m_WeightRatio[2][0];
3916 m_BCRegion.SetIndex(InputDimension-1,3);
3917 m_BC2Regions[3]=m_BCRegion;
3918 m_BC2Values[3]=-m_WeightRatio[2][0];
3919 m_BCRegion.SetIndex(InputDimension-1,4);
3920 m_BC2Regions[4]=m_BCRegion;
3921 m_BC2Values[4]=-m_WeightRatio[2][0];
3928 m_FirstRegion.SetSize(InputDimension-1,0);
3931 m_BCRegions.resize(2);
3932 m_BCValues.resize(2);
3933 m_BCRegion.SetIndex(InputDimension-1,2);
3934 m_BCRegions[0]=m_BCRegion;
3935 m_BCValues[0]=-m_WeightRatio[2][0];
3936 m_BCRegion.SetIndex(InputDimension-1,3);
3937 m_BCRegions[1]=m_BCRegion;
3938 m_BCValues[1]=-m_WeightRatio[2][0];
3941 m_SecondRegion.SetSize(InputDimension-1,1);
3942 m_SecondRegion.SetIndex(InputDimension-1,3);
3945 m_BC2Regions.resize(5);
3946 m_BC2Values.resize(5);
3947 m_BCRegion.SetIndex(InputDimension-1,0);
3948 m_BC2Regions[0]=m_BCRegion;
3949 m_BC2Values[0]=m_WeightRatio[2][0];
3950 m_BCRegion.SetIndex(InputDimension-1,1);
3951 m_BC2Regions[1]=m_BCRegion;
3953 m_BCRegion.SetIndex(InputDimension-1,2);
3954 m_BC2Regions[2]=m_BCRegion;
3955 m_BC2Values[2]=m_WeightRatio[2][0];
3956 m_BCRegion.SetIndex(InputDimension-1,3);
3957 m_BC2Regions[3]=m_BCRegion;
3958 m_BC2Values[3]=-m_WeightRatio[2][0];
3959 m_BCRegion.SetIndex(InputDimension-1,4);
3960 m_BC2Regions[4]=m_BCRegion;
3961 m_BC2Values[4]=-m_WeightRatio[2][0];
3964 m_ThirdRegion.SetSize(InputDimension-1,1);
3965 m_ThirdRegion.SetIndex(InputDimension-1,4);
3972 DD("supportindex > 3 ???");
3973 DD(m_SupportRegion.GetIndex(InputDimension-1));
3975 } // end switch index
3978 } // end case 7 shape
3982 // ----------------------------------------------------------------------
3983 // The diamond with 5 internal CP:
3984 // periodic, constrained to zero at the reference at position 3.5
3985 // Coeff row 0 1 2 BC1 3 4 BC2 5
3986 // PaddedCoeff R:0 1 2 3 4 5 6 7
3987 // BC1= -weights[3]/weights[1] ( R2+R5 ) - R4
3988 // BC2= weights[2]/weights[0] ( R0+R2-R5-R7 ) + R1
3989 // ---------------------------------------------------------------------
3990 switch(m_SupportRegion.GetIndex(InputDimension-1))
3995 m_FirstRegion.SetSize(InputDimension-1,3);
3996 m_FirstRegion.SetIndex(InputDimension-1,0);
3999 m_BCRegions.resize(3);
4000 m_BCValues.resize(3);
4001 m_BCRegion.SetIndex(InputDimension-1,2);
4002 m_BCRegions[0]=m_BCRegion;
4003 m_BCValues[0]=-m_WeightRatio[3][1];
4004 m_BCRegion.SetIndex(InputDimension-1,3);
4005 m_BCRegions[1]=m_BCRegion;
4007 m_BCRegion.SetIndex(InputDimension-1,4);
4008 m_BCRegions[2]=m_BCRegion;
4009 m_BCValues[2]=-m_WeightRatio[3][1];
4012 m_SecondRegion.SetSize(InputDimension-1,0);
4015 m_BC2Regions.resize(0);
4022 m_FirstRegion.SetSize(InputDimension-1,2);
4023 m_FirstRegion.SetIndex(InputDimension-1,1);
4026 m_BCRegions.resize(3);
4027 m_BCValues.resize(3);
4028 m_BCRegion.SetIndex(InputDimension-1,2);
4029 m_BCRegions[0]=m_BCRegion;
4030 m_BCValues[0]=-m_WeightRatio[3][1];
4031 m_BCRegion.SetIndex(InputDimension-1,3);
4032 m_BCRegions[1]=m_BCRegion;
4034 m_BCRegion.SetIndex(InputDimension-1,4);
4035 m_BCRegions[2]=m_BCRegion;
4036 m_BCValues[2]=-m_WeightRatio[3][1];
4039 m_SecondRegion.SetSize(InputDimension-1,1);
4040 m_SecondRegion.SetIndex(InputDimension-1,3);
4043 m_BC2Regions.resize(0);
4050 m_FirstRegion.SetSize(InputDimension-1,1);
4051 m_FirstRegion.SetIndex(InputDimension-1,2);
4054 m_BCRegions.resize(3);
4055 m_BCValues.resize(3);
4056 m_BCRegion.SetIndex(InputDimension-1,2);
4057 m_BCRegions[0]=m_BCRegion;
4058 m_BCValues[0]=-m_WeightRatio[3][1];
4059 m_BCRegion.SetIndex(InputDimension-1,3);
4060 m_BCRegions[1]=m_BCRegion;
4062 m_BCRegion.SetIndex(InputDimension-1,4);
4063 m_BCRegions[2]=m_BCRegion;
4064 m_BCValues[2]=-m_WeightRatio[3][1];
4067 m_SecondRegion.SetSize(InputDimension-1,2);
4068 m_SecondRegion.SetIndex(InputDimension-1,3);
4071 m_BC2Regions.resize(0);
4078 m_FirstRegion.SetSize(InputDimension-1,0);
4079 m_FirstRegion.SetIndex(InputDimension-1,0);
4082 m_BCRegions.resize(3);
4083 m_BCValues.resize(3);
4084 m_BCRegion.SetIndex(InputDimension-1,2);
4085 m_BCRegions[0]=m_BCRegion;
4086 m_BCValues[0]=-m_WeightRatio[3][1];
4087 m_BCRegion.SetIndex(InputDimension-1,3);
4088 m_BCRegions[1]=m_BCRegion;
4090 m_BCRegion.SetIndex(InputDimension-1,4);
4091 m_BCRegions[2]=m_BCRegion;
4092 m_BCValues[2]=-m_WeightRatio[3][1];
4095 m_SecondRegion.SetSize(InputDimension-1,2);
4096 m_SecondRegion.SetIndex(InputDimension-1,3);
4099 m_BC2Regions.resize(5);
4100 m_BC2Values.resize(5);
4101 m_BCRegion.SetIndex(InputDimension-1,0);
4102 m_BC2Regions[0]=m_BCRegion;
4103 m_BC2Values[0]=m_WeightRatio[2][0];
4104 m_BCRegion.SetIndex(InputDimension-1,1);
4105 m_BC2Regions[1]=m_BCRegion;
4107 m_BCRegion.SetIndex(InputDimension-1,2);
4108 m_BC2Regions[2]=m_BCRegion;
4109 m_BC2Values[2]=m_WeightRatio[2][0];
4110 m_BCRegion.SetIndex(InputDimension-1,4);
4111 m_BC2Regions[3]=m_BCRegion;
4112 m_BC2Values[3]=-m_WeightRatio[2][0];
4113 m_BCRegion.SetIndex(InputDimension-1,5);
4114 m_BC2Regions[4]=m_BCRegion;
4115 m_BC2Values[4]=-m_WeightRatio[2][0];
4123 m_FirstRegion.SetSize(InputDimension-1,2);
4124 m_FirstRegion.SetIndex(InputDimension-1,3);
4127 m_BCRegions.resize(5);
4128 m_BCValues.resize(5);
4129 m_BCRegion.SetIndex(InputDimension-1,0);
4130 m_BCRegions[0]=m_BCRegion;
4131 m_BCValues[0]=m_WeightRatio[2][0];
4132 m_BCRegion.SetIndex(InputDimension-1,1);
4133 m_BCRegions[1]=m_BCRegion;
4135 m_BCRegion.SetIndex(InputDimension-1,2);
4136 m_BCRegions[2]=m_BCRegion;
4137 m_BCValues[2]=m_WeightRatio[2][0];
4138 m_BCRegion.SetIndex(InputDimension-1,4);
4139 m_BCRegions[3]=m_BCRegion;
4140 m_BCValues[3]=-m_WeightRatio[2][0];
4141 m_BCRegion.SetIndex(InputDimension-1,5);
4142 m_BCRegions[4]=m_BCRegion;
4143 m_BCValues[4]=-m_WeightRatio[2][0];
4146 m_SecondRegion.SetSize(InputDimension-1,1);
4147 m_SecondRegion.SetIndex(InputDimension-1,5);
4150 m_BC2Regions.resize(0);
4157 DD("supportindex > 4 ???");
4158 DD(m_SupportRegion.GetIndex(InputDimension-1));
4160 } // end switch index
4163 } // end case 7 shape
4167 // ----------------------------------------------------------------------
4168 // The sputnik with 5 internal CP
4169 // Periodic, constrained to zero at the reference
4170 // at position 3 and one indepent extreme
4171 // Coeff row BC1 0 1 BC2 2 3 BC3 4
4172 // PaddedCoeff R: 0 1 2 3 4 5 6 7
4174 // BC2= -weights[3]/weights[1] ( R2+R5 ) - R4
4176 // ---------------------------------------------------------------------
4177 switch(m_SupportRegion.GetIndex(InputDimension-1))
4182 m_FirstRegion.SetSize(InputDimension-1,0);
4185 m_BCRegions.resize(3);
4186 m_BCValues.resize(3);
4187 m_BCRegion.SetIndex(InputDimension-1,1);
4188 m_BCRegions[0]=m_BCRegion;
4190 m_BCRegion.SetIndex(InputDimension-1,3);
4191 m_BCRegions[1]=m_BCRegion;
4193 m_BCRegion.SetIndex(InputDimension-1,4);
4194 m_BCRegions[2]=m_BCRegion;
4198 m_SecondRegion.SetSize(InputDimension-1,2);
4199 m_SecondRegion.SetIndex(InputDimension-1,0);
4202 m_BC2Regions.resize(3);
4203 m_BC2Values.resize(3);
4204 m_BCRegion.SetIndex(InputDimension-1,1);
4205 m_BC2Regions[0]=m_BCRegion;
4206 m_BC2Values[0]=-m_WeightRatio[3][1];
4207 m_BCRegion.SetIndex(InputDimension-1,2);
4208 m_BC2Regions[1]=m_BCRegion;
4210 m_BCRegion.SetIndex(InputDimension-1,3);
4211 m_BC2Regions[2]=m_BCRegion;
4212 m_BC2Values[2]=-m_WeightRatio[3][1];
4219 m_FirstRegion.SetSize(InputDimension-1,2);
4220 m_FirstRegion.SetIndex(InputDimension-1,0);
4223 m_BCRegions.resize(3);
4224 m_BCValues.resize(3);
4225 m_BCRegion.SetIndex(InputDimension-1,1);
4226 m_BCRegions[0]=m_BCRegion;
4227 m_BCValues[0]=-m_WeightRatio[3][1];
4228 m_BCRegion.SetIndex(InputDimension-1,2);
4229 m_BCRegions[1]=m_BCRegion;
4231 m_BCRegion.SetIndex(InputDimension-1,3);
4232 m_BCRegions[2]=m_BCRegion;
4233 m_BCValues[2]=-m_WeightRatio[3][1];
4236 m_SecondRegion.SetSize(InputDimension-1,1);
4237 m_SecondRegion.SetIndex(InputDimension-1,2);
4240 m_BC2Regions.resize(0);
4247 m_FirstRegion.SetSize(InputDimension-1,1);
4248 m_FirstRegion.SetIndex(InputDimension-1,1);
4251 m_BCRegions.resize(3);
4252 m_BCValues.resize(3);
4253 m_BCRegion.SetIndex(InputDimension-1,1);
4254 m_BCRegions[0]=m_BCRegion;
4255 m_BCValues[0]=-m_WeightRatio[3][1];
4256 m_BCRegion.SetIndex(InputDimension-1,2);
4257 m_BCRegions[1]=m_BCRegion;
4259 m_BCRegion.SetIndex(InputDimension-1,3);
4260 m_BCRegions[2]=m_BCRegion;
4261 m_BCValues[2]=-m_WeightRatio[3][1];
4264 m_SecondRegion.SetSize(InputDimension-1,2);
4265 m_SecondRegion.SetIndex(InputDimension-1,2);
4268 m_BC2Regions.resize(0);
4275 m_FirstRegion.SetSize(InputDimension-1,0);
4278 m_BCRegions.resize(3);
4279 m_BCValues.resize(3);
4280 m_BCRegion.SetIndex(InputDimension-1,1);
4281 m_BCRegions[0]=m_BCRegion;
4282 m_BCValues[0]=-m_WeightRatio[3][1];
4283 m_BCRegion.SetIndex(InputDimension-1,2);
4284 m_BCRegions[1]=m_BCRegion;
4286 m_BCRegion.SetIndex(InputDimension-1,3);
4287 m_BCRegions[2]=m_BCRegion;
4288 m_BCValues[2]=-m_WeightRatio[3][1];
4291 m_SecondRegion.SetSize(InputDimension-1,1);
4292 m_SecondRegion.SetIndex(InputDimension-1,2);
4295 m_BC2Regions.resize(1);
4296 m_BC2Values.resize(1);
4297 m_BCRegion.SetIndex(InputDimension-1,0);
4298 m_BC2Regions[0]=m_BCRegion;
4306 m_FirstRegion.SetSize(InputDimension-1,2);
4307 m_FirstRegion.SetIndex(InputDimension-1,2);
4310 m_BCRegions.resize(1);
4311 m_BCValues.resize(1);
4312 m_BCRegion.SetIndex(InputDimension-1,0);
4313 m_BCRegions[0]=m_BCRegion;
4317 m_SecondRegion.SetSize(InputDimension-1,1);
4318 m_SecondRegion.SetIndex(InputDimension-1,4);
4321 m_BC2Regions.resize(0);
4327 DD("supportindex > 4 ???");
4328 DD(m_SupportRegion.GetIndex(InputDimension-1));
4330 } // end switch index
4333 } // end case 9 shape
4338 DD ("Other shapes currently not implemented");
4341 } // end switch shape
4342 } // end wrap region
4346 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
4348 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions,NOutputDimensions>
4349 ::TransformPointToContinuousIndex( const InputPointType & point, ContinuousIndexType & index ) const
4353 itk::Vector<double, OutputDimension> tvector;
4355 for ( j = 0; j < OutputDimension; j++ )
4357 //JV find the index in the PADDED version
4358 tvector[j] = point[j] - this->m_PaddedGridOrigin[j];
4361 itk::Vector<double, OutputDimension> cvector;
4363 cvector = m_PointToIndex * tvector;
4365 for ( j = 0; j < OutputDimension; j++ )
4367 index[j] = static_cast< typename ContinuousIndexType::CoordRepType >( cvector[j] );