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

vigra/seededregiongrowing.hxx VIGRA

00001 /************************************************************************/
00002 /*                                                                      */
00003 /*         Copyright 1998-2010 by Ullrich Koethe, Hans Meine            */
00004 /*                                                                      */
00005 /*    This file is part of the VIGRA computer vision library.           */
00006 /*    The VIGRA Website is                                              */
00007 /*        http://hci.iwr.uni-heidelberg.de/vigra/                       */
00008 /*    Please direct questions, bug reports, and contributions to        */
00009 /*        ullrich.koethe@iwr.uni-heidelberg.de    or                    */
00010 /*        vigra@informatik.uni-hamburg.de                               */
00011 /*                                                                      */
00012 /*    Permission is hereby granted, free of charge, to any person       */
00013 /*    obtaining a copy of this software and associated documentation    */
00014 /*    files (the "Software"), to deal in the Software without           */
00015 /*    restriction, including without limitation the rights to use,      */
00016 /*    copy, modify, merge, publish, distribute, sublicense, and/or      */
00017 /*    sell copies of the Software, and to permit persons to whom the    */
00018 /*    Software is furnished to do so, subject to the following          */
00019 /*    conditions:                                                       */
00020 /*                                                                      */
00021 /*    The above copyright notice and this permission notice shall be    */
00022 /*    included in all copies or substantial portions of the             */
00023 /*    Software.                                                         */
00024 /*                                                                      */
00025 /*    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND    */
00026 /*    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   */
00027 /*    OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          */
00028 /*    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT       */
00029 /*    HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,      */
00030 /*    WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING      */
00031 /*    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR     */
00032 /*    OTHER DEALINGS IN THE SOFTWARE.                                   */
00033 /*                                                                      */
00034 /************************************************************************/
00035 
00036 #ifndef VIGRA_SEEDEDREGIONGROWING_HXX
00037 #define VIGRA_SEEDEDREGIONGROWING_HXX
00038 
00039 #include <vector>
00040 #include <stack>
00041 #include <queue>
00042 #include "utilities.hxx"
00043 #include "stdimage.hxx"
00044 #include "stdimagefunctions.hxx"
00045 #include "pixelneighborhood.hxx"
00046 #include "bucket_queue.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 /** \addtogroup SeededRegionGrowing Region Segmentation Algorithms
00165     Region growing, watersheds, and voronoi tesselation
00166 */
00167 //@{
00168 
00169 /********************************************************/
00170 /*                                                      */
00171 /*                    seededRegionGrowing               */
00172 /*                                                      */
00173 /********************************************************/
00174 
00175 /** Choose between different types of Region Growing */
00176 enum SRGType { 
00177     CompleteGrow = 0, 
00178     KeepContours = 1, 
00179     StopAtThreshold = 2, 
00180     SRGWatershedLabel = -1 
00181 };
00182 
00183 /** \brief Region Segmentation by means of Seeded Region Growing.
00184 
00185     This algorithm implements seeded region growing as described in
00186 
00187     R. Adams, L. Bischof: "<em> Seeded Region Growing</em>", IEEE Trans. on Pattern
00188     Analysis and Maschine Intelligence, vol 16, no 6, 1994, and
00189 
00190     Ullrich K&ouml;the:
00191     <em><a href="http://hci.iwr.uni-heidelberg.de/people/ukoethe/papers/index.php#cite_primary_segmentation">Primary Image Segmentation</a></em>,
00192     in: G. Sagerer, S.
00193     Posch, F. Kummert (eds.): Mustererkennung 1995, Proc. 17. DAGM-Symposium,
00194     Springer 1995
00195 
00196     The seed image is a partly segmented image which contains uniquely
00197     labeled regions (the seeds) and unlabeled pixels (the candidates, label 0).
00198     The highest seed label found in the seed image is returned by the algorithm.
00199     
00200     Seed regions can be as large as you wish and as small as one pixel. If
00201     there are no candidates, the algorithm will simply copy the seed image
00202     into the output image. Otherwise it will aggregate the candidates into
00203     the existing regions so that a cost function is minimized. 
00204     Candidates are taken from the neighborhood of the already assigned pixels, 
00205     where the type of neighborhood is determined by parameter <tt>neighborhood</tt>
00206     which can take the values <tt>FourNeighborCode()</tt> (the default) 
00207     or <tt>EightNeighborCode()</tt>. The algorithm basically works as follows 
00208     (illustrated for 4-neighborhood, but 8-neighborhood works in the same way):
00209 
00210     <ol>
00211 
00212     <li> Find all candidate pixels that are 4-adjacent to a seed region.
00213     Calculate the cost for aggregating each candidate into its adjacent region
00214     and put the candidates into a priority queue.
00215 
00216     <li> While( priority queue is not empty and termination criterion is not fulfilled)
00217 
00218         <ol>
00219 
00220         <li> Take the candidate with least cost from the queue. If it has not
00221         already been merged, merge it with it's adjacent region.
00222 
00223         <li> Put all candidates that are 4-adjacent to the pixel just processed
00224         into the priority queue.
00225 
00226         </ol>
00227 
00228     </ol>
00229 
00230     <tt>SRGType</tt> can take the following values:
00231     
00232     <DL>
00233     <DT><tt>CompleteGrow</tt> <DD> produce a complete tesselation of the volume (default).
00234     <DT><tt>KeepContours</tt> <DD> keep a 1-voxel wide unlabeled contour between all regions.
00235     <DT><tt>StopAtThreshold</tt> <DD> stop when the boundary indicator values exceed the 
00236                              threshold given by parameter <tt>max_cost</tt>.
00237     <DT><tt>KeepContours | StopAtThreshold</tt> <DD> keep 1-voxel wide contour and stop at given <tt>max_cost</tt>.
00238     </DL>
00239 
00240     The cost is determined jointly by the source image and the
00241     region statistics functor. The source image contains feature values for each
00242     pixel which will be used by the region statistics functor to calculate and
00243     update statistics for each region and to calculate the cost for each
00244     candidate. The <TT>RegionStatisticsArray</TT> must be compatible to the
00245     \ref ArrayOfRegionStatistics functor and contains an <em> array</em> of
00246     statistics objects for each region. The indices must correspond to the
00247     labels of the seed regions. The statistics for the initial regions must have
00248     been calculated prior to calling <TT>seededRegionGrowing()</TT> (for example by
00249     means of \ref inspectTwoImagesIf()).
00250 
00251     For each candidate
00252     <TT>x</TT> that is adjacent to region <TT>i</TT>, the algorithm will call
00253     <TT>stats[i].cost(as(x))</TT> to get the cost (where <TT>x</TT> is a <TT>SrcIterator</TT>
00254     and <TT>as</TT> is
00255     the SrcAccessor). When a candidate has been merged with a region, the
00256     statistics are updated by calling <TT>stats[i].operator()(as(x))</TT>. Since
00257     the <TT>RegionStatisticsArray</TT> is passed by reference, this will overwrite
00258     the original statistics.
00259 
00260     If a candidate could be merged into more than one regions with identical
00261     cost, the algorithm will favour the nearest region. If <tt>StopAtThreshold</tt> is active, 
00262     and the cost of the current candidate at any point in the algorithm exceeds the optional 
00263     <tt>max_cost</tt> value (which defaults to <tt>NumericTraits<double>::max()</tt>), 
00264     region growing is aborted, and all voxels not yet assigned to a region remain unlabeled.
00265 
00266     In some cases, the cost only depends on the feature value of the current
00267     pixel. Then the update operation will simply be a no-op, and the <TT>cost()</TT>
00268     function returns its argument. This behavior is implemented by the
00269     \ref SeedRgDirectValueFunctor. With <tt>SRGType == KeepContours</tt>,
00270     this is equivalent to the watershed algorithm.
00271 
00272     <b> Declarations:</b>
00273 
00274     pass arguments explicitly:
00275     \code
00276     namespace vigra {
00277         template <class SrcIterator, class SrcAccessor,
00278                   class SeedImageIterator, class SeedAccessor,
00279                   class DestIterator, class DestAccessor,
00280                   class RegionStatisticsArray, class Neighborhood>
00281         typename SeedAccessor::value_type 
00282         seededRegionGrowing(SrcIterator srcul, SrcIterator srclr, SrcAccessor as,
00283                             SeedImageIterator seedsul, SeedAccessor aseeds,
00284                             DestIterator destul, DestAccessor ad,
00285                             RegionStatisticsArray & stats,
00286                             SRGType srgType = CompleteGrow,
00287                             Neighborhood neighborhood = FourNeighborCode(),
00288                             double max_cost = NumericTraits<double>::max());
00289     }
00290     \endcode
00291 
00292     use argument objects in conjunction with \ref ArgumentObjectFactories :
00293     \code
00294     namespace vigra {
00295         template <class SrcIterator, class SrcAccessor,
00296                   class SeedImageIterator, class SeedAccessor,
00297                   class DestIterator, class DestAccessor,
00298                   class RegionStatisticsArray, class Neighborhood>
00299         typename SeedAccessor::value_type
00300         seededRegionGrowing(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00301                             pair<SeedImageIterator, SeedAccessor> seeds,
00302                             pair<DestIterator, DestAccessor> dest,
00303                             RegionStatisticsArray & stats,
00304                             SRGType srgType = CompleteGrow,
00305                             Neighborhood neighborhood = FourNeighborCode(),
00306                             double max_cost = NumericTraits<double>::max());
00307     }
00308     \endcode
00309 
00310     <b> Usage:</b>
00311 
00312     <b>\#include</b> <vigra/seededregiongrowing.hxx><br>
00313     Namespace: vigra
00314 
00315     Example: implementation of the voronoi tesselation
00316 
00317     \code
00318     vigra::BImage points(w,h);
00319     vigra::FImage dist(x,y);
00320 
00321     // empty edge image
00322     points = 0;
00323     dist = 0;
00324 
00325     int max_region_label = 100;
00326 
00327     // throw in some random points:
00328     for(int i = 1; i <= max_region_label; ++i)
00329            points(w * rand() / RAND_MAX , h * rand() / RAND_MAX) = i;
00330 
00331     // calculate Euclidean distance transform
00332     vigra::distanceTransform(srcImageRange(points), destImage(dist), 2);
00333 
00334     // init statistics functor
00335     vigra::ArrayOfRegionStatistics<vigra::SeedRgDirectValueFunctor<float> >
00336                                               stats(max_region_label);
00337 
00338     // find voronoi region of each point
00339     vigra:: seededRegionGrowing(srcImageRange(dist), srcImage(points),
00340                                destImage(points), stats);
00341     \endcode
00342 
00343     <b> Required Interface:</b>
00344 
00345     \code
00346     SrcIterator src_upperleft, src_lowerright;
00347     SeedImageIterator seed_upperleft;
00348     DestIterator dest_upperleft;
00349 
00350     SrcAccessor src_accessor;
00351     SeedAccessor seed_accessor;
00352     DestAccessor dest_accessor;
00353 
00354     RegionStatisticsArray stats;
00355 
00356     // calculate costs
00357     RegionStatisticsArray::value_type::cost_type cost =
00358         stats[seed_accessor(seed_upperleft)].cost(src_accessor(src_upperleft));
00359 
00360     // compare costs
00361     cost < cost;
00362 
00363     // update statistics
00364     stats[seed_accessor(seed_upperleft)](src_accessor(src_upperleft));
00365 
00366     // set result
00367     dest_accessor.set(seed_accessor(seed_upperleft), dest_upperleft);
00368     \endcode
00369 
00370     Further requirements are determined by the <TT>RegionStatisticsArray</TT>.
00371 */
00372 doxygen_overloaded_function(template <...> void seededRegionGrowing)
00373 
00374 template <class SrcIterator, class SrcAccessor,
00375           class SeedImageIterator, class SeedAccessor,
00376           class DestIterator, class DestAccessor,
00377           class RegionStatisticsArray, class Neighborhood>
00378 typename SeedAccessor::value_type
00379 seededRegionGrowing(SrcIterator srcul,
00380                     SrcIterator srclr, SrcAccessor as,
00381                     SeedImageIterator seedsul, SeedAccessor aseeds,
00382                     DestIterator destul, DestAccessor ad,
00383                     RegionStatisticsArray & stats,
00384                     SRGType srgType,
00385                     Neighborhood,
00386                     double max_cost)
00387 {
00388     int w = srclr.x - srcul.x;
00389     int h = srclr.y - srcul.y;
00390     int count = 0;
00391 
00392     SrcIterator isy = srcul, isx = srcul;  // iterators for the src image
00393 
00394     typedef typename SeedAccessor::value_type LabelType;
00395     typedef typename RegionStatisticsArray::value_type RegionStatistics;
00396     typedef typename RegionStatistics::cost_type CostType;
00397     typedef detail::SeedRgPixel<CostType> Pixel;
00398 
00399     typename Pixel::Allocator allocator;
00400 
00401     typedef std::priority_queue<Pixel *, std::vector<Pixel *>,
00402                                 typename Pixel::Compare>  SeedRgPixelHeap;
00403 
00404     // copy seed image in an image with border
00405     IImage regions(w+2, h+2);
00406     IImage::Iterator ir = regions.upperLeft() + Diff2D(1,1);
00407     IImage::Iterator iry, irx;
00408 
00409     initImageBorder(destImageRange(regions), 1, SRGWatershedLabel);
00410     copyImage(seedsul, seedsul+Diff2D(w,h), aseeds, ir, regions.accessor());
00411 
00412     // allocate and init memory for the results
00413 
00414     SeedRgPixelHeap pheap;
00415     int cneighbor, maxRegionLabel = 0;
00416     
00417     typedef typename Neighborhood::Direction Direction;
00418     int directionCount = Neighborhood::DirectionCount;
00419     
00420     Point2D pos(0,0);
00421     for(isy=srcul, iry=ir, pos.y=0; pos.y<h;
00422         ++pos.y, ++isy.y, ++iry.y)
00423     {
00424         for(isx=isy, irx=iry, pos.x=0; pos.x<w;
00425             ++pos.x, ++isx.x, ++irx.x)
00426         {
00427             if(*irx == 0)
00428             {
00429                 // find candidate pixels for growing and fill heap
00430                 for(int i=0; i<directionCount; i++)
00431                 {
00432                     // cneighbor = irx[dist[i]];
00433                     cneighbor = irx[Neighborhood::diff((Direction)i)];
00434                     if(cneighbor > 0)
00435                     {
00436                         CostType cost = stats[cneighbor].cost(as(isx));
00437 
00438                         Pixel * pixel =
00439                             allocator.create(pos, pos+Neighborhood::diff((Direction)i), cost, count++, cneighbor);
00440                         pheap.push(pixel);
00441                     }
00442                 }
00443             }
00444             else
00445             {
00446                 vigra_precondition((LabelType)*irx <= stats.maxRegionLabel(),
00447                     "seededRegionGrowing(): Largest label exceeds size of RegionStatisticsArray.");
00448                 if(maxRegionLabel < *irx)
00449                     maxRegionLabel = *irx;
00450             }
00451         }
00452     }
00453     
00454     // perform region growing
00455     while(pheap.size() != 0)
00456     {
00457         Pixel * pixel = pheap.top();
00458         pheap.pop();
00459 
00460         Point2D pos = pixel->location_;
00461         Point2D nearest = pixel->nearest_;
00462         int lab = pixel->label_;
00463         CostType cost = pixel->cost_;
00464 
00465         allocator.dismiss(pixel);
00466 
00467         if((srgType & StopAtThreshold) != 0 && cost > max_cost)
00468             break;
00469 
00470         irx = ir + pos;
00471         isx = srcul + pos;
00472 
00473         if(*irx) // already labelled region / watershed?
00474             continue;
00475 
00476         if((srgType & KeepContours) != 0)
00477         {
00478             for(int i=0; i<directionCount; i++)
00479             {
00480                 cneighbor = irx[Neighborhood::diff((Direction)i)];
00481                 if((cneighbor>0) && (cneighbor != lab))
00482                 {
00483                     lab = SRGWatershedLabel;
00484                     break;
00485                 }
00486             }
00487         }
00488 
00489         *irx = lab;
00490 
00491         if((srgType & KeepContours) == 0 || lab > 0)
00492         {
00493             // update statistics
00494             stats[*irx](as(isx));
00495 
00496             // search neighborhood
00497             // second pass: find new candidate pixels
00498             for(int i=0; i<directionCount; i++)
00499             {
00500                 if(irx[Neighborhood::diff((Direction)i)] == 0)
00501                 {
00502                     CostType cost = stats[lab].cost(as(isx, Neighborhood::diff((Direction)i)));
00503 
00504                     Pixel * new_pixel =
00505                         allocator.create(pos+Neighborhood::diff((Direction)i), nearest, cost, count++, lab);
00506                     pheap.push(new_pixel);
00507                 }
00508             }
00509         }
00510     }
00511     
00512     // free temporary memory
00513     while(pheap.size() != 0)
00514     {
00515         allocator.dismiss(pheap.top());
00516         pheap.pop();
00517     }
00518 
00519     // write result
00520     transformImage(ir, ir+Point2D(w,h), regions.accessor(), destul, ad,
00521                    detail::UnlabelWatersheds());
00522 
00523     return (LabelType)maxRegionLabel;
00524 }
00525 
00526 template <class SrcIterator, class SrcAccessor,
00527           class SeedImageIterator, class SeedAccessor,
00528           class DestIterator, class DestAccessor,
00529           class RegionStatisticsArray, class Neighborhood>
00530 inline typename SeedAccessor::value_type
00531 seededRegionGrowing(SrcIterator srcul,
00532                     SrcIterator srclr, SrcAccessor as,
00533                     SeedImageIterator seedsul, SeedAccessor aseeds,
00534                     DestIterator destul, DestAccessor ad,
00535                     RegionStatisticsArray & stats,
00536                     SRGType srgType,
00537                     Neighborhood n)
00538 {
00539     return seededRegionGrowing(srcul, srclr, as,
00540                                 seedsul, aseeds,
00541                                 destul, ad,
00542                                 stats, srgType, n, NumericTraits<double>::max());
00543 }
00544 
00545 
00546 
00547 template <class SrcIterator, class SrcAccessor,
00548           class SeedImageIterator, class SeedAccessor,
00549           class DestIterator, class DestAccessor,
00550           class RegionStatisticsArray>
00551 inline typename SeedAccessor::value_type
00552 seededRegionGrowing(SrcIterator srcul,
00553                     SrcIterator srclr, SrcAccessor as,
00554                     SeedImageIterator seedsul, SeedAccessor aseeds,
00555                     DestIterator destul, DestAccessor ad,
00556                     RegionStatisticsArray & stats,
00557                     SRGType srgType)
00558 {
00559     return seededRegionGrowing(srcul, srclr, as,
00560                                 seedsul, aseeds,
00561                                 destul, ad,
00562                                 stats, srgType, FourNeighborCode());
00563 }
00564 
00565 template <class SrcIterator, class SrcAccessor,
00566           class SeedImageIterator, class SeedAccessor,
00567           class DestIterator, class DestAccessor,
00568           class RegionStatisticsArray>
00569 inline typename SeedAccessor::value_type
00570 seededRegionGrowing(SrcIterator srcul,
00571                     SrcIterator srclr, SrcAccessor as,
00572                     SeedImageIterator seedsul, SeedAccessor aseeds,
00573                     DestIterator destul, DestAccessor ad,
00574                     RegionStatisticsArray & stats)
00575 {
00576     return seededRegionGrowing(srcul, srclr, as,
00577                                 seedsul, aseeds,
00578                                 destul, ad,
00579                                 stats, CompleteGrow);
00580 }
00581 
00582 template <class SrcIterator, class SrcAccessor,
00583           class SeedImageIterator, class SeedAccessor,
00584           class DestIterator, class DestAccessor,
00585           class RegionStatisticsArray, class Neighborhood>
00586 inline typename SeedAccessor::value_type
00587 seededRegionGrowing(triple<SrcIterator, SrcIterator, SrcAccessor> img1,
00588                     pair<SeedImageIterator, SeedAccessor> img3,
00589                     pair<DestIterator, DestAccessor> img4,
00590                     RegionStatisticsArray & stats,
00591                     SRGType srgType, 
00592                     Neighborhood n,
00593                     double max_cost)
00594 {
00595     return seededRegionGrowing(img1.first, img1.second, img1.third,
00596                                 img3.first, img3.second,
00597                                 img4.first, img4.second,
00598                                 stats, srgType, n, max_cost);
00599 }
00600 
00601 template <class SrcIterator, class SrcAccessor,
00602           class SeedImageIterator, class SeedAccessor,
00603           class DestIterator, class DestAccessor,
00604           class RegionStatisticsArray, class Neighborhood>
00605 inline typename SeedAccessor::value_type
00606 seededRegionGrowing(triple<SrcIterator, SrcIterator, SrcAccessor> img1,
00607                     pair<SeedImageIterator, SeedAccessor> img3,
00608                     pair<DestIterator, DestAccessor> img4,
00609                     RegionStatisticsArray & stats,
00610                     SRGType srgType, 
00611                     Neighborhood n)
00612 {
00613     return seededRegionGrowing(img1.first, img1.second, img1.third,
00614                                 img3.first, img3.second,
00615                                 img4.first, img4.second,
00616                                 stats, srgType, n, NumericTraits<double>::max());
00617 }
00618 
00619 template <class SrcIterator, class SrcAccessor,
00620           class SeedImageIterator, class SeedAccessor,
00621           class DestIterator, class DestAccessor,
00622           class RegionStatisticsArray>
00623 inline typename SeedAccessor::value_type
00624 seededRegionGrowing(triple<SrcIterator, SrcIterator, SrcAccessor> img1,
00625                     pair<SeedImageIterator, SeedAccessor> img3,
00626                     pair<DestIterator, DestAccessor> img4,
00627                     RegionStatisticsArray & stats,
00628                     SRGType srgType)
00629 {
00630     return seededRegionGrowing(img1.first, img1.second, img1.third,
00631                                 img3.first, img3.second,
00632                                 img4.first, img4.second,
00633                                 stats, srgType, FourNeighborCode());
00634 }
00635 
00636 template <class SrcIterator, class SrcAccessor,
00637           class SeedImageIterator, class SeedAccessor,
00638           class DestIterator, class DestAccessor,
00639           class RegionStatisticsArray>
00640 inline typename SeedAccessor::value_type
00641 seededRegionGrowing(triple<SrcIterator, SrcIterator, SrcAccessor> img1,
00642                     pair<SeedImageIterator, SeedAccessor> img3,
00643                     pair<DestIterator, DestAccessor> img4,
00644                     RegionStatisticsArray & stats)
00645 {
00646     return seededRegionGrowing(img1.first, img1.second, img1.third,
00647                             img3.first, img3.second,
00648                             img4.first, img4.second,
00649                             stats, CompleteGrow);
00650 }
00651 
00652 template <class SrcIterator, class SrcAccessor,
00653           class DestIterator, class DestAccessor,
00654           class RegionStatisticsArray, class Neighborhood>
00655 typename DestAccessor::value_type 
00656 fastSeededRegionGrowing(SrcIterator srcul, SrcIterator srclr, SrcAccessor as,
00657                         DestIterator destul, DestAccessor ad,
00658                         RegionStatisticsArray & stats,
00659                         SRGType srgType,
00660                         Neighborhood,
00661                         double max_cost,
00662                         std::ptrdiff_t bucket_count = 256)
00663 {
00664     typedef typename DestAccessor::value_type LabelType;
00665 
00666     vigra_precondition((srgType & KeepContours) == 0,
00667        "fastSeededRegionGrowing(): the turbo algorithm doesn't support 'KeepContours', sorry.");
00668     
00669     int w = srclr.x - srcul.x;
00670     int h = srclr.y - srcul.y;
00671 
00672     SrcIterator isy = srcul, isx = srcul;  // iterators for the src image
00673     DestIterator idy = destul, idx = destul;  // iterators for the dest image
00674 
00675     BucketQueue<Point2D, true> pqueue(bucket_count);
00676     LabelType maxRegionLabel = 0;
00677     
00678     Point2D pos(0,0);
00679     for(isy=srcul, idy = destul, pos.y=0; pos.y<h; ++pos.y, ++isy.y, ++idy.y)
00680     {
00681         for(isx=isy, idx=idy, pos.x=0; pos.x<w; ++pos.x, ++isx.x, ++idx.x)
00682         {
00683             LabelType label = ad(idx);
00684             if(label != 0)
00685             {
00686                 vigra_precondition(label <= stats.maxRegionLabel(),
00687                     "fastSeededRegionGrowing(): Largest label exceeds size of RegionStatisticsArray.");
00688 
00689                 if(maxRegionLabel < label)
00690                     maxRegionLabel = label;
00691                 
00692                 AtImageBorder atBorder = isAtImageBorder(pos.x, pos.y, w, h);
00693                 if(atBorder == NotAtBorder)
00694                 {
00695                     NeighborhoodCirculator<DestIterator, Neighborhood> c(idx), cend(c);
00696                     do
00697                     {
00698                         if(ad(c) == 0)
00699                         {
00700                             std::ptrdiff_t priority = (std::ptrdiff_t)stats[label].cost(as(isx));
00701                             pqueue.push(pos, priority);
00702                             break;
00703                         }
00704                     }
00705                     while(++c != cend);
00706                 }
00707                 else
00708                 {
00709                     RestrictedNeighborhoodCirculator<DestIterator, Neighborhood> 
00710                                                             c(idx, atBorder), cend(c);
00711                     do
00712                     {
00713                         if(ad(c) == 0)
00714                         {
00715                             std::ptrdiff_t priority = (std::ptrdiff_t)stats[label].cost(as(isx));
00716                             pqueue.push(pos, priority);
00717                             break;
00718                         }
00719                     }
00720                     while(++c != cend);
00721                 }
00722             }
00723         }
00724     }
00725     
00726     // perform region growing
00727     while(!pqueue.empty())
00728     {
00729         Point2D pos = pqueue.top();
00730         std::ptrdiff_t cost = pqueue.topPriority();
00731         pqueue.pop();
00732         
00733         if((srgType & StopAtThreshold) != 0 && cost > max_cost)
00734             break;
00735 
00736         idx = destul + pos;
00737         isx = srcul + pos;
00738         
00739         std::ptrdiff_t label = ad(idx);
00740 
00741         AtImageBorder atBorder = isAtImageBorder(pos.x, pos.y, w, h);
00742         if(atBorder == NotAtBorder)
00743         {
00744             NeighborhoodCirculator<DestIterator, Neighborhood> c(idx), cend(c);
00745             
00746             do
00747             {
00748                 std::ptrdiff_t nlabel = ad(c);
00749                 if(nlabel == 0)
00750                 {
00751                     ad.set(label, idx, c.diff());
00752                     std::ptrdiff_t priority = 
00753                            std::max((std::ptrdiff_t)stats[label].cost(as(isx, c.diff())), cost);
00754                     pqueue.push(pos+c.diff(), priority);
00755                 }
00756             }
00757             while(++c != cend);
00758         }
00759         else
00760         {
00761             RestrictedNeighborhoodCirculator<DestIterator, Neighborhood> 
00762                                                     c(idx, atBorder), cend(c);
00763             do
00764             {
00765                 std::ptrdiff_t nlabel = ad(c);
00766                 if(nlabel == 0)
00767                 {
00768                     ad.set(label, idx, c.diff());
00769                     std::ptrdiff_t priority = 
00770                            std::max((std::ptrdiff_t)stats[label].cost(as(isx, c.diff())), cost);
00771                     pqueue.push(pos+c.diff(), priority);
00772                 }
00773             }
00774             while(++c != cend);
00775         }
00776     }
00777     
00778     return maxRegionLabel;
00779 }
00780 
00781 template <class SrcIterator, class SrcAccessor,
00782           class DestIterator, class DestAccessor,
00783           class RegionStatisticsArray, class Neighborhood>
00784 inline typename DestAccessor::value_type 
00785 fastSeededRegionGrowing(SrcIterator srcul, SrcIterator srclr, SrcAccessor as,
00786                         DestIterator destul, DestAccessor ad,
00787                         RegionStatisticsArray & stats,
00788                         SRGType srgType,
00789                         Neighborhood n)
00790 {
00791     return fastSeededRegionGrowing(srcul, srclr, as,
00792                                     destul, ad,
00793                                     stats, srgType, n, NumericTraits<double>::max(), 256);
00794 }
00795 
00796 template <class SrcIterator, class SrcAccessor,
00797           class DestIterator, class DestAccessor,
00798           class RegionStatisticsArray>
00799 inline typename DestAccessor::value_type 
00800 fastSeededRegionGrowing(SrcIterator srcul, SrcIterator srclr, SrcAccessor as,
00801                         DestIterator destul, DestAccessor ad,
00802                         RegionStatisticsArray & stats,
00803                         SRGType srgType)
00804 {
00805     return fastSeededRegionGrowing(srcul, srclr, as,
00806                                     destul, ad,
00807                                     stats, srgType, FourNeighborCode());
00808 }
00809 
00810 template <class SrcIterator, class SrcAccessor,
00811           class DestIterator, class DestAccessor,
00812           class RegionStatisticsArray>
00813 inline typename DestAccessor::value_type 
00814 fastSeededRegionGrowing(SrcIterator srcul, SrcIterator srclr, SrcAccessor as,
00815                         DestIterator destul, DestAccessor ad,
00816                         RegionStatisticsArray & stats)
00817 {
00818     return fastSeededRegionGrowing(srcul, srclr, as,
00819                                     destul, ad,
00820                                     stats, CompleteGrow);
00821 }
00822 
00823 template <class SrcIterator, class SrcAccessor,
00824           class DestIterator, class DestAccessor,
00825           class RegionStatisticsArray, class Neighborhood>
00826 inline typename DestAccessor::value_type 
00827 fastSeededRegionGrowing(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00828                         pair<DestIterator, DestAccessor> dest,
00829                         RegionStatisticsArray & stats,
00830                         SRGType srgType, 
00831                         Neighborhood n,
00832                         double max_cost,
00833                         std::ptrdiff_t bucket_count = 256)
00834 {
00835     return fastSeededRegionGrowing(src.first, src.second, src.third,
00836                                     dest.first, dest.second,
00837                                     stats, srgType, n, max_cost, bucket_count);
00838 }
00839 
00840 /********************************************************/
00841 /*                                                      */
00842 /*               SeedRgDirectValueFunctor               */
00843 /*                                                      */
00844 /********************************************************/
00845 
00846 /** \brief Statistics functor to be used for seeded region growing.
00847 
00848     This functor can be used if the cost of a candidate during
00849     \ref seededRegionGrowing() is equal to the feature value of that
00850     candidate and does not depend on properties of the region it is going to
00851     be merged with.
00852 
00853     <b>\#include</b> <vigra/seededregiongrowing.hxx><br>
00854     Namespace: vigra
00855 
00856 
00857      <b> Required Interface:</b>
00858 
00859      no requirements
00860 */
00861 template <class Value>
00862 class SeedRgDirectValueFunctor
00863 {
00864   public:
00865         /** the functor's argument type
00866         */
00867     typedef Value argument_type;
00868 
00869         /** the functor's result type (unused, only necessary for
00870             use of SeedRgDirectValueFunctor in \ref vigra::ArrayOfRegionStatistics
00871         */
00872     typedef Value result_type;
00873 
00874         /** \deprecated use argument_type
00875         */
00876     typedef Value value_type;
00877 
00878         /** the return type of the cost() function
00879         */
00880     typedef Value cost_type;
00881 
00882         /** Do nothing (since we need not update region statistics).
00883         */
00884     void operator()(argument_type const &) const {}
00885 
00886         /** Return argument (since cost is identical to feature value)
00887         */
00888     cost_type const & cost(argument_type const & v) const
00889     {
00890         return v;
00891     }
00892 };
00893 
00894 //@}
00895 
00896 } // namespace vigra
00897 
00898 #endif // VIGRA_SEEDEDREGIONGROWING_HXX
00899 

© 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.8.0 (20 Sep 2011)