implemented mean shift segmentation with elimination of small segments, added tests

diff --git a/modules/gpu/include/opencv2/gpu/gpu.hpp b/modules/gpu/include/opencv2/gpu/gpu.hpp
index b8ee807..1af4088 100644 (file)
--- a/modules/gpu/include/opencv2/gpu/gpu.hpp
+++ b/modules/gpu/include/opencv2/gpu/gpu.hpp
@@ -473,6 +473,10 @@ namespace cv
         CV_EXPORTS void meanShiftProc(const GpuMat& src, GpuMat& dstr, GpuMat& dstsp, int sp, int sr, \r
             TermCriteria criteria = TermCriteria(TermCriteria::MAX_ITER + TermCriteria::EPS, 5, 1));\r
 \r
+        //! Does mean shift segmentation with elimiation of small regions.\r
+        CV_EXPORTS void meanShiftSegmentation(const GpuMat& src, Mat& dst, int sp, int sr, int minsize,\r
+            TermCriteria criteria = TermCriteria(TermCriteria::MAX_ITER + TermCriteria::EPS, 5, 1));\r
+\r
         //! Does coloring of disparity image: [0..ndisp) -> [0..240, 1, 1] in HSV.\r
         //! Supported types of input disparity: CV_8U, CV_16S.\r
         //! Output disparity has CV_8UC4 type in BGRA format (alpha = 255).\r

