[ VIGRA Homepage | Class Index | Function Index | File Index | Main Page ]
![]() |
vigra/splines.hxx | ![]() |
---|
00001 /************************************************************************/ 00002 /* */ 00003 /* Copyright 1998-2004 by Ullrich Koethe */ 00004 /* Cognitive Systems Group, University of Hamburg, Germany */ 00005 /* */ 00006 /* This file is part of the VIGRA computer vision library. */ 00007 /* ( Version 1.3.3, Aug 18 2005 ) */ 00008 /* You may use, modify, and distribute this software according */ 00009 /* to the terms stated in the LICENSE file included in */ 00010 /* the VIGRA distribution. */ 00011 /* */ 00012 /* The VIGRA Website is */ 00013 /* http://kogs-www.informatik.uni-hamburg.de/~koethe/vigra/ */ 00014 /* Please direct questions, bug reports, and contributions to */ 00015 /* koethe@informatik.uni-hamburg.de */ 00016 /* */ 00017 /* THIS SOFTWARE IS PROVIDED AS IS AND WITHOUT ANY EXPRESS OR */ 00018 /* IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED */ 00019 /* WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. */ 00020 /* */ 00021 /************************************************************************/ 00022 00023 #ifndef VIGRA_SPLINES_HXX 00024 #define VIGRA_SPLINES_HXX 00025 00026 #include <cmath> 00027 #include "vigra/config.hxx" 00028 #include "vigra/mathutil.hxx" 00029 #include "vigra/polynomial.hxx" 00030 #include "vigra/array_vector.hxx" 00031 #include "vigra/fixedpoint.hxx" 00032 00033 namespace vigra { 00034 00035 /** \addtogroup MathFunctions Mathematical Functions 00036 */ 00037 //@{ 00038 /* B-Splines of arbitrary order and interpolating Catmull/Rom splines. 00039 00040 <b>\#include</b> "<a href="splines_8hxx-source.html">vigra/splines.hxx</a>"<br> 00041 Namespace: vigra 00042 */ 00043 #ifndef NO_PARTIAL_TEMPLATE_SPECIALIZATION 00044 00045 /** Basic interface of the spline functors. 00046 00047 Implements the spline functions defined by the recursion 00048 00049 \f[ B_0(x) = \left\{ \begin{array}{ll} 00050 1 & -\frac{1}{2} \leq x < \frac{1}{2} \\ 00051 0 & \mbox{otherwise} 00052 \end{array}\right. 00053 \f] 00054 00055 and 00056 00057 \f[ B_n(x) = B_0(x) * B_{n-1}(x) 00058 \f] 00059 00060 where * denotes convolution, and <i>n</i> is the spline order given by the 00061 template parameter <tt>ORDER</tt>. These spline classes can be used as 00062 unary and binary functors, as kernels for \ref resamplingConvolveImage(), 00063 and as arguments for \ref vigra::SplineImageView. Note that the spline order 00064 is given as a template argument. 00065 00066 <b>\#include</b> "<a href="splines_8hxx-source.html">vigra/splines.hxx</a>"<br> 00067 Namespace: vigra 00068 */ 00069 template <int ORDER, class T = double> 00070 class BSplineBase 00071 { 00072 public: 00073 00074 /** the value type if used as a kernel in \ref resamplingConvolveImage(). 00075 */ 00076 typedef T value_type; 00077 /** the functor's unary argument type 00078 */ 00079 typedef T argument_type; 00080 /** the functor's first binary argument type 00081 */ 00082 typedef T first_argument_type; 00083 /** the functor's second binary argument type 00084 */ 00085 typedef unsigned int second_argument_type; 00086 /** the functor's result type (unary and binary) 00087 */ 00088 typedef T result_type; 00089 /** the spline order 00090 */ 00091 enum StaticOrder { order = ORDER }; 00092 00093 /** Create functor for gevine derivative of the spline. The spline's order 00094 is specified spline by the template argument <TT>ORDER</tt>. 00095 */ 00096 explicit BSplineBase(unsigned int derivativeOrder = 0) 00097 : s1_(derivativeOrder) 00098 {} 00099 00100 /** Unary function call. 00101 Returns the value of the spline with the derivative order given in the 00102 constructor. Note that only derivatives up to <tt>ORDER-1</tt> are 00103 continous, and derivatives above <tt>ORDER+1</tt> are zero. 00104 */ 00105 result_type operator()(argument_type x) const 00106 { 00107 return exec(x, derivativeOrder()); 00108 } 00109 00110 /** Binary function call. 00111 The given derivative order is added to the derivative order 00112 specified in the constructor. Note that only derivatives up to <tt>ORDER-1</tt> are 00113 continous, and derivatives above <tt>ORDER+1</tt> are zero. 00114 */ 00115 result_type operator()(first_argument_type x, second_argument_type derivative_order) const 00116 { 00117 return exec(x, derivativeOrder() + derivative_order); 00118 } 00119 00120 /** Index operator. Same as unary function call. 00121 */ 00122 value_type operator[](value_type x) const 00123 { return operator()(x); } 00124 00125 /** Get the required filter radius for a discrete approximation of the 00126 spline. Always equal to <tt>(ORDER + 1) / 2.0</tt>. 00127 */ 00128 double radius() const 00129 { return (ORDER + 1) * 0.5; } 00130 00131 /** Get the derivative order of the Gaussian. 00132 */ 00133 unsigned int derivativeOrder() const 00134 { return s1_.derivativeOrder(); } 00135 00136 /** Get the prefilter coefficients required for interpolation. 00137 To interpolate with a B-spline, \ref resamplingConvolveImage() 00138 can be used. However, the image to be interpolated must be 00139 pre-filtered using \ref recursiveFilterImage() with the filter 00140 coefficients given by this function. The length of the array 00141 corresponds to the number of times \ref recursiveFilterImage() 00142 has to be applied (zero length means no prefiltering necessary). 00143 */ 00144 ArrayVector<double> const & prefilterCoefficients() const 00145 { 00146 static ArrayVector<double> const & b = calculatePrefilterCoefficients(); 00147 return b; 00148 } 00149 00150 static ArrayVector<double> const & calculatePrefilterCoefficients(); 00151 00152 typedef T WeightMatrix[ORDER+1][ORDER+1]; 00153 00154 /** Get the coefficients to transform spline coefficients into 00155 the coefficients of the corresponding polynomial. 00156 Currently internally used in SplineImageView; needs more 00157 documentation ??? 00158 */ 00159 static WeightMatrix & weights() 00160 { 00161 static WeightMatrix & b = calculateWeightMatrix(); 00162 return b; 00163 } 00164 00165 static WeightMatrix & calculateWeightMatrix(); 00166 00167 protected: 00168 result_type exec(first_argument_type x, second_argument_type derivative_order) const; 00169 00170 BSplineBase<ORDER-1, T> s1_; 00171 }; 00172 00173 template <int ORDER, class T> 00174 typename BSplineBase<ORDER, T>::result_type 00175 BSplineBase<ORDER, T>::exec(first_argument_type x, second_argument_type derivative_order) const 00176 { 00177 if(derivative_order == 0) 00178 { 00179 T n12 = (ORDER + 1.0) / 2.0; 00180 return ((n12 + x) * s1_(x + 0.5) + (n12 - x) * s1_(x - 0.5)) / ORDER; 00181 } 00182 else 00183 { 00184 --derivative_order; 00185 return s1_(x + 0.5, derivative_order) - s1_(x - 0.5, derivative_order); 00186 } 00187 } 00188 00189 template <int ORDER, class T> 00190 ArrayVector<double> const & BSplineBase<ORDER, T>::calculatePrefilterCoefficients() 00191 { 00192 static ArrayVector<double> b; 00193 if(ORDER > 1) 00194 { 00195 static const int r = ORDER / 2; 00196 StaticPolynomial<2*r, double> p(2*r); 00197 BSplineBase spline; 00198 for(int i = 0; i <= 2*r; ++i) 00199 p[i] = spline(T(i-r)); 00200 ArrayVector<double> roots; 00201 polynomialRealRoots(p, roots); 00202 for(unsigned int i = 0; i < roots.size(); ++i) 00203 if(VIGRA_CSTD::fabs(roots[i]) < 1.0) 00204 b.push_back(roots[i]); 00205 } 00206 return b; 00207 } 00208 00209 template <int ORDER, class T> 00210 typename BSplineBase<ORDER, T>::WeightMatrix & 00211 BSplineBase<ORDER, T>::calculateWeightMatrix() 00212 { 00213 static WeightMatrix b; 00214 double faculty = 1.0; 00215 for(int d = 0; d <= ORDER; ++d) 00216 { 00217 if(d > 1) 00218 faculty *= d; 00219 double x = ORDER / 2; 00220 BSplineBase spline; 00221 for(int i = 0; i <= ORDER; ++i, --x) 00222 b[d][i] = spline(x, d) / faculty; 00223 } 00224 return b; 00225 } 00226 00227 /********************************************************/ 00228 /* */ 00229 /* BSpline<N, T> */ 00230 /* */ 00231 /********************************************************/ 00232 00233 /** Spline functors for arbitrary orders. 00234 00235 Provides the interface of \ref vigra::BSplineBase with a more convenient 00236 name -- see there for more documentation. 00237 */ 00238 template <int ORDER, class T = double> 00239 class BSpline 00240 : public BSplineBase<ORDER, T> 00241 { 00242 public: 00243 /** Constructor forwarded to the base class constructor.. 00244 */ 00245 explicit BSpline(unsigned int derivativeOrder = 0) 00246 : BSplineBase<ORDER, T>(derivativeOrder) 00247 {} 00248 }; 00249 00250 /********************************************************/ 00251 /* */ 00252 /* BSpline<0, T> */ 00253 /* */ 00254 /********************************************************/ 00255 00256 template <class T> 00257 class BSplineBase<0, T> 00258 { 00259 public: 00260 00261 typedef T value_type; 00262 typedef T argument_type; 00263 typedef T first_argument_type; 00264 typedef unsigned int second_argument_type; 00265 typedef T result_type; 00266 enum StaticOrder { order = 0 }; 00267 00268 explicit BSplineBase(unsigned int derivativeOrder = 0) 00269 : derivativeOrder_(derivativeOrder) 00270 {} 00271 00272 result_type operator()(argument_type x) const 00273 { 00274 return exec(x, derivativeOrder_); 00275 } 00276 00277 template <unsigned int IntBits, unsigned int FracBits> 00278 FixedPoint<IntBits, FracBits> operator()(FixedPoint<IntBits, FracBits> x) const 00279 { 00280 typedef FixedPoint<IntBits, FracBits> Value; 00281 return x.value < Value::ONE_HALF && -Value::ONE_HALF <= x.value 00282 ? Value(Value::ONE, FPNoShift) 00283 : Value(0, FPNoShift); 00284 } 00285 00286 result_type operator()(first_argument_type x, second_argument_type derivative_order) const 00287 { 00288 return exec(x, derivativeOrder_ + derivative_order); 00289 } 00290 00291 value_type operator[](value_type x) const 00292 { return operator()(x); } 00293 00294 double radius() const 00295 { return 0.5; } 00296 00297 unsigned int derivativeOrder() const 00298 { return derivativeOrder_; } 00299 00300 ArrayVector<double> const & prefilterCoefficients() const 00301 { 00302 static ArrayVector<double> b; 00303 return b; 00304 } 00305 00306 typedef T WeightMatrix[1][1]; 00307 static WeightMatrix & weights() 00308 { 00309 static T b[1][1] = {{ 1.0}}; 00310 return b; 00311 } 00312 00313 protected: 00314 result_type exec(first_argument_type x, second_argument_type derivative_order) const 00315 { 00316 if(derivative_order == 0) 00317 return x < 0.5 && -0.5 <= x ? 00318 1.0 00319 : 0.0; 00320 else 00321 return 0.0; 00322 } 00323 00324 unsigned int derivativeOrder_; 00325 }; 00326 00327 /********************************************************/ 00328 /* */ 00329 /* BSpline<1, T> */ 00330 /* */ 00331 /********************************************************/ 00332 00333 template <class T> 00334 class BSpline<1, T> 00335 { 00336 public: 00337 00338 typedef T value_type; 00339 typedef T argument_type; 00340 typedef T first_argument_type; 00341 typedef unsigned int second_argument_type; 00342 typedef T result_type; 00343 enum StaticOrder { order = 1 }; 00344 00345 explicit BSpline(unsigned int derivativeOrder = 0) 00346 : derivativeOrder_(derivativeOrder) 00347 {} 00348 00349 result_type operator()(argument_type x) const 00350 { 00351 return exec(x, derivativeOrder_); 00352 } 00353 00354 template <unsigned int IntBits, unsigned int FracBits> 00355 FixedPoint<IntBits, FracBits> operator()(FixedPoint<IntBits, FracBits> x) const 00356 { 00357 typedef FixedPoint<IntBits, FracBits> Value; 00358 int v = abs(x.value); 00359 return v < Value::ONE ? 00360 Value(Value::ONE - v, FPNoShift) 00361 : Value(0, FPNoShift); 00362 } 00363 00364 result_type operator()(first_argument_type x, second_argument_type derivative_order) const 00365 { 00366 return exec(x, derivativeOrder_ + derivative_order); 00367 } 00368 00369 value_type operator[](value_type x) const 00370 { return operator()(x); } 00371 00372 double radius() const 00373 { return 1.0; } 00374 00375 unsigned int derivativeOrder() const 00376 { return derivativeOrder_; } 00377 00378 ArrayVector<double> const & prefilterCoefficients() const 00379 { 00380 static ArrayVector<double> b; 00381 return b; 00382 } 00383 00384 typedef T WeightMatrix[2][2]; 00385 static WeightMatrix & weights() 00386 { 00387 static T b[2][2] = {{ 1.0, 0.0}, {-1.0, 1.0}}; 00388 return b; 00389 } 00390 00391 protected: 00392 T exec(T x, unsigned int derivative_order) const; 00393 00394 unsigned int derivativeOrder_; 00395 }; 00396 00397 template <class T> 00398 T BSpline<1, T>::exec(T x, unsigned int derivative_order) const 00399 { 00400 switch(derivative_order) 00401 { 00402 case 0: 00403 { 00404 x = VIGRA_CSTD::fabs(x); 00405 return x < 1.0 ? 00406 1.0 - x 00407 : 0.0; 00408 } 00409 case 1: 00410 { 00411 return x < 0.0 ? 00412 -1.0 <= x ? 00413 1.0 00414 : 0.0 00415 : x < 1.0 ? 00416 -1.0 00417 : 0.0; 00418 } 00419 default: 00420 return 0.0; 00421 } 00422 } 00423 00424 /********************************************************/ 00425 /* */ 00426 /* BSpline<2, T> */ 00427 /* */ 00428 /********************************************************/ 00429 00430 template <class T> 00431 class BSpline<2, T> 00432 { 00433 public: 00434 00435 typedef T value_type; 00436 typedef T argument_type; 00437 typedef T first_argument_type; 00438 typedef unsigned int second_argument_type; 00439 typedef T result_type; 00440 enum StaticOrder { order = 2 }; 00441 00442 explicit BSpline(unsigned int derivativeOrder = 0) 00443 : derivativeOrder_(derivativeOrder) 00444 {} 00445 00446 result_type operator()(argument_type x) const 00447 { 00448 return exec(x, derivativeOrder_); 00449 } 00450 00451 template <unsigned int IntBits, unsigned int FracBits> 00452 FixedPoint<IntBits, FracBits> operator()(FixedPoint<IntBits, FracBits> x) const 00453 { 00454 typedef FixedPoint<IntBits, FracBits> Value; 00455 enum { ONE_HALF = Value::ONE_HALF, THREE_HALVES = ONE_HALF * 3, THREE_QUARTERS = THREE_HALVES / 2, 00456 PREMULTIPLY_SHIFT1 = FracBits <= 16 ? 0 : FracBits - 16, 00457 PREMULTIPLY_SHIFT2 = FracBits - 1 <= 16 ? 0 : FracBits - 17, 00458 POSTMULTIPLY_SHIFT1 = FracBits - 2*PREMULTIPLY_SHIFT1, 00459 POSTMULTIPLY_SHIFT2 = FracBits - 2*PREMULTIPLY_SHIFT2 }; 00460 int v = abs(x.value); 00461 return v == ONE_HALF 00462 ? Value(ONE_HALF, FPNoShift) 00463 : v <= ONE_HALF 00464 ? Value(THREE_QUARTERS - 00465 (int)(sq((unsigned)v >> PREMULTIPLY_SHIFT2) >> POSTMULTIPLY_SHIFT2), FPNoShift) 00466 : v < THREE_HALVES 00467 ? Value((int)(sq((unsigned)(THREE_HALVES-v) >> PREMULTIPLY_SHIFT1) >> (POSTMULTIPLY_SHIFT1 + 1)), FPNoShift) 00468 : Value(0, FPNoShift); 00469 } 00470 00471 result_type operator()(first_argument_type x, second_argument_type derivative_order) const 00472 { 00473 return exec(x, derivativeOrder_ + derivative_order); 00474 } 00475 00476 value_type operator[](value_type x) const 00477 { return operator()(x); } 00478 00479 double radius() const 00480 { return 1.5; } 00481 00482 unsigned int derivativeOrder() const 00483 { return derivativeOrder_; } 00484 00485 ArrayVector<double> const & prefilterCoefficients() const 00486 { 00487 static ArrayVector<double> b(1, 2.0*M_SQRT2 - 3.0); 00488 return b; 00489 } 00490 00491 typedef T WeightMatrix[3][3]; 00492 static WeightMatrix & weights() 00493 { 00494 static T b[3][3] = {{ 0.125, 0.75, 0.125}, 00495 {-0.5, 0.0, 0.5}, 00496 { 0.5, -1.0, 0.5}}; 00497 return b; 00498 } 00499 00500 protected: 00501 result_type exec(first_argument_type x, second_argument_type derivative_order) const; 00502 00503 unsigned int derivativeOrder_; 00504 }; 00505 00506 template <class T> 00507 typename BSpline<2, T>::result_type 00508 BSpline<2, T>::exec(first_argument_type x, second_argument_type derivative_order) const 00509 { 00510 switch(derivative_order) 00511 { 00512 case 0: 00513 { 00514 x = VIGRA_CSTD::fabs(x); 00515 return x < 0.5 ? 00516 0.75 - x*x 00517 : x < 1.5 ? 00518 0.5 * sq(1.5 - x) 00519 : 0.0; 00520 } 00521 case 1: 00522 { 00523 return x >= -0.5 ? 00524 x <= 0.5 ? 00525 -2.0 * x 00526 : x < 1.5 ? 00527 x - 1.5 00528 : 0.0 00529 : x > -1.5 ? 00530 x + 1.5 00531 : 0.0; 00532 } 00533 case 2: 00534 { 00535 return x >= -0.5 ? 00536 x < 0.5 ? 00537 -2.0 00538 : x < 1.5 ? 00539 1.0 00540 : 0.0 00541 : x >= -1.5 ? 00542 1.0 00543 : 0.0; 00544 } 00545 default: 00546 return 0.0; 00547 } 00548 } 00549 00550 /********************************************************/ 00551 /* */ 00552 /* BSpline<3, T> */ 00553 /* */ 00554 /********************************************************/ 00555 00556 template <class T> 00557 class BSpline<3, T> 00558 { 00559 public: 00560 00561 typedef T value_type; 00562 typedef T argument_type; 00563 typedef T first_argument_type; 00564 typedef unsigned int second_argument_type; 00565 typedef T result_type; 00566 enum StaticOrder { order = 3 }; 00567 00568 explicit BSpline(unsigned int derivativeOrder = 0) 00569 : derivativeOrder_(derivativeOrder) 00570 {} 00571 00572 result_type operator()(argument_type x) const 00573 { 00574 return exec(x, derivativeOrder_); 00575 } 00576 00577 template <unsigned int IntBits, unsigned int FracBits> 00578 FixedPoint<IntBits, FracBits> operator()(FixedPoint<IntBits, FracBits> x) const 00579 { 00580 typedef FixedPoint<IntBits, FracBits> Value; 00581 enum { ONE = Value::ONE, TWO = 2 * ONE, TWO_THIRDS = TWO / 3, ONE_SIXTH = ONE / 6, 00582 PREMULTIPLY_SHIFT = FracBits <= 16 ? 0 : FracBits - 16, 00583 POSTMULTIPLY_SHIFT = FracBits - 2*PREMULTIPLY_SHIFT }; 00584 int v = abs(x.value); 00585 return v == ONE 00586 ? Value(ONE_SIXTH, FPNoShift) 00587 : v < ONE 00588 ? Value(TWO_THIRDS + 00589 (((int)(sq((unsigned)v >> PREMULTIPLY_SHIFT) >> (POSTMULTIPLY_SHIFT + PREMULTIPLY_SHIFT)) 00590 * (((v >> 1) - ONE) >> PREMULTIPLY_SHIFT)) >> POSTMULTIPLY_SHIFT), FPNoShift) 00591 : v < TWO 00592 ? Value((int)((sq((unsigned)(TWO-v) >> PREMULTIPLY_SHIFT) >> (POSTMULTIPLY_SHIFT + PREMULTIPLY_SHIFT)) 00593 * ((unsigned)(TWO-v) >> PREMULTIPLY_SHIFT) / 6) >> POSTMULTIPLY_SHIFT, FPNoShift) 00594 : Value(0, FPNoShift); 00595 } 00596 00597 result_type operator()(first_argument_type x, second_argument_type derivative_order) const 00598 { 00599 return exec(x, derivativeOrder_ + derivative_order); 00600 } 00601 00602 result_type dx(argument_type x) const 00603 { return operator()(x, 1); } 00604 00605 result_type dxx(argument_type x) const 00606 { return operator()(x, 2); } 00607 00608 value_type operator[](value_type x) const 00609 { return operator()(x); } 00610 00611 double radius() const 00612 { return 2.0; } 00613 00614 unsigned int derivativeOrder() const 00615 { return derivativeOrder_; } 00616 00617 ArrayVector<double> const & prefilterCoefficients() const 00618 { 00619 static ArrayVector<double> b(1, VIGRA_CSTD::sqrt(3.0) - 2.0); 00620 return b; 00621 } 00622 00623 typedef T WeightMatrix[4][4]; 00624 static WeightMatrix & weights() 00625 { 00626 static T b[4][4] = {{ 1.0 / 6.0, 2.0 / 3.0, 1.0 / 6.0, 0.0}, 00627 {-0.5, 0.0, 0.5, 0.0}, 00628 { 0.5, -1.0, 0.5, 0.0}, 00629 {-1.0 / 6.0, 0.5, -0.5, 1.0 / 6.0}}; 00630 return b; 00631 } 00632 00633 protected: 00634 result_type exec(first_argument_type x, second_argument_type derivative_order) const; 00635 00636 unsigned int derivativeOrder_; 00637 }; 00638 00639 template <class T> 00640 typename BSpline<3, T>::result_type 00641 BSpline<3, T>::exec(first_argument_type x, second_argument_type derivative_order) const 00642 { 00643 switch(derivative_order) 00644 { 00645 case 0: 00646 { 00647 x = VIGRA_CSTD::fabs(x); 00648 if(x < 1.0) 00649 { 00650 return 2.0/3.0 + x*x*(-1.0 + 0.5*x); 00651 } 00652 else if(x < 2.0) 00653 { 00654 x = 2.0 - x; 00655 return x*x*x/6.0; 00656 } 00657 else 00658 return 0.0; 00659 } 00660 case 1: 00661 { 00662 double s = x < 0.0 ? 00663 -1.0 00664 : 1.0; 00665 x = VIGRA_CSTD::fabs(x); 00666 return x < 1.0 ? 00667 s*x*(-2.0 + 1.5*x) 00668 : x < 2.0 ? 00669 -0.5*s*sq(2.0 - x) 00670 : 0.0; 00671 } 00672 case 2: 00673 { 00674 x = VIGRA_CSTD::fabs(x); 00675 return x < 1.0 ? 00676 3.0*x - 2.0 00677 : x < 2.0 ? 00678 2.0 - x 00679 : 0.0; 00680 } 00681 case 3: 00682 { 00683 return x < 0.0 ? 00684 x < -1.0 ? 00685 x < -2.0 ? 00686 0.0 00687 : 1.0 00688 : -3.0 00689 : x < 1.0 ? 00690 3.0 00691 : x < 2.0 ? 00692 -1.0 00693 : 0.0; 00694 } 00695 default: 00696 return 0.0; 00697 } 00698 } 00699 00700 typedef BSpline<3, double> CubicBSplineKernel; 00701 00702 /********************************************************/ 00703 /* */ 00704 /* BSpline<5, T> */ 00705 /* */ 00706 /********************************************************/ 00707 00708 template <class T> 00709 class BSpline<5, T> 00710 { 00711 public: 00712 00713 typedef T value_type; 00714 typedef T argument_type; 00715 typedef T first_argument_type; 00716 typedef unsigned int second_argument_type; 00717 typedef T result_type; 00718 enum StaticOrder { order = 5 }; 00719 00720 explicit BSpline(unsigned int derivativeOrder = 0) 00721 : derivativeOrder_(derivativeOrder) 00722 {} 00723 00724 result_type operator()(argument_type x) const 00725 { 00726 return exec(x, derivativeOrder_); 00727 } 00728 00729 result_type operator()(first_argument_type x, second_argument_type derivative_order) const 00730 { 00731 return exec(x, derivativeOrder_ + derivative_order); 00732 } 00733 00734 result_type dx(argument_type x) const 00735 { return operator()(x, 1); } 00736 00737 result_type dxx(argument_type x) const 00738 { return operator()(x, 2); } 00739 00740 result_type dx3(argument_type x) const 00741 { return operator()(x, 3); } 00742 00743 result_type dx4(argument_type x) const 00744 { return operator()(x, 4); } 00745 00746 value_type operator[](value_type x) const 00747 { return operator()(x); } 00748 00749 double radius() const 00750 { return 3.0; } 00751 00752 unsigned int derivativeOrder() const 00753 { return derivativeOrder_; } 00754 00755 ArrayVector<double> const & prefilterCoefficients() const 00756 { 00757 static ArrayVector<double> const & b = initPrefilterCoefficients(); 00758 return b; 00759 } 00760 00761 static ArrayVector<double> const & initPrefilterCoefficients() 00762 { 00763 static ArrayVector<double> b(2); 00764 b[0] = -0.43057534709997114; 00765 b[1] = -0.043096288203264652; 00766 return b; 00767 } 00768 00769 typedef T WeightMatrix[6][6]; 00770 static WeightMatrix & weights() 00771 { 00772 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}, 00773 {-1.0/24.0, -5.0/12.0, 0.0, 5.0/12.0, 1.0/24.0, 0.0}, 00774 { 1.0/12.0, 1.0/6.0, -0.5, 1.0/6.0, 1.0/12.0, 0.0}, 00775 {-1.0/12.0, 1.0/6.0, 0.0, -1.0/6.0, 1.0/12.0, 0.0}, 00776 { 1.0/24.0, -1.0/6.0, 0.25, -1.0/6.0, 1.0/24.0, 0.0}, 00777 {-1.0/120.0, 1.0/24.0, -1.0/12.0, 1.0/12.0, -1.0/24.0, 1.0/120.0}}; 00778 return b; 00779 } 00780 00781 protected: 00782 result_type exec(first_argument_type x, second_argument_type derivative_order) const; 00783 00784 unsigned int derivativeOrder_; 00785 }; 00786 00787 template <class T> 00788 typename BSpline<5, T>::result_type 00789 BSpline<5, T>::exec(first_argument_type x, second_argument_type derivative_order) const 00790 { 00791 switch(derivative_order) 00792 { 00793 case 0: 00794 { 00795 x = VIGRA_CSTD::fabs(x); 00796 if(x <= 1.0) 00797 { 00798 return 0.55 + x*x*(-0.5 + x*x*(0.25 - x/12.0)); 00799 } 00800 else if(x < 2.0) 00801 { 00802 return 17.0/40.0 + x*(0.625 + x*(-1.75 + x*(1.25 + x*(-0.375 + x/24.0)))); 00803 } 00804 else if(x < 3.0) 00805 { 00806 x = 3.0 - x; 00807 return x*sq(x*x) / 120.0; 00808 } 00809 else 00810 return 0.0; 00811 } 00812 case 1: 00813 { 00814 double s = x < 0.0 ? 00815 -1.0 : 00816 1.0; 00817 x = VIGRA_CSTD::fabs(x); 00818 if(x <= 1.0) 00819 { 00820 return s*x*(-1.0 + x*x*(1.0 - 5.0/12.0*x)); 00821 } 00822 else if(x < 2.0) 00823 { 00824 return s*(0.625 + x*(-3.5 + x*(3.75 + x*(-1.5 + 5.0/24.0*x)))); 00825 } 00826 else if(x < 3.0) 00827 { 00828 x = 3.0 - x; 00829 return s*sq(x*x) / -24.0; 00830 } 00831 else 00832 return 0.0; 00833 } 00834 case 2: 00835 { 00836 x = VIGRA_CSTD::fabs(x); 00837 if(x <= 1.0) 00838 { 00839 return -1.0 + x*x*(3.0 -5.0/3.0*x); 00840 } 00841 else if(x < 2.0) 00842 { 00843 return -3.5 + x*(7.5 + x*(-4.5 + 5.0/6.0*x)); 00844 } 00845 else if(x < 3.0) 00846 { 00847 x = 3.0 - x; 00848 return x*x*x / 6.0; 00849 } 00850 else 00851 return 0.0; 00852 } 00853 case 3: 00854 { 00855 double s = x < 0.0 ? 00856 -1.0 : 00857 1.0; 00858 x = VIGRA_CSTD::fabs(x); 00859 if(x <= 1.0) 00860 { 00861 return s*x*(6.0 - 5.0*x); 00862 } 00863 else if(x < 2.0) 00864 { 00865 return s*(7.5 + x*(-9.0 + 2.5*x)); 00866 } 00867 else if(x < 3.0) 00868 { 00869 x = 3.0 - x; 00870 return -0.5*s*x*x; 00871 } 00872 else 00873 return 0.0; 00874 } 00875 case 4: 00876 { 00877 x = VIGRA_CSTD::fabs(x); 00878 if(x <= 1.0) 00879 { 00880 return 6.0 - 10.0*x; 00881 } 00882 else if(x < 2.0) 00883 { 00884 return -9.0 + 5.0*x; 00885 } 00886 else if(x < 3.0) 00887 { 00888 return 3.0 - x; 00889 } 00890 else 00891 return 0.0; 00892 } 00893 case 5: 00894 { 00895 return x < 0.0 ? 00896 x < -2.0 ? 00897 x < -3.0 ? 00898 0.0 00899 : 1.0 00900 : x < -1.0 ? 00901 -5.0 00902 : 10.0 00903 : x < 2.0 ? 00904 x < 1.0 ? 00905 -10.0 00906 : 5.0 00907 : x < 3.0 ? 00908 -1.0 00909 : 0.0; 00910 } 00911 default: 00912 return 0.0; 00913 } 00914 } 00915 00916 typedef BSpline<5, double> QuinticBSplineKernel; 00917 00918 #endif // NO_PARTIAL_TEMPLATE_SPECIALIZATION 00919 00920 /********************************************************/ 00921 /* */ 00922 /* CatmullRomSpline */ 00923 /* */ 00924 /********************************************************/ 00925 00926 /** Interpolating 3-rd order splines. 00927 00928 Implements the Catmull/Rom cardinal function 00929 00930 \f[ f(x) = \left\{ \begin{array}{ll} 00931 \frac{3}{2}x^3 - \frac{5}{2}x^2 + 1 & |x| \leq 1 \\ 00932 -\frac{1}{2}x^3 + \frac{5}{2}x^2 -4x + 2 & |x| \leq 2 \\ 00933 0 & \mbox{otherwise} 00934 \end{array}\right. 00935 \f] 00936 00937 It can be used as a functor, and as a kernel for 00938 \ref resamplingConvolveImage() to create a differentiable interpolant 00939 of an image. However, it should be noted that a twice differentiable 00940 interpolant can be created with only slightly more effort by recursive 00941 prefiltering followed by convolution with a 3rd order B-spline. 00942 00943 <b>\#include</b> "<a href="splines_8hxx-source.html">vigra/splines.hxx</a>"<br> 00944 Namespace: vigra 00945 */ 00946 template <class T = double> 00947 class CatmullRomSpline 00948 { 00949 public: 00950 /** the kernel's value type 00951 */ 00952 typedef T value_type; 00953 /** the unary functor's argument type 00954 */ 00955 typedef T argument_type; 00956 /** the unary functor's result type 00957 */ 00958 typedef T result_type; 00959 /** the splines polynomial order 00960 */ 00961 enum StaticOrder { order = 3 }; 00962 00963 /** function (functor) call 00964 */ 00965 result_type operator()(argument_type x) const; 00966 00967 /** index operator -- same as operator() 00968 */ 00969 T operator[] (T x) const 00970 { return operator()(x); } 00971 00972 /** Radius of the function's support. 00973 Needed for \ref resamplingConvolveImage(), always 2. 00974 */ 00975 int radius() const 00976 {return 2;} 00977 00978 /** Derivative order of the function: always 0. 00979 */ 00980 unsigned int derivativeOrder() const 00981 { return 0; } 00982 00983 /** Prefilter coefficients for compatibility with \ref vigra::BSpline. 00984 (array has zero length, since prefiltering is not necessary). 00985 */ 00986 ArrayVector<double> const & prefilterCoefficients() const 00987 { 00988 static ArrayVector<double> b; 00989 return b; 00990 } 00991 }; 00992 00993 00994 template <class T> 00995 typename CatmullRomSpline<T>::result_type 00996 CatmullRomSpline<T>::operator()(argument_type x) const 00997 { 00998 x = VIGRA_CSTD::fabs(x); 00999 if (x <= 1.0) 01000 { 01001 return 1.0 + x * x * (-2.5 + 1.5 * x); 01002 } 01003 else if (x >= 2.0) 01004 { 01005 return 0.0; 01006 } 01007 else 01008 { 01009 return 2.0 + x * (-4.0 + x * (2.5 - 0.5 * x)); 01010 } 01011 } 01012 01013 01014 //@} 01015 01016 } // namespace vigra 01017 01018 01019 #endif /* VIGRA_SPLINES_HXX */
© Ullrich Köthe (koethe@informatik.uni-hamburg.de) |
html generated using doxygen and Python
|