[ 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> <<a href="splines_8hxx-source.html">vigra/splines.hxx</a>><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> <<a href="splines_8hxx-source.html">vigra/splines.hxx</a>><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 gevine 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             continous, 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             continous, 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> <<a href="splines_8hxx-source.html">vigra/splines.hxx</a>><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.7.0 (Thu Aug 25 2011)