Delaunay三角形分割

はじめに

今回は、Delaunay三角形分割とその双対構造であるVoronoi図について説明する。

Delaunay三角形分割とは

いま3次元空間内で定義される関数 $z=f(x,y)$ を考える。関数の形状は未知であり、観測点 $X_i=(x_i,y_i),i=1,\cdots,N$ における値のみが既知であるとする（下図参照：ここから引用した）。

3つの観測点 $X_i,X_j,X_k$ の近傍にある未観測点 $X_*$ における $f(X_*)$ の値を予測したいとき、以下の手順を考えることができる。

点 $X_*$ を囲む、3つの観測点から構成される三角形を見つける。
この三角形の頂点における関数の値は既知である。適当な加重平均を計算し、 $f(X_*)$ を求める。

上の手順を実現するには、2次元平面を観測点 $X_i,i=1,\cdots,N$ を頂点とする三角形に分割しなければならない（下図参照：ここから引用した）。

このときよく使われる手法がDelaunayの三角形分割である。この分割手法は有限要素法を行う際のメッシュ生成にも使われている。

Delaunay三角形分割の定義

平面上の点の集合（母点集合）を $S=\{X_1,\cdots,X_N\}$ とする。 $S$ から3つの母点 $X_i,X_j,X_k$ を取り出し、これら3点のみを円周上に持ち、この3点以外の母点を円内に持たないような円が存在するとき、この3点を頂点とする三角形を作る。これを繰り返したのものがDelaunay三角形分割である（下図参照）。

この分割により、2次元平面は「ふっくら」とした三角形に分割される。

双対構造：Voronoi図

Delaunay三角形分割が実現できれば、これと双対な構造であるVoronoi図を容易に作ることができる（下図参照：ここから引用した）。

上図の赤点（母点）を結ぶ構造がDelaunay三角形分割であり、母点を細胞核のように内部にもつ領域の集まりがVoronoi図である。Voronoi図の生成手順は以下の通りである。

三角形 $A$ の外接円 $C_A$ の中心を求める。
三角形 $A,B$ が隣接するとき $C_A,C_B$ の中心を結ぶ。
これを繰り返す。

例えば母点 $X_i$ を内部に持つVoronoi領域 $V_i$ を考える。この領域内の任意の点 $x$ と全母点との距離を考えたとき、 $x$ に最も近い母点は $X_i$ である。つまり、母点 $X_i$ が「支配」する領域が $V_i$ である。ひとつの適用例としてWorld Airports Voronoiなるものがある。これは世界中の空港の位置を母点としたVoronoi図である。

実装

Delaunay三角形分割とVoronoi図はともにOpenCVに実装されているので、これを利用する。サンプルコードは以下の通りである。

#include <iostream>
#include <random>
#include <vector>
#include <opencv2/opencv.hpp>
#include <cpplinq.hpp>
#include <boost/range/irange.hpp>

// size of a window
constexpr int IMAGE_SIZE{1000};

cv::Mat initialize_canvas(const std::vector<cv::Point2f>& points);
void calculate_voronoi(cv::Subdiv2D& subdiv, const std::vector<cv::Point2f>& points, cv::Mat& canvas);
void calculate_delaunay(cv::Subdiv2D& subdiv, const std::vector<cv::Point2f>& points, cv::Mat& canvas);

int main(int argc, const char * argv[]) {
    
    // make random points
    std::mt19937 mt{};
    std::uniform_int_distribution<int> distribution{0, IMAGE_SIZE - 1};
    std::vector<cv::Point2f> points{};
    constexpr int N_POINTS{100};
    points.reserve(N_POINTS);
    for ([[maybe_unused]] auto i : boost::irange(0, N_POINTS))
    {
        auto x = distribution(mt);
        auto y = distribution(mt);
        points.emplace_back(x, y);
    }
    
    // make an object
    cv::Subdiv2D subdiv{};
    subdiv.initDelaunay(cv::Rect(0, 0, IMAGE_SIZE, IMAGE_SIZE));
    subdiv.insert(points);
  
    // make a delaunay triangulation
    auto canvas = initialize_canvas(points);
    calculate_delaunay(subdiv, points, canvas);
    cv::imwrite("/Users/kumada/Projects/cct_blog/voronoi/opencv_voronoi/opencv_voronoi/delaunay.jpg", canvas);
    
    // make a voronoi diagram
    canvas = initialize_canvas(points);
    calculate_voronoi(subdiv, points, canvas);
    cv::imwrite("/Users/kumada/Projects/cct_blog/voronoi/opencv_voronoi/opencv_voronoi/voronoi.jpg", canvas);
    
    return 0;
}

cv::Mat initialize_canvas(const std::vector<cv::Point2f>& points)
{
    // display random points
    cv::Mat canvas = cv::Mat::zeros(IMAGE_SIZE, IMAGE_SIZE, CV_8UC3);
    for (const auto& point: points)
    {
        cv::circle(canvas, point, 4, cv::Scalar(0, 0, 255), -1);
    }
    return canvas;
}

void calculate_voronoi(cv::Subdiv2D& subdiv, const std::vector<cv::Point2f>& points, cv::Mat& canvas)
{

    // prepare output containers
    std::vector<int> idx{};
    std::vector<std::vector<cv::Point2f>> facet_lists{};
    std::vector<cv::Point2f> facet_cernters{};

    // calculcate voronoi
    subdiv.getVoronoiFacetList(idx, facet_lists, facet_cernters);

    // display voronoi structure
    for (const auto& list: facet_lists)
    {
        auto before = list.back();
        for (const auto& point: list)
        {
            cv::Point p1{static_cast<int>(before.x), static_cast<int>(before.y)};
            cv::Point p2{static_cast<int>(point.x), static_cast<int>(point.y)};
            cv::line(canvas, p1, p2, cv::Scalar(64, 255, 128));
            before = std::move(point);
        }
    }
}

