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>
387 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
388 ::GetNumberOfParameters(void) const
391 // The number of parameters equal SpaceDimension * number of
392 // of pixels in the grid region.
393 return ( static_cast<unsigned int>( SpaceDimension ) *
394 static_cast<unsigned int>( m_GridRegion.GetNumberOfPixels() ) );
399 // Get the padded number of parameters
400 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
402 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
403 ::GetPaddedNumberOfParameters(void) const
406 // The number of parameters equal SpaceDimension * number of
407 // of pixels in the grid region.
408 return ( static_cast<unsigned int>( SpaceDimension ) *
409 static_cast<unsigned int>( m_PaddedGridRegion.GetNumberOfPixels() ) );
415 // Get the number of parameters per dimension
416 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
418 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
419 ::GetNumberOfParametersPerDimension(void) const
421 // The number of parameters per dimension equal number of
422 // of pixels in the grid region.
423 return ( static_cast<unsigned int>( m_GridRegion.GetNumberOfPixels() ) );
428 // Set the grid region
429 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
431 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
432 ::SetGridRegion( const RegionType & region )
434 if ( m_GridRegion != region )
436 m_GridRegion = region;
437 m_PaddedGridRegion=region;
439 // JV set the padded region
440 typename CoefficientImageType::RegionType::SizeType paddedSize= region.GetSize();
442 //JV size dependes on shape
443 switch (m_TransformShape)
447 paddedSize[InputDimension-1]+=4;
451 paddedSize[InputDimension-1]+=4;
455 paddedSize[InputDimension-1]+=3;
459 paddedSize[InputDimension-1]+=2;
463 paddedSize[InputDimension-1]+=3;
466 paddedSize[InputDimension-1]=+1;
468 m_PaddedGridRegion.SetSize(paddedSize);
470 // Set regions for each coefficient and jacobian image
471 m_WrappedImage->SetRegions( m_GridRegion );
472 m_PaddedCoefficientImage->SetRegions( m_PaddedGridRegion );
473 m_PaddedCoefficientImage->Allocate();
474 for (unsigned int j=0; j <OutputDimension;j++)
476 m_JacobianImage[j]->SetRegions( m_GridRegion );
479 // JV used the padded version for the valid region
480 // Set the valid region
481 // If the grid spans the interval [start, last].
482 // The valid interval for evaluation is [start+offset, last-offset]
483 // when spline order is even.
484 // The valid interval for evaluation is [start+offset, last-offset)
485 // when spline order is odd.
486 // Where offset = vcl_floor(spline / 2 ).
487 // Note that the last pixel is not included in the valid region
488 // with odd spline orders.
489 typename RegionType::SizeType size = m_PaddedGridRegion.GetSize();
490 typename RegionType::IndexType index = m_PaddedGridRegion.GetIndex();
491 for ( unsigned int j = 0; j < NInputDimensions; j++ )
494 static_cast< typename RegionType::IndexValueType >( m_Offset[j] );
496 static_cast< typename RegionType::SizeValueType> ( 2 * m_Offset[j] );
497 m_ValidRegionLast[j] = index[j] +
498 static_cast< typename RegionType::IndexValueType >( size[j] ) - 1;
500 m_ValidRegion.SetSize( size );
501 m_ValidRegion.SetIndex( index );
503 // If we are using the default parameters, update their size and set to identity.
504 // Input parameters point to internal buffer => using default parameters.
505 if (m_InputParametersPointer == &m_InternalParametersBuffer)
507 // Check if we need to resize the default parameter buffer.
508 if ( m_InternalParametersBuffer.GetSize() != this->GetNumberOfParameters() )
510 m_InternalParametersBuffer.SetSize( this->GetNumberOfParameters() );
511 // Fill with zeros for identity.
512 m_InternalParametersBuffer.Fill( 0 );
521 // Set the grid spacing
522 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
524 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
525 ::SetGridSpacing( const SpacingType & spacing )
527 if ( m_GridSpacing != spacing )
529 m_GridSpacing = spacing;
531 // Set spacing for each coefficient and jacobian image
532 m_WrappedImage->SetSpacing( m_GridSpacing.GetDataPointer() );
533 m_PaddedCoefficientImage->SetSpacing( m_GridSpacing.GetDataPointer() );
534 for (unsigned int j=0; j <OutputDimension; j++) m_JacobianImage[j]->SetSpacing( m_GridSpacing.GetDataPointer() );
538 for( unsigned int i=0; i<OutputDimension; i++)
540 scale[i][i] = m_GridSpacing[i];
543 m_IndexToPoint = m_GridDirection * scale;
544 m_PointToIndex = m_IndexToPoint.GetInverse();
551 // Set the grid direction
552 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
554 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
555 ::SetGridDirection( const DirectionType & direction )
557 if ( m_GridDirection != direction )
559 m_GridDirection = direction;
561 // Set direction for each coefficient and jacobian image
562 m_WrappedImage->SetDirection( m_GridDirection );
563 m_PaddedCoefficientImage->SetDirection( m_GridDirection );
564 for (unsigned int j=0; j <OutputDimension; j++) m_JacobianImage[j]->SetDirection( m_GridDirection );
568 for( unsigned int i=0; i<OutputDimension; i++)
570 scale[i][i] = m_GridSpacing[i];
573 m_IndexToPoint = m_GridDirection * scale;
574 m_PointToIndex = m_IndexToPoint.GetInverse();
582 // Set the grid origin
583 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
585 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
586 ::SetGridOrigin( const OriginType& origin )
588 if( m_GridOrigin!=origin)
590 m_GridOrigin = origin;
592 // JV The origin depends on the shape
593 switch (m_TransformShape)
597 m_PaddedGridOrigin=origin;
598 m_PaddedGridOrigin[InputDimension-1]=origin[InputDimension-1]-2* m_GridSpacing[InputDimension-1];
602 m_PaddedGridOrigin=origin;
603 m_PaddedGridOrigin[InputDimension-1]=origin[InputDimension-1]-1* m_GridSpacing[InputDimension-1];
607 m_PaddedGridOrigin=origin;
608 m_PaddedGridOrigin[InputDimension-1]=origin[InputDimension-1]-1* m_GridSpacing[InputDimension-1];
612 m_PaddedGridOrigin=origin;
616 m_PaddedGridOrigin=origin;
617 m_PaddedGridOrigin[InputDimension-1]=origin[InputDimension-1]-1* m_GridSpacing[InputDimension-1];
620 m_PaddedGridOrigin=origin;
623 // Set origin for each coefficient and jacobianimage
624 m_WrappedImage->SetOrigin( m_GridOrigin.GetDataPointer() );
625 m_PaddedCoefficientImage->SetOrigin( m_PaddedGridOrigin.GetDataPointer() );
626 for (unsigned int j=0; j <OutputDimension; j++) m_JacobianImage[j]->SetOrigin( m_GridOrigin.GetDataPointer() );
633 // Set the parameters
634 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
636 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
639 if( m_InputParametersPointer )
641 ParametersType * parameters =
642 const_cast<ParametersType *>( m_InputParametersPointer );
643 parameters->Fill( 0.0 );
648 itkExceptionMacro( << "Input parameters for the spline haven't been set ! "
649 << "Set them using the SetParameters or SetCoefficientImage method first." );
654 // Set the parameters
655 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
657 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
658 ::SetParameters( const ParametersType & parameters )
661 // Check if the number of parameters match the
662 // Expected number of parameters
663 if ( parameters.Size() != this->GetNumberOfParameters() )
665 itkExceptionMacro(<<"Mismatched between parameters size "
667 << " and region size "
668 << m_GridRegion.GetNumberOfPixels() );
671 // Clean up buffered parameters
672 m_InternalParametersBuffer = ParametersType( 0 );
674 // Keep a reference to the input parameters
675 m_InputParametersPointer = ¶meters;
677 // Wrap flat array as images of coefficients
678 this->WrapAsImages();
680 //JV Set padded input to vector interpolator
681 m_VectorInterpolator->SetInputImage(this->GetPaddedCoefficientImage());
683 // Modified is always called since we just have a pointer to the
684 // parameters and cannot know if the parameters have changed.
689 // Set the Fixed Parameters
690 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
692 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
693 ::SetFixedParameters( const ParametersType & parameters )
696 // JV number should be exact, no defaults for spacing
697 if ( parameters.Size() != NInputDimensions * (5 + NInputDimensions)+3 )
699 itkExceptionMacro(<< "Mismatched between parameters size "
701 << " and number of fixed parameters "
702 << NInputDimensions * (5 + NInputDimensions)+3 );
704 /*********************************************************
705 Fixed Parameters store the following information:
710 // JV we add the splineOrders, LUTsamplingfactor, mask pointer and bulktransform pointer
718 The size of these is equal to the NInputDimensions
719 *********************************************************/
721 /** Set the Grid Parameters */
723 for (unsigned int i=0; i<NInputDimensions; i++)
725 gridSize[i] = static_cast<int> (parameters[i]);
727 RegionType bsplineRegion;
728 bsplineRegion.SetSize( gridSize );
730 /** Set the Origin Parameters */
732 for (unsigned int i=0; i<NInputDimensions; i++)
734 origin[i] = parameters[NInputDimensions+i];
737 /** Set the Spacing Parameters */
739 for (unsigned int i=0; i<NInputDimensions; i++)
741 spacing[i] = parameters[2*NInputDimensions+i];
744 /** Set the Direction Parameters */
745 DirectionType direction;
746 for (unsigned int di=0; di<NInputDimensions; di++)
748 for (unsigned int dj=0; dj<NInputDimensions; dj++)
750 direction[di][dj] = parameters[3*NInputDimensions+(di*NInputDimensions+dj)];
754 //JV add the spline orders
755 SizeType splineOrders;
756 for (unsigned int i=0; i<NInputDimensions; i++)
758 splineOrders[i]=parameters[(3+NInputDimensions)*NInputDimensions+i];
761 //JV add the LUT sampling factor
762 SizeType samplingFactors;
763 for (unsigned int i=0; i<NInputDimensions; i++)
765 samplingFactors[i]=parameters[(4+NInputDimensions)*NInputDimensions+i];
768 //JV add the MaskPointer
769 m_Mask=((MaskPointer)( (size_t)parameters[(5+NInputDimensions)*NInputDimensions]));
771 //JV add the MaskPointer
772 m_BulkTransform=((BulkTransformPointer)( (size_t)parameters[(5+NInputDimensions)*NInputDimensions+1]));
774 //JV add the TransformShape
775 m_TransformShape=((unsigned int)parameters[(5+NInputDimensions)*NInputDimensions+2]);
778 this->SetSplineOrders( splineOrders );
779 this->SetGridSpacing( spacing );
780 this->SetGridDirection( direction );
781 this->SetGridOrigin( origin );
782 this->SetGridRegion( bsplineRegion );
783 this->SetLUTSamplingFactors( samplingFactors );
788 // Wrap flat parameters as images
789 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
791 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
794 //JV Wrap parameter array in vectorial image, changed parameter order: A1x A1y A1z, A2x ....
795 PixelType * dataPointer =reinterpret_cast<PixelType *>( const_cast<double *>(m_InputParametersPointer->data_block() )) ;
796 unsigned int numberOfPixels = m_GridRegion.GetNumberOfPixels();
798 m_WrappedImage->GetPixelContainer()->SetImportPointer( dataPointer,numberOfPixels);//InputDimension
799 m_CoefficientImage = m_WrappedImage;
801 //=====================================
802 //JV Create padded structure adding BC
803 //=====================================
804 PadCoefficientImage();
806 //=====================================
807 //JV Wrap jacobian into OutputDimension X Vectorial images
808 //=====================================
809 #if ITK_VERSION_MAJOR >= 4
810 this->m_SharedDataBSplineJacobian.set_size( OutputDimension, this->GetNumberOfParameters() );
812 this->m_Jacobian.set_size( OutputDimension, this->GetNumberOfParameters() );
815 // Use memset to set the memory
816 // JV four rows of three comps of parameters
817 #if ITK_VERSION_MAJOR >= 4
818 JacobianPixelType * jacobianDataPointer = reinterpret_cast<JacobianPixelType *>(this->m_SharedDataBSplineJacobian.data_block());
820 JacobianPixelType * jacobianDataPointer = reinterpret_cast<JacobianPixelType *>(this->m_Jacobian.data_block());
822 memset(jacobianDataPointer, 0, OutputDimension*numberOfPixels*sizeof(JacobianPixelType));
824 for (unsigned int j=0; j<OutputDimension; j++)
826 m_JacobianImage[j]->GetPixelContainer()->
827 SetImportPointer( jacobianDataPointer, numberOfPixels );
828 jacobianDataPointer += numberOfPixels;
831 // Reset the J parameters
832 m_LastJacobianIndex = m_ValidRegion.GetIndex();
833 m_FirstRegion.SetSize(m_NullSize);
834 m_SecondRegion.SetSize(m_NullSize);
835 for ( unsigned int j = 0; j < SpaceDimension ; j++ )
837 m_FirstIterator[j]= IteratorType( m_JacobianImage[j], m_FirstRegion);
838 m_SecondIterator[j]= IteratorType( m_JacobianImage[j], m_SecondRegion);
841 m_BCValues.resize(0);
842 m_BCRegions.resize(0);
844 m_BC2Values.resize(0);
845 m_BC2Regions.resize(0);
847 m_BC3Values.resize(0);
848 m_BC3Regions.resize(0);
854 // Set the parameters by value
855 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
857 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
858 ::SetParametersByValue( const ParametersType & parameters )
861 // Check if the number of parameters match the
862 // Expected number of parameters
863 if ( parameters.Size() != this->GetNumberOfParameters() )
865 itkExceptionMacro(<<"Mismatched between parameters size "
867 << " and region size "
868 << m_GridRegion.GetNumberOfPixels() );
872 m_InternalParametersBuffer = parameters;
873 m_InputParametersPointer = &m_InternalParametersBuffer;
875 // Wrap flat array as images of coefficients
876 this->WrapAsImages();
878 //JV Set padded input to vector interpolator
879 m_VectorInterpolator->SetInputImage(this->GetPaddedCoefficientImage());
881 // Modified is always called since we just have a pointer to the
882 // Parameters and cannot know if the parameters have changed.
887 // Get the parameters
888 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
890 typename ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
892 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
893 ::GetParameters( void ) const
895 /** NOTE: For efficiency, this class does not keep a copy of the parameters -
896 * it just keeps pointer to input parameters.
898 if (NULL == m_InputParametersPointer)
900 itkExceptionMacro( <<"Cannot GetParameters() because m_InputParametersPointer is NULL. Perhaps SetCoefficientImage() has been called causing the NULL pointer." );
903 return (*m_InputParametersPointer);
907 // Get the parameters
908 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
910 typename ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
912 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
913 ::GetFixedParameters( void ) const
915 RegionType resRegion = this->GetGridRegion( );
917 for (unsigned int i=0; i<NInputDimensions; i++)
919 this->m_FixedParameters[i] = (resRegion.GetSize())[i];
921 for (unsigned int i=0; i<NInputDimensions; i++)
923 this->m_FixedParameters[NInputDimensions+i] = (this->GetGridOrigin())[i];
925 for (unsigned int i=0; i<NInputDimensions; i++)
927 this->m_FixedParameters[2*NInputDimensions+i] = (this->GetGridSpacing())[i];
929 for (unsigned int di=0; di<NInputDimensions; di++)
931 for (unsigned int dj=0; dj<NInputDimensions; dj++)
933 this->m_FixedParameters[3*NInputDimensions+(di*NInputDimensions+dj)] = (this->GetGridDirection())[di][dj];
937 //JV add splineOrders
938 for (unsigned int i=0; i<NInputDimensions; i++)
940 this->m_FixedParameters[(3+NInputDimensions)*NInputDimensions+i] = (this->GetSplineOrders())[i];
943 //JV add LUTsamplingFactor
944 for (unsigned int i=0; i<NInputDimensions; i++)
946 this->m_FixedParameters[(4+NInputDimensions)*NInputDimensions+i] = (this->GetLUTSamplingFactors())[i];
950 this->m_FixedParameters[(5+NInputDimensions)*NInputDimensions]=(double)((size_t) m_Mask);
952 //JV add the bulktransform pointer
953 this->m_FixedParameters[(5+NInputDimensions)*NInputDimensions+1]=(double)((size_t) m_BulkTransform);
955 //JV add the transform shape
956 this->m_FixedParameters[(5+NInputDimensions)*NInputDimensions+2]=(double)( m_TransformShape);
958 return (this->m_FixedParameters);
962 // Set the B-Spline coefficients using input images
963 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
965 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
966 ::SetCoefficientImage( CoefficientImagePointer image )
968 this->SetGridSpacing( image->GetSpacing() );
969 this->SetGridOrigin( image->GetOrigin() );
970 this->SetGridDirection( image->GetDirection() );
971 this->SetGridRegion( image->GetBufferedRegion() );
972 m_CoefficientImage = image;
975 // m_WrappedImage=m_CoefficientImage;
977 // Update the interpolator
978 this->PadCoefficientImage();
979 m_VectorInterpolator->SetInputImage(this->GetPaddedCoefficientImage());
981 // Clean up buffered parameters
982 m_InternalParametersBuffer = ParametersType( 0 );
983 m_InputParametersPointer = NULL;
988 // // Set the B-Spline coefficients using input images
989 // template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
991 // ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
992 // ::SetPaddedCoefficientImage( CoefficientImagePointer image )
994 // //JV modify the region
995 // typename CoefficientImageType::RegionType region=image->GetBufferedRegion();
996 // typename CoefficientImageType::RegionType::SizeType size=region.GetSize();
997 // size[InputDimension-1]-=2;
998 // region.SetSize(size);
1001 // this->SetGridRegion( region );
1002 // this->SetGridSpacing( image->GetSpacing() );
1003 // this->SetGridDirection( image->GetDirection() );
1004 // this->SetGridOrigin( image->GetOrigin() );
1005 // m_PaddedCoefficientImage = image;
1006 // this->ExtractCoefficientImage();
1007 // m_VectorInterpolator->SetInputImage(this->GetPaddedCoefficientImage());
1009 // // Clean up buffered parameters
1010 // m_InternalParametersBuffer = ParametersType( 0 );
1011 // m_InputParametersPointer = NULL;
1015 // Set the B-Spline coefficients using input images
1016 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
1017 typename ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>::CoefficientImageType::Pointer
1018 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
1019 ::ExtractTemporalRow(const typename CoefficientImageType::Pointer& coefficientImage, unsigned int temporalIndex)
1022 typename CoefficientImageType::RegionType sourceRegion=coefficientImage->GetLargestPossibleRegion();
1023 sourceRegion.SetSize(InputDimension-1, 1);
1024 sourceRegion.SetIndex(InputDimension-1, temporalIndex);
1027 typedef clitk::ExtractImageFilter<CoefficientImageType, CoefficientImageType> ExtractImageFilterType;
1028 typename ExtractImageFilterType::Pointer extract=ExtractImageFilterType::New();
1029 extract->SetInput(coefficientImage);
1030 extract->SetExtractionRegion(sourceRegion);
1032 typename CoefficientImageType::Pointer row= extract->GetOutput();
1034 // Set index to zero
1035 sourceRegion.SetIndex(InputDimension-1, 0);
1036 row->SetRegions(sourceRegion);
1041 // Set the B-Spline coefficients using input images
1042 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
1044 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
1045 ::PadCoefficientImage(void)
1048 // Define paste, extract and combine filters
1049 typedef itk::PasteImageFilter<CoefficientImageType, CoefficientImageType, CoefficientImageType> PasteImageFilterType;
1050 typedef clitk::ExtractImageFilter<CoefficientImageType, CoefficientImageType> ExtractImageFilterType;
1051 typedef clitk::LinearCombinationImageFilter<CoefficientImageType, CoefficientImageType> LinearCombinationFilterType;
1054 typename CoefficientImageType::RegionType sourceRegion=m_PaddedCoefficientImage->GetLargestPossibleRegion();
1055 typename CoefficientImageType::RegionType destinationRegion=m_PaddedCoefficientImage->GetLargestPossibleRegion();
1056 typename CoefficientImageType::RegionType::SizeType sourceSize=sourceRegion.GetSize();
1057 typename CoefficientImageType::RegionType::SizeType destinationSize=destinationRegion.GetSize();
1058 typename CoefficientImageType::IndexType sourceIndex=sourceRegion.GetIndex();
1059 typename CoefficientImageType::IndexType destinationIndex=destinationRegion.GetIndex();
1061 // JV Padding depends on the shape
1062 switch (m_TransformShape)
1067 2: rabbit 4 CP 3 DOF
1068 3: rabbit 5 CP 4 DOF
1069 4: sputnik 4 CP 4 DOF
1070 5: sputnik 5 CP 5 DOF
1071 6: diamond 6 CP 5 DOF
1072 7: diamond 7 CP 6 DOF
1077 // ----------------------------------------------------------------------
1078 // The egg with 4 internal CP (starting from inhale)
1079 // Periodic, constrained to zero at the reference
1080 // at position 3 and
1081 // Coeff row BC1 BC2 0 BC3 1 2 BC4
1082 // PaddedCoeff R: 0 1 2 3 4 5 6
1085 // BC3= -weights[2]/weights[0] ( R2+R4 )
1087 // ---------------------------------------------------------------------
1089 //---------------------------------
1090 // 1. First two temporal rows are identical: paste 1-2 to 0-1
1091 typename PasteImageFilterType::Pointer paster0=PasteImageFilterType::New();
1092 paster0->SetDestinationImage(m_PaddedCoefficientImage);
1093 paster0->SetSourceImage(m_CoefficientImage);
1094 sourceRegion.SetIndex(NInputDimensions-1,1);
1095 sourceRegion.SetSize(NInputDimensions-1,2);
1096 destinationIndex[InputDimension-1]=0;
1097 paster0->SetDestinationIndex(destinationIndex);
1098 paster0->SetSourceRegion(sourceRegion);
1100 //---------------------------------
1101 // 1. Next temporal row = paste 0 to 2
1102 typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1103 typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1104 paster1->SetDestinationImage(paster0->GetOutput());
1105 paster1->SetSourceImage(row0);
1106 destinationIndex[InputDimension-1]=2;
1107 paster1->SetDestinationIndex(destinationIndex);
1108 paster1->SetSourceRegion(row0->GetLargestPossibleRegion());
1110 //---------------------------------
1111 // 2. Middle row at index 3=BC
1112 // Extract row at index 1, 2 (of coeff image)
1113 typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1116 typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1117 combine1->SetInput(0,row0);
1118 combine1->SetInput(1,row1);
1119 combine1->SetA(-m_WeightRatio[2][0]);
1120 combine1->SetB(-m_WeightRatio[2][0]);
1122 typename CoefficientImageType::Pointer bc3Row=combine1->GetOutput();
1124 // Paste middleRow at index 3 (padded coeff)
1125 typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1126 paster2->SetDestinationImage(paster1->GetOutput());
1127 paster2->SetSourceImage(bc3Row);
1128 destinationIndex[InputDimension-1]=3;
1129 paster2->SetDestinationIndex(destinationIndex);
1130 paster2->SetSourceRegion(bc3Row->GetLargestPossibleRegion());
1132 //---------------------------------
1133 // 3. Next two temporal rows identical: paste 1,2 to 4,5
1134 typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1135 paster3->SetDestinationImage(paster2->GetOutput());
1136 paster3->SetSourceImage(m_CoefficientImage);
1137 sourceRegion.SetIndex(InputDimension-1,1);
1138 sourceRegion.SetSize(InputDimension-1,2);
1139 destinationIndex[InputDimension-1]=4;
1140 paster3->SetDestinationIndex(destinationIndex);
1141 paster3->SetSourceRegion(sourceRegion);
1143 // ---------------------------------
1144 // 4. Row at index 6=BC4= R2
1145 // Paste BC3 row at index 5
1146 typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
1147 paster4->SetDestinationImage(paster3->GetOutput());
1148 paster4->SetSourceImage(row0);
1149 destinationIndex[InputDimension-1]=6;
1150 paster4->SetDestinationIndex(destinationIndex);
1151 paster4->SetSourceRegion(row0->GetLargestPossibleRegion());
1153 // Update the chain!
1155 m_PaddedCoefficientImage= paster4->GetOutput();
1162 // ----------------------------------------------------------------------
1163 // The egg with 5 internal CP (starting from inhale)
1164 // Periodic, constrained to zero at the reference
1165 // at position 3 and
1166 // Coeff row BC1 BC2 0 BC3 1 2 3 BC4
1167 // PaddedCoeff R: 0 1 2 3 4 5 6 7
1170 // BC3= -weights[3]/weights[1] ( R2+R5 ) - R4
1172 // ---------------------------------------------------------------------
1173 //---------------------------------
1174 // 1. First two temporal rows are identical: paste 2-3 to 0-1
1175 typename PasteImageFilterType::Pointer paster0=PasteImageFilterType::New();
1176 paster0->SetDestinationImage(m_PaddedCoefficientImage);
1177 paster0->SetSourceImage(m_CoefficientImage);
1178 sourceRegion.SetIndex(NInputDimensions-1,2);
1179 sourceRegion.SetSize(NInputDimensions-1,2);
1180 destinationIndex[InputDimension-1]=0;
1181 paster0->SetDestinationIndex(destinationIndex);
1182 paster0->SetSourceRegion(sourceRegion);
1184 //---------------------------------
1185 // 1. Next temporal row = paste 0 to 2
1186 typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1187 typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1188 paster1->SetDestinationImage(paster0->GetOutput());
1189 paster1->SetSourceImage(row0);
1190 destinationIndex[InputDimension-1]=2;
1191 paster1->SetDestinationIndex(destinationIndex);
1192 paster1->SetSourceRegion(row0->GetLargestPossibleRegion());
1194 //---------------------------------
1195 // 2. Middle row at index 3=BC
1196 // Extract row at index 1, 2 (of coeff image)
1197 typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1198 typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
1201 typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1202 combine1->SetInput(0,row0);
1203 combine1->SetInput(1,row2);
1206 // combine1->Update();
1207 typename LinearCombinationFilterType::Pointer combine2=LinearCombinationFilterType::New();
1208 combine2->SetInput(0,row1);
1209 combine2->SetInput(1,combine1->GetOutput());
1210 combine2->SetA(-1.);
1211 combine2->SetB(-m_WeightRatio[3][1]);
1213 typename CoefficientImageType::Pointer bc3Row=combine2->GetOutput();
1215 // Paste middleRow at index 3 (padded coeff)
1216 typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1217 paster2->SetDestinationImage(paster1->GetOutput());
1218 paster2->SetSourceImage(bc3Row);
1219 destinationIndex[InputDimension-1]=3;
1220 paster2->SetDestinationIndex(destinationIndex);
1221 paster2->SetSourceRegion(bc3Row->GetLargestPossibleRegion());
1223 //---------------------------------
1224 // 3. Next three temporal rows identical: paste 1,2,3 to 4,5,6
1225 typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1226 paster3->SetDestinationImage(paster2->GetOutput());
1227 paster3->SetSourceImage(m_CoefficientImage);
1228 sourceRegion.SetIndex(InputDimension-1,1);
1229 sourceRegion.SetSize(InputDimension-1,3);
1230 destinationIndex[InputDimension-1]=4;
1231 paster3->SetDestinationIndex(destinationIndex);
1232 paster3->SetSourceRegion(sourceRegion);
1234 // ---------------------------------
1235 // 4. Row at index 7=BC4= R2
1236 // Paste BC3 row at index 5
1237 typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
1238 paster4->SetDestinationImage(paster3->GetOutput());
1239 paster4->SetSourceImage(row0);
1240 destinationIndex[InputDimension-1]=7;
1241 paster4->SetDestinationIndex(destinationIndex);
1242 paster4->SetSourceRegion(row0->GetLargestPossibleRegion());
1244 // Update the chain!
1246 m_PaddedCoefficientImage= paster4->GetOutput();
1252 // // ----------------------------------------------------------------------
1253 // // The egg with 5 internal CP:
1254 // // Periodic, C2 smooth everywhere and constrained to zero at the reference
1255 // // Coeff row R5 BC1 0 1 2 3 BC2 R2
1256 // // PaddedCoeff R: 0 1 2 3 4 5 6 7
1257 // // BC1= -weights[2]/weights[0] ( R2+R5)
1259 // // ---------------------------------------------------------------------
1261 // // Extract rows with index 0 and 3
1262 // typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1263 // typename CoefficientImageType::Pointer row3=this->ExtractTemporalRow(m_CoefficientImage, 3);
1265 // // Paste the first row
1266 // typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1267 // destinationIndex.Fill(0);
1268 // paster1->SetDestinationIndex(destinationIndex);
1269 // paster1->SetSourceRegion(row3->GetLargestPossibleRegion());
1270 // paster1->SetSourceImage(row3);
1271 // paster1->SetDestinationImage(m_PaddedCoefficientImage);
1273 // // Linearly Combine rows for BC1 and BC2
1274 // typedef clitk::LinearCombinationImageFilter<CoefficientImageType, CoefficientImageType> LinearCombinationFilterType;
1275 // typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1276 // combine1->SetFirstInput(row0);
1277 // combine1->SetSecondInput(row3);
1278 // combine1->SetA(-m_WeightRatio[2][0]);
1279 // combine1->SetB(-m_WeightRatio[2][0]);
1280 // combine1->Update();
1281 // typename CoefficientImageType::Pointer bcRow=combine1->GetOutput();
1283 // // Paste the second row
1284 // typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1285 // destinationIndex[InputDimension-1]=1;
1286 // paster2->SetDestinationIndex(destinationIndex);
1287 // paster2->SetSourceRegion(bcRow->GetLargestPossibleRegion());
1288 // paster2->SetSourceImage(bcRow);
1289 // paster2->SetDestinationImage(paster1->GetOutput());
1291 // // Paste the coefficientImage
1292 // typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1293 // destinationIndex[InputDimension-1]=2;
1294 // paster3->SetDestinationIndex(destinationIndex);
1295 // paster3->SetSourceRegion(m_CoefficientImage->GetLargestPossibleRegion());
1296 // paster3->SetSourceImage(m_CoefficientImage);
1297 // paster3->SetDestinationImage(paster2->GetOutput());
1299 // // Paste the last two rows
1300 // typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
1301 // destinationIndex[InputDimension-1]=6;
1302 // paster4->SetDestinationIndex(destinationIndex);
1303 // paster4->SetSourceRegion(bcRow->GetLargestPossibleRegion());
1304 // paster4->SetSourceImage(bcRow);
1305 // paster4->SetDestinationImage(paster3->GetOutput());
1307 // typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
1308 // destinationIndex[InputDimension-1]=7;
1309 // paster5->SetDestinationIndex(destinationIndex);
1310 // paster5->SetSourceRegion(row0->GetLargestPossibleRegion());
1311 // paster5->SetSourceImage(row0);
1312 // paster5->SetDestinationImage(paster4->GetOutput());
1314 // // Update the chain!
1315 // paster5->Update();
1316 // m_PaddedCoefficientImage= paster5->GetOutput();
1323 // ----------------------------------------------------------------------
1324 // The rabbit with 4 internal CP
1325 // Periodic, constrained to zero at the reference
1326 // at position 3 and the extremes fixed with anit-symm bc
1327 // Coeff row BC1 0 1 BC2 2 BC3 BC4
1328 // PaddedCoeff R: 0 1 2 3 4 5 6
1330 // BC2= -weights[2]/weights[0] ( R2+R4 )
1333 // ---------------------------------------------------------------------
1335 // ---------------------------------
1336 // 0. First Row =BC1
1337 typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1338 typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1339 typename LinearCombinationFilterType::Pointer combine0=LinearCombinationFilterType::New();
1340 combine0->SetInput(0,row0);
1341 combine0->SetInput(1,row1);
1343 combine0->SetB(-1.);
1345 typename CoefficientImageType::Pointer bc1Row=combine0->GetOutput();
1347 typename PasteImageFilterType::Pointer paster0=PasteImageFilterType::New();
1348 paster0->SetDestinationImage(m_PaddedCoefficientImage);
1349 paster0->SetSourceImage(bc1Row);
1350 destinationIndex[InputDimension-1]=0;
1351 paster0->SetDestinationIndex(destinationIndex);
1352 paster0->SetSourceRegion(bc1Row->GetLargestPossibleRegion());
1354 //---------------------------------
1355 // 1. Next two temporal rows are identical: paste 0-1 to 1-2
1356 typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1357 paster1->SetDestinationImage(paster0->GetOutput());
1358 paster1->SetSourceImage(m_CoefficientImage);
1359 sourceRegion.SetIndex(NInputDimensions-1,0);
1360 destinationIndex[InputDimension-1]=1;
1361 sourceRegion.SetSize(NInputDimensions-1,2);
1362 paster1->SetDestinationIndex(destinationIndex);
1363 paster1->SetSourceRegion(sourceRegion);
1365 //---------------------------------
1366 // 2. Middle row at index 3=BC
1367 // Extract row at index 1, 2 (of coeff image)
1368 typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
1371 typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1372 combine1->SetInput(0,row1);
1373 combine1->SetInput(1,row2);
1374 combine1->SetA(-m_WeightRatio[2][0]);
1375 combine1->SetB(-m_WeightRatio[2][0]);
1377 typename CoefficientImageType::Pointer bc2Row=combine1->GetOutput();
1379 // Paste middleRow at index 3 (padded coeff)
1380 typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1381 paster2->SetDestinationImage(paster1->GetOutput());
1382 paster2->SetSourceImage(bc2Row);
1383 destinationIndex[InputDimension-1]=3;
1384 paster2->SetDestinationIndex(destinationIndex);
1385 paster2->SetSourceRegion(bc2Row->GetLargestPossibleRegion());
1387 //---------------------------------
1388 // 3. Next temporal row is identical: paste 2 to 4
1389 typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1390 paster3->SetDestinationImage(paster2->GetOutput());
1391 paster3->SetSourceImage(row2);
1392 destinationIndex[InputDimension-1]=4;
1393 paster3->SetDestinationIndex(destinationIndex);
1394 paster3->SetSourceRegion(row2->GetLargestPossibleRegion());
1396 // ---------------------------------
1397 // 4. Row at index 5=BC (paddedcoeff image) R1
1398 // Paste BC3 row at index 5
1399 typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
1400 paster4->SetDestinationImage(paster3->GetOutput());
1401 paster4->SetSourceImage(row0);
1402 destinationIndex[InputDimension-1]=5;
1403 paster4->SetDestinationIndex(destinationIndex);
1404 paster4->SetSourceRegion(row0->GetLargestPossibleRegion());
1406 // ---------------------------------
1407 // 5. Paste BC4 row at index 6
1408 typename LinearCombinationFilterType::Pointer combine3=LinearCombinationFilterType::New();
1409 combine3->SetInput(0,row0);
1410 combine3->SetInput(1,row2);
1412 combine3->SetB(-1.);
1414 typename CoefficientImageType::Pointer bc4Row=combine3->GetOutput();
1415 typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
1416 paster5->SetDestinationImage(paster4->GetOutput());
1417 paster5->SetSourceImage(bc4Row);
1418 destinationIndex[InputDimension-1]=6;
1419 paster5->SetDestinationIndex(destinationIndex);
1420 paster5->SetSourceRegion(bc4Row->GetLargestPossibleRegion());
1422 // Update the chain!
1424 m_PaddedCoefficientImage= paster5->GetOutput();
1431 // ----------------------------------------------------------------------
1432 // The rabbit with 5 internal CP
1433 // Periodic, constrained to zero at the reference at position 3.5
1434 // and the extremes fixed with anti-symmetrical BC
1435 // Coeff row BC1 0 1 BC2 2 3 BC3 BC4
1436 // PaddedCoeff R: 0 1 2 3 4 5 6 7
1438 // BC2= -weights[3]/weights[1] ( R2+R5 ) - R4
1441 // ---------------------------------------------------------------------
1443 // ---------------------------------
1444 // 0. First Row =BC1
1445 typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1446 typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1447 typename LinearCombinationFilterType::Pointer combine0=LinearCombinationFilterType::New();
1448 combine0->SetInput(0,row0);
1449 combine0->SetInput(1,row1);
1451 combine0->SetB(-1.);
1453 typename CoefficientImageType::Pointer bc1Row=combine0->GetOutput();
1455 typename PasteImageFilterType::Pointer paster0=PasteImageFilterType::New();
1456 paster0->SetDestinationImage(m_PaddedCoefficientImage);
1457 paster0->SetSourceImage(bc1Row);
1458 destinationIndex[InputDimension-1]=0;
1459 paster0->SetDestinationIndex(destinationIndex);
1460 paster0->SetSourceRegion(bc1Row->GetLargestPossibleRegion());
1462 //---------------------------------
1463 // 1. Next two temporal rows are identical: paste 0-1 to 1-2
1464 typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1465 paster1->SetDestinationImage(paster0->GetOutput());
1466 paster1->SetSourceImage(m_CoefficientImage);
1467 sourceRegion.SetIndex(InputDimension-1,0);
1468 destinationIndex[InputDimension-1]=1;
1469 sourceRegion.SetSize(NInputDimensions-1,2);
1470 paster1->SetDestinationIndex(destinationIndex);
1471 paster1->SetSourceRegion(sourceRegion);
1473 //---------------------------------
1474 // 2. Middle row at index 3=BC
1475 // Extract row at index 1, 3 (of coeff image)
1476 //typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1477 typename CoefficientImageType::Pointer row3=this->ExtractTemporalRow(m_CoefficientImage, 3);
1480 typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1481 combine1->SetInput(0,row1);
1482 combine1->SetInput(1,row3);
1487 // Extract row at index 2 (coeff image)
1488 typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
1491 typename LinearCombinationFilterType::Pointer combine2=LinearCombinationFilterType::New();
1492 combine2->SetInput(0,combine1->GetOutput());
1493 combine2->SetInput(1,row2);
1494 combine2->SetA(-m_WeightRatio[3][1]);
1495 combine2->SetB(-1.);
1497 typename CoefficientImageType::Pointer bc2Row=combine2->GetOutput();
1499 // Paste middleRow at index 3 (padded coeff)
1500 typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1501 paster2->SetDestinationImage(paster1->GetOutput());
1502 paster2->SetSourceImage(bc2Row);
1503 destinationIndex[InputDimension-1]=3;
1504 paster2->SetDestinationIndex(destinationIndex);
1505 paster2->SetSourceRegion(bc2Row->GetLargestPossibleRegion());
1507 //---------------------------------
1508 // 3. Next two temporal rows are identical: paste 2,3 to 4,5
1509 typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1510 paster3->SetDestinationImage(paster2->GetOutput());
1511 paster3->SetSourceImage(m_CoefficientImage);
1512 sourceRegion.SetIndex(InputDimension-1,2);
1513 destinationIndex[InputDimension-1]=4;
1514 sourceRegion.SetSize(NInputDimensions-1,2);
1515 paster3->SetDestinationIndex(destinationIndex);
1516 paster3->SetSourceRegion(sourceRegion);
1518 // ---------------------------------
1519 // 4. Row at index 6=BC (paddedcoeff image)R1
1520 // Paste BC3 row at index 6
1521 typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
1522 paster4->SetDestinationImage(paster3->GetOutput());
1523 paster4->SetSourceImage(row0);
1524 destinationIndex[InputDimension-1]=6;
1525 paster4->SetDestinationIndex(destinationIndex);
1526 paster4->SetSourceRegion(row0->GetLargestPossibleRegion());
1528 // ---------------------------------
1529 // 5. Paste BC4 row at index 7
1530 typename LinearCombinationFilterType::Pointer combine3=LinearCombinationFilterType::New();
1531 combine3->SetInput(0,row0);
1532 combine3->SetInput(1,row3);
1534 combine3->SetB(-1.);
1536 typename CoefficientImageType::Pointer bc4Row=combine3->GetOutput();
1537 typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
1538 paster5->SetDestinationImage(paster4->GetOutput());
1539 paster5->SetSourceImage(bc4Row);
1540 destinationIndex[InputDimension-1]=7;
1541 paster5->SetDestinationIndex(destinationIndex);
1542 paster5->SetSourceRegion(bc4Row->GetLargestPossibleRegion());
1544 // Update the chain!
1546 m_PaddedCoefficientImage= paster5->GetOutput();
1553 // ----------------------------------------------------------------------
1554 // The sputnik with 4 internal CP
1555 // Periodic, constrained to zero at the reference
1556 // at position 3 and one indepent extremes copied
1557 // Coeff row BC1 0 1 BC2 2 BC2 3
1558 // PaddedCoeff R: 0 1 2 3 4 5 6
1560 // BC2= -weights[2]/weights[0] ( R2+R4 )
1561 // BC3= weights[2]/weights[0] ( R2-R4) + R1
1562 // ---------------------------------------------------------------------
1564 //---------------------------------
1565 // 1. First Row is equal to last row: paste 3 row to 0
1566 typename CoefficientImageType::Pointer row3=this->ExtractTemporalRow(m_CoefficientImage, 3);
1567 typename PasteImageFilterType::Pointer paster0=PasteImageFilterType::New();
1568 paster0->SetDestinationImage(m_PaddedCoefficientImage);
1569 paster0->SetSourceImage(row3);
1570 destinationIndex[InputDimension-1]=0;
1571 paster0->SetDestinationIndex(destinationIndex);
1572 paster0->SetSourceRegion(row3->GetLargestPossibleRegion());
1574 //---------------------------------
1575 // 1. Next two temporal rows are identical: paste 0-1 to 1-2
1576 typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1577 paster1->SetDestinationImage(paster0->GetOutput());
1578 paster1->SetSourceImage(m_CoefficientImage);
1579 sourceRegion.SetIndex(NInputDimensions-1,0);
1580 destinationIndex[InputDimension-1]=1;
1581 sourceRegion.SetSize(NInputDimensions-1,2);
1582 paster1->SetDestinationIndex(destinationIndex);
1583 paster1->SetSourceRegion(sourceRegion);
1585 //---------------------------------
1586 // 2. Middle row at index 3=BC
1587 // Extract row at index 1, 2 (of coeff image)
1588 typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1589 typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
1592 typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1593 combine1->SetInput(0,row1);
1594 combine1->SetInput(1,row2);
1595 combine1->SetA(-m_WeightRatio[2][0]);
1596 combine1->SetB(-m_WeightRatio[2][0]);
1598 typename CoefficientImageType::Pointer bc1Row=combine1->GetOutput();
1600 // Paste middleRow at index 3 (padded coeff)
1601 typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1602 paster2->SetDestinationImage(paster1->GetOutput());
1603 paster2->SetSourceImage(bc1Row);
1604 destinationIndex[InputDimension-1]=3;
1605 paster2->SetDestinationIndex(destinationIndex);
1606 paster2->SetSourceRegion(bc1Row->GetLargestPossibleRegion());
1608 //---------------------------------
1609 // 3. Next temporal row is identical: paste 2 to 4
1610 typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1611 paster3->SetDestinationImage(paster2->GetOutput());
1612 paster3->SetSourceImage(row2);
1613 destinationIndex[InputDimension-1]=4;
1614 paster3->SetDestinationIndex(destinationIndex);
1615 paster3->SetSourceRegion(row2->GetLargestPossibleRegion());
1617 //---------------------------------
1618 // 4. Final row at index 5=BC3 (paddedcoeff image)R1 R2 R4 corresponds to index in coeff image 0 1 2
1619 // Extract row at index 1, 2 (of coeff image)already done
1620 typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1623 typename LinearCombinationFilterType::Pointer combine3=LinearCombinationFilterType::New();
1624 combine3->SetInput(0,row1);
1625 combine3->SetInput(1,row2);
1627 combine3->SetB(-1.);
1630 typename LinearCombinationFilterType::Pointer combine4=LinearCombinationFilterType::New();
1631 combine4->SetInput(0,combine3->GetOutput());
1632 combine4->SetInput(1,row0);
1633 combine4->SetA(m_WeightRatio[2][0]);
1636 typename CoefficientImageType::Pointer bc2Row=combine4->GetOutput();
1638 // Paste BC row at index 5
1639 typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
1640 paster4->SetDestinationImage(paster3->GetOutput());
1641 paster4->SetSourceImage(bc2Row);
1642 destinationIndex[InputDimension-1]=5;
1643 paster4->SetDestinationIndex(destinationIndex);
1644 paster4->SetSourceRegion(bc2Row->GetLargestPossibleRegion());
1646 // Paste row 3 at index 6
1647 typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
1648 paster5->SetDestinationImage(paster4->GetOutput());
1649 paster5->SetSourceImage(row3);
1650 destinationIndex[InputDimension-1]=6;
1651 paster5->SetDestinationIndex(destinationIndex);
1652 paster5->SetSourceRegion(row3->GetLargestPossibleRegion());
1654 // Update the chain!
1656 m_PaddedCoefficientImage= paster5->GetOutput();
1664 // ----------------------------------------------------------------------
1665 // The sputnik with 5 internal CP
1666 // Periodic, constrained to zero at the reference
1667 // at position 3 and one indepent extreme
1668 // Coeff row BC1 0 1 BC2 2 3 BC3 4
1669 // PaddedCoeff R: 0 1 2 3 4 5 6 7
1670 // BC1= R2 + R5 - R7
1671 // BC2= -weights[3]/weights[1] ( R2+R5 ) - R4
1672 // BC3= R1 + 0.5 R2 - 0.5 R7
1673 // ----------------------------------------------------------------------
1674 //---------------------------------
1676 typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1677 typename CoefficientImageType::Pointer row3=this->ExtractTemporalRow(m_CoefficientImage, 3);
1678 typename CoefficientImageType::Pointer row4=this->ExtractTemporalRow(m_CoefficientImage, 4);
1681 typename LinearCombinationFilterType::Pointer combine0=LinearCombinationFilterType::New();
1682 combine0->SetInput(0,row1);
1683 combine0->SetInput(1,row3);
1686 //combine0->Update();
1687 typename LinearCombinationFilterType::Pointer combine0bis=LinearCombinationFilterType::New();
1688 combine0bis->SetInput(0,combine0->GetOutput());
1689 combine0bis->SetInput(1,row4);
1690 combine0bis->SetA(1.);
1691 combine0bis->SetB(-1.);
1692 combine0bis->Update();
1693 typename CoefficientImageType::Pointer bc1Row=combine0bis->GetOutput();
1695 typename PasteImageFilterType::Pointer paster0=PasteImageFilterType::New();
1696 paster0->SetDestinationImage(m_PaddedCoefficientImage);
1697 paster0->SetSourceImage(bc1Row);
1698 destinationIndex[InputDimension-1]=0;
1699 paster0->SetDestinationIndex(destinationIndex);
1700 paster0->SetSourceRegion(bc1Row->GetLargestPossibleRegion());
1702 //---------------------------------
1703 // 1. Next two temporal rows are identical: paste 0-1 to 1-2
1704 typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1705 paster1->SetDestinationImage(paster0->GetOutput());
1706 paster1->SetSourceImage(m_CoefficientImage);
1707 sourceRegion.SetIndex(NInputDimensions-1,0);
1708 destinationIndex[InputDimension-1]=1;
1709 sourceRegion.SetSize(NInputDimensions-1,2);
1710 paster1->SetDestinationIndex(destinationIndex);
1711 paster1->SetSourceRegion(sourceRegion);
1713 //---------------------------------
1714 // 2. Middle row at index 3=BC
1715 // Extract row at index 1, 2 (of coeff image)
1716 //typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1717 typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
1720 typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1721 combine1->SetInput(0,row1);
1722 combine1->SetInput(1,row3);
1727 typename LinearCombinationFilterType::Pointer combine2=LinearCombinationFilterType::New();
1728 combine2->SetInput(0,combine1->GetOutput());
1729 combine2->SetInput(1,row2);
1730 combine2->SetA(-m_WeightRatio[3][1]);
1731 combine2->SetB(-1.);
1733 typename CoefficientImageType::Pointer bc2Row=combine2->GetOutput();
1735 // Paste middleRow at index 3 (padded coeff)
1736 typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1737 paster2->SetDestinationImage(paster1->GetOutput());
1738 paster2->SetSourceImage(bc2Row);
1739 destinationIndex[InputDimension-1]=3;
1740 paster2->SetDestinationIndex(destinationIndex);
1741 paster2->SetSourceRegion(bc2Row->GetLargestPossibleRegion());
1743 //---------------------------------
1744 // 3. Next two temporal rows are identical: paste 2,3 to 4,5
1745 typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1746 paster3->SetDestinationImage(paster2->GetOutput());
1747 paster3->SetSourceImage(m_CoefficientImage);
1748 sourceRegion.SetIndex(InputDimension-1,2);
1749 destinationIndex[InputDimension-1]=4;
1750 sourceRegion.SetSize(NInputDimensions-1,2);
1751 paster3->SetDestinationIndex(destinationIndex);
1752 paster3->SetSourceRegion(sourceRegion);
1754 //---------------------------------
1755 // 4. Final row at index 6=BC3
1756 typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1759 typename LinearCombinationFilterType::Pointer combine3=LinearCombinationFilterType::New();
1760 combine3->SetInput(0,row1);
1761 combine3->SetInput(1,row4);
1763 combine3->SetB(-1.);
1764 typename LinearCombinationFilterType::Pointer combine4=LinearCombinationFilterType::New();
1765 combine4->SetInput(0,row0);
1766 combine4->SetInput(1,combine3->GetOutput());
1770 typename CoefficientImageType::Pointer bc3Row=combine4->GetOutput();
1772 // Paste BC row at index 6
1773 typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
1774 paster4->SetDestinationImage(paster3->GetOutput());
1775 paster4->SetSourceImage(bc3Row);
1776 destinationIndex[InputDimension-1]=6;
1777 paster4->SetDestinationIndex(destinationIndex);
1778 paster4->SetSourceRegion(bc3Row->GetLargestPossibleRegion());
1780 // Paste row 4 at index 7
1781 typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
1782 paster5->SetDestinationImage(paster4->GetOutput());
1783 paster5->SetSourceImage(row4);
1784 destinationIndex[InputDimension-1]=7;
1785 paster5->SetDestinationIndex(destinationIndex);
1786 paster5->SetSourceRegion(row4->GetLargestPossibleRegion());
1788 // Update the chain!
1790 m_PaddedCoefficientImage= paster5->GetOutput();
1797 // ----------------------------------------------------------------------
1798 // The diamond with 4 internal CP:
1799 // Periodic, constrained to zero at the reference at position 3
1800 // Coeff row 0 1 2 BC1 3 BC2 4
1801 // PaddedCoeff R:0 1 2 3 4 5 6
1802 // BC1= -weights[2]/weights[0] ( R2+R4 )
1803 // BC2= weights[2]/weights[0] ( R0+R2-R4-R6 ) + R1
1804 // ---------------------------------------------------------------------
1806 //---------------------------------
1807 // 1. First Three temporal rows are identical: paste 0-2 to 0-2
1808 typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1809 paster1->SetDestinationImage(m_PaddedCoefficientImage);
1810 paster1->SetSourceImage(m_CoefficientImage);
1811 sourceIndex.Fill(0);
1812 destinationIndex.Fill(0);
1813 sourceSize[NInputDimensions-1]=3;
1814 sourceRegion.SetSize(sourceSize);
1815 sourceRegion.SetIndex(sourceIndex);
1816 paster1->SetDestinationIndex(destinationIndex);
1817 paster1->SetSourceRegion(sourceRegion);
1819 //---------------------------------
1820 // 2. Middle row at index 3=BC
1821 // Extract row at index 0, 4 (of coeff image)
1822 typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
1823 typename CoefficientImageType::Pointer row3=this->ExtractTemporalRow(m_CoefficientImage, 3);
1826 typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1827 combine1->SetInput(0,row2);
1828 combine1->SetInput(1,row3);
1829 combine1->SetA(-m_WeightRatio[2][0]);
1830 combine1->SetB(-m_WeightRatio[2][0]);
1832 typename CoefficientImageType::Pointer bc1Row=combine1->GetOutput();
1834 // Paste middleRow at index 3 (padded coeff)
1835 typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1836 paster2->SetDestinationImage(paster1->GetOutput());
1837 paster2->SetSourceImage(bc1Row);
1838 destinationIndex[InputDimension-1]=3;
1839 paster2->SetDestinationIndex(destinationIndex);
1840 paster2->SetSourceRegion(bc1Row->GetLargestPossibleRegion());
1842 //---------------------------------
1843 // 3. Next row identical: paste 3 to 4
1844 typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1845 paster3->SetDestinationImage(paster2->GetOutput());
1846 paster3->SetSourceImage(row3);
1847 destinationIndex[InputDimension-1]=4;
1848 paster3->SetDestinationIndex(destinationIndex);
1849 paster3->SetSourceRegion(row3->GetLargestPossibleRegion());
1851 //---------------------------------
1852 // 4. Final row at index 6=BC (paddedcoeff image)R0 R2 R5 R7 R1 corresponds to index in coeff image 0 2 4 5 1
1853 // Extract row at index 0, 2 (of coeff image)already done
1854 typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1855 typename CoefficientImageType::Pointer row4=this->ExtractTemporalRow(m_CoefficientImage, 4);
1856 typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1859 typename LinearCombinationFilterType::Pointer combine3=LinearCombinationFilterType::New();
1860 combine3->SetInput(0,row0);
1861 combine3->SetInput(1,row2);
1866 typename LinearCombinationFilterType::Pointer combine4=LinearCombinationFilterType::New();
1867 combine4->SetInput(0,row3);
1868 combine4->SetInput(1,row4);
1873 typename LinearCombinationFilterType::Pointer combine5=LinearCombinationFilterType::New();
1874 combine5->SetInput(0,combine3->GetOutput());
1875 combine5->SetInput(1,combine4->GetOutput());
1877 combine5->SetB(-1.);
1880 typename LinearCombinationFilterType::Pointer combine6=LinearCombinationFilterType::New();
1881 combine6->SetInput(0,combine5->GetOutput());
1882 combine6->SetInput(1,row1);
1883 combine6->SetA(m_WeightRatio[2][0]);
1886 typename CoefficientImageType::Pointer bc2Row=combine6->GetOutput();
1888 // Paste BC row at index 5
1889 typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
1890 paster4->SetDestinationImage(paster3->GetOutput());
1891 paster4->SetSourceImage(bc2Row);
1892 destinationIndex[InputDimension-1]=5;
1893 paster4->SetDestinationIndex(destinationIndex);
1894 paster4->SetSourceRegion(bc2Row->GetLargestPossibleRegion());
1896 // Paste last row at index 6
1897 typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
1898 paster5->SetDestinationImage(paster4->GetOutput());
1899 paster5->SetSourceImage(row4);
1900 destinationIndex[InputDimension-1]=6;
1901 paster5->SetDestinationIndex(destinationIndex);
1902 paster5->SetSourceRegion(row4->GetLargestPossibleRegion());
1904 // Update the chain!
1906 m_PaddedCoefficientImage= paster5->GetOutput();
1912 // ----------------------------------------------------------------------
1913 // The diamond with 5 internal CP:
1914 // periodic, constrained to zero at the reference at position 3.5
1915 // Coeff row 0 1 2 BC1 3 4 BC2 5
1916 // PaddedCoeff R:0 1 2 3 4 5 6 7
1917 // BC1= -weights[3]/weights[1] ( R2+R5 ) - R4
1918 // BC2= weights[2]/weights[0] ( R0+R2-R5-R7 ) + R1
1919 // ---------------------------------------------------------------------
1921 //---------------------------------
1922 // 1. First Three temporal rows are identical: paste 0-2 to 0-2
1923 typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1924 paster1->SetDestinationImage(m_PaddedCoefficientImage);
1925 paster1->SetSourceImage(m_CoefficientImage);
1926 sourceIndex.Fill(0);
1927 destinationIndex.Fill(0);
1928 sourceSize[NInputDimensions-1]=3;
1929 sourceRegion.SetSize(sourceSize);
1930 sourceRegion.SetIndex(sourceIndex);
1931 paster1->SetDestinationIndex(destinationIndex);
1932 paster1->SetSourceRegion(sourceRegion);
1934 //---------------------------------
1935 // 2. Middle row at index 3=BC
1936 // Extract row at index 0, 4 (of coeff image)
1937 typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
1938 typename CoefficientImageType::Pointer row4=this->ExtractTemporalRow(m_CoefficientImage, 4);
1941 typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1942 combine1->SetInput(0,row2);
1943 combine1->SetInput(1,row4);
1948 // Extract row at index 3 (coeff image)
1949 typename CoefficientImageType::Pointer row3=this->ExtractTemporalRow(m_CoefficientImage, 3);
1952 typename LinearCombinationFilterType::Pointer combine2=LinearCombinationFilterType::New();
1953 combine2->SetInput(0,combine1->GetOutput());
1954 combine2->SetInput(1,row3);
1955 combine2->SetA(-m_WeightRatio[3][1] );
1956 combine2->SetB(-1.);
1958 typename CoefficientImageType::Pointer bc1Row=combine2->GetOutput();
1960 // Paste middleRow at index 3 (padded coeff)
1961 typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1962 paster2->SetDestinationImage(paster1->GetOutput());
1963 paster2->SetSourceImage(bc1Row);
1964 destinationIndex[InputDimension-1]=3;
1965 paster2->SetDestinationIndex(destinationIndex);
1966 paster2->SetSourceRegion(bc1Row->GetLargestPossibleRegion());
1968 //---------------------------------
1969 // 3. Next two temporal rows are identical: paste 3,4 to 4,5
1970 typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1971 paster3->SetDestinationImage(paster2->GetOutput());
1972 paster3->SetSourceImage(m_CoefficientImage);
1973 sourceIndex[InputDimension-1]=3;
1974 destinationIndex[InputDimension-1]=4;
1975 sourceSize[NInputDimensions-1]=2;
1976 sourceRegion.SetSize(sourceSize);
1977 sourceRegion.SetIndex(sourceIndex);
1978 paster3->SetDestinationIndex(destinationIndex);
1979 paster3->SetSourceRegion(sourceRegion);
1981 //---------------------------------
1982 // 4. Final row at index 6=BC (paddedcoeff image)R0 R2 R5 R7 R1 corresponds to index in coeff image 0 2 4 5 1
1983 // Extract row at index 0, 2 (of coeff image)already done
1984 typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1985 typename CoefficientImageType::Pointer row5=this->ExtractTemporalRow(m_CoefficientImage, 5);
1986 typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1989 typename LinearCombinationFilterType::Pointer combine3=LinearCombinationFilterType::New();
1990 combine3->SetInput(0,row0);
1991 combine3->SetInput(1,row2);
1996 typename LinearCombinationFilterType::Pointer combine4=LinearCombinationFilterType::New();
1997 combine4->SetInput(0,row4);
1998 combine4->SetInput(1,row5);
2003 typename LinearCombinationFilterType::Pointer combine5=LinearCombinationFilterType::New();
2004 combine5->SetInput(0,combine3->GetOutput());
2005 combine5->SetInput(1,combine4->GetOutput());
2007 combine5->SetB(-1.);
2010 typename LinearCombinationFilterType::Pointer combine6=LinearCombinationFilterType::New();
2011 combine6->SetInput(0,combine5->GetOutput());
2012 combine6->SetInput(1,row1);
2013 combine6->SetA(m_WeightRatio[2][0]);
2016 typename CoefficientImageType::Pointer bc2Row=combine6->GetOutput();
2018 // Paste BC row at index 6
2019 typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
2020 paster4->SetDestinationImage(paster3->GetOutput());
2021 paster4->SetSourceImage(bc2Row);
2022 destinationIndex[InputDimension-1]=6;
2023 paster4->SetDestinationIndex(destinationIndex);
2024 paster4->SetSourceRegion(bc2Row->GetLargestPossibleRegion());
2026 // Paste last row at index 7
2027 typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
2028 paster5->SetDestinationImage(paster4->GetOutput());
2029 paster5->SetSourceImage(row5);
2030 destinationIndex[InputDimension-1]=7;
2031 paster5->SetDestinationIndex(destinationIndex);
2032 paster5->SetSourceRegion(row5->GetLargestPossibleRegion());
2034 // Update the chain!
2036 m_PaddedCoefficientImage= paster5->GetOutput();
2043 // ----------------------------------------------------------------------
2044 // The sputnik with 5 internal CP T''(0)=T''(10)
2045 // Periodic, constrained to zero at the reference
2046 // at position 3 and one indepent extreme
2047 // Coeff row BC1 0 1 BC2 2 3 BC3 4
2048 // PaddedCoeff R: 0 1 2 3 4 5 6 7
2050 // BC2= -weights[3]/weights[1] ( R2+R5 ) - R4
2052 // ---------------------------------------------------------------------
2054 //---------------------------------
2056 typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
2057 typename CoefficientImageType::Pointer row3=this->ExtractTemporalRow(m_CoefficientImage, 3);
2058 typename CoefficientImageType::Pointer row4=this->ExtractTemporalRow(m_CoefficientImage, 4);
2061 typename LinearCombinationFilterType::Pointer combine0=LinearCombinationFilterType::New();
2062 combine0->SetInput(0,row3);
2063 combine0->SetInput(1,row4);
2066 typename LinearCombinationFilterType::Pointer combine0bis=LinearCombinationFilterType::New();
2067 combine0bis->SetInput(0,combine0->GetOutput());
2068 combine0bis->SetInput(1,row1);
2069 combine0bis->SetA(1.);
2070 combine0bis->SetB(-1.);
2071 combine0bis->Update();
2072 typename CoefficientImageType::Pointer bc1Row=combine0bis->GetOutput();
2074 typename PasteImageFilterType::Pointer paster0=PasteImageFilterType::New();
2075 paster0->SetDestinationImage(m_PaddedCoefficientImage);
2076 paster0->SetSourceImage(bc1Row);
2077 destinationIndex[InputDimension-1]=0;
2078 paster0->SetDestinationIndex(destinationIndex);
2079 paster0->SetSourceRegion(bc1Row->GetLargestPossibleRegion());
2081 //---------------------------------
2082 // 1. Next two temporal rows are identical: paste 0-1 to 1-2
2083 typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
2084 paster1->SetDestinationImage(paster0->GetOutput());
2085 paster1->SetSourceImage(m_CoefficientImage);
2086 sourceRegion.SetIndex(NInputDimensions-1,0);
2087 destinationIndex[InputDimension-1]=1;
2088 sourceRegion.SetSize(NInputDimensions-1,2);
2089 paster1->SetDestinationIndex(destinationIndex);
2090 paster1->SetSourceRegion(sourceRegion);
2092 //---------------------------------
2093 // 2. Middle row at index 3=BC
2094 // Extract row at index 1, 2 (of coeff image)
2095 // typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
2096 typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
2099 typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
2100 combine1->SetInput(0,row1);
2101 combine1->SetInput(1,row3);
2106 typename LinearCombinationFilterType::Pointer combine2=LinearCombinationFilterType::New();
2107 combine2->SetInput(0,combine1->GetOutput());
2108 combine2->SetInput(1,row2);
2109 combine2->SetA(-m_WeightRatio[3][1]);
2110 combine2->SetB(-1.);
2112 typename CoefficientImageType::Pointer bc2Row=combine2->GetOutput();
2114 // Paste middleRow at index 3 (padded coeff)
2115 typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
2116 paster2->SetDestinationImage(paster1->GetOutput());
2117 paster2->SetSourceImage(bc2Row);
2118 destinationIndex[InputDimension-1]=3;
2119 paster2->SetDestinationIndex(destinationIndex);
2120 paster2->SetSourceRegion(bc2Row->GetLargestPossibleRegion());
2122 //---------------------------------
2123 // 3. Next two temporal rows are identical: paste 2,3 to 4,5
2124 typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
2125 paster3->SetDestinationImage(paster2->GetOutput());
2126 paster3->SetSourceImage(m_CoefficientImage);
2127 sourceRegion.SetIndex(InputDimension-1,2);
2128 destinationIndex[InputDimension-1]=4;
2129 sourceRegion.SetSize(NInputDimensions-1,2);
2130 paster3->SetDestinationIndex(destinationIndex);
2131 paster3->SetSourceRegion(sourceRegion);
2133 //---------------------------------
2134 // 4. Final row at index 6=BC3
2135 typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
2138 // typename LinearCombinationFilterType::Pointer combine3=LinearCombinationFilterType::New();
2139 // combine3->SetInput(0,row0);
2140 // combine3->SetInput(1,row1);
2141 // combine3->SetA(1.);
2142 // combine3->SetB(0.5);
2143 // combine3->Update();
2144 // typename CoefficientImageType::Pointer bc3Row=combine3->GetOutput();
2146 // Paste BC row at index 6
2147 typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
2148 paster4->SetDestinationImage(paster3->GetOutput());
2149 paster4->SetSourceImage(row0);
2150 destinationIndex[InputDimension-1]=6;
2151 paster4->SetDestinationIndex(destinationIndex);
2152 paster4->SetSourceRegion(row0->GetLargestPossibleRegion());
2154 // Paste row 4 at index 7
2155 typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
2156 paster5->SetDestinationImage(paster4->GetOutput());
2157 paster5->SetSourceImage(row4);
2158 destinationIndex[InputDimension-1]=7;
2159 paster5->SetDestinationIndex(destinationIndex);
2160 paster5->SetSourceRegion(row4->GetLargestPossibleRegion());
2162 // Update the chain!
2164 m_PaddedCoefficientImage= paster5->GetOutput();
2170 DD ("Shape not available");
2176 // // Extract coefficients from padded version
2177 // template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
2179 // ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2180 // ::ExtractCoefficientImage( )
2182 // ////DD("extract coeff image");
2183 // typename ExtractImageFilterType::Pointer extractImageFilter = ExtractImageFilterType::New();
2184 // typename CoefficientImageType::RegionType extractionRegion=m_PaddedCoefficientImage->GetLargestPossibleRegion();
2185 // typename CoefficientImageType::RegionType::SizeType extractionSize=extractionRegion.GetSize();
2186 // typename CoefficientImageType::RegionType::IndexType extractionIndex = extractionRegion.GetIndex();
2187 // extractionSize[InputDimension-1]-=4;
2188 // extractionIndex[InputDimension-1]=2;
2189 // extractionRegion.SetSize(extractionSize);
2190 // extractionRegion.SetIndex(extractionIndex);
2191 // extractImageFilter->SetInput(m_PaddedCoefficientImage);
2192 // extractImageFilter->SetExtractionRegion(extractionRegion);
2193 // extractImageFilter->Update();
2194 // m_CoefficientImage=extractImageFilter->GetOutput();
2199 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
2201 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2202 ::PrintSelf(std::ostream &os, itk::Indent indent) const
2205 this->Superclass::PrintSelf(os, indent);
2207 os << indent << "GridRegion: " << m_GridRegion << std::endl;
2208 os << indent << "GridOrigin: " << m_GridOrigin << std::endl;
2209 os << indent << "GridSpacing: " << m_GridSpacing << std::endl;
2210 os << indent << "GridDirection: " << m_GridDirection << std::endl;
2211 os << indent << "IndexToPoint: " << m_IndexToPoint << std::endl;
2212 os << indent << "PointToIndex: " << m_PointToIndex << std::endl;
2214 os << indent << "CoefficientImage: [ ";
2215 os << m_CoefficientImage.GetPointer() << " ]" << std::endl;
2217 os << indent << "WrappedImage: [ ";
2218 os << m_WrappedImage.GetPointer() << " ]" << std::endl;
2220 os << indent << "InputParametersPointer: "
2221 << m_InputParametersPointer << std::endl;
2222 os << indent << "ValidRegion: " << m_ValidRegion << std::endl;
2223 os << indent << "LastJacobianIndex: " << m_LastJacobianIndex << std::endl;
2224 os << indent << "BulkTransform: ";
2225 os << m_BulkTransform << std::endl;
2227 if ( m_BulkTransform )
2229 os << indent << "BulkTransformType: "
2230 << m_BulkTransform->GetNameOfClass() << std::endl;
2232 os << indent << "VectorBSplineInterpolator: ";
2233 os << m_VectorInterpolator.GetPointer() << std::endl;
2234 os << indent << "Mask: ";
2235 os << m_Mask<< std::endl;
2240 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
2242 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2243 ::InsideValidRegion( const ContinuousIndexType& index ) const
2247 if ( !m_ValidRegion.IsInside( index ) )
2252 // JV verify all dimensions
2255 typedef typename ContinuousIndexType::ValueType ValueType;
2256 for( unsigned int j = 0; j < NInputDimensions; j++ )
2258 if (m_SplineOrderOdd[j])
2260 if ( index[j] >= static_cast<ValueType>( m_ValidRegionLast[j] ) )
2272 // Transform a point
2273 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
2274 typename ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2276 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2277 ::TransformPoint(const InputPointType &inputPoint) const
2280 InputPointType transformedPoint;
2281 OutputPointType outputPoint;
2284 if ( m_BulkTransform )
2286 transformedPoint = m_BulkTransform->TransformPoint( inputPoint );
2290 transformedPoint = inputPoint;
2293 // Deformable transform
2294 if ( m_PaddedCoefficientImage )
2296 // Check if inside mask
2297 if(m_Mask && !(m_Mask->IsInside(inputPoint) ) )
2299 // Outside: no (deformable) displacement
2300 return transformedPoint;
2303 // Check if inside valid region
2305 ContinuousIndexType index;
2306 this->TransformPointToContinuousIndex( inputPoint, index );
2307 inside = this->InsideValidRegion( index );
2310 // Outside: no (deformable) displacement
2311 DD("outside valid region");
2312 return transformedPoint;
2315 // Call the vector interpolator
2316 itk::Vector<TCoordRep,SpaceDimension> displacement=m_VectorInterpolator->EvaluateAtContinuousIndex(index);
2318 // JV add for the spatial dimensions
2319 outputPoint=transformedPoint;
2320 for (unsigned int i=0; i<NInputDimensions-1; i++)
2321 outputPoint[i] += displacement[i];
2327 itkWarningMacro( << "B-spline coefficients have not been set" );
2328 outputPoint = transformedPoint;
2336 //JV Deformably transform a point
2337 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
2338 typename ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2340 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2341 ::DeformablyTransformPoint(const InputPointType &inputPoint) const
2343 OutputPointType outputPoint;
2344 if ( m_PaddedCoefficientImage )
2347 // Check if inside mask
2348 if(m_Mask && !(m_Mask->IsInside(inputPoint) ) )
2350 // Outside: no (deformable) displacement
2356 ContinuousIndexType index;
2357 this->TransformPointToContinuousIndex( inputPoint, index );
2358 inside = this->InsideValidRegion( index );
2362 //outside: no (deformable) displacement
2363 outputPoint = inputPoint;
2367 // Call the vector interpolator
2368 itk::Vector<TCoordRep,SpaceDimension> displacement=m_VectorInterpolator->EvaluateAtContinuousIndex(index);
2370 // JV add for the spatial dimensions
2371 outputPoint=inputPoint;
2372 for (unsigned int i=0; i<NInputDimensions-1; i++)
2373 outputPoint[i] += displacement[i];
2376 // No coefficients available
2379 itkWarningMacro( << "B-spline coefficients have not been set" );
2380 outputPoint = inputPoint;
2387 // JV weights are identical as for transformpoint, could be done simultaneously in metric!!!!
2388 // Compute the Jacobian in one position
2389 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
2390 #if ITK_VERSION_MAJOR >= 4
2392 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2393 ::ComputeJacobianWithRespectToParameters( const InputPointType & point, JacobianType & jacobian) const
2396 typename ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2398 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2399 ::GetJacobian( const InputPointType & point ) const
2403 //========================================================
2404 // Zero all components of jacobian
2405 //========================================================
2406 // JV not thread safe (m_LastJacobianIndex), instantiate N transforms
2407 // NOTE: for efficiency, we only need to zero out the coefficients
2408 // that got fill last time this method was called.
2409 unsigned int j=0,b=0;
2411 //Set the previously-set to zero
2412 for ( j = 0; j < SpaceDimension; j++ )
2414 m_FirstIterator[j].GoToBegin();
2415 while ( !m_FirstIterator[j].IsAtEnd() )
2417 m_FirstIterator[j].Set( m_ZeroVector );
2418 ++(m_FirstIterator[j]);
2422 //Set the previously-set to zero
2423 for ( j = 0; j < SpaceDimension; j++ )
2425 m_SecondIterator[j].GoToBegin();
2426 while ( ! (m_SecondIterator[j]).IsAtEnd() )
2428 m_SecondIterator[j].Set( m_ZeroVector );
2429 ++(m_SecondIterator[j]);
2433 //Set the previously-set to zero
2435 for ( j = 0; j < SpaceDimension; j++ )
2437 m_ThirdIterator[j].GoToBegin();
2438 while ( ! (m_ThirdIterator[j]).IsAtEnd() )
2440 m_ThirdIterator[j].Set( m_ZeroVector );
2441 ++(m_ThirdIterator[j]);
2445 //Set the previously-set to zero
2447 for (b=0; b<m_BCSize;b++)
2448 if ( !( m_FirstRegion.IsInside(m_BCRegions[b]) | m_SecondRegion.IsInside(m_BCRegions[b]) ) )
2449 for ( j = 0; j < SpaceDimension; j++ )
2451 m_BCIterators[j][b].GoToBegin();
2452 while ( ! (m_BCIterators[j][b]).IsAtEnd() )
2454 m_BCIterators[j][b].Set( m_ZeroVector );
2455 ++(m_BCIterators[j][b]);
2459 //Set the previously-set to zero
2461 for (b=0; b<m_BC2Size;b++)
2462 if ( !( m_FirstRegion.IsInside(m_BC2Regions[b]) | m_SecondRegion.IsInside(m_BC2Regions[b]) ) )
2463 for ( j = 0; j < SpaceDimension; j++ )
2465 m_BC2Iterators[j][b].GoToBegin();
2466 while ( ! (m_BC2Iterators[j][b]).IsAtEnd() )
2468 m_BC2Iterators[j][b].Set( m_ZeroVector );
2469 ++(m_BC2Iterators[j][b]);
2473 //Set the previously-set to zero
2475 for (b=0; b<m_BC3Size;b++)
2476 if ( !( m_FirstRegion.IsInside(m_BC3Regions[b]) | m_SecondRegion.IsInside(m_BC3Regions[b]) ) )
2477 for ( j = 0; j < SpaceDimension; j++ )
2479 m_BC3Iterators[j][b].GoToBegin();
2480 while ( ! (m_BC3Iterators[j][b]).IsAtEnd() )
2482 m_BC3Iterators[j][b].Set( m_ZeroVector );
2483 ++(m_BC3Iterators[j][b]);
2488 //========================================================
2489 // For each dimension, copy the weight to the support region
2490 //========================================================
2492 // Check if inside mask
2493 if(m_Mask && !(m_Mask->IsInside(point) ) )
2495 // Outside: no (deformable) displacement
2496 #if ITK_VERSION_MAJOR >= 4
2497 jacobian = m_SharedDataBSplineJacobian;
2500 return this->m_Jacobian;
2505 this->TransformPointToContinuousIndex( point, m_Index );
2507 // NOTE: if the support region does not lie totally within the grid
2508 // we assume zero displacement and return the input point
2509 if ( !this->InsideValidRegion( m_Index ) )
2511 #if ITK_VERSION_MAJOR >= 4
2512 jacobian = m_SharedDataBSplineJacobian;
2515 return this->m_Jacobian;
2519 // Compute interpolation weights
2520 const WeightsDataType *weights=NULL;
2521 m_VectorInterpolator->EvaluateWeightsAtContinuousIndex( m_Index, &weights, m_LastJacobianIndex);
2524 m_SupportRegion.SetIndex( m_LastJacobianIndex );
2525 WrapRegion(m_SupportRegion, m_FirstRegion, m_SecondRegion, m_ThirdRegion, m_BCRegions, m_BCValues, m_BC2Regions, m_BC2Values, m_BC3Regions, m_BC3Values, m_InitialOffset);
2526 m_ThirdSize=m_ThirdRegion.GetSize()[InputDimension-1];
2527 m_BCSize=m_BCRegions.size();
2528 m_BC2Size=m_BC2Regions.size();
2529 m_BC3Size=m_BC3Regions.size();
2531 // Reset the iterators
2532 for ( j = 0; j < SpaceDimension ; j++ )
2534 m_FirstIterator[j] = IteratorType( m_JacobianImage[j], m_FirstRegion);
2535 m_SecondIterator[j] = IteratorType( m_JacobianImage[j], m_SecondRegion);
2536 if(m_ThirdSize) m_ThirdIterator[j] = IteratorType( m_JacobianImage[j], m_ThirdRegion);
2538 m_BCIterators[j].resize(m_BCSize);
2539 for (b=0; b<m_BCSize;b++)
2540 m_BCIterators[j][b]= IteratorType( m_JacobianImage[j], m_BCRegions[b]);
2541 m_BC2Iterators[j].resize(m_BC2Size);
2542 for (b=0; b<m_BC2Size;b++)
2543 m_BC2Iterators[j][b]= IteratorType( m_JacobianImage[j], m_BC2Regions[b]);
2544 m_BC3Iterators[j].resize(m_BC3Size);
2545 for (b=0; b<m_BC3Size;b++)
2546 m_BC3Iterators[j][b]= IteratorType( m_JacobianImage[j], m_BC3Regions[b]);
2549 // Skip if on a fixed condition
2552 if (m_BCSize) weights+=m_InitialOffset*m_BCRegions[0].GetNumberOfPixels();
2553 else std::cerr<<"InitialOffset without BCSize: Error!!!!!!!"<<std::endl;
2556 //copy weight to jacobian image
2557 for ( j = 0; j < SpaceDimension; j++ )
2559 // For each dimension, copy the weight to the support region
2560 while ( ! (m_FirstIterator[j]).IsAtEnd() )
2562 m_ZeroVector[j]=*weights;
2563 (m_FirstIterator[j]).Set( m_ZeroVector);
2564 ++(m_FirstIterator[j]);
2568 m_ZeroVector[j]=itk::NumericTraits<JacobianValueType>::Zero;
2569 if (j != SpaceDimension-1) weights-=m_FirstRegion.GetNumberOfPixels();
2572 // Skip BC1 and go to the second region
2573 if (m_BCSize) weights+=m_BCRegions[0].GetNumberOfPixels();
2575 // For each dimension, copy the weight to the support region
2576 //copy weight to jacobian image
2577 for ( j = 0; j < SpaceDimension; j++ )
2579 while ( ! (m_SecondIterator[j]).IsAtEnd() )
2581 m_ZeroVector[j]=*weights;
2582 (m_SecondIterator[j]).Set( m_ZeroVector);
2583 ++(m_SecondIterator[j]);
2586 weights-=m_SecondRegion.GetNumberOfPixels();
2587 m_ZeroVector[j]=itk::NumericTraits<JacobianValueType>::Zero;
2593 // Put pointer in correct position
2594 weights-=m_BCRegions[0].GetNumberOfPixels();
2596 for ( j = 0; j < SpaceDimension; j++ )
2598 for ( b=0; b < m_BCSize; b++ )
2600 while ( ! (m_BCIterators[j][b]).IsAtEnd() )
2602 //copy weight to jacobian image
2603 m_ZeroVector[j]=(*weights) * m_BCValues[b];
2604 (m_BCIterators[j][b]).Value()+= m_ZeroVector;
2605 ++(m_BCIterators[j][b]);
2608 weights-=m_BCRegions[b].GetNumberOfPixels();
2610 m_ZeroVector[j]=itk::NumericTraits<JacobianValueType>::Zero;
2612 weights+=m_BCRegions[m_BCSize-1].GetNumberOfPixels();
2615 // Add the BC2 to the weights
2618 // Move further in the weights pointer
2619 weights+=m_SecondRegion.GetNumberOfPixels();
2621 for ( j = 0; j < SpaceDimension; j++ )
2623 for ( b=0; b < m_BC2Size; b++ )
2625 while ( ! (m_BC2Iterators[j][b]).IsAtEnd() )
2627 //copy weight to jacobian image
2628 m_ZeroVector[j]=(*weights) * m_BC2Values[b];
2629 (m_BC2Iterators[j][b]).Value()+= m_ZeroVector;
2630 ++(m_BC2Iterators[j][b]);
2633 weights-=m_BC2Regions[b].GetNumberOfPixels();
2635 m_ZeroVector[j]=itk::NumericTraits<JacobianValueType>::Zero;
2637 // Move further in the weights pointer
2638 weights+=m_BC2Regions[m_BC2Size-1].GetNumberOfPixels();
2644 // For each dimension, copy the weight to the support region
2645 //copy weight to jacobian image
2646 for ( j = 0; j < SpaceDimension; j++ )
2648 while ( ! (m_ThirdIterator[j]).IsAtEnd() )
2650 m_ZeroVector[j]=*weights;
2651 (m_ThirdIterator[j]).Value()+= m_ZeroVector;
2652 ++(m_ThirdIterator[j]);
2656 // Move further in the weights pointer?
2657 if (j != SpaceDimension-1) weights-=m_ThirdRegion.GetNumberOfPixels();
2658 m_ZeroVector[j]=itk::NumericTraits<JacobianValueType>::Zero;
2662 // Add the BC3 to the weights
2665 for ( j = 0; j < SpaceDimension; j++ )
2667 for ( b=0; b < m_BC3Size; b++ )
2669 while ( ! (m_BC3Iterators[j][b]).IsAtEnd() )
2671 //copy weight to jacobian image
2672 m_ZeroVector[j]=(*weights) * m_BC3Values[b];
2673 (m_BC3Iterators[j][b]).Value()+= m_ZeroVector;
2674 ++(m_BC3Iterators[j][b]);
2677 weights-=m_BC3Regions[b].GetNumberOfPixels();
2679 m_ZeroVector[j]=itk::NumericTraits<JacobianValueType>::Zero;
2683 // Return the result
2684 #if ITK_VERSION_MAJOR >= 4
2685 jacobian = m_SharedDataBSplineJacobian;
2687 return this->m_Jacobian;
2692 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
2694 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2695 ::WrapRegion( const RegionType & m_SupportRegion,
2696 RegionType & m_FirstRegion,
2697 RegionType & m_SecondRegion,
2698 RegionType & m_ThirdRegion,
2699 std::vector<RegionType>& m_BCRegions,std::vector<double>& m_BCValues,
2700 std::vector<RegionType>& m_BC2Regions,std::vector<double>& m_BC2Values,
2701 std::vector<RegionType>& m_BC3Regions,std::vector<double>& m_BC3Values,
2702 unsigned int& m_InitialOffset ) const
2705 // Synchronize regions
2707 m_FirstRegion=m_SupportRegion;
2708 m_BCRegion=m_SupportRegion;
2709 m_BCRegion.SetSize(InputDimension-1,1);
2710 m_SecondRegion=m_SupportRegion;
2711 m_ThirdRegion=m_SupportRegion;
2712 m_ThirdRegion.SetSize(InputDimension-1,0);
2713 m_BC3Regions.resize(0);
2716 // BC depends on shape
2717 switch(m_TransformShape)
2722 2: rabbit 4 CP 3 DOF
2723 3: rabbit 5 CP 4 DOF
2724 4: sputnik 4 CP 4 DOF
2725 5: sputnik 5 CP 5 DOF
2726 6: diamond 6 CP 5 DOF
2727 7: diamond 7 CP 6 DOF
2732 // ----------------------------------------------------------------------
2733 // The egg with 4 internal CP (starting from inhale)
2734 // Periodic, constrained to zero at the reference
2735 // at position 3 and
2736 // Coeff row BC1 BC2 0 BC3 1 2 BC4
2737 // PaddedCoeff R: 0 1 2 3 4 5 6
2740 // BC3= -weights[2]/weights[0] ( R2+R4 )
2742 // ---------------------------------------------------------------------
2743 switch(m_SupportRegion.GetIndex(InputDimension-1))
2748 m_FirstRegion.SetSize(InputDimension-1,2);
2749 m_FirstRegion.SetIndex(InputDimension-1,1);
2752 m_BCRegions.resize(0);
2755 m_SecondRegion.SetSize(InputDimension-1,1);
2756 m_SecondRegion.SetIndex(InputDimension-1,0);
2759 m_BC2Regions.resize(2);
2760 m_BC2Values.resize(2);
2761 m_BCRegion.SetIndex(InputDimension-1,0);
2762 m_BC2Regions[0]=m_BCRegion;
2763 m_BC2Values[0]=-m_WeightRatio[2][0];
2764 m_BCRegion.SetIndex(InputDimension-1,1);
2765 m_BC2Regions[1]=m_BCRegion;
2766 m_BC2Values[1]= -m_WeightRatio[2][0];
2773 m_FirstRegion.SetSize(InputDimension-1,1);
2774 m_FirstRegion.SetIndex(InputDimension-1,2);
2777 m_BCRegions.resize(0);
2780 m_SecondRegion.SetSize(InputDimension-1,1);
2781 m_SecondRegion.SetIndex(InputDimension-1,0);
2784 m_BC2Regions.resize(2);
2785 m_BC2Values.resize(2);
2786 m_BCRegion.SetIndex(InputDimension-1,0);
2787 m_BC2Regions[0]=m_BCRegion;
2788 m_BC2Values[0]=-m_WeightRatio[2][0];
2789 m_BCRegion.SetIndex(InputDimension-1,1);
2790 m_BC2Regions[1]=m_BCRegion;
2791 m_BC2Values[1]= -m_WeightRatio[2][0];
2794 m_FirstRegion.SetSize(InputDimension-1,1);
2795 m_FirstRegion.SetIndex(InputDimension-1,1);
2802 m_FirstRegion.SetSize(InputDimension-1,1);
2803 m_FirstRegion.SetIndex(InputDimension-1,0);
2806 m_BCRegions.resize(2);
2807 m_BCValues.resize(2);
2808 m_BCRegion.SetIndex(InputDimension-1,0);
2809 m_BCRegions[0]=m_BCRegion;
2810 m_BCValues[0]=-m_WeightRatio[2][0];
2811 m_BCRegion.SetIndex(InputDimension-1,1);
2812 m_BCRegions[1]=m_BCRegion;
2813 m_BCValues[1]= -m_WeightRatio[2][0];
2816 m_SecondRegion.SetSize(InputDimension-1,2);
2817 m_SecondRegion.SetIndex(InputDimension-1,1);
2820 m_BC2Regions.resize(0);
2827 m_FirstRegion.SetSize(InputDimension-1,0);
2830 m_BCRegions.resize(2);
2831 m_BCValues.resize(2);
2832 m_BCRegion.SetIndex(InputDimension-1,0);
2833 m_BCRegions[0]=m_BCRegion;
2834 m_BCValues[0]=-m_WeightRatio[2][0];
2835 m_BCRegion.SetIndex(InputDimension-1,1);
2836 m_BCRegions[1]=m_BCRegion;
2837 m_BCValues[1]= -m_WeightRatio[2][0];
2840 m_SecondRegion.SetSize(InputDimension-1,2);
2841 m_SecondRegion.SetIndex(InputDimension-1,1);
2844 m_BC2Regions.resize(1);
2845 m_BC2Values.resize(1);
2846 m_BCRegion.SetIndex(InputDimension-1,0);
2847 m_BC2Regions[0]=m_BCRegion;
2855 DD("supportindex > 3 ???");
2856 DD(m_SupportRegion.GetIndex(InputDimension-1));
2858 } // end switch index
2865 // ----------------------------------------------------------------------
2866 // The egg with 5 internal CP (starting from inhale)
2867 // Periodic, constrained to zero at the reference
2868 // at position 3 and
2869 // Coeff row BC1 BC2 0 BC3 1 2 3 BC4
2870 // PaddedCoeff R: 0 1 2 3 4 5 6 7
2873 // BC3= -weights[3]/weights[1] ( R2+R5 ) - R4
2875 // ---------------------------------------------------------------------
2876 switch(m_SupportRegion.GetIndex(InputDimension-1))
2881 m_FirstRegion.SetSize(InputDimension-1,2);
2882 m_FirstRegion.SetIndex(InputDimension-1,2);
2885 m_BCRegions.resize(0);
2888 m_SecondRegion.SetSize(InputDimension-1,1);
2889 m_SecondRegion.SetIndex(InputDimension-1,0);
2892 m_BC2Regions.resize(3);
2893 m_BC2Values.resize(3);
2894 m_BCRegion.SetIndex(InputDimension-1,0);
2895 m_BC2Regions[0]=m_BCRegion;
2896 m_BC2Values[0]=-m_WeightRatio[3][1];
2897 m_BCRegion.SetIndex(InputDimension-1,1);
2898 m_BC2Regions[1]=m_BCRegion;
2900 m_BCRegion.SetIndex(InputDimension-1,2);
2901 m_BC2Regions[2]=m_BCRegion;
2902 m_BC2Values[2]=-m_WeightRatio[3][1];
2909 m_FirstRegion.SetSize(InputDimension-1,1);
2910 m_FirstRegion.SetIndex(InputDimension-1,3);
2913 m_BCRegions.resize(0);
2916 m_SecondRegion.SetSize(InputDimension-1,1);
2917 m_SecondRegion.SetIndex(InputDimension-1,0);
2920 m_BC2Regions.resize(3);
2921 m_BC2Values.resize(3);
2922 m_BCRegion.SetIndex(InputDimension-1,0);
2923 m_BC2Regions[0]=m_BCRegion;
2924 m_BC2Values[0]=-m_WeightRatio[3][1];
2925 m_BCRegion.SetIndex(InputDimension-1,1);
2926 m_BC2Regions[1]=m_BCRegion;
2928 m_BCRegion.SetIndex(InputDimension-1,2);
2929 m_BC2Regions[2]=m_BCRegion;
2930 m_BC2Values[2]=-m_WeightRatio[3][1];
2933 m_FirstRegion.SetSize(InputDimension-1,1);
2934 m_FirstRegion.SetIndex(InputDimension-1,1);
2941 m_FirstRegion.SetSize(InputDimension-1,1);
2942 m_FirstRegion.SetIndex(InputDimension-1,0);
2945 m_BCRegions.resize(3);
2946 m_BCValues.resize(3);
2947 m_BCRegion.SetIndex(InputDimension-1,0);
2948 m_BCRegions[0]=m_BCRegion;
2949 m_BCValues[0]=-m_WeightRatio[3][1];
2950 m_BCRegion.SetIndex(InputDimension-1,1);
2951 m_BCRegions[1]=m_BCRegion;
2953 m_BCRegion.SetIndex(InputDimension-1,2);
2954 m_BCRegions[2]=m_BCRegion;
2955 m_BCValues[2]=-m_WeightRatio[3][1];
2958 m_SecondRegion.SetSize(InputDimension-1,2);
2959 m_SecondRegion.SetIndex(InputDimension-1,1);
2962 m_BC2Regions.resize(0);
2969 m_FirstRegion.SetSize(InputDimension-1,0);
2972 m_BCRegions.resize(3);
2973 m_BCValues.resize(3);
2974 m_BCRegion.SetIndex(InputDimension-1,0);
2975 m_BCRegions[0]=m_BCRegion;
2976 m_BCValues[0]=-m_WeightRatio[3][1];
2977 m_BCRegion.SetIndex(InputDimension-1,1);
2978 m_BCRegions[1]=m_BCRegion;
2980 m_BCRegion.SetIndex(InputDimension-1,2);
2981 m_BCRegions[2]=m_BCRegion;
2982 m_BCValues[2]=-m_WeightRatio[3][1];
2985 m_SecondRegion.SetSize(InputDimension-1,3);
2986 m_SecondRegion.SetIndex(InputDimension-1,1);
2989 m_BC2Regions.resize(0);
2996 m_FirstRegion.SetSize(InputDimension-1,3);
2997 m_FirstRegion.SetIndex(InputDimension-1,1);
3000 m_BCRegions.resize(0);
3003 m_SecondRegion.SetSize(InputDimension-1,1);
3004 m_SecondRegion.SetIndex(InputDimension-1,0);
3007 m_BC2Regions.resize(0);
3014 DD("supportindex > 3 ???");
3015 DD(m_SupportRegion.GetIndex(InputDimension-1));
3017 } // end switch index
3021 // // ----------------------------------------------------------------------
3022 // // The egg with 5 internal CP:
3023 // // Periodic, C2 smooth everywhere and constrained to zero at the reference
3024 // // Coeff row R5 BC1 0 1 2 3 BC2 R2
3025 // // PaddedCoeff R: 0 1 2 3 4 5 6 7
3026 // // BC1= -weights[2]/weights[0] ( R2+R5)
3028 // // ---------------------------------------------------------------------
3029 // switch(m_SupportRegion.GetIndex(InputDimension-1))
3034 // m_FirstRegion.SetSize(InputDimension-1,1);
3035 // m_FirstRegion.SetIndex(InputDimension-1,3);
3038 // m_BCRegions.resize(2);
3039 // m_BCValues.resize(2);
3040 // m_BCRegion.SetIndex(InputDimension-1,0);
3041 // m_BCRegions[0]=m_BCRegion;
3042 // m_BCValues[0]=-m_WeightRatio[2][0];
3043 // m_BCRegion.SetIndex(InputDimension-1,3);
3044 // m_BCRegions[1]=m_BCRegion;
3045 // m_BCValues[1]=-m_WeightRatio[2][0];
3048 // m_SecondRegion.SetSize(InputDimension-1,2);
3049 // m_SecondRegion.SetIndex(InputDimension-1,0);
3052 // m_BC2Regions.resize(0);
3059 // m_FirstRegion.SetSize(InputDimension-1,0);
3062 // m_BCRegions.resize(2);
3063 // m_BCValues.resize(2);
3064 // m_BCRegion.SetIndex(InputDimension-1,0);
3065 // m_BCRegions[0]=m_BCRegion;
3066 // m_BCValues[0]=-m_WeightRatio[2][0];
3067 // m_BCRegion.SetIndex(InputDimension-1,3);
3068 // m_BCRegions[1]=m_BCRegion;
3069 // m_BCValues[1]=-m_WeightRatio[2][0];
3072 // m_SecondRegion.SetSize(InputDimension-1,3);
3073 // m_SecondRegion.SetIndex(InputDimension-1,0);
3076 // m_BC2Regions.resize(0);
3083 // m_FirstRegion.SetSize(InputDimension-1,4);
3084 // m_FirstRegion.SetIndex(InputDimension-1,0);
3087 // m_BCRegions.resize(0);
3090 // m_SecondRegion.SetSize(InputDimension-1,0);
3093 // m_BC2Regions.resize(0);
3100 // m_FirstRegion.SetSize(InputDimension-1,3);
3101 // m_FirstRegion.SetIndex(InputDimension-1,1);
3104 // m_BCRegions.resize(2);
3105 // m_BCValues.resize(2);
3106 // m_BCRegion.SetIndex(InputDimension-1,0);
3107 // m_BCRegions[0]=m_BCRegion;
3108 // m_BCValues[0]=-m_WeightRatio[2][0];
3109 // m_BCRegion.SetIndex(InputDimension-1,3);
3110 // m_BCRegions[1]=m_BCRegion;
3111 // m_BCValues[1]=-m_WeightRatio[2][0];
3114 // m_SecondRegion.SetSize(InputDimension-1,0);
3117 // m_BC2Regions.resize(0);
3125 // m_FirstRegion.SetSize(InputDimension-1,2);
3126 // m_FirstRegion.SetIndex(InputDimension-1,2);
3129 // m_BCRegions.resize(2);
3130 // m_BCValues.resize(2);
3131 // m_BCRegion.SetIndex(InputDimension-1,0);
3132 // m_BCRegions[0]=m_BCRegion;
3133 // m_BCValues[0]=-m_WeightRatio[2][0];
3134 // m_BCRegion.SetIndex(InputDimension-1,3);
3135 // m_BCRegions[1]=m_BCRegion;
3136 // m_BCValues[1]=-m_WeightRatio[2][0];
3139 // m_SecondRegion.SetSize(InputDimension-1,1);
3140 // m_SecondRegion.SetIndex(InputDimension-1,0);
3143 // m_BC2Regions.resize(0);
3150 // DD("supportindex > 4 ???");
3151 // DD(m_SupportRegion.GetIndex(InputDimension-1));
3152 // DD(m_TransformShape);
3154 // }// end swith index
3157 // } // end case 1 shape
3161 // ----------------------------------------------------------------------
3162 // The rabbit with 4 internal CP
3163 // Periodic, constrained to zero at the reference
3164 // at position 3 and the extremes fixed with anit-symm bc
3165 // Coeff row BC1 0 1 BC2 2 BC3 BC4
3166 // PaddedCoeff R: 0 1 2 3 4 5 6
3168 // BC2= -weights[2]/weights[0] ( R2+R4 )
3171 // ---------------------------------------------------------------------
3172 switch(m_SupportRegion.GetIndex(InputDimension-1))
3177 m_FirstRegion.SetSize(InputDimension-1,0);
3180 m_BCRegions.resize(2);
3181 m_BCValues.resize(2);
3182 m_BCRegion.SetIndex(InputDimension-1,0);
3183 m_BCRegions[0]=m_BCRegion;
3185 m_BCRegion.SetIndex(InputDimension-1,1);
3186 m_BCRegions[1]=m_BCRegion;
3190 m_SecondRegion.SetSize(InputDimension-1,2);
3191 m_SecondRegion.SetIndex(InputDimension-1,0);
3194 m_BC2Regions.resize(2);
3195 m_BC2Values.resize(2);
3196 m_BCRegion.SetIndex(InputDimension-1,1);
3197 m_BC2Regions[0]=m_BCRegion;
3198 m_BC2Values[0]=-m_WeightRatio[2][0];
3199 m_BCRegion.SetIndex(InputDimension-1,2);
3200 m_BC2Regions[1]=m_BCRegion;
3201 m_BC2Values[1]= -m_WeightRatio[2][0];
3208 m_FirstRegion.SetSize(InputDimension-1,2);
3209 m_FirstRegion.SetIndex(InputDimension-1,0);
3212 m_BCRegions.resize(2);
3213 m_BCValues.resize(2);
3214 m_BCRegion.SetIndex(InputDimension-1,1);
3215 m_BCRegions[0]=m_BCRegion;
3216 m_BCValues[0]=-m_WeightRatio[2][0];
3217 m_BCRegion.SetIndex(InputDimension-1,2);
3218 m_BCRegions[1]=m_BCRegion;
3219 m_BCValues[1]= -m_WeightRatio[2][0];
3222 m_SecondRegion.SetSize(InputDimension-1,1);
3223 m_SecondRegion.SetIndex(InputDimension-1,2);
3226 m_BC2Regions.resize(0);
3233 m_FirstRegion.SetSize(InputDimension-1,1);
3234 m_FirstRegion.SetIndex(InputDimension-1,1);
3237 m_BCRegions.resize(2);
3238 m_BCValues.resize(2);
3239 m_BCRegion.SetIndex(InputDimension-1,1);
3240 m_BCRegions[0]=m_BCRegion;
3241 m_BCValues[0]=-m_WeightRatio[2][0];
3242 m_BCRegion.SetIndex(InputDimension-1,2);
3243 m_BCRegions[1]=m_BCRegion;
3244 m_BCValues[1]= -m_WeightRatio[2][0];
3247 m_SecondRegion.SetSize(InputDimension-1,1);
3248 m_SecondRegion.SetIndex(InputDimension-1,2);
3251 m_BC2Regions.resize(1);
3252 m_BC2Values.resize(1);
3253 m_BCRegion.SetIndex(InputDimension-1,0);
3254 m_BC2Regions[0]=m_BCRegion;
3262 m_FirstRegion.SetSize(InputDimension-1,0);
3265 m_BCRegions.resize(2);
3266 m_BCValues.resize(2);
3267 m_BCRegion.SetIndex(InputDimension-1,1);
3268 m_BCRegions[0]=m_BCRegion;
3269 m_BCValues[0]=-m_WeightRatio[2][0];
3270 m_BCRegion.SetIndex(InputDimension-1,2);
3271 m_BCRegions[1]=m_BCRegion;
3272 m_BCValues[1]= -m_WeightRatio[2][0];
3275 m_SecondRegion.SetSize(InputDimension-1,1);
3276 m_SecondRegion.SetIndex(InputDimension-1,2);
3279 m_BC2Regions.resize(1);
3280 m_BC2Values.resize(1);
3281 m_BCRegion.SetIndex(InputDimension-1,0);
3282 m_BC2Regions[0]=m_BCRegion;
3286 m_BC3Regions.resize(2);
3287 m_BC3Values.resize(2);
3288 m_BCRegion.SetIndex(InputDimension-1,0);
3289 m_BC3Regions[0]=m_BCRegion;
3291 m_BCRegion.SetIndex(InputDimension-1,2);
3292 m_BC3Regions[1]=m_BCRegion;
3300 DD("supportindex > 3 ???");
3301 DD(m_SupportRegion.GetIndex(InputDimension-1));
3303 } // end switch index
3306 } // end case 2 shape
3309 // ----------------------------------------------------------------------
3310 // The rabbit with 5 internal CP
3311 // Periodic, constrained to zero at the reference at position 3.5
3312 // and the extremes fixed with anti-symmetrical BC
3313 // Coeff row BC1 0 1 BC2 2 3 BC3 BC4
3314 // PaddedCoeff R: 0 1 2 3 4 5 6 7
3316 // BC2= -weights[3]/weights[1] ( R2+R5 ) - R4
3319 // ---------------------------------------------------------------------
3320 switch(m_SupportRegion.GetIndex(InputDimension-1))
3325 m_FirstRegion.SetSize(InputDimension-1,0);
3328 m_BCRegions.resize(2);
3329 m_BCValues.resize(2);
3330 m_BCRegion.SetIndex(InputDimension-1,0);
3331 m_BCRegions[0]=m_BCRegion;
3333 m_BCRegion.SetIndex(InputDimension-1,1);
3334 m_BCRegions[1]=m_BCRegion;
3338 m_SecondRegion.SetSize(InputDimension-1,2);
3339 m_SecondRegion.SetIndex(InputDimension-1,0);
3342 m_BC2Regions.resize(3);
3343 m_BC2Values.resize(3);
3344 m_BCRegion.SetIndex(InputDimension-1,1);
3345 m_BC2Regions[0]=m_BCRegion;
3346 m_BC2Values[0]=-m_WeightRatio[3][1];
3347 m_BCRegion.SetIndex(InputDimension-1,2);
3348 m_BC2Regions[1]=m_BCRegion;
3350 m_BCRegion.SetIndex(InputDimension-1,3);
3351 m_BC2Regions[2]=m_BCRegion;
3352 m_BC2Values[2]=-m_WeightRatio[3][1];
3359 m_FirstRegion.SetSize(InputDimension-1,2);
3360 m_FirstRegion.SetIndex(InputDimension-1,0);
3363 m_BCRegions.resize(3);
3364 m_BCValues.resize(3);
3365 m_BCRegion.SetIndex(InputDimension-1,1);
3366 m_BCRegions[0]=m_BCRegion;
3367 m_BCValues[0]=-m_WeightRatio[3][1];
3368 m_BCRegion.SetIndex(InputDimension-1,2);
3369 m_BCRegions[1]=m_BCRegion;
3371 m_BCRegion.SetIndex(InputDimension-1,3);
3372 m_BCRegions[2]=m_BCRegion;
3373 m_BCValues[2]=-m_WeightRatio[3][1];
3376 m_SecondRegion.SetSize(InputDimension-1,1);
3377 m_SecondRegion.SetIndex(InputDimension-1,2);
3380 m_BC2Regions.resize(0);
3387 m_FirstRegion.SetSize(InputDimension-1,1);
3388 m_FirstRegion.SetIndex(InputDimension-1,1);
3391 m_BCRegions.resize(3);
3392 m_BCValues.resize(3);
3393 m_BCRegion.SetIndex(InputDimension-1,1);
3394 m_BCRegions[0]=m_BCRegion;
3395 m_BCValues[0]=-m_WeightRatio[3][1];
3396 m_BCRegion.SetIndex(InputDimension-1,2);
3397 m_BCRegions[1]=m_BCRegion;
3399 m_BCRegion.SetIndex(InputDimension-1,3);
3400 m_BCRegions[2]=m_BCRegion;
3401 m_BCValues[2]=-m_WeightRatio[3][1];
3404 m_SecondRegion.SetSize(InputDimension-1,2);
3405 m_SecondRegion.SetIndex(InputDimension-1,2);
3408 m_BC2Regions.resize(0);
3415 m_FirstRegion.SetSize(InputDimension-1,0);
3418 m_BCRegions.resize(3);
3419 m_BCValues.resize(3);
3420 m_BCRegion.SetIndex(InputDimension-1,1);
3421 m_BCRegions[0]=m_BCRegion;
3422 m_BCValues[0]=-m_WeightRatio[3][1];
3423 m_BCRegion.SetIndex(InputDimension-1,2);
3424 m_BCRegions[1]=m_BCRegion;
3426 m_BCRegion.SetIndex(InputDimension-1,3);
3427 m_BCRegions[2]=m_BCRegion;
3428 m_BCValues[2]=-m_WeightRatio[3][1];
3431 m_SecondRegion.SetSize(InputDimension-1,2);
3432 m_SecondRegion.SetIndex(InputDimension-1,2);
3435 m_BC2Regions.resize(1);
3436 m_BC2Values.resize(1);
3437 m_BCRegion.SetIndex(InputDimension-1,0);
3438 m_BC2Regions[0]=m_BCRegion;
3447 m_FirstRegion.SetSize(InputDimension-1,2);
3448 m_FirstRegion.SetIndex(InputDimension-1,2);
3451 m_BCRegions.resize(1);
3452 m_BCValues.resize(1);
3453 m_BCRegion.SetIndex(InputDimension-1,0);
3454 m_BCRegions[0]=m_BCRegion;
3458 m_SecondRegion.SetSize(InputDimension-1,0);
3461 m_BC2Regions.resize(2);
3462 m_BC2Values.resize(2);
3463 m_BCRegion.SetIndex(InputDimension-1,0);
3464 m_BC2Regions[0]=m_BCRegion;
3466 m_BCRegion.SetIndex(InputDimension-1,3);
3467 m_BC2Regions[1]=m_BCRegion;
3475 DD("supportindex > 4 ???");
3476 DD(m_SupportRegion.GetIndex(InputDimension-1));
3478 } // end switch index
3481 } // end case 3 shape
3485 // ----------------------------------------------------------------------
3486 // The sputnik with 4 internal CP
3487 // Periodic, constrained to zero at the reference
3488 // at position 3 and one indepent extremes copied
3489 // Coeff row BC1 0 1 BC2 2 BC2 3
3490 // PaddedCoeff R: 0 1 2 3 4 5 6
3492 // BC2= -weights[2]/weights[0] ( R2+R4 )
3493 // BC3= weights[2]/weights[0] ( R2-R4) + R1
3494 // ---------------------------------------------------------------------
3495 switch(m_SupportRegion.GetIndex(InputDimension-1))
3500 m_FirstRegion.SetSize(InputDimension-1,0);
3503 m_BCRegions.resize(1);
3504 m_BCValues.resize(1);
3505 m_BCRegion.SetIndex(InputDimension-1,3);
3506 m_BCRegions[0]=m_BCRegion;
3510 m_SecondRegion.SetSize(InputDimension-1,2);
3511 m_SecondRegion.SetIndex(InputDimension-1,0);
3514 m_BC2Regions.resize(2);
3515 m_BC2Values.resize(2);
3516 m_BCRegion.SetIndex(InputDimension-1,1);
3517 m_BC2Regions[0]=m_BCRegion;
3518 m_BC2Values[0]=-m_WeightRatio[2][0];
3519 m_BCRegion.SetIndex(InputDimension-1,2);
3520 m_BC2Regions[1]=m_BCRegion;
3521 m_BC2Values[1]= -m_WeightRatio[2][0];
3528 m_FirstRegion.SetSize(InputDimension-1,2);
3529 m_FirstRegion.SetIndex(InputDimension-1,0);
3532 m_BCRegions.resize(2);
3533 m_BCValues.resize(2);
3534 m_BCRegion.SetIndex(InputDimension-1,1);
3535 m_BCRegions[0]=m_BCRegion;
3536 m_BCValues[0]=-m_WeightRatio[2][0];
3537 m_BCRegion.SetIndex(InputDimension-1,2);
3538 m_BCRegions[1]=m_BCRegion;
3539 m_BCValues[1]= -m_WeightRatio[2][0];
3542 m_SecondRegion.SetSize(InputDimension-1,1);
3543 m_SecondRegion.SetIndex(InputDimension-1,2);
3546 m_BC2Regions.resize(0);
3553 m_FirstRegion.SetSize(InputDimension-1,1);
3554 m_FirstRegion.SetIndex(InputDimension-1,1);
3557 m_BCRegions.resize(2);
3558 m_BCValues.resize(2);
3559 m_BCRegion.SetIndex(InputDimension-1,1);
3560 m_BCRegions[0]=m_BCRegion;
3561 m_BCValues[0]=-m_WeightRatio[2][0];
3562 m_BCRegion.SetIndex(InputDimension-1,2);
3563 m_BCRegions[1]=m_BCRegion;
3564 m_BCValues[1]= -m_WeightRatio[2][0];
3567 m_SecondRegion.SetSize(InputDimension-1,1);
3568 m_SecondRegion.SetIndex(InputDimension-1,2);
3571 m_BC2Regions.resize(3);
3572 m_BC2Values.resize(3);
3573 m_BCRegion.SetIndex(InputDimension-1,0);
3574 m_BC2Regions[0]=m_BCRegion;
3576 m_BCRegion.SetIndex(InputDimension-1,1);
3577 m_BC2Regions[1]=m_BCRegion;
3578 m_BC2Values[1]=m_WeightRatio[2][0];
3579 m_BCRegion.SetIndex(InputDimension-1,2);
3580 m_BC2Regions[2]=m_BCRegion;
3581 m_BC2Values[2]=-m_WeightRatio[2][0];
3588 m_FirstRegion.SetSize(InputDimension-1,0);
3591 m_BCRegions.resize(2);
3592 m_BCValues.resize(2);
3593 m_BCRegion.SetIndex(InputDimension-1,1);
3594 m_BCRegions[0]=m_BCRegion;
3595 m_BCValues[0]=-m_WeightRatio[2][0];
3596 m_BCRegion.SetIndex(InputDimension-1,2);
3597 m_BCRegions[1]=m_BCRegion;
3598 m_BCValues[1]= -m_WeightRatio[2][0];
3601 m_SecondRegion.SetSize(InputDimension-1,1);
3602 m_SecondRegion.SetIndex(InputDimension-1,2);
3605 m_BC2Regions.resize(3);
3606 m_BC2Values.resize(3);
3607 m_BCRegion.SetIndex(InputDimension-1,0);
3608 m_BC2Regions[0]=m_BCRegion;
3610 m_BCRegion.SetIndex(InputDimension-1,1);
3611 m_BC2Regions[1]=m_BCRegion;
3612 m_BC2Values[1]=m_WeightRatio[2][0];
3613 m_BCRegion.SetIndex(InputDimension-1,2);
3614 m_BC2Regions[2]=m_BCRegion;
3615 m_BC2Values[2]=-m_WeightRatio[2][0];
3618 m_ThirdRegion.SetSize(InputDimension-1,1);
3619 m_ThirdRegion.SetIndex(InputDimension-1,3);
3626 DD("supportindex > 3 ???");
3627 DD(m_SupportRegion.GetIndex(InputDimension-1));
3629 } // end switch index
3632 } // end case 4 shape
3636 // ----------------------------------------------------------------------
3637 // The sputnik with 5 internal CP
3638 // Periodic, constrained to zero at the reference
3639 // at position 3 and one indepent extreme
3640 // Coeff row BC1 0 1 BC2 2 3 BC3 4
3641 // PaddedCoeff R: 0 1 2 3 4 5 6 7
3642 // BC1= R2 + R5 - R7
3643 // BC2= -weights[3]/weights[1] ( R2+R5 ) - R4
3644 // BC3= R1 + 0.5 R2 - 0.5 R7
3645 // ----------------------------------------------------------------------
3646 switch(m_SupportRegion.GetIndex(InputDimension-1))
3651 m_FirstRegion.SetSize(InputDimension-1,0);
3654 m_BCRegions.resize(3);
3655 m_BCValues.resize(3);
3656 m_BCRegion.SetIndex(InputDimension-1,1);
3657 m_BCRegions[0]=m_BCRegion;
3659 m_BCRegion.SetIndex(InputDimension-1,3);
3660 m_BCRegions[1]=m_BCRegion;
3662 m_BCRegion.SetIndex(InputDimension-1,4);
3663 m_BCRegions[2]=m_BCRegion;
3667 m_SecondRegion.SetSize(InputDimension-1,2);
3668 m_SecondRegion.SetIndex(InputDimension-1,0);
3671 m_BC2Regions.resize(3);
3672 m_BC2Values.resize(3);
3673 m_BCRegion.SetIndex(InputDimension-1,1);
3674 m_BC2Regions[0]=m_BCRegion;
3675 m_BC2Values[0]=-m_WeightRatio[3][1];
3676 m_BCRegion.SetIndex(InputDimension-1,2);
3677 m_BC2Regions[1]=m_BCRegion;
3679 m_BCRegion.SetIndex(InputDimension-1,3);
3680 m_BC2Regions[2]=m_BCRegion;
3681 m_BC2Values[2]=-m_WeightRatio[3][1];
3688 m_FirstRegion.SetSize(InputDimension-1,2);
3689 m_FirstRegion.SetIndex(InputDimension-1,0);
3692 m_BCRegions.resize(3);
3693 m_BCValues.resize(3);
3694 m_BCRegion.SetIndex(InputDimension-1,1);
3695 m_BCRegions[0]=m_BCRegion;
3696 m_BCValues[0]=-m_WeightRatio[3][1];
3697 m_BCRegion.SetIndex(InputDimension-1,2);
3698 m_BCRegions[1]=m_BCRegion;
3700 m_BCRegion.SetIndex(InputDimension-1,3);
3701 m_BCRegions[2]=m_BCRegion;
3702 m_BCValues[2]=-m_WeightRatio[3][1];
3705 m_SecondRegion.SetSize(InputDimension-1,1);
3706 m_SecondRegion.SetIndex(InputDimension-1,2);
3709 m_BC2Regions.resize(0);
3716 m_FirstRegion.SetSize(InputDimension-1,1);
3717 m_FirstRegion.SetIndex(InputDimension-1,1);
3720 m_BCRegions.resize(3);
3721 m_BCValues.resize(3);
3722 m_BCRegion.SetIndex(InputDimension-1,1);
3723 m_BCRegions[0]=m_BCRegion;
3724 m_BCValues[0]=-m_WeightRatio[3][1];
3725 m_BCRegion.SetIndex(InputDimension-1,2);
3726 m_BCRegions[1]=m_BCRegion;
3728 m_BCRegion.SetIndex(InputDimension-1,3);
3729 m_BCRegions[2]=m_BCRegion;
3730 m_BCValues[2]=-m_WeightRatio[3][1];
3733 m_SecondRegion.SetSize(InputDimension-1,2);
3734 m_SecondRegion.SetIndex(InputDimension-1,2);
3737 m_BC2Regions.resize(0);
3744 m_FirstRegion.SetSize(InputDimension-1,0);
3747 m_BCRegions.resize(3);
3748 m_BCValues.resize(3);
3749 m_BCRegion.SetIndex(InputDimension-1,1);
3750 m_BCRegions[0]=m_BCRegion;
3751 m_BCValues[0]=-m_WeightRatio[3][1];
3752 m_BCRegion.SetIndex(InputDimension-1,2);
3753 m_BCRegions[1]=m_BCRegion;
3755 m_BCRegion.SetIndex(InputDimension-1,3);
3756 m_BCRegions[2]=m_BCRegion;
3757 m_BCValues[2]=-m_WeightRatio[3][1];
3760 m_SecondRegion.SetSize(InputDimension-1,2);
3761 m_SecondRegion.SetIndex(InputDimension-1,2);
3764 m_BC2Regions.resize(3);
3765 m_BC2Values.resize(3);
3766 m_BCRegion.SetIndex(InputDimension-1,0);
3767 m_BC2Regions[0]=m_BCRegion;
3769 m_BCRegion.SetIndex(InputDimension-1,1);
3770 m_BC2Regions[1]=m_BCRegion;
3772 m_BCRegion.SetIndex(InputDimension-1,4);
3773 m_BC2Regions[2]=m_BCRegion;
3774 m_BC2Values[2]=-0.5;
3781 m_FirstRegion.SetSize(InputDimension-1,2);
3782 m_FirstRegion.SetIndex(InputDimension-1,2);
3785 m_BCRegions.resize(3);
3786 m_BCValues.resize(3);
3787 m_BCRegion.SetIndex(InputDimension-1,0);
3788 m_BCRegions[0]=m_BCRegion;
3790 m_BCRegion.SetIndex(InputDimension-1,1);
3791 m_BCRegions[1]=m_BCRegion;
3793 m_BCRegion.SetIndex(InputDimension-1,4);
3794 m_BCRegions[2]=m_BCRegion;
3798 m_SecondRegion.SetSize(InputDimension-1,1);
3799 m_SecondRegion.SetIndex(InputDimension-1,4);
3802 m_BC2Regions.resize(0);
3808 DD("supportindex > 4 ???");
3809 DD(m_SupportRegion.GetIndex(InputDimension-1));
3811 } // end switch index
3814 } // end case 5 shape
3821 // ----------------------------------------------------------------------
3822 // The diamond with 4 internal CP:
3823 // Periodic, constrained to zero at the reference at position 3
3824 // Coeff row 0 1 2 BC1 3 BC2 4
3825 // PaddedCoeff R:0 1 2 3 4 5 6
3826 // BC1= -weights[2]/weights[0] ( R2+R4 )
3827 // BC2= weights[2]/weights[0] ( R0+R2-R4-R6 ) + R1
3828 // ---------------------------------------------------------------------
3829 switch(m_SupportRegion.GetIndex(InputDimension-1))
3834 m_FirstRegion.SetSize(InputDimension-1,3);
3835 m_FirstRegion.SetIndex(InputDimension-1,0);
3838 m_BCRegions.resize(2);
3839 m_BCValues.resize(2);
3840 m_BCRegion.SetIndex(InputDimension-1,2);
3841 m_BCRegions[0]=m_BCRegion;
3842 m_BCValues[0]=-m_WeightRatio[2][0];
3843 m_BCRegion.SetIndex(InputDimension-1,3);
3844 m_BCRegions[1]=m_BCRegion;
3845 m_BCValues[1]=-m_WeightRatio[2][0];
3848 m_SecondRegion.SetSize(InputDimension-1,0);
3851 m_BC2Regions.resize(0);
3858 m_FirstRegion.SetSize(InputDimension-1,2);
3859 m_FirstRegion.SetIndex(InputDimension-1,1);
3862 m_BCRegions.resize(2);
3863 m_BCValues.resize(2);
3864 m_BCRegion.SetIndex(InputDimension-1,2);
3865 m_BCRegions[0]=m_BCRegion;
3866 m_BCValues[0]=-m_WeightRatio[2][0];
3867 m_BCRegion.SetIndex(InputDimension-1,3);
3868 m_BCRegions[1]=m_BCRegion;
3869 m_BCValues[1]=-m_WeightRatio[2][0];
3872 m_SecondRegion.SetSize(InputDimension-1,1);
3873 m_SecondRegion.SetIndex(InputDimension-1,3);
3876 m_BC2Regions.resize(0);
3883 m_FirstRegion.SetSize(InputDimension-1,1);
3884 m_FirstRegion.SetIndex(InputDimension-1,2);
3887 m_BCRegions.resize(2);
3888 m_BCValues.resize(2);
3889 m_BCRegion.SetIndex(InputDimension-1,2);
3890 m_BCRegions[0]=m_BCRegion;
3891 m_BCValues[0]=-m_WeightRatio[2][0];
3892 m_BCRegion.SetIndex(InputDimension-1,3);
3893 m_BCRegions[1]=m_BCRegion;
3894 m_BCValues[1]=-m_WeightRatio[2][0];
3897 m_SecondRegion.SetSize(InputDimension-1,1);
3898 m_SecondRegion.SetIndex(InputDimension-1,3);
3901 m_BC2Regions.resize(5);
3902 m_BC2Values.resize(5);
3903 m_BCRegion.SetIndex(InputDimension-1,0);
3904 m_BC2Regions[0]=m_BCRegion;
3905 m_BC2Values[0]=m_WeightRatio[2][0];
3906 m_BCRegion.SetIndex(InputDimension-1,1);
3907 m_BC2Regions[1]=m_BCRegion;
3909 m_BCRegion.SetIndex(InputDimension-1,2);
3910 m_BC2Regions[2]=m_BCRegion;
3911 m_BC2Values[2]=m_WeightRatio[2][0];
3912 m_BCRegion.SetIndex(InputDimension-1,3);
3913 m_BC2Regions[3]=m_BCRegion;
3914 m_BC2Values[3]=-m_WeightRatio[2][0];
3915 m_BCRegion.SetIndex(InputDimension-1,4);
3916 m_BC2Regions[4]=m_BCRegion;
3917 m_BC2Values[4]=-m_WeightRatio[2][0];
3924 m_FirstRegion.SetSize(InputDimension-1,0);
3927 m_BCRegions.resize(2);
3928 m_BCValues.resize(2);
3929 m_BCRegion.SetIndex(InputDimension-1,2);
3930 m_BCRegions[0]=m_BCRegion;
3931 m_BCValues[0]=-m_WeightRatio[2][0];
3932 m_BCRegion.SetIndex(InputDimension-1,3);
3933 m_BCRegions[1]=m_BCRegion;
3934 m_BCValues[1]=-m_WeightRatio[2][0];
3937 m_SecondRegion.SetSize(InputDimension-1,1);
3938 m_SecondRegion.SetIndex(InputDimension-1,3);
3941 m_BC2Regions.resize(5);
3942 m_BC2Values.resize(5);
3943 m_BCRegion.SetIndex(InputDimension-1,0);
3944 m_BC2Regions[0]=m_BCRegion;
3945 m_BC2Values[0]=m_WeightRatio[2][0];
3946 m_BCRegion.SetIndex(InputDimension-1,1);
3947 m_BC2Regions[1]=m_BCRegion;
3949 m_BCRegion.SetIndex(InputDimension-1,2);
3950 m_BC2Regions[2]=m_BCRegion;
3951 m_BC2Values[2]=m_WeightRatio[2][0];
3952 m_BCRegion.SetIndex(InputDimension-1,3);
3953 m_BC2Regions[3]=m_BCRegion;
3954 m_BC2Values[3]=-m_WeightRatio[2][0];
3955 m_BCRegion.SetIndex(InputDimension-1,4);
3956 m_BC2Regions[4]=m_BCRegion;
3957 m_BC2Values[4]=-m_WeightRatio[2][0];
3960 m_ThirdRegion.SetSize(InputDimension-1,1);
3961 m_ThirdRegion.SetIndex(InputDimension-1,4);
3968 DD("supportindex > 3 ???");
3969 DD(m_SupportRegion.GetIndex(InputDimension-1));
3971 } // end switch index
3974 } // end case 7 shape
3978 // ----------------------------------------------------------------------
3979 // The diamond with 5 internal CP:
3980 // periodic, constrained to zero at the reference at position 3.5
3981 // Coeff row 0 1 2 BC1 3 4 BC2 5
3982 // PaddedCoeff R:0 1 2 3 4 5 6 7
3983 // BC1= -weights[3]/weights[1] ( R2+R5 ) - R4
3984 // BC2= weights[2]/weights[0] ( R0+R2-R5-R7 ) + R1
3985 // ---------------------------------------------------------------------
3986 switch(m_SupportRegion.GetIndex(InputDimension-1))
3991 m_FirstRegion.SetSize(InputDimension-1,3);
3992 m_FirstRegion.SetIndex(InputDimension-1,0);
3995 m_BCRegions.resize(3);
3996 m_BCValues.resize(3);
3997 m_BCRegion.SetIndex(InputDimension-1,2);
3998 m_BCRegions[0]=m_BCRegion;
3999 m_BCValues[0]=-m_WeightRatio[3][1];
4000 m_BCRegion.SetIndex(InputDimension-1,3);
4001 m_BCRegions[1]=m_BCRegion;
4003 m_BCRegion.SetIndex(InputDimension-1,4);
4004 m_BCRegions[2]=m_BCRegion;
4005 m_BCValues[2]=-m_WeightRatio[3][1];
4008 m_SecondRegion.SetSize(InputDimension-1,0);
4011 m_BC2Regions.resize(0);
4018 m_FirstRegion.SetSize(InputDimension-1,2);
4019 m_FirstRegion.SetIndex(InputDimension-1,1);
4022 m_BCRegions.resize(3);
4023 m_BCValues.resize(3);
4024 m_BCRegion.SetIndex(InputDimension-1,2);
4025 m_BCRegions[0]=m_BCRegion;
4026 m_BCValues[0]=-m_WeightRatio[3][1];
4027 m_BCRegion.SetIndex(InputDimension-1,3);
4028 m_BCRegions[1]=m_BCRegion;
4030 m_BCRegion.SetIndex(InputDimension-1,4);
4031 m_BCRegions[2]=m_BCRegion;
4032 m_BCValues[2]=-m_WeightRatio[3][1];
4035 m_SecondRegion.SetSize(InputDimension-1,1);
4036 m_SecondRegion.SetIndex(InputDimension-1,3);
4039 m_BC2Regions.resize(0);
4046 m_FirstRegion.SetSize(InputDimension-1,1);
4047 m_FirstRegion.SetIndex(InputDimension-1,2);
4050 m_BCRegions.resize(3);
4051 m_BCValues.resize(3);
4052 m_BCRegion.SetIndex(InputDimension-1,2);
4053 m_BCRegions[0]=m_BCRegion;
4054 m_BCValues[0]=-m_WeightRatio[3][1];
4055 m_BCRegion.SetIndex(InputDimension-1,3);
4056 m_BCRegions[1]=m_BCRegion;
4058 m_BCRegion.SetIndex(InputDimension-1,4);
4059 m_BCRegions[2]=m_BCRegion;
4060 m_BCValues[2]=-m_WeightRatio[3][1];
4063 m_SecondRegion.SetSize(InputDimension-1,2);
4064 m_SecondRegion.SetIndex(InputDimension-1,3);
4067 m_BC2Regions.resize(0);
4074 m_FirstRegion.SetSize(InputDimension-1,0);
4075 m_FirstRegion.SetIndex(InputDimension-1,0);
4078 m_BCRegions.resize(3);
4079 m_BCValues.resize(3);
4080 m_BCRegion.SetIndex(InputDimension-1,2);
4081 m_BCRegions[0]=m_BCRegion;
4082 m_BCValues[0]=-m_WeightRatio[3][1];
4083 m_BCRegion.SetIndex(InputDimension-1,3);
4084 m_BCRegions[1]=m_BCRegion;
4086 m_BCRegion.SetIndex(InputDimension-1,4);
4087 m_BCRegions[2]=m_BCRegion;
4088 m_BCValues[2]=-m_WeightRatio[3][1];
4091 m_SecondRegion.SetSize(InputDimension-1,2);
4092 m_SecondRegion.SetIndex(InputDimension-1,3);
4095 m_BC2Regions.resize(5);
4096 m_BC2Values.resize(5);
4097 m_BCRegion.SetIndex(InputDimension-1,0);
4098 m_BC2Regions[0]=m_BCRegion;
4099 m_BC2Values[0]=m_WeightRatio[2][0];
4100 m_BCRegion.SetIndex(InputDimension-1,1);
4101 m_BC2Regions[1]=m_BCRegion;
4103 m_BCRegion.SetIndex(InputDimension-1,2);
4104 m_BC2Regions[2]=m_BCRegion;
4105 m_BC2Values[2]=m_WeightRatio[2][0];
4106 m_BCRegion.SetIndex(InputDimension-1,4);
4107 m_BC2Regions[3]=m_BCRegion;
4108 m_BC2Values[3]=-m_WeightRatio[2][0];
4109 m_BCRegion.SetIndex(InputDimension-1,5);
4110 m_BC2Regions[4]=m_BCRegion;
4111 m_BC2Values[4]=-m_WeightRatio[2][0];
4119 m_FirstRegion.SetSize(InputDimension-1,2);
4120 m_FirstRegion.SetIndex(InputDimension-1,3);
4123 m_BCRegions.resize(5);
4124 m_BCValues.resize(5);
4125 m_BCRegion.SetIndex(InputDimension-1,0);
4126 m_BCRegions[0]=m_BCRegion;
4127 m_BCValues[0]=m_WeightRatio[2][0];
4128 m_BCRegion.SetIndex(InputDimension-1,1);
4129 m_BCRegions[1]=m_BCRegion;
4131 m_BCRegion.SetIndex(InputDimension-1,2);
4132 m_BCRegions[2]=m_BCRegion;
4133 m_BCValues[2]=m_WeightRatio[2][0];
4134 m_BCRegion.SetIndex(InputDimension-1,4);
4135 m_BCRegions[3]=m_BCRegion;
4136 m_BCValues[3]=-m_WeightRatio[2][0];
4137 m_BCRegion.SetIndex(InputDimension-1,5);
4138 m_BCRegions[4]=m_BCRegion;
4139 m_BCValues[4]=-m_WeightRatio[2][0];
4142 m_SecondRegion.SetSize(InputDimension-1,1);
4143 m_SecondRegion.SetIndex(InputDimension-1,5);
4146 m_BC2Regions.resize(0);
4153 DD("supportindex > 4 ???");
4154 DD(m_SupportRegion.GetIndex(InputDimension-1));
4156 } // end switch index
4159 } // end case 7 shape
4163 // ----------------------------------------------------------------------
4164 // The sputnik with 5 internal CP
4165 // Periodic, constrained to zero at the reference
4166 // at position 3 and one indepent extreme
4167 // Coeff row BC1 0 1 BC2 2 3 BC3 4
4168 // PaddedCoeff R: 0 1 2 3 4 5 6 7
4170 // BC2= -weights[3]/weights[1] ( R2+R5 ) - R4
4172 // ---------------------------------------------------------------------
4173 switch(m_SupportRegion.GetIndex(InputDimension-1))
4178 m_FirstRegion.SetSize(InputDimension-1,0);
4181 m_BCRegions.resize(3);
4182 m_BCValues.resize(3);
4183 m_BCRegion.SetIndex(InputDimension-1,1);
4184 m_BCRegions[0]=m_BCRegion;
4186 m_BCRegion.SetIndex(InputDimension-1,3);
4187 m_BCRegions[1]=m_BCRegion;
4189 m_BCRegion.SetIndex(InputDimension-1,4);
4190 m_BCRegions[2]=m_BCRegion;
4194 m_SecondRegion.SetSize(InputDimension-1,2);
4195 m_SecondRegion.SetIndex(InputDimension-1,0);
4198 m_BC2Regions.resize(3);
4199 m_BC2Values.resize(3);
4200 m_BCRegion.SetIndex(InputDimension-1,1);
4201 m_BC2Regions[0]=m_BCRegion;
4202 m_BC2Values[0]=-m_WeightRatio[3][1];
4203 m_BCRegion.SetIndex(InputDimension-1,2);
4204 m_BC2Regions[1]=m_BCRegion;
4206 m_BCRegion.SetIndex(InputDimension-1,3);
4207 m_BC2Regions[2]=m_BCRegion;
4208 m_BC2Values[2]=-m_WeightRatio[3][1];
4215 m_FirstRegion.SetSize(InputDimension-1,2);
4216 m_FirstRegion.SetIndex(InputDimension-1,0);
4219 m_BCRegions.resize(3);
4220 m_BCValues.resize(3);
4221 m_BCRegion.SetIndex(InputDimension-1,1);
4222 m_BCRegions[0]=m_BCRegion;
4223 m_BCValues[0]=-m_WeightRatio[3][1];
4224 m_BCRegion.SetIndex(InputDimension-1,2);
4225 m_BCRegions[1]=m_BCRegion;
4227 m_BCRegion.SetIndex(InputDimension-1,3);
4228 m_BCRegions[2]=m_BCRegion;
4229 m_BCValues[2]=-m_WeightRatio[3][1];
4232 m_SecondRegion.SetSize(InputDimension-1,1);
4233 m_SecondRegion.SetIndex(InputDimension-1,2);
4236 m_BC2Regions.resize(0);
4243 m_FirstRegion.SetSize(InputDimension-1,1);
4244 m_FirstRegion.SetIndex(InputDimension-1,1);
4247 m_BCRegions.resize(3);
4248 m_BCValues.resize(3);
4249 m_BCRegion.SetIndex(InputDimension-1,1);
4250 m_BCRegions[0]=m_BCRegion;
4251 m_BCValues[0]=-m_WeightRatio[3][1];
4252 m_BCRegion.SetIndex(InputDimension-1,2);
4253 m_BCRegions[1]=m_BCRegion;
4255 m_BCRegion.SetIndex(InputDimension-1,3);
4256 m_BCRegions[2]=m_BCRegion;
4257 m_BCValues[2]=-m_WeightRatio[3][1];
4260 m_SecondRegion.SetSize(InputDimension-1,2);
4261 m_SecondRegion.SetIndex(InputDimension-1,2);
4264 m_BC2Regions.resize(0);
4271 m_FirstRegion.SetSize(InputDimension-1,0);
4274 m_BCRegions.resize(3);
4275 m_BCValues.resize(3);
4276 m_BCRegion.SetIndex(InputDimension-1,1);
4277 m_BCRegions[0]=m_BCRegion;
4278 m_BCValues[0]=-m_WeightRatio[3][1];
4279 m_BCRegion.SetIndex(InputDimension-1,2);
4280 m_BCRegions[1]=m_BCRegion;
4282 m_BCRegion.SetIndex(InputDimension-1,3);
4283 m_BCRegions[2]=m_BCRegion;
4284 m_BCValues[2]=-m_WeightRatio[3][1];
4287 m_SecondRegion.SetSize(InputDimension-1,1);
4288 m_SecondRegion.SetIndex(InputDimension-1,2);
4291 m_BC2Regions.resize(1);
4292 m_BC2Values.resize(1);
4293 m_BCRegion.SetIndex(InputDimension-1,0);
4294 m_BC2Regions[0]=m_BCRegion;
4302 m_FirstRegion.SetSize(InputDimension-1,2);
4303 m_FirstRegion.SetIndex(InputDimension-1,2);
4306 m_BCRegions.resize(1);
4307 m_BCValues.resize(1);
4308 m_BCRegion.SetIndex(InputDimension-1,0);
4309 m_BCRegions[0]=m_BCRegion;
4313 m_SecondRegion.SetSize(InputDimension-1,1);
4314 m_SecondRegion.SetIndex(InputDimension-1,4);
4317 m_BC2Regions.resize(0);
4323 DD("supportindex > 4 ???");
4324 DD(m_SupportRegion.GetIndex(InputDimension-1));
4326 } // end switch index
4329 } // end case 9 shape
4334 DD ("Other shapes currently not implemented");
4337 } // end switch shape
4338 } // end wrap region
4342 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
4344 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions,NOutputDimensions>
4345 ::TransformPointToContinuousIndex( const InputPointType & point, ContinuousIndexType & index ) const
4349 itk::Vector<double, OutputDimension> tvector;
4351 for ( j = 0; j < OutputDimension; j++ )
4353 //JV find the index in the PADDED version
4354 tvector[j] = point[j] - this->m_PaddedGridOrigin[j];
4357 itk::Vector<double, OutputDimension> cvector;
4359 cvector = m_PointToIndex * tvector;
4361 for ( j = 0; j < OutputDimension; j++ )
4363 index[j] = static_cast< typename ContinuousIndexType::CoordRepType >( cvector[j] );