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

vigra/stdconvolution.hxx |

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