[ 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.1 (3 Dec 2010)