]> Creatis software - clitk.git/blob - itk/clitkVectorBSplineInterpolateImageFunctionWithLUT.txx
Merge branch 'master' of /home/dsarrut/clitk3.server
[clitk.git] / itk / clitkVectorBSplineInterpolateImageFunctionWithLUT.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 _CLITKVECTORBSPLINEINTERPOLATEIMAGEFUNCTIONWITHLUT_TXX
19 #define _CLITKVECTORBSPLINEINTERPOLATEIMAGEFUNCTIONWITHLUT_TXX
20 /* =========================================================================
21
22 @file   itkBSplineInterpolateImageFunctionWithLUT.txx
23 @author jefvdmb@gmail.com
24
25 Copyright (c)
26 * CREATIS (Centre de Recherche et d'Applications en Traitement de l'Image).
27 All rights reserved. See Doc/License.txt or
28 http://www.creatis.insa-lyon.fr/Public/Gdcm/License.html for details.
29 * Léon Bérard cancer center, 28 rue Laënnec, 69373 Lyon cedex 08, France
30 * http://www.creatis.insa-lyon.fr/rio
31
32 This software is distributed WITHOUT ANY WARRANTY; without even the
33 implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
34 PURPOSE.  See the above copyright notices for more information.
35
36 ========================================================================= */
37 namespace clitk
38 {
39 //====================================================================
40 template <class TImageType, class TCoordRep, class TCoefficientType>
41 VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
42 VectorBSplineInterpolateImageFunctionWithLUT():Superclass()
43 {
44   // Set default values
45   SetLUTSamplingFactor(20);
46   SetSplineOrder(3);
47   mWeightsAreUpToDate = false;
48 }
49
50 //====================================================================
51
52 //====================================================================
53 template <class TImageType, class TCoordRep, class TCoefficientType>
54 void VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
55 SetLUTSamplingFactor(const int& s)
56 {
57   for(int i=0; i<TImageType::ImageDimension; i++) mSamplingFactors[i] = s;
58   mWeightsAreUpToDate = false;
59 }
60
61 //====================================================================
62
63 //====================================================================
64 template <class TImageType, class TCoordRep, class TCoefficientType>
65 void VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
66 SetLUTSamplingFactors(const SizeType& s)
67 {
68   for(int i=0; i<TImageType::ImageDimension; i++) mSamplingFactors[i] = s[i];
69   mWeightsAreUpToDate = false;
70 }
71
72 //====================================================================
73
74 //====================================================================
75 template <class TImageType, class TCoordRep, class TCoefficientType>
76 void VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
77 SetSplineOrder(const unsigned int& SplineOrder)
78 {
79   Superclass::SetSplineOrder(SplineOrder);
80   // Compute support and half size
81   static const int d = TImageType::ImageDimension;
82   for(int l=0; l<d; l++) {
83     mSplineOrders[l]= SplineOrder;
84     mSupport[l] = SplineOrder+1;
85     if (mSupport[l] % 2 == 0) { // support is even
86       mHalfSupport[l] = mSupport[l]/2-1;
87     } else mHalfSupport[l] = mSupport[l]/2; // support is odd (like cubic spline)
88   }
89   mSupportSize = 1;
90   for(int l=0; l<d; l++) {
91     mSupportSize *= mSupport[l];
92   }
93   mWeightsAreUpToDate = false;
94 }
95 //====================================================================
96
97 //JV
98 //====================================================================
99 template <class TImageType, class TCoordRep, class TCoefficientType>
100 void VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
101 SetSplineOrders(const SizeType& SplineOrders)
102 {
103   mSplineOrders=SplineOrders;
104
105   // Compute support and half size
106   static const int d = TImageType::ImageDimension;
107   for(int l=0; l<d; l++) {
108     mSupport[l] = mSplineOrders[l]+1;
109     if (mSupport[l] % 2 == 0) { // support is even
110       mHalfSupport[l] = mSupport[l]/2-1;
111     } else mHalfSupport[l] = mSupport[l]/2; // support is odd (like cubic spline)
112   }
113   mSupportSize = 1;
114   for(int l=0; l<d; l++) {
115     mSupportSize *= mSupport[l];
116   }
117   mWeightsAreUpToDate = false;
118 }
119 //====================================================================
120
121 //====================================================================
122 template <class TImageType, class TCoordRep, class TCoefficientType>
123 void VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
124 SetInputImage(const TImageType * inputData)
125 {
126   // JV  Call superclass (decomposition filter is executeed each time!)
127   // JV Should call clitkVectorBSplineDecompositionFilterWithOBD to allow different order by dimension
128   Superclass::SetInputImage(inputData);
129     
130   // Update the weightproperties
131   if (!inputData) return;
132   UpdateWeightsProperties();
133
134 }
135
136 //====================================================================
137 template <class TImageType, class TCoordRep, class TCoefficientType>
138 void VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
139 UpdateWeightsProperties()
140 {
141
142   // Compute Memory offset inside coefficients images (for looping over coefficients)
143   static const unsigned int d = TImageType::ImageDimension;
144   mInputMemoryOffset[0] = 1;
145   for(unsigned int l=1; l<d; l++) {
146     mInputMemoryOffset[l] =
147       mInputMemoryOffset[l-1]*this->m_Coefficients->GetLargestPossibleRegion().GetSize(l-1);
148   }
149
150     // Compute mSupportOffset according to input size
151     mSupportOffset.resize(mSupportSize);
152     mSupportIndex.resize(mSupportSize);
153     for(unsigned int l=0; l<d; l++) mSupportIndex[0][l] = 0;
154     for(unsigned int k=0; k<mSupportSize; k++) {
155       // Get memory offset
156       mSupportOffset[k] = itk::Index2Offset<TImageType::ImageDimension>(mSupportIndex[k], mInputMemoryOffset);
157       // next coefficient index
158       if (k != mSupportSize-1) {
159         for(unsigned int l=0; l<d; l++) mSupportIndex[k+1][l] = mSupportIndex[k][l];
160         int l=0;
161         bool stop = false;
162         while (!stop) {
163           mSupportIndex[k+1][l]++;
164           if (static_cast<unsigned int>(mSupportIndex[k+1][l]) == mSupport[l]) {
165             mSupportIndex[k+1][l] = 0;
166             l++;
167           }
168         else stop = true;
169         }
170       }
171     }
172
173     // Compute BSpline weights if not up to date
174     if (!mWeightsAreUpToDate)
175       UpdatePrecomputedWeights();
176   }
177 //====================================================================
178
179 //====================================================================
180 template <class TImageType, class TCoordRep, class TCoefficientType>
181 void VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
182 UpdatePrecomputedWeights()
183 {
184   mWeightsCalculator.SetSplineOrders(mSplineOrders);
185   mWeightsCalculator.SetSamplingFactors(mSamplingFactors);
186   mWeightsCalculator.ComputeTensorProducts();
187   mWeightsAreUpToDate = true;
188 }
189 //====================================================================
190
191 //====================================================================
192 template <class TImageType, class TCoordRep, class TCoefficientType>
193 typename VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::IndexType
194 VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
195 GetSampleIndexOfPixelPosition(const ContinuousIndexType & x, IndexType & EvaluateIndex) const
196 {
197
198   /*
199     WARNING : sometimes, a floating number x could not really be
200     represented in memory. In this case, the difference between x and
201     floor(x) can be almost 1 (instead of 0).  So we take into account
202     such special case, otherwise it could lead to annoying results.
203   */
204   //  static const TCoefficientType tiny = 1.0e-7;
205   IndexType index;
206
207   for(int l=0; l<TImageType::ImageDimension; l++) {
208     // Compute t1 = distance to floor
209     TCoefficientType t1 = x[l]- vcl_floor(x[l]);
210
211     // Compute index in precomputed weights table
212     TCoefficientType t2 = mSamplingFactors[l]*t1;
213     index[l] = (IndexValueType)lrint(t2);
214
215     // For even order : test if too close to 0.5 (but lower). In this
216     // case : take the next coefficient
217     if (!(mSplineOrders[l] & 1)) {
218       if (t1<0.5) {
219         if (mSamplingFactors[l] & 1) {
220           if (index[l] ==  (int) mSamplingFactors[l]/2+1) EvaluateIndex[l] = EvaluateIndex[l]+1;
221         }
222
223         else if (index[l] == (int) mSamplingFactors[l]/2) EvaluateIndex[l] = EvaluateIndex[l]+1;
224       }
225     }
226
227     // When to close to 1, take the next coefficient for odd order, but
228     // only change index for odd
229     if (index[l] == (int)mSamplingFactors[l]) {
230       index[l] = 0;
231       if (mSplineOrders[l] & 1) EvaluateIndex[l] = EvaluateIndex[l]+1;
232     }
233   }
234
235   // The end
236   return index;
237 }
238
239
240 //====================================================================
241
242 //====================================================================
243 template <class TImageType, class TCoordRep, class TCoefficientType>
244 typename VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::OutputType
245 VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
246 EvaluateAtContinuousIndex(const ContinuousIndexType & x) const
247 {
248
249     // JV Compute BSpline weights if not up to date! Problem const: pass image as last
250     //  if (!mWeightsAreUpToDate) UpdatePrecomputedWeights();
251   // For shorter coding
252   static const unsigned int d = TImageType::ImageDimension;
253
254   // Compute the index of the first interpolation coefficient in the coefficient image
255   IndexType evaluateIndex;
256   long indx;
257   for (unsigned int l=0; l<d; l++)  {
258     if (mSplineOrders[l] & 1) {  // Use this index calculation for odd splineOrder (like cubic)
259       indx = (long)vcl_floor(x[l]) - mSplineOrders[l] / 2 ; //this->m_SplineOrder / 2;
260       evaluateIndex[l] = indx;
261     } else { // Use this index calculation for even splineOrder
262       if (mSplineOrders[l] == 0) evaluateIndex[l] = (long)rint(x[l]);
263       else {
264         indx = (long)vcl_floor((x[l]+ 0.5)) - mSplineOrders[l] / 2; //this->m_SplineOrder / 2;
265         evaluateIndex[l] = indx;
266       }
267     }
268   }
269
270   // Compute index of precomputed weights and get pointer to first weights
271   const IndexType weightIndex = GetSampleIndexOfPixelPosition(x, evaluateIndex);
272   const TCoefficientType * pweights = mWeightsCalculator.GetFirstTensorProduct(weightIndex);
273
274   // Check boundaries
275   bool boundaryCase = false;
276   for (unsigned int l=0; l<d; l++) {
277     if ((evaluateIndex[l] < 0) ||
278         (evaluateIndex[l]+mSupport[l]) >= this->m_Coefficients->GetLargestPossibleRegion().GetSize(l)) {
279       boundaryCase = true;
280     }
281   }
282
283   // Pointer to support offset
284   const int * psupport;
285
286   // Special case for boundary (to be changed ....)
287   std::vector<int> correctedSupportOffset;
288   if (boundaryCase) {
289     std::vector<IndexType> correctedSupportIndex;//(mSupportSize);
290     correctedSupportIndex.resize(mSupportSize);
291     correctedSupportOffset.resize(mSupportSize);
292     for(unsigned int i=0; i<mSupportSize; i++) {
293       for (unsigned int l=0; l<d; l++) {
294         long a = mSupportIndex[i][l] + evaluateIndex[l];
295         long b = this->m_Coefficients->GetLargestPossibleRegion().GetSize(l);
296         long d2 = 2 * b - 2;
297         if (a < 0) {
298           correctedSupportIndex[i][l] = -a - d2*(-a/d2) - evaluateIndex[l];//mSupportIndex[i][l]-a;
299         } else {
300           if (a>=b) {
301             correctedSupportIndex[i][l] = d2 - a - evaluateIndex[l];
302           } else {
303             correctedSupportIndex[i][l] = mSupportIndex[i][l]; //a - d2*(a/d2) - EvaluateIndex[l];
304           }
305           /*
306             if (a>=b) {
307             correctedSupportIndex[i][l] = d2 - a - EvaluateIndex[l];//mSupportIndex[i][l] - (a-(b-1));
308             }
309             else {
310             correctedSupportIndex[i][l] = mSupportIndex[i][l];
311             }
312           */
313         }
314       }
315       // DD(correctedSupportIndex[i]);
316       correctedSupportOffset[i] = itk::Index2Offset<TImageType::ImageDimension>(correctedSupportIndex[i], mInputMemoryOffset);
317     }
318     psupport = &correctedSupportOffset[0];
319   } else {
320     psupport = &mSupportOffset[0];
321   }
322
323   // Get pointer to first coefficient. Even if in some boundary cases,
324   // EvaluateIndex is out of the coefficient image,
325   //JV
326   const CoefficientImagePixelType * pcoef = &(this->m_Coefficients->GetPixel(evaluateIndex));
327
328   // Main loop over BSpline support
329   CoefficientImagePixelType interpolated = itk::NumericTraits<CoefficientImagePixelType>::Zero;;
330   for (unsigned int p=0; p<mSupportSize; p++) {
331     interpolated += pcoef[*psupport] * (*pweights);
332     ++psupport;
333     ++pweights;
334   }
335
336   // Return interpolated value
337   return(interpolated);
338 }
339 //====================================================================
340
341
342 //====================================================================
343 template <class TImageType, class TCoordRep, class TCoefficientType>
344 void
345 VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
346 EvaluateWeightsAtContinuousIndex(const ContinuousIndexType&  x,  const TCoefficientType ** pweights, IndexType & evaluateIndex)const
347 {
348
349   // JV Compute BSpline weights if not up to date! Problem const: pass image as last
350   //  if (!mWeightsAreUpToDate) UpdatePrecomputedWeights();
351
352   // For shorter coding
353   static const unsigned int d = TImageType::ImageDimension;
354
355   // Compute the index of the first interpolation coefficient in the coefficient image
356   long indx;
357   for (unsigned int l=0; l<d; l++)  {
358     if (mSplineOrders[l] & 1) {  // Use this index calculation for odd splineOrder (like cubic)
359       indx = (long)vcl_floor(x[l]) - mSplineOrders[l] / 2 ; //this->m_SplineOrder / 2;
360       evaluateIndex[l] = indx;
361     } else { // Use this index calculation for even splineOrder
362       if (mSplineOrders[l] == 0) evaluateIndex[l] = (long)rint(x[l]);
363       else {
364         indx = (long)vcl_floor((x[l]+ 0.5)) - mSplineOrders[l] / 2; //this->m_SplineOrder / 2;
365         evaluateIndex[l] = indx;
366       }
367     }
368   }
369
370   // Compute index of precomputed weights and get pointer to first weights
371   const IndexType weightIndex = GetSampleIndexOfPixelPosition(x, evaluateIndex);
372   *pweights = mWeightsCalculator.GetFirstTensorProduct(weightIndex);
373
374 }
375 //====================================================================
376
377
378
379
380 } //end namespace
381
382 #endif //_CLITKVECTORBSPLINEINTERPOLATEIMAGEFUNCTIONWITHLUT_TXX