Having a matrix MxN of integers how to group them into polygons with boost geometry?

Question

We have a matrix of given integers (any from 1 to INT_MAX) like

1 2 3
1 3 3
1 3 3
100 2 1

We want to create polygons with same colors for each unique int in matrix so our polygons would have coords/groupings like shown here.

and we could generate image like this:

Which *(because of vectirisation that was performed would scale to such size like):

(sorry for crappy drawings)

Is it possible and how to do such thing with boost geometry?

Update:

So @sehe sad: I'd simply let Boost Geometry do most of the work. so I created this pixel by pixel class aeria grower using purely Boost.Geometry, compiles, runs but I need it to run on clustered data.. and I have 1000 by 1800 files of uchars (each unique uchar == data belongs to that claster). Problem with this code: on 18th line it gets SO WARY SLOW that each point creation starts to take more than one second=(

code:

//Boost
#include <boost/assign.hpp>
#include <boost/foreach.hpp>
#include <boost/shared_ptr.hpp>
#include <boost/algorithm/string.hpp>

#include <boost/geometry/geometry.hpp>
#include <boost/geometry/geometries/geometries.hpp>
#include <boost/geometry/multi/geometries/multi_polygon.hpp>
#include <boost/geometry/geometries/adapted/boost_tuple.hpp>

//and this is why we use Boost Geometry from Boost trunk 
//#include <boost/geometry/extensions/io/svg/svg_mapper.hpp>

BOOST_GEOMETRY_REGISTER_BOOST_TUPLE_CS(cs::cartesian)


void make_point(int x, int y,  boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<double> > & ring)
{
    using namespace boost::assign;

    boost::geometry::append(  ring,     boost::geometry::model::d2::point_xy<double>(x-1, y-1));
    boost::geometry::append(  ring,     boost::geometry::model::d2::point_xy<double>(x, y-1));
    boost::geometry::append(  ring,      boost::geometry::model::d2::point_xy<double>(x, y));
    boost::geometry::append(  ring,      boost::geometry::model::d2::point_xy<double>(x-1, y));
    boost::geometry::append(  ring,     boost::geometry::model::d2::point_xy<double>(x-1, y-1));
    boost::geometry::correct(ring);
}

void create_point(int x, int y, boost::geometry::model::multi_polygon< boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<double> > > & mp)
{
    boost::geometry::model::multi_polygon< boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<double> > > temp;
    boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<double> > ring;
    make_point(x, y, ring);
    boost::geometry::union_(mp, ring, temp);
    boost::geometry::correct(temp);
    mp=temp;
}

int main()
{
    using namespace boost::assign;


    boost::geometry::model::multi_polygon< boost::geometry::model::polygon < boost::geometry::model::d2::point_xy<double> > > pol, simpl;

    //read image
    std::ifstream in("1.mask", std::ios_base::in | std::ios_base::binary);
    int sx, sy;
    in.read(reinterpret_cast<char*>(&sy), sizeof(int));
    in.read(reinterpret_cast<char*>(&sx), sizeof(int));
    std::vector< std::vector<unsigned char> > image(sy);

    for(int i =1; i <= sy; i++)
    {
        std::vector<unsigned char> row(sx);
        in.read(reinterpret_cast<char*>(&row[0]), sx);
        image[i-1] = row;
    }

    //
    std::map<unsigned char,  boost::geometry::model::multi_polygon < boost::geometry::model::polygon < boost::geometry::model::d2::point_xy<double> > >  > layered_image;

    for(int y=1; y <= sy; y++)
    {
        for(int x=1; x <= sx; x++)
        {
            if (image[y-1][x-1] != 1)
            {
                create_point(x, y, layered_image[image[y-1][x-1]]);
                std::cout << x << " : " << y << std::endl;
            }
        }
    }
}

So as you can see my code suks.. so I decided to create a renderer for @sehe code:

#include <iostream>
#include <fstream>

#include <vector>
#include <string>
#include <set>

