mssegmentation.cpp

/*M///////////////////////////////////////////////////////////////////////////////////////
//
//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
//
//  By downloading, copying, installing or using the software you agree to this license.
//  If you do not agree to this license, do not download, install,
//  copy or use the software.
//
//
//                           License Agreement
//                For Open Source Computer Vision Library
//
// Copyright (C) 2010-2012, Institute Of Software Chinese Academy Of Science, all rights reserved.
// Copyright (C) 2010-2012, Advanced Micro Devices, Inc., all rights reserved.
// Third party copyrights are property of their respective owners.
//
// @Authors
//
// Redistribution and use in source and binary forms, with or without modification,
// are permitted provided that the following conditions are met:
//
//   * Redistribution's of source code must retain the above copyright notice,
//     this list of conditions and the following disclaimer.
//
//   * Redistribution's in binary form must reproduce the above copyright notice,
//     this list of conditions and the following disclaimer in the documentation
//     and/or other oclMaterials provided with the distribution.
//
//   * The name of the copyright holders may not be used to endorse or promote products
//     derived from this software without specific prior written permission.
//
// This software is provided by the copyright holders and contributors "as is" and
// any express or implied warranties, including, but not limited to, the implied
// warranties of merchantability and fitness for a particular purpose are disclaimed.
// In no event shall the Intel Corporation or contributors be liable for any direct,
// indirect, incidental, special, exemplary, or consequential damages
// (including, but not limited to, procurement of substitute goods or services;
// loss of use, data, or profits; or business interruption) however caused
// and on any theory of liability, whether in contract, strict liability,
// or tort (including negligence or otherwise) arising in any way out of
// the use of this software, even if advised of the possibility of such damage.
//
//M*/

#include "precomp.hpp"

using namespace std;

// Auxiliray stuff
namespace
{

    //
    // Declarations
    //

    class DjSets
    {
    public:
        DjSets(int n);
        int find(int elem);
        int merge(int set1, int set2);

        vector<int> parent;
        vector<int> rank;
        vector<int> size;
    private:
        DjSets(const DjSets &) {}
        DjSets operator =(const DjSets &);
    };

    template <typename T>
    struct GraphEdge
    {
        GraphEdge() {}
        GraphEdge(int to, int next, const T &val) : to(to), next(next), val(val) {}
        int to;
        int next;
        T val;
    };


    template <typename T>
    class Graph
    {
    public:
        typedef GraphEdge<T> Edge;

        Graph(int numv, int nume_max);

        void addEdge(int from, int to, const T &val = T());

        vector<int> start;
        vector<Edge> edges;

        int numv;
        int nume_max;
        int nume;
    private:
        Graph(const Graph &) {}
        Graph operator =(const Graph &) {}
    };


    struct SegmLinkVal
    {
        SegmLinkVal() {}
        SegmLinkVal(int dr, int dsp) : dr(dr), dsp(dsp) {}
        bool operator <(const SegmLinkVal &other) const
        {
            return dr + dsp < other.dr + other.dsp;
        }
        int dr;
        int dsp;
    };


    struct SegmLink
    {
        SegmLink() {}
        SegmLink(int from, int to, const SegmLinkVal &val)
            : from(from), to(to), val(val) {}
        bool operator <(const SegmLink &other) const
        {
            return val < other.val;
        }
        int from;
        int to;
        SegmLinkVal val;
    };

    //
    // Implementation
    //

    DjSets DjSets::operator = (const DjSets &/*obj*/)
    {
        //cout << "Invalid DjSets constructor\n";
        CV_Error(-1, "Invalid DjSets constructor\n");
        return *this;
    }

    DjSets::DjSets(int n) : parent(n), rank(n, 0), size(n, 1)
    {
        for (int i = 0; i < n; ++i)
            parent[i] = i;
    }


    inline int DjSets::find(int elem)
    {
        int set = elem;
        while (set != parent[set])
            set = parent[set];
        while (elem != parent[elem])
        {
            int next = parent[elem];
            parent[elem] = set;
            elem = next;
        }
        return set;
    }


