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>
ignoreBackgroundValue(const T & value)76 BlockwiseLabelOptions& ignoreBackgroundValue(const T& value)
77 {
78 LabelOptions::ignoreBackgroundValue(value);
79 return *this;
80 }
81
neighborhood(NeighborhoodType n)82 BlockwiseLabelOptions & neighborhood(NeighborhoodType n)
83 {
84 LabelOptions::neighborhood(n);
85 return *this;
86 }
87
blockShape(const Shape & shape)88 BlockwiseLabelOptions & blockShape(const Shape & shape)
89 {
90 BlockwiseOptions::blockShape(shape);
91 return *this;
92 }
93
94 template <class T, int N>
blockShape(const TinyVector<T,N> & shape)95 BlockwiseLabelOptions & blockShape(const TinyVector<T, N> & shape)
96 {
97 BlockwiseOptions::blockShape(shape);
98 return *this;
99 }
100
numThreads(const int n)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>
operator ()vigra::blockwise_labeling_detail::BorderVisitor120 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
blockwiseLabeling(DataBlocksIterator data_blocks_begin,DataBlocksIterator data_blocks_end,LabelBlocksIterator label_blocks_begin,LabelBlocksIterator label_blocks_end,BlockwiseLabelOptions const & options,Equal equal,Mapping & mapping)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 vigra_precondition(std::distance(data_blocks_begin,data_blocks_end) == std::distance(label_blocks_begin,label_blocks_end),
153 "the sizes of input ranges are different");
154
155 static const unsigned int Dimensions = DataBlocksIterator::dimension + 1;
156 MultiArray<Dimensions, Label> label_offsets(label_blocks_begin.shape());
157
158 bool has_background = options.hasBackgroundValue();
159
160 // mapping stage: label each block and save number of labels assigned in blocks before the current block in label_offsets
161 Label unmerged_label_number;
162 {
163 DataBlocksIterator data_blocks_it = data_blocks_begin;
164 LabelBlocksIterator label_blocks_it = label_blocks_begin;
165 typename MultiArray<Dimensions, Label>::iterator offsets_it = label_offsets.begin();
166 Label current_offset = 0;
167 // a la OPENMP_PRAGMA FOR
168
169 auto d = std::distance(data_blocks_begin, data_blocks_end);
170
171
172 std::vector<Label> nSeg(d);
173 //std::vector<int> ids(d);
174 //std::iota(ids.begin(), ids.end(), 0 );
175
176 parallel_foreach(options.getNumThreads(), d,
177 [&](const int /*threadId*/, const uint64_t i){
178 Label resVal = labelMultiArray(data_blocks_it[i], label_blocks_it[i],
179 options, equal);
180 if(has_background) // FIXME: reversed condition?
181 ++resVal;
182 nSeg[i] = resVal;
183 }
184 );
185
186 for(int i=0; i<d;++i){
187 offsets_it[i] = current_offset;
188 current_offset+=nSeg[i];
189 }
190
191
192 unmerged_label_number = current_offset;
193 if(!has_background)
194 ++unmerged_label_number;
195 }
196
197 // reduce stage: merge adjacent labels if the region overlaps
198 UnionFindArray<Label> global_unions(unmerged_label_number);
199 if(has_background)
200 {
201 // merge all labels that refer to background
202 for(typename MultiArray<Dimensions, Label>::iterator offsets_it = label_offsets.begin();
203 offsets_it != label_offsets.end();
204 ++offsets_it)
205 {
206 global_unions.makeUnion(0, *offsets_it);
207 }
208 }
209
210 typedef GridGraph<Dimensions, undirected_tag> Graph;
211 typedef typename Graph::edge_iterator EdgeIterator;
212 Graph blocks_graph(blocks_shape, options.getNeighborhood());
213 for(EdgeIterator it = blocks_graph.get_edge_iterator(); it != blocks_graph.get_edge_end_iterator(); ++it)
214 {
215 Shape u = blocks_graph.u(*it);
216 Shape v = blocks_graph.v(*it);
217 Shape difference = v - u;
218
219 BorderVisitor<Equal, Label> border_visitor;
220 border_visitor.u_label_offset = label_offsets[u];
221 border_visitor.v_label_offset = label_offsets[v];
222 border_visitor.global_unions = &global_unions;
223 border_visitor.equal = &equal;
224 visitBorder(data_blocks_begin[u], label_blocks_begin[u],
225 data_blocks_begin[v], label_blocks_begin[v],
226 difference, options.getNeighborhood(), border_visitor);
227 }
228
229 // fill mapping (local labels) -> (global labels)
230 Label last_label = global_unions.makeContiguous();
231 {
232 typename MultiArray<Dimensions, Label>::iterator offsets_it = label_offsets.begin();
233 Label offset = *offsets_it;
234 ++offsets_it;
235 typename Mapping::iterator mapping_it = mapping.begin();
236 for( ; offsets_it != label_offsets.end(); ++offsets_it, ++mapping_it)
237 {
238 mapping_it->clear();
239 Label next_offset = *offsets_it;
240 if(has_background)
241 {
242 for(Label current_label = offset; current_label != next_offset; ++current_label)
243 {
244 mapping_it->push_back(global_unions.findLabel(current_label));
245 }
246 }
247 else
248 {
249 mapping_it->push_back(0); // local labels start at 1
250 for(Label current_label = offset + 1; current_label != next_offset + 1; ++current_label)
251 {
252 mapping_it->push_back(global_unions.findLabel(current_label));
253 }
254 }
255
256 offset = next_offset;
257 }
258 // last block:
259 // instead of next_offset, use last_label+1
260 mapping_it->clear();
261 if(has_background)
262 {
263 for(Label current_label = offset; current_label != unmerged_label_number; ++current_label)
264 {
265 mapping_it->push_back(global_unions.findLabel(current_label));
266 }
267 }
268 else
269 {
270 mapping_it->push_back(0);
271 for(Label current_label = offset + 1; current_label != unmerged_label_number; ++current_label)
272 {
273 mapping_it->push_back(global_unions.findLabel(current_label));
274 }
275 }
276 }
277 return last_label;
278 }
279
280
281 template <class LabelBlocksIterator, class MappingIterator>
toGlobalLabels(LabelBlocksIterator label_blocks_begin,LabelBlocksIterator label_blocks_end,MappingIterator mapping_begin,MappingIterator mapping_end)282 void toGlobalLabels(LabelBlocksIterator label_blocks_begin, LabelBlocksIterator label_blocks_end,
283 MappingIterator mapping_begin, MappingIterator mapping_end)
284 {
285 typedef typename LabelBlocksIterator::value_type LabelBlock;
286 for( ; label_blocks_begin != label_blocks_end; ++label_blocks_begin, ++mapping_begin)
287 {
288 vigra_assert(mapping_begin != mapping_end, "");
289 for(typename LabelBlock::iterator labels_it = label_blocks_begin->begin();
290 labels_it != label_blocks_begin->end();
291 ++labels_it)
292 {
293 vigra_assert(*labels_it < mapping_begin->size(), "");
294 *labels_it = (*mapping_begin)[*labels_it];
295 }
296 }
297 }
298
299 } // namespace blockwise_labeling_detail
300
301 /*************************************************************/
302 /* */
303 /* labelMultiArrayBlockwise */
304 /* */
305 /*************************************************************/
306
307 /** \weakgroup ParallelProcessing
308 \sa labelMultiArrayBlockwise <B>(...)</B>
309 */
310
311 /** \brief Connected components labeling for MultiArrays and ChunkedArrays.
312
313 <b> Declarations:</b>
314
315 \code
316 namespace vigra {
317 // assign local labels and generate mapping (local labels) -> (global labels) for each chunk
318 template <unsigned int N, class Data, class S1,
319 class Label, class S2,
320 class Equal, class S3>
321 Label labelMultiArrayBlockwise(const MultiArrayView<N, Data, S1>& data,
322 MultiArrayView<N, Label, S2> labels,
323 const BlockwiseLabelOptions& options,
324 Equal equal,
325 MultiArrayView<N, std::vector<Label>, S3>& mapping);
326
327 // assign local labels and generate mapping (local labels) -> (global labels) for each chunk
328 template <unsigned int N, class T, class S1,
329 class Label, class S2,
330 class EqualityFunctor>
331 Label labelMultiArrayBlockwise(const ChunkedArray<N, Data>& data,
332 ChunkedArray<N, Label>& labels,
333 const BlockwiseLabelOptions& options,
334 Equal equal,
335 MultiArrayView<N, std::vector<Label>, S3>& mapping);
336
337 // assign global labels
338 template <unsigned int N, class Data, class S1,
339 class Label, class S2,
340 class Equal>
341 Label labelMultiArrayBlockwise(const MultiArrayView<N, Data, S1>& data,
342 MultiArrayView<N, Label, S2> labels,
343 const BlockwiseLabelOptions& options,
344 Equal equal);
345
346 // assign global labels
347 template <unsigned int N, class T, class S1,
348 class Label, class S2,
349 class EqualityFunctor = std::equal_to<T> >
350 Label labelMultiArrayBlockwise(const ChunkedArray<N, Data>& data,
351 ChunkedArray<N, Label>& labels,
352 const BlockwiseLabelOptions& options = BlockwiseLabelOptions(),
353 Equal equal = std::equal_to<T>());
354 }
355 \endcode
356
357 The resulting labeling is equivalent to a labeling by \ref labelMultiArray, that is, the connected components are the same but may have different ids.
358 \ref NeighborhoodType and background value (if any) can be specified with the LabelOptions object.
359 If the \a mapping parameter is provided, each chunk is labeled seperately and contiguously (starting at one, zero for background),
360 with \a mapping containing a mapping of local labels to global labels for each chunk.
361 Thus, the shape of 'mapping' has to be large enough to hold each chunk coordinate.
362
363 Return: the number of regions found (=largest global region label)
364
365 <b> Usage: </b>
366
367 <b>\#include </b> \<vigra/blockwise_labeling.hxx\><br>
368 Namespace: vigra
369
370 \code
371 Shape3 shape = Shape3(10);
372 Shape3 chunk_shape = Shape3(4);
373 ChunkedArrayLazy<3, int> data(shape, chunk_shape);
374 // fill data ...
375
376 ChunkedArrayLazy<3, size_t> labels(shape, chunk_shape);
377
378 MultiArray<3, std::vector<size_t> > mapping(Shape3(3)); // there are 3 chunks in each dimension
379
380 labelMultiArrayBlockwise(data, labels, LabelOptions().neighborhood(DirectNeighborhood).background(0),
381 std::equal_to<int>(), mapping);
382
383 // check out chunk in the middle
384 MultiArray<3, size_t> middle_chunk(Shape3(4));
385 labels.checkoutSubarray(Shape3(4), middle_chunk);
386
387 // print number of non-background labels assigned in middle_chunk
388 cout << mapping[Shape3(1)].size() << endl;
389
390 // get local label for voxel
391 // this may be the same value assigned to different component in another chunk
392 size_t local_label = middle_chunk[Shape3(2)];
393 // get global label for voxel
394 // if another voxel has the same label, they are in the same connected component albeit they may be in different chunks
395 size_t global_label = mapping[Shape3(1)][local_label
396 \endcode
397 */
doxygen_overloaded_function(template<...> unsigned int labelMultiArrayBlockwise)398 doxygen_overloaded_function(template <...> unsigned int labelMultiArrayBlockwise)
399
400 template <unsigned int N, class Data, class S1,
401 class Label, class S2,
402 class Equal, class S3>
403 Label labelMultiArrayBlockwise(const MultiArrayView<N, Data, S1>& data,
404 MultiArrayView<N, Label, S2> labels,
405 const BlockwiseLabelOptions& options,
406 Equal equal,
407 MultiArrayView<N, std::vector<Label>, S3>& mapping)
408 {
409 using namespace blockwise_labeling_detail;
410
411 typedef typename MultiArrayShape<N>::type Shape;
412 Shape block_shape(options.getBlockShapeN<N>());
413
414 MultiArray<N, MultiArrayView<N, Data, S1> > data_blocks = blockify(data, block_shape);
415 MultiArray<N, MultiArrayView<N, Label, S2> > label_blocks = blockify(labels, block_shape);
416 return blockwiseLabeling(data_blocks.begin(), data_blocks.end(),
417 label_blocks.begin(), label_blocks.end(),
418 options, equal, mapping);
419 }
420
421 template <unsigned int N, class Data, class S1,
422 class Label, class S2,
423 class Equal>
424 Label labelMultiArrayBlockwise(const MultiArrayView<N, Data, S1>& data,
425 MultiArrayView<N, Label, S2> labels,
426 const BlockwiseLabelOptions& options,
427 Equal equal)
428 {
429 using namespace blockwise_labeling_detail;
430
431 typedef typename MultiArrayShape<N>::type Shape;
432 Shape block_shape(options.getBlockShapeN<N>());
433
434 MultiArray<N, MultiArrayView<N, Data, S1> > data_blocks = blockify(data, block_shape);
435 MultiArray<N, MultiArrayView<N, Label, S2> > label_blocks = blockify(labels, block_shape);
436 MultiArray<N, std::vector<Label> > mapping(data_blocks.shape());
437 Label last_label = blockwiseLabeling(data_blocks.begin(), data_blocks.end(),
438 label_blocks.begin(), label_blocks.end(),
439 options, equal, mapping);
440
441 // replace local labels by global labels
442 toGlobalLabels(label_blocks.begin(), label_blocks.end(), mapping.begin(), mapping.end());
443 return last_label;
444 }
445
446 template <unsigned int N, class Data, class S1,
447 class Label, class S2>
448 Label labelMultiArrayBlockwise(const MultiArrayView<N, Data, S1>& data,
449 MultiArrayView<N, Label, S2> labels,
450 const BlockwiseLabelOptions& options = BlockwiseLabelOptions())
451 {
452 return labelMultiArrayBlockwise(data, labels, options, std::equal_to<Data>());
453 }
454
455
456 template <unsigned int N, class Data, class Label, class Equal, class S3>
457 Label labelMultiArrayBlockwise(const ChunkedArray<N, Data>& data,
458 ChunkedArray<N, Label>& labels,
459 const BlockwiseLabelOptions& options,
460 Equal equal,
461 MultiArrayView<N, std::vector<Label>, S3> mapping)
462 {
463 using namespace blockwise_labeling_detail;
464
465 vigra_precondition(options.getBlockShape().size() == 0,
466 "labelMultiArrayBlockwise(ChunkedArray, ...): custom block shapes not supported "
467 "(always uses the array's chunk shape).");
468
469 typedef typename ChunkedArray<N, Data>::shape_type Shape;
470
471 typedef typename ChunkedArray<N, Data>::chunk_const_iterator DataChunkIterator;
472 typedef typename ChunkedArray<N, Label>::chunk_iterator LabelChunkIterator;
473
474 DataChunkIterator data_chunks_begin = data.chunk_begin(Shape(0), data.shape());
475 LabelChunkIterator label_chunks_begin = labels.chunk_begin(Shape(0), labels.shape());
476
477 return blockwiseLabeling(data_chunks_begin, data_chunks_begin.getEndIterator(),
478 label_chunks_begin, label_chunks_begin.getEndIterator(),
479 options, equal, mapping);
480 }
481
482 template <unsigned int N, class Data, class Label, class Equal>
483 Label labelMultiArrayBlockwise(const ChunkedArray<N, Data>& data,
484 ChunkedArray<N, Label>& labels,
485 const BlockwiseLabelOptions& options,
486 Equal equal)
487 {
488 using namespace blockwise_labeling_detail;
489 MultiArray<N, std::vector<Label> > mapping(data.chunkArrayShape());
490 Label result = labelMultiArrayBlockwise(data, labels, options, equal, mapping);
491 typedef typename ChunkedArray<N, Data>::shape_type Shape;
492 toGlobalLabels(labels.chunk_begin(Shape(0), data.shape()), labels.chunk_end(Shape(0), data.shape()), mapping.begin(), mapping.end());
493 return result;
494 }
495
496 template <unsigned int N, class Data, class Label>
497 Label labelMultiArrayBlockwise(const ChunkedArray<N, Data>& data,
498 ChunkedArray<N, Label>& labels,
499 const BlockwiseLabelOptions& options = BlockwiseLabelOptions())
500 {
501 return labelMultiArrayBlockwise(data, labels, options, std::equal_to<Data>());
502 }
503
504 //@}
505
506 } // namespace vigra
507
508 #endif // VIGRA_BLOCKWISE_LABELING_HXX
509