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

slic.hxx VIGRA

1 /************************************************************************/
2 /* */
3 /* Copyright 2012-2013 by Ullrich Koethe and Thorsten Beier */
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_SLIC_HXX
37 #define VIGRA_SLIC_HXX
38 
39 #include "multi_array.hxx"
40 #include "multi_convolution.hxx"
41 #include "multi_labeling.hxx"
42 #include "numerictraits.hxx"
43 #include "accumulator.hxx"
44 #include "array_vector.hxx"
45 
46 namespace vigra {
47 
48 /** \addtogroup SeededRegionGrowing
49 */
50 //@{
51 
52 /********************************************************/
53 /* */
54 /* generateSlicSeeds */
55 /* */
56 /********************************************************/
57 
58 /** \brief Generate seeds for SLIC superpixel computation in arbitrary dimensions.
59 
60  The source array \a src must be a scalar boundary indicator such as the gradient
61  magnitude. Seeds are initially placed on a regular Cartesian grid with spacing
62  \a seedDist und then moved to the point with smallest boundary indicator within
63  a search region of radius \a searchRadius around the initial position. The resulting
64  points are then marked in the output array \a seeds by consecutive labels.
65 
66  The function returns the number of selected seeds, which equals the largest seed label
67  because labeling starts at 1.
68 
69  <b> Declaration:</b>
70 
71  use arbitrary-dimensional arrays:
72  \code
73  namespace vigra {
74  template <unsigned int N, class T, class S1,
75  class Label, class S2>
76  unsigned int
77  generateSlicSeeds(MultiArrayView<N, T, S1> const & src,
78  MultiArrayView<N, Label, S2> seeds,
79  unsigned int seedDist,
80  unsigned int searchRadius = 1);
81  }
82  \endcode
83 
84  <b> Usage:</b>
85 
86  <b>\#include</b> <vigra/slic.hxx><br>
87  Namespace: vigra
88 
89  \code
90  MultiArray<2, RGBValue<float> > src(Shape2(w, h));
91  ... // fill src image
92 
93  // transform image to Lab color space
94  transformImage(srcImageRange(src), destImage(src), RGBPrime2LabFunctor<float>());
95 
96  // compute image gradient magnitude at scale 1.0 as a boundary indicator
97  MultiArray<2, float> grad(src.shape());
98  gaussianGradientMagnitude(srcImageRange(src), destImage(grad), 1.0);
99 
100  MultiArray<2, unsigned int> seeds(src.shape());
101  int seedDistance = 15;
102 
103  // place seeds on a grid with distance 15, but then move it to the lowest gradient
104  // poistion in a 3x3 window
105  generateSlicSeeds(grad, seeds, seedDistance);
106  \endcode
107 
108  For more details and examples see slicSuperpixels().
109 */
110 doxygen_overloaded_function(template <...> unsigned int generateSlicSeeds)
111 
112 template <unsigned int N, class T, class S1,
113  class Label, class S2>
114 unsigned int
115 generateSlicSeeds(MultiArrayView<N, T, S1> const & boundaryIndicatorImage,
116  MultiArrayView<N, Label, S2> seeds,
117  unsigned int seedDist,
118  unsigned int searchRadius = 1)
119 {
120  typedef typename MultiArrayShape<N>::type Shape;
121 
122  seeds.init(0);
123  Shape shape(boundaryIndicatorImage.shape()),
124  seedShape(floor(shape / double(seedDist))),
125  offset((shape - (seedShape - Shape(1))*seedDist) / 2);
126 
127  unsigned int label = 0;
128  MultiCoordinateIterator<N> iter(seedShape),
129  end = iter.getEndIterator();
130  for(; iter != end; ++iter)
131  {
132  // define search window around current seed center
133  Shape center = (*iter)*seedDist + offset;
134  Shape startCoord = max(Shape(0), center-Shape(searchRadius));
135  Shape endCoord = min(center+Shape(searchRadius+1), shape);
136 
137  // find the coordinate of minimum boundary indicator in window
138  using namespace acc;
139  AccumulatorChain<CoupledArrays<N, T>,
140  Select<WeightArg<1>, Coord<ArgMinWeight> > > a;
141  extractFeatures(boundaryIndicatorImage.subarray(startCoord, endCoord), a);
142 
143  // add seed at minimum position, if not already occupied
144  Shape minCoord = get<Coord<ArgMinWeight> >(a) + startCoord;
145  if(seeds[minCoord] == 0)
146  seeds[minCoord] = ++label;
147  }
148  return label;
149 }
150 
151 /** \brief Options object for slicSuperpixels().
152 
153  <b> Usage:</b>
154 
155  see slicSuperpixels() for detailed examples.
156 */
158 {
159  /** \brief Create options object with default settings.
160 
161  Defaults are: perform 10 iterations, determine a size limit for superpixels automatically.
162  */
164  : iter(10),
165  sizeLimit(0)
166  {}
167 
168  /** \brief Number of iterations.
169 
170  Default: 10
171  */
172  SlicOptions & iterations(unsigned int i)
173  {
174  iter = i;
175  return *this;
176  }
177 
178  /** \brief Minimum superpixel size.
179 
180  If you set this to 1, no size filtering will be performed.
181 
182  Default: 0 (determine size limit automatically as <tt>average size / 4</tt>)
183  */
184  SlicOptions & minSize(unsigned int s)
185  {
186  sizeLimit = s;
187  return *this;
188  }
189 
190  unsigned int iter;
191  unsigned int sizeLimit;
192 };
193 
194 namespace detail {
195 
196 template <unsigned int N, class T, class Label>
197 class Slic
198 {
199  public:
200  //
201  typedef MultiArrayView<N, T> DataImageType;
202  typedef MultiArrayView<N, Label> LabelImageType;
203  typedef typename DataImageType::difference_type ShapeType;
204  typedef typename PromoteTraits<
205  typename NormTraits<T>::NormType,
206  typename NormTraits<MultiArrayIndex>::NormType
207  >::Promote DistanceType;
208 
209  Slic(DataImageType dataImage,
210  LabelImageType labelImage,
211  DistanceType intensityScaling,
212  int maxRadius,
213  SlicOptions const & options = SlicOptions());
214 
215  unsigned int execute();
216 
217  private:
218  void updateAssigments();
219  unsigned int postProcessing();
220 
221  typedef MultiArray<N,DistanceType> DistanceImageType;
222 
223  ShapeType shape_;
224  DataImageType dataImage_;
225  LabelImageType labelImage_;
226  DistanceImageType distance_;
227  int max_radius_;
228  DistanceType normalization_;
229  SlicOptions options_;
230 
231  typedef acc::Select<acc::DataArg<1>, acc::LabelArg<2>, acc::Mean, acc::RegionCenter> Statistics;
232  typedef acc::AccumulatorChainArray<CoupledArrays<N, T, Label>, Statistics> RegionFeatures;
233  RegionFeatures clusters_;
234 };
235 
236 
237 
238 template <unsigned int N, class T, class Label>
239 Slic<N, T, Label>::Slic(
240  DataImageType dataImage,
241  LabelImageType labelImage,
242  DistanceType intensityScaling,
243  int maxRadius,
244  SlicOptions const & options)
245 : shape_(dataImage.shape()),
246  dataImage_(dataImage),
247  labelImage_(labelImage),
248  distance_(shape_),
249  max_radius_(maxRadius),
250  normalization_(sq(intensityScaling) / sq(max_radius_)),
251  options_(options)
252 {
253  clusters_.ignoreLabel(0);
254 }
255 
256 template <unsigned int N, class T, class Label>
257 unsigned int Slic<N, T, Label>::execute()
258 {
259  // Do SLIC
260  for(size_t i=0; i<options_.iter; ++i)
261  {
262  // update mean for each cluster
263  clusters_.reset();
264  extractFeatures(dataImage_, labelImage_, clusters_);
265 
266  // update which pixels get assigned to which cluster
267  updateAssigments();
268  }
269 
270  return postProcessing();
271 }
272 
273 template <unsigned int N, class T, class Label>
274 void
275 Slic<N, T, Label>::updateAssigments()
276 {
277  using namespace acc;
278  distance_.init(NumericTraits<DistanceType>::max());
279  for(unsigned int c=1; c<=clusters_.maxRegionLabel(); ++c)
280  {
281  if(get<Count>(clusters_, c) == 0) // label doesn't exist
282  continue;
283 
284  typedef typename LookupTag<RegionCenter, RegionFeatures>::value_type CenterType;
285  CenterType center = get<RegionCenter>(clusters_, c);
286 
287  // get ROI limits around region center
288  ShapeType pixelCenter(round(center)),
289  startCoord(max(ShapeType(0), pixelCenter - ShapeType(max_radius_))),
290  endCoord(min(shape_, pixelCenter + ShapeType(max_radius_+1)));
291  center -= startCoord; // need center relative to ROI
292 
293  // setup iterators for ROI
294  typedef typename CoupledArrays<N, T, Label, DistanceType>::IteratorType Iterator;
295  Iterator iter = createCoupledIterator(dataImage_, labelImage_, distance_).
296  restrictToSubarray(startCoord, endCoord),
297  end = iter.getEndIterator();
298 
299  // only pixels within the ROI can be assigned to a cluster
300  for(; iter != end; ++iter)
301  {
302  // compute distance between cluster center and pixel
303  DistanceType spatialDist = squaredNorm(center-iter.point());
304  DistanceType colorDist = squaredNorm(get<Mean>(clusters_, c)-iter.template get<1>());
305  DistanceType dist = colorDist + normalization_*spatialDist;
306  // update label?
307  if(dist < iter.template get<3>())
308  {
309  iter.template get<2>() = static_cast<Label>(c);
310  iter.template get<3>() = dist;
311  }
312  }
313  }
314 }
315 
316 template <unsigned int N, class T, class Label>
317 unsigned int
318 Slic<N, T, Label>::postProcessing()
319 {
320  // get rid of regions below a size limit
321  MultiArray<N,Label> tmpLabelImage(labelImage_);
322  unsigned int maxLabel = labelMultiArray(tmpLabelImage, labelImage_, DirectNeighborhood);
323 
324  unsigned int sizeLimit = options_.sizeLimit == 0
325  ? (unsigned int)(0.25 * labelImage_.size() / maxLabel)
326  : options_.sizeLimit;
327  if(sizeLimit == 1)
328  return maxLabel;
329 
330  // determine region size
331  using namespace acc;
332  AccumulatorChainArray<CoupledArrays<N, Label>, Select<LabelArg<1>, Count> > sizes;
333  extractFeatures(labelImage_, sizes);
334 
335  typedef GridGraph<N, undirected_tag> Graph;
336  Graph graph(labelImage_.shape(), DirectNeighborhood);
337 
338  typedef typename Graph::NodeIt graph_scanner;
339  typedef typename Graph::OutBackArcIt neighbor_iterator;
340 
341  ArrayVector<Label> regions(maxLabel+1);
342 
343  // make sure that all regions exceed the sizeLimit
344  for (graph_scanner node(graph); node != lemon::INVALID; ++node)
345  {
346  Label label = labelImage_[*node];
347 
348  if(regions[label] > 0)
349  continue; // already processed
350 
351  regions[label] = label;
352 
353  if(get<Count>(sizes, label) < sizeLimit)
354  {
355  // region is too small => merge into an existing neighbor
356  for (neighbor_iterator arc(graph, node); arc != lemon::INVALID; ++arc)
357  {
358  regions[label] = regions[labelImage_[graph.target(*arc)]];
359  break;
360  }
361  }
362  }
363 
364  // make labels contiguous after possible merging
365  maxLabel = 0;
366  for(unsigned int i=1; i<=maxLabel; ++i)
367  {
368  if(regions[i] == i)
369  {
370  regions[i] = (Label)++maxLabel;
371  }
372  else
373  {
374  regions[i] = regions[regions[i]];
375  }
376  }
377 
378  // update labels
379  for (graph_scanner node(graph); node != lemon::INVALID; ++node)
380  {
381  labelImage_[*node] = regions[labelImage_[*node]];
382  }
383 
384  return maxLabel;
385 }
386 
387 } // namespace detail
388 
389 
390 /** \brief Compute SLIC superpixels in arbitrary dimensions.
391 
392  This function implements the algorithm described in
393 
394  R. Achanta et al.: <em>"SLIC Superpixels Compared to State-of-the-Art
395  Superpixel Methods"</em>, IEEE Trans. Patt. Analysis Mach. Intell. 34(11):2274-2281, 2012
396 
397  The value type <tt>T</tt> of the source array \a src must provide the necessary functionality
398  to compute averages and squared distances (i.e. it must fulfill the requirements of a
399  \ref LinearSpace and support squaredNorm(T const &)). This is true for all scalar types as well as
400  \ref vigra::TinyVector and \ref vigra::RGBValue. The output array \a labels will be filled
401  with labels designating membership of each point in one of the superpixel regions.
402 
403  The output array can optionally contain seeds (which will be overwritten by the output)
404  to give you full control over seed placement. If \a labels is empty, seeds will be created
405  automatically by an internal call to generateSlicSeeds().
406 
407  The parameter \a seedDistance specifies the radius of the window around each seed (or, more
408  precisely, around the present regions centers) where the algorithm looks for potential members
409  of the corresponding superpixel. It thus places an upper limit on the superpixel size. When seeds
410  are computed automatically, this parameter also determines the grid spacing for seed placement.
411 
412  The parameter \a intensityScaling is used to normalize (i.e. divide) the color/intensity difference
413  before it is compared with the spatial distance. This corresponds to parameter <i>m</i> in equation
414  (2) of the paper.
415 
416  The options object can be used to specify the number of iterations (<tt>SlicOptions::iterations()</tt>)
417  and an explicit minimal superpixel size (<tt>SlicOptions::minSize()</tt>). By default, the algorithm
418  merges all regions that are smaller than 1/4 of the average superpixel size.
419 
420  The function returns the number of superpixels, which equals the largest label
421  because labeling starts at 1.
422 
423  <b> Declaration:</b>
424 
425  use arbitrary-dimensional arrays:
426  \code
427  namespace vigra {
428  template <unsigned int N, class T, class S1,
429  class Label, class S2,
430  class DistanceType>
431  unsigned int
432  slicSuperpixels(MultiArrayView<N, T, S1> const & src,
433  MultiArrayView<N, Label, S2> labels,
434  DistanceType intensityScaling,
435  unsigned int seedDistance,
436  SlicOptions const & options = SlicOptions());
437  }
438  \endcode
439 
440  <b> Usage:</b>
441 
442  <b>\#include</b> <vigra/slic.hxx><br>
443  Namespace: vigra
444 
445  \code
446  MultiArray<2, RGBValue<float> > src(Shape2(w, h));
447  ... // fill src image
448 
449  // transform image to Lab color space
450  transformMultiArray(srcMultiArrayRange(src), destMultiArray(src), RGBPrime2LabFunctor<float>());
451 
452  MultiArray<2, unsigned int> labels(src.shape());
453  int seedDistance = 15;
454  double intensityScaling = 20.0;
455 
456  // compute seeds automatically, perform 40 iterations, and scale intensity differences
457  // down to 1/20 before comparing with spatial distances
458  slicSuperpixels(src, labels, intensityScaling, seedDistance, SlicOptions().iterations(40));
459  \endcode
460 
461  This works for arbitrary-dimensional arrays.
462 */
463 doxygen_overloaded_function(template <...> unsigned int slicSuperpixels)
464 
465 template <unsigned int N, class T, class S1,
466  class Label, class S2,
467  class DistanceType>
468 unsigned int
469 slicSuperpixels(MultiArrayView<N, T, S1> const & src,
470  MultiArrayView<N, Label, S2> labels,
471  DistanceType intensityScaling,
472  unsigned int seedDistance,
473  SlicOptions const & options = SlicOptions())
474 {
475  if(!labels.any())
476  {
477  typedef typename NormTraits<T>::NormType TmpType;
478  MultiArray<N, TmpType> grad(src.shape());
479  gaussianGradientMagnitude(src, grad, 1.0);
480  generateSlicSeeds(grad, labels, seedDistance);
481  }
482  return detail::Slic<N, T, Label>(src, labels, intensityScaling, seedDistance, options).execute();
483 }
484 
485 //@}
486 
487 } // namespace vigra
488 
489 #endif // VIGRA_SLIC_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 (Mon Nov 18 2013)