[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | 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.6.0, Aug 13 2008 )                                    */
    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 /*        ullrich.koethe@iwr.uni-heidelberg.de    or                    */
    00012 /*        vigra@informatik.uni-hamburg.de                               */
    00013 /*                                                                      */
    00014 /*    Permission is hereby granted, free of charge, to any person       */
    00015 /*    obtaining a copy of this software and associated documentation    */
    00016 /*    files (the "Software"), to deal in the Software without           */
    00017 /*    restriction, including without limitation the rights to use,      */
    00018 /*    copy, modify, merge, publish, distribute, sublicense, and/or      */
    00019 /*    sell copies of the Software, and to permit persons to whom the    */
    00020 /*    Software is furnished to do so, subject to the following          */
    00021 /*    conditions:                                                       */
    00022 /*                                                                      */
    00023 /*    The above copyright notice and this permission notice shall be    */
    00024 /*    included in all copies or substantial portions of the             */
    00025 /*    Software.                                                         */
    00026 /*                                                                      */
    00027 /*    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND    */
    00028 /*    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   */
    00029 /*    OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          */
    00030 /*    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT       */
    00031 /*    HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,      */
    00032 /*    WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING      */
    00033 /*    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR     */
    00034 /*    OTHER DEALINGS IN THE SOFTWARE.                                   */
    00035 /*                                                                      */
    00036 /************************************************************************/
    00037
    00038 #ifndef VIGRA_SPLINES_HXX
    00039 #define VIGRA_SPLINES_HXX
    00040
    00041 #include <cmath>
    00042 #include "config.hxx"
    00043 #include "mathutil.hxx"
    00044 #include "polynomial.hxx"
    00045 #include "array_vector.hxx"
    00046 #include "fixedpoint.hxx"
    00047
    00048 namespace vigra {
    00049
    00050 /** \addtogroup MathFunctions Mathematical Functions
    00051 */
    00052 //@{
    00053 /* B-Splines of arbitrary order and interpolating Catmull/Rom splines.
    00054
    00055     <b>\#include</b> <<a href="splines_8hxx-source.html">vigra/splines.hxx</a>><br>
    00056     Namespace: vigra
    00057 */
    00058 #ifndef NO_PARTIAL_TEMPLATE_SPECIALIZATION
    00059
    00060 /** Basic interface of the spline functors.
    00061
    00062     Implements the spline functions defined by the recursion
    00063
    00064     \f[ B_0(x) = \left\{ \begin{array}{ll}
    00065                                   1 & -\frac{1}{2} \leq x < \frac{1}{2} \\
    00066                                   0 & \mbox{otherwise}
    00067                         \end{array}\right.
    00068     \f]
    00069
    00070     and
    00071
    00072     \f[ B_n(x) = B_0(x) * B_{n-1}(x)
    00073     \f]
    00074
    00075     where * denotes convolution, and <i>n</i> is the spline order given by the
    00076     template parameter <tt>ORDER</tt>. These spline classes can be used as
    00077     unary and binary functors, as kernels for \ref resamplingConvolveImage(),
    00078     and as arguments for \ref vigra::SplineImageView. Note that the spline order
    00079     is given as a template argument.
    00080
    00081     <b>\#include</b> <<a href="splines_8hxx-source.html">vigra/splines.hxx</a>><br>
    00082     Namespace: vigra
    00083 */
    00084 template <int ORDER, class T = double>
    00085 class BSplineBase
    00086 {
    00087   public:
    00088
    00089         /** the value type if used as a kernel in \ref resamplingConvolveImage().
    00090         */
    00091     typedef T            value_type;
    00092         /** the functor's unary argument type
    00093         */
    00094     typedef T            argument_type;
    00095         /** the functor's first binary argument type
    00096         */
    00097     typedef T            first_argument_type;
    00098         /** the functor's second binary argument type
    00099         */
    00100     typedef unsigned int second_argument_type;
    00101         /** the functor's result type (unary and binary)
    00102         */
    00103     typedef T            result_type;
    00104         /** the spline order
    00105         */
    00106     enum StaticOrder { order = ORDER };
    00107
    00108         /** Create functor for gevine derivative of the spline. The spline's order
    00109             is specified spline by the template argument <TT>ORDER</tt>.
    00110         */
    00111     explicit BSplineBase(unsigned int derivativeOrder = 0)
    00112     : s1_(derivativeOrder)
    00113     {}
    00114
    00115         /** Unary function call.
    00116             Returns the value of the spline with the derivative order given in the
    00117             constructor. Note that only derivatives up to <tt>ORDER-1</tt> are
    00118             continous, and derivatives above <tt>ORDER+1</tt> are zero.
    00119         */
    00120     result_type operator()(argument_type x) const
    00121     {
    00122         return exec(x, derivativeOrder());
    00123     }
    00124
    00125         /** Binary function call.
    00126             The given derivative order is added to the derivative order
    00127             specified in the constructor. Note that only derivatives up to <tt>ORDER-1</tt> are
    00128             continous, and derivatives above <tt>ORDER+1</tt> are zero.
    00129         */
    00130     result_type operator()(first_argument_type x, second_argument_type derivative_order) const
    00131     {
    00132          return exec(x, derivativeOrder() + derivative_order);
    00133     }
    00134
    00135         /** Index operator. Same as unary function call.
    00136         */
    00137     value_type operator[](value_type x) const
    00138         { return operator()(x); }
    00139
    00140         /** Get the required filter radius for a discrete approximation of the
    00141             spline. Always equal to <tt>(ORDER + 1) / 2.0</tt>.
    00142         */
    00143     double radius() const
    00144         { return (ORDER + 1) * 0.5; }
    00145
    00146         /** Get the derivative order of the Gaussian.
    00147         */
    00148     unsigned int derivativeOrder() const
    00149         { return s1_.derivativeOrder(); }
    00150
    00151         /** Get the prefilter coefficients required for interpolation.
    00152             To interpolate with a B-spline, \ref resamplingConvolveImage()
    00153             can be used. However, the image to be interpolated must be
    00154             pre-filtered using \ref recursiveFilterX() and \ref recursiveFilterY()
    00155             with the filter coefficients given by this function. The length of the array
    00156             corresponds to how many times the above recursive filtering
    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; // (note: integer division)
    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     result_type dx(argument_type x) const
    00492         { return operator()(x, 1); }
    00493
    00494     value_type operator[](value_type x) const
    00495         { return operator()(x); }
    00496
    00497     double radius() const
    00498         { return 1.5; }
    00499
    00500     unsigned int derivativeOrder() const
    00501         { return derivativeOrder_; }
    00502
    00503     ArrayVector<double> const & prefilterCoefficients() const
    00504     {
    00505         static ArrayVector<double> b(1, 2.0*M_SQRT2 - 3.0);
    00506         return b;
    00507     }
    00508
    00509     typedef T WeightMatrix[3][3];
    00510     static WeightMatrix & weights()
    00511     {
    00512         static T b[3][3] = {{ 0.125, 0.75, 0.125},
    00513                             {-0.5, 0.0, 0.5},
    00514                             { 0.5, -1.0, 0.5}};
    00515         return b;
    00516     }
    00517
    00518   protected:
    00519     result_type exec(first_argument_type x, second_argument_type derivative_order) const;
    00520
    00521     unsigned int derivativeOrder_;
    00522 };
    00523
    00524 template <class T>
    00525 typename BSpline<2, T>::result_type
    00526 BSpline<2, T>::exec(first_argument_type x, second_argument_type derivative_order) const
    00527 {
    00528     switch(derivative_order)
    00529     {
    00530         case 0:
    00531         {
    00532             x = VIGRA_CSTD::fabs(x);
    00533             return x < 0.5 ?
    00534                     0.75 - x*x
    00535                     : x < 1.5 ?
    00536                         0.5 * sq(1.5 - x)
    00537                     : 0.0;
    00538         }
    00539         case 1:
    00540         {
    00541             return x >= -0.5 ?
    00542                      x <= 0.5 ?
    00543                        -2.0 * x
    00544                      : x < 1.5 ?
    00545                          x - 1.5
    00546                        : 0.0
    00547                    : x > -1.5 ?
    00548                        x + 1.5
    00549                      : 0.0;
    00550         }
    00551         case 2:
    00552         {
    00553             return x >= -0.5 ?
    00554                      x < 0.5 ?
    00555                          -2.0
    00556                      : x < 1.5 ?
    00557                          1.0
    00558                        : 0.0
    00559                    : x >= -1.5 ?
    00560                        1.0
    00561                      : 0.0;
    00562         }
    00563         default:
    00564             return 0.0;
    00565     }
    00566 }
    00567
    00568 /********************************************************/
    00569 /*                                                      */
    00570 /*                     BSpline<3, T>                    */
    00571 /*                                                      */
    00572 /********************************************************/
    00573
    00574 template <class T>
    00575 class BSpline<3, T>
    00576 {
    00577   public:
    00578
    00579     typedef T            value_type;
    00580     typedef T            argument_type;
    00581     typedef T            first_argument_type;
    00582     typedef unsigned int second_argument_type;
    00583     typedef T            result_type;
    00584     enum StaticOrder { order = 3 };
    00585
    00586     explicit BSpline(unsigned int derivativeOrder = 0)
    00587     : derivativeOrder_(derivativeOrder)
    00588     {}
    00589
    00590     result_type operator()(argument_type x) const
    00591     {
    00592         return exec(x, derivativeOrder_);
    00593     }
    00594
    00595     template <unsigned int IntBits, unsigned int FracBits>
    00596     FixedPoint<IntBits, FracBits> operator()(FixedPoint<IntBits, FracBits> x) const
    00597     {
    00598         typedef FixedPoint<IntBits, FracBits> Value;
    00599         enum { ONE = Value::ONE, TWO = 2 * ONE, TWO_THIRDS = TWO / 3, ONE_SIXTH = ONE / 6,
    00600                PREMULTIPLY_SHIFT = FracBits <= 16 ? 0 : FracBits - 16,
    00601                POSTMULTIPLY_SHIFT = FracBits - 2*PREMULTIPLY_SHIFT };
    00602         int v = abs(x.value);
    00603         return v == ONE
    00604                    ? Value(ONE_SIXTH, FPNoShift)
    00605                    : v < ONE
    00606                        ? Value(TWO_THIRDS +
    00607                                (((int)(sq((unsigned)v >> PREMULTIPLY_SHIFT) >> (POSTMULTIPLY_SHIFT + PREMULTIPLY_SHIFT))
    00608                                        * (((v >> 1) - ONE) >> PREMULTIPLY_SHIFT)) >> POSTMULTIPLY_SHIFT), FPNoShift)
    00609                        : v < TWO
    00610                             ? Value((int)((sq((unsigned)(TWO-v) >> PREMULTIPLY_SHIFT) >> (POSTMULTIPLY_SHIFT + PREMULTIPLY_SHIFT))
    00611                                       * ((unsigned)(TWO-v) >> PREMULTIPLY_SHIFT) / 6) >> POSTMULTIPLY_SHIFT, FPNoShift)
    00612                             : Value(0, FPNoShift);
    00613     }
    00614
    00615     result_type operator()(first_argument_type x, second_argument_type derivative_order) const
    00616     {
    00617          return exec(x, derivativeOrder_ + derivative_order);
    00618     }
    00619
    00620     result_type dx(argument_type x) const
    00621         { return operator()(x, 1); }
    00622
    00623     result_type dxx(argument_type x) const
    00624         { return operator()(x, 2); }
    00625
    00626     value_type operator[](value_type x) const
    00627         { return operator()(x); }
    00628
    00629     double radius() const
    00630         { return 2.0; }
    00631
    00632     unsigned int derivativeOrder() const
    00633         { return derivativeOrder_; }
    00634
    00635     ArrayVector<double> const & prefilterCoefficients() const
    00636     {
    00637         static ArrayVector<double> b(1, VIGRA_CSTD::sqrt(3.0) - 2.0);
    00638         return b;
    00639     }
    00640
    00641     typedef T WeightMatrix[4][4];
    00642     static WeightMatrix & weights()
    00643     {
    00644         static T b[4][4] = {{ 1.0 / 6.0, 2.0 / 3.0, 1.0 / 6.0, 0.0},
    00645                             {-0.5, 0.0, 0.5, 0.0},
    00646                             { 0.5, -1.0, 0.5, 0.0},
    00647                             {-1.0 / 6.0, 0.5, -0.5, 1.0 / 6.0}};
    00648         return b;
    00649     }
    00650
    00651   protected:
    00652     result_type exec(first_argument_type x, second_argument_type derivative_order) const;
    00653
    00654     unsigned int derivativeOrder_;
    00655 };
    00656
    00657 template <class T>
    00658 typename BSpline<3, T>::result_type
    00659 BSpline<3, T>::exec(first_argument_type x, second_argument_type derivative_order) const
    00660 {
    00661     switch(derivative_order)
    00662     {
    00663         case 0:
    00664         {
    00665             x = VIGRA_CSTD::fabs(x);
    00666             if(x < 1.0)
    00667             {
    00668                 return 2.0/3.0 + x*x*(-1.0 + 0.5*x);
    00669             }
    00670             else if(x < 2.0)
    00671             {
    00672                 x = 2.0 - x;
    00673                 return x*x*x/6.0;
    00674             }
    00675             else
    00676                 return 0.0;
    00677         }
    00678         case 1:
    00679         {
    00680             double s = x < 0.0 ?
    00681                          -1.0
    00682                        :  1.0;
    00683             x = VIGRA_CSTD::fabs(x);
    00684             return x < 1.0 ?
    00685                      s*x*(-2.0 + 1.5*x)
    00686                    : x < 2.0 ?
    00687                        -0.5*s*sq(2.0 - x)
    00688                      : 0.0;
    00689         }
    00690         case 2:
    00691         {
    00692             x = VIGRA_CSTD::fabs(x);
    00693             return x < 1.0 ?
    00694                      3.0*x - 2.0
    00695                    : x < 2.0 ?
    00696                        2.0 - x
    00697                      : 0.0;
    00698         }
    00699         case 3:
    00700         {
    00701             return x < 0.0 ?
    00702                      x < -1.0 ?
    00703                        x < -2.0 ?
    00704                          0.0
    00705                        : 1.0
    00706                      : -3.0
    00707                    : x < 1.0 ?
    00708                        3.0
    00709                      : x < 2.0 ?
    00710                          -1.0
    00711                        : 0.0;
    00712         }
    00713         default:
    00714             return 0.0;
    00715     }
    00716 }
    00717
    00718 typedef BSpline<3, double> CubicBSplineKernel;
    00719
    00720 /********************************************************/
    00721 /*                                                      */
    00722 /*                     BSpline<4, T>                    */
    00723 /*                                                      */
    00724 /********************************************************/
    00725
    00726 template <class T>
    00727 class BSpline<4, T>
    00728 {
    00729   public:
    00730
    00731     typedef T            value_type;
    00732     typedef T            argument_type;
    00733     typedef T            first_argument_type;
    00734     typedef unsigned int second_argument_type;
    00735     typedef T            result_type;
    00736     enum StaticOrder { order = 4 };
    00737
    00738     explicit BSpline(unsigned int derivativeOrder = 0)
    00739     : derivativeOrder_(derivativeOrder)
    00740     {}
    00741
    00742     result_type operator()(argument_type x) const
    00743     {
    00744         return exec(x, derivativeOrder_);
    00745     }
    00746
    00747     result_type operator()(first_argument_type x, second_argument_type derivative_order) const
    00748     {
    00749          return exec(x, derivativeOrder_ + derivative_order);
    00750     }
    00751
    00752     result_type dx(argument_type x) const
    00753         { return operator()(x, 1); }
    00754
    00755     result_type dxx(argument_type x) const
    00756         { return operator()(x, 2); }
    00757
    00758     result_type dx3(argument_type x) const
    00759         { return operator()(x, 3); }
    00760
    00761     value_type operator[](value_type x) const
    00762         { return operator()(x); }
    00763
    00764     double radius() const
    00765         { return 2.5; }
    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         // -19 + 4*sqrt(19) + 2*sqrt(2*(83 - 19*sqrt(19)))
    00780         b[0] = -0.361341225900220177092212841325;
    00781         // -19 - 4*sqrt(19) + 2*sqrt(2*(83 + 19*sqrt(19)))
    00782         b[1] = -0.013725429297339121360331226939;
    00783         return b;
    00784     }
    00785
    00786     typedef T WeightMatrix[5][5];
    00787     static WeightMatrix & weights()
    00788     {
    00789         static T b[5][5] = {{ 1.0/384.0, 19.0/96.0, 115.0/192.0, 19.0/96.0, 1.0/384.0},
    00790                             {-1.0/48.0, -11.0/24.0, 0.0, 11.0/24.0, 1.0/48.0},
    00791                             { 1.0/16.0, 1.0/4.0, -5.0/8.0, 1.0/4.0, 1.0/16.0},
    00792                             {-1.0/12.0, 1.0/6.0, 0.0, -1.0/6.0, 1.0/12.0},
    00793                             { 1.0/24.0, -1.0/6.0, 0.25, -1.0/6.0, 1.0/24.0}};
    00794         return b;
    00795     }
    00796
    00797   protected:
    00798     result_type exec(first_argument_type x, second_argument_type derivative_order) const;
    00799
    00800     unsigned int derivativeOrder_;
    00801 };
    00802
    00803 template <class T>
    00804 typename BSpline<4, T>::result_type
    00805 BSpline<4, T>::exec(first_argument_type x, second_argument_type derivative_order) const
    00806 {
    00807     switch(derivative_order)
    00808     {
    00809         case 0:
    00810         {
    00811             x = VIGRA_CSTD::fabs(x);
    00812             if(x <= 0.5)
    00813             {
    00814                 return 115.0/192.0 + x*x*(-0.625 + x*x*0.25);
    00815             }
    00816             else if(x < 1.5)
    00817             {
    00818                 return (55.0/16.0 + x*(1.25 + x*(-7.5 + x*(5.0 - x)))) / 6.0;
    00819             }
    00820             else if(x < 2.5)
    00821             {
    00822                 x = 2.5 - x;
    00823                 return sq(x*x) / 24.0;
    00824             }
    00825             else
    00826                 return 0.0;
    00827         }
    00828         case 1:
    00829         {
    00830             double s = x < 0.0 ?
    00831                           -1.0 :
    00832                            1.0;
    00833             x = VIGRA_CSTD::fabs(x);
    00834             if(x <= 0.5)
    00835             {
    00836                 return s*x*(-1.25 + x*x);
    00837             }
    00838             else if(x < 1.5)
    00839             {
    00840                 return s*(5.0 + x*(-60.0 + x*(60.0 - 16.0*x))) / 24.0;
    00841             }
    00842             else if(x < 2.5)
    00843             {
    00844                 x = 2.5 - x;
    00845                 return s*x*x*x / -6.0;
    00846             }
    00847             else
    00848                 return 0.0;
    00849         }
    00850         case 2:
    00851         {
    00852             x = VIGRA_CSTD::fabs(x);
    00853             if(x <= 0.5)
    00854             {
    00855                 return -1.25 + 3.0*x*x;
    00856             }
    00857             else if(x < 1.5)
    00858             {
    00859                 return -2.5 + x*(5.0 - 2.0*x);
    00860             }
    00861             else if(x < 2.5)
    00862             {
    00863                 x = 2.5 - x;
    00864                 return x*x / 2.0;
    00865             }
    00866             else
    00867                 return 0.0;
    00868         }
    00869         case 3:
    00870         {
    00871             double s = x < 0.0 ?
    00872                           -1.0 :
    00873                            1.0;
    00874             x = VIGRA_CSTD::fabs(x);
    00875             if(x <= 0.5)
    00876             {
    00877                 return s*x*6.0;
    00878             }
    00879             else if(x < 1.5)
    00880             {
    00881                 return s*(5.0 - 4.0*x);
    00882             }
    00883             else if(x < 2.5)
    00884             {
    00885                 return s*(x - 2.5);
    00886             }
    00887             else
    00888                 return 0.0;
    00889         }
    00890         case 4:
    00891         {
    00892             return x < 0.0
    00893                      ? x < -2.5
    00894                          ? 0.0
    00895                          : x < -1.5
    00896                              ? 1.0
    00897                              : x < -0.5
    00898                                  ? -4.0
    00899                                  : 6.0
    00900                      : x < 0.5
    00901                          ? 6.0
    00902                          : x < 1.5
    00903                              ? -4.0
    00904                              : x < 2.5
    00905                                  ? 1.0
    00906                                  : 0.0;
    00907         }
    00908         default:
    00909             return 0.0;
    00910     }
    00911 }
    00912
    00913 typedef BSpline<4, double> QuarticBSplineKernel;
    00914
    00915 /********************************************************/
    00916 /*                                                      */
    00917 /*                     BSpline<5, T>                    */
    00918 /*                                                      */
    00919 /********************************************************/
    00920
    00921 template <class T>
    00922 class BSpline<5, T>
    00923 {
    00924   public:
    00925
    00926     typedef T            value_type;
    00927     typedef T            argument_type;
    00928     typedef T            first_argument_type;
    00929     typedef unsigned int second_argument_type;
    00930     typedef T            result_type;
    00931     enum StaticOrder { order = 5 };
    00932
    00933     explicit BSpline(unsigned int derivativeOrder = 0)
    00934     : derivativeOrder_(derivativeOrder)
    00935     {}
    00936
    00937     result_type operator()(argument_type x) const
    00938     {
    00939         return exec(x, derivativeOrder_);
    00940     }
    00941
    00942     result_type operator()(first_argument_type x, second_argument_type derivative_order) const
    00943     {
    00944          return exec(x, derivativeOrder_ + derivative_order);
    00945     }
    00946
    00947     result_type dx(argument_type x) const
    00948         { return operator()(x, 1); }
    00949
    00950     result_type dxx(argument_type x) const
    00951         { return operator()(x, 2); }
    00952
    00953     result_type dx3(argument_type x) const
    00954         { return operator()(x, 3); }
    00955
    00956     result_type dx4(argument_type x) const
    00957         { return operator()(x, 4); }
    00958
    00959     value_type operator[](value_type x) const
    00960         { return operator()(x); }
    00961
    00962     double radius() const
    00963         { return 3.0; }
    00964
    00965     unsigned int derivativeOrder() const
    00966         { return derivativeOrder_; }
    00967
    00968     ArrayVector<double> const & prefilterCoefficients() const
    00969     {
    00970         static ArrayVector<double> const & b = initPrefilterCoefficients();
    00971         return b;
    00972     }
    00973
    00974     static ArrayVector<double> const & initPrefilterCoefficients()
    00975     {
    00976         static ArrayVector<double> b(2);
    00977         // -(13/2) + sqrt(105)/2 + sqrt(1/2*((135 - 13*sqrt(105))))
    00978         b[0] = -0.430575347099973791851434783493;
    00979         // (1/2)*((-13) - sqrt(105) + sqrt(2*((135 + 13*sqrt(105)))))
    00980         b[1] = -0.043096288203264653822712376822;
    00981         return b;
    00982     }
    00983
    00984     typedef T WeightMatrix[6][6];
    00985     static WeightMatrix & weights()
    00986     {
    00987         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},
    00988                             {-1.0/24.0, -5.0/12.0, 0.0, 5.0/12.0, 1.0/24.0, 0.0},
    00989                             { 1.0/12.0, 1.0/6.0, -0.5, 1.0/6.0, 1.0/12.0, 0.0},
    00990                             {-1.0/12.0, 1.0/6.0, 0.0, -1.0/6.0, 1.0/12.0, 0.0},
    00991                             { 1.0/24.0, -1.0/6.0, 0.25, -1.0/6.0, 1.0/24.0, 0.0},
    00992                             {-1.0/120.0, 1.0/24.0, -1.0/12.0, 1.0/12.0, -1.0/24.0, 1.0/120.0}};
    00993         return b;
    00994     }
    00995
    00996   protected:
    00997     result_type exec(first_argument_type x, second_argument_type derivative_order) const;
    00998
    00999     unsigned int derivativeOrder_;
    01000 };
    01001
    01002 template <class T>
    01003 typename BSpline<5, T>::result_type
    01004 BSpline<5, T>::exec(first_argument_type x, second_argument_type derivative_order) const
    01005 {
    01006     switch(derivative_order)
    01007     {
    01008         case 0:
    01009         {
    01010             x = VIGRA_CSTD::fabs(x);
    01011             if(x <= 1.0)
    01012             {
    01013                 return 0.55 + x*x*(-0.5 + x*x*(0.25 - x/12.0));
    01014             }
    01015             else if(x < 2.0)
    01016             {
    01017                 return 17.0/40.0 + x*(0.625 + x*(-1.75 + x*(1.25 + x*(-0.375 + x/24.0))));
    01018             }
    01019             else if(x < 3.0)
    01020             {
    01021                 x = 3.0 - x;
    01022                 return x*sq(x*x) / 120.0;
    01023             }
    01024             else
    01025                 return 0.0;
    01026         }
    01027         case 1:
    01028         {
    01029             double s = x < 0.0 ?
    01030                           -1.0 :
    01031                            1.0;
    01032             x = VIGRA_CSTD::fabs(x);
    01033             if(x <= 1.0)
    01034             {
    01035                 return s*x*(-1.0 + x*x*(1.0 - 5.0/12.0*x));
    01036             }
    01037             else if(x < 2.0)
    01038             {
    01039                 return s*(0.625 + x*(-3.5 + x*(3.75 + x*(-1.5 + 5.0/24.0*x))));
    01040             }
    01041             else if(x < 3.0)
    01042             {
    01043                 x = 3.0 - x;
    01044                 return s*sq(x*x) / -24.0;
    01045             }
    01046             else
    01047                 return 0.0;
    01048         }
    01049         case 2:
    01050         {
    01051             x = VIGRA_CSTD::fabs(x);
    01052             if(x <= 1.0)
    01053             {
    01054                 return -1.0 + x*x*(3.0 -5.0/3.0*x);
    01055             }
    01056             else if(x < 2.0)
    01057             {
    01058                 return -3.5 + x*(7.5 + x*(-4.5 + 5.0/6.0*x));
    01059             }
    01060             else if(x < 3.0)
    01061             {
    01062                 x = 3.0 - x;
    01063                 return x*x*x / 6.0;
    01064             }
    01065             else
    01066                 return 0.0;
    01067         }
    01068         case 3:
    01069         {
    01070             double s = x < 0.0 ?
    01071                           -1.0 :
    01072                            1.0;
    01073             x = VIGRA_CSTD::fabs(x);
    01074             if(x <= 1.0)
    01075             {
    01076                 return s*x*(6.0 - 5.0*x);
    01077             }
    01078             else if(x < 2.0)
    01079             {
    01080                 return s*(7.5 + x*(-9.0 + 2.5*x));
    01081             }
    01082             else if(x < 3.0)
    01083             {
    01084                 x = 3.0 - x;
    01085                 return -0.5*s*x*x;
    01086             }
    01087             else
    01088                 return 0.0;
    01089         }
    01090         case 4:
    01091         {
    01092             x = VIGRA_CSTD::fabs(x);
    01093             if(x <= 1.0)
    01094             {
    01095                 return 6.0 - 10.0*x;
    01096             }
    01097             else if(x < 2.0)
    01098             {
    01099                 return -9.0 + 5.0*x;
    01100             }
    01101             else if(x < 3.0)
    01102             {
    01103                 return 3.0 - x;
    01104             }
    01105             else
    01106                 return 0.0;
    01107         }
    01108         case 5:
    01109         {
    01110             return x < 0.0 ?
    01111                      x < -2.0 ?
    01112                        x < -3.0 ?
    01113                          0.0
    01114                        : 1.0
    01115                      : x < -1.0 ?
    01116                          -5.0
    01117                        : 10.0
    01118                    : x < 2.0 ?
    01119                        x < 1.0 ?
    01120                          -10.0
    01121                        : 5.0
    01122                      : x < 3.0 ?
    01123                          -1.0
    01124                        : 0.0;
    01125         }
    01126         default:
    01127             return 0.0;
    01128     }
    01129 }
    01130
    01131 typedef BSpline<5, double> QuinticBSplineKernel;
    01132
    01133 #endif // NO_PARTIAL_TEMPLATE_SPECIALIZATION
    01134
    01135 /********************************************************/
    01136 /*                                                      */
    01137 /*                      CatmullRomSpline                */
    01138 /*                                                      */
    01139 /********************************************************/
    01140
    01141 /** Interpolating 3-rd order splines.
    01142
    01143     Implements the Catmull/Rom cardinal function
    01144
    01145     \f[ f(x) = \left\{ \begin{array}{ll}
    01146                                   \frac{3}{2}x^3 - \frac{5}{2}x^2 + 1 & |x| \leq 1 \\
    01147                                   -\frac{1}{2}x^3 + \frac{5}{2}x^2 -4x + 2 & |x| \leq 2 \\
    01148                                   0 & \mbox{otherwise}
    01149                         \end{array}\right.
    01150     \f]
    01151
    01152     It can be used as a functor, and as a kernel for
    01153     \ref resamplingConvolveImage() to create a differentiable interpolant
    01154     of an image. However, it should be noted that a twice differentiable
    01155     interpolant can be created with only slightly more effort by recursive
    01156     prefiltering followed by convolution with a 3rd order B-spline.
    01157
    01158     <b>\#include</b> <<a href="splines_8hxx-source.html">vigra/splines.hxx</a>><br>
    01159     Namespace: vigra
    01160 */
    01161 template <class T = double>
    01162 class CatmullRomSpline
    01163 {
    01164 public:
    01165         /** the kernel's value type
    01166         */
    01167     typedef T value_type;
    01168         /** the unary functor's argument type
    01169         */
    01170     typedef T argument_type;
    01171         /** the unary functor's result type
    01172         */
    01173     typedef T result_type;
    01174         /** the splines polynomial order
    01175         */
    01176     enum StaticOrder { order = 3 };
    01177
    01178         /** function (functor) call
    01179         */
    01180     result_type operator()(argument_type x) const;
    01181
    01182         /** index operator -- same as operator()
    01183         */
    01184     T operator[] (T x) const
    01185         { return operator()(x); }
    01186
    01187         /** Radius of the function's support.
    01188             Needed for  \ref resamplingConvolveImage(), always 2.
    01189         */
    01190     int radius() const
    01191         {return 2;}
    01192
    01193         /** Derivative order of the function: always 0.
    01194         */
    01195     unsigned int derivativeOrder() const
    01196         { return 0; }
    01197
    01198         /** Prefilter coefficients for compatibility with \ref vigra::BSpline.
    01199             (array has zero length, since prefiltering is not necessary).
    01200         */
    01201     ArrayVector<double> const & prefilterCoefficients() const
    01202     {
    01203         static ArrayVector<double> b;
    01204         return b;
    01205     }
    01206 };
    01207
    01208
    01209 template <class T>
    01210 typename CatmullRomSpline<T>::result_type
    01211 CatmullRomSpline<T>::operator()(argument_type x) const
    01212 {
    01213     x = VIGRA_CSTD::fabs(x);
    01214     if (x <= 1.0)
    01215     {
    01216         return 1.0 + x * x * (-2.5 + 1.5 * x);
    01217     }
    01218     else if (x >= 2.0)
    01219     {
    01220         return 0.0;
    01221     }
    01222     else
    01223     {
    01224         return 2.0 + x * (-4.0 + x * (2.5 - 0.5 * x));
    01225     }
    01226 }
    01227
    01228
    01229 //@}
    01230
    01231 } // namespace vigra
    01232
    01233
    01234 #endif /* VIGRA_SPLINES_HXX */

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de) Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany | html generated using doxygen and Python VIGRA 1.6.0 (13 Aug 2008)