OpenCV C++/IFでボロノイ図・ドロネー図を描く

ボロノイ図・ドロネー図の作成については、大昔からOpenCV C/IFに有りました。C++には無いと思っていたのですが、少し前にいつの間にか有ることを知りました。

Cのボロノイ図まわりのAPIは、正直いまだにあまり理解していません。いまひとつ、つかみどころがない感じです。
http://opencv.jp/opencv-2svn/py/planar_subdivisions.html

しかしC++のI/Fはわりとわかりやすくなっています。これを使って、ボロノイ図・ドロネー図を描画するサンプルを残しておきます。

cv::Subdiv2Dの定義

publicのところだけ、opencv2/imgproc/imgproc.hppから拝借。マニュアルはまだ無いようです。

class CV_EXPORTS_W Subdiv2D
{
public:
    enum
    {
        PTLOC_ERROR = -2,
        PTLOC_OUTSIDE_RECT = -1,
        PTLOC_INSIDE = 0,
        PTLOC_VERTEX = 1,
        PTLOC_ON_EDGE = 2
    };

    enum
    {
        NEXT_AROUND_ORG   = 0x00,
        NEXT_AROUND_DST   = 0x22,
        PREV_AROUND_ORG   = 0x11,
        PREV_AROUND_DST   = 0x33,
        NEXT_AROUND_LEFT  = 0x13,
        NEXT_AROUND_RIGHT = 0x31,
        PREV_AROUND_LEFT  = 0x20,
        PREV_AROUND_RIGHT = 0x02
    };

    CV_WRAP Subdiv2D();
    CV_WRAP Subdiv2D(Rect rect);
    CV_WRAP void initDelaunay(Rect rect);

    CV_WRAP int insert(Point2f pt);
    CV_WRAP void insert(const vector<Point2f>& ptvec);
    CV_WRAP int locate(Point2f pt, CV_OUT int& edge, CV_OUT int& vertex);

    CV_WRAP int findNearest(Point2f pt, CV_OUT Point2f* nearestPt=0);
    CV_WRAP void getEdgeList(CV_OUT vector<Vec4f>& edgeList) const;
    CV_WRAP void getTriangleList(CV_OUT vector<Vec6f>& triangleList) const;
    CV_WRAP void getVoronoiFacetList(const vector<int>& idx, CV_OUT vector<vector<Point2f> >& facetList,
                                     CV_OUT vector<Point2f>& facetCenters);

    CV_WRAP Point2f getVertex(int vertex, CV_OUT int* firstEdge=0) const;

    CV_WRAP int getEdge( int edge, int nextEdgeType ) const;
    CV_WRAP int nextEdge(int edge) const;
    CV_WRAP int rotateEdge(int edge, int rotate) const;
    CV_WRAP int symEdge(int edge) const;
    CV_WRAP int edgeOrg(int edge, CV_OUT Point2f* orgpt=0) const;
    CV_WRAP int edgeDst(int edge, CV_OUT Point2f* dstpt=0) const;

protected:
    // ....
};

入力の用意

ここから実装開始。
例によってランダムに点を打ちまくります(0～600の範囲)。Subdiv2Dの都合で、cv::Point2fとして点を作っていきます。

std::mt19937 mt;
std::uniform_int_distribution<int> distribution(0, 599);
std::vector<cv::Point2f> points;
for(int i = 0; i < 100; i++)
{
    int x = distribution(mt);
    int y = distribution(mt);
    points.push_back(cv::Point2f(x, y));
}

乱数の発生には、最近の新しいC++標準で追加された<random>ヘッダにあるメルセンヌツイスタ(std::mt19937)を使用しています。OpenCVにも cv::RNG_MT19937 というクラスがあるので、そちらを使っても良いです。

Subdiv2D初期化

cv::Subdiv2D subdiv;
subdiv.initDelaunay(cv::Rect(0, 0, 600, 600));
subdiv.insert(points);

insert関数は、点を1つ入れるものと配列で一気に複数入れるものの2つのオーバーロードがあります。一気に入れた方が高速かなと思ったのですが、体感ではあまり違いはわかりませんでした。

ボロノイ図を取得

getVoronoiFacetListを使います。

std::vector<int> idx;
std::vector<std::vector<cv::Point2f>> facetLists;
std::vector<cv::Point2f> facetCenters;
subdiv.getVoronoiFacetList(idx, facetLists, facetCenters);

facetListsが「ボロノイ領域の頂点リストのリスト」、facetCentersがボロノイ領域の母点リストを表します。

idx（取りたいボロノイ領域のIDリスト）が困りました。そもそもIDって何だよ、どこからわかるんだよ、って感じでした。しかし、空のvectorを指定すれば全部取れるようです。それで全然オッケーです。

ボロノイ図の描画

ランダムに打った点とあわせて、ボロノイ図を描画します。

// ランダム点の描画
cv::Mat img = cv::Mat::zeros(600, 600, CV_8UC3);		
for(auto pt = points.begin(); pt != points.end(); pt++)
{
    cv::circle(img, *pt, 4, cv::Scalar(0,0,255), -1);			
}

// Voronoi図を描画
for(auto list = facetLists.begin(); list != facetLists.end(); list++)
{
    cv::Point2f before = list->back();
    for(auto pt = list->begin(); pt != list->end(); pt++)
    {
        cv::Point p1((int)before.x, (int)before.y);
	cv::Point p2((int)pt->x, (int)pt->y);
	cv::line(img, p1, p2, cv::Scalar(64,255,128));
	before = *pt;
    }
}

cv::namedWindow("voronoi");
cv::imshow("voronoi", img);
cv::waitKey();

ドロネー三角形分割の取得・描画

以下2つ、どちらでも描けます。

getEdgeList

ドロネー図を構成する辺のリストを得ることができます。辺は4次元のベクトル(Vec4f)として得られ、(x1, y1, x2, y2)の順で並んでいます。

// 辺のリストを取得
std::vector<cv::Vec4f> edgeList;
subdiv.getEdgeList(edgeList);

// 描画
for(auto edge = edgeList.begin(); edge != edgeList.end(); edge++)
{
    cv::Point p1(edge->val[0], edge->val[1]);
    cv::Point p2(edge->val[2], edge->val[3]);
    cv::line(imgEdges, p1, p2, cv::Scalar(48,128,48));
}

getTriangleList

こちらの関数では、三角形の3頂点がVec6fとしてまとまって得られます。getEdgeList同様に、それをばらして解釈していくことになります。若干コード量が多くなるのと、同じ線を2回描画することになりますが、基本は同じです。

// ドロネー三角形のリストを取得
std::vector<cv::Vec6f> triangles;
subdiv.getTriangleList(triangles);

// 描画
for(auto it = triangles.begin(); it != triangles.end(); it++)
{
    cv::Vec6f &vec = *it;
    cv::Point p1(vec[0], vec[1]);
    cv::Point p2(vec[2], vec[3]);
    cv::Point p3(vec[4], vec[5]);
    cv::line(img, p1, p2, cv::Scalar(64,255,128));
    cv::line(img, p2, p3, cv::Scalar(64,255,128));
    cv::line(img, p3, p1, cv::Scalar(64,255,128));
}

かなりわかりやすいですね。おかげさまでドロネー恐怖症がかなり直りました。