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