diff --git a/modules/gpu/src/mssegmentation.cpp b/modules/gpu/src/mssegmentation.cpp
new file mode 100644 (file)
index 0000000..0600ea9
--- /dev/null
+++ b/modules/gpu/src/mssegmentation.cpp
@@ -0,0 +1,475 @@
+/*M///////////////////////////////////////////////////////////////////////////////////////\r
+//\r
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.\r
+//\r
+//  By downloading, copying, installing or using the software you agree to this license.\r
+//  If you do not agree to this license, do not download, install,\r
+//  copy or use the software.\r
+//\r
+//\r
+//                           License Agreement\r
+//                For Open Source Computer Vision Library\r
+//\r
+// Copyright (C) 2000-2008, Intel Corporation, all rights reserved.\r
+// Copyright (C) 2009, Willow Garage Inc., all rights reserved.\r
+// Third party copyrights are property of their respective owners.\r
+//\r
+// Redistribution and use in source and binary forms, with or without modification,\r
+// are permitted provided that the following conditions are met:\r
+//\r
+//   * Redistribution's of source code must retain the above copyright notice,\r
+//     this list of conditions and the following disclaimer.\r
+//\r
+//   * Redistribution's in binary form must reproduce the above copyright notice,\r
+//     this list of conditions and the following disclaimer in the documentation\r
+//     and/or other materials provided with the distribution.\r
+//\r
+//   * The name of the copyright holders may not be used to endorse or promote products\r
+//     derived from this software without specific prior written permission.\r
+//\r
+// This software is provided by the copyright holders and contributors "as is" and\r
+// any express or implied warranties, including, but not limited to, the implied\r
+// warranties of merchantability and fitness for a particular purpose are disclaimed.\r
+// In no event shall the Intel Corporation or contributors be liable for any direct,\r
+// indirect, incidental, special, exemplary, or consequential damages\r
+// (including, but not limited to, procurement of substitute goods or services;\r
+// loss of use, data, or profits; or business interruption) however caused\r
+// and on any theory of liability, whether in contract, strict liability,\r
+// or tort (including negligence or otherwise) arising in any way out of\r
+// the use of this software, even if advised of the possibility of such damage.\r
+//\r
+//M*/\r
+\r
+#include <time.h>\r
+#include <vector>\r
+#include "precomp.hpp"\r
+\r
+#if !defined(HAVE_CUDA)\r
+\r
+namespace cv\r
+{\r
+namespace gpu\r
+{\r
+\r
+void meanShiftSegmentation(const GpuMat&, Mat&, int, int, int, TermCriteria) { throw_nogpu(); }\r
+\r
+} // namespace gpu\r
+} // namespace cv\r
+\r
+#else\r
+\r
+//#define _MSSEGMENTATION_DBG\r
+\r
+#ifdef _MSSEGMENTATION_DBG\r
+#include <iostream>\r
+#define LOG(s) std::cout << (s) << std::endl\r
+#define LOG2(s1, s2) std::cout << (s1) << (s2) << std::endl\r
+#define DBG(code) code\r
+#else\r
+#define LOG(s1)\r
+#define LOG2(s1, s2)\r
+#define DBG(code)\r
+#endif\r
+\r
+#define PIX(y, x) ((y) * ncols + (x))\r
+\r
+using namespace std;\r
+\r
+// Auxiliray stuff\r
+namespace\r
+{\r
+\r
+//\r
+// Declarations\r
+//\r
+\r
+class DjSets\r
+{\r
+public:\r
+    DjSets(int n);\r
+    ~DjSets();\r
+    int find(int elem) const;\r
+    int merge(int set1, int set2);\r
+\r
+    int* parent;\r
+    int* rank;\r
+    int* size;\r
+private:\r
+    DjSets(const DjSets&) {}\r
+    DjSets operator =(const DjSets&) {}\r
+};\r
+\r
+\r
+template <typename T>\r
+struct GraphEdge\r
+{\r
+    GraphEdge() {}\r
+    GraphEdge(int to, int next, const T& val) : to(to), next(next), val(val) {}\r
+    int to;\r
+    int next;\r
+    T val;\r
+};\r
+\r
+\r
+template <typename T>\r
+class Graph\r
+{\r
+public:\r
+    typedef GraphEdge<T> Edge;\r
+\r
+    Graph(int numv, int nume_max);\r
+    ~Graph();\r
+\r
+    void addEdge(int from, int to, const T& val=T());\r
+\r
+    int* start;\r
+    Edge* edges;\r
+\r
+    int numv;\r
+    int nume_max;\r
+    int nume;\r
+private:\r
+    Graph(const Graph&) {}\r
+    Graph operator =(const Graph&) {}\r
+};\r
+\r
+\r
+struct SegmLinkVal\r
+{\r
+    SegmLinkVal() {}\r
+    SegmLinkVal(int dr, int dsp) : dr(dr), dsp(dsp) {}\r
+    bool operator <(const SegmLinkVal& other) const\r
+    {\r
+        return dr + dsp < other.dr + other.dsp;\r
+    }\r
+    int dr;\r
+    int dsp;\r
+};\r
+\r
+\r
+struct SegmLink\r
+{\r
+    SegmLink() {}\r
+    SegmLink(int from, int to, const SegmLinkVal& val) \r
+        : from(from), to(to), val(val) {}\r
+    int from;\r
+    int to;\r
+    SegmLinkVal val;\r
+};\r
+\r
+\r
+struct SegmLinkCmp\r
+{\r
+    bool operator ()(const SegmLink& lhs, const SegmLink& rhs) const\r
+    {\r
+        return lhs.val < rhs.val;\r
+    }\r
+};\r
+\r
+//\r
+// Implementation\r
+//\r
+\r
+DjSets::DjSets(int n)\r
+{\r
+    parent = new int[n];\r
+    rank = new int[n];\r
+    size = new int[n];   \r
+    for (int i = 0; i < n; ++i)\r
+    {\r
+        parent[i] = i;\r
+        rank[i] = 0;\r
+        size[i] = 1;\r
+    }\r
+}\r
+\r
+\r
+DjSets::~DjSets()\r
+{\r
+    delete[] parent;\r
+    delete[] rank;\r
+    delete[] size;\r
+}\r
+\r
+\r
+inline int DjSets::find(int elem) const\r
+{\r
+    int set = elem;\r
+    while (set != parent[set])\r
+        set = parent[set];\r
+    while (elem != parent[elem])\r
+    {\r
+        int next = parent[elem];\r
+        parent[elem] = set;\r
+        elem = next;\r
+    }\r
+    return set;\r
+}\r
+\r
+\r
+inline int DjSets::merge(int set1, int set2)\r
+{\r
+    if (rank[set1] < rank[set2])\r
+    {\r
+        parent[set1] = set2;\r
+        size[set2] += size[set1];\r
+        return set2;\r
+    }\r
+    if (rank[set2] < rank[set1])\r
+    {\r
+        parent[set2] = set1;\r
+        size[set1] += size[set2];\r
+        return set1;\r
+    }\r
+    parent[set1] = set2;\r
+    rank[set2]++;\r
+    size[set2] += size[set1];\r
+    return set2;\r
+}\r
+\r
+\r
+template <typename T>\r
+Graph<T>::Graph(int numv, int nume_max)\r
+{\r
+    this->numv = numv;\r
+    this->nume_max = nume_max;\r
+    start = new int[numv];\r
+    for (int i = 0; i < numv; ++i)\r
+        start[i] = -1;\r
+    edges = new Edge[nume_max];\r
+    nume = 0;\r
+}\r
+\r
+\r
+template <typename T>\r
+Graph<T>::~Graph()\r
+{\r
+    delete[] start;\r
+    delete[] edges;\r
+}\r
+\r
+\r
+template <typename T>\r
+inline void Graph<T>::addEdge(int from, int to, const T& val)\r
+{\r
+    Edge* edge = edges + nume;\r
+    new (edge) SegmLink(to, start[from], val);\r
+    start[from] = nume;\r
+    nume++;\r
+}\r
+\r
+\r
+inline int sqr(int x)\r
+{\r
+    return x * x;\r
+}\r
+\r
+} // anonymous namespace\r
+\r
+\r
+namespace cv\r
+{\r
+namespace gpu\r
+{\r
+\r
+void meanShiftSegmentation(const GpuMat& src, Mat& dst, int sp, int sr, int minsize, TermCriteria criteria)\r
+{\r
+    CV_Assert(src.type() == CV_8UC4);\r
+    const int nrows = src.rows;\r
+    const int ncols = src.cols;\r
+    const int hr = sr;\r
+    const int hsp = sp;\r
+\r
+    DBG(clock_t start = clock());\r
+\r
+    // Perform mean shift procedure and obtain region and spatial maps\r
+    GpuMat h_rmap, h_spmap;\r
+    meanShiftProc(src, h_rmap, h_spmap, sp, sr, criteria);\r
+    Mat rmap = h_rmap;\r
+    Mat spmap = h_spmap;\r
+\r
+    LOG2("meanshift:", clock() - start);\r
+    DBG(start = clock());\r
+\r
+    Graph<SegmLinkVal> g(nrows * ncols, 4 * (nrows - 1) * (ncols - 1)\r
+                                        + (nrows - 1) + (ncols - 1));\r
+\r
+    LOG2("ragalloc:", clock() - start);\r
+    DBG(start = clock());\r
+\r
+    // Make region adjacent graph from image\r
+    // TODO: SSE?\r
+    Vec4b r1;\r
+    Vec4b r2[4];\r
+    Point_<short> sp1;\r
+    Point_<short> sp2[4];\r
+    int dr[4];\r
+    int dsp[4];\r
+    for (int y = 0; y < nrows - 1; ++y)\r
+    {\r
+        Vec4b* ry = rmap.ptr<Vec4b>(y);\r
+        Vec4b* ryp = rmap.ptr<Vec4b>(y + 1);\r
+        Point_<short>* spy = spmap.ptr<Point_<short> >(y);\r
+        Point_<short>* spyp = spmap.ptr<Point_<short> >(y + 1);\r
+        for (int x = 0; x < ncols - 1; ++x)\r
+        {\r
+            r1 = ry[x];\r
+            sp1 = spy[x];\r
+\r
+            r2[0] = ry[x + 1];\r
+            r2[1] = ryp[x];\r
+            r2[2] = ryp[x + 1];\r
+            r2[3] = ryp[x];\r
+\r
+            sp2[0] = spy[x + 1];\r
+            sp2[1] = spyp[x];\r
+            sp2[2] = spyp[x + 1];\r
+            sp2[3] = spyp[x];\r
+\r
+            dr[0] = sqr(r1[0] - r2[0][0]) + sqr(r1[1] - r2[0][1]) + sqr(r1[2] - r2[0][2]);\r
+            dr[1] = sqr(r1[0] - r2[1][0]) + sqr(r1[1] - r2[1][1]) + sqr(r1[2] - r2[1][2]);\r
+            dr[2] = sqr(r1[0] - r2[2][0]) + sqr(r1[1] - r2[2][1]) + sqr(r1[2] - r2[2][2]);\r
+            dsp[0] = sqr(sp1.x - sp2[0].x) + sqr(sp1.y - sp2[0].y);\r
+            dsp[1] = sqr(sp1.x - sp2[1].x) + sqr(sp1.y - sp2[1].y);\r
+            dsp[2] = sqr(sp1.x - sp2[2].x) + sqr(sp1.y - sp2[2].y); \r
+\r
+            r1 = ry[x + 1];\r
+            sp1 = spy[x + 1];\r
+\r
+            dr[3] = sqr(r1[0] - r2[3][0]) + sqr(r1[1] - r2[3][1]) + sqr(r1[2] - r2[3][2]);\r
+            dsp[3] = sqr(sp1.x - sp2[3].x) + sqr(sp1.y - sp2[3].y); \r
+\r
+            g.addEdge(PIX(y, x), PIX(y, x + 1), SegmLinkVal(dr[0], dsp[0]));\r
+            g.addEdge(PIX(y, x), PIX(y + 1, x), SegmLinkVal(dr[1], dsp[1]));\r
+            g.addEdge(PIX(y, x), PIX(y + 1, x + 1), SegmLinkVal(dr[2], dsp[2]));\r
+            g.addEdge(PIX(y, x + 1), PIX(y, x + 1), SegmLinkVal(dr[3], dsp[3]));\r
+        }\r
+    }\r
+    for (int y = 0; y < nrows - 1; ++y)\r
+    {\r
+        r1 = rmap.at<Vec4b>(y, ncols - 1);\r
+        r2[0] = rmap.at<Vec4b>(y + 1, ncols - 1);\r
+        sp1 = spmap.at<Point_<short> >(y, ncols - 1);\r
+        sp2[0] = spmap.at<Point_<short> >(y + 1, ncols - 1);\r
+        dr[0] = sqr(r1[0] - r2[0][0]) + sqr(r1[1] - r2[0][1]) + sqr(r1[2] - r2[0][2]);\r
+        dsp[0] = sqr(sp1.x - sp2[0].x) + sqr(sp1.y - sp2[0].y);   \r
+        g.addEdge(PIX(y, ncols - 1), PIX(y + 1, ncols - 1), SegmLinkVal(dr[0], dsp[0]));\r
+    }\r
+    for (int x = 0; x < ncols - 1; ++x)\r
+    {\r
+        r1 = rmap.at<Vec4b>(nrows - 1, x);\r
+        r2[0] = rmap.at<Vec4b>(nrows - 1, x + 1);\r
+        sp1 = spmap.at<Point_<short> >(nrows - 1, x);\r
+        sp2[0] = spmap.at<Point_<short> >(nrows - 1, x + 1);\r
+        dr[0] = sqr(r1[0] - r2[0][0]) + sqr(r1[1] - r2[0][1]) + sqr(r1[2] - r2[0][2]);\r
+        dsp[0] = sqr(sp1.x - sp2[0].x) + sqr(sp1.y - sp2[0].y);   \r
+        g.addEdge(PIX(nrows - 1, x), PIX(nrows - 1, x + 1), SegmLinkVal(dr[0], dsp[0]));\r
+    }\r
+\r
+    LOG2("raginit:", clock() - start);\r
+    DBG(start = clock());\r
+\r
+    DjSets comps(g.numv);\r
+\r
+    LOG2("djsetinit:", clock() - start);\r
+    DBG(start = clock());\r
+\r
+    // Find adjacent components\r
+    for (int v = 0; v < g.numv; ++v)\r
+    {\r
+        for (int e_it = g.start[v]; e_it != -1; e_it = g.edges[e_it].next)\r
+        {\r
+            int comp1 = comps.find(v);\r
+            int comp2 = comps.find(g.edges[e_it].to);\r
+            if (comp1 != comp2 && g.edges[e_it].val.dr < hr \r
+                               && g.edges[e_it].val.dsp < hsp)\r
+                comps.merge(comp1, comp2);\r
+        }\r
+    }\r
+\r
+    LOG2("findadjacent:", clock() - start);\r
+    DBG(start = clock());\r
+\r
+    vector<SegmLink> edges;\r
+    edges.reserve(g.numv);\r
+\r
+    LOG2("initedges:", clock() - start);\r
+    DBG(start = clock());\r
+\r
+    for (int v = 0; v < g.numv; ++v)\r
+    {\r
+        int comp1 = comps.find(v);\r
+        for (int e_it = g.start[v]; e_it != -1; e_it = g.edges[e_it].next)\r
+        {\r
+            int comp2 = comps.find(g.edges[e_it].to);\r
+            if (comp1 != comp2)\r
+                edges.push_back(SegmLink(comp1, comp2, g.edges[e_it].val));\r
+        }\r
+    }\r
+\r
+    LOG2("prepareforsort:", clock() - start);\r
+    DBG(start = clock());\r
+\r
+    // Sort all graph's edges connecting differnet components (in asceding order)\r
+    sort(edges.begin(), edges.end(), SegmLinkCmp());\r
+\r
+    LOG2("sortedges:", clock() - start);\r
+    DBG(start = clock());\r
+\r
+    // Exclude small components (starting from the nearest couple)\r
+    vector<SegmLink>::iterator e_it = edges.begin();\r
+    for (; e_it != edges.end(); ++e_it)\r
+    {\r
+        int comp1 = comps.find(e_it->from);\r
+        int comp2 = comps.find(e_it->to);\r
+        if (comp1 != comp2 && (comps.size[comp1] < minsize || comps.size[comp2] < minsize))\r
+            comps.merge(comp1, comp2);\r
+    }\r
+\r
+    LOG2("excludesmall:", clock() - start);\r
+    DBG(start = clock());\r
+\r
+    // Compute sum of the pixel's colors which are in the same segment\r
+    Mat h_src = src;\r
+    vector<Vec4i> sumcols(nrows * ncols, Vec4i(0, 0, 0, 0));\r
+    for (int y = 0; y < nrows; ++y)\r
+    {\r
+        Vec4b* h_srcy = h_src.ptr<Vec4b>(y);\r
+        for (int x = 0; x < ncols; ++x)\r
+        {\r
+            int parent = comps.find(PIX(y, x));\r
+            Vec4b col = h_srcy[x];\r
+            Vec4i& sumcol = sumcols[parent];\r
+            sumcol[0] += col[0];\r
+            sumcol[1] += col[1];\r
+            sumcol[2] += col[2];\r
+        }\r
+    }\r
+\r
+    LOG2("computesum:", clock() - start);\r
+    DBG(start = clock());\r
+\r
+    // Create final image, color of each segment is the average color of its pixels\r
+    dst.create(src.size(), src.type());\r
+\r
+    for (int y = 0; y < nrows; ++y)\r
+    {\r
+        Vec4b* dsty = dst.ptr<Vec4b>(y);\r
+        for (int x = 0; x < ncols; ++x)\r
+        {\r
+            int parent = comps.find(PIX(y, x));\r
+            const Vec4i& sumcol = sumcols[parent];\r
+            Vec4b& dstcol = dsty[x];\r
+            dstcol[0] = static_cast<uchar>(sumcol[0] / comps.size[parent]);\r
+            dstcol[1] = static_cast<uchar>(sumcol[1] / comps.size[parent]);\r
+            dstcol[2] = static_cast<uchar>(sumcol[2] / comps.size[parent]);\r
+        }\r
+    }\r
+\r
+    LOG2("createfinal:", clock() - start);\r
+}\r
+\r
+} // namespace gpu\r
+} // namespace cv\r
+\r
+#endif // #if !defined (HAVE_CUDA)
\ No newline at end of file

