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

vigra/splines.hxx |

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