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

watersheds.hxx
1 /************************************************************************/
2 /* */
3 /* Copyright 1998-2005 by Ullrich Koethe */
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_WATERSHEDS_HXX
37 #define VIGRA_WATERSHEDS_HXX
38 
39 #include <functional>
40 #include "mathutil.hxx"
41 #include "stdimage.hxx"
42 #include "pixelneighborhood.hxx"
43 #include "union_find.hxx"
44 
45 namespace vigra {
46 
47 template <class SrcIterator, class SrcAccessor,
48  class DestIterator, class DestAccessor,
49  class Neighborhood>
50 unsigned int watershedLabeling(SrcIterator upperlefts,
51  SrcIterator lowerrights, SrcAccessor sa,
52  DestIterator upperleftd, DestAccessor da,
53  Neighborhood)
54 {
55  typedef typename DestAccessor::value_type LabelType;
56 
57  int w = lowerrights.x - upperlefts.x;
58  int h = lowerrights.y - upperlefts.y;
59  int x,y;
60 
61  SrcIterator ys(upperlefts);
62  SrcIterator xs(ys);
63  DestIterator yd(upperleftd);
64  DestIterator xd(yd);
65 
66  // temporary image to store region labels
67  detail::UnionFindArray<LabelType> labels;
68 
69  // initialize the neighborhood circulators
70  NeighborOffsetCirculator<Neighborhood> ncstart(Neighborhood::CausalFirst);
71  NeighborOffsetCirculator<Neighborhood> ncstartBorder(Neighborhood::North);
72  NeighborOffsetCirculator<Neighborhood> ncend(Neighborhood::CausalLast);
73  ++ncend;
74  NeighborOffsetCirculator<Neighborhood> ncendBorder(Neighborhood::North);
75  ++ncendBorder;
76 
77  // pass 1: scan image from upper left to lower right
78  // to find connected components
79 
80  // Each component will be represented by a tree of pixels. Each
81  // pixel contains the scan order address of its parent in the
82  // tree. In order for pass 2 to work correctly, the parent must
83  // always have a smaller scan order address than the child.
84  // Therefore, we can merge trees only at their roots, because the
85  // root of the combined tree must have the smallest scan order
86  // address among all the tree's pixels/ nodes. The root of each
87  // tree is distinguished by pointing to itself (it contains its
88  // own scan order address). This condition is enforced whenever a
89  // new region is found or two regions are merged
90  da.set(labels.finalizeLabel(labels.nextFreeLabel()), xd);
91 
92  ++xs.x;
93  ++xd.x;
94  for(x = 1; x != w; ++x, ++xs.x, ++xd.x)
95  {
96  if((sa(xs) & Neighborhood::directionBit(Neighborhood::West)) ||
97  (sa(xs, Neighborhood::west()) & Neighborhood::directionBit(Neighborhood::East)))
98  {
99  da.set(da(xd, Neighborhood::west()), xd);
100  }
101  else
102  {
103  da.set(labels.finalizeLabel(labels.nextFreeLabel()), xd);
104  }
105  }
106 
107  ++ys.y;
108  ++yd.y;
109  for(y = 1; y != h; ++y, ++ys.y, ++yd.y)
110  {
111  xs = ys;
112  xd = yd;
113 
114  for(x = 0; x != w; ++x, ++xs.x, ++xd.x)
115  {
116  NeighborOffsetCirculator<Neighborhood> nc(x == w-1
117  ? ncstartBorder
118  : ncstart);
119  NeighborOffsetCirculator<Neighborhood> nce(x == 0
120  ? ncendBorder
121  : ncend);
122  LabelType currentLabel = labels.nextFreeLabel();
123  for(; nc != nce; ++nc)
124  {
125  if((sa(xs) & nc.directionBit()) || (sa(xs, *nc) & nc.oppositeDirectionBit()))
126  {
127  currentLabel = labels.makeUnion(da(xd,*nc), currentLabel);
128  }
129  }
130  da.set(labels.finalizeLabel(currentLabel), xd);
131  }
132  }
133 
134  unsigned int count = labels.makeContiguous();
135 
136  // pass 2: assign one label to each region (tree)
137  // so that labels form a consecutive sequence 1, 2, ...
138  yd = upperleftd;
139  for(y=0; y != h; ++y, ++yd.y)
140  {
141  DestIterator xd(yd);
142  for(x = 0; x != w; ++x, ++xd.x)
143  {
144  da.set(labels[da(xd)], xd);
145  }
146  }
147  return count;
148 }
149 
150 template <class SrcIterator, class SrcAccessor,
151  class DestIterator, class DestAccessor,
152  class Neighborhood>
153 unsigned int watershedLabeling(triple<SrcIterator, SrcIterator, SrcAccessor> src,
154  pair<DestIterator, DestAccessor> dest,
155  Neighborhood neighborhood)
156 {
157  return watershedLabeling(src.first, src.second, src.third,
158  dest.first, dest.second, neighborhood);
159 }
160 
161 
162 template <class SrcIterator, class SrcAccessor,
163  class DestIterator, class DestAccessor>
164 void prepareWatersheds(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
165  DestIterator upperleftd, DestAccessor da,
167 {
168  int w = lowerrights.x - upperlefts.x;
169  int h = lowerrights.y - upperlefts.y;
170  int x,y;
171 
172  SrcIterator ys(upperlefts);
173  SrcIterator xs(ys);
174 
175  DestIterator yd = upperleftd;
176 
177  for(y = 0; y != h; ++y, ++ys.y, ++yd.y)
178  {
179  xs = ys;
180  DestIterator xd = yd;
181 
182  for(x = 0; x != w; ++x, ++xs.x, ++xd.x)
183  {
184  AtImageBorder atBorder = isAtImageBorder(x,y,w,h);
185  typename SrcAccessor::value_type v = sa(xs);
186  // the following choice causes minima to point
187  // to their lowest neighbor -- would this be better???
188  // typename SrcAccessor::value_type v = NumericTraits<typename SrcAccessor::value_type>::max();
189  int o = 0; // means center is minimum
190  if(atBorder == NotAtBorder)
191  {
192  NeighborhoodCirculator<SrcIterator, FourNeighborCode> c(xs), cend(c);
193  do {
194  if(sa(c) <= v)
195  {
196  v = sa(c);
197  o = c.directionBit();
198  }
199  }
200  while(++c != cend);
201  }
202  else
203  {
204  RestrictedNeighborhoodCirculator<SrcIterator, FourNeighborCode> c(xs, atBorder), cend(c);
205  do {
206  if(sa(c) <= v)
207  {
208  v = sa(c);
209  o = c.directionBit();
210  }
211  }
212  while(++c != cend);
213  }
214  da.set(o, xd);
215  }
216  }
217 }
218 
219 template <class SrcIterator, class SrcAccessor,
220  class DestIterator, class DestAccessor>
221 void prepareWatersheds(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
222  DestIterator upperleftd, DestAccessor da,
224 {
225  int w = lowerrights.x - upperlefts.x;
226  int h = lowerrights.y - upperlefts.y;
227  int x,y;
228 
229  SrcIterator ys(upperlefts);
230  SrcIterator xs(ys);
231 
232  DestIterator yd = upperleftd;
233 
234  for(y = 0; y != h; ++y, ++ys.y, ++yd.y)
235  {
236  xs = ys;
237  DestIterator xd = yd;
238 
239  for(x = 0; x != w; ++x, ++xs.x, ++xd.x)
240  {
241  AtImageBorder atBorder = isAtImageBorder(x,y,w,h);
242  typename SrcAccessor::value_type v = sa(xs);
243  // the following choice causes minima to point
244  // to their lowest neighbor -- would this be better???
245  // typename SrcAccessor::value_type v = NumericTraits<typename SrcAccessor::value_type>::max();
246  int o = 0; // means center is minimum
247  if(atBorder == NotAtBorder)
248  {
249  // handle diagonal and principal neighbors separately
250  // so that principal neighbors are preferred when there are
251  // candidates with equal strength
252  NeighborhoodCirculator<SrcIterator, EightNeighborCode>
254  for(int i = 0; i < 4; ++i, c += 2)
255  {
256  if(sa(c) <= v)
257  {
258  v = sa(c);
259  o = c.directionBit();
260  }
261  }
262  --c;
263  for(int i = 0; i < 4; ++i, c += 2)
264  {
265  if(sa(c) <= v)
266  {
267  v = sa(c);
268  o = c.directionBit();
269  }
270  }
271  }
272  else
273  {
274  RestrictedNeighborhoodCirculator<SrcIterator, EightNeighborCode>
275  c(xs, atBorder), cend(c);
276  do
277  {
278  if(!c.isDiagonal())
279  continue;
280  if(sa(c) <= v)
281  {
282  v = sa(c);
283  o = c.directionBit();
284  }
285  }
286  while(++c != cend);
287  do
288  {
289  if(c.isDiagonal())
290  continue;
291  if(sa(c) <= v)
292  {
293  v = sa(c);
294  o = c.directionBit();
295  }
296  }
297  while(++c != cend);
298  }
299  da.set(o, xd);
300  }
301  }
302 }
303 
304 /** \addtogroup SeededRegionGrowing Region Segmentation Algorithms
305  Region growing, watersheds, and voronoi tesselation
306 */
307 //@{
308 
309 /********************************************************/
310 /* */
311 /* watersheds */
312 /* */
313 /********************************************************/
314 
315 /** \brief Region Segmentation by means of the watershed algorithm.
316 
317  This function implements the union-find version of the watershed algorithms
318  as described in
319 
320  J. Roerdink, R. Meijster: "<em>The watershed transform: definitions, algorithms,
321  and parallelization stretegies</em>", Fundamenta Informaticae, 41:187-228, 2000
322 
323  The source image is a boundary indicator such as the gradient magnitude
324  of the trace of the \ref boundaryTensor(). Local minima of the boundary indicator
325  are used as region seeds, and all other pixels are recursively assigned to the same
326  region as their lowest neighbor. Pass \ref vigra::EightNeighborCode or
327  \ref vigra::FourNeighborCode to determine the neighborhood where pixel values
328  are compared. The pixel type of the input image must be <tt>LessThanComparable</tt>.
329  The function uses accessors.
330 
331  Note that VIGRA provides an alternative implementaion of the watershed transform via
332  \ref seededRegionGrowing(). It is slower, but handles plateaus better
333  and allows to keep a one pixel wide boundary between regions.
334 
335  <b> Declarations:</b>
336 
337  pass arguments explicitly:
338  \code
339  namespace vigra {
340  template <class SrcIterator, class SrcAccessor,
341  class DestIterator, class DestAccessor,
342  class Neighborhood = EightNeighborCode>
343  unsigned int
344  watersheds(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
345  DestIterator upperleftd, DestAccessor da,
346  Neighborhood neighborhood = EightNeighborCode())
347  }
348  \endcode
349 
350  use argument objects in conjunction with \ref ArgumentObjectFactories :
351  \code
352  namespace vigra {
353  template <class SrcIterator, class SrcAccessor,
354  class DestIterator, class DestAccessor,
355  class Neighborhood = EightNeighborCode>
356  unsigned int
357  watersheds(triple<SrcIterator, SrcIterator, SrcAccessor> src,
358  pair<DestIterator, DestAccessor> dest,
359  Neighborhood neighborhood = EightNeighborCode())
360  }
361  \endcode
362 
363  <b> Usage:</b>
364 
365  <b>\#include</b> <<a href="watersheds_8hxx-source.html">vigra/watersheds.hxx</a>><br>
366  Namespace: vigra
367 
368  Example: watersheds of the gradient magnitude.
369 
370  \code
371  vigra::BImage in(w,h);
372  ... // read input data
373 
374  vigra::FImage gradx(x,y), grady(x,y), gradMag(x,y);
375  gaussianGradient(srcImageRange(src), destImage(gradx), destImage(grady), 3.0);
376  combineTwoImages(srcImageRange(gradx), srcImage(grady), destImage(gradMag),
377  vigra::MagnitudeFunctor<float>());
378 
379  // the pixel type of the destination image must be large enough to hold
380  // numbers up to 'max_region_label' to prevent overflow
381  vigra::IImage labeling(x,y);
382  int max_region_label = watersheds(srcImageRange(gradMag), destImage(labeling));
383 
384  \endcode
385 
386  <b> Required Interface:</b>
387 
388  \code
389  SrcImageIterator src_upperleft, src_lowerright;
390  DestImageIterator dest_upperleft;
391 
392  SrcAccessor src_accessor;
393  DestAccessor dest_accessor;
394 
395  // compare src values
396  src_accessor(src_upperleft) <= src_accessor(src_upperleft)
397 
398  // set result
399  int label;
400  dest_accessor.set(label, dest_upperleft);
401  \endcode
402 */
403 doxygen_overloaded_function(template <...> unsigned int watersheds)
404 
405 template <class SrcIterator, class SrcAccessor,
406  class DestIterator, class DestAccessor,
407  class Neighborhood>
408 unsigned int
409 watersheds(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
410  DestIterator upperleftd, DestAccessor da, Neighborhood neighborhood)
411 {
412  SImage orientationImage(lowerrights - upperlefts);
413  SImage::traverser yo = orientationImage.upperLeft();
414 
415  prepareWatersheds(upperlefts, lowerrights, sa,
416  orientationImage.upperLeft(), orientationImage.accessor(), neighborhood);
417  return watershedLabeling(orientationImage.upperLeft(), orientationImage.lowerRight(), orientationImage.accessor(),
418  upperleftd, da, neighborhood);
419 }
420 
421 template <class SrcIterator, class SrcAccessor,
422  class DestIterator, class DestAccessor>
423 inline unsigned int
424 watersheds(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
425  DestIterator upperleftd, DestAccessor da)
426 {
427  return watersheds(upperlefts, lowerrights, sa, upperleftd, da, EightNeighborCode());
428 }
429 
430 template <class SrcIterator, class SrcAccessor,
431  class DestIterator, class DestAccessor,
432  class Neighborhood>
433 inline unsigned int
434 watersheds(triple<SrcIterator, SrcIterator, SrcAccessor> src,
435  pair<DestIterator, DestAccessor> dest, Neighborhood neighborhood)
436 {
437  return watersheds(src.first, src.second, src.third, dest.first, dest.second, neighborhood);
438 }
439 
440 template <class SrcIterator, class SrcAccessor,
441  class DestIterator, class DestAccessor>
442 inline unsigned int
443 watersheds(triple<SrcIterator, SrcIterator, SrcAccessor> src,
444  pair<DestIterator, DestAccessor> dest)
445 {
446  return watersheds(src.first, src.second, src.third, dest.first, dest.second);
447 }
448 
449 //@}
450 
451 } // namespace vigra
452 
453 #endif // VIGRA_WATERSHEDS_HXX

© 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.7.1 (Tue Jul 10 2012)