[ 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 */
167 //@{
168 
169 /********************************************************/
170 /* */
171 /* seededRegionGrowing3D */
172 /* */
173 /********************************************************/
174 
175 /** \brief Three-dimensional Region Segmentation by means of Seeded Region Growing.
176 
177  This algorithm implements seeded region growing as described in
178 
179  The seed image is a partly segmented multi-dimensional array which contains uniquely
180  labeled regions (the seeds) and unlabeled voxels (the candidates, label 0).
181  Seed regions can be as large as you wish and as small as one voxel. If
182  there are no candidates, the algorithm will simply copy the seed array
183  into the output array. Otherwise it will aggregate the candidates into
184  the existing regions so that a cost function is minimized.
185  Candidates are taken from the neighborhood of the already assigned pixels,
186  where the type of neighborhood is determined by parameter <tt>neighborhood</tt>
187  which can take the values <tt>NeighborCode3DSix()</tt> (the default)
188  or <tt>NeighborCode3DTwentySix()</tt>. The algorithm basically works as follows
189  (illustrated for 6-neighborhood, but 26-neighborhood works in the same way):
190 
191  <ol>
192 
193  <li> Find all candidate pixels that are 6-adjacent to a seed region.
194  Calculate the cost for aggregating each candidate into its adjacent region
195  and put the candidates into a priority queue.
196 
197  <li> While( priority queue is not empty)
198 
199  <ol>
200 
201  <li> Take the candidate with least cost from the queue. If it has not
202  already been merged, merge it with it's adjacent region.
203 
204  <li> Put all candidates that are 4-adjacent to the pixel just processed
205  into the priority queue.
206 
207  </ol>
208 
209  </ol>
210 
211  <tt>SRGType</tt> can take the following values:
212 
213  <DL>
214  <DT><tt>CompleteGrow</tt> <DD> produce a complete tesselation of the volume (default).
215  <DT><tt>KeepContours</tt> <DD> keep a 1-voxel wide unlabeled contour between all regions.
216  <DT><tt>StopAtThreshold</tt> <DD> stop when the boundary indicator values exceed the
217  threshold given by parameter <tt>max_cost</tt>.
218  <DT><tt>KeepContours | StopAtThreshold</tt> <DD> keep 1-voxel wide contour and stop at given <tt>max_cost</tt>.
219  </DL>
220 
221  The cost is determined jointly by the source array and the
222  region statistics functor. The source array contains feature values for each
223  pixel which will be used by the region statistics functor to calculate and
224  update statistics for each region and to calculate the cost for each
225  candidate. The <TT>RegionStatisticsArray</TT> must be compatible to the
226  \ref ArrayOfRegionStatistics functor and contains an <em> array</em> of
227  statistics objects for each region. The indices must correspond to the
228  labels of the seed regions. The statistics for the initial regions must have
229  been calculated prior to calling <TT>seededRegionGrowing3D()</TT>
230 
231  For each candidate
232  <TT>x</TT> that is adjacent to region <TT>i</TT>, the algorithm will call
233  <TT>stats[i].cost(as(x))</TT> to get the cost (where <TT>x</TT> is a <TT>SrcImageIterator</TT>
234  and <TT>as</TT> is
235  the SrcAccessor). When a candidate has been merged with a region, the
236  statistics are updated by calling <TT>stats[i].operator()(as(x))</TT>. Since
237  the <TT>RegionStatisticsArray</TT> is passed by reference, this will overwrite
238  the original statistics.
239 
240  If a candidate could be merged into more than one regions with identical
241  cost, the algorithm will favour the nearest region. If <tt>StopAtThreshold</tt> is active,
242  and the cost of the current candidate at any point in the algorithm exceeds the optional
243  <tt>max_cost</tt> value (which defaults to <tt>NumericTraits<double>::max()</tt>),
244  region growing is aborted, and all voxels not yet assigned to a region remain unlabeled.
245 
246  In some cases, the cost only depends on the feature value of the current
247  voxel. Then the update operation will simply be a no-op, and the <TT>cost()</TT>
248  function returns its argument. This behavior is implemented by the
249  \ref SeedRgDirectValueFunctor.
250 
251  <b> Declarations:</b>
252 
253  pass arguments explicitly:
254  \code
255  namespace vigra {
256  template <class SrcImageIterator, class Shape, class SrcAccessor,
257  class SeedImageIterator, class SeedAccessor,
258  class DestImageIterator, class DestAccessor,
259  class RegionStatisticsArray, class Neighborhood>
260  void
261  seededRegionGrowing3D(SrcImageIterator srcul, Shape shape, SrcAccessor as,
262  SeedImageIterator seedsul, SeedAccessor aseeds,
263  DestImageIterator destul, DestAccessor ad,
264  RegionStatisticsArray & stats,
265  SRGType srgType = CompleteGrow,
266  Neighborhood neighborhood = NeighborCode3DSix(),
267  double max_cost = NumericTraits<double>::max());
268  }
269  \endcode
270 
271  use argument objects in conjunction with \ref ArgumentObjectFactories :
272  \code
273  namespace vigra {
274  template <class SrcImageIterator, class Shape, class SrcAccessor,
275  class SeedImageIterator, class SeedAccessor,
276  class DestImageIterator, class DestAccessor,
277  class RegionStatisticsArray, class Neighborhood>
278  void
279  seededRegionGrowing3D(triple<SrcImageIterator, Shape, SrcAccessor> src,
280  pair<SeedImageIterator, SeedAccessor> seeds,
281  pair<DestImageIterator, DestAccessor> dest,
282  RegionStatisticsArray & stats,
283  SRGType srgType = CompleteGrow,
284  Neighborhood neighborhood = NeighborCode3DSix(),
285  double max_cost = NumericTraits<double>::max());
286  }
287  \endcode
288 
289 */
290 doxygen_overloaded_function(template <...> void seededRegionGrowing3D)
291 
292 template <class SrcImageIterator, class Diff_type, class SrcAccessor,
293  class SeedImageIterator, class SeedAccessor,
294  class DestImageIterator, class DestAccessor,
295  class RegionStatisticsArray, class Neighborhood>
296 void
297 seededRegionGrowing3D(SrcImageIterator srcul, Diff_type shape, SrcAccessor as,
298  SeedImageIterator seedsul, SeedAccessor aseeds,
299  DestImageIterator destul, DestAccessor ad,
300  RegionStatisticsArray & stats,
301  SRGType srgType,
302  Neighborhood,
303  double max_cost)
304 {
305  //SrcImageIterator srclr = srcul + shape;
306  //int w = srclr.x - srcul.x;
307  int w = shape[0];
308  //int h = srclr.y - srcul.y;
309  int h = shape[1];
310  //int d = srclr.z - srcul.z;
311  int d = shape[2];
312  int count = 0;
313 
314  SrcImageIterator isy = srcul, isx = srcul, isz = srcul; // iterators for the src image
315 
316  typedef typename RegionStatisticsArray::value_type RegionStatistics;
317  typedef typename PromoteTraits<typename RegionStatistics::cost_type, double>::Promote CostType;
318  typedef detail::SeedRgVoxel<CostType, Diff_type> Voxel;
319 
320  typename Voxel::Allocator allocator;
321 
322  typedef std::priority_queue< Voxel *,
323  std::vector<Voxel *>,
324  typename Voxel::Compare > SeedRgVoxelHeap;
325  typedef MultiArray<3, int> IVolume;
326 
327  // copy seed image in an image with border
328  Diff_type regionshape = shape + Diff_type(2,2,2);
329  IVolume regions(regionshape);
330  MultiIterator<3,int> ir = regions.traverser_begin();
331  ir = ir + Diff_type(1,1,1);
332 
333  //IVolume::Iterator iry, irx, irz;
334  MultiIterator<3,int> iry, irx, irz;
335 
336  //initImageBorder(destImageRange(regions), 1, SRGWatershedLabel);
337  initMultiArrayBorder(destMultiArrayRange(regions), 1, SRGWatershedLabel);
338 
339  copyMultiArray(seedsul, Diff_type(w,h,d), aseeds, ir, AccessorTraits<int>::default_accessor());
340 
341  // allocate and init memory for the results
342 
343  SeedRgVoxelHeap pheap;
344  int cneighbor;
345 
346  #if 0
347  static const Diff_type dist[] = { Diff_type(-1, 0, 0), Diff_type( 0,-1, 0),
348  Diff_type( 1, 0, 0), Diff_type( 0, 1, 0),
349  Diff_type( 0, 0,-1), Diff_type( 0, 0, 1) };
350  #endif
351 
352  typedef typename Neighborhood::Direction Direction;
353  int directionCount = Neighborhood::DirectionCount;
354 
355  Diff_type pos(0,0,0);
356 
357  for(isz=srcul, irz=ir, pos[2]=0; pos[2]<d;
358  pos[2]++, isz.dim2()++, irz.dim2()++)
359  {
360  //std::cerr << "Z = " << pos[2] << std::endl;
361 
362  for(isy=isz, iry=irz, pos[1]=0; pos[1]<h;
363  pos[1]++, isy.dim1()++, iry.dim1()++)
364  {
365  //std::cerr << "Y = " << pos[1] << std::endl;
366 
367  for(isx=isy, irx=iry, pos[0]=0; pos[0]<w;
368  pos[0]++, isx.dim0()++, irx.dim0()++)
369  {
370  //std::cerr << "X = " << pos[0] << std::endl;
371 
372  if(*irx == 0)
373  {
374  // find candidate pixels for growing and fill heap
375  for(int i=0; i<directionCount; i++)
376  {
377  cneighbor = *(irx + Neighborhood::diff((Direction)i));
378  if(cneighbor > 0)
379  {
380  CostType cost = stats[cneighbor].cost(as(isx));
381 
382  Voxel * voxel =
383  allocator.create(pos, pos+Neighborhood::diff((Direction)i), cost, count++, cneighbor);
384  pheap.push(voxel);
385  }
386  }
387  }
388  }
389  }
390  }
391 
392  // perform region growing
393  while(pheap.size() != 0)
394  {
395  Voxel * voxel = pheap.top();
396  pheap.pop();
397 
398  Diff_type pos = voxel->location_;
399  Diff_type nearest = voxel->nearest_;
400  int lab = voxel->label_;
401  CostType cost = voxel->cost_;
402 
403  allocator.dismiss(voxel);
404 
405  if((srgType & StopAtThreshold) != 0 && cost > max_cost)
406  break;
407 
408  irx = ir + pos;
409  isx = srcul + pos;
410 
411  if(*irx) // already labelled region / watershed?
412  continue;
413 
414  if((srgType & KeepContours) != 0)
415  {
416  for(int i=0; i<directionCount; i++)
417  {
418  cneighbor = * (irx + Neighborhood::diff((Direction)i));
419  if((cneighbor>0) && (cneighbor != lab))
420  {
421  lab = SRGWatershedLabel;
422  break;
423  }
424  }
425  }
426 
427  *irx = lab;
428 
429  if((srgType & KeepContours) == 0 || lab > 0)
430  {
431  // update statistics
432  stats[*irx](as(isx));
433 
434  // search neighborhood
435  // second pass: find new candidate pixels
436  for(int i=0; i<directionCount; i++)
437  {
438  if(*(irx + Neighborhood::diff((Direction)i)) == 0)
439  {
440  CostType cost = stats[lab].cost(as(isx, Neighborhood::diff((Direction)i)));
441 
442  Voxel * new_voxel =
443  allocator.create(pos+Neighborhood::diff((Direction)i), nearest, cost, count++, lab);
444  pheap.push(new_voxel);
445  }
446  }
447  }
448  }
449 
450  // free temporary memory
451  while(pheap.size() != 0)
452  {
453  allocator.dismiss(pheap.top());
454  pheap.pop();
455  }
456 
457  // write result
458  transformMultiArray(ir, Diff_type(w,h,d), AccessorTraits<int>::default_accessor(),
459  destul, ad, detail::UnlabelWatersheds());
460 }
461 
462 template <class SrcImageIterator, class Diff_type, class SrcAccessor,
463  class SeedImageIterator, class SeedAccessor,
464  class DestImageIterator, class DestAccessor,
465  class RegionStatisticsArray, class Neighborhood >
466 inline void
467 seededRegionGrowing3D(SrcImageIterator srcul, Diff_type shape, SrcAccessor as,
468  SeedImageIterator seedsul, SeedAccessor aseeds,
469  DestImageIterator destul, DestAccessor ad,
470  RegionStatisticsArray & stats, SRGType srgType, Neighborhood n)
471 {
472  seededRegionGrowing3D( srcul, shape, as, seedsul, aseeds,
473  destul, ad, stats, srgType, n, NumericTraits<double>::max());
474 }
475 
476 template <class SrcImageIterator, class Diff_type, class SrcAccessor,
477  class SeedImageIterator, class SeedAccessor,
478  class DestImageIterator, class DestAccessor,
479  class RegionStatisticsArray >
480 inline void
481 seededRegionGrowing3D(SrcImageIterator srcul, Diff_type shape, SrcAccessor as,
482  SeedImageIterator seedsul, SeedAccessor aseeds,
483  DestImageIterator destul, DestAccessor ad,
484  RegionStatisticsArray & stats, SRGType srgType)
485 {
486  seededRegionGrowing3D( srcul, shape, as, seedsul, aseeds,
487  destul, ad, stats, srgType, NeighborCode3DSix());
488 }
489 
490 template <class SrcImageIterator, class Diff_type, class SrcAccessor,
491  class SeedImageIterator, class SeedAccessor,
492  class DestImageIterator, class DestAccessor,
493  class RegionStatisticsArray >
494 inline void
495 seededRegionGrowing3D(SrcImageIterator srcul, Diff_type shape, SrcAccessor as,
496  SeedImageIterator seedsul, SeedAccessor aseeds,
497  DestImageIterator destul, DestAccessor ad,
498  RegionStatisticsArray & stats)
499 {
500  seededRegionGrowing3D( srcul, shape, as, seedsul, aseeds, destul, ad,
501  stats, CompleteGrow);
502 }
503 
504 template <class SrcImageIterator, class Shape, class SrcAccessor,
505  class SeedImageIterator, class SeedAccessor,
506  class DestImageIterator, class DestAccessor,
507  class RegionStatisticsArray, class Neighborhood>
508 inline void
509 seededRegionGrowing3D(triple<SrcImageIterator, Shape, SrcAccessor> img1,
510  pair<SeedImageIterator, SeedAccessor> img3,
511  pair<DestImageIterator, DestAccessor> img4,
512  RegionStatisticsArray & stats,
513  SRGType srgType, Neighborhood n, double max_cost)
514 {
515  seededRegionGrowing3D(img1.first, img1.second, img1.third,
516  img3.first, img3.second,
517  img4.first, img4.second,
518  stats, srgType, n, max_cost);
519 }
520 
521 template <class SrcImageIterator, class Shape, class SrcAccessor,
522  class SeedImageIterator, class SeedAccessor,
523  class DestImageIterator, class DestAccessor,
524  class RegionStatisticsArray, class Neighborhood>
525 inline void
526 seededRegionGrowing3D(triple<SrcImageIterator, Shape, SrcAccessor> img1,
527  pair<SeedImageIterator, SeedAccessor> img3,
528  pair<DestImageIterator, DestAccessor> img4,
529  RegionStatisticsArray & stats,
530  SRGType srgType, Neighborhood n)
531 {
532  seededRegionGrowing3D(img1.first, img1.second, img1.third,
533  img3.first, img3.second,
534  img4.first, img4.second,
535  stats, srgType, n, NumericTraits<double>::max());
536 }
537 
538 template <class SrcImageIterator, class Shape, class SrcAccessor,
539  class SeedImageIterator, class SeedAccessor,
540  class DestImageIterator, class DestAccessor,
541  class RegionStatisticsArray>
542 inline void
543 seededRegionGrowing3D(triple<SrcImageIterator, Shape, SrcAccessor> img1,
544  pair<SeedImageIterator, SeedAccessor> img3,
545  pair<DestImageIterator, DestAccessor> img4,
546  RegionStatisticsArray & stats, SRGType srgType)
547 {
548  seededRegionGrowing3D(img1.first, img1.second, img1.third,
549  img3.first, img3.second,
550  img4.first, img4.second,
551  stats, srgType, NeighborCode3DSix());
552 }
553 
554 template <class SrcImageIterator, class Shape, class SrcAccessor,
555  class SeedImageIterator, class SeedAccessor,
556  class DestImageIterator, class DestAccessor,
557  class RegionStatisticsArray>
558 inline void
559 seededRegionGrowing3D(triple<SrcImageIterator, Shape, SrcAccessor> img1,
560  pair<SeedImageIterator, SeedAccessor> img3,
561  pair<DestImageIterator, DestAccessor> img4,
562  RegionStatisticsArray & stats)
563 {
564  seededRegionGrowing3D(img1.first, img1.second, img1.third,
565  img3.first, img3.second,
566  img4.first, img4.second,
567  stats);
568 }
569 
570 } // namespace vigra
571 
572 #endif // VIGRA_SEEDEDREGIONGROWING_HXX
573 

© 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.9.0 (Sat Oct 5 2013)