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

watersheds3d.hxx VIGRA

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  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 currentIndex = labels.nextFreeIndex(); // 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  currentIndex = labels.makeUnion(da(xd,*nc), currentIndex);
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  int dummy = x+(*nc)[0]; // prevents an apparently incorrect optimization in gcc 4.8
205  if (dummy<0)
206  {
207  std::cerr << "internal error " << dummy << std::endl;
208  }
209  // Direction of NTraversr Neighbor's direction bit is pointing
210  // = Direction of voxel towards us?
211  if((sa(xs) & nc.directionBit()) || (sa(xs,*nc) & nc.oppositeDirectionBit()))
212  {
213  currentIndex = labels.makeUnion(da(xd,*nc), currentIndex);
214  }
215  nc.turnTo(Neighborhood3D::nearBorderDirectionsCausal(atBorder,++j));
216  }
217  }
218  da.set(labels.finalizeIndex(currentIndex), xd);
219  }
220  }
221  }
222 
223  unsigned int count = labels.makeContiguous();
224 
225  // pass 2: assign one label to each region (tree)
226  // so that labels form a consecutive sequence 1, 2, ...
227  zd = d_Iter;
228  for(z=0; z != d; ++z, ++zd.dim2())
229  {
230  DestIterator yd(zd);
231 
232  for(y=0; y != h; ++y, ++yd.dim1())
233  {
234  DestIterator xd(yd);
235 
236  for(x = 0; x != w; ++x, ++xd.dim0())
237  {
238  da.set(labels.findLabel(da(xd)), xd);
239  }
240  }
241  }
242  return count;
243 }
244 
245 
246 /** \addtogroup SeededRegionGrowing
247 */
248 //@{
249 
250 /********************************************************/
251 /* */
252 /* watersheds3D */
253 /* */
254 /********************************************************/
255 
256 /** \brief Region Segmentation by means of the watershed algorithm.
257 
258  This function is deprecated, use \ref watershedsMultiArray() instead.
259 
260  <b> Declarations:</b>
261 
262  \deprecatedAPI{watersheds3D}
263  pass \ref MultiIteratorPage "MultiIterators" and \ref DataAccessors :
264  \code
265  namespace vigra {
266  template <class SrcIterator, class SrcAccessor,class SrcShape,
267  class DestIterator, class DestAccessor,
268  class Neighborhood3D>
269  unsigned int watersheds3D(SrcIterator s_Iter, SrcShape srcShape, SrcAccessor sa,
270  DestIterator d_Iter, DestAccessor da,
271  Neighborhood3D neighborhood3D);
272  }
273  \endcode
274  use argument objects in conjunction with \ref ArgumentObjectFactories :
275  \code
276  namespace vigra {
277  template <class SrcIterator, class SrcAccessor,class SrcShape,
278  class DestIterator, class DestAccessor,
279  class Neighborhood3D>
280  unsigned int watersheds3D(triple<SrcIterator, SrcShape, SrcAccessor> src,
281  pair<DestIterator, DestAccessor> dest,
282  Neighborhood3D neighborhood3D);
283  }
284  \endcode
285 
286  use with 3D-Six-Neighborhood:
287  \code
288  namespace vigra {
289 
290  template <class SrcIterator, class SrcAccessor,class SrcShape,
291  class DestIterator, class DestAccessor>
292  unsigned int watersheds3DSix(triple<SrcIterator, SrcShape, SrcAccessor> src,
293  pair<DestIterator, DestAccessor> dest);
294 
295  }
296  \endcode
297 
298  use with 3D-TwentySix-Neighborhood:
299  \code
300  namespace vigra {
301 
302  template <class SrcIterator, class SrcAccessor,class SrcShape,
303  class DestIterator, class DestAccessor>
304  unsigned int watersheds3DTwentySix(triple<SrcIterator, SrcShape, SrcAccessor> src,
305  pair<DestIterator, DestAccessor> dest);
306 
307  }
308  \endcode
309  \deprecatedEnd
310 
311  This function implements the union-find version of the watershed algorithms
312  as described in
313 
314  J. Roerdink, R. Meijster: <em>"The watershed transform: definitions, algorithms,
315  and parallelization strategies"</em>, Fundamenta Informaticae, 41:187-228, 2000
316 
317  The source volume is a boundary indicator such as the gradient magnitude
318  of the trace of the \ref boundaryTensor(). Local minima of the boundary indicator
319  are used as region seeds, and all other voxels are recursively assigned to the same
320  region as their lowest neighbor. Pass \ref vigra::NeighborCode3DSix or
321  \ref vigra::NeighborCode3DTwentySix to determine the neighborhood where voxel values
322  are compared. The voxel type of the input volume must be <tt>LessThanComparable</tt>.
323 
324  <b> Usage:</b>
325 
326  <b>\#include</b> <vigra/watersheds3D.hxx><br>
327  Namespace: vigra
328 
329  Example: watersheds3D of the gradient magnitude.
330 
331  \code
332  Shape3 shape(w, h, d);
333 
334  MultiArray<3, float> src(shape), grad(shape);
335  ...
336 
337  double scale = 1;
338  gaussianGradientMagnitude(src, grad, scale);
339 
340  MultiArray<3, int> labels(shape);
341 
342  // find 6-connected regions
343  int max_region_label = watersheds3DSix(grad, labels);
344 
345  // find 26-connected regions
346  max_region_label = watersheds3DTwentySix(grad, labels);
347  \endcode
348 
349  \deprecatedUsage{watersheds3D}
350  \code
351  typedef vigra::MultiArray<3,int> IntVolume;
352  typedef vigra::MultiArray<3,double> DVolume;
353  DVolume src(DVolume::difference_type(w,h,d));
354  IntVolume dest(IntVolume::difference_type(w,h,d));
355 
356  float gauss=1;
357 
358  vigra::MultiArray<3, vigra::TinyVector<float,3> > temp(IntVolume::difference_type(w,h,d));
359  vigra::gaussianGradientMultiArray(srcMultiArrayRange(vol),destMultiArray(temp),gauss);
360 
361  IntVolume::iterator temp_iter=temp.begin();
362  for(DVolume::iterator iter=src.begin(); iter!=src.end(); ++iter, ++temp_iter)
363  *iter = norm(*temp_iter);
364 
365  // find 6-connected regions
366  int max_region_label = vigra::watersheds3DSix(srcMultiArrayRange(src), destMultiArray(dest));
367 
368  // find 26-connected regions
369  max_region_label = vigra::watersheds3DTwentySix(srcMultiArrayRange(src), destMultiArray(dest));
370 
371  \endcode
372  <b> Required Interface:</b>
373  \code
374  SrcIterator src_begin;
375  SrcShape src_shape;
376  DestIterator dest_begin;
377 
378  SrcAccessor src_accessor;
379  DestAccessor dest_accessor;
380 
381  // compare src values
382  src_accessor(src_begin) <= src_accessor(src_begin)
383 
384  // set result
385  int label;
386  dest_accessor.set(label, dest_begin);
387  \endcode
388  \deprecatedEnd
389 */
390 doxygen_overloaded_function(template <...> unsigned int watersheds3D)
391 
392 template <class SrcIterator, class SrcAccessor, class SrcShape,
393  class DestIterator, class DestAccessor,
394  class Neighborhood3D>
395 unsigned int watersheds3D( SrcIterator s_Iter, SrcShape srcShape, SrcAccessor sa,
396  DestIterator d_Iter, DestAccessor da, Neighborhood3D neighborhood3D)
397 {
398  //create temporary volume to store the DAG of directions to minima
399  if ((int)Neighborhood3D::DirectionCount>7){ //If we have 3D-TwentySix Neighborhood
400 
401  vigra::MultiArray<3,int> orientationVolume(srcShape);
402 
403  preparewatersheds3D( s_Iter, srcShape, sa,
404  destMultiArray(orientationVolume).first, destMultiArray(orientationVolume).second,
405  neighborhood3D);
406 
407  return watershedLabeling3D( srcMultiArray(orientationVolume).first, srcShape, srcMultiArray(orientationVolume).second,
408  d_Iter, da,
409  neighborhood3D);
410  }
411  else{
412 
413  vigra::MultiArray<3,unsigned char> orientationVolume(srcShape);
414 
415  preparewatersheds3D( s_Iter, srcShape, sa,
416  destMultiArray(orientationVolume).first, destMultiArray(orientationVolume).second,
417  neighborhood3D);
418 
419  return watershedLabeling3D( srcMultiArray(orientationVolume).first, srcShape, srcMultiArray(orientationVolume).second,
420  d_Iter, da,
421  neighborhood3D);
422  }
423 }
424 
425 template <class SrcIterator, class SrcShape, class SrcAccessor,
426  class DestIterator, class DestAccessor>
427 inline unsigned int watersheds3DSix( triple<SrcIterator, SrcShape, SrcAccessor> src,
428  pair<DestIterator, DestAccessor> dest)
429 {
430  return watersheds3D(src.first, src.second, src.third, dest.first, dest.second, NeighborCode3DSix());
431 }
432 
433 template <class SrcIterator, class SrcShape, class SrcAccessor,
434  class DestIterator, class DestAccessor>
435 inline unsigned int watersheds3DTwentySix( triple<SrcIterator, SrcShape, SrcAccessor> src,
436  pair<DestIterator, DestAccessor> dest)
437 {
438  return watersheds3D(src.first, src.second, src.third, dest.first, dest.second, NeighborCode3DTwentySix());
439 }
440 
441 template <unsigned int N, class T1, class S1,
442  class T2, class S2>
443 inline unsigned int
444 watersheds3DSix(MultiArrayView<N, T1, S1> const & source,
445  MultiArrayView<N, T2, S2> dest)
446 {
447  vigra_precondition(source.shape() == dest.shape(),
448  "watersheds3DSix(): shape mismatch between input and output.");
449  return watersheds3DSix(srcMultiArrayRange(source), destMultiArray(dest));
450 }
451 
452 template <unsigned int N, class T1, class S1,
453  class T2, class S2>
454 inline unsigned int
455 watersheds3DTwentySix(MultiArrayView<N, T1, S1> const & source,
456  MultiArrayView<N, T2, S2> dest)
457 {
458  vigra_precondition(source.shape() == dest.shape(),
459  "watersheds3DTwentySix(): shape mismatch between input and output.");
460  return watersheds3DTwentySix(srcMultiArrayRange(source), destMultiArray(dest));
461 }
462 
463 }//namespace vigra
464 
465 #endif //VIGRA_watersheds3D_HXX
AtImageBorder AtVolumeBorder
Encode whether a voxel is near the volume border.
Definition: voxelneighborhood.hxx:72
unsigned int watersheds3D(...)
Region Segmentation by means of the watershed algorithm.
Neighborhood3DTwentySix::NeighborCode3D NeighborCode3DTwentySix
Definition: voxelneighborhood.hxx:1646
Main MultiArray class containing the memory management.
Definition: multi_array.hxx:2422
AtVolumeBorder isAtVolumeBorderCausal(int x, int y, int z, int width, int height, int)
Find out whether a voxel is at a scan-order relevant volume border. This function checks if x == 0 or...
Definition: voxelneighborhood.hxx:112
AtVolumeBorder isAtVolumeBorder(int x, int y, int z, int width, int height, int depth)
Find out whether a voxel is at the volume border.
Definition: voxelneighborhood.hxx:82
vigra::GridGraph< N, DirectedTag >::vertex_descriptor source(typename vigra::GridGraph< N, DirectedTag >::edge_descriptor const &e, vigra::GridGraph< N, DirectedTag > const &g)
Get a vertex descriptor for the start vertex of edge e in graph g (API: boost).
Definition: multi_gridgraph.hxx:2943
doxygen_overloaded_function(template<...> void separableConvolveBlockwise) template< unsigned int N
Separated convolution on ChunkedArrays.
Neighborhood3DSix::NeighborCode3D NeighborCode3DSix
Definition: voxelneighborhood.hxx:490
 
Definition: pixelneighborhood.hxx:70

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