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

vigra/separableconvolution.hxx

00001 /************************************************************************/
00002 /*                                                                      */
00003 /*               Copyright 1998-2002 by Ullrich Koethe                  */
00004 /*       Cognitive Systems Group, University of Hamburg, Germany        */
00005 /*                                                                      */
00006 /*    This file is part of the VIGRA computer vision library.           */
00007 /*    ( Version 1.6.0, Aug 13 2008 )                                    */
00008 /*    The VIGRA Website is                                              */
00009 /*        http://kogs-www.informatik.uni-hamburg.de/~koethe/vigra/      */
00010 /*    Please direct questions, bug reports, and contributions to        */
00011 /*        ullrich.koethe@iwr.uni-heidelberg.de    or                    */
00012 /*        vigra@informatik.uni-hamburg.de                               */
00013 /*                                                                      */
00014 /*    Permission is hereby granted, free of charge, to any person       */
00015 /*    obtaining a copy of this software and associated documentation    */
00016 /*    files (the "Software"), to deal in the Software without           */
00017 /*    restriction, including without limitation the rights to use,      */
00018 /*    copy, modify, merge, publish, distribute, sublicense, and/or      */
00019 /*    sell copies of the Software, and to permit persons to whom the    */
00020 /*    Software is furnished to do so, subject to the following          */
00021 /*    conditions:                                                       */
00022 /*                                                                      */
00023 /*    The above copyright notice and this permission notice shall be    */
00024 /*    included in all copies or substantial portions of the             */
00025 /*    Software.                                                         */
00026 /*                                                                      */
00027 /*    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND    */
00028 /*    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   */
00029 /*    OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          */
00030 /*    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT       */
00031 /*    HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,      */
00032 /*    WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING      */
00033 /*    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR     */
00034 /*    OTHER DEALINGS IN THE SOFTWARE.                                   */
00035 /*                                                                      */
00036 /************************************************************************/
00037 
00038 
00039 #ifndef VIGRA_SEPARABLECONVOLUTION_HXX
00040 #define VIGRA_SEPARABLECONVOLUTION_HXX
00041 
00042 #include <cmath>
00043 #include "utilities.hxx"
00044 #include "numerictraits.hxx"
00045 #include "imageiteratoradapter.hxx"
00046 #include "bordertreatment.hxx"
00047 #include "gaussians.hxx"
00048 #include "array_vector.hxx"
00049 
00050 namespace vigra {
00051 
00052 /********************************************************/
00053 /*                                                      */
00054 /*                internalConvolveLineWrap              */
00055 /*                                                      */
00056 /********************************************************/
00057 
00058 template <class SrcIterator, class SrcAccessor,
00059           class DestIterator, class DestAccessor,
00060           class KernelIterator, class KernelAccessor>
00061 void internalConvolveLineWrap(SrcIterator is, SrcIterator iend, SrcAccessor sa,
00062                               DestIterator id, DestAccessor da,
00063                               KernelIterator kernel, KernelAccessor ka,
00064                               int kleft, int kright)
00065 {
00066   //    int w = iend - is;
00067     int w = std::distance( is, iend );
00068 
00069     typedef typename NumericTraits<typename
00070                       SrcAccessor::value_type>::RealPromote SumType;
00071 
00072     SrcIterator ibegin = is;
00073 
00074     for(int x=0; x<w; ++x, ++is, ++id)
00075     {
00076         KernelIterator ik = kernel + kright;
00077         SumType sum = NumericTraits<SumType>::zero();
00078 
00079         if(x < kright)
00080         {
00081             int x0 = x - kright;
00082             SrcIterator iss = iend + x0;
00083 
00084             for(; x0; ++x0, --ik, ++iss)
00085             {
00086                 sum += ka(ik) * sa(iss);
00087             }
00088 
00089             iss = ibegin;
00090             SrcIterator isend = is + (1 - kleft);
00091             for(; iss != isend ; --ik, ++iss)
00092             {
00093                 sum += ka(ik) * sa(iss);
00094             }
00095         }
00096         else if(w-x <= -kleft)
00097         {
00098             SrcIterator iss = is + (-kright);
00099             SrcIterator isend = iend;
00100             for(; iss != isend ; --ik, ++iss)
00101             {
00102                 sum += ka(ik) * sa(iss);
00103             }
00104 
00105             int x0 = -kleft - w + x + 1;
00106             iss = ibegin;
00107 
00108             for(; x0; --x0, --ik, ++iss)
00109             {
00110                 sum += ka(ik) * sa(iss);
00111             }
00112         }
00113         else
00114         {
00115             SrcIterator iss = is - kright;
00116             SrcIterator isend = is + (1 - kleft);
00117             for(; iss != isend ; --ik, ++iss)
00118             {
00119                 sum += ka(ik) * sa(iss);
00120             }
00121         }
00122 
00123         da.set(NumericTraits<typename
00124                       DestAccessor::value_type>::fromRealPromote(sum), id);
00125     }
00126 }
00127 
00128 /********************************************************/
00129 /*                                                      */
00130 /*                internalConvolveLineClip              */
00131 /*                                                      */
00132 /********************************************************/
00133 
00134 template <class SrcIterator, class SrcAccessor,
00135           class DestIterator, class DestAccessor,
00136           class KernelIterator, class KernelAccessor,
00137           class Norm>
00138 void internalConvolveLineClip(SrcIterator is, SrcIterator iend, SrcAccessor sa,
00139                               DestIterator id, DestAccessor da,
00140                               KernelIterator kernel, KernelAccessor ka,
00141                               int kleft, int kright, Norm norm)
00142 {
00143   //    int w = iend - is;
00144     int w = std::distance( is, iend );
00145 
00146     typedef typename NumericTraits<typename
00147                       SrcAccessor::value_type>::RealPromote SumType;
00148 
00149     SrcIterator ibegin = is;
00150 
00151     for(int x=0; x<w; ++x, ++is, ++id)
00152     {
00153         KernelIterator ik = kernel + kright;
00154         SumType sum = NumericTraits<SumType>::zero();
00155 
00156         if(x < kright)
00157         {
00158             int x0 = x - kright;
00159             Norm clipped = NumericTraits<Norm>::zero();
00160 
00161             for(; x0; ++x0, --ik)
00162             {
00163                 clipped += ka(ik);
00164             }
00165 
00166             SrcIterator iss = ibegin;
00167             SrcIterator isend = is + (1 - kleft);
00168             for(; iss != isend ; --ik, ++iss)
00169             {
00170                 sum += ka(ik) * sa(iss);
00171             }
00172 
00173             sum = norm / (norm - clipped) * sum;
00174         }
00175         else if(w-x <= -kleft)
00176         {
00177             SrcIterator iss = is + (-kright);
00178             SrcIterator isend = iend;
00179             for(; iss != isend ; --ik, ++iss)
00180             {
00181                 sum += ka(ik) * sa(iss);
00182             }
00183 
00184             Norm clipped = NumericTraits<Norm>::zero();
00185 
00186             int x0 = -kleft - w + x + 1;
00187 
00188             for(; x0; --x0, --ik)
00189             {
00190                 clipped += ka(ik);
00191             }
00192 
00193             sum = norm / (norm - clipped) * sum;
00194         }
00195         else
00196         {
00197             SrcIterator iss = is + (-kright);
00198             SrcIterator isend = is + (1 - kleft);
00199             for(; iss != isend ; --ik, ++iss)
00200             {
00201                 sum += ka(ik) * sa(iss);
00202             }
00203         }
00204 
00205         da.set(NumericTraits<typename
00206                       DestAccessor::value_type>::fromRealPromote(sum), id);
00207     }
00208 }
00209 
00210 /********************************************************/
00211 /*                                                      */
00212 /*             internalConvolveLineReflect              */
00213 /*                                                      */
00214 /********************************************************/
00215 
00216 template <class SrcIterator, class SrcAccessor,
00217           class DestIterator, class DestAccessor,
00218           class KernelIterator, class KernelAccessor>
00219 void internalConvolveLineReflect(SrcIterator is, SrcIterator iend, SrcAccessor sa,
00220                               DestIterator id, DestAccessor da,
00221                               KernelIterator kernel, KernelAccessor ka,
00222                               int kleft, int kright)
00223 {
00224   //    int w = iend - is;
00225     int w = std::distance( is, iend );
00226 
00227     typedef typename NumericTraits<typename
00228                       SrcAccessor::value_type>::RealPromote SumType;
00229 
00230     SrcIterator ibegin = is;
00231 
00232     for(int x=0; x<w; ++x, ++is, ++id)
00233     {
00234         KernelIterator ik = kernel + kright;
00235         SumType sum = NumericTraits<SumType>::zero();
00236 
00237         if(x < kright)
00238         {
00239             int x0 = x - kright;
00240             SrcIterator iss = ibegin - x0;
00241 
00242             for(; x0; ++x0, --ik, --iss)
00243             {
00244                 sum += ka(ik) * sa(iss);
00245             }
00246 
00247             SrcIterator isend = is + (1 - kleft);
00248             for(; iss != isend ; --ik, ++iss)
00249             {
00250                 sum += ka(ik) * sa(iss);
00251             }
00252         }
00253         else if(w-x <= -kleft)
00254         {
00255             SrcIterator iss = is + (-kright);
00256             SrcIterator isend = iend;
00257             for(; iss != isend ; --ik, ++iss)
00258             {
00259                 sum += ka(ik) * sa(iss);
00260             }
00261 
00262             int x0 = -kleft - w + x + 1;
00263             iss = iend - 2;
00264 
00265             for(; x0; --x0, --ik, --iss)
00266             {
00267                 sum += ka(ik) * sa(iss);
00268             }
00269         }
00270         else
00271         {
00272             SrcIterator iss = is + (-kright);
00273             SrcIterator isend = is + (1 - kleft);
00274             for(; iss != isend ; --ik, ++iss)
00275             {
00276                 sum += ka(ik) * sa(iss);
00277             }
00278         }
00279 
00280         da.set(NumericTraits<typename
00281                       DestAccessor::value_type>::fromRealPromote(sum), id);
00282     }
00283 }
00284 
00285 /********************************************************/
00286 /*                                                      */
00287 /*             internalConvolveLineRepeat               */
00288 /*                                                      */
00289 /********************************************************/
00290 
00291 template <class SrcIterator, class SrcAccessor,
00292           class DestIterator, class DestAccessor,
00293           class KernelIterator, class KernelAccessor>
00294 void internalConvolveLineRepeat(SrcIterator is, SrcIterator iend, SrcAccessor sa,
00295                               DestIterator id, DestAccessor da,
00296                               KernelIterator kernel, KernelAccessor ka,
00297                               int kleft, int kright)
00298 {
00299   //    int w = iend - is;
00300     int w = std::distance( is, iend );
00301 
00302     typedef typename NumericTraits<typename
00303                       SrcAccessor::value_type>::RealPromote SumType;
00304 
00305     SrcIterator ibegin = is;
00306 
00307     for(int x=0; x<w; ++x, ++is, ++id)
00308     {
00309         KernelIterator ik = kernel + kright;
00310         SumType sum = NumericTraits<SumType>::zero();
00311 
00312         if(x < kright)
00313         {
00314             int x0 = x - kright;
00315             SrcIterator iss = ibegin;
00316 
00317             for(; x0; ++x0, --ik)
00318             {
00319                 sum += ka(ik) * sa(iss);
00320             }
00321 
00322             SrcIterator isend = is + (1 - kleft);
00323             for(; iss != isend ; --ik, ++iss)
00324             {
00325                 sum += ka(ik) * sa(iss);
00326             }
00327         }
00328         else if(w-x <= -kleft)
00329         {
00330             SrcIterator iss = is + (-kright);
00331             SrcIterator isend = iend;
00332             for(; iss != isend ; --ik, ++iss)
00333             {
00334                 sum += ka(ik) * sa(iss);
00335             }
00336 
00337             int x0 = -kleft - w + x + 1;
00338             iss = iend - 1;
00339 
00340             for(; x0; --x0, --ik)
00341             {
00342                 sum += ka(ik) * sa(iss);
00343             }
00344         }
00345         else
00346         {
00347             SrcIterator iss = is + (-kright);
00348             SrcIterator isend = is + (1 - kleft);
00349             for(; iss != isend ; --ik, ++iss)
00350             {
00351                 sum += ka(ik) * sa(iss);
00352             }
00353         }
00354 
00355         da.set(NumericTraits<typename
00356                       DestAccessor::value_type>::fromRealPromote(sum), id);
00357     }
00358 }
00359 
00360 /********************************************************/
00361 /*                                                      */
00362 /*              internalConvolveLineAvoid               */
00363 /*                                                      */
00364 /********************************************************/
00365 
00366 template <class SrcIterator, class SrcAccessor,
00367           class DestIterator, class DestAccessor,
00368           class KernelIterator, class KernelAccessor>
00369 void internalConvolveLineAvoid(SrcIterator is, SrcIterator iend, SrcAccessor sa,
00370                               DestIterator id, DestAccessor da,
00371                               KernelIterator kernel, KernelAccessor ka,
00372                               int kleft, int kright)
00373 {
00374   //    int w = iend - is;
00375     int w = std::distance( is, iend );
00376 
00377     typedef typename NumericTraits<typename
00378                       SrcAccessor::value_type>::RealPromote SumType;
00379 
00380     is += kright;
00381     id += kright;
00382 
00383     for(int x=kright; x<w+kleft; ++x, ++is, ++id)
00384     {
00385         KernelIterator ik = kernel + kright;
00386         SumType sum = NumericTraits<SumType>::zero();
00387 
00388         SrcIterator iss = is + (-kright);
00389         SrcIterator isend = is + (1 - kleft);
00390         for(; iss != isend ; --ik, ++iss)
00391         {
00392             sum += ka(ik) * sa(iss);
00393         }
00394 
00395         da.set(NumericTraits<typename
00396                       DestAccessor::value_type>::fromRealPromote(sum), id);
00397     }
00398 }
00399 
00400 /********************************************************/
00401 /*                                                      */
00402 /*         Separable convolution functions              */
00403 /*                                                      */
00404 /********************************************************/
00405 
00406 /** \addtogroup SeparableConvolution One-dimensional and separable convolution functions
00407 
00408     Perform 1D convolution and separable filtering in 2 dimensions.
00409 
00410     These generic convolution functions implement
00411     the standard convolution operation for a wide range of images and
00412     signals that fit into the required interface. They need a suitable
00413     kernel to operate.
00414 */
00415 //@{
00416 
00417 /** \brief Performs a 1-dimensional convolution of the source signal using the given
00418     kernel.
00419 
00420     The KernelIterator must point to the center iterator, and
00421     the kernel's size is given by its left (kleft <= 0) and right
00422     (kright >= 0) borders. The signal must always be larger than the kernel.
00423     At those positions where the kernel does not completely fit
00424     into the signal's range, the specified \ref BorderTreatmentMode is
00425     applied.
00426 
00427     The signal's value_type (SrcAccessor::value_type) must be a
00428     linear space over the kernel's value_type (KernelAccessor::value_type),
00429     i.e. addition of source values, multiplication with kernel values,
00430     and NumericTraits must be defined.
00431     The kernel's value_type must be an algebraic field,
00432     i.e. the arithmetic operations (+, -, *, /) and NumericTraits must
00433     be defined.
00434 
00435     <b> Declarations:</b>
00436 
00437     pass arguments explicitly:
00438     \code
00439     namespace vigra {
00440         template <class SrcIterator, class SrcAccessor,
00441                   class DestIterator, class DestAccessor,
00442                   class KernelIterator, class KernelAccessor>
00443         void convolveLine(SrcIterator is, SrcIterator isend, SrcAccessor sa,
00444                           DestIterator id, DestAccessor da,
00445                           KernelIterator ik, KernelAccessor ka,
00446                           int kleft, int kright, BorderTreatmentMode border)
00447     }
00448     \endcode
00449 
00450 
00451     use argument objects in conjunction with \ref ArgumentObjectFactories :
00452     \code
00453     namespace vigra {
00454         template <class SrcIterator, class SrcAccessor,
00455                   class DestIterator, class DestAccessor,
00456                   class KernelIterator, class KernelAccessor>
00457         void convolveLine(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00458                           pair<DestIterator, DestAccessor> dest,
00459                           tuple5<KernelIterator, KernelAccessor,
00460                                  int, int, BorderTreatmentMode> kernel)
00461     }
00462     \endcode
00463 
00464     <b> Usage:</b>
00465 
00466     <b>\#include</b> <<a href="separableconvolution_8hxx-source.html">vigra/separableconvolution.hxx</a>>
00467 
00468 
00469     \code
00470     std::vector<float> src, dest;
00471     ...
00472 
00473     // define binomial filter of size 5
00474     static float kernel[] =
00475            { 1.0/16.0, 4.0/16.0, 6.0/16.0, 4.0/16.0, 1.0/16.0};
00476 
00477     typedef vigra::StandardAccessor<float> FAccessor;
00478     typedef vigra::StandardAccessor<float> KernelAccessor;
00479 
00480 
00481     vigra::convolveLine(src.begin(), src.end(), FAccessor(), dest.begin(), FAccessor(),
00482              kernel+2, KernelAccessor(), -2, 2, BORDER_TREATMENT_REFLECT);
00483     //       ^^^^^^^^  this is the center of the kernel
00484 
00485     \endcode
00486 
00487     <b> Required Interface:</b>
00488 
00489     \code
00490     RandomAccessIterator is, isend;
00491     RandomAccessIterator id;
00492     RandomAccessIterator ik;
00493 
00494     SrcAccessor src_accessor;
00495     DestAccessor dest_accessor;
00496     KernelAccessor kernel_accessor;
00497 
00498     NumericTraits<SrcAccessor::value_type>::RealPromote s = src_accessor(is);
00499 
00500     s = s + s;
00501     s = kernel_accessor(ik) * s;
00502 
00503     dest_accessor.set(
00504         NumericTraits<DestAccessor::value_type>::fromRealPromote(s), id);
00505 
00506     \endcode
00507 
00508     If border == BORDER_TREATMENT_CLIP:
00509 
00510     \code
00511     NumericTraits<KernelAccessor::value_type>::RealPromote k = kernel_accessor(ik);
00512 
00513     k = k + k;
00514     k = k - k;
00515     k = k * k;
00516     k = k / k;
00517 
00518     \endcode
00519 
00520     <b> Preconditions:</b>
00521 
00522     \code
00523     kleft <= 0
00524     kright >= 0
00525     iend - is >= kright + kleft + 1
00526     \endcode
00527 
00528     If border == BORDER_TREATMENT_CLIP: Sum of kernel elements must be
00529     != 0.
00530 
00531 */
00532 doxygen_overloaded_function(template <...> void convolveLine)
00533 
00534 template <class SrcIterator, class SrcAccessor,
00535           class DestIterator, class DestAccessor,
00536           class KernelIterator, class KernelAccessor>
00537 void convolveLine(SrcIterator is, SrcIterator iend, SrcAccessor sa,
00538                   DestIterator id, DestAccessor da,
00539                   KernelIterator ik, KernelAccessor ka,
00540                   int kleft, int kright, BorderTreatmentMode border)
00541 {
00542     typedef typename KernelAccessor::value_type KernelValue;
00543 
00544     vigra_precondition(kleft <= 0,
00545                  "convolveLine(): kleft must be <= 0.\n");
00546     vigra_precondition(kright >= 0,
00547                  "convolveLine(): kright must be >= 0.\n");
00548 
00549     //    int w = iend - is;
00550     int w = std::distance( is, iend );
00551 
00552     vigra_precondition(w >= kright - kleft + 1,
00553                  "convolveLine(): kernel longer than line\n");
00554 
00555     switch(border)
00556     {
00557       case BORDER_TREATMENT_WRAP:
00558       {
00559         internalConvolveLineWrap(is, iend, sa, id, da, ik, ka, kleft, kright);
00560         break;
00561       }
00562       case BORDER_TREATMENT_AVOID:
00563       {
00564         internalConvolveLineAvoid(is, iend, sa, id, da, ik, ka, kleft, kright);
00565         break;
00566       }
00567       case BORDER_TREATMENT_REFLECT:
00568       {
00569         internalConvolveLineReflect(is, iend, sa, id, da, ik, ka, kleft, kright);
00570         break;
00571       }
00572       case BORDER_TREATMENT_REPEAT:
00573       {
00574         internalConvolveLineRepeat(is, iend, sa, id, da, ik, ka, kleft, kright);
00575         break;
00576       }
00577       case BORDER_TREATMENT_CLIP:
00578       {
00579         // find norm of kernel
00580         typedef typename KernelAccessor::value_type KT;
00581         KT norm = NumericTraits<KT>::zero();
00582         KernelIterator iik = ik + kleft;
00583         for(int i=kleft; i<=kright; ++i, ++iik) norm += ka(iik);
00584 
00585         vigra_precondition(norm != NumericTraits<KT>::zero(),
00586                      "convolveLine(): Norm of kernel must be != 0"
00587                      " in mode BORDER_TREATMENT_CLIP.\n");
00588 
00589         internalConvolveLineClip(is, iend, sa, id, da, ik, ka, kleft, kright, norm);
00590         break;
00591       }
00592       default:
00593       {
00594         vigra_precondition(0,
00595                      "convolveLine(): Unknown border treatment mode.\n");
00596       }
00597     }
00598 }
00599 
00600 template <class SrcIterator, class SrcAccessor,
00601           class DestIterator, class DestAccessor,
00602           class KernelIterator, class KernelAccessor>
00603 inline
00604 void convolveLine(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00605                   pair<DestIterator, DestAccessor> dest,
00606                   tuple5<KernelIterator, KernelAccessor,
00607                          int, int, BorderTreatmentMode> kernel)
00608 {
00609     convolveLine(src.first, src.second, src.third,
00610                  dest.first, dest.second,
00611                  kernel.first, kernel.second,
00612                  kernel.third, kernel.fourth, kernel.fifth);
00613 }
00614 
00615 /********************************************************/
00616 /*                                                      */
00617 /*                      separableConvolveX              */
00618 /*                                                      */
00619 /********************************************************/
00620 
00621 /** \brief Performs a 1 dimensional convolution in x direction.
00622 
00623     It calls \ref convolveLine() for every row of the image. See \ref convolveLine() 
00624     for more information about required interfaces and vigra_preconditions.
00625 
00626     <b> Declarations:</b>
00627 
00628     pass arguments explicitly:
00629     \code
00630     namespace vigra {
00631         template <class SrcImageIterator, class SrcAccessor,
00632                   class DestImageIterator, class DestAccessor,
00633                   class KernelIterator, class KernelAccessor>
00634         void separableConvolveX(SrcImageIterator supperleft,
00635                                 SrcImageIterator slowerright, SrcAccessor sa,
00636                                 DestImageIterator dupperleft, DestAccessor da,
00637                                 KernelIterator ik, KernelAccessor ka,
00638                                 int kleft, int kright, BorderTreatmentMode border)
00639     }
00640     \endcode
00641 
00642 
00643     use argument objects in conjunction with \ref ArgumentObjectFactories :
00644     \code
00645     namespace vigra {
00646         template <class SrcImageIterator, class SrcAccessor,
00647                   class DestImageIterator, class DestAccessor,
00648                   class KernelIterator, class KernelAccessor>
00649         void separableConvolveX(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
00650                                 pair<DestImageIterator, DestAccessor> dest,
00651                                 tuple5<KernelIterator, KernelAccessor,
00652                                              int, int, BorderTreatmentMode> kernel)
00653     }
00654     \endcode
00655 
00656     <b> Usage:</b>
00657 
00658     <b>\#include</b> <<a href="separableconvolution_8hxx-source.html">vigra/separableconvolution.hxx</a>>
00659 
00660 
00661     \code
00662     vigra::FImage src(w,h), dest(w,h);
00663     ...
00664 
00665     // define Gaussian kernel with std. deviation 3.0
00666     vigra::Kernel1D<double> kernel;
00667     kernel.initGaussian(3.0);
00668 
00669     vigra::separableConvolveX(srcImageRange(src), destImage(dest), kernel1d(kernel));
00670 
00671     \endcode
00672 
00673 */
00674 doxygen_overloaded_function(template <...> void separableConvolveX)
00675 
00676 template <class SrcIterator, class SrcAccessor,
00677           class DestIterator, class DestAccessor,
00678           class KernelIterator, class KernelAccessor>
00679 void separableConvolveX(SrcIterator supperleft,
00680                         SrcIterator slowerright, SrcAccessor sa,
00681                         DestIterator dupperleft, DestAccessor da,
00682                         KernelIterator ik, KernelAccessor ka,
00683                         int kleft, int kright, BorderTreatmentMode border)
00684 {
00685     typedef typename KernelAccessor::value_type KernelValue;
00686 
00687     vigra_precondition(kleft <= 0,
00688                  "separableConvolveX(): kleft must be <= 0.\n");
00689     vigra_precondition(kright >= 0,
00690                  "separableConvolveX(): kright must be >= 0.\n");
00691 
00692     int w = slowerright.x - supperleft.x;
00693     int h = slowerright.y - supperleft.y;
00694 
00695     vigra_precondition(w >= kright - kleft + 1,
00696                  "separableConvolveX(): kernel longer than line\n");
00697 
00698     int y;
00699 
00700     for(y=0; y<h; ++y, ++supperleft.y, ++dupperleft.y)
00701     {
00702         typename SrcIterator::row_iterator rs = supperleft.rowIterator();
00703         typename DestIterator::row_iterator rd = dupperleft.rowIterator();
00704 
00705         convolveLine(rs, rs+w, sa, rd, da,
00706                      ik, ka, kleft, kright, border);
00707     }
00708 }
00709 
00710 template <class SrcIterator, class SrcAccessor,
00711           class DestIterator, class DestAccessor,
00712           class KernelIterator, class KernelAccessor>
00713 inline void
00714 separableConvolveX(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00715                   pair<DestIterator, DestAccessor> dest,
00716                   tuple5<KernelIterator, KernelAccessor,
00717                          int, int, BorderTreatmentMode> kernel)
00718 {
00719     separableConvolveX(src.first, src.second, src.third,
00720                  dest.first, dest.second,
00721                  kernel.first, kernel.second,
00722                  kernel.third, kernel.fourth, kernel.fifth);
00723 }
00724 
00725 
00726 
00727 /********************************************************/
00728 /*                                                      */
00729 /*                      separableConvolveY              */
00730 /*                                                      */
00731 /********************************************************/
00732 
00733 /** \brief Performs a 1 dimensional convolution in y direction.
00734 
00735     It calls \ref convolveLine() for every column of the image. See \ref convolveLine() 
00736     for more information about required interfaces and vigra_preconditions.
00737 
00738     <b> Declarations:</b>
00739 
00740     pass arguments explicitly:
00741     \code
00742     namespace vigra {
00743         template <class SrcImageIterator, class SrcAccessor,
00744                   class DestImageIterator, class DestAccessor,
00745                   class KernelIterator, class KernelAccessor>
00746         void separableConvolveY(SrcImageIterator supperleft,
00747                                 SrcImageIterator slowerright, SrcAccessor sa,
00748                                 DestImageIterator dupperleft, DestAccessor da,
00749                                 KernelIterator ik, KernelAccessor ka,
00750                                 int kleft, int kright, BorderTreatmentMode border)
00751     }
00752     \endcode
00753 
00754 
00755     use argument objects in conjunction with \ref ArgumentObjectFactories :
00756     \code
00757     namespace vigra {
00758         template <class SrcImageIterator, class SrcAccessor,
00759                   class DestImageIterator, class DestAccessor,
00760                   class KernelIterator, class KernelAccessor>
00761         void separableConvolveY(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
00762                                 pair<DestImageIterator, DestAccessor> dest,
00763                                 tuple5<KernelIterator, KernelAccessor,
00764                                              int, int, BorderTreatmentMode> kernel)
00765     }
00766     \endcode
00767 
00768     <b> Usage:</b>
00769 
00770     <b>\#include</b> <<a href="separableconvolution_8hxx-source.html">vigra/separableconvolution.hxx</a>>
00771 
00772 
00773     \code
00774     vigra::FImage src(w,h), dest(w,h);
00775     ...
00776 
00777     // define Gaussian kernel with std. deviation 3.0
00778     vigra::Kernel1D kernel;
00779     kernel.initGaussian(3.0);
00780 
00781     vigra::separableConvolveY(srcImageRange(src), destImage(dest), kernel1d(kernel));
00782 
00783     \endcode
00784 
00785 */
00786 doxygen_overloaded_function(template <...> void separableConvolveY)
00787 
00788 template <class SrcIterator, class SrcAccessor,
00789           class DestIterator, class DestAccessor,
00790           class KernelIterator, class KernelAccessor>
00791 void separableConvolveY(SrcIterator supperleft,
00792                         SrcIterator slowerright, SrcAccessor sa,
00793                         DestIterator dupperleft, DestAccessor da,
00794                         KernelIterator ik, KernelAccessor ka,
00795                         int kleft, int kright, BorderTreatmentMode border)
00796 {
00797     typedef typename KernelAccessor::value_type KernelValue;
00798 
00799     vigra_precondition(kleft <= 0,
00800                  "separableConvolveY(): kleft must be <= 0.\n");
00801     vigra_precondition(kright >= 0,
00802                  "separableConvolveY(): kright must be >= 0.\n");
00803 
00804     int w = slowerright.x - supperleft.x;
00805     int h = slowerright.y - supperleft.y;
00806 
00807     vigra_precondition(h >= kright - kleft + 1,
00808                  "separableConvolveY(): kernel longer than line\n");
00809 
00810     int x;
00811 
00812     for(x=0; x<w; ++x, ++supperleft.x, ++dupperleft.x)
00813     {
00814         typename SrcIterator::column_iterator cs = supperleft.columnIterator();
00815         typename DestIterator::column_iterator cd = dupperleft.columnIterator();
00816 
00817         convolveLine(cs, cs+h, sa, cd, da,
00818                      ik, ka, kleft, kright, border);
00819     }
00820 }
00821 
00822 template <class SrcIterator, class SrcAccessor,
00823           class DestIterator, class DestAccessor,
00824           class KernelIterator, class KernelAccessor>
00825 inline void
00826 separableConvolveY(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00827                   pair<DestIterator, DestAccessor> dest,
00828                   tuple5<KernelIterator, KernelAccessor,
00829                          int, int, BorderTreatmentMode> kernel)
00830 {
00831     separableConvolveY(src.first, src.second, src.third,
00832                  dest.first, dest.second,
00833                  kernel.first, kernel.second,
00834                  kernel.third, kernel.fourth, kernel.fifth);
00835 }
00836 
00837 //@}
00838 
00839 /********************************************************/
00840 /*                                                      */
00841 /*                      Kernel1D                        */
00842 /*                                                      */
00843 /********************************************************/
00844 
00845 /** \brief Generic 1 dimensional convolution kernel.
00846 
00847     This kernel may be used for convolution of 1 dimensional signals or for
00848     separable convolution of multidimensional signals.
00849 
00850     Convlution functions access the kernel via a 1 dimensional random access
00851     iterator which they get by calling \ref center(). This iterator
00852     points to the center of the kernel. The kernel's size is given by its left() (<=0)
00853     and right() (>= 0) methods. The desired border treatment mode is
00854     returned by borderTreatment().
00855 
00856     The different init functions create a kernel with the specified
00857     properties. The kernel's value_type must be a linear space, i.e. it
00858     must define multiplication with doubles and NumericTraits.
00859 
00860 
00861     The kernel defines a factory function kernel1d() to create an argument object
00862     (see \ref KernelArgumentObjectFactories).
00863 
00864     <b> Usage:</b>
00865 
00866     <b>\#include</b> <<a href="stdconvolution_8hxx-source.html">vigra/stdconvolution.hxx</a>>
00867 
00868     \code
00869     vigra::FImage src(w,h), dest(w,h);
00870     ...
00871 
00872     // define Gaussian kernel with std. deviation 3.0
00873     vigra::Kernel1D kernel;
00874     kernel.initGaussian(3.0);
00875 
00876     vigra::separableConvolveX(srcImageRange(src), destImage(dest), kernel1d(kernel));
00877     \endcode
00878 
00879     <b> Required Interface:</b>
00880 
00881     \code
00882     value_type v = vigra::NumericTraits<value_type>::one(); // if norm is not
00883                                                             // given explicitly
00884     double d;
00885 
00886     v = d * v;
00887     \endcode
00888 */
00889 
00890 template <class ARITHTYPE>
00891 class Kernel1D
00892 {
00893   public:
00894     typedef ArrayVector<ARITHTYPE> InternalVector;
00895 
00896         /** the kernel's value type
00897         */
00898     typedef typename InternalVector::value_type value_type;
00899 
00900         /** the kernel's reference type
00901         */
00902     typedef typename InternalVector::reference reference;
00903 
00904         /** the kernel's const reference type
00905         */
00906     typedef typename InternalVector::const_reference const_reference;
00907 
00908         /** deprecated -- use Kernel1D::iterator
00909         */
00910     typedef typename InternalVector::iterator Iterator;
00911 
00912         /** 1D random access iterator over the kernel's values
00913         */
00914     typedef typename InternalVector::iterator iterator;
00915 
00916         /** const 1D random access iterator over the kernel's values
00917         */
00918     typedef typename InternalVector::const_iterator const_iterator;
00919 
00920         /** the kernel's accessor
00921         */
00922     typedef StandardAccessor<ARITHTYPE> Accessor;
00923 
00924         /** the kernel's const accessor
00925         */
00926     typedef StandardConstAccessor<ARITHTYPE> ConstAccessor;
00927 
00928     struct InitProxy
00929     {
00930         InitProxy(Iterator i, int count, value_type & norm)
00931         : iter_(i), base_(i),
00932           count_(count), sum_(count),
00933           norm_(norm)
00934         {}
00935 
00936         ~InitProxy()
00937         {
00938             vigra_precondition(count_ == 1 || count_ == sum_,
00939                   "Kernel1D::initExplicitly(): "
00940                   "Wrong number of init values.");
00941         }
00942 
00943         InitProxy & operator,(value_type const & v)
00944         {
00945             if(sum_ == count_) 
00946                 norm_ = *iter_;
00947 
00948             norm_ += v;
00949 
00950             --count_;
00951 
00952             if(count_ > 0)
00953             {
00954                 ++iter_;
00955                 *iter_ = v;
00956             }
00957 
00958             return *this;
00959         }
00960 
00961         Iterator iter_, base_;
00962         int count_, sum_;
00963         value_type & norm_;
00964     };
00965 
00966     static value_type one() { return NumericTraits<value_type>::one(); }
00967 
00968         /** Default constructor.
00969             Creates a kernel of size 1 which would copy the signal
00970             unchanged.
00971         */
00972     Kernel1D()
00973     : kernel_(),
00974       left_(0),
00975       right_(0),
00976       border_treatment_(BORDER_TREATMENT_REFLECT),
00977       norm_(one())
00978     {
00979         kernel_.push_back(norm_);
00980     }
00981 
00982         /** Copy constructor.
00983         */
00984     Kernel1D(Kernel1D const & k)
00985     : kernel_(k.kernel_),
00986       left_(k.left_),
00987       right_(k.right_),
00988       border_treatment_(k.border_treatment_),
00989       norm_(k.norm_)
00990     {}
00991 
00992         /** Copy assignment.
00993         */
00994     Kernel1D & operator=(Kernel1D const & k)
00995     {
00996         if(this != &k)
00997         {
00998             left_ = k.left_;
00999             right_ = k.right_;
01000             border_treatment_ = k.border_treatment_;
01001             norm_ = k.norm_;
01002             kernel_ = k.kernel_;
01003         }
01004         return *this;
01005     }
01006 
01007         /** Initialization.
01008             This initializes the kernel with the given constant. The norm becomes
01009             v*size().
01010 
01011             Instead of a single value an initializer list of length size()
01012             can be used like this:
01013 
01014             \code
01015             vigra::Kernel1D<float> roberts_gradient_x;
01016 
01017             roberts_gradient_x.initExplicitly(0, 1) = 1.0, -1.0;
01018             \endcode
01019 
01020             In this case, the norm will be set to the sum of the init values.
01021             An initializer list of wrong length will result in a run-time error.
01022         */
01023     InitProxy operator=(value_type const & v)
01024     {
01025         int size = right_ - left_ + 1;
01026         for(unsigned int i=0; i<kernel_.size(); ++i) kernel_[i] = v;
01027         norm_ = (double)size*v;
01028 
01029         return InitProxy(kernel_.begin(), size, norm_);
01030     }
01031 
01032         /** Destructor.
01033         */
01034     ~Kernel1D()
01035     {}
01036 
01037         /**
01038             Init as a sampled Gaussian function. The radius of the kernel is
01039             always 3*std_dev. '<tt>norm</tt>' denotes the sum of all bins of the kernel
01040             (i.e. the kernel is corrected for the normalization error introduced
01041              by windowing the Gaussian to a finite interval). However,
01042             if <tt>norm</tt> is 0.0, the kernel is normalized to 1 by the analytic
01043             expression for the Gaussian, and <b>no</b> correction for the windowing
01044             error is performed.
01045 
01046             Precondition:
01047             \code
01048             std_dev >= 0.0
01049             \endcode
01050 
01051             Postconditions:
01052             \code
01053             1. left()  == -(int)(3.0*std_dev + 0.5)
01054             2. right() ==  (int)(3.0*std_dev + 0.5)
01055             3. borderTreatment() == BORDER_TREATMENT_REFLECT
01056             4. norm() == norm
01057             \endcode
01058         */
01059     void initGaussian(double std_dev, value_type norm);
01060 
01061         /** Init as a Gaussian function with norm 1.
01062          */
01063     void initGaussian(double std_dev)
01064     {
01065         initGaussian(std_dev, one());
01066     }
01067 
01068 
01069         /**
01070             Init as Lindeberg's discrete analog of the Gaussian function. The radius of the kernel is
01071             always 3*std_dev. 'norm' denotes the sum of all bins of the kernel.
01072 
01073             Precondition:
01074             \code
01075             std_dev >= 0.0
01076             \endcode
01077 
01078             Postconditions:
01079             \code
01080             1. left()  == -(int)(3.0*std_dev + 0.5)
01081             2. right() ==  (int)(3.0*std_dev + 0.5)
01082             3. borderTreatment() == BORDER_TREATMENT_REFLECT
01083             4. norm() == norm
01084             \endcode
01085         */
01086     void initDiscreteGaussian(double std_dev, value_type norm);
01087 
01088         /** Init as a Lindeberg's discrete analog of the Gaussian function
01089             with norm 1.
01090          */
01091     void initDiscreteGaussian(double std_dev)
01092     {
01093         initDiscreteGaussian(std_dev, one());
01094     }
01095 
01096         /**
01097             Init as a Gaussian derivative of order '<tt>order</tt>'.
01098             The radius of the kernel is always <tt>3*std_dev + 0.5*order</tt>.
01099             '<tt>norm</tt>' denotes the norm of the kernel so that the
01100             following condition is fulfilled:
01101 
01102             \f[ \sum_{i=left()}^{right()}
01103                          \frac{(-i)^{order}kernel[i]}{order!} = norm
01104             \f]
01105 
01106             Thus, the kernel will be corrected for the error introduced
01107             by windowing the Gaussian to a finite interval. However,
01108             if <tt>norm</tt> is 0.0, the kernel is normalized to 1 by the analytic
01109             expression for the Gaussian derivative, and <b>no</b> correction for the
01110             windowing error is performed.
01111 
01112             Preconditions:
01113             \code
01114             1. std_dev >= 0.0
01115             2. order   >= 1
01116             \endcode
01117 
01118             Postconditions:
01119             \code
01120             1. left()  == -(int)(3.0*std_dev + 0.5*order + 0.5)
01121             2. right() ==  (int)(3.0*std_dev + 0.5*order + 0.5)
01122             3. borderTreatment() == BORDER_TREATMENT_REFLECT
01123             4. norm() == norm
01124             \endcode
01125         */
01126     void initGaussianDerivative(double std_dev, int order, value_type norm);
01127 
01128         /** Init as a Gaussian derivative with norm 1.
01129          */
01130     void initGaussianDerivative(double std_dev, int order)
01131     {
01132         initGaussianDerivative(std_dev, order, one());
01133     }
01134 
01135         /**
01136             Init an optimal 3-tap smoothing filter.
01137             The filter values are 
01138             
01139             \code
01140             [0.216, 0.568, 0.216]
01141             \endcode
01142             
01143             These values are optimal in the sense that the 3x3 filter obtained by separable application
01144             of this filter is the best possible 3x3 approximation to a Gaussian filter.
01145             The equivalent Gaussian has sigma = 0.680.
01146  
01147             Postconditions:
01148             \code
01149             1. left()  == -1
01150             2. right() ==  1
01151             3. borderTreatment() == BORDER_TREATMENT_REFLECT
01152             4. norm() == 1.0
01153             \endcode
01154         */
01155     void initOptimalSmoothing3()
01156     {
01157         this->initExplicitly(-1, 1) = 0.216, 0.568, 0.216;
01158         this->setBorderTreatment(BORDER_TREATMENT_REFLECT);
01159     }
01160 
01161         /**
01162             Init an optimal 3-tap smoothing filter to be used in the context of first derivative computation.
01163             This filter must be used in conjunction with the symmetric difference filter (see initSymmetricDifference()), 
01164             such that the difference filter is applied along one dimension, and the smoothing filter along the other.
01165             The filter values are 
01166             
01167             \code
01168             [0.224365, 0.55127, 0.224365]
01169             \endcode
01170             
01171             These values are optimal in the sense that the 3x3 filter obtained by combining 
01172             this filter with the symmetric difference is the best possible 3x3 approximation to a 
01173             Gaussian first derivative filter. The equivalent Gaussian has sigma = 0.675.
01174  
01175             Postconditions:
01176             \code
01177             1. left()  == -1
01178             2. right() ==  1
01179             3. borderTreatment() == BORDER_TREATMENT_REFLECT
01180             4. norm() == 1.0
01181             \endcode
01182         */
01183     void initOptimalFirstDerivativeSmoothing3()
01184     {
01185         this->initExplicitly(-1, 1) = 0.224365, 0.55127, 0.224365;
01186         this->setBorderTreatment(BORDER_TREATMENT_REFLECT);
01187     }
01188 
01189         /**
01190             Init an optimal 3-tap smoothing filter to be used in the context of second derivative computation.
01191             This filter must be used in conjunction with the 3-tap second difference filter (see initSecondDifference3()), 
01192             such that the difference filter is applied along one dimension, and the smoothing filter along the other.
01193             The filter values are 
01194             
01195             \code
01196             [0.13, 0.74, 0.13]
01197             \endcode
01198             
01199             These values are optimal in the sense that the 3x3 filter obtained by combining 
01200             this filter with the 3-tap second difference is the best possible 3x3 approximation to a 
01201             Gaussian second derivative filter. The equivalent Gaussian has sigma = 0.433.
01202  
01203             Postconditions:
01204             \code
01205             1. left()  == -1
01206             2. right() ==  1
01207             3. borderTreatment() == BORDER_TREATMENT_REFLECT
01208             4. norm() == 1.0
01209             \endcode
01210         */
01211     void initOptimalSecondDerivativeSmoothing3()
01212     {
01213         this->initExplicitly(-1, 1) = 0.13, 0.74, 0.13;
01214         this->setBorderTreatment(BORDER_TREATMENT_REFLECT);
01215     }
01216 
01217         /**
01218             Init an optimal 5-tap smoothing filter.
01219             The filter values are 
01220             
01221             \code
01222             [0.03134, 0.24, 0.45732, 0.24, 0.03134]
01223             \endcode
01224             
01225             These values are optimal in the sense that the 5x5 filter obtained by separable application
01226             of this filter is the best possible 5x5 approximation to a Gaussian filter.
01227             The equivalent Gaussian has sigma = 0.867.
01228  
01229             Postconditions:
01230             \code
01231             1. left()  == -2
01232             2. right() ==  2
01233             3. borderTreatment() == BORDER_TREATMENT_REFLECT
01234             4. norm() == 1.0
01235             \endcode
01236         */
01237     void initOptimalSmoothing5()
01238     {
01239         this->initExplicitly(-2, 2) = 0.03134, 0.24, 0.45732, 0.24, 0.03134;
01240         this->setBorderTreatment(BORDER_TREATMENT_REFLECT);
01241     }
01242 
01243         /**
01244             Init an optimal 5-tap smoothing filter to be used in the context of first derivative computation.
01245            This filter must be used in conjunction with the optimal 5-tap first derivative filter 
01246            (see initOptimalFirstDerivative5()),  such that the derivative filter is applied along one dimension, 
01247            and the smoothing filter along the other. The filter values are 
01248             
01249             \code
01250             [0.04255, 0.241, 0.4329, 0.241, 0.04255]
01251             \endcode
01252             
01253             These values are optimal in the sense that the 5x5 filter obtained by combining 
01254             this filter with the optimal 5-tap first derivative is the best possible 5x5 approximation to a 
01255             Gaussian first derivative filter. The equivalent Gaussian has sigma = 0.906.
01256  
01257             Postconditions:
01258             \code
01259             1. left()  == -2
01260             2. right() ==  2
01261             3. borderTreatment() == BORDER_TREATMENT_REFLECT
01262             4. norm() == 1.0
01263             \endcode
01264         */
01265     void initOptimalFirstDerivativeSmoothing5()
01266     {
01267         this->initExplicitly(-2, 2) = 0.04255, 0.241, 0.4329, 0.241, 0.04255;
01268         this->setBorderTreatment(BORDER_TREATMENT_REFLECT);
01269     }
01270 
01271         /**
01272             Init an optimal 5-tap smoothing filter to be used in the context of second derivative computation.
01273            This filter must be used in conjunction with the optimal 5-tap second derivative filter 
01274            (see initOptimalSecondDerivative5()), such that the derivative filter is applied along one dimension, 
01275            and the smoothing filter along the other. The filter values are 
01276             
01277             \code
01278             [0.0243, 0.23556, 0.48028, 0.23556, 0.0243]
01279             \endcode
01280             
01281             These values are optimal in the sense that the 5x5 filter obtained by combining 
01282             this filter with the optimal 5-tap second derivative is the best possible 5x5 approximation to a 
01283             Gaussian second derivative filter. The equivalent Gaussian has sigma = 0.817.
01284  
01285             Postconditions:
01286             \code
01287             1. left()  == -2
01288             2. right() ==  2
01289             3. borderTreatment() == BORDER_TREATMENT_REFLECT
01290             4. norm() == 1.0
01291             \endcode
01292         */
01293     void initOptimalSecondDerivativeSmoothing5()
01294     {
01295         this->initExplicitly(-2, 2) = 0.0243, 0.23556, 0.48028, 0.23556, 0.0243;
01296         this->setBorderTreatment(BORDER_TREATMENT_REFLECT);
01297     }
01298 
01299         /**
01300             Init a 5-tap filter as defined by Peter Burt in the context of pyramid creation.
01301             The filter values are
01302             
01303             \code
01304             [a, 0.25, 0.5-2*a, 0.25, a]
01305             \endcode
01306             
01307             The default <tt>a = 0.04785</tt> is optimal in the sense that it minimizes the difference
01308             to a true Gaussian filter (which would have sigma = 0.975). For other values of <tt>a</tt>, the scale 
01309             of the most similar Gaussian can be approximated by
01310             
01311             \code
01312             sigma = 5.1 * a + 0.731
01313             \endcode
01314  
01315             Preconditions:
01316             \code
01317             0 <= a <= 0.125
01318             \endcode
01319 
01320             Postconditions:
01321             \code
01322             1. left()  == -2
01323             2. right() ==  2
01324             3. borderTreatment() == BORDER_TREATMENT_REFLECT
01325             4. norm() == 1.0
01326             \endcode
01327         */
01328     void initBurtFilter(double a = 0.04785)
01329     {
01330         vigra_precondition(a >= 0.0 && a <= 0.125,
01331             "Kernel1D::initBurtFilter(): 0 <= a <= 0.125 required.");
01332         this->initExplicitly(-2, 2) = a, 0.25, 0.5 - 2.0*a, 0.25, a;
01333         this->setBorderTreatment(BORDER_TREATMENT_REFLECT);
01334     }
01335 
01336         /**
01337             Init as a Binomial filter. 'norm' denotes the sum of all bins
01338             of the kernel.
01339 
01340             Precondition:
01341             \code
01342             radius   >= 0
01343             \endcode
01344 
01345             Postconditions:
01346             \code
01347             1. left()  == -radius
01348             2. right() ==  radius
01349             3. borderTreatment() == BORDER_TREATMENT_REFLECT
01350             4. norm() == norm
01351             \endcode
01352         */
01353     void initBinomial(int radius, value_type norm);
01354 
01355         /** Init as a Binomial filter with norm 1.
01356          */
01357     void initBinomial(int radius)
01358     {
01359         initBinomial(radius, one());
01360     }
01361 
01362         /**
01363             Init as an Averaging filter. 'norm' denotes the sum of all bins
01364             of the kernel. The window size is (2*radius+1) * (2*radius+1)
01365 
01366             Precondition:
01367             \code
01368             radius   >= 0
01369             \endcode
01370 
01371             Postconditions:
01372             \code
01373             1. left()  == -radius
01374             2. right() ==  radius
01375             3. borderTreatment() == BORDER_TREATMENT_CLIP
01376             4. norm() == norm
01377             \endcode
01378         */
01379     void initAveraging(int radius, value_type norm);
01380 
01381         /** Init as an Averaging filter with norm 1.
01382          */
01383     void initAveraging(int radius)
01384     {
01385         initAveraging(radius, one());
01386     }
01387 
01388         /**
01389            Init as a symmetric gradient filter of the form
01390            <TT>[ 0.5 * norm, 0.0 * norm, -0.5 * norm]</TT>
01391            
01392            <b>Deprecated</b>. Use initSymmetricDifference() instead.
01393 
01394             Postconditions:
01395             \code
01396             1. left()  == -1
01397             2. right() ==  1
01398             3. borderTreatment() == BORDER_TREATMENT_REPEAT
01399             4. norm() == norm
01400             \endcode
01401         */
01402     void
01403     initSymmetricGradient(value_type norm )
01404     {
01405         initSymmetricDifference(norm);
01406         setBorderTreatment(BORDER_TREATMENT_REPEAT);
01407     }
01408 
01409         /** Init as a symmetric gradient filter with norm 1.
01410            
01411            <b>Deprecated</b>. Use initSymmetricDifference() instead.
01412          */
01413     void initSymmetricGradient()
01414     {
01415         initSymmetricGradient(one());
01416     }
01417 
01418     void
01419     initSymmetricDifference(value_type norm );
01420 
01421         /** Init as the 3-tap symmetric difference filter
01422              The filter values are
01423              
01424              \code
01425              [0.5, 0, -0.5]
01426              \endcode
01427 
01428             Postconditions:
01429             \code
01430             1. left()  == -1
01431             2. right() ==  1
01432             3. borderTreatment() == BORDER_TREATMENT_REFLECT
01433             4. norm() == 1.0
01434             \endcode
01435           */
01436     void initSymmetricDifference()
01437     {
01438         initSymmetricDifference(one());
01439     }
01440 
01441         /**
01442             Init the 3-tap second difference filter.
01443             The filter values are
01444              
01445              \code
01446              [1, -2, 1]
01447              \endcode
01448 
01449             Postconditions:
01450             \code
01451             1. left()  == -1
01452             2. right() ==  1
01453             3. borderTreatment() == BORDER_TREATMENT_REFLECT
01454             4. norm() == 1
01455             \endcode
01456         */
01457     void
01458     initSecondDifference3()
01459     {
01460         this->initExplicitly(-1, 1) = 1.0, -2.0, 1.0;
01461         this->setBorderTreatment(BORDER_TREATMENT_REFLECT);
01462     }
01463     
01464         /**
01465             Init an optimal 5-tap first derivative filter.
01466            This filter must be used in conjunction with the corresponding 5-tap smoothing filter 
01467            (see initOptimalFirstDerivativeSmoothing5()), such that the derivative filter is applied along one dimension,
01468             and the smoothing filter along the other.
01469             The filter values are 
01470             
01471             \code
01472             [0.1, 0.3, 0.0, -0.3, -0.1]
01473             \endcode
01474             
01475             These values are optimal in the sense that the 5x5 filter obtained by combining 
01476             this filter with the corresponding 5-tap smoothing filter is the best possible 5x5 approximation to a 
01477             Gaussian first derivative filter. The equivalent Gaussian has sigma = 0.906.
01478             
01479             If the filter is instead separably combined with itself, an almost optimal approximation of the
01480             mixed second Gaussian derivative at scale sigma = 0.899 results.
01481  
01482             Postconditions:
01483             \code
01484             1. left()  == -2
01485             2. right() ==  2
01486             3. borderTreatment() == BORDER_TREATMENT_REFLECT
01487             4. norm() == 1.0
01488             \endcode
01489         */
01490     void initOptimalFirstDerivative5()
01491     {
01492         this->initExplicitly(-2, 2) = 0.1, 0.3, 0.0, -0.3, -0.1;
01493         this->setBorderTreatment(BORDER_TREATMENT_REFLECT);
01494     }
01495     
01496         /**
01497             Init an optimal 5-tap second derivative filter.
01498            This filter must be used in conjunction with the corresponding 5-tap smoothing filter 
01499            (see initOptimalSecondDerivativeSmoothing5()), such that the derivative filter is applied along one dimension,
01500             and the smoothing filter along the other.
01501             The filter values are 
01502             
01503             \code
01504             [0.22075, 0.117, -0.6755, 0.117, 0.22075]
01505             \endcode
01506             
01507             These values are optimal in the sense that the 5x5 filter obtained by combining 
01508             this filter with the corresponding 5-tap smoothing filter is the best possible 5x5 approximation to a 
01509             Gaussian second derivative filter. The equivalent Gaussian has sigma = 0.817.
01510  
01511             Postconditions:
01512             \code
01513             1. left()  == -2
01514             2. right() ==  2
01515             3. borderTreatment() == BORDER_TREATMENT_REFLECT
01516             4. norm() == 1.0
01517             \endcode
01518         */
01519     void initOptimalSecondDerivative5()
01520     {
01521         this->initExplicitly(-2, 2) = 0.22075, 0.117, -0.6755, 0.117, 0.22075;
01522         this->setBorderTreatment(BORDER_TREATMENT_REFLECT);
01523     }
01524 
01525         /** Init the kernel by an explicit initializer list.
01526             The left and right boundaries of the kernel must be passed.
01527             A comma-separated initializer list is given after the assignment
01528             operator. This function is used like this:
01529 
01530             \code
01531             // define horizontal Roberts filter
01532             vigra::Kernel1D<float> roberts_gradient_x;
01533 
01534             roberts_gradient_x.initExplicitly(0, 1) = 1.0, -1.0;
01535             \endcode
01536 
01537             The norm is set to the sum of the initialzer values. If the wrong number of
01538             values is given, a run-time error results. It is, however, possible to give
01539             just one initializer. This creates an averaging filter with the given constant:
01540 
01541             \code
01542             vigra::Kernel1D<float> average5x1;
01543 
01544             average5x1.initExplicitly(-2, 2) = 1.0/5.0;
01545             \endcode
01546 
01547             Here, the norm is set to value*size().
01548 
01549             <b> Preconditions:</b>
01550 
01551             \code
01552 
01553             1. left <= 0
01554             2. right >= 0
01555             3. the number of values in the initializer list
01556                is 1 or equals the size of the kernel.
01557             \endcode
01558         */
01559     Kernel1D & initExplicitly(int left, int right)
01560     {
01561         vigra_precondition(left <= 0,
01562                      "Kernel1D::initExplicitly(): left border must be <= 0.");
01563         vigra_precondition(right >= 0,
01564                      "Kernel1D::initExplicitly(): right border must be >= 0.");
01565 
01566         right_ = right;
01567         left_ = left;
01568 
01569         kernel_.resize(right - left + 1);
01570 
01571         return *this;
01572     }
01573 
01574         /** Get iterator to center of kernel
01575 
01576             Postconditions:
01577             \code
01578 
01579             center()[left()] ... center()[right()] are valid kernel positions
01580             \endcode
01581         */
01582     iterator center()
01583     {
01584         return kernel_.begin() - left();
01585     }
01586 
01587     const_iterator center() const
01588     {
01589         return kernel_.begin() - left();
01590     }
01591 
01592         /** Access kernel value at specified location.
01593 
01594             Preconditions:
01595             \code
01596 
01597             left() <= location <= right()
01598             \endcode
01599         */
01600     reference operator[](int location)
01601     {
01602         return kernel_[location - left()];
01603     }
01604 
01605     const_reference operator[](int location) const
01606     {
01607         return kernel_[location - left()];
01608     }
01609 
01610         /** left border of kernel (inclusive), always <= 0
01611         */
01612     int left() const { return left_; }
01613 
01614         /** right border of kernel (inclusive), always >= 0
01615         */
01616     int right() const { return right_; }
01617 
01618         /** size of kernel (right() - left() + 1)
01619         */
01620     int size() const { return right_ - left_ + 1; }
01621 
01622         /** current border treatment mode
01623         */
01624     BorderTreatmentMode borderTreatment() const
01625     { return border_treatment_; }
01626 
01627         /** Set border treatment mode.
01628         */
01629     void setBorderTreatment( BorderTreatmentMode new_mode)
01630     { border_treatment_ = new_mode; }
01631 
01632         /** norm of kernel
01633         */
01634     value_type norm() const { return norm_; }
01635 
01636         /** set a new norm and normalize kernel, use the normalization formula
01637             for the given <tt>derivativeOrder</tt>.
01638         */
01639     void
01640     normalize(value_type norm, unsigned int derivativeOrder = 0, double offset = 0.0);
01641 
01642         /** normalize kernel to norm 1.
01643         */
01644     void
01645     normalize()
01646     {
01647         normalize(one());
01648     }
01649 
01650         /** get a const accessor
01651         */
01652     ConstAccessor accessor() const { return ConstAccessor(); }
01653 
01654         /** get an accessor
01655         */
01656     Accessor accessor() { return Accessor(); }
01657 
01658   private:
01659     InternalVector kernel_;
01660     int left_, right_;
01661     BorderTreatmentMode border_treatment_;
01662     value_type norm_;
01663 };
01664 
01665 template <class ARITHTYPE>
01666 void Kernel1D<ARITHTYPE>::normalize(value_type norm,
01667                           unsigned int derivativeOrder,
01668                           double offset)
01669 {
01670     typedef typename NumericTraits<value_type>::RealPromote TmpType;
01671 
01672     // find kernel sum
01673     Iterator k = kernel_.begin();
01674     TmpType sum = NumericTraits<TmpType>::zero();
01675 
01676     if(derivativeOrder == 0)
01677     {
01678         for(; k < kernel_.end(); ++k)
01679         {
01680             sum += *k;
01681         }
01682     }
01683     else
01684     {
01685         unsigned int faculty = 1;
01686         for(unsigned int i = 2; i <= derivativeOrder; ++i)
01687             faculty *= i;
01688         for(double x = left() + offset; k < kernel_.end(); ++x, ++k)
01689         {
01690             sum += *k * VIGRA_CSTD::pow(-x, int(derivativeOrder)) / faculty;
01691         }
01692     }
01693 
01694     vigra_precondition(sum != NumericTraits<value_type>::zero(),
01695                     "Kernel1D<ARITHTYPE>::normalize(): "
01696                     "Cannot normalize a kernel with sum = 0");
01697     // normalize
01698     sum = norm / sum;
01699     k = kernel_.begin();
01700     for(; k != kernel_.end(); ++k)
01701     {
01702         *k = *k * sum;
01703     }
01704 
01705     norm_ = norm;
01706 }
01707 
01708 /***********************************************************************/
01709 
01710 template <class ARITHTYPE>
01711 void Kernel1D<ARITHTYPE>::initGaussian(double std_dev,
01712                                        value_type norm)
01713 {
01714     vigra_precondition(std_dev >= 0.0,
01715               "Kernel1D::initGaussian(): Standard deviation must be >= 0.");
01716 
01717     if(std_dev > 0.0)
01718     {
01719         Gaussian<ARITHTYPE> gauss(std_dev);
01720 
01721         // first calculate required kernel sizes
01722         int radius = (int)(3.0 * std_dev + 0.5);
01723         if(radius == 0)
01724             radius = 1;
01725 
01726         // allocate the kernel
01727         kernel_.erase(kernel_.begin(), kernel_.end());
01728         kernel_.reserve(radius*2+1);
01729 
01730         for(ARITHTYPE x = -radius; x <= radius; ++x)
01731         {
01732             kernel_.push_back(gauss(x));
01733         }
01734         left_ = -radius;
01735         right_ = radius;
01736     }
01737     else
01738     {
01739         kernel_.erase(kernel_.begin(), kernel_.end());
01740         kernel_.push_back(1.0);
01741         left_ = 0;
01742         right_ = 0;
01743     }
01744 
01745     if(norm != 0.0)
01746         normalize(norm);
01747     else
01748         norm_ = 1.0;
01749 
01750     // best border treatment for Gaussians is BORDER_TREATMENT_REFLECT
01751     border_treatment_ = BORDER_TREATMENT_REFLECT;
01752 }
01753 
01754 /***********************************************************************/
01755 
01756 template <class ARITHTYPE>
01757 void Kernel1D<ARITHTYPE>::initDiscreteGaussian(double std_dev,
01758                                        value_type norm)
01759 {
01760     vigra_precondition(std_dev >= 0.0,
01761               "Kernel1D::initDiscreteGaussian(): Standard deviation must be >= 0.");
01762 
01763     if(std_dev > 0.0)
01764     {
01765         // first calculate required kernel sizes
01766         int radius = (int)(3.0*std_dev + 0.5);
01767         if(radius == 0)
01768             radius = 1;
01769 
01770         double f = 2.0 / std_dev / std_dev;
01771 
01772         // allocate the working array
01773         int maxIndex = (int)(2.0 * (radius + 5.0 * VIGRA_CSTD::sqrt((double)radius)) + 0.5);
01774         InternalVector warray(maxIndex+1);
01775         warray[maxIndex] = 0.0;
01776         warray[maxIndex-1] = 1.0;
01777 
01778         for(int i = maxIndex-2; i >= radius; --i)
01779         {
01780             warray[i] = warray[i+2] + f * (i+1) * warray[i+1];
01781             if(warray[i] > 1.0e40)
01782             {
01783                 warray[i+1] /= warray[i];
01784                 warray[i] = 1.0;
01785             }
01786         }
01787 
01788         // the following rescaling ensures that the numbers stay in a sensible range
01789         // during the rest of the iteration, so no other rescaling is needed
01790         double er = VIGRA_CSTD::exp(-radius*radius / (2.0*std_dev*std_dev));
01791         warray[radius+1] = er * warray[radius+1] / warray[radius];
01792         warray[radius] = er;
01793 
01794         for(int i = radius-1; i >= 0; --i)
01795         {
01796             warray[i] = warray[i+2] + f * (i+1) * warray[i+1];
01797             er += warray[i];
01798         }
01799 
01800         double scale = norm / (2*er - warray[0]);
01801 
01802         initExplicitly(-radius, radius);
01803         iterator c = center();
01804 
01805         for(int i=0; i<=radius; ++i)
01806         {
01807             c[i] = c[-i] = warray[i] * scale;
01808         }
01809     }
01810     else
01811     {
01812         kernel_.erase(kernel_.begin(), kernel_.end());
01813         kernel_.push_back(norm);
01814         left_ = 0;
01815         right_ = 0;
01816     }
01817 
01818     norm_ = norm;
01819 
01820     // best border treatment for Gaussians is BORDER_TREATMENT_REFLECT
01821     border_treatment_ = BORDER_TREATMENT_REFLECT;
01822 }
01823 
01824 /***********************************************************************/
01825 
01826 template <class ARITHTYPE>
01827 void
01828 Kernel1D<ARITHTYPE>::initGaussianDerivative(double std_dev,
01829                     int order,
01830                     value_type norm)
01831 {
01832     vigra_precondition(order >= 0,
01833               "Kernel1D::initGaussianDerivative(): Order must be >= 0.");
01834 
01835     if(order == 0)
01836     {
01837         initGaussian(std_dev, norm);
01838         return;
01839     }
01840 
01841     vigra_precondition(std_dev > 0.0,
01842               "Kernel1D::initGaussianDerivative(): "
01843               "Standard deviation must be > 0.");
01844 
01845     Gaussian<ARITHTYPE> gauss(std_dev, order);
01846 
01847     // first calculate required kernel sizes
01848     int radius = (int)(3.0 * std_dev + 0.5 * order + 0.5);
01849     if(radius == 0)
01850         radius = 1;
01851 
01852     // allocate the kernels
01853     kernel_.clear();
01854     kernel_.reserve(radius*2+1);
01855 
01856     // fill the kernel and calculate the DC component
01857     // introduced by truncation of the Gaussian
01858     ARITHTYPE dc = 0.0;
01859     for(ARITHTYPE x = -radius; x <= radius; ++x)
01860     {
01861         kernel_.push_back(gauss(x));
01862         dc += kernel_[kernel_.size()-1];
01863     }
01864     dc /= (2.0*radius + 1.0);
01865 
01866     // remove DC, but only if kernel correction is permitted by a non-zero
01867     // value for norm
01868     if(norm != 0.0)
01869     {
01870         for(unsigned int i=0; i < kernel_.size(); ++i)
01871         {
01872             kernel_[i] -= dc;
01873         }
01874     }
01875 
01876     left_ = -radius;
01877     right_ = radius;
01878 
01879     if(norm != 0.0)
01880         normalize(norm, order);
01881     else
01882         norm_ = 1.0;
01883 
01884     // best border treatment for Gaussian derivatives is
01885     // BORDER_TREATMENT_REFLECT
01886     border_treatment_ = BORDER_TREATMENT_REFLECT;
01887 }
01888 
01889 /***********************************************************************/
01890 
01891 template <class ARITHTYPE>
01892 void
01893 Kernel1D<ARITHTYPE>::initBinomial(int radius,
01894                                   value_type norm)
01895 {
01896     vigra_precondition(radius > 0,
01897               "Kernel1D::initBinomial(): Radius must be > 0.");
01898 
01899     // allocate the kernel
01900     InternalVector kernel(radius*2+1);
01901 
01902     int i,j;
01903     for(i=0; i<radius*2+1; ++i) kernel[i] = 0;
01904 
01905     // fill kernel
01906     typename InternalVector::iterator x = kernel.begin() + radius;
01907     x[radius] = 1.0;
01908 
01909     for(j=radius-1; j>=-radius; --j)
01910     {
01911         for(i=j; i<radius; ++i)
01912         {
01913             x[i] = (x[i] + x[i+1]) / 2.0;
01914         }
01915         x[radius] /= 2.0;
01916     }
01917 
01918     // normalize
01919     kernel_.erase(kernel_.begin(), kernel_.end());
01920     kernel_.reserve(radius*2+1);
01921 
01922     for(i=0; i<=radius*2+1; ++i)
01923     {
01924         kernel_.push_back(kernel[i] * norm);
01925     }
01926 
01927     left_ = -radius;
01928     right_ = radius;
01929     norm_ = norm;
01930 
01931     // best border treatment for Binomial is BORDER_TREATMENT_REFLECT
01932     border_treatment_ = BORDER_TREATMENT_REFLECT;
01933 }
01934 
01935 /***********************************************************************/
01936 
01937 template <class ARITHTYPE>
01938 void Kernel1D<ARITHTYPE>::initAveraging(int radius,
01939                                         value_type norm)
01940 {
01941     vigra_precondition(radius > 0,
01942               "Kernel1D::initAveraging(): Radius must be > 0.");
01943 
01944     // calculate scaling
01945     double scale = 1.0 / (radius * 2 + 1);
01946 
01947     // normalize
01948     kernel_.erase(kernel_.begin(), kernel_.end());
01949     kernel_.reserve(radius*2+1);
01950 
01951     for(int i=0; i<=radius*2+1; ++i)
01952     {
01953         kernel_.push_back(scale * norm);
01954     }
01955 
01956     left_ = -radius;
01957     right_ = radius;
01958     norm_ = norm;
01959 
01960     // best border treatment for Averaging is BORDER_TREATMENT_CLIP
01961     border_treatment_ = BORDER_TREATMENT_CLIP;
01962 }
01963 
01964 /***********************************************************************/
01965 
01966 template <class ARITHTYPE>
01967 void
01968 Kernel1D<ARITHTYPE>::initSymmetricDifference(value_type norm)
01969 {
01970     kernel_.erase(kernel_.begin(), kernel_.end());
01971     kernel_.reserve(3);
01972 
01973     kernel_.push_back(0.5 * norm);
01974     kernel_.push_back(0.0 * norm);
01975     kernel_.push_back(-0.5 * norm);
01976 
01977     left_ = -1;
01978     right_ = 1;
01979     norm_ = norm;
01980 
01981     // best border treatment for symmetric difference is
01982     // BORDER_TREATMENT_REFLECT
01983     border_treatment_ = BORDER_TREATMENT_REFLECT;
01984 }
01985 
01986 /**************************************************************/
01987 /*                                                            */
01988 /*         Argument object factories for Kernel1D             */
01989 /*                                                            */
01990 /*     (documentation: see vigra/convolution.hxx)             */
01991 /*                                                            */
01992 /**************************************************************/
01993 
01994 template <class KernelIterator, class KernelAccessor>
01995 inline
01996 tuple5<KernelIterator, KernelAccessor, int, int, BorderTreatmentMode>
01997 kernel1d(KernelIterator ik, KernelAccessor ka,
01998        int kleft, int kright, BorderTreatmentMode border)
01999 {
02000     return
02001       tuple5<KernelIterator, KernelAccessor, int, int, BorderTreatmentMode>(
02002                                                 ik, ka, kleft, kright, border);
02003 }
02004 
02005 template <class T>
02006 inline
02007 tuple5<typename Kernel1D<T>::const_iterator, typename Kernel1D<T>::ConstAccessor,
02008        int, int, BorderTreatmentMode>
02009 kernel1d(Kernel1D<T> const & k)
02010 
02011 {
02012     return
02013         tuple5<typename Kernel1D<T>::const_iterator, typename Kernel1D<T>::ConstAccessor,
02014                int, int, BorderTreatmentMode>(
02015                                      k.center(),
02016                                      k.accessor(),
02017                                      k.left(), k.right(),
02018                                      k.borderTreatment());
02019 }
02020 
02021 template <class T>
02022 inline
02023 tuple5<typename Kernel1D<T>::const_iterator, typename Kernel1D<T>::ConstAccessor,
02024        int, int, BorderTreatmentMode>
02025 kernel1d(Kernel1D<T> const & k, BorderTreatmentMode border)
02026 
02027 {
02028     return
02029         tuple5<typename Kernel1D<T>::const_iterator, typename Kernel1D<T>::ConstAccessor,
02030                int, int, BorderTreatmentMode>(
02031                                      k.center(),
02032                                      k.accessor(),
02033                                      k.left(), k.right(),
02034                                      border);
02035 }
02036 
02037 
02038 } // namespace vigra
02039 
02040 #endif // VIGRA_SEPARABLECONVOLUTION_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)