1 #ifndef _CLITKVECTORBSPLINEINTERPOLATEIMAGEFUNCTIONWITHLUT_TXX
2 #define _CLITKVECTORBSPLINEINTERPOLATEIMAGEFUNCTIONWITHLUT_TXX
4 /* =========================================================================
6 @file itkBSplineInterpolateImageFunctionWithLUT.txx
7 @author jefvdmb@gmail.com
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 VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
26 VectorBSplineInterpolateImageFunctionWithLUT():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 VectorBSplineInterpolateImageFunctionWithLUT<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 VectorBSplineInterpolateImageFunctionWithLUT<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 VectorBSplineInterpolateImageFunctionWithLUT<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 VectorBSplineInterpolateImageFunctionWithLUT<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 VectorBSplineInterpolateImageFunctionWithLUT<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 VectorBSplineInterpolateImageFunctionWithLUT<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);
144 if (!inputData) return;
145 UpdateWeightsProperties();
149 //====================================================================
150 template <class TImageType, class TCoordRep, class TCoefficientType>
151 void VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
152 UpdateWeightsProperties() {
154 // Compute Memory offset inside coefficients images (for looping over coefficients)
155 static const unsigned int d = TImageType::ImageDimension;
156 mInputMemoryOffset[0] = 1;
157 for(unsigned int l=1; l<d; l++) {
158 mInputMemoryOffset[l] =
159 mInputMemoryOffset[l-1]*this->m_Coefficients->GetLargestPossibleRegion().GetSize(l-1);
163 if (!mWeightsAreUpToDate)
165 // Compute mSupportOffset according to input size
166 mSupportOffset.resize(mSupportSize);
167 mSupportIndex.resize(mSupportSize);
168 for(unsigned int l=0; l<d; l++) mSupportIndex[0][l] = 0;
169 for(unsigned int k=0; k<mSupportSize; k++) {
171 mSupportOffset[k] = itk::Index2Offset<TImageType::ImageDimension>(mSupportIndex[k], mInputMemoryOffset);
172 // next coefficient index
173 if (k != mSupportSize-1) {
174 for(unsigned int l=0; l<d; l++) mSupportIndex[k+1][l] = mSupportIndex[k][l];
178 mSupportIndex[k+1][l]++;
179 if (static_cast<unsigned int>(mSupportIndex[k+1][l]) == mSupport[l]) {
180 mSupportIndex[k+1][l] = 0;
189 // for(unsigned int l=0; l<d; l++) {
190 // if (mOutputSpacing[l] == -1) {
191 // std::cerr << "Please use SetOutputSpacing before using BLUT interpolator" << std::endl;
196 // Compute BSpline weights if not up to date
197 //if (!mWeightsAreUpToDate)
198 UpdatePrecomputedWeights();
201 // *mIntrinsecErrorMax = *mIntrinsecError = 0.0;
202 // *mNumberOfError = 0;
204 //====================================================================
206 //====================================================================
207 template <class TImageType, class TCoordRep, class TCoefficientType>
208 void VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
209 UpdatePrecomputedWeights() {
210 // mLUTTimer.Reset();
211 // mLUTTimer.Start();
212 mWeightsCalculator.SetSplineOrders(mSplineOrders);
213 mWeightsCalculator.SetSamplingFactors(mSamplingFactors);
214 mWeightsCalculator.ComputeTensorProducts();
215 mWeightsAreUpToDate = true;
217 // coef = new TCoefficientType[mSupportSize];
218 // mCorrectedSupportIndex.resize(mSupportSize);
219 // mCorrectedSupportOffset.resize(mSupportSize);
221 // mLUTTimer.Print("LUT \t");
223 //====================================================================
225 //====================================================================
226 template <class TImageType, class TCoordRep, class TCoefficientType>
227 typename VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::IndexType
228 VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
229 GetSampleIndexOfPixelPosition(const ContinuousIndexType & x, IndexType & EvaluateIndex) const {
232 WARNING : sometimes, a floating number x could not really be
233 represented in memory. In this case, the difference between x and
234 floor(x) can be almost 1 (instead of 0). So we take into account
235 such special case, otherwise it could lead to annoying results.
237 // static const TCoefficientType tiny = 1.0e-7;
240 for(int l=0; l<TImageType::ImageDimension; l++)
242 // bool mChange = false;
244 // Compute t1 = distance to floor
245 TCoefficientType t1 = x[l]- vcl_floor(x[l]);
247 // Compute index in precomputed weights table
248 TCoefficientType t2 = mSamplingFactors[l]*t1;
249 index[l] = (IndexValueType)lrint(t2);
251 // For even order : test if too close to 0.5 (but lower). In this
252 // case : take the next coefficient
253 if (!(mSplineOrders[l] & 1)) {
255 if (mSamplingFactors[l] & 1) {
256 if (index[l] == (int) mSamplingFactors[l]/2+1) EvaluateIndex[l] = EvaluateIndex[l]+1;
259 else if (index[l] == (int) mSamplingFactors[l]/2) EvaluateIndex[l] = EvaluateIndex[l]+1;
263 // Statistics (to be removed)
265 *mIntrinsecError += fabs(index[l]-t2);
267 if (fabs(index[l]-t2)> *mIntrinsecErrorMax) *mIntrinsecErrorMax = fabs(index[l]-t2);
270 // When to close to 1, take the next coefficient for odd order, but
271 // only change index for odd
272 if (index[l] == (int)mSamplingFactors[l]) {
274 if (mSplineOrders[l] & 1) EvaluateIndex[l] = EvaluateIndex[l]+1;
283 //====================================================================
285 //====================================================================
286 template <class TImageType, class TCoordRep, class TCoefficientType>
287 typename VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::OutputType
288 VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
289 EvaluateAtContinuousIndex(const ContinuousIndexType & x) const {
291 // For shorter coding
292 static const unsigned int d = TImageType::ImageDimension;
294 // Compute the index of the first interpolation coefficient in the coefficient image
295 IndexType evaluateIndex;
297 for (unsigned int l=0; l<d; l++) {
298 if (mSplineOrders[l] & 1) { // Use this index calculation for odd splineOrder (like cubic)
299 indx = (long)vcl_floor(x[l]) - mSplineOrders[l] / 2 ; //this->m_SplineOrder / 2;
300 evaluateIndex[l] = indx;
302 else { // Use this index calculation for even splineOrder
303 if (mSplineOrders[l] == 0) evaluateIndex[l] = (long)rint(x[l]);
305 indx = (long)vcl_floor((x[l]+ 0.5)) - mSplineOrders[l] / 2; //this->m_SplineOrder / 2;
306 evaluateIndex[l] = indx;
311 // Compute index of precomputed weights and get pointer to first weights
312 const IndexType weightIndex = GetSampleIndexOfPixelPosition(x, evaluateIndex);
313 const TCoefficientType * pweights = mWeightsCalculator.GetFirstTensorProduct(weightIndex);
316 bool boundaryCase = false;
317 for (unsigned int l=0; l<d; l++) {
318 if ((evaluateIndex[l] < 0) ||
319 (evaluateIndex[l]+mSupport[l]) >= this->m_Coefficients->GetLargestPossibleRegion().GetSize(l)) {
324 // Pointer to support offset
325 const int * psupport;
327 // Special case for boundary (to be changed ....)
328 std::vector<int> correctedSupportOffset;
331 //std::vector<TCoefficientType> coef(mSupportSize);
332 // DD(EvaluateIndex);
333 //std::vector<int> CorrectedSupportOffset;//(mSupportSize);
334 std::vector<IndexType> correctedSupportIndex;//(mSupportSize);
335 correctedSupportIndex.resize(mSupportSize);
336 correctedSupportOffset.resize(mSupportSize);
337 for(unsigned int i=0; i<mSupportSize; i++) {
338 // DD(mSupportIndex[i]);
339 for (unsigned int l=0; l<d; l++) {
340 long a = mSupportIndex[i][l] + evaluateIndex[l];
341 long b = this->m_Coefficients->GetLargestPossibleRegion().GetSize(l);
346 correctedSupportIndex[i][l] = -a - d2*(-a/d2) - evaluateIndex[l];//mSupportIndex[i][l]-a;
350 correctedSupportIndex[i][l] = d2 - a - evaluateIndex[l];
353 correctedSupportIndex[i][l] = mSupportIndex[i][l]; //a - d2*(a/d2) - EvaluateIndex[l];
357 correctedSupportIndex[i][l] = d2 - a - EvaluateIndex[l];//mSupportIndex[i][l] - (a-(b-1));
360 correctedSupportIndex[i][l] = mSupportIndex[i][l];
365 // DD(correctedSupportIndex[i]);
366 correctedSupportOffset[i] = itk::Index2Offset<TImageType::ImageDimension>(correctedSupportIndex[i], mInputMemoryOffset);
368 // for (unsigned int l=0; l<d; l++) {
369 // EvaluateIndex[l] = EvaluateIndex[l] + correctedSupportIndex[0][l];
371 // DD(EvaluateIndex);
372 psupport = &correctedSupportOffset[0];
375 psupport = &mSupportOffset[0];
378 // Get pointer to first coefficient. Even if in some boundary cases,
379 // EvaluateIndex is out of the coefficient image,
381 const CoefficientImagePixelType * pcoef = &(this->m_Coefficients->GetPixel(evaluateIndex));
383 // Main loop over BSpline support
384 CoefficientImagePixelType interpolated = 0.0;
385 for (unsigned int p=0; p<mSupportSize; p++) {
386 interpolated += pcoef[*psupport] * (*pweights);
391 // Return interpolated value
392 return(interpolated);
394 //====================================================================
397 //====================================================================
398 template <class TImageType, class TCoordRep, class TCoefficientType>
400 VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
401 EvaluateWeightsAtContinuousIndex(const ContinuousIndexType& x, const TCoefficientType ** pweights, IndexType & evaluateIndex)const {
403 // JV Compute BSpline weights if not up to date! Problem const: pass image as last
404 // if (!mWeightsAreUpToDate) UpdatePrecomputedWeights();
406 // For shorter coding
407 static const unsigned int d = TImageType::ImageDimension;
409 // Compute the index of the first interpolation coefficient in the coefficient image
410 //IndexType evaluateIndex;
412 for (unsigned int l=0; l<d; l++) {
413 if (mSplineOrders[l] & 1) { // Use this index calculation for odd splineOrder (like cubic)
414 indx = (long)vcl_floor(x[l]) - mSplineOrders[l] / 2 ; //this->m_SplineOrder / 2;
415 evaluateIndex[l] = indx;
417 else { // Use this index calculation for even splineOrder
418 if (mSplineOrders[l] == 0) evaluateIndex[l] = (long)rint(x[l]);
420 indx = (long)vcl_floor((x[l]+ 0.5)) - mSplineOrders[l] / 2; //this->m_SplineOrder / 2;
421 evaluateIndex[l] = indx;
426 // Compute index of precomputed weights and get pointer to first weights
427 const IndexType weightIndex = GetSampleIndexOfPixelPosition(x, evaluateIndex);
428 *pweights = mWeightsCalculator.GetFirstTensorProduct(weightIndex);
431 //====================================================================
438 #endif //_CLITKVECTORBSPLINEINTERPOLATEIMAGEFUNCTIONWITHLUT_TXX