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://oncora1.lyon.fnclcc.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 _clitkVectorBSplineInterpolateImageFunction_txx
19 #define _clitkVectorBSplineInterpolateImageFunction_txx
20 #include "itkConfigure.h"
22 // Second, redirect to the optimized version if necessary
23 // #ifdef ITK_USE_OPTIMIZED_REGISTRATION_METHODS
24 // #include "itkOptVectorBSplineInterpolateImageFunction.txx"
27 #include "clitkVectorBSplineInterpolateImageFunction.h"
30 #include "itkImageLinearIteratorWithIndex.h"
31 #include "itkImageRegionConstIteratorWithIndex.h"
32 #include "itkImageRegionIterator.h"
34 #include "itkVector.h"
36 #include "itkMatrix.h"
44 template <class TImageType, class TCoordRep, class TCoefficientType>
45 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
46 ::VectorBSplineInterpolateImageFunction()
49 unsigned int SplineOrder = 3;
50 m_CoefficientFilter = CoefficientFilter::New();
51 // ***TODO: Should we store coefficients in a variable or retrieve from filter?
52 m_Coefficients = CoefficientImageType::New();
53 this->SetSplineOrder(SplineOrder);
54 this->m_UseImageDirection = false;
58 * Standard "PrintSelf" method
60 template <class TImageType, class TCoordRep, class TCoefficientType>
62 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
65 itk::Indent indent) const
67 Superclass::PrintSelf( os, indent );
68 os << indent << "Spline Order: " << m_SplineOrder << std::endl;
69 os << indent << "UseImageDirection = "
70 << (this->m_UseImageDirection ? "On" : "Off") << std::endl;
75 template <class TImageType, class TCoordRep, class TCoefficientType>
77 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
78 ::SetInputImage(const TImageType * inputData)
83 DD("calling decomposition filter");
84 m_CoefficientFilter->SetInput(inputData);
86 // the Coefficient Filter requires that the spline order and the input data be set.
87 // TODO: We need to ensure that this is only run once and only after both input and
88 // spline order have been set. Should we force an update after the
89 // splineOrder has been set also?
91 m_CoefficientFilter->Update();
92 m_Coefficients = m_CoefficientFilter->GetOutput();
95 // Call the Superclass implementation after, in case the filter
96 // pulls in more of the input image
97 Superclass::SetInputImage(inputData);
99 m_DataLength = inputData->GetBufferedRegion().GetSize();
104 m_Coefficients = NULL;
109 template <class TImageType, class TCoordRep, class TCoefficientType>
111 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
112 ::SetSplineOrder(unsigned int SplineOrder)
114 if (SplineOrder == m_SplineOrder)
118 m_SplineOrder = SplineOrder;
119 m_CoefficientFilter->SetSplineOrder( SplineOrder );
122 m_MaxNumberInterpolationPoints = 1;
123 for (unsigned int n=0; n < ImageDimension; n++)
125 m_MaxNumberInterpolationPoints *= ( m_SplineOrder + 1);
127 this->GeneratePointsToIndex( );
131 template <class TImageType, class TCoordRep, class TCoefficientType>
133 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
135 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
136 ::EvaluateAtContinuousIndex( const ContinuousIndexType & x ) const
138 vnl_matrix<long> EvaluateIndex(ImageDimension, ( m_SplineOrder + 1 ));
140 // compute the interpolation indexes
141 this->DetermineRegionOfSupport(EvaluateIndex, x, m_SplineOrder);
144 vnl_matrix<double> weights(ImageDimension, ( m_SplineOrder + 1 ));
145 SetInterpolationWeights( x, EvaluateIndex, weights, m_SplineOrder );
147 // Modify EvaluateIndex at the boundaries using mirror boundary conditions
148 this->ApplyMirrorBoundaryConditions(EvaluateIndex, m_SplineOrder);
150 // perform interpolation
152 itk::Vector<double, VectorDimension> interpolated;
153 for (unsigned int i=0; i< VectorDimension; i++) interpolated[i]=itk::NumericTraits<double>::Zero;
155 IndexType coefficientIndex;
156 // Step through eachpoint in the N-dimensional interpolation cube.
157 for (unsigned int p = 0; p < m_MaxNumberInterpolationPoints; p++)
159 // translate each step into the N-dimensional index.
160 // IndexType pointIndex = PointToIndex( p );
163 for (unsigned int n = 0; n < ImageDimension; n++ )
166 w *= weights[n][ m_PointsToIndex[p][n] ];
167 coefficientIndex[n] = EvaluateIndex[n][m_PointsToIndex[p][n]]; // Build up ND index for coefficients.
169 // Convert our step p to the appropriate point in ND space in the
170 // m_Coefficients cube.
171 //JV shouldn't be necessary
172 for (unsigned int i=0; i<VectorDimension; i++)
173 interpolated[i] += w * m_Coefficients->GetPixel(coefficientIndex)[i];
176 /* double interpolated = 0.0;
177 IndexType coefficientIndex;
178 // Step through eachpoint in the N-dimensional interpolation cube.
179 for (unsigned int sp = 0; sp <= m_SplineOrder; sp++)
181 for (unsigned int sp1=0; sp1 <= m_SplineOrder; sp1++)
185 for (unsigned int n1 = 0; n1 < ImageDimension; n1++ )
187 w *= weights[n1][ sp1 ];
188 coefficientIndex[n1] = EvaluateIndex[n1][sp]; // Build up ND index for coefficients.
191 interpolated += w * m_Coefficients->GetPixel(coefficientIndex);
195 return(interpolated);
200 template <class TImageType, class TCoordRep, class TCoefficientType>
202 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
203 :: CovariantVectorType
204 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
205 ::EvaluateDerivativeAtContinuousIndex( const ContinuousIndexType & x ) const
207 vnl_matrix<long> EvaluateIndex(ImageDimension, ( m_SplineOrder + 1 ));
209 // compute the interpolation indexes
210 // TODO: Do we need to revisit region of support for the derivatives?
211 this->DetermineRegionOfSupport(EvaluateIndex, x, m_SplineOrder);
214 vnl_matrix<double> weights(ImageDimension, ( m_SplineOrder + 1 ));
215 SetInterpolationWeights( x, EvaluateIndex, weights, m_SplineOrder );
217 vnl_matrix<double> weightsDerivative(ImageDimension, ( m_SplineOrder + 1));
218 SetDerivativeWeights( x, EvaluateIndex, weightsDerivative, ( m_SplineOrder ) );
220 // Modify EvaluateIndex at the boundaries using mirror boundary conditions
221 this->ApplyMirrorBoundaryConditions(EvaluateIndex, m_SplineOrder);
223 const InputImageType * inputImage = this->GetInputImage();
224 const typename InputImageType::SpacingType & spacing = inputImage->GetSpacing();
226 // Calculate derivative
227 CovariantVectorType derivativeValue;
229 IndexType coefficientIndex;
230 for (unsigned int n = 0; n < ImageDimension; n++)
232 derivativeValue[n] = 0.0;
233 for (unsigned int p = 0; p < m_MaxNumberInterpolationPoints; p++)
236 for (unsigned int n1 = 0; n1 < ImageDimension; n1++)
238 //coefficientIndex[n1] = EvaluateIndex[n1][sp];
239 coefficientIndex[n1] = EvaluateIndex[n1][m_PointsToIndex[p][n1]];
243 //w *= weights[n][ m_PointsToIndex[p][n] ];
244 tempValue *= weightsDerivative[n1][ m_PointsToIndex[p][n1] ];
248 tempValue *= weights[n1][ m_PointsToIndex[p][n1] ];
251 derivativeValue[n] += m_Coefficients->GetPixel(coefficientIndex) * tempValue ;
253 derivativeValue[n] /= spacing[n]; // take spacing into account
256 #ifdef ITK_USE_ORIENTED_IMAGE_DIRECTION
257 if( this->m_UseImageDirection )
259 CovariantVectorType orientedDerivative;
260 inputImage->TransformLocalVectorToPhysicalVector( derivativeValue, orientedDerivative );
261 return orientedDerivative;
265 return(derivativeValue);
269 template <class TImageType, class TCoordRep, class TCoefficientType>
271 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
272 ::SetInterpolationWeights( const ContinuousIndexType & x, const vnl_matrix<long> & EvaluateIndex,
273 vnl_matrix<double> & weights, unsigned int splineOrder ) const
275 // For speed improvements we could make each case a separate function and use
276 // function pointers to reference the correct weight order.
277 // Left as is for now for readability.
278 double w, w2, w4, t, t0, t1;
283 for (unsigned int n = 0; n < ImageDimension; n++)
285 w = x[n] - (double) EvaluateIndex[n][1];
286 weights[n][3] = (1.0 / 6.0) * w * w * w;
287 weights[n][0] = (1.0 / 6.0) + 0.5 * w * (w - 1.0) - weights[n][3];
288 weights[n][2] = w + weights[n][0] - 2.0 * weights[n][3];
289 weights[n][1] = 1.0 - weights[n][0] - weights[n][2] - weights[n][3];
293 for (unsigned int n = 0; n < ImageDimension; n++)
295 weights[n][0] = 1; // implements nearest neighbor
299 for (unsigned int n = 0; n < ImageDimension; n++)
301 w = x[n] - (double) EvaluateIndex[n][0];
303 weights[n][0] = 1.0 - w;
307 for (unsigned int n = 0; n < ImageDimension; n++)
310 w = x[n] - (double)EvaluateIndex[n][1];
311 weights[n][1] = 0.75 - w * w;
312 weights[n][2] = 0.5 * (w - weights[n][1] + 1.0);
313 weights[n][0] = 1.0 - weights[n][1] - weights[n][2];
317 for (unsigned int n = 0; n < ImageDimension; n++)
320 w = x[n] - (double)EvaluateIndex[n][2];
322 t = (1.0 / 6.0) * w2;
323 weights[n][0] = 0.5 - w;
324 weights[n][0] *= weights[n][0];
325 weights[n][0] *= (1.0 / 24.0) * weights[n][0];
326 t0 = w * (t - 11.0 / 24.0);
327 t1 = 19.0 / 96.0 + w2 * (0.25 - t);
328 weights[n][1] = t1 + t0;
329 weights[n][3] = t1 - t0;
330 weights[n][4] = weights[n][0] + t0 + 0.5 * w;
331 weights[n][2] = 1.0 - weights[n][0] - weights[n][1] - weights[n][3] - weights[n][4];
335 for (unsigned int n = 0; n < ImageDimension; n++)
338 w = x[n] - (double)EvaluateIndex[n][2];
340 weights[n][5] = (1.0 / 120.0) * w * w2 * w2;
345 weights[n][0] = (1.0 / 24.0) * (1.0 / 5.0 + w2 + w4) - weights[n][5];
346 t0 = (1.0 / 24.0) * (w2 * (w2 - 5.0) + 46.0 / 5.0);
347 t1 = (-1.0 / 12.0) * w * (t + 4.0);
348 weights[n][2] = t0 + t1;
349 weights[n][3] = t0 - t1;
350 t0 = (1.0 / 16.0) * (9.0 / 5.0 - t);
351 t1 = (1.0 / 24.0) * w * (w4 - w2 - 5.0);
352 weights[n][1] = t0 + t1;
353 weights[n][4] = t0 - t1;
357 // SplineOrder not implemented yet.
358 itk::ExceptionObject err(__FILE__, __LINE__);
359 err.SetLocation( ITK_LOCATION );
360 err.SetDescription( "SplineOrder must be between 0 and 5. Requested spline order has not been implemented yet." );
367 template <class TImageType, class TCoordRep, class TCoefficientType>
369 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
370 ::SetDerivativeWeights( const ContinuousIndexType & x, const vnl_matrix<long> & EvaluateIndex,
371 vnl_matrix<double> & weights, unsigned int splineOrder ) const
373 // For speed improvements we could make each case a separate function and use
374 // function pointers to reference the correct weight order.
375 // Another possiblity would be to loop inside the case statement (reducing the number
376 // of switch statement executions to one per routine call.
377 // Left as is for now for readability.
378 double w, w1, w2, w3, w4, w5, t, t0, t1, t2;
379 int derivativeSplineOrder = (int) splineOrder -1;
381 switch (derivativeSplineOrder)
384 // Calculates B(splineOrder) ( (x + 1/2) - xi) - B(splineOrder -1) ( (x - 1/2) - xi)
386 // Why would we want to do this?
387 for (unsigned int n = 0; n < ImageDimension; n++)
393 for (unsigned int n = 0; n < ImageDimension; n++)
395 weights[n][0] = -1.0;
400 for (unsigned int n = 0; n < ImageDimension; n++)
402 w = x[n] + 0.5 - (double)EvaluateIndex[n][1];
406 weights[n][0] = 0.0 - w1;
407 weights[n][1] = w1 - w;
413 for (unsigned int n = 0; n < ImageDimension; n++)
415 w = x[n] + .5 - (double)EvaluateIndex[n][2];
417 w3 = 0.5 * (w - w2 + 1.0);
420 weights[n][0] = 0.0 - w1;
421 weights[n][1] = w1 - w2;
422 weights[n][2] = w2 - w3;
428 for (unsigned int n = 0; n < ImageDimension; n++)
430 w = x[n] + 0.5 - (double)EvaluateIndex[n][2];
431 w4 = (1.0 / 6.0) * w * w * w;
432 w1 = (1.0 / 6.0) + 0.5 * w * (w - 1.0) - w4;
433 w3 = w + w1 - 2.0 * w4;
434 w2 = 1.0 - w1 - w3 - w4;
436 weights[n][0] = 0.0 - w1;
437 weights[n][1] = w1 - w2;
438 weights[n][2] = w2 - w3;
439 weights[n][3] = w3 - w4;
444 for (unsigned int n = 0; n < ImageDimension; n++)
446 w = x[n] + .5 - (double)EvaluateIndex[n][3];
448 t = (1.0 / 6.0) * t2;
451 w1 *= (1.0 / 24.0) * w1;
452 t0 = w * (t - 11.0 / 24.0);
453 t1 = 19.0 / 96.0 + t2 * (0.25 - t);
456 w5 = w1 + t0 + 0.5 * w;
457 w3 = 1.0 - w1 - w2 - w4 - w5;
459 weights[n][0] = 0.0 - w1;
460 weights[n][1] = w1 - w2;
461 weights[n][2] = w2 - w3;
462 weights[n][3] = w3 - w4;
463 weights[n][4] = w4 - w5;
469 // SplineOrder not implemented yet.
470 itk::ExceptionObject err(__FILE__, __LINE__);
471 err.SetLocation( ITK_LOCATION );
472 err.SetDescription( "SplineOrder (for derivatives) must be between 1 and 5. Requested spline order has not been implemented yet." );
480 // Generates m_PointsToIndex;
481 template <class TImageType, class TCoordRep, class TCoefficientType>
483 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
484 ::GeneratePointsToIndex( )
486 // m_PointsToIndex is used to convert a sequential location to an N-dimension
487 // index vector. This is precomputed to save time during the interpolation routine.
488 m_PointsToIndex.resize(m_MaxNumberInterpolationPoints);
489 for (unsigned int p = 0; p < m_MaxNumberInterpolationPoints; p++)
492 unsigned long indexFactor[ImageDimension];
494 for (int j=1; j< static_cast<int>(ImageDimension); j++)
496 indexFactor[j] = indexFactor[j-1] * ( m_SplineOrder + 1 );
498 for (int j = (static_cast<int>(ImageDimension) - 1); j >= 0; j--)
500 m_PointsToIndex[p][j] = pp / indexFactor[j];
501 pp = pp % indexFactor[j];
506 template <class TImageType, class TCoordRep, class TCoefficientType>
508 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
509 ::DetermineRegionOfSupport( vnl_matrix<long> & evaluateIndex,
510 const ContinuousIndexType & x,
511 unsigned int splineOrder ) const
515 // compute the interpolation indexes
516 for (unsigned int n = 0; n< ImageDimension; n++)
518 if (splineOrder & 1) // Use this index calculation for odd splineOrder
520 indx = (long)vcl_floor((float)x[n]) - splineOrder / 2;
521 for (unsigned int k = 0; k <= splineOrder; k++)
523 evaluateIndex[n][k] = indx++;
526 else // Use this index calculation for even splineOrder
528 indx = (long)vcl_floor((float)(x[n] + 0.5)) - splineOrder / 2;
529 for (unsigned int k = 0; k <= splineOrder; k++)
531 evaluateIndex[n][k] = indx++;
537 template <class TImageType, class TCoordRep, class TCoefficientType>
539 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
540 ::ApplyMirrorBoundaryConditions(vnl_matrix<long> & evaluateIndex,
541 unsigned int splineOrder) const
543 for (unsigned int n = 0; n < ImageDimension; n++)
545 long dataLength2 = 2 * m_DataLength[n] - 2;
547 // apply the mirror boundary conditions
548 // TODO: We could implement other boundary options beside mirror
549 if (m_DataLength[n] == 1)
551 for (unsigned int k = 0; k <= splineOrder; k++)
553 evaluateIndex[n][k] = 0;
558 for (unsigned int k = 0; k <= splineOrder; k++)
560 // btw - Think about this couldn't this be replaced with a more elagent modulus method?
561 evaluateIndex[n][k] = (evaluateIndex[n][k] < 0L) ? (-evaluateIndex[n][k] - dataLength2 * ((-evaluateIndex[n][k]) / dataLength2))
562 : (evaluateIndex[n][k] - dataLength2 * (evaluateIndex[n][k] / dataLength2));
563 if ((long) m_DataLength[n] <= evaluateIndex[n][k])
565 evaluateIndex[n][k] = dataLength2 - evaluateIndex[n][k];