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

vigra/splineimageview.hxx |

    00001 /************************************************************************/
    00002 /*                                                                      */
    00003 /*               Copyright 1998-2004 by Ullrich Koethe                  */
    00004 /*                                                                      */
    00005 /*    This file is part of the VIGRA computer vision library.           */
    00006 /*    The VIGRA Website is                                              */
    00007 /*        http://hci.iwr.uni-heidelberg.de/vigra/                       */
    00008 /*    Please direct questions, bug reports, and contributions to        */
    00009 /*        ullrich.koethe@iwr.uni-heidelberg.de    or                    */
    00010 /*        vigra@informatik.uni-hamburg.de                               */
    00011 /*                                                                      */
    00012 /*    Permission is hereby granted, free of charge, to any person       */
    00013 /*    obtaining a copy of this software and associated documentation    */
    00014 /*    files (the "Software"), to deal in the Software without           */
    00015 /*    restriction, including without limitation the rights to use,      */
    00016 /*    copy, modify, merge, publish, distribute, sublicense, and/or      */
    00017 /*    sell copies of the Software, and to permit persons to whom the    */
    00018 /*    Software is furnished to do so, subject to the following          */
    00019 /*    conditions:                                                       */
    00020 /*                                                                      */
    00021 /*    The above copyright notice and this permission notice shall be    */
    00022 /*    included in all copies or substantial portions of the             */
    00023 /*    Software.                                                         */
    00024 /*                                                                      */
    00025 /*    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND    */
    00026 /*    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   */
    00027 /*    OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          */
    00028 /*    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT       */
    00029 /*    HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,      */
    00030 /*    WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING      */
    00031 /*    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR     */
    00032 /*    OTHER DEALINGS IN THE SOFTWARE.                                   */
    00033 /*                                                                      */
    00034 /************************************************************************/
    00035
    00036 #ifndef VIGRA_SPLINEIMAGEVIEW_HXX
    00037 #define VIGRA_SPLINEIMAGEVIEW_HXX
    00038
    00039 #include "mathutil.hxx"
    00040 #include "recursiveconvolution.hxx"
    00041 #include "splines.hxx"
    00042 #include "array_vector.hxx"
    00043 #include "basicimage.hxx"
    00044 #include "copyimage.hxx"
    00045 #include "tinyvector.hxx"
    00046 #include "fixedpoint.hxx"
    00047 #include "multi_array.hxx"
    00048
    00049 namespace vigra {
    00050
    00051 /********************************************************/
    00052 /*                                                      */
    00053 /*                    SplineImageView                   */
    00054 /*                                                      */
    00055 /********************************************************/
    00056 /** \brief Create a continuous view onto a discrete image using splines.
    00057
    00058     This class is very useful if image values or derivatives at arbitrary
    00059     real-valued coordinates are needed. Access at such coordinates is implemented by
    00060     interpolating the given discrete image values with a spline of the
    00061     specified <tt>ORDER</TT>. Continuous derivatives are available up to
    00062     degree <tt>ORDER-1</TT>. If the requested coordinates are near the image border,
    00063     reflective boundary conditions are applied. In principle, this class can also be used
    00064     for image resizing, but here the functions from the <tt>resize...</tt> family are
    00065     more efficient, since they exploit the regularity of the sampling grid.
    00066
    00067     The <tt>SplineImageView</tt> template is explicitly specialized to make it as efficient as possible.
    00068     In particular, unnecessary copying of the image is avoided when the iterators passed
    00069     in the constructor originate from a \ref vigra::BasicImage. In addition, these specializations
    00070     provide function <tt>unchecked(...)</tt> that do not perform bounds checking. If the original image
    00071     is not a variant of \ref vigra::BasicImage, one can customize the internal representation by
    00072     using \ref vigra::SplineImageView0 or \ref vigra::SplineImageView1.
    00073
    00074     <b>Usage:</b>
    00075
    00076     <b>\#include</b> <vigra/splineimageview.hxx><br>
    00077     Namespace: vigra
    00078
    00079     \code
    00080     BImage img(w,h);
    00081     ... // fill img
    00082
    00083     // construct spline view for quadratic interpolation
    00084     SplineImageView<2, double> spi2(srcImageRange(img));
    00085
    00086     double x = ..., y = ...;
    00087     double v2 = spi2(x, y);
    00088
    00089     // construct spline view for linear interpolation
    00090     SplineImageView<1, UInt32> spi1(srcImageRange(img));
    00091
    00092     UInt32 v1 = spi1(x, y);
    00093
    00094     FixedPoint<16, 15> fx(...), fy(...);
    00095     UInt32 vf = spi1.unchecked(fx, fy); // caller is sure that (fx, fy) are valid coordinates
    00096     \endcode
    00097 */
    00098 template <int ORDER, class VALUETYPE>
    00099 class SplineImageView
    00100 {
    00101     typedef typename NumericTraits<VALUETYPE>::RealPromote InternalValue;
    00102
    00103   public:
    00104
    00105         /** The view's value type (return type of access and derivative functions).
    00106         */
    00107     typedef VALUETYPE value_type;
    00108
    00109         /** The view's squared norm type (return type of g2() etc.).
    00110         */
    00111     typedef typename NormTraits<VALUETYPE>::SquaredNormType SquaredNormType;
    00112
    00113         /** The view's size type.
    00114         */
    00115     typedef Size2D size_type;
    00116
    00117         /** The view's difference type.
    00118         */
    00119     typedef TinyVector<double, 2> difference_type;
    00120
    00121         /** The order of the spline used.
    00122         */
    00123     enum StaticOrder { order = ORDER };
    00124
    00125         /** The type of the internal image holding the spline coefficients.
    00126         */
    00127     typedef BasicImage<InternalValue> InternalImage;
    00128
    00129   private:
    00130     typedef typename InternalImage::traverser InternalTraverser;
    00131     typedef typename InternalTraverser::row_iterator InternalRowIterator;
    00132     typedef typename InternalTraverser::column_iterator InternalColumnIterator;
    00133     typedef BSpline<ORDER, double> Spline;
    00134
    00135     enum { ksize_ = ORDER + 1, kcenter_ = ORDER / 2 };
    00136
    00137   public:
    00138         /** Construct SplineImageView for the given image.
    00139
    00140             If <tt>skipPrefiltering = true</tt> (default: <tt>false</tt>), the recursive
    00141             prefilter of the cardinal spline function is not applied, resulting
    00142             in an approximating (smoothing) rather than interpolating spline. This is
    00143             especially useful if customized prefilters are to be applied.
    00144         */
    00145     template <class SrcIterator, class SrcAccessor>
    00146     SplineImageView(SrcIterator is, SrcIterator iend, SrcAccessor sa, bool skipPrefiltering = false)
    00147     : w_(iend.x - is.x), h_(iend.y - is.y), w1_(w_-1), h1_(h_-1),
    00148       x0_(kcenter_), x1_(w_ - kcenter_ - 2), y0_(kcenter_), y1_(h_ - kcenter_ - 2),
    00149       image_(w_, h_),
    00150       x_(-1.0), y_(-1.0),
    00151       u_(-1.0), v_(-1.0)
    00152     {
    00153         copyImage(srcIterRange(is, iend, sa), destImage(image_));
    00154         if(!skipPrefiltering)
    00155             init();
    00156     }
    00157
    00158         /** Construct SplineImageView for the given image.
    00159
    00160             If <tt>skipPrefiltering = true</tt> (default: <tt>false</tt>), the recursive
    00161             prefilter of the cardinal spline function is not applied, resulting
    00162             in an approximating (smoothing) rather than interpolating spline. This is
    00163             especially useful if customized prefilters are to be applied.
    00164         */
    00165     template <class SrcIterator, class SrcAccessor>
    00166     SplineImageView(triple<SrcIterator, SrcIterator, SrcAccessor> s, bool skipPrefiltering = false)
    00167     : w_(s.second.x - s.first.x), h_(s.second.y - s.first.y), w1_(w_-1), h1_(h_-1),
    00168       x0_(kcenter_), x1_(w_ - kcenter_ - 2), y0_(kcenter_), y1_(h_ - kcenter_ - 2),
    00169       image_(w_, h_),
    00170       x_(-1.0), y_(-1.0),
    00171       u_(-1.0), v_(-1.0)
    00172     {
    00173         copyImage(srcIterRange(s.first, s.second, s.third), destImage(image_));
    00174         if(!skipPrefiltering)
    00175             init();
    00176     }
    00177
    00178         /** Access interpolated function at real-valued coordinate <tt>(x, y)</tt>.
    00179             If <tt>(x, y)</tt> is near the image border or outside the image, the value
    00180             is calculated with reflective boundary conditions. An exception is thrown if the
    00181             coordinate is outside the first reflection.
    00182         */
    00183     value_type operator()(double x, double y) const;
    00184
    00185         /** Access derivative of order <tt>(dx, dy)</tt> at real-valued coordinate <tt>(x, y)</tt>.
    00186             If <tt>(x, y)</tt> is near the image border or outside the image, the value
    00187             is calculated with reflective boundary conditions. An exception is thrown if the
    00188             coordinate is outside the first reflection.
    00189         */
    00190     value_type operator()(double x, double y, unsigned int dx, unsigned int dy) const;
    00191
    00192         /** Access 1st derivative in x-direction at real-valued coordinate <tt>(x, y)</tt>.
    00193             Equivalent to <tt>splineView(x, y, 1, 0)</tt>.
    00194         */
    00195     value_type dx(double x, double y) const
    00196         { return operator()(x, y, 1, 0); }
    00197
    00198         /** Access 1st derivative in y-direction at real-valued coordinate <tt>(x, y)</tt>.
    00199             Equivalent to <tt>splineView(x, y, 0, 1)</tt>.
    00200         */
    00201     value_type dy(double x, double y) const
    00202         { return operator()(x, y, 0, 1); }
    00203
    00204         /** Access 2nd derivative in x-direction at real-valued coordinate <tt>(x, y)</tt>.
    00205             Equivalent to <tt>splineView(x, y, 2, 0)</tt>.
    00206         */
    00207     value_type dxx(double x, double y) const
    00208         { return operator()(x, y, 2, 0); }
    00209
    00210         /** Access mixed 2nd derivative at real-valued coordinate <tt>(x, y)</tt>.
    00211             Equivalent to <tt>splineView(x, y, 1, 1)</tt>.
    00212         */
    00213     value_type dxy(double x, double y) const
    00214         { return operator()(x, y, 1, 1); }
    00215
    00216         /** Access 2nd derivative in y-direction at real-valued coordinate <tt>(x, y)</tt>.
    00217             Equivalent to <tt>splineView(x, y, 0, 2)</tt>.
    00218         */
    00219     value_type dyy(double x, double y) const
    00220         { return operator()(x, y, 0, 2); }
    00221
    00222         /** Access 3rd derivative in x-direction at real-valued coordinate <tt>(x, y)</tt>.
    00223             Equivalent to <tt>splineView(x, y, 3, 0)</tt>.
    00224         */
    00225     value_type dx3(double x, double y) const
    00226         { return operator()(x, y, 3, 0); }
    00227
    00228         /** Access 3rd derivative in y-direction at real-valued coordinate <tt>(x, y)</tt>.
    00229             Equivalent to <tt>splineView(x, y, 0, 3)</tt>.
    00230         */
    00231     value_type dy3(double x, double y) const
    00232         { return operator()(x, y, 0, 3); }
    00233
    00234         /** Access mixed 3rd derivative dxxy at real-valued coordinate <tt>(x, y)</tt>.
    00235             Equivalent to <tt>splineView(x, y, 2, 1)</tt>.
    00236         */
    00237     value_type dxxy(double x, double y) const
    00238         { return operator()(x, y, 2, 1); }
    00239
    00240         /** Access mixed 3rd derivative dxyy at real-valued coordinate <tt>(x, y)</tt>.
    00241             Equivalent to <tt>splineView(x, y, 1, 2)</tt>.
    00242         */
    00243     value_type dxyy(double x, double y) const
    00244         { return operator()(x, y, 1, 2); }
    00245
    00246         /** Access interpolated function at real-valued coordinate <tt>d</tt>.
    00247             Equivalent to <tt>splineView(d[0], d[1])</tt>.
    00248         */
    00249     value_type operator()(difference_type const & d) const
    00250         { return operator()(d[0], d[1]); }
    00251
    00252         /** Access derivative of order <tt>(dx, dy)</tt> at real-valued coordinate <tt>d</tt>.
    00253             Equivalent to <tt>splineView(d[0], d[1], dx, dy)</tt>.
    00254         */
    00255     value_type operator()(difference_type const & d, unsigned int dx, unsigned int dy) const
    00256         { return operator()(d[0], d[1], dx, dy); }
    00257
    00258         /** Access 1st derivative in x-direction at real-valued coordinate <tt>d</tt>.
    00259             Equivalent to <tt>splineView.dx(d[0], d[1])</tt>.
    00260         */
    00261     value_type dx(difference_type const & d) const
    00262         { return dx(d[0], d[1]); }
    00263
    00264         /** Access 1st derivative in y-direction at real-valued coordinate <tt>d</tt>.
    00265             Equivalent to <tt>splineView.dy(d[0], d[1])</tt>.
    00266         */
    00267     value_type dy(difference_type const & d) const
    00268         { return dy(d[0], d[1]); }
    00269
    00270         /** Access 2nd derivative in x-direction at real-valued coordinate <tt>d</tt>.
    00271             Equivalent to <tt>splineView.dxx(d[0], d[1])</tt>.
    00272         */
    00273     value_type dxx(difference_type const & d) const
    00274         { return dxx(d[0], d[1]); }
    00275
    00276         /** Access mixed 2nd derivative at real-valued coordinate <tt>d</tt>.
    00277             Equivalent to <tt>splineView.dxy(d[0], d[1])</tt>.
    00278         */
    00279     value_type dxy(difference_type const & d) const
    00280         { return dxy(d[0], d[1]); }
    00281
    00282         /** Access 2nd derivative in y-direction at real-valued coordinate <tt>d</tt>.
    00283             Equivalent to <tt>splineView.dyy(d[0], d[1])</tt>.
    00284         */
    00285     value_type dyy(difference_type const & d) const
    00286         { return dyy(d[0], d[1]); }
    00287
    00288         /** Access 3rd derivative in x-direction at real-valued coordinate <tt>d</tt>.
    00289             Equivalent to <tt>splineView.dx3(d[0], d[1])</tt>.
    00290         */
    00291     value_type dx3(difference_type const & d) const
    00292         { return dx3(d[0], d[1]); }
    00293
    00294         /** Access 3rd derivative in y-direction at real-valued coordinate <tt>d</tt>.
    00295             Equivalent to <tt>splineView.dy3(d[0], d[1])</tt>.
    00296         */
    00297     value_type dy3(difference_type const & d) const
    00298         { return dy3(d[0], d[1]); }
    00299
    00300         /** Access mixed 3rd derivative dxxy at real-valued coordinate <tt>d</tt>.
    00301             Equivalent to <tt>splineView.dxxy(d[0], d[1])</tt>.
    00302         */
    00303     value_type dxxy(difference_type const & d) const
    00304         { return dxxy(d[0], d[1]); }
    00305
    00306         /** Access mixed 3rd derivative dxyy at real-valued coordinate <tt>d</tt>.
    00307             Equivalent to <tt>splineView.dxyy(d[0], d[1])</tt>.
    00308         */
    00309     value_type dxyy(difference_type const & d) const
    00310         { return dxyy(d[0], d[1]); }
    00311
    00312         /** Access gradient squared magnitude at real-valued coordinate <tt>(x, y)</tt>.
    00313         */
    00314     SquaredNormType g2(double x, double y) const;
    00315
    00316         /** Access 1st derivative in x-direction of gradient squared magnitude
    00317             at real-valued coordinate <tt>(x, y)</tt>.
    00318         */
    00319     SquaredNormType g2x(double x, double y) const;
    00320
    00321         /** Access 1st derivative in y-direction of gradient squared magnitude
    00322             at real-valued coordinate <tt>(x, y)</tt>.
    00323         */
    00324     SquaredNormType g2y(double x, double y) const;
    00325
    00326         /** Access 2nd derivative in x-direction of gradient squared magnitude
    00327             at real-valued coordinate <tt>(x, y)</tt>.
    00328         */
    00329     SquaredNormType g2xx(double x, double y) const;
    00330
    00331         /** Access mixed 2nd derivative of gradient squared magnitude
    00332             at real-valued coordinate <tt>(x, y)</tt>.
    00333         */
    00334     SquaredNormType g2xy(double x, double y) const;
    00335
    00336         /** Access 2nd derivative in y-direction of gradient squared magnitude
    00337             at real-valued coordinate <tt>(x, y)</tt>.
    00338         */
    00339     SquaredNormType g2yy(double x, double y) const;
    00340
    00341         /** Access gradient squared magnitude at real-valued coordinate <tt>d</tt>.
    00342         */
    00343     SquaredNormType g2(difference_type const & d) const
    00344         { return g2(d[0], d[1]); }
    00345
    00346         /** Access 1st derivative in x-direction of gradient squared magnitude
    00347             at real-valued coordinate <tt>d</tt>.
    00348         */
    00349     SquaredNormType g2x(difference_type const & d) const
    00350         { return g2x(d[0], d[1]); }
    00351
    00352         /** Access 1st derivative in y-direction of gradient squared magnitude
    00353             at real-valued coordinate <tt>d</tt>.
    00354         */
    00355     SquaredNormType g2y(difference_type const & d) const
    00356         { return g2y(d[0], d[1]); }
    00357
    00358         /** Access 2nd derivative in x-direction of gradient squared magnitude
    00359             at real-valued coordinate <tt>d</tt>.
    00360         */
    00361     SquaredNormType g2xx(difference_type const & d) const
    00362         { return g2xx(d[0], d[1]); }
    00363
    00364         /** Access mixed 2nd derivative of gradient squared magnitude
    00365             at real-valued coordinate <tt>d</tt>.
    00366         */
    00367     SquaredNormType g2xy(difference_type const & d) const
    00368         { return g2xy(d[0], d[1]); }
    00369
    00370         /** Access 2nd derivative in y-direction of gradient squared magnitude
    00371             at real-valued coordinate <tt>d</tt>.
    00372         */
    00373     SquaredNormType g2yy(difference_type const & d) const
    00374         { return g2yy(d[0], d[1]); }
    00375
    00376         /** The width of the image.
    00377             <tt>0 <= x <= width()-1</tt> is required for all access functions.
    00378         */
    00379     unsigned int width() const
    00380         { return w_; }
    00381
    00382         /** The height of the image.
    00383             <tt>0 <= y <= height()-1</tt> is required for all access functions.
    00384         */
    00385     unsigned int height() const
    00386         { return h_; }
    00387
    00388         /** The size of the image.
    00389             <tt>0 <= x <= size().x-1</tt> and <tt>0 <= y <= size().y-1</tt>
    00390             are required for all access functions.
    00391         */
    00392     size_type size() const
    00393         { return size_type(w_, h_); }
    00394
    00395         /** The shape of the image.
    00396             Same as size(), except for the return type.
    00397         */
    00398     TinyVector<unsigned int, 2> shape() const
    00399         { return TinyVector<unsigned int, 2>(w_, h_); }
    00400
    00401         /** The internal image holding the spline coefficients.
    00402         */
    00403     InternalImage const & image() const
    00404     {
    00405         return image_;
    00406     }
    00407
    00408         /** Get the array of polynomial coefficients for the facet containing
    00409             the point <tt>(x, y)</tt>. The array <tt>res</tt> must have
    00410             dimension <tt>(ORDER+1)x(ORDER+1)</tt>. From these coefficients, the
    00411             value of the interpolated function can be calculated by the following
    00412             algorithm
    00413
    00414             \code
    00415             SplineImageView<ORDER, float> view(...);
    00416             double x = ..., y = ...;
    00417             double dx, dy;
    00418
    00419             // calculate the local facet coordinates of x and y
    00420             if(ORDER % 2)
    00421             {
    00422                 // odd order => facet coordinates between 0 and 1
    00423                 dx = x - floor(x);
    00424                 dy = y - floor(y);
    00425             }
    00426             else
    00427             {
    00428                 // even order => facet coordinates between -0.5 and 0.5
    00429                 dx = x - floor(x + 0.5);
    00430                 dy = y - floor(y + 0.5);
    00431             }
    00432
    00433             BasicImage<float> coefficients;
    00434             view.coefficientArray(x, y, coefficients);
    00435
    00436             float f_x_y = 0.0;
    00437             for(int ny = 0; ny < ORDER + 1; ++ny)
    00438                 for(int nx = 0; nx < ORDER + 1; ++nx)
    00439                     f_x_y += pow(dx, nx) * pow(dy, ny) * coefficients(nx, ny);
    00440
    00441             assert(abs(f_x_y - view(x, y)) < 1e-6);
    00442             \endcode
    00443         */
    00444     template <class Array>
    00445     void coefficientArray(double x, double y, Array & res) const;
    00446
    00447         /** Check if x is in the original image range.
    00448             Equivalent to <tt>0 <= x <= width()-1</tt>.
    00449         */
    00450     bool isInsideX(double x) const
    00451     {
    00452         return x >= 0.0 && x <= width()-1.0;
    00453     }
    00454
    00455         /** Check if y is in the original image range.
    00456             Equivalent to <tt>0 <= y <= height()-1</tt>.
    00457         */
    00458     bool isInsideY(double y) const
    00459     {
    00460         return y >= 0.0 && y <= height()-1.0;
    00461     }
    00462
    00463         /** Check if x and y are in the original image range.
    00464             Equivalent to <tt>0 <= x <= width()-1</tt> and <tt>0 <= y <= height()-1</tt>.
    00465         */
    00466     bool isInside(double x, double y) const
    00467     {
    00468         return isInsideX(x) && isInsideY(y);
    00469     }
    00470
    00471         /** Check if x and y are in the valid range. Points outside the original image range are computed
    00472             by reflective boundary conditions, but only within the first reflection.
    00473             Equivalent to <tt>-width() + ORDER/2 + 2 < x < 2*width() - ORDER/2 - 2</tt> and
    00474             <tt>-height() + ORDER/2 + 2 < y < 2*height() - ORDER/2 - 2</tt>.
    00475         */
    00476     bool isValid(double x, double y) const
    00477     {
    00478         return x < w1_ + x1_ && x > -x1_ && y < h1_ + y1_ && y > -y1_;
    00479     }
    00480
    00481         /** Check whether the points <tt>(x0, y0)</tt> and <tt>(x1, y1)</tt> are in
    00482             the same spline facet. For odd order splines, facets span the range
    00483             <tt>(floor(x), floor(x)+1) x (floor(y), floor(y)+1)</tt> (i.e. we have
    00484             integer facet borders), whereas even order splines have facet between
    00485             half integer values
    00486             <tt>(floor(x)-0.5, floor(x)+0.5) x (floor(y)-0.5, floor(y)+0.5)</tt>.
    00487         */
    00488     bool sameFacet(double x0, double y0, double x1, double y1) const
    00489     {
    00490          x0 = VIGRA_CSTD::floor((ORDER % 2) ? x0 : x0 + 0.5);
    00491          y0 = VIGRA_CSTD::floor((ORDER % 2) ? y0 : y0 + 0.5);
    00492          x1 = VIGRA_CSTD::floor((ORDER % 2) ? x1 : x1 + 0.5);
    00493          y1 = VIGRA_CSTD::floor((ORDER % 2) ? y1 : y1 + 0.5);
    00494          return x0 == x1 && y0 == y1;
    00495     }
    00496
    00497   protected:
    00498
    00499     void init();
    00500     void calculateIndices(double x, double y) const;
    00501     void coefficients(double t, double * const & c) const;
    00502     void derivCoefficients(double t, unsigned int d, double * const & c) const;
    00503     value_type convolve() const;
    00504
    00505     unsigned int w_, h_;
    00506     int w1_, h1_;
    00507     double x0_, x1_, y0_, y1_;
    00508     InternalImage image_;
    00509     Spline k_;
    00510     mutable double x_, y_, u_, v_, kx_[ksize_], ky_[ksize_];
    00511     mutable int ix_[ksize_], iy_[ksize_];
    00512 };
    00513
    00514 template <int ORDER, class VALUETYPE>
    00515 void SplineImageView<ORDER, VALUETYPE>::init()
    00516 {
    00517     ArrayVector<double> const & b = k_.prefilterCoefficients();
    00518
    00519     for(unsigned int i=0; i<b.size(); ++i)
    00520     {
    00521         recursiveFilterX(srcImageRange(image_), destImage(image_), b[i], BORDER_TREATMENT_REFLECT);
    00522         recursiveFilterY(srcImageRange(image_), destImage(image_), b[i], BORDER_TREATMENT_REFLECT);
    00523     }
    00524 }
    00525
    00526 namespace detail
    00527 {
    00528
    00529 template <int i>
    00530 struct SplineImageViewUnrollLoop1
    00531 {
    00532     template <class Array>
    00533     static void exec(int c0, Array c)
    00534     {
    00535         SplineImageViewUnrollLoop1<i-1>::exec(c0, c);
    00536         c[i] = c0 + i;
    00537     }
    00538 };
    00539
    00540 template <>
    00541 struct SplineImageViewUnrollLoop1<0>
    00542 {
    00543     template <class Array>
    00544     static void exec(int c0, Array c)
    00545     {
    00546         c[0] = c0;
    00547     }
    00548 };
    00549
    00550 template <int i, class ValueType>
    00551 struct SplineImageViewUnrollLoop2
    00552 {
    00553     template <class Array1, class RowIterator, class Array2>
    00554     static ValueType
    00555     exec(Array1 k, RowIterator r, Array2 x)
    00556     {
    00557         return ValueType(k[i] * r[x[i]]) + SplineImageViewUnrollLoop2<i-1, ValueType>::exec(k, r, x);
    00558     }
    00559 };
    00560
    00561 template <class ValueType>
    00562 struct SplineImageViewUnrollLoop2<0, ValueType>
    00563 {
    00564     template <class Array1, class RowIterator, class Array2>
    00565     static ValueType
    00566     exec(Array1 k, RowIterator r, Array2 x)
    00567     {
    00568         return ValueType(k[0] * r[x[0]]);
    00569     }
    00570 };
    00571
    00572 } // namespace detail
    00573
    00574 template <int ORDER, class VALUETYPE>
    00575 void
    00576 SplineImageView<ORDER, VALUETYPE>::calculateIndices(double x, double y) const
    00577 {
    00578     if(x == x_ && y == y_)
    00579         return;   // still in cache
    00580
    00581     if(x > x0_ && x < x1_ && y > y0_ && y < y1_)
    00582     {
    00583         detail::SplineImageViewUnrollLoop1<ORDER>::exec(
    00584                                 (ORDER % 2) ? int(x - kcenter_) : int(x + 0.5 - kcenter_), ix_);
    00585         detail::SplineImageViewUnrollLoop1<ORDER>::exec(
    00586                                 (ORDER % 2) ? int(y - kcenter_) : int(y + 0.5 - kcenter_), iy_);
    00587
    00588         u_ = x - ix_[kcenter_];
    00589         v_ = y - iy_[kcenter_];
    00590     }
    00591     else
    00592     {
    00593         vigra_precondition(isValid(x,y),
    00594                     "SplineImageView::calculateIndices(): coordinates out of range.");
    00595
    00596         int xCenter = (ORDER % 2) ?
    00597                       (int)VIGRA_CSTD::floor(x) :
    00598                       (int)VIGRA_CSTD::floor(x + 0.5);
    00599         int yCenter = (ORDER % 2) ?
    00600                       (int)VIGRA_CSTD::floor(y) :
    00601                       (int)VIGRA_CSTD::floor(y + 0.5);
    00602
    00603         if(x >= x1_)
    00604         {
    00605             for(int i = 0; i < ksize_; ++i)
    00606                 ix_[i] = w1_ - vigra::abs(w1_ - xCenter - (i - kcenter_));
    00607         }
    00608         else
    00609         {
    00610             for(int i = 0; i < ksize_; ++i)
    00611                 ix_[i] = vigra::abs(xCenter - (kcenter_ - i));
    00612         }
    00613         if(y >= y1_)
    00614         {
    00615             for(int i = 0; i < ksize_; ++i)
    00616                 iy_[i] = h1_ - vigra::abs(h1_ - yCenter - (i - kcenter_));
    00617         }
    00618         else
    00619         {
    00620             for(int i = 0; i < ksize_; ++i)
    00621                 iy_[i] = vigra::abs(yCenter - (kcenter_ - i));
    00622         }
    00623         u_ = x - xCenter;
    00624         v_ = y - yCenter;
    00625     }
    00626     x_ = x;
    00627     y_ = y;
    00628 }
    00629
    00630 template <int ORDER, class VALUETYPE>
    00631 void SplineImageView<ORDER, VALUETYPE>::coefficients(double t, double * const & c) const
    00632 {
    00633     t += kcenter_;
    00634     for(int i = 0; i<ksize_; ++i)
    00635         c[i] = k_(t-i);
    00636 }
    00637
    00638 template <int ORDER, class VALUETYPE>
    00639 void SplineImageView<ORDER, VALUETYPE>::derivCoefficients(double t,
    00640                                                unsigned int d, double * const & c) const
    00641 {
    00642     t += kcenter_;
    00643     for(int i = 0; i<ksize_; ++i)
    00644         c[i] = k_(t-i, d);
    00645 }
    00646
    00647 template <int ORDER, class VALUETYPE>
    00648 VALUETYPE SplineImageView<ORDER, VALUETYPE>::convolve() const
    00649 {
    00650     typedef typename NumericTraits<VALUETYPE>::RealPromote RealPromote;
    00651     RealPromote sum;
    00652     sum = RealPromote(
    00653       ky_[0]*detail::SplineImageViewUnrollLoop2<ORDER, RealPromote>::exec(kx_, image_.rowBegin(iy_[0]), ix_));
    00654
    00655     for(int j=1; j<ksize_; ++j)
    00656     {
    00657         sum += RealPromote(
    00658           ky_[j]*detail::SplineImageViewUnrollLoop2<ORDER, RealPromote>::exec(kx_, image_.rowBegin(iy_[j]), ix_));
    00659     }
    00660     return detail::RequiresExplicitCast<VALUETYPE>::cast(sum);
    00661 }
    00662
    00663 template <int ORDER, class VALUETYPE>
    00664 template <class Array>
    00665 void
    00666 SplineImageView<ORDER, VALUETYPE>::coefficientArray(double x, double y, Array & res) const
    00667 {
    00668     typedef typename Array::value_type ResType;
    00669     typename Spline::WeightMatrix & weights = Spline::weights();
    00670     ResType tmp[ksize_][ksize_];
    00671
    00672     calculateIndices(x, y);
    00673     for(int j=0; j<ksize_; ++j)
    00674     {
    00675         for(int i=0; i<ksize_; ++i)
    00676         {
    00677             tmp[i][j] = ResType();
    00678             for(int k=0; k<ksize_; ++k)
    00679             {
    00680                 tmp[i][j] += weights[i][k]*image_(ix_[k], iy_[j]);
    00681             }
    00682        }
    00683     }
    00684     for(int j=0; j<ksize_; ++j)
    00685     {
    00686         for(int i=0; i<ksize_; ++i)
    00687         {
    00688             res(i,j) = ResType();
    00689             for(int k=0; k<ksize_; ++k)
    00690             {
    00691                 res(i,j) += weights[j][k]*tmp[i][k];
    00692             }
    00693         }
    00694     }
    00695 }
    00696
    00697 template <int ORDER, class VALUETYPE>
    00698 VALUETYPE SplineImageView<ORDER, VALUETYPE>::operator()(double x, double y) const
    00699 {
    00700     calculateIndices(x, y);
    00701     coefficients(u_, kx_);
    00702     coefficients(v_, ky_);
    00703     return convolve();
    00704 }
    00705
    00706 template <int ORDER, class VALUETYPE>
    00707 VALUETYPE SplineImageView<ORDER, VALUETYPE>::operator()(double x, double y,
    00708                                                  unsigned int dx, unsigned int dy) const
    00709 {
    00710     calculateIndices(x, y);
    00711     derivCoefficients(u_, dx, kx_);
    00712     derivCoefficients(v_, dy, ky_);
    00713     return convolve();
    00714 }
    00715
    00716 template <int ORDER, class VALUETYPE>
    00717 typename SplineImageView<ORDER, VALUETYPE>::SquaredNormType
    00718 SplineImageView<ORDER, VALUETYPE>::g2(double x, double y) const
    00719 {
    00720     return squaredNorm(dx(x,y)) + squaredNorm(dy(x,y));
    00721 }
    00722
    00723 template <int ORDER, class VALUETYPE>
    00724 typename SplineImageView<ORDER, VALUETYPE>::SquaredNormType
    00725 SplineImageView<ORDER, VALUETYPE>::g2x(double x, double y) const
    00726 {
    00727     return SquaredNormType(2.0)*(dot(dx(x,y), dxx(x,y)) + dot(dy(x,y), dxy(x,y)));
    00728 }
    00729
    00730 template <int ORDER, class VALUETYPE>
    00731 typename SplineImageView<ORDER, VALUETYPE>::SquaredNormType
    00732 SplineImageView<ORDER, VALUETYPE>::g2y(double x, double y) const
    00733 {
    00734     return SquaredNormType(2.0)*(dot(dx(x,y), dxy(x,y)) + dot(dy(x,y), dyy(x,y)));
    00735 }
    00736
    00737 template <int ORDER, class VALUETYPE>
    00738 typename SplineImageView<ORDER, VALUETYPE>::SquaredNormType
    00739 SplineImageView<ORDER, VALUETYPE>::g2xx(double x, double y) const
    00740 {
    00741     return SquaredNormType(2.0)*(squaredNorm(dxx(x,y)) + dot(dx(x,y), dx3(x,y)) +
    00742                                  squaredNorm(dxy(x,y)) + dot(dy(x,y), dxxy(x,y)));
    00743 }
    00744
    00745 template <int ORDER, class VALUETYPE>
    00746 typename SplineImageView<ORDER, VALUETYPE>::SquaredNormType
    00747 SplineImageView<ORDER, VALUETYPE>::g2yy(double x, double y) const
    00748 {
    00749     return SquaredNormType(2.0)*(squaredNorm(dxy(x,y)) + dot(dx(x,y), dxyy(x,y)) +
    00750                                  squaredNorm(dyy(x,y)) + dot(dy(x,y), dy3(x,y)));
    00751 }
    00752
    00753 template <int ORDER, class VALUETYPE>
    00754 typename SplineImageView<ORDER, VALUETYPE>::SquaredNormType
    00755 SplineImageView<ORDER, VALUETYPE>::g2xy(double x, double y) const
    00756 {
    00757     return SquaredNormType(2.0)*(dot(dx(x,y), dxxy(x,y)) + dot(dy(x,y), dxyy(x,y)) +
    00758                                  dot(dxy(x,y), dxx(x,y) + dyy(x,y)));
    00759 }
    00760
    00761 /********************************************************/
    00762 /*                                                      */
    00763 /*                    SplineImageView0                  */
    00764 /*                                                      */
    00765 /********************************************************/
    00766 template <class VALUETYPE, class INTERNAL_INDEXER>
    00767 class SplineImageView0Base
    00768 {
    00769     typedef typename INTERNAL_INDEXER::value_type InternalValue;
    00770   public:
    00771     typedef VALUETYPE value_type;
    00772     typedef typename NormTraits<VALUETYPE>::SquaredNormType SquaredNormType;
    00773     typedef Size2D size_type;
    00774     typedef TinyVector<double, 2> difference_type;
    00775     enum StaticOrder { order = 0 };
    00776
    00777   public:
    00778
    00779     SplineImageView0Base(unsigned int w, unsigned int h)
    00780     : w_(w), h_(h)
    00781     {}
    00782
    00783     SplineImageView0Base(int w, int h, INTERNAL_INDEXER i)
    00784     : w_(w), h_(h), internalIndexer_(i)
    00785     {}
    00786
    00787     template <unsigned IntBits1, unsigned FractionalBits1,
    00788               unsigned IntBits2, unsigned FractionalBits2>
    00789     value_type unchecked(FixedPoint<IntBits1, FractionalBits1> x,
    00790                          FixedPoint<IntBits2, FractionalBits2> y) const
    00791     {
    00792         return internalIndexer_(round(x), round(y));
    00793     }
    00794
    00795     template <unsigned IntBits1, unsigned FractionalBits1,
    00796               unsigned IntBits2, unsigned FractionalBits2>
    00797     value_type unchecked(FixedPoint<IntBits1, FractionalBits1> x,
    00798                          FixedPoint<IntBits2, FractionalBits2> y,
    00799                          unsigned int dx, unsigned int dy) const
    00800     {
    00801         if((dx != 0) || (dy != 0))
    00802             return NumericTraits<VALUETYPE>::zero();
    00803         return unchecked(x, y);
    00804     }
    00805
    00806     value_type unchecked(double x, double y) const
    00807     {
    00808         return internalIndexer_((int)(x + 0.5), (int)(y + 0.5));
    00809     }
    00810
    00811     value_type unchecked(double x, double y, unsigned int dx, unsigned int dy) const
    00812     {
    00813         if((dx != 0) || (dy != 0))
    00814             return NumericTraits<VALUETYPE>::zero();
    00815         return unchecked(x, y);
    00816     }
    00817
    00818     value_type operator()(double x, double y) const
    00819     {
    00820         int ix, iy;
    00821         if(x < 0.0)
    00822         {
    00823             ix = (int)(-x + 0.5);
    00824             vigra_precondition(ix <= (int)w_ - 1,
    00825                     "SplineImageView::operator(): coordinates out of range.");
    00826         }
    00827         else
    00828         {
    00829             ix = (int)(x + 0.5);
    00830             if(ix >= (int)w_)
    00831             {
    00832                 ix = 2*w_-2-ix;
    00833                 vigra_precondition(ix >= 0,
    00834                         "SplineImageView::operator(): coordinates out of range.");
    00835             }
    00836         }
    00837         if(y < 0.0)
    00838         {
    00839             iy = (int)(-y + 0.5);
    00840             vigra_precondition(iy <= (int)h_ - 1,
    00841                     "SplineImageView::operator(): coordinates out of range.");
    00842         }
    00843         else
    00844         {
    00845             iy = (int)(y + 0.5);
    00846             if(iy >= (int)h_)
    00847             {
    00848                 iy = 2*h_-2-iy;
    00849                 vigra_precondition(iy >= 0,
    00850                         "SplineImageView::operator(): coordinates out of range.");
    00851             }
    00852         }
    00853         return internalIndexer_(ix, iy);
    00854     }
    00855
    00856     value_type operator()(double x, double y, unsigned int dx, unsigned int dy) const
    00857     {
    00858         if((dx != 0) || (dy != 0))
    00859             return NumericTraits<VALUETYPE>::zero();
    00860         return operator()(x, y);
    00861     }
    00862
    00863     value_type dx(double x, double y) const
    00864         { return NumericTraits<VALUETYPE>::zero(); }
    00865
    00866     value_type dy(double x, double y) const
    00867         { return NumericTraits<VALUETYPE>::zero(); }
    00868
    00869     value_type dxx(double x, double y) const
    00870         { return NumericTraits<VALUETYPE>::zero(); }
    00871
    00872     value_type dxy(double x, double y) const
    00873         { return NumericTraits<VALUETYPE>::zero(); }
    00874
    00875     value_type dyy(double x, double y) const
    00876         { return NumericTraits<VALUETYPE>::zero(); }
    00877
    00878     value_type dx3(double x, double y) const
    00879         { return NumericTraits<VALUETYPE>::zero(); }
    00880
    00881     value_type dy3(double x, double y) const
    00882         { return NumericTraits<VALUETYPE>::zero(); }
    00883
    00884     value_type dxxy(double x, double y) const
    00885         { return NumericTraits<VALUETYPE>::zero(); }
    00886
    00887     value_type dxyy(double x, double y) const
    00888         { return NumericTraits<VALUETYPE>::zero(); }
    00889
    00890     value_type operator()(difference_type const & d) const
    00891         { return operator()(d[0], d[1]); }
    00892
    00893     value_type operator()(difference_type const & d, unsigned int dx, unsigned int dy) const
    00894         { return operator()(d[0], d[1], dx, dy); }
    00895
    00896     value_type dx(difference_type const & d) const
    00897         { return NumericTraits<VALUETYPE>::zero(); }
    00898
    00899     value_type dy(difference_type const & d) const
    00900         { return NumericTraits<VALUETYPE>::zero(); }
    00901
    00902     value_type dxx(difference_type const & d) const
    00903         { return NumericTraits<VALUETYPE>::zero(); }
    00904
    00905     value_type dxy(difference_type const & d) const
    00906         { return NumericTraits<VALUETYPE>::zero(); }
    00907
    00908     value_type dyy(difference_type const & d) const
    00909         { return NumericTraits<VALUETYPE>::zero(); }
    00910
    00911     value_type dx3(difference_type const & d) const
    00912         { return NumericTraits<VALUETYPE>::zero(); }
    00913
    00914     value_type dy3(difference_type const & d) const
    00915         { return NumericTraits<VALUETYPE>::zero(); }
    00916
    00917     value_type dxxy(difference_type const & d) const
    00918         { return NumericTraits<VALUETYPE>::zero(); }
    00919
    00920     value_type dxyy(difference_type const & d) const
    00921         { return NumericTraits<VALUETYPE>::zero(); }
    00922
    00923     SquaredNormType g2(double x, double y) const
    00924         { return NumericTraits<SquaredNormType>::zero(); }
    00925
    00926     SquaredNormType g2x(double x, double y) const
    00927         { return NumericTraits<SquaredNormType>::zero(); }
    00928
    00929     SquaredNormType g2y(double x, double y) const
    00930         { return NumericTraits<SquaredNormType>::zero(); }
    00931
    00932     SquaredNormType g2xx(double x, double y) const
    00933         { return NumericTraits<SquaredNormType>::zero(); }
    00934
    00935     SquaredNormType g2xy(double x, double y) const
    00936         { return NumericTraits<SquaredNormType>::zero(); }
    00937
    00938     SquaredNormType g2yy(double x, double y) const
    00939         { return NumericTraits<SquaredNormType>::zero(); }
    00940
    00941     SquaredNormType g2(difference_type const & d) const
    00942         { return NumericTraits<SquaredNormType>::zero(); }
    00943
    00944     SquaredNormType g2x(difference_type const & d) const
    00945         { return NumericTraits<SquaredNormType>::zero(); }
    00946
    00947     SquaredNormType g2y(difference_type const & d) const
    00948         { return NumericTraits<SquaredNormType>::zero(); }
    00949
    00950     SquaredNormType g2xx(difference_type const & d) const
    00951         { return NumericTraits<SquaredNormType>::zero(); }
    00952
    00953     SquaredNormType g2xy(difference_type const & d) const
    00954         { return NumericTraits<SquaredNormType>::zero(); }
    00955
    00956     SquaredNormType g2yy(difference_type const & d) const
    00957         { return NumericTraits<SquaredNormType>::zero(); }
    00958
    00959     unsigned int width() const
    00960         { return w_; }
    00961
    00962     unsigned int height() const
    00963         { return h_; }
    00964
    00965     size_type size() const
    00966         { return size_type(w_, h_); }
    00967
    00968     TinyVector<unsigned int, 2> shape() const
    00969         { return TinyVector<unsigned int, 2>(w_, h_); }
    00970
    00971     template <class Array>
    00972     void coefficientArray(double x, double y, Array & res) const
    00973     {
    00974         res(0, 0) = operator()(x,y);
    00975     }
    00976
    00977     bool isInsideX(double x) const
    00978     {
    00979         return x >= 0.0 && x <= width() - 1.0;
    00980     }
    00981
    00982     bool isInsideY(double y) const
    00983     {
    00984         return y >= 0.0 && y <= height() - 1.0;
    00985     }
    00986
    00987     bool isInside(double x, double y) const
    00988     {
    00989         return isInsideX(x) && isInsideY(y);
    00990     }
    00991
    00992     bool isValid(double x, double y) const
    00993     {
    00994         return x < 2.0*w_-2.0 && x > 1.0-w_ && y < 2.0*h_-2.0 && y > 1.0-h_;
    00995     }
    00996
    00997     bool sameFacet(double x0, double y0, double x1, double y1) const
    00998     {
    00999          x0 = VIGRA_CSTD::floor(x0 + 0.5);
    01000          y0 = VIGRA_CSTD::floor(y0 + 0.5);
    01001          x1 = VIGRA_CSTD::floor(x1 + 0.5);
    01002          y1 = VIGRA_CSTD::floor(y1 + 0.5);
    01003          return x0 == x1 && y0 == y1;
    01004     }
    01005
    01006   protected:
    01007     unsigned int w_, h_;
    01008     INTERNAL_INDEXER internalIndexer_;
    01009 };
    01010
    01011 /** \brief Create an image view for nearest-neighbor interpolation.
    01012
    01013     This class behaves like \ref vigra::SplineImageView&lt;0, ...&gt;, but one can pass
    01014     an additional template argument that determined the internal representation of the image.
    01015     If this is equal to the argument type passed in the constructor, the image is not copied.
    01016     By default, this works for \ref vigra::BasicImage, \ref vigra::BasicImageView,
    01017     \ref vigra::MultiArray&lt;2, ...&gt;, and \ref vigra::MultiArrayView&lt;2, ...&gt;.
    01018
    01019 */
    01020 template <class VALUETYPE, class INTERNAL_TRAVERSER = typename BasicImage<VALUETYPE>::const_traverser>
    01021 class SplineImageView0
    01022 : public SplineImageView0Base<VALUETYPE, INTERNAL_TRAVERSER>
    01023 {
    01024     typedef SplineImageView0Base<VALUETYPE, INTERNAL_TRAVERSER> Base;
    01025   public:
    01026     typedef typename Base::value_type value_type;
    01027     typedef typename Base::SquaredNormType SquaredNormType;
    01028     typedef typename Base::size_type size_type;
    01029     typedef typename Base::difference_type difference_type;
    01030     enum StaticOrder { order = Base::order };
    01031     typedef BasicImage<VALUETYPE> InternalImage;
    01032
    01033   protected:
    01034     typedef typename IteratorTraits<INTERNAL_TRAVERSER>::mutable_iterator InternalTraverser;
    01035     typedef typename IteratorTraits<InternalTraverser>::DefaultAccessor InternalAccessor;
    01036     typedef typename IteratorTraits<INTERNAL_TRAVERSER>::const_iterator InternalConstTraverser;
    01037     typedef typename IteratorTraits<InternalConstTraverser>::DefaultAccessor InternalConstAccessor;
    01038
    01039   public:
    01040
    01041         /* when traverser and accessor types passed to the constructor are the same as the corresponding
    01042            internal types, we need not copy the image (speed up)
    01043         */
    01044     SplineImageView0(InternalTraverser is, InternalTraverser iend, InternalAccessor sa)
    01045     : Base(iend.x - is.x, iend.y - is.y, is)
    01046     {}
    01047
    01048     SplineImageView0(triple<InternalTraverser, InternalTraverser, InternalAccessor> s)
    01049     : Base(s.second.x - s.first.x, s.second.y - s.first.y, s.first)
    01050     {}
    01051
    01052     SplineImageView0(InternalConstTraverser is, InternalConstTraverser iend, InternalConstAccessor sa)
    01053     : Base(iend.x - is.x, iend.y - is.y, is)
    01054     {}
    01055
    01056     SplineImageView0(triple<InternalConstTraverser, InternalConstTraverser, InternalConstAccessor> s)
    01057     : Base(s.second.x - s.first.x, s.second.y - s.first.y, s.first)
    01058     {}
    01059
    01060     template<class T, class SU>
    01061     SplineImageView0(MultiArrayView<2, T, SU> const & i)
    01062     : Base(i.shape(0), i.shape(1)),
    01063       image_(i.shape(0), i.shape(1))
    01064     {
    01065         for(unsigned int y=0; y<this->height(); ++y)
    01066             for(unsigned int x=0; x<this->width(); ++x)
    01067                 image_(x,y) = detail::RequiresExplicitCast<VALUETYPE>::cast(i(x,y));
    01068         this->internalIndexer_ = image_.upperLeft();
    01069     }
    01070
    01071     template <class SrcIterator, class SrcAccessor>
    01072     SplineImageView0(SrcIterator is, SrcIterator iend, SrcAccessor sa)
    01073     : Base(iend.x - is.x, iend.y - is.y),
    01074       image_(iend - is)
    01075     {
    01076         copyImage(srcIterRange(is, iend, sa), destImage(image_));
    01077         this->internalIndexer_ = image_.upperLeft();
    01078     }
    01079
    01080     template <class SrcIterator, class SrcAccessor>
    01081     SplineImageView0(triple<SrcIterator, SrcIterator, SrcAccessor> s)
    01082     : Base(s.second.x - s.first.x, s.second.y - s.first.y),
    01083       image_(s.second - s.first)
    01084     {
    01085         copyImage(s, destImage(image_));
    01086         this->internalIndexer_ = image_.upperLeft();
    01087     }
    01088
    01089     InternalImage const & image() const
    01090         { return image_; }
    01091
    01092   protected:
    01093     InternalImage image_;
    01094 };
    01095
    01096 template <class VALUETYPE, class StridedOrUnstrided>
    01097 class SplineImageView0<VALUETYPE, MultiArrayView<2, VALUETYPE, StridedOrUnstrided> >
    01098 : public SplineImageView0Base<VALUETYPE, MultiArrayView<2, VALUETYPE, StridedOrUnstrided> >
    01099 {
    01100     typedef SplineImageView0Base<VALUETYPE, MultiArrayView<2, VALUETYPE, StridedOrUnstrided> > Base;
    01101   public:
    01102     typedef typename Base::value_type value_type;
    01103     typedef typename Base::SquaredNormType SquaredNormType;
    01104     typedef typename Base::size_type size_type;
    01105     typedef typename Base::difference_type difference_type;
    01106     enum StaticOrder { order = Base::order };
    01107     typedef BasicImage<VALUETYPE> InternalImage;
    01108
    01109   protected:
    01110     typedef MultiArrayView<2, VALUETYPE, StridedOrUnstrided> InternalIndexer;
    01111
    01112   public:
    01113
    01114         /* when traverser and accessor types passed to the constructor are the same as the corresponding
    01115            internal types, we need not copy the image (speed up)
    01116         */
    01117     SplineImageView0(InternalIndexer const & i)
    01118     : Base(i.shape(0), i.shape(1), i)
    01119     {}
    01120
    01121     template<class T, class SU>
    01122     SplineImageView0(MultiArrayView<2, T, SU> const & i)
    01123     : Base(i.shape(0), i.shape(1)),
    01124       image_(i.shape(0), i.shape(1))
    01125     {
    01126         for(unsigned int y=0; y<this->height(); ++y)
    01127             for(unsigned int x=0; x<this->width(); ++x)
    01128                 image_(x,y) = detail::RequiresExplicitCast<VALUETYPE>::cast(i(x,y));
    01129         this->internalIndexer_ = InternalIndexer(typename InternalIndexer::difference_type(this->width(), this->height()),
    01130                                                  image_.data());
    01131     }
    01132
    01133     template <class SrcIterator, class SrcAccessor>
    01134     SplineImageView0(SrcIterator is, SrcIterator iend, SrcAccessor sa)
    01135     : Base(iend.x - is.x, iend.y - is.y),
    01136       image_(iend-is)
    01137     {
    01138         copyImage(srcIterRange(is, iend, sa), destImage(image_));
    01139         this->internalIndexer_ = InternalIndexer(typename InternalIndexer::difference_type(this->width(), this->height()),
    01140                                                  image_.data());
    01141     }
    01142
    01143     template <class SrcIterator, class SrcAccessor>
    01144     SplineImageView0(triple<SrcIterator, SrcIterator, SrcAccessor> s)
    01145     : Base(s.second.x - s.first.x, s.second.y - s.first.y),
    01146       image_(s.second - s.first)
    01147     {
    01148         copyImage(s, destImage(image_));
    01149         this->internalIndexer_ = InternalIndexer(typename InternalIndexer::difference_type(this->width(), this->height()),
    01150                                                  image_.data());
    01151     }
    01152
    01153     InternalImage const & image() const
    01154         { return image_; }
    01155
    01156   protected:
    01157     InternalImage image_;
    01158 };
    01159
    01160 template <class VALUETYPE>
    01161 class SplineImageView<0, VALUETYPE>
    01162 : public SplineImageView0<VALUETYPE>
    01163 {
    01164     typedef SplineImageView0<VALUETYPE> Base;
    01165   public:
    01166     typedef typename Base::value_type value_type;
    01167     typedef typename Base::SquaredNormType SquaredNormType;
    01168     typedef typename Base::size_type size_type;
    01169     typedef typename Base::difference_type difference_type;
    01170     enum StaticOrder { order = Base::order };
    01171     typedef typename Base::InternalImage InternalImage;
    01172
    01173   protected:
    01174     typedef typename Base::InternalTraverser InternalTraverser;
    01175     typedef typename Base::InternalAccessor InternalAccessor;
    01176     typedef typename Base::InternalConstTraverser InternalConstTraverser;
    01177     typedef typename Base::InternalConstAccessor InternalConstAccessor;
    01178
    01179 public:
    01180
    01181         /* when traverser and accessor types passed to the constructor are the same as the corresponding
    01182            internal types, we need not copy the image (speed up)
    01183         */
    01184     SplineImageView(InternalTraverser is, InternalTraverser iend, InternalAccessor sa, bool /* unused */ = false)
    01185     : Base(is, iend, sa)
    01186     {}
    01187
    01188     SplineImageView(triple<InternalTraverser, InternalTraverser, InternalAccessor> s, bool /* unused */ = false)
    01189     : Base(s)
    01190     {}
    01191
    01192     SplineImageView(InternalConstTraverser is, InternalConstTraverser iend, InternalConstAccessor sa, bool /* unused */ = false)
    01193     : Base(is, iend, sa)
    01194     {}
    01195
    01196     SplineImageView(triple<InternalConstTraverser, InternalConstTraverser, InternalConstAccessor> s, bool /* unused */ = false)
    01197     : Base(s)
    01198     {}
    01199
    01200     template <class SrcIterator, class SrcAccessor>
    01201     SplineImageView(SrcIterator is, SrcIterator iend, SrcAccessor sa, bool /* unused */ = false)
    01202     : Base(is, iend, sa)
    01203     {
    01204         copyImage(srcIterRange(is, iend, sa), destImage(this->image_));
    01205     }
    01206
    01207     template <class SrcIterator, class SrcAccessor>
    01208     SplineImageView(triple<SrcIterator, SrcIterator, SrcAccessor> s, bool /* unused */ = false)
    01209     : Base(s)
    01210     {
    01211         copyImage(s, destImage(this->image_));
    01212     }
    01213 };
    01214
    01215 /********************************************************/
    01216 /*                                                      */
    01217 /*                    SplineImageView1                  */
    01218 /*                                                      */
    01219 /********************************************************/
    01220 template <class VALUETYPE, class INTERNAL_INDEXER>
    01221 class SplineImageView1Base
    01222 {
    01223     typedef typename INTERNAL_INDEXER::value_type InternalValue;
    01224   public:
    01225     typedef VALUETYPE value_type;
    01226     typedef Size2D size_type;
    01227     typedef typename NormTraits<VALUETYPE>::SquaredNormType SquaredNormType;
    01228     typedef TinyVector<double, 2> difference_type;
    01229     enum StaticOrder { order = 1 };
    01230
    01231   public:
    01232
    01233     SplineImageView1Base(unsigned int w, unsigned int h)
    01234     : w_(w), h_(h)
    01235     {}
    01236
    01237     SplineImageView1Base(int w, int h, INTERNAL_INDEXER i)
    01238     : w_(w), h_(h), internalIndexer_(i)
    01239     {}
    01240
    01241     template <unsigned IntBits1, unsigned FractionalBits1,
    01242               unsigned IntBits2, unsigned FractionalBits2>
    01243     value_type unchecked(FixedPoint<IntBits1, FractionalBits1> x,
    01244                          FixedPoint<IntBits2, FractionalBits2> y) const
    01245     {
    01246         int ix = floor(x);
    01247         FixedPoint<0, FractionalBits1> tx = frac(x - FixedPoint<IntBits1, FractionalBits1>(ix));
    01248         FixedPoint<0, FractionalBits1> dtx = dual_frac(tx);
    01249         if(ix == (int)w_ - 1)
    01250         {
    01251             --ix;
    01252             tx.value = FixedPoint<0, FractionalBits1>::ONE;
    01253             dtx.value = 0;
    01254         }
    01255         int iy = floor(y);
    01256         FixedPoint<0, FractionalBits2> ty = frac(y - FixedPoint<IntBits2, FractionalBits2>(iy));
    01257         FixedPoint<0, FractionalBits2> dty = dual_frac(ty);
    01258         if(iy == (int)h_ - 1)
    01259         {
    01260             --iy;
    01261             ty.value = FixedPoint<0, FractionalBits2>::ONE;
    01262             dty.value = 0;
    01263         }
    01264         return fixed_point_cast<value_type>(
    01265                     dty*(dtx*fixedPoint(internalIndexer_(ix,iy)) +
    01266                                    tx*fixedPoint(internalIndexer_(ix+1,iy))) +
    01267                     ty *(dtx*fixedPoint(internalIndexer_(ix,iy+1)) +
    01268                                    tx*fixedPoint(internalIndexer_(ix+1,iy+1))));
    01269     }
    01270
    01271     template <unsigned IntBits1, unsigned FractionalBits1,
    01272               unsigned IntBits2, unsigned FractionalBits2>
    01273     value_type unchecked(FixedPoint<IntBits1, FractionalBits1> x,
    01274                          FixedPoint<IntBits2, FractionalBits2> y,
    01275                          unsigned int dx, unsigned int dy) const
    01276     {
    01277         int ix = floor(x);
    01278         FixedPoint<0, FractionalBits1> tx = frac(x - FixedPoint<IntBits1, FractionalBits1>(ix));
    01279         FixedPoint<0, FractionalBits1> dtx = dual_frac(tx);
    01280         if(ix == (int)w_ - 1)
    01281         {
    01282             --ix;
    01283             tx.value = FixedPoint<0, FractionalBits1>::ONE;
    01284             dtx.value = 0;
    01285         }
    01286         int iy = floor(y);
    01287         FixedPoint<0, FractionalBits2> ty = frac(y - FixedPoint<IntBits2, FractionalBits2>(iy));
    01288         FixedPoint<0, FractionalBits2> dty = dual_frac(ty);
    01289         if(iy == (int)h_ - 1)
    01290         {
    01291             --iy;
    01292             ty.value = FixedPoint<0, FractionalBits2>::ONE;
    01293             dty.value = 0;
    01294         }
    01295         switch(dx)
    01296         {
    01297           case 0:
    01298               switch(dy)
    01299               {
    01300                 case 0:
    01301                     return fixed_point_cast<value_type>(
    01302                                 dty*(dtx*fixedPoint(internalIndexer_(ix,iy)) +
    01303                                                tx*fixedPoint(internalIndexer_(ix+1,iy))) +
    01304                                 ty *(dtx*fixedPoint(internalIndexer_(ix,iy+1)) +
    01305                                                tx*fixedPoint(internalIndexer_(ix+1,iy+1))));
    01306                 case 1:
    01307                     return fixed_point_cast<value_type>(
    01308                            (dtx*fixedPoint(internalIndexer_(ix,iy+1)) + tx*fixedPoint(internalIndexer_(ix+1,iy+1))) -
    01309                            (dtx*fixedPoint(internalIndexer_(ix,iy)) + tx*fixedPoint(internalIndexer_(ix+1,iy))));
    01310                 default:
    01311                     return NumericTraits<VALUETYPE>::zero();
    01312               }
    01313           case 1:
    01314               switch(dy)
    01315               {
    01316                 case 0:
    01317                     return fixed_point_cast<value_type>(
    01318                                 dty*(fixedPoint(internalIndexer_(ix+1,iy)) - fixedPoint(internalIndexer_(ix,iy))) +
    01319                                 ty *(fixedPoint(internalIndexer_(ix+1,iy+1)) - fixedPoint(internalIndexer_(ix,iy+1))));
    01320                 case 1:
    01321                     return detail::RequiresExplicitCast<value_type>::cast(
    01322                                 (internalIndexer_(ix+1,iy+1) - internalIndexer_(ix,iy+1)) -
    01323                                 (internalIndexer_(ix+1,iy) - internalIndexer_(ix,iy)));
    01324                 default:
    01325                     return NumericTraits<VALUETYPE>::zero();
    01326               }
    01327           default:
    01328               return NumericTraits<VALUETYPE>::zero();
    01329         }
    01330     }
    01331
    01332     value_type unchecked(double x, double y) const
    01333     {
    01334         int ix = (int)std::floor(x);
    01335         if(ix == (int)w_ - 1)
    01336             --ix;
    01337         double tx = x - ix;
    01338         int iy = (int)std::floor(y);
    01339         if(iy == (int)h_ - 1)
    01340             --iy;
    01341         double ty = y - iy;
    01342         return NumericTraits<value_type>::fromRealPromote(
    01343                    (1.0-ty)*((1.0-tx)*internalIndexer_(ix,iy) + tx*internalIndexer_(ix+1,iy)) +
    01344                     ty *((1.0-tx)*internalIndexer_(ix,iy+1) + tx*internalIndexer_(ix+1,iy+1)));
    01345     }
    01346
    01347     value_type unchecked(double x, double y, unsigned int dx, unsigned int dy) const
    01348     {
    01349         int ix = (int)std::floor(x);
    01350         if(ix == (int)w_ - 1)
    01351             --ix;
    01352         double tx = x - ix;
    01353         int iy = (int)std::floor(y);
    01354         if(iy == (int)h_ - 1)
    01355             --iy;
    01356         double ty = y - iy;
    01357         switch(dx)
    01358         {
    01359           case 0:
    01360               switch(dy)
    01361               {
    01362                 case 0:
    01363                     return detail::RequiresExplicitCast<value_type>::cast(
    01364                                (1.0-ty)*((1.0-tx)*internalIndexer_(ix,iy) + tx*internalIndexer_(ix+1,iy)) +
    01365                                 ty *((1.0-tx)*internalIndexer_(ix,iy+1) + tx*internalIndexer_(ix+1,iy+1)));
    01366                 case 1:
    01367                     return detail::RequiresExplicitCast<value_type>::cast(
    01368                                ((1.0-tx)*internalIndexer_(ix,iy+1) + tx*internalIndexer_(ix+1,iy+1)) -
    01369                                ((1.0-tx)*internalIndexer_(ix,iy) + tx*internalIndexer_(ix+1,iy)));
    01370                 default:
    01371                     return NumericTraits<VALUETYPE>::zero();
    01372               }
    01373           case 1:
    01374               switch(dy)
    01375               {
    01376                 case 0:
    01377                     return detail::RequiresExplicitCast<value_type>::cast(
    01378                                (1.0-ty)*(internalIndexer_(ix+1,iy) - internalIndexer_(ix,iy)) +
    01379                                 ty *(internalIndexer_(ix+1,iy+1) - internalIndexer_(ix,iy+1)));
    01380                 case 1:
    01381                     return detail::RequiresExplicitCast<value_type>::cast(
    01382                               (internalIndexer_(ix+1,iy+1) - internalIndexer_(ix,iy+1)) -
    01383                               (internalIndexer_(ix+1,iy) - internalIndexer_(ix,iy)));
    01384                 default:
    01385                     return NumericTraits<VALUETYPE>::zero();
    01386               }
    01387           default:
    01388               return NumericTraits<VALUETYPE>::zero();
    01389         }
    01390     }
    01391
    01392     value_type operator()(double x, double y) const
    01393     {
    01394         return operator()(x, y, 0, 0);
    01395     }
    01396
    01397     value_type operator()(double x, double y, unsigned int dx, unsigned int dy) const
    01398     {
    01399         value_type mul = NumericTraits<value_type>::one();
    01400         if(x < 0.0)
    01401         {
    01402             x = -x;
    01403             vigra_precondition(x <= w_ - 1.0,
    01404                     "SplineImageView::operator(): coordinates out of range.");
    01405             if(dx % 2)
    01406                 mul = -mul;
    01407         }
    01408         else if(x > w_ - 1.0)
    01409         {
    01410             x = 2.0*w_-2.0-x;
    01411             vigra_precondition(x >= 0.0,
    01412                     "SplineImageView::operator(): coordinates out of range.");
    01413             if(dx % 2)
    01414                 mul = -mul;
    01415         }
    01416         if(y < 0.0)
    01417         {
    01418             y = -y;
    01419             vigra_precondition(y <= h_ - 1.0,
    01420                     "SplineImageView::operator(): coordinates out of range.");
    01421             if(dy % 2)
    01422                 mul = -mul;
    01423         }
    01424         else if(y > h_ - 1.0)
    01425         {
    01426             y = 2.0*h_-2.0-y;
    01427             vigra_precondition(y >= 0.0,
    01428                     "SplineImageView::operator(): coordinates out of range.");
    01429             if(dy % 2)
    01430                 mul = -mul;
    01431         }
    01432         return mul*unchecked(x, y, dx, dy);
    01433     }
    01434
    01435     value_type dx(double x, double y) const
    01436         { return operator()(x, y, 1, 0); }
    01437
    01438     value_type dy(double x, double y) const
    01439         { return operator()(x, y, 0, 1); }
    01440
    01441     value_type dxx(double x, double y) const
    01442         { return NumericTraits<VALUETYPE>::zero(); }
    01443
    01444     value_type dxy(double x, double y) const
    01445         { return operator()(x, y, 1, 1); }
    01446
    01447     value_type dyy(double x, double y) const
    01448         { return NumericTraits<VALUETYPE>::zero(); }
    01449
    01450     value_type dx3(double x, double y) const
    01451         { return NumericTraits<VALUETYPE>::zero(); }
    01452
    01453     value_type dy3(double x, double y) const
    01454         { return NumericTraits<VALUETYPE>::zero(); }
    01455
    01456     value_type dxxy(double x, double y) const
    01457         { return NumericTraits<VALUETYPE>::zero(); }
    01458
    01459     value_type dxyy(double x, double y) const
    01460         { return NumericTraits<VALUETYPE>::zero(); }
    01461
    01462     value_type operator()(difference_type const & d) const
    01463         { return operator()(d[0], d[1]); }
    01464
    01465     value_type operator()(difference_type const & d, unsigned int dx, unsigned int dy) const
    01466         { return operator()(d[0], d[1], dx, dy); }
    01467
    01468     value_type dx(difference_type const & d) const
    01469         { return operator()(d[0], d[1], 1, 0); }
    01470
    01471     value_type dy(difference_type const & d) const
    01472         { return operator()(d[0], d[1], 0, 1); }
    01473
    01474     value_type dxx(difference_type const & d) const
    01475         { return NumericTraits<VALUETYPE>::zero(); }
    01476
    01477     value_type dxy(difference_type const & d) const
    01478         { return operator()(d[0], d[1], 1, 1); }
    01479
    01480     value_type dyy(difference_type const & d) const
    01481         { return NumericTraits<VALUETYPE>::zero(); }
    01482
    01483     value_type dx3(difference_type const & d) const
    01484         { return NumericTraits<VALUETYPE>::zero(); }
    01485
    01486     value_type dy3(difference_type const & d) const
    01487         { return NumericTraits<VALUETYPE>::zero(); }
    01488
    01489     value_type dxxy(difference_type const & d) const
    01490         { return NumericTraits<VALUETYPE>::zero(); }
    01491
    01492     value_type dxyy(difference_type const & d) const
    01493         { return NumericTraits<VALUETYPE>::zero(); }
    01494
    01495     SquaredNormType g2(double x, double y) const
    01496         { return squaredNorm(dx(x,y)) + squaredNorm(dy(x,y)); }
    01497
    01498     SquaredNormType g2x(double x, double y) const
    01499         { return NumericTraits<SquaredNormType>::zero(); }
    01500
    01501     SquaredNormType g2y(double x, double y) const
    01502         { return NumericTraits<SquaredNormType>::zero(); }
    01503
    01504     SquaredNormType g2xx(double x, double y) const
    01505         { return NumericTraits<SquaredNormType>::zero(); }
    01506
    01507     SquaredNormType g2xy(double x, double y) const
    01508         { return NumericTraits<SquaredNormType>::zero(); }
    01509
    01510     SquaredNormType g2yy(double x, double y) const
    01511         { return NumericTraits<SquaredNormType>::zero(); }
    01512
    01513     SquaredNormType g2(difference_type const & d) const
    01514         { return g2(d[0], d[1]); }
    01515
    01516     SquaredNormType g2x(difference_type const & d) const
    01517         { return NumericTraits<SquaredNormType>::zero(); }
    01518
    01519     SquaredNormType g2y(difference_type const & d) const
    01520         { return NumericTraits<SquaredNormType>::zero(); }
    01521
    01522     SquaredNormType g2xx(difference_type const & d) const
    01523         { return NumericTraits<SquaredNormType>::zero(); }
    01524
    01525     SquaredNormType g2xy(difference_type const & d) const
    01526         { return NumericTraits<SquaredNormType>::zero(); }
    01527
    01528     SquaredNormType g2yy(difference_type const & d) const
    01529         { return NumericTraits<SquaredNormType>::zero(); }
    01530
    01531     unsigned int width() const
    01532         { return w_; }
    01533
    01534     unsigned int height() const
    01535         { return h_; }
    01536
    01537     size_type size() const
    01538         { return size_type(w_, h_); }
    01539
    01540     TinyVector<unsigned int, 2> shape() const
    01541         { return TinyVector<unsigned int, 2>(w_, h_); }
    01542
    01543     template <class Array>
    01544     void coefficientArray(double x, double y, Array & res) const;
    01545
    01546     void calculateIndices(double x, double y, int & ix, int & iy, int & ix1, int & iy1) const;
    01547
    01548     bool isInsideX(double x) const
    01549     {
    01550         return x >= 0.0 && x <= width() - 1.0;
    01551     }
    01552
    01553     bool isInsideY(double y) const
    01554     {
    01555         return y >= 0.0 && y <= height() - 1.0;
    01556     }
    01557
    01558     bool isInside(double x, double y) const
    01559     {
    01560         return isInsideX(x) && isInsideY(y);
    01561     }
    01562
    01563     bool isValid(double x, double y) const
    01564     {
    01565         return x < 2.0*w_-2.0 && x > 1.0-w_ && y < 2.0*h_-2.0 && y > 1.0-h_;
    01566     }
    01567
    01568     bool sameFacet(double x0, double y0, double x1, double y1) const
    01569     {
    01570          x0 = VIGRA_CSTD::floor(x0);
    01571          y0 = VIGRA_CSTD::floor(y0);
    01572          x1 = VIGRA_CSTD::floor(x1);
    01573          y1 = VIGRA_CSTD::floor(y1);
    01574          return x0 == x1 && y0 == y1;
    01575     }
    01576
    01577   protected:
    01578     unsigned int w_, h_;
    01579     INTERNAL_INDEXER internalIndexer_;
    01580 };
    01581
    01582 template <class VALUETYPE, class INTERNAL_INDEXER>
    01583 template <class Array>
    01584 void SplineImageView1Base<VALUETYPE, INTERNAL_INDEXER>::coefficientArray(double x, double y, Array & res) const
    01585 {
    01586     int ix, iy, ix1, iy1;
    01587     calculateIndices(x, y, ix, iy, ix1, iy1);
    01588     res(0,0) = internalIndexer_(ix,iy);
    01589     res(1,0) = internalIndexer_(ix1,iy) - internalIndexer_(ix,iy);
    01590     res(0,1) = internalIndexer_(ix,iy1) - internalIndexer_(ix,iy);
    01591     res(1,1) = internalIndexer_(ix,iy) - internalIndexer_(ix1,iy) -
    01592                internalIndexer_(ix,iy1) + internalIndexer_(ix1,iy1);
    01593 }
    01594
    01595 template <class VALUETYPE, class INTERNAL_INDEXER>
    01596 void SplineImageView1Base<VALUETYPE, INTERNAL_INDEXER>::calculateIndices(double x, double y, int & ix, int & iy, int & ix1, int & iy1) const
    01597 {
    01598     if(x < 0.0)
    01599     {
    01600         x = -x;
    01601         vigra_precondition(x <= w_ - 1.0,
    01602                 "SplineImageView::calculateIndices(): coordinates out of range.");
    01603         ix = (int)VIGRA_CSTD::ceil(x);
    01604         ix1 = ix - 1;
    01605     }
    01606     else if(x >= w_ - 1.0)
    01607     {
    01608         x = 2.0*w_-2.0-x;
    01609         vigra_precondition(x > 0.0,
    01610                 "SplineImageView::calculateIndices(): coordinates out of range.");
    01611         ix = (int)VIGRA_CSTD::ceil(x);
    01612         ix1 = ix - 1;
    01613     }
    01614     else
    01615     {
    01616         ix = (int)VIGRA_CSTD::floor(x);
    01617         ix1 = ix + 1;
    01618     }
    01619     if(y < 0.0)
    01620     {
    01621         y = -y;
    01622         vigra_precondition(y <= h_ - 1.0,
    01623                 "SplineImageView::calculateIndices(): coordinates out of range.");
    01624         iy = (int)VIGRA_CSTD::ceil(y);
    01625         iy1 = iy - 1;
    01626     }
    01627     else if(y >= h_ - 1.0)
    01628     {
    01629         y = 2.0*h_-2.0-y;
    01630         vigra_precondition(y > 0.0,
    01631                 "SplineImageView::calculateIndices(): coordinates out of range.");
    01632         iy = (int)VIGRA_CSTD::ceil(y);
    01633         iy1 = iy - 1;
    01634     }
    01635     else
    01636     {
    01637         iy = (int)VIGRA_CSTD::floor(y);
    01638         iy1 = iy + 1;
    01639     }
    01640 }
    01641
    01642 /** \brief Create an image view for bi-linear interpolation.
    01643
    01644     This class behaves like \ref vigra::SplineImageView&lt;1, ...&gt;, but one can pass
    01645     an additional template argument that determined the internal representation of the image.
    01646     If this is equal to the argument type passed in the constructor, the image is not copied.
    01647     By default, this works for \ref vigra::BasicImage, \ref vigra::BasicImageView,
    01648     \ref vigra::MultiArray&lt;2, ...&gt;, and \ref vigra::MultiArrayView&lt;2, ...&gt;.
    01649
    01650     In addition to the function provided by  \ref vigra::SplineImageView, there are functions
    01651     <tt>unchecked(x,y)</tt> and <tt>unchecked(x,y, xorder, yorder)</tt> which improve speed by
    01652     not applying bounds checking and reflective border treatment (<tt>isInside(x, y)</tt> must
    01653     be <tt>true</tt>), but otherwise behave identically to their checked counterparts.
    01654     In addition, <tt>x</tt> and <tt>y</tt> can have type \ref vigra::FixedPoint instead of
    01655     <tt>double</tt>.
    01656 */
    01657 template <class VALUETYPE, class INTERNAL_TRAVERSER = typename BasicImage<VALUETYPE>::const_traverser>
    01658 class SplineImageView1
    01659 : public SplineImageView1Base<VALUETYPE, INTERNAL_TRAVERSER>
    01660 {
    01661     typedef SplineImageView1Base<VALUETYPE, INTERNAL_TRAVERSER> Base;
    01662   public:
    01663     typedef typename Base::value_type value_type;
    01664     typedef typename Base::SquaredNormType SquaredNormType;
    01665     typedef typename Base::size_type size_type;
    01666     typedef typename Base::difference_type difference_type;
    01667     enum StaticOrder { order = Base::order };
    01668     typedef BasicImage<VALUETYPE> InternalImage;
    01669
    01670   protected:
    01671     typedef typename IteratorTraits<INTERNAL_TRAVERSER>::mutable_iterator InternalTraverser;
    01672     typedef typename IteratorTraits<InternalTraverser>::DefaultAccessor InternalAccessor;
    01673     typedef typename IteratorTraits<INTERNAL_TRAVERSER>::const_iterator InternalConstTraverser;
    01674     typedef typename IteratorTraits<InternalConstTraverser>::DefaultAccessor InternalConstAccessor;
    01675
    01676   public:
    01677
    01678         /* when traverser and accessor types passed to the constructor are the same as the corresponding
    01679            internal types, we need not copy the image (speed up)
    01680         */
    01681     SplineImageView1(InternalTraverser is, InternalTraverser iend, InternalAccessor sa)
    01682     : Base(iend.x - is.x, iend.y - is.y, is)
    01683     {}
    01684
    01685     SplineImageView1(triple<InternalTraverser, InternalTraverser, InternalAccessor> s)
    01686     : Base(s.second.x - s.first.x, s.second.y - s.first.y, s.first)
    01687     {}
    01688
    01689     SplineImageView1(InternalConstTraverser is, InternalConstTraverser iend, InternalConstAccessor sa)
    01690     : Base(iend.x - is.x, iend.y - is.y, is)
    01691     {}
    01692
    01693     SplineImageView1(triple<InternalConstTraverser, InternalConstTraverser, InternalConstAccessor> s)
    01694     : Base(s.second.x - s.first.x, s.second.y - s.first.y, s.first)
    01695     {}
    01696
    01697     template<class T, class SU>
    01698     SplineImageView1(MultiArrayView<2, T, SU> const & i)
    01699     : Base(i.shape(0), i.shape(1)),
    01700       image_(i.shape(0), i.shape(1))
    01701     {
    01702         for(unsigned int y=0; y<this->height(); ++y)
    01703             for(unsigned int x=0; x<this->width(); ++x)
    01704                 image_(x,y) = detail::RequiresExplicitCast<VALUETYPE>::cast(i(x,y));
    01705         this->internalIndexer_ = image_.upperLeft();
    01706     }
    01707
    01708     template <class SrcIterator, class SrcAccessor>
    01709     SplineImageView1(SrcIterator is, SrcIterator iend, SrcAccessor sa)
    01710     : Base(iend.x - is.x, iend.y - is.y),
    01711       image_(iend - is)
    01712     {
    01713         copyImage(srcIterRange(is, iend, sa), destImage(image_));
    01714         this->internalIndexer_ = image_.upperLeft();
    01715     }
    01716
    01717     template <class SrcIterator, class SrcAccessor>
    01718     SplineImageView1(triple<SrcIterator, SrcIterator, SrcAccessor> s)
    01719     : Base(s.second.x - s.first.x, s.second.y - s.first.y),
    01720       image_(s.second - s.first)
    01721     {
    01722         copyImage(s, destImage(image_));
    01723         this->internalIndexer_ = image_.upperLeft();
    01724     }
    01725
    01726     InternalImage const & image() const
    01727         { return image_; }
    01728
    01729   protected:
    01730     InternalImage image_;
    01731 };
    01732
    01733 template <class VALUETYPE, class StridedOrUnstrided>
    01734 class SplineImageView1<VALUETYPE, MultiArrayView<2, VALUETYPE, StridedOrUnstrided> >
    01735 : public SplineImageView1Base<VALUETYPE, MultiArrayView<2, VALUETYPE, StridedOrUnstrided> >
    01736 {
    01737     typedef SplineImageView1Base<VALUETYPE, MultiArrayView<2, VALUETYPE, StridedOrUnstrided> > Base;
    01738   public:
    01739     typedef typename Base::value_type value_type;
    01740     typedef typename Base::SquaredNormType SquaredNormType;
    01741     typedef typename Base::size_type size_type;
    01742     typedef typename Base::difference_type difference_type;
    01743     enum StaticOrder { order = Base::order };
    01744     typedef BasicImage<VALUETYPE> InternalImage;
    01745
    01746   protected:
    01747     typedef MultiArrayView<2, VALUETYPE, StridedOrUnstrided> InternalIndexer;
    01748
    01749   public:
    01750
    01751         /* when traverser and accessor types passed to the constructor are the same as the corresponding
    01752            internal types, we need not copy the image (speed up)
    01753         */
    01754     SplineImageView1(InternalIndexer const & i)
    01755     : Base(i.shape(0), i.shape(1), i)
    01756     {}
    01757
    01758     template<class T, class SU>
    01759     SplineImageView1(MultiArrayView<2, T, SU> const & i)
    01760     : Base(i.shape(0), i.shape(1)),
    01761       image_(i.shape(0), i.shape(1))
    01762     {
    01763         for(unsigned int y=0; y<this->height(); ++y)
    01764             for(unsigned int x=0; x<this->width(); ++x)
    01765                 image_(x,y) = detail::RequiresExplicitCast<VALUETYPE>::cast(i(x,y));
    01766         this->internalIndexer_ = InternalIndexer(typename InternalIndexer::difference_type(this->width(), this->height()),
    01767                                                  image_.data());
    01768     }
    01769
    01770     template <class SrcIterator, class SrcAccessor>
    01771     SplineImageView1(SrcIterator is, SrcIterator iend, SrcAccessor sa)
    01772     : Base(iend.x - is.x, iend.y - is.y),
    01773       image_(iend-is)
    01774     {
    01775         copyImage(srcIterRange(is, iend, sa), destImage(image_));
    01776         this->internalIndexer_ = InternalIndexer(typename InternalIndexer::difference_type(this->width(), this->height()),
    01777                                                  image_.data());
    01778     }
    01779
    01780     template <class SrcIterator, class SrcAccessor>
    01781     SplineImageView1(triple<SrcIterator, SrcIterator, SrcAccessor> s)
    01782     : Base(s.second.x - s.first.x, s.second.y - s.first.y),
    01783       image_(s.second - s.first)
    01784     {
    01785         copyImage(s, destImage(image_));
    01786         this->internalIndexer_ = InternalIndexer(typename InternalIndexer::difference_type(this->width(), this->height()),
    01787                                                  image_.data());
    01788     }
    01789
    01790     InternalImage const & image() const
    01791         { return image_; }
    01792
    01793   protected:
    01794     InternalImage image_;
    01795 };
    01796
    01797 template <class VALUETYPE>
    01798 class SplineImageView<1, VALUETYPE>
    01799 : public SplineImageView1<VALUETYPE>
    01800 {
    01801     typedef SplineImageView1<VALUETYPE> Base;
    01802   public:
    01803     typedef typename Base::value_type value_type;
    01804     typedef typename Base::SquaredNormType SquaredNormType;
    01805     typedef typename Base::size_type size_type;
    01806     typedef typename Base::difference_type difference_type;
    01807     enum StaticOrder { order = Base::order };
    01808     typedef typename Base::InternalImage InternalImage;
    01809
    01810   protected:
    01811     typedef typename Base::InternalTraverser InternalTraverser;
    01812     typedef typename Base::InternalAccessor InternalAccessor;
    01813     typedef typename Base::InternalConstTraverser InternalConstTraverser;
    01814     typedef typename Base::InternalConstAccessor InternalConstAccessor;
    01815
    01816 public:
    01817
    01818         /* when traverser and accessor types passed to the constructor are the same as the corresponding
    01819            internal types, we need not copy the image (speed up)
    01820         */
    01821     SplineImageView(InternalTraverser is, InternalTraverser iend, InternalAccessor sa, bool /* unused */ = false)
    01822     : Base(is, iend, sa)
    01823     {}
    01824
    01825     SplineImageView(triple<InternalTraverser, InternalTraverser, InternalAccessor> s, bool /* unused */ = false)
    01826     : Base(s)
    01827     {}
    01828
    01829     SplineImageView(InternalConstTraverser is, InternalConstTraverser iend, InternalConstAccessor sa, bool /* unused */ = false)
    01830     : Base(is, iend, sa)
    01831     {}
    01832
    01833     SplineImageView(triple<InternalConstTraverser, InternalConstTraverser, InternalConstAccessor> s, bool /* unused */ = false)
    01834     : Base(s)
    01835     {}
    01836
    01837     template <class SrcIterator, class SrcAccessor>
    01838     SplineImageView(SrcIterator is, SrcIterator iend, SrcAccessor sa, bool /* unused */ = false)
    01839     : Base(is, iend, sa)
    01840     {
    01841         copyImage(srcIterRange(is, iend, sa), destImage(this->image_));
    01842     }
    01843
    01844     template <class SrcIterator, class SrcAccessor>
    01845     SplineImageView(triple<SrcIterator, SrcIterator, SrcAccessor> s, bool /* unused */ = false)
    01846     : Base(s)
    01847     {
    01848         copyImage(s, destImage(this->image_));
    01849     }
    01850 };
    01851
    01852 } // namespace vigra
    01853
    01854
    01855 #endif /* VIGRA_SPLINEIMAGEVIEW_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.9.0 (Tue Nov 6 2012)