]> Creatis software - clitk.git/blob - itk/itkBSplineInterpolateImageFunctionWithLUT.txx.original
Put console back for windows
[clitk.git] / itk / itkBSplineInterpolateImageFunctionWithLUT.txx.original
1 #ifndef _ITKINTERPOLATEIMAGEFUNCTIONWITHLUT_
2 #define _ITKINTERPOLATEIMAGEFUNCTIONWITHLUT_
3
4 /* =========================================================================
5                                                                                 
6   @file   itkBSplineInterpolateImageFunctionWithLUT.txx
7   @author David Sarrut <david.sarrut@creatis.insa-lyon.fr>
8
9   Copyright (c) 
10   * CREATIS (Centre de Recherche et d'Applications en Traitement de l'Image). 
11     All rights reserved. See Doc/License.txt or
12     http://www.creatis.insa-lyon.fr/Public/Gdcm/License.html for details.
13   * L�on B�rard cancer center, 28 rue La�nnec, 69373 Lyon cedex 08, France
14   * http://www.creatis.insa-lyon.fr/rio
15                                                                                 
16   This software is distributed WITHOUT ANY WARRANTY; without even the
17   implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
18   PURPOSE.  See the above copyright notices for more information.
19                                                                              
20 ========================================================================= */
21 #include "itkCastImageFilter.h"
22
23 //====================================================================
24 template <class TImageType, class TCoordRep, class TCoefficientType>
25 BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
26 BSplineInterpolateImageFunctionWithLUT():Superclass() {  
27   // Set default values
28   for(int i=0; i<TImageType::ImageDimension; i++) mOutputSpacing[i] = -1;
29   SetLUTSamplingFactor(10);
30   SetSplineOrder(3);
31   mWeightsAreUpToDate = false;
32   // Following need to be pointer beacause values are updated into 
33   // "const" function
34   mIntrinsecError  = new double;
35   mNumberOfError = new long;
36   mIntrinsecErrorMax = new double;
37   mInputIsCoef = false;
38 }
39 //====================================================================
40
41 //====================================================================
42 template <class TImageType, class TCoordRep, class TCoefficientType>
43 void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
44 SetLUTSamplingFactor(int s) {  
45   for(int i=0; i<TImageType::ImageDimension; i++) mSamplingFactors[i] = s;
46   mWeightsAreUpToDate = false;
47 }
48
49 //====================================================================
50
51 //====================================================================
52 template <class TImageType, class TCoordRep, class TCoefficientType>
53 void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
54 SetOutputSpacing(const SpacingType & s) {
55   for(int i=0; i<TImageType::ImageDimension; i++) 
56     mOutputSpacing[i] = s[i];
57   mWeightsAreUpToDate = false;
58 }
59 //====================================================================
60
61 //====================================================================
62 template <class TImageType, class TCoordRep, class TCoefficientType>
63 void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
64 SetSplineOrder(unsigned int SplineOrder) {
65   Superclass::SetSplineOrder(SplineOrder);
66   // Compute support and half size
67   static const int d = TImageType::ImageDimension;
68   for(int l=0; l<d; l++) {
69     mSplineOrders[l]= SplineOrder;
70     mSupport[l] = SplineOrder+1;
71     if (mSupport[l] % 2 == 0) { // support is even
72       mHalfSupport[l] = mSupport[l]/2-1;
73     }
74     else mHalfSupport[l] = mSupport[l]/2; // support is odd (like cubic spline)
75   }
76   mSupportSize = 1;
77   for(int l=0; l<d; l++) {
78     mSupportSize *= mSupport[l];
79   }
80   mWeightsAreUpToDate = false;
81 }
82 //====================================================================
83
84 //JV
85 //====================================================================
86 template <class TImageType, class TCoordRep, class TCoefficientType>
87 void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
88 SetSplineOrders(SizeType SplineOrders) {
89   mSplineOrders=SplineOrders;
90   DD(mSplineOrders); 
91 // Compute support and half size
92   static const int d = TImageType::ImageDimension;
93   for(int l=0; l<d; l++) {
94     mSupport[l] = mSplineOrders[l]+1;
95     if (mSupport[l] % 2 == 0) { // support is even
96       mHalfSupport[l] = mSupport[l]/2-1;
97     }
98     else mHalfSupport[l] = mSupport[l]/2; // support is odd (like cubic spline)
99   }
100   mSupportSize = 1;
101   for(int l=0; l<d; l++) {
102     mSupportSize *= mSupport[l];
103   }
104   mWeightsAreUpToDate = false;
105 }
106 //====================================================================
107
108 //====================================================================
109 template <class TImageType, class TCoordRep, class TCoefficientType>
110 void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
111 SetInputImage(const TImageType * inputData) { 
112
113   // Call superclass SetInputImage and return if input is null
114   // clitk::Timer t;
115   //   t.Start();
116   //ds 
117   // this->InterpolateImageFunction<TImageType,TCoordRep>::SetInputImage(inputData);
118
119   //==============================
120   if (!mInputIsCoef) {
121     Superclass::SetInputImage(inputData);
122     //  t.Stop();
123   }
124   //==============================
125   else {
126     // If inputData is not float/double, we do not want to set
127     // mInputIsCoef to true, but the compiler goes here ...
128
129     if (typeid(typename TImageType::PixelType) != typeid(TCoefficientType)) {
130       std::cerr << "Input should be float or double to be set as 'coefficients' with SetInputImageIsCoefficient(true) and should be the same type than TCoefficientType" << std::endl;
131       exit(0);
132     }
133     else {
134       itkStaticConstMacro(ImageDimension, unsigned int, Superclass::ImageDimension);
135       typedef Image<TCoefficientType, ImageDimension> CoefficientImageType;
136       typedef itk::CastImageFilter< TImageType, CoefficientImageType> CastFilterType;
137       typename CastFilterType::Pointer castFilter=CastFilterType::New();
138       castFilter->SetInput(inputData);
139       castFilter->Update();
140       this->m_Coefficients=castFilter->GetOutput();
141       if ( this->m_Coefficients.IsNotNull())
142         this->m_DataLength = this->m_Coefficients->GetBufferedRegion().GetSize();
143     }
144
145 //     //JV input should be double or float. To keep it generic: cast to the CoefficientImageType
146 //     /** Dimension underlying input image. */
147 //     itkStaticConstMacro(ImageDimension, unsigned int, Superclass::ImageDimension);
148 //     /** Internal Coefficient typedef support */
149 //     DD( ImageDimension);
150 //     typedef Image<TCoefficientType,ImageDimension> CoefficientImageType;
151 //     typedef itk::CastImageFilter< TImageType, CoefficientImageType> CastFilterType;
152 //     typename CastFilterType::Pointer castFilter=CastFilterType::New();
153 //     castFilter->SetInput(inputData);
154 //     castFilter->Update();
155 //      DS this->InterpolateImageFunction<TImageType,TCoordRep>::m_Coefficients = inputData;
156
157 //     //JV cast to the coeff type, don't tell my mother
158 //     this->InterpolateImageFunction<TImageType,TCoordRep>::SetInputImage(inputData);
159 //     this->m_Coefficients=castFilter->GetOutput();
160 //     DS ??????    this->m_Coefficients = inputData;
161
162 //     //   //============================== old not generic solution
163 //     //      this->InterpolateImageFunction<TImageType,TCoordRep>::SetInputImage(inputData);
164 //          // DS this->InterpolateImageFunction<TImageType,TCoordRep>::m_Coefficients = inputData;
165 //     //     this->m_Coefficients = inputData;
166 //     //   //==============================    
167 //     if ( this->m_Coefficients.IsNotNull()) {
168 //       this->m_DataLength = this->m_Coefficients->GetBufferedRegion().GetSize();
169 //     }
170   }
171
172   //   t.Stop();  
173   if (!inputData) return;
174  //  mCoefficientTimer = t;
175
176   // Compute Memory offset inside coefficients images (for looping over coefficients)
177   static const int d = TImageType::ImageDimension;
178   mInputMemoryOffset[0] = 1;
179   for(int l=1; l<d; l++) {
180     mInputMemoryOffset[l] = 
181       mInputMemoryOffset[l-1]*this->m_Coefficients->GetLargestPossibleRegion().GetSize(l-1);
182   }
183   
184   // Compute mSupportOffset according to input size
185   mSupportOffset.resize(mSupportSize);
186   mSupportIndex.resize(mSupportSize);
187   for(int l=0; l<d; l++) mSupportIndex[0][l] = 0;
188   for(unsigned int k=0; k<mSupportSize; k++) {
189     // Get memory offset
190     mSupportOffset[k] = Index2Offset<TImageType::ImageDimension>(mSupportIndex[k], mInputMemoryOffset);
191     // next coefficient index
192     if (k != mSupportSize-1) {
193       for(int l=0; l<d; l++) mSupportIndex[k+1][l] = mSupportIndex[k][l];
194       int l=0;
195       bool stop = false;
196       while (!stop) {
197         mSupportIndex[k+1][l]++;
198         if (mSupportIndex[k+1][l] == (int)mSupport[l]) { //ds supportindex=uint et support=int
199           mSupportIndex[k+1][l] = 0;
200           l++;
201         }
202         else stop = true;
203       }
204     }
205   }
206
207   // Check  
208   for(int l=0; l<d; l++) {
209     if (mOutputSpacing[l] == -1) {
210       std::cerr << "Please use SetOutputSpacing before using bslut interpolator" << std::endl;
211       exit(0);
212     }
213   }
214
215   // Compute BSpline weights if not up to date
216   if (!mWeightsAreUpToDate) UpdatePrecomputedWeights();
217
218   // Initialisation
219   *mIntrinsecErrorMax = *mIntrinsecError = 0.0;
220   *mNumberOfError = 0;  
221 }
222 //====================================================================
223
224 //====================================================================
225 template <class TImageType, class TCoordRep, class TCoefficientType>
226 void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
227 UpdatePrecomputedWeights() {
228 //   mLUTTimer.Reset();
229 //   mLUTTimer.Start();
230   mWeightsCalculator.SetSplineOrders(mSplineOrders);
231   mWeightsCalculator.SetSamplingFactors(mSamplingFactors);
232   mWeightsCalculator.ComputeTensorProducts();
233   mWeightsAreUpToDate = true;
234   //DS
235   //   coef = new TCoefficientType[mSupportSize];   
236   //     mCorrectedSupportIndex.resize(mSupportSize);
237   //     mCorrectedSupportOffset.resize(mSupportSize);
238  //  mLUTTimer.Stop();
239 //   mLUTTimer.Print("LUT      \t");
240 }
241 //====================================================================
242
243 //====================================================================
244 template <class TImageType, class TCoordRep, class TCoefficientType>
245 typename BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::IndexType
246 BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
247 GetSampleIndexOfPixelPosition(const ContinuousIndexType & x, IndexType & EvaluateIndex) const {
248
249   /*
250     WARNING : sometimes, a floating number x could not really be
251     represented in memory. In this case, the difference between x and
252     floor(x) can be almost 1 (instead of 0).  So we take into account
253     such special case, otherwise it could lead to annoying results.
254   */
255   // static const TCoefficientType tiny = 1.0e-7;
256   IndexType index;
257
258   for(int l=0; l<TImageType::ImageDimension; l++) {
259   // bool mChange = false;
260
261   // Compute t1 = distance to floor 
262   TCoefficientType t1 = x[l]- vcl_floor(x[l]);
263     
264   // Compute index in precomputed weights table
265   TCoefficientType t2 = mSamplingFactors[l]*t1;
266   index[l] = (IndexValueType)lrint(t2);
267
268   // For even order : test if too close to 0.5 (but lower). In this
269   // case : take the next coefficient
270   if (!(mSplineOrders[l] & 1)) {
271     if (t1<0.5) {
272       if (mSamplingFactors[l] & 1) {
273         if (index[l] ==  (int) mSamplingFactors[l]/2+1) EvaluateIndex[l] = EvaluateIndex[l]+1;
274       }
275       else if (index[l] == (int) mSamplingFactors[l]/2) EvaluateIndex[l] = EvaluateIndex[l]+1;
276     }
277   }
278
279   // Statistics (to be removed)
280   /*
281   *mIntrinsecError += fabs(index[l]-t2);                    
282   (*mNumberOfError)++;
283   if (fabs(index[l]-t2)> *mIntrinsecErrorMax) *mIntrinsecErrorMax = fabs(index[l]-t2);
284   */
285
286   // When to close to 1, take the next coefficient for ordd order, but
287   // only change index for odd
288     if (index[l] == (int)mSamplingFactors[l]) {
289       index[l] = 0;
290       if (mSplineOrders[l] & 1) EvaluateIndex[l] = EvaluateIndex[l]+1;
291     }
292   }
293
294   // The end
295   return index;
296 }
297 //====================================================================
298
299 //====================================================================
300 template <class TImageType, class TCoordRep, class TCoefficientType>
301 typename BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::OutputType
302 BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
303 EvaluateAtContinuousIndex(const ContinuousIndexType & x) const {
304
305   // For shorter coding
306   static const unsigned int d = TImageType::ImageDimension;
307
308   // Compute the indexe of the first interpolation coefficient in the
309   // coefficient image
310   IndexType EvaluateIndex;
311   long indx;
312   for (unsigned int l=0; l<d; l++)  {
313     if (mSplineOrders[l] & 1) {  // Use this index calculation for odd splineOrder (like cubic)
314       indx = (long)vcl_floor(x[l]) - mSplineOrders[l] / 2 ; //this->m_SplineOrder / 2;
315       EvaluateIndex[l] = indx;
316     }
317     else { // Use this index calculation for even splineOrder
318       if (mSplineOrders[l] == 0) EvaluateIndex[l] = (long)rint(x[l]);
319       else {
320         indx = (long)vcl_floor((x[l]+ 0.5)) - mSplineOrders[l] / 2; //this->m_SplineOrder / 2;
321         EvaluateIndex[l] = indx;
322       }
323     }
324   }
325
326   // Compute index of precomputed weights and get pointer to first weights
327   const IndexType weightIndex = GetSampleIndexOfPixelPosition(x, EvaluateIndex);
328   const TCoefficientType * pweights = mWeightsCalculator.GetFirstTensorProduct(weightIndex);
329
330   // Check boundaries
331   bool mBoundaryCase = false;
332   for (unsigned int l=0; l<d; l++) {
333     if ((EvaluateIndex[l] < 0) ||
334         (EvaluateIndex[l]+mSupport[l]) >= this->m_Coefficients->GetLargestPossibleRegion().GetSize(l)) {
335       mBoundaryCase = true;
336     }
337   }
338
339   // Pointer to support offset
340   const int * psupport;
341     
342   // Special case for boundary (to be changed ....)
343   std::vector<int> mCorrectedSupportOffset;//(mSupportSize);
344   if (mBoundaryCase) {
345     //    return -1000;
346     //std::vector<TCoefficientType> coef(mSupportSize);    
347     // DD(EvaluateIndex);
348      std::vector<IndexType> mCorrectedSupportIndex;//(mSupportSize);
349     mCorrectedSupportIndex.resize(mSupportSize);
350     mCorrectedSupportOffset.resize(mSupportSize);
351     for(unsigned int i=0; i<mSupportSize; i++) {
352       // DD(mSupportIndex[i]);
353       for (unsigned int l=0; l<d; l++) {
354         long a = mSupportIndex[i][l] + EvaluateIndex[l];
355         long b = this->m_Coefficients->GetLargestPossibleRegion().GetSize(l);
356         // DD(a);
357         //      DD(b);
358         long d2 = 2 * b - 2;
359         if (a < 0) {
360           mCorrectedSupportIndex[i][l] = -a - d2*(-a/d2) - EvaluateIndex[l];//mSupportIndex[i][l]-a;
361         }
362         else {
363           if (a>=b) {
364             mCorrectedSupportIndex[i][l] = d2 - a - EvaluateIndex[l];
365           }
366           else {
367             mCorrectedSupportIndex[i][l] = mSupportIndex[i][l]; //a - d2*(a/d2) - EvaluateIndex[l];
368           }
369           /*
370           if (a>=b) {
371             mCorrectedSupportIndex[i][l] = d2 - a - EvaluateIndex[l];//mSupportIndex[i][l] - (a-(b-1));
372           }
373           else {
374             mCorrectedSupportIndex[i][l] = mSupportIndex[i][l];
375           }
376           */
377         }
378       }
379       // DD(mCorrectedSupportIndex[i]);
380       mCorrectedSupportOffset[i] = Index2Offset<TImageType::ImageDimension>(mCorrectedSupportIndex[i], mInputMemoryOffset);
381     }
382     // for (unsigned int l=0; l<d; l++) {
383     //       EvaluateIndex[l] = EvaluateIndex[l] + mCorrectedSupportIndex[0][l];
384     //     }
385     // DD(EvaluateIndex);
386     psupport = &mCorrectedSupportOffset[0];
387   }
388   else {
389     psupport = &mSupportOffset[0];
390   }
391
392   // Get pointer to first coefficient. Even if in some boundary cases,
393   // EvaluateIndex is out of the coefficient image, 
394   const TCoefficientType * pcoef = &(this->m_Coefficients->GetPixel(EvaluateIndex));
395
396   // Main loop over BSpline support
397   TCoefficientType interpolated = 0.0;
398   for (unsigned int p=0; p<mSupportSize; p++) {
399     interpolated += pcoef[*psupport] * (*pweights);
400     ++psupport;
401     ++pweights;
402   }
403
404   // Return interpolated value
405   return(interpolated);    
406 }
407 //====================================================================
408
409 #endif