[ 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 "array_vector.hxx"
43#include "fixedpoint.hxx"
44
45namespace vigra {
46
47namespace autodiff {
48
49template <class T, int N>
50class 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*/
88template <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
195template <int ORDER, class T>
196ArrayVector<double> BSplineBase<ORDER, T>::prefilterCoefficients_(getPrefilterCoefficients());
197
198template <int ORDER, class T>
199typename BSplineBase<ORDER, T>::WeightMatrix BSplineBase<ORDER, T>::weightMatrix_(calculateWeightMatrix());
200
201template <int ORDER, class T>
203BSplineBase<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
217template <int ORDER, class T>
219BSplineBase<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
244template <int ORDER, class T>
245typename BSplineBase<ORDER, T>::WeightMatrix
246BSplineBase<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*/
273template <int ORDER, class T = double>
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
291template <class T>
292class 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
371template <class T>
372ArrayVector<double> BSplineBase<0, T>::prefilterCoefficients_;
373
374template <class T>
375typename BSplineBase<0, T>::WeightMatrix BSplineBase<0, T>::weightMatrix_ = {{ 1.0 }};
376
377/********************************************************/
378/* */
379/* BSpline<1, T> */
380/* */
381/********************************************************/
382
383template <class T>
384class 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
457template <class T>
458ArrayVector<double> BSpline<1, T>::prefilterCoefficients_;
459
460template <class T>
461typename BSpline<1, T>::WeightMatrix BSpline<1, T>::weightMatrix_ = {{ 1.0, 0.0}, {-1.0, 1.0}};
462
463template <class T>
464T 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
496template <class T>
497class 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
585template <class T>
586ArrayVector<double> BSpline<2, T>::prefilterCoefficients_(1, 2.0*M_SQRT2 - 3.0);
587
588template <class T>
589typename 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
594template <class T>
596BSpline<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
644template <class T>
645class 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
742template <class T>
743ArrayVector<double> BSpline<3, T>::prefilterCoefficients_(1, VIGRA_CSTD::sqrt(3.0) - 2.0);
744
745template <class T>
746typename 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
752template <class T>
754BSpline<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
813typedef BSpline<3, double> CubicBSplineKernel;
814
815/********************************************************/
816/* */
817/* BSpline<4, T> */
818/* */
819/********************************************************/
820
821template <class T>
822class 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
906template <class T>
907ArrayVector<double> BSpline<4, T>::prefilterCoefficients_(BSplineBase<4, T>::getPrefilterCoefficients());
908
909template <class T>
910typename 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
917template <class T>
919BSpline<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
1027typedef BSpline<4, double> QuarticBSplineKernel;
1028
1029/********************************************************/
1030/* */
1031/* BSpline<5, T> */
1032/* */
1033/********************************************************/
1034
1035template <class T>
1036class 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
1123template <class T>
1124ArrayVector<double> BSpline<5, T>::prefilterCoefficients_(BSplineBase<5, T>::getPrefilterCoefficients());
1125
1126template <class T>
1127typename 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
1135template <class T>
1137BSpline<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
1264typedef 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*/
1294template <class T = double>
1296{
1297public:
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
1343template <class T>
1344ArrayVector<double> CatmullRomSpline<T>::prefilterCoefficients_;
1345
1346template <class T>
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 */
Definition array_vector.hxx:514
Definition splines.hxx:90
value_type operator[](value_type x) const
Definition splines.hxx:143
unsigned int second_argument_type
Definition splines.hxx:106
T value_type
Definition splines.hxx:97
double radius() const
Definition splines.hxx:149
ArrayVector< double > const & prefilterCoefficients() const
Definition splines.hxx:165
result_type operator()(argument_type x) const
Definition splines.hxx:126
T first_argument_type
Definition splines.hxx:103
StaticOrder
Definition splines.hxx:112
T result_type
Definition splines.hxx:109
BSplineBase(unsigned int derivativeOrder=0)
Definition splines.hxx:117
result_type operator()(first_argument_type x, second_argument_type derivative_order) const
Definition splines.hxx:136
T argument_type
Definition splines.hxx:100
static WeightMatrix const & weights()
Definition splines.hxx:177
unsigned int derivativeOrder() const
Definition splines.hxx:154
Definition splines.hxx:276
BSpline(unsigned int derivativeOrder=0)
Definition splines.hxx:280
Definition splines.hxx:1296
result_type operator()(argument_type x) const
Definition splines.hxx:1348
T value_type
Definition splines.hxx:1300
ArrayVector< double > const & prefilterCoefficients() const
Definition splines.hxx:1334
T operator[](T x) const
Definition splines.hxx:1317
StaticOrder
Definition splines.hxx:1309
int radius() const
Definition splines.hxx:1323
T result_type
Definition splines.hxx:1306
T argument_type
Definition splines.hxx:1303
unsigned int derivativeOrder() const
Definition splines.hxx:1328
Definition autodiff.hxx:88
NumericTraits< T >::Promote sq(T t)
The square function.
Definition mathutil.hxx:382
FFTWComplex< R >::NormType abs(const FFTWComplex< R > &a)
absolute value (= magnitude)
Definition fftw3.hxx:1002

© 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.12.1 (Thu Feb 27 2025)