    inline int DjSets::merge(int set1, int set2)
    {
        if (rank[set1] < rank[set2])
        {
            parent[set1] = set2;
            size[set2] += size[set1];
            return set2;
        }
        if (rank[set2] < rank[set1])
        {
            parent[set2] = set1;
            size[set1] += size[set2];
            return set1;
        }
        parent[set1] = set2;
        rank[set2]++;
        size[set2] += size[set1];
        return set2;
    }


    template <typename T>
    Graph<T>::Graph(int numv, int nume_max) : start(numv, -1), edges(nume_max)
    {
        this->numv = numv;
        this->nume_max = nume_max;
        nume = 0;
    }


    template <typename T>
    inline void Graph<T>::addEdge(int from, int to, const T &val)
    {
        edges[nume] = Edge(to, start[from], val);
        start[from] = nume;
        nume++;
    }


    inline int pix(int y, int x, int ncols)
    {
        return y * ncols + x;
    }


    inline int sqr(int x)
    {
        return x * x;
    }


    inline int dist2(const cv::Vec4b &lhs, const cv::Vec4b &rhs)
    {
        return sqr(lhs[0] - rhs[0]) + sqr(lhs[1] - rhs[1]) + sqr(lhs[2] - rhs[2]);
    }


    inline int dist2(const cv::Vec2s &lhs, const cv::Vec2s &rhs)
    {
        return sqr(lhs[0] - rhs[0]) + sqr(lhs[1] - rhs[1]);
    }

} // anonymous namespace

namespace cv
{
    namespace ocl
    {

