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

© Ullrich Köthe (koethe@informatik.uni-hamburg.de) Cognitive Systems Group, University of Hamburg, Germany | html generated using doxygen and Python VIGRA 1.3.2 (27 Jan 2005)