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"
56 #ifndef NO_PARTIAL_TEMPLATE_SPECIALIZATION
82 template <
int ORDER,
class T =
double>
142 {
return (ORDER + 1) * 0.5; }
147 {
return s1_.derivativeOrder(); }
165 typedef T WeightMatrix[ORDER+1][ORDER+1];
174 static WeightMatrix & b = calculateWeightMatrix();
178 static WeightMatrix & calculateWeightMatrix();
186 template <
int ORDER,
class T>
187 typename BSplineBase<ORDER, T>::result_type
188 BSplineBase<ORDER, T>::exec(first_argument_type x, second_argument_type derivative_order)
const
190 if(derivative_order == 0)
192 T n12 = (ORDER + 1.0) / 2.0;
193 return ((n12 + x) * s1_(x + 0.5) + (n12 - x) * s1_(x - 0.5)) / ORDER;
198 return s1_(x + 0.5, derivative_order) - s1_(x - 0.5, derivative_order);
202 template <
int ORDER,
class T>
203 ArrayVector<double>
const & BSplineBase<ORDER, T>::calculatePrefilterCoefficients()
205 static ArrayVector<double> b;
208 static const int r = ORDER / 2;
209 StaticPolynomial<2*r, double> p(2*r);
211 for(
int i = 0; i <= 2*r; ++i)
212 p[i] = spline(T(i-r));
213 ArrayVector<double> roots;
215 for(
unsigned int i = 0; i < roots.size(); ++i)
216 if(VIGRA_CSTD::fabs(roots[i]) < 1.0)
217 b.push_back(roots[i]);
222 template <
int ORDER,
class T>
223 typename BSplineBase<ORDER, T>::WeightMatrix &
224 BSplineBase<ORDER, T>::calculateWeightMatrix()
226 static WeightMatrix b;
227 double faculty = 1.0;
228 for(
int d = 0; d <= ORDER; ++d)
232 double x = ORDER / 2;
234 for(
int i = 0; i <= ORDER; ++i, --x)
235 b[d][i] = spline(x, d) / faculty;
251 template <
int ORDER,
class T =
double>
270 class BSplineBase<0, T>
287 return exec(x, derivativeOrder_);
290 template <
unsigned int IntBits,
unsigned int FracBits>
291 FixedPoint<IntBits, FracBits>
operator()(FixedPoint<IntBits, FracBits> x)
const
293 typedef FixedPoint<IntBits, FracBits> Value;
294 return x.value < Value::ONE_HALF && -Value::ONE_HALF <= x.value
295 ? Value(Value::ONE, FPNoShift)
296 : Value(0, FPNoShift);
301 return exec(x, derivativeOrder_ + derivative_order);
311 {
return derivativeOrder_; }
315 static ArrayVector<double> b;
319 typedef T WeightMatrix[1][1];
320 static WeightMatrix &
weights()
322 static T b[1][1] = {{ 1.0}};
329 if(derivative_order == 0)
330 return x < 0.5 && -0.5 <= x ?
337 unsigned int derivativeOrder_;
364 return exec(x, derivativeOrder_);
367 template <
unsigned int IntBits,
unsigned int FracBits>
368 FixedPoint<IntBits, FracBits>
operator()(FixedPoint<IntBits, FracBits> x)
const
370 typedef FixedPoint<IntBits, FracBits> Value;
371 int v =
abs(x.value);
372 return v < Value::ONE ?
373 Value(Value::ONE - v, FPNoShift)
374 : Value(0, FPNoShift);
379 return exec(x, derivativeOrder_ + derivative_order);
389 {
return derivativeOrder_; }
393 static ArrayVector<double> b;
397 typedef T WeightMatrix[2][2];
398 static WeightMatrix &
weights()
400 static T b[2][2] = {{ 1.0, 0.0}, {-1.0, 1.0}};
405 T exec(T x,
unsigned int derivative_order)
const;
407 unsigned int derivativeOrder_;
411 T BSpline<1, T>::exec(T x,
unsigned int derivative_order)
const
413 switch(derivative_order)
417 x = VIGRA_CSTD::fabs(x);
461 return exec(x, derivativeOrder_);
464 template <
unsigned int IntBits,
unsigned int FracBits>
465 FixedPoint<IntBits, FracBits>
operator()(FixedPoint<IntBits, FracBits> x)
const
467 typedef FixedPoint<IntBits, FracBits> Value;
468 enum { ONE_HALF = Value::ONE_HALF, THREE_HALVES = ONE_HALF * 3, THREE_QUARTERS = THREE_HALVES / 2,
469 PREMULTIPLY_SHIFT1 = FracBits <= 16 ? 0 : FracBits - 16,
470 PREMULTIPLY_SHIFT2 = FracBits - 1 <= 16 ? 0 : FracBits - 17,
471 POSTMULTIPLY_SHIFT1 = FracBits - 2*PREMULTIPLY_SHIFT1,
472 POSTMULTIPLY_SHIFT2 = FracBits - 2*PREMULTIPLY_SHIFT2 };
473 int v =
abs(x.value);
475 ? Value(ONE_HALF, FPNoShift)
477 ? Value(THREE_QUARTERS -
478 (int)(
sq((unsigned)v >> PREMULTIPLY_SHIFT2) >> POSTMULTIPLY_SHIFT2), FPNoShift)
480 ? Value((int)(
sq((unsigned)(THREE_HALVES-v) >> PREMULTIPLY_SHIFT1) >> (POSTMULTIPLY_SHIFT1 + 1)), FPNoShift)
481 : Value(0, FPNoShift);
486 return exec(x, derivativeOrder_ + derivative_order);
499 {
return derivativeOrder_; }
503 static ArrayVector<double> b(1, 2.0*M_SQRT2 - 3.0);
507 typedef T WeightMatrix[3][3];
508 static WeightMatrix &
weights()
510 static T b[3][3] = {{ 0.125, 0.75, 0.125},
519 unsigned int derivativeOrder_;
523 typename BSpline<2, T>::result_type
524 BSpline<2, T>::exec(first_argument_type x, second_argument_type derivative_order)
const
526 switch(derivative_order)
530 x = VIGRA_CSTD::fabs(x);
590 return exec(x, derivativeOrder_);
593 template <
unsigned int IntBits,
unsigned int FracBits>
594 FixedPoint<IntBits, FracBits>
operator()(FixedPoint<IntBits, FracBits> x)
const
596 typedef FixedPoint<IntBits, FracBits> Value;
597 enum { ONE = Value::ONE, TWO = 2 * ONE, TWO_THIRDS = TWO / 3, ONE_SIXTH = ONE / 6,
598 PREMULTIPLY_SHIFT = FracBits <= 16 ? 0 : FracBits - 16,
599 POSTMULTIPLY_SHIFT = FracBits - 2*PREMULTIPLY_SHIFT };
600 int v =
abs(x.value);
602 ? Value(ONE_SIXTH, FPNoShift)
605 (((int)(
sq((unsigned)v >> PREMULTIPLY_SHIFT) >> (POSTMULTIPLY_SHIFT + PREMULTIPLY_SHIFT))
606 * (((v >> 1) - ONE) >> PREMULTIPLY_SHIFT)) >> POSTMULTIPLY_SHIFT), FPNoShift)
608 ? Value((int)((
sq((unsigned)(TWO-v) >> PREMULTIPLY_SHIFT) >> (POSTMULTIPLY_SHIFT + PREMULTIPLY_SHIFT))
609 * ((unsigned)(TWO-v) >> PREMULTIPLY_SHIFT) / 6) >> POSTMULTIPLY_SHIFT, FPNoShift)
610 : Value(0, FPNoShift);
615 return exec(x, derivativeOrder_ + derivative_order);
631 {
return derivativeOrder_; }
639 typedef T WeightMatrix[4][4];
640 static WeightMatrix &
weights()
642 static T b[4][4] = {{ 1.0 / 6.0, 2.0 / 3.0, 1.0 / 6.0, 0.0},
643 {-0.5, 0.0, 0.5, 0.0},
644 { 0.5, -1.0, 0.5, 0.0},
645 {-1.0 / 6.0, 0.5, -0.5, 1.0 / 6.0}};
652 unsigned int derivativeOrder_;
656 typename BSpline<3, T>::result_type
657 BSpline<3, T>::exec(first_argument_type x, second_argument_type derivative_order)
const
659 switch(derivative_order)
663 x = VIGRA_CSTD::fabs(x);
666 return 2.0/3.0 + x*x*(-1.0 + 0.5*x);
681 x = VIGRA_CSTD::fabs(x);
690 x = VIGRA_CSTD::fabs(x);
716 typedef BSpline<3, double> CubicBSplineKernel;
742 return exec(x, derivativeOrder_);
747 return exec(x, derivativeOrder_ + derivative_order);
766 {
return derivativeOrder_; }
770 static ArrayVector<double>
const & b = initPrefilterCoefficients();
774 static ArrayVector<double>
const & initPrefilterCoefficients()
776 static ArrayVector<double> b(2);
778 b[0] = -0.361341225900220177092212841325;
780 b[1] = -0.013725429297339121360331226939;
784 typedef T WeightMatrix[5][5];
785 static WeightMatrix &
weights()
787 static T b[5][5] = {{ 1.0/384.0, 19.0/96.0, 115.0/192.0, 19.0/96.0, 1.0/384.0},
788 {-1.0/48.0, -11.0/24.0, 0.0, 11.0/24.0, 1.0/48.0},
789 { 1.0/16.0, 1.0/4.0, -5.0/8.0, 1.0/4.0, 1.0/16.0},
790 {-1.0/12.0, 1.0/6.0, 0.0, -1.0/6.0, 1.0/12.0},
791 { 1.0/24.0, -1.0/6.0, 0.25, -1.0/6.0, 1.0/24.0}};
798 unsigned int derivativeOrder_;
802 typename BSpline<4, T>::result_type
803 BSpline<4, T>::exec(first_argument_type x, second_argument_type derivative_order)
const
805 switch(derivative_order)
809 x = VIGRA_CSTD::fabs(x);
812 return 115.0/192.0 + x*x*(-0.625 + x*x*0.25);
816 return (55.0/16.0 + x*(1.25 + x*(-7.5 + x*(5.0 - x)))) / 6.0;
821 return sq(x*x) / 24.0;
831 x = VIGRA_CSTD::fabs(x);
834 return s*x*(-1.25 + x*x);
838 return s*(5.0 + x*(-60.0 + x*(60.0 - 16.0*x))) / 24.0;
843 return s*x*x*x / -6.0;
850 x = VIGRA_CSTD::fabs(x);
853 return -1.25 + 3.0*x*x;
857 return -2.5 + x*(5.0 - 2.0*x);
872 x = VIGRA_CSTD::fabs(x);
879 return s*(5.0 - 4.0*x);
911 typedef BSpline<4, double> QuarticBSplineKernel;
937 return exec(x, derivativeOrder_);
942 return exec(x, derivativeOrder_ + derivative_order);
964 {
return derivativeOrder_; }
968 static ArrayVector<double>
const & b = initPrefilterCoefficients();
972 static ArrayVector<double>
const & initPrefilterCoefficients()
974 static ArrayVector<double> b(2);
976 b[0] = -0.430575347099973791851434783493;
978 b[1] = -0.043096288203264653822712376822;
982 typedef T WeightMatrix[6][6];
983 static WeightMatrix &
weights()
985 static T b[6][6] = {{ 1.0/120.0, 13.0/60.0, 11.0/20.0, 13.0/60.0, 1.0/120.0, 0.0},
986 {-1.0/24.0, -5.0/12.0, 0.0, 5.0/12.0, 1.0/24.0, 0.0},
987 { 1.0/12.0, 1.0/6.0, -0.5, 1.0/6.0, 1.0/12.0, 0.0},
988 {-1.0/12.0, 1.0/6.0, 0.0, -1.0/6.0, 1.0/12.0, 0.0},
989 { 1.0/24.0, -1.0/6.0, 0.25, -1.0/6.0, 1.0/24.0, 0.0},
990 {-1.0/120.0, 1.0/24.0, -1.0/12.0, 1.0/12.0, -1.0/24.0, 1.0/120.0}};
997 unsigned int derivativeOrder_;
1001 typename BSpline<5, T>::result_type
1002 BSpline<5, T>::exec(first_argument_type x, second_argument_type derivative_order)
const
1004 switch(derivative_order)
1008 x = VIGRA_CSTD::fabs(x);
1011 return 0.55 + x*x*(-0.5 + x*x*(0.25 - x/12.0));
1015 return 17.0/40.0 + x*(0.625 + x*(-1.75 + x*(1.25 + x*(-0.375 + x/24.0))));
1020 return x*
sq(x*x) / 120.0;
1027 double s = x < 0.0 ?
1030 x = VIGRA_CSTD::fabs(x);
1033 return s*x*(-1.0 + x*x*(1.0 - 5.0/12.0*x));
1037 return s*(0.625 + x*(-3.5 + x*(3.75 + x*(-1.5 + 5.0/24.0*x))));
1042 return s*
sq(x*x) / -24.0;
1049 x = VIGRA_CSTD::fabs(x);
1052 return -1.0 + x*x*(3.0 -5.0/3.0*x);
1056 return -3.5 + x*(7.5 + x*(-4.5 + 5.0/6.0*x));
1068 double s = x < 0.0 ?
1071 x = VIGRA_CSTD::fabs(x);
1074 return s*x*(6.0 - 5.0*x);
1078 return s*(7.5 + x*(-9.0 + 2.5*x));
1090 x = VIGRA_CSTD::fabs(x);
1093 return 6.0 - 10.0*x;
1097 return -9.0 + 5.0*x;
1129 typedef BSpline<5, double> QuinticBSplineKernel;
1131 #endif // NO_PARTIAL_TEMPLATE_SPECIALIZATION
1159 template <
class T =
double>
1208 typename CatmullRomSpline<T>::result_type
1211 x = VIGRA_CSTD::fabs(x);
1214 return 1.0 + x * x * (-2.5 + 1.5 * x);
1222 return 2.0 + x * (-4.0 + x * (2.5 - 0.5 * x));