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

details vigra/stdconvolution.hxx VIGRA

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

© Ullrich Köthe (koethe@informatik.uni-hamburg.de)
Cognitive Systems Group, University of Hamburg, Germany

html generated using doxygen and Python
VIGRA 1.4.0 (21 Dec 2005)