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

vigra/resampling_convolution.hxx VIGRA

00001 /************************************************************************/
00002 /*                                                                      */
00003 /*               Copyright 1998-2004 by Ullrich Koethe                  */
00004 /*                                                                      */
00005 /*    This file is part of the VIGRA computer vision library.           */
00006 /*    The VIGRA Website is                                              */
00007 /*        http://hci.iwr.uni-heidelberg.de/vigra/                       */
00008 /*    Please direct questions, bug reports, and contributions to        */
00009 /*        ullrich.koethe@iwr.uni-heidelberg.de    or                    */
00010 /*        vigra@informatik.uni-hamburg.de                               */
00011 /*                                                                      */
00012 /*    Permission is hereby granted, free of charge, to any person       */
00013 /*    obtaining a copy of this software and associated documentation    */
00014 /*    files (the "Software"), to deal in the Software without           */
00015 /*    restriction, including without limitation the rights to use,      */
00016 /*    copy, modify, merge, publish, distribute, sublicense, and/or      */
00017 /*    sell copies of the Software, and to permit persons to whom the    */
00018 /*    Software is furnished to do so, subject to the following          */
00019 /*    conditions:                                                       */
00020 /*                                                                      */
00021 /*    The above copyright notice and this permission notice shall be    */
00022 /*    included in all copies or substantial portions of the             */
00023 /*    Software.                                                         */
00024 /*                                                                      */
00025 /*    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND    */
00026 /*    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   */
00027 /*    OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          */
00028 /*    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT       */
00029 /*    HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,      */
00030 /*    WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING      */
00031 /*    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR     */
00032 /*    OTHER DEALINGS IN THE SOFTWARE.                                   */
00033 /*                                                                      */
00034 /************************************************************************/
00035 
00036 #ifndef VIGRA_RESAMPLING_CONVOLUTION_HXX
00037 #define VIGRA_RESAMPLING_CONVOLUTION_HXX
00038 
00039 #include <cmath>
00040 #include "stdimage.hxx"
00041 #include "array_vector.hxx"
00042 #include "rational.hxx"
00043 #include "functortraits.hxx"
00044 #include "functorexpression.hxx"
00045 #include "transformimage.hxx"
00046 #include "imagecontainer.hxx"
00047 
00048 namespace vigra {
00049 
00050 namespace resampling_detail
00051 {
00052 
00053 struct MapTargetToSourceCoordinate
00054 {
00055     MapTargetToSourceCoordinate(Rational<int> const & samplingRatio,
00056                                 Rational<int> const & offset)
00057     : a(samplingRatio.denominator()*offset.denominator()),
00058       b(samplingRatio.numerator()*offset.numerator()),
00059       c(samplingRatio.numerator()*offset.denominator())
00060     {}
00061 
00062 //        the following functions are more efficient realizations of:
00063 //             rational_cast<T>(i / samplingRatio + offset);
00064 //        we need efficiency because this may be called in the inner loop
00065 
00066     int operator()(int i) const
00067     {
00068         return (i * a + b) / c;
00069     }
00070 
00071     double toDouble(int i) const
00072     {
00073         return double(i * a + b) / c;
00074     }
00075 
00076     Rational<int> toRational(int i) const
00077     {
00078         return Rational<int>(i * a + b, c);
00079     }
00080     
00081     bool isExpand2() const
00082     {
00083         return a == 1 && b == 0 && c == 2;
00084     }
00085     
00086     bool isReduce2() const
00087     {
00088         return a == 2 && b == 0 && c == 1;
00089     }
00090 
00091     int a, b, c;
00092 };
00093 
00094 } // namespace resampling_detail
00095 
00096 template <>
00097 class FunctorTraits<resampling_detail::MapTargetToSourceCoordinate>
00098 : public FunctorTraitsBase<resampling_detail::MapTargetToSourceCoordinate>
00099 {
00100   public:
00101     typedef VigraTrueType isUnaryFunctor;
00102 };
00103 
00104 template <class SrcIter, class SrcAcc,
00105           class DestIter, class DestAcc,
00106           class KernelArray>
00107 void
00108 resamplingExpandLine2(SrcIter s, SrcIter send, SrcAcc src,
00109                        DestIter d, DestIter dend, DestAcc dest,
00110                        KernelArray const & kernels)
00111 {
00112     typedef typename KernelArray::value_type Kernel;
00113     typedef typename KernelArray::const_reference KernelRef;
00114     typedef typename Kernel::const_iterator KernelIter;
00115 
00116     typedef typename
00117         PromoteTraits<typename SrcAcc::value_type, typename Kernel::value_type>::Promote
00118         TmpType;
00119 
00120     int wo = send - s;
00121     int wn = dend - d;
00122     int wo2 = 2*wo - 2;
00123     
00124     int ileft = std::max(kernels[0].right(), kernels[1].right());
00125     int iright = wo + std::min(kernels[0].left(), kernels[1].left()) - 1;
00126     for(int i = 0; i < wn; ++i, ++d)
00127     {
00128         int is = i / 2;
00129         KernelRef kernel = kernels[i & 1];
00130         KernelIter k = kernel.center() + kernel.right();
00131         TmpType sum = NumericTraits<TmpType>::zero();        
00132         if(is < ileft)
00133         {
00134             for(int m=is-kernel.right(); m <= is-kernel.left(); ++m, --k)
00135             {
00136                 int mm = (m < 0) 
00137                         ? -m 
00138                         : m;
00139                 sum += *k * src(s, mm);
00140             }        
00141         }
00142         else if(is > iright)
00143         {
00144             for(int m=is-kernel.right(); m <= is-kernel.left(); ++m, --k)
00145             {
00146                 int mm =  (m >= wo) 
00147                             ? wo2 - m
00148                             : m;
00149                 sum += *k * src(s, mm);
00150             }        
00151         }
00152         else
00153         {
00154             SrcIter ss = s + is - kernel.right();
00155             for(int m = 0; m < kernel.size(); ++m, --k, ++ss)
00156             {
00157                 sum += *k * src(ss);
00158             }        
00159         }
00160         dest.set(sum, d);
00161     }
00162 }
00163 
00164 template <class SrcIter, class SrcAcc,
00165           class DestIter, class DestAcc,
00166           class KernelArray>
00167 void
00168 resamplingReduceLine2(SrcIter s, SrcIter send, SrcAcc src,
00169                        DestIter d, DestIter dend, DestAcc dest,
00170                        KernelArray const & kernels)
00171 {
00172     typedef typename KernelArray::value_type Kernel;
00173     typedef typename KernelArray::const_reference KernelRef;
00174     typedef typename Kernel::const_iterator KernelIter;
00175 
00176     KernelRef kernel = kernels[0];
00177     KernelIter kbegin = kernel.center() + kernel.right();
00178 
00179     typedef typename
00180         PromoteTraits<typename SrcAcc::value_type, typename Kernel::value_type>::Promote
00181         TmpType;
00182 
00183     int wo = send - s;
00184     int wn = dend - d;
00185     int wo2 = 2*wo - 2;
00186     
00187     int ileft = kernel.right();
00188     int iright = wo + kernel.left() - 1;
00189     for(int i = 0; i < wn; ++i, ++d)
00190     {
00191         int is = 2 * i;
00192         KernelIter k = kbegin;
00193         TmpType sum = NumericTraits<TmpType>::zero();        
00194         if(is < ileft)
00195         {
00196             for(int m=is-kernel.right(); m <= is-kernel.left(); ++m, --k)
00197             {
00198                 int mm = (m < 0) 
00199                         ? -m 
00200                         : m;
00201                 sum += *k * src(s, mm);
00202             }        
00203         }
00204         else if(is > iright)
00205         {
00206             for(int m=is-kernel.right(); m <= is-kernel.left(); ++m, --k)
00207             {
00208                 int mm =  (m >= wo) 
00209                             ? wo2 - m
00210                             : m;
00211                 sum += *k * src(s, mm);
00212             }        
00213         }
00214         else
00215         {
00216             SrcIter ss = s + is - kernel.right();
00217             for(int m = 0; m < kernel.size(); ++m, --k, ++ss)
00218             {
00219                 sum += *k * src(ss);
00220             }        
00221         }
00222         dest.set(sum, d);
00223     }
00224 }
00225 
00226 /** \addtogroup ResamplingConvolutionFilters Resampling Convolution Filters
00227 
00228     These functions implement the convolution operation when the source and target images
00229     have different sizes. This is realized by accessing a continuous kernel at the
00230     appropriate non-integer positions. The technique is, for example, described in
00231     D. Schumacher: <i>General Filtered Image Rescaling</i>, in: Graphics Gems III,
00232     Academic Press, 1992.
00233 */
00234 //@{
00235 
00236 /********************************************************/
00237 /*                                                      */
00238 /*                resamplingConvolveLine                */
00239 /*                                                      */
00240 /********************************************************/
00241 
00242 /** \brief Performs a 1-dimensional resampling convolution of the source signal using the given
00243     set of kernels.
00244 
00245     This function is mainly used internally: It is called for each dimension of a 
00246     higher dimensional array in order to perform a separable resize operation.
00247 
00248     <b> Declaration:</b>
00249 
00250     <b>\#include</b> <vigra/resampling_convolution.hxx>
00251 
00252     \code
00253     namespace vigra {
00254         template <class SrcIter, class SrcAcc,
00255                   class DestIter, class DestAcc,
00256                   class KernelArray,
00257                   class Functor>
00258         void
00259         resamplingConvolveLine(SrcIter s, SrcIter send, SrcAcc src,
00260                                DestIter d, DestIter dend, DestAcc dest,
00261                                KernelArray const & kernels,
00262                                Functor mapTargetToSourceCoordinate)    
00263     }
00264     \endcode
00265 
00266 */
00267 doxygen_overloaded_function(template <...> void resamplingConvolveLine)
00268 
00269 template <class SrcIter, class SrcAcc,
00270           class DestIter, class DestAcc,
00271           class KernelArray,
00272           class Functor>
00273 void
00274 resamplingConvolveLine(SrcIter s, SrcIter send, SrcAcc src,
00275                        DestIter d, DestIter dend, DestAcc dest,
00276                        KernelArray const & kernels,
00277                        Functor mapTargetToSourceCoordinate)
00278 {
00279     if(mapTargetToSourceCoordinate.isExpand2())
00280     {
00281         resamplingExpandLine2(s, send, src, d, dend, dest, kernels);
00282         return;
00283     }
00284     if(mapTargetToSourceCoordinate.isReduce2())
00285     {
00286         resamplingReduceLine2(s, send, src, d, dend, dest, kernels);
00287         return;
00288     }
00289     
00290     typedef typename
00291         NumericTraits<typename SrcAcc::value_type>::RealPromote
00292         TmpType;
00293     typedef typename KernelArray::value_type Kernel;
00294     typedef typename Kernel::const_iterator KernelIter;
00295 
00296     int wo = send - s;
00297     int wn = dend - d;
00298     int wo2 = 2*wo - 2;
00299 
00300     int i;
00301     typename KernelArray::const_iterator kernel = kernels.begin();
00302     for(i=0; i<wn; ++i, ++d, ++kernel)
00303     {
00304         // use the kernels periodically
00305         if(kernel == kernels.end())
00306             kernel = kernels.begin();
00307 
00308         // calculate current target point into source location
00309         int is = mapTargetToSourceCoordinate(i);
00310 
00311         TmpType sum = NumericTraits<TmpType>::zero();
00312 
00313         int lbound = is - kernel->right(),
00314             hbound = is - kernel->left();
00315 
00316         KernelIter k = kernel->center() + kernel->right();
00317         if(lbound < 0 || hbound >= wo)
00318         {
00319             vigra_precondition(-lbound < wo && wo2 - hbound >= 0,
00320                 "resamplingConvolveLine(): kernel or offset larger than image.");
00321             for(int m=lbound; m <= hbound; ++m, --k)
00322             {
00323                 int mm = (m < 0) ?
00324                             -m :
00325                             (m >= wo) ?
00326                                 wo2 - m :
00327                                 m;
00328                 sum = TmpType(sum + *k * src(s, mm));
00329             }
00330         }
00331         else
00332         {
00333             SrcIter ss = s + lbound;
00334             SrcIter ssend = s + hbound;
00335 
00336             for(; ss <= ssend; ++ss, --k)
00337             {
00338                 sum = TmpType(sum + *k * src(ss));
00339             }
00340         }
00341 
00342         dest.set(sum, d);
00343     }
00344 }
00345 
00346 template <class Kernel, class MapCoordinate, class KernelArray>
00347 void
00348 createResamplingKernels(Kernel const & kernel,
00349              MapCoordinate const & mapCoordinate, KernelArray & kernels)
00350 {
00351     for(unsigned int idest = 0; idest < kernels.size(); ++idest)
00352     {
00353         int isrc = mapCoordinate(idest);
00354         double idsrc = mapCoordinate.toDouble(idest);
00355         double offset = idsrc - isrc;
00356         double radius = kernel.radius();
00357         int left = int(ceil(-radius - offset));
00358         int right = int(floor(radius - offset));
00359         kernels[idest].initExplicitly(left, right);
00360 
00361         double x = left + offset;
00362         for(int i = left; i <= right; ++i, ++x)
00363             kernels[idest][i] = kernel(x);
00364         kernels[idest].normalize(1.0, kernel.derivativeOrder(), offset);
00365     }
00366 }
00367 
00368 /** \brief Apply a resampling filter in the x-direction.
00369 
00370     This function implements a convolution operation in x-direction
00371     (i.e. applies a 1D filter to every row) where the width of the source
00372     and destination images differ. This is typically used to avoid aliasing if
00373     the image is scaled down, or to interpolate smoothly if the image is scaled up.
00374     The target coordinates are transformed into source coordinates by
00375 
00376     \code
00377     xsource = (xtarget - offset) / samplingRatio
00378     \endcode
00379 
00380     The <tt>samplingRatio</tt> and <tt>offset</tt> must be given as \ref vigra::Rational
00381     in order to avoid rounding errors in this transformation. It is required that for all
00382     pixels of the target image, <tt>xsource</tt> remains within the range of the source
00383     image (i.e. <tt>0 <= xsource <= sourceWidth-1</tt>. Since <tt>xsource</tt> is
00384     in general not an integer, the <tt>kernel</tt> must be a functor that can be accessed at
00385     arbitrary (<tt>double</tt>) coordinates. It must also provide a member function <tt>radius()</tt>
00386     which specifies the support (non-zero interval) of the kernel. VIGRA already
00387     provides a number of suitable functors, e.g. \ref vigra::Gaussian, \ref vigra::BSpline
00388     \ref vigra::CatmullRomSpline, and \ref vigra::CoscotFunction. The function
00389     \ref resizeImageSplineInterpolation() is implemented by means resamplingConvolveX() and
00390     resamplingConvolveY().
00391 
00392     <b> Declarations:</b>
00393 
00394     pass arguments explicitly:
00395     \code
00396     namespace vigra {
00397         template <class SrcIter, class SrcAcc,
00398                   class DestIter, class DestAcc,
00399                   class Kernel>
00400         void
00401         resamplingConvolveX(SrcIter sul, SrcIter slr, SrcAcc src,
00402                             DestIter dul, DestIter dlr, DestAcc dest,
00403                             Kernel const & kernel,
00404                             Rational<int> const & samplingRatio, Rational<int> const & offset);
00405     }
00406     \endcode
00407 
00408 
00409     use argument objects in conjunction with \ref ArgumentObjectFactories :
00410     \code
00411     namespace vigra {
00412         template <class SrcIter, class SrcAcc,
00413                   class DestIter, class DestAcc,
00414                   class Kernel>
00415         void
00416         resamplingConvolveX(triple<SrcIter, SrcIter, SrcAcc> src,
00417                             triple<DestIter, DestIter, DestAcc> dest,
00418                             Kernel const & kernel,
00419                             Rational<int> const & samplingRatio, Rational<int> const & offset);
00420     }
00421     \endcode
00422 
00423     <b> Usage:</b>
00424 
00425     <b>\#include</b> <vigra/resampling_convolution.hxx>
00426 
00427 
00428     \code
00429     Rational<int> ratio(2), offset(0);
00430 
00431     FImage src(w,h),
00432            dest(rational_cast<int>(ratio*w), h);
00433 
00434     float sigma = 2.0;
00435     Gaussian<float> smooth(sigma);
00436     ...
00437 
00438     // simultaneously enlarge and smooth source image
00439     resamplingConvolveX(srcImageRange(src), destImageRange(dest),
00440                         smooth, ratio, offset);
00441     \endcode
00442 
00443     <b> Required Interface:</b>
00444 
00445     \code
00446     Kernel kernel;
00447     int kernelRadius = kernel.radius();
00448     double x = ...;  // must be <= radius()
00449     double value = kernel(x);
00450     \endcode
00451 */
00452 doxygen_overloaded_function(template <...> void resamplingConvolveX)
00453 
00454 template <class SrcIter, class SrcAcc,
00455           class DestIter, class DestAcc,
00456           class Kernel>
00457 void
00458 resamplingConvolveX(SrcIter sul, SrcIter slr, SrcAcc src,
00459                     DestIter dul, DestIter dlr, DestAcc dest,
00460                     Kernel const & kernel,
00461                     Rational<int> const & samplingRatio, Rational<int> const & offset)
00462 {
00463     int wold = slr.x - sul.x;
00464     int wnew = dlr.x - dul.x;
00465 
00466     vigra_precondition(!samplingRatio.is_inf() && samplingRatio > 0,
00467                 "resamplingConvolveX(): sampling ratio must be > 0 and < infinity");
00468     vigra_precondition(!offset.is_inf(),
00469                 "resamplingConvolveX(): offset must be < infinity");
00470 
00471     int period = lcm(samplingRatio.numerator(), samplingRatio.denominator());
00472     resampling_detail::MapTargetToSourceCoordinate mapCoordinate(samplingRatio, offset);
00473 
00474     ArrayVector<Kernel1D<double> > kernels(period);
00475 
00476     createResamplingKernels(kernel, mapCoordinate, kernels);
00477 
00478     for(; sul.y < slr.y; ++sul.y, ++dul.y)
00479     {
00480         typename SrcIter::row_iterator sr = sul.rowIterator();
00481         typename DestIter::row_iterator dr = dul.rowIterator();
00482         resamplingConvolveLine(sr, sr+wold, src, dr, dr+wnew, dest,
00483                                kernels, mapCoordinate);
00484     }
00485 }
00486 
00487 template <class SrcIter, class SrcAcc,
00488           class DestIter, class DestAcc,
00489           class Kernel>
00490 inline void
00491 resamplingConvolveX(triple<SrcIter, SrcIter, SrcAcc> src,
00492                     triple<DestIter, DestIter, DestAcc> dest,
00493                     Kernel const & kernel,
00494                     Rational<int> const & samplingRatio, Rational<int> const & offset)
00495 {
00496     resamplingConvolveX(src.first, src.second, src.third,
00497                         dest.first, dest.second, dest.third,
00498                         kernel, samplingRatio, offset);
00499 }
00500 
00501 /********************************************************/
00502 /*                                                      */
00503 /*                  resamplingConvolveY                 */
00504 /*                                                      */
00505 /********************************************************/
00506 
00507 /** \brief Apply a resampling filter in the y-direction.
00508 
00509     This function implements a convolution operation in y-direction
00510     (i.e. applies a 1D filter to every column) where the height of the source
00511     and destination images differ. This is typically used to avoid aliasing if
00512     the image is scaled down, or to interpolate smoothly if the image is scaled up.
00513     The target coordinates are transformed into source coordinates by
00514 
00515     \code
00516     ysource = (ytarget - offset) / samplingRatio
00517     \endcode
00518 
00519     The <tt>samplingRatio</tt> and <tt>offset</tt> must be given as \ref vigra::Rational
00520     in order to avoid rounding errors in this transformation. It is required that for all
00521     pixels of the target image, <tt>ysource</tt> remains within the range of the source
00522     image (i.e. <tt>0 <= ysource <= sourceHeight-1</tt>. Since <tt>ysource</tt> is
00523     in general not an integer, the <tt>kernel</tt> must be a functor that can be accessed at
00524     arbitrary (<tt>double</tt>) coordinates. It must also provide a member function <tt>radius()</tt>
00525     which specifies the support (non-zero interval) of the kernel. VIGRA already
00526     provides a number of suitable functors, e.g. \ref vigra::Gaussian, \ref vigra::BSpline
00527     \ref vigra::CatmullRomSpline, and \ref vigra::CoscotFunction. The function
00528     \ref resizeImageSplineInterpolation() is implemented by means resamplingConvolveX() and
00529     resamplingConvolveY().
00530 
00531     <b> Declarations:</b>
00532 
00533     pass arguments explicitly:
00534     \code
00535     namespace vigra {
00536         template <class SrcIter, class SrcAcc,
00537                   class DestIter, class DestAcc,
00538                   class Kernel>
00539         void
00540         resamplingConvolveY(SrcIter sul, SrcIter slr, SrcAcc src,
00541                             DestIter dul, DestIter dlr, DestAcc dest,
00542                             Kernel const & kernel,
00543                             Rational<int> const & samplingRatio, Rational<int> const & offset);
00544     }
00545     \endcode
00546 
00547 
00548     use argument objects in conjunction with \ref ArgumentObjectFactories :
00549     \code
00550     namespace vigra {
00551         template <class SrcIter, class SrcAcc,
00552                   class DestIter, class DestAcc,
00553                   class Kernel>
00554         void
00555         resamplingConvolveY(triple<SrcIter, SrcIter, SrcAcc> src,
00556                             triple<DestIter, DestIter, DestAcc> dest,
00557                             Kernel const & kernel,
00558                             Rational<int> const & samplingRatio, Rational<int> const & offset);
00559     }
00560     \endcode
00561 
00562     <b> Usage:</b>
00563 
00564     <b>\#include</b> <vigra/resampling_convolution.hxx>
00565 
00566 
00567     \code
00568     Rational<int> ratio(2), offset(0);
00569 
00570     FImage src(w,h),
00571            dest(w, rational_cast<int>(ratio*h));
00572 
00573     float sigma = 2.0;
00574     Gaussian<float> smooth(sigma);
00575     ...
00576 
00577     // simultaneously enlarge and smooth source image
00578     resamplingConvolveY(srcImageRange(src), destImageRange(dest),
00579                         smooth, ratio, offset);
00580     \endcode
00581 
00582     <b> Required Interface:</b>
00583 
00584     \code
00585     Kernel kernel;
00586     int kernelRadius = kernel.radius();
00587     double y = ...;  // must be <= radius()
00588     double value = kernel(y);
00589     \endcode
00590 */
00591 doxygen_overloaded_function(template <...> void resamplingConvolveY)
00592 
00593 template <class SrcIter, class SrcAcc,
00594           class DestIter, class DestAcc,
00595           class Kernel>
00596 void
00597 resamplingConvolveY(SrcIter sul, SrcIter slr, SrcAcc src,
00598                     DestIter dul, DestIter dlr, DestAcc dest,
00599                     Kernel const & kernel,
00600                     Rational<int> const & samplingRatio, Rational<int> const & offset)
00601 {
00602     int hold = slr.y - sul.y;
00603     int hnew = dlr.y - dul.y;
00604 
00605     vigra_precondition(!samplingRatio.is_inf() && samplingRatio > 0,
00606                 "resamplingConvolveY(): sampling ratio must be > 0 and < infinity");
00607     vigra_precondition(!offset.is_inf(),
00608                 "resamplingConvolveY(): offset must be < infinity");
00609 
00610     int period = lcm(samplingRatio.numerator(), samplingRatio.denominator());
00611 
00612     resampling_detail::MapTargetToSourceCoordinate mapCoordinate(samplingRatio, offset);
00613 
00614     ArrayVector<Kernel1D<double> > kernels(period);
00615 
00616     createResamplingKernels(kernel, mapCoordinate, kernels);
00617 
00618     for(; sul.x < slr.x; ++sul.x, ++dul.x)
00619     {
00620         typename SrcIter::column_iterator sc = sul.columnIterator();
00621         typename DestIter::column_iterator dc = dul.columnIterator();
00622         resamplingConvolveLine(sc, sc+hold, src, dc, dc+hnew, dest,
00623                                kernels, mapCoordinate);
00624     }
00625 }
00626 
00627 template <class SrcIter, class SrcAcc,
00628           class DestIter, class DestAcc,
00629           class Kernel>
00630 inline void
00631 resamplingConvolveY(triple<SrcIter, SrcIter, SrcAcc> src,
00632                     triple<DestIter, DestIter, DestAcc> dest,
00633                     Kernel const & kernel,
00634                     Rational<int> const & samplingRatio, Rational<int> const & offset)
00635 {
00636     resamplingConvolveY(src.first, src.second, src.third,
00637                         dest.first, dest.second, dest.third,
00638                         kernel, samplingRatio, offset);
00639 }
00640 
00641 /********************************************************/
00642 /*                                                      */
00643 /*               resamplingConvolveImage                */
00644 /*                                                      */
00645 /********************************************************/
00646 
00647 /** \brief Apply two separable resampling filters successively, the first in x-direction,
00648            the second in y-direction.
00649 
00650     This function is a shorthand for the concatenation of a call to
00651     \ref resamplingConvolveX() and \ref resamplingConvolveY()
00652     with the given kernels. See there for detailed documentation.
00653 
00654     <b> Declarations:</b>
00655 
00656     pass arguments explicitly:
00657     \code
00658     namespace vigra {
00659         template <class SrcIterator, class SrcAccessor,
00660                   class DestIterator, class DestAccessor,
00661                   class KernelX, class KernelY>
00662         void resamplingConvolveImage(SrcIterator sul,SrcIterator slr, SrcAccessor src,
00663                            DestIterator dul, DestIterator dlr, DestAccessor dest,
00664                            KernelX const & kx,
00665                            Rational<int> const & samplingRatioX, Rational<int> const & offsetX,
00666                            KernelY const & ky,
00667                            Rational<int> const & samplingRatioY, Rational<int> const & offsetY);
00668     }
00669     \endcode
00670 
00671 
00672     use argument objects in conjunction with \ref ArgumentObjectFactories :
00673     \code
00674     namespace vigra {
00675         template <class SrcIterator, class SrcAccessor,
00676                   class DestIterator, class DestAccessor,
00677                   class KernelX, class KernelY>
00678         void
00679         resamplingConvolveImage(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00680                            triple<DestIterator, DestIterator, DestAccessor> dest,
00681                            KernelX const & kx,
00682                            Rational<int> const & samplingRatioX, Rational<int> const & offsetX,
00683                            KernelY const & ky,
00684                            Rational<int> const & samplingRatioY, Rational<int> const & offsetY);
00685     }
00686     \endcode
00687 
00688     <b> Usage:</b>
00689 
00690     <b>\#include</b> <vigra/resampling_convolution.hxx>
00691 
00692 
00693     \code
00694     Rational<int> xratio(2), yratio(3), offset(0);
00695 
00696     FImage src(w,h),
00697            dest(rational_cast<int>(xratio*w), rational_cast<int>(yratio*h));
00698 
00699     float sigma = 2.0;
00700     Gaussian<float> smooth(sigma);
00701     ...
00702 
00703     // simultaneously enlarge and smooth source image
00704     resamplingConvolveImage(srcImageRange(src), destImageRange(dest),
00705                             smooth, xratio, offset,
00706                             smooth, yratio, offset);
00707 
00708     \endcode
00709 
00710 */
00711 doxygen_overloaded_function(template <...> void resamplingConvolveImage)
00712 
00713 template <class SrcIterator, class SrcAccessor,
00714           class DestIterator, class DestAccessor,
00715           class KernelX, class KernelY>
00716 void resamplingConvolveImage(SrcIterator sul,SrcIterator slr, SrcAccessor src,
00717                    DestIterator dul, DestIterator dlr, DestAccessor dest,
00718                    KernelX const & kx,
00719                    Rational<int> const & samplingRatioX, Rational<int> const & offsetX,
00720                    KernelY const & ky,
00721                    Rational<int> const & samplingRatioY, Rational<int> const & offsetY)
00722 {
00723     typedef typename
00724         NumericTraits<typename SrcAccessor::value_type>::RealPromote
00725         TmpType;
00726 
00727     BasicImage<TmpType> tmp(dlr.x - dul.x, slr.y - sul.y);
00728 
00729     resamplingConvolveX(srcIterRange(sul, slr, src),
00730                         destImageRange(tmp),
00731                         kx, samplingRatioX, offsetX);
00732     resamplingConvolveY(srcImageRange(tmp),
00733                         destIterRange(dul, dlr, dest),
00734                         ky, samplingRatioY, offsetY);
00735 }
00736 
00737 template <class SrcIterator, class SrcAccessor,
00738           class DestIterator, class DestAccessor,
00739           class KernelX, class KernelY>
00740 inline void
00741 resamplingConvolveImage(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00742                    triple<DestIterator, DestIterator, DestAccessor> dest,
00743                    KernelX const & kx,
00744                    Rational<int> const & samplingRatioX, Rational<int> const & offsetX,
00745                    KernelY const & ky,
00746                    Rational<int> const & samplingRatioY, Rational<int> const & offsetY)
00747 {
00748     resamplingConvolveImage(src.first, src.second, src.third,
00749                             dest.first, dest.second, dest.third,
00750                             kx, samplingRatioX, offsetX,
00751                             ky, samplingRatioY, offsetY);
00752 }
00753 
00754 /** \brief Two-fold down-sampling for image pyramid construction.
00755 
00756     Sorry, no \ref detailedDocumentation() available yet.
00757 
00758     <b> Declarations:</b>
00759 
00760     <b>\#include</b> <vigra/resampling_convolution.hxx><br>
00761     Namespace: vigra
00762 
00763     pass arguments explicitly:
00764     \code
00765     namespace vigra {
00766         template <class SrcIterator, class SrcAccessor,
00767                   class DestIterator, class DestAccessor>
00768         void pyramidReduceBurtFilter(SrcIterator sul, SrcIterator slr, SrcAccessor src,
00769                                      DestIterator dul, DestIterator dlr, DestAccessor dest,
00770                                      double centerValue = 0.4);
00771     }
00772     \endcode
00773 
00774     use argument objects in conjunction with \ref ArgumentObjectFactories :
00775     \code
00776     namespace vigra {
00777         template <class SrcIterator, class SrcAccessor,
00778                   class DestIterator, class DestAccessor>
00779         void pyramidReduceBurtFilter(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00780                                      triple<DestIterator, DestIterator, DestAccessor> dest,
00781                                      double centerValue = 0.4);
00782     }
00783     \endcode
00784 
00785     use a \ref vigra::ImagePyramid :
00786     \code
00787     namespace vigra {
00788         template <class Image, class Alloc>
00789         void pyramidReduceBurtFilter(ImagePyramid<Image, Alloc> & pyramid, int fromLevel, int toLevel,
00790                                      double centerValue = 0.4);
00791     }
00792     \endcode
00793 */
00794 doxygen_overloaded_function(template <...> void pyramidReduceBurtFilter)
00795 
00796 template <class SrcIterator, class SrcAccessor,
00797           class DestIterator, class DestAccessor>
00798 void pyramidReduceBurtFilter(SrcIterator sul, SrcIterator slr, SrcAccessor src,
00799                              DestIterator dul, DestIterator dlr, DestAccessor dest,
00800                              double centerValue = 0.4)
00801 {
00802     vigra_precondition(0.25 <= centerValue && centerValue <= 0.5,
00803              "pyramidReduceBurtFilter(): centerValue must be between 0.25 and 0.5.");
00804              
00805     int wold = slr.x - sul.x;
00806     int wnew = dlr.x - dul.x;
00807     int hold = slr.y - sul.y;
00808     int hnew = dlr.y - dul.y;
00809     
00810     vigra_precondition(wnew == (wold + 1) / 2 && hnew == (hold + 1) / 2,
00811        "pyramidReduceBurtFilter(): oldSize = ceil(newSize / 2) required.");
00812     
00813     vigra_precondition(wnew == (wold + 1) / 2 && hnew == (hold + 1) / 2,
00814        "pyramidReduceBurtFilter(): oldSize = ceil(newSize / 2) required.");
00815     
00816     Rational<int> samplingRatio(1,2), offset(0);
00817     resampling_detail::MapTargetToSourceCoordinate mapCoordinate(samplingRatio, offset);
00818     
00819     ArrayVector<Kernel1D<double> > kernels(1);
00820     kernels[0].initExplicitly(-2, 2) = 0.25 - centerValue / 2.0, 0.25, centerValue, 0.25, 0.25 - centerValue / 2.0;
00821    
00822     typedef typename
00823         NumericTraits<typename SrcAccessor::value_type>::RealPromote
00824         TmpType;
00825     typedef BasicImage<TmpType> TmpImage;
00826     typedef typename TmpImage::traverser TmpIterator;
00827     
00828     BasicImage<TmpType> tmp(wnew, hold);
00829     
00830     TmpIterator tul = tmp.upperLeft();
00831 
00832     for(; sul.y < slr.y; ++sul.y, ++tul.y)
00833     {
00834         typename SrcIterator::row_iterator sr = sul.rowIterator();
00835         typename TmpIterator::row_iterator tr = tul.rowIterator();
00836         // FIXME: replace with reduceLineBurtFilter()
00837         resamplingConvolveLine(sr, sr+wold, src, tr, tr+wnew, tmp.accessor(),
00838                                kernels, mapCoordinate);
00839     }
00840     
00841     tul  = tmp.upperLeft();
00842 
00843     for(; dul.x < dlr.x; ++dul.x, ++tul.x)
00844     {
00845         typename DestIterator::column_iterator dc = dul.columnIterator();
00846         typename TmpIterator::column_iterator tc = tul.columnIterator();
00847         resamplingConvolveLine(tc, tc+hold, tmp.accessor(), dc, dc+hnew, dest,
00848                                kernels, mapCoordinate);
00849     }
00850 }
00851 
00852 template <class SrcIterator, class SrcAccessor,
00853           class DestIterator, class DestAccessor>
00854 inline
00855 void pyramidReduceBurtFilter(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00856                              triple<DestIterator, DestIterator, DestAccessor> dest,
00857                              double centerValue = 0.4)
00858 {
00859     pyramidReduceBurtFilter(src.first, src.second, src.third, 
00860                             dest.first, dest.second, dest.third, centerValue);
00861 }
00862 
00863 template <class Image, class Alloc>
00864 inline
00865 void pyramidReduceBurtFilter(ImagePyramid<Image, Alloc> & pyramid, int fromLevel, int toLevel,
00866                              double centerValue = 0.4)
00867 {
00868     vigra_precondition(fromLevel  < toLevel,
00869        "pyramidReduceBurtFilter(): fromLevel must be smaller than toLevel.");
00870     vigra_precondition(pyramid.lowestLevel() <= fromLevel && toLevel <= pyramid.highestLevel(),
00871        "pyramidReduceBurtFilter(): fromLevel and toLevel must be between the lowest and highest pyramid levels (inclusive).");
00872 
00873     for(int i=fromLevel+1; i <= toLevel; ++i)
00874         pyramidReduceBurtFilter(srcImageRange(pyramid[i-1]), destImageRange(pyramid[i]), centerValue);
00875 }
00876 
00877 /** \brief Two-fold up-sampling for image pyramid reconstruction.
00878 
00879     Sorry, no \ref detailedDocumentation() available yet.
00880 
00881     <b> Declarations:</b>
00882 
00883     <b>\#include</b> <vigra/resampling_convolution.hxx><br>
00884     Namespace: vigra
00885 
00886     pass arguments explicitly:
00887     \code
00888     namespace vigra {
00889         template <class SrcIterator, class SrcAccessor,
00890                   class DestIterator, class DestAccessor>
00891         void pyramidExpandBurtFilter(SrcIterator sul, SrcIterator slr, SrcAccessor src,
00892                                      DestIterator dul, DestIterator dlr, DestAccessor dest,
00893                                      double centerValue = 0.4);
00894     }
00895     \endcode
00896 
00897 
00898     use argument objects in conjunction with \ref ArgumentObjectFactories :
00899     \code
00900     namespace vigra {
00901         template <class SrcIterator, class SrcAccessor,
00902                   class DestIterator, class DestAccessor>
00903         void pyramidExpandBurtFilter(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00904                                      triple<DestIterator, DestIterator, DestAccessor> dest,
00905                                      double centerValue = 0.4);
00906     }
00907     \endcode
00908 
00909     use a \ref vigra::ImagePyramid :
00910     \code
00911     namespace vigra {
00912         template <class Image, class Alloc>
00913         void pyramidExpandBurtFilter(ImagePyramid<Image, Alloc> & pyramid, int fromLevel, int toLevel,
00914                                      double centerValue = 0.4);
00915     }
00916     \endcode
00917 */
00918 doxygen_overloaded_function(template <...> void pyramidExpandBurtFilter)
00919 
00920 template <class SrcIterator, class SrcAccessor,
00921           class DestIterator, class DestAccessor>
00922 void pyramidExpandBurtFilter(SrcIterator sul, SrcIterator slr, SrcAccessor src,
00923                              DestIterator dul, DestIterator dlr, DestAccessor dest,
00924                              double centerValue = 0.4)
00925 {
00926     vigra_precondition(0.25 <= centerValue && centerValue <= 0.5,
00927              "pyramidExpandBurtFilter(): centerValue must be between 0.25 and 0.5.");
00928              
00929     int wold = slr.x - sul.x;
00930     int wnew = dlr.x - dul.x;
00931     int hold = slr.y - sul.y;
00932     int hnew = dlr.y - dul.y;
00933     
00934     vigra_precondition(wold == (wnew + 1) / 2 && hold == (hnew + 1) / 2,
00935        "pyramidExpandBurtFilter(): oldSize = ceil(newSize / 2) required.");
00936     
00937     vigra_precondition(wold == (wnew + 1) / 2 && hold == (hnew + 1) / 2,
00938        "pyramidExpandBurtFilter(): oldSize = ceil(newSize / 2) required.");
00939     
00940     Rational<int> samplingRatio(2), offset(0);
00941     resampling_detail::MapTargetToSourceCoordinate mapCoordinate(samplingRatio, offset);
00942     
00943     ArrayVector<Kernel1D<double> > kernels(2);
00944     kernels[0].initExplicitly(-1, 1) = 0.5 - centerValue, 2.0*centerValue, 0.5 - centerValue;
00945     kernels[1].initExplicitly(-1, 0) = 0.5, 0.5;
00946    
00947     typedef typename
00948         NumericTraits<typename SrcAccessor::value_type>::RealPromote
00949         TmpType;
00950     typedef BasicImage<TmpType> TmpImage;
00951     typedef typename TmpImage::traverser TmpIterator;
00952     
00953     BasicImage<TmpType> tmp(wnew, hold);
00954     
00955     TmpIterator tul = tmp.upperLeft();
00956 
00957     for(; sul.y < slr.y; ++sul.y, ++tul.y)
00958     {
00959         typename SrcIterator::row_iterator sr = sul.rowIterator();
00960         typename TmpIterator::row_iterator tr = tul.rowIterator();
00961         // FIXME: replace with expandLineBurtFilter()
00962         resamplingConvolveLine(sr, sr+wold, src, tr, tr+wnew, tmp.accessor(),
00963                                kernels, mapCoordinate);
00964     }
00965     
00966     tul  = tmp.upperLeft();
00967 
00968     for(; dul.x < dlr.x; ++dul.x, ++tul.x)
00969     {
00970         typename DestIterator::column_iterator dc = dul.columnIterator();
00971         typename TmpIterator::column_iterator tc = tul.columnIterator();
00972         resamplingConvolveLine(tc, tc+hold, tmp.accessor(), dc, dc+hnew, dest,
00973                                kernels, mapCoordinate);
00974     }
00975 }
00976 
00977 template <class SrcIterator, class SrcAccessor,
00978           class DestIterator, class DestAccessor>
00979 inline
00980 void pyramidExpandBurtFilter(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00981                              triple<DestIterator, DestIterator, DestAccessor> dest,
00982                              double centerValue = 0.4)
00983 {
00984     pyramidExpandBurtFilter(src.first, src.second, src.third, 
00985                             dest.first, dest.second, dest.third, centerValue);
00986 }
00987 
00988 
00989 template <class Image, class Alloc>
00990 inline
00991 void pyramidExpandBurtFilter(ImagePyramid<Image, Alloc> & pyramid, int fromLevel, int toLevel,
00992                              double centerValue = 0.4)
00993 {
00994     vigra_precondition(fromLevel  > toLevel,
00995        "pyramidExpandBurtFilter(): fromLevel must be larger than toLevel.");
00996     vigra_precondition(pyramid.lowestLevel() <= toLevel && fromLevel <= pyramid.highestLevel(),
00997        "pyramidExpandBurtFilter(): fromLevel and toLevel must be between the lowest and highest pyramid levels (inclusive).");
00998 
00999     for(int i=fromLevel-1; i >= toLevel; --i)
01000         pyramidExpandBurtFilter(srcImageRange(pyramid[i+1]), destImageRange(pyramid[i]), centerValue);
01001 }
01002 
01003 /** \brief Create a Laplacian pyramid.
01004 
01005     Sorry, no \ref detailedDocumentation() available yet.
01006 
01007     <b>\#include</b> <vigra/resampling_convolution.hxx><br>
01008     Namespace: vigra
01009 */
01010 template <class Image, class Alloc>
01011 inline
01012 void pyramidReduceBurtLaplacian(ImagePyramid<Image, Alloc> & pyramid, int fromLevel, int toLevel,
01013                              double centerValue = 0.4)
01014 {
01015     using namespace functor;
01016     
01017     pyramidReduceBurtFilter(pyramid, fromLevel, toLevel, centerValue);
01018     for(int i=fromLevel; i < toLevel; ++i)
01019     {
01020         typename ImagePyramid<Image, Alloc>::value_type tmpImage(pyramid[i].size());
01021         pyramidExpandBurtFilter(srcImageRange(pyramid[i+1]), destImageRange(tmpImage), centerValue);
01022         combineTwoImages(srcImageRange(tmpImage), srcImage(pyramid[i]), destImage(pyramid[i]),
01023                        Arg1() - Arg2()); 
01024     }
01025 }
01026 
01027 /** \brief Reconstruct a Laplacian pyramid.
01028 
01029     Sorry, no \ref detailedDocumentation() available yet.
01030 
01031     <b>\#include</b> <vigra/resampling_convolution.hxx><br>
01032     Namespace: vigra
01033 */
01034 template <class Image, class Alloc>
01035 inline
01036 void pyramidExpandBurtLaplacian(ImagePyramid<Image, Alloc> & pyramid, int fromLevel, int toLevel,
01037                                 double centerValue = 0.4)
01038 {
01039     using namespace functor;
01040     
01041     vigra_precondition(fromLevel  > toLevel,
01042        "pyramidExpandBurtLaplacian(): fromLevel must be larger than toLevel.");
01043     vigra_precondition(pyramid.lowestLevel() <= toLevel && fromLevel <= pyramid.highestLevel(),
01044        "pyramidExpandBurtLaplacian(): fromLevel and toLevel must be between the lowest and highest pyramid levels (inclusive).");
01045 
01046     for(int i=fromLevel-1; i >= toLevel; --i)
01047     {
01048         typename ImagePyramid<Image, Alloc>::value_type tmpImage(pyramid[i].size());
01049         pyramidExpandBurtFilter(srcImageRange(pyramid[i+1]), destImageRange(tmpImage), centerValue);
01050         combineTwoImages(srcImageRange(tmpImage), srcImage(pyramid[i]), destImage(pyramid[i]),
01051                        Arg1() - Arg2()); 
01052     }
01053 }
01054 
01055 //@}
01056 
01057 } // namespace vigra
01058 
01059 
01060 #endif /* VIGRA_RESAMPLING_CONVOLUTION_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)