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>::ShapedBLUTSpatioTemporalDeformableTransform():Superclass(0)
38 //JV default spline order
39 for ( i=0;i<InputDimension; i++)
42 //JV default LUTSamplingfactor
43 for ( i=0;i<InputDimension; i++)
44 m_LUTSamplingFactors[i]=20;
46 for (i=0;i<InputDimension; i++)
47 m_SupportSize[i] = m_SplineOrders[i]+1;
49 //Instantiate a Vector Bspline Interpolator
50 m_VectorInterpolator =VectorInterpolatorType::New();
51 m_VectorInterpolator->SetLUTSamplingFactors(m_LUTSamplingFactors);
52 m_VectorInterpolator->SetSplineOrders(m_SplineOrders);
54 // Set Bulk transform to NULL (no bulk performed)
55 m_BulkTransform = NULL;
63 // Default grid size is zero
67 //JV region containing the parameters
68 m_GridRegion.SetSize( m_NullSize);
69 m_GridRegion.SetIndex( m_NullIndex );
71 //JV use second region over the images
72 m_PaddedGridRegion.SetSize(m_NullSize);
73 m_PaddedGridRegion.SetIndex(m_NullIndex);
75 //JV Maintain two origins
76 m_GridOrigin.Fill( 0.0 ); // default origin is all zeros
77 m_PaddedGridOrigin.Fill( 0.0 ); // default origin is all zeros
79 m_GridSpacing.Fill( 1.0 ); // default spacing is all ones
80 m_GridDirection.SetIdentity(); // default spacing is all ones
82 m_InternalParametersBuffer = ParametersType(0);
83 // Make sure the parameters pointer is not NULL after construction.
84 m_InputParametersPointer = &m_InternalParametersBuffer;
86 // Initialize coeffient images
87 m_WrappedImage = CoefficientImageType::New();
88 m_WrappedImage->SetRegions( m_GridRegion );
89 m_WrappedImage->SetOrigin( m_GridOrigin.GetDataPointer() );
90 m_WrappedImage->SetSpacing( m_GridSpacing.GetDataPointer() );
91 m_WrappedImage->SetDirection( m_GridDirection );
93 // JV Initialize the padded version
94 m_PaddedCoefficientImage = CoefficientImageType::New();
95 m_PaddedCoefficientImage->SetRegions( m_PaddedGridRegion );
96 m_PaddedCoefficientImage->SetOrigin( m_GridOrigin.GetDataPointer() );
97 m_PaddedCoefficientImage->SetSpacing( m_GridSpacing.GetDataPointer() );
98 m_PaddedCoefficientImage->SetDirection( m_GridDirection );
100 m_CoefficientImage = NULL;
102 // Variables for computing interpolation
103 for (i=0; i <InputDimension;i++)
105 m_Offset[i] = m_SplineOrders[i] / 2;
106 if ( m_SplineOrders[i] % 2 )
108 m_SplineOrderOdd[i] = true;
112 m_SplineOrderOdd[i] = false;
116 //JV padded should be used when checking
117 m_ValidRegion = m_PaddedGridRegion;
119 // Initialize jacobian images
120 for (i=0; i <OutputDimension;i++)
122 m_JacobianImage[i] = JacobianImageType::New();
123 m_JacobianImage[i]->SetRegions( m_GridRegion );
124 m_JacobianImage[i]->SetOrigin( m_GridOrigin.GetDataPointer() );
125 m_JacobianImage[i]->SetSpacing( m_GridSpacing.GetDataPointer() );
126 m_JacobianImage[i]->SetDirection( m_GridDirection );
129 /** Fixed Parameters store the following information:
134 //JV we add the splineOrders, LUTsamplingfactor, m_Mask and m_BulkTransform
140 The size of these is equal to the NInputDimensions
141 *********************************************************/
142 this->m_FixedParameters.SetSize ( NInputDimensions * (NInputDimensions + 5)+3 );
143 this->m_FixedParameters.Fill ( 0.0 );
144 for ( i=0; i<NInputDimensions; i++)
146 this->m_FixedParameters[2*NInputDimensions+i] = m_GridSpacing[i];
148 for (unsigned int di=0; di<NInputDimensions; di++)
150 for (unsigned int dj=0; dj<NInputDimensions; dj++)
152 this->m_FixedParameters[3*NInputDimensions+(di*NInputDimensions+dj)] = m_GridDirection[di][dj];
156 //JV add splineOrders
157 for ( i=0; i<NInputDimensions; i++)
159 this->m_FixedParameters[ ( (3+NInputDimensions)*NInputDimensions)+i] = (this->GetSplineOrders())[i];
162 //JV add LUTsamplingFactors
163 for ( i=0; i<NInputDimensions; i++)
165 this->m_FixedParameters[( (4+NInputDimensions)*NInputDimensions)+i ] = this->GetLUTSamplingFactors()[i];
168 // JV add the mask pointer
169 this->m_FixedParameters[( (5+NInputDimensions)*NInputDimensions)]=(double)((size_t)m_Mask);
171 // JV add the bulkTransform pointer
172 this->m_FixedParameters[( (5+NInputDimensions)*NInputDimensions) +1]=(double)((size_t)m_BulkTransform);
174 // JV add the Transform shape
175 this->m_FixedParameters[( (5+NInputDimensions)*NInputDimensions) +2]=(double)(m_TransformShape);
178 // Calculate the PointToIndex matrices
180 for( unsigned int i=0; i<OutputDimension; i++)
182 scale[i][i] = m_GridSpacing[i];
185 m_IndexToPoint = m_GridDirection * scale;
186 m_PointToIndex = m_IndexToPoint.GetInverse();
188 m_LastJacobianIndex = m_ValidRegion.GetIndex();
191 //Weights to be used for the BC
193 std::vector<double> weights(5);
202 m_Weights[0]=weights;
210 m_Weights[1]=weights;
218 m_Weights[2]=weights;
226 m_Weights[3]=weights;
228 // Update the WeightRatios
229 m_WeightRatio.resize(4);
230 for (unsigned int i=0; i<4; i++)
232 m_WeightRatio[i].resize(4);
233 for (unsigned int j=0; j<4; j++)
234 m_WeightRatio[i][j]=m_Weights[m_SplineOrders[InputDimension-1]][i]/ m_Weights[m_SplineOrders[InputDimension-1]][j];
237 //JV initialize some variables for jacobian calculation
238 m_SupportRegion.SetSize(m_SupportSize);
239 m_SupportIndex.Fill(0);
240 m_SupportRegion.SetIndex(m_SupportIndex);
241 for ( i = 0; i < SpaceDimension ; i++ )
242 m_ZeroVector[i]=itk::NumericTraits<JacobianValueType>::Zero;
245 m_FirstRegion.SetSize(m_NullSize);
246 m_FirstRegion.SetIndex(m_NullIndex);
247 m_SecondRegion.SetSize(m_NullSize);
248 m_SecondRegion.SetIndex(m_NullIndex);
249 m_ThirdRegion.SetSize(m_NullSize);
250 m_ThirdRegion.SetIndex(m_NullIndex);
252 m_BCValues.resize(0);
253 m_BCRegions.resize(0);
255 m_BC2Values.resize(0);
256 m_BC2Regions.resize(0);
258 m_BC3Values.resize(0);
259 m_BC3Regions.resize(0);
263 for ( i = 0; i < SpaceDimension ; i++ )
265 m_FirstIterator[i]= IteratorType( m_JacobianImage[i], m_FirstRegion);
266 m_SecondIterator[i]= IteratorType( m_JacobianImage[i], m_SecondRegion);
267 m_ThirdIterator[i]= IteratorType( m_JacobianImage[i], m_ThirdRegion);
268 m_BCIterators[i].resize(0);
269 m_BC2Iterators[i].resize(0);
270 m_BC3Iterators[i].resize(0);
279 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
280 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
281 ::~ShapedBLUTSpatioTemporalDeformableTransform()
287 // JV set Spline Order
288 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
290 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
291 ::SetSplineOrder(const unsigned int & splineOrder)
293 SizeType splineOrders;
294 for (unsigned int i=0;i<InputDimension;i++)splineOrders[i]=splineOrder;
296 this->SetSplineOrders(splineOrders);
300 // JV set Spline Orders
301 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
303 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
304 ::SetSplineOrders(const SizeType & splineOrders)
306 if(m_SplineOrders!=splineOrders)
308 m_SplineOrders=splineOrders;
310 //update the interpolation function
311 m_VectorInterpolator->SetSplineOrders(m_SplineOrders);
313 //update the varaibles for computing interpolation
314 for (unsigned int i=0; i <InputDimension;i++)
316 m_SupportSize[i] = m_SplineOrders[i]+1;
317 m_Offset[i] = m_SplineOrders[i] / 2;
319 if ( m_SplineOrders[i] % 2 )
321 m_SplineOrderOdd[i] = true;
325 m_SplineOrderOdd[i] = false;
329 //SupportSize is updated!, update regions
330 //JV initialize some variables for jacobian calculation
331 m_SupportRegion.SetSize(m_SupportSize);
332 m_SupportIndex.Fill(0);
333 m_SupportRegion.SetIndex(m_SupportIndex);
335 // Update the WeightRatios
336 for (unsigned int i=0; i<4; i++)
337 for (unsigned int j=0; j<4; j++)
338 m_WeightRatio[i][j]=m_Weights[m_SplineOrders[InputDimension-1]][i]/ m_Weights[m_SplineOrders[InputDimension-1]][j];
345 // JV set sampling factor
346 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
348 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
349 ::SetLUTSamplingFactor( const int & samplingFactor)
351 SizeType samplingFactors;
352 for (unsigned int i=0; i<NInputDimensions; i++)
353 samplingFactors[i]=samplingFactor;
355 //update the interpolation function
356 this->SetLUTSamplingFactors(samplingFactors);
360 // JV set sampling factors
361 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
363 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
364 ::SetLUTSamplingFactors( const SizeType & samplingFactors)
366 if(m_LUTSamplingFactors!=samplingFactors)
368 for (unsigned int i=0; i<NInputDimensions; i++)
369 m_LUTSamplingFactors[i]=samplingFactors[i];
371 //update the interpolation function
372 m_VectorInterpolator->SetLUTSamplingFactors(m_LUTSamplingFactors);
379 // Get the number of parameters
380 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
381 typename ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>::NumberOfParametersType
382 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
383 ::GetNumberOfParameters(void) const
386 // The number of parameters equal SpaceDimension * number of
387 // of pixels in the grid region.
388 return ( static_cast<unsigned int>( SpaceDimension ) *
389 static_cast<unsigned int>( m_GridRegion.GetNumberOfPixels() ) );
394 // Get the padded number of parameters
395 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
397 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
398 ::GetPaddedNumberOfParameters(void) const
401 // The number of parameters equal SpaceDimension * number of
402 // of pixels in the grid region.
403 return ( static_cast<unsigned int>( SpaceDimension ) *
404 static_cast<unsigned int>( m_PaddedGridRegion.GetNumberOfPixels() ) );
410 // Get the number of parameters per dimension
411 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
413 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
414 ::GetNumberOfParametersPerDimension(void) const
416 // The number of parameters per dimension equal number of
417 // of pixels in the grid region.
418 return ( static_cast<unsigned int>( m_GridRegion.GetNumberOfPixels() ) );
423 // Set the grid region
424 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
426 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
427 ::SetGridRegion( const RegionType & region )
429 if ( m_GridRegion != region )
431 m_GridRegion = region;
432 m_PaddedGridRegion=region;
434 // JV set the padded region
435 typename CoefficientImageType::RegionType::SizeType paddedSize= region.GetSize();
437 //JV size dependes on shape
438 switch (m_TransformShape)
442 paddedSize[InputDimension-1]+=4;
446 paddedSize[InputDimension-1]+=4;
450 paddedSize[InputDimension-1]+=3;
454 paddedSize[InputDimension-1]+=2;
458 paddedSize[InputDimension-1]+=3;
461 paddedSize[InputDimension-1]=+1;
463 m_PaddedGridRegion.SetSize(paddedSize);
465 // Set regions for each coefficient and jacobian image
466 m_WrappedImage->SetRegions( m_GridRegion );
467 m_PaddedCoefficientImage->SetRegions( m_PaddedGridRegion );
468 m_PaddedCoefficientImage->Allocate();
469 for (unsigned int j=0; j <OutputDimension;j++)
471 m_JacobianImage[j]->SetRegions( m_GridRegion );
474 // JV used the padded version for the valid region
475 // Set the valid region
476 // If the grid spans the interval [start, last].
477 // The valid interval for evaluation is [start+offset, last-offset]
478 // when spline order is even.
479 // The valid interval for evaluation is [start+offset, last-offset)
480 // when spline order is odd.
481 // Where offset = vcl_floor(spline / 2 ).
482 // Note that the last pixel is not included in the valid region
483 // with odd spline orders.
484 typename RegionType::SizeType size = m_PaddedGridRegion.GetSize();
485 typename RegionType::IndexType index = m_PaddedGridRegion.GetIndex();
486 for ( unsigned int j = 0; j < NInputDimensions; j++ )
489 static_cast< typename RegionType::IndexValueType >( m_Offset[j] );
491 static_cast< typename RegionType::SizeValueType> ( 2 * m_Offset[j] );
492 m_ValidRegionLast[j] = index[j] +
493 static_cast< typename RegionType::IndexValueType >( size[j] ) - 1;
495 m_ValidRegion.SetSize( size );
496 m_ValidRegion.SetIndex( index );
498 // If we are using the default parameters, update their size and set to identity.
499 // Input parameters point to internal buffer => using default parameters.
500 if (m_InputParametersPointer == &m_InternalParametersBuffer)
502 // Check if we need to resize the default parameter buffer.
503 if ( m_InternalParametersBuffer.GetSize() != this->GetNumberOfParameters() )
505 m_InternalParametersBuffer.SetSize( this->GetNumberOfParameters() );
506 // Fill with zeros for identity.
507 m_InternalParametersBuffer.Fill( 0 );
516 // Set the grid spacing
517 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
519 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
520 ::SetGridSpacing( const SpacingType & spacing )
522 if ( m_GridSpacing != spacing )
524 m_GridSpacing = spacing;
526 // Set spacing for each coefficient and jacobian image
527 m_WrappedImage->SetSpacing( m_GridSpacing.GetDataPointer() );
528 m_PaddedCoefficientImage->SetSpacing( m_GridSpacing.GetDataPointer() );
529 for (unsigned int j=0; j <OutputDimension; j++) m_JacobianImage[j]->SetSpacing( m_GridSpacing.GetDataPointer() );
533 for( unsigned int i=0; i<OutputDimension; i++)
535 scale[i][i] = m_GridSpacing[i];
538 m_IndexToPoint = m_GridDirection * scale;
539 m_PointToIndex = m_IndexToPoint.GetInverse();
546 // Set the grid direction
547 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
549 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
550 ::SetGridDirection( const DirectionType & direction )
552 if ( m_GridDirection != direction )
554 m_GridDirection = direction;
556 // Set direction for each coefficient and jacobian image
557 m_WrappedImage->SetDirection( m_GridDirection );
558 m_PaddedCoefficientImage->SetDirection( m_GridDirection );
559 for (unsigned int j=0; j <OutputDimension; j++) m_JacobianImage[j]->SetDirection( m_GridDirection );
563 for( unsigned int i=0; i<OutputDimension; i++)
565 scale[i][i] = m_GridSpacing[i];
568 m_IndexToPoint = m_GridDirection * scale;
569 m_PointToIndex = m_IndexToPoint.GetInverse();
577 // Set the grid origin
578 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
580 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
581 ::SetGridOrigin( const OriginType& origin )
583 if( m_GridOrigin!=origin)
585 m_GridOrigin = origin;
587 // JV The origin depends on the shape
588 switch (m_TransformShape)
592 m_PaddedGridOrigin=origin;
593 m_PaddedGridOrigin[InputDimension-1]=origin[InputDimension-1]-2* m_GridSpacing[InputDimension-1];
597 m_PaddedGridOrigin=origin;
598 m_PaddedGridOrigin[InputDimension-1]=origin[InputDimension-1]-1* 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;
611 m_PaddedGridOrigin=origin;
612 m_PaddedGridOrigin[InputDimension-1]=origin[InputDimension-1]-1* m_GridSpacing[InputDimension-1];
615 m_PaddedGridOrigin=origin;
618 // Set origin for each coefficient and jacobianimage
619 m_WrappedImage->SetOrigin( m_GridOrigin.GetDataPointer() );
620 m_PaddedCoefficientImage->SetOrigin( m_PaddedGridOrigin.GetDataPointer() );
621 for (unsigned int j=0; j <OutputDimension; j++) m_JacobianImage[j]->SetOrigin( m_GridOrigin.GetDataPointer() );
628 // Set the parameters
629 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
631 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
634 if( m_InputParametersPointer )
636 ParametersType * parameters =
637 const_cast<ParametersType *>( m_InputParametersPointer );
638 parameters->Fill( 0.0 );
643 itkExceptionMacro( << "Input parameters for the spline haven't been set ! "
644 << "Set them using the SetParameters or SetCoefficientImage method first." );
649 // Set the parameters
650 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
652 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
653 ::SetParameters( const ParametersType & parameters )
656 // Check if the number of parameters match the
657 // Expected number of parameters
658 if ( parameters.Size() != this->GetNumberOfParameters() )
660 itkExceptionMacro(<<"Mismatched between parameters size "
662 << " and region size "
663 << m_GridRegion.GetNumberOfPixels() );
666 // Clean up buffered parameters
667 m_InternalParametersBuffer = ParametersType( 0 );
669 // Keep a reference to the input parameters
670 m_InputParametersPointer = ¶meters;
672 // Wrap flat array as images of coefficients
673 this->WrapAsImages();
675 //JV Set padded input to vector interpolator
676 m_VectorInterpolator->SetInputImage(this->GetPaddedCoefficientImage());
678 // Modified is always called since we just have a pointer to the
679 // parameters and cannot know if the parameters have changed.
684 // Set the Fixed Parameters
685 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
687 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
688 ::SetFixedParameters( const ParametersType & parameters )
691 // JV number should be exact, no defaults for spacing
692 if ( parameters.Size() != NInputDimensions * (5 + NInputDimensions)+3 )
694 itkExceptionMacro(<< "Mismatched between parameters size "
696 << " and number of fixed parameters "
697 << NInputDimensions * (5 + NInputDimensions)+3 );
699 /*********************************************************
700 Fixed Parameters store the following information:
705 // JV we add the splineOrders, LUTsamplingfactor, mask pointer and bulktransform pointer
713 The size of these is equal to the NInputDimensions
714 *********************************************************/
716 /** Set the Grid Parameters */
718 for (unsigned int i=0; i<NInputDimensions; i++)
720 gridSize[i] = static_cast<int> (parameters[i]);
722 RegionType bsplineRegion;
723 bsplineRegion.SetSize( gridSize );
725 /** Set the Origin Parameters */
727 for (unsigned int i=0; i<NInputDimensions; i++)
729 origin[i] = parameters[NInputDimensions+i];
732 /** Set the Spacing Parameters */
734 for (unsigned int i=0; i<NInputDimensions; i++)
736 spacing[i] = parameters[2*NInputDimensions+i];
739 /** Set the Direction Parameters */
740 DirectionType direction;
741 for (unsigned int di=0; di<NInputDimensions; di++)
743 for (unsigned int dj=0; dj<NInputDimensions; dj++)
745 direction[di][dj] = parameters[3*NInputDimensions+(di*NInputDimensions+dj)];
749 //JV add the spline orders
750 SizeType splineOrders;
751 for (unsigned int i=0; i<NInputDimensions; i++)
753 splineOrders[i]=parameters[(3+NInputDimensions)*NInputDimensions+i];
756 //JV add the LUT sampling factor
757 SizeType samplingFactors;
758 for (unsigned int i=0; i<NInputDimensions; i++)
760 samplingFactors[i]=parameters[(4+NInputDimensions)*NInputDimensions+i];
763 //JV add the MaskPointer
764 m_Mask=((MaskPointer)( (size_t)parameters[(5+NInputDimensions)*NInputDimensions]));
766 //JV add the MaskPointer
767 m_BulkTransform=((BulkTransformPointer)( (size_t)parameters[(5+NInputDimensions)*NInputDimensions+1]));
769 //JV add the TransformShape
770 m_TransformShape=((unsigned int)parameters[(5+NInputDimensions)*NInputDimensions+2]);
773 this->SetSplineOrders( splineOrders );
774 this->SetGridSpacing( spacing );
775 this->SetGridDirection( direction );
776 this->SetGridOrigin( origin );
777 this->SetGridRegion( bsplineRegion );
778 this->SetLUTSamplingFactors( samplingFactors );
783 // Wrap flat parameters as images
784 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
786 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
789 //JV Wrap parameter array in vectorial image, changed parameter order: A1x A1y A1z, A2x ....
790 PixelType * dataPointer =reinterpret_cast<PixelType *>( const_cast<double *>(m_InputParametersPointer->data_block() )) ;
791 unsigned int numberOfPixels = m_GridRegion.GetNumberOfPixels();
793 m_WrappedImage->GetPixelContainer()->SetImportPointer( dataPointer,numberOfPixels);//InputDimension
794 m_CoefficientImage = m_WrappedImage;
796 //=====================================
797 //JV Create padded structure adding BC
798 //=====================================
799 PadCoefficientImage();
801 //=====================================
802 //JV Wrap jacobian into OutputDimension X Vectorial images
803 //=====================================
804 this->m_SharedDataBSplineJacobian.set_size( OutputDimension, this->GetNumberOfParameters() );
806 // Use memset to set the memory
807 // JV four rows of three comps of parameters
808 JacobianPixelType * jacobianDataPointer = reinterpret_cast<JacobianPixelType *>(this->m_SharedDataBSplineJacobian.data_block());
809 memset(jacobianDataPointer, 0, OutputDimension*numberOfPixels*sizeof(JacobianPixelType));
811 for (unsigned int j=0; j<OutputDimension; j++)
813 m_JacobianImage[j]->GetPixelContainer()->
814 SetImportPointer( jacobianDataPointer, numberOfPixels );
815 jacobianDataPointer += numberOfPixels;
818 // Reset the J parameters
819 m_LastJacobianIndex = m_ValidRegion.GetIndex();
820 m_FirstRegion.SetSize(m_NullSize);
821 m_SecondRegion.SetSize(m_NullSize);
822 for ( unsigned int j = 0; j < SpaceDimension ; j++ )
824 m_FirstIterator[j]= IteratorType( m_JacobianImage[j], m_FirstRegion);
825 m_SecondIterator[j]= IteratorType( m_JacobianImage[j], m_SecondRegion);
828 m_BCValues.resize(0);
829 m_BCRegions.resize(0);
831 m_BC2Values.resize(0);
832 m_BC2Regions.resize(0);
834 m_BC3Values.resize(0);
835 m_BC3Regions.resize(0);
841 // Set the parameters by value
842 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
844 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
845 ::SetParametersByValue( const ParametersType & parameters )
848 // Check if the number of parameters match the
849 // Expected number of parameters
850 if ( parameters.Size() != this->GetNumberOfParameters() )
852 itkExceptionMacro(<<"Mismatched between parameters size "
854 << " and region size "
855 << m_GridRegion.GetNumberOfPixels() );
859 m_InternalParametersBuffer = parameters;
860 m_InputParametersPointer = &m_InternalParametersBuffer;
862 // Wrap flat array as images of coefficients
863 this->WrapAsImages();
865 //JV Set padded input to vector interpolator
866 m_VectorInterpolator->SetInputImage(this->GetPaddedCoefficientImage());
868 // Modified is always called since we just have a pointer to the
869 // Parameters and cannot know if the parameters have changed.
874 // Get the parameters
875 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
877 typename ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
879 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
880 ::GetParameters( void ) const
882 /** NOTE: For efficiency, this class does not keep a copy of the parameters -
883 * it just keeps pointer to input parameters.
885 if (NULL == m_InputParametersPointer)
887 itkExceptionMacro( <<"Cannot GetParameters() because m_InputParametersPointer is NULL. Perhaps SetCoefficientImage() has been called causing the NULL pointer." );
890 return (*m_InputParametersPointer);
894 // Get the parameters
895 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
897 typename ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
899 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
900 ::GetFixedParameters( void ) const
902 RegionType resRegion = this->GetGridRegion( );
904 for (unsigned int i=0; i<NInputDimensions; i++)
906 this->m_FixedParameters[i] = (resRegion.GetSize())[i];
908 for (unsigned int i=0; i<NInputDimensions; i++)
910 this->m_FixedParameters[NInputDimensions+i] = (this->GetGridOrigin())[i];
912 for (unsigned int i=0; i<NInputDimensions; i++)
914 this->m_FixedParameters[2*NInputDimensions+i] = (this->GetGridSpacing())[i];
916 for (unsigned int di=0; di<NInputDimensions; di++)
918 for (unsigned int dj=0; dj<NInputDimensions; dj++)
920 this->m_FixedParameters[3*NInputDimensions+(di*NInputDimensions+dj)] = (this->GetGridDirection())[di][dj];
924 //JV add splineOrders
925 for (unsigned int i=0; i<NInputDimensions; i++)
927 this->m_FixedParameters[(3+NInputDimensions)*NInputDimensions+i] = (this->GetSplineOrders())[i];
930 //JV add LUTsamplingFactor
931 for (unsigned int i=0; i<NInputDimensions; i++)
933 this->m_FixedParameters[(4+NInputDimensions)*NInputDimensions+i] = (this->GetLUTSamplingFactors())[i];
937 this->m_FixedParameters[(5+NInputDimensions)*NInputDimensions]=(double)((size_t) m_Mask);
939 //JV add the bulktransform pointer
940 this->m_FixedParameters[(5+NInputDimensions)*NInputDimensions+1]=(double)((size_t) m_BulkTransform);
942 //JV add the transform shape
943 this->m_FixedParameters[(5+NInputDimensions)*NInputDimensions+2]=(double)( m_TransformShape);
945 return (this->m_FixedParameters);
949 // Set the B-Spline coefficients using input images
950 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
952 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
953 ::SetCoefficientImage( CoefficientImagePointer image )
955 this->SetGridSpacing( image->GetSpacing() );
956 this->SetGridOrigin( image->GetOrigin() );
957 this->SetGridDirection( image->GetDirection() );
958 this->SetGridRegion( image->GetBufferedRegion() );
959 m_CoefficientImage = image;
962 // m_WrappedImage=m_CoefficientImage;
964 // Update the interpolator
965 this->PadCoefficientImage();
966 m_VectorInterpolator->SetInputImage(this->GetPaddedCoefficientImage());
968 // Clean up buffered parameters
969 m_InternalParametersBuffer = ParametersType( 0 );
970 m_InputParametersPointer = NULL;
975 // // Set the B-Spline coefficients using input images
976 // template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
978 // ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
979 // ::SetPaddedCoefficientImage( CoefficientImagePointer image )
981 // //JV modify the region
982 // typename CoefficientImageType::RegionType region=image->GetBufferedRegion();
983 // typename CoefficientImageType::RegionType::SizeType size=region.GetSize();
984 // size[InputDimension-1]-=2;
985 // region.SetSize(size);
988 // this->SetGridRegion( region );
989 // this->SetGridSpacing( image->GetSpacing() );
990 // this->SetGridDirection( image->GetDirection() );
991 // this->SetGridOrigin( image->GetOrigin() );
992 // m_PaddedCoefficientImage = image;
993 // this->ExtractCoefficientImage();
994 // m_VectorInterpolator->SetInputImage(this->GetPaddedCoefficientImage());
996 // // Clean up buffered parameters
997 // m_InternalParametersBuffer = ParametersType( 0 );
998 // m_InputParametersPointer = NULL;
1002 // Set the B-Spline coefficients using input images
1003 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
1004 typename ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>::CoefficientImageType::Pointer
1005 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
1006 ::ExtractTemporalRow(const typename CoefficientImageType::Pointer& coefficientImage, unsigned int temporalIndex)
1009 typename CoefficientImageType::RegionType sourceRegion=coefficientImage->GetLargestPossibleRegion();
1010 sourceRegion.SetSize(InputDimension-1, 1);
1011 sourceRegion.SetIndex(InputDimension-1, temporalIndex);
1014 typedef clitk::ExtractImageFilter<CoefficientImageType, CoefficientImageType> ExtractImageFilterType;
1015 typename ExtractImageFilterType::Pointer extract=ExtractImageFilterType::New();
1016 extract->SetInput(coefficientImage);
1017 extract->SetExtractionRegion(sourceRegion);
1019 typename CoefficientImageType::Pointer row= extract->GetOutput();
1021 // Set index to zero
1022 sourceRegion.SetIndex(InputDimension-1, 0);
1023 row->SetRegions(sourceRegion);
1028 // Set the B-Spline coefficients using input images
1029 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
1031 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
1032 ::PadCoefficientImage(void)
1035 // Define paste, extract and combine filters
1036 typedef itk::PasteImageFilter<CoefficientImageType, CoefficientImageType, CoefficientImageType> PasteImageFilterType;
1037 typedef clitk::ExtractImageFilter<CoefficientImageType, CoefficientImageType> ExtractImageFilterType;
1038 typedef clitk::LinearCombinationImageFilter<CoefficientImageType, CoefficientImageType> LinearCombinationFilterType;
1041 typename CoefficientImageType::RegionType sourceRegion=m_PaddedCoefficientImage->GetLargestPossibleRegion();
1042 typename CoefficientImageType::RegionType destinationRegion=m_PaddedCoefficientImage->GetLargestPossibleRegion();
1043 typename CoefficientImageType::RegionType::SizeType sourceSize=sourceRegion.GetSize();
1044 typename CoefficientImageType::IndexType sourceIndex=sourceRegion.GetIndex();
1045 typename CoefficientImageType::IndexType destinationIndex=destinationRegion.GetIndex();
1047 // JV Padding depends on the shape
1048 switch (m_TransformShape)
1053 2: rabbit 4 CP 3 DOF
1054 3: rabbit 5 CP 4 DOF
1055 4: sputnik 4 CP 4 DOF
1056 5: sputnik 5 CP 5 DOF
1057 6: diamond 6 CP 5 DOF
1058 7: diamond 7 CP 6 DOF
1063 // ----------------------------------------------------------------------
1064 // The egg with 4 internal CP (starting from inhale)
1065 // Periodic, constrained to zero at the reference
1066 // at position 3 and
1067 // Coeff row BC1 BC2 0 BC3 1 2 BC4
1068 // PaddedCoeff R: 0 1 2 3 4 5 6
1071 // BC3= -weights[2]/weights[0] ( R2+R4 )
1073 // ---------------------------------------------------------------------
1075 //---------------------------------
1076 // 1. First two temporal rows are identical: paste 1-2 to 0-1
1077 typename PasteImageFilterType::Pointer paster0=PasteImageFilterType::New();
1078 paster0->SetDestinationImage(m_PaddedCoefficientImage);
1079 paster0->SetSourceImage(m_CoefficientImage);
1080 sourceRegion.SetIndex(NInputDimensions-1,1);
1081 sourceRegion.SetSize(NInputDimensions-1,2);
1082 destinationIndex[InputDimension-1]=0;
1083 paster0->SetDestinationIndex(destinationIndex);
1084 paster0->SetSourceRegion(sourceRegion);
1086 //---------------------------------
1087 // 1. Next temporal row = paste 0 to 2
1088 typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1089 typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1090 paster1->SetDestinationImage(paster0->GetOutput());
1091 paster1->SetSourceImage(row0);
1092 destinationIndex[InputDimension-1]=2;
1093 paster1->SetDestinationIndex(destinationIndex);
1094 paster1->SetSourceRegion(row0->GetLargestPossibleRegion());
1096 //---------------------------------
1097 // 2. Middle row at index 3=BC
1098 // Extract row at index 1, 2 (of coeff image)
1099 typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1102 typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1103 combine1->SetInput(0,row0);
1104 combine1->SetInput(1,row1);
1105 combine1->SetA(-m_WeightRatio[2][0]);
1106 combine1->SetB(-m_WeightRatio[2][0]);
1108 typename CoefficientImageType::Pointer bc3Row=combine1->GetOutput();
1110 // Paste middleRow at index 3 (padded coeff)
1111 typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1112 paster2->SetDestinationImage(paster1->GetOutput());
1113 paster2->SetSourceImage(bc3Row);
1114 destinationIndex[InputDimension-1]=3;
1115 paster2->SetDestinationIndex(destinationIndex);
1116 paster2->SetSourceRegion(bc3Row->GetLargestPossibleRegion());
1118 //---------------------------------
1119 // 3. Next two temporal rows identical: paste 1,2 to 4,5
1120 typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1121 paster3->SetDestinationImage(paster2->GetOutput());
1122 paster3->SetSourceImage(m_CoefficientImage);
1123 sourceRegion.SetIndex(InputDimension-1,1);
1124 sourceRegion.SetSize(InputDimension-1,2);
1125 destinationIndex[InputDimension-1]=4;
1126 paster3->SetDestinationIndex(destinationIndex);
1127 paster3->SetSourceRegion(sourceRegion);
1129 // ---------------------------------
1130 // 4. Row at index 6=BC4= R2
1131 // Paste BC3 row at index 5
1132 typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
1133 paster4->SetDestinationImage(paster3->GetOutput());
1134 paster4->SetSourceImage(row0);
1135 destinationIndex[InputDimension-1]=6;
1136 paster4->SetDestinationIndex(destinationIndex);
1137 paster4->SetSourceRegion(row0->GetLargestPossibleRegion());
1139 // Update the chain!
1141 m_PaddedCoefficientImage= paster4->GetOutput();
1148 // ----------------------------------------------------------------------
1149 // The egg with 5 internal CP (starting from inhale)
1150 // Periodic, constrained to zero at the reference
1151 // at position 3 and
1152 // Coeff row BC1 BC2 0 BC3 1 2 3 BC4
1153 // PaddedCoeff R: 0 1 2 3 4 5 6 7
1156 // BC3= -weights[3]/weights[1] ( R2+R5 ) - R4
1158 // ---------------------------------------------------------------------
1159 //---------------------------------
1160 // 1. First two temporal rows are identical: paste 2-3 to 0-1
1161 typename PasteImageFilterType::Pointer paster0=PasteImageFilterType::New();
1162 paster0->SetDestinationImage(m_PaddedCoefficientImage);
1163 paster0->SetSourceImage(m_CoefficientImage);
1164 sourceRegion.SetIndex(NInputDimensions-1,2);
1165 sourceRegion.SetSize(NInputDimensions-1,2);
1166 destinationIndex[InputDimension-1]=0;
1167 paster0->SetDestinationIndex(destinationIndex);
1168 paster0->SetSourceRegion(sourceRegion);
1170 //---------------------------------
1171 // 1. Next temporal row = paste 0 to 2
1172 typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1173 typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1174 paster1->SetDestinationImage(paster0->GetOutput());
1175 paster1->SetSourceImage(row0);
1176 destinationIndex[InputDimension-1]=2;
1177 paster1->SetDestinationIndex(destinationIndex);
1178 paster1->SetSourceRegion(row0->GetLargestPossibleRegion());
1180 //---------------------------------
1181 // 2. Middle row at index 3=BC
1182 // Extract row at index 1, 2 (of coeff image)
1183 typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1184 typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
1187 typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1188 combine1->SetInput(0,row0);
1189 combine1->SetInput(1,row2);
1192 // combine1->Update();
1193 typename LinearCombinationFilterType::Pointer combine2=LinearCombinationFilterType::New();
1194 combine2->SetInput(0,row1);
1195 combine2->SetInput(1,combine1->GetOutput());
1196 combine2->SetA(-1.);
1197 combine2->SetB(-m_WeightRatio[3][1]);
1199 typename CoefficientImageType::Pointer bc3Row=combine2->GetOutput();
1201 // Paste middleRow at index 3 (padded coeff)
1202 typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1203 paster2->SetDestinationImage(paster1->GetOutput());
1204 paster2->SetSourceImage(bc3Row);
1205 destinationIndex[InputDimension-1]=3;
1206 paster2->SetDestinationIndex(destinationIndex);
1207 paster2->SetSourceRegion(bc3Row->GetLargestPossibleRegion());
1209 //---------------------------------
1210 // 3. Next three temporal rows identical: paste 1,2,3 to 4,5,6
1211 typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1212 paster3->SetDestinationImage(paster2->GetOutput());
1213 paster3->SetSourceImage(m_CoefficientImage);
1214 sourceRegion.SetIndex(InputDimension-1,1);
1215 sourceRegion.SetSize(InputDimension-1,3);
1216 destinationIndex[InputDimension-1]=4;
1217 paster3->SetDestinationIndex(destinationIndex);
1218 paster3->SetSourceRegion(sourceRegion);
1220 // ---------------------------------
1221 // 4. Row at index 7=BC4= R2
1222 // Paste BC3 row at index 5
1223 typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
1224 paster4->SetDestinationImage(paster3->GetOutput());
1225 paster4->SetSourceImage(row0);
1226 destinationIndex[InputDimension-1]=7;
1227 paster4->SetDestinationIndex(destinationIndex);
1228 paster4->SetSourceRegion(row0->GetLargestPossibleRegion());
1230 // Update the chain!
1232 m_PaddedCoefficientImage= paster4->GetOutput();
1238 // // ----------------------------------------------------------------------
1239 // // The egg with 5 internal CP:
1240 // // Periodic, C2 smooth everywhere and constrained to zero at the reference
1241 // // Coeff row R5 BC1 0 1 2 3 BC2 R2
1242 // // PaddedCoeff R: 0 1 2 3 4 5 6 7
1243 // // BC1= -weights[2]/weights[0] ( R2+R5)
1245 // // ---------------------------------------------------------------------
1247 // // Extract rows with index 0 and 3
1248 // typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1249 // typename CoefficientImageType::Pointer row3=this->ExtractTemporalRow(m_CoefficientImage, 3);
1251 // // Paste the first row
1252 // typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1253 // destinationIndex.Fill(0);
1254 // paster1->SetDestinationIndex(destinationIndex);
1255 // paster1->SetSourceRegion(row3->GetLargestPossibleRegion());
1256 // paster1->SetSourceImage(row3);
1257 // paster1->SetDestinationImage(m_PaddedCoefficientImage);
1259 // // Linearly Combine rows for BC1 and BC2
1260 // typedef clitk::LinearCombinationImageFilter<CoefficientImageType, CoefficientImageType> LinearCombinationFilterType;
1261 // typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1262 // combine1->SetFirstInput(row0);
1263 // combine1->SetSecondInput(row3);
1264 // combine1->SetA(-m_WeightRatio[2][0]);
1265 // combine1->SetB(-m_WeightRatio[2][0]);
1266 // combine1->Update();
1267 // typename CoefficientImageType::Pointer bcRow=combine1->GetOutput();
1269 // // Paste the second row
1270 // typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1271 // destinationIndex[InputDimension-1]=1;
1272 // paster2->SetDestinationIndex(destinationIndex);
1273 // paster2->SetSourceRegion(bcRow->GetLargestPossibleRegion());
1274 // paster2->SetSourceImage(bcRow);
1275 // paster2->SetDestinationImage(paster1->GetOutput());
1277 // // Paste the coefficientImage
1278 // typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1279 // destinationIndex[InputDimension-1]=2;
1280 // paster3->SetDestinationIndex(destinationIndex);
1281 // paster3->SetSourceRegion(m_CoefficientImage->GetLargestPossibleRegion());
1282 // paster3->SetSourceImage(m_CoefficientImage);
1283 // paster3->SetDestinationImage(paster2->GetOutput());
1285 // // Paste the last two rows
1286 // typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
1287 // destinationIndex[InputDimension-1]=6;
1288 // paster4->SetDestinationIndex(destinationIndex);
1289 // paster4->SetSourceRegion(bcRow->GetLargestPossibleRegion());
1290 // paster4->SetSourceImage(bcRow);
1291 // paster4->SetDestinationImage(paster3->GetOutput());
1293 // typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
1294 // destinationIndex[InputDimension-1]=7;
1295 // paster5->SetDestinationIndex(destinationIndex);
1296 // paster5->SetSourceRegion(row0->GetLargestPossibleRegion());
1297 // paster5->SetSourceImage(row0);
1298 // paster5->SetDestinationImage(paster4->GetOutput());
1300 // // Update the chain!
1301 // paster5->Update();
1302 // m_PaddedCoefficientImage= paster5->GetOutput();
1309 // ----------------------------------------------------------------------
1310 // The rabbit with 4 internal CP
1311 // Periodic, constrained to zero at the reference
1312 // at position 3 and the extremes fixed with anit-symm bc
1313 // Coeff row BC1 0 1 BC2 2 BC3 BC4
1314 // PaddedCoeff R: 0 1 2 3 4 5 6
1316 // BC2= -weights[2]/weights[0] ( R2+R4 )
1319 // ---------------------------------------------------------------------
1321 // ---------------------------------
1322 // 0. First Row =BC1
1323 typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1324 typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1325 typename LinearCombinationFilterType::Pointer combine0=LinearCombinationFilterType::New();
1326 combine0->SetInput(0,row0);
1327 combine0->SetInput(1,row1);
1329 combine0->SetB(-1.);
1331 typename CoefficientImageType::Pointer bc1Row=combine0->GetOutput();
1333 typename PasteImageFilterType::Pointer paster0=PasteImageFilterType::New();
1334 paster0->SetDestinationImage(m_PaddedCoefficientImage);
1335 paster0->SetSourceImage(bc1Row);
1336 destinationIndex[InputDimension-1]=0;
1337 paster0->SetDestinationIndex(destinationIndex);
1338 paster0->SetSourceRegion(bc1Row->GetLargestPossibleRegion());
1340 //---------------------------------
1341 // 1. Next two temporal rows are identical: paste 0-1 to 1-2
1342 typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1343 paster1->SetDestinationImage(paster0->GetOutput());
1344 paster1->SetSourceImage(m_CoefficientImage);
1345 sourceRegion.SetIndex(NInputDimensions-1,0);
1346 destinationIndex[InputDimension-1]=1;
1347 sourceRegion.SetSize(NInputDimensions-1,2);
1348 paster1->SetDestinationIndex(destinationIndex);
1349 paster1->SetSourceRegion(sourceRegion);
1351 //---------------------------------
1352 // 2. Middle row at index 3=BC
1353 // Extract row at index 1, 2 (of coeff image)
1354 typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
1357 typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1358 combine1->SetInput(0,row1);
1359 combine1->SetInput(1,row2);
1360 combine1->SetA(-m_WeightRatio[2][0]);
1361 combine1->SetB(-m_WeightRatio[2][0]);
1363 typename CoefficientImageType::Pointer bc2Row=combine1->GetOutput();
1365 // Paste middleRow at index 3 (padded coeff)
1366 typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1367 paster2->SetDestinationImage(paster1->GetOutput());
1368 paster2->SetSourceImage(bc2Row);
1369 destinationIndex[InputDimension-1]=3;
1370 paster2->SetDestinationIndex(destinationIndex);
1371 paster2->SetSourceRegion(bc2Row->GetLargestPossibleRegion());
1373 //---------------------------------
1374 // 3. Next temporal row is identical: paste 2 to 4
1375 typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1376 paster3->SetDestinationImage(paster2->GetOutput());
1377 paster3->SetSourceImage(row2);
1378 destinationIndex[InputDimension-1]=4;
1379 paster3->SetDestinationIndex(destinationIndex);
1380 paster3->SetSourceRegion(row2->GetLargestPossibleRegion());
1382 // ---------------------------------
1383 // 4. Row at index 5=BC (paddedcoeff image) R1
1384 // Paste BC3 row at index 5
1385 typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
1386 paster4->SetDestinationImage(paster3->GetOutput());
1387 paster4->SetSourceImage(row0);
1388 destinationIndex[InputDimension-1]=5;
1389 paster4->SetDestinationIndex(destinationIndex);
1390 paster4->SetSourceRegion(row0->GetLargestPossibleRegion());
1392 // ---------------------------------
1393 // 5. Paste BC4 row at index 6
1394 typename LinearCombinationFilterType::Pointer combine3=LinearCombinationFilterType::New();
1395 combine3->SetInput(0,row0);
1396 combine3->SetInput(1,row2);
1398 combine3->SetB(-1.);
1400 typename CoefficientImageType::Pointer bc4Row=combine3->GetOutput();
1401 typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
1402 paster5->SetDestinationImage(paster4->GetOutput());
1403 paster5->SetSourceImage(bc4Row);
1404 destinationIndex[InputDimension-1]=6;
1405 paster5->SetDestinationIndex(destinationIndex);
1406 paster5->SetSourceRegion(bc4Row->GetLargestPossibleRegion());
1408 // Update the chain!
1410 m_PaddedCoefficientImage= paster5->GetOutput();
1417 // ----------------------------------------------------------------------
1418 // The rabbit with 5 internal CP
1419 // Periodic, constrained to zero at the reference at position 3.5
1420 // and the extremes fixed with anti-symmetrical BC
1421 // Coeff row BC1 0 1 BC2 2 3 BC3 BC4
1422 // PaddedCoeff R: 0 1 2 3 4 5 6 7
1424 // BC2= -weights[3]/weights[1] ( R2+R5 ) - R4
1427 // ---------------------------------------------------------------------
1429 // ---------------------------------
1430 // 0. First Row =BC1
1431 typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1432 typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1433 typename LinearCombinationFilterType::Pointer combine0=LinearCombinationFilterType::New();
1434 combine0->SetInput(0,row0);
1435 combine0->SetInput(1,row1);
1437 combine0->SetB(-1.);
1439 typename CoefficientImageType::Pointer bc1Row=combine0->GetOutput();
1441 typename PasteImageFilterType::Pointer paster0=PasteImageFilterType::New();
1442 paster0->SetDestinationImage(m_PaddedCoefficientImage);
1443 paster0->SetSourceImage(bc1Row);
1444 destinationIndex[InputDimension-1]=0;
1445 paster0->SetDestinationIndex(destinationIndex);
1446 paster0->SetSourceRegion(bc1Row->GetLargestPossibleRegion());
1448 //---------------------------------
1449 // 1. Next two temporal rows are identical: paste 0-1 to 1-2
1450 typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1451 paster1->SetDestinationImage(paster0->GetOutput());
1452 paster1->SetSourceImage(m_CoefficientImage);
1453 sourceRegion.SetIndex(InputDimension-1,0);
1454 destinationIndex[InputDimension-1]=1;
1455 sourceRegion.SetSize(NInputDimensions-1,2);
1456 paster1->SetDestinationIndex(destinationIndex);
1457 paster1->SetSourceRegion(sourceRegion);
1459 //---------------------------------
1460 // 2. Middle row at index 3=BC
1461 // Extract row at index 1, 3 (of coeff image)
1462 //typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1463 typename CoefficientImageType::Pointer row3=this->ExtractTemporalRow(m_CoefficientImage, 3);
1466 typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1467 combine1->SetInput(0,row1);
1468 combine1->SetInput(1,row3);
1473 // Extract row at index 2 (coeff image)
1474 typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
1477 typename LinearCombinationFilterType::Pointer combine2=LinearCombinationFilterType::New();
1478 combine2->SetInput(0,combine1->GetOutput());
1479 combine2->SetInput(1,row2);
1480 combine2->SetA(-m_WeightRatio[3][1]);
1481 combine2->SetB(-1.);
1483 typename CoefficientImageType::Pointer bc2Row=combine2->GetOutput();
1485 // Paste middleRow at index 3 (padded coeff)
1486 typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1487 paster2->SetDestinationImage(paster1->GetOutput());
1488 paster2->SetSourceImage(bc2Row);
1489 destinationIndex[InputDimension-1]=3;
1490 paster2->SetDestinationIndex(destinationIndex);
1491 paster2->SetSourceRegion(bc2Row->GetLargestPossibleRegion());
1493 //---------------------------------
1494 // 3. Next two temporal rows are identical: paste 2,3 to 4,5
1495 typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1496 paster3->SetDestinationImage(paster2->GetOutput());
1497 paster3->SetSourceImage(m_CoefficientImage);
1498 sourceRegion.SetIndex(InputDimension-1,2);
1499 destinationIndex[InputDimension-1]=4;
1500 sourceRegion.SetSize(NInputDimensions-1,2);
1501 paster3->SetDestinationIndex(destinationIndex);
1502 paster3->SetSourceRegion(sourceRegion);
1504 // ---------------------------------
1505 // 4. Row at index 6=BC (paddedcoeff image)R1
1506 // Paste BC3 row at index 6
1507 typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
1508 paster4->SetDestinationImage(paster3->GetOutput());
1509 paster4->SetSourceImage(row0);
1510 destinationIndex[InputDimension-1]=6;
1511 paster4->SetDestinationIndex(destinationIndex);
1512 paster4->SetSourceRegion(row0->GetLargestPossibleRegion());
1514 // ---------------------------------
1515 // 5. Paste BC4 row at index 7
1516 typename LinearCombinationFilterType::Pointer combine3=LinearCombinationFilterType::New();
1517 combine3->SetInput(0,row0);
1518 combine3->SetInput(1,row3);
1520 combine3->SetB(-1.);
1522 typename CoefficientImageType::Pointer bc4Row=combine3->GetOutput();
1523 typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
1524 paster5->SetDestinationImage(paster4->GetOutput());
1525 paster5->SetSourceImage(bc4Row);
1526 destinationIndex[InputDimension-1]=7;
1527 paster5->SetDestinationIndex(destinationIndex);
1528 paster5->SetSourceRegion(bc4Row->GetLargestPossibleRegion());
1530 // Update the chain!
1532 m_PaddedCoefficientImage= paster5->GetOutput();
1539 // ----------------------------------------------------------------------
1540 // The sputnik with 4 internal CP
1541 // Periodic, constrained to zero at the reference
1542 // at position 3 and one indepent extremes copied
1543 // Coeff row BC1 0 1 BC2 2 BC2 3
1544 // PaddedCoeff R: 0 1 2 3 4 5 6
1546 // BC2= -weights[2]/weights[0] ( R2+R4 )
1547 // BC3= weights[2]/weights[0] ( R2-R4) + R1
1548 // ---------------------------------------------------------------------
1550 //---------------------------------
1551 // 1. First Row is equal to last row: paste 3 row to 0
1552 typename CoefficientImageType::Pointer row3=this->ExtractTemporalRow(m_CoefficientImage, 3);
1553 typename PasteImageFilterType::Pointer paster0=PasteImageFilterType::New();
1554 paster0->SetDestinationImage(m_PaddedCoefficientImage);
1555 paster0->SetSourceImage(row3);
1556 destinationIndex[InputDimension-1]=0;
1557 paster0->SetDestinationIndex(destinationIndex);
1558 paster0->SetSourceRegion(row3->GetLargestPossibleRegion());
1560 //---------------------------------
1561 // 1. Next two temporal rows are identical: paste 0-1 to 1-2
1562 typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1563 paster1->SetDestinationImage(paster0->GetOutput());
1564 paster1->SetSourceImage(m_CoefficientImage);
1565 sourceRegion.SetIndex(NInputDimensions-1,0);
1566 destinationIndex[InputDimension-1]=1;
1567 sourceRegion.SetSize(NInputDimensions-1,2);
1568 paster1->SetDestinationIndex(destinationIndex);
1569 paster1->SetSourceRegion(sourceRegion);
1571 //---------------------------------
1572 // 2. Middle row at index 3=BC
1573 // Extract row at index 1, 2 (of coeff image)
1574 typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1575 typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
1578 typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1579 combine1->SetInput(0,row1);
1580 combine1->SetInput(1,row2);
1581 combine1->SetA(-m_WeightRatio[2][0]);
1582 combine1->SetB(-m_WeightRatio[2][0]);
1584 typename CoefficientImageType::Pointer bc1Row=combine1->GetOutput();
1586 // Paste middleRow at index 3 (padded coeff)
1587 typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1588 paster2->SetDestinationImage(paster1->GetOutput());
1589 paster2->SetSourceImage(bc1Row);
1590 destinationIndex[InputDimension-1]=3;
1591 paster2->SetDestinationIndex(destinationIndex);
1592 paster2->SetSourceRegion(bc1Row->GetLargestPossibleRegion());
1594 //---------------------------------
1595 // 3. Next temporal row is identical: paste 2 to 4
1596 typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1597 paster3->SetDestinationImage(paster2->GetOutput());
1598 paster3->SetSourceImage(row2);
1599 destinationIndex[InputDimension-1]=4;
1600 paster3->SetDestinationIndex(destinationIndex);
1601 paster3->SetSourceRegion(row2->GetLargestPossibleRegion());
1603 //---------------------------------
1604 // 4. Final row at index 5=BC3 (paddedcoeff image)R1 R2 R4 corresponds to index in coeff image 0 1 2
1605 // Extract row at index 1, 2 (of coeff image)already done
1606 typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1609 typename LinearCombinationFilterType::Pointer combine3=LinearCombinationFilterType::New();
1610 combine3->SetInput(0,row1);
1611 combine3->SetInput(1,row2);
1613 combine3->SetB(-1.);
1616 typename LinearCombinationFilterType::Pointer combine4=LinearCombinationFilterType::New();
1617 combine4->SetInput(0,combine3->GetOutput());
1618 combine4->SetInput(1,row0);
1619 combine4->SetA(m_WeightRatio[2][0]);
1622 typename CoefficientImageType::Pointer bc2Row=combine4->GetOutput();
1624 // Paste BC row at index 5
1625 typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
1626 paster4->SetDestinationImage(paster3->GetOutput());
1627 paster4->SetSourceImage(bc2Row);
1628 destinationIndex[InputDimension-1]=5;
1629 paster4->SetDestinationIndex(destinationIndex);
1630 paster4->SetSourceRegion(bc2Row->GetLargestPossibleRegion());
1632 // Paste row 3 at index 6
1633 typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
1634 paster5->SetDestinationImage(paster4->GetOutput());
1635 paster5->SetSourceImage(row3);
1636 destinationIndex[InputDimension-1]=6;
1637 paster5->SetDestinationIndex(destinationIndex);
1638 paster5->SetSourceRegion(row3->GetLargestPossibleRegion());
1640 // Update the chain!
1642 m_PaddedCoefficientImage= paster5->GetOutput();
1650 // ----------------------------------------------------------------------
1651 // The sputnik with 5 internal CP
1652 // Periodic, constrained to zero at the reference
1653 // at position 3 and one indepent extreme
1654 // Coeff row BC1 0 1 BC2 2 3 BC3 4
1655 // PaddedCoeff R: 0 1 2 3 4 5 6 7
1656 // BC1= R2 + R5 - R7
1657 // BC2= -weights[3]/weights[1] ( R2+R5 ) - R4
1658 // BC3= R1 + 0.5 R2 - 0.5 R7
1659 // ----------------------------------------------------------------------
1660 //---------------------------------
1662 typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1663 typename CoefficientImageType::Pointer row3=this->ExtractTemporalRow(m_CoefficientImage, 3);
1664 typename CoefficientImageType::Pointer row4=this->ExtractTemporalRow(m_CoefficientImage, 4);
1667 typename LinearCombinationFilterType::Pointer combine0=LinearCombinationFilterType::New();
1668 combine0->SetInput(0,row1);
1669 combine0->SetInput(1,row3);
1672 //combine0->Update();
1673 typename LinearCombinationFilterType::Pointer combine0bis=LinearCombinationFilterType::New();
1674 combine0bis->SetInput(0,combine0->GetOutput());
1675 combine0bis->SetInput(1,row4);
1676 combine0bis->SetA(1.);
1677 combine0bis->SetB(-1.);
1678 combine0bis->Update();
1679 typename CoefficientImageType::Pointer bc1Row=combine0bis->GetOutput();
1681 typename PasteImageFilterType::Pointer paster0=PasteImageFilterType::New();
1682 paster0->SetDestinationImage(m_PaddedCoefficientImage);
1683 paster0->SetSourceImage(bc1Row);
1684 destinationIndex[InputDimension-1]=0;
1685 paster0->SetDestinationIndex(destinationIndex);
1686 paster0->SetSourceRegion(bc1Row->GetLargestPossibleRegion());
1688 //---------------------------------
1689 // 1. Next two temporal rows are identical: paste 0-1 to 1-2
1690 typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1691 paster1->SetDestinationImage(paster0->GetOutput());
1692 paster1->SetSourceImage(m_CoefficientImage);
1693 sourceRegion.SetIndex(NInputDimensions-1,0);
1694 destinationIndex[InputDimension-1]=1;
1695 sourceRegion.SetSize(NInputDimensions-1,2);
1696 paster1->SetDestinationIndex(destinationIndex);
1697 paster1->SetSourceRegion(sourceRegion);
1699 //---------------------------------
1700 // 2. Middle row at index 3=BC
1701 // Extract row at index 1, 2 (of coeff image)
1702 //typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1703 typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
1706 typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1707 combine1->SetInput(0,row1);
1708 combine1->SetInput(1,row3);
1713 typename LinearCombinationFilterType::Pointer combine2=LinearCombinationFilterType::New();
1714 combine2->SetInput(0,combine1->GetOutput());
1715 combine2->SetInput(1,row2);
1716 combine2->SetA(-m_WeightRatio[3][1]);
1717 combine2->SetB(-1.);
1719 typename CoefficientImageType::Pointer bc2Row=combine2->GetOutput();
1721 // Paste middleRow at index 3 (padded coeff)
1722 typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1723 paster2->SetDestinationImage(paster1->GetOutput());
1724 paster2->SetSourceImage(bc2Row);
1725 destinationIndex[InputDimension-1]=3;
1726 paster2->SetDestinationIndex(destinationIndex);
1727 paster2->SetSourceRegion(bc2Row->GetLargestPossibleRegion());
1729 //---------------------------------
1730 // 3. Next two temporal rows are identical: paste 2,3 to 4,5
1731 typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1732 paster3->SetDestinationImage(paster2->GetOutput());
1733 paster3->SetSourceImage(m_CoefficientImage);
1734 sourceRegion.SetIndex(InputDimension-1,2);
1735 destinationIndex[InputDimension-1]=4;
1736 sourceRegion.SetSize(NInputDimensions-1,2);
1737 paster3->SetDestinationIndex(destinationIndex);
1738 paster3->SetSourceRegion(sourceRegion);
1740 //---------------------------------
1741 // 4. Final row at index 6=BC3
1742 typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1745 typename LinearCombinationFilterType::Pointer combine3=LinearCombinationFilterType::New();
1746 combine3->SetInput(0,row1);
1747 combine3->SetInput(1,row4);
1749 combine3->SetB(-1.);
1750 typename LinearCombinationFilterType::Pointer combine4=LinearCombinationFilterType::New();
1751 combine4->SetInput(0,row0);
1752 combine4->SetInput(1,combine3->GetOutput());
1756 typename CoefficientImageType::Pointer bc3Row=combine4->GetOutput();
1758 // Paste BC row at index 6
1759 typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
1760 paster4->SetDestinationImage(paster3->GetOutput());
1761 paster4->SetSourceImage(bc3Row);
1762 destinationIndex[InputDimension-1]=6;
1763 paster4->SetDestinationIndex(destinationIndex);
1764 paster4->SetSourceRegion(bc3Row->GetLargestPossibleRegion());
1766 // Paste row 4 at index 7
1767 typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
1768 paster5->SetDestinationImage(paster4->GetOutput());
1769 paster5->SetSourceImage(row4);
1770 destinationIndex[InputDimension-1]=7;
1771 paster5->SetDestinationIndex(destinationIndex);
1772 paster5->SetSourceRegion(row4->GetLargestPossibleRegion());
1774 // Update the chain!
1776 m_PaddedCoefficientImage= paster5->GetOutput();
1783 // ----------------------------------------------------------------------
1784 // The diamond with 4 internal CP:
1785 // Periodic, constrained to zero at the reference at position 3
1786 // Coeff row 0 1 2 BC1 3 BC2 4
1787 // PaddedCoeff R:0 1 2 3 4 5 6
1788 // BC1= -weights[2]/weights[0] ( R2+R4 )
1789 // BC2= weights[2]/weights[0] ( R0+R2-R4-R6 ) + R1
1790 // ---------------------------------------------------------------------
1792 //---------------------------------
1793 // 1. First Three temporal rows are identical: paste 0-2 to 0-2
1794 typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1795 paster1->SetDestinationImage(m_PaddedCoefficientImage);
1796 paster1->SetSourceImage(m_CoefficientImage);
1797 sourceIndex.Fill(0);
1798 destinationIndex.Fill(0);
1799 sourceSize[NInputDimensions-1]=3;
1800 sourceRegion.SetSize(sourceSize);
1801 sourceRegion.SetIndex(sourceIndex);
1802 paster1->SetDestinationIndex(destinationIndex);
1803 paster1->SetSourceRegion(sourceRegion);
1805 //---------------------------------
1806 // 2. Middle row at index 3=BC
1807 // Extract row at index 0, 4 (of coeff image)
1808 typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
1809 typename CoefficientImageType::Pointer row3=this->ExtractTemporalRow(m_CoefficientImage, 3);
1812 typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1813 combine1->SetInput(0,row2);
1814 combine1->SetInput(1,row3);
1815 combine1->SetA(-m_WeightRatio[2][0]);
1816 combine1->SetB(-m_WeightRatio[2][0]);
1818 typename CoefficientImageType::Pointer bc1Row=combine1->GetOutput();
1820 // Paste middleRow at index 3 (padded coeff)
1821 typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1822 paster2->SetDestinationImage(paster1->GetOutput());
1823 paster2->SetSourceImage(bc1Row);
1824 destinationIndex[InputDimension-1]=3;
1825 paster2->SetDestinationIndex(destinationIndex);
1826 paster2->SetSourceRegion(bc1Row->GetLargestPossibleRegion());
1828 //---------------------------------
1829 // 3. Next row identical: paste 3 to 4
1830 typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1831 paster3->SetDestinationImage(paster2->GetOutput());
1832 paster3->SetSourceImage(row3);
1833 destinationIndex[InputDimension-1]=4;
1834 paster3->SetDestinationIndex(destinationIndex);
1835 paster3->SetSourceRegion(row3->GetLargestPossibleRegion());
1837 //---------------------------------
1838 // 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
1839 // Extract row at index 0, 2 (of coeff image)already done
1840 typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1841 typename CoefficientImageType::Pointer row4=this->ExtractTemporalRow(m_CoefficientImage, 4);
1842 typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1845 typename LinearCombinationFilterType::Pointer combine3=LinearCombinationFilterType::New();
1846 combine3->SetInput(0,row0);
1847 combine3->SetInput(1,row2);
1852 typename LinearCombinationFilterType::Pointer combine4=LinearCombinationFilterType::New();
1853 combine4->SetInput(0,row3);
1854 combine4->SetInput(1,row4);
1859 typename LinearCombinationFilterType::Pointer combine5=LinearCombinationFilterType::New();
1860 combine5->SetInput(0,combine3->GetOutput());
1861 combine5->SetInput(1,combine4->GetOutput());
1863 combine5->SetB(-1.);
1866 typename LinearCombinationFilterType::Pointer combine6=LinearCombinationFilterType::New();
1867 combine6->SetInput(0,combine5->GetOutput());
1868 combine6->SetInput(1,row1);
1869 combine6->SetA(m_WeightRatio[2][0]);
1872 typename CoefficientImageType::Pointer bc2Row=combine6->GetOutput();
1874 // Paste BC row at index 5
1875 typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
1876 paster4->SetDestinationImage(paster3->GetOutput());
1877 paster4->SetSourceImage(bc2Row);
1878 destinationIndex[InputDimension-1]=5;
1879 paster4->SetDestinationIndex(destinationIndex);
1880 paster4->SetSourceRegion(bc2Row->GetLargestPossibleRegion());
1882 // Paste last row at index 6
1883 typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
1884 paster5->SetDestinationImage(paster4->GetOutput());
1885 paster5->SetSourceImage(row4);
1886 destinationIndex[InputDimension-1]=6;
1887 paster5->SetDestinationIndex(destinationIndex);
1888 paster5->SetSourceRegion(row4->GetLargestPossibleRegion());
1890 // Update the chain!
1892 m_PaddedCoefficientImage= paster5->GetOutput();
1898 // ----------------------------------------------------------------------
1899 // The diamond with 5 internal CP:
1900 // periodic, constrained to zero at the reference at position 3.5
1901 // Coeff row 0 1 2 BC1 3 4 BC2 5
1902 // PaddedCoeff R:0 1 2 3 4 5 6 7
1903 // BC1= -weights[3]/weights[1] ( R2+R5 ) - R4
1904 // BC2= weights[2]/weights[0] ( R0+R2-R5-R7 ) + R1
1905 // ---------------------------------------------------------------------
1907 //---------------------------------
1908 // 1. First Three temporal rows are identical: paste 0-2 to 0-2
1909 typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
1910 paster1->SetDestinationImage(m_PaddedCoefficientImage);
1911 paster1->SetSourceImage(m_CoefficientImage);
1912 sourceIndex.Fill(0);
1913 destinationIndex.Fill(0);
1914 sourceSize[NInputDimensions-1]=3;
1915 sourceRegion.SetSize(sourceSize);
1916 sourceRegion.SetIndex(sourceIndex);
1917 paster1->SetDestinationIndex(destinationIndex);
1918 paster1->SetSourceRegion(sourceRegion);
1920 //---------------------------------
1921 // 2. Middle row at index 3=BC
1922 // Extract row at index 0, 4 (of coeff image)
1923 typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
1924 typename CoefficientImageType::Pointer row4=this->ExtractTemporalRow(m_CoefficientImage, 4);
1927 typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
1928 combine1->SetInput(0,row2);
1929 combine1->SetInput(1,row4);
1934 // Extract row at index 3 (coeff image)
1935 typename CoefficientImageType::Pointer row3=this->ExtractTemporalRow(m_CoefficientImage, 3);
1938 typename LinearCombinationFilterType::Pointer combine2=LinearCombinationFilterType::New();
1939 combine2->SetInput(0,combine1->GetOutput());
1940 combine2->SetInput(1,row3);
1941 combine2->SetA(-m_WeightRatio[3][1] );
1942 combine2->SetB(-1.);
1944 typename CoefficientImageType::Pointer bc1Row=combine2->GetOutput();
1946 // Paste middleRow at index 3 (padded coeff)
1947 typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
1948 paster2->SetDestinationImage(paster1->GetOutput());
1949 paster2->SetSourceImage(bc1Row);
1950 destinationIndex[InputDimension-1]=3;
1951 paster2->SetDestinationIndex(destinationIndex);
1952 paster2->SetSourceRegion(bc1Row->GetLargestPossibleRegion());
1954 //---------------------------------
1955 // 3. Next two temporal rows are identical: paste 3,4 to 4,5
1956 typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
1957 paster3->SetDestinationImage(paster2->GetOutput());
1958 paster3->SetSourceImage(m_CoefficientImage);
1959 sourceIndex[InputDimension-1]=3;
1960 destinationIndex[InputDimension-1]=4;
1961 sourceSize[NInputDimensions-1]=2;
1962 sourceRegion.SetSize(sourceSize);
1963 sourceRegion.SetIndex(sourceIndex);
1964 paster3->SetDestinationIndex(destinationIndex);
1965 paster3->SetSourceRegion(sourceRegion);
1967 //---------------------------------
1968 // 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
1969 // Extract row at index 0, 2 (of coeff image)already done
1970 typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
1971 typename CoefficientImageType::Pointer row5=this->ExtractTemporalRow(m_CoefficientImage, 5);
1972 typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
1975 typename LinearCombinationFilterType::Pointer combine3=LinearCombinationFilterType::New();
1976 combine3->SetInput(0,row0);
1977 combine3->SetInput(1,row2);
1982 typename LinearCombinationFilterType::Pointer combine4=LinearCombinationFilterType::New();
1983 combine4->SetInput(0,row4);
1984 combine4->SetInput(1,row5);
1989 typename LinearCombinationFilterType::Pointer combine5=LinearCombinationFilterType::New();
1990 combine5->SetInput(0,combine3->GetOutput());
1991 combine5->SetInput(1,combine4->GetOutput());
1993 combine5->SetB(-1.);
1996 typename LinearCombinationFilterType::Pointer combine6=LinearCombinationFilterType::New();
1997 combine6->SetInput(0,combine5->GetOutput());
1998 combine6->SetInput(1,row1);
1999 combine6->SetA(m_WeightRatio[2][0]);
2002 typename CoefficientImageType::Pointer bc2Row=combine6->GetOutput();
2004 // Paste BC row at index 6
2005 typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
2006 paster4->SetDestinationImage(paster3->GetOutput());
2007 paster4->SetSourceImage(bc2Row);
2008 destinationIndex[InputDimension-1]=6;
2009 paster4->SetDestinationIndex(destinationIndex);
2010 paster4->SetSourceRegion(bc2Row->GetLargestPossibleRegion());
2012 // Paste last row at index 7
2013 typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
2014 paster5->SetDestinationImage(paster4->GetOutput());
2015 paster5->SetSourceImage(row5);
2016 destinationIndex[InputDimension-1]=7;
2017 paster5->SetDestinationIndex(destinationIndex);
2018 paster5->SetSourceRegion(row5->GetLargestPossibleRegion());
2020 // Update the chain!
2022 m_PaddedCoefficientImage= paster5->GetOutput();
2029 // ----------------------------------------------------------------------
2030 // The sputnik with 5 internal CP T''(0)=T''(10)
2031 // Periodic, constrained to zero at the reference
2032 // at position 3 and one indepent extreme
2033 // Coeff row BC1 0 1 BC2 2 3 BC3 4
2034 // PaddedCoeff R: 0 1 2 3 4 5 6 7
2036 // BC2= -weights[3]/weights[1] ( R2+R5 ) - R4
2038 // ---------------------------------------------------------------------
2040 //---------------------------------
2042 typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
2043 typename CoefficientImageType::Pointer row3=this->ExtractTemporalRow(m_CoefficientImage, 3);
2044 typename CoefficientImageType::Pointer row4=this->ExtractTemporalRow(m_CoefficientImage, 4);
2047 typename LinearCombinationFilterType::Pointer combine0=LinearCombinationFilterType::New();
2048 combine0->SetInput(0,row3);
2049 combine0->SetInput(1,row4);
2052 typename LinearCombinationFilterType::Pointer combine0bis=LinearCombinationFilterType::New();
2053 combine0bis->SetInput(0,combine0->GetOutput());
2054 combine0bis->SetInput(1,row1);
2055 combine0bis->SetA(1.);
2056 combine0bis->SetB(-1.);
2057 combine0bis->Update();
2058 typename CoefficientImageType::Pointer bc1Row=combine0bis->GetOutput();
2060 typename PasteImageFilterType::Pointer paster0=PasteImageFilterType::New();
2061 paster0->SetDestinationImage(m_PaddedCoefficientImage);
2062 paster0->SetSourceImage(bc1Row);
2063 destinationIndex[InputDimension-1]=0;
2064 paster0->SetDestinationIndex(destinationIndex);
2065 paster0->SetSourceRegion(bc1Row->GetLargestPossibleRegion());
2067 //---------------------------------
2068 // 1. Next two temporal rows are identical: paste 0-1 to 1-2
2069 typename PasteImageFilterType::Pointer paster1=PasteImageFilterType::New();
2070 paster1->SetDestinationImage(paster0->GetOutput());
2071 paster1->SetSourceImage(m_CoefficientImage);
2072 sourceRegion.SetIndex(NInputDimensions-1,0);
2073 destinationIndex[InputDimension-1]=1;
2074 sourceRegion.SetSize(NInputDimensions-1,2);
2075 paster1->SetDestinationIndex(destinationIndex);
2076 paster1->SetSourceRegion(sourceRegion);
2078 //---------------------------------
2079 // 2. Middle row at index 3=BC
2080 // Extract row at index 1, 2 (of coeff image)
2081 // typename CoefficientImageType::Pointer row1=this->ExtractTemporalRow(m_CoefficientImage, 1);
2082 typename CoefficientImageType::Pointer row2=this->ExtractTemporalRow(m_CoefficientImage, 2);
2085 typename LinearCombinationFilterType::Pointer combine1=LinearCombinationFilterType::New();
2086 combine1->SetInput(0,row1);
2087 combine1->SetInput(1,row3);
2092 typename LinearCombinationFilterType::Pointer combine2=LinearCombinationFilterType::New();
2093 combine2->SetInput(0,combine1->GetOutput());
2094 combine2->SetInput(1,row2);
2095 combine2->SetA(-m_WeightRatio[3][1]);
2096 combine2->SetB(-1.);
2098 typename CoefficientImageType::Pointer bc2Row=combine2->GetOutput();
2100 // Paste middleRow at index 3 (padded coeff)
2101 typename PasteImageFilterType::Pointer paster2=PasteImageFilterType::New();
2102 paster2->SetDestinationImage(paster1->GetOutput());
2103 paster2->SetSourceImage(bc2Row);
2104 destinationIndex[InputDimension-1]=3;
2105 paster2->SetDestinationIndex(destinationIndex);
2106 paster2->SetSourceRegion(bc2Row->GetLargestPossibleRegion());
2108 //---------------------------------
2109 // 3. Next two temporal rows are identical: paste 2,3 to 4,5
2110 typename PasteImageFilterType::Pointer paster3=PasteImageFilterType::New();
2111 paster3->SetDestinationImage(paster2->GetOutput());
2112 paster3->SetSourceImage(m_CoefficientImage);
2113 sourceRegion.SetIndex(InputDimension-1,2);
2114 destinationIndex[InputDimension-1]=4;
2115 sourceRegion.SetSize(NInputDimensions-1,2);
2116 paster3->SetDestinationIndex(destinationIndex);
2117 paster3->SetSourceRegion(sourceRegion);
2119 //---------------------------------
2120 // 4. Final row at index 6=BC3
2121 typename CoefficientImageType::Pointer row0=this->ExtractTemporalRow(m_CoefficientImage, 0);
2124 // typename LinearCombinationFilterType::Pointer combine3=LinearCombinationFilterType::New();
2125 // combine3->SetInput(0,row0);
2126 // combine3->SetInput(1,row1);
2127 // combine3->SetA(1.);
2128 // combine3->SetB(0.5);
2129 // combine3->Update();
2130 // typename CoefficientImageType::Pointer bc3Row=combine3->GetOutput();
2132 // Paste BC row at index 6
2133 typename PasteImageFilterType::Pointer paster4=PasteImageFilterType::New();
2134 paster4->SetDestinationImage(paster3->GetOutput());
2135 paster4->SetSourceImage(row0);
2136 destinationIndex[InputDimension-1]=6;
2137 paster4->SetDestinationIndex(destinationIndex);
2138 paster4->SetSourceRegion(row0->GetLargestPossibleRegion());
2140 // Paste row 4 at index 7
2141 typename PasteImageFilterType::Pointer paster5=PasteImageFilterType::New();
2142 paster5->SetDestinationImage(paster4->GetOutput());
2143 paster5->SetSourceImage(row4);
2144 destinationIndex[InputDimension-1]=7;
2145 paster5->SetDestinationIndex(destinationIndex);
2146 paster5->SetSourceRegion(row4->GetLargestPossibleRegion());
2148 // Update the chain!
2150 m_PaddedCoefficientImage= paster5->GetOutput();
2156 DD ("Shape not available");
2162 // // Extract coefficients from padded version
2163 // template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
2165 // ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2166 // ::ExtractCoefficientImage( )
2168 // ////DD("extract coeff image");
2169 // typename ExtractImageFilterType::Pointer extractImageFilter = ExtractImageFilterType::New();
2170 // typename CoefficientImageType::RegionType extractionRegion=m_PaddedCoefficientImage->GetLargestPossibleRegion();
2171 // typename CoefficientImageType::RegionType::SizeType extractionSize=extractionRegion.GetSize();
2172 // typename CoefficientImageType::RegionType::IndexType extractionIndex = extractionRegion.GetIndex();
2173 // extractionSize[InputDimension-1]-=4;
2174 // extractionIndex[InputDimension-1]=2;
2175 // extractionRegion.SetSize(extractionSize);
2176 // extractionRegion.SetIndex(extractionIndex);
2177 // extractImageFilter->SetInput(m_PaddedCoefficientImage);
2178 // extractImageFilter->SetExtractionRegion(extractionRegion);
2179 // extractImageFilter->Update();
2180 // m_CoefficientImage=extractImageFilter->GetOutput();
2185 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
2187 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2188 ::PrintSelf(std::ostream &os, itk::Indent indent) const
2191 this->Superclass::PrintSelf(os, indent);
2193 os << indent << "GridRegion: " << m_GridRegion << std::endl;
2194 os << indent << "GridOrigin: " << m_GridOrigin << std::endl;
2195 os << indent << "GridSpacing: " << m_GridSpacing << std::endl;
2196 os << indent << "GridDirection: " << m_GridDirection << std::endl;
2197 os << indent << "IndexToPoint: " << m_IndexToPoint << std::endl;
2198 os << indent << "PointToIndex: " << m_PointToIndex << std::endl;
2200 os << indent << "CoefficientImage: [ ";
2201 os << m_CoefficientImage.GetPointer() << " ]" << std::endl;
2203 os << indent << "WrappedImage: [ ";
2204 os << m_WrappedImage.GetPointer() << " ]" << std::endl;
2206 os << indent << "InputParametersPointer: "
2207 << m_InputParametersPointer << std::endl;
2208 os << indent << "ValidRegion: " << m_ValidRegion << std::endl;
2209 os << indent << "LastJacobianIndex: " << m_LastJacobianIndex << std::endl;
2210 os << indent << "BulkTransform: ";
2211 os << m_BulkTransform << std::endl;
2213 if ( m_BulkTransform )
2215 os << indent << "BulkTransformType: "
2216 << m_BulkTransform->GetNameOfClass() << std::endl;
2218 os << indent << "VectorBSplineInterpolator: ";
2219 os << m_VectorInterpolator.GetPointer() << std::endl;
2220 os << indent << "Mask: ";
2221 os << m_Mask<< std::endl;
2226 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
2228 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2229 ::InsideValidRegion( const ContinuousIndexType& index ) const
2233 if ( !m_ValidRegion.IsInside( index ) )
2238 // JV verify all dimensions
2241 typedef typename ContinuousIndexType::ValueType ValueType;
2242 for( unsigned int j = 0; j < NInputDimensions; j++ )
2244 if (m_SplineOrderOdd[j])
2246 if ( index[j] >= static_cast<ValueType>( m_ValidRegionLast[j] ) )
2258 // Transform a point
2259 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
2260 typename ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2262 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2263 ::TransformPoint(const InputPointType &inputPoint) const
2266 InputPointType transformedPoint;
2267 OutputPointType outputPoint;
2270 if ( m_BulkTransform )
2272 transformedPoint = m_BulkTransform->TransformPoint( inputPoint );
2276 transformedPoint = inputPoint;
2279 // Deformable transform
2280 if ( m_PaddedCoefficientImage )
2282 // Check if inside mask
2283 if(m_Mask && !(m_Mask->IsInside(inputPoint) ) )
2285 // Outside: no (deformable) displacement
2286 return transformedPoint;
2289 // Check if inside valid region
2291 ContinuousIndexType index;
2292 this->TransformPointToContinuousIndex( inputPoint, index );
2293 inside = this->InsideValidRegion( index );
2296 // Outside: no (deformable) displacement
2297 DD("outside valid region");
2298 return transformedPoint;
2301 // Call the vector interpolator
2302 itk::Vector<TCoordRep,SpaceDimension> displacement=m_VectorInterpolator->EvaluateAtContinuousIndex(index);
2304 // JV add for the spatial dimensions
2305 outputPoint=transformedPoint;
2306 for (unsigned int i=0; i<NInputDimensions-1; i++)
2307 outputPoint[i] += displacement[i];
2313 itkWarningMacro( << "B-spline coefficients have not been set" );
2314 outputPoint = transformedPoint;
2322 //JV Deformably transform a point
2323 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
2324 typename ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2326 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2327 ::DeformablyTransformPoint(const InputPointType &inputPoint) const
2329 OutputPointType outputPoint;
2330 if ( m_PaddedCoefficientImage )
2333 // Check if inside mask
2334 if(m_Mask && !(m_Mask->IsInside(inputPoint) ) )
2336 // Outside: no (deformable) displacement
2342 ContinuousIndexType index;
2343 this->TransformPointToContinuousIndex( inputPoint, index );
2344 inside = this->InsideValidRegion( index );
2348 //outside: no (deformable) displacement
2349 outputPoint = inputPoint;
2353 // Call the vector interpolator
2354 itk::Vector<TCoordRep,SpaceDimension> displacement=m_VectorInterpolator->EvaluateAtContinuousIndex(index);
2356 // JV add for the spatial dimensions
2357 outputPoint=inputPoint;
2358 for (unsigned int i=0; i<NInputDimensions-1; i++)
2359 outputPoint[i] += displacement[i];
2362 // No coefficients available
2365 itkWarningMacro( << "B-spline coefficients have not been set" );
2366 outputPoint = inputPoint;
2373 // JV weights are identical as for transformpoint, could be done simultaneously in metric!!!!
2374 // Compute the Jacobian in one position
2375 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
2377 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2378 ::ComputeJacobianWithRespectToParameters( const InputPointType & point, JacobianType & jacobian) const
2381 //========================================================
2382 // Zero all components of jacobian
2383 //========================================================
2384 // JV not thread safe (m_LastJacobianIndex), instantiate N transforms
2385 // NOTE: for efficiency, we only need to zero out the coefficients
2386 // that got fill last time this method was called.
2387 unsigned int j=0,b=0;
2389 //Set the previously-set to zero
2390 for ( j = 0; j < SpaceDimension; j++ )
2392 m_FirstIterator[j].GoToBegin();
2393 while ( !m_FirstIterator[j].IsAtEnd() )
2395 m_FirstIterator[j].Set( m_ZeroVector );
2396 ++(m_FirstIterator[j]);
2400 //Set the previously-set to zero
2401 for ( j = 0; j < SpaceDimension; j++ )
2403 m_SecondIterator[j].GoToBegin();
2404 while ( ! (m_SecondIterator[j]).IsAtEnd() )
2406 m_SecondIterator[j].Set( m_ZeroVector );
2407 ++(m_SecondIterator[j]);
2411 //Set the previously-set to zero
2413 for ( j = 0; j < SpaceDimension; j++ )
2415 m_ThirdIterator[j].GoToBegin();
2416 while ( ! (m_ThirdIterator[j]).IsAtEnd() )
2418 m_ThirdIterator[j].Set( m_ZeroVector );
2419 ++(m_ThirdIterator[j]);
2423 //Set the previously-set to zero
2425 for (b=0; b<m_BCSize;b++)
2426 if ( !( m_FirstRegion.IsInside(m_BCRegions[b]) | m_SecondRegion.IsInside(m_BCRegions[b]) ) )
2427 for ( j = 0; j < SpaceDimension; j++ )
2429 m_BCIterators[j][b].GoToBegin();
2430 while ( ! (m_BCIterators[j][b]).IsAtEnd() )
2432 m_BCIterators[j][b].Set( m_ZeroVector );
2433 ++(m_BCIterators[j][b]);
2437 //Set the previously-set to zero
2439 for (b=0; b<m_BC2Size;b++)
2440 if ( !( m_FirstRegion.IsInside(m_BC2Regions[b]) | m_SecondRegion.IsInside(m_BC2Regions[b]) ) )
2441 for ( j = 0; j < SpaceDimension; j++ )
2443 m_BC2Iterators[j][b].GoToBegin();
2444 while ( ! (m_BC2Iterators[j][b]).IsAtEnd() )
2446 m_BC2Iterators[j][b].Set( m_ZeroVector );
2447 ++(m_BC2Iterators[j][b]);
2451 //Set the previously-set to zero
2453 for (b=0; b<m_BC3Size;b++)
2454 if ( !( m_FirstRegion.IsInside(m_BC3Regions[b]) | m_SecondRegion.IsInside(m_BC3Regions[b]) ) )
2455 for ( j = 0; j < SpaceDimension; j++ )
2457 m_BC3Iterators[j][b].GoToBegin();
2458 while ( ! (m_BC3Iterators[j][b]).IsAtEnd() )
2460 m_BC3Iterators[j][b].Set( m_ZeroVector );
2461 ++(m_BC3Iterators[j][b]);
2466 //========================================================
2467 // For each dimension, copy the weight to the support region
2468 //========================================================
2470 // Check if inside mask
2471 if(m_Mask && !(m_Mask->IsInside(point) ) )
2473 // Outside: no (deformable) displacement
2474 jacobian = m_SharedDataBSplineJacobian;
2479 this->TransformPointToContinuousIndex( point, m_Index );
2481 // NOTE: if the support region does not lie totally within the grid
2482 // we assume zero displacement and return the input point
2483 if ( !this->InsideValidRegion( m_Index ) )
2485 jacobian = m_SharedDataBSplineJacobian;
2489 // Compute interpolation weights
2490 const WeightsDataType *weights=NULL;
2491 m_VectorInterpolator->EvaluateWeightsAtContinuousIndex( m_Index, &weights, m_LastJacobianIndex);
2494 m_SupportRegion.SetIndex( m_LastJacobianIndex );
2495 WrapRegion(m_SupportRegion, m_FirstRegion, m_SecondRegion, m_ThirdRegion, m_BCRegions, m_BCValues, m_BC2Regions, m_BC2Values, m_BC3Regions, m_BC3Values, m_InitialOffset);
2496 m_ThirdSize=m_ThirdRegion.GetSize()[InputDimension-1];
2497 m_BCSize=m_BCRegions.size();
2498 m_BC2Size=m_BC2Regions.size();
2499 m_BC3Size=m_BC3Regions.size();
2501 // Reset the iterators
2502 for ( j = 0; j < SpaceDimension ; j++ )
2504 m_FirstIterator[j] = IteratorType( m_JacobianImage[j], m_FirstRegion);
2505 m_SecondIterator[j] = IteratorType( m_JacobianImage[j], m_SecondRegion);
2506 if(m_ThirdSize) m_ThirdIterator[j] = IteratorType( m_JacobianImage[j], m_ThirdRegion);
2508 m_BCIterators[j].resize(m_BCSize);
2509 for (b=0; b<m_BCSize;b++)
2510 m_BCIterators[j][b]= IteratorType( m_JacobianImage[j], m_BCRegions[b]);
2511 m_BC2Iterators[j].resize(m_BC2Size);
2512 for (b=0; b<m_BC2Size;b++)
2513 m_BC2Iterators[j][b]= IteratorType( m_JacobianImage[j], m_BC2Regions[b]);
2514 m_BC3Iterators[j].resize(m_BC3Size);
2515 for (b=0; b<m_BC3Size;b++)
2516 m_BC3Iterators[j][b]= IteratorType( m_JacobianImage[j], m_BC3Regions[b]);
2519 // Skip if on a fixed condition
2522 if (m_BCSize) weights+=m_InitialOffset*m_BCRegions[0].GetNumberOfPixels();
2523 else std::cerr<<"InitialOffset without BCSize: Error!!!!!!!"<<std::endl;
2526 //copy weight to jacobian image
2527 for ( j = 0; j < SpaceDimension; j++ )
2529 // For each dimension, copy the weight to the support region
2530 while ( ! (m_FirstIterator[j]).IsAtEnd() )
2532 m_ZeroVector[j]=*weights;
2533 (m_FirstIterator[j]).Set( m_ZeroVector);
2534 ++(m_FirstIterator[j]);
2538 m_ZeroVector[j]=itk::NumericTraits<JacobianValueType>::Zero;
2539 if (j != SpaceDimension-1) weights-=m_FirstRegion.GetNumberOfPixels();
2542 // Skip BC1 and go to the second region
2543 if (m_BCSize) weights+=m_BCRegions[0].GetNumberOfPixels();
2545 // For each dimension, copy the weight to the support region
2546 //copy weight to jacobian image
2547 for ( j = 0; j < SpaceDimension; j++ )
2549 while ( ! (m_SecondIterator[j]).IsAtEnd() )
2551 m_ZeroVector[j]=*weights;
2552 (m_SecondIterator[j]).Set( m_ZeroVector);
2553 ++(m_SecondIterator[j]);
2556 weights-=m_SecondRegion.GetNumberOfPixels();
2557 m_ZeroVector[j]=itk::NumericTraits<JacobianValueType>::Zero;
2563 // Put pointer in correct position
2564 weights-=m_BCRegions[0].GetNumberOfPixels();
2566 for ( j = 0; j < SpaceDimension; j++ )
2568 for ( b=0; b < m_BCSize; b++ )
2570 while ( ! (m_BCIterators[j][b]).IsAtEnd() )
2572 //copy weight to jacobian image
2573 m_ZeroVector[j]=(*weights) * m_BCValues[b];
2574 (m_BCIterators[j][b]).Value()+= m_ZeroVector;
2575 ++(m_BCIterators[j][b]);
2578 weights-=m_BCRegions[b].GetNumberOfPixels();
2580 m_ZeroVector[j]=itk::NumericTraits<JacobianValueType>::Zero;
2582 weights+=m_BCRegions[m_BCSize-1].GetNumberOfPixels();
2585 // Add the BC2 to the weights
2588 // Move further in the weights pointer
2589 weights+=m_SecondRegion.GetNumberOfPixels();
2591 for ( j = 0; j < SpaceDimension; j++ )
2593 for ( b=0; b < m_BC2Size; b++ )
2595 while ( ! (m_BC2Iterators[j][b]).IsAtEnd() )
2597 //copy weight to jacobian image
2598 m_ZeroVector[j]=(*weights) * m_BC2Values[b];
2599 (m_BC2Iterators[j][b]).Value()+= m_ZeroVector;
2600 ++(m_BC2Iterators[j][b]);
2603 weights-=m_BC2Regions[b].GetNumberOfPixels();
2605 m_ZeroVector[j]=itk::NumericTraits<JacobianValueType>::Zero;
2607 // Move further in the weights pointer
2608 weights+=m_BC2Regions[m_BC2Size-1].GetNumberOfPixels();
2614 // For each dimension, copy the weight to the support region
2615 //copy weight to jacobian image
2616 for ( j = 0; j < SpaceDimension; j++ )
2618 while ( ! (m_ThirdIterator[j]).IsAtEnd() )
2620 m_ZeroVector[j]=*weights;
2621 (m_ThirdIterator[j]).Value()+= m_ZeroVector;
2622 ++(m_ThirdIterator[j]);
2626 // Move further in the weights pointer?
2627 if (j != SpaceDimension-1) weights-=m_ThirdRegion.GetNumberOfPixels();
2628 m_ZeroVector[j]=itk::NumericTraits<JacobianValueType>::Zero;
2632 // Add the BC3 to the weights
2635 for ( j = 0; j < SpaceDimension; j++ )
2637 for ( b=0; b < m_BC3Size; b++ )
2639 while ( ! (m_BC3Iterators[j][b]).IsAtEnd() )
2641 //copy weight to jacobian image
2642 m_ZeroVector[j]=(*weights) * m_BC3Values[b];
2643 (m_BC3Iterators[j][b]).Value()+= m_ZeroVector;
2644 ++(m_BC3Iterators[j][b]);
2647 weights-=m_BC3Regions[b].GetNumberOfPixels();
2649 m_ZeroVector[j]=itk::NumericTraits<JacobianValueType>::Zero;
2653 // Return the result
2654 jacobian = m_SharedDataBSplineJacobian;
2658 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
2660 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
2661 ::WrapRegion( const RegionType & m_SupportRegion,
2662 RegionType & m_FirstRegion,
2663 RegionType & m_SecondRegion,
2664 RegionType & m_ThirdRegion,
2665 std::vector<RegionType>& m_BCRegions,std::vector<double>& m_BCValues,
2666 std::vector<RegionType>& m_BC2Regions,std::vector<double>& m_BC2Values,
2667 std::vector<RegionType>& m_BC3Regions,std::vector<double>& m_BC3Values,
2668 unsigned int& m_InitialOffset ) const
2671 // Synchronize regions
2673 m_FirstRegion=m_SupportRegion;
2674 m_BCRegion=m_SupportRegion;
2675 m_BCRegion.SetSize(InputDimension-1,1);
2676 m_SecondRegion=m_SupportRegion;
2677 m_ThirdRegion=m_SupportRegion;
2678 m_ThirdRegion.SetSize(InputDimension-1,0);
2679 m_BC3Regions.resize(0);
2682 // BC depends on shape
2683 switch(m_TransformShape)
2688 2: rabbit 4 CP 3 DOF
2689 3: rabbit 5 CP 4 DOF
2690 4: sputnik 4 CP 4 DOF
2691 5: sputnik 5 CP 5 DOF
2692 6: diamond 6 CP 5 DOF
2693 7: diamond 7 CP 6 DOF
2698 // ----------------------------------------------------------------------
2699 // The egg with 4 internal CP (starting from inhale)
2700 // Periodic, constrained to zero at the reference
2701 // at position 3 and
2702 // Coeff row BC1 BC2 0 BC3 1 2 BC4
2703 // PaddedCoeff R: 0 1 2 3 4 5 6
2706 // BC3= -weights[2]/weights[0] ( R2+R4 )
2708 // ---------------------------------------------------------------------
2709 switch(m_SupportRegion.GetIndex(InputDimension-1))
2714 m_FirstRegion.SetSize(InputDimension-1,2);
2715 m_FirstRegion.SetIndex(InputDimension-1,1);
2718 m_BCRegions.resize(0);
2721 m_SecondRegion.SetSize(InputDimension-1,1);
2722 m_SecondRegion.SetIndex(InputDimension-1,0);
2725 m_BC2Regions.resize(2);
2726 m_BC2Values.resize(2);
2727 m_BCRegion.SetIndex(InputDimension-1,0);
2728 m_BC2Regions[0]=m_BCRegion;
2729 m_BC2Values[0]=-m_WeightRatio[2][0];
2730 m_BCRegion.SetIndex(InputDimension-1,1);
2731 m_BC2Regions[1]=m_BCRegion;
2732 m_BC2Values[1]= -m_WeightRatio[2][0];
2739 m_FirstRegion.SetSize(InputDimension-1,1);
2740 m_FirstRegion.SetIndex(InputDimension-1,2);
2743 m_BCRegions.resize(0);
2746 m_SecondRegion.SetSize(InputDimension-1,1);
2747 m_SecondRegion.SetIndex(InputDimension-1,0);
2750 m_BC2Regions.resize(2);
2751 m_BC2Values.resize(2);
2752 m_BCRegion.SetIndex(InputDimension-1,0);
2753 m_BC2Regions[0]=m_BCRegion;
2754 m_BC2Values[0]=-m_WeightRatio[2][0];
2755 m_BCRegion.SetIndex(InputDimension-1,1);
2756 m_BC2Regions[1]=m_BCRegion;
2757 m_BC2Values[1]= -m_WeightRatio[2][0];
2760 m_FirstRegion.SetSize(InputDimension-1,1);
2761 m_FirstRegion.SetIndex(InputDimension-1,1);
2768 m_FirstRegion.SetSize(InputDimension-1,1);
2769 m_FirstRegion.SetIndex(InputDimension-1,0);
2772 m_BCRegions.resize(2);
2773 m_BCValues.resize(2);
2774 m_BCRegion.SetIndex(InputDimension-1,0);
2775 m_BCRegions[0]=m_BCRegion;
2776 m_BCValues[0]=-m_WeightRatio[2][0];
2777 m_BCRegion.SetIndex(InputDimension-1,1);
2778 m_BCRegions[1]=m_BCRegion;
2779 m_BCValues[1]= -m_WeightRatio[2][0];
2782 m_SecondRegion.SetSize(InputDimension-1,2);
2783 m_SecondRegion.SetIndex(InputDimension-1,1);
2786 m_BC2Regions.resize(0);
2793 m_FirstRegion.SetSize(InputDimension-1,0);
2796 m_BCRegions.resize(2);
2797 m_BCValues.resize(2);
2798 m_BCRegion.SetIndex(InputDimension-1,0);
2799 m_BCRegions[0]=m_BCRegion;
2800 m_BCValues[0]=-m_WeightRatio[2][0];
2801 m_BCRegion.SetIndex(InputDimension-1,1);
2802 m_BCRegions[1]=m_BCRegion;
2803 m_BCValues[1]= -m_WeightRatio[2][0];
2806 m_SecondRegion.SetSize(InputDimension-1,2);
2807 m_SecondRegion.SetIndex(InputDimension-1,1);
2810 m_BC2Regions.resize(1);
2811 m_BC2Values.resize(1);
2812 m_BCRegion.SetIndex(InputDimension-1,0);
2813 m_BC2Regions[0]=m_BCRegion;
2821 DD("supportindex > 3 ???");
2822 DD(m_SupportRegion.GetIndex(InputDimension-1));
2824 } // end switch index
2831 // ----------------------------------------------------------------------
2832 // The egg with 5 internal CP (starting from inhale)
2833 // Periodic, constrained to zero at the reference
2834 // at position 3 and
2835 // Coeff row BC1 BC2 0 BC3 1 2 3 BC4
2836 // PaddedCoeff R: 0 1 2 3 4 5 6 7
2839 // BC3= -weights[3]/weights[1] ( R2+R5 ) - R4
2841 // ---------------------------------------------------------------------
2842 switch(m_SupportRegion.GetIndex(InputDimension-1))
2847 m_FirstRegion.SetSize(InputDimension-1,2);
2848 m_FirstRegion.SetIndex(InputDimension-1,2);
2851 m_BCRegions.resize(0);
2854 m_SecondRegion.SetSize(InputDimension-1,1);
2855 m_SecondRegion.SetIndex(InputDimension-1,0);
2858 m_BC2Regions.resize(3);
2859 m_BC2Values.resize(3);
2860 m_BCRegion.SetIndex(InputDimension-1,0);
2861 m_BC2Regions[0]=m_BCRegion;
2862 m_BC2Values[0]=-m_WeightRatio[3][1];
2863 m_BCRegion.SetIndex(InputDimension-1,1);
2864 m_BC2Regions[1]=m_BCRegion;
2866 m_BCRegion.SetIndex(InputDimension-1,2);
2867 m_BC2Regions[2]=m_BCRegion;
2868 m_BC2Values[2]=-m_WeightRatio[3][1];
2875 m_FirstRegion.SetSize(InputDimension-1,1);
2876 m_FirstRegion.SetIndex(InputDimension-1,3);
2879 m_BCRegions.resize(0);
2882 m_SecondRegion.SetSize(InputDimension-1,1);
2883 m_SecondRegion.SetIndex(InputDimension-1,0);
2886 m_BC2Regions.resize(3);
2887 m_BC2Values.resize(3);
2888 m_BCRegion.SetIndex(InputDimension-1,0);
2889 m_BC2Regions[0]=m_BCRegion;
2890 m_BC2Values[0]=-m_WeightRatio[3][1];
2891 m_BCRegion.SetIndex(InputDimension-1,1);
2892 m_BC2Regions[1]=m_BCRegion;
2894 m_BCRegion.SetIndex(InputDimension-1,2);
2895 m_BC2Regions[2]=m_BCRegion;
2896 m_BC2Values[2]=-m_WeightRatio[3][1];
2899 m_FirstRegion.SetSize(InputDimension-1,1);
2900 m_FirstRegion.SetIndex(InputDimension-1,1);
2907 m_FirstRegion.SetSize(InputDimension-1,1);
2908 m_FirstRegion.SetIndex(InputDimension-1,0);
2911 m_BCRegions.resize(3);
2912 m_BCValues.resize(3);
2913 m_BCRegion.SetIndex(InputDimension-1,0);
2914 m_BCRegions[0]=m_BCRegion;
2915 m_BCValues[0]=-m_WeightRatio[3][1];
2916 m_BCRegion.SetIndex(InputDimension-1,1);
2917 m_BCRegions[1]=m_BCRegion;
2919 m_BCRegion.SetIndex(InputDimension-1,2);
2920 m_BCRegions[2]=m_BCRegion;
2921 m_BCValues[2]=-m_WeightRatio[3][1];
2924 m_SecondRegion.SetSize(InputDimension-1,2);
2925 m_SecondRegion.SetIndex(InputDimension-1,1);
2928 m_BC2Regions.resize(0);
2935 m_FirstRegion.SetSize(InputDimension-1,0);
2938 m_BCRegions.resize(3);
2939 m_BCValues.resize(3);
2940 m_BCRegion.SetIndex(InputDimension-1,0);
2941 m_BCRegions[0]=m_BCRegion;
2942 m_BCValues[0]=-m_WeightRatio[3][1];
2943 m_BCRegion.SetIndex(InputDimension-1,1);
2944 m_BCRegions[1]=m_BCRegion;
2946 m_BCRegion.SetIndex(InputDimension-1,2);
2947 m_BCRegions[2]=m_BCRegion;
2948 m_BCValues[2]=-m_WeightRatio[3][1];
2951 m_SecondRegion.SetSize(InputDimension-1,3);
2952 m_SecondRegion.SetIndex(InputDimension-1,1);
2955 m_BC2Regions.resize(0);
2962 m_FirstRegion.SetSize(InputDimension-1,3);
2963 m_FirstRegion.SetIndex(InputDimension-1,1);
2966 m_BCRegions.resize(0);
2969 m_SecondRegion.SetSize(InputDimension-1,1);
2970 m_SecondRegion.SetIndex(InputDimension-1,0);
2973 m_BC2Regions.resize(0);
2980 DD("supportindex > 3 ???");
2981 DD(m_SupportRegion.GetIndex(InputDimension-1));
2983 } // end switch index
2987 // // ----------------------------------------------------------------------
2988 // // The egg with 5 internal CP:
2989 // // Periodic, C2 smooth everywhere and constrained to zero at the reference
2990 // // Coeff row R5 BC1 0 1 2 3 BC2 R2
2991 // // PaddedCoeff R: 0 1 2 3 4 5 6 7
2992 // // BC1= -weights[2]/weights[0] ( R2+R5)
2994 // // ---------------------------------------------------------------------
2995 // switch(m_SupportRegion.GetIndex(InputDimension-1))
3000 // m_FirstRegion.SetSize(InputDimension-1,1);
3001 // m_FirstRegion.SetIndex(InputDimension-1,3);
3004 // m_BCRegions.resize(2);
3005 // m_BCValues.resize(2);
3006 // m_BCRegion.SetIndex(InputDimension-1,0);
3007 // m_BCRegions[0]=m_BCRegion;
3008 // m_BCValues[0]=-m_WeightRatio[2][0];
3009 // m_BCRegion.SetIndex(InputDimension-1,3);
3010 // m_BCRegions[1]=m_BCRegion;
3011 // m_BCValues[1]=-m_WeightRatio[2][0];
3014 // m_SecondRegion.SetSize(InputDimension-1,2);
3015 // m_SecondRegion.SetIndex(InputDimension-1,0);
3018 // m_BC2Regions.resize(0);
3025 // m_FirstRegion.SetSize(InputDimension-1,0);
3028 // m_BCRegions.resize(2);
3029 // m_BCValues.resize(2);
3030 // m_BCRegion.SetIndex(InputDimension-1,0);
3031 // m_BCRegions[0]=m_BCRegion;
3032 // m_BCValues[0]=-m_WeightRatio[2][0];
3033 // m_BCRegion.SetIndex(InputDimension-1,3);
3034 // m_BCRegions[1]=m_BCRegion;
3035 // m_BCValues[1]=-m_WeightRatio[2][0];
3038 // m_SecondRegion.SetSize(InputDimension-1,3);
3039 // m_SecondRegion.SetIndex(InputDimension-1,0);
3042 // m_BC2Regions.resize(0);
3049 // m_FirstRegion.SetSize(InputDimension-1,4);
3050 // m_FirstRegion.SetIndex(InputDimension-1,0);
3053 // m_BCRegions.resize(0);
3056 // m_SecondRegion.SetSize(InputDimension-1,0);
3059 // m_BC2Regions.resize(0);
3066 // m_FirstRegion.SetSize(InputDimension-1,3);
3067 // m_FirstRegion.SetIndex(InputDimension-1,1);
3070 // m_BCRegions.resize(2);
3071 // m_BCValues.resize(2);
3072 // m_BCRegion.SetIndex(InputDimension-1,0);
3073 // m_BCRegions[0]=m_BCRegion;
3074 // m_BCValues[0]=-m_WeightRatio[2][0];
3075 // m_BCRegion.SetIndex(InputDimension-1,3);
3076 // m_BCRegions[1]=m_BCRegion;
3077 // m_BCValues[1]=-m_WeightRatio[2][0];
3080 // m_SecondRegion.SetSize(InputDimension-1,0);
3083 // m_BC2Regions.resize(0);
3091 // m_FirstRegion.SetSize(InputDimension-1,2);
3092 // m_FirstRegion.SetIndex(InputDimension-1,2);
3095 // m_BCRegions.resize(2);
3096 // m_BCValues.resize(2);
3097 // m_BCRegion.SetIndex(InputDimension-1,0);
3098 // m_BCRegions[0]=m_BCRegion;
3099 // m_BCValues[0]=-m_WeightRatio[2][0];
3100 // m_BCRegion.SetIndex(InputDimension-1,3);
3101 // m_BCRegions[1]=m_BCRegion;
3102 // m_BCValues[1]=-m_WeightRatio[2][0];
3105 // m_SecondRegion.SetSize(InputDimension-1,1);
3106 // m_SecondRegion.SetIndex(InputDimension-1,0);
3109 // m_BC2Regions.resize(0);
3116 // DD("supportindex > 4 ???");
3117 // DD(m_SupportRegion.GetIndex(InputDimension-1));
3118 // DD(m_TransformShape);
3120 // }// end swith index
3123 // } // end case 1 shape
3127 // ----------------------------------------------------------------------
3128 // The rabbit with 4 internal CP
3129 // Periodic, constrained to zero at the reference
3130 // at position 3 and the extremes fixed with anit-symm bc
3131 // Coeff row BC1 0 1 BC2 2 BC3 BC4
3132 // PaddedCoeff R: 0 1 2 3 4 5 6
3134 // BC2= -weights[2]/weights[0] ( R2+R4 )
3137 // ---------------------------------------------------------------------
3138 switch(m_SupportRegion.GetIndex(InputDimension-1))
3143 m_FirstRegion.SetSize(InputDimension-1,0);
3146 m_BCRegions.resize(2);
3147 m_BCValues.resize(2);
3148 m_BCRegion.SetIndex(InputDimension-1,0);
3149 m_BCRegions[0]=m_BCRegion;
3151 m_BCRegion.SetIndex(InputDimension-1,1);
3152 m_BCRegions[1]=m_BCRegion;
3156 m_SecondRegion.SetSize(InputDimension-1,2);
3157 m_SecondRegion.SetIndex(InputDimension-1,0);
3160 m_BC2Regions.resize(2);
3161 m_BC2Values.resize(2);
3162 m_BCRegion.SetIndex(InputDimension-1,1);
3163 m_BC2Regions[0]=m_BCRegion;
3164 m_BC2Values[0]=-m_WeightRatio[2][0];
3165 m_BCRegion.SetIndex(InputDimension-1,2);
3166 m_BC2Regions[1]=m_BCRegion;
3167 m_BC2Values[1]= -m_WeightRatio[2][0];
3174 m_FirstRegion.SetSize(InputDimension-1,2);
3175 m_FirstRegion.SetIndex(InputDimension-1,0);
3178 m_BCRegions.resize(2);
3179 m_BCValues.resize(2);
3180 m_BCRegion.SetIndex(InputDimension-1,1);
3181 m_BCRegions[0]=m_BCRegion;
3182 m_BCValues[0]=-m_WeightRatio[2][0];
3183 m_BCRegion.SetIndex(InputDimension-1,2);
3184 m_BCRegions[1]=m_BCRegion;
3185 m_BCValues[1]= -m_WeightRatio[2][0];
3188 m_SecondRegion.SetSize(InputDimension-1,1);
3189 m_SecondRegion.SetIndex(InputDimension-1,2);
3192 m_BC2Regions.resize(0);
3199 m_FirstRegion.SetSize(InputDimension-1,1);
3200 m_FirstRegion.SetIndex(InputDimension-1,1);
3203 m_BCRegions.resize(2);
3204 m_BCValues.resize(2);
3205 m_BCRegion.SetIndex(InputDimension-1,1);
3206 m_BCRegions[0]=m_BCRegion;
3207 m_BCValues[0]=-m_WeightRatio[2][0];
3208 m_BCRegion.SetIndex(InputDimension-1,2);
3209 m_BCRegions[1]=m_BCRegion;
3210 m_BCValues[1]= -m_WeightRatio[2][0];
3213 m_SecondRegion.SetSize(InputDimension-1,1);
3214 m_SecondRegion.SetIndex(InputDimension-1,2);
3217 m_BC2Regions.resize(1);
3218 m_BC2Values.resize(1);
3219 m_BCRegion.SetIndex(InputDimension-1,0);
3220 m_BC2Regions[0]=m_BCRegion;
3228 m_FirstRegion.SetSize(InputDimension-1,0);
3231 m_BCRegions.resize(2);
3232 m_BCValues.resize(2);
3233 m_BCRegion.SetIndex(InputDimension-1,1);
3234 m_BCRegions[0]=m_BCRegion;
3235 m_BCValues[0]=-m_WeightRatio[2][0];
3236 m_BCRegion.SetIndex(InputDimension-1,2);
3237 m_BCRegions[1]=m_BCRegion;
3238 m_BCValues[1]= -m_WeightRatio[2][0];
3241 m_SecondRegion.SetSize(InputDimension-1,1);
3242 m_SecondRegion.SetIndex(InputDimension-1,2);
3245 m_BC2Regions.resize(1);
3246 m_BC2Values.resize(1);
3247 m_BCRegion.SetIndex(InputDimension-1,0);
3248 m_BC2Regions[0]=m_BCRegion;
3252 m_BC3Regions.resize(2);
3253 m_BC3Values.resize(2);
3254 m_BCRegion.SetIndex(InputDimension-1,0);
3255 m_BC3Regions[0]=m_BCRegion;
3257 m_BCRegion.SetIndex(InputDimension-1,2);
3258 m_BC3Regions[1]=m_BCRegion;
3266 DD("supportindex > 3 ???");
3267 DD(m_SupportRegion.GetIndex(InputDimension-1));
3269 } // end switch index
3272 } // end case 2 shape
3275 // ----------------------------------------------------------------------
3276 // The rabbit with 5 internal CP
3277 // Periodic, constrained to zero at the reference at position 3.5
3278 // and the extremes fixed with anti-symmetrical BC
3279 // Coeff row BC1 0 1 BC2 2 3 BC3 BC4
3280 // PaddedCoeff R: 0 1 2 3 4 5 6 7
3282 // BC2= -weights[3]/weights[1] ( R2+R5 ) - R4
3285 // ---------------------------------------------------------------------
3286 switch(m_SupportRegion.GetIndex(InputDimension-1))
3291 m_FirstRegion.SetSize(InputDimension-1,0);
3294 m_BCRegions.resize(2);
3295 m_BCValues.resize(2);
3296 m_BCRegion.SetIndex(InputDimension-1,0);
3297 m_BCRegions[0]=m_BCRegion;
3299 m_BCRegion.SetIndex(InputDimension-1,1);
3300 m_BCRegions[1]=m_BCRegion;
3304 m_SecondRegion.SetSize(InputDimension-1,2);
3305 m_SecondRegion.SetIndex(InputDimension-1,0);
3308 m_BC2Regions.resize(3);
3309 m_BC2Values.resize(3);
3310 m_BCRegion.SetIndex(InputDimension-1,1);
3311 m_BC2Regions[0]=m_BCRegion;
3312 m_BC2Values[0]=-m_WeightRatio[3][1];
3313 m_BCRegion.SetIndex(InputDimension-1,2);
3314 m_BC2Regions[1]=m_BCRegion;
3316 m_BCRegion.SetIndex(InputDimension-1,3);
3317 m_BC2Regions[2]=m_BCRegion;
3318 m_BC2Values[2]=-m_WeightRatio[3][1];
3325 m_FirstRegion.SetSize(InputDimension-1,2);
3326 m_FirstRegion.SetIndex(InputDimension-1,0);
3329 m_BCRegions.resize(3);
3330 m_BCValues.resize(3);
3331 m_BCRegion.SetIndex(InputDimension-1,1);
3332 m_BCRegions[0]=m_BCRegion;
3333 m_BCValues[0]=-m_WeightRatio[3][1];
3334 m_BCRegion.SetIndex(InputDimension-1,2);
3335 m_BCRegions[1]=m_BCRegion;
3337 m_BCRegion.SetIndex(InputDimension-1,3);
3338 m_BCRegions[2]=m_BCRegion;
3339 m_BCValues[2]=-m_WeightRatio[3][1];
3342 m_SecondRegion.SetSize(InputDimension-1,1);
3343 m_SecondRegion.SetIndex(InputDimension-1,2);
3346 m_BC2Regions.resize(0);
3353 m_FirstRegion.SetSize(InputDimension-1,1);
3354 m_FirstRegion.SetIndex(InputDimension-1,1);
3357 m_BCRegions.resize(3);
3358 m_BCValues.resize(3);
3359 m_BCRegion.SetIndex(InputDimension-1,1);
3360 m_BCRegions[0]=m_BCRegion;
3361 m_BCValues[0]=-m_WeightRatio[3][1];
3362 m_BCRegion.SetIndex(InputDimension-1,2);
3363 m_BCRegions[1]=m_BCRegion;
3365 m_BCRegion.SetIndex(InputDimension-1,3);
3366 m_BCRegions[2]=m_BCRegion;
3367 m_BCValues[2]=-m_WeightRatio[3][1];
3370 m_SecondRegion.SetSize(InputDimension-1,2);
3371 m_SecondRegion.SetIndex(InputDimension-1,2);
3374 m_BC2Regions.resize(0);
3381 m_FirstRegion.SetSize(InputDimension-1,0);
3384 m_BCRegions.resize(3);
3385 m_BCValues.resize(3);
3386 m_BCRegion.SetIndex(InputDimension-1,1);
3387 m_BCRegions[0]=m_BCRegion;
3388 m_BCValues[0]=-m_WeightRatio[3][1];
3389 m_BCRegion.SetIndex(InputDimension-1,2);
3390 m_BCRegions[1]=m_BCRegion;
3392 m_BCRegion.SetIndex(InputDimension-1,3);
3393 m_BCRegions[2]=m_BCRegion;
3394 m_BCValues[2]=-m_WeightRatio[3][1];
3397 m_SecondRegion.SetSize(InputDimension-1,2);
3398 m_SecondRegion.SetIndex(InputDimension-1,2);
3401 m_BC2Regions.resize(1);
3402 m_BC2Values.resize(1);
3403 m_BCRegion.SetIndex(InputDimension-1,0);
3404 m_BC2Regions[0]=m_BCRegion;
3413 m_FirstRegion.SetSize(InputDimension-1,2);
3414 m_FirstRegion.SetIndex(InputDimension-1,2);
3417 m_BCRegions.resize(1);
3418 m_BCValues.resize(1);
3419 m_BCRegion.SetIndex(InputDimension-1,0);
3420 m_BCRegions[0]=m_BCRegion;
3424 m_SecondRegion.SetSize(InputDimension-1,0);
3427 m_BC2Regions.resize(2);
3428 m_BC2Values.resize(2);
3429 m_BCRegion.SetIndex(InputDimension-1,0);
3430 m_BC2Regions[0]=m_BCRegion;
3432 m_BCRegion.SetIndex(InputDimension-1,3);
3433 m_BC2Regions[1]=m_BCRegion;
3441 DD("supportindex > 4 ???");
3442 DD(m_SupportRegion.GetIndex(InputDimension-1));
3444 } // end switch index
3447 } // end case 3 shape
3451 // ----------------------------------------------------------------------
3452 // The sputnik with 4 internal CP
3453 // Periodic, constrained to zero at the reference
3454 // at position 3 and one indepent extremes copied
3455 // Coeff row BC1 0 1 BC2 2 BC2 3
3456 // PaddedCoeff R: 0 1 2 3 4 5 6
3458 // BC2= -weights[2]/weights[0] ( R2+R4 )
3459 // BC3= weights[2]/weights[0] ( R2-R4) + R1
3460 // ---------------------------------------------------------------------
3461 switch(m_SupportRegion.GetIndex(InputDimension-1))
3466 m_FirstRegion.SetSize(InputDimension-1,0);
3469 m_BCRegions.resize(1);
3470 m_BCValues.resize(1);
3471 m_BCRegion.SetIndex(InputDimension-1,3);
3472 m_BCRegions[0]=m_BCRegion;
3476 m_SecondRegion.SetSize(InputDimension-1,2);
3477 m_SecondRegion.SetIndex(InputDimension-1,0);
3480 m_BC2Regions.resize(2);
3481 m_BC2Values.resize(2);
3482 m_BCRegion.SetIndex(InputDimension-1,1);
3483 m_BC2Regions[0]=m_BCRegion;
3484 m_BC2Values[0]=-m_WeightRatio[2][0];
3485 m_BCRegion.SetIndex(InputDimension-1,2);
3486 m_BC2Regions[1]=m_BCRegion;
3487 m_BC2Values[1]= -m_WeightRatio[2][0];
3494 m_FirstRegion.SetSize(InputDimension-1,2);
3495 m_FirstRegion.SetIndex(InputDimension-1,0);
3498 m_BCRegions.resize(2);
3499 m_BCValues.resize(2);
3500 m_BCRegion.SetIndex(InputDimension-1,1);
3501 m_BCRegions[0]=m_BCRegion;
3502 m_BCValues[0]=-m_WeightRatio[2][0];
3503 m_BCRegion.SetIndex(InputDimension-1,2);
3504 m_BCRegions[1]=m_BCRegion;
3505 m_BCValues[1]= -m_WeightRatio[2][0];
3508 m_SecondRegion.SetSize(InputDimension-1,1);
3509 m_SecondRegion.SetIndex(InputDimension-1,2);
3512 m_BC2Regions.resize(0);
3519 m_FirstRegion.SetSize(InputDimension-1,1);
3520 m_FirstRegion.SetIndex(InputDimension-1,1);
3523 m_BCRegions.resize(2);
3524 m_BCValues.resize(2);
3525 m_BCRegion.SetIndex(InputDimension-1,1);
3526 m_BCRegions[0]=m_BCRegion;
3527 m_BCValues[0]=-m_WeightRatio[2][0];
3528 m_BCRegion.SetIndex(InputDimension-1,2);
3529 m_BCRegions[1]=m_BCRegion;
3530 m_BCValues[1]= -m_WeightRatio[2][0];
3533 m_SecondRegion.SetSize(InputDimension-1,1);
3534 m_SecondRegion.SetIndex(InputDimension-1,2);
3537 m_BC2Regions.resize(3);
3538 m_BC2Values.resize(3);
3539 m_BCRegion.SetIndex(InputDimension-1,0);
3540 m_BC2Regions[0]=m_BCRegion;
3542 m_BCRegion.SetIndex(InputDimension-1,1);
3543 m_BC2Regions[1]=m_BCRegion;
3544 m_BC2Values[1]=m_WeightRatio[2][0];
3545 m_BCRegion.SetIndex(InputDimension-1,2);
3546 m_BC2Regions[2]=m_BCRegion;
3547 m_BC2Values[2]=-m_WeightRatio[2][0];
3554 m_FirstRegion.SetSize(InputDimension-1,0);
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];
3584 m_ThirdRegion.SetSize(InputDimension-1,1);
3585 m_ThirdRegion.SetIndex(InputDimension-1,3);
3592 DD("supportindex > 3 ???");
3593 DD(m_SupportRegion.GetIndex(InputDimension-1));
3595 } // end switch index
3598 } // end case 4 shape
3602 // ----------------------------------------------------------------------
3603 // The sputnik with 5 internal CP
3604 // Periodic, constrained to zero at the reference
3605 // at position 3 and one indepent extreme
3606 // Coeff row BC1 0 1 BC2 2 3 BC3 4
3607 // PaddedCoeff R: 0 1 2 3 4 5 6 7
3608 // BC1= R2 + R5 - R7
3609 // BC2= -weights[3]/weights[1] ( R2+R5 ) - R4
3610 // BC3= R1 + 0.5 R2 - 0.5 R7
3611 // ----------------------------------------------------------------------
3612 switch(m_SupportRegion.GetIndex(InputDimension-1))
3617 m_FirstRegion.SetSize(InputDimension-1,0);
3620 m_BCRegions.resize(3);
3621 m_BCValues.resize(3);
3622 m_BCRegion.SetIndex(InputDimension-1,1);
3623 m_BCRegions[0]=m_BCRegion;
3625 m_BCRegion.SetIndex(InputDimension-1,3);
3626 m_BCRegions[1]=m_BCRegion;
3628 m_BCRegion.SetIndex(InputDimension-1,4);
3629 m_BCRegions[2]=m_BCRegion;
3633 m_SecondRegion.SetSize(InputDimension-1,2);
3634 m_SecondRegion.SetIndex(InputDimension-1,0);
3637 m_BC2Regions.resize(3);
3638 m_BC2Values.resize(3);
3639 m_BCRegion.SetIndex(InputDimension-1,1);
3640 m_BC2Regions[0]=m_BCRegion;
3641 m_BC2Values[0]=-m_WeightRatio[3][1];
3642 m_BCRegion.SetIndex(InputDimension-1,2);
3643 m_BC2Regions[1]=m_BCRegion;
3645 m_BCRegion.SetIndex(InputDimension-1,3);
3646 m_BC2Regions[2]=m_BCRegion;
3647 m_BC2Values[2]=-m_WeightRatio[3][1];
3654 m_FirstRegion.SetSize(InputDimension-1,2);
3655 m_FirstRegion.SetIndex(InputDimension-1,0);
3658 m_BCRegions.resize(3);
3659 m_BCValues.resize(3);
3660 m_BCRegion.SetIndex(InputDimension-1,1);
3661 m_BCRegions[0]=m_BCRegion;
3662 m_BCValues[0]=-m_WeightRatio[3][1];
3663 m_BCRegion.SetIndex(InputDimension-1,2);
3664 m_BCRegions[1]=m_BCRegion;
3666 m_BCRegion.SetIndex(InputDimension-1,3);
3667 m_BCRegions[2]=m_BCRegion;
3668 m_BCValues[2]=-m_WeightRatio[3][1];
3671 m_SecondRegion.SetSize(InputDimension-1,1);
3672 m_SecondRegion.SetIndex(InputDimension-1,2);
3675 m_BC2Regions.resize(0);
3682 m_FirstRegion.SetSize(InputDimension-1,1);
3683 m_FirstRegion.SetIndex(InputDimension-1,1);
3686 m_BCRegions.resize(3);
3687 m_BCValues.resize(3);
3688 m_BCRegion.SetIndex(InputDimension-1,1);
3689 m_BCRegions[0]=m_BCRegion;
3690 m_BCValues[0]=-m_WeightRatio[3][1];
3691 m_BCRegion.SetIndex(InputDimension-1,2);
3692 m_BCRegions[1]=m_BCRegion;
3694 m_BCRegion.SetIndex(InputDimension-1,3);
3695 m_BCRegions[2]=m_BCRegion;
3696 m_BCValues[2]=-m_WeightRatio[3][1];
3699 m_SecondRegion.SetSize(InputDimension-1,2);
3700 m_SecondRegion.SetIndex(InputDimension-1,2);
3703 m_BC2Regions.resize(0);
3710 m_FirstRegion.SetSize(InputDimension-1,0);
3713 m_BCRegions.resize(3);
3714 m_BCValues.resize(3);
3715 m_BCRegion.SetIndex(InputDimension-1,1);
3716 m_BCRegions[0]=m_BCRegion;
3717 m_BCValues[0]=-m_WeightRatio[3][1];
3718 m_BCRegion.SetIndex(InputDimension-1,2);
3719 m_BCRegions[1]=m_BCRegion;
3721 m_BCRegion.SetIndex(InputDimension-1,3);
3722 m_BCRegions[2]=m_BCRegion;
3723 m_BCValues[2]=-m_WeightRatio[3][1];
3726 m_SecondRegion.SetSize(InputDimension-1,2);
3727 m_SecondRegion.SetIndex(InputDimension-1,2);
3730 m_BC2Regions.resize(3);
3731 m_BC2Values.resize(3);
3732 m_BCRegion.SetIndex(InputDimension-1,0);
3733 m_BC2Regions[0]=m_BCRegion;
3735 m_BCRegion.SetIndex(InputDimension-1,1);
3736 m_BC2Regions[1]=m_BCRegion;
3738 m_BCRegion.SetIndex(InputDimension-1,4);
3739 m_BC2Regions[2]=m_BCRegion;
3740 m_BC2Values[2]=-0.5;
3747 m_FirstRegion.SetSize(InputDimension-1,2);
3748 m_FirstRegion.SetIndex(InputDimension-1,2);
3751 m_BCRegions.resize(3);
3752 m_BCValues.resize(3);
3753 m_BCRegion.SetIndex(InputDimension-1,0);
3754 m_BCRegions[0]=m_BCRegion;
3756 m_BCRegion.SetIndex(InputDimension-1,1);
3757 m_BCRegions[1]=m_BCRegion;
3759 m_BCRegion.SetIndex(InputDimension-1,4);
3760 m_BCRegions[2]=m_BCRegion;
3764 m_SecondRegion.SetSize(InputDimension-1,1);
3765 m_SecondRegion.SetIndex(InputDimension-1,4);
3768 m_BC2Regions.resize(0);
3774 DD("supportindex > 4 ???");
3775 DD(m_SupportRegion.GetIndex(InputDimension-1));
3777 } // end switch index
3780 } // end case 5 shape
3787 // ----------------------------------------------------------------------
3788 // The diamond with 4 internal CP:
3789 // Periodic, constrained to zero at the reference at position 3
3790 // Coeff row 0 1 2 BC1 3 BC2 4
3791 // PaddedCoeff R:0 1 2 3 4 5 6
3792 // BC1= -weights[2]/weights[0] ( R2+R4 )
3793 // BC2= weights[2]/weights[0] ( R0+R2-R4-R6 ) + R1
3794 // ---------------------------------------------------------------------
3795 switch(m_SupportRegion.GetIndex(InputDimension-1))
3800 m_FirstRegion.SetSize(InputDimension-1,3);
3801 m_FirstRegion.SetIndex(InputDimension-1,0);
3804 m_BCRegions.resize(2);
3805 m_BCValues.resize(2);
3806 m_BCRegion.SetIndex(InputDimension-1,2);
3807 m_BCRegions[0]=m_BCRegion;
3808 m_BCValues[0]=-m_WeightRatio[2][0];
3809 m_BCRegion.SetIndex(InputDimension-1,3);
3810 m_BCRegions[1]=m_BCRegion;
3811 m_BCValues[1]=-m_WeightRatio[2][0];
3814 m_SecondRegion.SetSize(InputDimension-1,0);
3817 m_BC2Regions.resize(0);
3824 m_FirstRegion.SetSize(InputDimension-1,2);
3825 m_FirstRegion.SetIndex(InputDimension-1,1);
3828 m_BCRegions.resize(2);
3829 m_BCValues.resize(2);
3830 m_BCRegion.SetIndex(InputDimension-1,2);
3831 m_BCRegions[0]=m_BCRegion;
3832 m_BCValues[0]=-m_WeightRatio[2][0];
3833 m_BCRegion.SetIndex(InputDimension-1,3);
3834 m_BCRegions[1]=m_BCRegion;
3835 m_BCValues[1]=-m_WeightRatio[2][0];
3838 m_SecondRegion.SetSize(InputDimension-1,1);
3839 m_SecondRegion.SetIndex(InputDimension-1,3);
3842 m_BC2Regions.resize(0);
3849 m_FirstRegion.SetSize(InputDimension-1,1);
3850 m_FirstRegion.SetIndex(InputDimension-1,2);
3853 m_BCRegions.resize(2);
3854 m_BCValues.resize(2);
3855 m_BCRegion.SetIndex(InputDimension-1,2);
3856 m_BCRegions[0]=m_BCRegion;
3857 m_BCValues[0]=-m_WeightRatio[2][0];
3858 m_BCRegion.SetIndex(InputDimension-1,3);
3859 m_BCRegions[1]=m_BCRegion;
3860 m_BCValues[1]=-m_WeightRatio[2][0];
3863 m_SecondRegion.SetSize(InputDimension-1,1);
3864 m_SecondRegion.SetIndex(InputDimension-1,3);
3867 m_BC2Regions.resize(5);
3868 m_BC2Values.resize(5);
3869 m_BCRegion.SetIndex(InputDimension-1,0);
3870 m_BC2Regions[0]=m_BCRegion;
3871 m_BC2Values[0]=m_WeightRatio[2][0];
3872 m_BCRegion.SetIndex(InputDimension-1,1);
3873 m_BC2Regions[1]=m_BCRegion;
3875 m_BCRegion.SetIndex(InputDimension-1,2);
3876 m_BC2Regions[2]=m_BCRegion;
3877 m_BC2Values[2]=m_WeightRatio[2][0];
3878 m_BCRegion.SetIndex(InputDimension-1,3);
3879 m_BC2Regions[3]=m_BCRegion;
3880 m_BC2Values[3]=-m_WeightRatio[2][0];
3881 m_BCRegion.SetIndex(InputDimension-1,4);
3882 m_BC2Regions[4]=m_BCRegion;
3883 m_BC2Values[4]=-m_WeightRatio[2][0];
3890 m_FirstRegion.SetSize(InputDimension-1,0);
3893 m_BCRegions.resize(2);
3894 m_BCValues.resize(2);
3895 m_BCRegion.SetIndex(InputDimension-1,2);
3896 m_BCRegions[0]=m_BCRegion;
3897 m_BCValues[0]=-m_WeightRatio[2][0];
3898 m_BCRegion.SetIndex(InputDimension-1,3);
3899 m_BCRegions[1]=m_BCRegion;
3900 m_BCValues[1]=-m_WeightRatio[2][0];
3903 m_SecondRegion.SetSize(InputDimension-1,1);
3904 m_SecondRegion.SetIndex(InputDimension-1,3);
3907 m_BC2Regions.resize(5);
3908 m_BC2Values.resize(5);
3909 m_BCRegion.SetIndex(InputDimension-1,0);
3910 m_BC2Regions[0]=m_BCRegion;
3911 m_BC2Values[0]=m_WeightRatio[2][0];
3912 m_BCRegion.SetIndex(InputDimension-1,1);
3913 m_BC2Regions[1]=m_BCRegion;
3915 m_BCRegion.SetIndex(InputDimension-1,2);
3916 m_BC2Regions[2]=m_BCRegion;
3917 m_BC2Values[2]=m_WeightRatio[2][0];
3918 m_BCRegion.SetIndex(InputDimension-1,3);
3919 m_BC2Regions[3]=m_BCRegion;
3920 m_BC2Values[3]=-m_WeightRatio[2][0];
3921 m_BCRegion.SetIndex(InputDimension-1,4);
3922 m_BC2Regions[4]=m_BCRegion;
3923 m_BC2Values[4]=-m_WeightRatio[2][0];
3926 m_ThirdRegion.SetSize(InputDimension-1,1);
3927 m_ThirdRegion.SetIndex(InputDimension-1,4);
3934 DD("supportindex > 3 ???");
3935 DD(m_SupportRegion.GetIndex(InputDimension-1));
3937 } // end switch index
3940 } // end case 7 shape
3944 // ----------------------------------------------------------------------
3945 // The diamond with 5 internal CP:
3946 // periodic, constrained to zero at the reference at position 3.5
3947 // Coeff row 0 1 2 BC1 3 4 BC2 5
3948 // PaddedCoeff R:0 1 2 3 4 5 6 7
3949 // BC1= -weights[3]/weights[1] ( R2+R5 ) - R4
3950 // BC2= weights[2]/weights[0] ( R0+R2-R5-R7 ) + R1
3951 // ---------------------------------------------------------------------
3952 switch(m_SupportRegion.GetIndex(InputDimension-1))
3957 m_FirstRegion.SetSize(InputDimension-1,3);
3958 m_FirstRegion.SetIndex(InputDimension-1,0);
3961 m_BCRegions.resize(3);
3962 m_BCValues.resize(3);
3963 m_BCRegion.SetIndex(InputDimension-1,2);
3964 m_BCRegions[0]=m_BCRegion;
3965 m_BCValues[0]=-m_WeightRatio[3][1];
3966 m_BCRegion.SetIndex(InputDimension-1,3);
3967 m_BCRegions[1]=m_BCRegion;
3969 m_BCRegion.SetIndex(InputDimension-1,4);
3970 m_BCRegions[2]=m_BCRegion;
3971 m_BCValues[2]=-m_WeightRatio[3][1];
3974 m_SecondRegion.SetSize(InputDimension-1,0);
3977 m_BC2Regions.resize(0);
3984 m_FirstRegion.SetSize(InputDimension-1,2);
3985 m_FirstRegion.SetIndex(InputDimension-1,1);
3988 m_BCRegions.resize(3);
3989 m_BCValues.resize(3);
3990 m_BCRegion.SetIndex(InputDimension-1,2);
3991 m_BCRegions[0]=m_BCRegion;
3992 m_BCValues[0]=-m_WeightRatio[3][1];
3993 m_BCRegion.SetIndex(InputDimension-1,3);
3994 m_BCRegions[1]=m_BCRegion;
3996 m_BCRegion.SetIndex(InputDimension-1,4);
3997 m_BCRegions[2]=m_BCRegion;
3998 m_BCValues[2]=-m_WeightRatio[3][1];
4001 m_SecondRegion.SetSize(InputDimension-1,1);
4002 m_SecondRegion.SetIndex(InputDimension-1,3);
4005 m_BC2Regions.resize(0);
4012 m_FirstRegion.SetSize(InputDimension-1,1);
4013 m_FirstRegion.SetIndex(InputDimension-1,2);
4016 m_BCRegions.resize(3);
4017 m_BCValues.resize(3);
4018 m_BCRegion.SetIndex(InputDimension-1,2);
4019 m_BCRegions[0]=m_BCRegion;
4020 m_BCValues[0]=-m_WeightRatio[3][1];
4021 m_BCRegion.SetIndex(InputDimension-1,3);
4022 m_BCRegions[1]=m_BCRegion;
4024 m_BCRegion.SetIndex(InputDimension-1,4);
4025 m_BCRegions[2]=m_BCRegion;
4026 m_BCValues[2]=-m_WeightRatio[3][1];
4029 m_SecondRegion.SetSize(InputDimension-1,2);
4030 m_SecondRegion.SetIndex(InputDimension-1,3);
4033 m_BC2Regions.resize(0);
4040 m_FirstRegion.SetSize(InputDimension-1,0);
4041 m_FirstRegion.SetIndex(InputDimension-1,0);
4044 m_BCRegions.resize(3);
4045 m_BCValues.resize(3);
4046 m_BCRegion.SetIndex(InputDimension-1,2);
4047 m_BCRegions[0]=m_BCRegion;
4048 m_BCValues[0]=-m_WeightRatio[3][1];
4049 m_BCRegion.SetIndex(InputDimension-1,3);
4050 m_BCRegions[1]=m_BCRegion;
4052 m_BCRegion.SetIndex(InputDimension-1,4);
4053 m_BCRegions[2]=m_BCRegion;
4054 m_BCValues[2]=-m_WeightRatio[3][1];
4057 m_SecondRegion.SetSize(InputDimension-1,2);
4058 m_SecondRegion.SetIndex(InputDimension-1,3);
4061 m_BC2Regions.resize(5);
4062 m_BC2Values.resize(5);
4063 m_BCRegion.SetIndex(InputDimension-1,0);
4064 m_BC2Regions[0]=m_BCRegion;
4065 m_BC2Values[0]=m_WeightRatio[2][0];
4066 m_BCRegion.SetIndex(InputDimension-1,1);
4067 m_BC2Regions[1]=m_BCRegion;
4069 m_BCRegion.SetIndex(InputDimension-1,2);
4070 m_BC2Regions[2]=m_BCRegion;
4071 m_BC2Values[2]=m_WeightRatio[2][0];
4072 m_BCRegion.SetIndex(InputDimension-1,4);
4073 m_BC2Regions[3]=m_BCRegion;
4074 m_BC2Values[3]=-m_WeightRatio[2][0];
4075 m_BCRegion.SetIndex(InputDimension-1,5);
4076 m_BC2Regions[4]=m_BCRegion;
4077 m_BC2Values[4]=-m_WeightRatio[2][0];
4085 m_FirstRegion.SetSize(InputDimension-1,2);
4086 m_FirstRegion.SetIndex(InputDimension-1,3);
4089 m_BCRegions.resize(5);
4090 m_BCValues.resize(5);
4091 m_BCRegion.SetIndex(InputDimension-1,0);
4092 m_BCRegions[0]=m_BCRegion;
4093 m_BCValues[0]=m_WeightRatio[2][0];
4094 m_BCRegion.SetIndex(InputDimension-1,1);
4095 m_BCRegions[1]=m_BCRegion;
4097 m_BCRegion.SetIndex(InputDimension-1,2);
4098 m_BCRegions[2]=m_BCRegion;
4099 m_BCValues[2]=m_WeightRatio[2][0];
4100 m_BCRegion.SetIndex(InputDimension-1,4);
4101 m_BCRegions[3]=m_BCRegion;
4102 m_BCValues[3]=-m_WeightRatio[2][0];
4103 m_BCRegion.SetIndex(InputDimension-1,5);
4104 m_BCRegions[4]=m_BCRegion;
4105 m_BCValues[4]=-m_WeightRatio[2][0];
4108 m_SecondRegion.SetSize(InputDimension-1,1);
4109 m_SecondRegion.SetIndex(InputDimension-1,5);
4112 m_BC2Regions.resize(0);
4119 DD("supportindex > 4 ???");
4120 DD(m_SupportRegion.GetIndex(InputDimension-1));
4122 } // end switch index
4125 } // end case 7 shape
4129 // ----------------------------------------------------------------------
4130 // The sputnik with 5 internal CP
4131 // Periodic, constrained to zero at the reference
4132 // at position 3 and one indepent extreme
4133 // Coeff row BC1 0 1 BC2 2 3 BC3 4
4134 // PaddedCoeff R: 0 1 2 3 4 5 6 7
4136 // BC2= -weights[3]/weights[1] ( R2+R5 ) - R4
4138 // ---------------------------------------------------------------------
4139 switch(m_SupportRegion.GetIndex(InputDimension-1))
4144 m_FirstRegion.SetSize(InputDimension-1,0);
4147 m_BCRegions.resize(3);
4148 m_BCValues.resize(3);
4149 m_BCRegion.SetIndex(InputDimension-1,1);
4150 m_BCRegions[0]=m_BCRegion;
4152 m_BCRegion.SetIndex(InputDimension-1,3);
4153 m_BCRegions[1]=m_BCRegion;
4155 m_BCRegion.SetIndex(InputDimension-1,4);
4156 m_BCRegions[2]=m_BCRegion;
4160 m_SecondRegion.SetSize(InputDimension-1,2);
4161 m_SecondRegion.SetIndex(InputDimension-1,0);
4164 m_BC2Regions.resize(3);
4165 m_BC2Values.resize(3);
4166 m_BCRegion.SetIndex(InputDimension-1,1);
4167 m_BC2Regions[0]=m_BCRegion;
4168 m_BC2Values[0]=-m_WeightRatio[3][1];
4169 m_BCRegion.SetIndex(InputDimension-1,2);
4170 m_BC2Regions[1]=m_BCRegion;
4172 m_BCRegion.SetIndex(InputDimension-1,3);
4173 m_BC2Regions[2]=m_BCRegion;
4174 m_BC2Values[2]=-m_WeightRatio[3][1];
4181 m_FirstRegion.SetSize(InputDimension-1,2);
4182 m_FirstRegion.SetIndex(InputDimension-1,0);
4185 m_BCRegions.resize(3);
4186 m_BCValues.resize(3);
4187 m_BCRegion.SetIndex(InputDimension-1,1);
4188 m_BCRegions[0]=m_BCRegion;
4189 m_BCValues[0]=-m_WeightRatio[3][1];
4190 m_BCRegion.SetIndex(InputDimension-1,2);
4191 m_BCRegions[1]=m_BCRegion;
4193 m_BCRegion.SetIndex(InputDimension-1,3);
4194 m_BCRegions[2]=m_BCRegion;
4195 m_BCValues[2]=-m_WeightRatio[3][1];
4198 m_SecondRegion.SetSize(InputDimension-1,1);
4199 m_SecondRegion.SetIndex(InputDimension-1,2);
4202 m_BC2Regions.resize(0);
4209 m_FirstRegion.SetSize(InputDimension-1,1);
4210 m_FirstRegion.SetIndex(InputDimension-1,1);
4213 m_BCRegions.resize(3);
4214 m_BCValues.resize(3);
4215 m_BCRegion.SetIndex(InputDimension-1,1);
4216 m_BCRegions[0]=m_BCRegion;
4217 m_BCValues[0]=-m_WeightRatio[3][1];
4218 m_BCRegion.SetIndex(InputDimension-1,2);
4219 m_BCRegions[1]=m_BCRegion;
4221 m_BCRegion.SetIndex(InputDimension-1,3);
4222 m_BCRegions[2]=m_BCRegion;
4223 m_BCValues[2]=-m_WeightRatio[3][1];
4226 m_SecondRegion.SetSize(InputDimension-1,2);
4227 m_SecondRegion.SetIndex(InputDimension-1,2);
4230 m_BC2Regions.resize(0);
4237 m_FirstRegion.SetSize(InputDimension-1,0);
4240 m_BCRegions.resize(3);
4241 m_BCValues.resize(3);
4242 m_BCRegion.SetIndex(InputDimension-1,1);
4243 m_BCRegions[0]=m_BCRegion;
4244 m_BCValues[0]=-m_WeightRatio[3][1];
4245 m_BCRegion.SetIndex(InputDimension-1,2);
4246 m_BCRegions[1]=m_BCRegion;
4248 m_BCRegion.SetIndex(InputDimension-1,3);
4249 m_BCRegions[2]=m_BCRegion;
4250 m_BCValues[2]=-m_WeightRatio[3][1];
4253 m_SecondRegion.SetSize(InputDimension-1,1);
4254 m_SecondRegion.SetIndex(InputDimension-1,2);
4257 m_BC2Regions.resize(1);
4258 m_BC2Values.resize(1);
4259 m_BCRegion.SetIndex(InputDimension-1,0);
4260 m_BC2Regions[0]=m_BCRegion;
4268 m_FirstRegion.SetSize(InputDimension-1,2);
4269 m_FirstRegion.SetIndex(InputDimension-1,2);
4272 m_BCRegions.resize(1);
4273 m_BCValues.resize(1);
4274 m_BCRegion.SetIndex(InputDimension-1,0);
4275 m_BCRegions[0]=m_BCRegion;
4279 m_SecondRegion.SetSize(InputDimension-1,1);
4280 m_SecondRegion.SetIndex(InputDimension-1,4);
4283 m_BC2Regions.resize(0);
4289 DD("supportindex > 4 ???");
4290 DD(m_SupportRegion.GetIndex(InputDimension-1));
4292 } // end switch index
4295 } // end case 9 shape
4300 DD ("Other shapes currently not implemented");
4303 } // end switch shape
4304 } // end wrap region
4308 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
4310 ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions,NOutputDimensions>
4311 ::TransformPointToContinuousIndex( const InputPointType & point, ContinuousIndexType & index ) const
4315 itk::Vector<double, OutputDimension> tvector;
4317 for ( j = 0; j < OutputDimension; j++ )
4319 //JV find the index in the PADDED version
4320 tvector[j] = point[j] - this->m_PaddedGridOrigin[j];
4323 itk::Vector<double, OutputDimension> cvector;
4325 cvector = m_PointToIndex * tvector;
4327 for ( j = 0; j < OutputDimension; j++ )
4329 index[j] = static_cast< typename ContinuousIndexType::CoordRepType >( cvector[j] );