ドロネー分割からのボロノイ分割
はじめに
ちょいとボロノイとドロネーに関して調べる機会があったので記事にしてみます。
以下結果です。
因みに、実行はRhino。コードはPythonです。
ですが、点群を取得、線の描画、以外はRhinoと依存関係を保有していませんので他の実行環境に応用することが可能だと思いますので、是非参考にしてみてください。
そもそもドロネー分割やボロノイ分割って？
ドロネーとボロノイをすごく雑にすごく手短に説明します。詳しくは先人たちが解説してくださっているはずなので検索してみてください。
イメージとしては、
任意の点集合を構成する点を母点とします。
・ドロネー分割は母点の隣接関係
・ボロノイ分割は母点の勢力範囲
と解釈できます。
ドロネー分割
今回使用したアルゴリズムの流れは以下のようになります。
- 点集合を内包する三角形を作成。
- 三角形の外接円内にある任意の点を選択。
- 選択した点を用い三角形を再分割。もとの三角形は削除。
- 2.で追加した任意の点を外接円で内包しているすべての三角形に対して3.の操作を行う。
- 3.と4.の操作で新たに作成された三角形の内、重複している三角形を削除。
- 2~5の流れを追加する点がなくなるまで繰り返す。
多分、イメージつきにくいと思いますので実際にコードと見比べてみてください。
以下、至らない点も多いとは思いますがソースコードです。
ファイルは実行ファイルのDelaunay.py モジュールのTriangle.py, Point.pyの３つで構成されています。
Delaunay.py
# coding: utf-8
import Rhino
import scriptcontext
import rhinoscriptsyntax as rs
import Triangle
import Point
points_obj_list = [] # Pointのオブジェクトを保持。
divide_triangle_list = [] # このリストに三角形分割を保持していく。[p1, p2, p3]
# Rhino上の点群を取得し座標値を取得。
points = rs.GetObjects('select points to use Delaunay', rs.filter.point)
for i in range(len(points)):
point = Point.Point(points[i])
point.system_guid_obj_to_coordinate()
points_obj_list.append(point)
# ドロネー分割の前準備。
# はじめの外接円を構成する点のインスタンス作成。この範囲内に点群を収めておく必要あり。
P1 = Point.Point([-100000, -100000, 0])
P1.change_coordinate()
P2 = Point.Point([100000, -100000, 0])
P2.change_coordinate()
P3 = Point.Point([0, 73000, 0])
P3.change_coordinate()
triangle = Triangle.Triangle(P1, P2, P3)
triangle.cul_center_coordinate_and_radius()
divide_triangle_list.append(triangle)
# ドロネー分割のメインアルゴリズム
for _ in range(len(points_obj_list)):
# 追加するポイントの選択。
select_point = points_obj_list.pop(0)
# 三角形の外接円内に追加したポイントが内包されていないか判定。内包されている場合はindexを保存。
temp_divide_triangle = []
count = 0
for i in range(len(divide_triangle_list)):
tri = divide_triangle_list[i + count]
check = tri.check_point_include_circumscribed_circle(select_point)
if check:
# 内包していた三角形の頂点を使用して再分割する。
new_triangle1 = Triangle.Triangle(tri.p1, tri.p2, select_point)
new_triangle2 = Triangle.Triangle(tri.p2, tri.p3, select_point)
new_triangle3 = Triangle.Triangle(tri.p1, tri.p3, select_point)
new_triangle1.cul_center_coordinate_and_radius()
new_triangle2.cul_center_coordinate_and_radius()
new_triangle3.cul_center_coordinate_and_radius()
temp_divide_triangle.append(new_triangle1)
temp_divide_triangle.append(new_triangle2)
temp_divide_triangle.append(new_triangle3)
del divide_triangle_list[i + count]
count = count - 1
else:
pass
# 重複している三角形を削除
_count = 0
for i in range(len(temp_divide_triangle)):
triangle_check1 = temp_divide_triangle[i + _count]
flag_del = False
count = 0
for j in range(len(temp_divide_triangle)):
if j <= i + _count:
continue
if len(temp_divide_triangle) - 1 < j + count:
break
triangle_check2 = temp_divide_triangle[j + count]
if triangle_check1.center_p == triangle_check2.center_p:
del temp_divide_triangle[j + count]
count = count - 1
flag_del = True
if flag_del:
del temp_divide_triangle[i + _count]
_count = _count - 1
if len(temp_divide_triangle) - 2 < i + _count:
break
# 重複三角形を削除した後にappend
divide_triangle_list.extend(temp_divide_triangle)
# ドロネー描画
for i in range(len(divide_triangle_list)):
divide_triangle_list[i].draw_divide_triangle()
# ボロノイ分割。
for i in range(len(divide_triangle_list)):
tri_1 = divide_triangle_list[i]
for j in range(len(divide_triangle_list)):
if i == j:
continue
tri_2 = divide_triangle_list[j]
count = 0
if tri_1.p1 == tri_2.p1 or tri_1.p1 == tri_2.p2 or tri_1.p1 == tri_2.p3:
count += 1
if tri_1.p2 == tri_2.p1 or tri_1.p2 == tri_2.p2 or tri_1.p2 == tri_2.p3:
count += 1
if tri_1.p3 == tri_2.p1 or tri_1.p3 == tri_2.p2 or tri_1.p3 == tri_2.p3:
count += 1
if count == 2:
line = Rhino.Geometry.Line(tri_1.center_p[0], tri_1.center_p[1], tri_1.center_p[2], tri_2.center_p[0],
tri_2.center_p[1], tri_2.center_p[2])
scriptcontext.doc.Objects.AddLine(line)
Point.py
# coding: utf-8
import rhinoscriptsyntax as rs
class Point:
def __init__(self, point):
self.point = point
self.x = None
self.y = None
self.z = None
self.coordinate = None
def system_guid_obj_to_coordinate(self):
'''System.GuidObjectをRhino.Objectに変換。クラス変数に座標値を保持させる'''
p = rs.coerce3dpoint(self.point)
self.x = p[0]
self.y = p[1]
self.z = p[2]
self.coordinate = [p[0], p[1], p[2]]
def change_coordinate(self):
'''はじめの外接円を構成する点を変換するメソッド'''
p = self.point
self.x = p[0]
self.y = p[1]
self.z = p[2]
self.coordinate = [p[0], p[1], p[2]]
Triangle.py
# coding: utf-8
import math
import Rhino
import scriptcontext
class Triangle:
def __init__(self, p1, p2, p3):
self.p1 = p1
self.p2 = p2
self.p3 = p3
self.radius = None
self.center_p = None
def cul_center_coordinate_and_radius(self):
'''任意の3点の外接円の中心と半径を取得。'''
# 外接円の中心点を計算
c = 2 * ((self.p2.x - self.p1.x)*(self.p3.y - self.p1.y) - (self.p2.y - self.p1.y)*(self.p3.x - self.p1.x))
x = ((self.p3.y - self.p1.y)*(pow(self.p2.x, 2) - pow(self.p1.x, 2) + pow(self.p2.y, 2) - pow(self.p1.y, 2)) +
(self.p1.y - self.p2.y)*(pow(self.p3.x, 2) - pow(self.p1.x, 2) + pow(self.p3.y, 2) - pow(self.p1.y, 2))) / c
y = ((self.p1.x - self.p3.x)*(pow(self.p2.x, 2) - pow(self.p1.x, 2) + pow(self.p2.y, 2) - pow(self.p1.y, 2)) +
(self.p2.x - self.p1.x)*(pow(self.p3.x, 2) - pow(self.p1.x, 2) + pow(self.p3.y, 2) - pow(self.p1.y, 2))) / c
self.center_p = [x, y, 0]
# 外接円の半径を計算
radius = math.sqrt(pow((self.p1.x - self.center_p[0]), 2) + pow((self.p1.y - self.center_p[1]), 2))
self.radius = radius
def draw_divide_triangle(self):
'''分割三角形をRhinoに描画するためのメソッド'''
line1 = Rhino.Geometry.Line(self.p1.x, self.p1.y, self.p1.z, self.p2.x, self.p2.y, self.p2.z)
line2 = Rhino.Geometry.Line(self.p1.x, self.p1.y, self.p1.z, self.p3.x, self.p3.y, self.p3.z)
line3 = Rhino.Geometry.Line(self.p2.x, self.p2.y, self.p2.z, self.p3.x, self.p3.y, self.p3.z)
scriptcontext.doc.Objects.AddLine(line1)
scriptcontext.doc.Objects.AddLine(line2)
scriptcontext.doc.Objects.AddLine(line3)
def check_point_include_circumscribed_circle(self, check_p):
'''外接円に任意の点が内包されているかどうか判定'''
distance = math.sqrt(pow((check_p.x - self.center_p[0]), 2) + pow((check_p.y - self.center_p[1]), 2))
if distance < self.radius:
return True
else:
return False
ボロノイ分割
ソースコードを読んでいただいた方はお気づきかもしれませんが、Delaunay.pyの末尾にボロノイ分割を実行するコードが紛れています。
ドロネー分割実行後であれば、Triangleインスタンスの情報を用いいて数行のコード追加で実行できます。
以下、抜き出したものです。
# ボロノイ分割。
for i in range(len(divide_triangle_list)):
tri_1 = divide_triangle_list[i]
for j in range(len(divide_triangle_list)):
if i == j:
continue
tri_2 = divide_triangle_list[j]
count = 0
if tri_1.p1 == tri_2.p1 or tri_1.p1 == tri_2.p2 or tri_1.p1 == tri_2.p3:
count += 1
if tri_1.p2 == tri_2.p1 or tri_1.p2 == tri_2.p2 or tri_1.p2 == tri_2.p3:
count += 1
if tri_1.p3 == tri_2.p1 or tri_1.p3 == tri_2.p2 or tri_1.p3 == tri_2.p3:
count += 1
if count == 2:
line = Rhino.Geometry.Line(tri_1.center_p[0], tri_1.center_p[1], tri_1.center_p[2], tri_2.center_p[0],
tri_2.center_p[1], tri_2.center_p[2])
scriptcontext.doc.Objects.AddLine(line)
ボロノイ分割のアルゴリズムは
・ドロネー分割で作成された共有する辺を持つ２つの三角形の外接円の中心点を直線で結ぶ。
以上です。（笑）
おわりに
単体で語られることが多い気がするボロノイとドロネーですが、コードで示されているように密接な関係があります。
隣接関係と勢力範囲を示す両者を組み合わせれば、様々なことに応用可能な気がします（気がするだけです笑）
最後までありがとうございました。