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

details vigra/resampling_convolution.hxx VIGRA

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.3.2, Jan 27 2005 )                                    */
00008 /*    You may use, modify, and distribute this software according       */
00009 /*    to the terms stated in the LICENSE file included in               */
00010 /*    the VIGRA distribution.                                           */
00011 /*                                                                      */
00012 /*    The VIGRA Website is                                              */
00013 /*        http://kogs-www.informatik.uni-hamburg.de/~koethe/vigra/      */
00014 /*    Please direct questions, bug reports, and contributions to        */
00015 /*        koethe@informatik.uni-hamburg.de                              */
00016 /*                                                                      */
00017 /*  THIS SOFTWARE IS PROVIDED AS IS AND WITHOUT ANY EXPRESS OR          */
00018 /*  IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED      */
00019 /*  WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. */
00020 /*                                                                      */
00021 /************************************************************************/
00022 
00023 #ifndef VIGRA_RESAMPLING_CONVOLUTION_HXX
00024 #define VIGRA_RESAMPLING_CONVOLUTION_HXX
00025 
00026 #include <cmath>
00027 #include "vigra/stdimage.hxx"
00028 #include "vigra/array_vector.hxx"
00029 #include "vigra/rational.hxx"
00030 #include "vigra/functortraits.hxx"
00031 
00032 namespace vigra {
00033 
00034 namespace resampling_detail
00035 {
00036 
00037 struct MapTargetToSourceCoordinate
00038 {
00039     MapTargetToSourceCoordinate(Rational<int> const & samplingRatio, 
00040                                 Rational<int> const & offset)
00041     : a(samplingRatio.denominator()*offset.denominator()),
00042       b(samplingRatio.numerator()*offset.numerator()),
00043       c(samplingRatio.numerator()*offset.denominator())
00044     {}
00045     
00046 //        the following funcions are more efficient realizations of:
00047 //             rational_cast<T>(i / samplingRatio + offset);
00048 //        we need efficiency because this may be called in the inner loop
00049 
00050     int operator()(int i) const
00051     {
00052         return (i * a + b) / c;
00053     }
00054     
00055     double toDouble(int i) const
00056     {
00057         return double(i * a + b) / c;
00058     }
00059 
00060     Rational<int> toRational(int i) const
00061     {
00062         return Rational<int>(i * a + b, c);
00063     }
00064 
00065     int a, b, c;
00066 };
00067 
00068 } // namespace resampling_detail
00069 
00070 template <>
00071 class FunctorTraits<resampling_detail::MapTargetToSourceCoordinate>
00072 : public FunctorTraitsBase<resampling_detail::MapTargetToSourceCoordinate>
00073 {
00074   public:
00075     typedef VigraTrueType isUnaryFunctor;
00076 };
00077 
00078 template <class SrcIter, class SrcAcc,
00079           class DestIter, class DestAcc,
00080           class KernelArray,
00081           class Functor>
00082 void 
00083 resamplingConvolveLine(SrcIter s, SrcIter send, SrcAcc src,
00084                        DestIter d, DestIter dend, DestAcc dest,
00085                        KernelArray const & kernels,
00086                        Functor mapTargetToSourceCoordinate)
00087 {
00088     typedef typename 
00089         NumericTraits<typename SrcAcc::value_type>::RealPromote
00090         TmpType;
00091     typedef typename KernelArray::value_type Kernel;
00092     typedef typename Kernel::const_iterator KernelIter;
00093     
00094     int wo = send - s;
00095     int wn = dend - d;
00096     int wo2 = 2*wo - 2;
00097     
00098     int i;
00099     typename KernelArray::const_iterator kernel = kernels.begin();
00100     for(i=0; i<wn; ++i, ++d, ++kernel)
00101     {
00102         // use the kernels periodically
00103         if(kernel == kernels.end())
00104             kernel = kernels.begin();
00105         
00106         // calculate current target point into source location
00107         int is = mapTargetToSourceCoordinate(i);
00108         
00109         TmpType sum = NumericTraits<TmpType>::zero();
00110 
00111         int lbound = is - kernel->right(), 
00112             hbound = is - kernel->left();
00113                     
00114         KernelIter k = kernel->center() + kernel->right();
00115         if(lbound < 0 || hbound >= wo)
00116         {    
00117             vigra_precondition(-lbound < wo && wo2 - hbound >= 0,
00118                 "resamplingConvolveLine(): kernel or offset larger than image.");
00119             for(int m=lbound; m <= hbound; ++m, --k)
00120             {
00121                 int mm = (m < 0) ?
00122                             -m :
00123                             (m >= wo) ?
00124                                 wo2 - m :
00125                                 m;
00126                 sum += *k * src(s, mm);
00127             }
00128         }
00129         else
00130         {
00131             SrcIter ss = s + lbound;
00132             SrcIter ssend = s + hbound;
00133             
00134             for(; ss <= ssend; ++ss, --k)
00135             {
00136                 sum += *k * src(ss);
00137             }
00138         }
00139         
00140         dest.set(sum, d);
00141     }
00142 }
00143 
00144 template <class Kernel, class MapCoordinate, class KernelArray>
00145 void
00146 createResamplingKernels(Kernel const & kernel, 
00147              MapCoordinate const & mapCoordinate, KernelArray & kernels)
00148 {
00149     for(unsigned int idest = 0; idest < kernels.size(); ++idest)
00150     {
00151         int isrc = mapCoordinate(idest);
00152         double idsrc = mapCoordinate.toDouble(idest);
00153         double offset = idsrc - isrc;
00154         double radius = kernel.radius();
00155         int left = int(ceil(-radius - offset));
00156         int right = int(floor(radius - offset));
00157         kernels[idest].initExplicitly(left, right);
00158         
00159         double x = left + offset;
00160         for(int i = left; i <= right; ++i, ++x)
00161             kernels[idest][i] = kernel(x);
00162         kernels[idest].normalize(1.0, kernel.derivativeOrder(), offset);
00163     }
00164 }
00165 
00166 /** \addtogroup ResamplingConvolutionFilters Resampling Convolution Filters
00167 
00168     These functions implement the convolution operation when the source and target images
00169     have different sizes. This is realized by accessing a continous kernel at the
00170     appropriate non-integer positions. The technique is, for example, described in
00171     D. Schumacher: <i>General Filtered Image Rescaling</i>, in: Graphics Gems III, 
00172     Academic Press, 1992.
00173 */
00174 //@{
00175 
00176 /********************************************************/
00177 /*                                                      */
00178 /*                  resamplingConvolveX                 */
00179 /*                                                      */
00180 /********************************************************/
00181 
00182 /** \brief Apply a resampling filter in the x-direction.
00183 
00184     This function implements a convolution operation in x-direction
00185     (i.e. applies a 1D filter to every row) where the width of the source
00186     and destination images differ. This is typically used to avoid aliasing if 
00187     the image is scaled down, or to interpolate smoothly if the image is scaled up.
00188     The target coordinates are transformed into source coordinates by
00189     
00190     \code
00191     xsource = (xtarget - offset) / samplingRatio
00192     \endcode
00193     
00194     The <tt>samplingRatio</tt> and <tt>offset</tt> must be given as \ref vigra::Rational 
00195     in order to avoid rounding errors in this transformation. It is required that for all
00196     pixels of the target image, <tt>xsource</tt> remains within the range of the source
00197     image (i.e. <tt>0 <= xsource <= sourceWidth-1</tt>. Since <tt>xsource</tt> is
00198     in general not an integer, the <tt>kernel</tt> must be a functor that can be accessed at 
00199     arbitrary (<tt>double</tt>) coordinates. It must also provide a member function <tt>radius()</tt>
00200     which specifies the support (non-zero interval) of the kernel. VIGRA already
00201     provides a number of suitable functors, e.g. \ref vigra::Gaussian, \ref vigra::BSpline
00202     \ref vigra::CatmullRomSpline, and \ref vigra::CoscotFunction. The function
00203     \ref resizeImageSplineInterpolation() is implemented by means resamplingConvolveX() and 
00204     resamplingConvolveY().
00205 
00206     <b> Declarations:</b>
00207 
00208     pass arguments explicitly:
00209     \code
00210     namespace vigra {
00211         template <class SrcIter, class SrcAcc,
00212                   class DestIter, class DestAcc,
00213                   class Kernel>
00214         void 
00215         resamplingConvolveX(SrcIter sul, SrcIter slr, SrcAcc src,
00216                             DestIter dul, DestIter dlr, DestAcc dest,
00217                             Kernel const & kernel,
00218                             Rational<int> const & samplingRatio, Rational<int> const & offset);
00219     }
00220     \endcode
00221 
00222 
00223     use argument objects in conjunction with \ref ArgumentObjectFactories:
00224     \code
00225     namespace vigra {
00226         template <class SrcIter, class SrcAcc,
00227                   class DestIter, class DestAcc,
00228                   class Kernel>
00229         void 
00230         resamplingConvolveX(triple<SrcIter, SrcIter, SrcAcc> src,
00231                             triple<DestIter, DestIter, DestAcc> dest,
00232                             Kernel const & kernel,
00233                             Rational<int> const & samplingRatio, Rational<int> const & offset);
00234     }
00235     \endcode
00236 
00237     <b> Usage:</b>
00238 
00239     <b>\#include</b> "<a href="resampling_convolution_8hxx-source.html">vigra/resampling_convolution.hxx</a>"
00240 
00241 
00242     \code
00243     Rational<int> ratio(2), offset(0);
00244 
00245     FImage src(w,h), 
00246            dest(rational_cast<int>(ratio*w), h);
00247 
00248     float sigma = 2.0;
00249     Gaussian<float> smooth(sigma);
00250     ...
00251 
00252     // simpultaneously enlarge and smooth source image
00253     resamplingConvolveX(srcImageRange(src), destImageRange(dest), 
00254                         smooth, ratio, offset);
00255     \endcode
00256 
00257     <b> Required Interface:</b>
00258     
00259     \code
00260     Kernel kernel;
00261     int kernelRadius = kernel.radius();
00262     double x = ...;  // must be <= radius()
00263     double value = kernel(x);
00264     \endcode
00265 */
00266 template <class SrcIter, class SrcAcc,
00267           class DestIter, class DestAcc,
00268           class Kernel>
00269 void 
00270 resamplingConvolveX(SrcIter sul, SrcIter slr, SrcAcc src,
00271                     DestIter dul, DestIter dlr, DestAcc dest,
00272                     Kernel const & kernel,
00273                     Rational<int> const & samplingRatio, Rational<int> const & offset)
00274 {
00275     int wold = slr.x - sul.x;
00276     int wnew = dlr.x - dul.x;
00277     
00278     vigra_precondition(!samplingRatio.is_inf() && samplingRatio > 0,
00279                 "resamplingConvolveX(): sampling ratio must be > 0 and < infinity");
00280     vigra_precondition(!offset.is_inf(),
00281                 "resamplingConvolveX(): offset must be < infinity");
00282 
00283     int period = lcm(samplingRatio.numerator(), samplingRatio.denominator());
00284     resampling_detail::MapTargetToSourceCoordinate mapCoordinate(samplingRatio, offset);
00285     
00286     ArrayVector<Kernel1D<double> > kernels(period);
00287     
00288     createResamplingKernels(kernel, mapCoordinate, kernels);
00289 
00290     for(; sul.y < slr.y; ++sul.y, ++dul.y)
00291     {
00292         typename SrcIter::row_iterator sr = sul.rowIterator();
00293         typename DestIter::row_iterator dr = dul.rowIterator();
00294         resamplingConvolveLine(sr, sr+wold, src, dr, dr+wnew, dest,
00295                                kernels, mapCoordinate);
00296     }
00297 }
00298 
00299 template <class SrcIter, class SrcAcc,
00300           class DestIter, class DestAcc,
00301           class Kernel>
00302 inline void 
00303 resamplingConvolveX(triple<SrcIter, SrcIter, SrcAcc> src,
00304                     triple<DestIter, DestIter, DestAcc> dest,
00305                     Kernel const & kernel,
00306                     Rational<int> const & samplingRatio, Rational<int> const & offset)
00307 {
00308     resamplingConvolveX(src.first, src.second, src.third,
00309                         dest.first, dest.second, dest.third,
00310                         kernel, samplingRatio, offset);
00311 }
00312 
00313 /********************************************************/
00314 /*                                                      */
00315 /*                  resamplingConvolveY                 */
00316 /*                                                      */
00317 /********************************************************/
00318 
00319 /** \brief Apply a resampling filter in the y-direction.
00320 
00321     This function implements a convolution operation in y-direction
00322     (i.e. applies a 1D filter to every column) where the height of the source
00323     and destination images differ. This is typically used to avoid aliasing if 
00324     the image is scaled down, or to interpolate smoothly if the image is scaled up.
00325     The target coordinates are transformed into source coordinates by
00326     
00327     \code
00328     ysource = (ytarget - offset) / samplingRatio
00329     \endcode
00330     
00331     The <tt>samplingRatio</tt> and <tt>offset</tt> must be given as \ref vigra::Rational 
00332     in order to avoid rounding errors in this transformation. It is required that for all
00333     pixels of the target image, <tt>ysource</tt> remains within the range of the source
00334     image (i.e. <tt>0 <= ysource <= sourceHeight-1</tt>. Since <tt>ysource</tt> is
00335     in general not an integer, the <tt>kernel</tt> must be a functor that can be accessed at 
00336     arbitrary (<tt>double</tt>) coordinates. It must also provide a member function <tt>radius()</tt>
00337     which specifies the support (non-zero interval) of the kernel. VIGRA already
00338     provides a number of suitable functors, e.g. \ref vigra::Gaussian, \ref vigra::BSpline
00339     \ref vigra::CatmullRomSpline, and \ref vigra::CoscotFunction. The function
00340     \ref resizeImageSplineInterpolation() is implemented by means resamplingConvolveX() and 
00341     resamplingConvolveY().
00342 
00343     <b> Declarations:</b>
00344 
00345     pass arguments explicitly:
00346     \code
00347     namespace vigra {
00348         template <class SrcIter, class SrcAcc,
00349                   class DestIter, class DestAcc,
00350                   class Kernel>
00351         void 
00352         resamplingConvolveY(SrcIter sul, SrcIter slr, SrcAcc src,
00353                             DestIter dul, DestIter dlr, DestAcc dest,
00354                             Kernel const & kernel,
00355                             Rational<int> const & samplingRatio, Rational<int> const & offset);
00356     }
00357     \endcode
00358 
00359 
00360     use argument objects in conjunction with \ref ArgumentObjectFactories:
00361     \code
00362     namespace vigra {
00363         template <class SrcIter, class SrcAcc,
00364                   class DestIter, class DestAcc,
00365                   class Kernel>
00366         void 
00367         resamplingConvolveY(triple<SrcIter, SrcIter, SrcAcc> src,
00368                             triple<DestIter, DestIter, DestAcc> dest,
00369                             Kernel const & kernel,
00370                             Rational<int> const & samplingRatio, Rational<int> const & offset);
00371     }
00372     \endcode
00373 
00374     <b> Usage:</b>
00375 
00376     <b>\#include</b> "<a href="resampling_convolution_8hxx-source.html">vigra/resampling_convolution.hxx</a>"
00377 
00378 
00379     \code
00380     Rational<int> ratio(2), offset(0);
00381 
00382     FImage src(w,h), 
00383            dest(w, rational_cast<int>(ratio*h));
00384 
00385     float sigma = 2.0;
00386     Gaussian<float> smooth(sigma);
00387     ...
00388 
00389     // simpultaneously enlarge and smooth source image
00390     resamplingConvolveY(srcImageRange(src), destImageRange(dest), 
00391                         smooth, ratio, offset);
00392     \endcode
00393 
00394     <b> Required Interface:</b>
00395     
00396     \code
00397     Kernel kernel;
00398     int kernelRadius = kernel.radius();
00399     double y = ...;  // must be <= radius()
00400     double value = kernel(y);
00401     \endcode
00402 */
00403 template <class SrcIter, class SrcAcc,
00404           class DestIter, class DestAcc,
00405           class Kernel>
00406 void 
00407 resamplingConvolveY(SrcIter sul, SrcIter slr, SrcAcc src,
00408                     DestIter dul, DestIter dlr, DestAcc dest,
00409                     Kernel const & kernel,
00410                     Rational<int> const & samplingRatio, Rational<int> const & offset)
00411 {
00412     int hold = slr.y - sul.y;
00413     int hnew = dlr.y - dul.y;
00414     
00415     vigra_precondition(!samplingRatio.is_inf() && samplingRatio > 0,
00416                 "resamplingConvolveY(): sampling ratio must be > 0 and < infinity");
00417     vigra_precondition(!offset.is_inf(),
00418                 "resamplingConvolveY(): offset must be < infinity");
00419 
00420     int period = lcm(samplingRatio.numerator(), samplingRatio.denominator());
00421     
00422     resampling_detail::MapTargetToSourceCoordinate mapCoordinate(samplingRatio, offset);
00423     
00424     ArrayVector<Kernel1D<double> > kernels(period);
00425     
00426     createResamplingKernels(kernel, mapCoordinate, kernels);
00427 
00428     for(; sul.x < slr.x; ++sul.x, ++dul.x)
00429     {
00430         typename SrcIter::column_iterator sc = sul.columnIterator();
00431         typename DestIter::column_iterator dc = dul.columnIterator();
00432         resamplingConvolveLine(sc, sc+hold, src, dc, dc+hnew, dest,
00433                                kernels, mapCoordinate);
00434     }
00435 }
00436 
00437 template <class SrcIter, class SrcAcc,
00438           class DestIter, class DestAcc,
00439           class Kernel>
00440 inline void 
00441 resamplingConvolveY(triple<SrcIter, SrcIter, SrcAcc> src,
00442                     triple<DestIter, DestIter, DestAcc> dest,
00443                     Kernel const & kernel,
00444                     Rational<int> const & samplingRatio, Rational<int> const & offset)
00445 {
00446     resamplingConvolveY(src.first, src.second, src.third,
00447                         dest.first, dest.second, dest.third,
00448                         kernel, samplingRatio, offset);
00449 }
00450 
00451 /********************************************************/
00452 /*                                                      */
00453 /*               resamplingConvolveImage                */
00454 /*                                                      */
00455 /********************************************************/
00456 
00457 /** \brief Apply two separable resampling filters successively, the first in x-direction, 
00458            the second in y-direction.
00459 
00460     This function is a shorthand for the concatenation of a call to
00461     \link ResamplingConvolutionFilters#resamplingConvolveX resamplingConvolveX\endlink()
00462     and \link ResamplingConvolutionFilters#resamplingConvolveY resamplingConvolveY\endlink() 
00463     with the given kernels. See there for detailed documentation.
00464 
00465     <b> Declarations:</b>
00466 
00467     pass arguments explicitly:
00468     \code
00469     namespace vigra {
00470         template <class SrcIterator, class SrcAccessor,
00471                   class DestIterator, class DestAccessor,
00472                   class KernelX, class KernelY>
00473         void resamplingConvolveImage(SrcIterator sul,SrcIterator slr, SrcAccessor src,
00474                            DestIterator dul, DestIterator dlr, DestAccessor dest,
00475                            KernelX const & kx, 
00476                            Rational<int> const & samplingRatioX, Rational<int> const & offsetX,
00477                            KernelY const & ky, 
00478                            Rational<int> const & samplingRatioY, Rational<int> const & offsetY);
00479     }
00480     \endcode
00481 
00482 
00483     use argument objects in conjunction with \ref ArgumentObjectFactories:
00484     \code
00485     namespace vigra {
00486         template <class SrcIterator, class SrcAccessor,
00487                   class DestIterator, class DestAccessor,
00488                   class KernelX, class KernelY>
00489         void
00490         resamplingConvolveImage(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00491                            triple<DestIterator, DestIterator, DestAccessor> dest,
00492                            KernelX const & kx, 
00493                            Rational<int> const & samplingRatioX, Rational<int> const & offsetX,
00494                            KernelY const & ky, 
00495                            Rational<int> const & samplingRatioY, Rational<int> const & offsetY);
00496     }
00497     \endcode
00498 
00499     <b> Usage:</b>
00500 
00501     <b>\#include</b> "<a href="resampling_convolution_8hxx-source.html">vigra/resampling_convolution.hxx</a>"
00502 
00503 
00504     \code
00505     Rational<int> xratio(2), yratio(3), offset(0);
00506 
00507     FImage src(w,h), 
00508            dest(rational_cast<int>(xratio*w), rational_cast<int>(yratio*h));
00509 
00510     float sigma = 2.0;
00511     Gaussian<float> smooth(sigma);
00512     ...
00513 
00514     // simpultaneously enlarge and smooth source image
00515     resamplingConvolveImage(srcImageRange(src), destImageRange(dest), 
00516                             smooth, xratio, offset,
00517                             smooth, yratio, offset);
00518 
00519     \endcode
00520 
00521 */
00522 template <class SrcIterator, class SrcAccessor,
00523           class DestIterator, class DestAccessor,
00524           class KernelX, class KernelY>
00525 void resamplingConvolveImage(SrcIterator sul,SrcIterator slr, SrcAccessor src,
00526                    DestIterator dul, DestIterator dlr, DestAccessor dest,
00527                    KernelX const & kx, 
00528                    Rational<int> const & samplingRatioX, Rational<int> const & offsetX,
00529                    KernelY const & ky, 
00530                    Rational<int> const & samplingRatioY, Rational<int> const & offsetY)
00531 {
00532     typedef typename
00533         NumericTraits<typename SrcAccessor::value_type>::RealPromote
00534         TmpType;
00535     
00536     BasicImage<TmpType> tmp(dlr.x - dul.x, slr.y - sul.y);
00537 
00538     resamplingConvolveX(srcIterRange(sul, slr, src),
00539                         destImageRange(tmp), 
00540                         kx, samplingRatioX, offsetX);
00541     resamplingConvolveY(srcImageRange(tmp),
00542                         destIterRange(dul, dlr, dest), 
00543                         ky, samplingRatioY, offsetY);
00544 }
00545 
00546 template <class SrcIterator, class SrcAccessor,
00547           class DestIterator, class DestAccessor,
00548           class KernelX, class KernelY>
00549 inline void
00550 resamplingConvolveImage(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00551                    triple<DestIterator, DestIterator, DestAccessor> dest,
00552                    KernelX const & kx, 
00553                    Rational<int> const & samplingRatioX, Rational<int> const & offsetX,
00554                    KernelY const & ky, 
00555                    Rational<int> const & samplingRatioY, Rational<int> const & offsetY)
00556 {
00557     resamplingConvolveImage(src.first, src.second, src.third,
00558                             dest.first, dest.second, dest.third,
00559                             kx, samplingRatioX, offsetX,
00560                             ky, samplingRatioY, offsetY);
00561 }
00562 
00563 } // namespace vigra 
00564 
00565 
00566 #endif /* VIGRA_RESAMPLING_CONVOLUTION_HXX */

© Ullrich Köthe (koethe@informatik.uni-hamburg.de)
Cognitive Systems Group, University of Hamburg, Germany

html generated using doxygen and Python
VIGRA 1.3.2 (27 Jan 2005)