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() {
44 //for(int i=0; i<TImageType::ImageDimension; i++) mOutputSpacing[i] = -1;
45 SetLUTSamplingFactor(20);
47 mWeightsAreUpToDate = false;
48 // Following need to be pointer beacause values are updated into
50 // mIntrinsecError = new double;
51 // mNumberOfError = new long;
52 // mIntrinsecErrorMax = new double;
53 // mInputIsCoef = false;
56 //====================================================================
58 //====================================================================
59 template <class TImageType, class TCoordRep, class TCoefficientType>
60 void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
61 SetLUTSamplingFactor(const int& s) {
62 for(int i=0; i<TImageType::ImageDimension; i++) mSamplingFactors[i] = s;
63 mWeightsAreUpToDate = false;
66 //====================================================================
68 //====================================================================
69 template <class TImageType, class TCoordRep, class TCoefficientType>
70 void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
71 SetLUTSamplingFactors(const SizeType& s) {
72 for(int i=0; i<TImageType::ImageDimension; i++) mSamplingFactors[i] = s[i];
73 mWeightsAreUpToDate = false;
76 //====================================================================
78 // //====================================================================
79 // template <class TImageType, class TCoordRep, class TCoefficientType>
80 // void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
81 // SetOutputSpacing(const SpacingType & s) {
82 // for(int i=0; i<TImageType::ImageDimension; i++)
83 // mOutputSpacing[i] = s[i];
84 // // mWeightsAreUpToDate = false;
86 //====================================================================
88 //====================================================================
89 template <class TImageType, class TCoordRep, class TCoefficientType>
90 void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
91 SetSplineOrder(const unsigned int& SplineOrder) {
92 Superclass::SetSplineOrder(SplineOrder);
93 // Compute support and half size
94 static const int d = TImageType::ImageDimension;
95 for(int l=0; l<d; l++) {
96 mSplineOrders[l]= SplineOrder;
97 mSupport[l] = SplineOrder+1;
98 if (mSupport[l] % 2 == 0) { // support is even
99 mHalfSupport[l] = mSupport[l]/2-1;
101 else mHalfSupport[l] = mSupport[l]/2; // support is odd (like cubic spline)
104 for(int l=0; l<d; l++) {
105 mSupportSize *= mSupport[l];
107 mWeightsAreUpToDate = false;
109 //====================================================================
112 //====================================================================
113 template <class TImageType, class TCoordRep, class TCoefficientType>
114 void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
115 SetSplineOrders(const SizeType& SplineOrders) {
116 mSplineOrders=SplineOrders;
118 // Compute support and half size
119 static const int d = TImageType::ImageDimension;
120 for(int l=0; l<d; l++) {
121 mSupport[l] = mSplineOrders[l]+1;
122 if (mSupport[l] % 2 == 0) { // support is even
123 mHalfSupport[l] = mSupport[l]/2-1;
125 else mHalfSupport[l] = mSupport[l]/2; // support is odd (like cubic spline)
128 for(int l=0; l<d; l++) {
129 mSupportSize *= mSupport[l];
131 mWeightsAreUpToDate = false;
133 //====================================================================
135 //====================================================================
136 template <class TImageType, class TCoordRep, class TCoefficientType>
137 void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
138 SetInputImage(const TImageType * inputData) {
140 //==============================
141 // if (!mInputIsCoef)
143 Superclass::SetInputImage(inputData);
146 //==============================
147 // //JV Don't call superclass (decomposition filter is executeed each time!)
150 // this->m_Coefficients = inputData;
151 // if ( this->m_Coefficients.IsNotNull())
153 // this->m_DataLength = this->m_Coefficients->GetBufferedRegion().GetSize();
156 // //Call super-superclass in case more input arrives
157 // itk::ImageFunction<TImageType, ITK_TYPENAME itk::NumericTraits<typename TImageType::PixelType>::RealType, TCoefficientType>::SetInputImage(inputData);
159 if (!inputData) return;
160 UpdateWeightsProperties();
164 //====================================================================
165 template <class TImageType, class TCoordRep, class TCoefficientType>
166 void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
167 UpdateWeightsProperties() {
169 // Compute Memory offset inside coefficients images (for looping over coefficients)
170 static const unsigned int d = TImageType::ImageDimension;
171 mInputMemoryOffset[0] = 1;
172 for(unsigned int l=1; l<d; l++) {
173 mInputMemoryOffset[l] =
174 mInputMemoryOffset[l-1]*this->m_Coefficients->GetLargestPossibleRegion().GetSize(l-1);
178 if (!mWeightsAreUpToDate)
180 // Compute mSupportOffset according to input size
181 mSupportOffset.resize(mSupportSize);
182 mSupportIndex.resize(mSupportSize);
183 for(unsigned int l=0; l<d; l++) mSupportIndex[0][l] = 0;
184 for(unsigned int k=0; k<mSupportSize; k++) {
186 mSupportOffset[k] = Index2Offset<TImageType::ImageDimension>(mSupportIndex[k], mInputMemoryOffset);
187 // next coefficient index
188 if (k != mSupportSize-1) {
189 for(unsigned int l=0; l<d; l++) mSupportIndex[k+1][l] = mSupportIndex[k][l];
193 mSupportIndex[k+1][l]++;
194 if (static_cast<unsigned int>(mSupportIndex[k+1][l]) == mSupport[l]) {
195 mSupportIndex[k+1][l] = 0;
204 // for(unsigned int l=0; l<d; l++) {
205 // if (mOutputSpacing[l] == -1) {
206 // std::cerr << "Please use SetOutputSpacing before using BLUT interpolator" << std::endl;
211 // Compute BSpline weights if not up to date
212 //if (!mWeightsAreUpToDate)
213 UpdatePrecomputedWeights();
216 // *mIntrinsecErrorMax = *mIntrinsecError = 0.0;
217 // *mNumberOfError = 0;
219 //====================================================================
221 //====================================================================
222 template <class TImageType, class TCoordRep, class TCoefficientType>
223 void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
224 UpdatePrecomputedWeights() {
225 // mLUTTimer.Reset();
226 // mLUTTimer.Start();
227 mWeightsCalculator.SetSplineOrders(mSplineOrders);
228 mWeightsCalculator.SetSamplingFactors(mSamplingFactors);
229 mWeightsCalculator.ComputeTensorProducts();
230 mWeightsAreUpToDate = true;
232 // coef = new TCoefficientType[mSupportSize];
233 // mCorrectedSupportIndex.resize(mSupportSize);
234 // mCorrectedSupportOffset.resize(mSupportSize);
236 // mLUTTimer.Print("LUT \t");
238 //====================================================================
240 //====================================================================
241 template <class TImageType, class TCoordRep, class TCoefficientType>
242 typename BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::IndexType
243 BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
244 GetSampleIndexOfPixelPosition(const ContinuousIndexType & x, IndexType & EvaluateIndex) const {
247 WARNING : sometimes, a floating number x could not really be
248 represented in memory. In this case, the difference between x and
249 floor(x) can be almost 1 (instead of 0). So we take into account
250 such special case, otherwise it could lead to annoying results.
252 // static const TCoefficientType tiny = 1.0e-7;
255 for(int l=0; l<TImageType::ImageDimension; l++)
257 // bool mChange = false;
259 // Compute t1 = distance to floor
260 TCoefficientType t1 = x[l]- vcl_floor(x[l]);
262 // Compute index in precomputed weights table
263 TCoefficientType t2 = mSamplingFactors[l]*t1;
264 index[l] = (IndexValueType)lrint(t2);
266 // For even order : test if too close to 0.5 (but lower). In this
267 // case : take the next coefficient
268 if (!(mSplineOrders[l] & 1)) {
270 if (mSamplingFactors[l] & 1) {
271 if (index[l] == (int) mSamplingFactors[l]/2+1) EvaluateIndex[l] = EvaluateIndex[l]+1;
274 else if (index[l] == (int) mSamplingFactors[l]/2) EvaluateIndex[l] = EvaluateIndex[l]+1;
278 // Statistics (to be removed)
280 *mIntrinsecError += fabs(index[l]-t2);
282 if (fabs(index[l]-t2)> *mIntrinsecErrorMax) *mIntrinsecErrorMax = fabs(index[l]-t2);
285 // When to close to 1, take the next coefficient for odd order, but
286 // only change index for odd
287 if (index[l] == (int)mSamplingFactors[l]) {
289 if (mSplineOrders[l] & 1) EvaluateIndex[l] = EvaluateIndex[l]+1;
298 //====================================================================
300 //====================================================================
301 template <class TImageType, class TCoordRep, class TCoefficientType>
302 typename BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::OutputType
303 BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
304 EvaluateAtContinuousIndex(const ContinuousIndexType & x) const {
306 // For shorter coding
307 static const unsigned int d = TImageType::ImageDimension;
309 // Compute the index of the first interpolation coefficient in the coefficient image
310 IndexType evaluateIndex;
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;
317 else { // Use this index calculation for even splineOrder
318 if (mSplineOrders[l] == 0) evaluateIndex[l] = (long)rint(x[l]);
320 indx = (long)vcl_floor((x[l]+ 0.5)) - mSplineOrders[l] / 2; //this->m_SplineOrder / 2;
321 evaluateIndex[l] = indx;
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);
331 bool boundaryCase = 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)) {
339 // Pointer to support offset
340 const int * psupport;
342 // Special case for boundary (to be changed ....)
343 std::vector<int> correctedSupportOffset;
346 //std::vector<TCoefficientType> coef(mSupportSize);
347 // DD(EvaluateIndex);
348 //std::vector<int> CorrectedSupportOffset;//(mSupportSize);
349 std::vector<IndexType> correctedSupportIndex;//(mSupportSize);
350 correctedSupportIndex.resize(mSupportSize);
351 correctedSupportOffset.resize(mSupportSize);
352 for(unsigned int i=0; i<mSupportSize; i++) {
353 // DD(mSupportIndex[i]);
354 for (unsigned int l=0; l<d; l++) {
355 long a = mSupportIndex[i][l] + evaluateIndex[l];
356 long b = this->m_Coefficients->GetLargestPossibleRegion().GetSize(l);
361 correctedSupportIndex[i][l] = -a - d2*(-a/d2) - evaluateIndex[l];//mSupportIndex[i][l]-a;
365 correctedSupportIndex[i][l] = d2 - a - evaluateIndex[l];
368 correctedSupportIndex[i][l] = mSupportIndex[i][l]; //a - d2*(a/d2) - EvaluateIndex[l];
372 correctedSupportIndex[i][l] = d2 - a - EvaluateIndex[l];//mSupportIndex[i][l] - (a-(b-1));
375 correctedSupportIndex[i][l] = mSupportIndex[i][l];
380 // DD(correctedSupportIndex[i]);
381 correctedSupportOffset[i] = itk::Index2Offset<TImageType::ImageDimension>(correctedSupportIndex[i], mInputMemoryOffset);
383 // for (unsigned int l=0; l<d; l++) {
384 // EvaluateIndex[l] = EvaluateIndex[l] + correctedSupportIndex[0][l];
386 // DD(EvaluateIndex);
387 psupport = &correctedSupportOffset[0];
390 psupport = &mSupportOffset[0];
393 // Get pointer to first coefficient. Even if in some boundary cases,
394 // EvaluateIndex is out of the coefficient image,
395 const TCoefficientType * pcoef = &(this->m_Coefficients->GetPixel(evaluateIndex));
397 // Main loop over BSpline support
398 TCoefficientType interpolated = 0.0;
399 for (unsigned int p=0; p<mSupportSize; p++) {
400 interpolated += pcoef[*psupport] * (*pweights);
405 // Return interpolated value
406 return(interpolated);
408 //====================================================================
411 #endif //_ITKINTERPOLATEIMAGEFUNCTIONWITHLUT_TXX