void calculate_delaunay(cv::Subdiv2D& subdiv, const std::vector<cv::Point2f>& points, cv::Mat& canvas)
{
    std::vector<cv::Vec4f> edge_list{};
    subdiv.getEdgeList(edge_list);
    
    for (const auto& edge: edge_list)
    {
        cv::Point p1{static_cast<int>(edge.val[0]), static_cast<int>(edge.val[1])};
        cv::Point p2{static_cast<int>(edge.val[2]), static_cast<int>(edge.val[3])};
        cv::line(canvas, p1, p2, cv::Scalar(48, 128, 48));
    }
}

#include <iostream>

#include <random>

#include <vector>

#include <opencv2/opencv.hpp>

#include <cpplinq.hpp>

#include <boost/range/irange.hpp>

// size of a window

constexpr int IMAGE_SIZE{1000};

cv::Mat initialize_canvas(const std::vector<cv::Point2f>& points);

void calculate_voronoi(cv::Subdiv2D& subdiv, const std::vector<cv::Point2f>& points, cv::Mat& canvas);

void calculate_delaunay(cv::Subdiv2D& subdiv, const std::vector<cv::Point2f>& points, cv::Mat& canvas);

int main(int argc, const char * argv[]) {

// make random points

std::mt19937 mt{};

std::uniform_int_distribution<int> distribution{0, IMAGE_SIZE - 1};

std::vector<cv::Point2f> points{};

constexpr int N_POINTS{100};

points.reserve(N_POINTS);

for ([[maybe_unused]] auto i : boost::irange(0, N_POINTS))

{

auto x = distribution(mt);

auto y = distribution(mt);

points.emplace_back(x, y);

}

// make an object

cv::Subdiv2D subdiv{};

subdiv.initDelaunay(cv::Rect(0, 0, IMAGE_SIZE, IMAGE_SIZE));

subdiv.insert(points);

// make a delaunay triangulation

auto canvas = initialize_canvas(points);

calculate_delaunay(subdiv, points, canvas);

cv::imwrite("/Users/kumada/Projects/cct_blog/voronoi/opencv_voronoi/opencv_voronoi/delaunay.jpg", canvas);

// make a voronoi diagram

canvas = initialize_canvas(points);

calculate_voronoi(subdiv, points, canvas);

cv::imwrite("/Users/kumada/Projects/cct_blog/voronoi/opencv_voronoi/opencv_voronoi/voronoi.jpg", canvas);

return 0;

}

cv::Mat initialize_canvas(const std::vector<cv::Point2f>& points)

{

// display random points

cv::Mat canvas = cv::Mat::zeros(IMAGE_SIZE, IMAGE_SIZE, CV_8UC3);

for (const auto& point: points)

{

cv::circle(canvas, point, 4, cv::Scalar(0, 0, 255), -1);

}

return canvas;

}

void calculate_voronoi(cv::Subdiv2D& subdiv, const std::vector<cv::Point2f>& points, cv::Mat& canvas)

{

// prepare output containers

std::vector<int> idx{};

std::vector<std::vector<cv::Point2f>> facet_lists{};

std::vector<cv::Point2f> facet_cernters{};

// calculcate voronoi

subdiv.getVoronoiFacetList(idx, facet_lists, facet_cernters);

// display voronoi structure

for (const auto& list: facet_lists)

{

auto before = list.back();

for (const auto& point: list)

{

cv::Point p1{static_cast<int>(before.x), static_cast<int>(before.y)};

cv::Point p2{static_cast<int>(point.x), static_cast<int>(point.y)};

cv::line(canvas, p1, p2, cv::Scalar(64, 255, 128));

before = std::move(point);

}

void calculate_delaunay(cv::Subdiv2D& subdiv, const std::vector<cv::Point2f>& points, cv::Mat& canvas)

{

std::vector<cv::Vec4f> edge_list{};

subdiv.getEdgeList(edge_list);

for (const auto& edge: edge_list)

{

cv::Point p1{static_cast<int>(edge.val[0]), static_cast<int>(edge.val[1])};

cv::Point p2{static_cast<int>(edge.val[2]), static_cast<int>(edge.val[3])};

cv::line(canvas, p1, p2, cv::Scalar(48, 128, 48));

}

これを実行すると以下の画像を得る。母点の数は1000である。最初の画像がDelaunay三角形分割、後の画像が同じ母点に対するVoronoi図である。

まとめ

今回は、Delaunay三角形分割とVoronoi図を紹介した。両者とも応用例の大変多い構造である。

参考文献

計算幾何プログラミング：アルゴリズムをプログラムに置き換えるとき、数式には出てこないさまざな例外を考える必要がある。その例外を回避する手法を詳細に論じた良書である。

コンピュータ・ジオメトリ：広く読まれているこの分野の教科書。

Kumada Seiya

仕事であろうとなかろうと勉強し続ける、その結果”中身”を知ったエンジニアになれる

Google Vision APIでOCR

新卒・中途の皆さん！魔法のワード、「ベンチャー」に惑わされるな

Delaunay三角形分割

はじめに

Delaunay三角形分割とは

Delaunay三角形分割の定義

双対構造：Voronoi図

実装

まとめ

参考文献

Kumada Seiya

最近の記事

DBSCANを用いたクラスタリング

レイアウト解析

ChatGPTを用いた画像異常検知

クォータニオンと回転行列の関係

OpenCVで深層学習

OpenCVのcv::Matに関する小ネタ集

3D CADデータと誤差

STLファイルのエラーチェック

3D CADデータの形状認識

STLファイルの体積計算

回転

Google Vision APIでOCR

アーカイブ

カテゴリー