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

details vigra/flatmorphology.hxx VIGRA

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

© Ullrich Köthe (koethe@informatik.uni-hamburg.de)
Cognitive Systems Group, University of Hamburg, Germany

html generated using doxygen and Python
VIGRA 1.4.0 (21 Dec 2005)