[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | 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 /*    The VIGRA Website is                                              */
00008 /*        http://kogs-www.informatik.uni-hamburg.de/~koethe/vigra/      */
00009 /*    Please direct questions, bug reports, and contributions to        */
00010 /*        ullrich.koethe@iwr.uni-heidelberg.de    or                    */
00011 /*        vigra@informatik.uni-hamburg.de                               */
00012 /*                                                                      */
00013 /*    Permission is hereby granted, free of charge, to any person       */
00014 /*    obtaining a copy of this software and associated documentation    */
00015 /*    files (the "Software"), to deal in the Software without           */
00016 /*    restriction, including without limitation the rights to use,      */
00017 /*    copy, modify, merge, publish, distribute, sublicense, and/or      */
00018 /*    sell copies of the Software, and to permit persons to whom the    */
00019 /*    Software is furnished to do so, subject to the following          */
00020 /*    conditions:                                                       */
00021 /*                                                                      */
00022 /*    The above copyright notice and this permission notice shall be    */
00023 /*    included in all copies or substantial portions of the             */
00024 /*    Software.                                                         */
00025 /*                                                                      */
00026 /*    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND    */
00027 /*    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   */
00028 /*    OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          */
00029 /*    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT       */
00030 /*    HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,      */
00031 /*    WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING      */
00032 /*    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR     */
00033 /*    OTHER DEALINGS IN THE SOFTWARE.                                   */
00034 /*                                                                      */
00035 /************************************************************************/
00036 
00037 
00038 #ifndef VIGRA_STDCONVOLUTION_HXX
00039 #define VIGRA_STDCONVOLUTION_HXX
00040 
00041 #include <cmath>
00042 #include "stdimage.hxx"
00043 #include "bordertreatment.hxx"
00044 #include "separableconvolution.hxx"
00045 #include "utilities.hxx"
00046 #include "sized_int.hxx"
00047 
00048 namespace vigra {
00049 
00050 template <class SrcIterator, class SrcAccessor,
00051           class DestIterator, class DestAccessor,
00052           class KernelIterator, class KernelAccessor,
00053           class KSumType>
00054 void internalPixelEvaluationByClip(int x, int y, int w, int h, SrcIterator xs,
00055                                    SrcAccessor src_acc, DestIterator xd, DestAccessor dest_acc,
00056                                    KernelIterator ki, Diff2D kul, Diff2D klr, KernelAccessor ak,
00057                                    KSumType norm)
00058 {
00059     typedef typename
00060         NumericTraits<typename SrcAccessor::value_type>::RealPromote SumType;
00061     typedef
00062         NumericTraits<typename DestAccessor::value_type> DestTraits;
00063 
00064     // calculate width and height of the kernel
00065     int kernel_width = klr.x - kul.x + 1;
00066     int kernel_height = klr.y - kul.y + 1;
00067 
00068     SumType sum = NumericTraits<SumType>::zero();
00069     int xx, yy;
00070     int x0, y0, x1, y1;
00071 
00072     y0 = (y<klr.y) ?  -y : -klr.y;
00073     y1 = (h-y-1<-kul.y) ? h-y-1 : -kul.y;
00074 
00075     x0 = (x<klr.x) ? -x : -klr.x;
00076     x1 = (w-x-1<-kul.x) ? w-x-1 : -kul.x;
00077 
00078     SrcIterator yys = xs + Diff2D(x0, y0);
00079     KernelIterator yk  = ki - Diff2D(x0, y0);
00080 
00081     KSumType ksum = NumericTraits<KSumType>::zero();
00082     kernel_width = x1 - x0 + 1;
00083     kernel_height = y1 - y0 + 1;
00084 
00085     //es wird zuerst abgeschnitten und dann gespigelt!
00086 
00087     for(yy=0; yy<kernel_height; ++yy, ++yys.y, --yk.y)
00088     {
00089         SrcIterator xxs = yys;
00090         KernelIterator xk  = yk;
00091 
00092         for(xx=0; xx<kernel_width; ++xx, ++xxs.x, --xk.x)
00093         {
00094             sum += ak(xk) * src_acc(xxs);
00095             ksum += ak(xk);
00096         }
00097     }
00098 
00099     //                      store average in destination pixel
00100     dest_acc.set(DestTraits::fromRealPromote((norm / ksum) * sum), xd);
00101 
00102 }
00103 
00104 
00105 #if 0
00106 
00107 template <class SrcIterator, class SrcAccessor,
00108           class DestIterator, class DestAccessor,
00109           class KernelIterator, class KernelAccessor>
00110 void internalPixelEvaluationByWrapReflectRepeat(int x, int y, int src_width, int src_height, SrcIterator xs,
00111                                                 SrcAccessor src_acc, DestIterator xd, DestAccessor dest_acc,
00112                                                 KernelIterator ki, Diff2D kul, Diff2D klr, KernelAccessor ak,
00113                                                 BorderTreatmentMode border)
00114 {
00115 
00116     typedef typename
00117         NumericTraits<typename SrcAccessor::value_type>::RealPromote SumType;
00118     typedef
00119         NumericTraits<typename DestAccessor::value_type> DestTraits;
00120 
00121     SumType sum = NumericTraits<SumType>::zero();
00122 
00123     SrcIterator src_ul = xs - Diff2D(x, y);
00124     SrcIterator src_lr = src_ul + Diff2D(src_width, src_height);
00125 
00126     SrcIterator yys = xs;
00127     KernelIterator yk  = ki;
00128 
00129     // calculate width and height of the kernel
00130     int kernel_width = klr.x - kul.x + 1;
00131     int kernel_height = klr.y - kul.y + 1;
00132 
00133     // where the kernel is beyond the borders:
00134     bool top_to_much = (y<klr.y) ? true : false;
00135     bool down_to_much = (src_height-y-1<-kul.y)? true : false;
00136     bool left_to_much = (x<klr.x)? true : false;
00137     bool right_to_much = (src_width-x-1<-kul.x)? true : false;
00138 
00139     // direction of iteration,
00140     // e.g. (-1, +1) for ll->ur or (-1, -1) for lr->ul
00141     Diff2D way_increment;
00142 
00143     /* Iteration is always done from valid to invalid range.
00144        The following tuple is composed as such:
00145        - If an invalid range is reached while iterating in X,
00146          a jump of border_increment.first is performed and
00147          border_increment.third is used for further iterating.
00148        - If an invalid range is reached while iterating in Y,
00149          a jump of border_increment.second is performed and
00150          border_increment.fourth is used for further iterating.
00151     */
00152     tuple4<int, int, int, int> border_increment;
00153     if (border == BORDER_TREATMENT_REPEAT){
00154         border_increment = tuple4<int, int, int, int>(1, 1, 0, 0);
00155     }else if (border == BORDER_TREATMENT_REFLECT){
00156         border_increment = tuple4<int, int, int, int>(2, 2, -1, -1);
00157     }else{ // BORDER_TREATMENT_WRAP
00158         border_increment = tuple4<int, int, int, int>(src_width, src_height, 1, 1);
00159     }
00160 
00161     pair<int, int> valid_step_count;
00162 
00163     if(left_to_much && !top_to_much && !down_to_much)
00164     {
00165         yys += klr;
00166         yk += kul;
00167         way_increment = Diff2D(-1, -1);
00168         border_increment.third = -border_increment.third;
00169         border_increment.fourth = -border_increment.fourth;
00170         valid_step_count = std::make_pair((yys - src_ul).x + 1, kernel_height);
00171     }
00172     else if(top_to_much && !left_to_much && !right_to_much)
00173     {
00174         yys += klr;
00175         yk += kul;
00176         way_increment = Diff2D(-1, -1);
00177         border_increment.third = -border_increment.third;
00178         border_increment.fourth = -border_increment.fourth;
00179         valid_step_count = std::make_pair(kernel_width, (yys - src_ul).y + 1);
00180     }
00181     else if(right_to_much && !top_to_much && !down_to_much)
00182     {
00183         yys += kul;
00184         yk += klr;
00185         way_increment = Diff2D(1, 1);
00186         border_increment.first = -border_increment.first;
00187         border_increment.second = -border_increment.second;
00188         valid_step_count = std::make_pair((src_lr - yys).x, kernel_height);
00189     }
00190     else if(down_to_much && !left_to_much && !right_to_much)
00191     {
00192         yys += kul;
00193         yk += klr;
00194         way_increment = Diff2D(1, 1);
00195         border_increment.first = -border_increment.first;
00196         border_increment.second = -border_increment.second;
00197         valid_step_count = std::make_pair(kernel_width, (src_lr - yys).y);
00198     }
00199     else if(down_to_much && left_to_much)
00200     {
00201         yys += kul + Diff2D(kernel_width - 1, 0);
00202         yk += kul + Diff2D(0, kernel_height - 1);
00203         way_increment = Diff2D(-1, 1);
00204         border_increment.second = -border_increment.second;
00205         border_increment.third = -border_increment.third;
00206         valid_step_count = std::make_pair((yys - src_ul).x + 1, (src_lr - yys).y);
00207     }
00208     else if(down_to_much && right_to_much)
00209     {
00210         yys += kul;
00211         yk += klr;
00212         way_increment = Diff2D(1, 1);
00213         border_increment.first = -border_increment.first;
00214         border_increment.second = -border_increment.second;
00215         valid_step_count = std::make_pair((src_lr - yys).x, (src_lr - yys).y);
00216     }
00217     else if(top_to_much && left_to_much)
00218     {
00219         yys += klr;
00220         yk += kul;
00221         way_increment = Diff2D(-1, -1);
00222         border_increment.third = -border_increment.third;
00223         border_increment.fourth = -border_increment.fourth;
00224         valid_step_count = std::make_pair((yys - src_ul).x + 1, (yys - src_ul).y + 1);
00225     }
00226     else
00227     { //top_to_much && right_to_much
00228         yys += kul + Diff2D(0, kernel_height - 1);
00229         yk += kul + Diff2D(kernel_width - 1, 0);
00230         way_increment = Diff2D(1, -1);
00231         border_increment.first = -border_increment.first;
00232         border_increment.fourth = -border_increment.fourth;
00233         valid_step_count = std::make_pair((src_lr - yys).x, (yys - src_ul).y + 1);
00234     }
00235 
00236     int yy = 0, xx;
00237 
00238     //laeuft den zulässigen Bereich in y-Richtung durch
00239     for(; yy < valid_step_count.second; ++yy, yys.y += way_increment.y, yk.y -= way_increment.y )
00240     {
00241         SrcIterator xxs = yys;
00242         KernelIterator xk  = yk;
00243 
00244         //laeuft den zulässigen Bereich in x-Richtung durch
00245         for(xx = 0; xx < valid_step_count.first; ++xx, xxs.x += way_increment.x, xk.x -= way_increment.x)
00246         {
00247             sum += ak(xk) * src_acc(xxs);
00248         }
00249 
00250         //Nächstes ++xxs.x wuerde in unzulässigen Bereich
00251         //bringen => Sprung in zulaessigen Bereich
00252         xxs.x += border_increment.first;
00253 
00254         for( ; xx < kernel_width; ++xx, xxs.x += border_increment.third, xk.x -= way_increment.x )
00255         {
00256             sum += ak(xk) * src_acc(xxs);
00257         }
00258     }
00259 
00260     //Nächstes ++yys.y wuerde in unzulässigen Bereich
00261     //bringen => Sprung in zulaessigen Bereich
00262     yys.y += border_increment.second;
00263 
00264     for( ; yy < kernel_height; ++yy, yys.y += border_increment.third, yk.y -= way_increment.y)
00265     {
00266         SrcIterator xxs = yys;
00267         KernelIterator xk  = yk;
00268 
00269         for(xx=0; xx < valid_step_count.first; ++xx, xxs.x += way_increment.x, xk.x -= way_increment.x)
00270         {
00271             sum += ak(xk) * src_acc(xxs);
00272         }
00273 
00274         //Sprung in den zulaessigen Bereich
00275         xxs.x += border_increment.first;
00276 
00277         for( ; xx < kernel_width; ++xx, xxs.x += border_increment.third, xk.x -= way_increment.x )
00278         {
00279             sum += ak(xk) * src_acc(xxs);
00280         }
00281     }
00282 
00283     // store average in destination pixel
00284     dest_acc.set(DestTraits::fromRealPromote(sum), xd);
00285 
00286 }// end of internalPixelEvaluationByWrapReflectRepeat
00287 #endif /* #if 0 */
00288 
00289 
00290 template <class SrcIterator, class SrcAccessor,
00291           class KernelIterator, class KernelAccessor,
00292           class SumType>
00293 void
00294 internalPixelEvaluationByWrapReflectRepeat(SrcIterator xs, SrcAccessor src_acc,
00295     KernelIterator xk, KernelAccessor ak,
00296     int left, int right, int kleft, int kright,
00297     int borderskipx, int borderinc, SumType & sum)
00298 {
00299     SrcIterator xxs = xs + left;
00300     KernelIterator xxk  = xk - left;
00301 
00302     for(int xx = left; xx <= right; ++xx, ++xxs, --xxk)
00303     {
00304         sum += ak(xxk) * src_acc(xxs);
00305     }
00306 
00307     xxs = xs + left - borderskipx;
00308     xxk = xk - left + 1;
00309     for(int xx = left - 1; xx >= -kright; --xx, xxs -= borderinc, ++xxk)
00310     {
00311         sum += ak(xxk) * src_acc(xxs);
00312     }
00313 
00314     xxs = xs + right + borderskipx;
00315     xxk = xk - right - 1;
00316     for(int xx = right + 1; xx <= -kleft; ++xx, xxs += borderinc, --xxk)
00317     {
00318         sum += ak(xxk) * src_acc(xxs);
00319     }
00320 }
00321 
00322 
00323 /** \addtogroup StandardConvolution Two-dimensional convolution functions
00324 
00325 Perform 2D non-separable convolution, with and without ROI mask.
00326 
00327 These generic convolution functions implement
00328 the standard 2D convolution operation for images that fit
00329 into the required interface. Arbitrary ROI's are supported
00330 by the mask version of the algorithm.
00331 The functions need a suitable 2D kernel to operate.
00332 */
00333 //@{
00334 
00335 /** \brief Performs a 2 dimensional convolution of the source image using the given
00336     kernel.
00337 
00338     The KernelIterator must point to the center of the kernel, and
00339     the kernel's size is given by its upper left (x and y of distance <= 0) and
00340     lower right (distance >= 0) corners. The image must always be larger than the
00341     kernel. At those positions where the kernel does not completely fit
00342     into the image, the specified \ref BorderTreatmentMode is
00343     applied. You can choice between following BorderTreatmentModes:
00344     <ul>
00345     <li>BORDER_TREATMENT_CLIP</li>
00346     <li>BORDER_TREATMENT_AVOID</li>
00347     <li>BORDER_TREATMENT_WRAP</li>
00348     <li>BORDER_TREATMENT_REFLECT</li>
00349     <li>BORDER_TREATMENT_REPEAT</li>
00350     </ul><br>
00351     The images's pixel type (SrcAccessor::value_type) must be a
00352     linear space over the kernel's value_type (KernelAccessor::value_type),
00353     i.e. addition of source values, multiplication with kernel values,
00354     and NumericTraits must be defined.
00355     The kernel's value_type must be an algebraic field,
00356     i.e. the arithmetic operations (+, -, *, /) and NumericTraits must
00357     be defined.
00358 
00359     <b> Declarations:</b>
00360 
00361     pass arguments explicitly:
00362     \code
00363     namespace vigra {
00364         template <class SrcIterator, class SrcAccessor,
00365                   class DestIterator, class DestAccessor,
00366                   class KernelIterator, class KernelAccessor>
00367         void convolveImage(SrcIterator src_ul, SrcIterator src_lr, SrcAccessor src_acc,
00368                            DestIterator dest_ul, DestAccessor dest_acc,
00369                            KernelIterator ki, KernelAccessor ak,
00370                            Diff2D kul, Diff2D klr, BorderTreatmentMode border);
00371     }
00372     \endcode
00373 
00374 
00375     use argument objects in conjunction with \ref ArgumentObjectFactories :
00376     \code
00377     namespace vigra {
00378         template <class SrcIterator, class SrcAccessor,
00379                   class DestIterator, class DestAccessor,
00380                   class KernelIterator, class KernelAccessor>
00381         void convolveImage(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00382                            pair<DestIterator, DestAccessor> dest,
00383                            tuple5<KernelIterator, KernelAccessor, Diff2D, Diff2D,
00384                            BorderTreatmentMode> kernel);
00385     }
00386     \endcode
00387 
00388     <b> Usage:</b>
00389 
00390     <b>\#include</b> <<a href="stdconvolution_8hxx-source.html">vigra/stdconvolution.hxx</a>><br>
00391     Namespace: vigra
00392 
00393 
00394     \code
00395     vigra::FImage src(w,h), dest(w,h);
00396     ...
00397 
00398     // define horizontal Sobel filter
00399     vigra::Kernel2D<float> sobel;
00400 
00401     sobel.initExplicitly(Diff2D(-1,-1), Diff2D(1,1)) =  // upper left and lower right
00402                          0.125, 0.0, -0.125,
00403                          0.25,  0.0, -0.25,
00404                          0.125, 0.0, -0.125;
00405     sobel.setBorderTreatment(vigra::BORDER_TREATMENT_REFLECT);
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     KernelSumType norm = NumericTraits<KernelSumType>::zero();
00502     if(border == BORDER_TREATMENT_CLIP)
00503     {
00504         // calculate the sum of the kernel elements for renormalization
00505         KernelIterator yk  = ki + klr;
00506 
00507         // determine sum within kernel (= norm)
00508         for(int y = 0; y < kernel_height; ++y, --yk.y)
00509         {
00510             KernelIterator xk  = yk;
00511             for(int x = 0; x < kernel_width; ++x, --xk.x)
00512             {
00513                 norm += ak(xk);
00514             }
00515         }
00516         vigra_precondition(norm != NumericTraits<KernelSumType>::zero(),
00517             "convolveImage(): Cannot use BORDER_TREATMENT_CLIP with a DC-free kernel");
00518     }
00519 
00520     // create iterators for the interior part of the image (where the kernel always fits into the image)
00521     DestIterator yd = dest_ul + Diff2D(klr.x, klr.y);
00522     SrcIterator ys = src_ul + Diff2D(klr.x, klr.y);
00523     SrcIterator send = src_lr + Diff2D(kul.x, kul.y);
00524 
00525     // iterate over the interior part
00526     for(; ys.y < send.y; ++ys.y, ++yd.y)
00527     {
00528         // create x iterators
00529         DestIterator xd(yd);
00530         SrcIterator xs(ys);
00531 
00532         for(; xs.x < send.x; ++xs.x, ++xd.x)
00533         {
00534             // init the sum
00535             SumType sum = NumericTraits<SumType>::zero();
00536 
00537             SrcIterator yys = xs - klr;
00538             SrcIterator yyend = xs - kul;
00539             KernelIterator yk  = ki + klr;
00540 
00541             for(; yys.y <= yyend.y; ++yys.y, --yk.y)
00542             {
00543                 typename SrcIterator::row_iterator xxs = yys.rowIterator();
00544                 typename SrcIterator::row_iterator xxe = xxs + kernel_width;
00545                 typename KernelIterator::row_iterator xk  = yk.rowIterator();
00546 
00547                 for(; xxs < xxe; ++xxs, --xk)
00548                 {
00549                     sum += ak(xk) * src_acc(xxs);
00550                 }
00551             }
00552 
00553             // store convolution result in destination pixel
00554             dest_acc.set(DestTraits::fromRealPromote(sum), xd);
00555         }
00556     }
00557 
00558     if(border == BORDER_TREATMENT_AVOID)
00559         return; // skip processing near the border
00560 
00561     int interiorskip = w + kul.x - klr.x - 1;
00562     int borderskipx;
00563     int borderskipy;
00564     int borderinc;
00565     if(border == BORDER_TREATMENT_REPEAT)
00566     {
00567         borderskipx = 0;
00568         borderskipy = 0;
00569         borderinc = 0;
00570     }
00571     else if(border == BORDER_TREATMENT_REFLECT)
00572     {
00573         borderskipx = -1;
00574         borderskipy = -1;
00575         borderinc = -1;
00576     }
00577     else if(border == BORDER_TREATMENT_WRAP)
00578     {
00579         borderskipx = -w+1;
00580         borderskipy = -h+1;
00581         borderinc = 1;
00582     }
00583 
00584     // create iterators for the entire image
00585     yd = dest_ul;
00586     ys = src_ul;
00587 
00588     // work on entire image (but skip the already computed points in the loop)
00589     for(int y = 0; y < h; ++y, ++ys.y, ++yd.y)
00590     {
00591         int top    = std::max(static_cast<IntBiggest>(-klr.y),
00592                               static_cast<IntBiggest>(src_ul.y - ys.y));
00593         int bottom = std::min(static_cast<IntBiggest>(-kul.y),
00594                               static_cast<IntBiggest>(src_lr.y - ys.y - 1));
00595 
00596         // create x iterators
00597         DestIterator xd(yd);
00598         SrcIterator xs(ys);
00599 
00600         for(int x = 0; x < w; ++x, ++xs.x, ++xd.x)
00601         {
00602             // check if we are away from the border
00603             if(y >= klr.y && y < h+kul.y && x == klr.x)
00604             {
00605                 // yes => skip the already computed points
00606                 x += interiorskip;
00607                 xs.x += interiorskip;
00608                 xd.x += interiorskip;
00609                 continue;
00610             }
00611             if (border == BORDER_TREATMENT_CLIP)
00612             {
00613                 internalPixelEvaluationByClip(x, y, w, h, xs, src_acc, xd, dest_acc, ki, kul, klr, ak, norm);
00614             }
00615             else
00616             {
00617                 int left  = std::max(-klr.x, src_ul.x - xs.x);
00618                 int right = std::min(-kul.x, src_lr.x - xs.x - 1);
00619 
00620                 // init the sum
00621                 SumType sum = NumericTraits<SumType>::zero();
00622 
00623                 // create iterators for the part of the kernel that fits into the image
00624                 SrcIterator yys = xs + Size2D(0, top);
00625                 KernelIterator yk  = ki - Size2D(0, top);
00626 
00627                 int yy;
00628                 for(yy = top; yy <= bottom; ++yy, ++yys.y, --yk.y)
00629                 {
00630                     internalPixelEvaluationByWrapReflectRepeat(yys.rowIterator(), src_acc, yk.rowIterator(), ak,
00631                          left, right, kul.x, klr.x, borderskipx, borderinc, sum);
00632                 }
00633                 yys = xs + Size2D(0, top - borderskipy);
00634                 yk  = ki - Size2D(0, top - 1);
00635                 for(yy = top - 1; yy >= -klr.y; --yy, yys.y -= borderinc, ++yk.y)
00636                 {
00637                     internalPixelEvaluationByWrapReflectRepeat(yys.rowIterator(), src_acc, yk.rowIterator(), ak,
00638                          left, right, kul.x, klr.x, borderskipx, borderinc, sum);
00639                 }
00640                 yys = xs + Size2D(0, bottom + borderskipy);
00641                 yk  = ki - Size2D(0, bottom + 1);
00642                 for(yy = bottom + 1; yy <= -kul.y; ++yy, yys.y += borderinc, --yk.y)
00643                 {
00644                     internalPixelEvaluationByWrapReflectRepeat(yys.rowIterator(), src_acc, yk.rowIterator(), ak,
00645                          left, right, kul.x, klr.x, borderskipx, borderinc, sum);
00646                 }
00647 
00648                 // store convolution result in destination pixel
00649                 dest_acc.set(DestTraits::fromRealPromote(sum), xd);
00650 
00651 //                internalPixelEvaluationByWrapReflectRepeat(x, y, w, h, xs, src_acc, xd, dest_acc, ki, kul, klr, ak, border);
00652             }
00653         }
00654     }
00655 }
00656 
00657 
00658 template <class SrcIterator, class SrcAccessor,
00659           class DestIterator, class DestAccessor,
00660           class KernelIterator, class KernelAccessor>
00661 inline
00662 void convolveImage(
00663                    triple<SrcIterator, SrcIterator, SrcAccessor> src,
00664                    pair<DestIterator, DestAccessor> dest,
00665                    tuple5<KernelIterator, KernelAccessor, Diff2D, Diff2D,
00666                    BorderTreatmentMode> kernel)
00667 {
00668     convolveImage(src.first, src.second, src.third,
00669                   dest.first, dest.second,
00670                   kernel.first, kernel.second, kernel.third,
00671                   kernel.fourth, kernel.fifth);
00672 }
00673 
00674 
00675 /** \brief Performs a 2-dimensional normalized convolution, i.e. convolution with a mask image.
00676 
00677     This functions computes
00678     <a href ="http://homepages.inf.ed.ac.uk/rbf/CVonline/LOCAL_COPIES/PIRODDI1/NormConv/NormConv.html">normalized
00679     convolution</a> as defined in
00680     Knutsson, H. and Westin, C-F.: <i>Normalized and differential convolution:
00681     Methods for Interpolation and Filtering of incomplete and uncertain data</i>.
00682     Proc. of the IEEE Conf. on Computer Vision and Pattern Recognition, 1993, 515-523.
00683 
00684     The mask image must be binary and encodes which pixels of the original image
00685     are valid. It is used as follows:
00686     Only pixel under the mask are used in the calculations. Whenever a part of the
00687     kernel lies outside the mask, it is ignored, and the kernel is renormalized to its
00688     original norm (analogous to the CLIP \ref BorderTreatmentMode). Thus, a useful convolution
00689     result is computed whenever <i>at least one valid pixel is within the current window</i>
00690     Thus, destination pixels not under the mask still receive a value if they are <i>near</i>
00691     the mask. Therefore, this algorithm is useful as an interpolator of sparse input data.
00692     If you are only interested in the destination values under the mask, you can perform
00693     a subsequent \ref copyImageIf().
00694 
00695     The KernelIterator must point to the center of the kernel, and
00696     the kernel's size is given by its upper left (x and y of distance <= 0) and
00697     lower right (distance >= 0) corners. The image must always be larger than the
00698     kernel. At those positions where the kernel does not completely fit
00699     into the image, the specified \ref BorderTreatmentMode is
00700     applied. Only BORDER_TREATMENT_CLIP and BORDER_TREATMENT_AVOID are currently
00701     supported.
00702 
00703     The images's pixel type (SrcAccessor::value_type) must be a
00704     linear space over the kernel's value_type (KernelAccessor::value_type),
00705     i.e. addition of source values, multiplication with kernel values,
00706     and NumericTraits must be defined.
00707     The kernel's value_type must be an algebraic field,
00708     i.e. the arithmetic operations (+, -, *, /) and NumericTraits must
00709     be defined.
00710 
00711     <b> Declarations:</b>
00712 
00713     pass arguments explicitly:
00714     \code
00715     namespace vigra {
00716         template <class SrcIterator, class SrcAccessor,
00717                   class MaskIterator, class MaskAccessor,
00718                   class DestIterator, class DestAccessor,
00719                   class KernelIterator, class KernelAccessor>
00720         void
00721         normalizedConvolveImage(SrcIterator src_ul, SrcIterator src_lr, SrcAccessor src_acc,
00722                                 MaskIterator mul, MaskAccessor am,
00723                                 DestIterator dest_ul, DestAccessor dest_acc,
00724                                 KernelIterator ki, KernelAccessor ak,
00725                                 Diff2D kul, Diff2D klr, BorderTreatmentMode border);
00726     }
00727     \endcode
00728 
00729 
00730     use argument objects in conjunction with \ref ArgumentObjectFactories :
00731     \code
00732     namespace vigra {
00733         template <class SrcIterator, class SrcAccessor,
00734                   class MaskIterator, class MaskAccessor,
00735                   class DestIterator, class DestAccessor,
00736                   class KernelIterator, class KernelAccessor>
00737         void normalizedConvolveImage(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00738                                      pair<MaskIterator, MaskAccessor> mask,
00739                                      pair<DestIterator, DestAccessor> dest,
00740                                      tuple5<KernelIterator, KernelAccessor, Diff2D, Diff2D,
00741                                      BorderTreatmentMode> kernel);
00742     }
00743     \endcode
00744 
00745     <b> Usage:</b>
00746 
00747     <b>\#include</b> <<a href="stdconvolution_8hxx-source.html">vigra/stdconvolution.hxx</a>><br>
00748     Namespace: vigra
00749 
00750 
00751     \code
00752     vigra::FImage src(w,h), dest(w,h);
00753     vigra::CImage mask(w,h);
00754     ...
00755 
00756     // define 3x3 binomial filter
00757     vigra::Kernel2D<float> binom;
00758 
00759     binom.initExplicitly(Diff2D(-1,-1), Diff2D(1,1)) =   // upper left and lower right
00760                          0.0625, 0.125, 0.0625,
00761                          0.125,  0.25,  0.125,
00762                          0.0625, 0.125, 0.0625;
00763 
00764     vigra::normalizedConvolveImage(srcImageRange(src), maskImage(mask), destImage(dest), kernel2d(binom));
00765     \endcode
00766 
00767     <b> Required Interface:</b>
00768 
00769     \code
00770     ImageIterator src_ul, src_lr;
00771     ImageIterator mul;
00772     ImageIterator dest_ul;
00773     ImageIterator ik;
00774 
00775     SrcAccessor src_accessor;
00776     MaskAccessor mask_accessor;
00777     DestAccessor dest_accessor;
00778     KernelAccessor kernel_accessor;
00779 
00780     NumericTraits<SrcAccessor::value_type>::RealPromote s = src_accessor(src_ul);
00781 
00782     s = s + s;
00783     s = kernel_accessor(ik) * s;
00784     s -= s;
00785 
00786     if(mask_accessor(mul)) ...;
00787 
00788     dest_accessor.set(
00789     NumericTraits<DestAccessor::value_type>::fromRealPromote(s), dest_ul);
00790 
00791     NumericTraits<KernelAccessor::value_type>::RealPromote k = kernel_accessor(ik);
00792 
00793     k += k;
00794     k -= k;
00795     k = k / k;
00796 
00797     \endcode
00798 
00799     <b> Preconditions:</b>
00800 
00801     \code
00802     kul.x <= 0
00803     kul.y <= 0
00804     klr.x >= 0
00805     klr.y >= 0
00806     src_lr.x - src_ul.x >= klr.x + kul.x + 1
00807     src_lr.y - src_ul.y >= klr.y + kul.y + 1
00808     border == BORDER_TREATMENT_CLIP || border == BORDER_TREATMENT_AVOID
00809     \endcode
00810 
00811     Sum of kernel elements must be != 0.
00812 
00813 */
00814 doxygen_overloaded_function(template <...> void normalizedConvolveImage)
00815 
00816 template <class SrcIterator, class SrcAccessor,
00817           class DestIterator, class DestAccessor,
00818           class MaskIterator, class MaskAccessor,
00819           class KernelIterator, class KernelAccessor>
00820 void
00821 normalizedConvolveImage(SrcIterator src_ul, SrcIterator src_lr, SrcAccessor src_acc,
00822                         MaskIterator mul, MaskAccessor am,
00823                         DestIterator dest_ul, DestAccessor dest_acc,
00824                         KernelIterator ki, KernelAccessor ak,
00825                         Diff2D kul, Diff2D klr, BorderTreatmentMode border)
00826 {
00827     vigra_precondition((border == BORDER_TREATMENT_CLIP  ||
00828                         border == BORDER_TREATMENT_AVOID),
00829                        "normalizedConvolveImage(): "
00830                        "Border treatment must be BORDER_TREATMENT_CLIP or BORDER_TREATMENT_AVOID.");
00831 
00832     vigra_precondition(kul.x <= 0 && kul.y <= 0,
00833                        "normalizedConvolveImage(): left borders must be <= 0.");
00834     vigra_precondition(klr.x >= 0 && klr.y >= 0,
00835                        "normalizedConvolveImage(): right borders must be >= 0.");
00836 
00837     // use traits to determine SumType as to prevent possible overflow
00838     typedef typename
00839         NumericTraits<typename SrcAccessor::value_type>::RealPromote SumType;
00840     typedef typename
00841         NumericTraits<typename KernelAccessor::value_type>::RealPromote KSumType;
00842     typedef
00843         NumericTraits<typename DestAccessor::value_type> DestTraits;
00844 
00845     // calculate width and height of the image
00846     int w = src_lr.x - src_ul.x;
00847     int h = src_lr.y - src_ul.y;
00848     int kernel_width = klr.x - kul.x + 1;
00849     int kernel_height = klr.y - kul.y + 1;
00850 
00851     int x,y;
00852     int ystart = (border == BORDER_TREATMENT_AVOID) ?  klr.y : 0;
00853     int yend   = (border == BORDER_TREATMENT_AVOID) ? h+kul.y : h;
00854     int xstart = (border == BORDER_TREATMENT_AVOID) ?  klr.x : 0;
00855     int xend   = (border == BORDER_TREATMENT_AVOID) ? w+kul.x : w;
00856 
00857     // create y iterators
00858     DestIterator yd = dest_ul + Diff2D(xstart, ystart);
00859     SrcIterator ys = src_ul + Diff2D(xstart, ystart);
00860     MaskIterator ym = mul + Diff2D(xstart, ystart);
00861 
00862     KSumType norm = ak(ki);
00863     int xx, yy;
00864     KernelIterator yk  = ki + klr;
00865     for(yy=0; yy<kernel_height; ++yy, --yk.y)
00866     {
00867         KernelIterator xk  = yk;
00868 
00869         for(xx=0; xx<kernel_width; ++xx, --xk.x)
00870         {
00871             norm += ak(xk);
00872         }
00873     }
00874     norm -= ak(ki);
00875 
00876 
00877     for(y=ystart; y < yend; ++y, ++ys.y, ++yd.y, ++ym.y)
00878     {
00879         // create x iterators
00880         DestIterator xd(yd);
00881         SrcIterator xs(ys);
00882         MaskIterator xm(ym);
00883 
00884         for(x=xstart; x < xend; ++x, ++xs.x, ++xd.x, ++xm.x)
00885         {
00886             // how much of the kernel fits into the image ?
00887             int x0, y0, x1, y1;
00888 
00889             y0 = (y<klr.y) ? -y : -klr.y;
00890             y1 = (h-y-1<-kul.y) ? h-y-1 : -kul.y;
00891             x0 = (x<klr.x) ? -x : -klr.x;
00892             x1 = (w-x-1<-kul.x) ? w-x-1 : -kul.x;
00893 
00894             bool first = true;
00895             // init the sum
00896             SumType sum;
00897             KSumType ksum;
00898 
00899             SrcIterator yys = xs + Diff2D(x0, y0);
00900             MaskIterator yym = xm + Diff2D(x0, y0);
00901             KernelIterator yk  = ki - Diff2D(x0, y0);
00902 
00903             int xx, kernel_width, kernel_height;
00904             kernel_width = x1 - x0 + 1;
00905             kernel_height = y1 - y0 + 1;
00906             for(yy=0; yy<kernel_height; ++yy, ++yys.y, --yk.y, ++yym.y)
00907             {
00908                 typename SrcIterator::row_iterator xxs = yys.rowIterator();
00909                 typename SrcIterator::row_iterator xxend = xxs + kernel_width;
00910                 typename MaskIterator::row_iterator xxm = yym.rowIterator();
00911                 typename KernelIterator::row_iterator xk  = yk.rowIterator();
00912 
00913                 for(xx=0; xxs < xxend; ++xxs, --xk, ++xxm)
00914                 {
00915                     if(!am(xxm)) continue;
00916 
00917                     if(first)
00918                     {
00919                         sum = ak(xk) * src_acc(xxs);
00920                         ksum = ak(xk);
00921                         first = false;
00922                     }
00923                     else
00924                     {
00925                         sum += ak(xk) * src_acc(xxs);
00926                         ksum += ak(xk);
00927                     }
00928                 }
00929             }
00930             // store average in destination pixel
00931             if(!first &&
00932                ksum != NumericTraits<KSumType>::zero())
00933             {
00934                 dest_acc.set(DestTraits::fromRealPromote((norm / ksum) * sum), xd);
00935             }
00936         }
00937     }
00938 }
00939 
00940 
00941 template <class SrcIterator, class SrcAccessor,
00942           class DestIterator, class DestAccessor,
00943           class MaskIterator, class MaskAccessor,
00944           class KernelIterator, class KernelAccessor>
00945 inline
00946 void normalizedConvolveImage(
00947                            triple<SrcIterator, SrcIterator, SrcAccessor> src,
00948                            pair<MaskIterator, MaskAccessor> mask,
00949                            pair<DestIterator, DestAccessor> dest,
00950                            tuple5<KernelIterator, KernelAccessor, Diff2D, Diff2D,
00951                            BorderTreatmentMode> kernel)
00952 {
00953     normalizedConvolveImage(src.first, src.second, src.third,
00954                             mask.first, mask.second,
00955                             dest.first, dest.second,
00956                             kernel.first, kernel.second, kernel.third,
00957                             kernel.fourth, kernel.fifth);
00958 }
00959 
00960 /** \brief Deprecated name of 2-dimensional normalized convolution, i.e. convolution with a mask image.
00961 
00962     See \ref normalizedConvolveImage() for documentation.
00963 
00964     <b> Declarations:</b>
00965 
00966     pass arguments explicitly:
00967     \code
00968     namespace vigra {
00969         template <class SrcIterator, class SrcAccessor,
00970                   class MaskIterator, class MaskAccessor,
00971                   class DestIterator, class DestAccessor,
00972                   class KernelIterator, class KernelAccessor>
00973         void
00974         convolveImageWithMask(SrcIterator src_ul, SrcIterator src_lr, SrcAccessor src_acc,
00975                               MaskIterator mul, MaskAccessor am,
00976                               DestIterator dest_ul, DestAccessor dest_acc,
00977                               KernelIterator ki, KernelAccessor ak,
00978                               Diff2D kul, Diff2D klr, BorderTreatmentMode border);
00979     }
00980     \endcode
00981 
00982 
00983     use argument objects in conjunction with \ref ArgumentObjectFactories :
00984     \code
00985     namespace vigra {
00986         template <class SrcIterator, class SrcAccessor,
00987                   class MaskIterator, class MaskAccessor,
00988                   class DestIterator, class DestAccessor,
00989                   class KernelIterator, class KernelAccessor>
00990         void convolveImageWithMask(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00991                                    pair<MaskIterator, MaskAccessor> mask,
00992                                    pair<DestIterator, DestAccessor> dest,
00993                                    tuple5<KernelIterator, KernelAccessor, Diff2D, Diff2D,
00994                                    BorderTreatmentMode> kernel);
00995     }
00996     \endcode
00997 */
00998 doxygen_overloaded_function(template <...> void convolveImageWithMask)
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 (ullrich.koethe@iwr.uni-heidelberg.de)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
VIGRA 1.6.0 (5 Nov 2009)