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

vigra/resizeimage.hxx

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