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

vigra/seededregiongrowing.hxx

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

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

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