//Boost
#include <boost/assign.hpp>
#include <boost/array.hpp>
#include <boost/algorithm/string.hpp>
#include <boost/geometry/geometry.hpp>
#include <boost/geometry/geometries/geometries.hpp>
#include <boost/geometry/multi/geometries/multi_polygon.hpp>
#include <boost/geometry/geometries/adapted/boost_tuple.hpp>
#include <boost/random/uniform_int.hpp>
#include <boost/random/mersenne_twister.hpp>


//and this is why we use Boost Geometry from Boost trunk 
#include <boost/geometry/extensions/io/svg/svg_mapper.hpp>

BOOST_GEOMETRY_REGISTER_BOOST_TUPLE_CS(cs::cartesian)

    namespace mxdetail
{
    typedef size_t cell_id; // row * COLS + col

    template <typename T> struct area
    {
        T value;
        typedef std::vector<cell_id> cells_t;
        cells_t cells;
    };

    template <typename T, size_t Rows, size_t Cols>
    std::vector<area<T> > getareas(const boost::array<boost::array<T, Cols>, Rows>& matrix)
    {
        typedef boost::array<boost::array<T, Cols>, Rows> mtx;
        std::vector<area<T> > areas;

        struct visitor_t
        {
            const mtx& matrix;
            std::set<cell_id> visited;

            visitor_t(const mtx& mtx) : matrix(mtx) { }

            area<T> start(const int row, const int col)
            {
                area<T> result;
                visit(row, col, result);
                return result;
            }

            void visit(const int row, const int col, area<T>& current)
            {
                const cell_id id = row*Cols+col;
                if (visited.end() != visited.find(id))
                    return;

                bool matches = current.cells.empty() || (matrix[row][col] == current.value);

                if (matches)
                {
                    visited.insert(id);
                    current.value = matrix[row][col];
                    current.cells.push_back(id);

                    // process neighbours
                    for (int nrow=std::max(0, row-1); nrow < std::min((int) Rows, row+2); nrow++)
                        for (int ncol=std::max(0, col-1); ncol < std::min((int) Cols, col+2); ncol++)
                            /* if (ncol!=col || nrow!=row) */
                            visit(nrow, ncol, current);
                }
            }
        } visitor(matrix);

        for (int r=0; r < (int) Rows; r++)
            for (int c=0; c < (int) Cols; c++)
            {
                mxdetail::area<int> area = visitor.start(r,c);
                if (!area.cells.empty()) // happens when startpoint already visited
                    areas.push_back(area);
            }

            return areas;
    }
}




typedef boost::array<int, 4> row;

template <typename T, size_t N>
boost::array<T, N> make_array(const T (&a)[N])
{
    boost::array<T, N> result;
    std::copy(a, a+N, result.begin());
    return result;
}


void make_point(int x, int y,  boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<double> > & ring)
{
    using namespace boost::assign;

    boost::geometry::append(  ring,     boost::geometry::model::d2::point_xy<double>(x-1, y-1));
    boost::geometry::append(  ring,     boost::geometry::model::d2::point_xy<double>(x, y-1));
    boost::geometry::append(  ring,      boost::geometry::model::d2::point_xy<double>(x, y));
    boost::geometry::append(  ring,      boost::geometry::model::d2::point_xy<double>(x-1, y));
    boost::geometry::append(  ring,     boost::geometry::model::d2::point_xy<double>(x-1, y-1));
    boost::geometry::correct(ring);
}

void create_point(int x, int y, boost::geometry::model::multi_polygon< boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<double> > > & mp)
{
    boost::geometry::model::multi_polygon< boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<double> > > temp;
    boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<double> > ring;
    make_point(x, y, ring);
    boost::geometry::union_(mp, ring, temp);
    boost::geometry::correct(temp);
    mp=temp;
}
boost::random::mt19937 rng;  
boost::random::uniform_int_distribution<> color(10,255);

