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

watersheds3d.hxx
1 /************************************************************************/
2 /* */
3 /* Copyright 2006-2007 by F. Heinrich, B. Seppke, 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_watersheds3D_HXX
37 #define VIGRA_watersheds3D_HXX
38 
39 #include "voxelneighborhood.hxx"
40 #include "multi_array.hxx"
41 #include "multi_localminmax.hxx"
42 #include "labelvolume.hxx"
43 #include "seededregiongrowing3d.hxx"
44 #include "watersheds.hxx"
45 
46 namespace vigra
47 {
48 
49 template <class SrcIterator, class SrcAccessor, class SrcShape,
50  class DestIterator, class DestAccessor, class Neighborhood3D>
51 int preparewatersheds3D( SrcIterator s_Iter, SrcShape srcShape, SrcAccessor sa,
52  DestIterator d_Iter, DestAccessor da, Neighborhood3D)
53 {
54  //basically needed for iteration and border-checks
55  int w = srcShape[0], h = srcShape[1], d = srcShape[2];
56  int x,y,z, local_min_count=0;
57 
58  //declare and define Iterators for all three dims at src
59  SrcIterator zs = s_Iter;
60  SrcIterator ys(zs);
61  SrcIterator xs(ys);
62 
63  //Declare Iterators for all three dims at dest
64  DestIterator zd = d_Iter;
65 
66  for(z = 0; z != d; ++z, ++zs.dim2(), ++zd.dim2())
67  {
68  ys = zs;
69  DestIterator yd(zd);
70 
71  for(y = 0; y != h; ++y, ++ys.dim1(), ++yd.dim1())
72  {
73  xs = ys;
74  DestIterator xd(yd);
75 
76  for(x = 0; x != w; ++x, ++xs.dim0(), ++xd.dim0())
77  {
78  AtVolumeBorder atBorder = isAtVolumeBorder(x,y,z,w,h,d);
79  typename SrcAccessor::value_type v = sa(xs);
80  // the following choice causes minima to point
81  // to their lowest neighbor -- would this be better???
82  // typename SrcAccessor::value_type v = NumericTraits<typename SrcAccessor::value_type>::max();
83  int o = 0; // means center is minimum
84  typename SrcAccessor::value_type my_v = v;
85  if(atBorder == NotAtBorder)
86  {
87  NeighborhoodCirculator<SrcIterator, Neighborhood3D> c(xs), cend(c);
88 
89  do {
90  if(sa(c) < v)
91  {
92  v = sa(c);
93  o = c.directionBit();
94  }
95  else if(sa(c) == my_v && my_v == v)
96  {
97  o = o | c.directionBit();
98  }
99  }
100  while(++c != cend);
101  }
102  else
103  {
104  RestrictedNeighborhoodCirculator<SrcIterator, Neighborhood3D> c(xs, atBorder), cend(c);
105  do {
106  if(sa(c) < v)
107  {
108  v = sa(c);
109  o = c.directionBit();
110  }
111  else if(sa(c) == my_v && my_v == v)
112  {
113  o = o | c.directionBit();
114  }
115  }
116  while(++c != cend);
117  }
118  if (o==0) local_min_count++;
119  da.set(o, xd);
120  }//end x-iteration
121  }//end y-iteration
122  }//end z-iteration
123  return local_min_count;
124 }
125 
126 template <class SrcIterator, class SrcAccessor,class SrcShape,
127  class DestIterator, class DestAccessor,
128  class Neighborhood3D>
129 unsigned int watershedLabeling3D( SrcIterator s_Iter, SrcShape srcShape, SrcAccessor sa,
130  DestIterator d_Iter, DestAccessor da,
131  Neighborhood3D)
132 {
133  typedef typename DestAccessor::value_type LabelType;
134 
135  //basically needed for iteration and border-checks
136  int w = srcShape[0], h = srcShape[1], d = srcShape[2];
137  int x,y,z;
138 
139  //declare and define Iterators for all three dims at src
140  SrcIterator zs = s_Iter;
141  DestIterator zd = d_Iter;
142 
143  // temporary image to store region labels
144  detail::UnionFindArray<LabelType> labels;
145 
146  // initialize the neighborhood traversers
147  NeighborOffsetCirculator<Neighborhood3D> nc(Neighborhood3D::CausalFirst);
148  NeighborOffsetCirculator<Neighborhood3D> nce(Neighborhood3D::CausalLast);
149  ++nce;
150  // pass 1: scan image from upper left front to lower right back
151  // to find connected components
152 
153  // Each component will be represented by a tree of pixels. Each
154  // pixel contains the scan order address of its parent in the
155  // tree. In order for pass 2 to work correctly, the parent must
156  // always have a smaller scan order address than the child.
157  // Therefore, we can merge trees only at their roots, because the
158  // root of the combined tree must have the smallest scan order
159  // address among all the tree's pixels/ nodes. The root of each
160  // tree is distinguished by pointing to itself (it contains its
161  // own scan order address). This condition is enforced whenever a
162  // new region is found or two regions are merged
163  for(z = 0; z != d; ++z, ++zs.dim2(), ++zd.dim2())
164  {
165  SrcIterator ys = zs;
166  DestIterator yd = zd;
167 
168  for(y = 0; y != h; ++y, ++ys.dim1(), ++yd.dim1())
169  {
170  SrcIterator xs = ys;
171  DestIterator xd = yd;
172 
173  for(x = 0; x != w; ++x, ++xs.dim0(), ++xd.dim0())
174  {
175  LabelType currentLabel = labels.nextFreeLabel(); // default: new region
176 
177  //check whether there is a special border treatment to be used or not
178  AtVolumeBorder atBorder = isAtVolumeBorderCausal(x,y,z,w,h,d);
179 
180  //We are not at the border!
181  if(atBorder == NotAtBorder)
182  {
183 
184  nc = NeighborOffsetCirculator<Neighborhood3D>(Neighborhood3D::CausalFirst);
185 
186  do
187  {
188  // Direction of NTraversr Neighbor's direction bit is pointing
189  // = Direction of voxel towards us?
190  if((sa(xs) & nc.directionBit()) || (sa(xs,*nc) & nc.oppositeDirectionBit()))
191  {
192  currentLabel = labels.makeUnion(da(xd,*nc), currentLabel);
193  }
194  ++nc;
195  }while(nc!=nce);
196  }
197  //we are at a border - handle this!!
198  else
199  {
200  nc = NeighborOffsetCirculator<Neighborhood3D>(Neighborhood3D::nearBorderDirectionsCausal(atBorder,0));
201  int j=0;
202  while(nc.direction() != Neighborhood3D::Error)
203  {
204  // Direction of NTraversr Neighbor's direction bit is pointing
205  // = Direction of voxel towards us?
206  if((sa(xs) & nc.directionBit()) || (sa(xs,*nc) & nc.oppositeDirectionBit()))
207  {
208  currentLabel = labels.makeUnion(da(xd,*nc), currentLabel);
209  }
210  nc.turnTo(Neighborhood3D::nearBorderDirectionsCausal(atBorder,++j));
211  }
212  }
213  da.set(labels.finalizeLabel(currentLabel), xd);
214  }
215  }
216  }
217 
218  unsigned int count = labels.makeContiguous();
219 
220  // pass 2: assign one label to each region (tree)
221  // so that labels form a consecutive sequence 1, 2, ...
222  zd = d_Iter;
223  for(z=0; z != d; ++z, ++zd.dim2())
224  {
225  DestIterator yd(zd);
226 
227  for(y=0; y != h; ++y, ++yd.dim1())
228  {
229  DestIterator xd(yd);
230 
231  for(x = 0; x != w; ++x, ++xd.dim0())
232  {
233  da.set(labels[da(xd)], xd);
234  }
235  }
236  }
237  return count;
238 }
239 
240 
241 /** \addtogroup SeededRegionGrowing
242 */
243 //@{
244 
245 /** \brief Generate seeds for watershed computation and seeded region growing.
246 
247  The source image is a boundary indicator such as the gradient magnitude
248  or the trace of the \ref boundaryTensor(). Seeds are generally generated
249  at locations where the boundaryness (i.e. the likelihood of the point being on the
250  boundary) is very small. In particular, seeds can be placed by either
251  looking for local minima (possibly including minimal plateaus) of the boundaryness,
252  of by looking at level sets (i.e. regions where the boundaryness is below a threshold).
253  Both methods can also be combined, so that only minima below a threshold are returned.
254  The particular seeding strategy is specified by the <tt>options</tt> object
255  (see \ref SeedOptions).
256 
257  The pixel type of the input image must be <tt>LessThanComparable</tt>.
258  The pixel type of the output image must be large enough to hold the labels for all seeds.
259  (typically, you will use <tt>UInt32</tt>). The function will label seeds by consecutive integers
260  (starting from 1) and returns the largest label it used.
261 
262  Pass \ref vigra::EightNeighborCode or \ref vigra::FourNeighborCode to determine the
263  neighborhood where pixel values are compared.
264 
265  The function uses accessors.
266 
267  <b> Declarations:</b>
268 
269  pass arguments explicitly:
270  \code
271  namespace vigra {
272  template <class SrcIterator, class SrcAccessor,
273  class DestIterator, class DestAccessor,
274  class Neighborhood = EightNeighborCode>
275  unsigned int
276  generateWatershedSeeds(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
277  DestIterator upperleftd, DestAccessor da,
278  Neighborhood neighborhood = EightNeighborCode(),
279  SeedOptions const & options = SeedOptions());
280  }
281  \endcode
282 
283  use argument objects in conjunction with \ref ArgumentObjectFactories :
284  \code
285  namespace vigra {
286  template <class SrcIterator, class SrcAccessor,
287  class DestIterator, class DestAccessor,
288  class Neighborhood = EightNeighborCode>
289  unsigned int
290  generateWatershedSeeds(triple<SrcIterator, SrcIterator, SrcAccessor> src,
291  pair<DestIterator, DestAccessor> dest,
292  Neighborhood neighborhood = EightNeighborCode(),
293  SeedOptions const & options = SeedOptions());
294  }
295  \endcode
296 
297  <b> Usage:</b>
298 
299  <b>\#include</b> <vigra/watersheds.hxx><br>
300  Namespace: vigra
301 
302  For detailed examples see watershedsRegionGrowing().
303 */
304 doxygen_overloaded_function(template <...> unsigned int generateWatershedSeeds3D)
305 
306 #if 0
307 template <unsigned int N, class T1, class C1, class T2, class C2>
308  class Neighborhood>
309 unsigned int
310 generateWatershedSeeds3D(MultiArrayView<N, T1, C1> in, MultiArrayView<N, T2, C2> out,
311  Neighborhood neighborhood,
312  SeedOptions const & options = SeedOptions())
313 {
314  using namespace functor;
315 
316  vigra_precondition(in.shape() == out.shape(),
317  "generateWatershedSeeds3D(): Shape mismatch between input and output.");
318 
319  vigra_precondition(options.mini != SeedOptions::LevelSets ||
320  options.thresholdIsValid<SrcType>(),
321  "generateWatershedSeeds3D(): SeedOptions.levelSets() must be specified with threshold.");
322 
323  MultiArray<N, UInt8> seeds(in.shape());
324 
325  if(options.mini == SeedOptions::LevelSets)
326  {
327  transformMultiArray(srcMultiArrayRange(in), destMultiArray(seeds),
328  ifThenElse(Arg1() <= Param(options.thresh), Param(1), Param(0)));
329  }
330  else
331  {
332  localMinima(in, seeds,
333  LocalMinmaxOptions().neighborhood(Neighborhood::DirectionCount)
334  .markWith(1.0)
335  .threshold(options.thresh)
336  .allowAtBorder()
337  .allowPlateaus(options.mini == SeedOptions::ExtendedMinima));
338  }
339 
340  return labelVolumeWithBackground(srcMultiArrayRange(seeds), destMultiArray(out),
341  neighborhood, 0);
342 }
343 
344 template <class SrcIterator, class SrcAccessor,
345  class DestIterator, class DestAccessor>
346 inline unsigned int
347 generateWatershedSeeds(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
348  DestIterator upperleftd, DestAccessor da,
349  SeedOptions const & options = SeedOptions())
350 {
351  return generateWatershedSeeds(upperlefts, lowerrights, sa, upperleftd, da,
352  EightNeighborCode(), options);
353 }
354 
355 template <class SrcIterator, class SrcAccessor,
356  class DestIterator, class DestAccessor,
357  class Neighborhood>
358 inline unsigned int
359 generateWatershedSeeds(triple<SrcIterator, SrcIterator, SrcAccessor> src,
360  pair<DestIterator, DestAccessor> dest,
361  Neighborhood neighborhood,
362  SeedOptions const & options = SeedOptions())
363 {
364  return generateWatershedSeeds(src.first, src.second, src.third,
365  dest.first, dest.second,
366  neighborhood, options);
367 }
368 
369 template <class SrcIterator, class SrcAccessor,
370  class DestIterator, class DestAccessor>
371 inline unsigned int
372 generateWatershedSeeds(triple<SrcIterator, SrcIterator, SrcAccessor> src,
373  pair<DestIterator, DestAccessor> dest,
374  SeedOptions const & options = SeedOptions())
375 {
376  return generateWatershedSeeds(src.first, src.second, src.third,
377  dest.first, dest.second,
378  EightNeighborCode(), options);
379 }
380 
381 #endif
382 
383 /********************************************************/
384 /* */
385 /* watersheds3D */
386 /* */
387 /********************************************************/
388 
389 /** \brief Region Segmentation by means of the watershed algorithm.
390 
391  <b> Declarations:</b>
392 
393  pass arguments explicitly:
394  \code
395  namespace vigra {
396  template <class SrcIterator, class SrcAccessor,class SrcShape,
397  class DestIterator, class DestAccessor,
398  class Neighborhood3D>
399  unsigned int watersheds3D(SrcIterator s_Iter, SrcShape srcShape, SrcAccessor sa,
400  DestIterator d_Iter, DestAccessor da,
401  Neighborhood3D neighborhood3D);
402  }
403  \endcode
404 
405  use argument objects in conjunction with \ref ArgumentObjectFactories :
406  \code
407  namespace vigra {
408  template <class SrcIterator, class SrcAccessor,class SrcShape,
409  class DestIterator, class DestAccessor,
410  class Neighborhood3D>
411  unsigned int watersheds3D(triple<SrcIterator, SrcShape, SrcAccessor> src,
412  pair<DestIterator, DestAccessor> dest,
413  Neighborhood3D neighborhood3D);
414  }
415  \endcode
416 
417  use with 3D-Six-Neighborhood:
418  \code
419  namespace vigra {
420 
421  template <class SrcIterator, class SrcAccessor,class SrcShape,
422  class DestIterator, class DestAccessor>
423  unsigned int watersheds3DSix(triple<SrcIterator, SrcShape, SrcAccessor> src,
424  pair<DestIterator, DestAccessor> dest);
425 
426  }
427  \endcode
428 
429  use with 3D-TwentySix-Neighborhood:
430  \code
431  namespace vigra {
432 
433  template <class SrcIterator, class SrcAccessor,class SrcShape,
434  class DestIterator, class DestAccessor>
435  unsigned int watersheds3DTwentySix(triple<SrcIterator, SrcShape, SrcAccessor> src,
436  pair<DestIterator, DestAccessor> dest);
437 
438  }
439  \endcode
440 
441  This function implements the union-find version of the watershed algorithms
442  as described in
443 
444  J. Roerdink, R. Meijster: "<em>The watershed transform: definitions, algorithms,
445  and parallelization strategies</em>", Fundamenta Informaticae, 41:187-228, 2000
446 
447  The source volume is a boundary indicator such as the gradient magnitude
448  of the trace of the \ref boundaryTensor(). Local minima of the boundary indicator
449  are used as region seeds, and all other voxels are recursively assigned to the same
450  region as their lowest neighbor. Pass \ref vigra::NeighborCode3DSix or
451  \ref vigra::NeighborCode3DTwentySix to determine the neighborhood where voxel values
452  are compared. The voxel type of the input volume must be <tt>LessThanComparable</tt>.
453  The function uses accessors.
454 
455  ...probably soon in VIGRA:
456  Note that VIGRA provides an alternative implementation of the watershed transform via
457  \ref seededRegionGrowing3D(). It is slower, but handles plateaus better
458  and allows to keep a one pixel wide boundary between regions.
459 
460  <b> Usage:</b>
461 
462  <b>\#include</b> <vigra/watersheds3D.hxx><br>
463  Namespace: vigra
464 
465  Example: watersheds3D of the gradient magnitude.
466 
467  \code
468  typedef vigra::MultiArray<3,int> IntVolume;
469  typedef vigra::MultiArray<3,double> DVolume;
470  DVolume src(DVolume::difference_type(w,h,d));
471  IntVolume dest(IntVolume::difference_type(w,h,d));
472 
473  float gauss=1;
474 
475  vigra::MultiArray<3, vigra::TinyVector<float,3> > temp(IntVolume::difference_type(w,h,d));
476  vigra::gaussianGradientMultiArray(srcMultiArrayRange(vol),destMultiArray(temp),gauss);
477 
478  IntVolume::iterator temp_iter=temp.begin();
479  for(DVolume::iterator iter=src.begin(); iter!=src.end(); ++iter, ++temp_iter)
480  *iter = norm(*temp_iter);
481 
482  // find 6-connected regions
483  int max_region_label = vigra::watersheds3DSix(srcMultiArrayRange(src), destMultiArray(dest));
484 
485  // find 26-connected regions
486  max_region_label = vigra::watersheds3DTwentySix(srcMultiArrayRange(src), destMultiArray(dest));
487 
488  \endcode
489 
490  <b> Required Interface:</b>
491 
492  \code
493  SrcIterator src_begin;
494  SrcShape src_shape;
495  DestIterator dest_begin;
496 
497  SrcAccessor src_accessor;
498  DestAccessor dest_accessor;
499 
500  // compare src values
501  src_accessor(src_begin) <= src_accessor(src_begin)
502 
503  // set result
504  int label;
505  dest_accessor.set(label, dest_begin);
506  \endcode
507 */
508 doxygen_overloaded_function(template <...> unsigned int watersheds3D)
509 
510 template <class SrcIterator, class SrcAccessor, class SrcShape,
511  class DestIterator, class DestAccessor,
512  class Neighborhood3D>
513 unsigned int watersheds3D( SrcIterator s_Iter, SrcShape srcShape, SrcAccessor sa,
514  DestIterator d_Iter, DestAccessor da, Neighborhood3D neighborhood3D)
515 {
516  //create temporary volume to store the DAG of directions to minima
517  if ((int)Neighborhood3D::DirectionCount>7){ //If we have 3D-TwentySix Neighborhood
518 
519  vigra::MultiArray<3,int> orientationVolume(srcShape);
520 
521  preparewatersheds3D( s_Iter, srcShape, sa,
522  destMultiArray(orientationVolume).first, destMultiArray(orientationVolume).second,
523  neighborhood3D);
524 
525  return watershedLabeling3D( srcMultiArray(orientationVolume).first, srcShape, srcMultiArray(orientationVolume).second,
526  d_Iter, da,
527  neighborhood3D);
528  }
529  else{
530 
531  vigra::MultiArray<3,unsigned char> orientationVolume(srcShape);
532 
533  preparewatersheds3D( s_Iter, srcShape, sa,
534  destMultiArray(orientationVolume).first, destMultiArray(orientationVolume).second,
535  neighborhood3D);
536 
537  return watershedLabeling3D( srcMultiArray(orientationVolume).first, srcShape, srcMultiArray(orientationVolume).second,
538  d_Iter, da,
539  neighborhood3D);
540  }
541 }
542 
543 template <class SrcIterator, class SrcShape, class SrcAccessor,
544  class DestIterator, class DestAccessor>
545 inline unsigned int watersheds3DSix( vigra::triple<SrcIterator, SrcShape, SrcAccessor> src,
546  vigra::pair<DestIterator, DestAccessor> dest)
547 {
548  return watersheds3D(src.first, src.second, src.third, dest.first, dest.second, NeighborCode3DSix());
549 }
550 
551 template <class SrcIterator, class SrcShape, class SrcAccessor,
552  class DestIterator, class DestAccessor>
553 inline unsigned int watersheds3DTwentySix( vigra::triple<SrcIterator, SrcShape, SrcAccessor> src,
554  vigra::pair<DestIterator, DestAccessor> dest)
555 {
556  return watersheds3D(src.first, src.second, src.third, dest.first, dest.second, NeighborCode3DTwentySix());
557 }
558 
559 }//namespace vigra
560 
561 #endif //VIGRA_watersheds3D_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.8.0 (Wed Sep 26 2012)