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

vigra/flatmorphology.hxx VIGRA

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