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 _CLITKVECTORBSPLINEINTERPOLATEIMAGEFUNCTIONWITHLUT_TXX
19 #define _CLITKVECTORBSPLINEINTERPOLATEIMAGEFUNCTIONWITHLUT_TXX
20 /* =========================================================================
22 @file itkBSplineInterpolateImageFunctionWithLUT.txx
23 @author jefvdmb@gmail.com
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 VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
42 VectorBSplineInterpolateImageFunctionWithLUT():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 VectorBSplineInterpolateImageFunctionWithLUT<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 VectorBSplineInterpolateImageFunctionWithLUT<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 VectorBSplineInterpolateImageFunctionWithLUT<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 VectorBSplineInterpolateImageFunctionWithLUT<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 VectorBSplineInterpolateImageFunctionWithLUT<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 VectorBSplineInterpolateImageFunctionWithLUT<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);
164 if (!inputData) return;
165 UpdateWeightsProperties();
169 //====================================================================
170 template <class TImageType, class TCoordRep, class TCoefficientType>
171 void VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
172 UpdateWeightsProperties()
175 // Compute Memory offset inside coefficients images (for looping over coefficients)
176 static const unsigned int d = TImageType::ImageDimension;
177 mInputMemoryOffset[0] = 1;
178 for(unsigned int l=1; l<d; l++) {
179 mInputMemoryOffset[l] =
180 mInputMemoryOffset[l-1]*this->m_Coefficients->GetLargestPossibleRegion().GetSize(l-1);
184 if (!mWeightsAreUpToDate) {
185 // Compute mSupportOffset according to input size
186 mSupportOffset.resize(mSupportSize);
187 mSupportIndex.resize(mSupportSize);
188 for(unsigned int l=0; l<d; l++) mSupportIndex[0][l] = 0;
189 for(unsigned int k=0; k<mSupportSize; k++) {
191 mSupportOffset[k] = itk::Index2Offset<TImageType::ImageDimension>(mSupportIndex[k], mInputMemoryOffset);
192 // next coefficient index
193 if (k != mSupportSize-1) {
194 for(unsigned int l=0; l<d; l++) mSupportIndex[k+1][l] = mSupportIndex[k][l];
198 mSupportIndex[k+1][l]++;
199 if (static_cast<unsigned int>(mSupportIndex[k+1][l]) == mSupport[l]) {
200 mSupportIndex[k+1][l] = 0;
208 // for(unsigned int l=0; l<d; l++) {
209 // if (mOutputSpacing[l] == -1) {
210 // std::cerr << "Please use SetOutputSpacing before using BLUT interpolator" << std::endl;
215 // Compute BSpline weights if not up to date
216 //if (!mWeightsAreUpToDate)
217 UpdatePrecomputedWeights();
220 // *mIntrinsecErrorMax = *mIntrinsecError = 0.0;
221 // *mNumberOfError = 0;
223 //====================================================================
225 //====================================================================
226 template <class TImageType, class TCoordRep, class TCoefficientType>
227 void VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
228 UpdatePrecomputedWeights()
230 // mLUTTimer.Reset();
231 // mLUTTimer.Start();
232 mWeightsCalculator.SetSplineOrders(mSplineOrders);
233 mWeightsCalculator.SetSamplingFactors(mSamplingFactors);
234 mWeightsCalculator.ComputeTensorProducts();
235 mWeightsAreUpToDate = true;
237 // coef = new TCoefficientType[mSupportSize];
238 // mCorrectedSupportIndex.resize(mSupportSize);
239 // mCorrectedSupportOffset.resize(mSupportSize);
241 // mLUTTimer.Print("LUT \t");
243 //====================================================================
245 //====================================================================
246 template <class TImageType, class TCoordRep, class TCoefficientType>
247 typename VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::IndexType
248 VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
249 GetSampleIndexOfPixelPosition(const ContinuousIndexType & x, IndexType & EvaluateIndex) const
253 WARNING : sometimes, a floating number x could not really be
254 represented in memory. In this case, the difference between x and
255 floor(x) can be almost 1 (instead of 0). So we take into account
256 such special case, otherwise it could lead to annoying results.
258 // static const TCoefficientType tiny = 1.0e-7;
261 for(int l=0; l<TImageType::ImageDimension; l++) {
262 // bool mChange = false;
264 // Compute t1 = distance to floor
265 TCoefficientType t1 = x[l]- vcl_floor(x[l]);
267 // Compute index in precomputed weights table
268 TCoefficientType t2 = mSamplingFactors[l]*t1;
269 index[l] = (IndexValueType)lrint(t2);
271 // For even order : test if too close to 0.5 (but lower). In this
272 // case : take the next coefficient
273 if (!(mSplineOrders[l] & 1)) {
275 if (mSamplingFactors[l] & 1) {
276 if (index[l] == (int) mSamplingFactors[l]/2+1) EvaluateIndex[l] = EvaluateIndex[l]+1;
279 else if (index[l] == (int) mSamplingFactors[l]/2) EvaluateIndex[l] = EvaluateIndex[l]+1;
283 // Statistics (to be removed)
285 *mIntrinsecError += fabs(index[l]-t2);
287 if (fabs(index[l]-t2)> *mIntrinsecErrorMax) *mIntrinsecErrorMax = fabs(index[l]-t2);
290 // When to close to 1, take the next coefficient for odd order, but
291 // only change index for odd
292 if (index[l] == (int)mSamplingFactors[l]) {
294 if (mSplineOrders[l] & 1) EvaluateIndex[l] = EvaluateIndex[l]+1;
303 //====================================================================
305 //====================================================================
306 template <class TImageType, class TCoordRep, class TCoefficientType>
307 typename VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::OutputType
308 VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
309 EvaluateAtContinuousIndex(const ContinuousIndexType & x) const
312 // For shorter coding
313 static const unsigned int d = TImageType::ImageDimension;
315 // Compute the index of the first interpolation coefficient in the coefficient image
316 IndexType evaluateIndex;
318 for (unsigned int l=0; l<d; l++) {
319 if (mSplineOrders[l] & 1) { // Use this index calculation for odd splineOrder (like cubic)
320 indx = (long)vcl_floor(x[l]) - mSplineOrders[l] / 2 ; //this->m_SplineOrder / 2;
321 evaluateIndex[l] = indx;
322 } else { // Use this index calculation for even splineOrder
323 if (mSplineOrders[l] == 0) evaluateIndex[l] = (long)rint(x[l]);
325 indx = (long)vcl_floor((x[l]+ 0.5)) - mSplineOrders[l] / 2; //this->m_SplineOrder / 2;
326 evaluateIndex[l] = indx;
331 // Compute index of precomputed weights and get pointer to first weights
332 const IndexType weightIndex = GetSampleIndexOfPixelPosition(x, evaluateIndex);
333 const TCoefficientType * pweights = mWeightsCalculator.GetFirstTensorProduct(weightIndex);
336 bool boundaryCase = false;
337 for (unsigned int l=0; l<d; l++) {
338 if ((evaluateIndex[l] < 0) ||
339 (evaluateIndex[l]+mSupport[l]) >= this->m_Coefficients->GetLargestPossibleRegion().GetSize(l)) {
344 // Pointer to support offset
345 const int * psupport;
347 // Special case for boundary (to be changed ....)
348 std::vector<int> correctedSupportOffset;
351 //std::vector<TCoefficientType> coef(mSupportSize);
352 // DD(EvaluateIndex);
353 //std::vector<int> CorrectedSupportOffset;//(mSupportSize);
354 std::vector<IndexType> correctedSupportIndex;//(mSupportSize);
355 correctedSupportIndex.resize(mSupportSize);
356 correctedSupportOffset.resize(mSupportSize);
357 for(unsigned int i=0; i<mSupportSize; i++) {
358 // DD(mSupportIndex[i]);
359 for (unsigned int l=0; l<d; l++) {
360 long a = mSupportIndex[i][l] + evaluateIndex[l];
361 long b = this->m_Coefficients->GetLargestPossibleRegion().GetSize(l);
366 correctedSupportIndex[i][l] = -a - d2*(-a/d2) - evaluateIndex[l];//mSupportIndex[i][l]-a;
369 correctedSupportIndex[i][l] = d2 - a - evaluateIndex[l];
371 correctedSupportIndex[i][l] = mSupportIndex[i][l]; //a - d2*(a/d2) - EvaluateIndex[l];
375 correctedSupportIndex[i][l] = d2 - a - EvaluateIndex[l];//mSupportIndex[i][l] - (a-(b-1));
378 correctedSupportIndex[i][l] = mSupportIndex[i][l];
383 // DD(correctedSupportIndex[i]);
384 correctedSupportOffset[i] = itk::Index2Offset<TImageType::ImageDimension>(correctedSupportIndex[i], mInputMemoryOffset);
386 // for (unsigned int l=0; l<d; l++) {
387 // EvaluateIndex[l] = EvaluateIndex[l] + correctedSupportIndex[0][l];
389 // DD(EvaluateIndex);
390 psupport = &correctedSupportOffset[0];
392 psupport = &mSupportOffset[0];
395 // Get pointer to first coefficient. Even if in some boundary cases,
396 // EvaluateIndex is out of the coefficient image,
398 const CoefficientImagePixelType * pcoef = &(this->m_Coefficients->GetPixel(evaluateIndex));
400 // Main loop over BSpline support
401 CoefficientImagePixelType interpolated = 0.0;
402 for (unsigned int p=0; p<mSupportSize; p++) {
403 interpolated += pcoef[*psupport] * (*pweights);
408 // Return interpolated value
409 return(interpolated);
411 //====================================================================
414 //====================================================================
415 template <class TImageType, class TCoordRep, class TCoefficientType>
417 VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
418 EvaluateWeightsAtContinuousIndex(const ContinuousIndexType& x, const TCoefficientType ** pweights, IndexType & evaluateIndex)const
421 // JV Compute BSpline weights if not up to date! Problem const: pass image as last
422 // if (!mWeightsAreUpToDate) UpdatePrecomputedWeights();
424 // For shorter coding
425 static const unsigned int d = TImageType::ImageDimension;
427 // Compute the index of the first interpolation coefficient in the coefficient image
428 //IndexType evaluateIndex;
430 for (unsigned int l=0; l<d; l++) {
431 if (mSplineOrders[l] & 1) { // Use this index calculation for odd splineOrder (like cubic)
432 indx = (long)vcl_floor(x[l]) - mSplineOrders[l] / 2 ; //this->m_SplineOrder / 2;
433 evaluateIndex[l] = indx;
434 } else { // Use this index calculation for even splineOrder
435 if (mSplineOrders[l] == 0) evaluateIndex[l] = (long)rint(x[l]);
437 indx = (long)vcl_floor((x[l]+ 0.5)) - mSplineOrders[l] / 2; //this->m_SplineOrder / 2;
438 evaluateIndex[l] = indx;
443 // Compute index of precomputed weights and get pointer to first weights
444 const IndexType weightIndex = GetSampleIndexOfPixelPosition(x, evaluateIndex);
445 *pweights = mWeightsCalculator.GetFirstTensorProduct(weightIndex);
448 //====================================================================
455 #endif //_CLITKVECTORBSPLINEINTERPOLATEIMAGEFUNCTIONWITHLUT_TXX