modules/gpu/src/mssegmentation.cpp

   1 /*M///////////////////////////////////////////////////////////////////////////////////////\r
   2 //\r
   3 //  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.\r
   4 //\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
   8 //\r
   9 //\r
  10 //                           License Agreement\r
  11 //                For Open Source Computer Vision Library\r
  12 //\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
  16 //\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
  19 //\r
  20 //   * Redistribution's of source code must retain the above copyright notice,\r
  21 //     this list of conditions and the following disclaimer.\r
  22 //\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
  26 //\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
  29 //\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
  40 //\r
  41 //M*/\r
  42 \r
  43 #include "precomp.hpp"\r
  44 \r
  45 #if !defined(HAVE_CUDA)\r
  46 \r
  47 void cv::gpu::meanShiftSegmentation(const GpuMat&, Mat&, int, int, int, TermCriteria) { throw_nogpu(); }\r
  48 \r
  49 #else\r
  50 \r
  51 using namespace std;\r
  52 \r
  53 // Auxiliray stuff\r
  54 namespace\r
  55 {\r
  56 \r
  57 //\r
  58 // Declarations\r
  59 //\r
  60 \r
  61 class DjSets\r
  62 {\r
  63 public:\r
  64     DjSets(int n);\r
  65     int find(int elem);\r
  66     int merge(int set1, int set2);\r
  67 \r
  68     vector<int> parent;\r
  69     vector<int> rank;\r
  70     vector<int> size;\r
  71 private:\r
  72     DjSets(const DjSets&);\r
  73     void operator =(const DjSets&);\r
  74 };\r
  75 \r
  76 \r
  77 template <typename T>\r
  78 struct GraphEdge\r
  79 {\r
  80     GraphEdge() {}\r
  81     GraphEdge(int to, int next, const T& val) : to(to), next(next), val(val) {}\r
  82     int to;\r
  83     int next;\r
  84     T val;\r
  85 };\r
  86 \r
  87 \r
  88 template <typename T>\r
  89 class Graph\r
  90 {\r
  91 public:\r
  92     typedef GraphEdge<T> Edge;\r
  93 \r
  94     Graph(int numv, int nume_max);\r
  95 \r
  96     void addEdge(int from, int to, const T& val=T());\r
  97 \r
  98     vector<int> start;\r
  99     vector<Edge> edges;\r
 100 \r
 101     int numv;\r
 102     int nume_max;\r
 103     int nume;\r
 104 private:\r
 105     Graph(const Graph&);\r
 106     void operator =(const Graph&);\r
 107 };\r
 108 \r
 109 \r
 110 struct SegmLinkVal\r
 111 {\r
 112     SegmLinkVal() {}\r
 113     SegmLinkVal(int dr, int dsp) : dr(dr), dsp(dsp) {}\r
 114     bool operator <(const SegmLinkVal& other) const\r
 115     {\r
 116         return dr + dsp < other.dr + other.dsp;\r
 117     }\r
 118     int dr;\r
 119     int dsp;\r
 120 };\r
 121 \r
 122 \r
 123 struct SegmLink\r
 124 {\r
 125     SegmLink() {}\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
 129     {\r
 130         return val < other.val;\r
 131     }\r
 132     int from;\r
 133     int to;\r
 134     SegmLinkVal val;\r
 135 };\r
 136 \r
 137 //\r
 138 // Implementation\r
 139 //\r
 140 \r
 141 DjSets::DjSets(int n) : parent(n), rank(n, 0), size(n, 1)\r
 142 {\r
 143     for (int i = 0; i < n; ++i)\r
 144         parent[i] = i;\r
 145 }\r
 146 \r
 147 \r
 148 inline int DjSets::find(int elem)\r
 149 {\r
 150     int set = elem;\r
 151     while (set != parent[set])\r
 152         set = parent[set];\r
 153     while (elem != parent[elem])\r
 154     {\r
 155         int next = parent[elem];\r
 156         parent[elem] = set;\r
 157         elem = next;\r
 158     }\r
 159     return set;\r
 160 }\r
 161 \r
 162 \r
 163 inline int DjSets::merge(int set1, int set2)\r
 164 {\r
 165     if (rank[set1] < rank[set2])\r
 166     {\r
 167         parent[set1] = set2;\r
 168         size[set2] += size[set1];\r
 169         return set2;\r
 170     }\r
 171     if (rank[set2] < rank[set1])\r
 172     {\r
 173         parent[set2] = set1;\r
 174         size[set1] += size[set2];\r
 175         return set1;\r
 176     }\r
 177     parent[set1] = set2;\r
 178     rank[set2]++;\r
 179     size[set2] += size[set1];\r
 180     return set2;\r
 181 }\r
 182 \r
 183 \r
 184 template <typename T>\r
 185 Graph<T>::Graph(int numv, int nume_max) : start(numv, -1), edges(nume_max)\r
 186 {\r
 187     this->numv = numv;\r
 188     this->nume_max = nume_max;\r
 189     nume = 0;\r
 190 }\r
 191 \r
 192 \r
 193 template <typename T>\r
 194 inline void Graph<T>::addEdge(int from, int to, const T& val)\r
 195 {\r
 196     edges[nume] = Edge(to, start[from], val);\r
 197     start[from] = nume;\r
 198     nume++;\r
 199 }\r
 200 \r
 201 \r
 202 inline int pix(int y, int x, int ncols)\r
 203 {\r
 204     return y * ncols + x;\r
 205 }\r
 206 \r
 207 \r
 208 inline int sqr(int x)\r
 209 {\r
 210     return x * x;\r
 211 }\r
 212 \r
 213 \r
 214 inline int dist2(const cv::Vec4b& lhs, const cv::Vec4b& rhs)\r
 215 {\r
 216     return sqr(lhs[0] - rhs[0]) + sqr(lhs[1] - rhs[1]) + sqr(lhs[2] - rhs[2]);\r
 217 }\r
 218 \r
 219 \r
 220 inline int dist2(const cv::Vec2s& lhs, const cv::Vec2s& rhs)\r
 221 {\r
 222     return sqr(lhs[0] - rhs[0]) + sqr(lhs[1] - rhs[1]);\r
 223 }\r
 224 \r
 225 } // anonymous namespace\r
 226 \r
 227 \r
 228 void cv::gpu::meanShiftSegmentation(const GpuMat& src, Mat& dst, int sp, int sr, int minsize, TermCriteria criteria)\r
 229 {\r
 230     CV_Assert(src.type() == CV_8UC4);\r
 231     const int nrows = src.rows;\r
 232     const int ncols = src.cols;\r
 233     const int hr = sr;\r
 234     const int hsp = sp;\r
 235 \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
 239     Mat rmap(d_rmap);\r
 240     Mat spmap(d_spmap);\r
 241 \r
 242     Graph<SegmLinkVal> g(nrows * ncols, 4 * (nrows - 1) * (ncols - 1)\r
 243                                         + (nrows - 1) + (ncols - 1));\r
 244 \r
 245     // Make region adjacent graph from image\r
 246     Vec4b r1;\r
 247     Vec4b r2[4];\r
 248     Vec2s sp1;\r
 249     Vec2s sp2[4];\r
 250     int dr[4];\r
 251     int dsp[4];\r
 252     for (int y = 0; y < nrows - 1; ++y)\r
 253     {\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
 259         {\r
 260             r1 = ry[x];\r
 261             sp1 = spy[x];\r
 262 \r
 263             r2[0] = ry[x + 1];\r
 264             r2[1] = ryp[x];\r
 265             r2[2] = ryp[x + 1];\r
 266             r2[3] = ryp[x];\r
 267 \r
 268             sp2[0] = spy[x + 1];\r
 269             sp2[1] = spyp[x];\r
 270             sp2[2] = spyp[x + 1];\r
 271             sp2[3] = spyp[x];\r
 272 \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
 279 \r
 280             r1 = ry[x + 1];\r
 281             sp1 = spy[x + 1];\r
 282 \r
 283             dr[3] = dist2(r1, r2[3]);\r
 284             dsp[3] = dist2(sp1, sp2[3]);\r
 285 \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
 290         }\r
 291     }\r
 292     for (int y = 0; y < nrows - 1; ++y)\r
 293     {\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
 301     }\r
 302     for (int x = 0; x < ncols - 1; ++x)\r
 303     {\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
 311     }\r
 312 \r
 313     DjSets comps(g.numv);\r
 314 \r
 315     // Find adjacent components\r
 316     for (int v = 0; v < g.numv; ++v)\r
 317     {\r
 318         for (int e_it = g.start[v]; e_it != -1; e_it = g.edges[e_it].next)\r
 319         {\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
 324         }\r
 325     }\r
 326 \r
 327     vector<SegmLink> edges;\r
 328     edges.reserve(g.numv);\r
 329 \r
 330     // Prepare edges connecting differnet components\r
 331     for (int v = 0; v < g.numv; ++v)\r
 332     {\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
 335         {\r
 336             int c2 = comps.find(g.edges[e_it].to);\r
 337             if (c1 != c2)\r
 338                 edges.push_back(SegmLink(c1, c2, g.edges[e_it].val));\r
 339         }\r
 340     }\r
 341 \r
 342     // Sort all graph's edges connecting differnet components (in asceding order)\r
 343     sort(edges.begin(), edges.end());\r
 344 \r
 345     // Exclude small components (starting from the nearest couple)\r
 346     for (size_t i = 0; i < edges.size(); ++i)\r
 347     {\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
 352     }\r
 353 \r
 354     // Compute sum of the pixel's colors which are in the same segment\r
 355     Mat h_src(src);\r
 356     vector<Vec4i> sumcols(nrows * ncols, Vec4i(0, 0, 0, 0));\r
 357     for (int y = 0; y < nrows; ++y)\r
 358     {\r
 359         Vec4b* h_srcy = h_src.ptr<Vec4b>(y);\r
 360         for (int x = 0; x < ncols; ++x)\r
 361         {\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
 368         }\r
 369     }\r
 370 \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
 373 \r
 374     for (int y = 0; y < nrows; ++y)\r
 375     {\r
 376         Vec4b* dsty = dst.ptr<Vec4b>(y);\r
 377         for (int x = 0; x < ncols; ++x)\r
 378         {\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
 385         }\r
 386     }\r
 387 }\r
 388 \r
 389 #endif // #if !defined (HAVE_CUDA)\r