]> Creatis software - clitk.git/blob - itk/itkBSplineInterpolateImageFunctionWithLUT.txx
Moved from repository clitk to clitk.private/tests_dav
[clitk.git] / itk / itkBSplineInterpolateImageFunctionWithLUT.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 _ITKINTERPOLATEIMAGEFUNCTIONWITHLUT_TXX
19 #define _ITKINTERPOLATEIMAGEFUNCTIONWITHLUT_TXX
20 #include <itkMath.h>
21
22 namespace itk
23 {
24   //====================================================================
25   template <class TImageType, class TCoordRep, class TCoefficientType>
26   BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
27   BSplineInterpolateImageFunctionWithLUT():Superclass() {  
28     
29     // Set default values
30     SetLUTSamplingFactor(20);
31     SetSplineOrder(3);
32     mWeightsAreUpToDate = false;
33   }
34
35   //====================================================================
36
37   //====================================================================
38   template <class TImageType, class TCoordRep, class TCoefficientType>
39   void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
40   SetLUTSamplingFactor(const int& s) {  
41     for(int i=0; i<TImageType::ImageDimension; i++) mSamplingFactors[i] = s;
42     mWeightsAreUpToDate = false;
43   }
44
45   //====================================================================
46
47   //====================================================================
48   template <class TImageType, class TCoordRep, class TCoefficientType>
49   void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
50   SetLUTSamplingFactors(const SizeType& s) {  
51     for(int i=0; i<TImageType::ImageDimension; i++) mSamplingFactors[i] = s[i];
52     mWeightsAreUpToDate = false;
53   }
54
55   //====================================================================
56
57   //====================================================================
58   template <class TImageType, class TCoordRep, class TCoefficientType>
59   void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
60   SetSplineOrder(const unsigned int& SplineOrder) {
61     Superclass::SetSplineOrder(SplineOrder);
62     // Compute support and half size
63     static const int d = TImageType::ImageDimension;
64     for(int l=0; l<d; l++) {
65       mSplineOrders[l]= SplineOrder;
66       mSupport[l] = SplineOrder+1;
67       if (mSupport[l] % 2 == 0) { // support is even
68         mHalfSupport[l] = mSupport[l]/2-1;
69       }
70       else mHalfSupport[l] = mSupport[l]/2; // support is odd (like cubic spline)
71     }
72     mSupportSize = 1;
73     for(int l=0; l<d; l++) {
74       mSupportSize *= mSupport[l];
75     }
76     mWeightsAreUpToDate = false;
77   }
78   //====================================================================
79
80   //JV
81   //====================================================================
82   template <class TImageType, class TCoordRep, class TCoefficientType>
83   void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
84   SetSplineOrders(const SizeType& SplineOrders) {
85     mSplineOrders=SplineOrders;
86
87     // Compute support and half size
88     static const int d = TImageType::ImageDimension;
89     for(int l=0; l<d; l++) {
90       mSupport[l] = mSplineOrders[l]+1;
91       if (mSupport[l] % 2 == 0) { // support is even
92         mHalfSupport[l] = mSupport[l]/2-1;
93       }
94       else mHalfSupport[l] = mSupport[l]/2; // support is odd (like cubic spline)
95     }
96     mSupportSize = 1;
97     for(int l=0; l<d; l++) {
98       mSupportSize *= mSupport[l];
99     }
100     mWeightsAreUpToDate = false;
101   }
102   //====================================================================
103
104   //====================================================================
105   template <class TImageType, class TCoordRep, class TCoefficientType>
106   void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
107   SetInputImage(const TImageType * inputData) { 
108
109     // JV  Call superclass (decomposition filter is executeed each time!)
110     // JV Should call itkBSplineDecompositionFilterWithOBD to allow different order by dimension
111     Superclass::SetInputImage(inputData);
112     
113     // Update the weightproperties
114     if (!inputData) return;
115     UpdateWeightsProperties();
116   }
117   
118   //====================================================================
119   template <class TImageType, class TCoordRep, class TCoefficientType>
120   void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
121   UpdateWeightsProperties() {
122     
123     // Compute Memory offset inside coefficients images (for looping over coefficients)
124     static const unsigned int d = TImageType::ImageDimension;
125     mInputMemoryOffset[0] = 1;
126     for(unsigned int l=1; l<d; l++) {
127       mInputMemoryOffset[l] = 
128         mInputMemoryOffset[l-1]*this->m_Coefficients->GetLargestPossibleRegion().GetSize(l-1);
129     }
130     
131     // Compute mSupportOffset according to input size
132     mSupportOffset.resize(mSupportSize);
133     mSupportIndex.resize(mSupportSize);
134     for(unsigned int l=0; l<d; l++) mSupportIndex[0][l] = 0;
135     for(unsigned int k=0; k<mSupportSize; k++) {
136       // Get memory offset
137       mSupportOffset[k] = Index2Offset<TImageType::ImageDimension>(mSupportIndex[k], mInputMemoryOffset);
138       // next coefficient index
139       if (k != mSupportSize-1) {
140         for(unsigned int l=0; l<d; l++) mSupportIndex[k+1][l] = mSupportIndex[k][l];
141         int l=0;
142         bool stop = false;
143         while (!stop) {
144           mSupportIndex[k+1][l]++;
145           if (static_cast<unsigned int>(mSupportIndex[k+1][l]) == mSupport[l]) {
146             mSupportIndex[k+1][l] = 0;
147             l++;
148           }
149               else stop = true;
150         }
151       }
152     }
153     
154     // Compute BSpline weights if not up to date
155     if (!mWeightsAreUpToDate) 
156       UpdatePrecomputedWeights();
157   }
158   //====================================================================
159
160   //====================================================================
161   template <class TImageType, class TCoordRep, class TCoefficientType>
162   void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
163   UpdatePrecomputedWeights() {
164
165     mWeightsCalculator.SetSplineOrders(mSplineOrders);
166     mWeightsCalculator.SetSamplingFactors(mSamplingFactors);
167     mWeightsCalculator.ComputeTensorProducts();
168     mWeightsAreUpToDate = true;
169   }
170   //====================================================================
171
172   //====================================================================
173   template <class TImageType, class TCoordRep, class TCoefficientType>
174   typename BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::IndexType
175   BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
176   GetSampleIndexOfPixelPosition(const ContinuousIndexType & x, IndexType & EvaluateIndex) const {
177
178     /*
179       WARNING : sometimes, a floating number x could not really be
180       represented in memory. In this case, the difference between x and
181       floor(x) can be almost 1 (instead of 0).  So we take into account
182       such special case, otherwise it could lead to annoying results.
183     */
184     //  static const TCoefficientType tiny = 1.0e-7;
185     IndexType index;
186
187     for(int l=0; l<TImageType::ImageDimension; l++) 
188       {
189         // Compute t1 = distance to floor 
190         TCoefficientType t1 = x[l]- vcl_floor(x[l]);
191     
192         // Compute index in precomputed weights table
193         TCoefficientType t2 = mSamplingFactors[l]*t1;
194         index[l] = (IndexValueType)Math::Round<TCoefficientType>(t2);
195
196         // For even order : test if too close to 0.5 (but lower). In this
197         // case : take the next coefficient
198         if (!(mSplineOrders[l] & 1)) {
199           if (t1<0.5) {
200             if (mSamplingFactors[l] & 1) {
201               if (index[l] ==  (int) mSamplingFactors[l]/2+1) EvaluateIndex[l] = EvaluateIndex[l]+1;
202             }
203
204             else if (index[l] == (int) mSamplingFactors[l]/2) EvaluateIndex[l] = EvaluateIndex[l]+1;
205           }
206         }
207
208         // When to close to 1, take the next coefficient for odd order, but
209         // only change index for odd
210         if (index[l] == (int)mSamplingFactors[l]) {
211           index[l] = 0;
212           if (mSplineOrders[l] & 1) EvaluateIndex[l] = EvaluateIndex[l]+1;
213         }
214       }
215
216     // The end
217     return index;
218   }
219
220
221   //====================================================================
222
223   //====================================================================
224   template <class TImageType, class TCoordRep, class TCoefficientType>
225   typename BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::OutputType
226   BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
227   EvaluateAtContinuousIndex(const ContinuousIndexType & x) const {
228
229     // JV Compute BSpline weights if not up to date! Problem const: pass image as last
230     //  if (!mWeightsAreUpToDate) UpdatePrecomputedWeights();
231
232     // For shorter coding
233     static const unsigned int d = TImageType::ImageDimension;
234
235     // Compute the index of the first interpolation coefficient in the coefficient image
236     IndexType evaluateIndex;
237     long indx;
238     for (unsigned int l=0; l<d; l++)  {
239       if (mSplineOrders[l] & 1) {  // Use this index calculation for odd splineOrder (like cubic)
240         indx = (long)vcl_floor(x[l]) - mSplineOrders[l] / 2 ; //this->m_SplineOrder / 2;
241         evaluateIndex[l] = indx;
242       }
243       else { // Use this index calculation for even splineOrder
244         if (mSplineOrders[l] == 0) evaluateIndex[l] = Math::Round<typename ContinuousIndexType::ValueType>(x[l]);
245         else {
246           indx = (long)vcl_floor((x[l]+ 0.5)) - mSplineOrders[l] / 2; //this->m_SplineOrder / 2;
247           evaluateIndex[l] = indx;
248         }
249       }
250     }
251   
252     // Compute index of precomputed weights and get pointer to first weights
253     const IndexType weightIndex = GetSampleIndexOfPixelPosition(x, evaluateIndex);
254     const TCoefficientType * pweights = mWeightsCalculator.GetFirstTensorProduct(weightIndex);
255
256     // Check boundaries
257     bool boundaryCase = false;
258     for (unsigned int l=0; l<d; l++) {
259       if ((evaluateIndex[l] < 0) ||
260           (evaluateIndex[l]+mSupport[l]) >= this->m_Coefficients->GetLargestPossibleRegion().GetSize(l)) {
261         boundaryCase = true;
262       }
263     }
264
265     // Pointer to support offset
266     const int * psupport;
267     
268     // Special case for boundary (to be changed ....)
269     std::vector<int> correctedSupportOffset;
270     if (boundaryCase) {
271     
272       std::vector<IndexType> correctedSupportIndex;//(mSupportSize);
273       correctedSupportIndex.resize(mSupportSize);
274       correctedSupportOffset.resize(mSupportSize);
275       for(unsigned int i=0; i<mSupportSize; i++) {
276         for (unsigned int l=0; l<d; l++) {
277           long a = mSupportIndex[i][l] + evaluateIndex[l];
278           long b = this->m_Coefficients->GetLargestPossibleRegion().GetSize(l);
279           long d2 = 2 * b - 2;
280           if (a < 0) {
281             correctedSupportIndex[i][l] = -a - d2*(-a/d2) - evaluateIndex[l];//mSupportIndex[i][l]-a;
282           }
283           else {
284             if (a>=b) {
285               correctedSupportIndex[i][l] = d2 - a - evaluateIndex[l];
286             }
287             else {
288               correctedSupportIndex[i][l] = mSupportIndex[i][l]; //a - d2*(a/d2) - EvaluateIndex[l];
289             }
290             /*
291               if (a>=b) {
292               correctedSupportIndex[i][l] = d2 - a - EvaluateIndex[l];//mSupportIndex[i][l] - (a-(b-1));
293               }
294               else {
295               correctedSupportIndex[i][l] = mSupportIndex[i][l];
296               }
297             */
298           }
299         }
300         correctedSupportOffset[i] = itk::Index2Offset<TImageType::ImageDimension>(correctedSupportIndex[i], mInputMemoryOffset);
301       }
302       psupport = &correctedSupportOffset[0];
303     }
304     else {
305       psupport = &mSupportOffset[0];
306     }
307
308     // Get pointer to first coefficient. Even if in some boundary cases,
309     // EvaluateIndex is out of the coefficient image, 
310     const TCoefficientType * pcoef = &(this->m_Coefficients->GetPixel(evaluateIndex));
311
312     // Main loop over BSpline support
313     TCoefficientType interpolated = itk::NumericTraits<TCoefficientType>::Zero;
314     for (unsigned int p=0; p<mSupportSize; p++) {
315       interpolated += pcoef[*psupport] * (*pweights);
316       ++psupport;
317       ++pweights;
318     }
319
320     // Return interpolated value
321     return(interpolated);    
322   }
323   //====================================================================
324
325 }
326 #endif //_ITKINTERPOLATEIMAGEFUNCTIONWITHLUT_TXX