[ 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.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)