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

vigra/flatmorphology.hxx

00001 /************************************************************************/
00002 /*                                                                      */
00003 /*               Copyright 1998-2002 by Ullrich Koethe                  */
00004 /*       Cognitive Systems Group, University of Hamburg, Germany        */
00005 /*                                                                      */
00006 /*    This file is part of the VIGRA computer vision library.           */
00007 /*    The VIGRA Website is                                              */
00008 /*        http://kogs-www.informatik.uni-hamburg.de/~koethe/vigra/      */
00009 /*    Please direct questions, bug reports, and contributions to        */
00010 /*        ullrich.koethe@iwr.uni-heidelberg.de    or                    */
00011 /*        vigra@informatik.uni-hamburg.de                               */
00012 /*                                                                      */
00013 /*    Permission is hereby granted, free of charge, to any person       */
00014 /*    obtaining a copy of this software and associated documentation    */
00015 /*    files (the "Software"), to deal in the Software without           */
00016 /*    restriction, including without limitation the rights to use,      */
00017 /*    copy, modify, merge, publish, distribute, sublicense, and/or      */
00018 /*    sell copies of the Software, and to permit persons to whom the    */
00019 /*    Software is furnished to do so, subject to the following          */
00020 /*    conditions:                                                       */
00021 /*                                                                      */
00022 /*    The above copyright notice and this permission notice shall be    */
00023 /*    included in all copies or substantial portions of the             */
00024 /*    Software.                                                         */
00025 /*                                                                      */
00026 /*    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND    */
00027 /*    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   */
00028 /*    OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          */
00029 /*    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT       */
00030 /*    HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,      */
00031 /*    WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING      */
00032 /*    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR     */
00033 /*    OTHER DEALINGS IN THE SOFTWARE.                                   */                
00034 /*                                                                      */
00035 /************************************************************************/
00036  
00037  
00038 #ifndef VIGRA_FLATMORPHOLOGY_HXX
00039 #define VIGRA_FLATMORPHOLOGY_HXX
00040 
00041 #include <cmath>
00042 #include <vector>
00043 #include "utilities.hxx"
00044 
00045 namespace vigra {
00046 
00047 /** \addtogroup Morphology Basic Morphological Operations
00048     Perform erosion, dilation, and median with disc structuring functions
00049     
00050     See also: \ref MultiArrayMorphology Separable morphology with parabola structuring functions in arbitrary dimensions
00051 */
00052 //@{
00053 
00054 /********************************************************/
00055 /*                                                      */
00056 /*                  discRankOrderFilter                 */
00057 /*                                                      */
00058 /********************************************************/
00059 
00060 /** \brief Apply rank order filter with disc structuring function to the image.
00061 
00062     The pixel values of the source image <b> must</b> be in the range
00063     0...255. Radius must be >= 0. Rank must be in the range 0.0 <= rank 
00064     <= 1.0. The filter acts as a minimum filter if rank = 0.0, 
00065     as a median if rank = 0.5, and as a maximum filter if rank = 1.0.
00066     Accessor are used to access the pixel data.
00067     
00068     <b> Declarations:</b>
00069     
00070     pass arguments explicitely:
00071     \code
00072     namespace vigra {
00073         template <class SrcIterator, class SrcAccessor,
00074                   class DestIterator, class DestAccessor>
00075         void
00076         discRankOrderFilter(SrcIterator upperleft1, 
00077                             SrcIterator lowerright1, SrcAccessor sa,
00078                             DestIterator upperleft2, DestAccessor da,
00079                             int radius, float rank)
00080     }
00081     \endcode
00082     
00083     
00084     use argument objects in conjunction with \ref ArgumentObjectFactories :
00085     \code
00086     namespace vigra {
00087         template <class SrcIterator, class SrcAccessor,
00088                   class DestIterator, class DestAccessor>
00089         void
00090         discRankOrderFilter(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00091                             pair<DestIterator, DestAccessor> dest,
00092                             int radius, float rank)
00093     }
00094     \endcode
00095     
00096     <b> Usage:</b>
00097     
00098         <b>\#include</b> <<a href="flatmorphology_8hxx-source.html">vigra/flatmorphology.hxx</a>><br>
00099     Namespace: vigra
00100     
00101     \code
00102     vigra::CImage src, dest;
00103     
00104     // do median filtering
00105     vigra::discRankOrderFilter(srcImageRange(src), destImage(dest), 10, 0.5);
00106     \endcode
00107 
00108     <b> Required Interface:</b>
00109     
00110     \code
00111     SrcIterator src_upperleft;
00112     DestIterator dest_upperleft;
00113     int x, y;
00114     unsigned char value;
00115     
00116     SrcAccessor src_accessor;
00117     DestAccessor dest_accessor;
00118     
00119     // value_type of accessor must be convertible to unsigned char
00120     value = src_accessor(src_upperleft, x, y); 
00121     
00122     dest_accessor.set(value, dest_upperleft, x, y);
00123     \endcode
00124     
00125     <b> Preconditions:</b>
00126     
00127     \code
00128     for all source pixels: 0 <= value <= 255
00129     
00130     (rank >= 0.0) && (rank <= 1.0)
00131     radius >= 0
00132     \endcode
00133     
00134 */
00135 doxygen_overloaded_function(template <...> void discRankOrderFilter)
00136 
00137 template <class SrcIterator, class SrcAccessor,
00138           class DestIterator, class DestAccessor>
00139 void
00140 discRankOrderFilter(SrcIterator upperleft1, 
00141                     SrcIterator lowerright1, SrcAccessor sa,
00142                     DestIterator upperleft2, DestAccessor da,
00143                     int radius, float rank)
00144 {
00145     vigra_precondition((rank >= 0.0) && (rank <= 1.0),
00146             "discRankOrderFilter(): Rank must be between 0 and 1"
00147             " (inclusive).");
00148     
00149     vigra_precondition(radius >= 0,
00150             "discRankOrderFilter(): Radius must be >= 0.");
00151     
00152     int i, x, y, xmax, ymax, xx, yy;
00153     int rankpos, winsize, leftsum;
00154     
00155     long hist[256];
00156     
00157     // prepare structuring function
00158     std::vector<int> struct_function(radius+1);
00159     struct_function[0] = radius;
00160     
00161     double r2 = (double)radius*radius;
00162     for(i=1; i<=radius; ++i)
00163     {
00164         double r = (double) i - 0.5;
00165         struct_function[i] = (int)(VIGRA_CSTD::sqrt(r2 - r*r) + 0.5);
00166     }
00167 
00168     int w = lowerright1.x - upperleft1.x;
00169     int h = lowerright1.y - upperleft1.y;
00170     
00171     SrcIterator ys(upperleft1);
00172     DestIterator yd(upperleft2);
00173 
00174     for(y=0; y<h; ++y, ++ys.y, ++yd.y)
00175     {
00176         SrcIterator xs(ys);
00177         DestIterator xd(yd);
00178         
00179         // first column
00180         int x0 = 0;
00181         int y0 = y;
00182         int x1 = w - 1;
00183         int y1 = h - y - 1;
00184 
00185         // clear histogram
00186         for(i=0; i<256; ++i) hist[i] = 0;
00187         winsize = 0;
00188         
00189         // init histogram
00190         ymax = (y1 < radius) ? y1 : radius;
00191         for(yy=0; yy<=ymax; ++yy)
00192         {
00193             xmax = (x1 < struct_function[yy]) ? x1 : struct_function[yy];
00194             for(xx=0; xx<=xmax; ++xx)
00195             {
00196                 hist[sa(xs, Diff2D(xx, yy))]++;
00197                 winsize++;
00198             }
00199         }
00200         
00201         ymax = (y0 < radius) ? y0 : radius;
00202         for(yy=1; yy<=ymax; ++yy)
00203         {
00204             xmax = (x1 < struct_function[yy]) ? x1 : struct_function[yy];
00205             for(xx=0; xx<=xmax; ++xx)
00206             {
00207                 hist[sa(xs, Diff2D(xx, -yy))]++;
00208                 winsize++;
00209             }
00210         }
00211     
00212         // find the desired histogramm bin 
00213         leftsum = 0;
00214         if(rank == 0.0)
00215         {
00216             for(i=0; i<256; i++)
00217             {
00218                 if(hist[i]) break;
00219             }
00220             rankpos = i;
00221         }
00222         else
00223         {
00224             for(i=0; i<256; i++)
00225             {
00226                 if((float)(hist[i]+leftsum) / winsize >= rank) break;
00227                 leftsum += hist[i];
00228             }
00229             rankpos = i;
00230         }
00231         
00232         da.set(rankpos, xd);
00233         
00234         ++xs.x;
00235         ++xd.x;
00236         
00237         // inner columns
00238         for(x=1; x<w; ++x, ++xs.x, ++xd.x)
00239         {
00240             x0 = x;
00241             y0 = y;
00242             x1 = w - x - 1;
00243             y1 = h - y - 1;
00244             
00245             // update histogramm 
00246             // remove pixels at left border 
00247             yy = (y1 < radius) ? y1 : radius;
00248             for(; yy>=0; yy--)
00249             {
00250                 unsigned char cur;
00251                 xx = struct_function[yy]+1;
00252                 if(xx > x0) break;
00253                 
00254                 cur = sa(xs, Diff2D(-xx, yy));
00255                 
00256                 hist[cur]--;
00257                 if(cur < rankpos) leftsum--;
00258                 winsize--;
00259             }
00260             yy = (y0 < radius) ? y0 : radius;
00261             for(; yy>=1; yy--)
00262             {
00263                 unsigned char cur;
00264                 xx = struct_function[yy]+1;
00265                 if(xx > x0) break;
00266                 
00267                 cur = sa(xs, Diff2D(-xx, -yy));
00268                 
00269                 hist[cur]--;
00270                 if(cur < rankpos) leftsum--;
00271                 winsize--;
00272             }
00273             
00274             // add pixels at right border 
00275             yy = (y1 < radius) ? y1 : radius;
00276             for(; yy>=0; yy--)
00277             {
00278                 unsigned char cur;
00279                 xx = struct_function[yy];
00280                 if(xx > x1) break;
00281                 
00282                 cur = sa(xs, Diff2D(xx, yy));
00283                 
00284                 hist[cur]++;
00285                 if(cur < rankpos) leftsum++;
00286                 winsize++;
00287             }
00288             yy = (y0 < radius) ? y0 : radius;
00289             for(; yy>=1; yy--)
00290             {
00291                 unsigned char cur;
00292                 xx = struct_function[yy];
00293                 if(xx > x1) break;
00294                 
00295                 cur = sa(xs, Diff2D(xx, -yy));
00296                 
00297                 hist[cur]++;
00298                 if(cur < rankpos) leftsum++;
00299                 winsize++;
00300             }
00301         
00302             // find the desired histogramm bin 
00303             if(rank == 0.0)
00304             {
00305                 if(leftsum == 0)
00306                 {
00307                     // search to the right 
00308                     for(i=rankpos; i<256; i++)
00309                     {
00310                         if(hist[i]) break;
00311                     }
00312                     rankpos = i;
00313                 }
00314                 else
00315                 {
00316                     // search to the left 
00317                     for(i=rankpos-1; i>=0; i--)
00318                     {
00319                         leftsum -= hist[i];
00320                         if(leftsum == 0) break;
00321                     }
00322                     rankpos = i;
00323                 }
00324             }
00325             else  // rank > 0.0 
00326             {
00327                 if((float)leftsum / winsize < rank)
00328                 {
00329                     // search to the right 
00330                     for(i=rankpos; i<256; i++)
00331                     {
00332                         if((float)(hist[i]+leftsum) / winsize >= rank) break;
00333                         leftsum+=hist[i];
00334                     }
00335                     rankpos = i;
00336                 }
00337                 else
00338                 {
00339                     /// search to the left 
00340                     for(i=rankpos-1; i>=0; i--)
00341                     {
00342                         leftsum-=hist[i];
00343                         if((float)leftsum / winsize < rank) break;
00344                     }
00345                     rankpos = i;
00346                 }
00347             }
00348                     
00349             da.set(rankpos, xd);
00350         }
00351     }
00352 }
00353 
00354 template <class SrcIterator, class SrcAccessor,
00355           class DestIterator, class DestAccessor>
00356 void
00357 discRankOrderFilter(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00358                     pair<DestIterator, DestAccessor> dest,
00359                     int radius, float rank)
00360 {
00361     discRankOrderFilter(src.first, src.second, src.third,
00362                         dest.first, dest.second,
00363                         radius, rank);
00364 }
00365 
00366 /********************************************************/
00367 /*                                                      */
00368 /*                      discErosion                     */
00369 /*                                                      */
00370 /********************************************************/
00371 
00372 /** \brief Apply erosion (minimum) filter with disc of given radius to image.
00373 
00374     This is an abbreviation vor the rank order filter with rank = 0.0.
00375     See \ref discRankOrderFilter() for more information.
00376     
00377     <b> Declarations:</b>
00378     
00379     pass arguments explicitely:
00380     \code
00381     namespace vigra {
00382         template <class SrcIterator, class SrcAccessor,
00383                   class DestIterator, class DestAccessor>
00384         void 
00385         discErosion(SrcIterator upperleft1, 
00386                     SrcIterator lowerright1, SrcAccessor sa,
00387                     DestIterator upperleft2, DestAccessor da,
00388                     int radius)
00389     }
00390     \endcode
00391     
00392     
00393     use argument objects in conjunction with \ref ArgumentObjectFactories :
00394     \code
00395     namespace vigra {
00396         template <class SrcIterator, class SrcAccessor,
00397                   class DestIterator, class DestAccessor>
00398         void
00399         discErosion(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00400                     pair<DestIterator, DestAccessor> dest,
00401                     int radius)
00402     }
00403     \endcode
00404 
00405 */
00406 doxygen_overloaded_function(template <...> void discErosion)
00407 
00408 template <class SrcIterator, class SrcAccessor,
00409           class DestIterator, class DestAccessor>
00410 inline void 
00411 discErosion(SrcIterator upperleft1, 
00412             SrcIterator lowerright1, SrcAccessor sa,
00413             DestIterator upperleft2, DestAccessor da,
00414             int radius)
00415 {
00416     vigra_precondition(radius >= 0, "discErosion(): Radius must be >= 0.");
00417     
00418     discRankOrderFilter(upperleft1, lowerright1, sa, 
00419                         upperleft2, da, radius, 0.0);
00420 }
00421 
00422 template <class SrcIterator, class SrcAccessor,
00423           class DestIterator, class DestAccessor>
00424 void
00425 discErosion(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00426             pair<DestIterator, DestAccessor> dest,
00427             int radius)
00428 {
00429     vigra_precondition(radius >= 0, "discErosion(): Radius must be >= 0.");
00430     
00431     discRankOrderFilter(src.first, src.second, src.third,
00432                         dest.first, dest.second,
00433                         radius, 0.0);
00434 }
00435 
00436 /********************************************************/
00437 /*                                                      */
00438 /*                     discDilation                     */
00439 /*                                                      */
00440 /********************************************************/
00441 
00442 /** \brief Apply dilation (maximum) filter with disc of given radius to image.
00443 
00444     This is an abbreviation vor the rank order filter with rank = 1.0.
00445     See \ref discRankOrderFilter() for more information.
00446     
00447     <b> Declarations:</b>
00448     
00449     pass arguments explicitely:
00450     \code
00451     namespace vigra {
00452         template <class SrcIterator, class SrcAccessor,
00453                   class DestIterator, class DestAccessor>
00454         void 
00455         discDilation(SrcIterator upperleft1, 
00456                     SrcIterator lowerright1, SrcAccessor sa,
00457                     DestIterator upperleft2, DestAccessor da,
00458                     int radius)
00459     }
00460     \endcode
00461     
00462     
00463     use argument objects in conjunction with \ref ArgumentObjectFactories :
00464     \code
00465     namespace vigra {
00466         template <class SrcIterator, class SrcAccessor,
00467                   class DestIterator, class DestAccessor>
00468         void
00469         discDilation(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00470                     pair<DestIterator, DestAccessor> dest,
00471                     int radius)
00472     }
00473     \endcode
00474 
00475 */
00476 doxygen_overloaded_function(template <...> void discDilation)
00477 
00478 template <class SrcIterator, class SrcAccessor,
00479           class DestIterator, class DestAccessor>
00480 inline void 
00481 discDilation(SrcIterator upperleft1, 
00482             SrcIterator lowerright1, SrcAccessor sa,
00483             DestIterator upperleft2, DestAccessor da,
00484             int radius)
00485 {
00486     vigra_precondition(radius >= 0, "discDilation(): Radius must be >= 0.");
00487     
00488     discRankOrderFilter(upperleft1, lowerright1, sa, 
00489                         upperleft2, da, radius, 1.0);
00490 }
00491 
00492 template <class SrcIterator, class SrcAccessor,
00493           class DestIterator, class DestAccessor>
00494 void
00495 discDilation(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00496             pair<DestIterator, DestAccessor> dest,
00497             int radius)
00498 {
00499     vigra_precondition(radius >= 0, "discDilation(): Radius must be >= 0.");
00500     
00501     discRankOrderFilter(src.first, src.second, src.third,
00502                         dest.first, dest.second,
00503                         radius, 1.0);
00504 }
00505 
00506 /********************************************************/
00507 /*                                                      */
00508 /*                      discMedian                      */
00509 /*                                                      */
00510 /********************************************************/
00511 
00512 /** \brief Apply median filter with disc of given radius to image.
00513 
00514     This is an abbreviation vor the rank order filter with rank = 0.5.
00515     See \ref discRankOrderFilter() for more information.
00516     
00517     <b> Declarations:</b>
00518     
00519     pass arguments explicitely:
00520     \code
00521     namespace vigra {
00522         template <class SrcIterator, class SrcAccessor,
00523                   class DestIterator, class DestAccessor>
00524         void 
00525         discMedian(SrcIterator upperleft1, 
00526                     SrcIterator lowerright1, SrcAccessor sa,
00527                     DestIterator upperleft2, DestAccessor da,
00528                     int radius)
00529     }
00530     \endcode
00531     
00532     
00533     use argument objects in conjunction with \ref ArgumentObjectFactories :
00534     \code
00535     namespace vigra {
00536         template <class SrcIterator, class SrcAccessor,
00537                   class DestIterator, class DestAccessor>
00538         void
00539         discMedian(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00540                     pair<DestIterator, DestAccessor> dest,
00541                     int radius)
00542     }
00543     \endcode
00544 
00545 */
00546 doxygen_overloaded_function(template <...> void discMedian)
00547 
00548 template <class SrcIterator, class SrcAccessor,
00549           class DestIterator, class DestAccessor>
00550 inline void 
00551 discMedian(SrcIterator upperleft1, 
00552             SrcIterator lowerright1, SrcAccessor sa,
00553             DestIterator upperleft2, DestAccessor da,
00554             int radius)
00555 {
00556     vigra_precondition(radius >= 0, "discMedian(): Radius must be >= 0.");
00557     
00558     discRankOrderFilter(upperleft1, lowerright1, sa, 
00559                         upperleft2, da, radius, 0.5);
00560 }
00561 
00562 template <class SrcIterator, class SrcAccessor,
00563           class DestIterator, class DestAccessor>
00564 void
00565 discMedian(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00566             pair<DestIterator, DestAccessor> dest,
00567             int radius)
00568 {
00569     vigra_precondition(radius >= 0, "discMedian(): Radius must be >= 0.");
00570     
00571     discRankOrderFilter(src.first, src.second, src.third,
00572                         dest.first, dest.second,
00573                         radius, 0.5);
00574 }
00575 
00576 /********************************************************/
00577 /*                                                      */
00578 /*            discRankOrderFilterWithMask               */
00579 /*                                                      */
00580 /********************************************************/
00581 
00582 /** \brief Apply rank order filter with disc structuring function to the image
00583     using a mask.
00584     
00585     The pixel values of the source image <b> must</b> be in the range
00586     0...255. Radius must be >= 0. Rank must be in the range 0.0 <= rank 
00587     <= 1.0. The filter acts as a minimum filter if rank = 0.0, 
00588     as a median if rank = 0.5, and as a maximum filter if rank = 1.0.
00589     Accessor are used to access the pixel data.
00590     
00591     The mask is only applied to th input image, i.e. the function
00592     generates an output wherever the current disc contains at least 
00593     one pixel with mask value 'true'. Source pixels with mask value
00594     'false' are ignored during the calculation of the rank order.
00595     
00596     <b> Declarations:</b>
00597     
00598     pass arguments explicitely:
00599     \code
00600     namespace vigra {
00601         template <class SrcIterator, class SrcAccessor,
00602                   class MaskIterator, class MaskAccessor,
00603                   class DestIterator, class DestAccessor>
00604         void
00605         discRankOrderFilterWithMask(SrcIterator upperleft1, 
00606                                     SrcIterator lowerright1, SrcAccessor sa,
00607                                     MaskIterator upperleftm, MaskAccessor mask,
00608                                     DestIterator upperleft2, DestAccessor da,
00609                                     int radius, float rank)
00610     }
00611     \endcode
00612     
00613     
00614     group arguments (use in conjunction with \ref ArgumentObjectFactories):
00615     \code
00616     namespace vigra {
00617         template <class SrcIterator, class SrcAccessor,
00618                   class MaskIterator, class MaskAccessor,
00619                   class DestIterator, class DestAccessor>
00620         void
00621         discRankOrderFilterWithMask(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00622                                     pair<MaskIterator, MaskAccessor> mask,
00623                                     pair<DestIterator, DestAccessor> dest,
00624                                     int radius, float rank)
00625     }
00626     \endcode
00627     
00628     <b> Usage:</b>
00629     
00630     <b>\#include</b> <<a href="flatmorphology_8hxx-source.html">vigra/flatmorphology.hxx</a>><br>
00631     Namespace: vigra
00632     
00633     \code
00634     vigra::CImage src, dest, mask;
00635     
00636     // do median filtering
00637     vigra::discRankOrderFilterWithMask(srcImageRange(src), 
00638                                 maskImage(mask), destImage(dest), 10, 0.5);
00639     \endcode
00640 
00641     <b> Required Interface:</b>
00642     
00643     \code
00644     SrcIterator src_upperleft;
00645     DestIterator dest_upperleft;
00646     MaskIterator mask_upperleft;
00647     int x, y;
00648     unsigned char value;
00649     
00650     SrcAccessor src_accessor;
00651     DestAccessor dest_accessor;
00652     MaskAccessor mask_accessor;
00653                      
00654     mask_accessor(mask_upperleft, x, y) // convertible to bool
00655     
00656     // value_type of accessor must be convertible to unsigned char
00657     value = src_accessor(src_upperleft, x, y); 
00658     
00659     dest_accessor.set(value, dest_upperleft, x, y);
00660     \endcode
00661     
00662     <b> Preconditions:</b>
00663     
00664     \code
00665     for all source pixels: 0 <= value <= 255
00666     
00667     (rank >= 0.0) && (rank <= 1.0)
00668     radius >= 0
00669     \endcode
00670     
00671 */
00672 doxygen_overloaded_function(template <...> void discRankOrderFilterWithMask)
00673 
00674 template <class SrcIterator, class SrcAccessor,
00675           class MaskIterator, class MaskAccessor,
00676           class DestIterator, class DestAccessor>
00677 void
00678 discRankOrderFilterWithMask(SrcIterator upperleft1, 
00679                             SrcIterator lowerright1, SrcAccessor sa,
00680                             MaskIterator upperleftm, MaskAccessor mask,
00681                             DestIterator upperleft2, DestAccessor da,
00682                             int radius, float rank)
00683 {
00684     vigra_precondition((rank >= 0.0) && (rank <= 1.0),
00685                  "discRankOrderFilter(): Rank must be between 0 and 1"
00686                  " (inclusive).");
00687     
00688     vigra_precondition(radius >= 0, "discRankOrderFilter(): Radius must be >= 0.");
00689     
00690     int i, x, y, xmax, ymax, xx, yy;
00691     int rankpos, winsize, leftsum;
00692     
00693     long hist[256];
00694     
00695     // prepare structuring function
00696     std::vector<int> struct_function(radius+1);
00697     struct_function[0] = radius;
00698     
00699     double r2 = (double)radius*radius;
00700     for(i=1; i<=radius; ++i)
00701     {
00702         double r = (double) i - 0.5;
00703         struct_function[i] = (int)(VIGRA_CSTD::sqrt(r2 - r*r) + 0.5);
00704     }
00705 
00706     int w = lowerright1.x - upperleft1.x;
00707     int h = lowerright1.y - upperleft1.y;
00708         
00709     SrcIterator ys(upperleft1);
00710     MaskIterator ym(upperleftm);
00711     DestIterator yd(upperleft2);
00712 
00713     for(y=0; y<h; ++y, ++ys.y, ++yd.y, ++ym.y)
00714     {
00715         SrcIterator xs(ys);
00716         MaskIterator xm(ym);
00717         DestIterator xd(yd);
00718         
00719         // first column
00720         int x0 = 0;
00721         int y0 = y;
00722         int x1 = w - 1;
00723         int y1 = h - y - 1;
00724 
00725         // clear histogram
00726         for(i=0; i<256; ++i) hist[i] = 0;
00727         winsize = 0;
00728         leftsum = 0;
00729         rankpos = 0;
00730         
00731         // init histogram
00732         ymax = (y1 < radius) ? y1 : radius;
00733         for(yy=0; yy<=ymax; ++yy)
00734         {
00735             xmax = (x1 < struct_function[yy]) ? x1 : struct_function[yy];
00736             for(xx=0; xx<=xmax; ++xx)
00737             {
00738                 Diff2D pos(xx, yy);
00739                 if(mask(xm, pos))
00740                 {
00741                     hist[sa(xs, pos)]++;
00742                     winsize++;
00743                 }
00744             }
00745         }
00746         
00747         ymax = (y0 < radius) ? y0 : radius;
00748         for(yy=1; yy<=ymax; ++yy)
00749         {
00750             xmax = (x1 < struct_function[yy]) ? x1 : struct_function[yy];
00751             for(xx=0; xx<=xmax; ++xx)
00752             {
00753                 Diff2D pos(xx, -yy);
00754                 if(mask(xm, pos))
00755                 {
00756                     hist[sa(xs, pos)]++;
00757                     winsize++;
00758                 }
00759             }
00760         }
00761     
00762         // find the desired histogramm bin 
00763         if(winsize) 
00764         {
00765             if(rank == 0.0)
00766             {
00767                 for(i=0; i<256; i++)
00768                 {
00769                     if(hist[i]) break;
00770                 }
00771                 rankpos = i;
00772             }
00773             else
00774             {
00775                 for(i=0; i<256; i++)
00776                 {
00777                     if((float)(hist[i]+leftsum) / winsize >= rank) break;
00778                     leftsum += hist[i];
00779                 }
00780                 rankpos = i;
00781             }
00782             
00783             da.set(rankpos, xd);
00784         }
00785             
00786         ++xs.x;
00787         ++xd.x;
00788         ++xm.x;
00789         
00790         // inner columns
00791         for(x=1; x<w; ++x, ++xs.x, ++xd.x, ++xm.x)
00792         {
00793             x0 = x;
00794             y0 = y;
00795             x1 = w - x - 1;
00796             y1 = h - y - 1;
00797             
00798             // update histogramm 
00799             // remove pixels at left border 
00800             yy = (y1 < radius) ? y1 : radius;
00801             for(; yy>=0; yy--)
00802             {
00803                 unsigned char cur;
00804                 xx = struct_function[yy]+1;
00805                 if(xx > x0) break;
00806                 
00807                 Diff2D pos(-xx, yy);
00808                 if(mask(xm, pos))
00809                 {
00810                     cur = sa(xs, pos);
00811                     
00812                     hist[cur]--;
00813                     if(cur < rankpos) leftsum--;
00814                     winsize--;
00815                 }
00816             }
00817             yy = (y0 < radius) ? y0 : radius;
00818             for(; yy>=1; yy--)
00819             {
00820                 unsigned char cur;
00821                 xx = struct_function[yy]+1;
00822                 if(xx > x0) break;
00823                 
00824                 Diff2D pos(-xx, -yy);
00825                 if(mask(xm, pos))
00826                 {
00827                     cur = sa(xs, pos);
00828                     
00829                     hist[cur]--;
00830                     if(cur < rankpos) leftsum--;
00831                     winsize--;
00832                 }
00833             }
00834             
00835             // add pixels at right border 
00836             yy = (y1 < radius) ? y1 : radius;
00837             for(; yy>=0; yy--)
00838             {
00839                 unsigned char cur;
00840                 xx = struct_function[yy];
00841                 if(xx > x1) break;
00842                 
00843                 Diff2D pos(xx, yy);
00844                 if(mask(xm, pos))
00845                 {
00846                     cur = sa(xs, pos);
00847                     
00848                     hist[cur]++;
00849                     if(cur < rankpos) leftsum++;
00850                     winsize++;
00851                 }
00852             }
00853             yy = (y0 < radius) ? y0 : radius;
00854             for(; yy>=1; yy--)
00855             {
00856                 unsigned char cur;
00857                 xx = struct_function[yy];
00858                 if(xx > x1) break;
00859                 
00860                 Diff2D pos(xx, -yy);
00861                 if(mask(xm, pos))
00862                 {
00863                     cur = sa(xs, pos);
00864                     
00865                     hist[cur]++;
00866                     if(cur < rankpos) leftsum++;
00867                     winsize++;
00868                 }
00869             }
00870         
00871             // find the desired histogramm bin 
00872             if(winsize) 
00873             {
00874                 if(rank == 0.0)
00875                 {
00876                     if(leftsum == 0)
00877                     {
00878                         // search to the right 
00879                         for(i=rankpos; i<256; i++)
00880                         {
00881                             if(hist[i]) break;
00882                         }
00883                         rankpos = i;
00884                     }
00885                     else
00886                     {
00887                         // search to the left 
00888                         for(i=rankpos-1; i>=0; i--)
00889                         {
00890                             leftsum -= hist[i];
00891                             if(leftsum == 0) break;
00892                         }
00893                         rankpos = i;
00894                     }
00895                 }
00896                 else  // rank > 0.0 
00897                 {
00898                     if((float)leftsum / winsize < rank)
00899                     {
00900                         // search to the right 
00901                         for(i=rankpos; i<256; i++)
00902                         {
00903                             if((float)(hist[i]+leftsum) / winsize >= rank) break;
00904                             leftsum+=hist[i];
00905                         }
00906                         rankpos = i;
00907                     }
00908                     else
00909                     {
00910                         /// search to the left 
00911                         for(i=rankpos-1; i>=0; i--)
00912                         {
00913                             leftsum-=hist[i];
00914                             if((float)leftsum / winsize < rank) break;
00915                         }
00916                         rankpos = i;
00917                     }
00918                 }
00919                         
00920                 da.set(rankpos, xd);
00921             }
00922             else
00923             {
00924                 leftsum = 0;
00925                 rankpos = 0;
00926             }
00927         }
00928     }
00929 }
00930 
00931 template <class SrcIterator, class SrcAccessor,
00932           class MaskIterator, class MaskAccessor,
00933           class DestIterator, class DestAccessor>
00934 void
00935 discRankOrderFilterWithMask(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00936                             pair<MaskIterator, MaskAccessor> mask,
00937                             pair<DestIterator, DestAccessor> dest,
00938                             int radius, float rank)
00939 {
00940     discRankOrderFilterWithMask(src.first, src.second, src.third,
00941                         mask.first, mask.second,
00942                         dest.first, dest.second,
00943                         radius, rank);
00944 }
00945 
00946 /********************************************************/
00947 /*                                                      */
00948 /*                 discErosionWithMask                  */
00949 /*                                                      */
00950 /********************************************************/
00951 
00952 /** \brief Apply erosion (minimum) filter with disc of given radius to image
00953     using a mask.
00954     
00955     This is an abbreviation vor the masked rank order filter with 
00956     rank = 0.0. See \ref discRankOrderFilterWithMask() for more information.
00957     
00958     <b> Declarations:</b>
00959     
00960     pass arguments explicitely:
00961     \code
00962     namespace vigra {
00963         template <class SrcIterator, class SrcAccessor,
00964                   class MaskIterator, class MaskAccessor,
00965                   class DestIterator, class DestAccessor>
00966         void 
00967         discErosionWithMask(SrcIterator upperleft1, 
00968                             SrcIterator lowerright1, SrcAccessor sa,
00969                             MaskIterator upperleftm, MaskAccessor mask,
00970                             DestIterator upperleft2, DestAccessor da,
00971                             int radius)
00972     }
00973     \endcode
00974     
00975     
00976     group arguments (use in conjunction with \ref ArgumentObjectFactories):
00977     \code
00978     namespace vigra {
00979         template <class SrcIterator, class SrcAccessor,
00980                   class MaskIterator, class MaskAccessor,
00981                   class DestIterator, class DestAccessor>
00982         void 
00983         discErosionWithMask(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00984                             pair<MaskIterator, MaskAccessor> mask,
00985                             pair<DestIterator, DestAccessor> dest,
00986                             int radius)
00987     }
00988     \endcode
00989 
00990 */
00991 doxygen_overloaded_function(template <...> void discErosionWithMask)
00992 
00993 template <class SrcIterator, class SrcAccessor,
00994           class MaskIterator, class MaskAccessor,
00995           class DestIterator, class DestAccessor>
00996 inline void 
00997 discErosionWithMask(SrcIterator upperleft1, 
00998                     SrcIterator lowerright1, SrcAccessor sa,
00999                     MaskIterator upperleftm, MaskAccessor mask,
01000                     DestIterator upperleft2, DestAccessor da,
01001                     int radius)
01002 {
01003     vigra_precondition(radius >= 0, "discErosionWithMask(): Radius must be >= 0.");
01004     
01005     discRankOrderFilterWithMask(upperleft1, lowerright1, sa, 
01006                                 upperleftm, mask,
01007                                 upperleft2, da,
01008                                 radius, 0.0);
01009 }
01010 
01011 template <class SrcIterator, class SrcAccessor,
01012           class MaskIterator, class MaskAccessor,
01013           class DestIterator, class DestAccessor>
01014 inline void 
01015 discErosionWithMask(triple<SrcIterator, SrcIterator, SrcAccessor> src,
01016                     pair<MaskIterator, MaskAccessor> mask,
01017                     pair<DestIterator, DestAccessor> dest,
01018                     int radius)
01019 {
01020     vigra_precondition(radius >= 0, "discErosionWithMask(): Radius must be >= 0.");
01021     
01022     discRankOrderFilterWithMask(src.first, src.second, src.third,
01023                         mask.first, mask.second,
01024                         dest.first, dest.second,
01025                         radius, 0.0);
01026 }
01027 
01028 /********************************************************/
01029 /*                                                      */
01030 /*                discDilationWithMask                  */
01031 /*                                                      */
01032 /********************************************************/
01033 
01034 /** \brief Apply dilation (maximum) filter with disc of given radius to image
01035     using a mask.
01036     
01037     This is an abbreviation vor the masked rank order filter with 
01038     rank = 1.0. See \ref discRankOrderFilterWithMask() for more information.
01039     
01040     <b> Declarations:</b>
01041     
01042     pass arguments explicitely:
01043     \code
01044     namespace vigra {
01045         template <class SrcIterator, class SrcAccessor,
01046                   class MaskIterator, class MaskAccessor,
01047                   class DestIterator, class DestAccessor>
01048         void 
01049         discDilationWithMask(SrcIterator upperleft1, 
01050                             SrcIterator lowerright1, SrcAccessor sa,
01051                             MaskIterator upperleftm, MaskAccessor mask,
01052                             DestIterator upperleft2, DestAccessor da,
01053                             int radius)
01054     }
01055     \endcode
01056     
01057     
01058     group arguments (use in conjunction with \ref ArgumentObjectFactories):
01059     \code
01060     namespace vigra {
01061         template <class SrcIterator, class SrcAccessor,
01062                   class MaskIterator, class MaskAccessor,
01063                   class DestIterator, class DestAccessor>
01064         void 
01065         discDilationWithMask(triple<SrcIterator, SrcIterator, SrcAccessor> src,
01066                             pair<MaskIterator, MaskAccessor> mask,
01067                             pair<DestIterator, DestAccessor> dest,
01068                             int radius)
01069     }
01070     \endcode
01071 
01072 */
01073 doxygen_overloaded_function(template <...> void discDilationWithMask)
01074 
01075 template <class SrcIterator, class SrcAccessor,
01076           class MaskIterator, class MaskAccessor,
01077           class DestIterator, class DestAccessor>
01078 inline void 
01079 discDilationWithMask(SrcIterator upperleft1, 
01080                     SrcIterator lowerright1, SrcAccessor sa,
01081                     MaskIterator upperleftm, MaskAccessor mask,
01082                     DestIterator upperleft2, DestAccessor da,
01083                     int radius)
01084 {
01085     vigra_precondition(radius >= 0, "discDilationWithMask(): Radius must be >= 0.");
01086     
01087     discRankOrderFilterWithMask(upperleft1, lowerright1, sa, 
01088                                 upperleftm, mask,
01089                                 upperleft2, da,
01090                                 radius, 1.0);
01091 }
01092 
01093 template <class SrcIterator, class SrcAccessor,
01094           class MaskIterator, class MaskAccessor,
01095           class DestIterator, class DestAccessor>
01096 inline void 
01097 discDilationWithMask(triple<SrcIterator, SrcIterator, SrcAccessor> src,
01098                     pair<MaskIterator, MaskAccessor> mask,
01099                     pair<DestIterator, DestAccessor> dest,
01100                     int radius)
01101 {
01102     vigra_precondition(radius >= 0, "discDilationWithMask(): Radius must be >= 0.");
01103     
01104     discRankOrderFilterWithMask(src.first, src.second, src.third,
01105                         mask.first, mask.second,
01106                         dest.first, dest.second,
01107                         radius, 1.0);
01108 }
01109 
01110 /********************************************************/
01111 /*                                                      */
01112 /*                 discMedianWithMask                   */
01113 /*                                                      */
01114 /********************************************************/
01115 
01116 /** \brief Apply median filter with disc of given radius to image
01117     using a mask.
01118     
01119     This is an abbreviation vor the masked rank order filter with 
01120     rank = 0.5. See \ref discRankOrderFilterWithMask() for more information.
01121     
01122     <b> Declarations:</b>
01123     
01124     pass arguments explicitely:
01125     \code
01126     namespace vigra {
01127         template <class SrcIterator, class SrcAccessor,
01128                   class MaskIterator, class MaskAccessor,
01129                   class DestIterator, class DestAccessor>
01130         void 
01131         discMedianWithMask(SrcIterator upperleft1, 
01132                             SrcIterator lowerright1, SrcAccessor sa,
01133                             MaskIterator upperleftm, MaskAccessor mask,
01134                             DestIterator upperleft2, DestAccessor da,
01135                             int radius)
01136     }
01137     \endcode
01138     
01139     
01140     group arguments (use in conjunction with \ref ArgumentObjectFactories):
01141     \code
01142     namespace vigra {
01143         template <class SrcIterator, class SrcAccessor,
01144                   class MaskIterator, class MaskAccessor,
01145                   class DestIterator, class DestAccessor>
01146         void 
01147         discMedianWithMask(triple<SrcIterator, SrcIterator, SrcAccessor> src,
01148                             pair<MaskIterator, MaskAccessor> mask,
01149                             pair<DestIterator, DestAccessor> dest,
01150                             int radius)
01151     }
01152     \endcode
01153 
01154 */
01155 doxygen_overloaded_function(template <...> void discMedianWithMask)
01156 
01157 template <class SrcIterator, class SrcAccessor,
01158           class MaskIterator, class MaskAccessor,
01159           class DestIterator, class DestAccessor>
01160 inline void 
01161 discMedianWithMask(SrcIterator upperleft1, 
01162                     SrcIterator lowerright1, SrcAccessor sa,
01163                     MaskIterator upperleftm, MaskAccessor mask,
01164                     DestIterator upperleft2, DestAccessor da,
01165                     int radius)
01166 {
01167     vigra_precondition(radius >= 0, "discMedianWithMask(): Radius must be >= 0.");
01168     
01169     discRankOrderFilterWithMask(upperleft1, lowerright1, sa, 
01170                                 upperleftm, mask,
01171                                 upperleft2, da,
01172                                 radius, 0.5);
01173 }
01174 
01175 template <class SrcIterator, class SrcAccessor,
01176           class MaskIterator, class MaskAccessor,
01177           class DestIterator, class DestAccessor>
01178 inline void 
01179 discMedianWithMask(triple<SrcIterator, SrcIterator, SrcAccessor> src,
01180                     pair<MaskIterator, MaskAccessor> mask,
01181                     pair<DestIterator, DestAccessor> dest,
01182                     int radius)
01183 {
01184     vigra_precondition(radius >= 0, "discMedianWithMask(): Radius must be >= 0.");
01185     
01186     discRankOrderFilterWithMask(src.first, src.second, src.third,
01187                         mask.first, mask.second,
01188                         dest.first, dest.second,
01189                         radius, 0.5);
01190 }
01191 
01192 //@}
01193 
01194 } // namespace vigra
01195 
01196 #endif // VIGRA_FLATMORPHOLOGY_HXX

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

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