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

vigra/splineimageview.hxx VIGRA

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)