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

# vigra/splineimageview.hxx

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