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

labelvolume.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_LABELVOLUME_HXX
37 #define VIGRA_LABELVOLUME_HXX
38 
39 
40 #include "voxelneighborhood.hxx"
41 #include "multi_array.hxx"
42 #include "union_find.hxx"
43 
44 namespace vigra{
45 
46 /** \addtogroup Labeling Connected Components Labeling
47  The 3-dimensional connected components algorithms may use either 6 or 26 connectivity.
48  By means of a functor the merge criterion can be defined arbitrarily.
49 */
50 //@{
51 
52 /********************************************************/
53 /* */
54 /* labelVolume */
55 /* */
56 /********************************************************/
57 
58 /** \brief Find the connected components of a segmented volume.
59 
60  <b> Declarations:</b>
61 
62  pass arguments explicitly:
63  \code
64  namespace vigra {
65 
66  template <class SrcIterator, class SrcAccessor,class SrcShape,
67  class DestIterator, class DestAccessor,
68  class Neighborhood3D>
69  unsigned int labelVolume(SrcIterator s_Iter, SrcShape srcShape, SrcAccessor sa,
70  DestIterator d_Iter, DestAccessor da,
71  Neighborhood3D neighborhood3D);
72 
73  template <class SrcIterator, class SrcAccessor,class SrcShape,
74  class DestIterator, class DestAccessor,
75  class Neighborhood3D, class EqualityFunctor>
76  unsigned int labelVolume(SrcIterator s_Iter, SrcShape srcShape, SrcAccessor sa,
77  DestIterator d_Iter, DestAccessor da,
78  Neighborhood3D neighborhood3D, EqualityFunctor equal);
79 
80  }
81  \endcode
82 
83  use argument objects in conjunction with \ref ArgumentObjectFactories :
84  \code
85  namespace vigra {
86 
87  template <class SrcIterator, class SrcAccessor,class SrcShape,
88  class DestIterator, class DestAccessor,
89  class Neighborhood3D>
90  unsigned int labelVolume(triple<SrcIterator, SrcShape, SrcAccessor> src,
91  pair<DestIterator, DestAccessor> dest,
92  Neighborhood3D neighborhood3D);
93 
94  template <class SrcIterator, class SrcAccessor,class SrcShape,
95  class DestIterator, class DestAccessor,
96  class Neighborhood3D, class EqualityFunctor>
97  unsigned int labelVolume(triple<SrcIterator, SrcShape, SrcAccessor> src,
98  pair<DestIterator, DestAccessor> dest,
99  Neighborhood3D neighborhood3D, EqualityFunctor equal);
100 
101  }
102  \endcode
103 
104  use with 3D-Six-Neighborhood:
105  \code
106  namespace vigra {
107 
108  template <class SrcIterator, class SrcAccessor,class SrcShape,
109  class DestIterator, class DestAccessor>
110  unsigned int labelVolumeSix(triple<SrcIterator, SrcShape, SrcAccessor> src,
111  pair<DestIterator, DestAccessor> dest);
112 
113  }
114  \endcode
115 
116  Connected components are defined as regions with uniform voxel
117  values. Thus, <TT>SrcAccessor::value_type</TT> either must be
118  equality comparable (first form), or an EqualityFunctor must be
119  provided that realizes the desired predicate (second form). The
120  destination's value type should be large enough to hold the labels
121  without overflow. Region numbers will be a consecutive sequence
122  starting with one and ending with the region number returned by
123  the function (inclusive).
124 
125  Return: the number of regions found (= largest region label)
126 
127  <b> Usage:</b>
128 
129  <b>\#include</b> <vigra/labelvolume.hxx><br>
130  Namespace: vigra
131 
132  \code
133  typedef vigra::MultiArray<3,int> IntVolume;
134  IntVolume src(IntVolume::difference_type(w,h,d));
135  IntVolume dest(IntVolume::difference_type(w,h,d));
136 
137  // find 6-connected regions
138  int max_region_label = vigra::labelVolumeSix(srcMultiArrayRange(src), destMultiArray(dest));
139 
140  // find 26-connected regions
141  int max_region_label = vigra::labelVolume(srcMultiArrayRange(src), destMultiArray(dest), NeighborCode3DTwentySix());
142  \endcode
143 
144  <b> Required Interface:</b>
145 
146  \code
147  SrcIterator src_begin;
148  SrcShape shape;
149  DestIterator dest_begin;
150 
151  SrcAccessor src_accessor;
152  DestAccessor dest_accessor;
153 
154  SrcAccessor::value_type u = src_accessor(src_begin);
155 
156  u == u // first form
157 
158  EqualityFunctor equal; // second form
159  equal(u, u) // second form
160 
161  int i;
162  dest_accessor.set(i, dest_begin);
163  \endcode
164 
165 */
166 doxygen_overloaded_function(template <...> unsigned int labelVolume)
167 
168 
169 template <class SrcIterator, class SrcAccessor,class SrcShape,
170  class DestIterator, class DestAccessor,
171  class Neighborhood3D>
172 unsigned int labelVolume(SrcIterator s_Iter, SrcShape srcShape, SrcAccessor sa,
173  DestIterator d_Iter, DestAccessor da,
174  Neighborhood3D neighborhood3D)
175 {
176  return labelVolume(s_Iter, srcShape, sa, d_Iter, da, neighborhood3D, std::equal_to<typename SrcAccessor::value_type>());
177 }
178 
179 template <class SrcIterator, class SrcAccessor,class SrcShape,
180  class DestIterator, class DestAccessor,
181  class Neighborhood3D>
182 unsigned int labelVolume(triple<SrcIterator, SrcShape, SrcAccessor> src,
183  pair<DestIterator, DestAccessor> dest,
184  Neighborhood3D neighborhood3D)
185 {
186  return labelVolume(src.first, src.second, src.third, dest.first, dest.second, neighborhood3D, std::equal_to<typename SrcAccessor::value_type>());
187 }
188 
189 template <class SrcIterator, class SrcAccessor,class SrcShape,
190  class DestIterator, class DestAccessor,
191  class Neighborhood3D, class EqualityFunctor>
192 unsigned int labelVolume(triple<SrcIterator, SrcShape, SrcAccessor> src,
193  pair<DestIterator, DestAccessor> dest,
194  Neighborhood3D neighborhood3D, EqualityFunctor equal)
195 {
196  return labelVolume(src.first, src.second, src.third, dest.first, dest.second, neighborhood3D, equal);
197 }
198 
199 template <class SrcIterator, class SrcAccessor,class SrcShape,
200  class DestIterator, class DestAccessor,
201  class Neighborhood3D, class EqualityFunctor>
202 unsigned int labelVolume(SrcIterator s_Iter, SrcShape srcShape, SrcAccessor sa,
203  DestIterator d_Iter, DestAccessor da,
204  Neighborhood3D, EqualityFunctor equal)
205 {
206  typedef typename DestAccessor::value_type LabelType;
207 
208  //basically needed for iteration and border-checks
209  int w = srcShape[0], h = srcShape[1], d = srcShape[2];
210  int x,y,z;
211 
212  // temporary image to store region labels
213  detail::UnionFindArray<LabelType> label;
214 
215  //Declare traversers for all three dims at target
216  SrcIterator zs = s_Iter;
217  DestIterator zd = d_Iter;
218 
219  // initialize the neighborhood traversers
220  NeighborOffsetCirculator<Neighborhood3D> nce(Neighborhood3D::CausalLast);
221  ++nce;
222  // pass 1: scan image from upper left front to lower right back
223  // to find connected components
224 
225  // Each component will be represented by a tree of pixels. Each
226  // pixel contains the scan order address of its parent in the
227  // tree. In order for pass 2 to work correctly, the parent must
228  // always have a smaller scan order address than the child.
229  // Therefore, we can merge trees only at their roots, because the
230  // root of the combined tree must have the smallest scan order
231  // address among all the tree's pixels/ nodes. The root of each
232  // tree is distinguished by pointing to itself (it contains its
233  // own scan order address). This condition is enforced whenever a
234  // new region is found or two regions are merged
235  for(z = 0; z != d; ++z, ++zs.dim2(), ++zd.dim2())
236  {
237  SrcIterator ys(zs);
238  DestIterator yd(zd);
239 
240  for(y = 0; y != h; ++y, ++ys.dim1(), ++yd.dim1())
241  {
242  SrcIterator xs(ys);
243  DestIterator xd(yd);
244 
245  for(x = 0; x != w; ++x, ++xs.dim0(), ++xd.dim0())
246  {
247  LabelType currentLabel = label.nextFreeLabel();
248 
249  //check whether there is a special border treatment to be used or not
250  AtVolumeBorder atBorder = isAtVolumeBorderCausal(x,y,z,w,h,d);
251 
252  //We are not at the border!
253  if(atBorder == NotAtBorder)
254  {
255  NeighborOffsetCirculator<Neighborhood3D> nc(Neighborhood3D::CausalFirst);
256 
257  do
258  {
259  // if colors are equal
260  if(equal(sa(xs), sa(xs, *nc)))
261  {
262  currentLabel = label.makeUnion(label[da(xd,*nc)], currentLabel);
263  }
264  ++nc;
265  }
266  while(nc!=nce);
267  }
268  else //we are at a border - handle this!!
269  {
270  NeighborOffsetCirculator<Neighborhood3D> nc(Neighborhood3D::nearBorderDirectionsCausal(atBorder,0));
271  int j=0;
272  while(nc.direction() != Neighborhood3D::Error)
273  {
274  /*
275  SrcShape s(x,y,z), sn = s + *nc;
276 
277  if (sn[0]<0 || sn[0]>=w || sn[1]<0 || sn[1]>=h || sn[2]<0 || sn[2]>=d)
278  {
279  std::cerr << "coordinate error at " << s << ", offset " << *nc << ", index " << (nc).direction() << " at border " <<
280  atBorder << std::endl;
281 
282  }
283  */
284  // colors equal???
285  if(equal(sa(xs), sa(xs, *nc)))
286  {
287  currentLabel = label.makeUnion(label[da(xd,*nc)], currentLabel);
288  }
289  nc.turnTo(Neighborhood3D::nearBorderDirectionsCausal(atBorder,++j));
290  }
291  }
292  da.set(label.finalizeLabel(currentLabel), xd);
293  }
294  }
295  }
296 
297  LabelType count = label.makeContiguous();
298 
299  // pass 2: assign one label to each region (tree)
300  // so that labels form a consecutive sequence 1, 2, ...
301  zd = d_Iter;
302  for(z=0; z != d; ++z, ++zd.dim2())
303  {
304  DestIterator yd(zd);
305 
306  for(y=0; y != h; ++y, ++yd.dim1())
307  {
308  DestIterator xd(yd);
309 
310  for(x = 0; x != w; ++x, ++xd.dim0())
311  {
312  da.set(label[da(xd)], xd);
313  }
314  }
315  }
316  return count;
317 }
318 
319 /********************************************************/
320 /* */
321 /* labelVolumeSix */
322 /* */
323 /********************************************************/
324 
325 /** \brief Find the connected components of a segmented volume
326  using the 6-neighborhood.
327 
328  See \ref labelVolume() for detailed documentation.
329 
330 */
331 template <class SrcIterator, class SrcAccessor,class SrcShape,
332  class DestIterator, class DestAccessor>
333 unsigned int labelVolumeSix(triple<SrcIterator, SrcShape, SrcAccessor> src,
334  pair<DestIterator, DestAccessor> dest)
335 {
336  return labelVolume(src.first, src.second, src.third, dest.first, dest.second, NeighborCode3DSix(), std::equal_to<typename SrcAccessor::value_type>());
337 }
338 
339 
340 
341 
342 /********************************************************/
343 /* */
344 /* labelVolumeWithBackground */
345 /* */
346 /********************************************************/
347 
348 /** \brief Find the connected components of a segmented volume,
349  excluding the background from labeling.
350 
351  <b> Declarations:</b>
352 
353  pass arguments explicitly:
354  \code
355  namespace vigra {
356 
357  template <class SrcIterator, class SrcAccessor,class SrcShape,
358  class DestIterator, class DestAccessor,
359  class Neighborhood3D, class ValueType>
360  unsigned int labelVolumeWithBackground( SrcIterator s_Iter, SrcShape srcShape, SrcAccessor sa,
361  DestIterator d_Iter, DestAccessor da,
362  Neighborhood3D neighborhood3D, ValueType background_value);
363 
364  template <class SrcIterator, class SrcAccessor,class SrcShape,
365  class DestIterator, class DestAccessor,
366  class Neighborhood3D, class ValueType, class EqualityFunctor>
367  unsigned int labelVolumeWithBackground( SrcIterator s_Iter, SrcShape srcShape, SrcAccessor sa,
368  DestIterator d_Iter, DestAccessor da,
369  Neighborhood3D neighborhood3D, ValueType background_value,
370  EqualityFunctor equal);
371 
372  }
373  \endcode
374 
375  use argument objects in conjunction with \ref ArgumentObjectFactories :
376  \code
377  namespace vigra {
378 
379  template <class SrcIterator, class SrcAccessor,class SrcShape,
380  class DestIterator, class DestAccessor,
381  class Neighborhood3D, class ValueType>
382  unsigned int labelVolumeWithBackground( triple<SrcIterator, SrcShape, SrcAccessor> src,
383  pair<DestIterator, DestAccessor> dest,
384  Neighborhood3D neighborhood3D, ValueType background_value);
385 
386  template <class SrcIterator, class SrcAccessor,class SrcShape,
387  class DestIterator, class DestAccessor,
388  class Neighborhood3D, class ValueType, class EqualityFunctor>
389  unsigned int labelVolumeWithBackground( triple<SrcIterator, SrcShape, SrcAccessor> src,
390  pair<DestIterator, DestAccessor> dest,
391  Neighborhood3D neighborhood3D, ValueType background_value,
392  EqualityFunctor equal);
393 
394  }
395  \endcode
396 
397  Connected components are defined as regions with uniform voxel
398  values. Thus, <TT>SrcAccessor::value_type</TT> either must be
399  equality comparable (first form), or an EqualityFunctor must be
400  provided that realizes the desired predicate (second form). All
401  voxel equal to the given '<TT>background_value</TT>' are ignored
402  when determining connected components and remain untouched in the
403  destination volume.
404 
405  The destination's value type should be large enough to hold the
406  labels without overflow. Region numbers will be a consecutive
407  sequence starting with one and ending with the region number
408  returned by the function (inclusive).
409 
410  Return: the number of regions found (= largest region label)
411 
412  <b> Usage:</b>
413 
414  <b>\#include</b> <vigra/labelvolume.hxx><br>
415  Namespace: vigra
416 
417  \code
418  typedef vigra::MultiArray<3,int> IntVolume;
419  IntVolume src(IntVolume::difference_type(w,h,d));
420  IntVolume dest(IntVolume::difference_type(w,h,d));
421 
422  // find 6-connected regions
423  int max_region_label = vigra::labelVolumeWithBackground(
424  srcMultiArrayRange(src), destMultiArray(dest), NeighborCode3DSix(), 0);
425  \endcode
426 
427  <b> Required Interface:</b>
428 
429  \code
430  SrcIterator src_begin;
431  SrcShape shape;
432  DestIterator dest_begin;
433 
434  SrcAccessor src_accessor;
435  DestAccessor dest_accessor;
436 
437  SrcAccessor::value_type u = src_accessor(src_begin);
438 
439  u == u // first form
440 
441  EqualityFunctor equal; // second form
442  equal(u, u) // second form
443 
444  int i;
445  dest_accessor.set(i, dest_begin);
446  \endcode
447 
448 */
449 doxygen_overloaded_function(template <...> unsigned int labelVolumeWithBackground)
450 
451 template <class SrcIterator, class SrcAccessor,class SrcShape,
452  class DestIterator, class DestAccessor,
453  class Neighborhood3D,
454  class ValueType>
455 unsigned int labelVolumeWithBackground(SrcIterator s_Iter, SrcShape srcShape, SrcAccessor sa,
456  DestIterator d_Iter, DestAccessor da,
457  Neighborhood3D neighborhood3D, ValueType backgroundValue)
458 {
459  return labelVolumeWithBackground(s_Iter, srcShape, sa, d_Iter, da, neighborhood3D, backgroundValue, std::equal_to<typename SrcAccessor::value_type>());
460 }
461 
462 template <class SrcIterator, class SrcAccessor,class SrcShape,
463  class DestIterator, class DestAccessor,
464  class Neighborhood3D,
465  class ValueType>
466 unsigned int labelVolumeWithBackground(triple<SrcIterator, SrcShape, SrcAccessor> src,
467  pair<DestIterator, DestAccessor> dest,
468  Neighborhood3D neighborhood3D, ValueType backgroundValue)
469 {
470  return labelVolumeWithBackground(src.first, src.second, src.third, dest.first, dest.second, neighborhood3D, backgroundValue, std::equal_to<typename SrcAccessor::value_type>());
471 }
472 
473 template <class SrcIterator, class SrcAccessor,class SrcShape,
474  class DestIterator, class DestAccessor,
475  class Neighborhood3D,
476  class ValueType, class EqualityFunctor>
477 unsigned int labelVolumeWithBackground(triple<SrcIterator, SrcShape, SrcAccessor> src,
478  pair<DestIterator, DestAccessor> dest,
479  Neighborhood3D neighborhood3D, ValueType backgroundValue, EqualityFunctor equal)
480 {
481  return labelVolumeWithBackground(src.first, src.second, src.third, dest.first, dest.second, neighborhood3D, backgroundValue, equal);
482 }
483 
484 template <class SrcIterator, class SrcAccessor,class SrcShape,
485  class DestIterator, class DestAccessor,
486  class Neighborhood3D,
487  class ValueType, class EqualityFunctor>
488 unsigned int labelVolumeWithBackground(SrcIterator s_Iter, SrcShape srcShape, SrcAccessor sa,
489  DestIterator d_Iter, DestAccessor da,
490  Neighborhood3D,
491  ValueType backgroundValue, EqualityFunctor equal)
492 {
493  typedef typename DestAccessor::value_type LabelType;
494 
495  //basically needed for iteration and border-checks
496  int w = srcShape[0], h = srcShape[1], d = srcShape[2];
497  int x,y,z;
498 
499  // temporary image to store region labels
500  detail::UnionFindArray<LabelType> label;
501 
502  //Declare traversers for all three dims at target
503  SrcIterator zs = s_Iter;
504  DestIterator zd = d_Iter;
505 
506  // initialize the neighborhood traversers
507  NeighborOffsetCirculator<Neighborhood3D> nce(Neighborhood3D::CausalLast);
508  ++nce;
509  // pass 1: scan image from upper left front to lower right back
510  // to find connected components
511 
512  // Each component will be represented by a tree of pixels. Each
513  // pixel contains the scan order address of its parent in the
514  // tree. In order for pass 2 to work correctly, the parent must
515  // always have a smaller scan order address than the child.
516  // Therefore, we can merge trees only at their roots, because the
517  // root of the combined tree must have the smallest scan order
518  // address among all the tree's pixels/ nodes. The root of each
519  // tree is distinguished by pointing to itself (it contains its
520  // own scan order address). This condition is enforced whenever a
521  // new region is found or two regions are merged
522  for(z = 0; z != d; ++z, ++zs.dim2(), ++zd.dim2())
523  {
524  SrcIterator ys(zs);
525  DestIterator yd(zd);
526 
527  for(y = 0; y != h; ++y, ++ys.dim1(), ++yd.dim1())
528  {
529  SrcIterator xs(ys);
530  DestIterator xd(yd);
531 
532  for(x = 0; x != w; ++x, ++xs.dim0(), ++xd.dim0())
533  {
534  if(equal(sa(xs), backgroundValue))
535  {
536  da.set(label[0], xd);
537  continue;
538  }
539 
540  LabelType currentLabel = label.nextFreeLabel();
541 
542  //check whether there is a special border treatment to be used or not
543  AtVolumeBorder atBorder = isAtVolumeBorderCausal(x,y,z,w,h,d);
544 
545  //We are not at the border!
546  if(atBorder == NotAtBorder)
547  {
548  NeighborOffsetCirculator<Neighborhood3D> nc(Neighborhood3D::CausalFirst);
549 
550  do
551  {
552  // if colors are equal
553  if(equal(sa(xs), sa(xs, *nc)))
554  {
555  currentLabel = label.makeUnion(label[da(xd,*nc)], currentLabel);
556  }
557  ++nc;
558  }
559  while(nc!=nce);
560  }
561  else //we are at a border - handle this!!
562  {
563  NeighborOffsetCirculator<Neighborhood3D> nc(Neighborhood3D::nearBorderDirectionsCausal(atBorder,0));
564  int j=0;
565  while(nc.direction() != Neighborhood3D::Error)
566  {
567  // colors equal???
568  if(equal(sa(xs), sa(xs, *nc)))
569  {
570  currentLabel = label.makeUnion(label[da(xd,*nc)], currentLabel);
571  }
572  nc.turnTo(Neighborhood3D::nearBorderDirectionsCausal(atBorder,++j));
573  }
574  }
575  da.set(label.finalizeLabel(currentLabel), xd);
576  }
577  }
578  }
579 
580  LabelType count = label.makeContiguous();
581 
582  // pass 2: assign one label to each region (tree)
583  // so that labels form a consecutive sequence 1, 2, ...
584  zd = d_Iter;
585  for(z=0; z != d; ++z, ++zd.dim2())
586  {
587  DestIterator yd(zd);
588 
589  for(y=0; y != h; ++y, ++yd.dim1())
590  {
591  DestIterator xd(yd);
592 
593  for(x = 0; x != w; ++x, ++xd.dim0())
594  {
595  da.set(label[da(xd)], xd);
596  }
597  }
598  }
599  return count;
600 }
601 
602 //@}
603 
604 } //end of namespace vigra
605 
606 #endif //VIGRA_LABELVOLUME_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.9.0 (Tue Oct 22 2013)