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

details vigra/separableconvolution.hxx VIGRA

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