[ 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.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)