std::string fill_rule()
{ 
    int red, green, blue;
    red = color(rng);
    green = color(rng);
    blue = color(rng);

    std::ostringstream rule;
    rule << "fill-rule:nonzero;fill-opacity:0.5;fill:rgb("
        << red  << ","  << green << "," << blue
        << ");stroke:rgb("
        << (red - 5) << "," << (green - 5) << "," << (blue -5) 
        <<  ");stroke-width:2";
    return rule.str();
}

int main()
{
    int sx = 4;
    int sy = 5;

    int row0[] = { 1  , 2, 3, 3, };
    int row1[] = { 1  , 3, 3, 3,};
    int row2[] = { 1  , 3, 3, 3, };
    int row3[] = { 2  , 2, 1, 2, };
    int row4[] = { 100, 2, 2, 2, };

    boost::array<row, 5> matrix;
    matrix[0] = make_array(row0);
    matrix[1] = make_array(row1);
    matrix[2] = make_array(row2);
    matrix[3] = make_array(row3);
    matrix[4] = make_array(row4);

    typedef std::vector<mxdetail::area<int> > areas_t;
    typedef areas_t::value_type::cells_t cells_t; 

    areas_t areas = mxdetail::getareas(matrix);

    using namespace boost::assign;

    typedef boost::geometry::model::polygon
        <
        boost::geometry::model::d2::point_xy<double>
        > polygon;

    typedef boost::geometry::model::multi_polygon<polygon> mp;

    typedef boost::geometry::point_type<mp>::type point_type;

    std::string filename = "draw.svg";
    std::ofstream svg(filename.c_str());


    boost::geometry::svg_mapper<point_type> mapper(svg, 400, 400);
    for (areas_t::const_iterator it=areas.begin(); it!=areas.end(); ++it)
    {
        mp pol;
        std::cout << "area of " << it->value << ": ";
        for (cells_t::const_iterator pit=it->cells.begin(); pit!=it->cells.end(); ++pit)
        {
            int row = *pit / 3, col = *pit % 3;
            std::cout << "(" << row << "," << col << "), ";
            create_point( (row+1), (col+1), pol);
        }

        std::cout << std::endl;
        mapper.add(pol);
        mapper.map(pol, fill_rule());
    }
    std::cout << "areas detected: " << areas.size() << std::endl;
    std::cin.get();
}

this code is compilable but it sucks (seems I did not get how to work with arrays after all...):

Answer 1

In short, if I got the question right, I'd simply let Boost Geometry do most of the work.

For a sample matrix of NxM, create NxM 'flyweight' rectangle polygons to correspond to each matrix cell visually.

Now, using an iterative deepening algorithm, find all groups:

* for each _unvisited_ cell in matrix
  * start a new group
  * [visit:] 
     - mark _visited_
     - for each neighbour with equal value: 
          - add to curent group and
          - recurse [visit:]

Note that the result of this algorithm could be distinct groups with the same values (representing disjunct polygons). E.g. the value 2 from the sample in the OP would result in two groups.

Now for each group, simply call Boost Geometry's Union_ algorithm to find the consolidated polygon to represent that group.

Sample implementation

Update Here is a non-optimized implementation in C++11:

Edit See here for C++03 version (using Boost)

The sample data used in the test corresponds to the matrix from the question.

#include <iostream>
#include <array>
#include <vector>
#include <set>

namespace mxdetail
{
    typedef size_t cell_id; // row * COLS + col

    template <typename T> struct area
    {
        T value;
        std::vector<cell_id> cells;
    };

