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

details vigra/stdconvolution.hxx VIGRA

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)