[ 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 /*    ( Version 1.6.0, Aug 13 2008 )                                    */
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 /*        ullrich.koethe@iwr.uni-heidelberg.de    or                    */
00012 /*        vigra@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         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 doxygen_overloaded_function(template <...> void seededRegionGrowing)
00349 
00350 template <class SrcImageIterator, class SrcAccessor,
00351           class SeedImageIterator, class SeedAccessor,
00352           class DestImageIterator, class DestAccessor,
00353           class RegionStatisticsArray>
00354 void seededRegionGrowing(SrcImageIterator srcul,
00355                          SrcImageIterator srclr, SrcAccessor as,
00356                          SeedImageIterator seedsul, SeedAccessor aseeds,
00357                          DestImageIterator destul, DestAccessor ad,
00358                          RegionStatisticsArray & stats,
00359                          const SRGType srgType)
00360 {
00361     int w = srclr.x - srcul.x;
00362     int h = srclr.y - srcul.y;
00363     int count = 0;
00364 
00365     SrcImageIterator isy = srcul, isx = srcul;  // iterators for the src image
00366 
00367     typedef typename RegionStatisticsArray::value_type RegionStatistics;
00368     typedef typename RegionStatistics::cost_type CostType;
00369     typedef detail::SeedRgPixel<CostType> Pixel;
00370 
00371     typename Pixel::Allocator allocator;
00372 
00373     typedef std::priority_queue<Pixel *, std::vector<Pixel *>,
00374                                 typename Pixel::Compare>  SeedRgPixelHeap;
00375 
00376     // copy seed image in an image with border
00377     IImage regions(w+2, h+2);
00378     IImage::Iterator ir = regions.upperLeft() + Diff2D(1,1);
00379     IImage::Iterator iry, irx;
00380 
00381     initImageBorder(destImageRange(regions), 1, SRGWatershedLabel);
00382     copyImage(seedsul, seedsul+Diff2D(w,h), aseeds, ir, regions.accessor());
00383 
00384     // allocate and init memory for the results
00385 
00386     SeedRgPixelHeap pheap;
00387     int cneighbor;
00388 
00389     static const Diff2D dist[] = { Diff2D(-1,0), Diff2D(0,-1),
00390                                    Diff2D(1,0),  Diff2D(0,1) };
00391 
00392     Point2D pos(0,0);
00393     for(isy=srcul, iry=ir, pos.y=0; pos.y<h;
00394         ++pos.y, ++isy.y, ++iry.y)
00395     {
00396         for(isx=isy, irx=iry, pos.x=0; pos.x<w;
00397             ++pos.x, ++isx.x, ++irx.x)
00398         {
00399             if(*irx == 0)
00400             {
00401                 // find candidate pixels for growing and fill heap
00402                 for(int i=0; i<4; i++)
00403                 {
00404                     cneighbor = irx[dist[i]];
00405                     if(cneighbor > 0)
00406                     {
00407                         CostType cost = stats[cneighbor].cost(as(isx));
00408 
00409                         Pixel * pixel =
00410                             allocator.create(pos, pos+dist[i], cost, count++, cneighbor);
00411                         pheap.push(pixel);
00412                     }
00413                 }
00414             }
00415         }
00416     }
00417 
00418     // perform region growing
00419     while(pheap.size() != 0)
00420     {
00421         Pixel * pixel = pheap.top();
00422         pheap.pop();
00423 
00424         Point2D pos = pixel->location_;
00425         Point2D nearest = pixel->nearest_;
00426         int lab = pixel->label_;
00427 
00428         allocator.dismiss(pixel);
00429 
00430         irx = ir + pos;
00431         isx = srcul + pos;
00432 
00433         if(*irx) // already labelled region / watershed?
00434             continue;
00435 
00436         if(srgType == KeepContours)
00437         {
00438             for(int i=0; i<4; i++)
00439             {
00440                 cneighbor = irx[dist[i]];
00441                 if((cneighbor>0) && (cneighbor != lab))
00442                 {
00443                     lab = SRGWatershedLabel;
00444                     break;
00445                 }
00446             }
00447         }
00448 
00449         *irx = lab;
00450 
00451         if((srgType != KeepContours) || (lab > 0))
00452         {
00453             // update statistics
00454             stats[*irx](as(isx));
00455 
00456             // search neighborhood
00457             // second pass: find new candidate pixels
00458             for(int i=0; i<4; i++)
00459             {
00460                 if(irx[dist[i]] == 0)
00461                 {
00462                     CostType cost = stats[lab].cost(as(isx, dist[i]));
00463 
00464                     Pixel * new_pixel =
00465                         allocator.create(pos+dist[i], nearest, cost, count++, lab);
00466                     pheap.push(new_pixel);
00467                 }
00468             }
00469         }
00470     }
00471 
00472     // write result
00473     transformImage(ir, ir+Point2D(w,h), regions.accessor(), destul, ad,
00474                    detail::UnlabelWatersheds());
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(SrcImageIterator srcul,
00483                     SrcImageIterator srclr, SrcAccessor as,
00484                     SeedImageIterator seedsul, SeedAccessor aseeds,
00485                     DestImageIterator destul, DestAccessor ad,
00486                     RegionStatisticsArray & stats)
00487 {
00488     seededRegionGrowing(srcul, srclr, as,
00489                         seedsul, aseeds,
00490                         destul, ad,
00491                         stats, CompleteGrow);
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                     SRGType srgType)
00504 {
00505     seededRegionGrowing(img1.first, img1.second, img1.third,
00506                         img3.first, img3.second,
00507                         img4.first, img4.second,
00508                         stats, srgType);
00509 }
00510 
00511 template <class SrcImageIterator, class SrcAccessor,
00512           class SeedImageIterator, class SeedAccessor,
00513           class DestImageIterator, class DestAccessor,
00514           class RegionStatisticsArray>
00515 inline void
00516 seededRegionGrowing(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> img1,
00517                     pair<SeedImageIterator, SeedAccessor> img3,
00518                     pair<DestImageIterator, DestAccessor> img4,
00519                     RegionStatisticsArray & stats)
00520 {
00521     seededRegionGrowing(img1.first, img1.second, img1.third,
00522                         img3.first, img3.second,
00523                         img4.first, img4.second,
00524                         stats, CompleteGrow);
00525 }
00526 
00527 /********************************************************/
00528 /*                                                      */
00529 /*               SeedRgDirectValueFunctor               */
00530 /*                                                      */
00531 /********************************************************/
00532 
00533 /** \brief Statistics functor to be used for seeded region growing.
00534 
00535     This functor can be used if the cost of a candidate during
00536     \ref seededRegionGrowing() is equal to the feature value of that
00537     candidate and does not depend on properties of the region it is going to
00538     be merged with.
00539 
00540     <b>\#include</b> <<a href="seededregiongrowing_8hxx-source.html">vigra/seededregiongrowing.hxx</a>><br>
00541     Namespace: vigra
00542 
00543 
00544      <b> Required Interface:</b>
00545 
00546      no requirements
00547 */
00548 template <class Value>
00549 class SeedRgDirectValueFunctor
00550 {
00551   public:
00552         /** the functor's argument type
00553         */
00554     typedef Value argument_type;
00555 
00556         /** the functor's result type (unused, only necessary for
00557             use of SeedRgDirectValueFunctor in \ref vigra::ArrayOfRegionStatistics
00558         */
00559     typedef Value result_type;
00560 
00561         /** \deprecated use argument_type
00562         */
00563     typedef Value value_type;
00564 
00565         /** the return type of the cost() function
00566         */
00567     typedef Value cost_type;
00568 
00569         /** Do nothing (since we need not update region statistics).
00570         */
00571     void operator()(argument_type const &) const {}
00572 
00573         /** Return argument (since cost is identical to feature value)
00574         */
00575     cost_type const & cost(argument_type const & v) const
00576     {
00577         return v;
00578     }
00579 };
00580 
00581 //@}
00582 
00583 } // namespace vigra
00584 
00585 #endif // VIGRA_SEEDEDREGIONGROWING_HXX
00586 

© 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 (13 Aug 2008)