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

blockwise_labeling.hxx |

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 */
66 class BlockwiseLabelOptions
67 : public LabelOptions
68 , public BlockwiseOptions
69 {
70 public:
71 typedef BlockwiseOptions::Shape Shape;
72
73 // reimplement setter functions to allow chaining
74
75 template <class T>
76 BlockwiseLabelOptions& ignoreBackgroundValue(const T& value)
77 {
78 LabelOptions::ignoreBackgroundValue(value);
79 return *this;
80 }
81
82 BlockwiseLabelOptions & neighborhood(NeighborhoodType n)
83 {
84 LabelOptions::neighborhood(n);
85 return *this;
86 }
87
88 BlockwiseLabelOptions & blockShape(const Shape & shape)
89 {
90 BlockwiseOptions::blockShape(shape);
91 return *this;
92 }
93
94 template <class T, int N>
95 BlockwiseLabelOptions & blockShape(const TinyVector<T, N> & shape)
96 {
97 BlockwiseOptions::blockShape(shape);
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
208 typedef GridGraph<Dimensions, undirected_tag> Graph;
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>
397 Label labelMultiArrayBlockwise(const MultiArrayView<N, Data, S1>& data,
398 MultiArrayView<N, Label, S2> labels,
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>
418 Label labelMultiArrayBlockwise(const MultiArrayView<N, Data, S1>& data,
419 MultiArrayView<N, Label, S2> labels,
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>
442 Label labelMultiArrayBlockwise(const MultiArrayView<N, Data, S1>& data,
443 MultiArrayView<N, Label, S2> labels,
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>
451 Label labelMultiArrayBlockwise(const ChunkedArray<N, Data>& data,
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>
477 Label labelMultiArrayBlockwise(const ChunkedArray<N, Data>& data,
478 ChunkedArray<N, Label>& labels,
479 const BlockwiseLabelOptions& options,
480 Equal equal)
481 {
482 using namespace blockwise_labeling_detail;
483 MultiArray<N, std::vector<Label> > mapping(data.chunkArrayShape());
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>
491 Label labelMultiArrayBlockwise(const ChunkedArray<N, Data>& data,
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
vigra::ChunkedArray::chunk_end
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
vigra::GridGraph
Define a grid graph in arbitrary dimensions.
Definition: multi_fwd.hxx:217
vigra::StridedScanOrderIterator
Sequential iterator for MultiArrayView.
Definition: multi_fwd.hxx:161
vigra::LabelOptions
Option object for labelMultiArray().
Definition: multi_labeling.hxx:309
vigra::BlockwiseOptions::getBlockShapeN
TinyVector< MultiArrayIndex, N > getBlockShapeN() const
Definition: multi_blockwise.hxx:89
vigra::ChunkedArray::chunk_begin
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
vigra::LabelOptions::getNeighborhood
NeighborhoodType getNeighborhood() const
Query the neighborhood type (direct or indirect).
Definition: multi_labeling.hxx:334
vigra::MultiArray
Main MultiArray class containing the memory management.
Definition: multi_array.hxx:2422
vigra::ChunkedArray
Interface and base class for chunked arrays.
Definition: multi_array_chunked.hxx:470
vigra
Definition: accessor.hxx:43
vigra::LabelOptions::hasBackgroundValue
bool hasBackgroundValue() const
Check if some background value is to be ignored.
Definition: multi_labeling.hxx:359
vigra::doxygen_overloaded_function
doxygen_overloaded_function(template<... > void separableConvolveBlockwise) template< unsigned int N
Separated convolution on ChunkedArrays.
vigra::ArrayVector< MultiArrayIndex >
vigra::ChunkedArray::chunkArrayShape
virtual shape_type chunkArrayShape() const
Number of chunks along each coordinate direction.
Definition: multi_array_chunked.hxx:1691
vigra::BlockwiseLabelOptions
Definition: blockwise_labeling.hxx:66
vigra::BlockwiseOptions
Definition: multi_blockwise.hxx:54
vigra::labelMultiArrayBlockwise
unsigned int labelMultiArrayBlockwise(...)
Connected components labeling for MultiArrays and ChunkedArrays.
vigra::ChunkedArray::shape
shape_type const & shape() const
Return the shape in this array.
vigra::BlockwiseOptions::getBlockShape
Shape const & getBlockShape() const
Definition: multi_blockwise.hxx:71
vigra::parallel_foreach
void parallel_foreach(...)
Apply a functor to all items in a range in parallel.
vigra::TinyVector< T, N >
vigra::LabelOptions::neighborhood
LabelOptions & neighborhood(NeighborhoodType n)
Choose direct or indirect neighborhood.
Definition: multi_labeling.hxx:326
vigra::MultiArrayView
Base class for, and view to, vigra::MultiArray.
Definition: multi_array.hxx:652
vigra::LabelOptions::ignoreBackgroundValue
LabelOptions & ignoreBackgroundValue(T const &t)
Set background value.
Definition: multi_labeling.hxx:351
vigra::ArrayVectorView::size
size_type size() const
Definition: array_vector.hxx:358
vigra::GridGraph::u
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
vigra::BlockwiseOptions::blockShape
BlockwiseOptions & blockShape(const Shape &blockShape)
Definition: multi_blockwise.hxx:114
vigra::labelMultiArray
unsigned int labelMultiArray(...)
Find the connected components of a MultiArray with arbitrary many dimensions.
vigra::NeighborhoodType
NeighborhoodType
Choose the neighborhood system in a dimension-independent way.
Definition: multi_fwd.hxx:186
vigra::ParallelOptions::getNumThreads
int getNumThreads() const
Get desired number of threads.
Definition: threadpool.hxx:85

© 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