[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]

splines.hxx VIGRA

1 /************************************************************************/
2 /* */
3 /* Copyright 1998-2004 by Ullrich Koethe */
4 /* */
5 /* This file is part of the VIGRA computer vision library. */
6 /* The VIGRA Website is */
7 /* http://hci.iwr.uni-heidelberg.de/vigra/ */
8 /* Please direct questions, bug reports, and contributions to */
9 /* ullrich.koethe@iwr.uni-heidelberg.de or */
10 /* vigra@informatik.uni-hamburg.de */
11 /* */
12 /* Permission is hereby granted, free of charge, to any person */
13 /* obtaining a copy of this software and associated documentation */
14 /* files (the "Software"), to deal in the Software without */
15 /* restriction, including without limitation the rights to use, */
16 /* copy, modify, merge, publish, distribute, sublicense, and/or */
17 /* sell copies of the Software, and to permit persons to whom the */
18 /* Software is furnished to do so, subject to the following */
19 /* conditions: */
20 /* */
21 /* The above copyright notice and this permission notice shall be */
22 /* included in all copies or substantial portions of the */
23 /* Software. */
24 /* */
25 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */
26 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */
27 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */
28 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */
29 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */
30 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */
31 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */
32 /* OTHER DEALINGS IN THE SOFTWARE. */
33 /* */
34 /************************************************************************/
35 
36 #ifndef VIGRA_SPLINES_HXX
37 #define VIGRA_SPLINES_HXX
38 
39 #include <cmath>
40 #include "config.hxx"
41 #include "mathutil.hxx"
42 #include "array_vector.hxx"
43 #include "fixedpoint.hxx"
44 
45 namespace vigra {
46 
47 namespace autodiff {
48 
49 template <class T, int N>
50 class DualVector;
51 
52 } // namespace autodiff
53 
54 /** \addtogroup MathFunctions Mathematical Functions
55 */
56 //@{
57 /* B-Splines of arbitrary order and interpolating Catmull/Rom splines.
58 
59  <b>\#include</b> <vigra/splines.hxx><br>
60  Namespace: vigra
61 */
62 #ifndef NO_PARTIAL_TEMPLATE_SPECIALIZATION
63 
64 /** Basic interface of the spline functors.
65 
66  Implements the spline functions defined by the recursion
67 
68  \f[ B_0(x) = \left\{ \begin{array}{ll}
69  1 & -\frac{1}{2} \leq x < \frac{1}{2} \\
70  0 & \mbox{otherwise}
71  \end{array}\right.
72  \f]
73 
74  and
75 
76  \f[ B_n(x) = B_0(x) * B_{n-1}(x)
77  \f]
78 
79  where * denotes convolution, and <i>n</i> is the spline order given by the template
80  parameter <tt>ORDER</tt> with <tt>ORDER < 18</tt>. These spline classes can be used as
81  unary and binary functors, as kernels for \ref resamplingConvolveImage(),
82  and as arguments for \ref vigra::SplineImageView. Note that the spline order
83  is given as a template argument.
84 
85  <b>\#include</b> <vigra/splines.hxx><br>
86  Namespace: vigra
87 */
88 template <int ORDER, class T = double>
90 {
91  public:
92 
93  static_assert (ORDER < 18 , "BSpline: ORDER must be less than 18." );
94 
95  /** the value type if used as a kernel in \ref resamplingConvolveImage().
96  */
97  typedef T value_type;
98  /** the functor's unary argument type
99  */
100  typedef T argument_type;
101  /** the functor's first binary argument type
102  */
104  /** the functor's second binary argument type
105  */
106  typedef unsigned int second_argument_type;
107  /** the functor's result type (unary and binary)
108  */
109  typedef T result_type;
110  /** the spline order
111  */
112  enum StaticOrder { order = ORDER };
113 
114  /** Create functor for given derivative of the spline. The spline's order
115  is specified spline by the template argument <TT>ORDER</tt>.
116  */
117  explicit BSplineBase(unsigned int derivativeOrder = 0)
118  : s1_(derivativeOrder)
119  {}
120 
121  /** Unary function call.
122  Returns the value of the spline with the derivative order given in the
123  constructor. Note that only derivatives up to <tt>ORDER-1</tt> are
124  continuous, and derivatives above <tt>ORDER+1</tt> are zero.
125  */
127  {
128  return exec(x, derivativeOrder());
129  }
130 
131  /** Binary function call.
132  The given derivative order is added to the derivative order
133  specified in the constructor. Note that only derivatives up to <tt>ORDER-1</tt> are
134  continuous, and derivatives above <tt>ORDER+1</tt> are zero.
135  */
137  {
138  return exec(x, derivativeOrder() + derivative_order);
139  }
140 
141  /** Index operator. Same as unary function call.
142  */
144  { return operator()(x); }
145 
146  /** Get the required filter radius for a discrete approximation of the
147  spline. Always equal to <tt>(ORDER + 1) / 2.0</tt>.
148  */
149  double radius() const
150  { return (ORDER + 1) * 0.5; }
151 
152  /** Get the derivative order of the Gaussian.
153  */
154  unsigned int derivativeOrder() const
155  { return s1_.derivativeOrder(); }
156 
157  /** Get the prefilter coefficients required for interpolation.
158  To interpolate with a B-spline, \ref resamplingConvolveImage()
159  can be used. However, the image to be interpolated must be
160  pre-filtered using \ref recursiveFilterX() and \ref recursiveFilterY()
161  with the filter coefficients given by this function. The length of the array
162  corresponds to how many times the above recursive filtering
163  has to be applied (zero length means no prefiltering necessary).
164  */
166  {
167  return prefilterCoefficients_;
168  }
169 
170  typedef ArrayVector<ArrayVector<T> > WeightMatrix;
171 
172  /** Get the coefficients to transform spline coefficients into
173  the coefficients of the corresponding polynomial.
174  Currently internally used in SplineImageView; needs more
175  documentation ???
176  */
177  static WeightMatrix const & weights()
178  {
179  return weightMatrix_;
180  }
181 
182  static ArrayVector<double> getPrefilterCoefficients();
183 
184  protected:
185  result_type exec(first_argument_type x, second_argument_type derivative_order) const;
186 
187  // factory function for the weight matrix
188  static WeightMatrix calculateWeightMatrix();
189 
190  BSplineBase<ORDER-1, T> s1_;
191  static ArrayVector<double> prefilterCoefficients_;
192  static WeightMatrix weightMatrix_;
193 };
194 
195 template <int ORDER, class T>
196 ArrayVector<double> BSplineBase<ORDER, T>::prefilterCoefficients_(getPrefilterCoefficients());
197 
198 template <int ORDER, class T>
199 typename BSplineBase<ORDER, T>::WeightMatrix BSplineBase<ORDER, T>::weightMatrix_(calculateWeightMatrix());
200 
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
204 {
205  if(derivative_order == 0)
206  {
207  T n12 = (ORDER + 1.0) / 2.0;
208  return ((n12 + x) * s1_(x + 0.5) + (n12 - x) * s1_(x - 0.5)) / ORDER;
209  }
210  else
211  {
212  --derivative_order;
213  return s1_(x + 0.5, derivative_order) - s1_(x - 0.5, derivative_order);
214  }
215 }
216 
217 template <int ORDER, class T>
218 ArrayVector<double>
219 BSplineBase<ORDER, T>::getPrefilterCoefficients()
220 {
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 }
240  };
241  return ArrayVector<double>(coeffs[ORDER], coeffs[ORDER]+ORDER/2);
242 }
243 
244 template <int ORDER, class T>
245 typename BSplineBase<ORDER, T>::WeightMatrix
246 BSplineBase<ORDER, T>::calculateWeightMatrix()
247 {
248  WeightMatrix res(ORDER+1, ArrayVector<T>(ORDER+1));
249  double faculty = 1.0;
250  for(int d = 0; d <= ORDER; ++d)
251  {
252  if(d > 1)
253  faculty *= d;
254  double x = ORDER / 2; // (note: integer division)
255  BSplineBase spline;
256  for(int i = 0; i <= ORDER; ++i, --x)
257  res[d][i] = spline(x, d) / faculty;
258  }
259  return res;
260 }
261 
262 /********************************************************/
263 /* */
264 /* BSpline<N, T> */
265 /* */
266 /********************************************************/
267 
268 /** Spline functors for arbitrary orders.
269 
270  Provides the interface of \ref vigra::BSplineBase with a more convenient
271  name -- see there for more documentation.
272 */
273 template <int ORDER, class T = double>
274 class BSpline
275 : public BSplineBase<ORDER, T>
276 {
277  public:
278  /** Constructor forwarded to the base class constructor..
279  */
280  explicit BSpline(unsigned int derivativeOrder = 0)
281  : BSplineBase<ORDER, T>(derivativeOrder)
282  {}
283 };
284 
285 /********************************************************/
286 /* */
287 /* BSpline<0, T> */
288 /* */
289 /********************************************************/
290 
291 template <class T>
292 class BSplineBase<0, T>
293 {
294  public:
295 
296  typedef T value_type;
297  typedef T argument_type;
298  typedef T first_argument_type;
299  typedef unsigned int second_argument_type;
300  typedef T result_type;
301  enum StaticOrder { order = 0 };
302 
303  explicit BSplineBase(unsigned int derivativeOrder = 0)
304  : derivativeOrder_(derivativeOrder)
305  {}
306 
308  {
309  return exec(x, derivativeOrder_);
310  }
311 
312  template <unsigned int IntBits, unsigned int FracBits>
313  FixedPoint<IntBits, FracBits> operator()(FixedPoint<IntBits, FracBits> x) const
314  {
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);
319  }
320 
321  template <class U, int N>
322  autodiff::DualVector<U, N> operator()(autodiff::DualVector<U, N> const & x) const
323  {
324  return x < 0.5 && -0.5 <= x
325  ? autodiff::DualVector<U, N>(1.0)
326  : autodiff::DualVector<U, N>(0.0);
327  }
328 
330  {
331  return exec(x, derivativeOrder_ + derivative_order);
332  }
333 
335  { return operator()(x); }
336 
337  double radius() const
338  { return 0.5; }
339 
340  unsigned int derivativeOrder() const
341  { return derivativeOrder_; }
342 
343  ArrayVector<double> const & prefilterCoefficients() const
344  {
345  return prefilterCoefficients_;
346  }
347 
348  typedef T WeightMatrix[1][1];
349 
350  static WeightMatrix const & weights()
351  {
352  return weightMatrix_;
353  }
354 
355  protected:
356  result_type exec(first_argument_type x, second_argument_type derivative_order) const
357  {
358  if(derivative_order == 0)
359  return x < 0.5 && -0.5 <= x ?
360  1.0
361  : 0.0;
362  else
363  return 0.0;
364  }
365 
366  unsigned int derivativeOrder_;
367  static ArrayVector<double> prefilterCoefficients_;
368  static WeightMatrix weightMatrix_;
369 };
370 
371 template <class T>
372 ArrayVector<double> BSplineBase<0, T>::prefilterCoefficients_;
373 
374 template <class T>
375 typename BSplineBase<0, T>::WeightMatrix BSplineBase<0, T>::weightMatrix_ = {{ 1.0 }};
376 
377 /********************************************************/
378 /* */
379 /* BSpline<1, T> */
380 /* */
381 /********************************************************/
382 
383 template <class T>
384 class BSpline<1, T>
385 {
386  public:
387 
388  typedef T value_type;
389  typedef T argument_type;
390  typedef T first_argument_type;
391  typedef unsigned int second_argument_type;
392  typedef T result_type;
393  enum StaticOrder { order = 1 };
394 
395  explicit BSpline(unsigned int derivativeOrder = 0)
396  : derivativeOrder_(derivativeOrder)
397  {}
398 
400  {
401  return exec(x, derivativeOrder_);
402  }
403 
404  template <unsigned int IntBits, unsigned int FracBits>
405  FixedPoint<IntBits, FracBits> operator()(FixedPoint<IntBits, FracBits> x) const
406  {
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);
412  }
413 
414  template <class U, int N>
415  autodiff::DualVector<U, N> operator()(autodiff::DualVector<U, N> x) const
416  {
417  x = abs(x);
418  return x < 1.0
419  ? 1.0 - x
420  : autodiff::DualVector<U, N>(0.0);
421  }
422 
424  {
425  return exec(x, derivativeOrder_ + derivative_order);
426  }
427 
429  { return operator()(x); }
430 
431  double radius() const
432  { return 1.0; }
433 
434  unsigned int derivativeOrder() const
435  { return derivativeOrder_; }
436 
437  ArrayVector<double> const & prefilterCoefficients() const
438  {
439  return prefilterCoefficients_;
440  }
441 
442  typedef T WeightMatrix[2][2];
443 
444  static WeightMatrix const & weights()
445  {
446  return weightMatrix_;
447  }
448 
449  protected:
450  T exec(T x, unsigned int derivative_order) const;
451 
452  unsigned int derivativeOrder_;
453  static ArrayVector<double> prefilterCoefficients_;
454  static WeightMatrix weightMatrix_;
455 };
456 
457 template <class T>
458 ArrayVector<double> BSpline<1, T>::prefilterCoefficients_;
459 
460 template <class T>
461 typename BSpline<1, T>::WeightMatrix BSpline<1, T>::weightMatrix_ = {{ 1.0, 0.0}, {-1.0, 1.0}};
462 
463 template <class T>
464 T BSpline<1, T>::exec(T x, unsigned int derivative_order) const
465 {
466  switch(derivative_order)
467  {
468  case 0:
469  {
470  x = VIGRA_CSTD::fabs(x);
471  return x < 1.0 ?
472  1.0 - x
473  : 0.0;
474  }
475  case 1:
476  {
477  return x < 0.0 ?
478  -1.0 <= x ?
479  1.0
480  : 0.0
481  : x < 1.0 ?
482  -1.0
483  : 0.0;
484  }
485  default:
486  return 0.0;
487  }
488 }
489 
490 /********************************************************/
491 /* */
492 /* BSpline<2, T> */
493 /* */
494 /********************************************************/
495 
496 template <class T>
497 class BSpline<2, T>
498 {
499  public:
500 
501  typedef T value_type;
502  typedef T argument_type;
503  typedef T first_argument_type;
504  typedef unsigned int second_argument_type;
505  typedef T result_type;
506  enum StaticOrder { order = 2 };
507 
508  explicit BSpline(unsigned int derivativeOrder = 0)
509  : derivativeOrder_(derivativeOrder)
510  {}
511 
513  {
514  return exec(x, derivativeOrder_);
515  }
516 
517  template <unsigned int IntBits, unsigned int FracBits>
518  FixedPoint<IntBits, FracBits> operator()(FixedPoint<IntBits, FracBits> x) const
519  {
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);
527  return v == ONE_HALF
528  ? Value(ONE_HALF, FPNoShift)
529  : v <= ONE_HALF
530  ? Value(THREE_QUARTERS -
531  (int)(sq((unsigned)v >> PREMULTIPLY_SHIFT2) >> POSTMULTIPLY_SHIFT2), FPNoShift)
532  : v < THREE_HALVES
533  ? Value((int)(sq((unsigned)(THREE_HALVES-v) >> PREMULTIPLY_SHIFT1) >> (POSTMULTIPLY_SHIFT1 + 1)), FPNoShift)
534  : Value(0, FPNoShift);
535  }
536 
537  template <class U, int N>
538  autodiff::DualVector<U, N> operator()(autodiff::DualVector<U, N> x) const
539  {
540  x = abs(x);
541  return x < 0.5
542  ? 0.75 - x*x
543  : x < 1.5
544  ? 0.5 * sq(1.5 - x)
545  : autodiff::DualVector<U, N>(0.0);
546  }
547 
549  {
550  return exec(x, derivativeOrder_ + derivative_order);
551  }
552 
553  result_type dx(argument_type x) const
554  { return operator()(x, 1); }
555 
557  { return operator()(x); }
558 
559  double radius() const
560  { return 1.5; }
561 
562  unsigned int derivativeOrder() const
563  { return derivativeOrder_; }
564 
565  ArrayVector<double> const & prefilterCoefficients() const
566  {
567  return prefilterCoefficients_;
568  }
569 
570  typedef T WeightMatrix[3][3];
571 
572  static WeightMatrix const & weights()
573  {
574  return weightMatrix_;
575  }
576 
577  protected:
578  result_type exec(first_argument_type x, second_argument_type derivative_order) const;
579 
580  unsigned int derivativeOrder_;
581  static ArrayVector<double> prefilterCoefficients_;
582  static WeightMatrix weightMatrix_;
583 };
584 
585 template <class T>
586 ArrayVector<double> BSpline<2, T>::prefilterCoefficients_(1, 2.0*M_SQRT2 - 3.0);
587 
588 template <class T>
589 typename BSpline<2, T>::WeightMatrix BSpline<2, T>::weightMatrix_ =
590  {{ 0.125, 0.75, 0.125},
591  {-0.5, 0.0, 0.5},
592  { 0.5, -1.0, 0.5}};
593 
594 template <class T>
595 typename BSpline<2, T>::result_type
596 BSpline<2, T>::exec(first_argument_type x, second_argument_type derivative_order) const
597 {
598  switch(derivative_order)
599  {
600  case 0:
601  {
602  x = VIGRA_CSTD::fabs(x);
603  return x < 0.5 ?
604  0.75 - x*x
605  : x < 1.5 ?
606  0.5 * sq(1.5 - x)
607  : 0.0;
608  }
609  case 1:
610  {
611  return x >= -0.5 ?
612  x <= 0.5 ?
613  -2.0 * x
614  : x < 1.5 ?
615  x - 1.5
616  : 0.0
617  : x > -1.5 ?
618  x + 1.5
619  : 0.0;
620  }
621  case 2:
622  {
623  return x >= -0.5 ?
624  x < 0.5 ?
625  -2.0
626  : x < 1.5 ?
627  1.0
628  : 0.0
629  : x >= -1.5 ?
630  1.0
631  : 0.0;
632  }
633  default:
634  return 0.0;
635  }
636 }
637 
638 /********************************************************/
639 /* */
640 /* BSpline<3, T> */
641 /* */
642 /********************************************************/
643 
644 template <class T>
645 class BSpline<3, T>
646 {
647  public:
648 
649  typedef T value_type;
650  typedef T argument_type;
651  typedef T first_argument_type;
652  typedef unsigned int second_argument_type;
653  typedef T result_type;
654  enum StaticOrder { order = 3 };
655 
656  explicit BSpline(unsigned int derivativeOrder = 0)
657  : derivativeOrder_(derivativeOrder)
658  {}
659 
661  {
662  return exec(x, derivativeOrder_);
663  }
664 
665  template <unsigned int IntBits, unsigned int FracBits>
666  FixedPoint<IntBits, FracBits> operator()(FixedPoint<IntBits, FracBits> x) const
667  {
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);
673  return v == ONE
674  ? Value(ONE_SIXTH, FPNoShift)
675  : v < ONE
676  ? Value(TWO_THIRDS +
677  (((int)(sq((unsigned)v >> PREMULTIPLY_SHIFT) >> (POSTMULTIPLY_SHIFT + PREMULTIPLY_SHIFT))
678  * (((v >> 1) - ONE) >> PREMULTIPLY_SHIFT)) >> POSTMULTIPLY_SHIFT), FPNoShift)
679  : v < TWO
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);
683  }
684 
685  template <class U, int N>
686  autodiff::DualVector<U, N> operator()(autodiff::DualVector<U, N> x) const
687  {
688  x = abs(x);
689  if(x < 1.0)
690  {
691  return 2.0/3.0 + x*x*(-1.0 + 0.5*x);
692  }
693  else if(x < 2.0)
694  {
695  x = 2.0 - x;
696  return x*x*x/6.0;
697  }
698  else
699  return autodiff::DualVector<U, N>(0.0);
700  }
701 
703  {
704  return exec(x, derivativeOrder_ + derivative_order);
705  }
706 
707  result_type dx(argument_type x) const
708  { return operator()(x, 1); }
709 
710  result_type dxx(argument_type x) const
711  { return operator()(x, 2); }
712 
714  { return operator()(x); }
715 
716  double radius() const
717  { return 2.0; }
718 
719  unsigned int derivativeOrder() const
720  { return derivativeOrder_; }
721 
722  ArrayVector<double> const & prefilterCoefficients() const
723  {
724  return prefilterCoefficients_;
725  }
726 
727  typedef T WeightMatrix[4][4];
728 
729  static WeightMatrix const & weights()
730  {
731  return weightMatrix_;
732  }
733 
734  protected:
735  result_type exec(first_argument_type x, second_argument_type derivative_order) const;
736 
737  unsigned int derivativeOrder_;
738  static ArrayVector<double> prefilterCoefficients_;
739  static WeightMatrix weightMatrix_;
740 };
741 
742 template <class T>
743 ArrayVector<double> BSpline<3, T>::prefilterCoefficients_(1, VIGRA_CSTD::sqrt(3.0) - 2.0);
744 
745 template <class T>
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}};
751 
752 template <class T>
753 typename BSpline<3, T>::result_type
754 BSpline<3, T>::exec(first_argument_type x, second_argument_type derivative_order) const
755 {
756  switch(derivative_order)
757  {
758  case 0:
759  {
760  x = VIGRA_CSTD::fabs(x);
761  if(x < 1.0)
762  {
763  return 2.0/3.0 + x*x*(-1.0 + 0.5*x);
764  }
765  else if(x < 2.0)
766  {
767  x = 2.0 - x;
768  return x*x*x/6.0;
769  }
770  else
771  return 0.0;
772  }
773  case 1:
774  {
775  double s = x < 0.0 ?
776  -1.0
777  : 1.0;
778  x = VIGRA_CSTD::fabs(x);
779  return x < 1.0 ?
780  s*x*(-2.0 + 1.5*x)
781  : x < 2.0 ?
782  -0.5*s*sq(2.0 - x)
783  : 0.0;
784  }
785  case 2:
786  {
787  x = VIGRA_CSTD::fabs(x);
788  return x < 1.0 ?
789  3.0*x - 2.0
790  : x < 2.0 ?
791  2.0 - x
792  : 0.0;
793  }
794  case 3:
795  {
796  return x < 0.0 ?
797  x < -1.0 ?
798  x < -2.0 ?
799  0.0
800  : 1.0
801  : -3.0
802  : x < 1.0 ?
803  3.0
804  : x < 2.0 ?
805  -1.0
806  : 0.0;
807  }
808  default:
809  return 0.0;
810  }
811 }
812 
813 typedef BSpline<3, double> CubicBSplineKernel;
814 
815 /********************************************************/
816 /* */
817 /* BSpline<4, T> */
818 /* */
819 /********************************************************/
820 
821 template <class T>
822 class BSpline<4, T>
823 {
824  public:
825 
826  typedef T value_type;
827  typedef T argument_type;
828  typedef T first_argument_type;
829  typedef unsigned int second_argument_type;
830  typedef T result_type;
831  enum StaticOrder { order = 4 };
832 
833  explicit BSpline(unsigned int derivativeOrder = 0)
834  : derivativeOrder_(derivativeOrder)
835  {}
836 
838  {
839  return exec(x, derivativeOrder_);
840  }
841 
843  {
844  return exec(x, derivativeOrder_ + derivative_order);
845  }
846 
847  template <class U, int N>
848  autodiff::DualVector<U, N> operator()(autodiff::DualVector<U, N> x) const
849  {
850  x = abs(x);
851  if(x <= 0.5)
852  {
853  return 115.0/192.0 + x*x*(-0.625 + x*x*0.25);
854  }
855  else if(x < 1.5)
856  {
857  return (55.0/16.0 + x*(1.25 + x*(-7.5 + x*(5.0 - x)))) / 6.0;
858  }
859  else if(x < 2.5)
860  {
861  x = 2.5 - x;
862  return sq(x*x) / 24.0;
863  }
864  else
865  return autodiff::DualVector<U, N>(0.0);
866  }
867 
868  result_type dx(argument_type x) const
869  { return operator()(x, 1); }
870 
871  result_type dxx(argument_type x) const
872  { return operator()(x, 2); }
873 
874  result_type dx3(argument_type x) const
875  { return operator()(x, 3); }
876 
878  { return operator()(x); }
879 
880  double radius() const
881  { return 2.5; }
882 
883  unsigned int derivativeOrder() const
884  { return derivativeOrder_; }
885 
886  ArrayVector<double> const & prefilterCoefficients() const
887  {
888  return prefilterCoefficients_;
889  }
890 
891  typedef T WeightMatrix[5][5];
892 
893  static WeightMatrix const & weights()
894  {
895  return weightMatrix_;
896  }
897 
898  protected:
899  result_type exec(first_argument_type x, second_argument_type derivative_order) const;
900 
901  unsigned int derivativeOrder_;
902  static ArrayVector<double> prefilterCoefficients_;
903  static WeightMatrix weightMatrix_;
904 };
905 
906 template <class T>
907 ArrayVector<double> BSpline<4, T>::prefilterCoefficients_(BSplineBase<4, T>::getPrefilterCoefficients());
908 
909 template <class T>
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}};
916 
917 template <class T>
918 typename BSpline<4, T>::result_type
919 BSpline<4, T>::exec(first_argument_type x, second_argument_type derivative_order) const
920 {
921  switch(derivative_order)
922  {
923  case 0:
924  {
925  x = VIGRA_CSTD::fabs(x);
926  if(x <= 0.5)
927  {
928  return 115.0/192.0 + x*x*(-0.625 + x*x*0.25);
929  }
930  else if(x < 1.5)
931  {
932  return (55.0/16.0 + x*(1.25 + x*(-7.5 + x*(5.0 - x)))) / 6.0;
933  }
934  else if(x < 2.5)
935  {
936  x = 2.5 - x;
937  return sq(x*x) / 24.0;
938  }
939  else
940  return 0.0;
941  }
942  case 1:
943  {
944  double s = x < 0.0 ?
945  -1.0 :
946  1.0;
947  x = VIGRA_CSTD::fabs(x);
948  if(x <= 0.5)
949  {
950  return s*x*(-1.25 + x*x);
951  }
952  else if(x < 1.5)
953  {
954  return s*(5.0 + x*(-60.0 + x*(60.0 - 16.0*x))) / 24.0;
955  }
956  else if(x < 2.5)
957  {
958  x = 2.5 - x;
959  return s*x*x*x / -6.0;
960  }
961  else
962  return 0.0;
963  }
964  case 2:
965  {
966  x = VIGRA_CSTD::fabs(x);
967  if(x <= 0.5)
968  {
969  return -1.25 + 3.0*x*x;
970  }
971  else if(x < 1.5)
972  {
973  return -2.5 + x*(5.0 - 2.0*x);
974  }
975  else if(x < 2.5)
976  {
977  x = 2.5 - x;
978  return x*x / 2.0;
979  }
980  else
981  return 0.0;
982  }
983  case 3:
984  {
985  double s = x < 0.0 ?
986  -1.0 :
987  1.0;
988  x = VIGRA_CSTD::fabs(x);
989  if(x <= 0.5)
990  {
991  return s*x*6.0;
992  }
993  else if(x < 1.5)
994  {
995  return s*(5.0 - 4.0*x);
996  }
997  else if(x < 2.5)
998  {
999  return s*(x - 2.5);
1000  }
1001  else
1002  return 0.0;
1003  }
1004  case 4:
1005  {
1006  return x < 0.0
1007  ? x < -2.5
1008  ? 0.0
1009  : x < -1.5
1010  ? 1.0
1011  : x < -0.5
1012  ? -4.0
1013  : 6.0
1014  : x < 0.5
1015  ? 6.0
1016  : x < 1.5
1017  ? -4.0
1018  : x < 2.5
1019  ? 1.0
1020  : 0.0;
1021  }
1022  default:
1023  return 0.0;
1024  }
1025 }
1026 
1027 typedef BSpline<4, double> QuarticBSplineKernel;
1028 
1029 /********************************************************/
1030 /* */
1031 /* BSpline<5, T> */
1032 /* */
1033 /********************************************************/
1034 
1035 template <class T>
1036 class BSpline<5, T>
1037 {
1038  public:
1039 
1040  typedef T value_type;
1041  typedef T argument_type;
1042  typedef T first_argument_type;
1043  typedef unsigned int second_argument_type;
1044  typedef T result_type;
1045  enum StaticOrder { order = 5 };
1046 
1047  explicit BSpline(unsigned int derivativeOrder = 0)
1048  : derivativeOrder_(derivativeOrder)
1049  {}
1050 
1052  {
1053  return exec(x, derivativeOrder_);
1054  }
1055 
1057  {
1058  return exec(x, derivativeOrder_ + derivative_order);
1059  }
1060 
1061  template <class U, int N>
1062  autodiff::DualVector<U, N> operator()(autodiff::DualVector<U, N> x) const
1063  {
1064  x = abs(x);
1065  if(x <= 1.0)
1066  {
1067  return 0.55 + x*x*(-0.5 + x*x*(0.25 - x/12.0));
1068  }
1069  else if(x < 2.0)
1070  {
1071  return 17.0/40.0 + x*(0.625 + x*(-1.75 + x*(1.25 + x*(-0.375 + x/24.0))));
1072  }
1073  else if(x < 3.0)
1074  {
1075  x = 3.0 - x;
1076  return x*sq(x*x) / 120.0;
1077  }
1078  else
1079  return autodiff::DualVector<U, N>(0.0);
1080  }
1081 
1082  result_type dx(argument_type x) const
1083  { return operator()(x, 1); }
1084 
1085  result_type dxx(argument_type x) const
1086  { return operator()(x, 2); }
1087 
1088  result_type dx3(argument_type x) const
1089  { return operator()(x, 3); }
1090 
1091  result_type dx4(argument_type x) const
1092  { return operator()(x, 4); }
1093 
1095  { return operator()(x); }
1096 
1097  double radius() const
1098  { return 3.0; }
1099 
1100  unsigned int derivativeOrder() const
1101  { return derivativeOrder_; }
1102 
1103  ArrayVector<double> const & prefilterCoefficients() const
1104  {
1105  return prefilterCoefficients_;
1106  }
1107 
1108  typedef T WeightMatrix[6][6];
1109 
1110  static WeightMatrix const & weights()
1111  {
1112  return weightMatrix_;
1113  }
1114 
1115  protected:
1116  result_type exec(first_argument_type x, second_argument_type derivative_order) const;
1117 
1118  unsigned int derivativeOrder_;
1119  static ArrayVector<double> prefilterCoefficients_;
1120  static WeightMatrix weightMatrix_;
1121 };
1122 
1123 template <class T>
1124 ArrayVector<double> BSpline<5, T>::prefilterCoefficients_(BSplineBase<5, T>::getPrefilterCoefficients());
1125 
1126 template <class T>
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}};
1134 
1135 template <class T>
1136 typename BSpline<5, T>::result_type
1137 BSpline<5, T>::exec(first_argument_type x, second_argument_type derivative_order) const
1138 {
1139  switch(derivative_order)
1140  {
1141  case 0:
1142  {
1143  x = VIGRA_CSTD::fabs(x);
1144  if(x <= 1.0)
1145  {
1146  return 0.55 + x*x*(-0.5 + x*x*(0.25 - x/12.0));
1147  }
1148  else if(x < 2.0)
1149  {
1150  return 17.0/40.0 + x*(0.625 + x*(-1.75 + x*(1.25 + x*(-0.375 + x/24.0))));
1151  }
1152  else if(x < 3.0)
1153  {
1154  x = 3.0 - x;
1155  return x*sq(x*x) / 120.0;
1156  }
1157  else
1158  return 0.0;
1159  }
1160  case 1:
1161  {
1162  double s = x < 0.0 ?
1163  -1.0 :
1164  1.0;
1165  x = VIGRA_CSTD::fabs(x);
1166  if(x <= 1.0)
1167  {
1168  return s*x*(-1.0 + x*x*(1.0 - 5.0/12.0*x));
1169  }
1170  else if(x < 2.0)
1171  {
1172  return s*(0.625 + x*(-3.5 + x*(3.75 + x*(-1.5 + 5.0/24.0*x))));
1173  }
1174  else if(x < 3.0)
1175  {
1176  x = 3.0 - x;
1177  return s*sq(x*x) / -24.0;
1178  }
1179  else
1180  return 0.0;
1181  }
1182  case 2:
1183  {
1184  x = VIGRA_CSTD::fabs(x);
1185  if(x <= 1.0)
1186  {
1187  return -1.0 + x*x*(3.0 -5.0/3.0*x);
1188  }
1189  else if(x < 2.0)
1190  {
1191  return -3.5 + x*(7.5 + x*(-4.5 + 5.0/6.0*x));
1192  }
1193  else if(x < 3.0)
1194  {
1195  x = 3.0 - x;
1196  return x*x*x / 6.0;
1197  }
1198  else
1199  return 0.0;
1200  }
1201  case 3:
1202  {
1203  double s = x < 0.0 ?
1204  -1.0 :
1205  1.0;
1206  x = VIGRA_CSTD::fabs(x);
1207  if(x <= 1.0)
1208  {
1209  return s*x*(6.0 - 5.0*x);
1210  }
1211  else if(x < 2.0)
1212  {
1213  return s*(7.5 + x*(-9.0 + 2.5*x));
1214  }
1215  else if(x < 3.0)
1216  {
1217  x = 3.0 - x;
1218  return -0.5*s*x*x;
1219  }
1220  else
1221  return 0.0;
1222  }
1223  case 4:
1224  {
1225  x = VIGRA_CSTD::fabs(x);
1226  if(x <= 1.0)
1227  {
1228  return 6.0 - 10.0*x;
1229  }
1230  else if(x < 2.0)
1231  {
1232  return -9.0 + 5.0*x;
1233  }
1234  else if(x < 3.0)
1235  {
1236  return 3.0 - x;
1237  }
1238  else
1239  return 0.0;
1240  }
1241  case 5:
1242  {
1243  return x < 0.0 ?
1244  x < -2.0 ?
1245  x < -3.0 ?
1246  0.0
1247  : 1.0
1248  : x < -1.0 ?
1249  -5.0
1250  : 10.0
1251  : x < 2.0 ?
1252  x < 1.0 ?
1253  -10.0
1254  : 5.0
1255  : x < 3.0 ?
1256  -1.0
1257  : 0.0;
1258  }
1259  default:
1260  return 0.0;
1261  }
1262 }
1263 
1264 typedef BSpline<5, double> QuinticBSplineKernel;
1265 
1266 #endif // NO_PARTIAL_TEMPLATE_SPECIALIZATION
1267 
1268 /********************************************************/
1269 /* */
1270 /* CatmullRomSpline */
1271 /* */
1272 /********************************************************/
1273 
1274 /** Interpolating 3-rd order splines.
1275 
1276  Implements the Catmull/Rom cardinal function
1277 
1278  \f[ f(x) = \left\{ \begin{array}{ll}
1279  \frac{3}{2}x^3 - \frac{5}{2}x^2 + 1 & |x| \leq 1 \\
1280  -\frac{1}{2}x^3 + \frac{5}{2}x^2 -4x + 2 & |x| \leq 2 \\
1281  0 & \mbox{otherwise}
1282  \end{array}\right.
1283  \f]
1284 
1285  It can be used as a functor, and as a kernel for
1286  \ref resamplingConvolveImage() to create a differentiable interpolant
1287  of an image. However, it should be noted that a twice differentiable
1288  interpolant can be created with only slightly more effort by recursive
1289  prefiltering followed by convolution with a 3rd order B-spline.
1290 
1291  <b>\#include</b> <vigra/splines.hxx><br>
1292  Namespace: vigra
1293 */
1294 template <class T = double>
1296 {
1297 public:
1298  /** the kernel's value type
1299  */
1300  typedef T value_type;
1301  /** the unary functor's argument type
1302  */
1303  typedef T argument_type;
1304  /** the unary functor's result type
1305  */
1306  typedef T result_type;
1307  /** the splines polynomial order
1308  */
1309  enum StaticOrder { order = 3 };
1310 
1311  /** function (functor) call
1312  */
1314 
1315  /** index operator -- same as operator()
1316  */
1317  T operator[] (T x) const
1318  { return operator()(x); }
1319 
1320  /** Radius of the function's support.
1321  Needed for \ref resamplingConvolveImage(), always 2.
1322  */
1323  int radius() const
1324  {return 2;}
1325 
1326  /** Derivative order of the function: always 0.
1327  */
1328  unsigned int derivativeOrder() const
1329  { return 0; }
1330 
1331  /** Prefilter coefficients for compatibility with \ref vigra::BSpline.
1332  (array has zero length, since prefiltering is not necessary).
1333  */
1335  {
1336  return prefilterCoefficients_;
1337  }
1338 
1339  protected:
1340  static ArrayVector<double> prefilterCoefficients_;
1341 };
1342 
1343 template <class T>
1344 ArrayVector<double> CatmullRomSpline<T>::prefilterCoefficients_;
1345 
1346 template <class T>
1347 typename CatmullRomSpline<T>::result_type
1349 {
1350  x = VIGRA_CSTD::fabs(x);
1351  if (x <= 1.0)
1352  {
1353  return 1.0 + x * x * (-2.5 + 1.5 * x);
1354  }
1355  else if (x >= 2.0)
1356  {
1357  return 0.0;
1358  }
1359  else
1360  {
1361  return 2.0 + x * (-4.0 + x * (2.5 - 0.5 * x));
1362  }
1363 }
1364 
1365 
1366 //@}
1367 
1368 } // namespace vigra
1369 
1370 
1371 #endif /* VIGRA_SPLINES_HXX */
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

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
vigra 1.11.1 (Fri May 19 2017)