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

details vigra/splineimageview.hxx VIGRA

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)