        void meanShiftSegmentation(const oclMat &src, Mat &dst, int sp, int sr, int minsize, TermCriteria criteria)
        {
            CV_Assert(src.type() == CV_8UC4);
            const int nrows = src.rows;
            const int ncols = src.cols;
            const int hr = sr;
            const int hsp = sp;

            // Perform mean shift procedure and obtain region and spatial maps
            oclMat h_rmap, h_spmap;
            meanShiftProc(src, h_rmap, h_spmap, sp, sr, criteria);
            Mat rmap = h_rmap;
            Mat spmap = h_spmap;

            Graph<SegmLinkVal> g(nrows * ncols, 4 * (nrows - 1) * (ncols - 1)
                                 + (nrows - 1) + (ncols - 1));

            // Make region adjacent graph from image
            Vec4b r1;
            Vec4b r2[4];
            Vec2s sp1;
            Vec2s sp2[4];
            int dr[4];
            int dsp[4];
            for (int y = 0; y < nrows - 1; ++y)
            {
                Vec4b *ry = rmap.ptr<Vec4b>(y);
                Vec4b *ryp = rmap.ptr<Vec4b>(y + 1);
                Vec2s *spy = spmap.ptr<Vec2s>(y);
                Vec2s *spyp = spmap.ptr<Vec2s>(y + 1);
                for (int x = 0; x < ncols - 1; ++x)
                {
                    r1 = ry[x];
                    sp1 = spy[x];

                    r2[0] = ry[x + 1];
                    r2[1] = ryp[x];
                    r2[2] = ryp[x + 1];
                    r2[3] = ryp[x];

                    sp2[0] = spy[x + 1];
                    sp2[1] = spyp[x];
                    sp2[2] = spyp[x + 1];
                    sp2[3] = spyp[x];

                    dr[0] = dist2(r1, r2[0]);
                    dr[1] = dist2(r1, r2[1]);
                    dr[2] = dist2(r1, r2[2]);
                    dsp[0] = dist2(sp1, sp2[0]);
                    dsp[1] = dist2(sp1, sp2[1]);
                    dsp[2] = dist2(sp1, sp2[2]);

                    r1 = ry[x + 1];
                    sp1 = spy[x + 1];

                    dr[3] = dist2(r1, r2[3]);
                    dsp[3] = dist2(sp1, sp2[3]);

                    g.addEdge(pix(y, x, ncols), pix(y, x + 1, ncols), SegmLinkVal(dr[0], dsp[0]));
                    g.addEdge(pix(y, x, ncols), pix(y + 1, x, ncols), SegmLinkVal(dr[1], dsp[1]));
                    g.addEdge(pix(y, x, ncols), pix(y + 1, x + 1, ncols), SegmLinkVal(dr[2], dsp[2]));
                    g.addEdge(pix(y, x + 1, ncols), pix(y + 1, x, ncols), SegmLinkVal(dr[3], dsp[3]));
                }
            }
            for (int y = 0; y < nrows - 1; ++y)
            {
                r1 = rmap.at<Vec4b>(y, ncols - 1);
                r2[0] = rmap.at<Vec4b>(y + 1, ncols - 1);
                sp1 = spmap.at<Vec2s>(y, ncols - 1);
                sp2[0] = spmap.at<Vec2s>(y + 1, ncols - 1);
                dr[0] = dist2(r1, r2[0]);
                dsp[0] = dist2(sp1, sp2[0]);
                g.addEdge(pix(y, ncols - 1, ncols), pix(y + 1, ncols - 1, ncols), SegmLinkVal(dr[0], dsp[0]));
            }
            for (int x = 0; x < ncols - 1; ++x)
            {
                r1 = rmap.at<Vec4b>(nrows - 1, x);
                r2[0] = rmap.at<Vec4b>(nrows - 1, x + 1);
                sp1 = spmap.at<Vec2s>(nrows - 1, x);
                sp2[0] = spmap.at<Vec2s>(nrows - 1, x + 1);
                dr[0] = dist2(r1, r2[0]);
                dsp[0] = dist2(sp1, sp2[0]);
                g.addEdge(pix(nrows - 1, x, ncols), pix(nrows - 1, x + 1, ncols), SegmLinkVal(dr[0], dsp[0]));
            }

            DjSets comps(g.numv);

            // Find adjacent components
            for (int v = 0; v < g.numv; ++v)
            {
                for (int e_it = g.start[v]; e_it != -1; e_it = g.edges[e_it].next)
                {
                    int c1 = comps.find(v);
                    int c2 = comps.find(g.edges[e_it].to);
                    if (c1 != c2 && g.edges[e_it].val.dr < hr && g.edges[e_it].val.dsp < hsp)
                        comps.merge(c1, c2);
                }
            }

            vector<SegmLink> edges;
            edges.reserve(g.numv);

            // Prepare edges connecting differnet components
            for (int v = 0; v < g.numv; ++v)
            {
                int c1 = comps.find(v);
                for (int e_it = g.start[v]; e_it != -1; e_it = g.edges[e_it].next)
                {
                    int c2 = comps.find(g.edges[e_it].to);
                    if (c1 != c2)
                        edges.push_back(SegmLink(c1, c2, g.edges[e_it].val));
                }
            }

            // Sort all graph's edges connecting differnet components (in asceding order)
            sort(edges.begin(), edges.end());

            // Exclude small components (starting from the nearest couple)
            for (size_t i = 0; i < edges.size(); ++i)
            {
                int c1 = comps.find(edges[i].from);
                int c2 = comps.find(edges[i].to);
                if (c1 != c2 && (comps.size[c1] < minsize || comps.size[c2] < minsize))
                    comps.merge(c1, c2);
            }

            // Compute sum of the pixel's colors which are in the same segment
            Mat h_src = src;
            vector<Vec4i> sumcols(nrows * ncols, Vec4i(0, 0, 0, 0));
            for (int y = 0; y < nrows; ++y)
            {
                Vec4b *h_srcy = h_src.ptr<Vec4b>(y);
                for (int x = 0; x < ncols; ++x)
                {
                    int parent = comps.find(pix(y, x, ncols));
                    Vec4b col = h_srcy[x];
                    Vec4i &sumcol = sumcols[parent];
                    sumcol[0] += col[0];
                    sumcol[1] += col[1];
                    sumcol[2] += col[2];
                }
            }

            // Create final image, color of each segment is the average color of its pixels
            dst.create(src.size(), src.type());

            for (int y = 0; y < nrows; ++y)
            {
                Vec4b *dsty = dst.ptr<Vec4b>(y);
                for (int x = 0; x < ncols; ++x)
                {
                    int parent = comps.find(pix(y, x, ncols));
                    const Vec4i &sumcol = sumcols[parent];
                    Vec4b &dstcol = dsty[x];
                    dstcol[0] = static_cast<uchar>(sumcol[0] / comps.size[parent]);
                    dstcol[1] = static_cast<uchar>(sumcol[1] / comps.size[parent]);
                    dstcol[2] = static_cast<uchar>(sumcol[2] / comps.size[parent]);
                }
            }
        }

    }
}