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

| vigra/splineimageview.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.3.2, Jan 27 2005 )                                    */
    00008 /*    You may use, modify, and distribute this software according       */
    00009 /*    to the terms stated in the LICENSE file included in               */
    00010 /*    the VIGRA distribution.                                           */
    00011 /*                                                                      */
    00012 /*    The VIGRA Website is                                              */
    00013 /*        http://kogs-www.informatik.uni-hamburg.de/~koethe/vigra/      */
    00014 /*    Please direct questions, bug reports, and contributions to        */
    00015 /*        koethe@informatik.uni-hamburg.de                              */
    00016 /*                                                                      */
    00017 /*  THIS SOFTWARE IS PROVIDED AS IS AND WITHOUT ANY EXPRESS OR          */
    00018 /*  IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED      */
    00019 /*  WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. */
    00020 /*                                                                      */
    00021 /************************************************************************/
    00022
    00023 #ifndef VIGRA_SPLINEIMAGEVIEW_HXX
    00024 #define VIGRA_SPLINEIMAGEVIEW_HXX
    00025
    00026 #include "vigra/mathutil.hxx"
    00027 #include "vigra/recursiveconvolution.hxx"
    00028 #include "vigra/splines.hxx"
    00029 #include "vigra/array_vector.hxx"
    00030
    00031 namespace vigra {
    00032
    00033 /********************************************************/
    00034 /*                                                      */
    00035 /*                    SplineImageView                   */
    00036 /*                                                      */
    00037 /********************************************************/
    00038 /** \brief Create a continuous view onto a discrete image using splines.
    00039
    00040     This class is very useful if image values or derivatives at arbitrary
    00041     real-valued coordinates are needed. Access at such coordinates is implemented by
    00042     interpolating the given discrete image values with a spline of the
    00043     specified <tt>ORDER</TT>. Continuous derivatives are available up to
    00044     degree <tt>ORDER-1</TT>. If the requested coordinates are near the image border,
    00045     reflective boundary conditions are applied. In principle, this class can also be used
    00046     for image resizing, but here the functions from the <tt>resize...</tt> family are
    00047     more efficient.
    00048
    00049 */
    00050 template <int ORDER, class VALUETYPE>
    00051 class SplineImageView
    00052 {
    00053     typedef typename NumericTraits<VALUETYPE>::RealPromote InternalValue;
    00054
    00055   public:
    00056
    00057         /** The view's value type (return type of access and derivative functions).
    00058         */
    00059     typedef VALUETYPE value_type;
    00060
    00061         /** The view's size type.
    00062         */
    00063     typedef Size2D size_type;
    00064
    00065         /** The view's difference type.
    00066         */
    00067     typedef TinyVector<double, 2> difference_type;
    00068
    00069         /** The order of the spline used.
    00070         */
    00071     enum StaticOrder { order = ORDER };
    00072
    00073         /** The type of the internal image holding the spline coefficients.
    00074         */
    00075     typedef BasicImage<InternalValue> InternalImage;
    00076
    00077   private:
    00078     typedef typename InternalImage::traverser InternalTraverser;
    00079     typedef typename InternalTraverser::row_iterator InternalRowIterator;
    00080     typedef typename InternalTraverser::column_iterator InternalColumnIterator;
    00081     typedef BSpline<ORDER, double> Spline;
    00082
    00083     enum { ksize_ = ORDER + 1, kcenter_ = ORDER / 2 };
    00084
    00085   public:
    00086         /** Construct SplineImageView for the given image.
    00087
    00088             If <tt>skipPrefiltering = true</tt> (default: <tt>false</tt>), the recursive
    00089             prefilter of the cardinal spline function is not applied, resulting
    00090             in an approximating (smoothing) rather than interpolating spline. This is
    00091             especially useful if customized prefilters are to be applied.
    00092         */
    00093     template <class SrcIterator, class SrcAccessor>
    00094     SplineImageView(SrcIterator is, SrcIterator iend, SrcAccessor sa, bool skipPrefiltering = false)
    00095     : w_(iend.x - is.x), h_(iend.y - is.y), w1_(w_-1), h1_(h_-1),
    00096       x0_(kcenter_), x1_(w_ - kcenter_ - 2), y0_(kcenter_), y1_(h_ - kcenter_ - 2),
    00097       image_(w_, h_),
    00098       x_(-1.0), y_(-1.0),
    00099       u_(-1.0), v_(-1.0)
    00100     {
    00101         copyImage(srcIterRange(is, iend, sa), destImage(image_));
    00102         if(!skipPrefiltering)
    00103             init();
    00104     }
    00105
    00106         /** Construct SplineImageView for the given image.
    00107
    00108             If <tt>skipPrefiltering = true</tt> (default: <tt>false</tt>), the recursive
    00109             prefilter of the cardinal spline function is not applied, resulting
    00110             in an approximating (smoothing) rather than interpolating spline. This is
    00111             especially useful if customized prefilters are to be applied.
    00112         */
    00113     template <class SrcIterator, class SrcAccessor>
    00114     SplineImageView(triple<SrcIterator, SrcIterator, SrcAccessor> s, bool skipPrefiltering = false)
    00115     : w_(s.second.x - s.first.x), h_(s.second.y - s.first.y), w1_(w_-1), h1_(h_-1),
    00116       x0_(kcenter_), x1_(w_ - kcenter_ - 2), y0_(kcenter_), y1_(h_ - kcenter_ - 2),
    00117       image_(w_, h_),
    00118       x_(-1.0), y_(-1.0),
    00119       u_(-1.0), v_(-1.0)
    00120     {
    00121         copyImage(srcIterRange(s.first, s.second, s.third), destImage(image_));
    00122         if(!skipPrefiltering)
    00123             init();
    00124     }
    00125
    00126         /** Access interpolated function at real-valued coordinate <tt>(x, y)</tt>.
    00127         */
    00128     value_type operator()(double x, double y) const;
    00129
    00130         /** Access derivative of order <tt>(dx, dy)</tt> at real-valued coordinate <tt>(x, y)</tt>.
    00131         */
    00132     value_type operator()(double x, double y, unsigned int dx, unsigned int dy) const;
    00133
    00134         /** Access 1st derivative in x-direction at real-valued coordinate <tt>(x, y)</tt>.
    00135             Equivalent to <tt>splineView(x, y, 1, 0)</tt>.
    00136         */
    00137     value_type dx(double x, double y) const
    00138         { return operator()(x, y, 1, 0); }
    00139
    00140         /** Access 1st derivative in y-direction at real-valued coordinate <tt>(x, y)</tt>.
    00141             Equivalent to <tt>splineView(x, y, 0, 1)</tt>.
    00142         */
    00143     value_type dy(double x, double y) const
    00144         { return operator()(x, y, 0, 1); }
    00145
    00146         /** Access 2nd derivative in x-direction at real-valued coordinate <tt>(x, y)</tt>.
    00147             Equivalent to <tt>splineView(x, y, 2, 0)</tt>.
    00148         */
    00149     value_type dxx(double x, double y) const
    00150         { return operator()(x, y, 2, 0); }
    00151
    00152         /** Access mixed 2nd derivative at real-valued coordinate <tt>(x, y)</tt>.
    00153             Equivalent to <tt>splineView(x, y, 1, 1)</tt>.
    00154         */
    00155     value_type dxy(double x, double y) const
    00156         { return operator()(x, y, 1, 1); }
    00157
    00158         /** Access 2nd derivative in y-direction at real-valued coordinate <tt>(x, y)</tt>.
    00159             Equivalent to <tt>splineView(x, y, 0, 2)</tt>.
    00160         */
    00161     value_type dyy(double x, double y) const
    00162         { return operator()(x, y, 0, 2); }
    00163
    00164         /** Access 3rd derivative in x-direction at real-valued coordinate <tt>(x, y)</tt>.
    00165             Equivalent to <tt>splineView(x, y, 3, 0)</tt>.
    00166         */
    00167     value_type dx3(double x, double y) const
    00168         { return operator()(x, y, 3, 0); }
    00169
    00170         /** Access 3rd derivative in y-direction at real-valued coordinate <tt>(x, y)</tt>.
    00171             Equivalent to <tt>splineView(x, y, 0, 3)</tt>.
    00172         */
    00173     value_type dy3(double x, double y) const
    00174         { return operator()(x, y, 0, 3); }
    00175
    00176         /** Access mixed 3rd derivative dxxy at real-valued coordinate <tt>(x, y)</tt>.
    00177             Equivalent to <tt>splineView(x, y, 2, 1)</tt>.
    00178         */
    00179     value_type dxxy(double x, double y) const
    00180         { return operator()(x, y, 2, 1); }
    00181
    00182         /** Access mixed 3rd derivative dxyy at real-valued coordinate <tt>(x, y)</tt>.
    00183             Equivalent to <tt>splineView(x, y, 1, 2)</tt>.
    00184         */
    00185     value_type dxyy(double x, double y) const
    00186         { return operator()(x, y, 1, 2); }
    00187
    00188         /** Access interpolated function at real-valued coordinate <tt>d</tt>.
    00189             Equivalent to <tt>splineView(d[0], d[1])</tt>.
    00190         */
    00191     value_type operator()(difference_type const & d) const
    00192         { return operator()(d[0], d[1]); }
    00193
    00194         /** Access derivative of order <tt>(dx, dy)</tt> at real-valued coordinate <tt>d</tt>.
    00195             Equivalent to <tt>splineView(d[0], d[1], dx, dy)</tt>.
    00196         */
    00197     value_type operator()(difference_type const & d, unsigned int dx, unsigned int dy) const
    00198         { return operator()(d[0], d[1], dx, dy); }
    00199
    00200         /** Access 1st derivative in x-direction at real-valued coordinate <tt>d</tt>.
    00201             Equivalent to <tt>splineView.dx(d[0], d[1])</tt>.
    00202         */
    00203     value_type dx(difference_type const & d) const
    00204         { return dx(d[0], d[1]); }
    00205
    00206         /** Access 1st derivative in y-direction at real-valued coordinate <tt>d</tt>.
    00207             Equivalent to <tt>splineView.dy(d[0], d[1])</tt>.
    00208         */
    00209     value_type dy(difference_type const & d) const
    00210         { return dy(d[0], d[1]); }
    00211
    00212         /** Access 2nd derivative in x-direction at real-valued coordinate <tt>d</tt>.
    00213             Equivalent to <tt>splineView.dxx(d[0], d[1])</tt>.
    00214         */
    00215     value_type dxx(difference_type const & d) const
    00216         { return dxx(d[0], d[1]); }
    00217
    00218         /** Access mixed 2nd derivative at real-valued coordinate <tt>d</tt>.
    00219             Equivalent to <tt>splineView.dxy(d[0], d[1])</tt>.
    00220         */
    00221     value_type dxy(difference_type const & d) const
    00222         { return dxy(d[0], d[1]); }
    00223
    00224         /** Access 2nd derivative in y-direction at real-valued coordinate <tt>d</tt>.
    00225             Equivalent to <tt>splineView.dyy(d[0], d[1])</tt>.
    00226         */
    00227     value_type dyy(difference_type const & d) const
    00228         { return dyy(d[0], d[1]); }
    00229
    00230         /** Access 3rd derivative in x-direction at real-valued coordinate <tt>d</tt>.
    00231             Equivalent to <tt>splineView.dx3(d[0], d[1])</tt>.
    00232         */
    00233     value_type dx3(difference_type const & d) const
    00234         { return dx3(d[0], d[1]); }
    00235
    00236         /** Access 3rd derivative in y-direction at real-valued coordinate <tt>d</tt>.
    00237             Equivalent to <tt>splineView.dy3(d[0], d[1])</tt>.
    00238         */
    00239     value_type dy3(difference_type const & d) const
    00240         { return dy3(d[0], d[1]); }
    00241
    00242         /** Access mixed 3rd derivative dxxy at real-valued coordinate <tt>d</tt>.
    00243             Equivalent to <tt>splineView.dxxy(d[0], d[1])</tt>.
    00244         */
    00245     value_type dxxy(difference_type const & d) const
    00246         { return dxxy(d[0], d[1]); }
    00247
    00248         /** Access mixed 3rd derivative dxyy at real-valued coordinate <tt>d</tt>.
    00249             Equivalent to <tt>splineView.dxyy(d[0], d[1])</tt>.
    00250         */
    00251     value_type dxyy(difference_type const & d) const
    00252         { return dxyy(d[0], d[1]); }
    00253
    00254         /** Access gradient squared magnitude at real-valued coordinate <tt>(x, y)</tt>.
    00255         */
    00256     value_type g2(double x, double y) const;
    00257
    00258         /** Access 1st derivative in x-direction of gradient squared magnitude
    00259             at real-valued coordinate <tt>(x, y)</tt>.
    00260         */
    00261     value_type g2x(double x, double y) const;
    00262
    00263         /** Access 1st derivative in y-direction of gradient squared magnitude
    00264             at real-valued coordinate <tt>(x, y)</tt>.
    00265         */
    00266     value_type g2y(double x, double y) const;
    00267
    00268         /** Access 2nd derivative in x-direction of gradient squared magnitude
    00269             at real-valued coordinate <tt>(x, y)</tt>.
    00270         */
    00271     value_type g2xx(double x, double y) const;
    00272
    00273         /** Access mixed 2nd derivative of gradient squared magnitude
    00274             at real-valued coordinate <tt>(x, y)</tt>.
    00275         */
    00276     value_type g2xy(double x, double y) const;
    00277
    00278         /** Access 2nd derivative in y-direction of gradient squared magnitude
    00279             at real-valued coordinate <tt>(x, y)</tt>.
    00280         */
    00281     value_type g2yy(double x, double y) const;
    00282
    00283         /** Access gradient squared magnitude at real-valued coordinate <tt>d</tt>.
    00284         */
    00285     value_type g2(difference_type const & d) const
    00286         { return g2(d[0], d[1]); }
    00287
    00288         /** Access 1st derivative in x-direction of gradient squared magnitude
    00289             at real-valued coordinate <tt>d</tt>.
    00290         */
    00291     value_type g2x(difference_type const & d) const
    00292         { return g2x(d[0], d[1]); }
    00293
    00294         /** Access 1st derivative in y-direction of gradient squared magnitude
    00295             at real-valued coordinate <tt>d</tt>.
    00296         */
    00297     value_type g2y(difference_type const & d) const
    00298         { return g2y(d[0], d[1]); }
    00299
    00300         /** Access 2nd derivative in x-direction of gradient squared magnitude
    00301             at real-valued coordinate <tt>d</tt>.
    00302         */
    00303     value_type g2xx(difference_type const & d) const
    00304         { return g2xx(d[0], d[1]); }
    00305
    00306         /** Access mixed 2nd derivative of gradient squared magnitude
    00307             at real-valued coordinate <tt>d</tt>.
    00308         */
    00309     value_type g2xy(difference_type const & d) const
    00310         { return g2xy(d[0], d[1]); }
    00311
    00312         /** Access 2nd derivative in y-direction of gradient squared magnitude
    00313             at real-valued coordinate <tt>d</tt>.
    00314         */
    00315     value_type g2yy(difference_type const & d) const
    00316         { return g2yy(d[0], d[1]); }
    00317
    00318         /** The width of the image.
    00319             <tt>0 &lt;= x &lt;= width()-1</tt> is required for all access functions.
    00320         */
    00321     unsigned int width() const
    00322         { return w_; }
    00323
    00324         /** The height of the image.
    00325             <tt>0 &lt;= y &lt;= height()-1</tt> is required for all access functions.
    00326         */
    00327     unsigned int height() const
    00328         { return h_; }
    00329
    00330         /** The size of the image.
    00331             <tt>0 &lt;= x &lt;= size().x-1</tt> and <tt>0 &lt;= y &lt;= size().y-1</tt>
    00332             are required for all access functions.
    00333         */
    00334     size_type size() const
    00335         { return size_type(w_, h_); }
    00336
    00337         /** The internal image holding the spline coefficients.
    00338         */
    00339     InternalImage const & image() const
    00340     {
    00341         return image_;
    00342     }
    00343
    00344         /** Get the array of polynomial coefficients for the facet containing
    00345             the point <tt>(x, y)</tt>. The array <tt>res</tt> will be resized to
    00346             dimension <tt>(ORDER+1)x(ORDER+1)</tt>. From these coefficients, the
    00347             value of the interpolated function can be calculated by the following
    00348             algorithm
    00349
    00350             \code
    00351             SplineImageView<ORDER, float> view(...);
    00352             double x = ..., y = ...;
    00353             double dx, dy;
    00354
    00355             // calculate the local facet coordinates of x and y
    00356             if(ORDER % 2)
    00357             {
    00358                 // odd order => facet coordinates between 0 and 1
    00359                 dx = x - floor(x);
    00360                 dy = y - floor(y);
    00361             }
    00362             else
    00363             {
    00364                 // even order => facet coordinates between -0.5 and 0.5
    00365                 dx = x - floor(x + 0.5);
    00366                 dy = y - floor(y + 0.5);
    00367             }
    00368
    00369             BasicImage<float> coefficients;
    00370             view.coefficientArray(x, y, coefficients);
    00371
    00372             float f_x_y = 0.0;
    00373             for(int ny = 0; ny < ORDER + 1; ++ny)
    00374                 for(int nx = 0; nx < ORDER + 1; ++nx)
    00375                     f_x_y += pow(dx, nx) * pow(dy, ny) * coefficients(nx, ny);
    00376
    00377             assert(abs(f_x_y - view(x, y)) < 1e-6);
    00378             \endcode
    00379         */
    00380     template <class Array>
    00381     void coefficientArray(double x, double y, Array & res) const;
    00382
    00383         /** Check if x is in the valid range.
    00384             Equivalent to <tt>0 &lt;= x &lt;= width()-1</tt>.
    00385         */
    00386     bool isInsideX(double x) const
    00387     {
    00388         return x >= 0.0 && x <= width()-1.0;
    00389     }
    00390
    00391         /** Check if y is in the valid range.
    00392             Equivalent to <tt>0 &lt;= y &lt;= height()-1</tt>.
    00393         */
    00394     bool isInsideY(double y) const
    00395     {
    00396         return y >= 0.0 && y <= height()-1.0;
    00397     }
    00398
    00399         /** Check if x and y are in the valid range.
    00400             Equivalent to <tt>0 &lt;= x &lt;= width()-1</tt> and <tt>0 &lt;= y &lt;= height()-1</tt>.
    00401         */
    00402     bool isInside(double x, double y) const
    00403     {
    00404         return isInsideX(x) && isInsideY(y);
    00405     }
    00406
    00407         /** Check whether the points <tt>(x0, y0)</tt> and <tt>(x1, y1)</tt> are in
    00408             the same spline facet. For odd order splines, facets span the range
    00409             <tt>(floor(x), floor(x)+1) x (floor(y), floor(y)+1)</tt> (i.e. we have
    00410             integer facet borders), whereas even order splines have facet between
    00411             half integer values
    00412             <tt>(floor(x)-0.5, floor(x)+0.5) x (floor(y)-0.5, floor(y)+0.5)</tt>.
    00413         */
    00414     bool sameFacet(double x0, double y0, double x1, double y1) const
    00415     {
    00416          x0 = VIGRA_CSTD::floor((ORDER % 2) ? x0 : x0 + 0.5);
    00417          y0 = VIGRA_CSTD::floor((ORDER % 2) ? y0 : y0 + 0.5);
    00418          x1 = VIGRA_CSTD::floor((ORDER % 2) ? x1 : x1 + 0.5);
    00419          y1 = VIGRA_CSTD::floor((ORDER % 2) ? y1 : y1 + 0.5);
    00420          return x0 == x1 && y0 == y1;
    00421     }
    00422
    00423   protected:
    00424
    00425     void init();
    00426     void calculateIndices(double x, double y) const;
    00427     void coefficients(double t, double * const & c) const;
    00428     void derivCoefficients(double t, unsigned int d, double * const & c) const;
    00429     value_type convolve() const;
    00430
    00431     unsigned int w_, h_;
    00432     int w1_, h1_;
    00433     double x0_, x1_, y0_, y1_;
    00434     InternalImage image_;
    00435     Spline k_;
    00436     mutable double x_, y_, u_, v_, kx_[ksize_], ky_[ksize_];
    00437     mutable int ix_[ksize_], iy_[ksize_];
    00438 };
    00439
    00440 template <int ORDER, class VALUETYPE>
    00441 void SplineImageView<ORDER, VALUETYPE>::init()
    00442 {
    00443     ArrayVector<double> const & b = k_.prefilterCoefficients();
    00444
    00445     for(unsigned int i=0; i<b.size(); ++i)
    00446     {
    00447         recursiveFilterX(srcImageRange(image_), destImage(image_), b[i], BORDER_TREATMENT_REFLECT);
    00448         recursiveFilterY(srcImageRange(image_), destImage(image_), b[i], BORDER_TREATMENT_REFLECT);
    00449     }
    00450 }
    00451
    00452 namespace detail
    00453 {
    00454
    00455 template <int i>
    00456 struct SplineImageViewUnrollLoop1
    00457 {
    00458     template <class Array>
    00459     static void exec(int c0, Array c)
    00460     {
    00461         SplineImageViewUnrollLoop1<i-1>::exec(c0, c);
    00462         c[i] = c0 + i;
    00463     }
    00464 };
    00465
    00466 template <>
    00467 struct SplineImageViewUnrollLoop1<0>
    00468 {
    00469     template <class Array>
    00470     static void exec(int c0, Array c)
    00471     {
    00472         c[0] = c0;
    00473     }
    00474 };
    00475
    00476 template <int i>
    00477 struct SplineImageViewUnrollLoop2
    00478 {
    00479     template <class Array1, class Image, class Array2>
    00480     static typename Image::value_type
    00481     exec(Array1 k, Image const & img, Array2 x, int y)
    00482     {
    00483         return k[i] * img(x[i], y) + SplineImageViewUnrollLoop2<i-1>::exec(k, img, x, y);
    00484     }
    00485 };
    00486
    00487 template <>
    00488 struct SplineImageViewUnrollLoop2<0>
    00489 {
    00490     template <class Array1, class Image, class Array2>
    00491     static typename Image::value_type
    00492     exec(Array1 k, Image const & img, Array2 x, int y)
    00493     {
    00494         return k[0] * img(x[0], y);
    00495     }
    00496 };
    00497
    00498 } // namespace detail
    00499
    00500 template <int ORDER, class VALUETYPE>
    00501 void
    00502 SplineImageView<ORDER, VALUETYPE>::calculateIndices(double x, double y) const
    00503 {
    00504     if(x == x_ && y == y_)
    00505         return;   // still in cache
    00506
    00507     if(x > x0_ && x < x1_ && y > y0_ && y < y1_)
    00508     {
    00509         detail::SplineImageViewUnrollLoop1<ORDER>::exec(
    00510                                 (ORDER % 2) ? int(x - kcenter_) : int(x + 0.5 - kcenter_), ix_);
    00511         detail::SplineImageViewUnrollLoop1<ORDER>::exec(
    00512                                 (ORDER % 2) ? int(y - kcenter_) : int(y + 0.5 - kcenter_), iy_);
    00513     }
    00514     else
    00515     {
    00516         vigra_precondition(isInside(x, y),
    00517              "SplineImageView<ORDER, VALUETYPE>::calculateIndices(): index out of bounds.");
    00518
    00519         ix_[kcenter_] = (ORDER % 2) ?
    00520                              (int)x :
    00521                              (int)(x + 0.5);
    00522         iy_[kcenter_] = (ORDER % 2) ?
    00523                              (int)y :
    00524                              (int)(y + 0.5);
    00525
    00526         for(int i=0; i<kcenter_; ++i)
    00527         {
    00528             ix_[i] = vigra::abs(ix_[kcenter_] - (kcenter_ - i));
    00529             iy_[i] = vigra::abs(iy_[kcenter_] - (kcenter_ - i));
    00530         }
    00531         for(int i=kcenter_+1; i<ksize_; ++i)
    00532         {
    00533             ix_[i] = w1_ - vigra::abs(w1_ - ix_[kcenter_] - (i - kcenter_));
    00534             iy_[i] = h1_ - vigra::abs(h1_ - iy_[kcenter_] - (i - kcenter_));
    00535         }
    00536     }
    00537     x_ = x;
    00538     y_ = y;
    00539     u_ = x - ix_[kcenter_];
    00540     v_ = y - iy_[kcenter_];
    00541 }
    00542
    00543 template <int ORDER, class VALUETYPE>
    00544 void SplineImageView<ORDER, VALUETYPE>::coefficients(double t, double * const & c) const
    00545 {
    00546     t += kcenter_;
    00547     for(int i = 0; i<ksize_; ++i)
    00548         c[i] = k_(t-i);
    00549 }
    00550
    00551 template <int ORDER, class VALUETYPE>
    00552 void SplineImageView<ORDER, VALUETYPE>::derivCoefficients(double t,
    00553                                                unsigned int d, double * const & c) const
    00554 {
    00555     t += kcenter_;
    00556     for(int i = 0; i<ksize_; ++i)
    00557         c[i] = k_(t-i, d);
    00558 }
    00559
    00560
    00561 template <int ORDER, class VALUETYPE>
    00562 VALUETYPE SplineImageView<ORDER, VALUETYPE>::convolve() const
    00563 {
    00564     InternalValue sum;
    00565     sum = ky_[0]*detail::SplineImageViewUnrollLoop2<ORDER>::exec(kx_, image_, ix_, iy_[0]);
    00566
    00567     for(int j=1; j<ksize_; ++j)
    00568     {
    00569         sum += ky_[j]*detail::SplineImageViewUnrollLoop2<ORDER>::exec(kx_, image_, ix_, iy_[j]);
    00570     }
    00571     return NumericTraits<VALUETYPE>::fromRealPromote(sum);
    00572 }
    00573
    00574 template <int ORDER, class VALUETYPE>
    00575 template <class Array>
    00576 void
    00577 SplineImageView<ORDER, VALUETYPE>::coefficientArray(double x, double y, Array & res) const
    00578 {
    00579     typename Spline::WeightMatrix & weights = Spline::weights();
    00580     InternalValue tmp[ksize_][ksize_];
    00581
    00582     calculateIndices(x, y);
    00583     for(int j=0; j<ksize_; ++j)
    00584     {
    00585         for(int i=0; i<ksize_; ++i)
    00586         {
    00587             tmp[i][j] = 0.0;
    00588             for(int k=0; k<ksize_; ++k)
    00589             {
    00590                 tmp[i][j] += weights[i][k]*image_(ix_[k], iy_[j]);
    00591             }
    00592        }
    00593     }
    00594     res.resize(ksize_, ksize_);
    00595     for(int j=0; j<ksize_; ++j)
    00596     {
    00597         for(int i=0; i<ksize_; ++i)
    00598         {
    00599             res(i,j) = 0.0;
    00600             for(int k=0; k<ksize_; ++k)
    00601             {
    00602                 res(i,j) += weights[j][k]*tmp[i][k];
    00603             }
    00604         }
    00605     }
    00606 }
    00607
    00608 template <int ORDER, class VALUETYPE>
    00609 VALUETYPE SplineImageView<ORDER, VALUETYPE>::operator()(double x, double y) const
    00610 {
    00611     calculateIndices(x, y);
    00612     coefficients(u_, kx_);
    00613     coefficients(v_, ky_);
    00614     return convolve();
    00615 }
    00616
    00617 template <int ORDER, class VALUETYPE>
    00618 VALUETYPE SplineImageView<ORDER, VALUETYPE>::operator()(double x, double y,
    00619                                                  unsigned int dx, unsigned int dy) const
    00620 {
    00621     calculateIndices(x, y);
    00622     derivCoefficients(u_, dx, kx_);
    00623     derivCoefficients(v_, dy, ky_);
    00624     return convolve();
    00625 }
    00626
    00627 template <int ORDER, class VALUETYPE>
    00628 VALUETYPE SplineImageView<ORDER, VALUETYPE>::g2(double x, double y) const
    00629 {
    00630     return sq(dx(x,y)) + sq(dy(x,y));
    00631 }
    00632
    00633 template <int ORDER, class VALUETYPE>
    00634 VALUETYPE SplineImageView<ORDER, VALUETYPE>::g2x(double x, double y) const
    00635 {
    00636     return 2.0*(dx(x,y) * dxx(x,y) + dy(x,y) * dxy(x,y));
    00637 }
    00638
    00639 template <int ORDER, class VALUETYPE>
    00640 VALUETYPE SplineImageView<ORDER, VALUETYPE>::g2y(double x, double y) const
    00641 {
    00642     return 2.0*(dx(x,y) * dxy(x,y) + dy(x,y) * dyy(x,y));
    00643 }
    00644
    00645 template <int ORDER, class VALUETYPE>
    00646 VALUETYPE SplineImageView<ORDER, VALUETYPE>::g2xx(double x, double y) const
    00647 {
    00648     return 2.0*(sq(dxx(x,y)) + dx(x,y) * dx3(x,y) + sq(dxy(x,y)) + dy(x,y) * dxxy(x,y));
    00649 }
    00650
    00651 template <int ORDER, class VALUETYPE>
    00652 VALUETYPE SplineImageView<ORDER, VALUETYPE>::g2yy(double x, double y) const
    00653 {
    00654     return 2.0*(sq(dxy(x,y)) + dx(x,y) * dxyy(x,y) + sq(dyy(x,y)) + dy(x,y) * dy3(x,y));
    00655 }
    00656
    00657 template <int ORDER, class VALUETYPE>
    00658 VALUETYPE SplineImageView<ORDER, VALUETYPE>::g2xy(double x, double y) const
    00659 {
    00660     return 2.0*(dx(x,y) * dxxy(x,y) + dy(x,y) * dxyy(x,y) + dxy(x,y) * (dxx(x,y) + dyy(x,y)));
    00661 }
    00662
    00663 } // namespace vigra
    00664
    00665
    00666 #endif /* VIGRA_SPLINEIMAGEVIEW_HXX */

© Ullrich Köthe (koethe@informatik.uni-hamburg.de) Cognitive Systems Group, University of Hamburg, Germany | html generated using doxygen and Python VIGRA 1.3.2 (27 Jan 2005)