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
23 //====================================================================
24 template <class TImageType, class TCoordRep, class TCoefficientType>
25 BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
26 BSplineInterpolateImageFunctionWithLUT():Superclass() {
29 SetLUTSamplingFactor(20);
31 mWeightsAreUpToDate = false;
34 //====================================================================
36 //====================================================================
37 template <class TImageType, class TCoordRep, class TCoefficientType>
38 void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
39 SetLUTSamplingFactor(const int& s) {
40 for(int i=0; i<TImageType::ImageDimension; i++) mSamplingFactors[i] = s;
41 mWeightsAreUpToDate = false;
44 //====================================================================
46 //====================================================================
47 template <class TImageType, class TCoordRep, class TCoefficientType>
48 void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
49 SetLUTSamplingFactors(const SizeType& s) {
50 for(int i=0; i<TImageType::ImageDimension; i++) mSamplingFactors[i] = s[i];
51 mWeightsAreUpToDate = false;
54 //====================================================================
56 //====================================================================
57 template <class TImageType, class TCoordRep, class TCoefficientType>
58 void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
59 SetSplineOrder(const unsigned int& SplineOrder) {
60 Superclass::SetSplineOrder(SplineOrder);
61 // Compute support and half size
62 static const int d = TImageType::ImageDimension;
63 for(int l=0; l<d; l++) {
64 mSplineOrders[l]= SplineOrder;
65 mSupport[l] = SplineOrder+1;
66 if (mSupport[l] % 2 == 0) { // support is even
67 mHalfSupport[l] = mSupport[l]/2-1;
69 else mHalfSupport[l] = mSupport[l]/2; // support is odd (like cubic spline)
72 for(int l=0; l<d; l++) {
73 mSupportSize *= mSupport[l];
75 mWeightsAreUpToDate = false;
77 //====================================================================
80 //====================================================================
81 template <class TImageType, class TCoordRep, class TCoefficientType>
82 void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
83 SetSplineOrders(const SizeType& SplineOrders) {
84 mSplineOrders=SplineOrders;
86 // Compute support and half size
87 static const int d = TImageType::ImageDimension;
88 for(int l=0; l<d; l++) {
89 mSupport[l] = mSplineOrders[l]+1;
90 if (mSupport[l] % 2 == 0) { // support is even
91 mHalfSupport[l] = mSupport[l]/2-1;
93 else mHalfSupport[l] = mSupport[l]/2; // support is odd (like cubic spline)
96 for(int l=0; l<d; l++) {
97 mSupportSize *= mSupport[l];
99 mWeightsAreUpToDate = false;
101 //====================================================================
103 //====================================================================
104 template <class TImageType, class TCoordRep, class TCoefficientType>
105 void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
106 SetInputImage(const TImageType * inputData) {
108 // JV Call superclass (decomposition filter is executeed each time!)
109 // JV Should call itkBSplineDecompositionFilterWithOBD to allow different order by dimension
110 Superclass::SetInputImage(inputData);
112 // Update the weightproperties
113 if (!inputData) return;
114 UpdateWeightsProperties();
117 //====================================================================
118 template <class TImageType, class TCoordRep, class TCoefficientType>
119 void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
120 UpdateWeightsProperties() {
122 // Compute Memory offset inside coefficients images (for looping over coefficients)
123 static const unsigned int d = TImageType::ImageDimension;
124 mInputMemoryOffset[0] = 1;
125 for(unsigned int l=1; l<d; l++) {
126 mInputMemoryOffset[l] =
127 mInputMemoryOffset[l-1]*this->m_Coefficients->GetLargestPossibleRegion().GetSize(l-1);
130 // Compute mSupportOffset according to input size
131 mSupportOffset.resize(mSupportSize);
132 mSupportIndex.resize(mSupportSize);
133 for(unsigned int l=0; l<d; l++) mSupportIndex[0][l] = 0;
134 for(unsigned int k=0; k<mSupportSize; k++) {
136 mSupportOffset[k] = Index2Offset<TImageType::ImageDimension>(mSupportIndex[k], mInputMemoryOffset);
137 // next coefficient index
138 if (k != mSupportSize-1) {
139 for(unsigned int l=0; l<d; l++) mSupportIndex[k+1][l] = mSupportIndex[k][l];
143 mSupportIndex[k+1][l]++;
144 if (static_cast<unsigned int>(mSupportIndex[k+1][l]) == mSupport[l]) {
145 mSupportIndex[k+1][l] = 0;
153 // Compute BSpline weights if not up to date
154 if (!mWeightsAreUpToDate)
155 UpdatePrecomputedWeights();
157 //====================================================================
159 //====================================================================
160 template <class TImageType, class TCoordRep, class TCoefficientType>
161 void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
162 UpdatePrecomputedWeights() {
164 mWeightsCalculator.SetSplineOrders(mSplineOrders);
165 mWeightsCalculator.SetSamplingFactors(mSamplingFactors);
166 mWeightsCalculator.ComputeTensorProducts();
167 mWeightsAreUpToDate = true;
169 //====================================================================
171 //====================================================================
172 template <class TImageType, class TCoordRep, class TCoefficientType>
173 typename BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::IndexType
174 BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
175 GetSampleIndexOfPixelPosition(const ContinuousIndexType & x, IndexType & EvaluateIndex) const {
178 WARNING : sometimes, a floating number x could not really be
179 represented in memory. In this case, the difference between x and
180 floor(x) can be almost 1 (instead of 0). So we take into account
181 such special case, otherwise it could lead to annoying results.
183 // static const TCoefficientType tiny = 1.0e-7;
186 for(int l=0; l<TImageType::ImageDimension; l++)
188 // Compute t1 = distance to floor
189 TCoefficientType t1 = x[l]- vcl_floor(x[l]);
191 // Compute index in precomputed weights table
192 TCoefficientType t2 = mSamplingFactors[l]*t1;
193 index[l] = (IndexValueType)lrint(t2);
195 // For even order : test if too close to 0.5 (but lower). In this
196 // case : take the next coefficient
197 if (!(mSplineOrders[l] & 1)) {
199 if (mSamplingFactors[l] & 1) {
200 if (index[l] == (int) mSamplingFactors[l]/2+1) EvaluateIndex[l] = EvaluateIndex[l]+1;
203 else if (index[l] == (int) mSamplingFactors[l]/2) EvaluateIndex[l] = EvaluateIndex[l]+1;
207 // When to close to 1, take the next coefficient for odd order, but
208 // only change index for odd
209 if (index[l] == (int)mSamplingFactors[l]) {
211 if (mSplineOrders[l] & 1) EvaluateIndex[l] = EvaluateIndex[l]+1;
220 //====================================================================
222 //====================================================================
223 template <class TImageType, class TCoordRep, class TCoefficientType>
224 typename BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::OutputType
225 BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
226 EvaluateAtContinuousIndex(const ContinuousIndexType & x) const {
228 // JV Compute BSpline weights if not up to date! Problem const: pass image as last
229 // if (!mWeightsAreUpToDate) UpdatePrecomputedWeights();
231 // For shorter coding
232 static const unsigned int d = TImageType::ImageDimension;
234 // Compute the index of the first interpolation coefficient in the coefficient image
235 IndexType evaluateIndex;
237 for (unsigned int l=0; l<d; l++) {
238 if (mSplineOrders[l] & 1) { // Use this index calculation for odd splineOrder (like cubic)
239 indx = (long)vcl_floor(x[l]) - mSplineOrders[l] / 2 ; //this->m_SplineOrder / 2;
240 evaluateIndex[l] = indx;
242 else { // Use this index calculation for even splineOrder
243 if (mSplineOrders[l] == 0) evaluateIndex[l] = (long)rint(x[l]);
245 indx = (long)vcl_floor((x[l]+ 0.5)) - mSplineOrders[l] / 2; //this->m_SplineOrder / 2;
246 evaluateIndex[l] = indx;
251 // Compute index of precomputed weights and get pointer to first weights
252 const IndexType weightIndex = GetSampleIndexOfPixelPosition(x, evaluateIndex);
253 const TCoefficientType * pweights = mWeightsCalculator.GetFirstTensorProduct(weightIndex);
256 bool boundaryCase = false;
257 for (unsigned int l=0; l<d; l++) {
258 if ((evaluateIndex[l] < 0) ||
259 (evaluateIndex[l]+mSupport[l]) >= this->m_Coefficients->GetLargestPossibleRegion().GetSize(l)) {
264 // Pointer to support offset
265 const int * psupport;
267 // Special case for boundary (to be changed ....)
268 std::vector<int> correctedSupportOffset;
271 std::vector<IndexType> correctedSupportIndex;//(mSupportSize);
272 correctedSupportIndex.resize(mSupportSize);
273 correctedSupportOffset.resize(mSupportSize);
274 for(unsigned int i=0; i<mSupportSize; i++) {
275 for (unsigned int l=0; l<d; l++) {
276 long a = mSupportIndex[i][l] + evaluateIndex[l];
277 long b = this->m_Coefficients->GetLargestPossibleRegion().GetSize(l);
280 correctedSupportIndex[i][l] = -a - d2*(-a/d2) - evaluateIndex[l];//mSupportIndex[i][l]-a;
284 correctedSupportIndex[i][l] = d2 - a - evaluateIndex[l];
287 correctedSupportIndex[i][l] = mSupportIndex[i][l]; //a - d2*(a/d2) - EvaluateIndex[l];
291 correctedSupportIndex[i][l] = d2 - a - EvaluateIndex[l];//mSupportIndex[i][l] - (a-(b-1));
294 correctedSupportIndex[i][l] = mSupportIndex[i][l];
299 correctedSupportOffset[i] = itk::Index2Offset<TImageType::ImageDimension>(correctedSupportIndex[i], mInputMemoryOffset);
301 psupport = &correctedSupportOffset[0];
304 psupport = &mSupportOffset[0];
307 // Get pointer to first coefficient. Even if in some boundary cases,
308 // EvaluateIndex is out of the coefficient image,
309 const TCoefficientType * pcoef = &(this->m_Coefficients->GetPixel(evaluateIndex));
311 // Main loop over BSpline support
312 TCoefficientType interpolated = itk::NumericTraits<TCoefficientType>::Zero;
313 for (unsigned int p=0; p<mSupportSize; p++) {
314 interpolated += pcoef[*psupport] * (*pweights);
319 // Return interpolated value
320 return(interpolated);
322 //====================================================================
325 #endif //_ITKINTERPOLATEIMAGEFUNCTIONWITHLUT_TXX