[ VIGRA Homepage | Class Index | Function Index | File Index | Main Page ]
![]() |
vigra/splines.hxx | ![]() |
---|
00001 /************************************************************************/ 00002 /* */ 00003 /* Copyright 1998-2004 by Ullrich Koethe */ 00004 /* Cognitive Systems Group, University of Hamburg, Germany */ 00005 /* */ 00006 /* This file is part of the VIGRA computer vision library. */ 00007 /* ( Version 1.5.0, Dec 07 2006 ) */ 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 /* koethe@informatik.uni-hamburg.de or */ 00012 /* vigra@kogs1.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 recursiveFilterImage() with the filter 00155 coefficients given by this function. The length of the array 00156 corresponds to the number of times \ref recursiveFilterImage() 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; 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 value_type operator[](value_type x) const 00492 { return operator()(x); } 00493 00494 double radius() const 00495 { return 1.5; } 00496 00497 unsigned int derivativeOrder() const 00498 { return derivativeOrder_; } 00499 00500 ArrayVector<double> const & prefilterCoefficients() const 00501 { 00502 static ArrayVector<double> b(1, 2.0*M_SQRT2 - 3.0); 00503 return b; 00504 } 00505 00506 typedef T WeightMatrix[3][3]; 00507 static WeightMatrix & weights() 00508 { 00509 static T b[3][3] = {{ 0.125, 0.75, 0.125}, 00510 {-0.5, 0.0, 0.5}, 00511 { 0.5, -1.0, 0.5}}; 00512 return b; 00513 } 00514 00515 protected: 00516 result_type exec(first_argument_type x, second_argument_type derivative_order) const; 00517 00518 unsigned int derivativeOrder_; 00519 }; 00520 00521 template <class T> 00522 typename BSpline<2, T>::result_type 00523 BSpline<2, T>::exec(first_argument_type x, second_argument_type derivative_order) const 00524 { 00525 switch(derivative_order) 00526 { 00527 case 0: 00528 { 00529 x = VIGRA_CSTD::fabs(x); 00530 return x < 0.5 ? 00531 0.75 - x*x 00532 : x < 1.5 ? 00533 0.5 * sq(1.5 - x) 00534 : 0.0; 00535 } 00536 case 1: 00537 { 00538 return x >= -0.5 ? 00539 x <= 0.5 ? 00540 -2.0 * x 00541 : x < 1.5 ? 00542 x - 1.5 00543 : 0.0 00544 : x > -1.5 ? 00545 x + 1.5 00546 : 0.0; 00547 } 00548 case 2: 00549 { 00550 return x >= -0.5 ? 00551 x < 0.5 ? 00552 -2.0 00553 : x < 1.5 ? 00554 1.0 00555 : 0.0 00556 : x >= -1.5 ? 00557 1.0 00558 : 0.0; 00559 } 00560 default: 00561 return 0.0; 00562 } 00563 } 00564 00565 /********************************************************/ 00566 /* */ 00567 /* BSpline<3, T> */ 00568 /* */ 00569 /********************************************************/ 00570 00571 template <class T> 00572 class BSpline<3, T> 00573 { 00574 public: 00575 00576 typedef T value_type; 00577 typedef T argument_type; 00578 typedef T first_argument_type; 00579 typedef unsigned int second_argument_type; 00580 typedef T result_type; 00581 enum StaticOrder { order = 3 }; 00582 00583 explicit BSpline(unsigned int derivativeOrder = 0) 00584 : derivativeOrder_(derivativeOrder) 00585 {} 00586 00587 result_type operator()(argument_type x) const 00588 { 00589 return exec(x, derivativeOrder_); 00590 } 00591 00592 template <unsigned int IntBits, unsigned int FracBits> 00593 FixedPoint<IntBits, FracBits> operator()(FixedPoint<IntBits, FracBits> x) const 00594 { 00595 typedef FixedPoint<IntBits, FracBits> Value; 00596 enum { ONE = Value::ONE, TWO = 2 * ONE, TWO_THIRDS = TWO / 3, ONE_SIXTH = ONE / 6, 00597 PREMULTIPLY_SHIFT = FracBits <= 16 ? 0 : FracBits - 16, 00598 POSTMULTIPLY_SHIFT = FracBits - 2*PREMULTIPLY_SHIFT }; 00599 int v = abs(x.value); 00600 return v == ONE 00601 ? Value(ONE_SIXTH, FPNoShift) 00602 : v < ONE 00603 ? Value(TWO_THIRDS + 00604 (((int)(sq((unsigned)v >> PREMULTIPLY_SHIFT) >> (POSTMULTIPLY_SHIFT + PREMULTIPLY_SHIFT)) 00605 * (((v >> 1) - ONE) >> PREMULTIPLY_SHIFT)) >> POSTMULTIPLY_SHIFT), FPNoShift) 00606 : v < TWO 00607 ? Value((int)((sq((unsigned)(TWO-v) >> PREMULTIPLY_SHIFT) >> (POSTMULTIPLY_SHIFT + PREMULTIPLY_SHIFT)) 00608 * ((unsigned)(TWO-v) >> PREMULTIPLY_SHIFT) / 6) >> POSTMULTIPLY_SHIFT, FPNoShift) 00609 : Value(0, FPNoShift); 00610 } 00611 00612 result_type operator()(first_argument_type x, second_argument_type derivative_order) const 00613 { 00614 return exec(x, derivativeOrder_ + derivative_order); 00615 } 00616 00617 result_type dx(argument_type x) const 00618 { return operator()(x, 1); } 00619 00620 result_type dxx(argument_type x) const 00621 { return operator()(x, 2); } 00622 00623 value_type operator[](value_type x) const 00624 { return operator()(x); } 00625 00626 double radius() const 00627 { return 2.0; } 00628 00629 unsigned int derivativeOrder() const 00630 { return derivativeOrder_; } 00631 00632 ArrayVector<double> const & prefilterCoefficients() const 00633 { 00634 static ArrayVector<double> b(1, VIGRA_CSTD::sqrt(3.0) - 2.0); 00635 return b; 00636 } 00637 00638 typedef T WeightMatrix[4][4]; 00639 static WeightMatrix & weights() 00640 { 00641 static T b[4][4] = {{ 1.0 / 6.0, 2.0 / 3.0, 1.0 / 6.0, 0.0}, 00642 {-0.5, 0.0, 0.5, 0.0}, 00643 { 0.5, -1.0, 0.5, 0.0}, 00644 {-1.0 / 6.0, 0.5, -0.5, 1.0 / 6.0}}; 00645 return b; 00646 } 00647 00648 protected: 00649 result_type exec(first_argument_type x, second_argument_type derivative_order) const; 00650 00651 unsigned int derivativeOrder_; 00652 }; 00653 00654 template <class T> 00655 typename BSpline<3, T>::result_type 00656 BSpline<3, T>::exec(first_argument_type x, second_argument_type derivative_order) const 00657 { 00658 switch(derivative_order) 00659 { 00660 case 0: 00661 { 00662 x = VIGRA_CSTD::fabs(x); 00663 if(x < 1.0) 00664 { 00665 return 2.0/3.0 + x*x*(-1.0 + 0.5*x); 00666 } 00667 else if(x < 2.0) 00668 { 00669 x = 2.0 - x; 00670 return x*x*x/6.0; 00671 } 00672 else 00673 return 0.0; 00674 } 00675 case 1: 00676 { 00677 double s = x < 0.0 ? 00678 -1.0 00679 : 1.0; 00680 x = VIGRA_CSTD::fabs(x); 00681 return x < 1.0 ? 00682 s*x*(-2.0 + 1.5*x) 00683 : x < 2.0 ? 00684 -0.5*s*sq(2.0 - x) 00685 : 0.0; 00686 } 00687 case 2: 00688 { 00689 x = VIGRA_CSTD::fabs(x); 00690 return x < 1.0 ? 00691 3.0*x - 2.0 00692 : x < 2.0 ? 00693 2.0 - x 00694 : 0.0; 00695 } 00696 case 3: 00697 { 00698 return x < 0.0 ? 00699 x < -1.0 ? 00700 x < -2.0 ? 00701 0.0 00702 : 1.0 00703 : -3.0 00704 : x < 1.0 ? 00705 3.0 00706 : x < 2.0 ? 00707 -1.0 00708 : 0.0; 00709 } 00710 default: 00711 return 0.0; 00712 } 00713 } 00714 00715 typedef BSpline<3, double> CubicBSplineKernel; 00716 00717 /********************************************************/ 00718 /* */ 00719 /* BSpline<5, T> */ 00720 /* */ 00721 /********************************************************/ 00722 00723 template <class T> 00724 class BSpline<5, T> 00725 { 00726 public: 00727 00728 typedef T value_type; 00729 typedef T argument_type; 00730 typedef T first_argument_type; 00731 typedef unsigned int second_argument_type; 00732 typedef T result_type; 00733 enum StaticOrder { order = 5 }; 00734 00735 explicit BSpline(unsigned int derivativeOrder = 0) 00736 : derivativeOrder_(derivativeOrder) 00737 {} 00738 00739 result_type operator()(argument_type x) const 00740 { 00741 return exec(x, derivativeOrder_); 00742 } 00743 00744 result_type operator()(first_argument_type x, second_argument_type derivative_order) const 00745 { 00746 return exec(x, derivativeOrder_ + derivative_order); 00747 } 00748 00749 result_type dx(argument_type x) const 00750 { return operator()(x, 1); } 00751 00752 result_type dxx(argument_type x) const 00753 { return operator()(x, 2); } 00754 00755 result_type dx3(argument_type x) const 00756 { return operator()(x, 3); } 00757 00758 result_type dx4(argument_type x) const 00759 { return operator()(x, 4); } 00760 00761 value_type operator[](value_type x) const 00762 { return operator()(x); } 00763 00764 double radius() const 00765 { return 3.0; } 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 b[0] = -0.43057534709997114; 00780 b[1] = -0.043096288203264652; 00781 return b; 00782 } 00783 00784 typedef T WeightMatrix[6][6]; 00785 static WeightMatrix & weights() 00786 { 00787 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}, 00788 {-1.0/24.0, -5.0/12.0, 0.0, 5.0/12.0, 1.0/24.0, 0.0}, 00789 { 1.0/12.0, 1.0/6.0, -0.5, 1.0/6.0, 1.0/12.0, 0.0}, 00790 {-1.0/12.0, 1.0/6.0, 0.0, -1.0/6.0, 1.0/12.0, 0.0}, 00791 { 1.0/24.0, -1.0/6.0, 0.25, -1.0/6.0, 1.0/24.0, 0.0}, 00792 {-1.0/120.0, 1.0/24.0, -1.0/12.0, 1.0/12.0, -1.0/24.0, 1.0/120.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<5, T>::result_type 00804 BSpline<5, 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 <= 1.0) 00812 { 00813 return 0.55 + x*x*(-0.5 + x*x*(0.25 - x/12.0)); 00814 } 00815 else if(x < 2.0) 00816 { 00817 return 17.0/40.0 + x*(0.625 + x*(-1.75 + x*(1.25 + x*(-0.375 + x/24.0)))); 00818 } 00819 else if(x < 3.0) 00820 { 00821 x = 3.0 - x; 00822 return x*sq(x*x) / 120.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 <= 1.0) 00834 { 00835 return s*x*(-1.0 + x*x*(1.0 - 5.0/12.0*x)); 00836 } 00837 else if(x < 2.0) 00838 { 00839 return s*(0.625 + x*(-3.5 + x*(3.75 + x*(-1.5 + 5.0/24.0*x)))); 00840 } 00841 else if(x < 3.0) 00842 { 00843 x = 3.0 - x; 00844 return s*sq(x*x) / -24.0; 00845 } 00846 else 00847 return 0.0; 00848 } 00849 case 2: 00850 { 00851 x = VIGRA_CSTD::fabs(x); 00852 if(x <= 1.0) 00853 { 00854 return -1.0 + x*x*(3.0 -5.0/3.0*x); 00855 } 00856 else if(x < 2.0) 00857 { 00858 return -3.5 + x*(7.5 + x*(-4.5 + 5.0/6.0*x)); 00859 } 00860 else if(x < 3.0) 00861 { 00862 x = 3.0 - x; 00863 return x*x*x / 6.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 <= 1.0) 00875 { 00876 return s*x*(6.0 - 5.0*x); 00877 } 00878 else if(x < 2.0) 00879 { 00880 return s*(7.5 + x*(-9.0 + 2.5*x)); 00881 } 00882 else if(x < 3.0) 00883 { 00884 x = 3.0 - x; 00885 return -0.5*s*x*x; 00886 } 00887 else 00888 return 0.0; 00889 } 00890 case 4: 00891 { 00892 x = VIGRA_CSTD::fabs(x); 00893 if(x <= 1.0) 00894 { 00895 return 6.0 - 10.0*x; 00896 } 00897 else if(x < 2.0) 00898 { 00899 return -9.0 + 5.0*x; 00900 } 00901 else if(x < 3.0) 00902 { 00903 return 3.0 - x; 00904 } 00905 else 00906 return 0.0; 00907 } 00908 case 5: 00909 { 00910 return x < 0.0 ? 00911 x < -2.0 ? 00912 x < -3.0 ? 00913 0.0 00914 : 1.0 00915 : x < -1.0 ? 00916 -5.0 00917 : 10.0 00918 : x < 2.0 ? 00919 x < 1.0 ? 00920 -10.0 00921 : 5.0 00922 : x < 3.0 ? 00923 -1.0 00924 : 0.0; 00925 } 00926 default: 00927 return 0.0; 00928 } 00929 } 00930 00931 typedef BSpline<5, double> QuinticBSplineKernel; 00932 00933 #endif // NO_PARTIAL_TEMPLATE_SPECIALIZATION 00934 00935 /********************************************************/ 00936 /* */ 00937 /* CatmullRomSpline */ 00938 /* */ 00939 /********************************************************/ 00940 00941 /** Interpolating 3-rd order splines. 00942 00943 Implements the Catmull/Rom cardinal function 00944 00945 \f[ f(x) = \left\{ \begin{array}{ll} 00946 \frac{3}{2}x^3 - \frac{5}{2}x^2 + 1 & |x| \leq 1 \\ 00947 -\frac{1}{2}x^3 + \frac{5}{2}x^2 -4x + 2 & |x| \leq 2 \\ 00948 0 & \mbox{otherwise} 00949 \end{array}\right. 00950 \f] 00951 00952 It can be used as a functor, and as a kernel for 00953 \ref resamplingConvolveImage() to create a differentiable interpolant 00954 of an image. However, it should be noted that a twice differentiable 00955 interpolant can be created with only slightly more effort by recursive 00956 prefiltering followed by convolution with a 3rd order B-spline. 00957 00958 <b>\#include</b> "<a href="splines_8hxx-source.html">vigra/splines.hxx</a>"<br> 00959 Namespace: vigra 00960 */ 00961 template <class T = double> 00962 class CatmullRomSpline 00963 { 00964 public: 00965 /** the kernel's value type 00966 */ 00967 typedef T value_type; 00968 /** the unary functor's argument type 00969 */ 00970 typedef T argument_type; 00971 /** the unary functor's result type 00972 */ 00973 typedef T result_type; 00974 /** the splines polynomial order 00975 */ 00976 enum StaticOrder { order = 3 }; 00977 00978 /** function (functor) call 00979 */ 00980 result_type operator()(argument_type x) const; 00981 00982 /** index operator -- same as operator() 00983 */ 00984 T operator[] (T x) const 00985 { return operator()(x); } 00986 00987 /** Radius of the function's support. 00988 Needed for \ref resamplingConvolveImage(), always 2. 00989 */ 00990 int radius() const 00991 {return 2;} 00992 00993 /** Derivative order of the function: always 0. 00994 */ 00995 unsigned int derivativeOrder() const 00996 { return 0; } 00997 00998 /** Prefilter coefficients for compatibility with \ref vigra::BSpline. 00999 (array has zero length, since prefiltering is not necessary). 01000 */ 01001 ArrayVector<double> const & prefilterCoefficients() const 01002 { 01003 static ArrayVector<double> b; 01004 return b; 01005 } 01006 }; 01007 01008 01009 template <class T> 01010 typename CatmullRomSpline<T>::result_type 01011 CatmullRomSpline<T>::operator()(argument_type x) const 01012 { 01013 x = VIGRA_CSTD::fabs(x); 01014 if (x <= 1.0) 01015 { 01016 return 1.0 + x * x * (-2.5 + 1.5 * x); 01017 } 01018 else if (x >= 2.0) 01019 { 01020 return 0.0; 01021 } 01022 else 01023 { 01024 return 2.0 + x * (-4.0 + x * (2.5 - 0.5 * x)); 01025 } 01026 } 01027 01028 01029 //@} 01030 01031 } // namespace vigra 01032 01033 01034 #endif /* VIGRA_SPLINES_HXX */
© Ullrich Köthe (koethe@informatik.uni-hamburg.de) |
html generated using doxygen and Python
|