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

seededregiongrowing.hxx VIGRA

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

© 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.10.0