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

blockwise_labeling.hxx VIGRA

1 /************************************************************************/
2 /* */
3 /* Copyright 2013-2014 by Martin Bidlingmaier and 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_BLOCKWISE_LABELING_HXX
37 #define VIGRA_BLOCKWISE_LABELING_HXX
38 
39 #include <algorithm>
40 
41 #include "threadpool.hxx"
42 #include "counting_iterator.hxx"
43 #include "multi_gridgraph.hxx"
44 #include "multi_labeling.hxx"
45 #include "multi_blockwise.hxx"
46 #include "union_find.hxx"
47 #include "multi_array_chunked.hxx"
48 #include "metaprogramming.hxx"
49 
50 #include "visit_border.hxx"
51 #include "blockify.hxx"
52 
53 namespace vigra
54 {
55 
56 /** \addtogroup Labeling
57 */
58 //@{
59 
60  /** Options object for labelMultiArrayBlockwise().
61 
62  It is simply a subclass of both \ref vigra::LabelOptions
63  and \ref vigra::BlockwiseOptions. See there for
64  detailed documentation.
65  */
67 : public LabelOptions
68 , public BlockwiseOptions
69 {
70 public:
72 
73  // reimplement setter functions to allow chaining
74 
75  template <class T>
76  BlockwiseLabelOptions& ignoreBackgroundValue(const T& value)
77  {
79  return *this;
80  }
81 
82  BlockwiseLabelOptions & neighborhood(NeighborhoodType n)
83  {
85  return *this;
86  }
87 
88  BlockwiseLabelOptions & blockShape(const Shape & shape)
89  {
91  return *this;
92  }
93 
94  template <class T, int N>
95  BlockwiseLabelOptions & blockShape(const TinyVector<T, N> & shape)
96  {
98  return *this;
99  }
100 
101  BlockwiseLabelOptions & numThreads(const int n)
102  {
103  BlockwiseOptions::numThreads(n);
104  return *this;
105  }
106 };
107 
108 namespace blockwise_labeling_detail
109 {
110 
111 template <class Equal, class Label>
112 struct BorderVisitor
113 {
114  Label u_label_offset;
115  Label v_label_offset;
116  UnionFindArray<Label>* global_unions;
117  Equal* equal;
118 
119  template <class Data, class Shape>
120  void operator()(const Data& u_data, Label& u_label, const Data& v_data, Label& v_label, const Shape& diff)
121  {
122  if(labeling_equality::callEqual(*equal, u_data, v_data, diff))
123  {
124  global_unions->makeUnion(u_label + u_label_offset, v_label + v_label_offset);
125  }
126  }
127 };
128 
129 // needed by MSVC
130 template <class LabelBlocksIterator>
131 struct BlockwiseLabelingResult
132 {
133  typedef typename LabelBlocksIterator::value_type::value_type type;
134 };
135 
136 template <class DataBlocksIterator, class LabelBlocksIterator,
137  class Equal, class Mapping>
138 typename BlockwiseLabelingResult<LabelBlocksIterator>::type
139 blockwiseLabeling(DataBlocksIterator data_blocks_begin, DataBlocksIterator data_blocks_end,
140  LabelBlocksIterator label_blocks_begin, LabelBlocksIterator label_blocks_end,
141  BlockwiseLabelOptions const & options,
142  Equal equal,
143  Mapping& mapping)
144 {
145  typedef typename LabelBlocksIterator::value_type::value_type Label;
146  typedef typename DataBlocksIterator::shape_type Shape;
147 
148  Shape blocks_shape = data_blocks_begin.shape();
149  vigra_precondition(blocks_shape == label_blocks_begin.shape() &&
150  blocks_shape == mapping.shape(),
151  "shapes of blocks of blocks do not match");
152 
153  static const unsigned int Dimensions = DataBlocksIterator::dimension + 1;
154  MultiArray<Dimensions, Label> label_offsets(label_blocks_begin.shape());
155 
156  bool has_background = options.hasBackgroundValue();
157 
158  // mapping stage: label each block and save number of labels assigned in blocks before the current block in label_offsets
159  Label unmerged_label_number;
160  {
161  DataBlocksIterator data_blocks_it = data_blocks_begin;
162  LabelBlocksIterator label_blocks_it = label_blocks_begin;
163  typename MultiArray<Dimensions, Label>::iterator offsets_it = label_offsets.begin();
164  Label current_offset = 0;
165  // a la OPENMP_PRAGMA FOR
166 
167  auto d = std::distance(data_blocks_begin, data_blocks_end);
168 
169 
170  std::vector<Label> nSeg(d);
171  //std::vector<int> ids(d);
172  //std::iota(ids.begin(), ids.end(), 0 );
173 
174  parallel_foreach(options.getNumThreads(), d,
175  [&](const int threadId, const uint64_t i){
176  Label resVal = labelMultiArray(data_blocks_it[i], label_blocks_it[i],
177  options, equal);
178  if(has_background) // FIXME: reversed condition?
179  ++resVal;
180  nSeg[i] = resVal;
181  }
182  );
183 
184  for(int i=0; i<d;++i){
185  offsets_it[i] = current_offset;
186  current_offset+=nSeg[i];
187  }
188 
189 
190  unmerged_label_number = current_offset;
191  if(!has_background)
192  ++unmerged_label_number;
193  }
194 
195  // reduce stage: merge adjacent labels if the region overlaps
196  UnionFindArray<Label> global_unions(unmerged_label_number);
197  if(has_background)
198  {
199  // merge all labels that refer to background
200  for(typename MultiArray<Dimensions, Label>::iterator offsets_it = label_offsets.begin();
201  offsets_it != label_offsets.end();
202  ++offsets_it)
203  {
204  global_unions.makeUnion(0, *offsets_it);
205  }
206  }
207 
209  typedef typename Graph::edge_iterator EdgeIterator;
210  Graph blocks_graph(blocks_shape, options.getNeighborhood());
211  for(EdgeIterator it = blocks_graph.get_edge_iterator(); it != blocks_graph.get_edge_end_iterator(); ++it)
212  {
213  Shape u = blocks_graph.u(*it);
214  Shape v = blocks_graph.v(*it);
215  Shape difference = v - u;
216 
217  BorderVisitor<Equal, Label> border_visitor;
218  border_visitor.u_label_offset = label_offsets[u];
219  border_visitor.v_label_offset = label_offsets[v];
220  border_visitor.global_unions = &global_unions;
221  border_visitor.equal = &equal;
222  visitBorder(data_blocks_begin[u], label_blocks_begin[u],
223  data_blocks_begin[v], label_blocks_begin[v],
224  difference, options.getNeighborhood(), border_visitor);
225  }
226 
227  // fill mapping (local labels) -> (global labels)
228  Label last_label = global_unions.makeContiguous();
229  {
230  typename MultiArray<Dimensions, Label>::iterator offsets_it = label_offsets.begin();
231  Label offset = *offsets_it;
232  ++offsets_it;
233  typename Mapping::iterator mapping_it = mapping.begin();
234  for( ; offsets_it != label_offsets.end(); ++offsets_it, ++mapping_it)
235  {
236  mapping_it->clear();
237  Label next_offset = *offsets_it;
238  if(has_background)
239  {
240  for(Label current_label = offset; current_label != next_offset; ++current_label)
241  {
242  mapping_it->push_back(global_unions.findLabel(current_label));
243  }
244  }
245  else
246  {
247  mapping_it->push_back(0); // local labels start at 1
248  for(Label current_label = offset + 1; current_label != next_offset + 1; ++current_label)
249  {
250  mapping_it->push_back(global_unions.findLabel(current_label));
251  }
252  }
253 
254  offset = next_offset;
255  }
256  // last block:
257  // instead of next_offset, use last_label+1
258  mapping_it->clear();
259  if(has_background)
260  {
261  for(Label current_label = offset; current_label != unmerged_label_number; ++current_label)
262  {
263  mapping_it->push_back(global_unions.findLabel(current_label));
264  }
265  }
266  else
267  {
268  mapping_it->push_back(0);
269  for(Label current_label = offset + 1; current_label != unmerged_label_number; ++current_label)
270  {
271  mapping_it->push_back(global_unions.findLabel(current_label));
272  }
273  }
274  }
275  return last_label;
276 }
277 
278 
279 template <class LabelBlocksIterator, class MappingIterator>
280 void toGlobalLabels(LabelBlocksIterator label_blocks_begin, LabelBlocksIterator label_blocks_end,
281  MappingIterator mapping_begin, MappingIterator mapping_end)
282 {
283  typedef typename LabelBlocksIterator::value_type LabelBlock;
284  for( ; label_blocks_begin != label_blocks_end; ++label_blocks_begin, ++mapping_begin)
285  {
286  vigra_assert(mapping_begin != mapping_end, "");
287  for(typename LabelBlock::iterator labels_it = label_blocks_begin->begin();
288  labels_it != label_blocks_begin->end();
289  ++labels_it)
290  {
291  vigra_assert(*labels_it < mapping_begin->size(), "");
292  *labels_it = (*mapping_begin)[*labels_it];
293  }
294  }
295 }
296 
297 } // namespace blockwise_labeling_detail
298 
299 /*************************************************************/
300 /* */
301 /* labelMultiArrayBlockwise */
302 /* */
303 /*************************************************************/
304 
305 /** \brief Connected components labeling for MultiArrays and ChunkedArrays.
306 
307  <b> Declarations:</b>
308 
309  \code
310  namespace vigra {
311  // assign local labels and generate mapping (local labels) -> (global labels) for each chunk
312  template <unsigned int N, class Data, class S1,
313  class Label, class S2,
314  class Equal, class S3>
315  Label labelMultiArrayBlockwise(const MultiArrayView<N, Data, S1>& data,
316  MultiArrayView<N, Label, S2> labels,
317  const BlockwiseLabelOptions& options,
318  Equal equal,
319  MultiArrayView<N, std::vector<Label>, S3>& mapping);
320 
321  // assign local labels and generate mapping (local labels) -> (global labels) for each chunk
322  template <unsigned int N, class T, class S1,
323  class Label, class S2,
324  class EqualityFunctor>
325  Label labelMultiArrayBlockwise(const ChunkedArray<N, Data>& data,
326  ChunkedArray<N, Label>& labels,
327  const BlockwiseLabelOptions& options,
328  Equal equal,
329  MultiArrayView<N, std::vector<Label>, S3>& mapping);
330 
331  // assign global labels
332  template <unsigned int N, class Data, class S1,
333  class Label, class S2,
334  class Equal>
335  Label labelMultiArrayBlockwise(const MultiArrayView<N, Data, S1>& data,
336  MultiArrayView<N, Label, S2> labels,
337  const BlockwiseLabelOptions& options,
338  Equal equal);
339 
340  // assign global labels
341  template <unsigned int N, class T, class S1,
342  class Label, class S2,
343  class EqualityFunctor = std::equal_to<T> >
344  Label labelMultiArrayBlockwise(const ChunkedArray<N, Data>& data,
345  ChunkedArray<N, Label>& labels,
346  const BlockwiseLabelOptions& options = BlockwiseLabelOptions(),
347  Equal equal = std::equal_to<T>());
348  }
349  \endcode
350 
351  The resulting labeling is equivalent to a labeling by \ref labelMultiArray, that is, the connected components are the same but may have different ids.
352  \ref NeighborhoodType and background value (if any) can be specified with the LabelOptions object.
353  If the \a mapping parameter is provided, each chunk is labeled seperately and contiguously (starting at one, zero for background),
354  with \a mapping containing a mapping of local labels to global labels for each chunk.
355  Thus, the shape of 'mapping' has to be large enough to hold each chunk coordinate.
356 
357  Return: the number of regions found (=largest global region label)
358 
359  <b> Usage: </b>
360 
361  <b>\#include </b> <vigra/blockwise_labeling.hxx><br>
362  Namespace: vigra
363 
364  \code
365  Shape3 shape = Shape3(10);
366  Shape3 chunk_shape = Shape3(4);
367  ChunkedArrayLazy<3, int> data(shape, chunk_shape);
368  // fill data ...
369 
370  ChunkedArrayLazy<3, size_t> labels(shape, chunk_shape);
371 
372  MultiArray<3, std::vector<size_t> > mapping(Shape3(3)); // there are 3 chunks in each dimension
373 
374  labelMultiArrayBlockwise(data, labels, LabelOptions().neighborhood(DirectNeighborhood).background(0),
375  std::equal_to<int>(), mapping);
376 
377  // check out chunk in the middle
378  MultiArray<3, size_t> middle_chunk(Shape3(4));
379  labels.checkoutSubarray(Shape3(4), middle_chunk);
380 
381  // print number of non-background labels assigned in middle_chunk
382  cout << mapping[Shape3(1)].size() << endl;
383 
384  // get local label for voxel
385  // this may be the same value assigned to different component in another chunk
386  size_t local_label = middle_chunk[Shape3(2)];
387  // get global label for voxel
388  // if another voxel has the same label, they are in the same connected component albeit they may be in different chunks
389  size_t global_label = mapping[Shape3(1)][local_label
390  \endcode
391  */
392 doxygen_overloaded_function(template <...> unsigned int labelMultiArrayBlockwise)
393 
394 template <unsigned int N, class Data, class S1,
395  class Label, class S2,
396  class Equal, class S3>
399  const BlockwiseLabelOptions& options,
400  Equal equal,
401  MultiArrayView<N, std::vector<Label>, S3>& mapping)
402 {
403  using namespace blockwise_labeling_detail;
404 
405  typedef typename MultiArrayShape<N>::type Shape;
406  Shape block_shape(options.getBlockShapeN<N>());
407 
408  MultiArray<N, MultiArrayView<N, Data, S1> > data_blocks = blockify(data, block_shape);
409  MultiArray<N, MultiArrayView<N, Label, S2> > label_blocks = blockify(labels, block_shape);
410  return blockwiseLabeling(data_blocks.begin(), data_blocks.end(),
411  label_blocks.begin(), label_blocks.end(),
412  options, equal, mapping);
413 }
414 
415 template <unsigned int N, class Data, class S1,
416  class Label, class S2,
417  class Equal>
420  const BlockwiseLabelOptions& options,
421  Equal equal)
422 {
423  using namespace blockwise_labeling_detail;
424 
425  typedef typename MultiArrayShape<N>::type Shape;
426  Shape block_shape(options.getBlockShapeN<N>());
427 
428  MultiArray<N, MultiArrayView<N, Data, S1> > data_blocks = blockify(data, block_shape);
429  MultiArray<N, MultiArrayView<N, Label, S2> > label_blocks = blockify(labels, block_shape);
430  MultiArray<N, std::vector<Label> > mapping(data_blocks.shape());
431  Label last_label = blockwiseLabeling(data_blocks.begin(), data_blocks.end(),
432  label_blocks.begin(), label_blocks.end(),
433  options, equal, mapping);
434 
435  // replace local labels by global labels
436  toGlobalLabels(label_blocks.begin(), label_blocks.end(), mapping.begin(), mapping.end());
437  return last_label;
438 }
439 
440 template <unsigned int N, class Data, class S1,
441  class Label, class S2>
444  const BlockwiseLabelOptions& options = BlockwiseLabelOptions())
445 {
446  return labelMultiArrayBlockwise(data, labels, options, std::equal_to<Data>());
447 }
448 
449 
450 template <unsigned int N, class Data, class Label, class Equal, class S3>
452  ChunkedArray<N, Label>& labels,
453  const BlockwiseLabelOptions& options,
454  Equal equal,
455  MultiArrayView<N, std::vector<Label>, S3> mapping)
456 {
457  using namespace blockwise_labeling_detail;
458 
459  vigra_precondition(options.getBlockShape().size() == 0,
460  "labelMultiArrayBlockwise(ChunkedArray, ...): custom block shapes not supported "
461  "(always uses the array's chunk shape).");
462 
463  typedef typename ChunkedArray<N, Data>::shape_type Shape;
464 
465  typedef typename ChunkedArray<N, Data>::chunk_const_iterator DataChunkIterator;
466  typedef typename ChunkedArray<N, Label>::chunk_iterator LabelChunkIterator;
467 
468  DataChunkIterator data_chunks_begin = data.chunk_begin(Shape(0), data.shape());
469  LabelChunkIterator label_chunks_begin = labels.chunk_begin(Shape(0), labels.shape());
470 
471  return blockwiseLabeling(data_chunks_begin, data_chunks_begin.getEndIterator(),
472  label_chunks_begin, label_chunks_begin.getEndIterator(),
473  options, equal, mapping);
474 }
475 
476 template <unsigned int N, class Data, class Label, class Equal>
478  ChunkedArray<N, Label>& labels,
479  const BlockwiseLabelOptions& options,
480  Equal equal)
481 {
482  using namespace blockwise_labeling_detail;
484  Label result = labelMultiArrayBlockwise(data, labels, options, equal, mapping);
485  typedef typename ChunkedArray<N, Data>::shape_type Shape;
486  toGlobalLabels(labels.chunk_begin(Shape(0), data.shape()), labels.chunk_end(Shape(0), data.shape()), mapping.begin(), mapping.end());
487  return result;
488 }
489 
490 template <unsigned int N, class Data, class Label>
492  ChunkedArray<N, Label>& labels,
493  const BlockwiseLabelOptions& options = BlockwiseLabelOptions())
494 {
495  return labelMultiArrayBlockwise(data, labels, options, std::equal_to<Data>());
496 }
497 
498 //@}
499 
500 } // namespace vigra
501 
502 #endif // VIGRA_BLOCKWISE_LABELING_HXX
chunk_iterator chunk_end(shape_type const &start, shape_type const &stop)
Create the end iterator for iteration over all chunks intersected by the given ROI.
Definition: multi_array_chunked.hxx:2439
Define a grid graph in arbitrary dimensions.
Definition: multi_fwd.hxx:217
Sequential iterator for MultiArrayView.
Definition: multi_fwd.hxx:161
Option object for labelMultiArray().
Definition: multi_labeling.hxx:309
chunk_iterator chunk_begin(shape_type const &start, shape_type const &stop)
Create an iterator over all chunks intersected by the given ROI.
Definition: multi_array_chunked.hxx:2430
int getNumThreads() const
Get desired number of threads.
Definition: threadpool.hxx:90
Main MultiArray class containing the memory management.
Definition: multi_array.hxx:2422
shape_type const & shape() const
Return the shape in this array.
Interface and base class for chunked arrays.
Definition: multi_array_chunked.hxx:470
Definition: accessor.hxx:43
doxygen_overloaded_function(template<... > void separableConvolveBlockwise) template< unsigned int N
Separated convolution on ChunkedArrays.
bool hasBackgroundValue() const
Check if some background value is to be ignored.
Definition: multi_labeling.hxx:359
Definition: blockwise_labeling.hxx:66
TinyVector< MultiArrayIndex, N > getBlockShapeN() const
Definition: multi_blockwise.hxx:89
Definition: multi_blockwise.hxx:54
unsigned int labelMultiArrayBlockwise(...)
Connected components labeling for MultiArrays and ChunkedArrays.
Node u(Edge const &e) const
Get the start node of the given edge e (API: LEMON, the boost::graph API provides the free function ...
Definition: multi_gridgraph.hxx:2382
void parallel_foreach(...)
Apply a functor to all items in a range in parallel.
LabelOptions & neighborhood(NeighborhoodType n)
Choose direct or indirect neighborhood.
Definition: multi_labeling.hxx:326
NeighborhoodType getNeighborhood() const
Query the neighborhood type (direct or indirect).
Definition: multi_labeling.hxx:334
Shape const & getBlockShape() const
Definition: multi_blockwise.hxx:71
virtual shape_type chunkArrayShape() const
Number of chunks along each coordinate direction.
Definition: multi_array_chunked.hxx:1691
Base class for, and view to, vigra::MultiArray.
Definition: multi_array.hxx:652
LabelOptions & ignoreBackgroundValue(T const &t)
Set background value.
Definition: multi_labeling.hxx:351
size_type size() const
Definition: array_vector.hxx:358
BlockwiseOptions & blockShape(const Shape &blockShape)
Definition: multi_blockwise.hxx:114
unsigned int labelMultiArray(...)
Find the connected components of a MultiArray with arbitrary many dimensions.
NeighborhoodType
Choose the neighborhood system in a dimension-independent way.
Definition: multi_fwd.hxx:186

© 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