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