diff --git a/tests/gpu/src/mssegmentation.cpp b/tests/gpu/src/mssegmentation.cpp
new file mode 100644 (file)
index 0000000..478de07
--- /dev/null
+++ b/tests/gpu/src/mssegmentation.cpp
@@ -0,0 +1,103 @@
+/*M///////////////////////////////////////////////////////////////////////////////////////\r
+//\r
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.\r
+//\r
+//  By downloading, copying, installing or using the software you agree to this license.\r
+//  If you do not agree to this license, do not download, install,\r
+//  copy or use the software.\r
+//\r
+//\r
+//                        Intel License Agreement\r
+//                For Open Source Computer Vision Library\r
+//\r
+// Copyright (C) 2000, Intel Corporation, all rights reserved.\r
+// Third party copyrights are property of their respective owners.\r
+//\r
+// Redistribution and use in source and binary forms, with or without modification,\r
+// are permitted provided that the following conditions are met:\r
+//\r
+//   * Redistribution's of source code must retain the above copyright notice,\r
+//     this list of conditions and the following disclaimer.\r
+//\r
+//   * Redistribution's in binary form must reproduce the above copyright notice,\r
+//     this list of conditions and the following disclaimer in the documentation\r
+//     and/or other materials provided with the distribution.\r
+//\r
+//   * The name of Intel Corporation may not be used to endorse or promote products\r
+//     derived from this software without specific prior written permission.\r
+//\r
+// This software is provided by the copyright holders and contributors "as is" and\r
+// any express or implied warranties, including, but not limited to, the implied\r
+// warranties of merchantability and fitness for a particular purpose are disclaimed.\r
+// In no event shall the Intel Corporation or contributors be liable for any direct,\r
+// indirect, incidental, special, exemplary, or consequential damages\r
+// (including, but not limited to, procurement of substitute goods or services;\r
+// loss of use, data, or profits; or business interruption) however caused\r
+// and on any theory of liability, whether in contract, strict liability,\r
+// or tort (including negligence or otherwise) arising in any way out of\r
+// the use of this software, even if advised of the possibility of such damage.\r
+//\r
+//M*/\r
+\r
+#include <opencv2/opencv.hpp>\r
+#include <opencv2/gpu/gpu.hpp>\r
+#include <iostream>\r
+#include <string>\r
+#include <iosfwd>\r
+#include "gputest.hpp"\r
+using namespace cv;\r
+using namespace cv::gpu;\r
+using namespace std;\r
+\r
+struct CV_GpuMeanShiftSegmentationTest : public CvTest {\r
+    CV_GpuMeanShiftSegmentationTest() : CvTest( "GPU-MeanShiftSegmentation", "MeanShiftSegmentation" ) {}\r
+\r
+    void run(int) \r
+    {\r
+        try \r
+        {\r
+            Mat img_rgb = imread(string(ts->get_data_path()) + "meanshift/cones.png");\r
+            if (img_rgb.empty())\r
+            {\r
+                ts->set_failed_test_info(CvTS::FAIL_MISSING_TEST_DATA);\r
+                return;\r
+            }\r
+\r
+            Mat img;\r
+            cvtColor(img_rgb, img, CV_BGR2BGRA);\r
+\r
+            for (int minsize = 0; minsize < 2000; minsize = (minsize + 1) * 2) \r
+            {\r
+                stringstream path;\r
+                path << ts->get_data_path() << "meanshift/cones_segmented_sp10_sr10_minsize" << minsize << ".png";\r
+\r
+                Mat dst;\r
+                meanShiftSegmentation((GpuMat)img, dst, 10, 10, minsize);\r
+                Mat dst_rgb;\r
+                cvtColor(dst, dst_rgb, CV_BGRA2BGR);\r
+\r
+                //imwrite(path.str(), dst_rgb);\r
+                Mat dst_ref = imread(path.str());\r
+                if (dst_ref.empty()) \r
+                {\r
+                    ts->set_failed_test_info(CvTS::FAIL_MISSING_TEST_DATA);\r
+                    return;\r
+                }\r
+                if (abs(cv::norm(dst_rgb - dst_ref, NORM_INF)) > 1e-3) \r
+                {\r
+                    ts->printf(CvTS::LOG, "\ndiffers from image *minsize%d.png\n", minsize);\r
+                    ts->set_failed_test_info(CvTS::FAIL_BAD_ACCURACY);\r
+                    return;\r
+                }\r
+            }\r
+        }\r
+        catch (const cv::Exception& e) \r
+        {\r
+            if (!check_and_treat_gpu_exception(e, ts))\r
+                throw;\r
+            return;\r
+        }\r
+\r
+        ts->set_failed_test_info(CvTS::OK);\r
+    }    \r
+} ms_segm_test;
\ No newline at end of file