[ VIGRA Homepage | Class Index | Function Index | File Index | Main Page ]

details vigra/splines.hxx VIGRA

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

© Ullrich Köthe (koethe@informatik.uni-hamburg.de)
Cognitive Systems Group, University of Hamburg, Germany

html generated using doxygen and Python
VIGRA 1.3.3 (18 Aug 2005)