vigra - vigra/separableconvolution.hxx Source File

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

vigra/separableconvolution.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_SEPARABLECONVOLUTION_HXX
    00038 #define VIGRA_SEPARABLECONVOLUTION_HXX
    00039
    00040 #include <cmath>
    00041 #include "utilities.hxx"
    00042 #include "numerictraits.hxx"
    00043 #include "imageiteratoradapter.hxx"
    00044 #include "bordertreatment.hxx"
    00045 #include "gaussians.hxx"
    00046 #include "array_vector.hxx"
    00047
    00048 namespace vigra {
    00049
    00050 /********************************************************/
    00051 /*                                                      */
    00052 /*            internalConvolveLineOptimistic            */
    00053 /*                                                      */
    00054 /********************************************************/
    00055
    00056 // This function assumes that the input array is actually larger than
    00057 // the range [is, iend), so that it can safely access values outside
    00058 // this range. This is useful if (1) we work on a small ROI, or
    00059 // (2) we enlarge the input by copying with border treatment.
    00060 template <class SrcIterator, class SrcAccessor,
    00061           class DestIterator, class DestAccessor,
    00062           class KernelIterator, class KernelAccessor>
    00063 void internalConvolveLineOptimistic(SrcIterator is, SrcIterator iend, SrcAccessor sa,
    00064                                     DestIterator id, DestAccessor da,
    00065                                     KernelIterator kernel, KernelAccessor ka,
    00066                                     int kleft, int kright)
    00067 {
    00068     typedef typename PromoteTraits<
    00069             typename SrcAccessor::value_type,
    00070             typename KernelAccessor::value_type>::Promote SumType;
    00071
    00072     int w = std::distance( is, iend );
    00073     int kw = kright - kleft + 1;
    00074     for(int x=0; x<w; ++x, ++is, ++id)
    00075     {
    00076         SrcIterator iss = is + (-kright);
    00077         KernelIterator ik = kernel + kright;
    00078         SumType sum = NumericTraits<SumType>::zero();
    00079
    00080         for(int k = 0; k < kw; ++k, --ik, ++iss)
    00081         {
    00082             sum += ka(ik) * sa(iss);
    00083         }
    00084
    00085         da.set(detail::RequiresExplicitCast<typename
    00086                       DestAccessor::value_type>::cast(sum), id);
    00087     }
    00088 }
    00089
    00090 namespace detail {
    00091
    00092 // dest array must have size = stop - start + kright - kleft
    00093 template <class SrcIterator, class SrcAccessor,
    00094           class DestIterator, class DestAccessor>
    00095 void
    00096 copyLineWithBorderTreatment(SrcIterator is, SrcIterator iend, SrcAccessor sa,
    00097                             DestIterator id, DestAccessor da,
    00098                             int start, int stop,
    00099                             int kleft, int kright,
    00100                             BorderTreatmentMode borderTreatment)
    00101 {
    00102     int w = std::distance( is, iend );
    00103     int leftBorder = start - kright;
    00104     int rightBorder = stop - kleft;
    00105     int copyEnd = std::min(w, rightBorder);
    00106
    00107     if(leftBorder < 0)
    00108     {
    00109         switch(borderTreatment)
    00110         {
    00111             case BORDER_TREATMENT_WRAP:
    00112             {
    00113                 for(; leftBorder<0; ++leftBorder, ++id)
    00114                     da.set(sa(iend, leftBorder), id);
    00115                 break;
    00116             }
    00117             case BORDER_TREATMENT_AVOID:
    00118             {
    00119                 // nothing to do
    00120                 break;
    00121             }
    00122             case BORDER_TREATMENT_REFLECT:
    00123             {
    00124                 for(; leftBorder<0; ++leftBorder, ++id)
    00125                     da.set(sa(is, -leftBorder), id);
    00126                 break;
    00127             }
    00128             case BORDER_TREATMENT_REPEAT:
    00129             {
    00130                 for(; leftBorder<0; ++leftBorder, ++id)
    00131                     da.set(sa(is), id);
    00132                 break;
    00133             }
    00134             case BORDER_TREATMENT_CLIP:
    00135             {
    00136                 vigra_precondition(false,
    00137                              "copyLineWithBorderTreatment() internal error: not applicable to BORDER_TREATMENT_CLIP.");
    00138                 break;
    00139             }
    00140             case BORDER_TREATMENT_ZEROPAD:
    00141             {
    00142                 for(; leftBorder<0; ++leftBorder, ++id)
    00143                     da.set(NumericTraits<typename DestAccessor::value_type>::zero(), id);
    00144                 break;
    00145             }
    00146             default:
    00147             {
    00148                 vigra_precondition(false,
    00149                              "copyLineWithBorderTreatment(): Unknown border treatment mode.");
    00150             }
    00151         }
    00152     }
    00153
    00154     SrcIterator iss = is + leftBorder;
    00155     vigra_invariant( leftBorder < copyEnd,
    00156         "copyLineWithBorderTreatment(): assertion failed.");
    00157     for(; leftBorder<copyEnd; ++leftBorder, ++id, ++iss)
    00158         da.set(sa(iss), id);
    00159
    00160     if(copyEnd < rightBorder)
    00161     {
    00162         switch(borderTreatment)
    00163         {
    00164             case BORDER_TREATMENT_WRAP:
    00165             {
    00166                 for(; copyEnd<rightBorder; ++copyEnd, ++id, ++is)
    00167                     da.set(sa(is), id);
    00168                 break;
    00169             }
    00170             case BORDER_TREATMENT_AVOID:
    00171             {
    00172                 // nothing to do
    00173                 break;
    00174             }
    00175             case BORDER_TREATMENT_REFLECT:
    00176             {
    00177                 iss -= 2;
    00178                 for(; copyEnd<rightBorder; ++copyEnd, ++id, --iss)
    00179                     da.set(sa(iss), id);
    00180                 break;
    00181             }
    00182             case BORDER_TREATMENT_REPEAT:
    00183             {
    00184                 --iss;
    00185                 for(; copyEnd<rightBorder; ++copyEnd, ++id)
    00186                     da.set(sa(iss), id);
    00187                 break;
    00188             }
    00189             case BORDER_TREATMENT_CLIP:
    00190             {
    00191                 vigra_precondition(false,
    00192                              "copyLineWithBorderTreatment() internal error: not applicable to BORDER_TREATMENT_CLIP.");
    00193                 break;
    00194             }
    00195             case BORDER_TREATMENT_ZEROPAD:
    00196             {
    00197                 for(; copyEnd<rightBorder; ++copyEnd, ++id)
    00198                     da.set(NumericTraits<typename DestAccessor::value_type>::zero(), id);
    00199                 break;
    00200             }
    00201             default:
    00202             {
    00203                 vigra_precondition(false,
    00204                              "copyLineWithBorderTreatment(): Unknown border treatment mode.");
    00205             }
    00206         }
    00207     }
    00208 }
    00209
    00210 } // namespace detail
    00211
    00212 /********************************************************/
    00213 /*                                                      */
    00214 /*                internalConvolveLineWrap              */
    00215 /*                                                      */
    00216 /********************************************************/
    00217
    00218 template <class SrcIterator, class SrcAccessor,
    00219           class DestIterator, class DestAccessor,
    00220           class KernelIterator, class KernelAccessor>
    00221 void internalConvolveLineWrap(SrcIterator is, SrcIterator iend, SrcAccessor sa,
    00222                               DestIterator id, DestAccessor da,
    00223                               KernelIterator kernel, KernelAccessor ka,
    00224                               int kleft, int kright,
    00225                               int start = 0, int stop = 0)
    00226 {
    00227     int w = std::distance( is, iend );
    00228
    00229     typedef typename PromoteTraits<
    00230             typename SrcAccessor::value_type,
    00231             typename KernelAccessor::value_type>::Promote SumType;
    00232
    00233     SrcIterator ibegin = is;
    00234
    00235     if(stop == 0)
    00236         stop = w;
    00237     is += start;
    00238
    00239     for(int x=start; x<stop; ++x, ++is, ++id)
    00240     {
    00241         KernelIterator ik = kernel + kright;
    00242         SumType sum = NumericTraits<SumType>::zero();
    00243
    00244         if(x < kright)
    00245         {
    00246             int x0 = x - kright;
    00247             SrcIterator iss = iend + x0;
    00248
    00249             for(; x0; ++x0, --ik, ++iss)
    00250             {
    00251                 sum += ka(ik) * sa(iss);
    00252             }
    00253
    00254             iss = ibegin;
    00255             if(w-x <= -kleft)
    00256             {
    00257                 SrcIterator isend = iend;
    00258                 for(; iss != isend ; --ik, ++iss)
    00259                 {
    00260                     sum += ka(ik) * sa(iss);
    00261                 }
    00262
    00263                 int x0 = -kleft - w + x + 1;
    00264                 iss = ibegin;
    00265
    00266                 for(; x0; --x0, --ik, ++iss)
    00267                 {
    00268                     sum += ka(ik) * sa(iss);
    00269                 }
    00270             }
    00271             else
    00272             {
    00273                 SrcIterator isend = is + (1 - kleft);
    00274                 for(; iss != isend ; --ik, ++iss)
    00275                 {
    00276                     sum += ka(ik) * sa(iss);
    00277                 }
    00278             }
    00279         }
    00280         else if(w-x <= -kleft)
    00281         {
    00282             SrcIterator iss = is + (-kright);
    00283             SrcIterator isend = iend;
    00284             for(; iss != isend ; --ik, ++iss)
    00285             {
    00286                 sum += ka(ik) * sa(iss);
    00287             }
    00288
    00289             int x0 = -kleft - w + x + 1;
    00290             iss = ibegin;
    00291
    00292             for(; x0; --x0, --ik, ++iss)
    00293             {
    00294                 sum += ka(ik) * sa(iss);
    00295             }
    00296         }
    00297         else
    00298         {
    00299             SrcIterator iss = is - kright;
    00300             SrcIterator isend = is + (1 - kleft);
    00301             for(; iss != isend ; --ik, ++iss)
    00302             {
    00303                 sum += ka(ik) * sa(iss);
    00304             }
    00305         }
    00306
    00307         da.set(detail::RequiresExplicitCast<typename
    00308                       DestAccessor::value_type>::cast(sum), id);
    00309     }
    00310 }
    00311
    00312 /********************************************************/
    00313 /*                                                      */
    00314 /*                internalConvolveLineClip              */
    00315 /*                                                      */
    00316 /********************************************************/
    00317
    00318 template <class SrcIterator, class SrcAccessor,
    00319           class DestIterator, class DestAccessor,
    00320           class KernelIterator, class KernelAccessor,
    00321           class Norm>
    00322 void internalConvolveLineClip(SrcIterator is, SrcIterator iend, SrcAccessor sa,
    00323                               DestIterator id, DestAccessor da,
    00324                               KernelIterator kernel, KernelAccessor ka,
    00325                               int kleft, int kright, Norm norm,
    00326                               int start = 0, int stop = 0)
    00327 {
    00328     int w = std::distance( is, iend );
    00329
    00330     typedef typename PromoteTraits<
    00331             typename SrcAccessor::value_type,
    00332             typename KernelAccessor::value_type>::Promote SumType;
    00333
    00334     SrcIterator ibegin = is;
    00335
    00336     if(stop == 0)
    00337         stop = w;
    00338     is += start;
    00339
    00340     for(int x=start; x<stop; ++x, ++is, ++id)
    00341     {
    00342         KernelIterator ik = kernel + kright;
    00343         SumType sum = NumericTraits<SumType>::zero();
    00344
    00345         if(x < kright)
    00346         {
    00347             int x0 = x - kright;
    00348             Norm clipped = NumericTraits<Norm>::zero();
    00349
    00350             for(; x0; ++x0, --ik)
    00351             {
    00352                 clipped += ka(ik);
    00353             }
    00354
    00355             SrcIterator iss = ibegin;
    00356             if(w-x <= -kleft)
    00357             {
    00358                 SrcIterator isend = iend;
    00359                 for(; iss != isend ; --ik, ++iss)
    00360                 {
    00361                     sum += ka(ik) * sa(iss);
    00362                 }
    00363
    00364                 int x0 = -kleft - w + x + 1;
    00365
    00366                 for(; x0; --x0, --ik)
    00367                 {
    00368                     clipped += ka(ik);
    00369                 }
    00370             }
    00371             else
    00372             {
    00373                 SrcIterator isend = is + (1 - kleft);
    00374                 for(; iss != isend ; --ik, ++iss)
    00375                 {
    00376                     sum += ka(ik) * sa(iss);
    00377                 }
    00378             }
    00379
    00380             sum = norm / (norm - clipped) * sum;
    00381         }
    00382         else if(w-x <= -kleft)
    00383         {
    00384             SrcIterator iss = is + (-kright);
    00385             SrcIterator isend = iend;
    00386             for(; iss != isend ; --ik, ++iss)
    00387             {
    00388                 sum += ka(ik) * sa(iss);
    00389             }
    00390
    00391             Norm clipped = NumericTraits<Norm>::zero();
    00392
    00393             int x0 = -kleft - w + x + 1;
    00394
    00395             for(; x0; --x0, --ik)
    00396             {
    00397                 clipped += ka(ik);
    00398             }
    00399
    00400             sum = norm / (norm - clipped) * sum;
    00401         }
    00402         else
    00403         {
    00404             SrcIterator iss = is + (-kright);
    00405             SrcIterator isend = is + (1 - kleft);
    00406             for(; iss != isend ; --ik, ++iss)
    00407             {
    00408                 sum += ka(ik) * sa(iss);
    00409             }
    00410         }
    00411
    00412         da.set(detail::RequiresExplicitCast<typename
    00413                       DestAccessor::value_type>::cast(sum), id);
    00414     }
    00415 }
    00416
    00417 /********************************************************/
    00418 /*                                                      */
    00419 /*             internalConvolveLineZeropad              */
    00420 /*                                                      */
    00421 /********************************************************/
    00422
    00423 template <class SrcIterator, class SrcAccessor,
    00424           class DestIterator, class DestAccessor,
    00425           class KernelIterator, class KernelAccessor>
    00426 void internalConvolveLineZeropad(SrcIterator is, SrcIterator iend, SrcAccessor sa,
    00427                                  DestIterator id, DestAccessor da,
    00428                                  KernelIterator kernel, KernelAccessor ka,
    00429                                  int kleft, int kright,
    00430                                  int start = 0, int stop = 0)
    00431 {
    00432     int w = std::distance( is, iend );
    00433
    00434     typedef typename PromoteTraits<
    00435             typename SrcAccessor::value_type,
    00436             typename KernelAccessor::value_type>::Promote SumType;
    00437
    00438     SrcIterator ibegin = is;
    00439
    00440     if(stop == 0)
    00441         stop = w;
    00442     is += start;
    00443
    00444     for(int x=start; x<stop; ++x, ++is, ++id)
    00445     {
    00446         SumType sum = NumericTraits<SumType>::zero();
    00447
    00448         if(x < kright)
    00449         {
    00450             KernelIterator ik = kernel + x;
    00451             SrcIterator iss = ibegin;
    00452
    00453             if(w-x <= -kleft)
    00454             {
    00455                 SrcIterator isend = iend;
    00456                 for(; iss != isend ; --ik, ++iss)
    00457                 {
    00458                     sum += ka(ik) * sa(iss);
    00459                 }
    00460             }
    00461             else
    00462             {
    00463                 SrcIterator isend = is + (1 - kleft);
    00464                 for(; iss != isend ; --ik, ++iss)
    00465                 {
    00466                     sum += ka(ik) * sa(iss);
    00467                 }
    00468             }
    00469         }
    00470         else if(w-x <= -kleft)
    00471         {
    00472             KernelIterator ik = kernel + kright;
    00473             SrcIterator iss = is + (-kright);
    00474             SrcIterator isend = iend;
    00475             for(; iss != isend ; --ik, ++iss)
    00476             {
    00477                 sum += ka(ik) * sa(iss);
    00478             }
    00479         }
    00480         else
    00481         {
    00482             KernelIterator ik = kernel + kright;
    00483             SrcIterator iss = is + (-kright);
    00484             SrcIterator isend = is + (1 - kleft);
    00485             for(; iss != isend ; --ik, ++iss)
    00486             {
    00487                 sum += ka(ik) * sa(iss);
    00488             }
    00489         }
    00490
    00491         da.set(detail::RequiresExplicitCast<typename
    00492                       DestAccessor::value_type>::cast(sum), id);
    00493     }
    00494 }
    00495
    00496 /********************************************************/
    00497 /*                                                      */
    00498 /*             internalConvolveLineReflect              */
    00499 /*                                                      */
    00500 /********************************************************/
    00501
    00502 template <class SrcIterator, class SrcAccessor,
    00503           class DestIterator, class DestAccessor,
    00504           class KernelIterator, class KernelAccessor>
    00505 void internalConvolveLineReflect(SrcIterator is, SrcIterator iend, SrcAccessor sa,
    00506                               DestIterator id, DestAccessor da,
    00507                               KernelIterator kernel, KernelAccessor ka,
    00508                               int kleft, int kright,
    00509                               int start = 0, int stop = 0)
    00510 {
    00511     int w = std::distance( is, iend );
    00512
    00513     typedef typename PromoteTraits<
    00514             typename SrcAccessor::value_type,
    00515             typename KernelAccessor::value_type>::Promote SumType;
    00516
    00517     SrcIterator ibegin = is;
    00518
    00519     if(stop == 0)
    00520         stop = w;
    00521     is += start;
    00522
    00523     for(int x=start; x<stop; ++x, ++is, ++id)
    00524     {
    00525         KernelIterator ik = kernel + kright;
    00526         SumType sum = NumericTraits<SumType>::zero();
    00527
    00528         if(x < kright)
    00529         {
    00530             int x0 = x - kright;
    00531             SrcIterator iss = ibegin - x0;
    00532
    00533             for(; x0; ++x0, --ik, --iss)
    00534             {
    00535                 sum += ka(ik) * sa(iss);
    00536             }
    00537
    00538             if(w-x <= -kleft)
    00539             {
    00540                 SrcIterator isend = iend;
    00541                 for(; iss != isend ; --ik, ++iss)
    00542                 {
    00543                     sum += ka(ik) * sa(iss);
    00544                 }
    00545
    00546                 int x0 = -kleft - w + x + 1;
    00547                 iss = iend - 2;
    00548
    00549                 for(; x0; --x0, --ik, --iss)
    00550                 {
    00551                     sum += ka(ik) * sa(iss);
    00552                 }
    00553             }
    00554             else
    00555             {
    00556                 SrcIterator isend = is + (1 - kleft);
    00557                 for(; iss != isend ; --ik, ++iss)
    00558                 {
    00559                     sum += ka(ik) * sa(iss);
    00560                 }
    00561             }
    00562         }
    00563         else if(w-x <= -kleft)
    00564         {
    00565             SrcIterator iss = is + (-kright);
    00566             SrcIterator isend = iend;
    00567             for(; iss != isend ; --ik, ++iss)
    00568             {
    00569                 sum += ka(ik) * sa(iss);
    00570             }
    00571
    00572             int x0 = -kleft - w + x + 1;
    00573             iss = iend - 2;
    00574
    00575             for(; x0; --x0, --ik, --iss)
    00576             {
    00577                 sum += ka(ik) * sa(iss);
    00578             }
    00579         }
    00580         else
    00581         {
    00582             SrcIterator iss = is + (-kright);
    00583             SrcIterator isend = is + (1 - kleft);
    00584             for(; iss != isend ; --ik, ++iss)
    00585             {
    00586                 sum += ka(ik) * sa(iss);
    00587             }
    00588         }
    00589
    00590         da.set(detail::RequiresExplicitCast<typename
    00591                       DestAccessor::value_type>::cast(sum), id);
    00592     }
    00593 }
    00594
    00595 /********************************************************/
    00596 /*                                                      */
    00597 /*             internalConvolveLineRepeat               */
    00598 /*                                                      */
    00599 /********************************************************/
    00600
    00601 template <class SrcIterator, class SrcAccessor,
    00602           class DestIterator, class DestAccessor,
    00603           class KernelIterator, class KernelAccessor>
    00604 void internalConvolveLineRepeat(SrcIterator is, SrcIterator iend, SrcAccessor sa,
    00605                               DestIterator id, DestAccessor da,
    00606                               KernelIterator kernel, KernelAccessor ka,
    00607                               int kleft, int kright,
    00608                               int start = 0, int stop = 0)
    00609 {
    00610     int w = std::distance( is, iend );
    00611
    00612     typedef typename PromoteTraits<
    00613             typename SrcAccessor::value_type,
    00614             typename KernelAccessor::value_type>::Promote SumType;
    00615
    00616     SrcIterator ibegin = is;
    00617
    00618     if(stop == 0)
    00619         stop = w;
    00620     is += start;
    00621
    00622     for(int x=start; x<stop; ++x, ++is, ++id)
    00623     {
    00624         KernelIterator ik = kernel + kright;
    00625         SumType sum = NumericTraits<SumType>::zero();
    00626
    00627         if(x < kright)
    00628         {
    00629             int x0 = x - kright;
    00630             SrcIterator iss = ibegin;
    00631
    00632             for(; x0; ++x0, --ik)
    00633             {
    00634                 sum += ka(ik) * sa(iss);
    00635             }
    00636
    00637             if(w-x <= -kleft)
    00638             {
    00639                 SrcIterator isend = iend;
    00640                 for(; iss != isend ; --ik, ++iss)
    00641                 {
    00642                     sum += ka(ik) * sa(iss);
    00643                 }
    00644
    00645                 int x0 = -kleft - w + x + 1;
    00646                 iss = iend - 1;
    00647
    00648                 for(; x0; --x0, --ik)
    00649                 {
    00650                     sum += ka(ik) * sa(iss);
    00651                 }
    00652             }
    00653             else
    00654             {
    00655                 SrcIterator isend = is + (1 - kleft);
    00656                 for(; iss != isend ; --ik, ++iss)
    00657                 {
    00658                     sum += ka(ik) * sa(iss);
    00659                 }
    00660             }
    00661         }
    00662         else if(w-x <= -kleft)
    00663         {
    00664             SrcIterator iss = is + (-kright);
    00665             SrcIterator isend = iend;
    00666             for(; iss != isend ; --ik, ++iss)
    00667             {
    00668                 sum += ka(ik) * sa(iss);
    00669             }
    00670
    00671             int x0 = -kleft - w + x + 1;
    00672             iss = iend - 1;
    00673
    00674             for(; x0; --x0, --ik)
    00675             {
    00676                 sum += ka(ik) * sa(iss);
    00677             }
    00678         }
    00679         else
    00680         {
    00681             SrcIterator iss = is + (-kright);
    00682             SrcIterator isend = is + (1 - kleft);
    00683             for(; iss != isend ; --ik, ++iss)
    00684             {
    00685                 sum += ka(ik) * sa(iss);
    00686             }
    00687         }
    00688
    00689         da.set(detail::RequiresExplicitCast<typename
    00690                       DestAccessor::value_type>::cast(sum), id);
    00691     }
    00692 }
    00693
    00694 /********************************************************/
    00695 /*                                                      */
    00696 /*              internalConvolveLineAvoid               */
    00697 /*                                                      */
    00698 /********************************************************/
    00699
    00700 template <class SrcIterator, class SrcAccessor,
    00701           class DestIterator, class DestAccessor,
    00702           class KernelIterator, class KernelAccessor>
    00703 void internalConvolveLineAvoid(SrcIterator is, SrcIterator iend, SrcAccessor sa,
    00704                               DestIterator id, DestAccessor da,
    00705                               KernelIterator kernel, KernelAccessor ka,
    00706                               int kleft, int kright,
    00707                               int start = 0, int stop = 0)
    00708 {
    00709     int w = std::distance( is, iend );
    00710     if(start < stop) // we got a valid subrange
    00711     {
    00712         if(w + kleft < stop)
    00713             stop = w + kleft;
    00714         if(start < kright)
    00715         {
    00716             id += kright - start;
    00717             start = kright;
    00718         }
    00719     }
    00720     else
    00721     {
    00722         id += kright;
    00723         start = kright;
    00724         stop = w + kleft;
    00725     }
    00726
    00727     typedef typename PromoteTraits<
    00728             typename SrcAccessor::value_type,
    00729             typename KernelAccessor::value_type>::Promote SumType;
    00730
    00731     is += start;
    00732
    00733     for(int x=start; x<stop; ++x, ++is, ++id)
    00734     {
    00735         KernelIterator ik = kernel + kright;
    00736         SumType sum = NumericTraits<SumType>::zero();
    00737
    00738         SrcIterator iss = is + (-kright);
    00739         SrcIterator isend = is + (1 - kleft);
    00740         for(; iss != isend ; --ik, ++iss)
    00741         {
    00742             sum += ka(ik) * sa(iss);
    00743         }
    00744
    00745         da.set(detail::RequiresExplicitCast<typename
    00746                       DestAccessor::value_type>::cast(sum), id);
    00747     }
    00748 }
    00749
    00750 /********************************************************/
    00751 /*                                                      */
    00752 /*         Separable convolution functions              */
    00753 /*                                                      */
    00754 /********************************************************/
    00755
    00756 /** \addtogroup SeparableConvolution One-dimensional and separable convolution functions
    00757
    00758     Perform 1D convolution and separable filtering in 2 dimensions.
    00759
    00760     These generic convolution functions implement
    00761     the standard convolution operation for a wide range of images and
    00762     signals that fit into the required interface. They need a suitable
    00763     kernel to operate.
    00764 */
    00765 //@{
    00766
    00767 /** \brief Performs a 1-dimensional convolution of the source signal using the given
    00768     kernel.
    00769
    00770     The KernelIterator must point to the center iterator, and
    00771     the kernel's size is given by its left (kleft <= 0) and right
    00772     (kright >= 0) borders. The signal must always be larger than the kernel.
    00773     At those positions where the kernel does not completely fit
    00774     into the signal's range, the specified \ref BorderTreatmentMode is
    00775     applied.
    00776
    00777     The signal's value_type (SrcAccessor::value_type) must be a
    00778     linear space over the kernel's value_type (KernelAccessor::value_type),
    00779     i.e. addition of source values, multiplication with kernel values,
    00780     and NumericTraits must be defined.
    00781     The kernel's value_type must be an algebraic field,
    00782     i.e. the arithmetic operations (+, -, *, /) and NumericTraits must
    00783     be defined.
    00784
    00785     If <tt>start</tt> and <tt>stop</tt> are non-zero, the relation
    00786     <tt>0 <= start < stop <= width</tt> must hold (where <tt>width</tt>
    00787     is the length of the input array). The convolution is then restricted to that
    00788     subrange, and it is assumed that the output array only refers to that
    00789     subrange (i.e. <tt>id</tt> points to the element corresponding to
    00790     <tt>start</tt>). If <tt>start</tt> and <tt>stop</tt> are both zero
    00791     (the default), the entire array is convolved.
    00792
    00793     <b> Declarations:</b>
    00794
    00795     pass arguments explicitly:
    00796     \code
    00797     namespace vigra {
    00798         template <class SrcIterator, class SrcAccessor,
    00799                   class DestIterator, class DestAccessor,
    00800                   class KernelIterator, class KernelAccessor>
    00801         void convolveLine(SrcIterator is, SrcIterator isend, SrcAccessor sa,
    00802                           DestIterator id, DestAccessor da,
    00803                           KernelIterator ik, KernelAccessor ka,
    00804                           int kleft, int kright, BorderTreatmentMode border,
    00805                           int start = 0, int stop = 0 )
    00806     }
    00807     \endcode
    00808
    00809
    00810     use argument objects in conjunction with \ref ArgumentObjectFactories :
    00811     \code
    00812     namespace vigra {
    00813         template <class SrcIterator, class SrcAccessor,
    00814                   class DestIterator, class DestAccessor,
    00815                   class KernelIterator, class KernelAccessor>
    00816         void convolveLine(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    00817                           pair<DestIterator, DestAccessor> dest,
    00818                           tuple5<KernelIterator, KernelAccessor,
    00819                                  int, int, BorderTreatmentMode> kernel,
    00820                            int start = 0, int stop = 0)
    00821     }
    00822     \endcode
    00823
    00824     <b> Usage:</b>
    00825
    00826     <b>\#include</b> <vigra/separableconvolution.hxx>
    00827
    00828
    00829     \code
    00830     std::vector<float> src, dest;
    00831     ...
    00832
    00833     // define binomial filter of size 5
    00834     static float kernel[] =
    00835            { 1.0/16.0, 4.0/16.0, 6.0/16.0, 4.0/16.0, 1.0/16.0};
    00836
    00837     typedef vigra::StandardAccessor<float> FAccessor;
    00838     typedef vigra::StandardAccessor<float> KernelAccessor;
    00839
    00840
    00841     vigra::convolveLine(src.begin(), src.end(), FAccessor(), dest.begin(), FAccessor(),
    00842              kernel+2, KernelAccessor(), -2, 2, BORDER_TREATMENT_REFLECT);
    00843     //       ^^^^^^^^  this is the center of the kernel
    00844
    00845     \endcode
    00846
    00847     <b> Required Interface:</b>
    00848
    00849     \code
    00850     RandomAccessIterator is, isend;
    00851     RandomAccessIterator id;
    00852     RandomAccessIterator ik;
    00853
    00854     SrcAccessor src_accessor;
    00855     DestAccessor dest_accessor;
    00856     KernelAccessor kernel_accessor;
    00857
    00858     NumericTraits<SrcAccessor::value_type>::RealPromote s = src_accessor(is);
    00859
    00860     s = s + s;
    00861     s = kernel_accessor(ik) * s;
    00862
    00863     dest_accessor.set(
    00864         NumericTraits<DestAccessor::value_type>::fromRealPromote(s), id);
    00865
    00866     \endcode
    00867
    00868     If border == BORDER_TREATMENT_CLIP:
    00869
    00870     \code
    00871     NumericTraits<KernelAccessor::value_type>::RealPromote k = kernel_accessor(ik);
    00872
    00873     k = k + k;
    00874     k = k - k;
    00875     k = k * k;
    00876     k = k / k;
    00877
    00878     \endcode
    00879
    00880     <b> Preconditions:</b>
    00881
    00882     \code
    00883     kleft <= 0
    00884     kright >= 0
    00885     iend - is >= kright + kleft + 1
    00886     \endcode
    00887
    00888     If border == BORDER_TREATMENT_CLIP: Sum of kernel elements must be
    00889     != 0.
    00890
    00891 */
    00892 doxygen_overloaded_function(template <...> void convolveLine)
    00893
    00894 template <class SrcIterator, class SrcAccessor,
    00895           class DestIterator, class DestAccessor,
    00896           class KernelIterator, class KernelAccessor>
    00897 void convolveLine(SrcIterator is, SrcIterator iend, SrcAccessor sa,
    00898                   DestIterator id, DestAccessor da,
    00899                   KernelIterator ik, KernelAccessor ka,
    00900                   int kleft, int kright, BorderTreatmentMode border,
    00901                   int start = 0, int stop = 0)
    00902 {
    00903     typedef typename KernelAccessor::value_type KernelValue;
    00904
    00905     vigra_precondition(kleft <= 0,
    00906                  "convolveLine(): kleft must be <= 0.\n");
    00907     vigra_precondition(kright >= 0,
    00908                  "convolveLine(): kright must be >= 0.\n");
    00909
    00910     //    int w = iend - is;
    00911     int w = std::distance( is, iend );
    00912
    00913     vigra_precondition(w >= std::max(kright, -kleft) + 1,
    00914                  "convolveLine(): kernel longer than line.\n");
    00915
    00916     if(stop != 0)
    00917         vigra_precondition(0 <= start && start < stop && stop <= w,
    00918                         "convolveLine(): invalid subrange (start, stop).\n");
    00919
    00920     typedef typename PromoteTraits<
    00921             typename SrcAccessor::value_type,
    00922             typename KernelAccessor::value_type>::Promote SumType;
    00923     ArrayVector<SumType> a(iend - is);
    00924     switch(border)
    00925     {
    00926       case BORDER_TREATMENT_WRAP:
    00927       {
    00928         internalConvolveLineWrap(is, iend, sa, id, da, ik, ka, kleft, kright, start, stop);
    00929         break;
    00930       }
    00931       case BORDER_TREATMENT_AVOID:
    00932       {
    00933         internalConvolveLineAvoid(is, iend, sa, id, da, ik, ka, kleft, kright, start, stop);
    00934         break;
    00935       }
    00936       case BORDER_TREATMENT_REFLECT:
    00937       {
    00938         internalConvolveLineReflect(is, iend, sa, id, da, ik, ka, kleft, kright, start, stop);
    00939         break;
    00940       }
    00941       case BORDER_TREATMENT_REPEAT:
    00942       {
    00943         internalConvolveLineRepeat(is, iend, sa, id, da, ik, ka, kleft, kright, start, stop);
    00944         break;
    00945       }
    00946       case BORDER_TREATMENT_CLIP:
    00947       {
    00948         // find norm of kernel
    00949         typedef typename KernelAccessor::value_type KT;
    00950         KT norm = NumericTraits<KT>::zero();
    00951         KernelIterator iik = ik + kleft;
    00952         for(int i=kleft; i<=kright; ++i, ++iik)
    00953             norm += ka(iik);
    00954
    00955         vigra_precondition(norm != NumericTraits<KT>::zero(),
    00956                      "convolveLine(): Norm of kernel must be != 0"
    00957                      " in mode BORDER_TREATMENT_CLIP.\n");
    00958
    00959         internalConvolveLineClip(is, iend, sa, id, da, ik, ka, kleft, kright, norm, start, stop);
    00960         break;
    00961       }
    00962       case BORDER_TREATMENT_ZEROPAD:
    00963       {
    00964         internalConvolveLineZeropad(is, iend, sa, id, da, ik, ka, kleft, kright, start, stop);
    00965         break;
    00966       }
    00967       default:
    00968       {
    00969         vigra_precondition(0,
    00970                      "convolveLine(): Unknown border treatment mode.\n");
    00971       }
    00972     }
    00973 }
    00974
    00975 template <class SrcIterator, class SrcAccessor,
    00976           class DestIterator, class DestAccessor,
    00977           class KernelIterator, class KernelAccessor>
    00978 inline
    00979 void convolveLine(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    00980                   pair<DestIterator, DestAccessor> dest,
    00981                   tuple5<KernelIterator, KernelAccessor,
    00982                          int, int, BorderTreatmentMode> kernel,
    00983                   int start = 0, int stop = 0)
    00984 {
    00985     convolveLine(src.first, src.second, src.third,
    00986                  dest.first, dest.second,
    00987                  kernel.first, kernel.second,
    00988                  kernel.third, kernel.fourth, kernel.fifth, start, stop);
    00989 }
    00990
    00991 /********************************************************/
    00992 /*                                                      */
    00993 /*                      separableConvolveX              */
    00994 /*                                                      */
    00995 /********************************************************/
    00996
    00997 /** \brief Performs a 1 dimensional convolution in x direction.
    00998
    00999     It calls \ref convolveLine() for every row of the image. See \ref convolveLine()
    01000     for more information about required interfaces and vigra_preconditions.
    01001
    01002     <b> Declarations:</b>
    01003
    01004     pass arguments explicitly:
    01005     \code
    01006     namespace vigra {
    01007         template <class SrcImageIterator, class SrcAccessor,
    01008                   class DestImageIterator, class DestAccessor,
    01009                   class KernelIterator, class KernelAccessor>
    01010         void separableConvolveX(SrcImageIterator supperleft,
    01011                                 SrcImageIterator slowerright, SrcAccessor sa,
    01012                                 DestImageIterator dupperleft, DestAccessor da,
    01013                                 KernelIterator ik, KernelAccessor ka,
    01014                                 int kleft, int kright, BorderTreatmentMode border)
    01015     }
    01016     \endcode
    01017
    01018
    01019     use argument objects in conjunction with \ref ArgumentObjectFactories :
    01020     \code
    01021     namespace vigra {
    01022         template <class SrcImageIterator, class SrcAccessor,
    01023                   class DestImageIterator, class DestAccessor,
    01024                   class KernelIterator, class KernelAccessor>
    01025         void separableConvolveX(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
    01026                                 pair<DestImageIterator, DestAccessor> dest,
    01027                                 tuple5<KernelIterator, KernelAccessor,
    01028                                              int, int, BorderTreatmentMode> kernel)
    01029     }
    01030     \endcode
    01031
    01032     <b> Usage:</b>
    01033
    01034     <b>\#include</b> <vigra/separableconvolution.hxx>
    01035
    01036
    01037     \code
    01038     vigra::FImage src(w,h), dest(w,h);
    01039     ...
    01040
    01041     // define Gaussian kernel with std. deviation 3.0
    01042     vigra::Kernel1D<double> kernel;
    01043     kernel.initGaussian(3.0);
    01044
    01045     vigra::separableConvolveX(srcImageRange(src), destImage(dest), kernel1d(kernel));
    01046
    01047     \endcode
    01048
    01049 */
    01050 doxygen_overloaded_function(template <...> void separableConvolveX)
    01051
    01052 template <class SrcIterator, class SrcAccessor,
    01053           class DestIterator, class DestAccessor,
    01054           class KernelIterator, class KernelAccessor>
    01055 void separableConvolveX(SrcIterator supperleft,
    01056                         SrcIterator slowerright, SrcAccessor sa,
    01057                         DestIterator dupperleft, DestAccessor da,
    01058                         KernelIterator ik, KernelAccessor ka,
    01059                         int kleft, int kright, BorderTreatmentMode border)
    01060 {
    01061     typedef typename KernelAccessor::value_type KernelValue;
    01062
    01063     vigra_precondition(kleft <= 0,
    01064                  "separableConvolveX(): kleft must be <= 0.\n");
    01065     vigra_precondition(kright >= 0,
    01066                  "separableConvolveX(): kright must be >= 0.\n");
    01067
    01068     int w = slowerright.x - supperleft.x;
    01069     int h = slowerright.y - supperleft.y;
    01070
    01071     vigra_precondition(w >= std::max(kright, -kleft) + 1,
    01072                  "separableConvolveX(): kernel longer than line\n");
    01073
    01074     int y;
    01075
    01076     for(y=0; y<h; ++y, ++supperleft.y, ++dupperleft.y)
    01077     {
    01078         typename SrcIterator::row_iterator rs = supperleft.rowIterator();
    01079         typename DestIterator::row_iterator rd = dupperleft.rowIterator();
    01080
    01081         convolveLine(rs, rs+w, sa, rd, da,
    01082                      ik, ka, kleft, kright, border);
    01083     }
    01084 }
    01085
    01086 template <class SrcIterator, class SrcAccessor,
    01087           class DestIterator, class DestAccessor,
    01088           class KernelIterator, class KernelAccessor>
    01089 inline void
    01090 separableConvolveX(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    01091                   pair<DestIterator, DestAccessor> dest,
    01092                   tuple5<KernelIterator, KernelAccessor,
    01093                          int, int, BorderTreatmentMode> kernel)
    01094 {
    01095     separableConvolveX(src.first, src.second, src.third,
    01096                  dest.first, dest.second,
    01097                  kernel.first, kernel.second,
    01098                  kernel.third, kernel.fourth, kernel.fifth);
    01099 }
    01100
    01101
    01102
    01103 /********************************************************/
    01104 /*                                                      */
    01105 /*                      separableConvolveY              */
    01106 /*                                                      */
    01107 /********************************************************/
    01108
    01109 /** \brief Performs a 1 dimensional convolution in y direction.
    01110
    01111     It calls \ref convolveLine() for every column of the image. See \ref convolveLine()
    01112     for more information about required interfaces and vigra_preconditions.
    01113
    01114     <b> Declarations:</b>
    01115
    01116     pass arguments explicitly:
    01117     \code
    01118     namespace vigra {
    01119         template <class SrcImageIterator, class SrcAccessor,
    01120                   class DestImageIterator, class DestAccessor,
    01121                   class KernelIterator, class KernelAccessor>
    01122         void separableConvolveY(SrcImageIterator supperleft,
    01123                                 SrcImageIterator slowerright, SrcAccessor sa,
    01124                                 DestImageIterator dupperleft, DestAccessor da,
    01125                                 KernelIterator ik, KernelAccessor ka,
    01126                                 int kleft, int kright, BorderTreatmentMode border)
    01127     }
    01128     \endcode
    01129
    01130
    01131     use argument objects in conjunction with \ref ArgumentObjectFactories :
    01132     \code
    01133     namespace vigra {
    01134         template <class SrcImageIterator, class SrcAccessor,
    01135                   class DestImageIterator, class DestAccessor,
    01136                   class KernelIterator, class KernelAccessor>
    01137         void separableConvolveY(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
    01138                                 pair<DestImageIterator, DestAccessor> dest,
    01139                                 tuple5<KernelIterator, KernelAccessor,
    01140                                              int, int, BorderTreatmentMode> kernel)
    01141     }
    01142     \endcode
    01143
    01144     <b> Usage:</b>
    01145
    01146     <b>\#include</b> <vigra/separableconvolution.hxx>
    01147
    01148
    01149     \code
    01150     vigra::FImage src(w,h), dest(w,h);
    01151     ...
    01152
    01153     // define Gaussian kernel with std. deviation 3.0
    01154     vigra::Kernel1D kernel;
    01155     kernel.initGaussian(3.0);
    01156
    01157     vigra::separableConvolveY(srcImageRange(src), destImage(dest), kernel1d(kernel));
    01158
    01159     \endcode
    01160
    01161 */
    01162 doxygen_overloaded_function(template <...> void separableConvolveY)
    01163
    01164 template <class SrcIterator, class SrcAccessor,
    01165           class DestIterator, class DestAccessor,
    01166           class KernelIterator, class KernelAccessor>
    01167 void separableConvolveY(SrcIterator supperleft,
    01168                         SrcIterator slowerright, SrcAccessor sa,
    01169                         DestIterator dupperleft, DestAccessor da,
    01170                         KernelIterator ik, KernelAccessor ka,
    01171                         int kleft, int kright, BorderTreatmentMode border)
    01172 {
    01173     typedef typename KernelAccessor::value_type KernelValue;
    01174
    01175     vigra_precondition(kleft <= 0,
    01176                  "separableConvolveY(): kleft must be <= 0.\n");
    01177     vigra_precondition(kright >= 0,
    01178                  "separableConvolveY(): kright must be >= 0.\n");
    01179
    01180     int w = slowerright.x - supperleft.x;
    01181     int h = slowerright.y - supperleft.y;
    01182
    01183     vigra_precondition(h >= std::max(kright, -kleft) + 1,
    01184                  "separableConvolveY(): kernel longer than line\n");
    01185
    01186     int x;
    01187
    01188     for(x=0; x<w; ++x, ++supperleft.x, ++dupperleft.x)
    01189     {
    01190         typename SrcIterator::column_iterator cs = supperleft.columnIterator();
    01191         typename DestIterator::column_iterator cd = dupperleft.columnIterator();
    01192
    01193         convolveLine(cs, cs+h, sa, cd, da,
    01194                      ik, ka, kleft, kright, border);
    01195     }
    01196 }
    01197
    01198 template <class SrcIterator, class SrcAccessor,
    01199           class DestIterator, class DestAccessor,
    01200           class KernelIterator, class KernelAccessor>
    01201 inline void
    01202 separableConvolveY(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    01203                   pair<DestIterator, DestAccessor> dest,
    01204                   tuple5<KernelIterator, KernelAccessor,
    01205                          int, int, BorderTreatmentMode> kernel)
    01206 {
    01207     separableConvolveY(src.first, src.second, src.third,
    01208                  dest.first, dest.second,
    01209                  kernel.first, kernel.second,
    01210                  kernel.third, kernel.fourth, kernel.fifth);
    01211 }
    01212
    01213 //@}
    01214
    01215 /********************************************************/
    01216 /*                                                      */
    01217 /*                      Kernel1D                        */
    01218 /*                                                      */
    01219 /********************************************************/
    01220
    01221 /** \brief Generic 1 dimensional convolution kernel.
    01222
    01223     This kernel may be used for convolution of 1 dimensional signals or for
    01224     separable convolution of multidimensional signals.
    01225
    01226     Convolution functions access the kernel via a 1 dimensional random access
    01227     iterator which they get by calling \ref center(). This iterator
    01228     points to the center of the kernel. The kernel's size is given by its left() (<=0)
    01229     and right() (>= 0) methods. The desired border treatment mode is
    01230     returned by borderTreatment().
    01231
    01232     The different init functions create a kernel with the specified
    01233     properties. The kernel's value_type must be a linear space, i.e. it
    01234     must define multiplication with doubles and NumericTraits.
    01235
    01236
    01237     The kernel defines a factory function kernel1d() to create an argument object
    01238     (see \ref KernelArgumentObjectFactories).
    01239
    01240     <b> Usage:</b>
    01241
    01242     <b>\#include</b> <vigra/stdconvolution.hxx>
    01243
    01244     \code
    01245     vigra::FImage src(w,h), dest(w,h);
    01246     ...
    01247
    01248     // define Gaussian kernel with std. deviation 3.0
    01249     vigra::Kernel1D kernel;
    01250     kernel.initGaussian(3.0);
    01251
    01252     vigra::separableConvolveX(srcImageRange(src), destImage(dest), kernel1d(kernel));
    01253     \endcode
    01254
    01255     <b> Required Interface:</b>
    01256
    01257     \code
    01258     value_type v = vigra::NumericTraits<value_type>::one(); // if norm is not
    01259                                                             // given explicitly
    01260     double d;
    01261
    01262     v = d * v;
    01263     \endcode
    01264 */
    01265
    01266 template <class ARITHTYPE>
    01267 class Kernel1D
    01268 {
    01269   public:
    01270     typedef ArrayVector<ARITHTYPE> InternalVector;
    01271
    01272         /** the kernel's value type
    01273         */
    01274     typedef typename InternalVector::value_type value_type;
    01275
    01276         /** the kernel's reference type
    01277         */
    01278     typedef typename InternalVector::reference reference;
    01279
    01280         /** the kernel's const reference type
    01281         */
    01282     typedef typename InternalVector::const_reference const_reference;
    01283
    01284         /** deprecated -- use Kernel1D::iterator
    01285         */
    01286     typedef typename InternalVector::iterator Iterator;
    01287
    01288         /** 1D random access iterator over the kernel's values
    01289         */
    01290     typedef typename InternalVector::iterator iterator;
    01291
    01292         /** const 1D random access iterator over the kernel's values
    01293         */
    01294     typedef typename InternalVector::const_iterator const_iterator;
    01295
    01296         /** the kernel's accessor
    01297         */
    01298     typedef StandardAccessor<ARITHTYPE> Accessor;
    01299
    01300         /** the kernel's const accessor
    01301         */
    01302     typedef StandardConstAccessor<ARITHTYPE> ConstAccessor;
    01303
    01304     struct InitProxy
    01305     {
    01306         InitProxy(Iterator i, int count, value_type & norm)
    01307         : iter_(i), base_(i),
    01308           count_(count), sum_(count),
    01309           norm_(norm)
    01310         {}
    01311
    01312         ~InitProxy()
    01313 #ifndef _MSC_VER
    01314              throw(PreconditionViolation)
    01315 #endif
    01316         {
    01317             vigra_precondition(count_ == 1 || count_ == sum_,
    01318                   "Kernel1D::initExplicitly(): "
    01319                   "Wrong number of init values.");
    01320         }
    01321
    01322         InitProxy & operator,(value_type const & v)
    01323         {
    01324             if(sum_ == count_)
    01325                 norm_ = *iter_;
    01326
    01327             norm_ += v;
    01328
    01329             --count_;
    01330
    01331             if(count_ > 0)
    01332             {
    01333                 ++iter_;
    01334                 *iter_ = v;
    01335             }
    01336             return *this;
    01337         }
    01338
    01339         Iterator iter_, base_;
    01340         int count_, sum_;
    01341         value_type & norm_;
    01342     };
    01343
    01344     static value_type one() { return NumericTraits<value_type>::one(); }
    01345
    01346         /** Default constructor.
    01347             Creates a kernel of size 1 which would copy the signal
    01348             unchanged.
    01349         */
    01350     Kernel1D()
    01351     : kernel_(),
    01352       left_(0),
    01353       right_(0),
    01354       border_treatment_(BORDER_TREATMENT_REFLECT),
    01355       norm_(one())
    01356     {
    01357         kernel_.push_back(norm_);
    01358     }
    01359
    01360         /** Copy constructor.
    01361         */
    01362     Kernel1D(Kernel1D const & k)
    01363     : kernel_(k.kernel_),
    01364       left_(k.left_),
    01365       right_(k.right_),
    01366       border_treatment_(k.border_treatment_),
    01367       norm_(k.norm_)
    01368     {}
    01369
    01370         /** Construct from kernel with different element type, e.g. double => FixedPoint16.
    01371         */
    01372     template <class U>
    01373     Kernel1D(Kernel1D<U> const & k)
    01374     : kernel_(k.center()+k.left(), k.center()+k.right()+1),
    01375       left_(k.left()),
    01376       right_(k.right()),
    01377       border_treatment_(k.borderTreatment()),
    01378       norm_(k.norm())
    01379     {}
    01380
    01381         /** Copy assignment.
    01382         */
    01383     Kernel1D & operator=(Kernel1D const & k)
    01384     {
    01385         if(this != &k)
    01386         {
    01387             left_ = k.left_;
    01388             right_ = k.right_;
    01389             border_treatment_ = k.border_treatment_;
    01390             norm_ = k.norm_;
    01391             kernel_ = k.kernel_;
    01392         }
    01393         return *this;
    01394     }
    01395
    01396         /** Initialization.
    01397             This initializes the kernel with the given constant. The norm becomes
    01398             v*size().
    01399
    01400             Instead of a single value an initializer list of length size()
    01401             can be used like this:
    01402
    01403             \code
    01404             vigra::Kernel1D<float> roberts_gradient_x;
    01405
    01406             roberts_gradient_x.initExplicitly(0, 1) = 1.0, -1.0;
    01407             \endcode
    01408
    01409             In this case, the norm will be set to the sum of the init values.
    01410             An initializer list of wrong length will result in a run-time error.
    01411         */
    01412     InitProxy operator=(value_type const & v)
    01413     {
    01414         int size = right_ - left_ + 1;
    01415         for(unsigned int i=0; i<kernel_.size(); ++i) kernel_[i] = v;
    01416         norm_ = (double)size*v;
    01417
    01418         return InitProxy(kernel_.begin(), size, norm_);
    01419     }
    01420
    01421         /** Destructor.
    01422         */
    01423     ~Kernel1D()
    01424     {}
    01425
    01426         /**
    01427             Init as a sampled Gaussian function. The radius of the kernel is
    01428             always 3*std_dev. '<tt>norm</tt>' denotes the sum of all bins of the kernel
    01429             (i.e. the kernel is corrected for the normalization error introduced
    01430              by windowing the Gaussian to a finite interval). However,
    01431             if <tt>norm</tt> is 0.0, the kernel is normalized to 1 by the analytic
    01432             expression for the Gaussian, and <b>no</b> correction for the windowing
    01433             error is performed. If <tt>windowRatio = 0.0</tt>, the radius of the filter
    01434             window is <tt>radius = round(3.0 * std_dev)</tt>, otherwise it is
    01435             <tt>radius = round(windowRatio * std_dev)</tt> (where <tt>windowRatio > 0.0</tt>
    01436             is required).
    01437
    01438             Precondition:
    01439             \code
    01440             std_dev >= 0.0
    01441             \endcode
    01442
    01443             Postconditions:
    01444             \code
    01445             1. left()  == -(int)(3.0*std_dev + 0.5)
    01446             2. right() ==  (int)(3.0*std_dev + 0.5)
    01447             3. borderTreatment() == BORDER_TREATMENT_REFLECT
    01448             4. norm() == norm
    01449             \endcode
    01450         */
    01451     void initGaussian(double std_dev, value_type norm, double windowRatio = 0.0);
    01452
    01453         /** Init as a Gaussian function with norm 1.
    01454          */
    01455     void initGaussian(double std_dev)
    01456     {
    01457         initGaussian(std_dev, one());
    01458     }
    01459
    01460
    01461         /**
    01462             Init as Lindeberg's discrete analog of the Gaussian function. The radius of the kernel is
    01463             always 3*std_dev. 'norm' denotes the sum of all bins of the kernel.
    01464
    01465             Precondition:
    01466             \code
    01467             std_dev >= 0.0
    01468             \endcode
    01469
    01470             Postconditions:
    01471             \code
    01472             1. left()  == -(int)(3.0*std_dev + 0.5)
    01473             2. right() ==  (int)(3.0*std_dev + 0.5)
    01474             3. borderTreatment() == BORDER_TREATMENT_REFLECT
    01475             4. norm() == norm
    01476             \endcode
    01477         */
    01478     void initDiscreteGaussian(double std_dev, value_type norm);
    01479
    01480         /** Init as a Lindeberg's discrete analog of the Gaussian function
    01481             with norm 1.
    01482          */
    01483     void initDiscreteGaussian(double std_dev)
    01484     {
    01485         initDiscreteGaussian(std_dev, one());
    01486     }
    01487
    01488         /**
    01489             Init as a Gaussian derivative of order '<tt>order</tt>'.
    01490             The radius of the kernel is always <tt>3*std_dev + 0.5*order</tt>.
    01491             '<tt>norm</tt>' denotes the norm of the kernel so that the
    01492             following condition is fulfilled:
    01493
    01494             \f[ \sum_{i=left()}^{right()}
    01495                          \frac{(-i)^{order}kernel[i]}{order!} = norm
    01496             \f]
    01497
    01498             Thus, the kernel will be corrected for the error introduced
    01499             by windowing the Gaussian to a finite interval. However,
    01500             if <tt>norm</tt> is 0.0, the kernel is normalized to 1 by the analytic
    01501             expression for the Gaussian derivative, and <b>no</b> correction for the
    01502             windowing error is performed. If <tt>windowRatio = 0.0</tt>, the radius
    01503             of the filter window is <tt>radius = round(3.0 * std_dev + 0.5 * order)</tt>,
    01504             otherwise it is <tt>radius = round(windowRatio * std_dev)</tt> (where
    01505             <tt>windowRatio > 0.0</tt> is required).
    01506
    01507             Preconditions:
    01508             \code
    01509             1. std_dev >= 0.0
    01510             2. order   >= 1
    01511             \endcode
    01512
    01513             Postconditions:
    01514             \code
    01515             1. left()  == -(int)(3.0*std_dev + 0.5*order + 0.5)
    01516             2. right() ==  (int)(3.0*std_dev + 0.5*order + 0.5)
    01517             3. borderTreatment() == BORDER_TREATMENT_REFLECT
    01518             4. norm() == norm
    01519             \endcode
    01520         */
    01521     void initGaussianDerivative(double std_dev, int order, value_type norm, double windowRatio = 0.0);
    01522
    01523         /** Init as a Gaussian derivative with norm 1.
    01524          */
    01525     void initGaussianDerivative(double std_dev, int order)
    01526     {
    01527         initGaussianDerivative(std_dev, order, one());
    01528     }
    01529
    01530         /**
    01531             Init an optimal 3-tap smoothing filter.
    01532             The filter values are
    01533
    01534             \code
    01535             [0.216, 0.568, 0.216]
    01536             \endcode
    01537
    01538             These values are optimal in the sense that the 3x3 filter obtained by separable application
    01539             of this filter is the best possible 3x3 approximation to a Gaussian filter.
    01540             The equivalent Gaussian has sigma = 0.680.
    01541
    01542             Postconditions:
    01543             \code
    01544             1. left()  == -1
    01545             2. right() ==  1
    01546             3. borderTreatment() == BORDER_TREATMENT_REFLECT
    01547             4. norm() == 1.0
    01548             \endcode
    01549         */
    01550     void initOptimalSmoothing3()
    01551     {
    01552         this->initExplicitly(-1, 1) = 0.216, 0.568, 0.216;
    01553         this->setBorderTreatment(BORDER_TREATMENT_REFLECT);
    01554     }
    01555
    01556         /**
    01557             Init an optimal 3-tap smoothing filter to be used in the context of first derivative computation.
    01558             This filter must be used in conjunction with the symmetric difference filter (see initSymmetricDifference()),
    01559             such that the difference filter is applied along one dimension, and the smoothing filter along the other.
    01560             The filter values are
    01561
    01562             \code
    01563             [0.224365, 0.55127, 0.224365]
    01564             \endcode
    01565
    01566             These values are optimal in the sense that the 3x3 filter obtained by combining
    01567             this filter with the symmetric difference is the best possible 3x3 approximation to a
    01568             Gaussian first derivative filter. The equivalent Gaussian has sigma = 0.675.
    01569
    01570             Postconditions:
    01571             \code
    01572             1. left()  == -1
    01573             2. right() ==  1
    01574             3. borderTreatment() == BORDER_TREATMENT_REFLECT
    01575             4. norm() == 1.0
    01576             \endcode
    01577         */
    01578     void initOptimalFirstDerivativeSmoothing3()
    01579     {
    01580         this->initExplicitly(-1, 1) = 0.224365, 0.55127, 0.224365;
    01581         this->setBorderTreatment(BORDER_TREATMENT_REFLECT);
    01582     }
    01583
    01584         /**
    01585             Init an optimal 3-tap smoothing filter to be used in the context of second derivative computation.
    01586             This filter must be used in conjunction with the 3-tap second difference filter (see initSecondDifference3()),
    01587             such that the difference filter is applied along one dimension, and the smoothing filter along the other.
    01588             The filter values are
    01589
    01590             \code
    01591             [0.13, 0.74, 0.13]
    01592             \endcode
    01593
    01594             These values are optimal in the sense that the 3x3 filter obtained by combining
    01595             this filter with the 3-tap second difference is the best possible 3x3 approximation to a
    01596             Gaussian second derivative filter. The equivalent Gaussian has sigma = 0.433.
    01597
    01598             Postconditions:
    01599             \code
    01600             1. left()  == -1
    01601             2. right() ==  1
    01602             3. borderTreatment() == BORDER_TREATMENT_REFLECT
    01603             4. norm() == 1.0
    01604             \endcode
    01605         */
    01606     void initOptimalSecondDerivativeSmoothing3()
    01607     {
    01608         this->initExplicitly(-1, 1) = 0.13, 0.74, 0.13;
    01609         this->setBorderTreatment(BORDER_TREATMENT_REFLECT);
    01610     }
    01611
    01612         /**
    01613             Init an optimal 5-tap smoothing filter.
    01614             The filter values are
    01615
    01616             \code
    01617             [0.03134, 0.24, 0.45732, 0.24, 0.03134]
    01618             \endcode
    01619
    01620             These values are optimal in the sense that the 5x5 filter obtained by separable application
    01621             of this filter is the best possible 5x5 approximation to a Gaussian filter.
    01622             The equivalent Gaussian has sigma = 0.867.
    01623
    01624             Postconditions:
    01625             \code
    01626             1. left()  == -2
    01627             2. right() ==  2
    01628             3. borderTreatment() == BORDER_TREATMENT_REFLECT
    01629             4. norm() == 1.0
    01630             \endcode
    01631         */
    01632     void initOptimalSmoothing5()
    01633     {
    01634         this->initExplicitly(-2, 2) = 0.03134, 0.24, 0.45732, 0.24, 0.03134;
    01635         this->setBorderTreatment(BORDER_TREATMENT_REFLECT);
    01636     }
    01637
    01638         /**
    01639             Init an optimal 5-tap smoothing filter to be used in the context of first derivative computation.
    01640            This filter must be used in conjunction with the optimal 5-tap first derivative filter
    01641            (see initOptimalFirstDerivative5()),  such that the derivative filter is applied along one dimension,
    01642            and the smoothing filter along the other. The filter values are
    01643
    01644             \code
    01645             [0.04255, 0.241, 0.4329, 0.241, 0.04255]
    01646             \endcode
    01647
    01648             These values are optimal in the sense that the 5x5 filter obtained by combining
    01649             this filter with the optimal 5-tap first derivative is the best possible 5x5 approximation to a
    01650             Gaussian first derivative filter. The equivalent Gaussian has sigma = 0.906.
    01651
    01652             Postconditions:
    01653             \code
    01654             1. left()  == -2
    01655             2. right() ==  2
    01656             3. borderTreatment() == BORDER_TREATMENT_REFLECT
    01657             4. norm() == 1.0
    01658             \endcode
    01659         */
    01660     void initOptimalFirstDerivativeSmoothing5()
    01661     {
    01662         this->initExplicitly(-2, 2) = 0.04255, 0.241, 0.4329, 0.241, 0.04255;
    01663         this->setBorderTreatment(BORDER_TREATMENT_REFLECT);
    01664     }
    01665
    01666         /**
    01667             Init an optimal 5-tap smoothing filter to be used in the context of second derivative computation.
    01668            This filter must be used in conjunction with the optimal 5-tap second derivative filter
    01669            (see initOptimalSecondDerivative5()), such that the derivative filter is applied along one dimension,
    01670            and the smoothing filter along the other. The filter values are
    01671
    01672             \code
    01673             [0.0243, 0.23556, 0.48028, 0.23556, 0.0243]
    01674             \endcode
    01675
    01676             These values are optimal in the sense that the 5x5 filter obtained by combining
    01677             this filter with the optimal 5-tap second derivative is the best possible 5x5 approximation to a
    01678             Gaussian second derivative filter. The equivalent Gaussian has sigma = 0.817.
    01679
    01680             Postconditions:
    01681             \code
    01682             1. left()  == -2
    01683             2. right() ==  2
    01684             3. borderTreatment() == BORDER_TREATMENT_REFLECT
    01685             4. norm() == 1.0
    01686             \endcode
    01687         */
    01688     void initOptimalSecondDerivativeSmoothing5()
    01689     {
    01690         this->initExplicitly(-2, 2) = 0.0243, 0.23556, 0.48028, 0.23556, 0.0243;
    01691         this->setBorderTreatment(BORDER_TREATMENT_REFLECT);
    01692     }
    01693
    01694         /**
    01695             Init a 5-tap filter as defined by Peter Burt in the context of pyramid creation.
    01696             The filter values are
    01697
    01698             \code
    01699             [a, 0.25, 0.5-2*a, 0.25, a]
    01700             \endcode
    01701
    01702             The default <tt>a = 0.04785</tt> is optimal in the sense that it minimizes the difference
    01703             to a true Gaussian filter (which would have sigma = 0.975). For other values of <tt>a</tt>, the scale
    01704             of the most similar Gaussian can be approximated by
    01705
    01706             \code
    01707             sigma = 5.1 * a + 0.731
    01708             \endcode
    01709
    01710             Preconditions:
    01711             \code
    01712             0 <= a <= 0.125
    01713             \endcode
    01714
    01715             Postconditions:
    01716             \code
    01717             1. left()  == -2
    01718             2. right() ==  2
    01719             3. borderTreatment() == BORDER_TREATMENT_REFLECT
    01720             4. norm() == 1.0
    01721             \endcode
    01722         */
    01723     void initBurtFilter(double a = 0.04785)
    01724     {
    01725         vigra_precondition(a >= 0.0 && a <= 0.125,
    01726             "Kernel1D::initBurtFilter(): 0 <= a <= 0.125 required.");
    01727         this->initExplicitly(-2, 2) = a, 0.25, 0.5 - 2.0*a, 0.25, a;
    01728         this->setBorderTreatment(BORDER_TREATMENT_REFLECT);
    01729     }
    01730
    01731         /**
    01732             Init as a Binomial filter. 'norm' denotes the sum of all bins
    01733             of the kernel.
    01734
    01735             Precondition:
    01736             \code
    01737             radius   >= 0
    01738             \endcode
    01739
    01740             Postconditions:
    01741             \code
    01742             1. left()  == -radius
    01743             2. right() ==  radius
    01744             3. borderTreatment() == BORDER_TREATMENT_REFLECT
    01745             4. norm() == norm
    01746             \endcode
    01747         */
    01748     void initBinomial(int radius, value_type norm);
    01749
    01750         /** Init as a Binomial filter with norm 1.
    01751          */
    01752     void initBinomial(int radius)
    01753     {
    01754         initBinomial(radius, one());
    01755     }
    01756
    01757         /**
    01758             Init as an Averaging filter. 'norm' denotes the sum of all bins
    01759             of the kernel. The window size is (2*radius+1) * (2*radius+1)
    01760
    01761             Precondition:
    01762             \code
    01763             radius   >= 0
    01764             \endcode
    01765
    01766             Postconditions:
    01767             \code
    01768             1. left()  == -radius
    01769             2. right() ==  radius
    01770             3. borderTreatment() == BORDER_TREATMENT_CLIP
    01771             4. norm() == norm
    01772             \endcode
    01773         */
    01774     void initAveraging(int radius, value_type norm);
    01775
    01776         /** Init as an Averaging filter with norm 1.
    01777          */
    01778     void initAveraging(int radius)
    01779     {
    01780         initAveraging(radius, one());
    01781     }
    01782
    01783         /**
    01784            Init as a symmetric gradient filter of the form
    01785            <TT>[ 0.5 * norm, 0.0 * norm, -0.5 * norm]</TT>
    01786
    01787            <b>Deprecated</b>. Use initSymmetricDifference() instead.
    01788
    01789             Postconditions:
    01790             \code
    01791             1. left()  == -1
    01792             2. right() ==  1
    01793             3. borderTreatment() == BORDER_TREATMENT_REPEAT
    01794             4. norm() == norm
    01795             \endcode
    01796         */
    01797     void
    01798     initSymmetricGradient(value_type norm )
    01799     {
    01800         initSymmetricDifference(norm);
    01801         setBorderTreatment(BORDER_TREATMENT_REPEAT);
    01802     }
    01803
    01804         /** Init as a symmetric gradient filter with norm 1.
    01805
    01806            <b>Deprecated</b>. Use initSymmetricDifference() instead.
    01807          */
    01808     void initSymmetricGradient()
    01809     {
    01810         initSymmetricGradient(one());
    01811     }
    01812
    01813         /** Init as the 2-tap forward difference filter.
    01814              The filter values are
    01815
    01816             \code
    01817             [1.0, -1.0]
    01818             \endcode
    01819
    01820             (note that filters are reflected by the convolution algorithm,
    01821              and we get a forward difference after reflection).
    01822
    01823             Postconditions:
    01824             \code
    01825             1. left()  == -1
    01826             2. right() ==  0
    01827             3. borderTreatment() == BORDER_TREATMENT_REFLECT
    01828             4. norm() == 1.0
    01829             \endcode
    01830           */
    01831     void initForwardDifference()
    01832     {
    01833         this->initExplicitly(-1, 0) = 1.0, -1.0;
    01834         this->setBorderTreatment(BORDER_TREATMENT_REFLECT);
    01835     }
    01836
    01837         /** Init as the 2-tap backward difference filter.
    01838             The filter values are
    01839
    01840             \code
    01841             [1.0, -1.0]
    01842             \endcode
    01843
    01844             (note that filters are reflected by the convolution algorithm,
    01845              and we get a forward difference after reflection).
    01846
    01847             Postconditions:
    01848             \code
    01849             1. left()  == 0
    01850             2. right() ==  1
    01851             3. borderTreatment() == BORDER_TREATMENT_REFLECT
    01852             4. norm() == 1.0
    01853             \endcode
    01854           */
    01855     void initBackwardDifference()
    01856     {
    01857         this->initExplicitly(0, 1) = 1.0, -1.0;
    01858         this->setBorderTreatment(BORDER_TREATMENT_REFLECT);
    01859     }
    01860
    01861     void
    01862     initSymmetricDifference(value_type norm );
    01863
    01864         /** Init as the 3-tap symmetric difference filter
    01865             The filter values are
    01866
    01867             \code
    01868             [0.5, 0, -0.5]
    01869             \endcode
    01870
    01871             Postconditions:
    01872             \code
    01873             1. left()  == -1
    01874             2. right() ==  1
    01875             3. borderTreatment() == BORDER_TREATMENT_REFLECT
    01876             4. norm() == 1.0
    01877             \endcode
    01878           */
    01879     void initSymmetricDifference()
    01880     {
    01881         initSymmetricDifference(one());
    01882     }
    01883
    01884         /**
    01885             Init the 3-tap second difference filter.
    01886             The filter values are
    01887
    01888             \code
    01889             [1, -2, 1]
    01890             \endcode
    01891
    01892             Postconditions:
    01893             \code
    01894             1. left()  == -1
    01895             2. right() ==  1
    01896             3. borderTreatment() == BORDER_TREATMENT_REFLECT
    01897             4. norm() == 1
    01898             \endcode
    01899         */
    01900     void
    01901     initSecondDifference3()
    01902     {
    01903         this->initExplicitly(-1, 1) = 1.0, -2.0, 1.0;
    01904         this->setBorderTreatment(BORDER_TREATMENT_REFLECT);
    01905     }
    01906
    01907         /**
    01908             Init an optimal 5-tap first derivative filter.
    01909            This filter must be used in conjunction with the corresponding 5-tap smoothing filter
    01910            (see initOptimalFirstDerivativeSmoothing5()), such that the derivative filter is applied along one dimension,
    01911             and the smoothing filter along the other.
    01912             The filter values are
    01913
    01914             \code
    01915             [0.1, 0.3, 0.0, -0.3, -0.1]
    01916             \endcode
    01917
    01918             These values are optimal in the sense that the 5x5 filter obtained by combining
    01919             this filter with the corresponding 5-tap smoothing filter is the best possible 5x5 approximation to a
    01920             Gaussian first derivative filter. The equivalent Gaussian has sigma = 0.906.
    01921
    01922             If the filter is instead separably combined with itself, an almost optimal approximation of the
    01923             mixed second Gaussian derivative at scale sigma = 0.899 results.
    01924
    01925             Postconditions:
    01926             \code
    01927             1. left()  == -2
    01928             2. right() ==  2
    01929             3. borderTreatment() == BORDER_TREATMENT_REFLECT
    01930             4. norm() == 1.0
    01931             \endcode
    01932         */
    01933     void initOptimalFirstDerivative5()
    01934     {
    01935         this->initExplicitly(-2, 2) = 0.1, 0.3, 0.0, -0.3, -0.1;
    01936         this->setBorderTreatment(BORDER_TREATMENT_REFLECT);
    01937     }
    01938
    01939         /**
    01940             Init an optimal 5-tap second derivative filter.
    01941            This filter must be used in conjunction with the corresponding 5-tap smoothing filter
    01942            (see initOptimalSecondDerivativeSmoothing5()), such that the derivative filter is applied along one dimension,
    01943             and the smoothing filter along the other.
    01944             The filter values are
    01945
    01946             \code
    01947             [0.22075, 0.117, -0.6755, 0.117, 0.22075]
    01948             \endcode
    01949
    01950             These values are optimal in the sense that the 5x5 filter obtained by combining
    01951             this filter with the corresponding 5-tap smoothing filter is the best possible 5x5 approximation to a
    01952             Gaussian second derivative filter. The equivalent Gaussian has sigma = 0.817.
    01953
    01954             Postconditions:
    01955             \code
    01956             1. left()  == -2
    01957             2. right() ==  2
    01958             3. borderTreatment() == BORDER_TREATMENT_REFLECT
    01959             4. norm() == 1.0
    01960             \endcode
    01961         */
    01962     void initOptimalSecondDerivative5()
    01963     {
    01964         this->initExplicitly(-2, 2) = 0.22075, 0.117, -0.6755, 0.117, 0.22075;
    01965         this->setBorderTreatment(BORDER_TREATMENT_REFLECT);
    01966     }
    01967
    01968         /** Init the kernel by an explicit initializer list.
    01969             The left and right boundaries of the kernel must be passed.
    01970             A comma-separated initializer list is given after the assignment
    01971             operator. This function is used like this:
    01972
    01973             \code
    01974             // define horizontal Roberts filter
    01975             vigra::Kernel1D<float> roberts_gradient_x;
    01976
    01977             roberts_gradient_x.initExplicitly(0, 1) = 1.0, -1.0;
    01978             \endcode
    01979
    01980             The norm is set to the sum of the initializer values. If the wrong number of
    01981             values is given, a run-time error results. It is, however, possible to give
    01982             just one initializer. This creates an averaging filter with the given constant:
    01983
    01984             \code
    01985             vigra::Kernel1D<float> average5x1;
    01986
    01987             average5x1.initExplicitly(-2, 2) = 1.0/5.0;
    01988             \endcode
    01989
    01990             Here, the norm is set to value*size().
    01991
    01992             <b> Preconditions:</b>
    01993
    01994             \code
    01995
    01996             1. left <= 0
    01997             2. right >= 0
    01998             3. the number of values in the initializer list
    01999                is 1 or equals the size of the kernel.
    02000             \endcode
    02001         */
    02002     Kernel1D & initExplicitly(int left, int right)
    02003     {
    02004         vigra_precondition(left <= 0,
    02005                      "Kernel1D::initExplicitly(): left border must be <= 0.");
    02006         vigra_precondition(right >= 0,
    02007                      "Kernel1D::initExplicitly(): right border must be >= 0.");
    02008
    02009         right_ = right;
    02010         left_ = left;
    02011
    02012         kernel_.resize(right - left + 1);
    02013
    02014         return *this;
    02015     }
    02016
    02017         /** Get iterator to center of kernel
    02018
    02019             Postconditions:
    02020             \code
    02021
    02022             center()[left()] ... center()[right()] are valid kernel positions
    02023             \endcode
    02024         */
    02025     iterator center()
    02026     {
    02027         return kernel_.begin() - left();
    02028     }
    02029
    02030     const_iterator center() const
    02031     {
    02032         return kernel_.begin() - left();
    02033     }
    02034
    02035         /** Access kernel value at specified location.
    02036
    02037             Preconditions:
    02038             \code
    02039
    02040             left() <= location <= right()
    02041             \endcode
    02042         */
    02043     reference operator[](int location)
    02044     {
    02045         return kernel_[location - left()];
    02046     }
    02047
    02048     const_reference operator[](int location) const
    02049     {
    02050         return kernel_[location - left()];
    02051     }
    02052
    02053         /** left border of kernel (inclusive), always <= 0
    02054         */
    02055     int left() const { return left_; }
    02056
    02057         /** right border of kernel (inclusive), always >= 0
    02058         */
    02059     int right() const { return right_; }
    02060
    02061         /** size of kernel (right() - left() + 1)
    02062         */
    02063     int size() const { return right_ - left_ + 1; }
    02064
    02065         /** current border treatment mode
    02066         */
    02067     BorderTreatmentMode borderTreatment() const
    02068     { return border_treatment_; }
    02069
    02070         /** Set border treatment mode.
    02071         */
    02072     void setBorderTreatment( BorderTreatmentMode new_mode)
    02073     { border_treatment_ = new_mode; }
    02074
    02075         /** norm of kernel
    02076         */
    02077     value_type norm() const { return norm_; }
    02078
    02079         /** set a new norm and normalize kernel, use the normalization formula
    02080             for the given <tt>derivativeOrder</tt>.
    02081         */
    02082     void
    02083     normalize(value_type norm, unsigned int derivativeOrder = 0, double offset = 0.0);
    02084
    02085         /** normalize kernel to norm 1.
    02086         */
    02087     void
    02088     normalize()
    02089     {
    02090         normalize(one());
    02091     }
    02092
    02093         /** get a const accessor
    02094         */
    02095     ConstAccessor accessor() const { return ConstAccessor(); }
    02096
    02097         /** get an accessor
    02098         */
    02099     Accessor accessor() { return Accessor(); }
    02100
    02101   private:
    02102     InternalVector kernel_;
    02103     int left_, right_;
    02104     BorderTreatmentMode border_treatment_;
    02105     value_type norm_;
    02106 };
    02107
    02108 template <class ARITHTYPE>
    02109 void Kernel1D<ARITHTYPE>::normalize(value_type norm,
    02110                           unsigned int derivativeOrder,
    02111                           double offset)
    02112 {
    02113     typedef typename NumericTraits<value_type>::RealPromote TmpType;
    02114
    02115     // find kernel sum
    02116     Iterator k = kernel_.begin();
    02117     TmpType sum = NumericTraits<TmpType>::zero();
    02118
    02119     if(derivativeOrder == 0)
    02120     {
    02121         for(; k < kernel_.end(); ++k)
    02122         {
    02123             sum += *k;
    02124         }
    02125     }
    02126     else
    02127     {
    02128         unsigned int faculty = 1;
    02129         for(unsigned int i = 2; i <= derivativeOrder; ++i)
    02130             faculty *= i;
    02131         for(double x = left() + offset; k < kernel_.end(); ++x, ++k)
    02132         {
    02133             sum = TmpType(sum + *k * VIGRA_CSTD::pow(-x, int(derivativeOrder)) / faculty);
    02134         }
    02135     }
    02136
    02137     vigra_precondition(sum != NumericTraits<value_type>::zero(),
    02138                     "Kernel1D<ARITHTYPE>::normalize(): "
    02139                     "Cannot normalize a kernel with sum = 0");
    02140     // normalize
    02141     sum = norm / sum;
    02142     k = kernel_.begin();
    02143     for(; k != kernel_.end(); ++k)
    02144     {
    02145         *k = *k * sum;
    02146     }
    02147
    02148     norm_ = norm;
    02149 }
    02150
    02151 /***********************************************************************/
    02152
    02153 template <class ARITHTYPE>
    02154 void Kernel1D<ARITHTYPE>::initGaussian(double std_dev,
    02155                                        value_type norm,
    02156                                        double windowRatio)
    02157 {
    02158     vigra_precondition(std_dev >= 0.0,
    02159               "Kernel1D::initGaussian(): Standard deviation must be >= 0.");
    02160     vigra_precondition(windowRatio >= 0.0,
    02161               "Kernel1D::initGaussian(): windowRatio must be >= 0.");
    02162
    02163     if(std_dev > 0.0)
    02164     {
    02165         Gaussian<ARITHTYPE> gauss((ARITHTYPE)std_dev);
    02166
    02167         // first calculate required kernel sizes
    02168         int radius;
    02169         if (windowRatio == 0.0)
    02170             radius = (int)(3.0 * std_dev + 0.5);
    02171         else
    02172             radius = (int)(windowRatio * std_dev + 0.5);
    02173         if(radius == 0)
    02174             radius = 1;
    02175
    02176         // allocate the kernel
    02177         kernel_.erase(kernel_.begin(), kernel_.end());
    02178         kernel_.reserve(radius*2+1);
    02179
    02180         for(ARITHTYPE x = -(ARITHTYPE)radius; x <= (ARITHTYPE)radius; ++x)
    02181         {
    02182             kernel_.push_back(gauss(x));
    02183         }
    02184         left_ = -radius;
    02185         right_ = radius;
    02186     }
    02187     else
    02188     {
    02189         kernel_.erase(kernel_.begin(), kernel_.end());
    02190         kernel_.push_back(1.0);
    02191         left_ = 0;
    02192         right_ = 0;
    02193     }
    02194
    02195     if(norm != 0.0)
    02196         normalize(norm);
    02197     else
    02198         norm_ = 1.0;
    02199
    02200     // best border treatment for Gaussians is BORDER_TREATMENT_REFLECT
    02201     border_treatment_ = BORDER_TREATMENT_REFLECT;
    02202 }
    02203
    02204 /***********************************************************************/
    02205
    02206 template <class ARITHTYPE>
    02207 void Kernel1D<ARITHTYPE>::initDiscreteGaussian(double std_dev,
    02208                                        value_type norm)
    02209 {
    02210     vigra_precondition(std_dev >= 0.0,
    02211               "Kernel1D::initDiscreteGaussian(): Standard deviation must be >= 0.");
    02212
    02213     if(std_dev > 0.0)
    02214     {
    02215         // first calculate required kernel sizes
    02216         int radius = (int)(3.0*std_dev + 0.5);
    02217         if(radius == 0)
    02218             radius = 1;
    02219
    02220         double f = 2.0 / std_dev / std_dev;
    02221
    02222         // allocate the working array
    02223         int maxIndex = (int)(2.0 * (radius + 5.0 * VIGRA_CSTD::sqrt((double)radius)) + 0.5);
    02224         InternalVector warray(maxIndex+1);
    02225         warray[maxIndex] = 0.0;
    02226         warray[maxIndex-1] = 1.0;
    02227
    02228         for(int i = maxIndex-2; i >= radius; --i)
    02229         {
    02230             warray[i] = warray[i+2] + f * (i+1) * warray[i+1];
    02231             if(warray[i] > 1.0e40)
    02232             {
    02233                 warray[i+1] /= warray[i];
    02234                 warray[i] = 1.0;
    02235             }
    02236         }
    02237
    02238         // the following rescaling ensures that the numbers stay in a sensible range
    02239         // during the rest of the iteration, so no other rescaling is needed
    02240         double er = VIGRA_CSTD::exp(-radius*radius / (2.0*std_dev*std_dev));
    02241         warray[radius+1] = er * warray[radius+1] / warray[radius];
    02242         warray[radius] = er;
    02243
    02244         for(int i = radius-1; i >= 0; --i)
    02245         {
    02246             warray[i] = warray[i+2] + f * (i+1) * warray[i+1];
    02247             er += warray[i];
    02248         }
    02249
    02250         double scale = norm / (2*er - warray[0]);
    02251
    02252         initExplicitly(-radius, radius);
    02253         iterator c = center();
    02254
    02255         for(int i=0; i<=radius; ++i)
    02256         {
    02257             c[i] = c[-i] = warray[i] * scale;
    02258         }
    02259     }
    02260     else
    02261     {
    02262         kernel_.erase(kernel_.begin(), kernel_.end());
    02263         kernel_.push_back(norm);
    02264         left_ = 0;
    02265         right_ = 0;
    02266     }
    02267
    02268     norm_ = norm;
    02269
    02270     // best border treatment for Gaussians is BORDER_TREATMENT_REFLECT
    02271     border_treatment_ = BORDER_TREATMENT_REFLECT;
    02272 }
    02273
    02274 /***********************************************************************/
    02275
    02276 template <class ARITHTYPE>
    02277 void
    02278 Kernel1D<ARITHTYPE>::initGaussianDerivative(double std_dev,
    02279                                             int order,
    02280                                             value_type norm,
    02281                                             double windowRatio)
    02282 {
    02283     vigra_precondition(order >= 0,
    02284               "Kernel1D::initGaussianDerivative(): Order must be >= 0.");
    02285
    02286     if(order == 0)
    02287     {
    02288         initGaussian(std_dev, norm, windowRatio);
    02289         return;
    02290     }
    02291
    02292     vigra_precondition(std_dev > 0.0,
    02293               "Kernel1D::initGaussianDerivative(): "
    02294               "Standard deviation must be > 0.");
    02295     vigra_precondition(windowRatio >= 0.0,
    02296               "Kernel1D::initGaussianDerivative(): windowRatio must be >= 0.");
    02297
    02298     Gaussian<ARITHTYPE> gauss((ARITHTYPE)std_dev, order);
    02299
    02300     // first calculate required kernel sizes
    02301     int radius;
    02302     if(windowRatio == 0.0)
    02303         radius = (int)(3.0 * std_dev + 0.5 * order + 0.5);
    02304     else
    02305         radius = (int)(windowRatio * std_dev + 0.5);
    02306     if(radius == 0)
    02307         radius = 1;
    02308
    02309     // allocate the kernels
    02310     kernel_.clear();
    02311     kernel_.reserve(radius*2+1);
    02312
    02313     // fill the kernel and calculate the DC component
    02314     // introduced by truncation of the Gaussian
    02315     ARITHTYPE dc = 0.0;
    02316     for(ARITHTYPE x = -(ARITHTYPE)radius; x <= (ARITHTYPE)radius; ++x)
    02317     {
    02318         kernel_.push_back(gauss(x));
    02319         dc += kernel_[kernel_.size()-1];
    02320     }
    02321     dc = ARITHTYPE(dc / (2.0*radius + 1.0));
    02322
    02323     // remove DC, but only if kernel correction is permitted by a non-zero
    02324     // value for norm
    02325     if(norm != 0.0)
    02326     {
    02327         for(unsigned int i=0; i < kernel_.size(); ++i)
    02328         {
    02329             kernel_[i] -= dc;
    02330         }
    02331     }
    02332
    02333     left_ = -radius;
    02334     right_ = radius;
    02335
    02336     if(norm != 0.0)
    02337         normalize(norm, order);
    02338     else
    02339         norm_ = 1.0;
    02340
    02341     // best border treatment for Gaussian derivatives is
    02342     // BORDER_TREATMENT_REFLECT
    02343     border_treatment_ = BORDER_TREATMENT_REFLECT;
    02344 }
    02345
    02346 /***********************************************************************/
    02347
    02348 template <class ARITHTYPE>
    02349 void
    02350 Kernel1D<ARITHTYPE>::initBinomial(int radius,
    02351                                   value_type norm)
    02352 {
    02353     vigra_precondition(radius > 0,
    02354               "Kernel1D::initBinomial(): Radius must be > 0.");
    02355
    02356     // allocate the kernel
    02357     InternalVector(radius*2+1).swap(kernel_);
    02358     typename InternalVector::iterator x = kernel_.begin() + radius;
    02359
    02360     // fill kernel
    02361     x[radius] = norm;
    02362     for(int j=radius-1; j>=-radius; --j)
    02363     {
    02364         x[j] = 0.5 * x[j+1];
    02365         for(int i=j+1; i<radius; ++i)
    02366         {
    02367             x[i] = 0.5 * (x[i] + x[i+1]);
    02368         }
    02369         x[radius] *= 0.5;
    02370     }
    02371
    02372     left_ = -radius;
    02373     right_ = radius;
    02374     norm_ = norm;
    02375
    02376     // best border treatment for Binomial is BORDER_TREATMENT_REFLECT
    02377     border_treatment_ = BORDER_TREATMENT_REFLECT;
    02378 }
    02379
    02380 /***********************************************************************/
    02381
    02382 template <class ARITHTYPE>
    02383 void Kernel1D<ARITHTYPE>::initAveraging(int radius,
    02384                                         value_type norm)
    02385 {
    02386     vigra_precondition(radius > 0,
    02387               "Kernel1D::initAveraging(): Radius must be > 0.");
    02388
    02389     // calculate scaling
    02390     double scale = 1.0 / (radius * 2 + 1);
    02391
    02392     // normalize
    02393     kernel_.erase(kernel_.begin(), kernel_.end());
    02394     kernel_.reserve(radius*2+1);
    02395
    02396     for(int i=0; i<=radius*2+1; ++i)
    02397     {
    02398         kernel_.push_back(scale * norm);
    02399     }
    02400
    02401     left_ = -radius;
    02402     right_ = radius;
    02403     norm_ = norm;
    02404
    02405     // best border treatment for Averaging is BORDER_TREATMENT_CLIP
    02406     border_treatment_ = BORDER_TREATMENT_CLIP;
    02407 }
    02408
    02409 /***********************************************************************/
    02410
    02411 template <class ARITHTYPE>
    02412 void
    02413 Kernel1D<ARITHTYPE>::initSymmetricDifference(value_type norm)
    02414 {
    02415     kernel_.erase(kernel_.begin(), kernel_.end());
    02416     kernel_.reserve(3);
    02417
    02418     kernel_.push_back(ARITHTYPE(0.5 * norm));
    02419     kernel_.push_back(ARITHTYPE(0.0 * norm));
    02420     kernel_.push_back(ARITHTYPE(-0.5 * norm));
    02421
    02422     left_ = -1;
    02423     right_ = 1;
    02424     norm_ = norm;
    02425
    02426     // best border treatment for symmetric difference is
    02427     // BORDER_TREATMENT_REFLECT
    02428     border_treatment_ = BORDER_TREATMENT_REFLECT;
    02429 }
    02430
    02431 /**************************************************************/
    02432 /*                                                            */
    02433 /*         Argument object factories for Kernel1D             */
    02434 /*                                                            */
    02435 /*     (documentation: see vigra/convolution.hxx)             */
    02436 /*                                                            */
    02437 /**************************************************************/
    02438
    02439 template <class KernelIterator, class KernelAccessor>
    02440 inline
    02441 tuple5<KernelIterator, KernelAccessor, int, int, BorderTreatmentMode>
    02442 kernel1d(KernelIterator ik, KernelAccessor ka,
    02443        int kleft, int kright, BorderTreatmentMode border)
    02444 {
    02445     return
    02446       tuple5<KernelIterator, KernelAccessor, int, int, BorderTreatmentMode>(
    02447                                                 ik, ka, kleft, kright, border);
    02448 }
    02449
    02450 template <class T>
    02451 inline
    02452 tuple5<typename Kernel1D<T>::const_iterator, typename Kernel1D<T>::ConstAccessor,
    02453        int, int, BorderTreatmentMode>
    02454 kernel1d(Kernel1D<T> const & k)
    02455
    02456 {
    02457     return
    02458         tuple5<typename Kernel1D<T>::const_iterator, typename Kernel1D<T>::ConstAccessor,
    02459                int, int, BorderTreatmentMode>(
    02460                                      k.center(),
    02461                                      k.accessor(),
    02462                                      k.left(), k.right(),
    02463                                      k.borderTreatment());
    02464 }
    02465
    02466 template <class T>
    02467 inline
    02468 tuple5<typename Kernel1D<T>::const_iterator, typename Kernel1D<T>::ConstAccessor,
    02469        int, int, BorderTreatmentMode>
    02470 kernel1d(Kernel1D<T> const & k, BorderTreatmentMode border)
    02471
    02472 {
    02473     return
    02474         tuple5<typename Kernel1D<T>::const_iterator, typename Kernel1D<T>::ConstAccessor,
    02475                int, int, BorderTreatmentMode>(
    02476                                      k.center(),
    02477                                      k.accessor(),
    02478                                      k.left(), k.right(),
    02479                                      border);
    02480 }
    02481
    02482
    02483 } // namespace vigra
    02484
    02485 #endif // VIGRA_SEPARABLECONVOLUTION_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.9.0 (Tue Nov 6 2012)