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