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

vigra/separableconvolution.hxx VIGRA

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