1 /*M///////////////////////////////////////////////////////////////////////////////////////
\r
3 // IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
\r
5 // By downloading, copying, installing or using the software you agree to this license.
\r
6 // If you do not agree to this license, do not download, install,
\r
7 // copy or use the software.
\r
10 // License Agreement
\r
11 // For Open Source Computer Vision Library
\r
13 // Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
\r
14 // Copyright (C) 2009, Willow Garage Inc., all rights reserved.
\r
15 // Third party copyrights are property of their respective owners.
\r
17 // Redistribution and use in source and binary forms, with or without modification,
\r
18 // are permitted provided that the following conditions are met:
\r
20 // * Redistribution's of source code must retain the above copyright notice,
\r
21 // this list of conditions and the following disclaimer.
\r
23 // * Redistribution's in binary form must reproduce the above copyright notice,
\r
24 // this list of conditions and the following disclaimer in the documentation
\r
25 // and/or other materials provided with the distribution.
\r
27 // * The name of the copyright holders may not be used to endorse or promote products
\r
28 // derived from this software without specific prior written permission.
\r
30 // This software is provided by the copyright holders and contributors "as is" and
\r
31 // any express or implied warranties, including, but not limited to, the implied
\r
32 // warranties of merchantability and fitness for a particular purpose are disclaimed.
\r
33 // In no event shall the Intel Corporation or contributors be liable for any direct,
\r
34 // indirect, incidental, special, exemplary, or consequential damages
\r
35 // (including, but not limited to, procurement of substitute goods or services;
\r
36 // loss of use, data, or profits; or business interruption) however caused
\r
37 // and on any theory of liability, whether in contract, strict liability,
\r
38 // or tort (including negligence or otherwise) arising in any way out of
\r
39 // the use of this software, even if advised of the possibility of such damage.
\r
43 #include "precomp.hpp"
\r
45 #if !defined(HAVE_CUDA)
\r
47 void cv::gpu::meanShiftSegmentation(const GpuMat&, Mat&, int, int, int, TermCriteria) { throw_nogpu(); }
\r
51 using namespace std;
\r
66 int merge(int set1, int set2);
\r
72 DjSets(const DjSets&);
\r
73 void operator =(const DjSets&);
\r
77 template <typename T>
\r
81 GraphEdge(int to, int next, const T& val) : to(to), next(next), val(val) {}
\r
88 template <typename T>
\r
92 typedef GraphEdge<T> Edge;
\r
94 Graph(int numv, int nume_max);
\r
96 void addEdge(int from, int to, const T& val=T());
\r
105 Graph(const Graph&);
\r
106 void operator =(const Graph&);
\r
113 SegmLinkVal(int dr, int dsp) : dr(dr), dsp(dsp) {}
\r
114 bool operator <(const SegmLinkVal& other) const
\r
116 return dr + dsp < other.dr + other.dsp;
\r
126 SegmLink(int from, int to, const SegmLinkVal& val)
\r
127 : from(from), to(to), val(val) {}
\r
128 bool operator <(const SegmLink& other) const
\r
130 return val < other.val;
\r
141 DjSets::DjSets(int n) : parent(n), rank(n, 0), size(n, 1)
\r
143 for (int i = 0; i < n; ++i)
\r
148 inline int DjSets::find(int elem)
\r
151 while (set != parent[set])
\r
153 while (elem != parent[elem])
\r
155 int next = parent[elem];
\r
156 parent[elem] = set;
\r
163 inline int DjSets::merge(int set1, int set2)
\r
165 if (rank[set1] < rank[set2])
\r
167 parent[set1] = set2;
\r
168 size[set2] += size[set1];
\r
171 if (rank[set2] < rank[set1])
\r
173 parent[set2] = set1;
\r
174 size[set1] += size[set2];
\r
177 parent[set1] = set2;
\r
179 size[set2] += size[set1];
\r
184 template <typename T>
\r
185 Graph<T>::Graph(int numv, int nume_max) : start(numv, -1), edges(nume_max)
\r
188 this->nume_max = nume_max;
\r
193 template <typename T>
\r
194 inline void Graph<T>::addEdge(int from, int to, const T& val)
\r
196 edges[nume] = Edge(to, start[from], val);
\r
197 start[from] = nume;
\r
202 inline int pix(int y, int x, int ncols)
\r
204 return y * ncols + x;
\r
208 inline int sqr(int x)
\r
214 inline int dist2(const cv::Vec4b& lhs, const cv::Vec4b& rhs)
\r
216 return sqr(lhs[0] - rhs[0]) + sqr(lhs[1] - rhs[1]) + sqr(lhs[2] - rhs[2]);
\r
220 inline int dist2(const cv::Vec2s& lhs, const cv::Vec2s& rhs)
\r
222 return sqr(lhs[0] - rhs[0]) + sqr(lhs[1] - rhs[1]);
\r
225 } // anonymous namespace
\r
228 void cv::gpu::meanShiftSegmentation(const GpuMat& src, Mat& dst, int sp, int sr, int minsize, TermCriteria criteria)
\r
230 CV_Assert(src.type() == CV_8UC4);
\r
231 const int nrows = src.rows;
\r
232 const int ncols = src.cols;
\r
234 const int hsp = sp;
\r
236 // Perform mean shift procedure and obtain region and spatial maps
\r
237 GpuMat d_rmap, d_spmap;
\r
238 meanShiftProc(src, d_rmap, d_spmap, sp, sr, criteria);
\r
240 Mat spmap(d_spmap);
\r
242 Graph<SegmLinkVal> g(nrows * ncols, 4 * (nrows - 1) * (ncols - 1)
\r
243 + (nrows - 1) + (ncols - 1));
\r
245 // Make region adjacent graph from image
\r
252 for (int y = 0; y < nrows - 1; ++y)
\r
254 Vec4b* ry = rmap.ptr<Vec4b>(y);
\r
255 Vec4b* ryp = rmap.ptr<Vec4b>(y + 1);
\r
256 Vec2s* spy = spmap.ptr<Vec2s>(y);
\r
257 Vec2s* spyp = spmap.ptr<Vec2s>(y + 1);
\r
258 for (int x = 0; x < ncols - 1; ++x)
\r
265 r2[2] = ryp[x + 1];
\r
268 sp2[0] = spy[x + 1];
\r
270 sp2[2] = spyp[x + 1];
\r
273 dr[0] = dist2(r1, r2[0]);
\r
274 dr[1] = dist2(r1, r2[1]);
\r
275 dr[2] = dist2(r1, r2[2]);
\r
276 dsp[0] = dist2(sp1, sp2[0]);
\r
277 dsp[1] = dist2(sp1, sp2[1]);
\r
278 dsp[2] = dist2(sp1, sp2[2]);
\r
283 dr[3] = dist2(r1, r2[3]);
\r
284 dsp[3] = dist2(sp1, sp2[3]);
\r
286 g.addEdge(pix(y, x, ncols), pix(y, x + 1, ncols), SegmLinkVal(dr[0], dsp[0]));
\r
287 g.addEdge(pix(y, x, ncols), pix(y + 1, x, ncols), SegmLinkVal(dr[1], dsp[1]));
\r
288 g.addEdge(pix(y, x, ncols), pix(y + 1, x + 1, ncols), SegmLinkVal(dr[2], dsp[2]));
\r
289 g.addEdge(pix(y, x + 1, ncols), pix(y + 1, x, ncols), SegmLinkVal(dr[3], dsp[3]));
\r
292 for (int y = 0; y < nrows - 1; ++y)
\r
294 r1 = rmap.at<Vec4b>(y, ncols - 1);
\r
295 r2[0] = rmap.at<Vec4b>(y + 1, ncols - 1);
\r
296 sp1 = spmap.at<Vec2s>(y, ncols - 1);
\r
297 sp2[0] = spmap.at<Vec2s>(y + 1, ncols - 1);
\r
298 dr[0] = dist2(r1, r2[0]);
\r
299 dsp[0] = dist2(sp1, sp2[0]);
\r
300 g.addEdge(pix(y, ncols - 1, ncols), pix(y + 1, ncols - 1, ncols), SegmLinkVal(dr[0], dsp[0]));
\r
302 for (int x = 0; x < ncols - 1; ++x)
\r
304 r1 = rmap.at<Vec4b>(nrows - 1, x);
\r
305 r2[0] = rmap.at<Vec4b>(nrows - 1, x + 1);
\r
306 sp1 = spmap.at<Vec2s>(nrows - 1, x);
\r
307 sp2[0] = spmap.at<Vec2s>(nrows - 1, x + 1);
\r
308 dr[0] = dist2(r1, r2[0]);
\r
309 dsp[0] = dist2(sp1, sp2[0]);
\r
310 g.addEdge(pix(nrows - 1, x, ncols), pix(nrows - 1, x + 1, ncols), SegmLinkVal(dr[0], dsp[0]));
\r
313 DjSets comps(g.numv);
\r
315 // Find adjacent components
\r
316 for (int v = 0; v < g.numv; ++v)
\r
318 for (int e_it = g.start[v]; e_it != -1; e_it = g.edges[e_it].next)
\r
320 int c1 = comps.find(v);
\r
321 int c2 = comps.find(g.edges[e_it].to);
\r
322 if (c1 != c2 && g.edges[e_it].val.dr < hr && g.edges[e_it].val.dsp < hsp)
\r
323 comps.merge(c1, c2);
\r
327 vector<SegmLink> edges;
\r
328 edges.reserve(g.numv);
\r
330 // Prepare edges connecting differnet components
\r
331 for (int v = 0; v < g.numv; ++v)
\r
333 int c1 = comps.find(v);
\r
334 for (int e_it = g.start[v]; e_it != -1; e_it = g.edges[e_it].next)
\r
336 int c2 = comps.find(g.edges[e_it].to);
\r
338 edges.push_back(SegmLink(c1, c2, g.edges[e_it].val));
\r
342 // Sort all graph's edges connecting differnet components (in asceding order)
\r
343 sort(edges.begin(), edges.end());
\r
345 // Exclude small components (starting from the nearest couple)
\r
346 for (size_t i = 0; i < edges.size(); ++i)
\r
348 int c1 = comps.find(edges[i].from);
\r
349 int c2 = comps.find(edges[i].to);
\r
350 if (c1 != c2 && (comps.size[c1] < minsize || comps.size[c2] < minsize))
\r
351 comps.merge(c1, c2);
\r
354 // Compute sum of the pixel's colors which are in the same segment
\r
356 vector<Vec4i> sumcols(nrows * ncols, Vec4i(0, 0, 0, 0));
\r
357 for (int y = 0; y < nrows; ++y)
\r
359 Vec4b* h_srcy = h_src.ptr<Vec4b>(y);
\r
360 for (int x = 0; x < ncols; ++x)
\r
362 int parent = comps.find(pix(y, x, ncols));
\r
363 Vec4b col = h_srcy[x];
\r
364 Vec4i& sumcol = sumcols[parent];
\r
365 sumcol[0] += col[0];
\r
366 sumcol[1] += col[1];
\r
367 sumcol[2] += col[2];
\r
371 // Create final image, color of each segment is the average color of its pixels
\r
372 dst.create(src.size(), src.type());
\r
374 for (int y = 0; y < nrows; ++y)
\r
376 Vec4b* dsty = dst.ptr<Vec4b>(y);
\r
377 for (int x = 0; x < ncols; ++x)
\r
379 int parent = comps.find(pix(y, x, ncols));
\r
380 const Vec4i& sumcol = sumcols[parent];
\r
381 Vec4b& dstcol = dsty[x];
\r
382 dstcol[0] = static_cast<uchar>(sumcol[0] / comps.size[parent]);
\r
383 dstcol[1] = static_cast<uchar>(sumcol[1] / comps.size[parent]);
\r
384 dstcol[2] = static_cast<uchar>(sumcol[2] / comps.size[parent]);
\r
389 #endif // #if !defined (HAVE_CUDA)
\r