1 #ifndef _ITKINTERPOLATEIMAGEFUNCTIONWITHLUT_TXX
2 #define _ITKINTERPOLATEIMAGEFUNCTIONWITHLUT_TXX
4 /* =========================================================================
6 @file itkBSplineInterpolateImageFunctionWithLUT.txx
7 @author David Sarrut <david.sarrut@creatis.insa-lyon.fr>
10 * CREATIS (Centre de Recherche et d'Applications en Traitement de l'Image).
11 All rights reserved. See Doc/License.txt or
12 http://www.creatis.insa-lyon.fr/Public/Gdcm/License.html for details.
13 * Léon Bérard cancer center, 28 rue Laënnec, 69373 Lyon cedex 08, France
14 * http://www.creatis.insa-lyon.fr/rio
16 This software is distributed WITHOUT ANY WARRANTY; without even the
17 implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
18 PURPOSE. See the above copyright notices for more information.
20 ========================================================================= */
23 //====================================================================
24 template <class TImageType, class TCoordRep, class TCoefficientType>
25 BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
26 BSplineInterpolateImageFunctionWithLUT():Superclass() {
28 //for(int i=0; i<TImageType::ImageDimension; i++) mOutputSpacing[i] = -1;
29 SetLUTSamplingFactor(20);
31 mWeightsAreUpToDate = false;
32 // Following need to be pointer beacause values are updated into
34 // mIntrinsecError = new double;
35 // mNumberOfError = new long;
36 // mIntrinsecErrorMax = new double;
37 // mInputIsCoef = false;
40 //====================================================================
42 //====================================================================
43 template <class TImageType, class TCoordRep, class TCoefficientType>
44 void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
45 SetLUTSamplingFactor(const int& s) {
46 for(int i=0; i<TImageType::ImageDimension; i++) mSamplingFactors[i] = s;
47 mWeightsAreUpToDate = false;
50 //====================================================================
52 //====================================================================
53 template <class TImageType, class TCoordRep, class TCoefficientType>
54 void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
55 SetLUTSamplingFactors(const SizeType& s) {
56 for(int i=0; i<TImageType::ImageDimension; i++) mSamplingFactors[i] = s[i];
57 mWeightsAreUpToDate = false;
60 //====================================================================
62 // //====================================================================
63 // template <class TImageType, class TCoordRep, class TCoefficientType>
64 // void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
65 // SetOutputSpacing(const SpacingType & s) {
66 // for(int i=0; i<TImageType::ImageDimension; i++)
67 // mOutputSpacing[i] = s[i];
68 // // mWeightsAreUpToDate = false;
70 //====================================================================
72 //====================================================================
73 template <class TImageType, class TCoordRep, class TCoefficientType>
74 void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
75 SetSplineOrder(const unsigned int& SplineOrder) {
76 Superclass::SetSplineOrder(SplineOrder);
77 // Compute support and half size
78 static const int d = TImageType::ImageDimension;
79 for(int l=0; l<d; l++) {
80 mSplineOrders[l]= SplineOrder;
81 mSupport[l] = SplineOrder+1;
82 if (mSupport[l] % 2 == 0) { // support is even
83 mHalfSupport[l] = mSupport[l]/2-1;
85 else mHalfSupport[l] = mSupport[l]/2; // support is odd (like cubic spline)
88 for(int l=0; l<d; l++) {
89 mSupportSize *= mSupport[l];
91 mWeightsAreUpToDate = false;
93 //====================================================================
96 //====================================================================
97 template <class TImageType, class TCoordRep, class TCoefficientType>
98 void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
99 SetSplineOrders(const SizeType& SplineOrders) {
100 mSplineOrders=SplineOrders;
102 // Compute support and half size
103 static const int d = TImageType::ImageDimension;
104 for(int l=0; l<d; l++) {
105 mSupport[l] = mSplineOrders[l]+1;
106 if (mSupport[l] % 2 == 0) { // support is even
107 mHalfSupport[l] = mSupport[l]/2-1;
109 else mHalfSupport[l] = mSupport[l]/2; // support is odd (like cubic spline)
112 for(int l=0; l<d; l++) {
113 mSupportSize *= mSupport[l];
115 mWeightsAreUpToDate = false;
117 //====================================================================
119 //====================================================================
120 template <class TImageType, class TCoordRep, class TCoefficientType>
121 void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
122 SetInputImage(const TImageType * inputData) {
124 //==============================
125 // if (!mInputIsCoef)
127 Superclass::SetInputImage(inputData);
130 //==============================
131 // //JV Don't call superclass (decomposition filter is executeed each time!)
134 // this->m_Coefficients = inputData;
135 // if ( this->m_Coefficients.IsNotNull())
137 // this->m_DataLength = this->m_Coefficients->GetBufferedRegion().GetSize();
140 // //Call super-superclass in case more input arrives
141 // itk::ImageFunction<TImageType, ITK_TYPENAME itk::NumericTraits<typename TImageType::PixelType>::RealType, TCoefficientType>::SetInputImage(inputData);
143 if (!inputData) return;
144 UpdateWeightsProperties();
148 //====================================================================
149 template <class TImageType, class TCoordRep, class TCoefficientType>
150 void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
151 UpdateWeightsProperties() {
153 // Compute Memory offset inside coefficients images (for looping over coefficients)
154 static const unsigned int d = TImageType::ImageDimension;
155 mInputMemoryOffset[0] = 1;
156 for(unsigned int l=1; l<d; l++) {
157 mInputMemoryOffset[l] =
158 mInputMemoryOffset[l-1]*this->m_Coefficients->GetLargestPossibleRegion().GetSize(l-1);
162 if (!mWeightsAreUpToDate)
164 // Compute mSupportOffset according to input size
165 mSupportOffset.resize(mSupportSize);
166 mSupportIndex.resize(mSupportSize);
167 for(unsigned int l=0; l<d; l++) mSupportIndex[0][l] = 0;
168 for(unsigned int k=0; k<mSupportSize; k++) {
170 mSupportOffset[k] = Index2Offset<TImageType::ImageDimension>(mSupportIndex[k], mInputMemoryOffset);
171 // next coefficient index
172 if (k != mSupportSize-1) {
173 for(unsigned int l=0; l<d; l++) mSupportIndex[k+1][l] = mSupportIndex[k][l];
177 mSupportIndex[k+1][l]++;
178 if (static_cast<unsigned int>(mSupportIndex[k+1][l]) == mSupport[l]) {
179 mSupportIndex[k+1][l] = 0;
188 // for(unsigned int l=0; l<d; l++) {
189 // if (mOutputSpacing[l] == -1) {
190 // std::cerr << "Please use SetOutputSpacing before using BLUT interpolator" << std::endl;
195 // Compute BSpline weights if not up to date
196 //if (!mWeightsAreUpToDate)
197 UpdatePrecomputedWeights();
200 // *mIntrinsecErrorMax = *mIntrinsecError = 0.0;
201 // *mNumberOfError = 0;
203 //====================================================================
205 //====================================================================
206 template <class TImageType, class TCoordRep, class TCoefficientType>
207 void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
208 UpdatePrecomputedWeights() {
209 // mLUTTimer.Reset();
210 // mLUTTimer.Start();
211 mWeightsCalculator.SetSplineOrders(mSplineOrders);
212 mWeightsCalculator.SetSamplingFactors(mSamplingFactors);
213 mWeightsCalculator.ComputeTensorProducts();
214 mWeightsAreUpToDate = true;
216 // coef = new TCoefficientType[mSupportSize];
217 // mCorrectedSupportIndex.resize(mSupportSize);
218 // mCorrectedSupportOffset.resize(mSupportSize);
220 // mLUTTimer.Print("LUT \t");
222 //====================================================================
224 //====================================================================
225 template <class TImageType, class TCoordRep, class TCoefficientType>
226 typename BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::IndexType
227 BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
228 GetSampleIndexOfPixelPosition(const ContinuousIndexType & x, IndexType & EvaluateIndex) const {
231 WARNING : sometimes, a floating number x could not really be
232 represented in memory. In this case, the difference between x and
233 floor(x) can be almost 1 (instead of 0). So we take into account
234 such special case, otherwise it could lead to annoying results.
236 // static const TCoefficientType tiny = 1.0e-7;
239 for(int l=0; l<TImageType::ImageDimension; l++)
241 // bool mChange = false;
243 // Compute t1 = distance to floor
244 TCoefficientType t1 = x[l]- vcl_floor(x[l]);
246 // Compute index in precomputed weights table
247 TCoefficientType t2 = mSamplingFactors[l]*t1;
248 index[l] = (IndexValueType)lrint(t2);
250 // For even order : test if too close to 0.5 (but lower). In this
251 // case : take the next coefficient
252 if (!(mSplineOrders[l] & 1)) {
254 if (mSamplingFactors[l] & 1) {
255 if (index[l] == (int) mSamplingFactors[l]/2+1) EvaluateIndex[l] = EvaluateIndex[l]+1;
258 else if (index[l] == (int) mSamplingFactors[l]/2) EvaluateIndex[l] = EvaluateIndex[l]+1;
262 // Statistics (to be removed)
264 *mIntrinsecError += fabs(index[l]-t2);
266 if (fabs(index[l]-t2)> *mIntrinsecErrorMax) *mIntrinsecErrorMax = fabs(index[l]-t2);
269 // When to close to 1, take the next coefficient for odd order, but
270 // only change index for odd
271 if (index[l] == (int)mSamplingFactors[l]) {
273 if (mSplineOrders[l] & 1) EvaluateIndex[l] = EvaluateIndex[l]+1;
282 //====================================================================
284 //====================================================================
285 template <class TImageType, class TCoordRep, class TCoefficientType>
286 typename BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::OutputType
287 BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
288 EvaluateAtContinuousIndex(const ContinuousIndexType & x) const {
290 // For shorter coding
291 static const unsigned int d = TImageType::ImageDimension;
293 // Compute the index of the first interpolation coefficient in the coefficient image
294 IndexType evaluateIndex;
296 for (unsigned int l=0; l<d; l++) {
297 if (mSplineOrders[l] & 1) { // Use this index calculation for odd splineOrder (like cubic)
298 indx = (long)vcl_floor(x[l]) - mSplineOrders[l] / 2 ; //this->m_SplineOrder / 2;
299 evaluateIndex[l] = indx;
301 else { // Use this index calculation for even splineOrder
302 if (mSplineOrders[l] == 0) evaluateIndex[l] = (long)rint(x[l]);
304 indx = (long)vcl_floor((x[l]+ 0.5)) - mSplineOrders[l] / 2; //this->m_SplineOrder / 2;
305 evaluateIndex[l] = indx;
310 // Compute index of precomputed weights and get pointer to first weights
311 const IndexType weightIndex = GetSampleIndexOfPixelPosition(x, evaluateIndex);
312 const TCoefficientType * pweights = mWeightsCalculator.GetFirstTensorProduct(weightIndex);
315 bool boundaryCase = false;
316 for (unsigned int l=0; l<d; l++) {
317 if ((evaluateIndex[l] < 0) ||
318 (evaluateIndex[l]+mSupport[l]) >= this->m_Coefficients->GetLargestPossibleRegion().GetSize(l)) {
323 // Pointer to support offset
324 const int * psupport;
326 // Special case for boundary (to be changed ....)
327 std::vector<int> correctedSupportOffset;
330 //std::vector<TCoefficientType> coef(mSupportSize);
331 // DD(EvaluateIndex);
332 //std::vector<int> CorrectedSupportOffset;//(mSupportSize);
333 std::vector<IndexType> correctedSupportIndex;//(mSupportSize);
334 correctedSupportIndex.resize(mSupportSize);
335 correctedSupportOffset.resize(mSupportSize);
336 for(unsigned int i=0; i<mSupportSize; i++) {
337 // DD(mSupportIndex[i]);
338 for (unsigned int l=0; l<d; l++) {
339 long a = mSupportIndex[i][l] + evaluateIndex[l];
340 long b = this->m_Coefficients->GetLargestPossibleRegion().GetSize(l);
345 correctedSupportIndex[i][l] = -a - d2*(-a/d2) - evaluateIndex[l];//mSupportIndex[i][l]-a;
349 correctedSupportIndex[i][l] = d2 - a - evaluateIndex[l];
352 correctedSupportIndex[i][l] = mSupportIndex[i][l]; //a - d2*(a/d2) - EvaluateIndex[l];
356 correctedSupportIndex[i][l] = d2 - a - EvaluateIndex[l];//mSupportIndex[i][l] - (a-(b-1));
359 correctedSupportIndex[i][l] = mSupportIndex[i][l];
364 // DD(correctedSupportIndex[i]);
365 correctedSupportOffset[i] = itk::Index2Offset<TImageType::ImageDimension>(correctedSupportIndex[i], mInputMemoryOffset);
367 // for (unsigned int l=0; l<d; l++) {
368 // EvaluateIndex[l] = EvaluateIndex[l] + correctedSupportIndex[0][l];
370 // DD(EvaluateIndex);
371 psupport = &correctedSupportOffset[0];
374 psupport = &mSupportOffset[0];
377 // Get pointer to first coefficient. Even if in some boundary cases,
378 // EvaluateIndex is out of the coefficient image,
379 const TCoefficientType * pcoef = &(this->m_Coefficients->GetPixel(evaluateIndex));
381 // Main loop over BSpline support
382 TCoefficientType interpolated = 0.0;
383 for (unsigned int p=0; p<mSupportSize; p++) {
384 interpolated += pcoef[*psupport] * (*pweights);
389 // Return interpolated value
390 return(interpolated);
392 //====================================================================
395 #endif //_ITKINTERPOLATEIMAGEFUNCTIONWITHLUT_TXX