Delaunay triangulation by Python
Pythonのmatplotlibを使ってドロネー三角分割をやってみた。
1 Introduction
ランダムな点の集合から適当に線を引っ張って格子(メッシュ)を作りたいことがある。これを実現する方法として，delaunay(ドロネー，デローニ)三角分割というアルゴリズムがある。簡単にできる方法がないか調べると，Pythonでこのアルゴリズムを使えるモジュールとして以下がみつかった。
- scipy
- matplotlib
- triangle
- pyhull
この中で最も手軽なのはmatplotlibだろう。scipyはインストールが若干面倒。matplotlibだと三角を作りながらグラフも書けるのでこれがよい。
matplotlibの中でもdelaunay三角分割を行えるモジュールには以下の二つがある。
- matplotlib.delaunay
- matplotlib.tri
matplotlib.triはmatplotlibで三角形を取り扱うモジュールとなっている。matplotlib.delaunayはmatplotlibの中でdelaunay三角分割を行うためのモジュールとなっている。
迷うところだが，基本的にはmatplotlib.triを使えばよい。理由は大きく以下2点である。
- matplotlib.triを読み込むとmatplotlib.delaunayも自動で読み込んでいること。
- site-packages/matplotlib.tri/triangulation.pyを開くと確認可能。
- matplotlib.triのほうが公式のマニュアルの説明が充実していること。
- 参考：triangular grids — Matplotlib 1.4.0 documentation http://matplotlib.org/api/tri_api.htm
以下のサイトでは，matplotlib.delaunayを使ったサンプルを示しており参考にした。
参考：ドロネー三角形 matplotlibで | BoxHeadRoom http://boxheadroom.com/2010/09/18/matplotlib_delaunay
2 Sample
実際にmatplotlib.triを使ってドロネー三角分割で適当な点の集合から三角形メッシュを作ってプロットした。以下がそのコードだ。
#!/usr/bin/env python2.7
# coding: utf-8
# (File name: delaunay-test.py)
# Author: SENOO, Ken
# (Last Update: 2014-09-20T19:31+09:00)
# License: MIT
import matplotlib.pyplot as plt
import matplotlib.tri as tri
import numpy as np
x=np.random.randn(10)
y=np.random.randn(10)
triangle=tri.Triangulation(x,y)
plt.clf()
## plot triangle mesh
# plt.triplot(x,y,triangle.triangles, “ro-“)
plt.triplot(triangle, “ro-“)
## plot node label
for i in range(len(x)):
plt.text(x[i], y[i], i)
plt.text(x[i], y[i], i)
## center of triangle
cx=[]; cy=[]
for tria in triangle.triangles:
cx.append(np.mean([x[i] for i in tria], dtype=np.float32))
cy.append(np.mean([y[i] for i in tria], dtype=np.float32))
## plot triangle center label
plt.plot(cx, cy, “^b”)
for tria in range(len(triangle.triangles)):
plt.text(cx[tria], cy[tria], tria)
plt.savefig(“tri.png”, bbox_inches=”tight”)
これを実行すると例えば以下のような図ができる。np.random.randnによりランダムに点を発生させているので同じ図はできない。
3 Description
簡単に上記のコードの要点を説明する。
最初にモジュールを読み込み，ランダムな点を生成している。以下のコードでドロネー三角分割を行っている。
triangle=tri.Triangulation(x,y)
これにより，triangle変数にmatplotlib.tri.triangulation.Triangulationオブジェクトを格納している。このオブジェクトに三角形の座標や位相の情報が入っている。主な属性は以下の表にまとめた。
表 6.3: Triangulationオブジェクトの主な属性
属性
|
説明
|
x, y
|
ノードの座標値。
|
edge
|
点同士をどう結ぶかの情報が入っている。[[2, 3], …, [3,5]]のような(N, 2)の配列。
|
triangles
|
三角形の構成ノード番号(反時計回り)の配列。(N, 3)の配列。
|
neighbors
|
ノード番号を共有している三角形番号の配列。(N, 3)の配列。共有していなければ-1が入る。
|
この他にもmaskなどある。情報は上記の表の属性にアクセスできればどうにかなるだろう。
続いて，以下のコードでプロットしている。
## plot triangle mesh
# plt.triplot(x,y,triangle.triangles, “ro-“)
plt.triplot(triangle, “ro-“)
plt.triplot()という関数は三角形のメッシュをプロットする関数となっている。Triangulationオブジェクトを想定して作られており，このオブジェクトを渡せば簡単にプロットしてくれる。その他にも，三角形の構成ノードリストを渡せばプロットしてくれる。
最後に，ノードと三角形の対応を見るためにラベルをプロットしている。三角形の中心座標を取得するために三角形の3点の平均をとって新しく配列を作りなおしている。
## center of triangle
cx=[]; cy=[]
for tria in triangle.triangles:
cx.append(np.mean([x[i] for i in tria], dtype=np.float32))
cy.append(np.mean([y[i] for i in tria], dtype=np.float32))
もう少し配列操作を使った一括演算によるうまい方法があるような気もする。
4 Conclusion
matplotlib.triを使ってドロネー三角分割を行った。中心座標を算出してノードと三角形ごとに番号のラベルを振り，グラフにプロットした。今後何かで三角形メッシュを作る必要が合った時に役に立つだろう。