    template <typename T, size_t Rows, size_t Cols>
        std::vector<area<T> > getareas(const std::array<std::array<T, Cols>, Rows>& matrix)
    {
        typedef std::array<std::array<T, Cols>, Rows> mtx;
        std::vector<area<T> > areas;

        struct visitor_t
        {
            const mtx& matrix;
            std::set<cell_id> visited;

            visitor_t(const mtx& mtx) : matrix(mtx) { }

            area<T> start(const int row, const int col)
            {
                area<T> result;
                visit(row, col, result);
                return result;
            }

            void visit(const int row, const int col, area<T>& current)
            {
                const cell_id id = row*Cols+col;
                if (visited.end() != visited.find(id))
                    return;

                bool matches = current.cells.empty() || (matrix[row][col] == current.value);

                if (matches)
                {
                    visited.insert(id);
                    current.value = matrix[row][col];
                    current.cells.push_back(id);

                    // process neighbours
                    for (int nrow=std::max(0, row-1); nrow < std::min((int) Rows, row+2); nrow++)
                    for (int ncol=std::max(0, col-1); ncol < std::min((int) Cols, col+2); ncol++)
                        /* if (ncol!=col || nrow!=row) */
                            visit(nrow, ncol, current);
                }
            }
        } visitor(matrix);

        for (int r=0; r < Rows; r++)
            for (int c=0; c < Cols; c++)
            {
                auto area = visitor.start(r,c);
                if (!area.cells.empty()) // happens when startpoint already visited
                    areas.push_back(area);
            }

        return areas;
    }
}

int main()
{
    typedef std::array<int, 3> row;
    std::array<row, 4> matrix = { 
        row { 1  , 2, 3, },
        row { 1  , 3, 3, },
        row { 1  , 3, 3, },
        row { 100, 2, 1, },
    };

    auto areas = mxdetail::getareas(matrix);

    std::cout << "areas detected: " << areas.size() << std::endl;
    for (const auto& area : areas)
    {
        std::cout << "area of " << area.value << ": ";
        for (auto pt : area.cells)
        {
            int row = pt / 3, col = pt % 3;
            std::cout << "(" << row << "," << col << "), ";
        }
        std::cout << std::endl;
    }
}

Compiled with gcc-4.6 -std=c++0x the output is:

areas detected: 6
area of 1: (0,0), (1,0), (2,0), 
area of 2: (0,1), 
area of 3: (0,2), (1,1), (1,2), (2,1), (2,2), 
area of 100: (3,0), 
area of 2: (3,1), 
area of 1: (3,2), 

Answer 2

When number of points is big (say, more than 1000x1000), the solution above would gobble a lot of memory. And this is exactly what happened to the topic-starter.

Below I show more scalable approach.

I would separate two problems here: one is to find the areas, another is to convert them into polygons.

The first problem is actually equivalent to finding the connected components of the grid graph where neighbors has edges if and only if they have equal "colors" attached to it. One can use a grid graph from boost-graph.

#include <boost/graph/grid_graph.hpp>
// Define dimension lengths, a MxN in this case
boost::array<int, 2> lengths = { { M, N } };

// Create a MxN two-dimensional, unwrapped grid graph 
grid_graph<2> graph(lengths);

Next, we should convert a given matrix M into an edge filter: grid edges are present iff the "color" of the neighbors are the same.

template <class Matrix>
struct GridEdgeFilter
{  
    typedef grid_graph<2> grid;
    GridEdgeFilter(const Matrix & m, const grid&g):_m(&m),_g(&g){}

    /// \return true iff edge is present in the graph
    bool operator()(grid::edge_descriptor e) const
    {
        grid::vertex_descriptor src = source(e,*_g), tgt = target(e,*_g);
        //src[0] is x-coord of src, etc. The value (*m)[x,y] is the color of the point (x,y).
        //Edge is preserved iff matrix values are equal
        return (*_m)[src[0],src[1]] == (*_m)[tgt[0],tgt[1]];
    }

    const Matrix * _m;
    const grid* _g;
};

Finally, we define a boost::filtered_graph of grid and EdgeFilter and call Boost.Graph algorithm for connected components.

Each connected component represents a set of points of a single color i.e. exactly the area we want to transform into a polygon.

Here we have another issue. Boost.Geometry only allows to merge polygons one by one. Hence it becomes very slow when number of polygons is big.

The better way is to use Boost.Polygon, namely its Property Merge functionality. One starts with empty property_merge object, and goes on by inserting rectangles of given color (you can set color as a property). Then one calls the method merge and gets a polygon_set as the output.