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

| vigra/stdconvolution.hxx |
--- | --- | ---

    00001 /************************************************************************/
    00002 /*                                                                      */
    00003 /*               Copyright 1998-2002 by Ullrich Koethe                  */
    00004 /*       Cognitive Systems Group, University of Hamburg, Germany        */
    00005 /*                                                                      */
    00006 /*    This file is part of the VIGRA computer vision library.           */
    00007 /*    ( Version 1.3.2, Jan 27 2005 )                                    */
    00008 /*    You may use, modify, and distribute this software according       */
    00009 /*    to the terms stated in the LICENSE file included in               */
    00010 /*    the VIGRA distribution.                                           */
    00011 /*                                                                      */
    00012 /*    The VIGRA Website is                                              */
    00013 /*        http://kogs-www.informatik.uni-hamburg.de/~koethe/vigra/      */
    00014 /*    Please direct questions, bug reports, and contributions to        */
    00015 /*        koethe@informatik.uni-hamburg.de                              */
    00016 /*                                                                      */
    00017 /*  THIS SOFTWARE IS PROVIDED AS IS AND WITHOUT ANY EXPRESS OR          */
    00018 /*  IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED      */
    00019 /*  WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. */
    00020 /*                                                                      */
    00021 /************************************************************************/
    00022
    00023
    00024 #ifndef VIGRA_STDCONVOLUTION_HXX
    00025 #define VIGRA_STDCONVOLUTION_HXX
    00026
    00027 #include <cmath>
    00028 #include "vigra/stdimage.hxx"
    00029 #include "vigra/bordertreatment.hxx"
    00030 #include "vigra/separableconvolution.hxx"
    00031 #include "vigra/utilities.hxx"
    00032
    00033 namespace vigra {
    00034
    00035 template <class SrcIterator, class SrcAccessor,
    00036           class DestIterator, class DestAccessor,
    00037           class KernelIterator, class KernelAccessor,
    00038           class KSumType>
    00039 void internalPixelEvaluationByClip(int x, int y, int w, int h, SrcIterator xs,
    00040                                    SrcAccessor src_acc, DestIterator xd, DestAccessor dest_acc,
    00041                                    KernelIterator ki, Diff2D kul, Diff2D klr, KernelAccessor ak,
    00042                                    KSumType norm)
    00043 {
    00044     typedef typename
    00045         NumericTraits<typename SrcAccessor::value_type>::RealPromote SumType;
    00046     typedef
    00047         NumericTraits<typename DestAccessor::value_type> DestTraits;
    00048
    00049     // calculate width and height of the kernel
    00050     int kernel_width = klr.x - kul.x + 1;
    00051     int kernel_height = klr.y - kul.y + 1;
    00052
    00053     SumType sum = NumericTraits<SumType>::zero();
    00054     int xx, yy;
    00055     int x0, y0, x1, y1;
    00056
    00057     y0 = (y<klr.y) ?  -y : -klr.y;
    00058     y1 = (h-y-1<-kul.y) ? h-y-1 : -kul.y;
    00059
    00060     x0 = (x<klr.x) ? -x : -klr.x;
    00061     x1 = (w-x-1<-kul.x) ? w-x-1 : -kul.x;
    00062
    00063     SrcIterator yys = xs + Diff2D(x0, y0);
    00064     KernelIterator yk  = ki - Diff2D(x0, y0);
    00065
    00066     KSumType ksum = NumericTraits<KSumType>::zero();
    00067     kernel_width = x1 - x0 + 1;
    00068     kernel_height = y1 - y0 + 1;
    00069
    00070     //es wird zuerst abgeschnitten und dann gespigelt!
    00071
    00072     for(yy=0; yy<kernel_height; ++yy, ++yys.y, --yk.y)
    00073     {
    00074         SrcIterator xxs = yys;
    00075         KernelIterator xk  = yk;
    00076
    00077         for(xx=0; xx<kernel_width; ++xx, ++xxs.x, --xk.x)
    00078         {
    00079             sum += ak(xk) * src_acc(xxs);
    00080             ksum += ak(xk);
    00081         }
    00082     }
    00083
    00084     //                      store average in destination pixel
    00085     dest_acc.set(DestTraits::fromRealPromote((norm / ksum) * sum), xd);
    00086
    00087 }
    00088
    00089
    00090 #if 0
    00091
    00092 template <class SrcIterator, class SrcAccessor,
    00093           class DestIterator, class DestAccessor,
    00094           class KernelIterator, class KernelAccessor>
    00095 void internalPixelEvaluationByWrapReflectRepeat(int x, int y, int src_width, int src_height, SrcIterator xs,
    00096                                                 SrcAccessor src_acc, DestIterator xd, DestAccessor dest_acc,
    00097                                                 KernelIterator ki, Diff2D kul, Diff2D klr, KernelAccessor ak,
    00098                                                 BorderTreatmentMode border)
    00099 {
    00100
    00101     typedef typename
    00102         NumericTraits<typename SrcAccessor::value_type>::RealPromote SumType;
    00103     typedef
    00104         NumericTraits<typename DestAccessor::value_type> DestTraits;
    00105
    00106     SumType sum = NumericTraits<SumType>::zero();
    00107
    00108     SrcIterator src_ul = xs - Diff2D(x, y);
    00109     SrcIterator src_lr = src_ul + Diff2D(src_width, src_height);
    00110
    00111     SrcIterator yys = xs;
    00112     KernelIterator yk  = ki;
    00113
    00114     // calculate width and height of the kernel
    00115     int kernel_width = klr.x - kul.x + 1;
    00116     int kernel_height = klr.y - kul.y + 1;
    00117
    00118     //Zeigt an wo der Kernel ’ber die Grenzen hinausgeht
    00119     bool top_to_much = (y<klr.y) ? true : false;
    00120     bool down_to_much = (src_height-y-1<-kul.y)? true : false;
    00121     bool left_to_much = (x<klr.x)? true : false;
    00122     bool right_to_much = (src_width-x-1<-kul.x)? true : false;
    00123
    00124     //Die Richtung x und y !!!
    00125     //in der bei der Iteration ’ber das aktuelle Bereich im Bild
    00126     //iteriert wird. Also wenn von ur->ll dann (-1, +1) und wenn lr->ul
    00127     //dann (-1, -1).
    00128     Diff2D way_increment;
    00129
    00130     /* iteriert wird immer aus dem g’ltigen in den ung’ltigen
    00131        Bereich! dieser Tupel setzt sich wie folgt zusammen:
    00132        1. Wird bei der Iteration in X-Richtung ung’ltiger Bereich
    00133        erreicht so wird mit border_increment.first gesprungen und
    00134        mit border_increment.third weiter iteriert.
    00135        2. Wird bei der Iteration in Y-Richtung ung’ltiger Bereich
    00136        erreicht so wird mit border_increment.second gesprungen und
    00137        mit border_increment.fourth weiter iteriert.
    00138     */
    00139     tuple4<int, int, int, int> border_increment;
    00140     if (border == BORDER_TREATMENT_REPEAT){
    00141         border_increment = tuple4<int, int, int, int>(1, 1, 0, 0);
    00142     }else if (border == BORDER_TREATMENT_REFLECT){
    00143         border_increment = tuple4<int, int, int, int>(2, 2, -1, -1);
    00144     }else{ // BORDER_TREATMENT_WRAP
    00145         border_increment = tuple4<int, int, int, int>(src_width, src_height, 1, 1);
    00146     }
    00147
    00148     pair<int, int> valid_step_count;
    00149
    00150     if(left_to_much && !top_to_much && !down_to_much)
    00151     {
    00152         yys += klr;
    00153         yk += kul;
    00154         way_increment = Diff2D(-1, -1);
    00155         border_increment.third = -border_increment.third;
    00156         border_increment.fourth = -border_increment.fourth;
    00157         valid_step_count = std::make_pair((yys - src_ul).x + 1, kernel_height);
    00158     }
    00159     else if(top_to_much && !left_to_much && !right_to_much)
    00160     {
    00161         yys += klr;
    00162         yk += kul;
    00163         way_increment = Diff2D(-1, -1);
    00164         border_increment.third = -border_increment.third;
    00165         border_increment.fourth = -border_increment.fourth;
    00166         valid_step_count = std::make_pair(kernel_width, (yys - src_ul).y + 1);
    00167     }
    00168     else if(right_to_much && !top_to_much && !down_to_much)
    00169     {
    00170         yys += kul;
    00171         yk += klr;
    00172         way_increment = Diff2D(1, 1);
    00173         border_increment.first = -border_increment.first;
    00174         border_increment.second = -border_increment.second;
    00175         valid_step_count = std::make_pair((src_lr - yys).x, kernel_height);
    00176     }
    00177     else if(down_to_much && !left_to_much && !right_to_much)
    00178     {
    00179         yys += kul;
    00180         yk += klr;
    00181         way_increment = Diff2D(1, 1);
    00182         border_increment.first = -border_increment.first;
    00183         border_increment.second = -border_increment.second;
    00184         valid_step_count = std::make_pair(kernel_width, (src_lr - yys).y);
    00185     }
    00186     else if(down_to_much && left_to_much)
    00187     {
    00188         yys += kul + Diff2D(kernel_width - 1, 0);
    00189         yk += kul + Diff2D(0, kernel_height - 1);
    00190         way_increment = Diff2D(-1, 1);
    00191         border_increment.second = -border_increment.second;
    00192         border_increment.third = -border_increment.third;
    00193         valid_step_count = std::make_pair((yys - src_ul).x + 1, (src_lr - yys).y);
    00194     }
    00195     else if(down_to_much && right_to_much)
    00196     {
    00197         yys += kul;
    00198         yk += klr;
    00199         way_increment = Diff2D(1, 1);
    00200         border_increment.first = -border_increment.first;
    00201         border_increment.second = -border_increment.second;
    00202         valid_step_count = std::make_pair((src_lr - yys).x, (src_lr - yys).y);
    00203     }
    00204     else if(top_to_much && left_to_much)
    00205     {
    00206         yys += klr;
    00207         yk += kul;
    00208         way_increment = Diff2D(-1, -1);
    00209         border_increment.third = -border_increment.third;
    00210         border_increment.fourth = -border_increment.fourth;
    00211         valid_step_count = std::make_pair((yys - src_ul).x + 1, (yys - src_ul).y + 1);
    00212     }
    00213     else
    00214     { //top_to_much && right_to_much
    00215         yys += kul + Diff2D(0, kernel_height - 1);
    00216         yk += kul + Diff2D(kernel_width - 1, 0);
    00217         way_increment = Diff2D(1, -1);
    00218         border_increment.first = -border_increment.first;
    00219         border_increment.fourth = -border_increment.fourth;
    00220         valid_step_count = std::make_pair((src_lr - yys).x, (yys - src_ul).y + 1);
    00221     }
    00222
    00223     int yy = 0, xx;
    00224
    00225     //laeuft den zul„ssigen Bereich in y-Richtung durch
    00226     for(; yy < valid_step_count.second; ++yy, yys.y += way_increment.y, yk.y -= way_increment.y )
    00227     {
    00228         SrcIterator xxs = yys;
    00229         KernelIterator xk  = yk;
    00230
    00231         //laeuft den zul„ssigen Bereich in x-Richtung durch
    00232         for(xx = 0; xx < valid_step_count.first; ++xx, xxs.x += way_increment.x, xk.x -= way_increment.x)
    00233         {
    00234             sum += ak(xk) * src_acc(xxs);
    00235         }
    00236
    00237         //N„chstes ++xxs.x wuerde in unzul„ssigen Bereich
    00238         //bringen => Sprung in zulaessigen Bereich
    00239         xxs.x += border_increment.first;
    00240
    00241         for( ; xx < kernel_width; ++xx, xxs.x += border_increment.third, xk.x -= way_increment.x )
    00242         {
    00243             sum += ak(xk) * src_acc(xxs);
    00244         }
    00245     }
    00246
    00247     //N„chstes ++yys.y wuerde in unzul„ssigen Bereich
    00248     //bringen => Sprung in zulaessigen Bereich
    00249     yys.y += border_increment.second;
    00250
    00251     for( ; yy < kernel_height; ++yy, yys.y += border_increment.third, yk.y -= way_increment.y)
    00252     {
    00253         SrcIterator xxs = yys;
    00254         KernelIterator xk  = yk;
    00255
    00256         for(xx=0; xx < valid_step_count.first; ++xx, xxs.x += way_increment.x, xk.x -= way_increment.x)
    00257         {
    00258             sum += ak(xk) * src_acc(xxs);
    00259         }
    00260
    00261         //Sprung in den zulaessigen Bereich
    00262         xxs.x += border_increment.first;
    00263
    00264         for( ; xx < kernel_width; ++xx, xxs.x += border_increment.third, xk.x -= way_increment.x )
    00265         {
    00266             sum += ak(xk) * src_acc(xxs);
    00267         }
    00268     }
    00269
    00270     // store average in destination pixel
    00271     dest_acc.set(DestTraits::fromRealPromote(sum), xd);
    00272
    00273 }// end of internalPixelEvaluationByWrapReflectRepeat
    00274 #endif /* #if 0 */
    00275
    00276
    00277 template <class SrcIterator, class SrcAccessor,
    00278           class KernelIterator, class KernelAccessor,
    00279           class SumType>
    00280 void
    00281 internalPixelEvaluationByWrapReflectRepeat(SrcIterator xs, SrcAccessor src_acc,
    00282     KernelIterator xk, KernelAccessor ak,
    00283     int left, int right, int kleft, int kright,
    00284     int borderskipx, int borderinc, SumType & sum)
    00285 {
    00286     SrcIterator xxs = xs + left;
    00287     KernelIterator xxk  = xk - left;
    00288
    00289     for(int xx = left; xx <= right; ++xx, ++xxs, --xxk)
    00290     {
    00291         sum += ak(xxk) * src_acc(xxs);
    00292     }
    00293
    00294     xxs = xs + left - borderskipx;
    00295     xxk = xk - left + 1;
    00296     for(int xx = left - 1; xx >= -kright; --xx, xxs -= borderinc, ++xxk)
    00297     {
    00298         sum += ak(xxk) * src_acc(xxs);
    00299     }
    00300
    00301     xxs = xs + right + borderskipx;
    00302     xxk = xk - right - 1;
    00303     for(int xx = right + 1; xx <= -kleft; ++xx, xxs += borderinc, --xxk)
    00304     {
    00305         sum += ak(xxk) * src_acc(xxs);
    00306     }
    00307 }
    00308
    00309
    00310 /** \addtogroup StandardConvolution Two-dimensional convolution functions
    00311
    00312 Perform 2D non-separable convolution, with and without ROI mask.
    00313
    00314 These generic convolution functions implement
    00315 the standard 2D convolution operation for images that fit
    00316 into the required interface. Arbitrary ROI's are supported
    00317 by the mask version of the algorithm.
    00318 The functions need a suitable 2D kernel to operate.
    00319 */
    00320 //@{
    00321
    00322 /** \brief Performs a 2 dimensional convolution of the source image using the given
    00323     kernel.
    00324
    00325     The KernelIterator must point to the center of the kernel, and
    00326     the kernel's size is given by its upper left (x and y of distance <= 0) and
    00327     lower right (distance >= 0) corners. The image must always be larger than the
    00328     kernel. At those positions where the kernel does not completely fit
    00329     into the image, the specified \ref BorderTreatmentMode is
    00330     applied. You can choice between following BorderTreatmentModes:
    00331     <ul>
    00332     <li>BORDER_TREATMENT_CLIP</li>
    00333     <li>BORDER_TREATMENT_AVOID</li>
    00334     <li>BORDER_TREATMENT_WRAP</li>
    00335     <li>BORDER_TREATMENT_REFLECT</li>
    00336     <li>BORDER_TREATMENT_REPEAT</li>
    00337     </ul><br>
    00338     The images's pixel type (SrcAccessor::value_type) must be a
    00339     linear space over the kernel's value_type (KernelAccessor::value_type),
    00340     i.e. addition of source values, multiplication with kernel values,
    00341     and NumericTraits must be defined.
    00342     The kernel's value_type must be an algebraic field,
    00343     i.e. the arithmetic operations (+, -, *, /) and NumericTraits must
    00344     be defined.
    00345
    00346     <b> Declarations:</b>
    00347
    00348     pass arguments explicitly:
    00349     \code
    00350     namespace vigra {
    00351         template <class SrcIterator, class SrcAccessor,
    00352                   class DestIterator, class DestAccessor,
    00353                   class KernelIterator, class KernelAccessor>
    00354         void convolveImage(SrcIterator src_ul, SrcIterator src_lr, SrcAccessor src_acc,
    00355                            DestIterator dest_ul, DestAccessor dest_acc,
    00356                            KernelIterator ki, KernelAccessor ak,
    00357                            Diff2D kul, Diff2D klr, BorderTreatmentMode border);
    00358     }
    00359     \endcode
    00360
    00361
    00362     use argument objects in conjunction with \ref ArgumentObjectFactories:
    00363     \code
    00364     namespace vigra {
    00365         template <class SrcIterator, class SrcAccessor,
    00366                   class DestIterator, class DestAccessor,
    00367                   class KernelIterator, class KernelAccessor>
    00368         void convolveImage(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    00369                            pair<DestIterator, DestAccessor> dest,
    00370                            tuple5<KernelIterator, KernelAccessor, Diff2D, Diff2D,
    00371                            BorderTreatmentMode> kernel);
    00372     }
    00373     \endcode
    00374
    00375     <b> Usage:</b>
    00376
    00377     <b>\#include</b> "<a href="stdconvolution_8hxx-source.html">vigra/stdconvolution.hxx</a>"<br>
    00378     Namespace: vigra
    00379
    00380
    00381     \code
    00382     vigra::FImage src(w,h), dest(w,h);
    00383     ...
    00384
    00385     // define horizontal Sobel filter
    00386     vigra::Kernel2D<float> sobel;
    00387
    00388     sobel.initExplicitly(Diff2D(-1,-1), Diff2D(1,1)) =  // upper left and lower right
    00389                          0.125, 0.0, -0.125,
    00390                          0.25,  0.0, -0.25,
    00391                          0.125, 0.0, -0.125;
    00392
    00393     vigra::convolveImage(srcImageRange(src), destImage(dest), kernel2d(sobel));
    00394     \endcode
    00395
    00396     <b> Required Interface:</b>
    00397
    00398     \code
    00399     ImageIterator src_ul, src_lr;
    00400     ImageIterator dest_ul;
    00401     ImageIterator ik;
    00402
    00403     SrcAccessor src_accessor;
    00404     DestAccessor dest_accessor;
    00405     KernelAccessor kernel_accessor;
    00406
    00407     NumericTraits<SrcAccessor::value_type>::RealPromote s = src_accessor(src_ul);
    00408
    00409     s = s + s;
    00410     s = kernel_accessor(ik) * s;
    00411     s -= s;
    00412
    00413     dest_accessor.set(
    00414     NumericTraits<DestAccessor::value_type>::fromRealPromote(s), dest_ul);
    00415
    00416     NumericTraits<KernelAccessor::value_type>::RealPromote k = kernel_accessor(ik);
    00417
    00418     k += k;
    00419     k -= k;
    00420     k = k / k;
    00421
    00422     \endcode
    00423
    00424     <b> Preconditions:</b>
    00425
    00426     \code
    00427     kul.x <= 0
    00428     kul.y <= 0
    00429     klr.x >= 0
    00430     klr.y >= 0
    00431     src_lr.x - src_ul.x >= klr.x + kul.x + 1
    00432     src_lr.y - src_ul.y >= klr.y + kul.y + 1
    00433     \endcode
    00434
    00435     If border == BORDER_TREATMENT_CLIP: Sum of kernel elements must be
    00436     != 0.
    00437
    00438 */
    00439 template <class SrcIterator, class SrcAccessor,
    00440           class DestIterator, class DestAccessor,
    00441           class KernelIterator, class KernelAccessor>
    00442 void convolveImage(SrcIterator src_ul, SrcIterator src_lr, SrcAccessor src_acc,
    00443                    DestIterator dest_ul, DestAccessor dest_acc,
    00444                    KernelIterator ki, KernelAccessor ak,
    00445                    Diff2D kul, Diff2D klr, BorderTreatmentMode border)
    00446 {
    00447     vigra_precondition((border == BORDER_TREATMENT_CLIP    ||
    00448                         border == BORDER_TREATMENT_AVOID   ||
    00449                         border == BORDER_TREATMENT_REFLECT ||
    00450                         border == BORDER_TREATMENT_REPEAT  ||
    00451                         border == BORDER_TREATMENT_WRAP),
    00452                        "convolveImage():\n"
    00453                        "  Border treatment must be one of follow treatments:\n"
    00454                        "  - BORDER_TREATMENT_CLIP\n"
    00455                        "  - BORDER_TREATMENT_AVOID\n"
    00456                        "  - BORDER_TREATMENT_REFLECT\n"
    00457                        "  - BORDER_TREATMENT_REPEAT\n"
    00458                        "  - BORDER_TREATMENT_WRAP\n");
    00459
    00460     vigra_precondition(kul.x <= 0 && kul.y <= 0,
    00461                        "convolveImage(): coordinates of "
    00462                        "kernel's upper left must be <= 0.");
    00463     vigra_precondition(klr.x >= 0 && klr.y >= 0,
    00464                        "convolveImage(): coordinates of "
    00465                        "kernel's lower right must be >= 0.");
    00466
    00467     // use traits to determine SumType as to prevent possible overflow
    00468     typedef typename
    00469         NumericTraits<typename SrcAccessor::value_type>::RealPromote SumType;
    00470     typedef typename
    00471         NumericTraits<typename KernelAccessor::value_type>::RealPromote KernelSumType;
    00472     typedef
    00473         NumericTraits<typename DestAccessor::value_type> DestTraits;
    00474
    00475     // calculate width and height of the image
    00476     int w = src_lr.x - src_ul.x;
    00477     int h = src_lr.y - src_ul.y;
    00478
    00479     // calculate width and height of the kernel
    00480     int kernel_width = klr.x - kul.x + 1;
    00481     int kernel_height = klr.y - kul.y + 1;
    00482
    00483     vigra_precondition(w >= kernel_width && h >= kernel_height,
    00484                        "convolveImage(): kernel larger than image.");
    00485
    00486     int x,y;
    00487
    00488     KernelSumType norm = NumericTraits<KernelSumType>::zero();
    00489     if(border == BORDER_TREATMENT_CLIP)
    00490     {
    00491         // caluclate the sum of the kernel elements for renormalization
    00492         KernelIterator yk  = ki + klr;
    00493
    00494         //Die Summe der Punkte im Kernel wird ermittelt (= norm)
    00495         for(y=0; y<kernel_height; ++y, --yk.y)
    00496         {
    00497             KernelIterator xk  = yk;
    00498             for(x=0; x<kernel_width; ++x, --xk.x)
    00499             {
    00500                 norm += ak(xk);
    00501             }
    00502         }
    00503         vigra_precondition(norm != NumericTraits<KernelSumType>::zero(),
    00504             "convolveImage(): Cannot use BORDER_TREATMENT_CLIP with a DC-free kernel");
    00505     }
    00506
    00507     // create iterators for the interior part of the image (where the kernel always fits into the image)
    00508     DestIterator yd = dest_ul + Diff2D(klr.x, klr.y);
    00509     SrcIterator ys = src_ul + Diff2D(klr.x, klr.y);
    00510     SrcIterator send = src_lr + Diff2D(kul.x, kul.y);
    00511
    00512     // iterate over the interior part
    00513     for(; ys.y < send.y; ++ys.y, ++yd.y)
    00514     {
    00515         // create x iterators
    00516         DestIterator xd(yd);
    00517         SrcIterator xs(ys);
    00518
    00519         for(; xs.x < send.x; ++x, ++xs.x, ++xd.x)
    00520         {
    00521             // init the sum
    00522             SumType sum = NumericTraits<SumType>::zero();
    00523
    00524             SrcIterator yys = xs - klr;
    00525             SrcIterator yyend = xs - kul;
    00526             KernelIterator yk  = ki + klr;
    00527
    00528             for(; yys.y <= yyend.y; ++yys.y, --yk.y)
    00529             {
    00530                 typename SrcIterator::row_iterator xxs = yys.rowIterator();
    00531                 typename SrcIterator::row_iterator xxe = xxs + kernel_width;
    00532                 typename KernelIterator::row_iterator xk  = yk.rowIterator();
    00533
    00534                 for(; xxs < xxe; ++xxs, --xk)
    00535                 {
    00536                     sum += ak(xk) * src_acc(xxs);
    00537                 }
    00538             }
    00539
    00540             // store convolution result in destination pixel
    00541             dest_acc.set(DestTraits::fromRealPromote(sum), xd);
    00542         }
    00543     }
    00544
    00545     if(border == BORDER_TREATMENT_AVOID)
    00546         return; // skip processing near the border
    00547
    00548     int interiorskip = w + kul.x - klr.x - 1;
    00549     int borderskipx;
    00550     int borderskipy;
    00551     int borderinc;
    00552     if(border == BORDER_TREATMENT_REPEAT)
    00553     {
    00554         borderskipx = 0;
    00555         borderskipy = 0;
    00556         borderinc = 0;
    00557     }
    00558     else if(border == BORDER_TREATMENT_REFLECT)
    00559     {
    00560         borderskipx = -1;
    00561         borderskipy = -1;
    00562         borderinc = -1;
    00563     }
    00564     else if(border == BORDER_TREATMENT_WRAP)
    00565     {
    00566         borderskipx = -w+1;
    00567         borderskipy = -h+1;
    00568         borderinc = 1;
    00569     }
    00570
    00571     // create iterators for the entire image
    00572     yd = dest_ul;
    00573     ys = src_ul;
    00574
    00575     // go over the entire image (but skip the already computed points in the loop)
    00576     for(y=0; y < h; ++y, ++ys.y, ++yd.y)
    00577     {
    00578         int top    = std::max(-klr.y, src_ul.y - ys.y);
    00579         int bottom = std::min(-kul.y, src_lr.y - ys.y - 1);
    00580
    00581         // create x iterators
    00582         DestIterator xd(yd);
    00583         SrcIterator xs(ys);
    00584
    00585         for(x=0; x < w; ++x, ++xs.x, ++xd.x)
    00586         {
    00587             // check if we are away from the border
    00588             if(y >= klr.y && y < h+kul.y && x == klr.x)
    00589             {
    00590                 // yes => skip the already computed points
    00591                 x += interiorskip;
    00592                 xs.x += interiorskip;
    00593                 xd.x += interiorskip;
    00594                 continue;
    00595             }
    00596             if (border == BORDER_TREATMENT_CLIP)
    00597             {
    00598                 internalPixelEvaluationByClip(x, y, w, h, xs, src_acc, xd, dest_acc, ki, kul, klr, ak, norm);
    00599             }
    00600             else
    00601             {
    00602                 int left   = std::max(-klr.x, src_ul.x - xs.x);
    00603                 int right  = std::min(-kul.x, src_lr.x - xs.x - 1);
    00604
    00605                 // init the sum
    00606                 SumType sum = NumericTraits<SumType>::zero();
    00607
    00608                 // create iterators for the part of the kernel that fits into the image
    00609                 SrcIterator yys = xs + Size2D(0, top);
    00610                 KernelIterator yk  = ki - Size2D(0, top);
    00611
    00612                 int yy;
    00613                 for(yy = top; yy <= bottom; ++yy, ++yys.y, --yk.y)
    00614                 {
    00615                     internalPixelEvaluationByWrapReflectRepeat(yys.rowIterator(), src_acc, yk.rowIterator(), ak,
    00616                          left, right, kul.x, klr.x, borderskipx, borderinc, sum);
    00617                 }
    00618                 yys = xs + Size2D(0, top - borderskipy);
    00619                 yk  = ki - Size2D(0, top - 1);
    00620                 for(yy = top - 1; yy >= -klr.y; --yy, yys.y -= borderinc, ++yk.y)
    00621                 {
    00622                     internalPixelEvaluationByWrapReflectRepeat(yys.rowIterator(), src_acc, yk.rowIterator(), ak,
    00623                          left, right, kul.x, klr.x, borderskipx, borderinc, sum);
    00624                 }
    00625                 yys = xs + Size2D(0, bottom + borderskipy);
    00626                 yk  = ki - Size2D(0, bottom + 1);
    00627                 for(yy = bottom + 1; yy <= -kul.y; ++yy, yys.y += borderinc, --yk.y)
    00628                 {
    00629                     internalPixelEvaluationByWrapReflectRepeat(yys.rowIterator(), src_acc, yk.rowIterator(), ak,
    00630                          left, right, kul.x, klr.x, borderskipx, borderinc, sum);
    00631                 }
    00632
    00633                 // store convolution result in destination pixel
    00634                 dest_acc.set(DestTraits::fromRealPromote(sum), xd);
    00635
    00636 //                internalPixelEvaluationByWrapReflectRepeat(x, y, w, h, xs, src_acc, xd, dest_acc, ki, kul, klr, ak, border);
    00637             }
    00638         }
    00639     }
    00640 }
    00641
    00642
    00643 template <class SrcIterator, class SrcAccessor,
    00644           class DestIterator, class DestAccessor,
    00645           class KernelIterator, class KernelAccessor>
    00646 inline
    00647 void convolveImage(
    00648                    triple<SrcIterator, SrcIterator, SrcAccessor> src,
    00649                    pair<DestIterator, DestAccessor> dest,
    00650                    tuple5<KernelIterator, KernelAccessor, Diff2D, Diff2D,
    00651                    BorderTreatmentMode> kernel)
    00652 {
    00653     convolveImage(src.first, src.second, src.third,
    00654                   dest.first, dest.second,
    00655                   kernel.first, kernel.second, kernel.third,
    00656                   kernel.fourth, kernel.fifth);
    00657 }
    00658
    00659
    00660 /** \brief Performs a 2-dimensional normalized convolution, i.e. convolution with a mask image.
    00661
    00662     This functions computes
    00663     <a href ="http://homepages.inf.ed.ac.uk/rbf/CVonline/LOCAL_COPIES/PIRODDI1/NormConv/NormConv.html">normalized
    00664     convolution</a> as defined in
    00665     Knutsson, H. and Westin, C-F.: <i>Normalized and differential convolution:
    00666     Methods for Interpolation and Filtering of incomplete and uncertain data</i>.
    00667     Proc. of the IEEE Conf. on Computer Vision and Pattern Recognition, 1993, 515-523.
    00668
    00669     The mask image must be binary and encodes which pixels of the original image
    00670     are valid. It is used as follows:
    00671     Only pixel under the mask are used in the calculations. Whenever a part of the
    00672     kernel lies outside the mask, it is ignored, and the kernel is renormalized to its
    00673     original norm (analogous to the CLIP \ref BorderTreatmentMode). Thus, a useful convolution
    00674     result is computed whenever <i>at least one valid pixel is within the current window</i>
    00675     Thus, destination pixels not under the mask still receive a value if they are <i>near</i>
    00676     the mask. Therefore, this algorithm is useful as an interpolator of sparse input data.
    00677     If you are only interested in the destination values under the mask, you can perform
    00678     a subsequent \ref copyImageIf().
    00679
    00680     The KernelIterator must point to the center of the kernel, and
    00681     the kernel's size is given by its upper left (x and y of distance <= 0) and
    00682     lower right (distance >= 0) corners. The image must always be larger than the
    00683     kernel. At those positions where the kernel does not completely fit
    00684     into the image, the specified \ref BorderTreatmentMode is
    00685     applied. Only BORDER_TREATMENT_CLIP and BORDER_TREATMENT_AVOID are currently
    00686     supported.
    00687
    00688     The images's pixel type (SrcAccessor::value_type) must be a
    00689     linear space over the kernel's value_type (KernelAccessor::value_type),
    00690     i.e. addition of source values, multiplication with kernel values,
    00691     and NumericTraits must be defined.
    00692     The kernel's value_type must be an algebraic field,
    00693     i.e. the arithmetic operations (+, -, *, /) and NumericTraits must
    00694     be defined.
    00695
    00696     <b> Declarations:</b>
    00697
    00698     pass arguments explicitly:
    00699     \code
    00700     namespace vigra {
    00701         template <class SrcIterator, class SrcAccessor,
    00702                   class MaskIterator, class MaskAccessor,
    00703                   class DestIterator, class DestAccessor,
    00704                   class KernelIterator, class KernelAccessor>
    00705         void
    00706         normalizedConvolveImage(SrcIterator src_ul, SrcIterator src_lr, SrcAccessor src_acc,
    00707                                 MaskIterator mul, MaskAccessor am,
    00708                                 DestIterator dest_ul, DestAccessor dest_acc,
    00709                                 KernelIterator ki, KernelAccessor ak,
    00710                                 Diff2D kul, Diff2D klr, BorderTreatmentMode border);
    00711     }
    00712     \endcode
    00713
    00714
    00715     use argument objects in conjunction with \ref ArgumentObjectFactories:
    00716     \code
    00717     namespace vigra {
    00718         template <class SrcIterator, class SrcAccessor,
    00719                   class MaskIterator, class MaskAccessor,
    00720                   class DestIterator, class DestAccessor,
    00721                   class KernelIterator, class KernelAccessor>
    00722         inline
    00723         void normalizedConvolveImage(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    00724                                      pair<MaskIterator, MaskAccessor> mask,
    00725                                      pair<DestIterator, DestAccessor> dest,
    00726                                      tuple5<KernelIterator, KernelAccessor, Diff2D, Diff2D,
    00727                                      BorderTreatmentMode> kernel);
    00728     }
    00729     \endcode
    00730
    00731     <b> Usage:</b>
    00732
    00733     <b>\#include</b> "<a href="stdconvolution_8hxx-source.html">vigra/stdconvolution.hxx</a>"<br>
    00734     Namespace: vigra
    00735
    00736
    00737     \code
    00738     vigra::FImage src(w,h), dest(w,h);
    00739     vigra::CImage mask(w,h);
    00740     ...
    00741
    00742     // define 3x3 binomial filter
    00743     vigra::Kernel2D<float> binom;
    00744
    00745     binom.initExplicitly(Diff2D(-1,-1), Diff2D(1,1)) =   // upper left and lower right
    00746                          0.0625, 0.125, 0.0625,
    00747                          0.125,  0.25,  0.125,
    00748                          0.0625, 0.125, 0.0625;
    00749
    00750     vigra::normalizedConvolveImage(srcImageRange(src), maskImage(mask), destImage(dest), kernel2d(binom));
    00751     \endcode
    00752
    00753     <b> Required Interface:</b>
    00754
    00755     \code
    00756     ImageIterator src_ul, src_lr;
    00757     ImageIterator mul;
    00758     ImageIterator dest_ul;
    00759     ImageIterator ik;
    00760
    00761     SrcAccessor src_accessor;
    00762     MaskAccessor mask_accessor;
    00763     DestAccessor dest_accessor;
    00764     KernelAccessor kernel_accessor;
    00765
    00766     NumericTraits<SrcAccessor::value_type>::RealPromote s = src_accessor(src_ul);
    00767
    00768     s = s + s;
    00769     s = kernel_accessor(ik) * s;
    00770     s -= s;
    00771
    00772     if(mask_accessor(mul)) ...;
    00773
    00774     dest_accessor.set(
    00775     NumericTraits<DestAccessor::value_type>::fromRealPromote(s), dest_ul);
    00776
    00777     NumericTraits<KernelAccessor::value_type>::RealPromote k = kernel_accessor(ik);
    00778
    00779     k += k;
    00780     k -= k;
    00781     k = k / k;
    00782
    00783     \endcode
    00784
    00785     <b> Preconditions:</b>
    00786
    00787     \code
    00788     kul.x <= 0
    00789     kul.y <= 0
    00790     klr.x >= 0
    00791     klr.y >= 0
    00792     src_lr.x - src_ul.x >= klr.x + kul.x + 1
    00793     src_lr.y - src_ul.y >= klr.y + kul.y + 1
    00794     border == BORDER_TREATMENT_CLIP || border == BORDER_TREATMENT_AVOID
    00795     \endcode
    00796
    00797     Sum of kernel elements must be != 0.
    00798
    00799 */
    00800 template <class SrcIterator, class SrcAccessor,
    00801           class DestIterator, class DestAccessor,
    00802           class MaskIterator, class MaskAccessor,
    00803           class KernelIterator, class KernelAccessor>
    00804 void
    00805 normalizedConvolveImage(SrcIterator src_ul, SrcIterator src_lr, SrcAccessor src_acc,
    00806                         MaskIterator mul, MaskAccessor am,
    00807                         DestIterator dest_ul, DestAccessor dest_acc,
    00808                         KernelIterator ki, KernelAccessor ak,
    00809                         Diff2D kul, Diff2D klr, BorderTreatmentMode border)
    00810 {
    00811     vigra_precondition((border == BORDER_TREATMENT_CLIP  ||
    00812                         border == BORDER_TREATMENT_AVOID),
    00813                        "normalizedConvolveImage(): "
    00814                        "Border treatment must be BORDER_TREATMENT_CLIP or BORDER_TREATMENT_AVOID.");
    00815
    00816     vigra_precondition(kul.x <= 0 && kul.y <= 0,
    00817                        "normalizedConvolveImage(): left borders must be <= 0.");
    00818     vigra_precondition(klr.x >= 0 && klr.y >= 0,
    00819                        "normalizedConvolveImage(): right borders must be >= 0.");
    00820
    00821     // use traits to determine SumType as to prevent possible overflow
    00822     typedef typename
    00823         NumericTraits<typename SrcAccessor::value_type>::RealPromote SumType;
    00824     typedef typename
    00825         NumericTraits<typename KernelAccessor::value_type>::RealPromote KSumType;
    00826     typedef
    00827         NumericTraits<typename DestAccessor::value_type> DestTraits;
    00828
    00829     // calculate width and height of the image
    00830     int w = src_lr.x - src_ul.x;
    00831     int h = src_lr.y - src_ul.y;
    00832     int kernel_width = klr.x - kul.x + 1;
    00833     int kernel_height = klr.y - kul.y + 1;
    00834
    00835     int x,y;
    00836     int ystart = (border == BORDER_TREATMENT_AVOID) ?  klr.y : 0;
    00837     int yend   = (border == BORDER_TREATMENT_AVOID) ? h+kul.y : h;
    00838     int xstart = (border == BORDER_TREATMENT_AVOID) ?  klr.x : 0;
    00839     int xend   = (border == BORDER_TREATMENT_AVOID) ? w+kul.x : w;
    00840
    00841     // create y iterators
    00842     DestIterator yd = dest_ul + Diff2D(xstart, ystart);
    00843     SrcIterator ys = src_ul + Diff2D(xstart, ystart);
    00844     MaskIterator ym = mul + Diff2D(xstart, ystart);
    00845
    00846     KSumType norm = ak(ki);
    00847     int xx, yy;
    00848     KernelIterator yk  = ki + klr;
    00849     for(yy=0; yy<kernel_height; ++yy, --yk.y)
    00850     {
    00851         KernelIterator xk  = yk;
    00852
    00853         for(xx=0; xx<kernel_width; ++xx, --xk.x)
    00854         {
    00855             norm += ak(xk);
    00856         }
    00857     }
    00858     norm -= ak(ki);
    00859
    00860
    00861     for(y=ystart; y < yend; ++y, ++ys.y, ++yd.y, ++ym.y)
    00862     {
    00863         // create x iterators
    00864         DestIterator xd(yd);
    00865         SrcIterator xs(ys);
    00866         MaskIterator xm(ym);
    00867
    00868         for(x=xstart; x < xend; ++x, ++xs.x, ++xd.x, ++xm.x)
    00869         {
    00870             // how much of the kernel fits into the image ?
    00871             int x0, y0, x1, y1;
    00872
    00873             y0 = (y<klr.y) ? -y : -klr.y;
    00874             y1 = (h-y-1<-kul.y) ? h-y-1 : -kul.y;
    00875             x0 = (x<klr.x) ? -x : -klr.x;
    00876             x1 = (w-x-1<-kul.x) ? w-x-1 : -kul.x;
    00877
    00878             bool first = true;
    00879             // init the sum
    00880             SumType sum;
    00881             KSumType ksum;
    00882
    00883             SrcIterator yys = xs + Diff2D(x0, y0);
    00884             MaskIterator yym = xm + Diff2D(x0, y0);
    00885             KernelIterator yk  = ki - Diff2D(x0, y0);
    00886
    00887             int xx, kernel_width, kernel_height;
    00888             kernel_width = x1 - x0 + 1;
    00889             kernel_height = y1 - y0 + 1;
    00890             for(yy=0; yy<kernel_height; ++yy, ++yys.y, --yk.y, ++yym.y)
    00891             {
    00892                 typename SrcIterator::row_iterator xxs = yys.rowIterator();
    00893                 typename SrcIterator::row_iterator xxend = xxs + kernel_width;
    00894                 typename MaskIterator::row_iterator xxm = yym.rowIterator();
    00895                 typename KernelIterator::row_iterator xk  = yk.rowIterator();
    00896
    00897                 for(xx=0; xxs < xxend; ++xxs.x, --xk.x, ++xxm.x)
    00898                 {
    00899                     if(!am(xxm)) continue;
    00900
    00901                     if(first)
    00902                     {
    00903                         sum = ak(xk) * src_acc(xxs);
    00904                         ksum = ak(xk);
    00905                         first = false;
    00906                     }
    00907                     else
    00908                     {
    00909                         sum += ak(xk) * src_acc(xxs);
    00910                         ksum += ak(xk);
    00911                     }
    00912                 }
    00913             }
    00914             // store average in destination pixel
    00915             if(!first &&
    00916                ksum != NumericTraits<KSumType>::zero())
    00917             {
    00918                 dest_acc.set(DestTraits::fromRealPromote((norm / ksum) * sum), xd);
    00919             }
    00920         }
    00921     }
    00922 }
    00923
    00924
    00925 template <class SrcIterator, class SrcAccessor,
    00926           class DestIterator, class DestAccessor,
    00927           class MaskIterator, class MaskAccessor,
    00928           class KernelIterator, class KernelAccessor>
    00929 inline
    00930 void normalizedConvolveImage(
    00931                            triple<SrcIterator, SrcIterator, SrcAccessor> src,
    00932                            pair<MaskIterator, MaskAccessor> mask,
    00933                            pair<DestIterator, DestAccessor> dest,
    00934                            tuple5<KernelIterator, KernelAccessor, Diff2D, Diff2D,
    00935                            BorderTreatmentMode> kernel)
    00936 {
    00937     normalizedConvolveImage(src.first, src.second, src.third,
    00938                             mask.first, mask.second,
    00939                             dest.first, dest.second,
    00940                             kernel.first, kernel.second, kernel.third,
    00941                             kernel.fourth, kernel.fifth);
    00942 }
    00943
    00944 /** \brief Deprecated name of 2-dimensional normalized convolution, i.e. convolution with a mask image.
    00945
    00946     See \ref normalizedConvolveImage() for documentation.
    00947
    00948     <b> Declarations:</b>
    00949
    00950     pass arguments explicitly:
    00951     \code
    00952     namespace vigra {
    00953         template <class SrcIterator, class SrcAccessor,
    00954                   class MaskIterator, class MaskAccessor,
    00955                   class DestIterator, class DestAccessor,
    00956                   class KernelIterator, class KernelAccessor>
    00957         void
    00958         convolveImageWithMask(SrcIterator src_ul, SrcIterator src_lr, SrcAccessor src_acc,
    00959                               MaskIterator mul, MaskAccessor am,
    00960                               DestIterator dest_ul, DestAccessor dest_acc,
    00961                               KernelIterator ki, KernelAccessor ak,
    00962                               Diff2D kul, Diff2D klr, BorderTreatmentMode border);
    00963     }
    00964     \endcode
    00965
    00966
    00967     use argument objects in conjunction with \ref ArgumentObjectFactories:
    00968     \code
    00969     namespace vigra {
    00970         template <class SrcIterator, class SrcAccessor,
    00971                   class MaskIterator, class MaskAccessor,
    00972                   class DestIterator, class DestAccessor,
    00973                   class KernelIterator, class KernelAccessor>
    00974         inline
    00975         void convolveImageWithMask(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    00976                                    pair<MaskIterator, MaskAccessor> mask,
    00977                                    pair<DestIterator, DestAccessor> dest,
    00978                                    tuple5<KernelIterator, KernelAccessor, Diff2D, Diff2D,
    00979                                    BorderTreatmentMode> kernel);
    00980     }
    00981     \endcode
    00982 */
    00983 template <class SrcIterator, class SrcAccessor,
    00984           class DestIterator, class DestAccessor,
    00985           class MaskIterator, class MaskAccessor,
    00986           class KernelIterator, class KernelAccessor>
    00987 inline void
    00988 convolveImageWithMask(SrcIterator src_ul, SrcIterator src_lr, SrcAccessor src_acc,
    00989                       MaskIterator mul, MaskAccessor am,
    00990                       DestIterator dest_ul, DestAccessor dest_acc,
    00991                       KernelIterator ki, KernelAccessor ak,
    00992                       Diff2D kul, Diff2D klr, BorderTreatmentMode border)
    00993 {
    00994     normalizedConvolveImage(src_ul, src_lr, src_acc,
    00995                             mul, am,
    00996                             dest_ul, dest_acc,
    00997                             ki, ak, kul, klr, border);
    00998 }
    00999
    01000 template <class SrcIterator, class SrcAccessor,
    01001           class DestIterator, class DestAccessor,
    01002           class MaskIterator, class MaskAccessor,
    01003           class KernelIterator, class KernelAccessor>
    01004 inline
    01005 void convolveImageWithMask(
    01006                            triple<SrcIterator, SrcIterator, SrcAccessor> src,
    01007                            pair<MaskIterator, MaskAccessor> mask,
    01008                            pair<DestIterator, DestAccessor> dest,
    01009                            tuple5<KernelIterator, KernelAccessor, Diff2D, Diff2D,
    01010                            BorderTreatmentMode> kernel)
    01011 {
    01012     normalizedConvolveImage(src.first, src.second, src.third,
    01013                             mask.first, mask.second,
    01014                             dest.first, dest.second,
    01015                             kernel.first, kernel.second, kernel.third,
    01016                             kernel.fourth, kernel.fifth);
    01017 }
    01018
    01019 //@}
    01020
    01021 /********************************************************/
    01022 /*                                                      */
    01023 /*                      Kernel2D                        */
    01024 /*                                                      */
    01025 /********************************************************/
    01026
    01027 /** \brief Generic 2 dimensional convolution kernel.
    01028
    01029     This kernel may be used for convolution of 2 dimensional signals.
    01030
    01031     Convolution functions access the kernel via an ImageIterator
    01032     which they get by calling \ref center(). This iterator
    01033     points to the center of the kernel. The kernel's size is given by its upperLeft()
    01034     (upperLeft().x <= 0, upperLeft().y <= 0)
    01035     and lowerRight() (lowerRight().x >= 0, lowerRight().y >= 0) methods.
    01036     The desired border treatment mode is returned by borderTreatment().
    01037     (Note that the \ref StandardConvolution "2D convolution functions" don't currently
    01038     support all modes.)
    01039
    01040     The different init functions create a kernel with the specified
    01041     properties. The requirements for the kernel's value_type depend
    01042     on the init function used. At least NumericTraits must be defined.
    01043
    01044     The kernel defines a factory function kernel2d() to create an argument object
    01045     (see \ref KernelArgumentObjectFactories).
    01046
    01047     <b> Usage:</b>
    01048
    01049     <b>\#include</b> "<a href="stdconvolution_8hxx-source.html">vigra/stdconvolution.hxx</a>"<br>
    01050     Namespace: vigra
    01051
    01052     \code
    01053     vigra::FImage src(w,h), dest(w,h);
    01054     ...
    01055
    01056     // define horizontal Sobel filter
    01057     vigra::Kernel2D<float> sobel;
    01058
    01059     sobel.initExplicitly(Diff2D(-1,-1), Diff2D(1,1)) =  // upper left and lower right
    01060                          0.125, 0.0, -0.125,
    01061                          0.25,  0.0, -0.25,
    01062                          0.125, 0.0, -0.125;
    01063
    01064     vigra::convolveImage(srcImageRange(src), destImage(dest), kernel2d(sobel));
    01065     \endcode
    01066
    01067     <b> Required Interface:</b>
    01068
    01069     \code
    01070     value_type v = NumericTraits<value_type>::one();
    01071     \endcode
    01072
    01073     See also the init functions.
    01074 */
    01075 template <class ARITHTYPE>
    01076 class Kernel2D
    01077 {
    01078 public:
    01079         /** the kernel's value type
    01080          */
    01081     typedef ARITHTYPE value_type;
    01082
    01083         /** 2D random access iterator over the kernel's values
    01084          */
    01085     typedef typename BasicImage<value_type>::traverser Iterator;
    01086
    01087         /** const 2D random access iterator over the kernel's values
    01088          */
    01089     typedef typename BasicImage<value_type>::const_traverser ConstIterator;
    01090
    01091         /** the kernel's accessor
    01092          */
    01093     typedef typename BasicImage<value_type>::Accessor Accessor;
    01094
    01095         /** the kernel's const accessor
    01096          */
    01097     typedef typename BasicImage<value_type>::ConstAccessor ConstAccessor;
    01098
    01099     struct InitProxy
    01100     {
    01101         typedef typename
    01102         BasicImage<value_type>::ScanOrderIterator Iterator;
    01103
    01104         InitProxy(Iterator i, int count, value_type & norm)
    01105             : iter_(i), base_(i),
    01106               count_(count), sum_(count),
    01107               norm_(norm)
    01108         {}
    01109
    01110         ~InitProxy()
    01111         {
    01112             vigra_precondition(count_ == 1 || count_ == sum_,
    01113                                "Kernel2D::initExplicitly(): "
    01114                                "Too few init values.");
    01115         }
    01116
    01117         InitProxy & operator,(value_type const & v)
    01118         {
    01119             if(count_ == sum_)  norm_ = *iter_;
    01120
    01121             --count_;
    01122             vigra_precondition(count_ > 0,
    01123                                "Kernel2D::initExplicitly(): "
    01124                                "Too many init values.");
    01125
    01126             norm_ += v;
    01127
    01128             ++iter_;
    01129             *iter_ = v;
    01130
    01131             return *this;
    01132         }
    01133
    01134         Iterator iter_, base_;
    01135         int count_, sum_;
    01136         value_type & norm_;
    01137     };
    01138
    01139     static value_type one() { return NumericTraits<value_type>::one(); }
    01140
    01141         /** Default constructor.
    01142             Creates a kernel of size 1x1 which would copy the signal
    01143             unchanged.
    01144         */
    01145     Kernel2D()
    01146         : kernel_(1, 1, one()),
    01147           left_(0, 0),
    01148           right_(0, 0),
    01149       norm_(one()),
    01150           border_treatment_(BORDER_TREATMENT_CLIP)
    01151     {}
    01152
    01153         /** Copy constructor.
    01154          */
    01155     Kernel2D(Kernel2D const & k)
    01156         : kernel_(k.kernel_),
    01157           left_(k.left_),
    01158           right_(k.right_),
    01159           norm_(k.norm_),
    01160           border_treatment_(k.border_treatment_)
    01161     {}
    01162
    01163         /** Copy assignment.
    01164          */
    01165     Kernel2D & operator=(Kernel2D const & k)
    01166     {
    01167         if(this != &k)
    01168         {
    01169         kernel_ = k.kernel_;
    01170             left_ = k.left_;
    01171             right_ = k.right_;
    01172             norm_ = k.norm_;
    01173         border_treatment_ = k.border_treatment_;
    01174         }
    01175         return *this;
    01176     }
    01177
    01178         /** Initialization.
    01179             This initializes the kernel with the given constant. The norm becomes
    01180             v*width()*height().
    01181
    01182             Instead of a single value an initializer list of length width()*height()
    01183             can be used like this:
    01184
    01185             \code
    01186             vigra::Kernel2D<float> binom;
    01187
    01188             binom.initExplicitly(Diff2D(-1,-1), Diff2D(1,1)) =
    01189             0.0625, 0.125, 0.0625,
    01190             0.125,  0.25,  0.125,
    01191             0.0625, 0.125, 0.0625;
    01192             \endcode
    01193
    01194             In this case, the norm will be set to the sum of the init values.
    01195             An initializer list of wrong length will result in a run-time error.
    01196         */
    01197     InitProxy operator=(value_type const & v)
    01198     {
    01199         int size = (right_.x - left_.x + 1) *
    01200                    (right_.y - left_.y + 1);
    01201         kernel_ = v;
    01202         norm_ = (double)size*v;
    01203
    01204         return InitProxy(kernel_.begin(), size, norm_);
    01205     }
    01206
    01207         /** Destructor.
    01208          */
    01209     ~Kernel2D()
    01210     {}
    01211
    01212         /** Init the 2D kernel as the cartesian product of two 1D kernels
    01213             of type \ref Kernel1D. The norm becomes the product of the two original
    01214             norms.
    01215
    01216             <b> Required Interface:</b>
    01217
    01218             The kernel's value_type must be a linear algebra.
    01219
    01220             \code
    01221             vigra::Kernel2D<...>::value_type v;
    01222             v = v * v;
    01223             \endcode
    01224         */
    01225     void initSeparable(Kernel1D<value_type> & kx,
    01226                        Kernel1D<value_type> & ky)
    01227     {
    01228         left_ = Diff2D(kx.left(), ky.left());
    01229         right_ = Diff2D(kx.right(), ky.right());
    01230         int w = right_.x - left_.x + 1;
    01231         int h = right_.y - left_.y + 1;
    01232         kernel_.resize(w, h);
    01233
    01234         norm_ = kx.norm() * ky.norm();
    01235
    01236         typedef typename Kernel1D<value_type>::Iterator KIter;
    01237         typename Kernel1D<value_type>::Accessor ka;
    01238
    01239         KIter kiy = ky.center() + left_.y;
    01240         Iterator iy = center() + left_;
    01241
    01242         for(int y=left_.y; y<=right_.y; ++y, ++kiy, ++iy.y)
    01243         {
    01244             KIter kix = kx.center() + left_.x;
    01245             Iterator ix = iy;
    01246             for(int x=left_.x; x<=right_.x; ++x, ++kix, ++ix.x)
    01247             {
    01248                 *ix = ka(kix) * ka(kiy);
    01249             }
    01250         }
    01251     }
    01252
    01253         /** Init the 2D kernel as the cartesian product of two 1D kernels
    01254             given explicitly by iterators and sizes. The norm becomes the
    01255             sum of the resulting kernel values.
    01256
    01257             <b> Required Interface:</b>
    01258
    01259             The kernel's value_type must be a linear algebra.
    01260
    01261             \code
    01262             vigra::Kernel2D<...>::value_type v;
    01263             v = v * v;
    01264             v += v;
    01265             \endcode
    01266
    01267             <b> Preconditions:</b>
    01268
    01269             \code
    01270             xleft <= 0;
    01271             xright >= 0;
    01272             yleft <= 0;
    01273             yright >= 0;
    01274             \endcode
    01275         */
    01276     template <class KernelIterator>
    01277     void initSeparable(KernelIterator kxcenter, int xleft, int xright,
    01278                        KernelIterator kycenter, int yleft, int yright)
    01279     {
    01280         vigra_precondition(xleft <= 0 && yleft <= 0,
    01281                            "Kernel2D::initSeparable(): left borders must be <= 0.");
    01282         vigra_precondition(xright >= 0 && yright >= 0,
    01283                            "Kernel2D::initSeparable(): right borders must be >= 0.");
    01284
    01285         left_ = Point2D(xleft, yleft);
    01286         right_ = Point2D(xright, yright);
    01287
    01288         int w = right_.x - left_.x + 1;
    01289         int h = right_.y - left_.y + 1;
    01290         kernel_.resize(w, h);
    01291
    01292         KernelIterator kiy = kycenter + left_.y;
    01293         Iterator iy = center() + left_;
    01294
    01295         for(int y=left_.y; y<=right_.y; ++y, ++kiy, ++iy.y)
    01296         {
    01297             KernelIterator kix = kxcenter + left_.x;
    01298             Iterator ix = iy;
    01299             for(int x=left_.x; x<=right_.x; ++x, ++kix, ++ix.x)
    01300             {
    01301                 *ix = *kix * *kiy;
    01302             }
    01303         }
    01304
    01305         typename BasicImage<value_type>::iterator i = kernel_.begin();
    01306         typename BasicImage<value_type>::iterator iend = kernel_.end();
    01307         norm_ = *i;
    01308         ++i;
    01309
    01310         for(; i!= iend; ++i)
    01311         {
    01312             norm_ += *i;
    01313         }
    01314     }
    01315
    01316         /** Init the 2D kernel as a circular averaging filter. The norm will be
    01317             calculated as
    01318             <TT>NumericTraits<value_type>::one() / (number of non-zero kernel values)</TT>.
    01319             The kernel's value_type must be a linear space.
    01320
    01321             <b> Required Interface:</b>
    01322
    01323             \code
    01324             value_type v = vigra::NumericTraits<value_type>::one();
    01325
    01326             double d;
    01327             v = d * v;
    01328             \endcode
    01329
    01330             <b> Precondition:</b>
    01331
    01332             \code
    01333             radius > 0;
    01334             \endcode
    01335         */
    01336     void initDisk(int radius)
    01337     {
    01338         vigra_precondition(radius > 0,
    01339                            "Kernel2D::initDisk(): radius must be > 0.");
    01340
    01341         left_ = Point2D(-radius, -radius);
    01342         right_ = Point2D(radius, radius);
    01343         int w = right_.x - left_.x + 1;
    01344         int h = right_.y - left_.y + 1;
    01345         kernel_.resize(w, h);
    01346         norm_ = NumericTraits<value_type>::one();
    01347
    01348         kernel_ = NumericTraits<value_type>::zero();
    01349         double count = 0.0;
    01350
    01351         Iterator k = center();
    01352         double r2 = (double)radius*radius;
    01353
    01354         int i;
    01355         for(i=0; i<= radius; ++i)
    01356         {
    01357             double r = (double) i - 0.5;
    01358             int w = (int)(VIGRA_CSTD::sqrt(r2 - r*r) + 0.5);
    01359             for(int j=-w; j<=w; ++j)
    01360             {
    01361                 k(j, i) = NumericTraits<value_type>::one();
    01362                 k(j, -i) = NumericTraits<value_type>::one();
    01363                 count += (i != 0) ? 2.0 : 1.0;
    01364             }
    01365         }
    01366
    01367         count = 1.0 / count;
    01368
    01369         for(int y=-radius; y<=radius; ++y)
    01370         {
    01371             for(int x=-radius; x<=radius; ++x)
    01372             {
    01373                 k(x,y) = count * k(x,y);
    01374             }
    01375         }
    01376     }
    01377
    01378         /** Init the kernel by an explicit initializer list.
    01379             The upper left and lower right corners of the kernel must be passed.
    01380             A comma-separated initializer list is given after the assignment operator.
    01381             This function is used like this:
    01382
    01383             \code
    01384             // define horizontal Sobel filter
    01385             vigra::Kernel2D<float> sobel;
    01386
    01387             sobel.initExplicitly(Diff2D(-1,-1), Diff2D(1,1)) =
    01388             0.125, 0.0, -0.125,
    01389             0.25,  0.0, -0.25,
    01390             0.125, 0.0, -0.125;
    01391             \endcode
    01392
    01393             The norm is set to the sum of the initialzer values. If the wrong number of
    01394             values is given, a run-time error results. It is, however, possible to give
    01395             just one initializer. This creates an averaging filter with the given constant:
    01396
    01397             \code
    01398             vigra::Kernel2D<float> average3x3;
    01399
    01400             average3x3.initExplicitly(Diff2D(-1,-1), Diff2D(1,1)) = 1.0/9.0;
    01401             \endcode
    01402
    01403             Here, the norm is set to value*width()*height().
    01404
    01405             <b> Preconditions:</b>
    01406
    01407             \code
    01408             1. upperleft.x <= 0;
    01409             2. upperleft.y <= 0;
    01410             3. lowerright.x >= 0;
    01411             4. lowerright.y >= 0;
    01412             5. the number of values in the initializer list
    01413             is 1 or equals the size of the kernel.
    01414             \endcode
    01415         */
    01416     Kernel2D & initExplicitly(Diff2D upperleft, Diff2D lowerright)
    01417     {
    01418         vigra_precondition(upperleft.x <= 0 && upperleft.y <= 0,
    01419                            "Kernel2D::initExplicitly(): left borders must be <= 0.");
    01420         vigra_precondition(lowerright.x >= 0 && lowerright.y >= 0,
    01421                            "Kernel2D::initExplicitly(): right borders must be >= 0.");
    01422
    01423         left_ = Point2D(upperleft);
    01424         right_ = Point2D(lowerright);
    01425
    01426         int w = right_.x - left_.x + 1;
    01427         int h = right_.y - left_.y + 1;
    01428         kernel_.resize(w, h);
    01429
    01430         return *this;
    01431     }
    01432
    01433         /** Coordinates of the upper left corner of the kernel.
    01434          */
    01435     Point2D upperLeft() const { return left_; }
    01436
    01437         /** Coordinates of the lower right corner of the kernel.
    01438          */
    01439     Point2D lowerRight() const { return right_; }
    01440
    01441         /** Width of the kernel.
    01442          */
    01443     int width() const { return right_.x - left_.x + 1; }
    01444
    01445         /** Height of the kernel.
    01446          */
    01447     int height() const { return right_.y - left_.y + 1; }
    01448
    01449         /** ImageIterator that points to the center of the kernel (coordinate (0,0)).
    01450          */
    01451     Iterator center() { return kernel_.upperLeft() - left_; }
    01452
    01453         /** ImageIterator that points to the center of the kernel (coordinate (0,0)).
    01454          */
    01455     ConstIterator center() const { return kernel_.upperLeft() - left_; }
    01456
    01457         /** Access kernel entry at given position.
    01458          */
    01459     value_type & operator()(int x, int y)
    01460     { return kernel_[Diff2D(x,y) - left_]; }
    01461
    01462         /** Read kernel entry at given position.
    01463          */
    01464     value_type operator()(int x, int y) const
    01465     { return kernel_[Diff2D(x,y) - left_]; }
    01466
    01467         /** Access kernel entry at given position.
    01468          */
    01469     value_type & operator[](Diff2D const & d)
    01470     { return kernel_[d - left_]; }
    01471
    01472         /** Read kernel entry at given position.
    01473          */
    01474     value_type operator[](Diff2D const & d) const
    01475     { return kernel_[d - left_]; }
    01476
    01477         /** Norm of the kernel (i.e. sum of its elements).
    01478          */
    01479     value_type norm() const { return norm_; }
    01480
    01481         /** The kernels default accessor.
    01482          */
    01483     Accessor accessor() { return Accessor(); }
    01484
    01485         /** The kernels default const accessor.
    01486          */
    01487     ConstAccessor accessor() const { return ConstAccessor(); }
    01488
    01489         /** Normalize the kernel to the given value. (The norm is the sum of all kernel
    01490             elements.) The kernel's value_type must be a division algebra or
    01491             algebraic field.
    01492
    01493             <b> Required Interface:</b>
    01494
    01495             \code
    01496             value_type v = vigra::NumericTraits<value_type>::one(); // if norm is not
    01497                                                                     // given explicitly
    01498
    01499             v += v;
    01500             v = v * v;
    01501             v = v / v;
    01502             \endcode
    01503         */
    01504     void normalize(value_type norm)
    01505     {
    01506         typename BasicImage<value_type>::iterator i = kernel_.begin();
    01507         typename BasicImage<value_type>::iterator iend = kernel_.end();
    01508         typename NumericTraits<value_type>::RealPromote sum = *i;
    01509         ++i;
    01510
    01511         for(; i!= iend; ++i)
    01512         {
    01513             sum += *i;
    01514         }
    01515
    01516         sum = norm / sum;
    01517         i = kernel_.begin();
    01518         for(; i != iend; ++i)
    01519         {
    01520             *i = *i * sum;
    01521         }
    01522
    01523         norm_ = norm;
    01524     }
    01525
    01526         /** Normalize the kernel to norm 1.
    01527          */
    01528     void normalize()
    01529     {
    01530         normalize(one());
    01531     }
    01532
    01533         /** current border treatment mode
    01534          */
    01535     BorderTreatmentMode borderTreatment() const
    01536     { return border_treatment_; }
    01537
    01538         /** Set border treatment mode.
    01539             Only <TT>BORDER_TREATMENT_CLIP</TT> and <TT>BORDER_TREATMENT_AVOID</TT> are currently
    01540             allowed.
    01541         */
    01542     void setBorderTreatment( BorderTreatmentMode new_mode)
    01543     {
    01544         vigra_precondition((new_mode == BORDER_TREATMENT_CLIP    ||
    01545                             new_mode == BORDER_TREATMENT_AVOID   ||
    01546                             new_mode == BORDER_TREATMENT_REFLECT ||
    01547                             new_mode == BORDER_TREATMENT_REPEAT  ||
    01548                             new_mode == BORDER_TREATMENT_WRAP),
    01549                            "convolveImage():\n"
    01550                            "  Border treatment must be one of follow treatments:\n"
    01551                            "  - BORDER_TREATMENT_CLIP\n"
    01552                            "  - BORDER_TREATMENT_AVOID\n"
    01553                            "  - BORDER_TREATMENT_REFLECT\n"
    01554                            "  - BORDER_TREATMENT_REPEAT\n"
    01555                            "  - BORDER_TREATMENT_WRAP\n");
    01556
    01557         border_treatment_ = new_mode;
    01558     }
    01559
    01560
    01561 private:
    01562     BasicImage<value_type> kernel_;
    01563     Point2D left_, right_;
    01564     value_type norm_;
    01565     BorderTreatmentMode border_treatment_;
    01566 };
    01567
    01568 /**************************************************************/
    01569 /*                                                            */
    01570 /*         Argument object factories for Kernel2D             */
    01571 /*                                                            */
    01572 /*     (documentation: see vigra/convolution.hxx)             */
    01573 /*                                                            */
    01574 /**************************************************************/
    01575
    01576 template <class KernelIterator, class KernelAccessor>
    01577 inline
    01578 tuple5<KernelIterator, KernelAccessor, Diff2D, Diff2D, BorderTreatmentMode>
    01579 kernel2d(KernelIterator ik, KernelAccessor ak, Diff2D kul, Diff2D klr,
    01580          BorderTreatmentMode border)
    01581
    01582 {
    01583     return
    01584         tuple5<KernelIterator, KernelAccessor, Diff2D, Diff2D, BorderTreatmentMode> (
    01585                                                              ik, ak, kul, klr, border);
    01586 }
    01587
    01588 template <class T>
    01589 inline
    01590 tuple5<typename Kernel2D<T>::ConstIterator,
    01591        typename Kernel2D<T>::ConstAccessor,
    01592        Diff2D, Diff2D, BorderTreatmentMode>
    01593 kernel2d(Kernel2D<T> const & k)
    01594
    01595 {
    01596     return
    01597         tuple5<typename Kernel2D<T>::ConstIterator,
    01598                typename Kernel2D<T>::ConstAccessor,
    01599                Diff2D, Diff2D, BorderTreatmentMode>(
    01600             k.center(),
    01601             k.accessor(),
    01602             k.upperLeft(), k.lowerRight(),
    01603             k.borderTreatment());
    01604 }
    01605
    01606 template <class T>
    01607 inline
    01608 tuple5<typename Kernel2D<T>::ConstIterator,
    01609        typename Kernel2D<T>::ConstAccessor,
    01610        Diff2D, Diff2D, BorderTreatmentMode>
    01611 kernel2d(Kernel2D<T> const & k, BorderTreatmentMode border)
    01612
    01613 {
    01614     return
    01615         tuple5<typename Kernel2D<T>::ConstIterator,
    01616                typename Kernel2D<T>::ConstAccessor,
    01617                Diff2D, Diff2D, BorderTreatmentMode>(
    01618             k.center(),
    01619             k.accessor(),
    01620             k.upperLeft(), k.lowerRight(),
    01621             border);
    01622 }
    01623
    01624
    01625 } // namespace vigra
    01626
    01627 #endif // VIGRA_STDCONVOLUTION_HXX

© Ullrich Köthe (koethe@informatik.uni-hamburg.de) Cognitive Systems Group, University of Hamburg, Germany | html generated using doxygen and Python VIGRA 1.3.2 (27 Jan 2005)