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