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

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

© 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.9.0 (Sun Aug 10 2014)