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 Superpixels
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 */
doxygen_overloaded_function(template<...> unsigned int generateSlicSeeds)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 */
157 struct SlicOptions
158 {
159 /** \brief Create options object with default settings.
160
161 Defaults are: perform 10 iterations, determine a size limit for superpixels automatically.
162 */
SlicOptionsvigra::SlicOptions163 SlicOptions()
164 : iter(10),
165 sizeLimit(0)
166 {}
167
168 /** \brief Number of iterations.
169
170 Default: 10
171 */
iterationsvigra::SlicOptions172 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 */
minSizevigra::SlicOptions184 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>
Slic(DataImageType dataImage,LabelImageType labelImage,DistanceType intensityScaling,int maxRadius,SlicOptions const & options)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>
execute()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
updateAssigments()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
postProcessing()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 vigra::UnionFindArray<Label> regions(maxLabel+1);
342 ArrayVector<unsigned char> done(maxLabel+1, false);
343
344 // make sure that all regions exceed the sizeLimit
345 for (graph_scanner node(graph); node != lemon::INVALID; ++node)
346 {
347 Label label = labelImage_[*node];
348
349 if(done[label])
350 continue; // already processed
351
352 if(get<Count>(sizes, label) < sizeLimit)
353 {
354 // region is too small => merge into a neighbor
355 for (neighbor_iterator arc(graph, node); arc != lemon::INVALID; ++arc)
356 {
357 Label other = labelImage_[graph.target(*arc)];
358 if(label != other)
359 {
360 regions.makeUnion(label, other);
361 done[label] = true;
362 break;
363 }
364 }
365 }
366 else
367 {
368 done[label] = true;
369 }
370 }
371
372 // make labels contiguous after possible merging
373 Label newMaxLabel = regions.makeContiguous();
374 for (graph_scanner node(graph); node != lemon::INVALID; ++node)
375 {
376 labelImage_[*node] = regions.findLabel(labelImage_[*node]);
377 }
378
379 return (unsigned int)newMaxLabel;
380 }
381
382 } // namespace detail
383
384
385 /** \brief Compute SLIC superpixels in arbitrary dimensions.
386
387 This function implements the algorithm described in
388
389 R. Achanta et al.: <em>"SLIC Superpixels Compared to State-of-the-Art
390 Superpixel Methods"</em>, IEEE Trans. Patt. Analysis Mach. Intell. 34(11):2274-2281, 2012
391
392 The value type <tt>T</tt> of the source array \a src must provide the necessary functionality
393 to compute averages and squared distances (i.e. it must fulfill the requirements of a
394 \ref LinearSpace and support squaredNorm(T const &)). This is true for all scalar types as well as
395 \ref vigra::TinyVector and \ref vigra::RGBValue. The output array \a labels will be filled
396 with labels designating membership of each point in one of the superpixel regions.
397
398 The output array can optionally contain seeds (which will be overwritten by the output)
399 to give you full control over seed placement. If \a labels is empty, seeds will be created
400 automatically by an internal call to generateSlicSeeds().
401
402 The parameter \a seedDistance specifies the radius of the window around each seed (or, more
403 precisely, around the present regions centers) where the algorithm looks for potential members
404 of the corresponding superpixel. It thus places an upper limit on the superpixel size. When seeds
405 are computed automatically, this parameter also determines the grid spacing for seed placement.
406
407 The parameter \a intensityScaling is used to normalize (i.e. divide) the color/intensity difference
408 before it is compared with the spatial distance. This corresponds to parameter <i>m</i> in equation
409 (2) of the paper.
410
411 The options object can be used to specify the number of iterations (<tt>SlicOptions::iterations()</tt>)
412 and an explicit minimal superpixel size (<tt>SlicOptions::minSize()</tt>). By default, the algorithm
413 merges all regions that are smaller than 1/4 of the average superpixel size.
414
415 The function returns the number of superpixels, which equals the largest label
416 because labeling starts at 1.
417
418 <b> Declaration:</b>
419
420 use arbitrary-dimensional arrays:
421 \code
422 namespace vigra {
423 template <unsigned int N, class T, class S1,
424 class Label, class S2,
425 class DistanceType>
426 unsigned int
427 slicSuperpixels(MultiArrayView<N, T, S1> const & src,
428 MultiArrayView<N, Label, S2> labels,
429 DistanceType intensityScaling,
430 unsigned int seedDistance,
431 SlicOptions const & options = SlicOptions());
432 }
433 \endcode
434
435 <b> Usage:</b>
436
437 <b>\#include</b> \<vigra/slic.hxx\><br>
438 Namespace: vigra
439
440 \code
441 MultiArray<2, RGBValue<float> > src(Shape2(w, h));
442 ... // fill src image
443
444 // transform image to Lab color space
445 transformMultiArray(srcMultiArrayRange(src), destMultiArray(src), RGBPrime2LabFunctor<float>());
446
447 MultiArray<2, unsigned int> labels(src.shape());
448 int seedDistance = 15;
449 double intensityScaling = 20.0;
450
451 // compute seeds automatically, perform 40 iterations, and scale intensity differences
452 // down to 1/20 before comparing with spatial distances
453 slicSuperpixels(src, labels, intensityScaling, seedDistance, SlicOptions().iterations(40));
454 \endcode
455
456 This works for arbitrary-dimensional arrays.
457 */
doxygen_overloaded_function(template<...> unsigned int slicSuperpixels)458 doxygen_overloaded_function(template <...> unsigned int slicSuperpixels)
459
460 template <unsigned int N, class T, class S1,
461 class Label, class S2,
462 class DistanceType>
463 unsigned int
464 slicSuperpixels(MultiArrayView<N, T, S1> const & src,
465 MultiArrayView<N, Label, S2> labels,
466 DistanceType intensityScaling,
467 unsigned int seedDistance,
468 SlicOptions const & options = SlicOptions())
469 {
470 if(!labels.any())
471 {
472 typedef typename NormTraits<T>::NormType TmpType;
473 MultiArray<N, TmpType> grad(src.shape());
474 gaussianGradientMagnitude(src, grad, 1.0);
475 generateSlicSeeds(grad, labels, seedDistance);
476 }
477 return detail::Slic<N, T, Label>(src, labels, intensityScaling, seedDistance, options).execute();
478 }
479
480 //@}
481
482 } // namespace vigra
483
484 #endif // VIGRA_SLIC_HXX
485