modules/cudaimgproc/src/mssegmentation.cpp

   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) { 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)
 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);
 239     Mat rmap(d_rmap);
 240     Mat spmap(d_spmap);
 241
 242     Graph<SegmLinkVal> g(nrows * ncols, 4 * (nrows - 1) * (ncols - 1)
 243                                         + (nrows - 1) + (ncols - 1));
 244
 245     // Make region adjacent graph from image
 246     Vec4b r1;
 247     Vec4b r2[4];
 248     Vec2s sp1;
 249     Vec2s sp2[4];
 250     int dr[4];
 251     int dsp[4];
 252     for (int y = 0; y < nrows - 1; ++y)
 253     {
 254         Vec4b* ry = rmap.ptr<Vec4b>(y);
 255         Vec4b* ryp = rmap.ptr<Vec4b>(y + 1);
 256         Vec2s* spy = spmap.ptr<Vec2s>(y);
 257         Vec2s* spyp = spmap.ptr<Vec2s>(y + 1);
 258         for (int x = 0; x < ncols - 1; ++x)
 259         {
 260             r1 = ry[x];
 261             sp1 = spy[x];
 262
 263             r2[0] = ry[x + 1];
 264             r2[1] = ryp[x];
 265             r2[2] = ryp[x + 1];
 266             r2[3] = ryp[x];
 267
 268             sp2[0] = spy[x + 1];
 269             sp2[1] = spyp[x];
 270             sp2[2] = spyp[x + 1];
 271             sp2[3] = spyp[x];
 272
 273             dr[0] = dist2(r1, r2[0]);
 274             dr[1] = dist2(r1, r2[1]);
 275             dr[2] = dist2(r1, r2[2]);
 276             dsp[0] = dist2(sp1, sp2[0]);
 277             dsp[1] = dist2(sp1, sp2[1]);
 278             dsp[2] = dist2(sp1, sp2[2]);
 279
 280             r1 = ry[x + 1];
 281             sp1 = spy[x + 1];
 282
 283             dr[3] = dist2(r1, r2[3]);
 284             dsp[3] = dist2(sp1, sp2[3]);
 285
 286             g.addEdge(pix(y, x, ncols), pix(y, x + 1, ncols), SegmLinkVal(dr[0], dsp[0]));
 287             g.addEdge(pix(y, x, ncols), pix(y + 1, x, ncols), SegmLinkVal(dr[1], dsp[1]));
 288             g.addEdge(pix(y, x, ncols), pix(y + 1, x + 1, ncols), SegmLinkVal(dr[2], dsp[2]));
 289             g.addEdge(pix(y, x + 1, ncols), pix(y + 1, x, ncols), SegmLinkVal(dr[3], dsp[3]));
 290         }
 291     }
 292     for (int y = 0; y < nrows - 1; ++y)
 293     {
 294         r1 = rmap.at<Vec4b>(y, ncols - 1);
 295         r2[0] = rmap.at<Vec4b>(y + 1, ncols - 1);
 296         sp1 = spmap.at<Vec2s>(y, ncols - 1);
 297         sp2[0] = spmap.at<Vec2s>(y + 1, ncols - 1);
 298         dr[0] = dist2(r1, r2[0]);
 299         dsp[0] = dist2(sp1, sp2[0]);
 300         g.addEdge(pix(y, ncols - 1, ncols), pix(y + 1, ncols - 1, ncols), SegmLinkVal(dr[0], dsp[0]));
 301     }
 302     for (int x = 0; x < ncols - 1; ++x)
 303     {
 304         r1 = rmap.at<Vec4b>(nrows - 1, x);
 305         r2[0] = rmap.at<Vec4b>(nrows - 1, x + 1);
 306         sp1 = spmap.at<Vec2s>(nrows - 1, x);
 307         sp2[0] = spmap.at<Vec2s>(nrows - 1, x + 1);
 308         dr[0] = dist2(r1, r2[0]);
 309         dsp[0] = dist2(sp1, sp2[0]);
 310         g.addEdge(pix(nrows - 1, x, ncols), pix(nrows - 1, x + 1, ncols), SegmLinkVal(dr[0], dsp[0]));
 311     }
 312
 313     DjSets comps(g.numv);
 314
 315     // Find adjacent components
 316     for (int v = 0; v < g.numv; ++v)
 317     {
 318         for (int e_it = g.start[v]; e_it != -1; e_it = g.edges[e_it].next)
 319         {
 320             int c1 = comps.find(v);
 321             int c2 = comps.find(g.edges[e_it].to);
 322             if (c1 != c2 && g.edges[e_it].val.dr < hr && g.edges[e_it].val.dsp < hsp)
 323                 comps.merge(c1, c2);
 324         }
 325     }
 326
 327     std::vector<SegmLink> edges;
 328     edges.reserve(g.numv);
 329
 330     // Prepare edges connecting differnet components
 331     for (int v = 0; v < g.numv; ++v)
 332     {
 333         int c1 = comps.find(v);
 334         for (int e_it = g.start[v]; e_it != -1; e_it = g.edges[e_it].next)
 335         {
 336             int c2 = comps.find(g.edges[e_it].to);
 337             if (c1 != c2)
 338                 edges.push_back(SegmLink(c1, c2, g.edges[e_it].val));
 339         }
 340     }
 341
 342     // Sort all graph's edges connecting differnet components (in asceding order)
 343     std::sort(edges.begin(), edges.end());
 344
 345     // Exclude small components (starting from the nearest couple)
 346     for (size_t i = 0; i < edges.size(); ++i)
 347     {
 348         int c1 = comps.find(edges[i].from);
 349         int c2 = comps.find(edges[i].to);
 350         if (c1 != c2 && (comps.size[c1] < minsize || comps.size[c2] < minsize))
 351             comps.merge(c1, c2);
 352     }
 353
 354     // Compute sum of the pixel's colors which are in the same segment
 355     Mat h_src(src);
 356     std::vector<Vec4i> sumcols(nrows * ncols, Vec4i(0, 0, 0, 0));
 357     for (int y = 0; y < nrows; ++y)
 358     {
 359         Vec4b* h_srcy = h_src.ptr<Vec4b>(y);
 360         for (int x = 0; x < ncols; ++x)
 361         {
 362             int parent = comps.find(pix(y, x, ncols));
 363             Vec4b col = h_srcy[x];
 364             Vec4i& sumcol = sumcols[parent];
 365             sumcol[0] += col[0];
 366             sumcol[1] += col[1];
 367             sumcol[2] += col[2];
 368         }
 369     }
 370
 371     // Create final image, color of each segment is the average color of its pixels
 372     _dst.create(src.size(), src.type());
 373     Mat dst = _dst.getMat();
 374
 375     for (int y = 0; y < nrows; ++y)
 376     {
 377         Vec4b* dsty = dst.ptr<Vec4b>(y);
 378         for (int x = 0; x < ncols; ++x)
 379         {
 380             int parent = comps.find(pix(y, x, ncols));
 381             const Vec4i& sumcol = sumcols[parent];
 382             Vec4b& dstcol = dsty[x];
 383             dstcol[0] = static_cast<uchar>(sumcol[0] / comps.size[parent]);
 384             dstcol[1] = static_cast<uchar>(sumcol[1] / comps.size[parent]);
 385             dstcol[2] = static_cast<uchar>(sumcol[2] / comps.size[parent]);
 386             dstcol[3] = 255;
 387         }
 388     }
 389 }
 390
 391 #endif // #if !defined (HAVE_CUDA) || defined (CUDA_DISABLER)