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