]> Creatis software - clitk.git/blob - itk/clitkVectorBSplineInterpolateImageFunction.txx
Initial revision
[clitk.git] / itk / clitkVectorBSplineInterpolateImageFunction.txx
1 #ifndef _clitkVectorBSplineInterpolateImageFunction_txx
2 #define _clitkVectorBSplineInterpolateImageFunction_txx
3
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"
8
9 // Second, redirect to the optimized version if necessary
10 // #ifdef ITK_USE_OPTIMIZED_REGISTRATION_METHODS
11 // #include "itkOptVectorBSplineInterpolateImageFunction.txx"
12 // #else
13
14 #include "clitkVectorBSplineInterpolateImageFunction.h"
15
16
17 #include "itkImageLinearIteratorWithIndex.h"
18 #include "itkImageRegionConstIteratorWithIndex.h"
19 #include "itkImageRegionIterator.h"
20
21 #include "itkVector.h"
22
23 #include "itkMatrix.h"
24
25 namespace clitk
26 {
27
28 /**
29  * Constructor
30  */
31 template <class TImageType, class TCoordRep, class TCoefficientType>
32 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
33 ::VectorBSplineInterpolateImageFunction()
34 {
35   m_SplineOrder = 0;
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;
42 }
43
44 /**
45  * Standard "PrintSelf" method
46  */
47 template <class TImageType, class TCoordRep, class TCoefficientType>
48 void
49 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
50 ::PrintSelf(
51   std::ostream& os, 
52   itk::Indent indent) const
53 {
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;
58
59 }
60
61
62 template <class TImageType, class TCoordRep, class TCoefficientType>
63 void 
64 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
65 ::SetInputImage(const TImageType * inputData)
66 {
67   if ( inputData )
68     {
69         
70       DD("calling decomposition filter");
71       m_CoefficientFilter->SetInput(inputData);
72       
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?
77           
78           m_CoefficientFilter->Update();
79           m_Coefficients = m_CoefficientFilter->GetOutput();
80           
81        
82       // Call the Superclass implementation after, in case the filter
83       // pulls in  more of the input image
84       Superclass::SetInputImage(inputData);
85       
86       m_DataLength = inputData->GetBufferedRegion().GetSize();
87      
88     }
89   else
90    {
91    m_Coefficients = NULL;
92    }
93 }
94
95
96 template <class TImageType, class TCoordRep, class TCoefficientType>
97 void 
98 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
99 ::SetSplineOrder(unsigned int SplineOrder)
100 {
101   if (SplineOrder == m_SplineOrder)
102     {
103     return;
104     }
105   m_SplineOrder = SplineOrder;
106   m_CoefficientFilter->SetSplineOrder( SplineOrder );
107
108   //this->SetPoles();/
109   m_MaxNumberInterpolationPoints = 1;
110   for (unsigned int n=0; n < ImageDimension; n++)
111     {
112     m_MaxNumberInterpolationPoints *= ( m_SplineOrder + 1);
113     } 
114   this->GeneratePointsToIndex( );
115 }
116
117
118 template <class TImageType, class TCoordRep, class TCoefficientType>
119 typename 
120 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
121 ::OutputType
122 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
123 ::EvaluateAtContinuousIndex( const ContinuousIndexType & x ) const
124 {
125   vnl_matrix<long>        EvaluateIndex(ImageDimension, ( m_SplineOrder + 1 ));
126
127   // compute the interpolation indexes 
128   this->DetermineRegionOfSupport(EvaluateIndex, x, m_SplineOrder); 
129
130   // Determine weights
131   vnl_matrix<double>        weights(ImageDimension, ( m_SplineOrder + 1 ));
132   SetInterpolationWeights( x, EvaluateIndex, weights, m_SplineOrder );
133
134   // Modify EvaluateIndex at the boundaries using mirror boundary conditions
135   this->ApplyMirrorBoundaryConditions(EvaluateIndex, m_SplineOrder);
136   
137   // perform interpolation
138   //JV 
139   itk::Vector<double, VectorDimension> interpolated; 
140     for (unsigned int i=0; i< VectorDimension; i++) interpolated[i]=itk::NumericTraits<double>::Zero;
141
142   IndexType coefficientIndex;
143   // Step through eachpoint in the N-dimensional interpolation cube.
144   for (unsigned int p = 0; p < m_MaxNumberInterpolationPoints; p++)
145     {
146     // translate each step into the N-dimensional index.
147     //      IndexType pointIndex = PointToIndex( p );
148
149     double w = 1.0;
150     for (unsigned int n = 0; n < ImageDimension; n++ )
151       {
152
153         w *= weights[n][ m_PointsToIndex[p][n] ];
154         coefficientIndex[n] = EvaluateIndex[n][m_PointsToIndex[p][n]];  // Build up ND index for coefficients.
155       }
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];
161     }
162     
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++)
167     {
168     for (unsigned int sp1=0; sp1 <= m_SplineOrder; sp1++)
169       {
170
171       double w = 1.0;
172       for (unsigned int n1 = 0; n1 < ImageDimension; n1++ )
173         {
174         w *= weights[n1][ sp1 ];
175         coefficientIndex[n1] = EvaluateIndex[n1][sp];  // Build up ND index for coefficients.
176         }
177
178         interpolated += w * m_Coefficients->GetPixel(coefficientIndex);
179       }  
180     }
181 */
182   return(interpolated);
183     
184 }
185
186
187 template <class TImageType, class TCoordRep, class TCoefficientType>
188 typename 
189 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
190 :: CovariantVectorType
191 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
192 ::EvaluateDerivativeAtContinuousIndex( const ContinuousIndexType & x ) const
193 {
194   vnl_matrix<long>        EvaluateIndex(ImageDimension, ( m_SplineOrder + 1 ));
195
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);
199
200   // Determine weights
201   vnl_matrix<double>        weights(ImageDimension, ( m_SplineOrder + 1 ));
202   SetInterpolationWeights( x, EvaluateIndex, weights, m_SplineOrder );
203
204   vnl_matrix<double>        weightsDerivative(ImageDimension, ( m_SplineOrder + 1));
205   SetDerivativeWeights( x, EvaluateIndex, weightsDerivative, ( m_SplineOrder ) );
206
207   // Modify EvaluateIndex at the boundaries using mirror boundary conditions
208   this->ApplyMirrorBoundaryConditions(EvaluateIndex, m_SplineOrder);
209   
210   const InputImageType * inputImage = this->GetInputImage();
211   const typename InputImageType::SpacingType & spacing = inputImage->GetSpacing();
212
213   // Calculate derivative
214   CovariantVectorType derivativeValue;
215   double tempValue;
216   IndexType coefficientIndex;
217   for (unsigned int n = 0; n < ImageDimension; n++)
218     {
219     derivativeValue[n] = 0.0;
220     for (unsigned int p = 0; p < m_MaxNumberInterpolationPoints; p++)
221       {
222       tempValue = 1.0 ; 
223       for (unsigned int n1 = 0; n1 < ImageDimension; n1++)
224         {
225         //coefficientIndex[n1] = EvaluateIndex[n1][sp];
226         coefficientIndex[n1] = EvaluateIndex[n1][m_PointsToIndex[p][n1]];
227
228         if (n1 == n)
229           {
230           //w *= weights[n][ m_PointsToIndex[p][n] ];
231           tempValue *= weightsDerivative[n1][ m_PointsToIndex[p][n1] ];
232           }
233         else
234           {
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     {
246     CovariantVectorType orientedDerivative;
247     inputImage->TransformLocalVectorToPhysicalVector( derivativeValue, orientedDerivative );
248     return orientedDerivative;
249     }
250 #endif
251
252   return(derivativeValue);
253 }
254
255
256 template <class TImageType, class TCoordRep, class TCoefficientType>
257 void 
258 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
259 ::SetInterpolationWeights( const ContinuousIndexType & x, const vnl_matrix<long> & EvaluateIndex, 
260                            vnl_matrix<double> & weights, unsigned int splineOrder ) const
261 {
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;
266   
267   switch (splineOrder)
268     {
269     case 3:
270       for (unsigned int n = 0; n < ImageDimension; n++)
271         {
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];
277         }
278       break;
279     case 0:
280       for (unsigned int n = 0; n < ImageDimension; n++)
281         {
282         weights[n][0] = 1; // implements nearest neighbor
283         }
284       break;
285     case 1:
286       for (unsigned int n = 0; n < ImageDimension; n++)
287         {
288         w = x[n] - (double) EvaluateIndex[n][0];
289         weights[n][1] = w;
290         weights[n][0] = 1.0 - w;
291         }
292       break;
293     case 2:
294       for (unsigned int n = 0; n < ImageDimension; n++)
295         {
296         /* x */
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];
301         }
302       break;
303     case 4:
304       for (unsigned int n = 0; n < ImageDimension; n++)
305         {
306         /* x */
307         w = x[n] - (double)EvaluateIndex[n][2];
308         w2 = w * w;
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];
319         }
320       break;
321     case 5:
322       for (unsigned int n = 0; n < ImageDimension; n++)
323         {
324         /* x */
325         w = x[n] - (double)EvaluateIndex[n][2];
326         w2 = w * w;
327         weights[n][5] = (1.0 / 120.0) * w * w2 * w2;
328         w2 -= w;
329         w4 = w2 * w2;
330         w -= 0.5;
331         t = w2 * (w2 - 3.0);
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;
341         }
342       break;
343     default:
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." );
348       throw err;
349       break;
350     }
351     
352 }
353
354 template <class TImageType, class TCoordRep, class TCoefficientType>
355 void 
356 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
357 ::SetDerivativeWeights( const ContinuousIndexType & x, const vnl_matrix<long> & EvaluateIndex, 
358                         vnl_matrix<double> & weights, unsigned int splineOrder ) const
359 {
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;
367   
368   switch (derivativeSplineOrder)
369     {
370     
371     // Calculates B(splineOrder) ( (x + 1/2) - xi) - B(splineOrder -1) ( (x - 1/2) - xi)
372     case -1:
373       // Why would we want to do this?
374       for (unsigned int n = 0; n < ImageDimension; n++)
375         {
376         weights[n][0] = 0.0;
377         }
378       break;
379     case 0:
380       for (unsigned int n = 0; n < ImageDimension; n++)
381         {
382         weights[n][0] = -1.0;
383         weights[n][1] =  1.0;
384         }
385       break;
386     case 1:
387       for (unsigned int n = 0; n < ImageDimension; n++)
388         {
389         w = x[n] + 0.5 - (double)EvaluateIndex[n][1];
390         // w2 = w;
391         w1 = 1.0 - w;
392
393         weights[n][0] = 0.0 - w1;
394         weights[n][1] = w1 - w;
395         weights[n][2] = w; 
396         }
397       break;
398     case 2:
399       
400       for (unsigned int n = 0; n < ImageDimension; n++)
401         {
402         w = x[n] + .5 - (double)EvaluateIndex[n][2];
403         w2 = 0.75 - w * w;
404         w3 = 0.5 * (w - w2 + 1.0);
405         w1 = 1.0 - w2 - w3;
406
407         weights[n][0] = 0.0 - w1;
408         weights[n][1] = w1 - w2;
409         weights[n][2] = w2 - w3;
410         weights[n][3] = w3; 
411         }
412       break;
413     case 3:
414       
415       for (unsigned int n = 0; n < ImageDimension; n++)
416         {
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;
422
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;
427         weights[n][4] = w4;
428         }
429       break;
430     case 4:
431       for (unsigned int n = 0; n < ImageDimension; n++)
432         {
433         w = x[n] + .5 - (double)EvaluateIndex[n][3];
434         t2 = w * w;
435         t = (1.0 / 6.0) * t2;
436         w1 = 0.5 - w;
437         w1 *= w1;
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);
441         w2 = t1 + t0;
442         w4 = t1 - t0;
443         w5 = w1 + t0 + 0.5 * w;
444         w3 = 1.0 - w1 - w2 - w4 - w5;
445
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;
451         weights[n][5] = w5;
452         }
453       break;
454       
455     default:
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." );
460       throw err;
461       break;
462     }
463     
464 }
465
466
467 // Generates m_PointsToIndex;
468 template <class TImageType, class TCoordRep, class TCoefficientType>
469 void
470 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
471 ::GeneratePointsToIndex( )
472 {
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++)
477     {
478     int pp = p;
479     unsigned long indexFactor[ImageDimension];
480     indexFactor[0] = 1;
481     for (int j=1; j< static_cast<int>(ImageDimension); j++)
482       {
483       indexFactor[j] = indexFactor[j-1] * ( m_SplineOrder + 1 );
484       }
485     for (int j = (static_cast<int>(ImageDimension) - 1); j >= 0; j--)
486       {
487       m_PointsToIndex[p][j] = pp / indexFactor[j];
488       pp = pp % indexFactor[j];
489       }
490     }
491 }
492
493 template <class TImageType, class TCoordRep, class TCoefficientType>
494 void
495 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
496 ::DetermineRegionOfSupport( vnl_matrix<long> & evaluateIndex, 
497                             const ContinuousIndexType & x, 
498                             unsigned int splineOrder ) const
499
500   long indx;
501
502 // compute the interpolation indexes 
503   for (unsigned int n = 0; n< ImageDimension; n++)
504     {
505     if (splineOrder & 1)     // Use this index calculation for odd splineOrder
506       {
507       indx = (long)vcl_floor((float)x[n]) - splineOrder / 2;
508       for (unsigned int k = 0; k <= splineOrder; k++)
509         {
510         evaluateIndex[n][k] = indx++;
511         }
512       }
513     else                       // Use this index calculation for even splineOrder
514       { 
515       indx = (long)vcl_floor((float)(x[n] + 0.5)) - splineOrder / 2;
516       for (unsigned int k = 0; k <= splineOrder; k++)
517         {
518         evaluateIndex[n][k] = indx++;
519         }
520       }
521     }
522 }
523
524 template <class TImageType, class TCoordRep, class TCoefficientType>
525 void
526 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
527 ::ApplyMirrorBoundaryConditions(vnl_matrix<long> & evaluateIndex, 
528                                 unsigned int splineOrder) const
529 {
530   for (unsigned int n = 0; n < ImageDimension; n++)
531     {
532     long dataLength2 = 2 * m_DataLength[n] - 2;
533
534     // apply the mirror boundary conditions 
535     // TODO:  We could implement other boundary options beside mirror
536     if (m_DataLength[n] == 1)
537       {
538       for (unsigned int k = 0; k <= splineOrder; k++)
539         {
540         evaluateIndex[n][k] = 0;
541         }
542       }
543     else
544       {
545       for (unsigned int k = 0; k <= splineOrder; k++)
546         {
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])
551           {
552           evaluateIndex[n][k] = dataLength2 - evaluateIndex[n][k];
553           }
554         }
555       }
556     }
557 }
558
559 } // namespace itk
560
561 #endif
562
563 //#endif