[ 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 "multi_gridgraph.hxx"
40 #include "multi_labeling.hxx"
41 #include "union_find.hxx"
42 #include "multi_array_chunked.hxx"
43 #include "metaprogramming.hxx"
44 
45 #include "visit_border.hxx"
46 #include "blockify.hxx"
47 
48 namespace vigra
49 {
50 
51 /** \addtogroup Labeling
52 */
53 //@{
54 
55 class LabelOptions;
56 
57 namespace blockwise_labeling_detail
58 {
59 
60 template <class Equal, class Label>
61 struct BorderVisitor
62 {
63  Label u_label_offset;
64  Label v_label_offset;
65  UnionFindArray<Label>* global_unions;
66  Equal* equal;
67 
68  template <class Data, class Shape>
69  void operator()(const Data& u_data, Label& u_label, const Data& v_data, Label& v_label, const Shape& diff)
70  {
71  if(labeling_equality::callEqual(*equal, u_data, v_data, diff))
72  {
73  global_unions->makeUnion(u_label + u_label_offset, v_label + v_label_offset);
74  }
75  }
76 };
77 
78 // needed by MSVC
79 template <class LabelBlocksIterator>
80 struct BlockwiseLabelingResult
81 {
82  typedef typename LabelBlocksIterator::value_type::value_type type;
83 };
84 
85 template <class DataBlocksIterator, class LabelBlocksIterator, class Equal, class Value, class Mapping>
86 typename BlockwiseLabelingResult<LabelBlocksIterator>::type
87 blockwiseLabeling(DataBlocksIterator data_blocks_begin, DataBlocksIterator data_blocks_end,
88  LabelBlocksIterator label_blocks_begin, LabelBlocksIterator label_blocks_end,
89  NeighborhoodType neighborhood, Equal equal,
90  const Value* background_value,
91  Mapping& mapping)
92 {
93  typedef typename LabelBlocksIterator::value_type::value_type Label;
94  typedef typename DataBlocksIterator::shape_type Shape;
95 
96  Shape blocks_shape = data_blocks_begin.shape();
97  vigra_precondition(blocks_shape == label_blocks_begin.shape() &&
98  blocks_shape == mapping.shape(),
99  "shapes of blocks of blocks do not match");
100 
101  static const unsigned int Dimensions = DataBlocksIterator::dimension + 1;
102  MultiArray<Dimensions, Label> label_offsets(label_blocks_begin.shape());
103 
104  // mapping stage: label each block and save number of labels assigned in blocks before the current block in label_offsets
105  Label unmerged_label_number;
106  {
107  DataBlocksIterator data_blocks_it = data_blocks_begin;
108  LabelBlocksIterator label_blocks_it = label_blocks_begin;
109  typename MultiArray<Dimensions, Label>::iterator offsets_it = label_offsets.begin();
110  Label current_offset = 0;
111  for( ; data_blocks_it != data_blocks_end; ++data_blocks_it, ++label_blocks_it, ++offsets_it)
112  {
113  vigra_assert(label_blocks_it != label_blocks_end && offsets_it != label_offsets.end(), "");
114  *offsets_it = current_offset;
115  if(background_value)
116  {
117  current_offset += 1 + labelMultiArrayWithBackground(*data_blocks_it, *label_blocks_it,
118  neighborhood, *background_value, equal);
119  }
120  else
121  {
122  current_offset += labelMultiArray(*data_blocks_it, *label_blocks_it,
123  neighborhood, equal);
124  }
125  }
126  unmerged_label_number = current_offset;
127  if(!background_value)
128  ++unmerged_label_number;
129  }
130 
131  // reduce stage: merge adjacent labels if the region overlaps
132  UnionFindArray<Label> global_unions(unmerged_label_number);
133  if(background_value)
134  {
135  // merge all labels that refer to background
136  for(typename MultiArray<Dimensions, Label>::iterator offsets_it = label_offsets.begin();
137  offsets_it != label_offsets.end();
138  ++offsets_it)
139  {
140  global_unions.makeUnion(0, *offsets_it);
141  }
142  }
143 
144  typedef GridGraph<Dimensions, undirected_tag> Graph;
145  typedef typename Graph::edge_iterator EdgeIterator;
146  Graph blocks_graph(blocks_shape, neighborhood);
147  for(EdgeIterator it = blocks_graph.get_edge_iterator(); it != blocks_graph.get_edge_end_iterator(); ++it)
148  {
149  Shape u = blocks_graph.u(*it);
150  Shape v = blocks_graph.v(*it);
151  Shape difference = v - u;
152 
153  BorderVisitor<Equal, Label> border_visitor;
154  border_visitor.u_label_offset = label_offsets[u];
155  border_visitor.v_label_offset = label_offsets[v];
156  border_visitor.global_unions = &global_unions;
157  border_visitor.equal = &equal;
158  visitBorder(data_blocks_begin[u], label_blocks_begin[u],
159  data_blocks_begin[v], label_blocks_begin[v],
160  difference, neighborhood, border_visitor);
161  }
162 
163  // fill mapping (local labels) -> (global labels)
164  Label last_label = global_unions.makeContiguous();
165  {
166  typename MultiArray<Dimensions, Label>::iterator offsets_it = label_offsets.begin();
167  Label offset = *offsets_it;
168  ++offsets_it;
169  typename Mapping::iterator mapping_it = mapping.begin();
170  for( ; offsets_it != label_offsets.end(); ++offsets_it, ++mapping_it)
171  {
172  mapping_it->clear();
173  Label next_offset = *offsets_it;
174  if(background_value)
175  {
176  for(Label current_label = offset; current_label != next_offset; ++current_label)
177  {
178  mapping_it->push_back(global_unions.findLabel(current_label));
179  }
180  }
181  else
182  {
183  mapping_it->push_back(0); // local labels start at 1
184  for(Label current_label = offset + 1; current_label != next_offset + 1; ++current_label)
185  {
186  mapping_it->push_back(global_unions.findLabel(current_label));
187  }
188  }
189 
190  offset = next_offset;
191  }
192  // last block:
193  // instead of next_offset, use last_label+1
194  mapping_it->clear();
195  if(background_value)
196  {
197  for(Label current_label = offset; current_label != unmerged_label_number; ++current_label)
198  {
199  mapping_it->push_back(global_unions.findLabel(current_label));
200  }
201  }
202  else
203  {
204  mapping_it->push_back(0);
205  for(Label current_label = offset + 1; current_label != unmerged_label_number; ++current_label)
206  {
207  mapping_it->push_back(global_unions.findLabel(current_label));
208  }
209  }
210  }
211  return last_label;
212 }
213 
214 
215 template <class LabelBlocksIterator, class MappingIterator>
216 void toGlobalLabels(LabelBlocksIterator label_blocks_begin, LabelBlocksIterator label_blocks_end,
217  MappingIterator mapping_begin, MappingIterator mapping_end)
218 {
219  typedef typename LabelBlocksIterator::value_type LabelBlock;
220  for( ; label_blocks_begin != label_blocks_end; ++label_blocks_begin, ++mapping_begin)
221  {
222  vigra_assert(mapping_begin != mapping_end, "");
223  for(typename LabelBlock::iterator labels_it = label_blocks_begin->begin();
224  labels_it != label_blocks_begin->end();
225  ++labels_it)
226  {
227  vigra_assert(*labels_it < mapping_begin->size(), "");
228  *labels_it = (*mapping_begin)[*labels_it];
229  }
230  }
231 }
232 
233 
234 static const MultiArrayIndex default_block_side_length = 128;
235 
236 template <class T>
237 const T* getBackground(const LabelOptions& options);
238 template <class T>
239 const T* getBlockShape(const LabelOptions& options);
240 NeighborhoodType getNeighborhood(const LabelOptions& options);
241 
242 } // namespace blockwise_labeling_detail
243 
244 class LabelOptions
245 {
246 private:
247  struct type_erasure_base
248  {
249  virtual ~type_erasure_base()
250  {}
251  };
252  template <class T>
253  struct type_erasure
254  : type_erasure_base
255  {
256  type_erasure(const T& contained_obj)
257  : obj(contained_obj)
258  {}
259  T obj;
260  };
261 
262  VIGRA_UNIQUE_PTR<type_erasure_base> background_value_;
263  VIGRA_UNIQUE_PTR<type_erasure_base> block_shape_;
264  NeighborhoodType neighborhood_;
265 public:
266  LabelOptions()
267  : neighborhood_(DirectNeighborhood)
268  {}
269 private:
270  LabelOptions(const LabelOptions&); // deleted
271  LabelOptions& operator=(const LabelOptions&); // deleted
272 public:
273  template <class T>
274  LabelOptions& background(const T& background_value)
275  {
276  vigra_precondition(background_value_.get() == 0, "background set twice");
277  background_value_ = VIGRA_UNIQUE_PTR<type_erasure_base>(new type_erasure<T>(background_value));
278  return *this;
279  }
280  template <int N>
281  LabelOptions& blockShape(const TinyVector<MultiArrayIndex, N>& block_shape)
282  {
283  vigra_precondition(block_shape_.get() == 0, "block shape set twice");
284  block_shape_ = VIGRA_UNIQUE_PTR<type_erasure_base>(new type_erasure<TinyVector<MultiArrayIndex, N> >(block_shape));
285  return *this;
286  }
287 
288  LabelOptions& neighborhood(NeighborhoodType n)
289  {
290  neighborhood_ = n;
291  return *this;
292  }
293 
294 
295  template <class T>
296  friend const T* blockwise_labeling_detail::getBackground(const LabelOptions& options);
297  template <class T>
298  friend const T* blockwise_labeling_detail::getBlockShape(const LabelOptions& options);
299  friend NeighborhoodType blockwise_labeling_detail::getNeighborhood(const LabelOptions& options);
300 };
301 
302 namespace blockwise_labeling_detail
303 {
304 
305 template <class T>
306 const T* getBackground(const LabelOptions& options)
307 {
308  if(options.background_value_.get() == 0)
309  return 0;
310 
311  LabelOptions::type_erasure<T>* background = dynamic_cast<LabelOptions::type_erasure<T>*>(options.background_value_.get());
312  vigra_precondition(background != 0, "background value type and data type do not match");
313  return &background->obj;
314 }
315 template <class T>
316 const T* getBlockShape(const LabelOptions& options)
317 {
318  if(options.block_shape_.get() == 0)
319  return 0;
320 
321  LabelOptions::type_erasure<T>* block_shape = dynamic_cast<LabelOptions::type_erasure<T>*>(options.block_shape_.get());
322  vigra_precondition(block_shape != 0, "shape type has invalid dimension");
323  return &block_shape->obj;
324 }
325 
326 inline NeighborhoodType getNeighborhood(const LabelOptions& options)
327 {
328  return options.neighborhood_;
329 }
330 
331 }
332 
333 
334 
335 template <unsigned int N, class Data, class S1,
336  class Label, class S2,
337  class Equal, class S3>
338 Label labelMultiArrayBlockwise(const MultiArrayView<N, Data, S1>& data,
339  MultiArrayView<N, Label, S2> labels,
340  const LabelOptions& options,
341  Equal equal, MultiArrayView<N, std::vector<Label>, S3>& mapping)
342 {
343  using namespace blockwise_labeling_detail;
344  TinyVector<MultiArrayIndex, N> block_shape(default_block_side_length);
345  const TinyVector<MultiArrayIndex, N>* options_block_shape = getBlockShape<TinyVector<MultiArrayIndex, N> >(options);
346  if(options_block_shape)
347  block_shape = *options_block_shape;
348  const Data* background_value = getBackground<Data>(options);
349  NeighborhoodType neighborhood = getNeighborhood(options);
350 
351  MultiArray<N, MultiArrayView<N, Data, S1> > data_blocks = blockify(data, block_shape);
352  MultiArray<N, MultiArrayView<N, Label, S2> > label_blocks = blockify(labels, block_shape);
353  return blockwiseLabeling(data_blocks.begin(), data_blocks.end(),
354  label_blocks.begin(), label_blocks.end(),
355  neighborhood, equal, background_value, mapping);
356 }
357 template <unsigned int N, class Data, class S1,
358  class Label, class S2,
359  class Equal>
360 Label labelMultiArrayBlockwise(const MultiArrayView<N, Data, S1>& data,
361  MultiArrayView<N, Label, S2> labels, const LabelOptions& options, Equal equal) {
362  using namespace blockwise_labeling_detail;
363  TinyVector<MultiArrayIndex, N> block_shape(default_block_side_length);
364  const TinyVector<MultiArrayIndex, N>* options_block_shape = getBlockShape<TinyVector<MultiArrayIndex, N> >(options);
365  if(options_block_shape)
366  block_shape = *options_block_shape;
367  const Data* background_value = getBackground<Data>(options);
368  NeighborhoodType neighborhood = getNeighborhood(options);
369 
370  MultiArray<N, MultiArrayView<N, Data, S1> > data_blocks = blockify(data, block_shape);
371  MultiArray<N, MultiArrayView<N, Label, S2> > label_blocks = blockify(labels, block_shape);
372  MultiArray<N, std::vector<Label> > mapping(data_blocks.shape());
373  Label last_label = blockwiseLabeling(data_blocks.begin(), data_blocks.end(),
374  label_blocks.begin(), label_blocks.end(),
375  neighborhood, equal, background_value, mapping);
376 
377  // replace local labels by global labels
378  toGlobalLabels(label_blocks.begin(), label_blocks.end(), mapping.begin(), mapping.end());
379  return last_label;
380 }
381 template <unsigned int N, class Data, class S1,
382  class Label, class S2>
383 Label labelMultiArrayBlockwise(const MultiArrayView<N, Data, S1>& data,
384  MultiArrayView<N, Label, S2> labels,
385  const LabelOptions& options = LabelOptions())
386 {
387  return labelMultiArrayBlockwise(data, labels, options, std::equal_to<Data>());
388 }
389 
390 
391 
392 /*************************************************************/
393 /* */
394 /* labelMultiArrayBlockwise */
395 /* */
396 /*************************************************************/
397 
398 /** \brief Connected components labeling for ChunkedArrays.
399 
400  <b> Declarations:</b>
401 
402  \code
403  namespace vigra {
404  // assign local labels and generate mapping (local labels) -> (global labels) for each chunk
405  template <unsigned int N, class T, class S1,
406  class Label, class S2,
407  class EqualityFunctor>
408  Label labelMultiArrayBlockwise(const ChunkedArray<N, Data>& data,
409  ChunkedArray<N, Label>& labels,
410  const LabelOptions& options,
411  Equal equal,
412  MultiArrayView<N, std::vector<Label>, S3>& mapping);
413  // assign global labels
414  template <unsigned int N, class T, class S1,
415  class Label, class S2,
416  class EqualityFunctor = std::equal_to<T> >
417  Label labelMultiArrayBlockwise(const ChunkedArray<N, Data>& data,
418  ChunkedArray<N, Label>& labels,
419  const LabelOptions& options = LabelOptions(),
420  Equal equal = std::equal_to<T>());
421  }
422  \endcode
423 
424  The resulting labeling is equivalent to a labeling by \ref labelMultiArray, that is, the connected components are the same but may have different ids.
425  \ref NeighborhoodType and background value (if any) can be specified with the LabelOptions object.
426  If the \a mapping parameter is provided, each chunk is labeled seperately and contiguously (starting at one, zero for background),
427  with \a mapping containing a mapping of local labels to global labels for each chunk.
428  Thus, the shape of 'mapping' has to be large enough to hold each chunk coordinate.
429 
430  Return: the number of regions found (=largest global region label)
431 
432  <b> Usage: </b>
433 
434  <b>\#include </b> <vigra/blockwise_labeling.hxx><br>
435  Namespace: vigra
436 
437  \code
438  Shape3 shape = Shape3(10);
439  Shape3 chunk_shape = Shape3(4);
440  ChunkedArrayLazy<3, int> data(shape, chunk_shape);
441  // fill data ...
442 
443  ChunkedArrayLazy<3, size_t> labels(shape, chunk_shape);
444 
445  MultiArray<3, std::vector<size_t> > mapping(Shape3(3)); // there are 3 chunks in each dimension
446 
447  labelMultiArrayBlockwise(data, labels, LabelOptions().neighborhood(DirectNeighborhood).background(0),
448  std::equal_to<int>(), mapping);
449 
450  // check out chunk in the middle
451  MultiArray<3, size_t> middle_chunk(Shape3(4));
452  labels.checkoutSubarray(Shape3(4), middle_chunk);
453 
454  // print number of non-background labels assigned in middle_chunk
455  cout << mapping[Shape3(1)].size() << endl;
456 
457  // get local label for voxel
458  // this may be the same value assigned to different component in another chunk
459  size_t local_label = middle_chunk[Shape3(2)];
460  // get global label for voxel
461  // if another voxel has the same label, they are in the same connected component albeit they may be in different chunks
462  size_t global_label = mapping[Shape3(1)][local_label
463  \endcode
464  */
465 doxygen_overloaded_function(template <...> unsigned int labelMultiArrayBlockwise)
466 
467 template <unsigned int N, class Data, class Label, class Equal, class S3>
468 Label labelMultiArrayBlockwise(const ChunkedArray<N, Data>& data,
469  ChunkedArray<N, Label>& labels,
470  const LabelOptions& options,
471  Equal equal, MultiArrayView<N, std::vector<Label>, S3> mapping)
472 {
473  using namespace blockwise_labeling_detail;
474  const TinyVector<MultiArrayIndex, N>* options_block_shape = getBlockShape<TinyVector<MultiArrayIndex, N> >(options);
475  vigra_precondition(options_block_shape == 0, "block shape not supported for chunked arrays, uses chunk size per default");
476 
477  typedef typename ChunkedArray<N, Data>::shape_type Shape;
478 
479  const Data* background_value = getBackground<Data>(options);
480  NeighborhoodType neighborhood = getNeighborhood(options);
481 
482  typedef typename ChunkedArray<N, Data>::chunk_const_iterator DataChunkIterator;
483  typedef typename ChunkedArray<N, Label>::chunk_iterator LabelChunkIterator;
484 
485  DataChunkIterator data_chunks_begin = data.chunk_begin(Shape(0), data.shape());
486  LabelChunkIterator label_chunks_begin = labels.chunk_begin(Shape(0), labels.shape());
487 
488  return blockwiseLabeling(data_chunks_begin, data_chunks_begin.getEndIterator(),
489  label_chunks_begin, label_chunks_begin.getEndIterator(),
490  neighborhood, equal, background_value, mapping);
491 }
492 template <unsigned int N, class Data, class Label, class Equal>
493 Label labelMultiArrayBlockwise(const ChunkedArray<N, Data>& data,
494  ChunkedArray<N, Label>& labels,
495  const LabelOptions& options, Equal equal)
496 {
497  using namespace blockwise_labeling_detail;
498  MultiArray<N, std::vector<Label> > mapping(data.chunkArrayShape());
499  Label result = labelMultiArrayBlockwise(data, labels, options, equal, mapping);
500  typedef typename ChunkedArray<N, Data>::shape_type Shape;
501  toGlobalLabels(labels.chunk_begin(Shape(0), data.shape()), labels.chunk_end(Shape(0), data.shape()), mapping.begin(), mapping.end());
502  return result;
503 }
504 template <unsigned int N, class Data, class Label>
505 Label labelMultiArrayBlockwise(const ChunkedArray<N, Data>& data,
506  ChunkedArray<N, Label>& labels,
507  const LabelOptions& options = LabelOptions())
508 {
509  return labelMultiArrayBlockwise(data, labels, options, std::equal_to<Data>());
510 }
511 
512 //@}
513 
514 } // namespace vigra
515 
516 #endif // VIGRA_BLOCKWISE_LABELING_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.10.0 (Thu Jan 8 2015)