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

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

html generated using doxygen and Python
VIGRA 1.2.0 (7 Aug 2003)