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

blockwise_watersheds.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_WATERSHEDS_HXX
37 #define VIGRA_BLOCKWISE_WATERSHEDS_HXX
38 
39 #include "multi_array.hxx"
40 #include "multi_gridgraph.hxx"
41 #include "blockify.hxx"
42 #include "blockwise_labeling.hxx"
43 #include "overlapped_blocks.hxx"
44 
45 #include <limits>
46 
47 namespace vigra
48 {
49 
50 /** \addtogroup SeededRegionGrowing
51 */
52 //@{
53 
54 namespace blockwise_watersheds_detail
55 {
56 
57 template <class DataArray, class DirectionsBlocksIterator>
58 void prepareBlockwiseWatersheds(const Overlaps<DataArray>& overlaps,
59  DirectionsBlocksIterator directions_blocks_begin,
60  NeighborhoodType neighborhood)
61 {
62  static const unsigned int N = DataArray::actual_dimension;
63  typedef typename MultiArrayShape<N>::type Shape;
64  typedef typename DirectionsBlocksIterator::value_type DirectionsBlock;
65  Shape shape = overlaps.shape();
66  vigra_assert(shape == directions_blocks_begin.shape(), "");
67 
68  MultiCoordinateIterator<N> it(shape);
69  MultiCoordinateIterator<N> end = it.getEndIterator();
70  for( ; it != end; ++it)
71  {
72  DirectionsBlock directions_block = directions_blocks_begin[*it];
73  OverlappingBlock<DataArray> data_block = overlaps[*it];
74 
75  typedef GridGraph<N, undirected_tag> Graph;
76  typedef typename Graph::NodeIt GraphScanner;
77  typedef typename Graph::OutArcIt NeighborIterator;
78 
79  Graph graph(data_block.block.shape(), neighborhood);
80  for(GraphScanner node(graph); node != lemon::INVALID; ++node)
81  {
82  if(within(*node, data_block.inner_bounds))
83  {
84  typedef typename DataArray::value_type Data;
85  Data lowest_neighbor = data_block.block[*node];
86 
87  typedef typename DirectionsBlock::value_type Direction;
88  Direction lowest_neighbor_direction = std::numeric_limits<unsigned short>::max();
89 
90  for(NeighborIterator arc(graph, *node); arc != lemon::INVALID; ++arc)
91  {
92  Shape neighbor_coordinates = graph.target(*arc);
93  Data neighbor_data = data_block.block[neighbor_coordinates];
94  if(neighbor_data < lowest_neighbor)
95  {
96  lowest_neighbor = neighbor_data;
97  lowest_neighbor_direction = arc.neighborIndex();
98  }
99  }
100  directions_block[*node - data_block.inner_bounds.first] = lowest_neighbor_direction;
101  }
102  }
103  }
104 }
105 
106 template <unsigned int N>
107 struct UnionFindWatershedsEquality
108 {
109  // FIXME: this graph object shouldn't be necessary, most functions (and state) of graph are not used
110  // this probably needs some refactoring in GridGraph
111  GridGraph<N, undirected_tag>* graph;
112 
113  template <class Shape>
114  bool operator()(unsigned short u, const unsigned short v, const Shape& diff) const
115  {
116  static const unsigned short plateau_id = std::numeric_limits<unsigned short>::max();
117  return (u == plateau_id && v == plateau_id) ||
118  (u != plateau_id && graph->neighborOffset(u) == diff) ||
119  (v != plateau_id && graph->neighborOffset(graph->oppositeIndex(v)) == diff);
120  }
121 
122  struct WithDiffTag
123  {};
124 };
125 
126 } // namespace blockwise_watersheds_detail
127 
128 template <unsigned int N, class Data, class S1,
129  class Label, class S2>
130 Label unionFindWatershedsBlockwise(MultiArrayView<N, Data, S1> data,
131  MultiArrayView<N, Label, S2> labels,
132  NeighborhoodType neighborhood = DirectNeighborhood,
133  const typename MultiArrayView<N, Data, S1>::difference_type& block_shape =
134  typename MultiArrayView<N, Data, S1>::difference_type(128))
135 {
136  using namespace blockwise_watersheds_detail;
137 
138  typedef typename MultiArrayView<N, Data, S1>::difference_type Shape;
139  Shape shape = data.shape();
140  vigra_precondition(shape == labels.shape(), "shapes of data and labels do not match");
141 
142  MultiArray<N, unsigned short> directions(shape);
143 
144  MultiArray<N, MultiArrayView<N, unsigned short> > directions_blocks = blockify(directions, block_shape);
145 
146  Overlaps<MultiArrayView<N, Data, S1> > overlaps(data, block_shape, Shape(1), Shape(1));
147  prepareBlockwiseWatersheds(overlaps, directions_blocks.begin(), neighborhood);
148  GridGraph<N, undirected_tag> graph(data.shape(), neighborhood);
149  UnionFindWatershedsEquality<N> equal = {&graph};
150  return labelMultiArrayBlockwise(directions, labels, LabelOptions().neighborhood(neighborhood).blockShape(block_shape), equal);
151 }
152 
153 /*************************************************************/
154 /* */
155 /* unionFindWatershedsBlockwise */
156 /* */
157 /*************************************************************/
158 
159 /** \brief Blockise union-find watersheds transform for ChunkedArrays.
160 
161  <b> Declaration:</b>
162 
163  \code
164  namespace vigra {
165  template <unsigned int N, class Data, class Label>
166  Label unionFindWatershedsBlockwise(const ChunkedArray<N, Data>& data,
167  ChunkedArray<N, Label>& labels,
168  NeighborhoodType neighborhood = DirectNeighborhood);
169 
170  // provide temporary directions storage
171  template <unsigned int N, class Data, class Label>
172  Label unionFindWatershedsBlockwise(const ChunkedArray<N, Data>& data,
173  ChunkedArray<N, Label>& labels,
174  NeighborhoodType neighborhood,
175  ChunkedArray<N, unsigned short>& temporary_storage);
176  }
177  \endcode
178 
179  The resulting labeling is equivalent to a labeling by \ref watershedsUnionFind, that is,
180  the components are the same but may have different ids.
181  If \a temporary_storage is provided, this array is used for intermediate result storage.
182  Otherwise, a newly created \ref vigra::ChunkedArrayLazy is used.
183 
184  Return: the number of labels assigned (=largest label, because labels start at one)
185 
186  <b> Usage: </b>
187 
188  <b>\#include </b> <vigra/blockwise_watersheds.hxx><br>
189  Namespace: vigra
190 
191  \code
192  Shape3 shape = Shape3(10);
193  Shape3 chunk_shape = Shape3(4);
194  ChunkedArrayLazy<3, int> data(shape, chunk_shape);
195  // fill data ...
196 
197  ChunkedArrayLazy<3, size_t> labels(shape, chunk_shape);
198 
199  unionFindWatershedsBlockwise(data, labels, IndirectNeighborhood);
200  \endcode
201  */
202 doxygen_overloaded_function(template <...> unsigned int unionFindWatershedsBlockwise)
203 
204 template <unsigned int N, class Data, class Label>
205 Label unionFindWatershedsBlockwise(const ChunkedArray<N, Data>& data,
206  ChunkedArray<N, Label>& labels,
207  NeighborhoodType neighborhood,
208  ChunkedArray<N, unsigned short>& directions)
209 {
210  using namespace blockwise_watersheds_detail;
211 
212  typedef typename ChunkedArray<N, Data>::shape_type Shape;
213  Shape shape = data.shape();
214  vigra_precondition(shape == labels.shape() && shape == directions.shape(), "shapes of data and labels do not match");
215  Shape chunk_shape = data.chunkShape();
216  vigra_precondition(chunk_shape == labels.chunkShape() && chunk_shape == directions.chunkShape(), "chunk shapes do not match");
217 
218  Overlaps<ChunkedArray<N, Data> > overlaps(data, data.chunkShape(), Shape(1), Shape(1));
219 
220  prepareBlockwiseWatersheds(overlaps, directions.chunk_begin(Shape(0), shape), neighborhood);
221 
222  GridGraph<N, undirected_tag> graph(shape, neighborhood);
223  UnionFindWatershedsEquality<N> equal = {&graph};
224  return labelMultiArrayBlockwise(directions, labels, LabelOptions().neighborhood(neighborhood), equal);
225 }
226 
227 template <unsigned int N, class Data,
228  class Label>
229 inline Label
230 unionFindWatershedsBlockwise(const ChunkedArray<N, Data>& data,
231  ChunkedArray<N, Label>& labels,
232  NeighborhoodType neighborhood = DirectNeighborhood)
233 {
234  ChunkedArrayLazy<N, unsigned short> directions(data.shape(), data.chunkShape());
235  return unionFindWatershedsBlockwise(data, labels, neighborhood, directions);
236 }
237 
238 //@}
239 
240 } // namespace vigra
241 
242 #endif // VIGRA_BLOCKWISE_WATERSHEDS_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)