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() {
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 VectorBSplineInterpolateImageFunctionWithLUT<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 VectorBSplineInterpolateImageFunctionWithLUT<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 VectorBSplineInterpolateImageFunctionWithLUT<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 VectorBSplineInterpolateImageFunctionWithLUT<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 VectorBSplineInterpolateImageFunctionWithLUT<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 VectorBSplineInterpolateImageFunctionWithLUT<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);
160 if (!inputData) return;
161 UpdateWeightsProperties();
165 //====================================================================
166 template <class TImageType, class TCoordRep, class TCoefficientType>
167 void VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
168 UpdateWeightsProperties() {
170 // Compute Memory offset inside coefficients images (for looping over coefficients)
171 static const unsigned int d = TImageType::ImageDimension;
172 mInputMemoryOffset[0] = 1;
173 for(unsigned int l=1; l<d; l++) {
174 mInputMemoryOffset[l] =
175 mInputMemoryOffset[l-1]*this->m_Coefficients->GetLargestPossibleRegion().GetSize(l-1);
179 if (!mWeightsAreUpToDate)
181 // Compute mSupportOffset according to input size
182 mSupportOffset.resize(mSupportSize);
183 mSupportIndex.resize(mSupportSize);
184 for(unsigned int l=0; l<d; l++) mSupportIndex[0][l] = 0;
185 for(unsigned int k=0; k<mSupportSize; k++) {
187 mSupportOffset[k] = itk::Index2Offset<TImageType::ImageDimension>(mSupportIndex[k], mInputMemoryOffset);
188 // next coefficient index
189 if (k != mSupportSize-1) {
190 for(unsigned int l=0; l<d; l++) mSupportIndex[k+1][l] = mSupportIndex[k][l];
194 mSupportIndex[k+1][l]++;
195 if (static_cast<unsigned int>(mSupportIndex[k+1][l]) == mSupport[l]) {
196 mSupportIndex[k+1][l] = 0;
205 // for(unsigned int l=0; l<d; l++) {
206 // if (mOutputSpacing[l] == -1) {
207 // std::cerr << "Please use SetOutputSpacing before using BLUT interpolator" << std::endl;
212 // Compute BSpline weights if not up to date
213 //if (!mWeightsAreUpToDate)
214 UpdatePrecomputedWeights();
217 // *mIntrinsecErrorMax = *mIntrinsecError = 0.0;
218 // *mNumberOfError = 0;
220 //====================================================================
222 //====================================================================
223 template <class TImageType, class TCoordRep, class TCoefficientType>
224 void VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
225 UpdatePrecomputedWeights() {
226 // mLUTTimer.Reset();
227 // mLUTTimer.Start();
228 mWeightsCalculator.SetSplineOrders(mSplineOrders);
229 mWeightsCalculator.SetSamplingFactors(mSamplingFactors);
230 mWeightsCalculator.ComputeTensorProducts();
231 mWeightsAreUpToDate = true;
233 // coef = new TCoefficientType[mSupportSize];
234 // mCorrectedSupportIndex.resize(mSupportSize);
235 // mCorrectedSupportOffset.resize(mSupportSize);
237 // mLUTTimer.Print("LUT \t");
239 //====================================================================
241 //====================================================================
242 template <class TImageType, class TCoordRep, class TCoefficientType>
243 typename VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::IndexType
244 VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
245 GetSampleIndexOfPixelPosition(const ContinuousIndexType & x, IndexType & EvaluateIndex) const {
248 WARNING : sometimes, a floating number x could not really be
249 represented in memory. In this case, the difference between x and
250 floor(x) can be almost 1 (instead of 0). So we take into account
251 such special case, otherwise it could lead to annoying results.
253 // static const TCoefficientType tiny = 1.0e-7;
256 for(int l=0; l<TImageType::ImageDimension; l++)
258 // bool mChange = false;
260 // Compute t1 = distance to floor
261 TCoefficientType t1 = x[l]- vcl_floor(x[l]);
263 // Compute index in precomputed weights table
264 TCoefficientType t2 = mSamplingFactors[l]*t1;
265 index[l] = (IndexValueType)lrint(t2);
267 // For even order : test if too close to 0.5 (but lower). In this
268 // case : take the next coefficient
269 if (!(mSplineOrders[l] & 1)) {
271 if (mSamplingFactors[l] & 1) {
272 if (index[l] == (int) mSamplingFactors[l]/2+1) EvaluateIndex[l] = EvaluateIndex[l]+1;
275 else if (index[l] == (int) mSamplingFactors[l]/2) EvaluateIndex[l] = EvaluateIndex[l]+1;
279 // Statistics (to be removed)
281 *mIntrinsecError += fabs(index[l]-t2);
283 if (fabs(index[l]-t2)> *mIntrinsecErrorMax) *mIntrinsecErrorMax = fabs(index[l]-t2);
286 // When to close to 1, take the next coefficient for odd order, but
287 // only change index for odd
288 if (index[l] == (int)mSamplingFactors[l]) {
290 if (mSplineOrders[l] & 1) EvaluateIndex[l] = EvaluateIndex[l]+1;
299 //====================================================================
301 //====================================================================
302 template <class TImageType, class TCoordRep, class TCoefficientType>
303 typename VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::OutputType
304 VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
305 EvaluateAtContinuousIndex(const ContinuousIndexType & x) const {
307 // For shorter coding
308 static const unsigned int d = TImageType::ImageDimension;
310 // Compute the index of the first interpolation coefficient in the coefficient image
311 IndexType evaluateIndex;
313 for (unsigned int l=0; l<d; l++) {
314 if (mSplineOrders[l] & 1) { // Use this index calculation for odd splineOrder (like cubic)
315 indx = (long)vcl_floor(x[l]) - mSplineOrders[l] / 2 ; //this->m_SplineOrder / 2;
316 evaluateIndex[l] = indx;
318 else { // Use this index calculation for even splineOrder
319 if (mSplineOrders[l] == 0) evaluateIndex[l] = (long)rint(x[l]);
321 indx = (long)vcl_floor((x[l]+ 0.5)) - mSplineOrders[l] / 2; //this->m_SplineOrder / 2;
322 evaluateIndex[l] = indx;
327 // Compute index of precomputed weights and get pointer to first weights
328 const IndexType weightIndex = GetSampleIndexOfPixelPosition(x, evaluateIndex);
329 const TCoefficientType * pweights = mWeightsCalculator.GetFirstTensorProduct(weightIndex);
332 bool boundaryCase = false;
333 for (unsigned int l=0; l<d; l++) {
334 if ((evaluateIndex[l] < 0) ||
335 (evaluateIndex[l]+mSupport[l]) >= this->m_Coefficients->GetLargestPossibleRegion().GetSize(l)) {
340 // Pointer to support offset
341 const int * psupport;
343 // Special case for boundary (to be changed ....)
344 std::vector<int> correctedSupportOffset;
347 //std::vector<TCoefficientType> coef(mSupportSize);
348 // DD(EvaluateIndex);
349 //std::vector<int> CorrectedSupportOffset;//(mSupportSize);
350 std::vector<IndexType> correctedSupportIndex;//(mSupportSize);
351 correctedSupportIndex.resize(mSupportSize);
352 correctedSupportOffset.resize(mSupportSize);
353 for(unsigned int i=0; i<mSupportSize; i++) {
354 // DD(mSupportIndex[i]);
355 for (unsigned int l=0; l<d; l++) {
356 long a = mSupportIndex[i][l] + evaluateIndex[l];
357 long b = this->m_Coefficients->GetLargestPossibleRegion().GetSize(l);
362 correctedSupportIndex[i][l] = -a - d2*(-a/d2) - evaluateIndex[l];//mSupportIndex[i][l]-a;
366 correctedSupportIndex[i][l] = d2 - a - evaluateIndex[l];
369 correctedSupportIndex[i][l] = mSupportIndex[i][l]; //a - d2*(a/d2) - EvaluateIndex[l];
373 correctedSupportIndex[i][l] = d2 - a - EvaluateIndex[l];//mSupportIndex[i][l] - (a-(b-1));
376 correctedSupportIndex[i][l] = mSupportIndex[i][l];
381 // DD(correctedSupportIndex[i]);
382 correctedSupportOffset[i] = itk::Index2Offset<TImageType::ImageDimension>(correctedSupportIndex[i], mInputMemoryOffset);
384 // for (unsigned int l=0; l<d; l++) {
385 // EvaluateIndex[l] = EvaluateIndex[l] + correctedSupportIndex[0][l];
387 // DD(EvaluateIndex);
388 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,
397 const CoefficientImagePixelType * pcoef = &(this->m_Coefficients->GetPixel(evaluateIndex));
399 // Main loop over BSpline support
400 CoefficientImagePixelType interpolated = 0.0;
401 for (unsigned int p=0; p<mSupportSize; p++) {
402 interpolated += pcoef[*psupport] * (*pweights);
407 // Return interpolated value
408 return(interpolated);
410 //====================================================================
413 //====================================================================
414 template <class TImageType, class TCoordRep, class TCoefficientType>
416 VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
417 EvaluateWeightsAtContinuousIndex(const ContinuousIndexType& x, const TCoefficientType ** pweights, IndexType & evaluateIndex)const {
419 // JV Compute BSpline weights if not up to date! Problem const: pass image as last
420 // if (!mWeightsAreUpToDate) UpdatePrecomputedWeights();
422 // For shorter coding
423 static const unsigned int d = TImageType::ImageDimension;
425 // Compute the index of the first interpolation coefficient in the coefficient image
426 //IndexType evaluateIndex;
428 for (unsigned int l=0; l<d; l++) {
429 if (mSplineOrders[l] & 1) { // Use this index calculation for odd splineOrder (like cubic)
430 indx = (long)vcl_floor(x[l]) - mSplineOrders[l] / 2 ; //this->m_SplineOrder / 2;
431 evaluateIndex[l] = indx;
433 else { // Use this index calculation for even splineOrder
434 if (mSplineOrders[l] == 0) evaluateIndex[l] = (long)rint(x[l]);
436 indx = (long)vcl_floor((x[l]+ 0.5)) - mSplineOrders[l] / 2; //this->m_SplineOrder / 2;
437 evaluateIndex[l] = indx;
442 // Compute index of precomputed weights and get pointer to first weights
443 const IndexType weightIndex = GetSampleIndexOfPixelPosition(x, evaluateIndex);
444 *pweights = mWeightsCalculator.GetFirstTensorProduct(weightIndex);
447 //====================================================================
454 #endif //_CLITKVECTORBSPLINEINTERPOLATEIMAGEFUNCTIONWITHLUT_TXX