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