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

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