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

watersheds3d.hxx |

1 /************************************************************************/
2 /* */
3 /* Copyright 2006-2007 by F. Heinrich, B. Seppke, 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_watersheds3D_HXX
37 #define VIGRA_watersheds3D_HXX
38
39 #include "voxelneighborhood.hxx"
40 #include "multi_array.hxx"
41 #include "multi_localminmax.hxx"
42 #include "labelvolume.hxx"
43 #include "seededregiongrowing3d.hxx"
44 #include "watersheds.hxx"
45
46 namespace vigra
47 {
48
49 template <class SrcIterator, class SrcAccessor, class SrcShape,
50 class DestIterator, class DestAccessor, class Neighborhood3D>
51 int preparewatersheds3D( SrcIterator s_Iter, SrcShape srcShape, SrcAccessor sa,
52 DestIterator d_Iter, DestAccessor da, Neighborhood3D)
53 {
54 //basically needed for iteration and border-checks
55 int w = srcShape[0], h = srcShape[1], d = srcShape[2];
56 int x,y,z, local_min_count=0;
57
58 //declare and define Iterators for all three dims at src
59 SrcIterator zs = s_Iter;
60 SrcIterator ys(zs);
61 SrcIterator xs(ys);
62
63 //Declare Iterators for all three dims at dest
64 DestIterator zd = d_Iter;
65
66 for(z = 0; z != d; ++z, ++zs.dim2(), ++zd.dim2())
67 {
68 ys = zs;
69 DestIterator yd(zd);
70
71 for(y = 0; y != h; ++y, ++ys.dim1(), ++yd.dim1())
72 {
73 xs = ys;
74 DestIterator xd(yd);
75
76 for(x = 0; x != w; ++x, ++xs.dim0(), ++xd.dim0())
77 {
78 AtVolumeBorder atBorder = isAtVolumeBorder(x,y,z,w,h,d);
79 typename SrcAccessor::value_type v = sa(xs);
80 // the following choice causes minima to point
81 // to their lowest neighbor -- would this be better???
82 // typename SrcAccessor::value_type v = NumericTraits<typename SrcAccessor::value_type>::max();
83 int o = 0; // means center is minimum
84 typename SrcAccessor::value_type my_v = v;
85 if(atBorder == NotAtBorder)
86 {
87 NeighborhoodCirculator<SrcIterator, Neighborhood3D> c(xs), cend(c);
88
89 do {
90 if(sa(c) < v)
91 {
92 v = sa(c);
93 o = c.directionBit();
94 }
95 else if(sa(c) == my_v && my_v == v)
96 {
97 o = o | c.directionBit();
98 }
99 }
100 while(++c != cend);
101 }
102 else
103 {
104 RestrictedNeighborhoodCirculator<SrcIterator, Neighborhood3D> c(xs, atBorder), cend(c);
105 do {
106 if(sa(c) < v)
107 {
108 v = sa(c);
109 o = c.directionBit();
110 }
111 else if(sa(c) == my_v && my_v == v)
112 {
113 o = o | c.directionBit();
114 }
115 }
116 while(++c != cend);
117 }
118 if (o==0) local_min_count++;
119 da.set(o, xd);
120 }//end x-iteration
121 }//end y-iteration
122 }//end z-iteration
123 return local_min_count;
124 }
125
126 template <class SrcIterator, class SrcAccessor,class SrcShape,
127 class DestIterator, class DestAccessor,
128 class Neighborhood3D>
129 unsigned int watershedLabeling3D( SrcIterator s_Iter, SrcShape srcShape, SrcAccessor sa,
130 DestIterator d_Iter, DestAccessor da,
131 Neighborhood3D)
132 {
133 typedef typename DestAccessor::value_type LabelType;
134
135 //basically needed for iteration and border-checks
136 int w = srcShape[0], h = srcShape[1], d = srcShape[2];
137 int x,y,z;
138
139 //declare and define Iterators for all three dims at src
140 SrcIterator zs = s_Iter;
141 DestIterator zd = d_Iter;
142
143 // temporary image to store region labels
144 UnionFindArray<LabelType> labels;
145
146 // initialize the neighborhood traversers
147 NeighborOffsetCirculator<Neighborhood3D> nc(Neighborhood3D::CausalFirst);
148 NeighborOffsetCirculator<Neighborhood3D> nce(Neighborhood3D::CausalLast);
149 ++nce;
150 // pass 1: scan image from upper left front to lower right back
151 // to find connected components
152
153 // Each component will be represented by a tree of pixels. Each
154 // pixel contains the scan order address of its parent in the
155 // tree. In order for pass 2 to work correctly, the parent must
156 // always have a smaller scan order address than the child.
157 // Therefore, we can merge trees only at their roots, because the
158 // root of the combined tree must have the smallest scan order
159 // address among all the tree's pixels/ nodes. The root of each
160 // tree is distinguished by pointing to itself (it contains its
161 // own scan order address). This condition is enforced whenever a
162 // new region is found or two regions are merged
163 for(z = 0; z != d; ++z, ++zs.dim2(), ++zd.dim2())
164 {
165 SrcIterator ys = zs;
166 DestIterator yd = zd;
167
168 for(y = 0; y != h; ++y, ++ys.dim1(), ++yd.dim1())
169 {
170 SrcIterator xs = ys;
171 DestIterator xd = yd;
172
173 for(x = 0; x != w; ++x, ++xs.dim0(), ++xd.dim0())
174 {
175 LabelType currentIndex = labels.nextFreeIndex(); // default: new region
176
177 //check whether there is a special border treatment to be used or not
178 AtVolumeBorder atBorder = isAtVolumeBorderCausal(x,y,z,w,h,d);
179
180 //We are not at the border!
181 if(atBorder == NotAtBorder)
182 {
183
184 nc = NeighborOffsetCirculator<Neighborhood3D>(Neighborhood3D::CausalFirst);
185
186 do
187 {
188 // Direction of NTraversr Neighbor's direction bit is pointing
189 // = Direction of voxel towards us?
190 if((sa(xs) & nc.directionBit()) || (sa(xs,*nc) & nc.oppositeDirectionBit()))
191 {
192 currentIndex = labels.makeUnion(da(xd,*nc), currentIndex);
193 }
194 ++nc;
195 }while(nc!=nce);
196 }
197 //we are at a border - handle this!!
198 else
199 {
200 nc = NeighborOffsetCirculator<Neighborhood3D>(Neighborhood3D::nearBorderDirectionsCausal(atBorder,0));
201 int j=0;
202 while(nc.direction() != Neighborhood3D::Error)
203 {
204 int dummy = x+(*nc)[0]; // prevents an apparently incorrect optimization in gcc 4.8
205 if (dummy<0)
206 {
207 std::cerr << "internal error " << dummy << std::endl;
208 }
209 // Direction of NTraversr Neighbor's direction bit is pointing
210 // = Direction of voxel towards us?
211 if((sa(xs) & nc.directionBit()) || (sa(xs,*nc) & nc.oppositeDirectionBit()))
212 {
213 currentIndex = labels.makeUnion(da(xd,*nc), currentIndex);
214 }
215 nc.turnTo(Neighborhood3D::nearBorderDirectionsCausal(atBorder,++j));
216 }
217 }
218 da.set(labels.finalizeIndex(currentIndex), xd);
219 }
220 }
221 }
222
223 unsigned int count = labels.makeContiguous();
224
225 // pass 2: assign one label to each region (tree)
226 // so that labels form a consecutive sequence 1, 2, ...
227 zd = d_Iter;
228 for(z=0; z != d; ++z, ++zd.dim2())
229 {
230 DestIterator yd(zd);
231
232 for(y=0; y != h; ++y, ++yd.dim1())
233 {
234 DestIterator xd(yd);
235
236 for(x = 0; x != w; ++x, ++xd.dim0())
237 {
238 da.set(labels.findLabel(da(xd)), xd);
239 }
240 }
241 }
242 return count;
243 }
244
245
246 /** \addtogroup SeededRegionGrowing
247 */
248 //@{
249
250 /********************************************************/
251 /* */
252 /* watersheds3D */
253 /* */
254 /********************************************************/
255
256 /** \brief Region Segmentation by means of the watershed algorithm.
257
258 This function is deprecated, use \ref watershedsMultiArray() instead.
259
260 <b> Declarations:</b>
261
262 \deprecatedAPI{watersheds3D}
263 pass \ref MultiIteratorPage "MultiIterators" and \ref DataAccessors :
264 \code
265 namespace vigra {
266 template <class SrcIterator, class SrcAccessor,class SrcShape,
267 class DestIterator, class DestAccessor,
268 class Neighborhood3D>
269 unsigned int watersheds3D(SrcIterator s_Iter, SrcShape srcShape, SrcAccessor sa,
270 DestIterator d_Iter, DestAccessor da,
271 Neighborhood3D neighborhood3D);
272 }
273 \endcode
274 use argument objects in conjunction with \ref ArgumentObjectFactories :
275 \code
276 namespace vigra {
277 template <class SrcIterator, class SrcAccessor,class SrcShape,
278 class DestIterator, class DestAccessor,
279 class Neighborhood3D>
280 unsigned int watersheds3D(triple<SrcIterator, SrcShape, SrcAccessor> src,
281 pair<DestIterator, DestAccessor> dest,
282 Neighborhood3D neighborhood3D);
283 }
284 \endcode
285
286 use with 3D-Six-Neighborhood:
287 \code
288 namespace vigra {
289
290 template <class SrcIterator, class SrcAccessor,class SrcShape,
291 class DestIterator, class DestAccessor>
292 unsigned int watersheds3DSix(triple<SrcIterator, SrcShape, SrcAccessor> src,
293 pair<DestIterator, DestAccessor> dest);
294
295 }
296 \endcode
297
298 use with 3D-TwentySix-Neighborhood:
299 \code
300 namespace vigra {
301
302 template <class SrcIterator, class SrcAccessor,class SrcShape,
303 class DestIterator, class DestAccessor>
304 unsigned int watersheds3DTwentySix(triple<SrcIterator, SrcShape, SrcAccessor> src,
305 pair<DestIterator, DestAccessor> dest);
306
307 }
308 \endcode
309 \deprecatedEnd
310
311 This function implements the union-find version of the watershed algorithms
312 as described in
313
314 J. Roerdink, R. Meijster: <em>"The watershed transform: definitions, algorithms,
315 and parallelization strategies"</em>, Fundamenta Informaticae, 41:187-228, 2000
316
317 The source volume is a boundary indicator such as the gradient magnitude
318 of the trace of the \ref boundaryTensor(). Local minima of the boundary indicator
319 are used as region seeds, and all other voxels are recursively assigned to the same
320 region as their lowest neighbor. Pass \ref vigra::NeighborCode3DSix or
321 \ref vigra::NeighborCode3DTwentySix to determine the neighborhood where voxel values
322 are compared. The voxel type of the input volume must be <tt>LessThanComparable</tt>.
323
324 <b> Usage:</b>
325
326 <b>\#include</b> <vigra/watersheds3D.hxx><br>
327 Namespace: vigra
328
329 Example: watersheds3D of the gradient magnitude.
330
331 \code
332 Shape3 shape(w, h, d);
333
334 MultiArray<3, float> src(shape), grad(shape);
335 ...
336
337 double scale = 1;
338 gaussianGradientMagnitude(src, grad, scale);
339
340 MultiArray<3, int> labels(shape);
341
342 // find 6-connected regions
343 int max_region_label = watersheds3DSix(grad, labels);
344
345 // find 26-connected regions
346 max_region_label = watersheds3DTwentySix(grad, labels);
347 \endcode
348
349 \deprecatedUsage{watersheds3D}
350 \code
351 typedef vigra::MultiArray<3,int> IntVolume;
352 typedef vigra::MultiArray<3,double> DVolume;
353 DVolume src(DVolume::difference_type(w,h,d));
354 IntVolume dest(IntVolume::difference_type(w,h,d));
355
356 float gauss=1;
357
358 vigra::MultiArray<3, vigra::TinyVector<float,3> > temp(IntVolume::difference_type(w,h,d));
359 vigra::gaussianGradientMultiArray(srcMultiArrayRange(vol),destMultiArray(temp),gauss);
360
361 IntVolume::iterator temp_iter=temp.begin();
362 for(DVolume::iterator iter=src.begin(); iter!=src.end(); ++iter, ++temp_iter)
363 *iter = norm(*temp_iter);
364
365 // find 6-connected regions
366 int max_region_label = vigra::watersheds3DSix(srcMultiArrayRange(src), destMultiArray(dest));
367
368 // find 26-connected regions
369 max_region_label = vigra::watersheds3DTwentySix(srcMultiArrayRange(src), destMultiArray(dest));
370
371 \endcode
372 <b> Required Interface:</b>
373 \code
374 SrcIterator src_begin;
375 SrcShape src_shape;
376 DestIterator dest_begin;
377
378 SrcAccessor src_accessor;
379 DestAccessor dest_accessor;
380
381 // compare src values
382 src_accessor(src_begin) <= src_accessor(src_begin)
383
384 // set result
385 int label;
386 dest_accessor.set(label, dest_begin);
387 \endcode
388 \deprecatedEnd
389 */
390 doxygen_overloaded_function(template <...> unsigned int watersheds3D)
391
392 template <class SrcIterator, class SrcAccessor, class SrcShape,
393 class DestIterator, class DestAccessor,
394 class Neighborhood3D>
395 unsigned int watersheds3D( SrcIterator s_Iter, SrcShape srcShape, SrcAccessor sa,
396 DestIterator d_Iter, DestAccessor da, Neighborhood3D neighborhood3D)
397 {
398 //create temporary volume to store the DAG of directions to minima
399 if ((int)Neighborhood3D::DirectionCount>7){ //If we have 3D-TwentySix Neighborhood
400
401 vigra::MultiArray<3,int> orientationVolume(srcShape);
402
403 preparewatersheds3D( s_Iter, srcShape, sa,
404 destMultiArray(orientationVolume).first, destMultiArray(orientationVolume).second,
405 neighborhood3D);
406
407 return watershedLabeling3D( srcMultiArray(orientationVolume).first, srcShape, srcMultiArray(orientationVolume).second,
408 d_Iter, da,
409 neighborhood3D);
410 }
411 else{
412
413 vigra::MultiArray<3,unsigned char> orientationVolume(srcShape);
414
415 preparewatersheds3D( s_Iter, srcShape, sa,
416 destMultiArray(orientationVolume).first, destMultiArray(orientationVolume).second,
417 neighborhood3D);
418
419 return watershedLabeling3D( srcMultiArray(orientationVolume).first, srcShape, srcMultiArray(orientationVolume).second,
420 d_Iter, da,
421 neighborhood3D);
422 }
423 }
424
425 template <class SrcIterator, class SrcShape, class SrcAccessor,
426 class DestIterator, class DestAccessor>
427 inline unsigned int watersheds3DSix( triple<SrcIterator, SrcShape, SrcAccessor> src,
428 pair<DestIterator, DestAccessor> dest)
429 {
430 return watersheds3D(src.first, src.second, src.third, dest.first, dest.second, NeighborCode3DSix());
431 }
432
433 template <class SrcIterator, class SrcShape, class SrcAccessor,
434 class DestIterator, class DestAccessor>
435 inline unsigned int watersheds3DTwentySix( triple<SrcIterator, SrcShape, SrcAccessor> src,
436 pair<DestIterator, DestAccessor> dest)
437 {
438 return watersheds3D(src.first, src.second, src.third, dest.first, dest.second, NeighborCode3DTwentySix());
439 }
440
441 template <unsigned int N, class T1, class S1,
442 class T2, class S2>
443 inline unsigned int
444 watersheds3DSix(MultiArrayView<N, T1, S1> const & source,
445 MultiArrayView<N, T2, S2> dest)
446 {
447 vigra_precondition(source.shape() == dest.shape(),
448 "watersheds3DSix(): shape mismatch between input and output.");
449 return watersheds3DSix(srcMultiArrayRange(source), destMultiArray(dest));
450 }
451
452 template <unsigned int N, class T1, class S1,
453 class T2, class S2>
454 inline unsigned int
455 watersheds3DTwentySix(MultiArrayView<N, T1, S1> const & source,
456 MultiArrayView<N, T2, S2> dest)
457 {
458 vigra_precondition(source.shape() == dest.shape(),
459 "watersheds3DTwentySix(): shape mismatch between input and output.");
460 return watersheds3DTwentySix(srcMultiArrayRange(source), destMultiArray(dest));
461 }
462
463 }//namespace vigra
464
465 #endif //VIGRA_watersheds3D_HXX
vigra::AtVolumeBorder
AtImageBorder AtVolumeBorder
Encode whether a voxel is near the volume border.
Definition: voxelneighborhood.hxx:72
vigra::watersheds3D
unsigned int watersheds3D(...)
Region Segmentation by means of the watershed algorithm.
vigra::NeighborCode3DTwentySix
Neighborhood3DTwentySix::NeighborCode3D NeighborCode3DTwentySix
Definition: voxelneighborhood.hxx:1646
vigra::MultiArray
Main MultiArray class containing the memory management.
Definition: multi_array.hxx:2422
vigra::isAtVolumeBorderCausal
AtVolumeBorder isAtVolumeBorderCausal(int x, int y, int z, int width, int height, int)
Find out whether a voxel is at a scan-order relevant volume border. This function checks if x == 0 or...
Definition: voxelneighborhood.hxx:112
vigra::isAtVolumeBorder
AtVolumeBorder isAtVolumeBorder(int x, int y, int z, int width, int height, int depth)
Find out whether a voxel is at the volume border.
Definition: voxelneighborhood.hxx:82
vigra
Definition: accessor.hxx:43
vigra::doxygen_overloaded_function
doxygen_overloaded_function(template<... > void separableConvolveBlockwise) template< unsigned int N
Separated convolution on ChunkedArrays.
vigra::boost_graph::source
vigra::GridGraph< N, DirectedTag >::vertex_descriptor source(typename vigra::GridGraph< N, DirectedTag >::edge_descriptor const &e, vigra::GridGraph< N, DirectedTag > const &g)
Get a vertex descriptor for the start vertex of edge e in graph g (API: boost).
Definition: multi_gridgraph.hxx:2943
vigra::NeighborCode3DSix
Neighborhood3DSix::NeighborCode3D NeighborCode3DSix
Definition: voxelneighborhood.hxx:490
vigra::NotAtBorder
&#160;
Definition: pixelneighborhood.hxx:70

© 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