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

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

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
vigra 1.7.0 (Thu Aug 25 2011)