1 /*M///////////////////////////////////////////////////////////////////////////////////////
3 // IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
5 // By downloading, copying, installing or using the software you agree to this license.
6 // If you do not agree to this license, do not download, install,
7 // copy or use the software.
11 // For Open Source Computer Vision Library
13 // Copyright (C) 2010-2012, Institute Of Software Chinese Academy Of Science, all rights reserved.
14 // Copyright (C) 2010-2012, Advanced Micro Devices, Inc., all rights reserved.
15 // Third party copyrights are property of their respective owners.
19 // Redistribution and use in source and binary forms, with or without modification,
20 // are permitted provided that the following conditions are met:
22 // * Redistribution's of source code must retain the above copyright notice,
23 // this list of conditions and the following disclaimer.
25 // * Redistribution's in binary form must reproduce the above copyright notice,
26 // this list of conditions and the following disclaimer in the documentation
27 // and/or other oclMaterials provided with the distribution.
29 // * The name of the copyright holders may not be used to endorse or promote products
30 // derived from this software without specific prior written permission.
32 // This software is provided by the copyright holders and contributors "as is" and
33 // any express or implied warranties, including, but not limited to, the implied
34 // warranties of merchantability and fitness for a particular purpose are disclaimed.
35 // In no event shall the Intel Corporation or contributors be liable for any direct,
36 // indirect, incidental, special, exemplary, or consequential damages
37 // (including, but not limited to, procurement of substitute goods or services;
38 // loss of use, data, or profits; or business interruption) however caused
39 // and on any theory of liability, whether in contract, strict liability,
40 // or tort (including negligence or otherwise) arising in any way out of
41 // the use of this software, even if advised of the possibility of such damage.
45 #include "precomp.hpp"
46 #include "opencl_kernels.hpp"
49 using namespace cv::ocl;
64 int merge(int set1, int set2);
66 std::vector<int> parent;
67 std::vector<int> rank;
68 std::vector<int> size;
70 DjSets(const DjSets &) {}
71 DjSets operator =(const DjSets &);
78 GraphEdge(int to, int next, const T &val) : to(to), next(next), val(val) {}
89 typedef GraphEdge<T> Edge;
91 Graph(int numv, int nume_max);
93 void addEdge(int from, int to, const T &val = T());
95 std::vector<int> start;
96 std::vector<Edge> edges;
102 Graph(const Graph &) {}
103 Graph operator =(const Graph &) {}
110 SegmLinkVal(int dr, int dsp) : dr(dr), dsp(dsp) {}
111 bool operator <(const SegmLinkVal &other) const
113 return dr + dsp < other.dr + other.dsp;
123 SegmLink(int from, int to, const SegmLinkVal &val)
124 : from(from), to(to), val(val) {}
125 bool operator <(const SegmLink &other) const
127 return val < other.val;
138 DjSets DjSets::operator = (const DjSets &/*obj*/)
140 //cout << "Invalid DjSets constructor\n";
141 CV_Error(-1, "Invalid DjSets constructor\n");
145 DjSets::DjSets(int n) : parent(n), rank(n, 0), size(n, 1)
147 for (int i = 0; i < n; ++i)
152 inline int DjSets::find(int elem)
155 while (set != parent[set])
157 while (elem != parent[elem])
159 int next = parent[elem];
167 inline int DjSets::merge(int set1, int set2)
169 if (rank[set1] < rank[set2])
172 size[set2] += size[set1];
175 if (rank[set2] < rank[set1])
178 size[set1] += size[set2];
183 size[set2] += size[set1];
188 template <typename T>
189 Graph<T>::Graph(int numv, int nume_max) : start(numv, -1), edges(nume_max)
192 this->nume_max = nume_max;
197 template <typename T>
198 inline void Graph<T>::addEdge(int from, int to, const T &val)
200 edges[nume] = Edge(to, start[from], val);
206 inline int pix(int y, int x, int ncols)
208 return y * ncols + x;
212 inline int sqr(int x)
218 inline int dist2(const cv::Vec4b &lhs, const cv::Vec4b &rhs)
220 return sqr(lhs[0] - rhs[0]) + sqr(lhs[1] - rhs[1]) + sqr(lhs[2] - rhs[2]);
224 inline int dist2(const cv::Vec2s &lhs, const cv::Vec2s &rhs)
226 return sqr(lhs[0] - rhs[0]) + sqr(lhs[1] - rhs[1]);
229 } // anonymous namespace
236 void meanShiftSegmentation(const oclMat &src, Mat &dst, int sp, int sr, int minsize, TermCriteria criteria)
238 CV_Assert(src.type() == CV_8UC4);
239 const int nrows = src.rows;
240 const int ncols = src.cols;
244 // Perform mean shift procedure and obtain region and spatial maps
245 oclMat h_rmap, h_spmap;
246 meanShiftProc(src, h_rmap, h_spmap, sp, sr, criteria);
250 Graph<SegmLinkVal> g(nrows * ncols, 4 * (nrows - 1) * (ncols - 1)
251 + (nrows - 1) + (ncols - 1));
253 // Make region adjacent graph from image
260 for (int y = 0; y < nrows - 1; ++y)
262 Vec4b *ry = rmap.ptr<Vec4b>(y);
263 Vec4b *ryp = rmap.ptr<Vec4b>(y + 1);
264 Vec2s *spy = spmap.ptr<Vec2s>(y);
265 Vec2s *spyp = spmap.ptr<Vec2s>(y + 1);
266 for (int x = 0; x < ncols - 1; ++x)
278 sp2[2] = spyp[x + 1];
281 dr[0] = dist2(r1, r2[0]);
282 dr[1] = dist2(r1, r2[1]);
283 dr[2] = dist2(r1, r2[2]);
284 dsp[0] = dist2(sp1, sp2[0]);
285 dsp[1] = dist2(sp1, sp2[1]);
286 dsp[2] = dist2(sp1, sp2[2]);
291 dr[3] = dist2(r1, r2[3]);
292 dsp[3] = dist2(sp1, sp2[3]);
294 g.addEdge(pix(y, x, ncols), pix(y, x + 1, ncols), SegmLinkVal(dr[0], dsp[0]));
295 g.addEdge(pix(y, x, ncols), pix(y + 1, x, ncols), SegmLinkVal(dr[1], dsp[1]));
296 g.addEdge(pix(y, x, ncols), pix(y + 1, x + 1, ncols), SegmLinkVal(dr[2], dsp[2]));
297 g.addEdge(pix(y, x + 1, ncols), pix(y + 1, x, ncols), SegmLinkVal(dr[3], dsp[3]));
300 for (int y = 0; y < nrows - 1; ++y)
302 r1 = rmap.at<Vec4b>(y, ncols - 1);
303 r2[0] = rmap.at<Vec4b>(y + 1, ncols - 1);
304 sp1 = spmap.at<Vec2s>(y, ncols - 1);
305 sp2[0] = spmap.at<Vec2s>(y + 1, ncols - 1);
306 dr[0] = dist2(r1, r2[0]);
307 dsp[0] = dist2(sp1, sp2[0]);
308 g.addEdge(pix(y, ncols - 1, ncols), pix(y + 1, ncols - 1, ncols), SegmLinkVal(dr[0], dsp[0]));
310 for (int x = 0; x < ncols - 1; ++x)
312 r1 = rmap.at<Vec4b>(nrows - 1, x);
313 r2[0] = rmap.at<Vec4b>(nrows - 1, x + 1);
314 sp1 = spmap.at<Vec2s>(nrows - 1, x);
315 sp2[0] = spmap.at<Vec2s>(nrows - 1, x + 1);
316 dr[0] = dist2(r1, r2[0]);
317 dsp[0] = dist2(sp1, sp2[0]);
318 g.addEdge(pix(nrows - 1, x, ncols), pix(nrows - 1, x + 1, ncols), SegmLinkVal(dr[0], dsp[0]));
321 DjSets comps(g.numv);
323 // Find adjacent components
324 for (int v = 0; v < g.numv; ++v)
326 for (int e_it = g.start[v]; e_it != -1; e_it = g.edges[e_it].next)
328 int c1 = comps.find(v);
329 int c2 = comps.find(g.edges[e_it].to);
330 if (c1 != c2 && g.edges[e_it].val.dr < hr && g.edges[e_it].val.dsp < hsp)
335 std::vector<SegmLink> edges;
336 edges.reserve(g.numv);
338 // Prepare edges connecting differnet components
339 for (int v = 0; v < g.numv; ++v)
341 int c1 = comps.find(v);
342 for (int e_it = g.start[v]; e_it != -1; e_it = g.edges[e_it].next)
344 int c2 = comps.find(g.edges[e_it].to);
346 edges.push_back(SegmLink(c1, c2, g.edges[e_it].val));
350 // Sort all graph's edges connecting differnet components (in asceding order)
351 sort(edges.begin(), edges.end());
353 // Exclude small components (starting from the nearest couple)
354 for (size_t i = 0; i < edges.size(); ++i)
356 int c1 = comps.find(edges[i].from);
357 int c2 = comps.find(edges[i].to);
358 if (c1 != c2 && (comps.size[c1] < minsize || comps.size[c2] < minsize))
362 // Compute sum of the pixel's colors which are in the same segment
364 std::vector<Vec4i> sumcols(nrows * ncols, Vec4i(0, 0, 0, 0));
365 for (int y = 0; y < nrows; ++y)
367 Vec4b *h_srcy = h_src.ptr<Vec4b>(y);
368 for (int x = 0; x < ncols; ++x)
370 int parent = comps.find(pix(y, x, ncols));
371 Vec4b col = h_srcy[x];
372 Vec4i &sumcol = sumcols[parent];
379 // Create final image, color of each segment is the average color of its pixels
380 dst.create(src.size(), src.type());
382 for (int y = 0; y < nrows; ++y)
384 Vec4b *dsty = dst.ptr<Vec4b>(y);
385 for (int x = 0; x < ncols; ++x)
387 int parent = comps.find(pix(y, x, ncols));
388 const Vec4i &sumcol = sumcols[parent];
389 Vec4b &dstcol = dsty[x];
390 dstcol[0] = static_cast<uchar>(sumcol[0] / comps.size[parent]);
391 dstcol[1] = static_cast<uchar>(sumcol[1] / comps.size[parent]);
392 dstcol[2] = static_cast<uchar>(sumcol[2] / comps.size[parent]);