[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]
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.6.0, Aug 13 2008 ) */ 00008 /* The VIGRA Website is */ 00009 /* http://kogs-www.informatik.uni-hamburg.de/~koethe/vigra/ */ 00010 /* Please direct questions, bug reports, and contributions to */ 00011 /* ullrich.koethe@iwr.uni-heidelberg.de or */ 00012 /* vigra@informatik.uni-hamburg.de */ 00013 /* */ 00014 /* Permission is hereby granted, free of charge, to any person */ 00015 /* obtaining a copy of this software and associated documentation */ 00016 /* files (the "Software"), to deal in the Software without */ 00017 /* restriction, including without limitation the rights to use, */ 00018 /* copy, modify, merge, publish, distribute, sublicense, and/or */ 00019 /* sell copies of the Software, and to permit persons to whom the */ 00020 /* Software is furnished to do so, subject to the following */ 00021 /* conditions: */ 00022 /* */ 00023 /* The above copyright notice and this permission notice shall be */ 00024 /* included in all copies or substantial portions of the */ 00025 /* Software. */ 00026 /* */ 00027 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */ 00028 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */ 00029 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */ 00030 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */ 00031 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */ 00032 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */ 00033 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */ 00034 /* OTHER DEALINGS IN THE SOFTWARE. */ 00035 /* */ 00036 /************************************************************************/ 00037 00038 #ifndef VIGRA_SPLINES_HXX 00039 #define VIGRA_SPLINES_HXX 00040 00041 #include <cmath> 00042 #include "config.hxx" 00043 #include "mathutil.hxx" 00044 #include "polynomial.hxx" 00045 #include "array_vector.hxx" 00046 #include "fixedpoint.hxx" 00047 00048 namespace vigra { 00049 00050 /** \addtogroup MathFunctions Mathematical Functions 00051 */ 00052 //@{ 00053 /* B-Splines of arbitrary order and interpolating Catmull/Rom splines. 00054 00055 <b>\#include</b> <<a href="splines_8hxx-source.html">vigra/splines.hxx</a>><br> 00056 Namespace: vigra 00057 */ 00058 #ifndef NO_PARTIAL_TEMPLATE_SPECIALIZATION 00059 00060 /** Basic interface of the spline functors. 00061 00062 Implements the spline functions defined by the recursion 00063 00064 \f[ B_0(x) = \left\{ \begin{array}{ll} 00065 1 & -\frac{1}{2} \leq x < \frac{1}{2} \\ 00066 0 & \mbox{otherwise} 00067 \end{array}\right. 00068 \f] 00069 00070 and 00071 00072 \f[ B_n(x) = B_0(x) * B_{n-1}(x) 00073 \f] 00074 00075 where * denotes convolution, and <i>n</i> is the spline order given by the 00076 template parameter <tt>ORDER</tt>. These spline classes can be used as 00077 unary and binary functors, as kernels for \ref resamplingConvolveImage(), 00078 and as arguments for \ref vigra::SplineImageView. Note that the spline order 00079 is given as a template argument. 00080 00081 <b>\#include</b> <<a href="splines_8hxx-source.html">vigra/splines.hxx</a>><br> 00082 Namespace: vigra 00083 */ 00084 template <int ORDER, class T = double> 00085 class BSplineBase 00086 { 00087 public: 00088 00089 /** the value type if used as a kernel in \ref resamplingConvolveImage(). 00090 */ 00091 typedef T value_type; 00092 /** the functor's unary argument type 00093 */ 00094 typedef T argument_type; 00095 /** the functor's first binary argument type 00096 */ 00097 typedef T first_argument_type; 00098 /** the functor's second binary argument type 00099 */ 00100 typedef unsigned int second_argument_type; 00101 /** the functor's result type (unary and binary) 00102 */ 00103 typedef T result_type; 00104 /** the spline order 00105 */ 00106 enum StaticOrder { order = ORDER }; 00107 00108 /** Create functor for gevine derivative of the spline. The spline's order 00109 is specified spline by the template argument <TT>ORDER</tt>. 00110 */ 00111 explicit BSplineBase(unsigned int derivativeOrder = 0) 00112 : s1_(derivativeOrder) 00113 {} 00114 00115 /** Unary function call. 00116 Returns the value of the spline with the derivative order given in the 00117 constructor. Note that only derivatives up to <tt>ORDER-1</tt> are 00118 continous, and derivatives above <tt>ORDER+1</tt> are zero. 00119 */ 00120 result_type operator()(argument_type x) const 00121 { 00122 return exec(x, derivativeOrder()); 00123 } 00124 00125 /** Binary function call. 00126 The given derivative order is added to the derivative order 00127 specified in the constructor. Note that only derivatives up to <tt>ORDER-1</tt> are 00128 continous, and derivatives above <tt>ORDER+1</tt> are zero. 00129 */ 00130 result_type operator()(first_argument_type x, second_argument_type derivative_order) const 00131 { 00132 return exec(x, derivativeOrder() + derivative_order); 00133 } 00134 00135 /** Index operator. Same as unary function call. 00136 */ 00137 value_type operator[](value_type x) const 00138 { return operator()(x); } 00139 00140 /** Get the required filter radius for a discrete approximation of the 00141 spline. Always equal to <tt>(ORDER + 1) / 2.0</tt>. 00142 */ 00143 double radius() const 00144 { return (ORDER + 1) * 0.5; } 00145 00146 /** Get the derivative order of the Gaussian. 00147 */ 00148 unsigned int derivativeOrder() const 00149 { return s1_.derivativeOrder(); } 00150 00151 /** Get the prefilter coefficients required for interpolation. 00152 To interpolate with a B-spline, \ref resamplingConvolveImage() 00153 can be used. However, the image to be interpolated must be 00154 pre-filtered using \ref recursiveFilterX() and \ref recursiveFilterY() 00155 with the filter coefficients given by this function. The length of the array 00156 corresponds to how many times the above recursive filtering 00157 has to be applied (zero length means no prefiltering necessary). 00158 */ 00159 ArrayVector<double> const & prefilterCoefficients() const 00160 { 00161 static ArrayVector<double> const & b = calculatePrefilterCoefficients(); 00162 return b; 00163 } 00164 00165 static ArrayVector<double> const & calculatePrefilterCoefficients(); 00166 00167 typedef T WeightMatrix[ORDER+1][ORDER+1]; 00168 00169 /** Get the coefficients to transform spline coefficients into 00170 the coefficients of the corresponding polynomial. 00171 Currently internally used in SplineImageView; needs more 00172 documentation ??? 00173 */ 00174 static WeightMatrix & weights() 00175 { 00176 static WeightMatrix & b = calculateWeightMatrix(); 00177 return b; 00178 } 00179 00180 static WeightMatrix & calculateWeightMatrix(); 00181 00182 protected: 00183 result_type exec(first_argument_type x, second_argument_type derivative_order) const; 00184 00185 BSplineBase<ORDER-1, T> s1_; 00186 }; 00187 00188 template <int ORDER, class T> 00189 typename BSplineBase<ORDER, T>::result_type 00190 BSplineBase<ORDER, T>::exec(first_argument_type x, second_argument_type derivative_order) const 00191 { 00192 if(derivative_order == 0) 00193 { 00194 T n12 = (ORDER + 1.0) / 2.0; 00195 return ((n12 + x) * s1_(x + 0.5) + (n12 - x) * s1_(x - 0.5)) / ORDER; 00196 } 00197 else 00198 { 00199 --derivative_order; 00200 return s1_(x + 0.5, derivative_order) - s1_(x - 0.5, derivative_order); 00201 } 00202 } 00203 00204 template <int ORDER, class T> 00205 ArrayVector<double> const & BSplineBase<ORDER, T>::calculatePrefilterCoefficients() 00206 { 00207 static ArrayVector<double> b; 00208 if(ORDER > 1) 00209 { 00210 static const int r = ORDER / 2; 00211 StaticPolynomial<2*r, double> p(2*r); 00212 BSplineBase spline; 00213 for(int i = 0; i <= 2*r; ++i) 00214 p[i] = spline(T(i-r)); 00215 ArrayVector<double> roots; 00216 polynomialRealRoots(p, roots); 00217 for(unsigned int i = 0; i < roots.size(); ++i) 00218 if(VIGRA_CSTD::fabs(roots[i]) < 1.0) 00219 b.push_back(roots[i]); 00220 } 00221 return b; 00222 } 00223 00224 template <int ORDER, class T> 00225 typename BSplineBase<ORDER, T>::WeightMatrix & 00226 BSplineBase<ORDER, T>::calculateWeightMatrix() 00227 { 00228 static WeightMatrix b; 00229 double faculty = 1.0; 00230 for(int d = 0; d <= ORDER; ++d) 00231 { 00232 if(d > 1) 00233 faculty *= d; 00234 double x = ORDER / 2; // (note: integer division) 00235 BSplineBase spline; 00236 for(int i = 0; i <= ORDER; ++i, --x) 00237 b[d][i] = spline(x, d) / faculty; 00238 } 00239 return b; 00240 } 00241 00242 /********************************************************/ 00243 /* */ 00244 /* BSpline<N, T> */ 00245 /* */ 00246 /********************************************************/ 00247 00248 /** Spline functors for arbitrary orders. 00249 00250 Provides the interface of \ref vigra::BSplineBase with a more convenient 00251 name -- see there for more documentation. 00252 */ 00253 template <int ORDER, class T = double> 00254 class BSpline 00255 : public BSplineBase<ORDER, T> 00256 { 00257 public: 00258 /** Constructor forwarded to the base class constructor.. 00259 */ 00260 explicit BSpline(unsigned int derivativeOrder = 0) 00261 : BSplineBase<ORDER, T>(derivativeOrder) 00262 {} 00263 }; 00264 00265 /********************************************************/ 00266 /* */ 00267 /* BSpline<0, T> */ 00268 /* */ 00269 /********************************************************/ 00270 00271 template <class T> 00272 class BSplineBase<0, T> 00273 { 00274 public: 00275 00276 typedef T value_type; 00277 typedef T argument_type; 00278 typedef T first_argument_type; 00279 typedef unsigned int second_argument_type; 00280 typedef T result_type; 00281 enum StaticOrder { order = 0 }; 00282 00283 explicit BSplineBase(unsigned int derivativeOrder = 0) 00284 : derivativeOrder_(derivativeOrder) 00285 {} 00286 00287 result_type operator()(argument_type x) const 00288 { 00289 return exec(x, derivativeOrder_); 00290 } 00291 00292 template <unsigned int IntBits, unsigned int FracBits> 00293 FixedPoint<IntBits, FracBits> operator()(FixedPoint<IntBits, FracBits> x) const 00294 { 00295 typedef FixedPoint<IntBits, FracBits> Value; 00296 return x.value < Value::ONE_HALF && -Value::ONE_HALF <= x.value 00297 ? Value(Value::ONE, FPNoShift) 00298 : Value(0, FPNoShift); 00299 } 00300 00301 result_type operator()(first_argument_type x, second_argument_type derivative_order) const 00302 { 00303 return exec(x, derivativeOrder_ + derivative_order); 00304 } 00305 00306 value_type operator[](value_type x) const 00307 { return operator()(x); } 00308 00309 double radius() const 00310 { return 0.5; } 00311 00312 unsigned int derivativeOrder() const 00313 { return derivativeOrder_; } 00314 00315 ArrayVector<double> const & prefilterCoefficients() const 00316 { 00317 static ArrayVector<double> b; 00318 return b; 00319 } 00320 00321 typedef T WeightMatrix[1][1]; 00322 static WeightMatrix & weights() 00323 { 00324 static T b[1][1] = {{ 1.0}}; 00325 return b; 00326 } 00327 00328 protected: 00329 result_type exec(first_argument_type x, second_argument_type derivative_order) const 00330 { 00331 if(derivative_order == 0) 00332 return x < 0.5 && -0.5 <= x ? 00333 1.0 00334 : 0.0; 00335 else 00336 return 0.0; 00337 } 00338 00339 unsigned int derivativeOrder_; 00340 }; 00341 00342 /********************************************************/ 00343 /* */ 00344 /* BSpline<1, T> */ 00345 /* */ 00346 /********************************************************/ 00347 00348 template <class T> 00349 class BSpline<1, T> 00350 { 00351 public: 00352 00353 typedef T value_type; 00354 typedef T argument_type; 00355 typedef T first_argument_type; 00356 typedef unsigned int second_argument_type; 00357 typedef T result_type; 00358 enum StaticOrder { order = 1 }; 00359 00360 explicit BSpline(unsigned int derivativeOrder = 0) 00361 : derivativeOrder_(derivativeOrder) 00362 {} 00363 00364 result_type operator()(argument_type x) const 00365 { 00366 return exec(x, derivativeOrder_); 00367 } 00368 00369 template <unsigned int IntBits, unsigned int FracBits> 00370 FixedPoint<IntBits, FracBits> operator()(FixedPoint<IntBits, FracBits> x) const 00371 { 00372 typedef FixedPoint<IntBits, FracBits> Value; 00373 int v = abs(x.value); 00374 return v < Value::ONE ? 00375 Value(Value::ONE - v, FPNoShift) 00376 : Value(0, FPNoShift); 00377 } 00378 00379 result_type operator()(first_argument_type x, second_argument_type derivative_order) const 00380 { 00381 return exec(x, derivativeOrder_ + derivative_order); 00382 } 00383 00384 value_type operator[](value_type x) const 00385 { return operator()(x); } 00386 00387 double radius() const 00388 { return 1.0; } 00389 00390 unsigned int derivativeOrder() const 00391 { return derivativeOrder_; } 00392 00393 ArrayVector<double> const & prefilterCoefficients() const 00394 { 00395 static ArrayVector<double> b; 00396 return b; 00397 } 00398 00399 typedef T WeightMatrix[2][2]; 00400 static WeightMatrix & weights() 00401 { 00402 static T b[2][2] = {{ 1.0, 0.0}, {-1.0, 1.0}}; 00403 return b; 00404 } 00405 00406 protected: 00407 T exec(T x, unsigned int derivative_order) const; 00408 00409 unsigned int derivativeOrder_; 00410 }; 00411 00412 template <class T> 00413 T BSpline<1, T>::exec(T x, unsigned int derivative_order) const 00414 { 00415 switch(derivative_order) 00416 { 00417 case 0: 00418 { 00419 x = VIGRA_CSTD::fabs(x); 00420 return x < 1.0 ? 00421 1.0 - x 00422 : 0.0; 00423 } 00424 case 1: 00425 { 00426 return x < 0.0 ? 00427 -1.0 <= x ? 00428 1.0 00429 : 0.0 00430 : x < 1.0 ? 00431 -1.0 00432 : 0.0; 00433 } 00434 default: 00435 return 0.0; 00436 } 00437 } 00438 00439 /********************************************************/ 00440 /* */ 00441 /* BSpline<2, T> */ 00442 /* */ 00443 /********************************************************/ 00444 00445 template <class T> 00446 class BSpline<2, T> 00447 { 00448 public: 00449 00450 typedef T value_type; 00451 typedef T argument_type; 00452 typedef T first_argument_type; 00453 typedef unsigned int second_argument_type; 00454 typedef T result_type; 00455 enum StaticOrder { order = 2 }; 00456 00457 explicit BSpline(unsigned int derivativeOrder = 0) 00458 : derivativeOrder_(derivativeOrder) 00459 {} 00460 00461 result_type operator()(argument_type x) const 00462 { 00463 return exec(x, derivativeOrder_); 00464 } 00465 00466 template <unsigned int IntBits, unsigned int FracBits> 00467 FixedPoint<IntBits, FracBits> operator()(FixedPoint<IntBits, FracBits> x) const 00468 { 00469 typedef FixedPoint<IntBits, FracBits> Value; 00470 enum { ONE_HALF = Value::ONE_HALF, THREE_HALVES = ONE_HALF * 3, THREE_QUARTERS = THREE_HALVES / 2, 00471 PREMULTIPLY_SHIFT1 = FracBits <= 16 ? 0 : FracBits - 16, 00472 PREMULTIPLY_SHIFT2 = FracBits - 1 <= 16 ? 0 : FracBits - 17, 00473 POSTMULTIPLY_SHIFT1 = FracBits - 2*PREMULTIPLY_SHIFT1, 00474 POSTMULTIPLY_SHIFT2 = FracBits - 2*PREMULTIPLY_SHIFT2 }; 00475 int v = abs(x.value); 00476 return v == ONE_HALF 00477 ? Value(ONE_HALF, FPNoShift) 00478 : v <= ONE_HALF 00479 ? Value(THREE_QUARTERS - 00480 (int)(sq((unsigned)v >> PREMULTIPLY_SHIFT2) >> POSTMULTIPLY_SHIFT2), FPNoShift) 00481 : v < THREE_HALVES 00482 ? Value((int)(sq((unsigned)(THREE_HALVES-v) >> PREMULTIPLY_SHIFT1) >> (POSTMULTIPLY_SHIFT1 + 1)), FPNoShift) 00483 : Value(0, FPNoShift); 00484 } 00485 00486 result_type operator()(first_argument_type x, second_argument_type derivative_order) const 00487 { 00488 return exec(x, derivativeOrder_ + derivative_order); 00489 } 00490 00491 result_type dx(argument_type x) const 00492 { return operator()(x, 1); } 00493 00494 value_type operator[](value_type x) const 00495 { return operator()(x); } 00496 00497 double radius() const 00498 { return 1.5; } 00499 00500 unsigned int derivativeOrder() const 00501 { return derivativeOrder_; } 00502 00503 ArrayVector<double> const & prefilterCoefficients() const 00504 { 00505 static ArrayVector<double> b(1, 2.0*M_SQRT2 - 3.0); 00506 return b; 00507 } 00508 00509 typedef T WeightMatrix[3][3]; 00510 static WeightMatrix & weights() 00511 { 00512 static T b[3][3] = {{ 0.125, 0.75, 0.125}, 00513 {-0.5, 0.0, 0.5}, 00514 { 0.5, -1.0, 0.5}}; 00515 return b; 00516 } 00517 00518 protected: 00519 result_type exec(first_argument_type x, second_argument_type derivative_order) const; 00520 00521 unsigned int derivativeOrder_; 00522 }; 00523 00524 template <class T> 00525 typename BSpline<2, T>::result_type 00526 BSpline<2, T>::exec(first_argument_type x, second_argument_type derivative_order) const 00527 { 00528 switch(derivative_order) 00529 { 00530 case 0: 00531 { 00532 x = VIGRA_CSTD::fabs(x); 00533 return x < 0.5 ? 00534 0.75 - x*x 00535 : x < 1.5 ? 00536 0.5 * sq(1.5 - x) 00537 : 0.0; 00538 } 00539 case 1: 00540 { 00541 return x >= -0.5 ? 00542 x <= 0.5 ? 00543 -2.0 * x 00544 : x < 1.5 ? 00545 x - 1.5 00546 : 0.0 00547 : x > -1.5 ? 00548 x + 1.5 00549 : 0.0; 00550 } 00551 case 2: 00552 { 00553 return x >= -0.5 ? 00554 x < 0.5 ? 00555 -2.0 00556 : x < 1.5 ? 00557 1.0 00558 : 0.0 00559 : x >= -1.5 ? 00560 1.0 00561 : 0.0; 00562 } 00563 default: 00564 return 0.0; 00565 } 00566 } 00567 00568 /********************************************************/ 00569 /* */ 00570 /* BSpline<3, T> */ 00571 /* */ 00572 /********************************************************/ 00573 00574 template <class T> 00575 class BSpline<3, T> 00576 { 00577 public: 00578 00579 typedef T value_type; 00580 typedef T argument_type; 00581 typedef T first_argument_type; 00582 typedef unsigned int second_argument_type; 00583 typedef T result_type; 00584 enum StaticOrder { order = 3 }; 00585 00586 explicit BSpline(unsigned int derivativeOrder = 0) 00587 : derivativeOrder_(derivativeOrder) 00588 {} 00589 00590 result_type operator()(argument_type x) const 00591 { 00592 return exec(x, derivativeOrder_); 00593 } 00594 00595 template <unsigned int IntBits, unsigned int FracBits> 00596 FixedPoint<IntBits, FracBits> operator()(FixedPoint<IntBits, FracBits> x) const 00597 { 00598 typedef FixedPoint<IntBits, FracBits> Value; 00599 enum { ONE = Value::ONE, TWO = 2 * ONE, TWO_THIRDS = TWO / 3, ONE_SIXTH = ONE / 6, 00600 PREMULTIPLY_SHIFT = FracBits <= 16 ? 0 : FracBits - 16, 00601 POSTMULTIPLY_SHIFT = FracBits - 2*PREMULTIPLY_SHIFT }; 00602 int v = abs(x.value); 00603 return v == ONE 00604 ? Value(ONE_SIXTH, FPNoShift) 00605 : v < ONE 00606 ? Value(TWO_THIRDS + 00607 (((int)(sq((unsigned)v >> PREMULTIPLY_SHIFT) >> (POSTMULTIPLY_SHIFT + PREMULTIPLY_SHIFT)) 00608 * (((v >> 1) - ONE) >> PREMULTIPLY_SHIFT)) >> POSTMULTIPLY_SHIFT), FPNoShift) 00609 : v < TWO 00610 ? Value((int)((sq((unsigned)(TWO-v) >> PREMULTIPLY_SHIFT) >> (POSTMULTIPLY_SHIFT + PREMULTIPLY_SHIFT)) 00611 * ((unsigned)(TWO-v) >> PREMULTIPLY_SHIFT) / 6) >> POSTMULTIPLY_SHIFT, FPNoShift) 00612 : Value(0, FPNoShift); 00613 } 00614 00615 result_type operator()(first_argument_type x, second_argument_type derivative_order) const 00616 { 00617 return exec(x, derivativeOrder_ + derivative_order); 00618 } 00619 00620 result_type dx(argument_type x) const 00621 { return operator()(x, 1); } 00622 00623 result_type dxx(argument_type x) const 00624 { return operator()(x, 2); } 00625 00626 value_type operator[](value_type x) const 00627 { return operator()(x); } 00628 00629 double radius() const 00630 { return 2.0; } 00631 00632 unsigned int derivativeOrder() const 00633 { return derivativeOrder_; } 00634 00635 ArrayVector<double> const & prefilterCoefficients() const 00636 { 00637 static ArrayVector<double> b(1, VIGRA_CSTD::sqrt(3.0) - 2.0); 00638 return b; 00639 } 00640 00641 typedef T WeightMatrix[4][4]; 00642 static WeightMatrix & weights() 00643 { 00644 static T b[4][4] = {{ 1.0 / 6.0, 2.0 / 3.0, 1.0 / 6.0, 0.0}, 00645 {-0.5, 0.0, 0.5, 0.0}, 00646 { 0.5, -1.0, 0.5, 0.0}, 00647 {-1.0 / 6.0, 0.5, -0.5, 1.0 / 6.0}}; 00648 return b; 00649 } 00650 00651 protected: 00652 result_type exec(first_argument_type x, second_argument_type derivative_order) const; 00653 00654 unsigned int derivativeOrder_; 00655 }; 00656 00657 template <class T> 00658 typename BSpline<3, T>::result_type 00659 BSpline<3, T>::exec(first_argument_type x, second_argument_type derivative_order) const 00660 { 00661 switch(derivative_order) 00662 { 00663 case 0: 00664 { 00665 x = VIGRA_CSTD::fabs(x); 00666 if(x < 1.0) 00667 { 00668 return 2.0/3.0 + x*x*(-1.0 + 0.5*x); 00669 } 00670 else if(x < 2.0) 00671 { 00672 x = 2.0 - x; 00673 return x*x*x/6.0; 00674 } 00675 else 00676 return 0.0; 00677 } 00678 case 1: 00679 { 00680 double s = x < 0.0 ? 00681 -1.0 00682 : 1.0; 00683 x = VIGRA_CSTD::fabs(x); 00684 return x < 1.0 ? 00685 s*x*(-2.0 + 1.5*x) 00686 : x < 2.0 ? 00687 -0.5*s*sq(2.0 - x) 00688 : 0.0; 00689 } 00690 case 2: 00691 { 00692 x = VIGRA_CSTD::fabs(x); 00693 return x < 1.0 ? 00694 3.0*x - 2.0 00695 : x < 2.0 ? 00696 2.0 - x 00697 : 0.0; 00698 } 00699 case 3: 00700 { 00701 return x < 0.0 ? 00702 x < -1.0 ? 00703 x < -2.0 ? 00704 0.0 00705 : 1.0 00706 : -3.0 00707 : x < 1.0 ? 00708 3.0 00709 : x < 2.0 ? 00710 -1.0 00711 : 0.0; 00712 } 00713 default: 00714 return 0.0; 00715 } 00716 } 00717 00718 typedef BSpline<3, double> CubicBSplineKernel; 00719 00720 /********************************************************/ 00721 /* */ 00722 /* BSpline<4, T> */ 00723 /* */ 00724 /********************************************************/ 00725 00726 template <class T> 00727 class BSpline<4, T> 00728 { 00729 public: 00730 00731 typedef T value_type; 00732 typedef T argument_type; 00733 typedef T first_argument_type; 00734 typedef unsigned int second_argument_type; 00735 typedef T result_type; 00736 enum StaticOrder { order = 4 }; 00737 00738 explicit BSpline(unsigned int derivativeOrder = 0) 00739 : derivativeOrder_(derivativeOrder) 00740 {} 00741 00742 result_type operator()(argument_type x) const 00743 { 00744 return exec(x, derivativeOrder_); 00745 } 00746 00747 result_type operator()(first_argument_type x, second_argument_type derivative_order) const 00748 { 00749 return exec(x, derivativeOrder_ + derivative_order); 00750 } 00751 00752 result_type dx(argument_type x) const 00753 { return operator()(x, 1); } 00754 00755 result_type dxx(argument_type x) const 00756 { return operator()(x, 2); } 00757 00758 result_type dx3(argument_type x) const 00759 { return operator()(x, 3); } 00760 00761 value_type operator[](value_type x) const 00762 { return operator()(x); } 00763 00764 double radius() const 00765 { return 2.5; } 00766 00767 unsigned int derivativeOrder() const 00768 { return derivativeOrder_; } 00769 00770 ArrayVector<double> const & prefilterCoefficients() const 00771 { 00772 static ArrayVector<double> const & b = initPrefilterCoefficients(); 00773 return b; 00774 } 00775 00776 static ArrayVector<double> const & initPrefilterCoefficients() 00777 { 00778 static ArrayVector<double> b(2); 00779 // -19 + 4*sqrt(19) + 2*sqrt(2*(83 - 19*sqrt(19))) 00780 b[0] = -0.361341225900220177092212841325; 00781 // -19 - 4*sqrt(19) + 2*sqrt(2*(83 + 19*sqrt(19))) 00782 b[1] = -0.013725429297339121360331226939; 00783 return b; 00784 } 00785 00786 typedef T WeightMatrix[5][5]; 00787 static WeightMatrix & weights() 00788 { 00789 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}, 00790 {-1.0/48.0, -11.0/24.0, 0.0, 11.0/24.0, 1.0/48.0}, 00791 { 1.0/16.0, 1.0/4.0, -5.0/8.0, 1.0/4.0, 1.0/16.0}, 00792 {-1.0/12.0, 1.0/6.0, 0.0, -1.0/6.0, 1.0/12.0}, 00793 { 1.0/24.0, -1.0/6.0, 0.25, -1.0/6.0, 1.0/24.0}}; 00794 return b; 00795 } 00796 00797 protected: 00798 result_type exec(first_argument_type x, second_argument_type derivative_order) const; 00799 00800 unsigned int derivativeOrder_; 00801 }; 00802 00803 template <class T> 00804 typename BSpline<4, T>::result_type 00805 BSpline<4, T>::exec(first_argument_type x, second_argument_type derivative_order) const 00806 { 00807 switch(derivative_order) 00808 { 00809 case 0: 00810 { 00811 x = VIGRA_CSTD::fabs(x); 00812 if(x <= 0.5) 00813 { 00814 return 115.0/192.0 + x*x*(-0.625 + x*x*0.25); 00815 } 00816 else if(x < 1.5) 00817 { 00818 return (55.0/16.0 + x*(1.25 + x*(-7.5 + x*(5.0 - x)))) / 6.0; 00819 } 00820 else if(x < 2.5) 00821 { 00822 x = 2.5 - x; 00823 return sq(x*x) / 24.0; 00824 } 00825 else 00826 return 0.0; 00827 } 00828 case 1: 00829 { 00830 double s = x < 0.0 ? 00831 -1.0 : 00832 1.0; 00833 x = VIGRA_CSTD::fabs(x); 00834 if(x <= 0.5) 00835 { 00836 return s*x*(-1.25 + x*x); 00837 } 00838 else if(x < 1.5) 00839 { 00840 return s*(5.0 + x*(-60.0 + x*(60.0 - 16.0*x))) / 24.0; 00841 } 00842 else if(x < 2.5) 00843 { 00844 x = 2.5 - x; 00845 return s*x*x*x / -6.0; 00846 } 00847 else 00848 return 0.0; 00849 } 00850 case 2: 00851 { 00852 x = VIGRA_CSTD::fabs(x); 00853 if(x <= 0.5) 00854 { 00855 return -1.25 + 3.0*x*x; 00856 } 00857 else if(x < 1.5) 00858 { 00859 return -2.5 + x*(5.0 - 2.0*x); 00860 } 00861 else if(x < 2.5) 00862 { 00863 x = 2.5 - x; 00864 return x*x / 2.0; 00865 } 00866 else 00867 return 0.0; 00868 } 00869 case 3: 00870 { 00871 double s = x < 0.0 ? 00872 -1.0 : 00873 1.0; 00874 x = VIGRA_CSTD::fabs(x); 00875 if(x <= 0.5) 00876 { 00877 return s*x*6.0; 00878 } 00879 else if(x < 1.5) 00880 { 00881 return s*(5.0 - 4.0*x); 00882 } 00883 else if(x < 2.5) 00884 { 00885 return s*(x - 2.5); 00886 } 00887 else 00888 return 0.0; 00889 } 00890 case 4: 00891 { 00892 return x < 0.0 00893 ? x < -2.5 00894 ? 0.0 00895 : x < -1.5 00896 ? 1.0 00897 : x < -0.5 00898 ? -4.0 00899 : 6.0 00900 : x < 0.5 00901 ? 6.0 00902 : x < 1.5 00903 ? -4.0 00904 : x < 2.5 00905 ? 1.0 00906 : 0.0; 00907 } 00908 default: 00909 return 0.0; 00910 } 00911 } 00912 00913 typedef BSpline<4, double> QuarticBSplineKernel; 00914 00915 /********************************************************/ 00916 /* */ 00917 /* BSpline<5, T> */ 00918 /* */ 00919 /********************************************************/ 00920 00921 template <class T> 00922 class BSpline<5, T> 00923 { 00924 public: 00925 00926 typedef T value_type; 00927 typedef T argument_type; 00928 typedef T first_argument_type; 00929 typedef unsigned int second_argument_type; 00930 typedef T result_type; 00931 enum StaticOrder { order = 5 }; 00932 00933 explicit BSpline(unsigned int derivativeOrder = 0) 00934 : derivativeOrder_(derivativeOrder) 00935 {} 00936 00937 result_type operator()(argument_type x) const 00938 { 00939 return exec(x, derivativeOrder_); 00940 } 00941 00942 result_type operator()(first_argument_type x, second_argument_type derivative_order) const 00943 { 00944 return exec(x, derivativeOrder_ + derivative_order); 00945 } 00946 00947 result_type dx(argument_type x) const 00948 { return operator()(x, 1); } 00949 00950 result_type dxx(argument_type x) const 00951 { return operator()(x, 2); } 00952 00953 result_type dx3(argument_type x) const 00954 { return operator()(x, 3); } 00955 00956 result_type dx4(argument_type x) const 00957 { return operator()(x, 4); } 00958 00959 value_type operator[](value_type x) const 00960 { return operator()(x); } 00961 00962 double radius() const 00963 { return 3.0; } 00964 00965 unsigned int derivativeOrder() const 00966 { return derivativeOrder_; } 00967 00968 ArrayVector<double> const & prefilterCoefficients() const 00969 { 00970 static ArrayVector<double> const & b = initPrefilterCoefficients(); 00971 return b; 00972 } 00973 00974 static ArrayVector<double> const & initPrefilterCoefficients() 00975 { 00976 static ArrayVector<double> b(2); 00977 // -(13/2) + sqrt(105)/2 + sqrt(1/2*((135 - 13*sqrt(105)))) 00978 b[0] = -0.430575347099973791851434783493; 00979 // (1/2)*((-13) - sqrt(105) + sqrt(2*((135 + 13*sqrt(105))))) 00980 b[1] = -0.043096288203264653822712376822; 00981 return b; 00982 } 00983 00984 typedef T WeightMatrix[6][6]; 00985 static WeightMatrix & weights() 00986 { 00987 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}, 00988 {-1.0/24.0, -5.0/12.0, 0.0, 5.0/12.0, 1.0/24.0, 0.0}, 00989 { 1.0/12.0, 1.0/6.0, -0.5, 1.0/6.0, 1.0/12.0, 0.0}, 00990 {-1.0/12.0, 1.0/6.0, 0.0, -1.0/6.0, 1.0/12.0, 0.0}, 00991 { 1.0/24.0, -1.0/6.0, 0.25, -1.0/6.0, 1.0/24.0, 0.0}, 00992 {-1.0/120.0, 1.0/24.0, -1.0/12.0, 1.0/12.0, -1.0/24.0, 1.0/120.0}}; 00993 return b; 00994 } 00995 00996 protected: 00997 result_type exec(first_argument_type x, second_argument_type derivative_order) const; 00998 00999 unsigned int derivativeOrder_; 01000 }; 01001 01002 template <class T> 01003 typename BSpline<5, T>::result_type 01004 BSpline<5, T>::exec(first_argument_type x, second_argument_type derivative_order) const 01005 { 01006 switch(derivative_order) 01007 { 01008 case 0: 01009 { 01010 x = VIGRA_CSTD::fabs(x); 01011 if(x <= 1.0) 01012 { 01013 return 0.55 + x*x*(-0.5 + x*x*(0.25 - x/12.0)); 01014 } 01015 else if(x < 2.0) 01016 { 01017 return 17.0/40.0 + x*(0.625 + x*(-1.75 + x*(1.25 + x*(-0.375 + x/24.0)))); 01018 } 01019 else if(x < 3.0) 01020 { 01021 x = 3.0 - x; 01022 return x*sq(x*x) / 120.0; 01023 } 01024 else 01025 return 0.0; 01026 } 01027 case 1: 01028 { 01029 double s = x < 0.0 ? 01030 -1.0 : 01031 1.0; 01032 x = VIGRA_CSTD::fabs(x); 01033 if(x <= 1.0) 01034 { 01035 return s*x*(-1.0 + x*x*(1.0 - 5.0/12.0*x)); 01036 } 01037 else if(x < 2.0) 01038 { 01039 return s*(0.625 + x*(-3.5 + x*(3.75 + x*(-1.5 + 5.0/24.0*x)))); 01040 } 01041 else if(x < 3.0) 01042 { 01043 x = 3.0 - x; 01044 return s*sq(x*x) / -24.0; 01045 } 01046 else 01047 return 0.0; 01048 } 01049 case 2: 01050 { 01051 x = VIGRA_CSTD::fabs(x); 01052 if(x <= 1.0) 01053 { 01054 return -1.0 + x*x*(3.0 -5.0/3.0*x); 01055 } 01056 else if(x < 2.0) 01057 { 01058 return -3.5 + x*(7.5 + x*(-4.5 + 5.0/6.0*x)); 01059 } 01060 else if(x < 3.0) 01061 { 01062 x = 3.0 - x; 01063 return x*x*x / 6.0; 01064 } 01065 else 01066 return 0.0; 01067 } 01068 case 3: 01069 { 01070 double s = x < 0.0 ? 01071 -1.0 : 01072 1.0; 01073 x = VIGRA_CSTD::fabs(x); 01074 if(x <= 1.0) 01075 { 01076 return s*x*(6.0 - 5.0*x); 01077 } 01078 else if(x < 2.0) 01079 { 01080 return s*(7.5 + x*(-9.0 + 2.5*x)); 01081 } 01082 else if(x < 3.0) 01083 { 01084 x = 3.0 - x; 01085 return -0.5*s*x*x; 01086 } 01087 else 01088 return 0.0; 01089 } 01090 case 4: 01091 { 01092 x = VIGRA_CSTD::fabs(x); 01093 if(x <= 1.0) 01094 { 01095 return 6.0 - 10.0*x; 01096 } 01097 else if(x < 2.0) 01098 { 01099 return -9.0 + 5.0*x; 01100 } 01101 else if(x < 3.0) 01102 { 01103 return 3.0 - x; 01104 } 01105 else 01106 return 0.0; 01107 } 01108 case 5: 01109 { 01110 return x < 0.0 ? 01111 x < -2.0 ? 01112 x < -3.0 ? 01113 0.0 01114 : 1.0 01115 : x < -1.0 ? 01116 -5.0 01117 : 10.0 01118 : x < 2.0 ? 01119 x < 1.0 ? 01120 -10.0 01121 : 5.0 01122 : x < 3.0 ? 01123 -1.0 01124 : 0.0; 01125 } 01126 default: 01127 return 0.0; 01128 } 01129 } 01130 01131 typedef BSpline<5, double> QuinticBSplineKernel; 01132 01133 #endif // NO_PARTIAL_TEMPLATE_SPECIALIZATION 01134 01135 /********************************************************/ 01136 /* */ 01137 /* CatmullRomSpline */ 01138 /* */ 01139 /********************************************************/ 01140 01141 /** Interpolating 3-rd order splines. 01142 01143 Implements the Catmull/Rom cardinal function 01144 01145 \f[ f(x) = \left\{ \begin{array}{ll} 01146 \frac{3}{2}x^3 - \frac{5}{2}x^2 + 1 & |x| \leq 1 \\ 01147 -\frac{1}{2}x^3 + \frac{5}{2}x^2 -4x + 2 & |x| \leq 2 \\ 01148 0 & \mbox{otherwise} 01149 \end{array}\right. 01150 \f] 01151 01152 It can be used as a functor, and as a kernel for 01153 \ref resamplingConvolveImage() to create a differentiable interpolant 01154 of an image. However, it should be noted that a twice differentiable 01155 interpolant can be created with only slightly more effort by recursive 01156 prefiltering followed by convolution with a 3rd order B-spline. 01157 01158 <b>\#include</b> <<a href="splines_8hxx-source.html">vigra/splines.hxx</a>><br> 01159 Namespace: vigra 01160 */ 01161 template <class T = double> 01162 class CatmullRomSpline 01163 { 01164 public: 01165 /** the kernel's value type 01166 */ 01167 typedef T value_type; 01168 /** the unary functor's argument type 01169 */ 01170 typedef T argument_type; 01171 /** the unary functor's result type 01172 */ 01173 typedef T result_type; 01174 /** the splines polynomial order 01175 */ 01176 enum StaticOrder { order = 3 }; 01177 01178 /** function (functor) call 01179 */ 01180 result_type operator()(argument_type x) const; 01181 01182 /** index operator -- same as operator() 01183 */ 01184 T operator[] (T x) const 01185 { return operator()(x); } 01186 01187 /** Radius of the function's support. 01188 Needed for \ref resamplingConvolveImage(), always 2. 01189 */ 01190 int radius() const 01191 {return 2;} 01192 01193 /** Derivative order of the function: always 0. 01194 */ 01195 unsigned int derivativeOrder() const 01196 { return 0; } 01197 01198 /** Prefilter coefficients for compatibility with \ref vigra::BSpline. 01199 (array has zero length, since prefiltering is not necessary). 01200 */ 01201 ArrayVector<double> const & prefilterCoefficients() const 01202 { 01203 static ArrayVector<double> b; 01204 return b; 01205 } 01206 }; 01207 01208 01209 template <class T> 01210 typename CatmullRomSpline<T>::result_type 01211 CatmullRomSpline<T>::operator()(argument_type x) const 01212 { 01213 x = VIGRA_CSTD::fabs(x); 01214 if (x <= 1.0) 01215 { 01216 return 1.0 + x * x * (-2.5 + 1.5 * x); 01217 } 01218 else if (x >= 2.0) 01219 { 01220 return 0.0; 01221 } 01222 else 01223 { 01224 return 2.0 + x * (-4.0 + x * (2.5 - 0.5 * x)); 01225 } 01226 } 01227 01228 01229 //@} 01230 01231 } // namespace vigra 01232 01233 01234 #endif /* VIGRA_SPLINES_HXX */
© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de) |
html generated using doxygen and Python
|