[ VIGRA Homepage | Class Index | Function Index | File Index | 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 /*    ( Version 1.3.2, Jan 27 2005 )                                    */
    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 conjunction 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(srcul, srclr, as,
    00472                         seedsul, aseeds,
    00473                         destul, ad,
    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.3.2 (27 Jan 2005)