simple_delaunay_refinement.py

#!/usr/bin/env python
# coding:utf-8

from __future__ import print_function
import os
import math
import numpy as np
import scipy.spatial
import matplotlib.pyplot as plt
from random import random
from datetime import datetime

# 単位円内にランダムに置いた点でドロネー分割をする。
# その後、点の位置をいい感じに調整する。

# 点の数
NUM_POINTS = 300
# 初期値での点の間の最短距離
MIN_DIST_BETWEEN_POINTS = 0.01

# 調整の繰り返し最大数
MAX_ITERATION = 100
# 各点の調整後の位置の誤差(の２乗)の最大値
MAX_ERROR_SQUARED = 1e-5

def plot(points, tri):
    u"""
    Args:
        points (numpy.ndarray, shape=(x, 2)):
        tri (scipy.spatial.qhull.Delaunay):
    """
    fig = plt.figure(figsize=(8,8), dpi=80)
    ax = fig.add_subplot(1,1,1)
    lim = 1.2
    ax.axis([-lim, lim, -lim, lim])

    plt.triplot(points[:,0], points[:,1], tri.simplices.copy(),
            color='black', linewidth=1)
    plt.plot(points[:,0], points[:,1], 'o', color='black', markersize=4)

    plt.show()

def generatePoints(num_points, min_dist=0.01):
    u"""点の初期値を作成。(-1 < x < 1, -1 < y < 1)
    Args:
        num_points (int): 点の数
        min_dist (float): 点と点の間の最短距離
    Returns:
        numpy.ndarray: shape=(num_points, 2)
    """
    points = []
    wpi = math.pi * 2.0
    min_dist2 = min_dist**2
    while True:
        t = wpi * random()
        m = math.sqrt(random())
        p = np.array([math.cos(t), math.sin(t)]) * m
        ng = False
        for q in points:
            if np.sum((p - q)**2) < min_dist2:
                ng = True
                break
        if ng:
            continue
        points.append(p)
        if len(points) >= num_points:
            break
    return np.array(points)

def main():
    # 点の初期値を作成(単位円の内部に)
    points = generatePoints(NUM_POINTS, MIN_DIST_BETWEEN_POINTS)
    print("points generated")

    # 初期状態をプロット
    D = scipy.spatial.Delaunay(points)
    plot(D.points, D)

    time_start = datetime.now()

    for n in range(MAX_ITERATION):
        if n % 5 == 0:
            print(n, end=" ")

        n_points = []  # 調整後の点のリスト
        ng = False

        # 各点(i=index)について以下を行う
        for i in range(len(D.points)):
            # 点i
            point = D.points[i]
            # 点iを含む三角形を抽出
            idxs = D.simplices[(D.simplices == i).any(axis=1)]
            # indexの重複を取り除き、各indexの出現回数を取得
            idxs, counts = np.unique(idxs, return_counts=1)
            # 点iは以下の計算から除外(...しなくてもよい？)
            idxs = idxs[idxs != i]
            # indexに対応する座標を取得
            ps = D.points.take(idxs, axis=0)
            # 座標の平均値を取得。調整後の点iの位置とする
            avg = np.mean(ps, axis=0)
            
            # 点iが外周上の場合、円周上まで移動させる
            # (出現回数1の点があれば、点iは三角形で囲まれていない)
            if (counts == 1).any():
                avg = avg / np.sqrt(np.sum(avg**2))
            # 元の座標との距離(の２乗)が大きければ要再処理とする
            if np.sum((avg - point)**2) > MAX_ERROR_SQUARED:
                ng = True
            n_points.append(avg)
            
        points = np.array(n_points)
        D = scipy.spatial.Delaunay(points)
        
        if not ng:
            break
 
    print("\ntime=", datetime.now() - time_start)
    plot(D.points, D)

if __name__=="__main__":
    main()