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

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

html generated using doxygen and Python
VIGRA 1.4.0 (21 Dec 2005)