36 #ifndef VIGRA_SPLINES_HXX
37 #define VIGRA_SPLINES_HXX
41 #include "mathutil.hxx"
42 #include "array_vector.hxx"
43 #include "fixedpoint.hxx"
49 template <
class T,
int N>
62 #ifndef NO_PARTIAL_TEMPLATE_SPECIALIZATION
88 template <
int ORDER,
class T =
double>
93 static_assert (ORDER < 18 ,
"BSpline: ORDER must be less than 18." );
150 {
return (ORDER + 1) * 0.5; }
155 {
return s1_.derivativeOrder(); }
167 return prefilterCoefficients_;
179 return weightMatrix_;
188 static WeightMatrix calculateWeightMatrix();
192 static WeightMatrix weightMatrix_;
195 template <
int ORDER,
class T>
196 ArrayVector<double> BSplineBase<ORDER, T>::prefilterCoefficients_(getPrefilterCoefficients());
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>::getPrefilterCoefficients()
221 static const double coeffs[18][8] = {
222 { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
223 { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
224 { -0.17157287525380971, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
225 { -0.26794919243112281, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
226 { -0.36134122590022018, -0.01372542929733912, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
227 { -0.43057534709997379, -0.04309628820326465, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
228 { -0.48829458930304398, -0.081679271076237972, -0.0014141518083258175, 0.0, 0.0, 0.0, 0.0, 0.0 },
229 { -0.53528043079643672, -0.1225546151923274, -0.0091486948096082786, 0.0, 0.0, 0.0, 0.0, 0.0 },
230 { -0.57468690924876631, -0.16303526929728085, -0.023632294694844857, -0.00015382131064169087, 0.0, 0.0, 0.0, 0.0 },
231 { -0.60799738916862989, -0.20175052019315337, -0.043222608540481752, -0.0021213069031808186, 0.0, 0.0, 0.0, 0.0 },
232 { -0.63655066396942439, -0.2381827983775629, -0.065727033228308585, -0.0075281946755486927, -1.6982762823274658e-5, 0.0, 0.0, 0.0 },
233 { -0.66126606890072925, -0.27218034929478602, -0.089759599793713341, -0.016669627366234657, -0.00051055753444650205, 0.0, 0.0, 0.0 },
234 { -0.68286488419772362, -0.30378079328825425, -0.11435052002713579, -0.028836190198663809, -0.0025161662172613372, -1.8833056450639017e-6, 0.0, 0.0 },
235 { -0.70189425181681642, -0.33310723293062366, -0.13890111319431958, -0.043213866740363663, -0.0067380314152449142, -0.00012510011321441875, 0.0, 0.0 },
236 { -0.71878378723997516, -0.3603190719169625, -0.1630335147992984, -0.059089482194831018, -0.013246756734847919, -0.00086402404095333829, -2.0913096775275374e-7, 0.0 },
237 { -0.73387257168487741, -0.38558573427843323, -0.18652010845096478, -0.075907592047668185, -0.02175206579654047, -0.0028011514820764556, -3.093568045147443e-5, 0.0 },
238 { -0.747432387772212103, -0.409073604757528353, -0.29228719338953817, -9.32547189803214355e-2 -3.18677061204386616e-2, -6.25840678512839046e-3, -3.01565363306955866e-4, -2.32324863642097035e-8 },
239 { -0.75968322407189071, -0.43093965318039579, -0.23108984359927232, -0.1108289933162471, -0.043213911456684129, -0.011258183689471605, -0.0011859331251521767, -7.6875625812546846e-6 }
241 return ArrayVector<double>(coeffs[ORDER], coeffs[ORDER]+ORDER/2);
244 template <
int ORDER,
class T>
245 typename BSplineBase<ORDER, T>::WeightMatrix
246 BSplineBase<ORDER, T>::calculateWeightMatrix()
248 WeightMatrix res(ORDER+1, ArrayVector<T>(ORDER+1));
249 double faculty = 1.0;
250 for(
int d = 0; d <= ORDER; ++d)
254 double x = ORDER / 2;
256 for(
int i = 0; i <= ORDER; ++i, --x)
257 res[d][i] = spline(x, d) / faculty;
273 template <
int ORDER,
class T =
double>
292 class BSplineBase<0, T>
309 return exec(x, derivativeOrder_);
312 template <
unsigned int IntBits,
unsigned int FracBits>
313 FixedPoint<IntBits, FracBits>
operator()(FixedPoint<IntBits, FracBits> x)
const
315 typedef FixedPoint<IntBits, FracBits> Value;
316 return x.value < Value::ONE_HALF && -Value::ONE_HALF <= x.value
317 ? Value(Value::ONE, FPNoShift)
318 : Value(0, FPNoShift);
321 template <
class U,
int N>
322 autodiff::DualVector<U, N>
operator()(autodiff::DualVector<U, N>
const & x)
const
324 return x < 0.5 && -0.5 <= x
325 ? autodiff::DualVector<U, N>(1.0)
326 : autodiff::DualVector<U, N>(0.0);
331 return exec(x, derivativeOrder_ + derivative_order);
341 {
return derivativeOrder_; }
345 return prefilterCoefficients_;
348 typedef T WeightMatrix[1][1];
350 static WeightMatrix
const &
weights()
352 return weightMatrix_;
358 if(derivative_order == 0)
359 return x < 0.5 && -0.5 <= x ?
366 unsigned int derivativeOrder_;
367 static ArrayVector<double> prefilterCoefficients_;
368 static WeightMatrix weightMatrix_;
372 ArrayVector<double> BSplineBase<0, T>::prefilterCoefficients_;
375 typename BSplineBase<0, T>::WeightMatrix BSplineBase<0, T>::weightMatrix_ = {{ 1.0 }};
401 return exec(x, derivativeOrder_);
404 template <
unsigned int IntBits,
unsigned int FracBits>
405 FixedPoint<IntBits, FracBits>
operator()(FixedPoint<IntBits, FracBits> x)
const
407 typedef FixedPoint<IntBits, FracBits> Value;
408 int v =
abs(x.value);
409 return v < Value::ONE ?
410 Value(Value::ONE - v, FPNoShift)
411 : Value(0, FPNoShift);
414 template <
class U,
int N>
415 autodiff::DualVector<U, N>
operator()(autodiff::DualVector<U, N> x)
const
420 : autodiff::DualVector<U, N>(0.0);
425 return exec(x, derivativeOrder_ + derivative_order);
435 {
return derivativeOrder_; }
439 return prefilterCoefficients_;
442 typedef T WeightMatrix[2][2];
444 static WeightMatrix
const &
weights()
446 return weightMatrix_;
450 T exec(T x,
unsigned int derivative_order)
const;
452 unsigned int derivativeOrder_;
453 static ArrayVector<double> prefilterCoefficients_;
454 static WeightMatrix weightMatrix_;
458 ArrayVector<double> BSpline<1, T>::prefilterCoefficients_;
461 typename BSpline<1, T>::WeightMatrix BSpline<1, T>::weightMatrix_ = {{ 1.0, 0.0}, {-1.0, 1.0}};
464 T BSpline<1, T>::exec(T x,
unsigned int derivative_order)
const
466 switch(derivative_order)
470 x = VIGRA_CSTD::fabs(x);
514 return exec(x, derivativeOrder_);
517 template <
unsigned int IntBits,
unsigned int FracBits>
518 FixedPoint<IntBits, FracBits>
operator()(FixedPoint<IntBits, FracBits> x)
const
520 typedef FixedPoint<IntBits, FracBits> Value;
521 enum { ONE_HALF = Value::ONE_HALF, THREE_HALVES = ONE_HALF * 3, THREE_QUARTERS = THREE_HALVES / 2,
522 PREMULTIPLY_SHIFT1 = FracBits <= 16 ? 0 : FracBits - 16,
523 PREMULTIPLY_SHIFT2 = FracBits - 1 <= 16 ? 0 : FracBits - 17,
524 POSTMULTIPLY_SHIFT1 = FracBits - 2*PREMULTIPLY_SHIFT1,
525 POSTMULTIPLY_SHIFT2 = FracBits - 2*PREMULTIPLY_SHIFT2 };
526 int v =
abs(x.value);
528 ? Value(ONE_HALF, FPNoShift)
530 ? Value(THREE_QUARTERS -
531 (int)(
sq((unsigned)v >> PREMULTIPLY_SHIFT2) >> POSTMULTIPLY_SHIFT2), FPNoShift)
533 ? Value((int)(
sq((unsigned)(THREE_HALVES-v) >> PREMULTIPLY_SHIFT1) >> (POSTMULTIPLY_SHIFT1 + 1)), FPNoShift)
534 : Value(0, FPNoShift);
537 template <
class U,
int N>
538 autodiff::DualVector<U, N>
operator()(autodiff::DualVector<U, N> x)
const
545 : autodiff::DualVector<U, N>(0.0);
550 return exec(x, derivativeOrder_ + derivative_order);
563 {
return derivativeOrder_; }
567 return prefilterCoefficients_;
570 typedef T WeightMatrix[3][3];
572 static WeightMatrix
const &
weights()
574 return weightMatrix_;
580 unsigned int derivativeOrder_;
581 static ArrayVector<double> prefilterCoefficients_;
582 static WeightMatrix weightMatrix_;
586 ArrayVector<double> BSpline<2, T>::prefilterCoefficients_(1, 2.0*M_SQRT2 - 3.0);
589 typename BSpline<2, T>::WeightMatrix BSpline<2, T>::weightMatrix_ =
590 {{ 0.125, 0.75, 0.125},
595 typename BSpline<2, T>::result_type
596 BSpline<2, T>::exec(first_argument_type x, second_argument_type derivative_order)
const
598 switch(derivative_order)
602 x = VIGRA_CSTD::fabs(x);
662 return exec(x, derivativeOrder_);
665 template <
unsigned int IntBits,
unsigned int FracBits>
666 FixedPoint<IntBits, FracBits>
operator()(FixedPoint<IntBits, FracBits> x)
const
668 typedef FixedPoint<IntBits, FracBits> Value;
669 enum { ONE = Value::ONE, TWO = 2 * ONE, TWO_THIRDS = TWO / 3, ONE_SIXTH = ONE / 6,
670 PREMULTIPLY_SHIFT = FracBits <= 16 ? 0 : FracBits - 16,
671 POSTMULTIPLY_SHIFT = FracBits - 2*PREMULTIPLY_SHIFT };
672 int v =
abs(x.value);
674 ? Value(ONE_SIXTH, FPNoShift)
677 (((int)(
sq((unsigned)v >> PREMULTIPLY_SHIFT) >> (POSTMULTIPLY_SHIFT + PREMULTIPLY_SHIFT))
678 * (((v >> 1) - ONE) >> PREMULTIPLY_SHIFT)) >> POSTMULTIPLY_SHIFT), FPNoShift)
680 ? Value((int)((
sq((unsigned)(TWO-v) >> PREMULTIPLY_SHIFT) >> (POSTMULTIPLY_SHIFT + PREMULTIPLY_SHIFT))
681 * ((unsigned)(TWO-v) >> PREMULTIPLY_SHIFT) / 6) >> POSTMULTIPLY_SHIFT, FPNoShift)
682 : Value(0, FPNoShift);
685 template <
class U,
int N>
686 autodiff::DualVector<U, N>
operator()(autodiff::DualVector<U, N> x)
const
691 return 2.0/3.0 + x*x*(-1.0 + 0.5*x);
699 return autodiff::DualVector<U, N>(0.0);
704 return exec(x, derivativeOrder_ + derivative_order);
720 {
return derivativeOrder_; }
724 return prefilterCoefficients_;
727 typedef T WeightMatrix[4][4];
729 static WeightMatrix
const &
weights()
731 return weightMatrix_;
737 unsigned int derivativeOrder_;
738 static ArrayVector<double> prefilterCoefficients_;
739 static WeightMatrix weightMatrix_;
743 ArrayVector<double> BSpline<3, T>::prefilterCoefficients_(1,
VIGRA_CSTD::sqrt(3.0) - 2.0);
746 typename BSpline<3, T>::WeightMatrix BSpline<3, T>::weightMatrix_ =
747 {{ 1.0 / 6.0, 2.0 / 3.0, 1.0 / 6.0, 0.0},
748 {-0.5, 0.0, 0.5, 0.0},
749 { 0.5, -1.0, 0.5, 0.0},
750 {-1.0 / 6.0, 0.5, -0.5, 1.0 / 6.0}};
753 typename BSpline<3, T>::result_type
754 BSpline<3, T>::exec(first_argument_type x, second_argument_type derivative_order)
const
756 switch(derivative_order)
760 x = VIGRA_CSTD::fabs(x);
763 return 2.0/3.0 + x*x*(-1.0 + 0.5*x);
778 x = VIGRA_CSTD::fabs(x);
787 x = VIGRA_CSTD::fabs(x);
813 typedef BSpline<3, double> CubicBSplineKernel;
839 return exec(x, derivativeOrder_);
844 return exec(x, derivativeOrder_ + derivative_order);
847 template <
class U,
int N>
848 autodiff::DualVector<U, N>
operator()(autodiff::DualVector<U, N> x)
const
853 return 115.0/192.0 + x*x*(-0.625 + x*x*0.25);
857 return (55.0/16.0 + x*(1.25 + x*(-7.5 + x*(5.0 - x)))) / 6.0;
862 return sq(x*x) / 24.0;
865 return autodiff::DualVector<U, N>(0.0);
884 {
return derivativeOrder_; }
888 return prefilterCoefficients_;
891 typedef T WeightMatrix[5][5];
893 static WeightMatrix
const &
weights()
895 return weightMatrix_;
901 unsigned int derivativeOrder_;
902 static ArrayVector<double> prefilterCoefficients_;
903 static WeightMatrix weightMatrix_;
907 ArrayVector<double> BSpline<4, T>::prefilterCoefficients_(BSplineBase<4, T>::getPrefilterCoefficients());
910 typename BSpline<4, T>::WeightMatrix BSpline<4, T>::weightMatrix_ =
911 {{ 1.0/384.0, 19.0/96.0, 115.0/192.0, 19.0/96.0, 1.0/384.0},
912 {-1.0/48.0, -11.0/24.0, 0.0, 11.0/24.0, 1.0/48.0},
913 { 1.0/16.0, 1.0/4.0, -5.0/8.0, 1.0/4.0, 1.0/16.0},
914 {-1.0/12.0, 1.0/6.0, 0.0, -1.0/6.0, 1.0/12.0},
915 { 1.0/24.0, -1.0/6.0, 0.25, -1.0/6.0, 1.0/24.0}};
918 typename BSpline<4, T>::result_type
919 BSpline<4, T>::exec(first_argument_type x, second_argument_type derivative_order)
const
921 switch(derivative_order)
925 x = VIGRA_CSTD::fabs(x);
928 return 115.0/192.0 + x*x*(-0.625 + x*x*0.25);
932 return (55.0/16.0 + x*(1.25 + x*(-7.5 + x*(5.0 - x)))) / 6.0;
937 return sq(x*x) / 24.0;
947 x = VIGRA_CSTD::fabs(x);
950 return s*x*(-1.25 + x*x);
954 return s*(5.0 + x*(-60.0 + x*(60.0 - 16.0*x))) / 24.0;
959 return s*x*x*x / -6.0;
966 x = VIGRA_CSTD::fabs(x);
969 return -1.25 + 3.0*x*x;
973 return -2.5 + x*(5.0 - 2.0*x);
988 x = VIGRA_CSTD::fabs(x);
995 return s*(5.0 - 4.0*x);
1027 typedef BSpline<4, double> QuarticBSplineKernel;
1053 return exec(x, derivativeOrder_);
1058 return exec(x, derivativeOrder_ + derivative_order);
1061 template <
class U,
int N>
1062 autodiff::DualVector<U, N>
operator()(autodiff::DualVector<U, N> x)
const
1067 return 0.55 + x*x*(-0.5 + x*x*(0.25 - x/12.0));
1071 return 17.0/40.0 + x*(0.625 + x*(-1.75 + x*(1.25 + x*(-0.375 + x/24.0))));
1076 return x*
sq(x*x) / 120.0;
1079 return autodiff::DualVector<U, N>(0.0);
1101 {
return derivativeOrder_; }
1105 return prefilterCoefficients_;
1108 typedef T WeightMatrix[6][6];
1110 static WeightMatrix
const &
weights()
1112 return weightMatrix_;
1118 unsigned int derivativeOrder_;
1119 static ArrayVector<double> prefilterCoefficients_;
1120 static WeightMatrix weightMatrix_;
1124 ArrayVector<double> BSpline<5, T>::prefilterCoefficients_(BSplineBase<5, T>::getPrefilterCoefficients());
1127 typename BSpline<5, T>::WeightMatrix BSpline<5, T>::weightMatrix_ =
1128 {{ 1.0/120.0, 13.0/60.0, 11.0/20.0, 13.0/60.0, 1.0/120.0, 0.0},
1129 {-1.0/24.0, -5.0/12.0, 0.0, 5.0/12.0, 1.0/24.0, 0.0},
1130 { 1.0/12.0, 1.0/6.0, -0.5, 1.0/6.0, 1.0/12.0, 0.0},
1131 {-1.0/12.0, 1.0/6.0, 0.0, -1.0/6.0, 1.0/12.0, 0.0},
1132 { 1.0/24.0, -1.0/6.0, 0.25, -1.0/6.0, 1.0/24.0, 0.0},
1133 {-1.0/120.0, 1.0/24.0, -1.0/12.0, 1.0/12.0, -1.0/24.0, 1.0/120.0}};
1136 typename BSpline<5, T>::result_type
1137 BSpline<5, T>::exec(first_argument_type x, second_argument_type derivative_order)
const
1139 switch(derivative_order)
1143 x = VIGRA_CSTD::fabs(x);
1146 return 0.55 + x*x*(-0.5 + x*x*(0.25 - x/12.0));
1150 return 17.0/40.0 + x*(0.625 + x*(-1.75 + x*(1.25 + x*(-0.375 + x/24.0))));
1155 return x*
sq(x*x) / 120.0;
1162 double s = x < 0.0 ?
1165 x = VIGRA_CSTD::fabs(x);
1168 return s*x*(-1.0 + x*x*(1.0 - 5.0/12.0*x));
1172 return s*(0.625 + x*(-3.5 + x*(3.75 + x*(-1.5 + 5.0/24.0*x))));
1177 return s*
sq(x*x) / -24.0;
1184 x = VIGRA_CSTD::fabs(x);
1187 return -1.0 + x*x*(3.0 -5.0/3.0*x);
1191 return -3.5 + x*(7.5 + x*(-4.5 + 5.0/6.0*x));
1203 double s = x < 0.0 ?
1206 x = VIGRA_CSTD::fabs(x);
1209 return s*x*(6.0 - 5.0*x);
1213 return s*(7.5 + x*(-9.0 + 2.5*x));
1225 x = VIGRA_CSTD::fabs(x);
1228 return 6.0 - 10.0*x;
1232 return -9.0 + 5.0*x;
1264 typedef BSpline<5, double> QuinticBSplineKernel;
1266 #endif // NO_PARTIAL_TEMPLATE_SPECIALIZATION
1294 template <
class T =
double>
1336 return prefilterCoefficients_;
1344 ArrayVector<double> CatmullRomSpline<T>::prefilterCoefficients_;
1347 typename CatmullRomSpline<T>::result_type
1350 x = VIGRA_CSTD::fabs(x);
1353 return 1.0 + x * x * (-2.5 + 1.5 * x);
1361 return 2.0 + x * (-4.0 + x * (2.5 - 0.5 * x));
T operator[](T x) const
Definition: splines.hxx:1317
StaticOrder
Definition: splines.hxx:1309
double radius() const
Definition: splines.hxx:149
value_type operator[](value_type x) const
Definition: splines.hxx:143
T result_type
Definition: splines.hxx:1306
T argument_type
Definition: splines.hxx:1303
ArrayVector< double > const & prefilterCoefficients() const
Definition: splines.hxx:165
static WeightMatrix const & weights()
Definition: splines.hxx:177
Definition: splines.hxx:274
StaticOrder
Definition: splines.hxx:112
Definition: splines.hxx:89
T argument_type
Definition: splines.hxx:100
T first_argument_type
Definition: splines.hxx:103
NumericTraits< T >::Promote sq(T t)
The square function.
Definition: mathutil.hxx:382
result_type operator()(argument_type x) const
Definition: splines.hxx:126
int radius() const
Definition: splines.hxx:1323
unsigned int derivativeOrder() const
Definition: splines.hxx:1328
T value_type
Definition: splines.hxx:1300
unsigned int derivativeOrder() const
Definition: splines.hxx:154
BSplineBase(unsigned int derivativeOrder=0)
Definition: splines.hxx:117
unsigned int second_argument_type
Definition: splines.hxx:106
ArrayVector< double > const & prefilterCoefficients() const
Definition: splines.hxx:1334
T result_type
Definition: splines.hxx:109
FFTWComplex< R >::NormType abs(const FFTWComplex< R > &a)
absolute value (= magnitude)
Definition: fftw3.hxx:1002
BSpline(unsigned int derivativeOrder=0)
Definition: splines.hxx:280
result_type operator()(first_argument_type x, second_argument_type derivative_order) const
Definition: splines.hxx:136
T value_type
Definition: splines.hxx:93
result_type operator()(argument_type x) const
Definition: splines.hxx:1348
SquareRootTraits< FixedPoint< IntBits, FracBits > >::SquareRootResult sqrt(FixedPoint< IntBits, FracBits > v)
square root.
Definition: fixedpoint.hxx:616
Definition: splines.hxx:1295