]> Creatis software - clitk.git/blob - itk/clitkVectorBSplineInterpolateImageFunction.txx
Merge branch 'master' of git.creatis.insa-lyon.fr:clitk
[clitk.git] / itk / clitkVectorBSplineInterpolateImageFunction.txx
1 /*=========================================================================
2   Program:   vv                     http://www.creatis.insa-lyon.fr/rio/vv
3
4   Authors belong to:
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
8
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.
12
13   It is distributed under dual licence
14
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"
21
22 // Second, redirect to the optimized version if necessary
23 // #ifdef ITK_USE_OPTIMIZED_REGISTRATION_METHODS
24 // #include "itkOptVectorBSplineInterpolateImageFunction.txx"
25 // #else
26
27 #include "clitkVectorBSplineInterpolateImageFunction.h"
28
29
30 #include "itkImageLinearIteratorWithIndex.h"
31 #include "itkImageRegionConstIteratorWithIndex.h"
32 #include "itkImageRegionIterator.h"
33
34 #include "itkVector.h"
35
36 #include "itkMatrix.h"
37
38 namespace clitk
39 {
40
41 /**
42  * Constructor
43  */
44 template <class TImageType, class TCoordRep, class TCoefficientType>
45 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
46 ::VectorBSplineInterpolateImageFunction()
47 {
48   m_SplineOrder = 0;
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;
55 }
56
57 /**
58  * Standard "PrintSelf" method
59  */
60 template <class TImageType, class TCoordRep, class TCoefficientType>
61 void
62 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
63 ::PrintSelf(
64   std::ostream& os,
65   itk::Indent indent) const
66 {
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;
71
72 }
73
74
75 template <class TImageType, class TCoordRep, class TCoefficientType>
76 void
77 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
78 ::SetInputImage(const TImageType * inputData)
79 {
80   if ( inputData ) {
81
82     DD("calling decomposition filter");
83     m_CoefficientFilter->SetInput(inputData);
84
85     // the Coefficient Filter requires that the spline order and the input data be set.
86     // TODO:  We need to ensure that this is only run once and only after both input and
87     //        spline order have been set. Should we force an update after the
88     //        splineOrder has been set also?
89
90     m_CoefficientFilter->Update();
91     m_Coefficients = m_CoefficientFilter->GetOutput();
92
93
94     // Call the Superclass implementation after, in case the filter
95     // pulls in  more of the input image
96     Superclass::SetInputImage(inputData);
97
98     m_DataLength = inputData->GetBufferedRegion().GetSize();
99
100   } else {
101     m_Coefficients = NULL;
102   }
103 }
104
105
106 template <class TImageType, class TCoordRep, class TCoefficientType>
107 void
108 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
109 ::SetSplineOrder(unsigned int SplineOrder)
110 {
111   if (SplineOrder == m_SplineOrder) {
112     return;
113   }
114   m_SplineOrder = SplineOrder;
115   m_CoefficientFilter->SetSplineOrder( SplineOrder );
116
117   //this->SetPoles();/
118   m_MaxNumberInterpolationPoints = 1;
119   for (unsigned int n=0; n < ImageDimension; n++) {
120     m_MaxNumberInterpolationPoints *= ( m_SplineOrder + 1);
121   }
122   this->GeneratePointsToIndex( );
123 }
124
125
126 template <class TImageType, class TCoordRep, class TCoefficientType>
127 typename
128 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
129 ::OutputType
130 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
131 ::EvaluateAtContinuousIndex( const ContinuousIndexType & x ) const
132 {
133   vnl_matrix<long>        EvaluateIndex(ImageDimension, ( m_SplineOrder + 1 ));
134
135   // compute the interpolation indexes
136   this->DetermineRegionOfSupport(EvaluateIndex, x, m_SplineOrder);
137
138   // Determine weights
139   vnl_matrix<double>        weights(ImageDimension, ( m_SplineOrder + 1 ));
140   SetInterpolationWeights( x, EvaluateIndex, weights, m_SplineOrder );
141
142   // Modify EvaluateIndex at the boundaries using mirror boundary conditions
143   this->ApplyMirrorBoundaryConditions(EvaluateIndex, m_SplineOrder);
144
145   // perform interpolation
146   //JV
147   itk::Vector<double, VectorDimension> interpolated;
148   for (unsigned int i=0; i< VectorDimension; i++) interpolated[i]=itk::NumericTraits<double>::Zero;
149
150   IndexType coefficientIndex;
151   // Step through eachpoint in the N-dimensional interpolation cube.
152   for (unsigned int p = 0; p < m_MaxNumberInterpolationPoints; p++) {
153     // translate each step into the N-dimensional index.
154     //      IndexType pointIndex = PointToIndex( p );
155
156     double w = 1.0;
157     for (unsigned int n = 0; n < ImageDimension; n++ ) {
158
159       w *= weights[n][ m_PointsToIndex[p][n] ];
160       coefficientIndex[n] = EvaluateIndex[n][m_PointsToIndex[p][n]];  // Build up ND index for coefficients.
161     }
162     // Convert our step p to the appropriate point in ND space in the
163     // m_Coefficients cube.
164     //JV shouldn't be necessary
165     for (unsigned int i=0; i<VectorDimension; i++)
166       interpolated[i] += w * m_Coefficients->GetPixel(coefficientIndex)[i];
167   }
168
169   /*  double interpolated = 0.0;
170     IndexType coefficientIndex;
171     // Step through eachpoint in the N-dimensional interpolation cube.
172     for (unsigned int sp = 0; sp <= m_SplineOrder; sp++)
173       {
174       for (unsigned int sp1=0; sp1 <= m_SplineOrder; sp1++)
175         {
176
177         double w = 1.0;
178         for (unsigned int n1 = 0; n1 < ImageDimension; n1++ )
179           {
180           w *= weights[n1][ sp1 ];
181           coefficientIndex[n1] = EvaluateIndex[n1][sp];  // Build up ND index for coefficients.
182           }
183
184           interpolated += w * m_Coefficients->GetPixel(coefficientIndex);
185         }
186       }
187   */
188   return(interpolated);
189
190 }
191
192
193 template <class TImageType, class TCoordRep, class TCoefficientType>
194 typename
195 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
196 :: CovariantVectorType
197 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
198 ::EvaluateDerivativeAtContinuousIndex( const ContinuousIndexType & x ) const
199 {
200   vnl_matrix<long>        EvaluateIndex(ImageDimension, ( m_SplineOrder + 1 ));
201
202   // compute the interpolation indexes
203   // TODO: Do we need to revisit region of support for the derivatives?
204   this->DetermineRegionOfSupport(EvaluateIndex, x, m_SplineOrder);
205
206   // Determine weights
207   vnl_matrix<double>        weights(ImageDimension, ( m_SplineOrder + 1 ));
208   SetInterpolationWeights( x, EvaluateIndex, weights, m_SplineOrder );
209
210   vnl_matrix<double>        weightsDerivative(ImageDimension, ( m_SplineOrder + 1));
211   SetDerivativeWeights( x, EvaluateIndex, weightsDerivative, ( m_SplineOrder ) );
212
213   // Modify EvaluateIndex at the boundaries using mirror boundary conditions
214   this->ApplyMirrorBoundaryConditions(EvaluateIndex, m_SplineOrder);
215
216   const InputImageType * inputImage = this->GetInputImage();
217   const typename InputImageType::SpacingType & spacing = inputImage->GetSpacing();
218
219   // Calculate derivative
220   CovariantVectorType derivativeValue;
221   double tempValue;
222   IndexType coefficientIndex;
223   for (unsigned int n = 0; n < ImageDimension; n++) {
224     derivativeValue[n] = 0.0;
225     for (unsigned int p = 0; p < m_MaxNumberInterpolationPoints; p++) {
226       tempValue = 1.0 ;
227       for (unsigned int n1 = 0; n1 < ImageDimension; n1++) {
228         //coefficientIndex[n1] = EvaluateIndex[n1][sp];
229         coefficientIndex[n1] = EvaluateIndex[n1][m_PointsToIndex[p][n1]];
230
231         if (n1 == n) {
232           //w *= weights[n][ m_PointsToIndex[p][n] ];
233           tempValue *= weightsDerivative[n1][ m_PointsToIndex[p][n1] ];
234         } else {
235           tempValue *= weights[n1][ m_PointsToIndex[p][n1] ];
236         }
237       }
238       derivativeValue[n] += m_Coefficients->GetPixel(coefficientIndex) * tempValue ;
239     }
240     derivativeValue[n] /= spacing[n];   // take spacing into account
241   }
242
243 #ifdef ITK_USE_ORIENTED_IMAGE_DIRECTION
244   if( this->m_UseImageDirection ) {
245     CovariantVectorType orientedDerivative;
246     inputImage->TransformLocalVectorToPhysicalVector( derivativeValue, orientedDerivative );
247     return orientedDerivative;
248   }
249 #endif
250
251   return(derivativeValue);
252 }
253
254
255 template <class TImageType, class TCoordRep, class TCoefficientType>
256 void
257 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
258 ::SetInterpolationWeights( const ContinuousIndexType & x, const vnl_matrix<long> & EvaluateIndex,
259                            vnl_matrix<double> & weights, unsigned int splineOrder ) const
260 {
261   // For speed improvements we could make each case a separate function and use
262   // function pointers to reference the correct weight order.
263   // Left as is for now for readability.
264   double w, w2, w4, t, t0, t1;
265
266   switch (splineOrder) {
267   case 3:
268     for (unsigned int n = 0; n < ImageDimension; n++) {
269       w = x[n] - (double) EvaluateIndex[n][1];
270       weights[n][3] = (1.0 / 6.0) * w * w * w;
271       weights[n][0] = (1.0 / 6.0) + 0.5 * w * (w - 1.0) - weights[n][3];
272       weights[n][2] = w + weights[n][0] - 2.0 * weights[n][3];
273       weights[n][1] = 1.0 - weights[n][0] - weights[n][2] - weights[n][3];
274     }
275     break;
276   case 0:
277     for (unsigned int n = 0; n < ImageDimension; n++) {
278       weights[n][0] = 1; // implements nearest neighbor
279     }
280     break;
281   case 1:
282     for (unsigned int n = 0; n < ImageDimension; n++) {
283       w = x[n] - (double) EvaluateIndex[n][0];
284       weights[n][1] = w;
285       weights[n][0] = 1.0 - w;
286     }
287     break;
288   case 2:
289     for (unsigned int n = 0; n < ImageDimension; n++) {
290       /* x */
291       w = x[n] - (double)EvaluateIndex[n][1];
292       weights[n][1] = 0.75 - w * w;
293       weights[n][2] = 0.5 * (w - weights[n][1] + 1.0);
294       weights[n][0] = 1.0 - weights[n][1] - weights[n][2];
295     }
296     break;
297   case 4:
298     for (unsigned int n = 0; n < ImageDimension; n++) {
299       /* x */
300       w = x[n] - (double)EvaluateIndex[n][2];
301       w2 = w * w;
302       t = (1.0 / 6.0) * w2;
303       weights[n][0] = 0.5 - w;
304       weights[n][0] *= weights[n][0];
305       weights[n][0] *= (1.0 / 24.0) * weights[n][0];
306       t0 = w * (t - 11.0 / 24.0);
307       t1 = 19.0 / 96.0 + w2 * (0.25 - t);
308       weights[n][1] = t1 + t0;
309       weights[n][3] = t1 - t0;
310       weights[n][4] = weights[n][0] + t0 + 0.5 * w;
311       weights[n][2] = 1.0 - weights[n][0] - weights[n][1] - weights[n][3] - weights[n][4];
312     }
313     break;
314   case 5:
315     for (unsigned int n = 0; n < ImageDimension; n++) {
316       /* x */
317       w = x[n] - (double)EvaluateIndex[n][2];
318       w2 = w * w;
319       weights[n][5] = (1.0 / 120.0) * w * w2 * w2;
320       w2 -= w;
321       w4 = w2 * w2;
322       w -= 0.5;
323       t = w2 * (w2 - 3.0);
324       weights[n][0] = (1.0 / 24.0) * (1.0 / 5.0 + w2 + w4) - weights[n][5];
325       t0 = (1.0 / 24.0) * (w2 * (w2 - 5.0) + 46.0 / 5.0);
326       t1 = (-1.0 / 12.0) * w * (t + 4.0);
327       weights[n][2] = t0 + t1;
328       weights[n][3] = t0 - t1;
329       t0 = (1.0 / 16.0) * (9.0 / 5.0 - t);
330       t1 = (1.0 / 24.0) * w * (w4 - w2 - 5.0);
331       weights[n][1] = t0 + t1;
332       weights[n][4] = t0 - t1;
333     }
334     break;
335   default:
336     // SplineOrder not implemented yet.
337     itk::ExceptionObject err(__FILE__, __LINE__);
338     err.SetLocation( ITK_LOCATION );
339     err.SetDescription( "SplineOrder must be between 0 and 5. Requested spline order has not been implemented yet." );
340     throw err;
341     break;
342   }
343
344 }
345
346 template <class TImageType, class TCoordRep, class TCoefficientType>
347 void
348 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
349 ::SetDerivativeWeights( const ContinuousIndexType & x, const vnl_matrix<long> & EvaluateIndex,
350                         vnl_matrix<double> & weights, unsigned int splineOrder ) const
351 {
352   // For speed improvements we could make each case a separate function and use
353   // function pointers to reference the correct weight order.
354   // Another possiblity would be to loop inside the case statement (reducing the number
355   // of switch statement executions to one per routine call.
356   // Left as is for now for readability.
357   double w, w1, w2, w3, w4, w5, t, t0, t1, t2;
358   int derivativeSplineOrder = (int) splineOrder -1;
359
360   switch (derivativeSplineOrder) {
361
362     // Calculates B(splineOrder) ( (x + 1/2) - xi) - B(splineOrder -1) ( (x - 1/2) - xi)
363   case -1:
364     // Why would we want to do this?
365     for (unsigned int n = 0; n < ImageDimension; n++) {
366       weights[n][0] = 0.0;
367     }
368     break;
369   case 0:
370     for (unsigned int n = 0; n < ImageDimension; n++) {
371       weights[n][0] = -1.0;
372       weights[n][1] =  1.0;
373     }
374     break;
375   case 1:
376     for (unsigned int n = 0; n < ImageDimension; n++) {
377       w = x[n] + 0.5 - (double)EvaluateIndex[n][1];
378       // w2 = w;
379       w1 = 1.0 - w;
380
381       weights[n][0] = 0.0 - w1;
382       weights[n][1] = w1 - w;
383       weights[n][2] = w;
384     }
385     break;
386   case 2:
387
388     for (unsigned int n = 0; n < ImageDimension; n++) {
389       w = x[n] + .5 - (double)EvaluateIndex[n][2];
390       w2 = 0.75 - w * w;
391       w3 = 0.5 * (w - w2 + 1.0);
392       w1 = 1.0 - w2 - w3;
393
394       weights[n][0] = 0.0 - w1;
395       weights[n][1] = w1 - w2;
396       weights[n][2] = w2 - w3;
397       weights[n][3] = w3;
398     }
399     break;
400   case 3:
401
402     for (unsigned int n = 0; n < ImageDimension; n++) {
403       w = x[n] + 0.5 - (double)EvaluateIndex[n][2];
404       w4 = (1.0 / 6.0) * w * w * w;
405       w1 = (1.0 / 6.0) + 0.5 * w * (w - 1.0) - w4;
406       w3 = w + w1 - 2.0 * w4;
407       w2 = 1.0 - w1 - w3 - w4;
408
409       weights[n][0] = 0.0 - w1;
410       weights[n][1] = w1 - w2;
411       weights[n][2] = w2 - w3;
412       weights[n][3] = w3 - w4;
413       weights[n][4] = w4;
414     }
415     break;
416   case 4:
417     for (unsigned int n = 0; n < ImageDimension; n++) {
418       w = x[n] + .5 - (double)EvaluateIndex[n][3];
419       t2 = w * w;
420       t = (1.0 / 6.0) * t2;
421       w1 = 0.5 - w;
422       w1 *= w1;
423       w1 *= (1.0 / 24.0) * w1;
424       t0 = w * (t - 11.0 / 24.0);
425       t1 = 19.0 / 96.0 + t2 * (0.25 - t);
426       w2 = t1 + t0;
427       w4 = t1 - t0;
428       w5 = w1 + t0 + 0.5 * w;
429       w3 = 1.0 - w1 - w2 - w4 - w5;
430
431       weights[n][0] = 0.0 - w1;
432       weights[n][1] = w1 - w2;
433       weights[n][2] = w2 - w3;
434       weights[n][3] = w3 - w4;
435       weights[n][4] = w4 - w5;
436       weights[n][5] = w5;
437     }
438     break;
439
440   default:
441     // SplineOrder not implemented yet.
442     itk::ExceptionObject err(__FILE__, __LINE__);
443     err.SetLocation( ITK_LOCATION );
444     err.SetDescription( "SplineOrder (for derivatives) must be between 1 and 5. Requested spline order has not been implemented yet." );
445     throw err;
446     break;
447   }
448
449 }
450
451
452 // Generates m_PointsToIndex;
453 template <class TImageType, class TCoordRep, class TCoefficientType>
454 void
455 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
456 ::GeneratePointsToIndex( )
457 {
458   // m_PointsToIndex is used to convert a sequential location to an N-dimension
459   // index vector.  This is precomputed to save time during the interpolation routine.
460   m_PointsToIndex.resize(m_MaxNumberInterpolationPoints);
461   for (unsigned int p = 0; p < m_MaxNumberInterpolationPoints; p++) {
462     int pp = p;
463     unsigned long indexFactor[ImageDimension];
464     indexFactor[0] = 1;
465     for (int j=1; j< static_cast<int>(ImageDimension); j++) {
466       indexFactor[j] = indexFactor[j-1] * ( m_SplineOrder + 1 );
467     }
468     for (int j = (static_cast<int>(ImageDimension) - 1); j >= 0; j--) {
469       m_PointsToIndex[p][j] = pp / indexFactor[j];
470       pp = pp % indexFactor[j];
471     }
472   }
473 }
474
475 template <class TImageType, class TCoordRep, class TCoefficientType>
476 void
477 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
478 ::DetermineRegionOfSupport( vnl_matrix<long> & evaluateIndex,
479                             const ContinuousIndexType & x,
480                             unsigned int splineOrder ) const
481 {
482   long indx;
483
484 // compute the interpolation indexes
485   for (unsigned int n = 0; n< ImageDimension; n++) {
486     if (splineOrder & 1) {   // Use this index calculation for odd splineOrder
487       indx = (long)vcl_floor((float)x[n]) - splineOrder / 2;
488       for (unsigned int k = 0; k <= splineOrder; k++) {
489         evaluateIndex[n][k] = indx++;
490       }
491     } else {                   // Use this index calculation for even splineOrder
492       indx = (long)vcl_floor((float)(x[n] + 0.5)) - splineOrder / 2;
493       for (unsigned int k = 0; k <= splineOrder; k++) {
494         evaluateIndex[n][k] = indx++;
495       }
496     }
497   }
498 }
499
500 template <class TImageType, class TCoordRep, class TCoefficientType>
501 void
502 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
503 ::ApplyMirrorBoundaryConditions(vnl_matrix<long> & evaluateIndex,
504                                 unsigned int splineOrder) const
505 {
506   for (unsigned int n = 0; n < ImageDimension; n++) {
507     long dataLength2 = 2 * m_DataLength[n] - 2;
508
509     // apply the mirror boundary conditions
510     // TODO:  We could implement other boundary options beside mirror
511     if (m_DataLength[n] == 1) {
512       for (unsigned int k = 0; k <= splineOrder; k++) {
513         evaluateIndex[n][k] = 0;
514       }
515     } else {
516       for (unsigned int k = 0; k <= splineOrder; k++) {
517         // btw - Think about this couldn't this be replaced with a more elagent modulus method?
518         evaluateIndex[n][k] = (evaluateIndex[n][k] < 0L) ? (-evaluateIndex[n][k] - dataLength2 * ((-evaluateIndex[n][k]) / dataLength2))
519                               : (evaluateIndex[n][k] - dataLength2 * (evaluateIndex[n][k] / dataLength2));
520         if ((long) m_DataLength[n] <= evaluateIndex[n][k]) {
521           evaluateIndex[n][k] = dataLength2 - evaluateIndex[n][k];
522         }
523       }
524     }
525   }
526 }
527
528 } // namespace itk
529
530 #endif
531
532 //#endif