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

seededregiongrowing.hxx
1 /************************************************************************/
2 /* */
3 /* Copyright 1998-2010 by Ullrich Koethe, Hans Meine */
4 /* */
5 /* This file is part of the VIGRA computer vision library. */
6 /* The VIGRA Website is */
7 /* http://hci.iwr.uni-heidelberg.de/vigra/ */
8 /* Please direct questions, bug reports, and contributions to */
9 /* ullrich.koethe@iwr.uni-heidelberg.de or */
10 /* vigra@informatik.uni-hamburg.de */
11 /* */
12 /* Permission is hereby granted, free of charge, to any person */
13 /* obtaining a copy of this software and associated documentation */
14 /* files (the "Software"), to deal in the Software without */
15 /* restriction, including without limitation the rights to use, */
16 /* copy, modify, merge, publish, distribute, sublicense, and/or */
17 /* sell copies of the Software, and to permit persons to whom the */
18 /* Software is furnished to do so, subject to the following */
19 /* conditions: */
20 /* */
21 /* The above copyright notice and this permission notice shall be */
22 /* included in all copies or substantial portions of the */
23 /* Software. */
24 /* */
25 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */
26 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */
27 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */
28 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */
29 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */
30 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */
31 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */
32 /* OTHER DEALINGS IN THE SOFTWARE. */
33 /* */
34 /************************************************************************/
35 
36 #ifndef VIGRA_SEEDEDREGIONGROWING_HXX
37 #define VIGRA_SEEDEDREGIONGROWING_HXX
38 
39 #include <vector>
40 #include <stack>
41 #include <queue>
42 #include "utilities.hxx"
43 #include "stdimage.hxx"
44 #include "stdimagefunctions.hxx"
45 #include "pixelneighborhood.hxx"
46 #include "bucket_queue.hxx"
47 
48 namespace vigra {
49 
50 namespace detail {
51 
52 template <class COST>
53 class SeedRgPixel
54 {
55 public:
56  Point2D location_, nearest_;
57  COST cost_;
58  int count_;
59  int label_;
60  int dist_;
61 
62  SeedRgPixel()
63  : location_(0,0), nearest_(0,0), cost_(0), count_(0), label_(0)
64  {}
65 
66  SeedRgPixel(Point2D const & location, Point2D const & nearest,
67  COST const & cost, int const & count, int const & label)
68  : location_(location), nearest_(nearest),
69  cost_(cost), count_(count), label_(label)
70  {
71  int dx = location_.x - nearest_.x;
72  int dy = location_.y - nearest_.y;
73  dist_ = dx * dx + dy * dy;
74  }
75 
76  void set(Point2D const & location, Point2D const & nearest,
77  COST const & cost, int const & count, int const & label)
78  {
79  location_ = location;
80  nearest_ = nearest;
81  cost_ = cost;
82  count_ = count;
83  label_ = label;
84 
85  int dx = location_.x - nearest_.x;
86  int dy = location_.y - nearest_.y;
87  dist_ = dx * dx + dy * dy;
88  }
89 
90  struct Compare
91  {
92  // must implement > since priority_queue looks for largest element
93  bool operator()(SeedRgPixel const & l,
94  SeedRgPixel const & r) const
95  {
96  if(r.cost_ == l.cost_)
97  {
98  if(r.dist_ == l.dist_) return r.count_ < l.count_;
99 
100  return r.dist_ < l.dist_;
101  }
102 
103  return r.cost_ < l.cost_;
104  }
105  bool operator()(SeedRgPixel const * l,
106  SeedRgPixel const * r) const
107  {
108  if(r->cost_ == l->cost_)
109  {
110  if(r->dist_ == l->dist_) return r->count_ < l->count_;
111 
112  return r->dist_ < l->dist_;
113  }
114 
115  return r->cost_ < l->cost_;
116  }
117  };
118 
119  struct Allocator
120  {
121  ~Allocator()
122  {
123  while(!freelist_.empty())
124  {
125  delete freelist_.top();
126  freelist_.pop();
127  }
128  }
129 
130  SeedRgPixel *
131  create(Point2D const & location, Point2D const & nearest,
132  COST const & cost, int const & count, int const & label)
133  {
134  if(!freelist_.empty())
135  {
136  SeedRgPixel * res = freelist_.top();
137  freelist_.pop();
138  res->set(location, nearest, cost, count, label);
139  return res;
140  }
141 
142  return new SeedRgPixel(location, nearest, cost, count, label);
143  }
144 
145  void dismiss(SeedRgPixel * p)
146  {
147  freelist_.push(p);
148  }
149 
150  std::stack<SeedRgPixel<COST> *> freelist_;
151  };
152 };
153 
154 struct UnlabelWatersheds
155 {
156  int operator()(int label) const
157  {
158  return label < 0 ? 0 : label;
159  }
160 };
161 
162 } // namespace detail
163 
164 /** \addtogroup SeededRegionGrowing Region Segmentation Algorithms
165  Region growing, watersheds, and voronoi tesselation
166 */
167 //@{
168 
169 /********************************************************/
170 /* */
171 /* seededRegionGrowing */
172 /* */
173 /********************************************************/
174 
175 /** Choose between different types of Region Growing */
176 enum SRGType {
177  CompleteGrow = 0,
178  KeepContours = 1,
179  StopAtThreshold = 2,
180  SRGWatershedLabel = -1
181 };
182 
183 /** \brief Region Segmentation by means of Seeded Region Growing.
184 
185  This algorithm implements seeded region growing as described in
186 
187  R. Adams, L. Bischof: "<em> Seeded Region Growing</em>", IEEE Trans. on Pattern
188  Analysis and Maschine Intelligence, vol 16, no 6, 1994, and
189 
190  Ullrich K&ouml;the:
191  <em><a href="http://hci.iwr.uni-heidelberg.de/people/ukoethe/papers/index.php#cite_primary_segmentation">Primary Image Segmentation</a></em>,
192  in: G. Sagerer, S.
193  Posch, F. Kummert (eds.): Mustererkennung 1995, Proc. 17. DAGM-Symposium,
194  Springer 1995
195 
196  The seed image is a partly segmented image which contains uniquely
197  labeled regions (the seeds) and unlabeled pixels (the candidates, label 0).
198  The highest seed label found in the seed image is returned by the algorithm.
199 
200  Seed regions can be as large as you wish and as small as one pixel. If
201  there are no candidates, the algorithm will simply copy the seed image
202  into the output image. Otherwise it will aggregate the candidates into
203  the existing regions so that a cost function is minimized.
204  Candidates are taken from the neighborhood of the already assigned pixels,
205  where the type of neighborhood is determined by parameter <tt>neighborhood</tt>
206  which can take the values <tt>FourNeighborCode()</tt> (the default)
207  or <tt>EightNeighborCode()</tt>. The algorithm basically works as follows
208  (illustrated for 4-neighborhood, but 8-neighborhood works in the same way):
209 
210  <ol>
211 
212  <li> Find all candidate pixels that are 4-adjacent to a seed region.
213  Calculate the cost for aggregating each candidate into its adjacent region
214  and put the candidates into a priority queue.
215 
216  <li> While( priority queue is not empty and termination criterion is not fulfilled)
217 
218  <ol>
219 
220  <li> Take the candidate with least cost from the queue. If it has not
221  already been merged, merge it with it's adjacent region.
222 
223  <li> Put all candidates that are 4-adjacent to the pixel just processed
224  into the priority queue.
225 
226  </ol>
227 
228  </ol>
229 
230  <tt>SRGType</tt> can take the following values:
231 
232  <DL>
233  <DT><tt>CompleteGrow</tt> <DD> produce a complete tesselation of the volume (default).
234  <DT><tt>KeepContours</tt> <DD> keep a 1-voxel wide unlabeled contour between all regions.
235  <DT><tt>StopAtThreshold</tt> <DD> stop when the boundary indicator values exceed the
236  threshold given by parameter <tt>max_cost</tt>.
237  <DT><tt>KeepContours | StopAtThreshold</tt> <DD> keep 1-voxel wide contour and stop at given <tt>max_cost</tt>.
238  </DL>
239 
240  The cost is determined jointly by the source image and the
241  region statistics functor. The source image contains feature values for each
242  pixel which will be used by the region statistics functor to calculate and
243  update statistics for each region and to calculate the cost for each
244  candidate. The <TT>RegionStatisticsArray</TT> must be compatible to the
245  \ref ArrayOfRegionStatistics functor and contains an <em> array</em> of
246  statistics objects for each region. The indices must correspond to the
247  labels of the seed regions. The statistics for the initial regions must have
248  been calculated prior to calling <TT>seededRegionGrowing()</TT> (for example by
249  means of \ref inspectTwoImagesIf()).
250 
251  For each candidate
252  <TT>x</TT> that is adjacent to region <TT>i</TT>, the algorithm will call
253  <TT>stats[i].cost(as(x))</TT> to get the cost (where <TT>x</TT> is a <TT>SrcIterator</TT>
254  and <TT>as</TT> is
255  the SrcAccessor). When a candidate has been merged with a region, the
256  statistics are updated by calling <TT>stats[i].operator()(as(x))</TT>. Since
257  the <TT>RegionStatisticsArray</TT> is passed by reference, this will overwrite
258  the original statistics.
259 
260  If a candidate could be merged into more than one regions with identical
261  cost, the algorithm will favour the nearest region. If <tt>StopAtThreshold</tt> is active,
262  and the cost of the current candidate at any point in the algorithm exceeds the optional
263  <tt>max_cost</tt> value (which defaults to <tt>NumericTraits<double>::max()</tt>),
264  region growing is aborted, and all voxels not yet assigned to a region remain unlabeled.
265 
266  In some cases, the cost only depends on the feature value of the current
267  pixel. Then the update operation will simply be a no-op, and the <TT>cost()</TT>
268  function returns its argument. This behavior is implemented by the
269  \ref SeedRgDirectValueFunctor. With <tt>SRGType == KeepContours</tt>,
270  this is equivalent to the watershed algorithm.
271 
272  <b> Declarations:</b>
273 
274  pass arguments explicitly:
275  \code
276  namespace vigra {
277  template <class SrcIterator, class SrcAccessor,
278  class SeedImageIterator, class SeedAccessor,
279  class DestIterator, class DestAccessor,
280  class RegionStatisticsArray, class Neighborhood>
281  typename SeedAccessor::value_type
282  seededRegionGrowing(SrcIterator srcul, SrcIterator srclr, SrcAccessor as,
283  SeedImageIterator seedsul, SeedAccessor aseeds,
284  DestIterator destul, DestAccessor ad,
285  RegionStatisticsArray & stats,
286  SRGType srgType = CompleteGrow,
287  Neighborhood neighborhood = FourNeighborCode(),
288  double max_cost = NumericTraits<double>::max());
289  }
290  \endcode
291 
292  use argument objects in conjunction with \ref ArgumentObjectFactories :
293  \code
294  namespace vigra {
295  template <class SrcIterator, class SrcAccessor,
296  class SeedImageIterator, class SeedAccessor,
297  class DestIterator, class DestAccessor,
298  class RegionStatisticsArray, class Neighborhood>
299  typename SeedAccessor::value_type
300  seededRegionGrowing(triple<SrcIterator, SrcIterator, SrcAccessor> src,
301  pair<SeedImageIterator, SeedAccessor> seeds,
302  pair<DestIterator, DestAccessor> dest,
303  RegionStatisticsArray & stats,
304  SRGType srgType = CompleteGrow,
305  Neighborhood neighborhood = FourNeighborCode(),
306  double max_cost = NumericTraits<double>::max());
307  }
308  \endcode
309 
310  <b> Usage:</b>
311 
312  <b>\#include</b> <vigra/seededregiongrowing.hxx><br>
313  Namespace: vigra
314 
315  Example: implementation of the voronoi tesselation
316 
317  \code
318  vigra::BImage points(w,h);
319  vigra::FImage dist(x,y);
320 
321  // empty edge image
322  points = 0;
323  dist = 0;
324 
325  int max_region_label = 100;
326 
327  // throw in some random points:
328  for(int i = 1; i <= max_region_label; ++i)
329  points(w * rand() / RAND_MAX , h * rand() / RAND_MAX) = i;
330 
331  // calculate Euclidean distance transform
332  vigra::distanceTransform(srcImageRange(points), destImage(dist), 2);
333 
334  // init statistics functor
335  vigra::ArrayOfRegionStatistics<vigra::SeedRgDirectValueFunctor<float> >
336  stats(max_region_label);
337 
338  // find voronoi region of each point
339  vigra:: seededRegionGrowing(srcImageRange(dist), srcImage(points),
340  destImage(points), stats);
341  \endcode
342 
343  <b> Required Interface:</b>
344 
345  \code
346  SrcIterator src_upperleft, src_lowerright;
347  SeedImageIterator seed_upperleft;
348  DestIterator dest_upperleft;
349 
350  SrcAccessor src_accessor;
351  SeedAccessor seed_accessor;
352  DestAccessor dest_accessor;
353 
354  RegionStatisticsArray stats;
355 
356  // calculate costs
357  RegionStatisticsArray::value_type::cost_type cost =
358  stats[seed_accessor(seed_upperleft)].cost(src_accessor(src_upperleft));
359 
360  // compare costs
361  cost < cost;
362 
363  // update statistics
364  stats[seed_accessor(seed_upperleft)](src_accessor(src_upperleft));
365 
366  // set result
367  dest_accessor.set(seed_accessor(seed_upperleft), dest_upperleft);
368  \endcode
369 
370  Further requirements are determined by the <TT>RegionStatisticsArray</TT>.
371 */
372 doxygen_overloaded_function(template <...> void seededRegionGrowing)
373 
374 template <class SrcIterator, class SrcAccessor,
375  class SeedImageIterator, class SeedAccessor,
376  class DestIterator, class DestAccessor,
377  class RegionStatisticsArray, class Neighborhood>
378 typename SeedAccessor::value_type
379 seededRegionGrowing(SrcIterator srcul,
380  SrcIterator srclr, SrcAccessor as,
381  SeedImageIterator seedsul, SeedAccessor aseeds,
382  DestIterator destul, DestAccessor ad,
383  RegionStatisticsArray & stats,
384  SRGType srgType,
385  Neighborhood,
386  double max_cost)
387 {
388  int w = srclr.x - srcul.x;
389  int h = srclr.y - srcul.y;
390  int count = 0;
391 
392  SrcIterator isy = srcul, isx = srcul; // iterators for the src image
393 
394  typedef typename SeedAccessor::value_type LabelType;
395  typedef typename RegionStatisticsArray::value_type RegionStatistics;
396  typedef typename RegionStatistics::cost_type CostType;
397  typedef detail::SeedRgPixel<CostType> Pixel;
398 
399  typename Pixel::Allocator allocator;
400 
401  typedef std::priority_queue<Pixel *, std::vector<Pixel *>,
402  typename Pixel::Compare> SeedRgPixelHeap;
403 
404  // copy seed image in an image with border
405  IImage regions(w+2, h+2);
406  IImage::Iterator ir = regions.upperLeft() + Diff2D(1,1);
407  IImage::Iterator iry, irx;
408 
409  initImageBorder(destImageRange(regions), 1, SRGWatershedLabel);
410  copyImage(seedsul, seedsul+Diff2D(w,h), aseeds, ir, regions.accessor());
411 
412  // allocate and init memory for the results
413 
414  SeedRgPixelHeap pheap;
415  int cneighbor, maxRegionLabel = 0;
416 
417  typedef typename Neighborhood::Direction Direction;
418  int directionCount = Neighborhood::DirectionCount;
419 
420  Point2D pos(0,0);
421  for(isy=srcul, iry=ir, pos.y=0; pos.y<h;
422  ++pos.y, ++isy.y, ++iry.y)
423  {
424  for(isx=isy, irx=iry, pos.x=0; pos.x<w;
425  ++pos.x, ++isx.x, ++irx.x)
426  {
427  if(*irx == 0)
428  {
429  // find candidate pixels for growing and fill heap
430  for(int i=0; i<directionCount; i++)
431  {
432  // cneighbor = irx[dist[i]];
433  cneighbor = irx[Neighborhood::diff((Direction)i)];
434  if(cneighbor > 0)
435  {
436  CostType cost = stats[cneighbor].cost(as(isx));
437 
438  Pixel * pixel =
439  allocator.create(pos, pos+Neighborhood::diff((Direction)i), cost, count++, cneighbor);
440  pheap.push(pixel);
441  }
442  }
443  }
444  else
445  {
446  vigra_precondition((LabelType)*irx <= stats.maxRegionLabel(),
447  "seededRegionGrowing(): Largest label exceeds size of RegionStatisticsArray.");
448  if(maxRegionLabel < *irx)
449  maxRegionLabel = *irx;
450  }
451  }
452  }
453 
454  // perform region growing
455  while(pheap.size() != 0)
456  {
457  Pixel * pixel = pheap.top();
458  pheap.pop();
459 
460  Point2D pos = pixel->location_;
461  Point2D nearest = pixel->nearest_;
462  int lab = pixel->label_;
463  CostType cost = pixel->cost_;
464 
465  allocator.dismiss(pixel);
466 
467  if((srgType & StopAtThreshold) != 0 && cost > max_cost)
468  break;
469 
470  irx = ir + pos;
471  isx = srcul + pos;
472 
473  if(*irx) // already labelled region / watershed?
474  continue;
475 
476  if((srgType & KeepContours) != 0)
477  {
478  for(int i=0; i<directionCount; i++)
479  {
480  cneighbor = irx[Neighborhood::diff((Direction)i)];
481  if((cneighbor>0) && (cneighbor != lab))
482  {
483  lab = SRGWatershedLabel;
484  break;
485  }
486  }
487  }
488 
489  *irx = lab;
490 
491  if((srgType & KeepContours) == 0 || lab > 0)
492  {
493  // update statistics
494  stats[*irx](as(isx));
495 
496  // search neighborhood
497  // second pass: find new candidate pixels
498  for(int i=0; i<directionCount; i++)
499  {
500  if(irx[Neighborhood::diff((Direction)i)] == 0)
501  {
502  CostType cost = stats[lab].cost(as(isx, Neighborhood::diff((Direction)i)));
503 
504  Pixel * new_pixel =
505  allocator.create(pos+Neighborhood::diff((Direction)i), nearest, cost, count++, lab);
506  pheap.push(new_pixel);
507  }
508  }
509  }
510  }
511 
512  // free temporary memory
513  while(pheap.size() != 0)
514  {
515  allocator.dismiss(pheap.top());
516  pheap.pop();
517  }
518 
519  // write result
520  transformImage(ir, ir+Point2D(w,h), regions.accessor(), destul, ad,
521  detail::UnlabelWatersheds());
522 
523  return (LabelType)maxRegionLabel;
524 }
525 
526 template <class SrcIterator, class SrcAccessor,
527  class SeedImageIterator, class SeedAccessor,
528  class DestIterator, class DestAccessor,
529  class RegionStatisticsArray, class Neighborhood>
530 inline typename SeedAccessor::value_type
531 seededRegionGrowing(SrcIterator srcul,
532  SrcIterator srclr, SrcAccessor as,
533  SeedImageIterator seedsul, SeedAccessor aseeds,
534  DestIterator destul, DestAccessor ad,
535  RegionStatisticsArray & stats,
536  SRGType srgType,
537  Neighborhood n)
538 {
539  return seededRegionGrowing(srcul, srclr, as,
540  seedsul, aseeds,
541  destul, ad,
542  stats, srgType, n, NumericTraits<double>::max());
543 }
544 
545 
546 
547 template <class SrcIterator, class SrcAccessor,
548  class SeedImageIterator, class SeedAccessor,
549  class DestIterator, class DestAccessor,
550  class RegionStatisticsArray>
551 inline typename SeedAccessor::value_type
552 seededRegionGrowing(SrcIterator srcul,
553  SrcIterator srclr, SrcAccessor as,
554  SeedImageIterator seedsul, SeedAccessor aseeds,
555  DestIterator destul, DestAccessor ad,
556  RegionStatisticsArray & stats,
557  SRGType srgType)
558 {
559  return seededRegionGrowing(srcul, srclr, as,
560  seedsul, aseeds,
561  destul, ad,
562  stats, srgType, FourNeighborCode());
563 }
564 
565 template <class SrcIterator, class SrcAccessor,
566  class SeedImageIterator, class SeedAccessor,
567  class DestIterator, class DestAccessor,
568  class RegionStatisticsArray>
569 inline typename SeedAccessor::value_type
570 seededRegionGrowing(SrcIterator srcul,
571  SrcIterator srclr, SrcAccessor as,
572  SeedImageIterator seedsul, SeedAccessor aseeds,
573  DestIterator destul, DestAccessor ad,
574  RegionStatisticsArray & stats)
575 {
576  return seededRegionGrowing(srcul, srclr, as,
577  seedsul, aseeds,
578  destul, ad,
579  stats, CompleteGrow);
580 }
581 
582 template <class SrcIterator, class SrcAccessor,
583  class SeedImageIterator, class SeedAccessor,
584  class DestIterator, class DestAccessor,
585  class RegionStatisticsArray, class Neighborhood>
586 inline typename SeedAccessor::value_type
587 seededRegionGrowing(triple<SrcIterator, SrcIterator, SrcAccessor> img1,
588  pair<SeedImageIterator, SeedAccessor> img3,
589  pair<DestIterator, DestAccessor> img4,
590  RegionStatisticsArray & stats,
591  SRGType srgType,
592  Neighborhood n,
593  double max_cost)
594 {
595  return seededRegionGrowing(img1.first, img1.second, img1.third,
596  img3.first, img3.second,
597  img4.first, img4.second,
598  stats, srgType, n, max_cost);
599 }
600 
601 template <class SrcIterator, class SrcAccessor,
602  class SeedImageIterator, class SeedAccessor,
603  class DestIterator, class DestAccessor,
604  class RegionStatisticsArray, class Neighborhood>
605 inline typename SeedAccessor::value_type
606 seededRegionGrowing(triple<SrcIterator, SrcIterator, SrcAccessor> img1,
607  pair<SeedImageIterator, SeedAccessor> img3,
608  pair<DestIterator, DestAccessor> img4,
609  RegionStatisticsArray & stats,
610  SRGType srgType,
611  Neighborhood n)
612 {
613  return seededRegionGrowing(img1.first, img1.second, img1.third,
614  img3.first, img3.second,
615  img4.first, img4.second,
616  stats, srgType, n, NumericTraits<double>::max());
617 }
618 
619 template <class SrcIterator, class SrcAccessor,
620  class SeedImageIterator, class SeedAccessor,
621  class DestIterator, class DestAccessor,
622  class RegionStatisticsArray>
623 inline typename SeedAccessor::value_type
624 seededRegionGrowing(triple<SrcIterator, SrcIterator, SrcAccessor> img1,
625  pair<SeedImageIterator, SeedAccessor> img3,
626  pair<DestIterator, DestAccessor> img4,
627  RegionStatisticsArray & stats,
628  SRGType srgType)
629 {
630  return seededRegionGrowing(img1.first, img1.second, img1.third,
631  img3.first, img3.second,
632  img4.first, img4.second,
633  stats, srgType, FourNeighborCode());
634 }
635 
636 template <class SrcIterator, class SrcAccessor,
637  class SeedImageIterator, class SeedAccessor,
638  class DestIterator, class DestAccessor,
639  class RegionStatisticsArray>
640 inline typename SeedAccessor::value_type
641 seededRegionGrowing(triple<SrcIterator, SrcIterator, SrcAccessor> img1,
642  pair<SeedImageIterator, SeedAccessor> img3,
643  pair<DestIterator, DestAccessor> img4,
644  RegionStatisticsArray & stats)
645 {
646  return seededRegionGrowing(img1.first, img1.second, img1.third,
647  img3.first, img3.second,
648  img4.first, img4.second,
649  stats, CompleteGrow);
650 }
651 
652 template <class SrcIterator, class SrcAccessor,
653  class DestIterator, class DestAccessor,
654  class RegionStatisticsArray, class Neighborhood>
655 typename DestAccessor::value_type
656 fastSeededRegionGrowing(SrcIterator srcul, SrcIterator srclr, SrcAccessor as,
657  DestIterator destul, DestAccessor ad,
658  RegionStatisticsArray & stats,
659  SRGType srgType,
660  Neighborhood,
661  double max_cost,
662  std::ptrdiff_t bucket_count = 256)
663 {
664  typedef typename DestAccessor::value_type LabelType;
665 
666  vigra_precondition((srgType & KeepContours) == 0,
667  "fastSeededRegionGrowing(): the turbo algorithm doesn't support 'KeepContours', sorry.");
668 
669  int w = srclr.x - srcul.x;
670  int h = srclr.y - srcul.y;
671 
672  SrcIterator isy = srcul, isx = srcul; // iterators for the src image
673  DestIterator idy = destul, idx = destul; // iterators for the dest image
674 
675  BucketQueue<Point2D, true> pqueue(bucket_count);
676  LabelType maxRegionLabel = 0;
677 
678  Point2D pos(0,0);
679  for(isy=srcul, idy = destul, pos.y=0; pos.y<h; ++pos.y, ++isy.y, ++idy.y)
680  {
681  for(isx=isy, idx=idy, pos.x=0; pos.x<w; ++pos.x, ++isx.x, ++idx.x)
682  {
683  LabelType label = ad(idx);
684  if(label != 0)
685  {
686  vigra_precondition(label <= stats.maxRegionLabel(),
687  "fastSeededRegionGrowing(): Largest label exceeds size of RegionStatisticsArray.");
688 
689  if(maxRegionLabel < label)
690  maxRegionLabel = label;
691 
692  AtImageBorder atBorder = isAtImageBorder(pos.x, pos.y, w, h);
693  if(atBorder == NotAtBorder)
694  {
695  NeighborhoodCirculator<DestIterator, Neighborhood> c(idx), cend(c);
696  do
697  {
698  if(ad(c) == 0)
699  {
700  std::ptrdiff_t priority = (std::ptrdiff_t)stats[label].cost(as(isx));
701  pqueue.push(pos, priority);
702  break;
703  }
704  }
705  while(++c != cend);
706  }
707  else
708  {
709  RestrictedNeighborhoodCirculator<DestIterator, Neighborhood>
710  c(idx, atBorder), cend(c);
711  do
712  {
713  if(ad(c) == 0)
714  {
715  std::ptrdiff_t priority = (std::ptrdiff_t)stats[label].cost(as(isx));
716  pqueue.push(pos, priority);
717  break;
718  }
719  }
720  while(++c != cend);
721  }
722  }
723  }
724  }
725 
726  // perform region growing
727  while(!pqueue.empty())
728  {
729  Point2D pos = pqueue.top();
730  std::ptrdiff_t cost = pqueue.topPriority();
731  pqueue.pop();
732 
733  if((srgType & StopAtThreshold) != 0 && cost > max_cost)
734  break;
735 
736  idx = destul + pos;
737  isx = srcul + pos;
738 
739  std::ptrdiff_t label = ad(idx);
740 
741  AtImageBorder atBorder = isAtImageBorder(pos.x, pos.y, w, h);
742  if(atBorder == NotAtBorder)
743  {
744  NeighborhoodCirculator<DestIterator, Neighborhood> c(idx), cend(c);
745 
746  do
747  {
748  std::ptrdiff_t nlabel = ad(c);
749  if(nlabel == 0)
750  {
751  ad.set(label, idx, c.diff());
752  std::ptrdiff_t priority =
753  std::max((std::ptrdiff_t)stats[label].cost(as(isx, c.diff())), cost);
754  pqueue.push(pos+c.diff(), priority);
755  }
756  }
757  while(++c != cend);
758  }
759  else
760  {
761  RestrictedNeighborhoodCirculator<DestIterator, Neighborhood>
762  c(idx, atBorder), cend(c);
763  do
764  {
765  std::ptrdiff_t nlabel = ad(c);
766  if(nlabel == 0)
767  {
768  ad.set(label, idx, c.diff());
769  std::ptrdiff_t priority =
770  std::max((std::ptrdiff_t)stats[label].cost(as(isx, c.diff())), cost);
771  pqueue.push(pos+c.diff(), priority);
772  }
773  }
774  while(++c != cend);
775  }
776  }
777 
778  return maxRegionLabel;
779 }
780 
781 template <class SrcIterator, class SrcAccessor,
782  class DestIterator, class DestAccessor,
783  class RegionStatisticsArray, class Neighborhood>
784 inline typename DestAccessor::value_type
785 fastSeededRegionGrowing(SrcIterator srcul, SrcIterator srclr, SrcAccessor as,
786  DestIterator destul, DestAccessor ad,
787  RegionStatisticsArray & stats,
788  SRGType srgType,
789  Neighborhood n)
790 {
791  return fastSeededRegionGrowing(srcul, srclr, as,
792  destul, ad,
793  stats, srgType, n, NumericTraits<double>::max(), 256);
794 }
795 
796 template <class SrcIterator, class SrcAccessor,
797  class DestIterator, class DestAccessor,
798  class RegionStatisticsArray>
799 inline typename DestAccessor::value_type
800 fastSeededRegionGrowing(SrcIterator srcul, SrcIterator srclr, SrcAccessor as,
801  DestIterator destul, DestAccessor ad,
802  RegionStatisticsArray & stats,
803  SRGType srgType)
804 {
805  return fastSeededRegionGrowing(srcul, srclr, as,
806  destul, ad,
807  stats, srgType, FourNeighborCode());
808 }
809 
810 template <class SrcIterator, class SrcAccessor,
811  class DestIterator, class DestAccessor,
812  class RegionStatisticsArray>
813 inline typename DestAccessor::value_type
814 fastSeededRegionGrowing(SrcIterator srcul, SrcIterator srclr, SrcAccessor as,
815  DestIterator destul, DestAccessor ad,
816  RegionStatisticsArray & stats)
817 {
818  return fastSeededRegionGrowing(srcul, srclr, as,
819  destul, ad,
820  stats, CompleteGrow);
821 }
822 
823 template <class SrcIterator, class SrcAccessor,
824  class DestIterator, class DestAccessor,
825  class RegionStatisticsArray, class Neighborhood>
826 inline typename DestAccessor::value_type
827 fastSeededRegionGrowing(triple<SrcIterator, SrcIterator, SrcAccessor> src,
828  pair<DestIterator, DestAccessor> dest,
829  RegionStatisticsArray & stats,
830  SRGType srgType,
831  Neighborhood n,
832  double max_cost,
833  std::ptrdiff_t bucket_count = 256)
834 {
835  return fastSeededRegionGrowing(src.first, src.second, src.third,
836  dest.first, dest.second,
837  stats, srgType, n, max_cost, bucket_count);
838 }
839 
840 /********************************************************/
841 /* */
842 /* SeedRgDirectValueFunctor */
843 /* */
844 /********************************************************/
845 
846 /** \brief Statistics functor to be used for seeded region growing.
847 
848  This functor can be used if the cost of a candidate during
849  \ref seededRegionGrowing() is equal to the feature value of that
850  candidate and does not depend on properties of the region it is going to
851  be merged with.
852 
853  <b>\#include</b> <vigra/seededregiongrowing.hxx><br>
854  Namespace: vigra
855 
856 
857  <b> Required Interface:</b>
858 
859  no requirements
860 */
861 template <class Value>
863 {
864  public:
865  /** the functor's argument type
866  */
867  typedef Value argument_type;
868 
869  /** the functor's result type (unused, only necessary for
870  use of SeedRgDirectValueFunctor in \ref vigra::ArrayOfRegionStatistics
871  */
872  typedef Value result_type;
873 
874  /** \deprecated use argument_type
875  */
876  typedef Value value_type;
877 
878  /** the return type of the cost() function
879  */
880  typedef Value cost_type;
881 
882  /** Do nothing (since we need not update region statistics).
883  */
884  void operator()(argument_type const &) const {}
885 
886  /** Return argument (since cost is identical to feature value)
887  */
888  cost_type const & cost(argument_type const & v) const
889  {
890  return v;
891  }
892 };
893 
894 //@}
895 
896 } // namespace vigra
897 
898 #endif // VIGRA_SEEDEDREGIONGROWING_HXX
899 

© 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.9.0 (Tue Sep 24 2013)