36 #ifndef VIGRA_SPLINES_HXX
37 #define VIGRA_SPLINES_HXX
41 #include "mathutil.hxx"
42 #include "polynomial.hxx"
43 #include "array_vector.hxx"
44 #include "fixedpoint.hxx"
50 template <
class T,
int N>
63 #ifndef NO_PARTIAL_TEMPLATE_SPECIALIZATION
89 template <
int ORDER,
class T =
double>
135 result_type
operator()(first_argument_type x, second_argument_type derivative_order)
const
149 {
return (ORDER + 1) * 0.5; }
154 {
return s1_.derivativeOrder(); }
166 return prefilterCoefficients_;
178 return weightMatrix_;
182 result_type exec(first_argument_type x, second_argument_type derivative_order)
const;
188 static WeightMatrix calculateWeightMatrix();
192 static WeightMatrix weightMatrix_;
195 template <
int ORDER,
class T>
196 ArrayVector<double> BSplineBase<ORDER, T>::prefilterCoefficients_(BSplineBase<ORDER, T>::calculatePrefilterCoefficients());
198 template <
int ORDER,
class T>
199 typename BSplineBase<ORDER, T>::WeightMatrix BSplineBase<ORDER, T>::weightMatrix_(calculateWeightMatrix());
201 template <
int ORDER,
class T>
202 typename BSplineBase<ORDER, T>::result_type
203 BSplineBase<ORDER, T>::exec(first_argument_type x, second_argument_type derivative_order)
const
205 if(derivative_order == 0)
207 T n12 = (ORDER + 1.0) / 2.0;
208 return ((n12 + x) * s1_(x + 0.5) + (n12 - x) * s1_(x - 0.5)) / ORDER;
213 return s1_(x + 0.5, derivative_order) - s1_(x - 0.5, derivative_order);
217 template <
int ORDER,
class T>
219 BSplineBase<ORDER, T>::calculatePrefilterCoefficients()
221 ArrayVector<double> res;
224 const int r = ORDER / 2;
225 StaticPolynomial<2*r, double> p(2*r);
227 for(
int i = 0; i <= 2*r; ++i)
228 p[i] = spline(T(i-r));
229 ArrayVector<double> roots;
231 for(
unsigned int i = 0; i < roots.size(); ++i)
232 if(VIGRA_CSTD::fabs(roots[i]) < 1.0)
233 res.push_back(roots[i]);
238 template <
int ORDER,
class T>
239 typename BSplineBase<ORDER, T>::WeightMatrix
240 BSplineBase<ORDER, T>::calculateWeightMatrix()
242 WeightMatrix res(ORDER+1, ArrayVector<T>(ORDER+1));
243 double faculty = 1.0;
244 for(
int d = 0; d <= ORDER; ++d)
248 double x = ORDER / 2;
250 for(
int i = 0; i <= ORDER; ++i, --x)
251 res[d][i] = spline(x, d) / faculty;
267 template <
int ORDER,
class T =
double>
286 class BSplineBase<0, T>
303 return exec(x, derivativeOrder_);
306 template <
unsigned int IntBits,
unsigned int FracBits>
307 FixedPoint<IntBits, FracBits>
operator()(FixedPoint<IntBits, FracBits> x)
const
309 typedef FixedPoint<IntBits, FracBits> Value;
310 return x.value < Value::ONE_HALF && -Value::ONE_HALF <= x.value
311 ? Value(Value::ONE, FPNoShift)
312 : Value(0, FPNoShift);
315 template <
class U,
int N>
316 autodiff::DualVector<U, N>
operator()(autodiff::DualVector<U, N>
const & x)
const
318 return x < 0.5 && -0.5 <= x
319 ? autodiff::DualVector<U, N>(1.0)
320 : autodiff::DualVector<U, N>(0.0);
323 result_type
operator()(first_argument_type x, second_argument_type derivative_order)
const
325 return exec(x, derivativeOrder_ + derivative_order);
335 {
return derivativeOrder_; }
339 return prefilterCoefficients_;
342 typedef T WeightMatrix[1][1];
344 static WeightMatrix
const &
weights()
346 return weightMatrix_;
350 result_type exec(first_argument_type x, second_argument_type derivative_order)
const
352 if(derivative_order == 0)
353 return x < 0.5 && -0.5 <= x ?
360 unsigned int derivativeOrder_;
361 static ArrayVector<double> prefilterCoefficients_;
362 static WeightMatrix weightMatrix_;
366 ArrayVector<double> BSplineBase<0, T>::prefilterCoefficients_;
369 typename BSplineBase<0, T>::WeightMatrix BSplineBase<0, T>::weightMatrix_ = {{ 1.0 }};
395 return exec(x, derivativeOrder_);
398 template <
unsigned int IntBits,
unsigned int FracBits>
399 FixedPoint<IntBits, FracBits>
operator()(FixedPoint<IntBits, FracBits> x)
const
401 typedef FixedPoint<IntBits, FracBits> Value;
402 int v =
abs(x.value);
403 return v < Value::ONE ?
404 Value(Value::ONE - v, FPNoShift)
405 : Value(0, FPNoShift);
408 template <
class U,
int N>
409 autodiff::DualVector<U, N>
operator()(autodiff::DualVector<U, N> x)
const
414 : autodiff::DualVector<U, N>(0.0);
417 result_type
operator()(first_argument_type x, second_argument_type derivative_order)
const
419 return exec(x, derivativeOrder_ + derivative_order);
429 {
return derivativeOrder_; }
433 return prefilterCoefficients_;
436 typedef T WeightMatrix[2][2];
438 static WeightMatrix
const &
weights()
440 return weightMatrix_;
444 T exec(T x,
unsigned int derivative_order)
const;
446 unsigned int derivativeOrder_;
447 static ArrayVector<double> prefilterCoefficients_;
448 static WeightMatrix weightMatrix_;
452 ArrayVector<double> BSpline<1, T>::prefilterCoefficients_;
455 typename BSpline<1, T>::WeightMatrix BSpline<1, T>::weightMatrix_ = {{ 1.0, 0.0}, {-1.0, 1.0}};
458 T BSpline<1, T>::exec(T x,
unsigned int derivative_order)
const
460 switch(derivative_order)
464 x = VIGRA_CSTD::fabs(x);
508 return exec(x, derivativeOrder_);
511 template <
unsigned int IntBits,
unsigned int FracBits>
512 FixedPoint<IntBits, FracBits>
operator()(FixedPoint<IntBits, FracBits> x)
const
514 typedef FixedPoint<IntBits, FracBits> Value;
515 enum { ONE_HALF = Value::ONE_HALF, THREE_HALVES = ONE_HALF * 3, THREE_QUARTERS = THREE_HALVES / 2,
516 PREMULTIPLY_SHIFT1 = FracBits <= 16 ? 0 : FracBits - 16,
517 PREMULTIPLY_SHIFT2 = FracBits - 1 <= 16 ? 0 : FracBits - 17,
518 POSTMULTIPLY_SHIFT1 = FracBits - 2*PREMULTIPLY_SHIFT1,
519 POSTMULTIPLY_SHIFT2 = FracBits - 2*PREMULTIPLY_SHIFT2 };
520 int v =
abs(x.value);
522 ? Value(ONE_HALF, FPNoShift)
524 ? Value(THREE_QUARTERS -
525 (int)(
sq((unsigned)v >> PREMULTIPLY_SHIFT2) >> POSTMULTIPLY_SHIFT2), FPNoShift)
527 ? Value((int)(
sq((unsigned)(THREE_HALVES-v) >> PREMULTIPLY_SHIFT1) >> (POSTMULTIPLY_SHIFT1 + 1)), FPNoShift)
528 : Value(0, FPNoShift);
531 template <
class U,
int N>
532 autodiff::DualVector<U, N>
operator()(autodiff::DualVector<U, N> x)
const
539 : autodiff::DualVector<U, N>(0.0);
542 result_type
operator()(first_argument_type x, second_argument_type derivative_order)
const
544 return exec(x, derivativeOrder_ + derivative_order);
547 result_type dx(argument_type x)
const
557 {
return derivativeOrder_; }
561 return prefilterCoefficients_;
564 typedef T WeightMatrix[3][3];
566 static WeightMatrix
const &
weights()
568 return weightMatrix_;
572 result_type exec(first_argument_type x, second_argument_type derivative_order)
const;
574 unsigned int derivativeOrder_;
575 static ArrayVector<double> prefilterCoefficients_;
576 static WeightMatrix weightMatrix_;
580 ArrayVector<double> BSpline<2, T>::prefilterCoefficients_(1, 2.0*M_SQRT2 - 3.0);
583 typename BSpline<2, T>::WeightMatrix BSpline<2, T>::weightMatrix_ =
584 {{ 0.125, 0.75, 0.125},
589 typename BSpline<2, T>::result_type
590 BSpline<2, T>::exec(first_argument_type x, second_argument_type derivative_order)
const
592 switch(derivative_order)
596 x = VIGRA_CSTD::fabs(x);
656 return exec(x, derivativeOrder_);
659 template <
unsigned int IntBits,
unsigned int FracBits>
660 FixedPoint<IntBits, FracBits>
operator()(FixedPoint<IntBits, FracBits> x)
const
662 typedef FixedPoint<IntBits, FracBits> Value;
663 enum { ONE = Value::ONE, TWO = 2 * ONE, TWO_THIRDS = TWO / 3, ONE_SIXTH = ONE / 6,
664 PREMULTIPLY_SHIFT = FracBits <= 16 ? 0 : FracBits - 16,
665 POSTMULTIPLY_SHIFT = FracBits - 2*PREMULTIPLY_SHIFT };
666 int v =
abs(x.value);
668 ? Value(ONE_SIXTH, FPNoShift)
671 (((int)(
sq((unsigned)v >> PREMULTIPLY_SHIFT) >> (POSTMULTIPLY_SHIFT + PREMULTIPLY_SHIFT))
672 * (((v >> 1) - ONE) >> PREMULTIPLY_SHIFT)) >> POSTMULTIPLY_SHIFT), FPNoShift)
674 ? Value((int)((
sq((unsigned)(TWO-v) >> PREMULTIPLY_SHIFT) >> (POSTMULTIPLY_SHIFT + PREMULTIPLY_SHIFT))
675 * ((unsigned)(TWO-v) >> PREMULTIPLY_SHIFT) / 6) >> POSTMULTIPLY_SHIFT, FPNoShift)
676 : Value(0, FPNoShift);
679 template <
class U,
int N>
680 autodiff::DualVector<U, N>
operator()(autodiff::DualVector<U, N> x)
const
685 return 2.0/3.0 + x*x*(-1.0 + 0.5*x);
693 return autodiff::DualVector<U, N>(0.0);
696 result_type
operator()(first_argument_type x, second_argument_type derivative_order)
const
698 return exec(x, derivativeOrder_ + derivative_order);
701 result_type dx(argument_type x)
const
704 result_type dxx(argument_type x)
const
714 {
return derivativeOrder_; }
718 return prefilterCoefficients_;
721 typedef T WeightMatrix[4][4];
723 static WeightMatrix
const &
weights()
725 return weightMatrix_;
729 result_type exec(first_argument_type x, second_argument_type derivative_order)
const;
731 unsigned int derivativeOrder_;
732 static ArrayVector<double> prefilterCoefficients_;
733 static WeightMatrix weightMatrix_;
737 ArrayVector<double> BSpline<3, T>::prefilterCoefficients_(1,
VIGRA_CSTD::sqrt(3.0) - 2.0);
740 typename BSpline<3, T>::WeightMatrix BSpline<3, T>::weightMatrix_ =
741 {{ 1.0 / 6.0, 2.0 / 3.0, 1.0 / 6.0, 0.0},
742 {-0.5, 0.0, 0.5, 0.0},
743 { 0.5, -1.0, 0.5, 0.0},
744 {-1.0 / 6.0, 0.5, -0.5, 1.0 / 6.0}};
747 typename BSpline<3, T>::result_type
748 BSpline<3, T>::exec(first_argument_type x, second_argument_type derivative_order)
const
750 switch(derivative_order)
754 x = VIGRA_CSTD::fabs(x);
757 return 2.0/3.0 + x*x*(-1.0 + 0.5*x);
772 x = VIGRA_CSTD::fabs(x);
781 x = VIGRA_CSTD::fabs(x);
807 typedef BSpline<3, double> CubicBSplineKernel;
833 return exec(x, derivativeOrder_);
836 result_type
operator()(first_argument_type x, second_argument_type derivative_order)
const
838 return exec(x, derivativeOrder_ + derivative_order);
841 template <
class U,
int N>
842 autodiff::DualVector<U, N>
operator()(autodiff::DualVector<U, N> x)
const
847 return 115.0/192.0 + x*x*(-0.625 + x*x*0.25);
851 return (55.0/16.0 + x*(1.25 + x*(-7.5 + x*(5.0 - x)))) / 6.0;
856 return sq(x*x) / 24.0;
859 return autodiff::DualVector<U, N>(0.0);
862 result_type dx(argument_type x)
const
865 result_type dxx(argument_type x)
const
868 result_type dx3(argument_type x)
const
878 {
return derivativeOrder_; }
882 return prefilterCoefficients_;
885 typedef T WeightMatrix[5][5];
887 static WeightMatrix
const &
weights()
889 return weightMatrix_;
893 result_type exec(first_argument_type x, second_argument_type derivative_order)
const;
895 static ArrayVector<double> calculatePrefilterCoefficients()
897 ArrayVector<double> b(2);
899 b[0] = -0.361341225900220177092212841325;
901 b[1] = -0.013725429297339121360331226939;
905 unsigned int derivativeOrder_;
906 static ArrayVector<double> prefilterCoefficients_;
907 static WeightMatrix weightMatrix_;
911 ArrayVector<double> BSpline<4, T>::prefilterCoefficients_(calculatePrefilterCoefficients());
914 typename BSpline<4, T>::WeightMatrix BSpline<4, T>::weightMatrix_ =
915 {{ 1.0/384.0, 19.0/96.0, 115.0/192.0, 19.0/96.0, 1.0/384.0},
916 {-1.0/48.0, -11.0/24.0, 0.0, 11.0/24.0, 1.0/48.0},
917 { 1.0/16.0, 1.0/4.0, -5.0/8.0, 1.0/4.0, 1.0/16.0},
918 {-1.0/12.0, 1.0/6.0, 0.0, -1.0/6.0, 1.0/12.0},
919 { 1.0/24.0, -1.0/6.0, 0.25, -1.0/6.0, 1.0/24.0}};
922 typename BSpline<4, T>::result_type
923 BSpline<4, T>::exec(first_argument_type x, second_argument_type derivative_order)
const
925 switch(derivative_order)
929 x = VIGRA_CSTD::fabs(x);
932 return 115.0/192.0 + x*x*(-0.625 + x*x*0.25);
936 return (55.0/16.0 + x*(1.25 + x*(-7.5 + x*(5.0 - x)))) / 6.0;
941 return sq(x*x) / 24.0;
951 x = VIGRA_CSTD::fabs(x);
954 return s*x*(-1.25 + x*x);
958 return s*(5.0 + x*(-60.0 + x*(60.0 - 16.0*x))) / 24.0;
963 return s*x*x*x / -6.0;
970 x = VIGRA_CSTD::fabs(x);
973 return -1.25 + 3.0*x*x;
977 return -2.5 + x*(5.0 - 2.0*x);
992 x = VIGRA_CSTD::fabs(x);
999 return s*(5.0 - 4.0*x);
1031 typedef BSpline<4, double> QuarticBSplineKernel;
1055 result_type
operator()(argument_type x)
const
1057 return exec(x, derivativeOrder_);
1060 result_type
operator()(first_argument_type x, second_argument_type derivative_order)
const
1062 return exec(x, derivativeOrder_ + derivative_order);
1065 template <
class U,
int N>
1066 autodiff::DualVector<U, N>
operator()(autodiff::DualVector<U, N> x)
const
1071 return 0.55 + x*x*(-0.5 + x*x*(0.25 - x/12.0));
1075 return 17.0/40.0 + x*(0.625 + x*(-1.75 + x*(1.25 + x*(-0.375 + x/24.0))));
1080 return x*
sq(x*x) / 120.0;
1083 return autodiff::DualVector<U, N>(0.0);
1086 result_type dx(argument_type x)
const
1089 result_type dxx(argument_type x)
const
1092 result_type dx3(argument_type x)
const
1095 result_type dx4(argument_type x)
const
1105 {
return derivativeOrder_; }
1109 return prefilterCoefficients_;
1112 typedef T WeightMatrix[6][6];
1114 static WeightMatrix
const &
weights()
1116 return weightMatrix_;
1120 result_type exec(first_argument_type x, second_argument_type derivative_order)
const;
1122 static ArrayVector<double> calculatePrefilterCoefficients()
1124 ArrayVector<double> b(2);
1126 b[0] = -0.430575347099973791851434783493;
1128 b[1] = -0.043096288203264653822712376822;
1132 unsigned int derivativeOrder_;
1133 static ArrayVector<double> prefilterCoefficients_;
1134 static WeightMatrix weightMatrix_;
1138 ArrayVector<double> BSpline<5, T>::prefilterCoefficients_(calculatePrefilterCoefficients());
1141 typename BSpline<5, T>::WeightMatrix BSpline<5, T>::weightMatrix_ =
1142 {{ 1.0/120.0, 13.0/60.0, 11.0/20.0, 13.0/60.0, 1.0/120.0, 0.0},
1143 {-1.0/24.0, -5.0/12.0, 0.0, 5.0/12.0, 1.0/24.0, 0.0},
1144 { 1.0/12.0, 1.0/6.0, -0.5, 1.0/6.0, 1.0/12.0, 0.0},
1145 {-1.0/12.0, 1.0/6.0, 0.0, -1.0/6.0, 1.0/12.0, 0.0},
1146 { 1.0/24.0, -1.0/6.0, 0.25, -1.0/6.0, 1.0/24.0, 0.0},
1147 {-1.0/120.0, 1.0/24.0, -1.0/12.0, 1.0/12.0, -1.0/24.0, 1.0/120.0}};
1150 typename BSpline<5, T>::result_type
1151 BSpline<5, T>::exec(first_argument_type x, second_argument_type derivative_order)
const
1153 switch(derivative_order)
1157 x = VIGRA_CSTD::fabs(x);
1160 return 0.55 + x*x*(-0.5 + x*x*(0.25 - x/12.0));
1164 return 17.0/40.0 + x*(0.625 + x*(-1.75 + x*(1.25 + x*(-0.375 + x/24.0))));
1169 return x*
sq(x*x) / 120.0;
1176 double s = x < 0.0 ?
1179 x = VIGRA_CSTD::fabs(x);
1182 return s*x*(-1.0 + x*x*(1.0 - 5.0/12.0*x));
1186 return s*(0.625 + x*(-3.5 + x*(3.75 + x*(-1.5 + 5.0/24.0*x))));
1191 return s*
sq(x*x) / -24.0;
1198 x = VIGRA_CSTD::fabs(x);
1201 return -1.0 + x*x*(3.0 -5.0/3.0*x);
1205 return -3.5 + x*(7.5 + x*(-4.5 + 5.0/6.0*x));
1217 double s = x < 0.0 ?
1220 x = VIGRA_CSTD::fabs(x);
1223 return s*x*(6.0 - 5.0*x);
1227 return s*(7.5 + x*(-9.0 + 2.5*x));
1239 x = VIGRA_CSTD::fabs(x);
1242 return 6.0 - 10.0*x;
1246 return -9.0 + 5.0*x;
1278 typedef BSpline<5, double> QuinticBSplineKernel;
1280 #endif // NO_PARTIAL_TEMPLATE_SPECIALIZATION
1308 template <
class T =
double>
1327 result_type
operator()(argument_type x)
const;
1350 return prefilterCoefficients_;
1358 ArrayVector<double> CatmullRomSpline<T>::prefilterCoefficients_;
1361 typename CatmullRomSpline<T>::result_type
1364 x = VIGRA_CSTD::fabs(x);
1367 return 1.0 + x * x * (-2.5 + 1.5 * x);
1375 return 2.0 + x * (-4.0 + x * (2.5 - 0.5 * x));
T operator[](T x) const
Definition: splines.hxx:1331
double radius() const
Definition: splines.hxx:148
value_type operator[](value_type x) const
Definition: splines.hxx:142
T result_type
Definition: splines.hxx:1320
T argument_type
Definition: splines.hxx:1317
ArrayVector< double > const & prefilterCoefficients() const
Definition: splines.hxx:164
Definition: accessor.hxx:43
static WeightMatrix const & weights()
Definition: splines.hxx:176
Definition: splines.hxx:268
Definition: splines.hxx:90
T argument_type
Definition: splines.hxx:99
T first_argument_type
Definition: splines.hxx:102
NumericTraits< T >::Promote sq(T t)
The square function.
Definition: mathutil.hxx:344
result_type operator()(argument_type x) const
Definition: splines.hxx:125
int radius() const
Definition: splines.hxx:1337
unsigned int derivativeOrder() const
Definition: splines.hxx:1342
T value_type
Definition: splines.hxx:1314
bool polynomialRealRoots(POLYNOMIAL const &p, VECTOR &roots, bool polishRoots)
Definition: polynomial.hxx:1076
StaticOrder
Definition: splines.hxx:111
unsigned int derivativeOrder() const
Definition: splines.hxx:153
BSplineBase(unsigned int derivativeOrder=0)
Definition: splines.hxx:116
unsigned int second_argument_type
Definition: splines.hxx:105
ArrayVector< double > const & prefilterCoefficients() const
Definition: splines.hxx:1348
T result_type
Definition: splines.hxx:108
FFTWComplex< R >::NormType abs(const FFTWComplex< R > &a)
absolute value (= magnitude)
Definition: fftw3.hxx:1002
BSpline(unsigned int derivativeOrder=0)
Definition: splines.hxx:274
result_type operator()(first_argument_type x, second_argument_type derivative_order) const
Definition: splines.hxx:135
T value_type
Definition: splines.hxx:96
result_type operator()(argument_type x) const
Definition: splines.hxx:1362
SquareRootTraits< FixedPoint< IntBits, FracBits > >::SquareRootResult sqrt(FixedPoint< IntBits, FracBits > v)
square root.
Definition: fixedpoint.hxx:616
Definition: splines.hxx:1309