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

details vigra/seededregiongrowing.hxx VIGRA

00001 /************************************************************************/
00002 /*                                                                      */
00003 /*         Copyright 1998-2003 by Ullrich Koethe, Hans Meine            */
00004 /*       Cognitive Systems Group, University of Hamburg, Germany        */
00005 /*                                                                      */
00006 /*    This file is part of the VIGRA computer vision library.           */
00007 /*    ( Version 1.5.0, Dec 07 2006 )                                    */
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 #ifndef VIGRA_SEEDEDREGIONGROWING_HXX
00039 #define VIGRA_SEEDEDREGIONGROWING_HXX
00040 
00041 #include <vector>
00042 #include <stack>
00043 #include <queue>
00044 #include "utilities.hxx"
00045 #include "stdimage.hxx"
00046 #include "stdimagefunctions.hxx"
00047 
00048 namespace vigra {
00049 
00050 namespace detail {
00051 
00052 template <class COST>
00053 class SeedRgPixel
00054 {
00055 public:
00056     Point2D location_, nearest_;
00057     COST cost_;
00058     int count_;
00059     int label_;
00060     int dist_;
00061 
00062     SeedRgPixel()
00063     : location_(0,0), nearest_(0,0), cost_(0), count_(0), label_(0)
00064     {}
00065 
00066     SeedRgPixel(Point2D const & location, Point2D const & nearest,
00067                 COST const & cost, int const & count, int const & label)
00068     : location_(location), nearest_(nearest),
00069       cost_(cost), count_(count), label_(label)
00070     {
00071         int dx = location_.x - nearest_.x;
00072         int dy = location_.y - nearest_.y;
00073         dist_ = dx * dx + dy * dy;
00074     }
00075 
00076     void set(Point2D const & location, Point2D const & nearest,
00077              COST const & cost, int const & count, int const & label)
00078     {
00079         location_ = location;
00080         nearest_ = nearest;
00081         cost_ = cost;
00082         count_ = count;
00083         label_ = label;
00084 
00085         int dx = location_.x - nearest_.x;
00086         int dy = location_.y - nearest_.y;
00087         dist_ = dx * dx + dy * dy;
00088     }
00089 
00090     struct Compare
00091     {
00092         // must implement > since priority_queue looks for largest element
00093         bool operator()(SeedRgPixel const & l,
00094                         SeedRgPixel const & r) const
00095         {
00096             if(r.cost_ == l.cost_)
00097             {
00098                 if(r.dist_ == l.dist_) return r.count_ < l.count_;
00099 
00100                 return r.dist_ < l.dist_;
00101             }
00102 
00103             return r.cost_ < l.cost_;
00104         }
00105         bool operator()(SeedRgPixel const * l,
00106                         SeedRgPixel const * r) const
00107         {
00108             if(r->cost_ == l->cost_)
00109             {
00110                 if(r->dist_ == l->dist_) return r->count_ < l->count_;
00111 
00112                 return r->dist_ < l->dist_;
00113             }
00114 
00115             return r->cost_ < l->cost_;
00116         }
00117     };
00118 
00119     struct Allocator
00120     {
00121         ~Allocator()
00122         {
00123             while(!freelist_.empty())
00124             {
00125                 delete freelist_.top();
00126                 freelist_.pop();
00127             }
00128         }
00129 
00130         SeedRgPixel *
00131         create(Point2D const & location, Point2D const & nearest,
00132                COST const & cost, int const & count, int const & label)
00133         {
00134             if(!freelist_.empty())
00135             {
00136                 SeedRgPixel * res = freelist_.top();
00137                 freelist_.pop();
00138                 res->set(location, nearest, cost, count, label);
00139                 return res;
00140             }
00141 
00142             return new SeedRgPixel(location, nearest, cost, count, label);
00143         }
00144 
00145         void dismiss(SeedRgPixel * p)
00146         {
00147             freelist_.push(p);
00148         }
00149 
00150         std::stack<SeedRgPixel<COST> *> freelist_;
00151     };
00152 };
00153 
00154 struct UnlabelWatersheds
00155 {
00156     int operator()(int label) const
00157     {
00158         return label < 0 ? 0 : label;
00159     }
00160 };
00161 
00162 } // namespace detail
00163 
00164 enum SRGType { KeepContours, CompleteGrow, SRGWatershedLabel = -1 };
00165 
00166 /** \addtogroup SeededRegionGrowing Region Segmentation Algorithms
00167     Region growing, watersheds, and voronoi tesselation
00168 */
00169 //@{
00170 
00171 /********************************************************/
00172 /*                                                      */
00173 /*                    seededRegionGrowing               */
00174 /*                                                      */
00175 /********************************************************/
00176 
00177 /** \brief Region Segmentation by means of Seeded Region Growing.
00178 
00179     This algorithm implements seeded region growing as described in
00180 
00181     R. Adams, L. Bischof: "<em> Seeded Region Growing</em>", IEEE Trans. on Pattern
00182     Analysis and Maschine Intelligence, vol 16, no 6, 1994, and
00183 
00184     Ullrich K&ouml;the:
00185     <em> "<a href="http://kogs-www.informatik.uni-hamburg.de/~koethe/papers/#primary">Primary Image Segmentation</a>"</em>,
00186     in: G. Sagerer, S.
00187     Posch, F. Kummert (eds.): Mustererkennung 1995, Proc. 17. DAGM-Symposium,
00188     Springer 1995
00189 
00190     The seed image is a partly segmented image which contains uniquely
00191     labeled regions (the seeds) and unlabeled pixels (the candidates, label 0).
00192     Seed regions can be as large as you wish and as small as one pixel. If
00193     there are no candidates, the algorithm will simply copy the seed image
00194     into the output image. Otherwise it will aggregate the candidates into
00195     the existing regions so that a cost function is minimized. This
00196     works as follows:
00197 
00198     <ol>
00199 
00200     <li> Find all candidate pixels that are 4-adjacent to a seed region.
00201     Calculate the cost for aggregating each candidate into its adajacent region
00202     and put the candidates into a priority queue.
00203 
00204     <li> While( priority queue is not empty)
00205 
00206         <ol>
00207 
00208         <li> Take the candidate with least cost from the queue. If it has not
00209         already been merged, merge it with it's adjacent region.
00210 
00211         <li> Put all candidates that are 4-adjacent to the pixel just processed
00212         into the priority queue.
00213 
00214         </ol>
00215 
00216     </ol>
00217 
00218     If <tt>SRGType == CompleteGrow</tt> (the default), this algorithm
00219     will produce a complete 4-connected tesselation of the image.
00220     If <tt>SRGType == KeepContours</tt>, a one-pixel-wide border will be left
00221     between the regions. The border pixels get label 0 (zero).
00222 
00223     The cost is determined jointly by the source image and the
00224     region statistics functor. The source image contains feature values for each
00225     pixel which will be used by the region statistics functor to calculate and
00226     update statistics for each region and to calculate the cost for each
00227     candidate. The <TT>RegionStatisticsArray</TT> must be compatible to the
00228     \ref ArrayOfRegionStatistics functor and contains an <em> array</em> of
00229     statistics objects for each region. The indices must correspond to the
00230     labels of the seed regions. The statistics for the initial regions must have
00231     been calculated prior to calling <TT>seededRegionGrowing()</TT> (for example by
00232     means of \ref inspectTwoImagesIf()).
00233 
00234     For each candidate
00235     <TT>x</TT> that is adjacent to region <TT>i</TT>, the algorithm will call
00236     <TT>stats[i].cost(as(x))</TT> to get the cost (where <TT>x</TT> is a <TT>SrcImageIterator</TT>
00237     and <TT>as</TT> is
00238     the SrcAccessor). When a candidate has been merged with a region, the
00239     statistics are updated by calling <TT>stats[i].operator()(as(x))</TT>. Since
00240     the <TT>RegionStatisticsArray</TT> is passed by reference, this will overwrite
00241     the original statistics.
00242 
00243     If a candidate could be merged into more than one regions with identical
00244     cost, the algorithm will favour the nearest region.
00245 
00246     In some cases, the cost only depends on the feature value of the current
00247     pixel. Then the update operation will simply be a no-op, and the <TT>cost()</TT>
00248     function returns its argument. This behavior is implemented by the
00249     \ref SeedRgDirectValueFunctor. With <tt>SRGType == KeepContours</tt>,
00250     this is equivalent to the watershed algorithm.
00251 
00252     <b> Declarations:</b>
00253 
00254     pass arguments explicitly:
00255     \code
00256     namespace vigra {
00257         template <class SrcImageIterator, class SrcAccessor,
00258                   class SeedImageIterator, class SeedAccessor,
00259                   class DestImageIterator, class DestAccessor,
00260                   class RegionStatisticsArray>
00261         void seededRegionGrowing(SrcImageIterator srcul,
00262                                  SrcImageIterator srclr, SrcAccessor as,
00263                                  SeedImageIterator seedsul, SeedAccessor aseeds,
00264                                  DestImageIterator destul, DestAccessor ad,
00265                                  RegionStatisticsArray & stats,
00266                                  SRGType srgType = CompleteGrow);
00267     }
00268     \endcode
00269 
00270     use argument objects in conjunction with \ref ArgumentObjectFactories:
00271     \code
00272     namespace vigra {
00273         template <class SrcImageIterator, class SrcAccessor,
00274                   class SeedImageIterator, class SeedAccessor,
00275                   class DestImageIterator, class DestAccessor,
00276                   class RegionStatisticsArray>
00277         inline void
00278         seededRegionGrowing(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> img1,
00279                             pair<SeedImageIterator, SeedAccessor> img3,
00280                             pair<DestImageIterator, DestAccessor> img4,
00281                             RegionStatisticsArray & stats,
00282                             SRGType srgType = CompleteGrow);
00283     }
00284     \endcode
00285 
00286     <b> Usage:</b>
00287 
00288     <b>\#include</b> "<a href="seededregiongrowing_8hxx-source.html">vigra/seededregiongrowing.hxx</a>"<br>
00289     Namespace: vigra
00290 
00291     Example: implementation of the voronoi tesselation
00292 
00293     \code
00294     vigra::BImage points(w,h);
00295     vigra::FImage dist(x,y);
00296 
00297     // empty edge image
00298     points = 0;
00299     dist = 0;
00300 
00301     int max_region_label = 100;
00302 
00303     // throw in some random points:
00304     for(int i = 1; i <= max_region_label; ++i)
00305            points(w * rand() / RAND_MAX , h * rand() / RAND_MAX) = i;
00306 
00307     // calculate Euclidean distance transform
00308     vigra::distanceTransform(srcImageRange(points), destImage(dist), 2);
00309 
00310     // init statistics functor
00311     vigra::ArrayOfRegionStatistics<vigra::SeedRgDirectValueFunctor<float> >
00312                                               stats(max_region_label);
00313 
00314     // find voronoi region of each point
00315    vigra:: seededRegionGrowing(srcImageRange(dist), srcImage(points),
00316                                destImage(points), stats);
00317     \endcode
00318 
00319     <b> Required Interface:</b>
00320 
00321     \code
00322     SrcImageIterator src_upperleft, src_lowerright;
00323     SeedImageIterator seed_upperleft;
00324     DestImageIterator dest_upperleft;
00325 
00326     SrcAccessor src_accessor;
00327     SeedAccessor seed_accessor;
00328     DestAccessor dest_accessor;
00329 
00330     RegionStatisticsArray stats;
00331 
00332     // calculate costs
00333     RegionStatisticsArray::value_type::cost_type cost =
00334         stats[seed_accessor(seed_upperleft)].cost(src_accessor(src_upperleft));
00335 
00336     // compare costs
00337     cost < cost;
00338 
00339     // update statistics
00340     stats[seed_accessor(seed_upperleft)](src_accessor(src_upperleft));
00341 
00342     // set result
00343     dest_accessor.set(seed_accessor(seed_upperleft), dest_upperleft);
00344     \endcode
00345 
00346     Further requirements are determined by the <TT>RegionStatisticsArray</TT>.
00347 */
00348 template <class SrcImageIterator, class SrcAccessor,
00349           class SeedImageIterator, class SeedAccessor,
00350           class DestImageIterator, class DestAccessor,
00351           class RegionStatisticsArray>
00352 void seededRegionGrowing(SrcImageIterator srcul,
00353                          SrcImageIterator srclr, SrcAccessor as,
00354                          SeedImageIterator seedsul, SeedAccessor aseeds,
00355                          DestImageIterator destul, DestAccessor ad,
00356                          RegionStatisticsArray & stats,
00357                          const SRGType srgType)
00358 {
00359     int w = srclr.x - srcul.x;
00360     int h = srclr.y - srcul.y;
00361     int count = 0;
00362 
00363     SrcImageIterator isy = srcul, isx = srcul;  // iterators for the src image
00364 
00365     typedef typename RegionStatisticsArray::value_type RegionStatistics;
00366     typedef typename RegionStatistics::cost_type CostType;
00367     typedef detail::SeedRgPixel<CostType> Pixel;
00368 
00369     typename Pixel::Allocator allocator;
00370 
00371     typedef std::priority_queue<Pixel *, std::vector<Pixel *>,
00372                                 typename Pixel::Compare>  SeedRgPixelHeap;
00373 
00374     // copy seed image in an image with border
00375     IImage regions(w+2, h+2);
00376     IImage::Iterator ir = regions.upperLeft() + Diff2D(1,1);
00377     IImage::Iterator iry, irx;
00378 
00379     initImageBorder(destImageRange(regions), 1, SRGWatershedLabel);
00380     copyImage(seedsul, seedsul+Diff2D(w,h), aseeds, ir, regions.accessor());
00381 
00382     // allocate and init memory for the results
00383 
00384     SeedRgPixelHeap pheap;
00385     int cneighbor;
00386 
00387     static const Diff2D dist[] = { Diff2D(-1,0), Diff2D(0,-1),
00388                                    Diff2D(1,0),  Diff2D(0,1) };
00389 
00390     Point2D pos(0,0);
00391     for(isy=srcul, iry=ir, pos.y=0; pos.y<h;
00392         ++pos.y, ++isy.y, ++iry.y)
00393     {
00394         for(isx=isy, irx=iry, pos.x=0; pos.x<w;
00395             ++pos.x, ++isx.x, ++irx.x)
00396         {
00397             if(*irx == 0)
00398             {
00399                 // find candidate pixels for growing and fill heap
00400                 for(int i=0; i<4; i++)
00401                 {
00402                     cneighbor = irx[dist[i]];
00403                     if(cneighbor > 0)
00404                     {
00405                         CostType cost = stats[cneighbor].cost(as(isx));
00406 
00407                         Pixel * pixel =
00408                             allocator.create(pos, pos+dist[i], cost, count++, cneighbor);
00409                         pheap.push(pixel);
00410                     }
00411                 }
00412             }
00413         }
00414     }
00415 
00416     // perform region growing
00417     while(pheap.size() != 0)
00418     {
00419         Pixel * pixel = pheap.top();
00420         pheap.pop();
00421 
00422         Point2D pos = pixel->location_;
00423         Point2D nearest = pixel->nearest_;
00424         int lab = pixel->label_;
00425 
00426         allocator.dismiss(pixel);
00427 
00428         irx = ir + pos;
00429         isx = srcul + pos;
00430 
00431         if(*irx) // already labelled region / watershed?
00432             continue;
00433 
00434         if(srgType == KeepContours)
00435         {
00436             for(int i=0; i<4; i++)
00437             {
00438                 cneighbor = irx[dist[i]];
00439                 if((cneighbor>0) && (cneighbor != lab))
00440                 {
00441                     lab = SRGWatershedLabel;
00442                     break;
00443                 }
00444             }
00445         }
00446 
00447         *irx = lab;
00448 
00449         if((srgType != KeepContours) || (lab > 0))
00450         {
00451             // update statistics
00452             stats[*irx](as(isx));
00453 
00454             // search neighborhood
00455             // second pass: find new candidate pixels
00456             for(int i=0; i<4; i++)
00457             {
00458                 if(irx[dist[i]] == 0)
00459                 {
00460                     CostType cost = stats[lab].cost(as(isx, dist[i]));
00461 
00462                     Pixel * new_pixel =
00463                         allocator.create(pos+dist[i], nearest, cost, count++, lab);
00464                     pheap.push(new_pixel);
00465                 }
00466             }
00467         }
00468     }
00469 
00470     // write result
00471     transformImage(ir, ir+Point2D(w,h), regions.accessor(), destul, ad,
00472                    detail::UnlabelWatersheds());
00473 }
00474 
00475 template <class SrcImageIterator, class SrcAccessor,
00476           class SeedImageIterator, class SeedAccessor,
00477           class DestImageIterator, class DestAccessor,
00478           class RegionStatisticsArray>
00479 inline void
00480 seededRegionGrowing(SrcImageIterator srcul,
00481                     SrcImageIterator srclr, SrcAccessor as,
00482                     SeedImageIterator seedsul, SeedAccessor aseeds,
00483                     DestImageIterator destul, DestAccessor ad,
00484                     RegionStatisticsArray & stats)
00485 {
00486     seededRegionGrowing(srcul, srclr, as,
00487                         seedsul, aseeds,
00488                         destul, ad,
00489                         stats, CompleteGrow);
00490 }
00491 
00492 template <class SrcImageIterator, class SrcAccessor,
00493           class SeedImageIterator, class SeedAccessor,
00494           class DestImageIterator, class DestAccessor,
00495           class RegionStatisticsArray>
00496 inline void
00497 seededRegionGrowing(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> img1,
00498                     pair<SeedImageIterator, SeedAccessor> img3,
00499                     pair<DestImageIterator, DestAccessor> img4,
00500                     RegionStatisticsArray & stats,
00501                     SRGType srgType)
00502 {
00503     seededRegionGrowing(img1.first, img1.second, img1.third,
00504                         img3.first, img3.second,
00505                         img4.first, img4.second,
00506                         stats, srgType);
00507 }
00508 
00509 template <class SrcImageIterator, class SrcAccessor,
00510           class SeedImageIterator, class SeedAccessor,
00511           class DestImageIterator, class DestAccessor,
00512           class RegionStatisticsArray>
00513 inline void
00514 seededRegionGrowing(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> img1,
00515                     pair<SeedImageIterator, SeedAccessor> img3,
00516                     pair<DestImageIterator, DestAccessor> img4,
00517                     RegionStatisticsArray & stats)
00518 {
00519     seededRegionGrowing(img1.first, img1.second, img1.third,
00520                         img3.first, img3.second,
00521                         img4.first, img4.second,
00522                         stats, CompleteGrow);
00523 }
00524 
00525 /********************************************************/
00526 /*                                                      */
00527 /*               SeedRgDirectValueFunctor               */
00528 /*                                                      */
00529 /********************************************************/
00530 
00531 /** \brief Statistics functor to be used for seeded region growing.
00532 
00533     This functor can be used if the cost of a candidate during
00534     \ref seededRegionGrowing() is equal to the feature value of that
00535     candidate and does not depend on properties of the region it is going to
00536     be merged with.
00537 
00538     <b>\#include</b> "<a href="seededregiongrowing_8hxx-source.html">vigra/seededregiongrowing.hxx</a>"<br>
00539     Namespace: vigra
00540 
00541 
00542      <b> Required Interface:</b>
00543 
00544      no requirements
00545 */
00546 template <class Value>
00547 class SeedRgDirectValueFunctor
00548 {
00549   public:
00550         /** the functor's argument type
00551         */
00552     typedef Value argument_type;
00553 
00554         /** the functor's result type (unused, only necessary for
00555             use of SeedRgDirectValueFunctor in \ref vigra::ArrayOfRegionStatistics
00556         */
00557     typedef Value result_type;
00558 
00559         /** \deprecated use argument_type
00560         */
00561     typedef Value value_type;
00562 
00563         /** the return type of the cost() function
00564         */
00565     typedef Value cost_type;
00566 
00567         /** Do nothing (since we need not update region statistics).
00568         */
00569     void operator()(argument_type const &) const {}
00570 
00571         /** Return argument (since cost is identical to feature value)
00572         */
00573     cost_type const & cost(argument_type const & v) const
00574     {
00575         return v;
00576     }
00577 };
00578 
00579 //@}
00580 
00581 } // namespace vigra
00582 
00583 #endif // VIGRA_SEEDEDREGIONGROWING_HXX
00584 

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

html generated using doxygen and Python
VIGRA 1.5.0 (7 Dec 2006)