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