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

vigra/splines.hxx

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

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
VIGRA 1.6.0 (13 Aug 2008)