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

vigra/separableconvolution.hxx VIGRA

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)