- //====================================================================
-
- //====================================================================
- template <class TImageType, class TCoordRep, class TCoefficientType>
- void VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
- UpdatePrecomputedWeights() {
- // mLUTTimer.Reset();
- // mLUTTimer.Start();
- mWeightsCalculator.SetSplineOrders(mSplineOrders);
- mWeightsCalculator.SetSamplingFactors(mSamplingFactors);
- mWeightsCalculator.ComputeTensorProducts();
- mWeightsAreUpToDate = true;
- //DS
- // coef = new TCoefficientType[mSupportSize];
- // mCorrectedSupportIndex.resize(mSupportSize);
- // mCorrectedSupportOffset.resize(mSupportSize);
- // mLUTTimer.Stop();
- // mLUTTimer.Print("LUT \t");
- }
- //====================================================================
-
- //====================================================================
- template <class TImageType, class TCoordRep, class TCoefficientType>
- typename VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::IndexType
- VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
- GetSampleIndexOfPixelPosition(const ContinuousIndexType & x, IndexType & EvaluateIndex) const {
-
- /*
- WARNING : sometimes, a floating number x could not really be
- represented in memory. In this case, the difference between x and
- floor(x) can be almost 1 (instead of 0). So we take into account
- such special case, otherwise it could lead to annoying results.
- */
- // static const TCoefficientType tiny = 1.0e-7;
- IndexType index;
-
- for(int l=0; l<TImageType::ImageDimension; l++)
- {
- // bool mChange = false;
-
- // Compute t1 = distance to floor
- TCoefficientType t1 = x[l]- vcl_floor(x[l]);
-
- // Compute index in precomputed weights table
- TCoefficientType t2 = mSamplingFactors[l]*t1;
- index[l] = (IndexValueType)lrint(t2);
-
- // For even order : test if too close to 0.5 (but lower). In this
- // case : take the next coefficient
- if (!(mSplineOrders[l] & 1)) {
- if (t1<0.5) {
- if (mSamplingFactors[l] & 1) {
- if (index[l] == (int) mSamplingFactors[l]/2+1) EvaluateIndex[l] = EvaluateIndex[l]+1;
- }
-
- else if (index[l] == (int) mSamplingFactors[l]/2) EvaluateIndex[l] = EvaluateIndex[l]+1;
- }
- }
-
- // Statistics (to be removed)
- /*
- *mIntrinsecError += fabs(index[l]-t2);
- (*mNumberOfError)++;
- if (fabs(index[l]-t2)> *mIntrinsecErrorMax) *mIntrinsecErrorMax = fabs(index[l]-t2);
- */
-
- // When to close to 1, take the next coefficient for odd order, but
- // only change index for odd
- if (index[l] == (int)mSamplingFactors[l]) {
- index[l] = 0;
- if (mSplineOrders[l] & 1) EvaluateIndex[l] = EvaluateIndex[l]+1;
- }
+//====================================================================
+
+//====================================================================
+template <class TImageType, class TCoordRep, class TCoefficientType>
+void VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
+UpdatePrecomputedWeights()
+{
+ mWeightsCalculator.SetSplineOrders(mSplineOrders);
+ mWeightsCalculator.SetSamplingFactors(mSamplingFactors);
+ mWeightsCalculator.ComputeTensorProducts();
+ mWeightsAreUpToDate = true;
+}
+//====================================================================
+
+//====================================================================
+template <class TImageType, class TCoordRep, class TCoefficientType>
+typename VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::IndexType
+VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
+GetSampleIndexOfPixelPosition(const ContinuousIndexType & x, IndexType & EvaluateIndex) const
+{
+
+ /*
+ WARNING : sometimes, a floating number x could not really be
+ represented in memory. In this case, the difference between x and
+ floor(x) can be almost 1 (instead of 0). So we take into account
+ such special case, otherwise it could lead to annoying results.
+ */
+ // static const TCoefficientType tiny = 1.0e-7;
+ IndexType index;
+
+ for(int l=0; l<TImageType::ImageDimension; l++) {
+ // Compute t1 = distance to floor
+ TCoefficientType t1 = x[l]- vcl_floor(x[l]);
+
+ // Compute index in precomputed weights table
+ TCoefficientType t2 = mSamplingFactors[l]*t1;
+ index[l] = (IndexValueType)itk::Math::Floor<IndexValueType,TCoefficientType>(t2);
+
+ // For even order : test if too close to 0.5 (but lower). In this
+ // case : take the next coefficient
+ if (!(mSplineOrders[l] & 1)) {
+ if (t1<0.5) {
+ if (mSamplingFactors[l] & 1) {
+ if (index[l] == (int) mSamplingFactors[l]/2+1) EvaluateIndex[l] = EvaluateIndex[l]+1;
+ }
+
+ else if (index[l] == (int) mSamplingFactors[l]/2) EvaluateIndex[l] = EvaluateIndex[l]+1;