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

seededregiongrowing3d.hxx
1 /************************************************************************/
2 /* */
3 /* Copyright 2003-2007 by Kasim Terzic, Christian-Dennis Rahn */
4 /* and Ullrich Koethe */
5 /* */
6 /* This file is part of the VIGRA computer vision library. */
7 /* The VIGRA Website is */
8 /* http://hci.iwr.uni-heidelberg.de/vigra/ */
9 /* Please direct questions, bug reports, and contributions to */
10 /* ullrich.koethe@iwr.uni-heidelberg.de or */
11 /* vigra@informatik.uni-hamburg.de */
12 /* */
13 /* Permission is hereby granted, free of charge, to any person */
14 /* obtaining a copy of this software and associated documentation */
15 /* files (the "Software"), to deal in the Software without */
16 /* restriction, including without limitation the rights to use, */
17 /* copy, modify, merge, publish, distribute, sublicense, and/or */
18 /* sell copies of the Software, and to permit persons to whom the */
19 /* Software is furnished to do so, subject to the following */
20 /* conditions: */
21 /* */
22 /* The above copyright notice and this permission notice shall be */
23 /* included in all copies or substantial portions of the */
24 /* Software. */
25 /* */
26 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */
27 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */
28 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */
29 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */
30 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */
31 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */
32 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */
33 /* OTHER DEALINGS IN THE SOFTWARE. */
34 /* */
35 /************************************************************************/
36 
37 #ifndef VIGRA_SEEDEDREGIONGROWING_3D_HXX
38 #define VIGRA_SEEDEDREGIONGROWING_3D_HXX
39 
40 #include <vector>
41 #include <stack>
42 #include <queue>
43 #include "utilities.hxx"
44 #include "stdimage.hxx"
45 #include "stdimagefunctions.hxx"
46 #include "seededregiongrowing.hxx"
47 #include "multi_pointoperators.hxx"
48 #include "voxelneighborhood.hxx"
49 
50 namespace vigra {
51 
52 namespace detail {
53 
54 template <class COST, class Diff_type>
55 class SeedRgVoxel
56 {
57 public:
58  Diff_type location_, nearest_;
59  COST cost_;
60  int count_;
61  int label_;
62  int dist_;
63 
64  SeedRgVoxel()
65  //: location_(0,0,0), nearest_(0,0,0), cost_(0), count_(0), label_(0)
66  {
67  location_ = Diff_type(0,0,0);
68  nearest_ = Diff_type(0,0,0);
69  cost_ = 0;
70  count_ = 0;
71  label_ = 0;
72  }
73 
74  SeedRgVoxel(Diff_type const & location, Diff_type const & nearest,
75  COST const & cost, int const & count, int const & label)
76  : location_(location), nearest_(nearest),
77  cost_(cost), count_(count), label_(label)
78  {
79  int dx = location_[0] - nearest_[0];
80  int dy = location_[1] - nearest_[1];
81  int dz = location_[2] - nearest_[2];
82  dist_ = dx * dx + dy * dy + dz * dz;
83  }
84 
85  void set(Diff_type const & location, Diff_type const & nearest,
86  COST const & cost, int const & count, int const & label)
87  {
88  location_ = location;
89  nearest_ = nearest;
90  cost_ = cost;
91  count_ = count;
92  label_ = label;
93 
94  int dx = location_[0] - nearest_[0];
95  int dy = location_[1] - nearest_[1];
96  int dz = location_[2] - nearest_[2];
97  dist_ = dx * dx + dy * dy + dz * dz;
98  }
99 
100  struct Compare
101  {
102  // must implement > since priority_queue looks for largest element
103  bool operator()(SeedRgVoxel const & l,
104  SeedRgVoxel const & r) const
105  {
106  if(r.cost_ == l.cost_)
107  {
108  if(r.dist_ == l.dist_) return r.count_ < l.count_;
109 
110  return r.dist_ < l.dist_;
111  }
112 
113  return r.cost_ < l.cost_;
114  }
115  bool operator()(SeedRgVoxel const * l,
116  SeedRgVoxel const * r) const
117  {
118  if(r->cost_ == l->cost_)
119  {
120  if(r->dist_ == l->dist_) return r->count_ < l->count_;
121 
122  return r->dist_ < l->dist_;
123  }
124 
125  return r->cost_ < l->cost_;
126  }
127  };
128 
129  struct Allocator
130  {
131  ~Allocator()
132  {
133  while(!freelist_.empty())
134  {
135  delete freelist_.top();
136  freelist_.pop();
137  }
138  }
139 
140  SeedRgVoxel * create(Diff_type const & location, Diff_type const & nearest,
141  COST const & cost, int const & count, int const & label)
142  {
143  if(!freelist_.empty())
144  {
145  SeedRgVoxel * res = freelist_.top();
146  freelist_.pop();
147  res->set(location, nearest, cost, count, label);
148  return res;
149  }
150 
151  return new SeedRgVoxel(location, nearest, cost, count, label);
152  }
153 
154  void dismiss(SeedRgVoxel * p)
155  {
156  freelist_.push(p);
157  }
158 
159  std::stack<SeedRgVoxel<COST,Diff_type> *> freelist_;
160  };
161 };
162 
163 } // namespace detail
164 
165 /** \addtogroup SeededRegionGrowing
166  Region segmentation and voronoi tesselation
167 */
168 //@{
169 
170 /********************************************************/
171 /* */
172 /* seededRegionGrowing3D */
173 /* */
174 /********************************************************/
175 
176 /** \brief Three-dimensional Region Segmentation by means of Seeded Region Growing.
177 
178  This algorithm implements seeded region growing as described in
179 
180  The seed image is a partly segmented multi-dimensional array which contains uniquely
181  labeled regions (the seeds) and unlabeled voxels (the candidates, label 0).
182  Seed regions can be as large as you wish and as small as one voxel. If
183  there are no candidates, the algorithm will simply copy the seed array
184  into the output array. Otherwise it will aggregate the candidates into
185  the existing regions so that a cost function is minimized.
186  Candidates are taken from the neighborhood of the already assigned pixels,
187  where the type of neighborhood is determined by parameter <tt>neighborhood</tt>
188  which can take the values <tt>NeighborCode3DSix()</tt> (the default)
189  or <tt>NeighborCode3DTwentySix()</tt>. The algorithm basically works as follows
190  (illustrated for 6-neighborhood, but 26-neighborhood works in the same way):
191 
192  <ol>
193 
194  <li> Find all candidate pixels that are 6-adjacent to a seed region.
195  Calculate the cost for aggregating each candidate into its adajacent region
196  and put the candidates into a priority queue.
197 
198  <li> While( priority queue is not empty)
199 
200  <ol>
201 
202  <li> Take the candidate with least cost from the queue. If it has not
203  already been merged, merge it with it's adjacent region.
204 
205  <li> Put all candidates that are 4-adjacent to the pixel just processed
206  into the priority queue.
207 
208  </ol>
209 
210  </ol>
211 
212  <tt>SRGType</tt> can take the following values:
213 
214  <DL>
215  <DT><tt>CompleteGrow</tt> <DD> produce a complete tesselation of the volume (default).
216  <DT><tt>KeepContours</tt> <DD> keep a 1-voxel wide unlabeled contour between all regions.
217  <DT><tt>StopAtThreshold</tt> <DD> stop when the boundary indicator values exceed the
218  threshold given by parameter <tt>max_cost</tt>.
219  <DT><tt>KeepContours | StopAtThreshold</tt> <DD> keep 1-voxel wide contour and stop at given <tt>max_cost</tt>.
220  </DL>
221 
222  The cost is determined jointly by the source array and the
223  region statistics functor. The source array contains feature values for each
224  pixel which will be used by the region statistics functor to calculate and
225  update statistics for each region and to calculate the cost for each
226  candidate. The <TT>RegionStatisticsArray</TT> must be compatible to the
227  \ref ArrayOfRegionStatistics functor and contains an <em> array</em> of
228  statistics objects for each region. The indices must correspond to the
229  labels of the seed regions. The statistics for the initial regions must have
230  been calculated prior to calling <TT>seededRegionGrowing3D()</TT>
231 
232  For each candidate
233  <TT>x</TT> that is adjacent to region <TT>i</TT>, the algorithm will call
234  <TT>stats[i].cost(as(x))</TT> to get the cost (where <TT>x</TT> is a <TT>SrcImageIterator</TT>
235  and <TT>as</TT> is
236  the SrcAccessor). When a candidate has been merged with a region, the
237  statistics are updated by calling <TT>stats[i].operator()(as(x))</TT>. Since
238  the <TT>RegionStatisticsArray</TT> is passed by reference, this will overwrite
239  the original statistics.
240 
241  If a candidate could be merged into more than one regions with identical
242  cost, the algorithm will favour the nearest region. If <tt>StopAtThreshold</tt> is active,
243  and the cost of the current candidate at any point in the algorithm exceeds the optional
244  <tt>max_cost</tt> value (which defaults to <tt>NumericTraits<double>::max()</tt>),
245  region growing is aborted, and all voxels not yet assigned to a region remain unlabeled.
246 
247  In some cases, the cost only depends on the feature value of the current
248  voxel. Then the update operation will simply be a no-op, and the <TT>cost()</TT>
249  function returns its argument. This behavior is implemented by the
250  \ref SeedRgDirectValueFunctor.
251 
252  <b> Declarations:</b>
253 
254  pass arguments explicitly:
255  \code
256  namespace vigra {
257  template <class SrcImageIterator, class Shape, class SrcAccessor,
258  class SeedImageIterator, class SeedAccessor,
259  class DestImageIterator, class DestAccessor,
260  class RegionStatisticsArray, class Neighborhood>
261  void
262  seededRegionGrowing3D(SrcImageIterator srcul, Shape shape, SrcAccessor as,
263  SeedImageIterator seedsul, SeedAccessor aseeds,
264  DestImageIterator destul, DestAccessor ad,
265  RegionStatisticsArray & stats,
266  SRGType srgType = CompleteGrow,
267  Neighborhood neighborhood = NeighborCode3DSix(),
268  double max_cost = NumericTraits<double>::max());
269  }
270  \endcode
271 
272  use argument objects in conjunction with \ref ArgumentObjectFactories :
273  \code
274  namespace vigra {
275  template <class SrcImageIterator, class Shape, class SrcAccessor,
276  class SeedImageIterator, class SeedAccessor,
277  class DestImageIterator, class DestAccessor,
278  class RegionStatisticsArray, class Neighborhood>
279  void
280  seededRegionGrowing3D(triple<SrcImageIterator, Shape, SrcAccessor> src,
281  pair<SeedImageIterator, SeedAccessor> seeds,
282  pair<DestImageIterator, DestAccessor> dest,
283  RegionStatisticsArray & stats,
284  SRGType srgType = CompleteGrow,
285  Neighborhood neighborhood = NeighborCode3DSix(),
286  double max_cost = NumericTraits<double>::max());
287  }
288  \endcode
289 
290 */
292 
293 template <class SrcImageIterator, class Diff_type, class SrcAccessor,
294  class SeedImageIterator, class SeedAccessor,
295  class DestImageIterator, class DestAccessor,
296  class RegionStatisticsArray, class Neighborhood>
297 void
298 seededRegionGrowing3D(SrcImageIterator srcul, Diff_type shape, SrcAccessor as,
299  SeedImageIterator seedsul, SeedAccessor aseeds,
300  DestImageIterator destul, DestAccessor ad,
301  RegionStatisticsArray & stats,
302  SRGType srgType,
303  Neighborhood,
304  double max_cost)
305 {
306  SrcImageIterator srclr = srcul + shape;
307  //int w = srclr.x - srcul.x;
308  int w = shape[0];
309  //int h = srclr.y - srcul.y;
310  int h = shape[1];
311  //int d = srclr.z - srcul.z;
312  int d = shape[2];
313  int count = 0;
314 
315  SrcImageIterator isy = srcul, isx = srcul, isz = srcul; // iterators for the src image
316 
317  typedef typename RegionStatisticsArray::value_type RegionStatistics;
318  typedef typename PromoteTraits<typename RegionStatistics::cost_type, double>::Promote CostType;
319  typedef detail::SeedRgVoxel<CostType, Diff_type> Voxel;
320 
321  typename Voxel::Allocator allocator;
322 
323  typedef std::priority_queue< Voxel *,
324  std::vector<Voxel *>,
325  typename Voxel::Compare > SeedRgVoxelHeap;
326  typedef MultiArray<3, int> IVolume;
327 
328  // copy seed image in an image with border
329  Diff_type regionshape = shape + Diff_type(2,2,2);
330  IVolume regions(regionshape);
331  MultiIterator<3,int> ir = regions.traverser_begin();
332  ir = ir + Diff_type(1,1,1);
333 
334  //IVolume::Iterator iry, irx, irz;
335  MultiIterator<3,int> iry, irx, irz;
336 
337  //initImageBorder(destImageRange(regions), 1, SRGWatershedLabel);
338  initMultiArrayBorder(destMultiArrayRange(regions), 1, SRGWatershedLabel);
339 
340  copyMultiArray(seedsul, Diff_type(w,h,d), aseeds, ir, AccessorTraits<int>::default_accessor());
341 
342  // allocate and init memory for the results
343 
344  SeedRgVoxelHeap pheap;
345  int cneighbor;
346 
347  #if 0
348  static const Diff_type dist[] = { Diff_type(-1, 0, 0), Diff_type( 0,-1, 0),
349  Diff_type( 1, 0, 0), Diff_type( 0, 1, 0),
350  Diff_type( 0, 0,-1), Diff_type( 0, 0, 1) };
351  #endif
352 
353  typedef typename Neighborhood::Direction Direction;
354  int directionCount = Neighborhood::DirectionCount;
355 
356  Diff_type pos(0,0,0);
357 
358  for(isz=srcul, irz=ir, pos[2]=0; pos[2]<d;
359  pos[2]++, isz.dim2()++, irz.dim2()++)
360  {
361  //std::cerr << "Z = " << pos[2] << std::endl;
362 
363  for(isy=isz, iry=irz, pos[1]=0; pos[1]<h;
364  pos[1]++, isy.dim1()++, iry.dim1()++)
365  {
366  //std::cerr << "Y = " << pos[1] << std::endl;
367 
368  for(isx=isy, irx=iry, pos[0]=0; pos[0]<w;
369  pos[0]++, isx.dim0()++, irx.dim0()++)
370  {
371  //std::cerr << "X = " << pos[0] << std::endl;
372 
373  if(*irx == 0)
374  {
375  // find candidate pixels for growing and fill heap
376  for(int i=0; i<directionCount; i++)
377  {
378  cneighbor = *(irx + Neighborhood::diff((Direction)i));
379  if(cneighbor > 0)
380  {
381  CostType cost = stats[cneighbor].cost(as(isx));
382 
383  Voxel * voxel =
384  allocator.create(pos, pos+Neighborhood::diff((Direction)i), cost, count++, cneighbor);
385  pheap.push(voxel);
386  }
387  }
388  }
389  }
390  }
391  }
392 
393  // perform region growing
394  while(pheap.size() != 0)
395  {
396  Voxel * voxel = pheap.top();
397  pheap.pop();
398 
399  Diff_type pos = voxel->location_;
400  Diff_type nearest = voxel->nearest_;
401  int lab = voxel->label_;
402  CostType cost = voxel->cost_;
403 
404  allocator.dismiss(voxel);
405 
406  if((srgType & StopAtThreshold) != 0 && cost > max_cost)
407  break;
408 
409  irx = ir + pos;
410  isx = srcul + pos;
411 
412  if(*irx) // already labelled region / watershed?
413  continue;
414 
415  if((srgType & KeepContours) != 0)
416  {
417  for(int i=0; i<directionCount; i++)
418  {
419  cneighbor = * (irx + Neighborhood::diff((Direction)i));
420  if((cneighbor>0) && (cneighbor != lab))
421  {
422  lab = SRGWatershedLabel;
423  break;
424  }
425  }
426  }
427 
428  *irx = lab;
429 
430  if((srgType & KeepContours) == 0 || lab > 0)
431  {
432  // update statistics
433  stats[*irx](as(isx));
434 
435  // search neighborhood
436  // second pass: find new candidate pixels
437  for(int i=0; i<directionCount; i++)
438  {
439  if(*(irx + Neighborhood::diff((Direction)i)) == 0)
440  {
441  CostType cost = stats[lab].cost(as(isx, Neighborhood::diff((Direction)i)));
442 
443  Voxel * new_voxel =
444  allocator.create(pos+Neighborhood::diff((Direction)i), nearest, cost, count++, lab);
445  pheap.push(new_voxel);
446  }
447  }
448  }
449  }
450 
451  // free temporary memory
452  while(pheap.size() != 0)
453  {
454  allocator.dismiss(pheap.top());
455  pheap.pop();
456  }
457 
458  // write result
459  transformMultiArray(ir, Diff_type(w,h,d), AccessorTraits<int>::default_accessor(),
460  destul, ad, detail::UnlabelWatersheds());
461 }
462 
463 template <class SrcImageIterator, class Diff_type, class SrcAccessor,
464  class SeedImageIterator, class SeedAccessor,
465  class DestImageIterator, class DestAccessor,
466  class RegionStatisticsArray, class Neighborhood >
467 inline void
468 seededRegionGrowing3D(SrcImageIterator srcul, Diff_type shape, SrcAccessor as,
469  SeedImageIterator seedsul, SeedAccessor aseeds,
470  DestImageIterator destul, DestAccessor ad,
471  RegionStatisticsArray & stats, SRGType srgType, Neighborhood n)
472 {
473  seededRegionGrowing3D( srcul, shape, as, seedsul, aseeds,
474  destul, ad, stats, srgType, n, NumericTraits<double>::max());
475 }
476 
477 template <class SrcImageIterator, class Diff_type, class SrcAccessor,
478  class SeedImageIterator, class SeedAccessor,
479  class DestImageIterator, class DestAccessor,
480  class RegionStatisticsArray >
481 inline void
482 seededRegionGrowing3D(SrcImageIterator srcul, Diff_type shape, SrcAccessor as,
483  SeedImageIterator seedsul, SeedAccessor aseeds,
484  DestImageIterator destul, DestAccessor ad,
485  RegionStatisticsArray & stats, SRGType srgType)
486 {
487  seededRegionGrowing3D( srcul, shape, as, seedsul, aseeds,
488  destul, ad, stats, srgType, NeighborCode3DSix());
489 }
490 
491 template <class SrcImageIterator, class Diff_type, class SrcAccessor,
492  class SeedImageIterator, class SeedAccessor,
493  class DestImageIterator, class DestAccessor,
494  class RegionStatisticsArray >
495 inline void
496 seededRegionGrowing3D(SrcImageIterator srcul, Diff_type shape, SrcAccessor as,
497  SeedImageIterator seedsul, SeedAccessor aseeds,
498  DestImageIterator destul, DestAccessor ad,
499  RegionStatisticsArray & stats)
500 {
501  seededRegionGrowing3D( srcul, shape, as, seedsul, aseeds, destul, ad,
502  stats, CompleteGrow);
503 }
504 
505 template <class SrcImageIterator, class Shape, class SrcAccessor,
506  class SeedImageIterator, class SeedAccessor,
507  class DestImageIterator, class DestAccessor,
508  class RegionStatisticsArray, class Neighborhood>
509 inline void
510 seededRegionGrowing3D(triple<SrcImageIterator, Shape, SrcAccessor> img1,
511  pair<SeedImageIterator, SeedAccessor> img3,
512  pair<DestImageIterator, DestAccessor> img4,
513  RegionStatisticsArray & stats,
514  SRGType srgType, Neighborhood n, double max_cost)
515 {
516  seededRegionGrowing3D(img1.first, img1.second, img1.third,
517  img3.first, img3.second,
518  img4.first, img4.second,
519  stats, srgType, n, max_cost);
520 }
521 
522 template <class SrcImageIterator, class Shape, class SrcAccessor,
523  class SeedImageIterator, class SeedAccessor,
524  class DestImageIterator, class DestAccessor,
525  class RegionStatisticsArray, class Neighborhood>
526 inline void
527 seededRegionGrowing3D(triple<SrcImageIterator, Shape, SrcAccessor> img1,
528  pair<SeedImageIterator, SeedAccessor> img3,
529  pair<DestImageIterator, DestAccessor> img4,
530  RegionStatisticsArray & stats,
531  SRGType srgType, Neighborhood n)
532 {
533  seededRegionGrowing3D(img1.first, img1.second, img1.third,
534  img3.first, img3.second,
535  img4.first, img4.second,
536  stats, srgType, n, NumericTraits<double>::max());
537 }
538 
539 template <class SrcImageIterator, class Shape, class SrcAccessor,
540  class SeedImageIterator, class SeedAccessor,
541  class DestImageIterator, class DestAccessor,
542  class RegionStatisticsArray>
543 inline void
544 seededRegionGrowing3D(triple<SrcImageIterator, Shape, SrcAccessor> img1,
545  pair<SeedImageIterator, SeedAccessor> img3,
546  pair<DestImageIterator, DestAccessor> img4,
547  RegionStatisticsArray & stats, SRGType srgType)
548 {
549  seededRegionGrowing3D(img1.first, img1.second, img1.third,
550  img3.first, img3.second,
551  img4.first, img4.second,
552  stats, srgType, NeighborCode3DSix());
553 }
554 
555 template <class SrcImageIterator, class Shape, class SrcAccessor,
556  class SeedImageIterator, class SeedAccessor,
557  class DestImageIterator, class DestAccessor,
558  class RegionStatisticsArray>
559 inline void
560 seededRegionGrowing3D(triple<SrcImageIterator, Shape, SrcAccessor> img1,
561  pair<SeedImageIterator, SeedAccessor> img3,
562  pair<DestImageIterator, DestAccessor> img4,
563  RegionStatisticsArray & stats)
564 {
565  seededRegionGrowing3D(img1.first, img1.second, img1.third,
566  img3.first, img3.second,
567  img4.first, img4.second,
568  stats);
569 }
570 
571 } // namespace vigra
572 
573 #endif // VIGRA_SEEDEDREGIONGROWING_HXX
574 

© 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.7.1 (Wed Mar 12 2014)