1 #ifndef _clitkVectorBSplineInterpolateImageFunction_txx
2 #define _clitkVectorBSplineInterpolateImageFunction_txx
4 // First, make sure that we include the configuration file.
5 // This line may be removed once the ThreadSafeTransform gets
6 // integrated into ITK.
7 #include "itkConfigure.h"
9 // Second, redirect to the optimized version if necessary
10 // #ifdef ITK_USE_OPTIMIZED_REGISTRATION_METHODS
11 // #include "itkOptVectorBSplineInterpolateImageFunction.txx"
14 #include "clitkVectorBSplineInterpolateImageFunction.h"
17 #include "itkImageLinearIteratorWithIndex.h"
18 #include "itkImageRegionConstIteratorWithIndex.h"
19 #include "itkImageRegionIterator.h"
21 #include "itkVector.h"
23 #include "itkMatrix.h"
31 template <class TImageType, class TCoordRep, class TCoefficientType>
32 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
33 ::VectorBSplineInterpolateImageFunction()
36 unsigned int SplineOrder = 3;
37 m_CoefficientFilter = CoefficientFilter::New();
38 // ***TODO: Should we store coefficients in a variable or retrieve from filter?
39 m_Coefficients = CoefficientImageType::New();
40 this->SetSplineOrder(SplineOrder);
41 this->m_UseImageDirection = false;
45 * Standard "PrintSelf" method
47 template <class TImageType, class TCoordRep, class TCoefficientType>
49 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
52 itk::Indent indent) const
54 Superclass::PrintSelf( os, indent );
55 os << indent << "Spline Order: " << m_SplineOrder << std::endl;
56 os << indent << "UseImageDirection = "
57 << (this->m_UseImageDirection ? "On" : "Off") << std::endl;
62 template <class TImageType, class TCoordRep, class TCoefficientType>
64 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
65 ::SetInputImage(const TImageType * inputData)
70 DD("calling decomposition filter");
71 m_CoefficientFilter->SetInput(inputData);
73 // the Coefficient Filter requires that the spline order and the input data be set.
74 // TODO: We need to ensure that this is only run once and only after both input and
75 // spline order have been set. Should we force an update after the
76 // splineOrder has been set also?
78 m_CoefficientFilter->Update();
79 m_Coefficients = m_CoefficientFilter->GetOutput();
82 // Call the Superclass implementation after, in case the filter
83 // pulls in more of the input image
84 Superclass::SetInputImage(inputData);
86 m_DataLength = inputData->GetBufferedRegion().GetSize();
91 m_Coefficients = NULL;
96 template <class TImageType, class TCoordRep, class TCoefficientType>
98 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
99 ::SetSplineOrder(unsigned int SplineOrder)
101 if (SplineOrder == m_SplineOrder)
105 m_SplineOrder = SplineOrder;
106 m_CoefficientFilter->SetSplineOrder( SplineOrder );
109 m_MaxNumberInterpolationPoints = 1;
110 for (unsigned int n=0; n < ImageDimension; n++)
112 m_MaxNumberInterpolationPoints *= ( m_SplineOrder + 1);
114 this->GeneratePointsToIndex( );
118 template <class TImageType, class TCoordRep, class TCoefficientType>
120 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
122 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
123 ::EvaluateAtContinuousIndex( const ContinuousIndexType & x ) const
125 vnl_matrix<long> EvaluateIndex(ImageDimension, ( m_SplineOrder + 1 ));
127 // compute the interpolation indexes
128 this->DetermineRegionOfSupport(EvaluateIndex, x, m_SplineOrder);
131 vnl_matrix<double> weights(ImageDimension, ( m_SplineOrder + 1 ));
132 SetInterpolationWeights( x, EvaluateIndex, weights, m_SplineOrder );
134 // Modify EvaluateIndex at the boundaries using mirror boundary conditions
135 this->ApplyMirrorBoundaryConditions(EvaluateIndex, m_SplineOrder);
137 // perform interpolation
139 itk::Vector<double, VectorDimension> interpolated;
140 for (unsigned int i=0; i< VectorDimension; i++) interpolated[i]=itk::NumericTraits<double>::Zero;
142 IndexType coefficientIndex;
143 // Step through eachpoint in the N-dimensional interpolation cube.
144 for (unsigned int p = 0; p < m_MaxNumberInterpolationPoints; p++)
146 // translate each step into the N-dimensional index.
147 // IndexType pointIndex = PointToIndex( p );
150 for (unsigned int n = 0; n < ImageDimension; n++ )
153 w *= weights[n][ m_PointsToIndex[p][n] ];
154 coefficientIndex[n] = EvaluateIndex[n][m_PointsToIndex[p][n]]; // Build up ND index for coefficients.
156 // Convert our step p to the appropriate point in ND space in the
157 // m_Coefficients cube.
158 //JV shouldn't be necessary
159 for (unsigned int i=0; i<VectorDimension; i++)
160 interpolated[i] += w * m_Coefficients->GetPixel(coefficientIndex)[i];
163 /* double interpolated = 0.0;
164 IndexType coefficientIndex;
165 // Step through eachpoint in the N-dimensional interpolation cube.
166 for (unsigned int sp = 0; sp <= m_SplineOrder; sp++)
168 for (unsigned int sp1=0; sp1 <= m_SplineOrder; sp1++)
172 for (unsigned int n1 = 0; n1 < ImageDimension; n1++ )
174 w *= weights[n1][ sp1 ];
175 coefficientIndex[n1] = EvaluateIndex[n1][sp]; // Build up ND index for coefficients.
178 interpolated += w * m_Coefficients->GetPixel(coefficientIndex);
182 return(interpolated);
187 template <class TImageType, class TCoordRep, class TCoefficientType>
189 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
190 :: CovariantVectorType
191 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
192 ::EvaluateDerivativeAtContinuousIndex( const ContinuousIndexType & x ) const
194 vnl_matrix<long> EvaluateIndex(ImageDimension, ( m_SplineOrder + 1 ));
196 // compute the interpolation indexes
197 // TODO: Do we need to revisit region of support for the derivatives?
198 this->DetermineRegionOfSupport(EvaluateIndex, x, m_SplineOrder);
201 vnl_matrix<double> weights(ImageDimension, ( m_SplineOrder + 1 ));
202 SetInterpolationWeights( x, EvaluateIndex, weights, m_SplineOrder );
204 vnl_matrix<double> weightsDerivative(ImageDimension, ( m_SplineOrder + 1));
205 SetDerivativeWeights( x, EvaluateIndex, weightsDerivative, ( m_SplineOrder ) );
207 // Modify EvaluateIndex at the boundaries using mirror boundary conditions
208 this->ApplyMirrorBoundaryConditions(EvaluateIndex, m_SplineOrder);
210 const InputImageType * inputImage = this->GetInputImage();
211 const typename InputImageType::SpacingType & spacing = inputImage->GetSpacing();
213 // Calculate derivative
214 CovariantVectorType derivativeValue;
216 IndexType coefficientIndex;
217 for (unsigned int n = 0; n < ImageDimension; n++)
219 derivativeValue[n] = 0.0;
220 for (unsigned int p = 0; p < m_MaxNumberInterpolationPoints; p++)
223 for (unsigned int n1 = 0; n1 < ImageDimension; n1++)
225 //coefficientIndex[n1] = EvaluateIndex[n1][sp];
226 coefficientIndex[n1] = EvaluateIndex[n1][m_PointsToIndex[p][n1]];
230 //w *= weights[n][ m_PointsToIndex[p][n] ];
231 tempValue *= weightsDerivative[n1][ m_PointsToIndex[p][n1] ];
235 tempValue *= weights[n1][ m_PointsToIndex[p][n1] ];
238 derivativeValue[n] += m_Coefficients->GetPixel(coefficientIndex) * tempValue ;
240 derivativeValue[n] /= spacing[n]; // take spacing into account
243 #ifdef ITK_USE_ORIENTED_IMAGE_DIRECTION
244 if( this->m_UseImageDirection )
246 CovariantVectorType orientedDerivative;
247 inputImage->TransformLocalVectorToPhysicalVector( derivativeValue, orientedDerivative );
248 return orientedDerivative;
252 return(derivativeValue);
256 template <class TImageType, class TCoordRep, class TCoefficientType>
258 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
259 ::SetInterpolationWeights( const ContinuousIndexType & x, const vnl_matrix<long> & EvaluateIndex,
260 vnl_matrix<double> & weights, unsigned int splineOrder ) const
262 // For speed improvements we could make each case a separate function and use
263 // function pointers to reference the correct weight order.
264 // Left as is for now for readability.
265 double w, w2, w4, t, t0, t1;
270 for (unsigned int n = 0; n < ImageDimension; n++)
272 w = x[n] - (double) EvaluateIndex[n][1];
273 weights[n][3] = (1.0 / 6.0) * w * w * w;
274 weights[n][0] = (1.0 / 6.0) + 0.5 * w * (w - 1.0) - weights[n][3];
275 weights[n][2] = w + weights[n][0] - 2.0 * weights[n][3];
276 weights[n][1] = 1.0 - weights[n][0] - weights[n][2] - weights[n][3];
280 for (unsigned int n = 0; n < ImageDimension; n++)
282 weights[n][0] = 1; // implements nearest neighbor
286 for (unsigned int n = 0; n < ImageDimension; n++)
288 w = x[n] - (double) EvaluateIndex[n][0];
290 weights[n][0] = 1.0 - w;
294 for (unsigned int n = 0; n < ImageDimension; n++)
297 w = x[n] - (double)EvaluateIndex[n][1];
298 weights[n][1] = 0.75 - w * w;
299 weights[n][2] = 0.5 * (w - weights[n][1] + 1.0);
300 weights[n][0] = 1.0 - weights[n][1] - weights[n][2];
304 for (unsigned int n = 0; n < ImageDimension; n++)
307 w = x[n] - (double)EvaluateIndex[n][2];
309 t = (1.0 / 6.0) * w2;
310 weights[n][0] = 0.5 - w;
311 weights[n][0] *= weights[n][0];
312 weights[n][0] *= (1.0 / 24.0) * weights[n][0];
313 t0 = w * (t - 11.0 / 24.0);
314 t1 = 19.0 / 96.0 + w2 * (0.25 - t);
315 weights[n][1] = t1 + t0;
316 weights[n][3] = t1 - t0;
317 weights[n][4] = weights[n][0] + t0 + 0.5 * w;
318 weights[n][2] = 1.0 - weights[n][0] - weights[n][1] - weights[n][3] - weights[n][4];
322 for (unsigned int n = 0; n < ImageDimension; n++)
325 w = x[n] - (double)EvaluateIndex[n][2];
327 weights[n][5] = (1.0 / 120.0) * w * w2 * w2;
332 weights[n][0] = (1.0 / 24.0) * (1.0 / 5.0 + w2 + w4) - weights[n][5];
333 t0 = (1.0 / 24.0) * (w2 * (w2 - 5.0) + 46.0 / 5.0);
334 t1 = (-1.0 / 12.0) * w * (t + 4.0);
335 weights[n][2] = t0 + t1;
336 weights[n][3] = t0 - t1;
337 t0 = (1.0 / 16.0) * (9.0 / 5.0 - t);
338 t1 = (1.0 / 24.0) * w * (w4 - w2 - 5.0);
339 weights[n][1] = t0 + t1;
340 weights[n][4] = t0 - t1;
344 // SplineOrder not implemented yet.
345 itk::ExceptionObject err(__FILE__, __LINE__);
346 err.SetLocation( ITK_LOCATION );
347 err.SetDescription( "SplineOrder must be between 0 and 5. Requested spline order has not been implemented yet." );
354 template <class TImageType, class TCoordRep, class TCoefficientType>
356 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
357 ::SetDerivativeWeights( const ContinuousIndexType & x, const vnl_matrix<long> & EvaluateIndex,
358 vnl_matrix<double> & weights, unsigned int splineOrder ) const
360 // For speed improvements we could make each case a separate function and use
361 // function pointers to reference the correct weight order.
362 // Another possiblity would be to loop inside the case statement (reducing the number
363 // of switch statement executions to one per routine call.
364 // Left as is for now for readability.
365 double w, w1, w2, w3, w4, w5, t, t0, t1, t2;
366 int derivativeSplineOrder = (int) splineOrder -1;
368 switch (derivativeSplineOrder)
371 // Calculates B(splineOrder) ( (x + 1/2) - xi) - B(splineOrder -1) ( (x - 1/2) - xi)
373 // Why would we want to do this?
374 for (unsigned int n = 0; n < ImageDimension; n++)
380 for (unsigned int n = 0; n < ImageDimension; n++)
382 weights[n][0] = -1.0;
387 for (unsigned int n = 0; n < ImageDimension; n++)
389 w = x[n] + 0.5 - (double)EvaluateIndex[n][1];
393 weights[n][0] = 0.0 - w1;
394 weights[n][1] = w1 - w;
400 for (unsigned int n = 0; n < ImageDimension; n++)
402 w = x[n] + .5 - (double)EvaluateIndex[n][2];
404 w3 = 0.5 * (w - w2 + 1.0);
407 weights[n][0] = 0.0 - w1;
408 weights[n][1] = w1 - w2;
409 weights[n][2] = w2 - w3;
415 for (unsigned int n = 0; n < ImageDimension; n++)
417 w = x[n] + 0.5 - (double)EvaluateIndex[n][2];
418 w4 = (1.0 / 6.0) * w * w * w;
419 w1 = (1.0 / 6.0) + 0.5 * w * (w - 1.0) - w4;
420 w3 = w + w1 - 2.0 * w4;
421 w2 = 1.0 - w1 - w3 - w4;
423 weights[n][0] = 0.0 - w1;
424 weights[n][1] = w1 - w2;
425 weights[n][2] = w2 - w3;
426 weights[n][3] = w3 - w4;
431 for (unsigned int n = 0; n < ImageDimension; n++)
433 w = x[n] + .5 - (double)EvaluateIndex[n][3];
435 t = (1.0 / 6.0) * t2;
438 w1 *= (1.0 / 24.0) * w1;
439 t0 = w * (t - 11.0 / 24.0);
440 t1 = 19.0 / 96.0 + t2 * (0.25 - t);
443 w5 = w1 + t0 + 0.5 * w;
444 w3 = 1.0 - w1 - w2 - w4 - w5;
446 weights[n][0] = 0.0 - w1;
447 weights[n][1] = w1 - w2;
448 weights[n][2] = w2 - w3;
449 weights[n][3] = w3 - w4;
450 weights[n][4] = w4 - w5;
456 // SplineOrder not implemented yet.
457 itk::ExceptionObject err(__FILE__, __LINE__);
458 err.SetLocation( ITK_LOCATION );
459 err.SetDescription( "SplineOrder (for derivatives) must be between 1 and 5. Requested spline order has not been implemented yet." );
467 // Generates m_PointsToIndex;
468 template <class TImageType, class TCoordRep, class TCoefficientType>
470 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
471 ::GeneratePointsToIndex( )
473 // m_PointsToIndex is used to convert a sequential location to an N-dimension
474 // index vector. This is precomputed to save time during the interpolation routine.
475 m_PointsToIndex.resize(m_MaxNumberInterpolationPoints);
476 for (unsigned int p = 0; p < m_MaxNumberInterpolationPoints; p++)
479 unsigned long indexFactor[ImageDimension];
481 for (int j=1; j< static_cast<int>(ImageDimension); j++)
483 indexFactor[j] = indexFactor[j-1] * ( m_SplineOrder + 1 );
485 for (int j = (static_cast<int>(ImageDimension) - 1); j >= 0; j--)
487 m_PointsToIndex[p][j] = pp / indexFactor[j];
488 pp = pp % indexFactor[j];
493 template <class TImageType, class TCoordRep, class TCoefficientType>
495 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
496 ::DetermineRegionOfSupport( vnl_matrix<long> & evaluateIndex,
497 const ContinuousIndexType & x,
498 unsigned int splineOrder ) const
502 // compute the interpolation indexes
503 for (unsigned int n = 0; n< ImageDimension; n++)
505 if (splineOrder & 1) // Use this index calculation for odd splineOrder
507 indx = (long)vcl_floor((float)x[n]) - splineOrder / 2;
508 for (unsigned int k = 0; k <= splineOrder; k++)
510 evaluateIndex[n][k] = indx++;
513 else // Use this index calculation for even splineOrder
515 indx = (long)vcl_floor((float)(x[n] + 0.5)) - splineOrder / 2;
516 for (unsigned int k = 0; k <= splineOrder; k++)
518 evaluateIndex[n][k] = indx++;
524 template <class TImageType, class TCoordRep, class TCoefficientType>
526 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
527 ::ApplyMirrorBoundaryConditions(vnl_matrix<long> & evaluateIndex,
528 unsigned int splineOrder) const
530 for (unsigned int n = 0; n < ImageDimension; n++)
532 long dataLength2 = 2 * m_DataLength[n] - 2;
534 // apply the mirror boundary conditions
535 // TODO: We could implement other boundary options beside mirror
536 if (m_DataLength[n] == 1)
538 for (unsigned int k = 0; k <= splineOrder; k++)
540 evaluateIndex[n][k] = 0;
545 for (unsigned int k = 0; k <= splineOrder; k++)
547 // btw - Think about this couldn't this be replaced with a more elagent modulus method?
548 evaluateIndex[n][k] = (evaluateIndex[n][k] < 0L) ? (-evaluateIndex[n][k] - dataLength2 * ((-evaluateIndex[n][k]) / dataLength2))
549 : (evaluateIndex[n][k] - dataLength2 * (evaluateIndex[n][k] / dataLength2));
550 if ((long) m_DataLength[n] <= evaluateIndex[n][k])
552 evaluateIndex[n][k] = dataLength2 - evaluateIndex[n][k];