[ 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 "polynomial.hxx"
43 #include "array_vector.hxx"
44 #include "fixedpoint.hxx"
45 
46 namespace vigra {
47 
48 namespace autodiff {
49 
50 template <class T, int N>
51 class DualVector;
52 
53 } // namespace autodiff
54 
55 /** \addtogroup MathFunctions Mathematical Functions
56 */
57 //@{
58 /* B-Splines of arbitrary order and interpolating Catmull/Rom splines.
59 
60  <b>\#include</b> <vigra/splines.hxx><br>
61  Namespace: vigra
62 */
63 #ifndef NO_PARTIAL_TEMPLATE_SPECIALIZATION
64 
65 /** Basic interface of the spline functors.
66 
67  Implements the spline functions defined by the recursion
68 
69  \f[ B_0(x) = \left\{ \begin{array}{ll}
70  1 & -\frac{1}{2} \leq x < \frac{1}{2} \\
71  0 & \mbox{otherwise}
72  \end{array}\right.
73  \f]
74 
75  and
76 
77  \f[ B_n(x) = B_0(x) * B_{n-1}(x)
78  \f]
79 
80  where * denotes convolution, and <i>n</i> is the spline order given by the
81  template parameter <tt>ORDER</tt>. These spline classes can be used as
82  unary and binary functors, as kernels for \ref resamplingConvolveImage(),
83  and as arguments for \ref vigra::SplineImageView. Note that the spline order
84  is given as a template argument.
85 
86  <b>\#include</b> <vigra/splines.hxx><br>
87  Namespace: vigra
88 */
89 template <int ORDER, class T = double>
91 {
92  public:
93 
94  /** the value type if used as a kernel in \ref resamplingConvolveImage().
95  */
96  typedef T value_type;
97  /** the functor's unary argument type
98  */
99  typedef T argument_type;
100  /** the functor's first binary argument type
101  */
103  /** the functor's second binary argument type
104  */
105  typedef unsigned int second_argument_type;
106  /** the functor's result type (unary and binary)
107  */
108  typedef T result_type;
109  /** the spline order
110  */
111  enum StaticOrder { order = ORDER };
112 
113  /** Create functor for given derivative of the spline. The spline's order
114  is specified spline by the template argument <TT>ORDER</tt>.
115  */
116  explicit BSplineBase(unsigned int derivativeOrder = 0)
117  : s1_(derivativeOrder)
118  {}
119 
120  /** Unary function call.
121  Returns the value of the spline with the derivative order given in the
122  constructor. Note that only derivatives up to <tt>ORDER-1</tt> are
123  continuous, and derivatives above <tt>ORDER+1</tt> are zero.
124  */
126  {
127  return exec(x, derivativeOrder());
128  }
129 
130  /** Binary function call.
131  The given derivative order is added to the derivative order
132  specified in the constructor. Note that only derivatives up to <tt>ORDER-1</tt> are
133  continuous, and derivatives above <tt>ORDER+1</tt> are zero.
134  */
136  {
137  return exec(x, derivativeOrder() + derivative_order);
138  }
139 
140  /** Index operator. Same as unary function call.
141  */
143  { return operator()(x); }
144 
145  /** Get the required filter radius for a discrete approximation of the
146  spline. Always equal to <tt>(ORDER + 1) / 2.0</tt>.
147  */
148  double radius() const
149  { return (ORDER + 1) * 0.5; }
150 
151  /** Get the derivative order of the Gaussian.
152  */
153  unsigned int derivativeOrder() const
154  { return s1_.derivativeOrder(); }
155 
156  /** Get the prefilter coefficients required for interpolation.
157  To interpolate with a B-spline, \ref resamplingConvolveImage()
158  can be used. However, the image to be interpolated must be
159  pre-filtered using \ref recursiveFilterX() and \ref recursiveFilterY()
160  with the filter coefficients given by this function. The length of the array
161  corresponds to how many times the above recursive filtering
162  has to be applied (zero length means no prefiltering necessary).
163  */
165  {
166  return prefilterCoefficients_;
167  }
168 
169  typedef ArrayVector<ArrayVector<T> > WeightMatrix;
170 
171  /** Get the coefficients to transform spline coefficients into
172  the coefficients of the corresponding polynomial.
173  Currently internally used in SplineImageView; needs more
174  documentation ???
175  */
176  static WeightMatrix const & weights()
177  {
178  return weightMatrix_;
179  }
180 
181  protected:
182  result_type exec(first_argument_type x, second_argument_type derivative_order) const;
183 
184  // factory function for the prefilter coefficients array
185  static ArrayVector<double> calculatePrefilterCoefficients();
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_(BSplineBase<ORDER, T>::calculatePrefilterCoefficients());
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>::calculatePrefilterCoefficients()
220 {
221  ArrayVector<double> res;
222  if(ORDER > 1)
223  {
224  const int r = ORDER / 2;
225  StaticPolynomial<2*r, double> p(2*r);
226  BSplineBase spline;
227  for(int i = 0; i <= 2*r; ++i)
228  p[i] = spline(T(i-r));
229  ArrayVector<double> roots;
230  polynomialRealRoots(p, 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]);
234  }
235  return res;
236 }
237 
238 template <int ORDER, class T>
239 typename BSplineBase<ORDER, T>::WeightMatrix
240 BSplineBase<ORDER, T>::calculateWeightMatrix()
241 {
242  WeightMatrix res(ORDER+1, ArrayVector<T>(ORDER+1));
243  double faculty = 1.0;
244  for(int d = 0; d <= ORDER; ++d)
245  {
246  if(d > 1)
247  faculty *= d;
248  double x = ORDER / 2; // (note: integer division)
249  BSplineBase spline;
250  for(int i = 0; i <= ORDER; ++i, --x)
251  res[d][i] = spline(x, d) / faculty;
252  }
253  return res;
254 }
255 
256 /********************************************************/
257 /* */
258 /* BSpline<N, T> */
259 /* */
260 /********************************************************/
261 
262 /** Spline functors for arbitrary orders.
263 
264  Provides the interface of \ref vigra::BSplineBase with a more convenient
265  name -- see there for more documentation.
266 */
267 template <int ORDER, class T = double>
268 class BSpline
269 : public BSplineBase<ORDER, T>
270 {
271  public:
272  /** Constructor forwarded to the base class constructor..
273  */
274  explicit BSpline(unsigned int derivativeOrder = 0)
275  : BSplineBase<ORDER, T>(derivativeOrder)
276  {}
277 };
278 
279 /********************************************************/
280 /* */
281 /* BSpline<0, T> */
282 /* */
283 /********************************************************/
284 
285 template <class T>
286 class BSplineBase<0, T>
287 {
288  public:
289 
290  typedef T value_type;
291  typedef T argument_type;
292  typedef T first_argument_type;
293  typedef unsigned int second_argument_type;
294  typedef T result_type;
295  enum StaticOrder { order = 0 };
296 
297  explicit BSplineBase(unsigned int derivativeOrder = 0)
298  : derivativeOrder_(derivativeOrder)
299  {}
300 
302  {
303  return exec(x, derivativeOrder_);
304  }
305 
306  template <unsigned int IntBits, unsigned int FracBits>
307  FixedPoint<IntBits, FracBits> operator()(FixedPoint<IntBits, FracBits> x) const
308  {
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);
313  }
314 
315  template <class U, int N>
316  autodiff::DualVector<U, N> operator()(autodiff::DualVector<U, N> const & x) const
317  {
318  return x < 0.5 && -0.5 <= x
319  ? autodiff::DualVector<U, N>(1.0)
320  : autodiff::DualVector<U, N>(0.0);
321  }
322 
324  {
325  return exec(x, derivativeOrder_ + derivative_order);
326  }
327 
329  { return operator()(x); }
330 
331  double radius() const
332  { return 0.5; }
333 
334  unsigned int derivativeOrder() const
335  { return derivativeOrder_; }
336 
337  ArrayVector<double> const & prefilterCoefficients() const
338  {
339  return prefilterCoefficients_;
340  }
341 
342  typedef T WeightMatrix[1][1];
343 
344  static WeightMatrix const & weights()
345  {
346  return weightMatrix_;
347  }
348 
349  protected:
350  result_type exec(first_argument_type x, second_argument_type derivative_order) const
351  {
352  if(derivative_order == 0)
353  return x < 0.5 && -0.5 <= x ?
354  1.0
355  : 0.0;
356  else
357  return 0.0;
358  }
359 
360  unsigned int derivativeOrder_;
361  static ArrayVector<double> prefilterCoefficients_;
362  static WeightMatrix weightMatrix_;
363 };
364 
365 template <class T>
366 ArrayVector<double> BSplineBase<0, T>::prefilterCoefficients_;
367 
368 template <class T>
369 typename BSplineBase<0, T>::WeightMatrix BSplineBase<0, T>::weightMatrix_ = {{ 1.0 }};
370 
371 /********************************************************/
372 /* */
373 /* BSpline<1, T> */
374 /* */
375 /********************************************************/
376 
377 template <class T>
378 class BSpline<1, T>
379 {
380  public:
381 
382  typedef T value_type;
383  typedef T argument_type;
384  typedef T first_argument_type;
385  typedef unsigned int second_argument_type;
386  typedef T result_type;
387  enum StaticOrder { order = 1 };
388 
389  explicit BSpline(unsigned int derivativeOrder = 0)
390  : derivativeOrder_(derivativeOrder)
391  {}
392 
394  {
395  return exec(x, derivativeOrder_);
396  }
397 
398  template <unsigned int IntBits, unsigned int FracBits>
399  FixedPoint<IntBits, FracBits> operator()(FixedPoint<IntBits, FracBits> x) const
400  {
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);
406  }
407 
408  template <class U, int N>
409  autodiff::DualVector<U, N> operator()(autodiff::DualVector<U, N> x) const
410  {
411  x = abs(x);
412  return x < 1.0
413  ? 1.0 - x
414  : autodiff::DualVector<U, N>(0.0);
415  }
416 
418  {
419  return exec(x, derivativeOrder_ + derivative_order);
420  }
421 
423  { return operator()(x); }
424 
425  double radius() const
426  { return 1.0; }
427 
428  unsigned int derivativeOrder() const
429  { return derivativeOrder_; }
430 
431  ArrayVector<double> const & prefilterCoefficients() const
432  {
433  return prefilterCoefficients_;
434  }
435 
436  typedef T WeightMatrix[2][2];
437 
438  static WeightMatrix const & weights()
439  {
440  return weightMatrix_;
441  }
442 
443  protected:
444  T exec(T x, unsigned int derivative_order) const;
445 
446  unsigned int derivativeOrder_;
447  static ArrayVector<double> prefilterCoefficients_;
448  static WeightMatrix weightMatrix_;
449 };
450 
451 template <class T>
452 ArrayVector<double> BSpline<1, T>::prefilterCoefficients_;
453 
454 template <class T>
455 typename BSpline<1, T>::WeightMatrix BSpline<1, T>::weightMatrix_ = {{ 1.0, 0.0}, {-1.0, 1.0}};
456 
457 template <class T>
458 T BSpline<1, T>::exec(T x, unsigned int derivative_order) const
459 {
460  switch(derivative_order)
461  {
462  case 0:
463  {
464  x = VIGRA_CSTD::fabs(x);
465  return x < 1.0 ?
466  1.0 - x
467  : 0.0;
468  }
469  case 1:
470  {
471  return x < 0.0 ?
472  -1.0 <= x ?
473  1.0
474  : 0.0
475  : x < 1.0 ?
476  -1.0
477  : 0.0;
478  }
479  default:
480  return 0.0;
481  }
482 }
483 
484 /********************************************************/
485 /* */
486 /* BSpline<2, T> */
487 /* */
488 /********************************************************/
489 
490 template <class T>
491 class BSpline<2, T>
492 {
493  public:
494 
495  typedef T value_type;
496  typedef T argument_type;
497  typedef T first_argument_type;
498  typedef unsigned int second_argument_type;
499  typedef T result_type;
500  enum StaticOrder { order = 2 };
501 
502  explicit BSpline(unsigned int derivativeOrder = 0)
503  : derivativeOrder_(derivativeOrder)
504  {}
505 
507  {
508  return exec(x, derivativeOrder_);
509  }
510 
511  template <unsigned int IntBits, unsigned int FracBits>
512  FixedPoint<IntBits, FracBits> operator()(FixedPoint<IntBits, FracBits> x) const
513  {
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);
521  return v == ONE_HALF
522  ? Value(ONE_HALF, FPNoShift)
523  : v <= ONE_HALF
524  ? Value(THREE_QUARTERS -
525  (int)(sq((unsigned)v >> PREMULTIPLY_SHIFT2) >> POSTMULTIPLY_SHIFT2), FPNoShift)
526  : v < THREE_HALVES
527  ? Value((int)(sq((unsigned)(THREE_HALVES-v) >> PREMULTIPLY_SHIFT1) >> (POSTMULTIPLY_SHIFT1 + 1)), FPNoShift)
528  : Value(0, FPNoShift);
529  }
530 
531  template <class U, int N>
532  autodiff::DualVector<U, N> operator()(autodiff::DualVector<U, N> x) const
533  {
534  x = abs(x);
535  return x < 0.5
536  ? 0.75 - x*x
537  : x < 1.5
538  ? 0.5 * sq(1.5 - x)
539  : autodiff::DualVector<U, N>(0.0);
540  }
541 
543  {
544  return exec(x, derivativeOrder_ + derivative_order);
545  }
546 
547  result_type dx(argument_type x) const
548  { return operator()(x, 1); }
549 
551  { return operator()(x); }
552 
553  double radius() const
554  { return 1.5; }
555 
556  unsigned int derivativeOrder() const
557  { return derivativeOrder_; }
558 
559  ArrayVector<double> const & prefilterCoefficients() const
560  {
561  return prefilterCoefficients_;
562  }
563 
564  typedef T WeightMatrix[3][3];
565 
566  static WeightMatrix const & weights()
567  {
568  return weightMatrix_;
569  }
570 
571  protected:
572  result_type exec(first_argument_type x, second_argument_type derivative_order) const;
573 
574  unsigned int derivativeOrder_;
575  static ArrayVector<double> prefilterCoefficients_;
576  static WeightMatrix weightMatrix_;
577 };
578 
579 template <class T>
580 ArrayVector<double> BSpline<2, T>::prefilterCoefficients_(1, 2.0*M_SQRT2 - 3.0);
581 
582 template <class T>
583 typename BSpline<2, T>::WeightMatrix BSpline<2, T>::weightMatrix_ =
584  {{ 0.125, 0.75, 0.125},
585  {-0.5, 0.0, 0.5},
586  { 0.5, -1.0, 0.5}};
587 
588 template <class T>
589 typename BSpline<2, T>::result_type
590 BSpline<2, T>::exec(first_argument_type x, second_argument_type derivative_order) const
591 {
592  switch(derivative_order)
593  {
594  case 0:
595  {
596  x = VIGRA_CSTD::fabs(x);
597  return x < 0.5 ?
598  0.75 - x*x
599  : x < 1.5 ?
600  0.5 * sq(1.5 - x)
601  : 0.0;
602  }
603  case 1:
604  {
605  return x >= -0.5 ?
606  x <= 0.5 ?
607  -2.0 * x
608  : x < 1.5 ?
609  x - 1.5
610  : 0.0
611  : x > -1.5 ?
612  x + 1.5
613  : 0.0;
614  }
615  case 2:
616  {
617  return x >= -0.5 ?
618  x < 0.5 ?
619  -2.0
620  : x < 1.5 ?
621  1.0
622  : 0.0
623  : x >= -1.5 ?
624  1.0
625  : 0.0;
626  }
627  default:
628  return 0.0;
629  }
630 }
631 
632 /********************************************************/
633 /* */
634 /* BSpline<3, T> */
635 /* */
636 /********************************************************/
637 
638 template <class T>
639 class BSpline<3, T>
640 {
641  public:
642 
643  typedef T value_type;
644  typedef T argument_type;
645  typedef T first_argument_type;
646  typedef unsigned int second_argument_type;
647  typedef T result_type;
648  enum StaticOrder { order = 3 };
649 
650  explicit BSpline(unsigned int derivativeOrder = 0)
651  : derivativeOrder_(derivativeOrder)
652  {}
653 
655  {
656  return exec(x, derivativeOrder_);
657  }
658 
659  template <unsigned int IntBits, unsigned int FracBits>
660  FixedPoint<IntBits, FracBits> operator()(FixedPoint<IntBits, FracBits> x) const
661  {
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);
667  return v == ONE
668  ? Value(ONE_SIXTH, FPNoShift)
669  : v < ONE
670  ? Value(TWO_THIRDS +
671  (((int)(sq((unsigned)v >> PREMULTIPLY_SHIFT) >> (POSTMULTIPLY_SHIFT + PREMULTIPLY_SHIFT))
672  * (((v >> 1) - ONE) >> PREMULTIPLY_SHIFT)) >> POSTMULTIPLY_SHIFT), FPNoShift)
673  : v < TWO
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);
677  }
678 
679  template <class U, int N>
680  autodiff::DualVector<U, N> operator()(autodiff::DualVector<U, N> x) const
681  {
682  x = abs(x);
683  if(x < 1.0)
684  {
685  return 2.0/3.0 + x*x*(-1.0 + 0.5*x);
686  }
687  else if(x < 2.0)
688  {
689  x = 2.0 - x;
690  return x*x*x/6.0;
691  }
692  else
693  return autodiff::DualVector<U, N>(0.0);
694  }
695 
697  {
698  return exec(x, derivativeOrder_ + derivative_order);
699  }
700 
701  result_type dx(argument_type x) const
702  { return operator()(x, 1); }
703 
704  result_type dxx(argument_type x) const
705  { return operator()(x, 2); }
706 
708  { return operator()(x); }
709 
710  double radius() const
711  { return 2.0; }
712 
713  unsigned int derivativeOrder() const
714  { return derivativeOrder_; }
715 
716  ArrayVector<double> const & prefilterCoefficients() const
717  {
718  return prefilterCoefficients_;
719  }
720 
721  typedef T WeightMatrix[4][4];
722 
723  static WeightMatrix const & weights()
724  {
725  return weightMatrix_;
726  }
727 
728  protected:
729  result_type exec(first_argument_type x, second_argument_type derivative_order) const;
730 
731  unsigned int derivativeOrder_;
732  static ArrayVector<double> prefilterCoefficients_;
733  static WeightMatrix weightMatrix_;
734 };
735 
736 template <class T>
737 ArrayVector<double> BSpline<3, T>::prefilterCoefficients_(1, VIGRA_CSTD::sqrt(3.0) - 2.0);
738 
739 template <class T>
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}};
745 
746 template <class T>
747 typename BSpline<3, T>::result_type
748 BSpline<3, T>::exec(first_argument_type x, second_argument_type derivative_order) const
749 {
750  switch(derivative_order)
751  {
752  case 0:
753  {
754  x = VIGRA_CSTD::fabs(x);
755  if(x < 1.0)
756  {
757  return 2.0/3.0 + x*x*(-1.0 + 0.5*x);
758  }
759  else if(x < 2.0)
760  {
761  x = 2.0 - x;
762  return x*x*x/6.0;
763  }
764  else
765  return 0.0;
766  }
767  case 1:
768  {
769  double s = x < 0.0 ?
770  -1.0
771  : 1.0;
772  x = VIGRA_CSTD::fabs(x);
773  return x < 1.0 ?
774  s*x*(-2.0 + 1.5*x)
775  : x < 2.0 ?
776  -0.5*s*sq(2.0 - x)
777  : 0.0;
778  }
779  case 2:
780  {
781  x = VIGRA_CSTD::fabs(x);
782  return x < 1.0 ?
783  3.0*x - 2.0
784  : x < 2.0 ?
785  2.0 - x
786  : 0.0;
787  }
788  case 3:
789  {
790  return x < 0.0 ?
791  x < -1.0 ?
792  x < -2.0 ?
793  0.0
794  : 1.0
795  : -3.0
796  : x < 1.0 ?
797  3.0
798  : x < 2.0 ?
799  -1.0
800  : 0.0;
801  }
802  default:
803  return 0.0;
804  }
805 }
806 
807 typedef BSpline<3, double> CubicBSplineKernel;
808 
809 /********************************************************/
810 /* */
811 /* BSpline<4, T> */
812 /* */
813 /********************************************************/
814 
815 template <class T>
816 class BSpline<4, T>
817 {
818  public:
819 
820  typedef T value_type;
821  typedef T argument_type;
822  typedef T first_argument_type;
823  typedef unsigned int second_argument_type;
824  typedef T result_type;
825  enum StaticOrder { order = 4 };
826 
827  explicit BSpline(unsigned int derivativeOrder = 0)
828  : derivativeOrder_(derivativeOrder)
829  {}
830 
832  {
833  return exec(x, derivativeOrder_);
834  }
835 
837  {
838  return exec(x, derivativeOrder_ + derivative_order);
839  }
840 
841  template <class U, int N>
842  autodiff::DualVector<U, N> operator()(autodiff::DualVector<U, N> x) const
843  {
844  x = abs(x);
845  if(x <= 0.5)
846  {
847  return 115.0/192.0 + x*x*(-0.625 + x*x*0.25);
848  }
849  else if(x < 1.5)
850  {
851  return (55.0/16.0 + x*(1.25 + x*(-7.5 + x*(5.0 - x)))) / 6.0;
852  }
853  else if(x < 2.5)
854  {
855  x = 2.5 - x;
856  return sq(x*x) / 24.0;
857  }
858  else
859  return autodiff::DualVector<U, N>(0.0);
860  }
861 
862  result_type dx(argument_type x) const
863  { return operator()(x, 1); }
864 
865  result_type dxx(argument_type x) const
866  { return operator()(x, 2); }
867 
868  result_type dx3(argument_type x) const
869  { return operator()(x, 3); }
870 
872  { return operator()(x); }
873 
874  double radius() const
875  { return 2.5; }
876 
877  unsigned int derivativeOrder() const
878  { return derivativeOrder_; }
879 
880  ArrayVector<double> const & prefilterCoefficients() const
881  {
882  return prefilterCoefficients_;
883  }
884 
885  typedef T WeightMatrix[5][5];
886 
887  static WeightMatrix const & weights()
888  {
889  return weightMatrix_;
890  }
891 
892  protected:
893  result_type exec(first_argument_type x, second_argument_type derivative_order) const;
894 
895  static ArrayVector<double> calculatePrefilterCoefficients()
896  {
897  ArrayVector<double> b(2);
898  // -19 + 4*sqrt(19) + 2*sqrt(2*(83 - 19*sqrt(19)))
899  b[0] = -0.361341225900220177092212841325;
900  // -19 - 4*sqrt(19) + 2*sqrt(2*(83 + 19*sqrt(19)))
901  b[1] = -0.013725429297339121360331226939;
902  return b;
903  }
904 
905  unsigned int derivativeOrder_;
906  static ArrayVector<double> prefilterCoefficients_;
907  static WeightMatrix weightMatrix_;
908 };
909 
910 template <class T>
911 ArrayVector<double> BSpline<4, T>::prefilterCoefficients_(calculatePrefilterCoefficients());
912 
913 template <class T>
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}};
920 
921 template <class T>
922 typename BSpline<4, T>::result_type
923 BSpline<4, T>::exec(first_argument_type x, second_argument_type derivative_order) const
924 {
925  switch(derivative_order)
926  {
927  case 0:
928  {
929  x = VIGRA_CSTD::fabs(x);
930  if(x <= 0.5)
931  {
932  return 115.0/192.0 + x*x*(-0.625 + x*x*0.25);
933  }
934  else if(x < 1.5)
935  {
936  return (55.0/16.0 + x*(1.25 + x*(-7.5 + x*(5.0 - x)))) / 6.0;
937  }
938  else if(x < 2.5)
939  {
940  x = 2.5 - x;
941  return sq(x*x) / 24.0;
942  }
943  else
944  return 0.0;
945  }
946  case 1:
947  {
948  double s = x < 0.0 ?
949  -1.0 :
950  1.0;
951  x = VIGRA_CSTD::fabs(x);
952  if(x <= 0.5)
953  {
954  return s*x*(-1.25 + x*x);
955  }
956  else if(x < 1.5)
957  {
958  return s*(5.0 + x*(-60.0 + x*(60.0 - 16.0*x))) / 24.0;
959  }
960  else if(x < 2.5)
961  {
962  x = 2.5 - x;
963  return s*x*x*x / -6.0;
964  }
965  else
966  return 0.0;
967  }
968  case 2:
969  {
970  x = VIGRA_CSTD::fabs(x);
971  if(x <= 0.5)
972  {
973  return -1.25 + 3.0*x*x;
974  }
975  else if(x < 1.5)
976  {
977  return -2.5 + x*(5.0 - 2.0*x);
978  }
979  else if(x < 2.5)
980  {
981  x = 2.5 - x;
982  return x*x / 2.0;
983  }
984  else
985  return 0.0;
986  }
987  case 3:
988  {
989  double s = x < 0.0 ?
990  -1.0 :
991  1.0;
992  x = VIGRA_CSTD::fabs(x);
993  if(x <= 0.5)
994  {
995  return s*x*6.0;
996  }
997  else if(x < 1.5)
998  {
999  return s*(5.0 - 4.0*x);
1000  }
1001  else if(x < 2.5)
1002  {
1003  return s*(x - 2.5);
1004  }
1005  else
1006  return 0.0;
1007  }
1008  case 4:
1009  {
1010  return x < 0.0
1011  ? x < -2.5
1012  ? 0.0
1013  : x < -1.5
1014  ? 1.0
1015  : x < -0.5
1016  ? -4.0
1017  : 6.0
1018  : x < 0.5
1019  ? 6.0
1020  : x < 1.5
1021  ? -4.0
1022  : x < 2.5
1023  ? 1.0
1024  : 0.0;
1025  }
1026  default:
1027  return 0.0;
1028  }
1029 }
1030 
1031 typedef BSpline<4, double> QuarticBSplineKernel;
1032 
1033 /********************************************************/
1034 /* */
1035 /* BSpline<5, T> */
1036 /* */
1037 /********************************************************/
1038 
1039 template <class T>
1040 class BSpline<5, T>
1041 {
1042  public:
1043 
1044  typedef T value_type;
1045  typedef T argument_type;
1046  typedef T first_argument_type;
1047  typedef unsigned int second_argument_type;
1048  typedef T result_type;
1049  enum StaticOrder { order = 5 };
1050 
1051  explicit BSpline(unsigned int derivativeOrder = 0)
1052  : derivativeOrder_(derivativeOrder)
1053  {}
1054 
1056  {
1057  return exec(x, derivativeOrder_);
1058  }
1059 
1061  {
1062  return exec(x, derivativeOrder_ + derivative_order);
1063  }
1064 
1065  template <class U, int N>
1066  autodiff::DualVector<U, N> operator()(autodiff::DualVector<U, N> x) const
1067  {
1068  x = abs(x);
1069  if(x <= 1.0)
1070  {
1071  return 0.55 + x*x*(-0.5 + x*x*(0.25 - x/12.0));
1072  }
1073  else if(x < 2.0)
1074  {
1075  return 17.0/40.0 + x*(0.625 + x*(-1.75 + x*(1.25 + x*(-0.375 + x/24.0))));
1076  }
1077  else if(x < 3.0)
1078  {
1079  x = 3.0 - x;
1080  return x*sq(x*x) / 120.0;
1081  }
1082  else
1083  return autodiff::DualVector<U, N>(0.0);
1084  }
1085 
1086  result_type dx(argument_type x) const
1087  { return operator()(x, 1); }
1088 
1089  result_type dxx(argument_type x) const
1090  { return operator()(x, 2); }
1091 
1092  result_type dx3(argument_type x) const
1093  { return operator()(x, 3); }
1094 
1095  result_type dx4(argument_type x) const
1096  { return operator()(x, 4); }
1097 
1099  { return operator()(x); }
1100 
1101  double radius() const
1102  { return 3.0; }
1103 
1104  unsigned int derivativeOrder() const
1105  { return derivativeOrder_; }
1106 
1107  ArrayVector<double> const & prefilterCoefficients() const
1108  {
1109  return prefilterCoefficients_;
1110  }
1111 
1112  typedef T WeightMatrix[6][6];
1113 
1114  static WeightMatrix const & weights()
1115  {
1116  return weightMatrix_;
1117  }
1118 
1119  protected:
1120  result_type exec(first_argument_type x, second_argument_type derivative_order) const;
1121 
1122  static ArrayVector<double> calculatePrefilterCoefficients()
1123  {
1124  ArrayVector<double> b(2);
1125  // -(13/2) + sqrt(105)/2 + sqrt(1/2*((135 - 13*sqrt(105))))
1126  b[0] = -0.430575347099973791851434783493;
1127  // (1/2)*((-13) - sqrt(105) + sqrt(2*((135 + 13*sqrt(105)))))
1128  b[1] = -0.043096288203264653822712376822;
1129  return b;
1130  }
1131 
1132  unsigned int derivativeOrder_;
1133  static ArrayVector<double> prefilterCoefficients_;
1134  static WeightMatrix weightMatrix_;
1135 };
1136 
1137 template <class T>
1138 ArrayVector<double> BSpline<5, T>::prefilterCoefficients_(calculatePrefilterCoefficients());
1139 
1140 template <class T>
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}};
1148 
1149 template <class T>
1150 typename BSpline<5, T>::result_type
1151 BSpline<5, T>::exec(first_argument_type x, second_argument_type derivative_order) const
1152 {
1153  switch(derivative_order)
1154  {
1155  case 0:
1156  {
1157  x = VIGRA_CSTD::fabs(x);
1158  if(x <= 1.0)
1159  {
1160  return 0.55 + x*x*(-0.5 + x*x*(0.25 - x/12.0));
1161  }
1162  else if(x < 2.0)
1163  {
1164  return 17.0/40.0 + x*(0.625 + x*(-1.75 + x*(1.25 + x*(-0.375 + x/24.0))));
1165  }
1166  else if(x < 3.0)
1167  {
1168  x = 3.0 - x;
1169  return x*sq(x*x) / 120.0;
1170  }
1171  else
1172  return 0.0;
1173  }
1174  case 1:
1175  {
1176  double s = x < 0.0 ?
1177  -1.0 :
1178  1.0;
1179  x = VIGRA_CSTD::fabs(x);
1180  if(x <= 1.0)
1181  {
1182  return s*x*(-1.0 + x*x*(1.0 - 5.0/12.0*x));
1183  }
1184  else if(x < 2.0)
1185  {
1186  return s*(0.625 + x*(-3.5 + x*(3.75 + x*(-1.5 + 5.0/24.0*x))));
1187  }
1188  else if(x < 3.0)
1189  {
1190  x = 3.0 - x;
1191  return s*sq(x*x) / -24.0;
1192  }
1193  else
1194  return 0.0;
1195  }
1196  case 2:
1197  {
1198  x = VIGRA_CSTD::fabs(x);
1199  if(x <= 1.0)
1200  {
1201  return -1.0 + x*x*(3.0 -5.0/3.0*x);
1202  }
1203  else if(x < 2.0)
1204  {
1205  return -3.5 + x*(7.5 + x*(-4.5 + 5.0/6.0*x));
1206  }
1207  else if(x < 3.0)
1208  {
1209  x = 3.0 - x;
1210  return x*x*x / 6.0;
1211  }
1212  else
1213  return 0.0;
1214  }
1215  case 3:
1216  {
1217  double s = x < 0.0 ?
1218  -1.0 :
1219  1.0;
1220  x = VIGRA_CSTD::fabs(x);
1221  if(x <= 1.0)
1222  {
1223  return s*x*(6.0 - 5.0*x);
1224  }
1225  else if(x < 2.0)
1226  {
1227  return s*(7.5 + x*(-9.0 + 2.5*x));
1228  }
1229  else if(x < 3.0)
1230  {
1231  x = 3.0 - x;
1232  return -0.5*s*x*x;
1233  }
1234  else
1235  return 0.0;
1236  }
1237  case 4:
1238  {
1239  x = VIGRA_CSTD::fabs(x);
1240  if(x <= 1.0)
1241  {
1242  return 6.0 - 10.0*x;
1243  }
1244  else if(x < 2.0)
1245  {
1246  return -9.0 + 5.0*x;
1247  }
1248  else if(x < 3.0)
1249  {
1250  return 3.0 - x;
1251  }
1252  else
1253  return 0.0;
1254  }
1255  case 5:
1256  {
1257  return x < 0.0 ?
1258  x < -2.0 ?
1259  x < -3.0 ?
1260  0.0
1261  : 1.0
1262  : x < -1.0 ?
1263  -5.0
1264  : 10.0
1265  : x < 2.0 ?
1266  x < 1.0 ?
1267  -10.0
1268  : 5.0
1269  : x < 3.0 ?
1270  -1.0
1271  : 0.0;
1272  }
1273  default:
1274  return 0.0;
1275  }
1276 }
1277 
1278 typedef BSpline<5, double> QuinticBSplineKernel;
1279 
1280 #endif // NO_PARTIAL_TEMPLATE_SPECIALIZATION
1281 
1282 /********************************************************/
1283 /* */
1284 /* CatmullRomSpline */
1285 /* */
1286 /********************************************************/
1287 
1288 /** Interpolating 3-rd order splines.
1289 
1290  Implements the Catmull/Rom cardinal function
1291 
1292  \f[ f(x) = \left\{ \begin{array}{ll}
1293  \frac{3}{2}x^3 - \frac{5}{2}x^2 + 1 & |x| \leq 1 \\
1294  -\frac{1}{2}x^3 + \frac{5}{2}x^2 -4x + 2 & |x| \leq 2 \\
1295  0 & \mbox{otherwise}
1296  \end{array}\right.
1297  \f]
1298 
1299  It can be used as a functor, and as a kernel for
1300  \ref resamplingConvolveImage() to create a differentiable interpolant
1301  of an image. However, it should be noted that a twice differentiable
1302  interpolant can be created with only slightly more effort by recursive
1303  prefiltering followed by convolution with a 3rd order B-spline.
1304 
1305  <b>\#include</b> <vigra/splines.hxx><br>
1306  Namespace: vigra
1307 */
1308 template <class T = double>
1310 {
1311 public:
1312  /** the kernel's value type
1313  */
1314  typedef T value_type;
1315  /** the unary functor's argument type
1316  */
1317  typedef T argument_type;
1318  /** the unary functor's result type
1319  */
1320  typedef T result_type;
1321  /** the splines polynomial order
1322  */
1323  enum StaticOrder { order = 3 };
1324 
1325  /** function (functor) call
1326  */
1328 
1329  /** index operator -- same as operator()
1330  */
1331  T operator[] (T x) const
1332  { return operator()(x); }
1333 
1334  /** Radius of the function's support.
1335  Needed for \ref resamplingConvolveImage(), always 2.
1336  */
1337  int radius() const
1338  {return 2;}
1339 
1340  /** Derivative order of the function: always 0.
1341  */
1342  unsigned int derivativeOrder() const
1343  { return 0; }
1344 
1345  /** Prefilter coefficients for compatibility with \ref vigra::BSpline.
1346  (array has zero length, since prefiltering is not necessary).
1347  */
1349  {
1350  return prefilterCoefficients_;
1351  }
1352 
1353  protected:
1354  static ArrayVector<double> prefilterCoefficients_;
1355 };
1356 
1357 template <class T>
1358 ArrayVector<double> CatmullRomSpline<T>::prefilterCoefficients_;
1359 
1360 template <class T>
1361 typename CatmullRomSpline<T>::result_type
1363 {
1364  x = VIGRA_CSTD::fabs(x);
1365  if (x <= 1.0)
1366  {
1367  return 1.0 + x * x * (-2.5 + 1.5 * x);
1368  }
1369  else if (x >= 2.0)
1370  {
1371  return 0.0;
1372  }
1373  else
1374  {
1375  return 2.0 + x * (-4.0 + x * (2.5 - 0.5 * x));
1376  }
1377 }
1378 
1379 
1380 //@}
1381 
1382 } // namespace vigra
1383 
1384 
1385 #endif /* VIGRA_SPLINES_HXX */
T operator[](T x) const
Definition: splines.hxx:1331
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
ArrayVector< double > const & prefilterCoefficients() const
Definition: splines.hxx:164
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: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
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

© 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.10.0 (Fri Feb 21 2014)