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