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