1 /*M/////////////////////////////////////////////////////////////////////////////////////// 2 // 3 // IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING. 4 // 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. 8 // 9 // 10 // License Agreement 11 // For Open Source Computer Vision Library 12 // 13 // Copyright (C) 2000-2008, Intel Corporation, all rights reserved. 14 // Copyright (C) 2009, Willow Garage Inc., all rights reserved. 15 // Third party copyrights are property of their respective owners. 16 // 17 // Redistribution and use in source and binary forms, with or without modification, 18 // are permitted provided that the following conditions are met: 19 // 20 // * Redistribution's of source code must retain the above copyright notice, 21 // this list of conditions and the following disclaimer. 22 // 23 // * Redistribution's in binary form must reproduce the above copyright notice, 24 // this list of conditions and the following disclaimer in the documentation 25 // and/or other materials provided with the distribution. 26 // 27 // * The name of the copyright holders may not be used to endorse or promote products 28 // derived from this software without specific prior written permission. 29 // 30 // This software is provided by the copyright holders and contributors "as is" and 31 // any express or implied warranties, including, but not limited to, the implied 32 // warranties of merchantability and fitness for a particular purpose are disclaimed. 33 // In no event shall the Intel Corporation or contributors be liable for any direct, 34 // indirect, incidental, special, exemplary, or consequential damages 35 // (including, but not limited to, procurement of substitute goods or services; 36 // loss of use, data, or profits; or business interruption) however caused 37 // and on any theory of liability, whether in contract, strict liability, 38 // or tort (including negligence or otherwise) arising in any way out of 39 // the use of this software, even if advised of the possibility of such damage. 40 // 41 //M*/ 42 #include "precomp.hpp" 43 44 #if !defined HAVE_CUDA || defined(CUDA_DISABLER) 45 46 void cv::cuda::meanShiftSegmentation(InputArray, OutputArray, int, int, int, TermCriteria, Stream&) { throw_no_cuda(); } 47 48 #else 49 50 // Auxiliray stuff 51 namespace 52 { 53 54 // 55 // Declarations 56 // 57 58 class DjSets 59 { 60 public: 61 DjSets(int n); 62 int find(int elem); 63 int merge(int set1, int set2); 64 65 std::vector<int> parent; 66 std::vector<int> rank; 67 std::vector<int> size; 68 private: 69 DjSets(const DjSets&); 70 void operator =(const DjSets&); 71 }; 72 73 74 template <typename T> 75 struct GraphEdge 76 { 77 GraphEdge() {} 78 GraphEdge(int to_, int next_, const T& val_) : to(to_), next(next_), val(val_) {} 79 int to; 80 int next; 81 T val; 82 }; 83 84 85 template <typename T> 86 class Graph 87 { 88 public: 89 typedef GraphEdge<T> Edge; 90 91 Graph(int numv, int nume_max); 92 93 void addEdge(int from, int to, const T& val=T()); 94 95 std::vector<int> start; 96 std::vector<Edge> edges; 97 98 int numv; 99 int nume_max; 100 int nume; 101 private: 102 Graph(const Graph&); 103 void operator =(const Graph&); 104 }; 105 106 107 struct SegmLinkVal 108 { 109 SegmLinkVal() {} 110 SegmLinkVal(int dr_, int dsp_) : dr(dr_), dsp(dsp_) {} 111 bool operator <(const SegmLinkVal& other) const 112 { 113 return dr + dsp < other.dr + other.dsp; 114 } 115 int dr; 116 int dsp; 117 }; 118 119 120 struct SegmLink 121 { 122 SegmLink() {} 123 SegmLink(int from_, int to_, const SegmLinkVal& val_) 124 : from(from_), to(to_), val(val_) {} 125 bool operator <(const SegmLink& other) const 126 { 127 return val < other.val; 128 } 129 int from; 130 int to; 131 SegmLinkVal val; 132 }; 133 134 // 135 // Implementation 136 // 137 138 DjSets::DjSets(int n) : parent(n), rank(n, 0), size(n, 1) 139 { 140 for (int i = 0; i < n; ++i) 141 parent[i] = i; 142 } 143 144 145 inline int DjSets::find(int elem) 146 { 147 int set = elem; 148 while (set != parent[set]) 149 set = parent[set]; 150 while (elem != parent[elem]) 151 { 152 int next = parent[elem]; 153 parent[elem] = set; 154 elem = next; 155 } 156 return set; 157 } 158 159 160 inline int DjSets::merge(int set1, int set2) 161 { 162 if (rank[set1] < rank[set2]) 163 { 164 parent[set1] = set2; 165 size[set2] += size[set1]; 166 return set2; 167 } 168 if (rank[set2] < rank[set1]) 169 { 170 parent[set2] = set1; 171 size[set1] += size[set2]; 172 return set1; 173 } 174 parent[set1] = set2; 175 rank[set2]++; 176 size[set2] += size[set1]; 177 return set2; 178 } 179 180 181 template <typename T> 182 Graph<T>::Graph(int numv_, int nume_max_) : start(numv_, -1), edges(nume_max_) 183 { 184 this->numv = numv_; 185 this->nume_max = nume_max_; 186 nume = 0; 187 } 188 189 190 template <typename T> 191 inline void Graph<T>::addEdge(int from, int to, const T& val) 192 { 193 edges[nume] = Edge(to, start[from], val); 194 start[from] = nume; 195 nume++; 196 } 197 198 199 inline int pix(int y, int x, int ncols) 200 { 201 return y * ncols + x; 202 } 203 204 205 inline int sqr(int x) 206 { 207 return x * x; 208 } 209 210 211 inline int dist2(const cv::Vec4b& lhs, const cv::Vec4b& rhs) 212 { 213 return sqr(lhs[0] - rhs[0]) + sqr(lhs[1] - rhs[1]) + sqr(lhs[2] - rhs[2]); 214 } 215 216 217 inline int dist2(const cv::Vec2s& lhs, const cv::Vec2s& rhs) 218 { 219 return sqr(lhs[0] - rhs[0]) + sqr(lhs[1] - rhs[1]); 220 } 221 222 } // anonymous namespace 223 224 225 void cv::cuda::meanShiftSegmentation(InputArray _src, OutputArray _dst, int sp, int sr, int minsize, TermCriteria criteria, Stream& stream) 226 { 227 GpuMat src = _src.getGpuMat(); 228 229 CV_Assert( src.type() == CV_8UC4 ); 230 231 const int nrows = src.rows; 232 const int ncols = src.cols; 233 const int hr = sr; 234 const int hsp = sp; 235 236 // Perform mean shift procedure and obtain region and spatial maps 237 GpuMat d_rmap, d_spmap; 238 cuda::meanShiftProc(src, d_rmap, d_spmap, sp, sr, criteria, stream); 239 240 stream.waitForCompletion(); 241 242 Mat rmap(d_rmap); 243 Mat spmap(d_spmap); 244 245 Graph<SegmLinkVal> g(nrows * ncols, 4 * (nrows - 1) * (ncols - 1) 246 + (nrows - 1) + (ncols - 1)); 247 248 // Make region adjacent graph from image 249 Vec4b r1; 250 Vec4b r2[4]; 251 Vec2s sp1; 252 Vec2s sp2[4]; 253 int dr[4]; 254 int dsp[4]; 255 for (int y = 0; y < nrows - 1; ++y) 256 { 257 Vec4b* ry = rmap.ptr<Vec4b>(y); 258 Vec4b* ryp = rmap.ptr<Vec4b>(y + 1); 259 Vec2s* spy = spmap.ptr<Vec2s>(y); 260 Vec2s* spyp = spmap.ptr<Vec2s>(y + 1); 261 for (int x = 0; x < ncols - 1; ++x) 262 { 263 r1 = ry[x]; 264 sp1 = spy[x]; 265 266 r2[0] = ry[x + 1]; 267 r2[1] = ryp[x]; 268 r2[2] = ryp[x + 1]; 269 r2[3] = ryp[x]; 270 271 sp2[0] = spy[x + 1]; 272 sp2[1] = spyp[x]; 273 sp2[2] = spyp[x + 1]; 274 sp2[3] = spyp[x]; 275 276 dr[0] = dist2(r1, r2[0]); 277 dr[1] = dist2(r1, r2[1]); 278 dr[2] = dist2(r1, r2[2]); 279 dsp[0] = dist2(sp1, sp2[0]); 280 dsp[1] = dist2(sp1, sp2[1]); 281 dsp[2] = dist2(sp1, sp2[2]); 282 283 r1 = ry[x + 1]; 284 sp1 = spy[x + 1]; 285 286 dr[3] = dist2(r1, r2[3]); 287 dsp[3] = dist2(sp1, sp2[3]); 288 289 g.addEdge(pix(y, x, ncols), pix(y, x + 1, ncols), SegmLinkVal(dr[0], dsp[0])); 290 g.addEdge(pix(y, x, ncols), pix(y + 1, x, ncols), SegmLinkVal(dr[1], dsp[1])); 291 g.addEdge(pix(y, x, ncols), pix(y + 1, x + 1, ncols), SegmLinkVal(dr[2], dsp[2])); 292 g.addEdge(pix(y, x + 1, ncols), pix(y + 1, x, ncols), SegmLinkVal(dr[3], dsp[3])); 293 } 294 } 295 for (int y = 0; y < nrows - 1; ++y) 296 { 297 r1 = rmap.at<Vec4b>(y, ncols - 1); 298 r2[0] = rmap.at<Vec4b>(y + 1, ncols - 1); 299 sp1 = spmap.at<Vec2s>(y, ncols - 1); 300 sp2[0] = spmap.at<Vec2s>(y + 1, ncols - 1); 301 dr[0] = dist2(r1, r2[0]); 302 dsp[0] = dist2(sp1, sp2[0]); 303 g.addEdge(pix(y, ncols - 1, ncols), pix(y + 1, ncols - 1, ncols), SegmLinkVal(dr[0], dsp[0])); 304 } 305 for (int x = 0; x < ncols - 1; ++x) 306 { 307 r1 = rmap.at<Vec4b>(nrows - 1, x); 308 r2[0] = rmap.at<Vec4b>(nrows - 1, x + 1); 309 sp1 = spmap.at<Vec2s>(nrows - 1, x); 310 sp2[0] = spmap.at<Vec2s>(nrows - 1, x + 1); 311 dr[0] = dist2(r1, r2[0]); 312 dsp[0] = dist2(sp1, sp2[0]); 313 g.addEdge(pix(nrows - 1, x, ncols), pix(nrows - 1, x + 1, ncols), SegmLinkVal(dr[0], dsp[0])); 314 } 315 316 DjSets comps(g.numv); 317 318 // Find adjacent components 319 for (int v = 0; v < g.numv; ++v) 320 { 321 for (int e_it = g.start[v]; e_it != -1; e_it = g.edges[e_it].next) 322 { 323 int c1 = comps.find(v); 324 int c2 = comps.find(g.edges[e_it].to); 325 if (c1 != c2 && g.edges[e_it].val.dr < hr && g.edges[e_it].val.dsp < hsp) 326 comps.merge(c1, c2); 327 } 328 } 329 330 std::vector<SegmLink> edges; 331 edges.reserve(g.numv); 332 333 // Prepare edges connecting differnet components 334 for (int v = 0; v < g.numv; ++v) 335 { 336 int c1 = comps.find(v); 337 for (int e_it = g.start[v]; e_it != -1; e_it = g.edges[e_it].next) 338 { 339 int c2 = comps.find(g.edges[e_it].to); 340 if (c1 != c2) 341 edges.push_back(SegmLink(c1, c2, g.edges[e_it].val)); 342 } 343 } 344 345 // Sort all graph's edges connecting differnet components (in asceding order) 346 std::sort(edges.begin(), edges.end()); 347 348 // Exclude small components (starting from the nearest couple) 349 for (size_t i = 0; i < edges.size(); ++i) 350 { 351 int c1 = comps.find(edges[i].from); 352 int c2 = comps.find(edges[i].to); 353 if (c1 != c2 && (comps.size[c1] < minsize || comps.size[c2] < minsize)) 354 comps.merge(c1, c2); 355 } 356 357 // Compute sum of the pixel's colors which are in the same segment 358 Mat h_src(src); 359 std::vector<Vec4i> sumcols(nrows * ncols, Vec4i(0, 0, 0, 0)); 360 for (int y = 0; y < nrows; ++y) 361 { 362 Vec4b* h_srcy = h_src.ptr<Vec4b>(y); 363 for (int x = 0; x < ncols; ++x) 364 { 365 int parent = comps.find(pix(y, x, ncols)); 366 Vec4b col = h_srcy[x]; 367 Vec4i& sumcol = sumcols[parent]; 368 sumcol[0] += col[0]; 369 sumcol[1] += col[1]; 370 sumcol[2] += col[2]; 371 } 372 } 373 374 // Create final image, color of each segment is the average color of its pixels 375 _dst.create(src.size(), src.type()); 376 Mat dst = _dst.getMat(); 377 378 for (int y = 0; y < nrows; ++y) 379 { 380 Vec4b* dsty = dst.ptr<Vec4b>(y); 381 for (int x = 0; x < ncols; ++x) 382 { 383 int parent = comps.find(pix(y, x, ncols)); 384 const Vec4i& sumcol = sumcols[parent]; 385 Vec4b& dstcol = dsty[x]; 386 dstcol[0] = static_cast<uchar>(sumcol[0] / comps.size[parent]); 387 dstcol[1] = static_cast<uchar>(sumcol[1] / comps.size[parent]); 388 dstcol[2] = static_cast<uchar>(sumcol[2] / comps.size[parent]); 389 dstcol[3] = 255; 390 } 391 } 392 } 393 394 #endif // #if !defined (HAVE_CUDA) || defined (CUDA_DISABLER) 395