1 /*=========================================================================
2 Program: vv http://www.creatis.insa-lyon.fr/rio/vv
5 - University of LYON http://www.universite-lyon.fr/
6 - Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr
7 - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr
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.
13 It is distributed under dual licence
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 /* =========================================================================
22 @file itkBSplineInterpolateImageFunctionWithLUT.txx
23 @author David Sarrut <david.sarrut@creatis.insa-lyon.fr>
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
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.
36 ========================================================================= */
39 //====================================================================
40 template <class TImageType, class TCoordRep, class TCoefficientType>
41 BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
42 BSplineInterpolateImageFunctionWithLUT():Superclass()
45 //for(int i=0; i<TImageType::ImageDimension; i++) mOutputSpacing[i] = -1;
46 SetLUTSamplingFactor(20);
48 mWeightsAreUpToDate = false;
49 // Following need to be pointer beacause values are updated into
51 // mIntrinsecError = new double;
52 // mNumberOfError = new long;
53 // mIntrinsecErrorMax = new double;
54 // mInputIsCoef = false;
57 //====================================================================
59 //====================================================================
60 template <class TImageType, class TCoordRep, class TCoefficientType>
61 void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
62 SetLUTSamplingFactor(const int& s)
64 for(int i=0; i<TImageType::ImageDimension; i++) mSamplingFactors[i] = s;
65 mWeightsAreUpToDate = false;
68 //====================================================================
70 //====================================================================
71 template <class TImageType, class TCoordRep, class TCoefficientType>
72 void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
73 SetLUTSamplingFactors(const SizeType& s)
75 for(int i=0; i<TImageType::ImageDimension; i++) mSamplingFactors[i] = s[i];
76 mWeightsAreUpToDate = false;
79 //====================================================================
81 // //====================================================================
82 // template <class TImageType, class TCoordRep, class TCoefficientType>
83 // void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
84 // SetOutputSpacing(const SpacingType & s) {
85 // for(int i=0; i<TImageType::ImageDimension; i++)
86 // mOutputSpacing[i] = s[i];
87 // // mWeightsAreUpToDate = false;
89 //====================================================================
91 //====================================================================
92 template <class TImageType, class TCoordRep, class TCoefficientType>
93 void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
94 SetSplineOrder(const unsigned int& SplineOrder)
96 Superclass::SetSplineOrder(SplineOrder);
97 // Compute support and half size
98 static const int d = TImageType::ImageDimension;
99 for(int l=0; l<d; l++) {
100 mSplineOrders[l]= SplineOrder;
101 mSupport[l] = SplineOrder+1;
102 if (mSupport[l] % 2 == 0) { // support is even
103 mHalfSupport[l] = mSupport[l]/2-1;
104 } else mHalfSupport[l] = mSupport[l]/2; // support is odd (like cubic spline)
107 for(int l=0; l<d; l++) {
108 mSupportSize *= mSupport[l];
110 mWeightsAreUpToDate = false;
112 //====================================================================
115 //====================================================================
116 template <class TImageType, class TCoordRep, class TCoefficientType>
117 void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
118 SetSplineOrders(const SizeType& SplineOrders)
120 mSplineOrders=SplineOrders;
122 // Compute support and half size
123 static const int d = TImageType::ImageDimension;
124 for(int l=0; l<d; l++) {
125 mSupport[l] = mSplineOrders[l]+1;
126 if (mSupport[l] % 2 == 0) { // support is even
127 mHalfSupport[l] = mSupport[l]/2-1;
128 } else mHalfSupport[l] = mSupport[l]/2; // support is odd (like cubic spline)
131 for(int l=0; l<d; l++) {
132 mSupportSize *= mSupport[l];
134 mWeightsAreUpToDate = false;
136 //====================================================================
138 //====================================================================
139 template <class TImageType, class TCoordRep, class TCoefficientType>
140 void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
141 SetInputImage(const TImageType * inputData)
144 //==============================
145 // if (!mInputIsCoef)
147 Superclass::SetInputImage(inputData);
150 //==============================
151 // //JV Don't call superclass (decomposition filter is executeed each time!)
154 // this->m_Coefficients = inputData;
155 // if ( this->m_Coefficients.IsNotNull())
157 // this->m_DataLength = this->m_Coefficients->GetBufferedRegion().GetSize();
160 // //Call super-superclass in case more input arrives
161 // itk::ImageFunction<TImageType, ITK_TYPENAME itk::NumericTraits<typename TImageType::PixelType>::RealType, TCoefficientType>::SetInputImage(inputData);
163 if (!inputData) return;
164 UpdateWeightsProperties();
168 //====================================================================
169 template <class TImageType, class TCoordRep, class TCoefficientType>
170 void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
171 UpdateWeightsProperties()
174 // Compute Memory offset inside coefficients images (for looping over coefficients)
175 static const unsigned int d = TImageType::ImageDimension;
176 mInputMemoryOffset[0] = 1;
177 for(unsigned int l=1; l<d; l++) {
178 mInputMemoryOffset[l] =
179 mInputMemoryOffset[l-1]*this->m_Coefficients->GetLargestPossibleRegion().GetSize(l-1);
183 if (!mWeightsAreUpToDate) {
184 // Compute mSupportOffset according to input size
185 mSupportOffset.resize(mSupportSize);
186 mSupportIndex.resize(mSupportSize);
187 for(unsigned int l=0; l<d; l++) mSupportIndex[0][l] = 0;
188 for(unsigned int k=0; k<mSupportSize; k++) {
190 mSupportOffset[k] = Index2Offset<TImageType::ImageDimension>(mSupportIndex[k], mInputMemoryOffset);
191 // next coefficient index
192 if (k != mSupportSize-1) {
193 for(unsigned int l=0; l<d; l++) mSupportIndex[k+1][l] = mSupportIndex[k][l];
197 mSupportIndex[k+1][l]++;
198 if (static_cast<unsigned int>(mSupportIndex[k+1][l]) == mSupport[l]) {
199 mSupportIndex[k+1][l] = 0;
207 // for(unsigned int l=0; l<d; l++) {
208 // if (mOutputSpacing[l] == -1) {
209 // std::cerr << "Please use SetOutputSpacing before using BLUT interpolator" << std::endl;
214 // Compute BSpline weights if not up to date
215 //if (!mWeightsAreUpToDate)
216 UpdatePrecomputedWeights();
219 // *mIntrinsecErrorMax = *mIntrinsecError = 0.0;
220 // *mNumberOfError = 0;
222 //====================================================================
224 //====================================================================
225 template <class TImageType, class TCoordRep, class TCoefficientType>
226 void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
227 UpdatePrecomputedWeights()
229 // mLUTTimer.Reset();
230 // mLUTTimer.Start();
231 mWeightsCalculator.SetSplineOrders(mSplineOrders);
232 mWeightsCalculator.SetSamplingFactors(mSamplingFactors);
233 mWeightsCalculator.ComputeTensorProducts();
234 mWeightsAreUpToDate = true;
236 // coef = new TCoefficientType[mSupportSize];
237 // mCorrectedSupportIndex.resize(mSupportSize);
238 // mCorrectedSupportOffset.resize(mSupportSize);
240 // mLUTTimer.Print("LUT \t");
242 //====================================================================
244 //====================================================================
245 template <class TImageType, class TCoordRep, class TCoefficientType>
246 typename BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::IndexType
247 BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
248 GetSampleIndexOfPixelPosition(const ContinuousIndexType & x, IndexType & EvaluateIndex) const
252 WARNING : sometimes, a floating number x could not really be
253 represented in memory. In this case, the difference between x and
254 floor(x) can be almost 1 (instead of 0). So we take into account
255 such special case, otherwise it could lead to annoying results.
257 // static const TCoefficientType tiny = 1.0e-7;
260 for(int l=0; l<TImageType::ImageDimension; l++) {
261 // bool mChange = false;
263 // Compute t1 = distance to floor
264 TCoefficientType t1 = x[l]- vcl_floor(x[l]);
266 // Compute index in precomputed weights table
267 TCoefficientType t2 = mSamplingFactors[l]*t1;
268 index[l] = (IndexValueType)lrint(t2);
270 // For even order : test if too close to 0.5 (but lower). In this
271 // case : take the next coefficient
272 if (!(mSplineOrders[l] & 1)) {
274 if (mSamplingFactors[l] & 1) {
275 if (index[l] == (int) mSamplingFactors[l]/2+1) EvaluateIndex[l] = EvaluateIndex[l]+1;
278 else if (index[l] == (int) mSamplingFactors[l]/2) EvaluateIndex[l] = EvaluateIndex[l]+1;
282 // Statistics (to be removed)
284 *mIntrinsecError += fabs(index[l]-t2);
286 if (fabs(index[l]-t2)> *mIntrinsecErrorMax) *mIntrinsecErrorMax = fabs(index[l]-t2);
289 // When to close to 1, take the next coefficient for odd order, but
290 // only change index for odd
291 if (index[l] == (int)mSamplingFactors[l]) {
293 if (mSplineOrders[l] & 1) EvaluateIndex[l] = EvaluateIndex[l]+1;
302 //====================================================================
304 //====================================================================
305 template <class TImageType, class TCoordRep, class TCoefficientType>
306 typename BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::OutputType
307 BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
308 EvaluateAtContinuousIndex(const ContinuousIndexType & x) const
311 // For shorter coding
312 static const unsigned int d = TImageType::ImageDimension;
314 // Compute the index of the first interpolation coefficient in the coefficient image
315 IndexType evaluateIndex;
317 for (unsigned int l=0; l<d; l++) {
318 if (mSplineOrders[l] & 1) { // Use this index calculation for odd splineOrder (like cubic)
319 indx = (long)vcl_floor(x[l]) - mSplineOrders[l] / 2 ; //this->m_SplineOrder / 2;
320 evaluateIndex[l] = indx;
321 } else { // Use this index calculation for even splineOrder
322 if (mSplineOrders[l] == 0) evaluateIndex[l] = (long)rint(x[l]);
324 indx = (long)vcl_floor((x[l]+ 0.5)) - mSplineOrders[l] / 2; //this->m_SplineOrder / 2;
325 evaluateIndex[l] = indx;
330 // Compute index of precomputed weights and get pointer to first weights
331 const IndexType weightIndex = GetSampleIndexOfPixelPosition(x, evaluateIndex);
332 const TCoefficientType * pweights = mWeightsCalculator.GetFirstTensorProduct(weightIndex);
335 bool boundaryCase = false;
336 for (unsigned int l=0; l<d; l++) {
337 if ((evaluateIndex[l] < 0) ||
338 (evaluateIndex[l]+mSupport[l]) >= this->m_Coefficients->GetLargestPossibleRegion().GetSize(l)) {
343 // Pointer to support offset
344 const int * psupport;
346 // Special case for boundary (to be changed ....)
347 std::vector<int> correctedSupportOffset;
350 //std::vector<TCoefficientType> coef(mSupportSize);
351 // DD(EvaluateIndex);
352 //std::vector<int> CorrectedSupportOffset;//(mSupportSize);
353 std::vector<IndexType> correctedSupportIndex;//(mSupportSize);
354 correctedSupportIndex.resize(mSupportSize);
355 correctedSupportOffset.resize(mSupportSize);
356 for(unsigned int i=0; i<mSupportSize; i++) {
357 // DD(mSupportIndex[i]);
358 for (unsigned int l=0; l<d; l++) {
359 long a = mSupportIndex[i][l] + evaluateIndex[l];
360 long b = this->m_Coefficients->GetLargestPossibleRegion().GetSize(l);
365 correctedSupportIndex[i][l] = -a - d2*(-a/d2) - evaluateIndex[l];//mSupportIndex[i][l]-a;
368 correctedSupportIndex[i][l] = d2 - a - evaluateIndex[l];
370 correctedSupportIndex[i][l] = mSupportIndex[i][l]; //a - d2*(a/d2) - EvaluateIndex[l];
374 correctedSupportIndex[i][l] = d2 - a - EvaluateIndex[l];//mSupportIndex[i][l] - (a-(b-1));
377 correctedSupportIndex[i][l] = mSupportIndex[i][l];
382 // DD(correctedSupportIndex[i]);
383 correctedSupportOffset[i] = itk::Index2Offset<TImageType::ImageDimension>(correctedSupportIndex[i], mInputMemoryOffset);
385 // for (unsigned int l=0; l<d; l++) {
386 // EvaluateIndex[l] = EvaluateIndex[l] + correctedSupportIndex[0][l];
388 // DD(EvaluateIndex);
389 psupport = &correctedSupportOffset[0];
391 psupport = &mSupportOffset[0];
394 // Get pointer to first coefficient. Even if in some boundary cases,
395 // EvaluateIndex is out of the coefficient image,
396 const TCoefficientType * pcoef = &(this->m_Coefficients->GetPixel(evaluateIndex));
398 // Main loop over BSpline support
399 TCoefficientType interpolated = 0.0;
400 for (unsigned int p=0; p<mSupportSize; p++) {
401 interpolated += pcoef[*psupport] * (*pweights);
406 // Return interpolated value
407 return(interpolated);
409 //====================================================================
412 #endif //_ITKINTERPOLATEIMAGEFUNCTIONWITHLUT_TXX