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