vigra - vigra/resizeimage.hxx Source File

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

vigra/resizeimage.hxx |

    00001 /************************************************************************/
    00002 /*                                                                      */
    00003 /*               Copyright 1998-2004 by Ullrich Koethe                  */
    00004 /*                                                                      */
    00005 /*    This file is part of the VIGRA computer vision library.           */
    00006 /*    The VIGRA Website is                                              */
    00007 /*        http://hci.iwr.uni-heidelberg.de/vigra/                       */
    00008 /*    Please direct questions, bug reports, and contributions to        */
    00009 /*        ullrich.koethe@iwr.uni-heidelberg.de    or                    */
    00010 /*        vigra@informatik.uni-hamburg.de                               */
    00011 /*                                                                      */
    00012 /*    Permission is hereby granted, free of charge, to any person       */
    00013 /*    obtaining a copy of this software and associated documentation    */
    00014 /*    files (the "Software"), to deal in the Software without           */
    00015 /*    restriction, including without limitation the rights to use,      */
    00016 /*    copy, modify, merge, publish, distribute, sublicense, and/or      */
    00017 /*    sell copies of the Software, and to permit persons to whom the    */
    00018 /*    Software is furnished to do so, subject to the following          */
    00019 /*    conditions:                                                       */
    00020 /*                                                                      */
    00021 /*    The above copyright notice and this permission notice shall be    */
    00022 /*    included in all copies or substantial portions of the             */
    00023 /*    Software.                                                         */
    00024 /*                                                                      */
    00025 /*    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND    */
    00026 /*    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   */
    00027 /*    OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          */
    00028 /*    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT       */
    00029 /*    HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,      */
    00030 /*    WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING      */
    00031 /*    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR     */
    00032 /*    OTHER DEALINGS IN THE SOFTWARE.                                   */
    00033 /*                                                                      */
    00034 /************************************************************************/
    00035
    00036
    00037 #ifndef VIGRA_RESIZEIMAGE_HXX
    00038 #define VIGRA_RESIZEIMAGE_HXX
    00039
    00040 #include <vector>
    00041 #include "utilities.hxx"
    00042 #include "numerictraits.hxx"
    00043 #include "stdimage.hxx"
    00044 #include "recursiveconvolution.hxx"
    00045 #include "separableconvolution.hxx"
    00046 #include "resampling_convolution.hxx"
    00047 #include "splines.hxx"
    00048
    00049 namespace vigra {
    00050
    00051 /*****************************************************************/
    00052 /*                                                               */
    00053 /*                         CoscotFunction                        */
    00054 /*                                                               */
    00055 /*****************************************************************/
    00056
    00057 /*! The Coscot interpolation function.
    00058
    00059     Implements the Coscot interpolation function proposed by Maria Magnusson Seger
    00060     (maria@isy.liu.se) in the context of tomographic reconstruction. It provides a fast
    00061     transition between the pass- and stop-bands and minimal ripple outside the transition
    00062     region. Both properties are important for this application and can be tuned by the parameters
    00063     <i>m</i> and <i>h</i> (with defaults 3 and 0.5). The function is defined by
    00064
    00065     \f[ f_{m,h}(x) = \left\{ \begin{array}{ll}
    00066                                    \frac{1}{2m}\sin(\pi x)\cot(\pi x / (2 m))(h + (1-h)\cos(\pi x/m)) & |x| \leq m \\
    00067                                   0 & \mbox{otherwise}
    00068                         \end{array}\right.
    00069     \f]
    00070
    00071     It can be used as a functor, and as a kernel for
    00072     \ref resamplingConvolveImage() to create a differentiable interpolant
    00073     of an image.
    00074
    00075     <b>\#include</b> <vigra/resizeimage.hxx><br>
    00076     Namespace: vigra
    00077
    00078     \ingroup MathFunctions
    00079 */
    00080 template <class T>
    00081 class CoscotFunction
    00082 {
    00083   public:
    00084
    00085         /** the kernel's value type
    00086         */
    00087     typedef T            value_type;
    00088         /** the unary functor's argument type
    00089         */
    00090     typedef T            argument_type;
    00091         /** the splines polynomial order
    00092         */
    00093     typedef T            result_type;
    00094
    00095     CoscotFunction(unsigned int m = 3, double h = 0.5)
    00096     : m_(m),
    00097       h_(h)
    00098     {}
    00099
    00100         /** function (functor) call
    00101         */
    00102     result_type operator()(argument_type x) const
    00103     {
    00104         return x == 0.0 ?
    00105                     1.0
    00106                   : abs(x) < m_ ?
    00107                         VIGRA_CSTD::sin(M_PI*x) / VIGRA_CSTD::tan(M_PI * x / 2.0 / m_) *
    00108                              (h_ + (1.0 - h_) * VIGRA_CSTD::cos(M_PI * x / m_)) / 2.0 / m_
    00109                       : 0.0;
    00110     }
    00111
    00112         /** index operator -- same as operator()
    00113         */
    00114     value_type operator[](value_type x) const
    00115         { return operator()(x); }
    00116
    00117         /** Radius of the function's support.
    00118             Needed for  \ref resamplingConvolveImage(), equals m.
    00119         */
    00120     double radius() const
    00121         { return m_; }
    00122
    00123         /** Derivative order of the function: always 0.
    00124         */
    00125     unsigned int derivativeOrder() const
    00126         { return 0; }
    00127
    00128         /** Prefilter coefficients for compatibility with \ref vigra::BSpline.
    00129             (array has zero length, since prefiltering is not necessary).
    00130         */
    00131     ArrayVector<double> const & prefilterCoefficients() const
    00132     {
    00133         static ArrayVector<double> b;
    00134         return b;
    00135     }
    00136
    00137   protected:
    00138
    00139     unsigned int m_;
    00140     double h_;
    00141 };
    00142
    00143 /** \addtogroup GeometricTransformations Geometric Transformations
    00144     Zoom up and down by repeating pixels, or using various interpolation schemes.
    00145
    00146     See also: \ref resamplingConvolveImage(), \ref resampleImage(), \ref resizeMultiArraySplineInterpolation()
    00147
    00148     <b>\#include</b> <vigra/stdimagefunctions.hxx><br>
    00149     <b>or</b><br>
    00150     <b>\#include</b> <vigra/resizeimage.hxx><br>
    00151 */
    00152 //@{
    00153
    00154 /********************************************************/
    00155 /*                                                      */
    00156 /*               resizeLineNoInterpolation              */
    00157 /*                                                      */
    00158 /********************************************************/
    00159
    00160 template <class SrcIterator, class SrcAccessor,
    00161           class DestIterator, class DestAccessor>
    00162 void
    00163 resizeLineNoInterpolation(SrcIterator i1, SrcIterator iend, SrcAccessor as,
    00164                            DestIterator id, DestIterator idend, DestAccessor ad)
    00165 {
    00166     int wold = iend - i1;
    00167     int wnew = idend - id;
    00168
    00169     if(wnew == 1)
    00170     {
    00171         ad.set(as(i1), id);
    00172         return;
    00173     }
    00174
    00175     double dx = (double)(wold - 1) / (wnew - 1);
    00176     double x = 0.5;
    00177     for(; id != idend; ++id, x += dx)
    00178     {
    00179         int ix = (int)x;
    00180         ad.set(as(i1, ix), id);
    00181     }
    00182 }
    00183
    00184 /********************************************************/
    00185 /*                                                      */
    00186 /*              resizeImageNoInterpolation              */
    00187 /*                                                      */
    00188 /********************************************************/
    00189
    00190 /** \brief Resize image by repeating the nearest pixel values.
    00191
    00192     This algorithm is very fast and does not require any arithmetic on
    00193     the pixel types.
    00194
    00195     The range of both the input and output images (resp. regions) must
    00196     be given. Both images must have a size of at least 2x2 pixels. The
    00197     scaling factors are then calculated accordingly. Destination
    00198     pixels are directly copied from the appropriate source pixels.
    00199
    00200     The function uses accessors.
    00201
    00202     <b> Declarations:</b>
    00203
    00204     pass arguments explicitly:
    00205     \code
    00206     namespace vigra {
    00207         template <class SrcImageIterator, class SrcAccessor,
    00208                   class DestImageIterator, class DestAccessor>
    00209         void
    00210         resizeImageNoInterpolation(
    00211               SrcImageIterator is, SrcImageIterator iend, SrcAccessor sa,
    00212               DestImageIterator id, DestImageIterator idend, DestAccessor da)
    00213     }
    00214     \endcode
    00215
    00216
    00217     use argument objects in conjunction with \ref ArgumentObjectFactories :
    00218     \code
    00219     namespace vigra {
    00220         template <class SrcImageIterator, class SrcAccessor,
    00221                   class DestImageIterator, class DestAccessor>
    00222         void
    00223         resizeImageNoInterpolation(
    00224               triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
    00225               triple<DestImageIterator, DestImageIterator, DestAccessor> dest)
    00226     }
    00227     \endcode
    00228
    00229     <b> Usage:</b>
    00230
    00231         <b>\#include</b> <vigra/resizeimage.hxx><br>
    00232         Namespace: vigra
    00233
    00234     \code
    00235     vigra::resizeImageNoInterpolation(
    00236                src.upperLeft(), src.lowerRight(), src.accessor(),
    00237                dest.upperLeft(), dest.lowerRight(), dest.accessor());
    00238
    00239     \endcode
    00240
    00241     <b> Required Interface:</b>
    00242
    00243     \code
    00244     SrcImageIterator src_upperleft, src_lowerright;
    00245     DestImageIterator dest_upperleft, src_lowerright;
    00246
    00247     SrcAccessor src_accessor;
    00248     DestAccessor dest_accessor;
    00249
    00250     dest_accessor.set(src_accessor(src_upperleft), dest_upperleft);
    00251
    00252     \endcode
    00253
    00254     <b> Preconditions:</b>
    00255
    00256     \code
    00257     src_lowerright.x - src_upperleft.x > 1
    00258     src_lowerright.y - src_upperleft.y > 1
    00259     dest_lowerright.x - dest_upperleft.x > 1
    00260     dest_lowerright.y - dest_upperleft.y > 1
    00261     \endcode
    00262
    00263 */
    00264 doxygen_overloaded_function(template <...> void resizeImageNoInterpolation)
    00265
    00266 template <class SrcIterator, class SrcAccessor,
    00267           class DestIterator, class DestAccessor>
    00268 void
    00269 resizeImageNoInterpolation(SrcIterator is, SrcIterator iend, SrcAccessor sa,
    00270                       DestIterator id, DestIterator idend, DestAccessor da)
    00271 {
    00272     int w = iend.x - is.x;
    00273     int h = iend.y - is.y;
    00274
    00275     int wnew = idend.x - id.x;
    00276     int hnew = idend.y - id.y;
    00277
    00278     vigra_precondition((w > 1) && (h > 1),
    00279                  "resizeImageNoInterpolation(): "
    00280                  "Source image too small.\n");
    00281     vigra_precondition((wnew > 1) && (hnew > 1),
    00282                  "resizeImageNoInterpolation(): "
    00283                  "Destination image too small.\n");
    00284
    00285     typedef BasicImage<typename SrcAccessor::value_type> TmpImage;
    00286     typedef typename TmpImage::traverser TmpImageIterator;
    00287
    00288     TmpImage tmp(w, hnew);
    00289
    00290     TmpImageIterator yt = tmp.upperLeft();
    00291
    00292     for(int x=0; x<w; ++x, ++is.x, ++yt.x)
    00293     {
    00294         typename SrcIterator::column_iterator c1 = is.columnIterator();
    00295         typename TmpImageIterator::column_iterator ct = yt.columnIterator();
    00296
    00297         resizeLineNoInterpolation(c1, c1 + h, sa, ct, ct + hnew, tmp.accessor());
    00298     }
    00299
    00300     yt = tmp.upperLeft();
    00301
    00302     for(int y=0; y < hnew; ++y, ++yt.y, ++id.y)
    00303     {
    00304         typename DestIterator::row_iterator rd = id.rowIterator();
    00305         typename TmpImageIterator::row_iterator rt = yt.rowIterator();
    00306
    00307         resizeLineNoInterpolation(rt, rt + w, tmp.accessor(), rd, rd + wnew, da);
    00308     }
    00309 }
    00310
    00311 template <class SrcIterator, class SrcAccessor,
    00312           class DestIterator, class DestAccessor>
    00313 inline
    00314 void
    00315 resizeImageNoInterpolation(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    00316                            triple<DestIterator, DestIterator, DestAccessor> dest)
    00317 {
    00318     resizeImageNoInterpolation(src.first, src.second, src.third,
    00319                                    dest.first, dest.second, dest.third);
    00320 }
    00321
    00322 /********************************************************/
    00323 /*                                                      */
    00324 /*             resizeLineLinearInterpolation            */
    00325 /*                                                      */
    00326 /********************************************************/
    00327
    00328 template <class SrcIterator, class SrcAccessor,
    00329           class DestIterator, class DestAccessor>
    00330 void
    00331 resizeLineLinearInterpolation(SrcIterator i1, SrcIterator iend, SrcAccessor as,
    00332                            DestIterator id, DestIterator idend, DestAccessor ad)
    00333 {
    00334     int wold = iend - i1;
    00335     int wnew = idend - id;
    00336
    00337     if((wold <= 1) || (wnew <= 1)) return; // oder error ?
    00338
    00339     typedef
    00340         NumericTraits<typename DestAccessor::value_type> DestTraits;
    00341     typedef typename DestTraits::RealPromote RealPromote;
    00342
    00343     ad.set(DestTraits::fromRealPromote(as(i1)), id);
    00344     ++id;
    00345
    00346     --iend, --idend;
    00347     ad.set(DestTraits::fromRealPromote(as(iend)), idend);
    00348
    00349     double dx = (double)(wold - 1) / (wnew - 1);
    00350     double x = dx;
    00351
    00352     for(; id != idend; ++id, x += dx)
    00353     {
    00354         if(x >= 1.0)
    00355         {
    00356             int xx = (int)x;
    00357             i1 += xx;
    00358             x -= (double)xx;
    00359         }
    00360         double x1 = 1.0 - x;
    00361
    00362         ad.set(DestTraits::fromRealPromote(RealPromote(x1 * as(i1) + x * as(i1, 1))), id);
    00363     }
    00364 }
    00365
    00366 /********************************************************/
    00367 /*                                                      */
    00368 /*           resizeImageLinearInterpolation             */
    00369 /*                                                      */
    00370 /********************************************************/
    00371
    00372 /** \brief Resize image using linear interpolation.
    00373
    00374     The function uses the standard separable bilinear interpolation algorithm to
    00375     obtain a good compromise between quality and speed.
    00376
    00377     The range must of both the input and output images (resp. regions)
    00378     must be given. Both images must have a size of at
    00379     least 2x2. The scaling factors are then calculated
    00380     accordingly. If the source image is larger than the destination, it
    00381     is smoothed (band limited) using a recursive
    00382     exponential filter. The source value_type (SrcAccessor::value_type) must
    00383     be a linear space, i.e. it must support addition, multiplication
    00384     with a scalar real number and \ref NumericTraits "NumericTraits".
    00385     The function uses accessors.
    00386
    00387     <b> Declarations:</b>
    00388
    00389     pass arguments explicitly:
    00390     \code
    00391     namespace vigra {
    00392         template <class SrcImageIterator, class SrcAccessor,
    00393                   class DestImageIterator, class DestAccessor>
    00394         void
    00395         resizeImageLinearInterpolation(
    00396               SrcImageIterator is, SrcImageIterator iend, SrcAccessor sa,
    00397               DestImageIterator id, DestImageIterator idend, DestAccessor da)
    00398     }
    00399     \endcode
    00400
    00401
    00402     use argument objects in conjunction with \ref ArgumentObjectFactories :
    00403     \code
    00404     namespace vigra {
    00405         template <class SrcImageIterator, class SrcAccessor,
    00406                   class DestImageIterator, class DestAccessor>
    00407         void
    00408         resizeImageLinearInterpolation(
    00409               triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
    00410               triple<DestImageIterator, DestImageIterator, DestAccessor> dest)
    00411     }
    00412     \endcode
    00413
    00414     <b> Usage:</b>
    00415
    00416         <b>\#include</b> <vigra/resizeimage.hxx><br>
    00417         Namespace: vigra
    00418
    00419     \code
    00420     vigra::resizeImageLinearInterpolation(
    00421                src.upperLeft(), src.lowerRight(), src.accessor(),
    00422                dest.upperLeft(), dest.lowerRight(), dest.accessor());
    00423
    00424     \endcode
    00425
    00426     <b> Required Interface:</b>
    00427
    00428     \code
    00429     SrcImageIterator src_upperleft, src_lowerright;
    00430     DestImageIterator dest_upperleft, src_lowerright;
    00431
    00432     SrcAccessor src_accessor;
    00433     DestAccessor dest_accessor;
    00434
    00435     NumericTraits<SrcAccessor::value_type>::RealPromote
    00436                              u = src_accessor(src_upperleft),
    00437                  v = src_accessor(src_upperleft, 1);
    00438     double d;
    00439
    00440     u = d * v;
    00441     u = u + v;
    00442
    00443     dest_accessor.set(
    00444         NumericTraits<DestAccessor::value_type>::fromRealPromote(u),
    00445     dest_upperleft);
    00446
    00447     \endcode
    00448
    00449     <b> Preconditions:</b>
    00450
    00451     \code
    00452     src_lowerright.x - src_upperleft.x > 1
    00453     src_lowerright.y - src_upperleft.y > 1
    00454     dest_lowerright.x - dest_upperleft.x > 1
    00455     dest_lowerright.y - dest_upperleft.y > 1
    00456     \endcode
    00457
    00458 */
    00459 doxygen_overloaded_function(template <...> void resizeImageLinearInterpolation)
    00460
    00461 template <class SrcIterator, class SrcAccessor,
    00462           class DestIterator, class DestAccessor>
    00463 void
    00464 resizeImageLinearInterpolation(SrcIterator is, SrcIterator iend, SrcAccessor sa,
    00465                       DestIterator id, DestIterator idend, DestAccessor da)
    00466 {
    00467     int w = iend.x - is.x;
    00468     int h = iend.y - is.y;
    00469
    00470     int wnew = idend.x - id.x;
    00471     int hnew = idend.y - id.y;
    00472
    00473     vigra_precondition((w > 1) && (h > 1),
    00474                  "resizeImageLinearInterpolation(): "
    00475                  "Source image too small.\n");
    00476     vigra_precondition((wnew > 1) && (hnew > 1),
    00477                  "resizeImageLinearInterpolation(): "
    00478                  "Destination image too small.\n");
    00479
    00480     double const scale = 2.0;
    00481
    00482     typedef typename SrcAccessor::value_type SRCVT;
    00483     typedef typename NumericTraits<SRCVT>::RealPromote TMPTYPE;
    00484     typedef BasicImage<TMPTYPE> TmpImage;
    00485     typedef typename TmpImage::traverser TmpImageIterator;
    00486
    00487     BasicImage<TMPTYPE> tmp(w, hnew);
    00488     BasicImage<TMPTYPE> line((h > w) ? h : w, 1);
    00489
    00490     int x,y;
    00491
    00492     typename BasicImage<TMPTYPE>::Iterator yt = tmp.upperLeft();
    00493     typename TmpImageIterator::row_iterator lt = line.upperLeft().rowIterator();
    00494
    00495     for(x=0; x<w; ++x, ++is.x, ++yt.x)
    00496     {
    00497         typename SrcIterator::column_iterator c1 = is.columnIterator();
    00498         typename TmpImageIterator::column_iterator ct = yt.columnIterator();
    00499
    00500         if(hnew < h)
    00501         {
    00502             recursiveSmoothLine(c1, c1 + h, sa,
    00503                  lt, line.accessor(), (double)h/hnew/scale);
    00504
    00505             resizeLineLinearInterpolation(lt, lt + h, line.accessor(),
    00506                                           ct, ct + hnew, tmp.accessor());
    00507         }
    00508         else
    00509         {
    00510             resizeLineLinearInterpolation(c1, c1 + h, sa,
    00511                                           ct, ct + hnew, tmp.accessor());
    00512         }
    00513     }
    00514
    00515     yt = tmp.upperLeft();
    00516
    00517     for(y=0; y < hnew; ++y, ++yt.y, ++id.y)
    00518     {
    00519         typename DestIterator::row_iterator rd = id.rowIterator();
    00520         typename TmpImageIterator::row_iterator rt = yt.rowIterator();
    00521
    00522         if(wnew < w)
    00523         {
    00524             recursiveSmoothLine(rt, rt + w, tmp.accessor(),
    00525                               lt, line.accessor(), (double)w/wnew/scale);
    00526
    00527             resizeLineLinearInterpolation(lt, lt + w, line.accessor(),
    00528                                           rd, rd + wnew, da);
    00529         }
    00530         else
    00531         {
    00532             resizeLineLinearInterpolation(rt, rt + w, tmp.accessor(),
    00533                                           rd, rd + wnew, da);
    00534         }
    00535     }
    00536 }
    00537
    00538 template <class SrcIterator, class SrcAccessor,
    00539           class DestIterator, class DestAccessor>
    00540 inline
    00541 void
    00542 resizeImageLinearInterpolation(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    00543                                triple<DestIterator, DestIterator, DestAccessor> dest)
    00544 {
    00545     resizeImageLinearInterpolation(src.first, src.second, src.third,
    00546                                    dest.first, dest.second, dest.third);
    00547 }
    00548
    00549 /***************************************************************/
    00550 /*                                                             */
    00551 /*                resizeImageSplineInterpolation               */
    00552 /*                                                             */
    00553 /***************************************************************/
    00554
    00555 /** \brief Resize image using B-spline interpolation.
    00556
    00557     The function implements separable spline interpolation algorithm described in
    00558
    00559     M. Unser, A. Aldroubi, M. Eden, <i>"B-Spline Signal Processing"</i>
    00560     IEEE Transactions on Signal Processing, vol. 41, no. 2, pp. 821-833 (part I),
    00561     pp. 834-848 (part II), 1993.
    00562
    00563     to obtain optimal interpolation quality and speed. You may pass the function
    00564     a spline of arbitrary order (e.g. <TT>BSpline<ORDER, double></tt> or
    00565     <TT>CatmullRomSpline<double></tt>). The default is a third order spline
    00566     which gives a twice continuously differentiable interpolant.
    00567     The implementation ensures that image values are interpolated rather
    00568     than smoothed by first calling a recursive (sharpening) prefilter as
    00569     described in the above paper. Then the actual interpolation is done
    00570     using \ref resamplingConvolveLine().
    00571
    00572     The range of both the input and output images (resp. regions)
    00573     must be given. The input image must have a size of at
    00574     least 4x4, the destination of at least 2x2. The scaling factors are then calculated
    00575     accordingly. If the source image is larger than the destination, it
    00576     is smoothed (band limited) using a recursive
    00577     exponential filter. The source value_type (SrcAccessor::value_type) must
    00578     be a linear algebra, i.e. it must support addition, subtraction,
    00579     and multiplication (+, -, *), multiplication with a scalar
    00580     real number and \ref NumericTraits "NumericTraits".
    00581     The function uses accessors.
    00582
    00583     <b> Declarations:</b>
    00584
    00585     pass arguments explicitly:
    00586     \code
    00587     namespace vigra {
    00588         template <class SrcImageIterator, class SrcAccessor,
    00589                   class DestImageIterator, class DestAccessor,
    00590                   class SPLINE>
    00591         void
    00592         resizeImageSplineInterpolation(
    00593               SrcImageIterator is, SrcImageIterator iend, SrcAccessor sa,
    00594               DestImageIterator id, DestImageIterator idend, DestAccessor da,
    00595               SPLINE spline = BSpline<3, double>())
    00596     }
    00597     \endcode
    00598
    00599
    00600     use argument objects in conjunction with \ref ArgumentObjectFactories :
    00601     \code
    00602     namespace vigra {
    00603         template <class SrcImageIterator, class SrcAccessor,
    00604                   class DestImageIterator, class DestAccessor,
    00605                   class SPLINE>
    00606         void
    00607         resizeImageSplineInterpolation(
    00608               triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
    00609               triple<DestImageIterator, DestImageIterator, DestAccessor> dest,
    00610               SPLINE spline = BSpline<3, double>())
    00611     }
    00612     \endcode
    00613
    00614     <b> Usage:</b>
    00615
    00616         <b>\#include</b> <vigra/resizeimage.hxx><br>
    00617         Namespace: vigra
    00618
    00619     \code
    00620     vigra::resizeImageSplineInterpolation(
    00621                src.upperLeft(), src.lowerRight(), src.accessor(),
    00622                dest.upperLeft(), dest.lowerRight(), dest.accessor());
    00623
    00624     \endcode
    00625
    00626     <b> Required Interface:</b>
    00627
    00628     \code
    00629     SrcImageIterator src_upperleft, src_lowerright;
    00630     DestImageIterator dest_upperleft, src_lowerright;
    00631
    00632     SrcAccessor src_accessor;
    00633     DestAccessor dest_accessor;
    00634
    00635     NumericTraits<SrcAccessor::value_type>::RealPromote
    00636                              u = src_accessor(src_upperleft),
    00637                  v = src_accessor(src_upperleft, 1);
    00638     double d;
    00639
    00640     u = d * v;
    00641     u = u + v;
    00642     u = u - v;
    00643     u = u * v;
    00644     u += v;
    00645     u -= v;
    00646
    00647     dest_accessor.set(
    00648         NumericTraits<DestAccessor::value_type>::fromRealPromote(u),
    00649     dest_upperleft);
    00650
    00651     \endcode
    00652
    00653     <b> Preconditions:</b>
    00654
    00655     \code
    00656     src_lowerright.x - src_upperleft.x > 3
    00657     src_lowerright.y - src_upperleft.y > 3
    00658     dest_lowerright.x - dest_upperleft.x > 1
    00659     dest_lowerright.y - dest_upperleft.y > 1
    00660     \endcode
    00661
    00662 */
    00663 doxygen_overloaded_function(template <...> void resizeImageSplineInterpolation)
    00664
    00665 template <class SrcIterator, class SrcAccessor,
    00666           class DestIterator, class DestAccessor,
    00667           class SPLINE>
    00668 void
    00669 resizeImageSplineInterpolation(
    00670     SrcIterator src_iter, SrcIterator src_iter_end, SrcAccessor src_acc,
    00671     DestIterator dest_iter, DestIterator dest_iter_end, DestAccessor dest_acc,
    00672     SPLINE const & spline)
    00673 {
    00674
    00675     int width_old = src_iter_end.x - src_iter.x;
    00676     int height_old = src_iter_end.y - src_iter.y;
    00677
    00678     int width_new = dest_iter_end.x - dest_iter.x;
    00679     int height_new = dest_iter_end.y - dest_iter.y;
    00680
    00681     vigra_precondition((width_old > 1) && (height_old > 1),
    00682                  "resizeImageSplineInterpolation(): "
    00683                  "Source image too small.\n");
    00684
    00685     vigra_precondition((width_new > 1) && (height_new > 1),
    00686                  "resizeImageSplineInterpolation(): "
    00687                  "Destination image too small.\n");
    00688
    00689     Rational<int> xratio(width_new - 1, width_old - 1);
    00690     Rational<int> yratio(height_new - 1, height_old - 1);
    00691     Rational<int> offset(0);
    00692     resampling_detail::MapTargetToSourceCoordinate xmapCoordinate(xratio, offset);
    00693     resampling_detail::MapTargetToSourceCoordinate ymapCoordinate(yratio, offset);
    00694     int xperiod = lcm(xratio.numerator(), xratio.denominator());
    00695     int yperiod = lcm(yratio.numerator(), yratio.denominator());
    00696
    00697     double const scale = 2.0;
    00698
    00699     typedef typename SrcAccessor::value_type SRCVT;
    00700     typedef typename NumericTraits<SRCVT>::RealPromote TMPTYPE;
    00701     typedef BasicImage<TMPTYPE> TmpImage;
    00702     typedef typename TmpImage::traverser TmpImageIterator;
    00703
    00704     BasicImage<TMPTYPE> tmp(width_old, height_new);
    00705
    00706     BasicImage<TMPTYPE> line((height_old > width_old) ? height_old : width_old, 1);
    00707     typename BasicImage<TMPTYPE>::Accessor tmp_acc = tmp.accessor();
    00708     ArrayVector<double> const & prefilterCoeffs = spline.prefilterCoefficients();
    00709
    00710     int x,y;
    00711
    00712     ArrayVector<Kernel1D<double> > kernels(yperiod);
    00713     createResamplingKernels(spline, ymapCoordinate, kernels);
    00714
    00715     typename BasicImage<TMPTYPE>::Iterator y_tmp = tmp.upperLeft();
    00716     typename TmpImageIterator::row_iterator line_tmp = line.upperLeft().rowIterator();
    00717
    00718     for(x=0; x<width_old; ++x, ++src_iter.x, ++y_tmp.x)
    00719     {
    00720
    00721         typename SrcIterator::column_iterator c_src = src_iter.columnIterator();
    00722         typename TmpImageIterator::column_iterator c_tmp = y_tmp.columnIterator();
    00723
    00724         if(prefilterCoeffs.size() == 0)
    00725         {
    00726             if(height_new >= height_old)
    00727             {
    00728                 resamplingConvolveLine(c_src, c_src + height_old, src_acc,
    00729                                        c_tmp, c_tmp + height_new, tmp_acc,
    00730                                        kernels, ymapCoordinate);
    00731             }
    00732             else
    00733             {
    00734                 recursiveSmoothLine(c_src, c_src + height_old, src_acc,
    00735                      line_tmp, line.accessor(), (double)height_old/height_new/scale);
    00736                 resamplingConvolveLine(line_tmp, line_tmp + height_old, line.accessor(),
    00737                                        c_tmp, c_tmp + height_new, tmp_acc,
    00738                                        kernels, ymapCoordinate);
    00739             }
    00740         }
    00741         else
    00742         {
    00743             recursiveFilterLine(c_src, c_src + height_old, src_acc,
    00744                                 line_tmp, line.accessor(),
    00745                                 prefilterCoeffs[0], BORDER_TREATMENT_REFLECT);
    00746             for(unsigned int b = 1; b < prefilterCoeffs.size(); ++b)
    00747             {
    00748                 recursiveFilterLine(line_tmp, line_tmp + height_old, line.accessor(),
    00749                                     line_tmp, line.accessor(),
    00750                                     prefilterCoeffs[b], BORDER_TREATMENT_REFLECT);
    00751             }
    00752             if(height_new < height_old)
    00753             {
    00754                 recursiveSmoothLine(line_tmp, line_tmp + height_old, line.accessor(),
    00755                      line_tmp, line.accessor(), (double)height_old/height_new/scale);
    00756             }
    00757             resamplingConvolveLine(line_tmp, line_tmp + height_old, line.accessor(),
    00758                                    c_tmp, c_tmp + height_new, tmp_acc,
    00759                                    kernels, ymapCoordinate);
    00760         }
    00761     }
    00762
    00763     y_tmp = tmp.upperLeft();
    00764
    00765     kernels.resize(xperiod);
    00766     createResamplingKernels(spline, xmapCoordinate, kernels);
    00767
    00768     for(y=0; y < height_new; ++y, ++y_tmp.y, ++dest_iter.y)
    00769     {
    00770         typename DestIterator::row_iterator r_dest = dest_iter.rowIterator();
    00771         typename TmpImageIterator::row_iterator r_tmp = y_tmp.rowIterator();
    00772
    00773         if(prefilterCoeffs.size() == 0)
    00774         {
    00775             if(width_new >= width_old)
    00776             {
    00777                 resamplingConvolveLine(r_tmp, r_tmp + width_old, tmp.accessor(),
    00778                                        r_dest, r_dest + width_new, dest_acc,
    00779                                        kernels, xmapCoordinate);
    00780             }
    00781             else
    00782             {
    00783                 recursiveSmoothLine(r_tmp, r_tmp + width_old, tmp.accessor(),
    00784                                   line_tmp, line.accessor(), (double)width_old/width_new/scale);
    00785                 resamplingConvolveLine(line_tmp, line_tmp + width_old, line.accessor(),
    00786                                        r_dest, r_dest + width_new, dest_acc,
    00787                                        kernels, xmapCoordinate);
    00788             }
    00789         }
    00790         else
    00791         {
    00792             recursiveFilterLine(r_tmp, r_tmp + width_old, tmp.accessor(),
    00793                                 line_tmp, line.accessor(),
    00794                                 prefilterCoeffs[0], BORDER_TREATMENT_REFLECT);
    00795             for(unsigned int b = 1; b < prefilterCoeffs.size(); ++b)
    00796             {
    00797                 recursiveFilterLine(line_tmp, line_tmp + width_old, line.accessor(),
    00798                                     line_tmp, line.accessor(),
    00799                                     prefilterCoeffs[b], BORDER_TREATMENT_REFLECT);
    00800             }
    00801             if(width_new < width_old)
    00802             {
    00803                 recursiveSmoothLine(line_tmp, line_tmp + width_old, line.accessor(),
    00804                                     line_tmp, line.accessor(), (double)width_old/width_new/scale);
    00805             }
    00806             resamplingConvolveLine(line_tmp, line_tmp + width_old, line.accessor(),
    00807                                    r_dest, r_dest + width_new, dest_acc,
    00808                                    kernels, xmapCoordinate);
    00809         }
    00810     }
    00811 }
    00812
    00813 template <class SrcIterator, class SrcAccessor,
    00814           class DestIterator, class DestAccessor,
    00815           class SPLINE>
    00816 inline
    00817 void
    00818 resizeImageSplineInterpolation(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    00819                       triple<DestIterator, DestIterator, DestAccessor> dest,
    00820                       SPLINE const & spline)
    00821 {
    00822     resizeImageSplineInterpolation(src.first, src.second, src.third,
    00823                                    dest.first, dest.second, dest.third, spline);
    00824 }
    00825
    00826 template <class SrcIterator, class SrcAccessor,
    00827           class DestIterator, class DestAccessor>
    00828 void
    00829 resizeImageSplineInterpolation(SrcIterator is, SrcIterator iend, SrcAccessor sa,
    00830                       DestIterator id, DestIterator idend, DestAccessor da)
    00831 {
    00832     resizeImageSplineInterpolation(is, iend, sa, id, idend, da, BSpline<3, double>());
    00833 }
    00834
    00835 template <class SrcIterator, class SrcAccessor,
    00836           class DestIterator, class DestAccessor>
    00837 inline
    00838 void
    00839 resizeImageSplineInterpolation(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    00840                       triple<DestIterator, DestIterator, DestAccessor> dest)
    00841 {
    00842     resizeImageSplineInterpolation(src.first, src.second, src.third,
    00843                                    dest.first, dest.second, dest.third);
    00844 }
    00845
    00846 /*****************************************************************/
    00847 /*                                                               */
    00848 /*              resizeImageCatmullRomInterpolation               */
    00849 /*                                                               */
    00850 /*****************************************************************/
    00851
    00852 /** \brief Resize image using the Catmull/Rom interpolation function.
    00853
    00854     The function calls like \ref resizeImageSplineInterpolation() with
    00855     \ref vigra::CatmullRomSpline as an interpolation kernel.
    00856     The interpolated function has one continuous derivative.
    00857     (See \ref resizeImageSplineInterpolation() for more documentation)
    00858
    00859     <b> Declarations:</b>
    00860
    00861     pass arguments explicitly:
    00862     \code
    00863     namespace vigra {
    00864         template <class SrcIterator, class SrcAccessor,
    00865                   class DestIterator, class DestAccessor>
    00866         void
    00867         resizeImageCatmullRomInterpolation(SrcIterator src_iter, SrcIterator src_iter_end, SrcAccessor src_acc,
    00868                               DestIterator dest_iter, DestIterator dest_iter_end, DestAccessor dest_acc);
    00869     }
    00870     \endcode
    00871
    00872
    00873     use argument objects in conjunction with \ref ArgumentObjectFactories :
    00874     \code
    00875     namespace vigra {
    00876         template <class SrcIterator, class SrcAccessor,
    00877                   class DestIterator, class DestAccessor>
    00878         void
    00879         resizeImageCatmullRomInterpolation(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    00880                               triple<DestIterator, DestIterator, DestAccessor> dest);
    00881     }
    00882     \endcode
    00883
    00884
    00885     <b>\#include</b> <vigra/resizeimage.hxx><br>
    00886     Namespace: vigra
    00887
    00888 */
    00889 doxygen_overloaded_function(template <...> void resizeImageCatmullRomInterpolation)
    00890
    00891 template <class SrcIterator, class SrcAccessor,
    00892           class DestIterator, class DestAccessor>
    00893 inline void
    00894 resizeImageCatmullRomInterpolation(SrcIterator src_iter, SrcIterator src_iter_end, SrcAccessor src_acc,
    00895                       DestIterator dest_iter, DestIterator dest_iter_end, DestAccessor dest_acc)
    00896 {
    00897     resizeImageSplineInterpolation(src_iter, src_iter_end, src_acc, dest_iter, dest_iter_end, dest_acc,
    00898                                   CatmullRomSpline<double>());
    00899 }
    00900
    00901 template <class SrcIterator, class SrcAccessor,
    00902           class DestIterator, class DestAccessor>
    00903 inline
    00904 void
    00905 resizeImageCatmullRomInterpolation(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    00906                       triple<DestIterator, DestIterator, DestAccessor> dest)
    00907 {
    00908     resizeImageCatmullRomInterpolation(src.first, src.second, src.third,
    00909                                      dest.first, dest.second, dest.third);
    00910 }
    00911
    00912 #if 0
    00913 /*****************************************************************/
    00914 /*                                                               */
    00915 /*              resizeImageCubicInterpolation                    */
    00916 /*                                                               */
    00917 /*****************************************************************/
    00918
    00919 /** \brief Resize image using the cardinal B-spline interpolation function.
    00920
    00921     The function calls like \ref resizeImageSplineInterpolation() with
    00922     \ref vigra::BSpline<3, double> and prefiltering as an interpolation kernel.
    00923     The interpolated function has two continuous derivatives.
    00924     (See \ref resizeImageSplineInterpolation() for more documentation)
    00925
    00926     <b>\#include</b> <vigra/resizeimage.hxx><br>
    00927     Namespace: vigra
    00928
    00929 */
    00930 template <class SrcIterator, class SrcAccessor,
    00931           class DestIterator, class DestAccessor>
    00932 void
    00933 resizeImageCubicInterpolation(SrcIterator src_iter, SrcIterator src_iter_end, SrcAccessor src_acc,
    00934                       DestIterator dest_iter, DestIterator dest_iter_end, DestAccessor dest_acc)
    00935 {
    00936     resizeImageSplineInterpolation(src_iter, src_iter_end, src_acc, dest_iter, dest_iter_end, dest_acc,
    00937                                   BSpline<3, double>());
    00938 }
    00939
    00940 template <class SrcIterator, class SrcAccessor,
    00941           class DestIterator, class DestAccessor>
    00942 inline
    00943 void
    00944 resizeImageCubicInterpolation(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    00945                       triple<DestIterator, DestIterator, DestAccessor> dest)
    00946 {
    00947     resizeImageCubicInterpolation(src.first, src.second, src.third,
    00948                                    dest.first, dest.second, dest.third);
    00949 }
    00950 #endif
    00951
    00952 /*****************************************************************/
    00953 /*                                                               */
    00954 /*              resizeImageCoscotInterpolation                   */
    00955 /*                                                               */
    00956 /*****************************************************************/
    00957
    00958 /** \brief Resize image using the Coscot interpolation function.
    00959
    00960     The function calls \ref resizeImageSplineInterpolation() with
    00961     \ref vigra::CoscotFunction as an interpolation kernel.
    00962     The interpolated function has one continuous derivative.
    00963     (See \ref resizeImageSplineInterpolation() for more documentation)
    00964
    00965     <b> Declarations:</b>
    00966
    00967     pass arguments explicitly:
    00968     \code
    00969     namespace vigra {
    00970         template <class SrcIterator, class SrcAccessor,
    00971                   class DestIterator, class DestAccessor>
    00972         void
    00973         resizeImageCoscotInterpolation(SrcIterator src_iter, SrcIterator src_iter_end, SrcAccessor src_acc,
    00974                               DestIterator dest_iter, DestIterator dest_iter_end, DestAccessor dest_acc);
    00975     }
    00976     \endcode
    00977
    00978
    00979     use argument objects in conjunction with \ref ArgumentObjectFactories :
    00980     \code
    00981     namespace vigra {
    00982         template <class SrcIterator, class SrcAccessor,
    00983                   class DestIterator, class DestAccessor>
    00984         void
    00985         resizeImageCoscotInterpolation(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    00986                               triple<DestIterator, DestIterator, DestAccessor> dest);
    00987     }
    00988     \endcode
    00989
    00990
    00991     <b>\#include</b> <vigra/resizeimage.hxx><br>
    00992     Namespace: vigra
    00993
    00994 */
    00995 doxygen_overloaded_function(template <...> void resizeImageCoscotInterpolation)
    00996
    00997 template <class SrcIterator, class SrcAccessor,
    00998           class DestIterator, class DestAccessor>
    00999 void
    01000 resizeImageCoscotInterpolation(SrcIterator src_iter, SrcIterator src_iter_end, SrcAccessor src_acc,
    01001                       DestIterator dest_iter, DestIterator dest_iter_end, DestAccessor dest_acc)
    01002 {
    01003     resizeImageSplineInterpolation(src_iter, src_iter_end, src_acc, dest_iter, dest_iter_end, dest_acc,
    01004                                    CoscotFunction<double>());
    01005 }
    01006
    01007 template <class SrcIterator, class SrcAccessor,
    01008           class DestIterator, class DestAccessor>
    01009 inline
    01010 void
    01011 resizeImageCoscotInterpolation(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    01012                       triple<DestIterator, DestIterator, DestAccessor> dest)
    01013 {
    01014     resizeImageCoscotInterpolation(src.first, src.second, src.third,
    01015                                    dest.first, dest.second, dest.third);
    01016 }
    01017
    01018
    01019 #if 0 // old version of the spline interpolation algorithm
    01020
    01021 /********************************************************/
    01022 /*                                                      */
    01023 /*           resizeCalculateSplineCoefficients          */
    01024 /*         (internally used by resize functions)        */
    01025 /*                                                      */
    01026 /********************************************************/
    01027
    01028 template <class SrcIterator, class SrcAccessor, class VALUETYPE>
    01029 void
    01030 resizeCalculateSplineCoefficients(SrcIterator i1, SrcIterator iend,
    01031                 SrcAccessor a, VALUETYPE * i2)
    01032 {
    01033     int n = iend - i1;
    01034
    01035     if(n <= 0) return;
    01036
    01037     VALUETYPE zero = NumericTraits<VALUETYPE>::zero();
    01038     VALUETYPE two = 2.0 * NumericTraits<VALUETYPE>::one();
    01039     VALUETYPE half = 0.5 * NumericTraits<VALUETYPE>::one();
    01040
    01041     *i2 = zero;
    01042     if(n == 1) return;
    01043
    01044     std::vector<VALUETYPE> vec(n);
    01045     typename std::vector<VALUETYPE>::iterator u = vec.begin();
    01046
    01047     *u = zero;
    01048
    01049     for(++i1, ++i2, ++u, --iend; i1 != iend; ++i1, ++i2, ++u)
    01050     {
    01051         VALUETYPE p = 0.5 * i2[-1] + two;
    01052         *i2 = half / p;
    01053         *u = 3.0 *(a(i1,1) - 2.0 * a(i1) + a(i1, -1)) - 0.5 * u[-1] / p;
    01054     }
    01055
    01056     *i2 = zero;
    01057
    01058     for(--i2, --u; u != vec; --u, --i2)
    01059     {
    01060         *i2 = *i2 * i2[1] + *u;
    01061     }
    01062 }
    01063
    01064 /********************************************************/
    01065 /*                                                      */
    01066 /*         resizeImageInternalSplineGradient            */
    01067 /*                                                      */
    01068 /********************************************************/
    01069
    01070 template <class SrcIterator, class SrcAccessor,
    01071           class DoubleIterator, class TempIterator, class DestIterator>
    01072 void
    01073 resizeImageInternalSplineGradient(SrcIterator in, SrcIterator inend, SrcAccessor sa,
    01074                          DoubleIterator tmp, TempIterator r, DestIterator id)
    01075 {
    01076     int w = inend - in;
    01077
    01078     int x;
    01079
    01080     typedef typename SrcAccessor::value_type SRCVT;
    01081     typedef typename NumericTraits<SRCVT>::RealPromote TMPTYPE;
    01082
    01083     // calculate border derivatives
    01084     SrcIterator xs = in;
    01085     TMPTYPE p0 = -11.0/6.0 * sa(xs);  ++xs;
    01086             p0 += 3.0 * sa(xs);  ++xs;
    01087             p0 += -1.5 * sa(xs);  ++xs;
    01088             p0 += 1.0/3.0 * sa(xs);
    01089
    01090     xs = in + w-1;
    01091     TMPTYPE pw = 11.0/6.0 * sa(xs);  --xs;
    01092             pw += -3.0 * sa(xs);  --xs;
    01093             pw +=  1.5 * sa(xs);  --xs;
    01094             pw += -1.0/3.0 * sa(xs);
    01095
    01096     xs = in + 2;
    01097     SrcIterator xs1 = in;
    01098
    01099     for(x=1; x<w-1; ++x, ++xs, ++xs1)
    01100     {
    01101         r[x] = 3.0 * (sa(xs) - sa(xs1));
    01102     }
    01103
    01104     r[1] -= p0;
    01105     r[w-2] -= pw;
    01106
    01107     double q = 0.25;
    01108
    01109     id[0] = p0;
    01110     id[w-1] = pw;
    01111     id[1] = 0.25 * r[1];
    01112
    01113     for(x=2; x<w-1; ++x)
    01114     {
    01115         tmp[x] = q;
    01116         q = 1.0 / (4.0 - q);
    01117         id[x] = q * (r[x] - id[x-1]);
    01118     }
    01119
    01120     for(x=w-3; x>=1; --x)
    01121     {
    01122         id[x] -= tmp[x+1]*id[x+1];
    01123     }
    01124 }
    01125
    01126 /********************************************************/
    01127 /*                                                      */
    01128 /*         resizeImageInternalSplineInterpolation       */
    01129 /*                                                      */
    01130 /********************************************************/
    01131
    01132 template <class SrcIterator, class SrcAccessor,
    01133           class DestIterator, class DestAccessor>
    01134 void
    01135 resizeImageInternalSplineInterpolation(SrcIterator is, SrcIterator iend, SrcAccessor sa,
    01136                       DestIterator id, DestIterator idend, DestAccessor da)
    01137 {
    01138     int w = iend.x - is.x;
    01139     int h = iend.y - is.y;
    01140
    01141     int wnew = idend.x - id.x;
    01142     int hnew = idend.y - id.y;
    01143
    01144     typedef typename SrcAccessor::value_type SRCVT;
    01145     typedef typename NumericTraits<SRCVT>::RealPromote TMPTYPE;
    01146     typedef typename BasicImage<TMPTYPE>::Iterator TMPITER;
    01147     typedef
    01148         NumericTraits<typename DestAccessor::value_type> DestTraits;
    01149
    01150     BasicImage<TMPTYPE> dx(w,h);
    01151     BasicImage<TMPTYPE> dy(w,h);
    01152     BasicImage<TMPTYPE> dxy(w,h);
    01153     BasicImage<TMPTYPE> W(4,4), W1(4,4);
    01154     std::vector<TMPTYPE> R(w > h ? w : h);
    01155     std::vector<double> tmp(w > h ? w : h);
    01156
    01157     typename BasicImage<TMPTYPE>::Accessor ta;
    01158
    01159     SrcIterator in = is;
    01160
    01161     TMPITER idx = dx.upperLeft();
    01162     TMPITER idy = dy.upperLeft();
    01163     TMPITER idxy = dxy.upperLeft();
    01164     typename std::vector<TMPTYPE>::iterator r = R.begin();
    01165     typename std::vector<double>::iterator it = tmp.begin();
    01166
    01167     double ig[] = { 1.0, 0.0, -3.0,  2.0,
    01168                     0.0, 1.0, -2.0,  1.0,
    01169                     0.0, 0.0,  3.0, -2.0,
    01170                     0.0, 0.0, -1.0,  1.0 };
    01171
    01172     int x, y, i, j, k;
    01173
    01174
    01175     // calculate x derivatives
    01176     for(y=0; y<h; ++y, ++in.y, ++idx.y)
    01177     {
    01178         typename SrcIterator::row_iterator sr = in.rowIterator();
    01179         typename TMPITER::row_iterator dr = idx.rowIterator();
    01180         resizeImageInternalSplineGradient(sr, sr+w, sa,
    01181                                           it, r, dr);
    01182     }
    01183
    01184     in = is;
    01185
    01186     // calculate y derivatives
    01187     for(x=0; x<w; ++x, ++in.x, ++idy.x)
    01188     {
    01189         typename SrcIterator::column_iterator sc = in.columnIterator();
    01190         typename TMPITER::column_iterator dc = idy.columnIterator();
    01191         resizeImageInternalSplineGradient(sc, sc+h, sa,
    01192                                           it, r, dc);
    01193     }
    01194
    01195     in = is;
    01196     idy = dy.upperLeft();
    01197
    01198     // calculate mixed derivatives
    01199     for(y=0; y<h; ++y, ++idy.y, ++idxy.y)
    01200     {
    01201         typename TMPITER::row_iterator sr = idy.rowIterator();
    01202         typename TMPITER::row_iterator dr = idxy.rowIterator();
    01203         resizeImageInternalSplineGradient(sr, sr+w, ta,
    01204                                           it, r, dr);
    01205     }
    01206
    01207     double du = (double)(w-1) / (wnew-1);
    01208     double dv = (double)(h-1) / (hnew-1);
    01209     double ov = 0.0;
    01210     int oy = 0;
    01211     int yy = oy;
    01212
    01213     DestIterator xxd = id, yyd = id;
    01214
    01215     static Diff2D down(0,1), right(1,0), downright(1,1);
    01216
    01217     for(y=0; y<h-1; ++y, ++in.y, ov -= 1.0)
    01218     {
    01219         if(y < h-2 && ov >= 1.0) continue;
    01220         int y1 = y+1;
    01221         double v = ov;
    01222         double ou = 0.0;
    01223         int ox = 0;
    01224         int xx = ox;
    01225
    01226         SrcIterator xs = in;
    01227         for(x=0; x<w-1; ++x, ++xs.x, ou -= 1.0)
    01228         {
    01229             if(x < w-2 && ou >= 1.0) continue;
    01230             int x1 = x+1;
    01231             double u = ou;
    01232
    01233             DestIterator xd = id + Diff2D(ox,oy);
    01234             W[0][0] = sa(xs);
    01235             W[0][1] = dy(x, y);
    01236             W[0][2] = sa(xs, down);
    01237             W[0][3] = dy(x, y1);
    01238             W[1][0] = dx(x, y);
    01239             W[1][1] = dxy(x, y);
    01240             W[1][2] = dx(x, y1);
    01241             W[1][3] = dxy(x, y1);
    01242             W[2][0] = sa(xs, right);
    01243             W[2][1] = dy(x1,y);
    01244             W[2][2] = sa(xs, downright);
    01245             W[2][3] = dy(x1, y1);
    01246             W[3][0] = dx(x1, y);
    01247             W[3][1] = dxy(x1, y);
    01248             W[3][2] = dx(x1, y1);
    01249             W[3][3] = dxy(x1, y1);
    01250
    01251             for(i=0; i<4; ++i)
    01252             {
    01253                 for(j=0; j<4; ++j)
    01254                 {
    01255                     W1[j][i] = ig[j] * W[0][i];
    01256                     for(k=1; k<4; ++k)
    01257                     {
    01258                         W1[j][i] += ig[j+4*k] * W[k][i];
    01259                     }
    01260                 }
    01261             }
    01262             for(i=0; i<4; ++i)
    01263             {
    01264                 for(j=0; j<4; ++j)
    01265                 {
    01266                     W[j][i] = ig[i] * W1[j][0];
    01267                     for(k=1; k<4; ++k)
    01268                     {
    01269                        W[j][i] += ig[4*k+i] * W1[j][k];
    01270                     }
    01271                 }
    01272             }
    01273
    01274             TMPTYPE a1,a2,a3,a4;
    01275
    01276             yyd = xd;
    01277             for(v=ov, yy=oy; v<1.0; v+=dv, ++yyd.y, ++yy)
    01278             {
    01279                 a1 = W[0][0] + v * (W[0][1] +
    01280                                v * (W[0][2] + v * W[0][3]));
    01281                 a2 = W[1][0] + v * (W[1][1] +
    01282                                v * (W[1][2] + v * W[1][3]));
    01283                 a3 = W[2][0] + v * (W[2][1] +
    01284                                v * (W[2][2] + v * W[2][3]));
    01285                 a4 = W[3][0] + v * (W[3][1] +
    01286                                v * (W[3][2] + v * W[3][3]));
    01287
    01288                 xxd = yyd;
    01289                 for(u=ou, xx=ox; u<1.0; u+=du, ++xxd.x, ++xx)
    01290                 {
    01291                     da.set(DestTraits::fromRealPromote(a1 + u * (a2 + u * (a3 + u * a4))), xxd);
    01292                 }
    01293
    01294                 if(xx == wnew-1)
    01295                 {
    01296                     da.set(DestTraits::fromRealPromote(a1 + a2 + a3 + a4), xxd);
    01297                 }
    01298             }
    01299
    01300             if(yy == hnew-1)
    01301             {
    01302                 a1 = W[0][0] + W[0][1] + W[0][2] + W[0][3];
    01303                 a2 = W[1][0] + W[1][1] + W[1][2] + W[1][3];
    01304                 a3 = W[2][0] + W[2][1] + W[2][2] + W[2][3];
    01305                 a4 = W[3][0] + W[3][1] + W[3][2] + W[3][3];
    01306
    01307                 DestIterator xxd = yyd;
    01308                 for(u=ou, xx=ox; u<1.0; u+=du, ++xxd.x, ++xx)
    01309                 {
    01310                     da.set(DestTraits::fromRealPromote(a1 + u * (a2 + u * (a3 + u * a4))), xxd);
    01311                 }
    01312
    01313                 if(xx == wnew-1)
    01314                 {
    01315                     da.set(DestTraits::fromRealPromote(a1 + a2 + a3 + a4), xxd);
    01316                 }
    01317             }
    01318
    01319             ou = u;
    01320             ox = xx;
    01321         }
    01322         ov = v;
    01323         oy = yy;
    01324     }
    01325 }
    01326
    01327 template <class SrcIterator, class SrcAccessor,
    01328           class DestIterator, class DestAccessor>
    01329 void
    01330 resizeImageSplineInterpolation(SrcIterator is, SrcIterator iend, SrcAccessor sa,
    01331                       DestIterator id, DestIterator idend, DestAccessor da)
    01332 {
    01333     int w = iend.x - is.x;
    01334     int h = iend.y - is.y;
    01335
    01336     int wnew = idend.x - id.x;
    01337     int hnew = idend.y - id.y;
    01338
    01339     vigra_precondition((w > 3) && (h > 3),
    01340                  "resizeImageSplineInterpolation(): "
    01341                  "Source image too small.\n");
    01342     vigra_precondition((wnew > 1) && (hnew > 1),
    01343                  "resizeImageSplineInterpolation(): "
    01344                  "Destination image too small.\n");
    01345
    01346     double scale = 2.0;
    01347
    01348     if(wnew < w || hnew < h)
    01349     {
    01350         typedef typename SrcAccessor::value_type SRCVT;
    01351         typedef typename NumericTraits<SRCVT>::RealPromote TMPTYPE;
    01352         typedef typename BasicImage<TMPTYPE>::Iterator TMPITER;
    01353
    01354         BasicImage<TMPTYPE> t(w,h);
    01355         TMPITER it = t.upperLeft();
    01356
    01357         if(wnew < w)
    01358         {
    01359             recursiveSmoothX(is, iend, sa,
    01360                     it, t.accessor(), (double)w/wnew/scale);
    01361
    01362             if(hnew < h)
    01363             {
    01364                recursiveSmoothY(it, t.lowerRight(), t.accessor(),
    01365                     it, t.accessor(), (double)h/hnew/scale);
    01366             }
    01367         }
    01368         else
    01369         {
    01370            recursiveSmoothY(is, iend, sa,
    01371                     it, t.accessor(), (double)h/hnew/scale);
    01372         }
    01373
    01374         resizeImageInternalSplineInterpolation(it, t.lowerRight(), t.accessor(),
    01375                                                id, idend, da);
    01376     }
    01377     else
    01378     {
    01379         resizeImageInternalSplineInterpolation(is, iend, sa, id, idend, da);
    01380     }
    01381 }
    01382
    01383 template <class SrcIterator, class SrcAccessor,
    01384           class DestIterator, class DestAccessor>
    01385 inline
    01386 void
    01387 resizeImageSplineInterpolation(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    01388                       triple<DestIterator, DestIterator, DestAccessor> dest)
    01389 {
    01390     resizeImageSplineInterpolation(src.first, src.second, src.third,
    01391                                    dest.first, dest.second, dest.third);
    01392 }
    01393 #endif  // old algorithm version
    01394
    01395 //@}
    01396
    01397 } // namespace vigra
    01398
    01399 #endif // VIGRA_RESIZEIMAGE_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)