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>
117 : s1_(derivativeOrder)
127 return exec(x, derivativeOrder());
135 result_type
operator()(first_argument_type x, second_argument_type derivative_order)
const 137 return exec(x, derivativeOrder() + derivative_order);
143 {
return operator()(x); }
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>
198 template <
int ORDER,
class T>
201 template <
int ORDER,
class T>
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>
224 const int r = ORDER / 2;
227 for(
int i = 0; i <= 2*r; ++i)
228 p[i] = spline(T(i-r));
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>
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>
274 explicit BSpline(
unsigned int derivativeOrder = 0)
297 explicit BSplineBase(
unsigned int derivativeOrder = 0)
298 : derivativeOrder_(derivativeOrder)
301 result_type operator()(argument_type x)
const 303 return exec(x, derivativeOrder_);
306 template <
unsigned int IntBits,
unsigned int FracBits>
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>
318 return x < 0.5 && -0.5 <= x
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);
328 value_type operator[](value_type x)
const 329 {
return operator()(x); }
331 double radius()
const 334 unsigned int derivativeOrder()
const 335 {
return derivativeOrder_; }
339 return prefilterCoefficients_;
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_;
362 static WeightMatrix weightMatrix_;
389 explicit BSpline(
unsigned int derivativeOrder = 0)
390 : derivativeOrder_(derivativeOrder)
393 result_type operator()(argument_type x)
const 395 return exec(x, derivativeOrder_);
398 template <
unsigned int IntBits,
unsigned int FracBits>
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>
417 result_type operator()(first_argument_type x, second_argument_type derivative_order)
const 419 return exec(x, derivativeOrder_ + derivative_order);
422 value_type operator[](value_type x)
const 423 {
return operator()(x); }
425 double radius()
const 428 unsigned int derivativeOrder()
const 429 {
return derivativeOrder_; }
433 return prefilterCoefficients_;
438 static WeightMatrix
const & weights()
440 return weightMatrix_;
444 T exec(T x,
unsigned int derivative_order)
const;
446 unsigned int derivativeOrder_;
448 static WeightMatrix weightMatrix_;
460 switch(derivative_order)
464 x = VIGRA_CSTD::fabs(x);
502 explicit BSpline(
unsigned int derivativeOrder = 0)
503 : derivativeOrder_(derivativeOrder)
506 result_type operator()(argument_type x)
const 508 return exec(x, derivativeOrder_);
511 template <
unsigned int IntBits,
unsigned int FracBits>
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>
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 548 {
return operator()(x, 1); }
550 value_type operator[](value_type x)
const 551 {
return operator()(x); }
553 double radius()
const 556 unsigned int derivativeOrder()
const 557 {
return derivativeOrder_; }
561 return prefilterCoefficients_;
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_;
576 static WeightMatrix weightMatrix_;
584 {{ 0.125, 0.75, 0.125},
592 switch(derivative_order)
596 x = VIGRA_CSTD::fabs(x);
650 explicit BSpline(
unsigned int derivativeOrder = 0)
651 : derivativeOrder_(derivativeOrder)
654 result_type operator()(argument_type x)
const 656 return exec(x, derivativeOrder_);
659 template <
unsigned int IntBits,
unsigned int FracBits>
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>
685 return 2.0/3.0 + x*x*(-1.0 + 0.5*x);
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 702 {
return operator()(x, 1); }
704 result_type dxx(argument_type x)
const 705 {
return operator()(x, 2); }
707 value_type operator[](value_type x)
const 708 {
return operator()(x); }
710 double radius()
const 713 unsigned int derivativeOrder()
const 714 {
return derivativeOrder_; }
718 return prefilterCoefficients_;
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_;
733 static WeightMatrix 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}};
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);
827 explicit BSpline(
unsigned int derivativeOrder = 0)
828 : derivativeOrder_(derivativeOrder)
831 result_type operator()(argument_type x)
const 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>
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;
862 result_type dx(argument_type x)
const 863 {
return operator()(x, 1); }
865 result_type dxx(argument_type x)
const 866 {
return operator()(x, 2); }
868 result_type dx3(argument_type x)
const 869 {
return operator()(x, 3); }
871 value_type operator[](value_type x)
const 872 {
return operator()(x); }
874 double radius()
const 877 unsigned int derivativeOrder()
const 878 {
return derivativeOrder_; }
882 return prefilterCoefficients_;
887 static WeightMatrix
const & weights()
889 return weightMatrix_;
893 result_type exec(first_argument_type x, second_argument_type derivative_order)
const;
899 b[0] = -0.361341225900220177092212841325;
901 b[1] = -0.013725429297339121360331226939;
905 unsigned int derivativeOrder_;
907 static WeightMatrix 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}};
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);
1051 explicit BSpline(
unsigned int derivativeOrder = 0)
1052 : derivativeOrder_(derivativeOrder)
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>
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;
1086 result_type dx(argument_type x)
const 1087 {
return operator()(x, 1); }
1089 result_type dxx(argument_type x)
const 1090 {
return operator()(x, 2); }
1092 result_type dx3(argument_type x)
const 1093 {
return operator()(x, 3); }
1095 result_type dx4(argument_type x)
const 1096 {
return operator()(x, 4); }
1098 value_type operator[](value_type x)
const 1099 {
return operator()(x); }
1101 double radius()
const 1104 unsigned int derivativeOrder()
const 1105 {
return derivativeOrder_; }
1109 return prefilterCoefficients_;
1114 static WeightMatrix
const & weights()
1116 return weightMatrix_;
1120 result_type exec(first_argument_type x, second_argument_type derivative_order)
const;
1126 b[0] = -0.430575347099973791851434783493;
1128 b[1] = -0.043096288203264653822712376822;
1132 unsigned int derivativeOrder_;
1134 static WeightMatrix 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}};
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;
1280 #endif // NO_PARTIAL_TEMPLATE_SPECIALIZATION 1308 template <
class T =
double>
1327 result_type operator()(argument_type x)
const;
1331 T operator[] (T x)
const 1332 {
return operator()(x); }
1350 return prefilterCoefficients_;
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));
StaticOrder
Definition: splines.hxx:1323
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
Definition: autodiff.hxx:87
ArrayVector< double > const & prefilterCoefficients() const
Definition: splines.hxx:164
Definition: polynomial.hxx:52
Definition: accessor.hxx:43
static WeightMatrix const & weights()
Definition: splines.hxx:176
Definition: splines.hxx:268
StaticOrder
Definition: splines.hxx:111
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:365
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
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
size_type size() const
Definition: array_vector.hxx:358
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
Definition: fixedpoint.hxx:47
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