Kei Minagawa's Blog

皆川圭(@keimina)のブログ、Pythonで試したことを書いていきます

ボロノイ図を書いて、マウスクリックした時に一番近いデータポイントをを表示するプログラムを実装する

ボロノイ図を書くプログラムと、マウスクリックした時に一番近いデータポイントをを表示するプログラムを実装します。

1 ボロノイ図とは

ボロノイ図の詳細は Wikipedia を参照してください。

ja.wikipedia.org

簡単な例をあげると、xy平面上に点Aと点Bがあったとして、点Xが点Aに近ければ赤、点Bに近ければ青というように、点Xに色を塗るとします。xy平面上にある点全ての点Xに前者の操作を繰り返していくと、全ての領域が赤または青に塗られると思います。この例では、点Aと点Bの二つしかありませんでしたが、多数の点でこのようなことを行うと、色のついた領域が四角形や六角形などの多角形になり、独特な図になります。このような図形をボロノイ図と呼ぶと認識しています。

2 ボロノイ図を書くプログラム実装する

2.1 コード

ボロノイ図を書くコードを実装します。コードは以下の通りです。

# -*-coding:utf-8 -*-
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

np.random.seed(0)

# グリッドに関する設定
GRID_N = 100
GRID_M= 100

STEP_x = 1.0/GRID_N
STEP_y = 1.0/GRID_M
mx, my = np.meshgrid(np.linspace(0, 1, GRID_N), np.linspace(0, 1, GRID_M))

# データに関する設定
DATA_N = 10
df = pd.DataFrame(np.random.random((DATA_N, 2)), columns=['x', 'y'])

# データをブロードキャスト(データを3階のテンソルへ変換)
x_broad = np.broadcast_to(df.x[:, np.newaxis, np.newaxis], (DATA_N, GRID_N, GRID_N))
y_broad = np.broadcast_to(df.y[:, np.newaxis, np.newaxis], (DATA_N, GRID_M, GRID_N))

# 2乗距離が最小の平面のインデックスを求める
dx = x_broad - mx
dy = y_broad - my

r = dx**2 + dy**2
z = r.argmin(axis=0)

# 描画
fig, ax = plt.subplots()
ax.plot(df.x, df.y, lw=0, marker='o')
ax.pcolormesh(mx, my, z)

# アスペクト比を 1:1 にする
ax.set_aspect('equal')
plt.show(True)

2.2 コードの説明

コードの説明は以下の通りです。

説明
21行目
x_broad = np.broadcast_to(df.x[:, np.newaxis, np.newaxis], (DATA_N, GRID_N, GRID_N))
np.broadcast_to 関数を使用して、配列の 1D -> 3D へのブロードキャストを行っています。np.broadcast_to 関数については、以下を参照してください。

keimina.hatenablog.jp

25
dx = x_broad - mx
np.meshgrid で作成した行列 x からデータのポイントをブロードキャストして作成した配列を引き算しています。 3次元と2次元の計算ですので、これもブロードキャストの演算が行われます。イメージ的には0軸が上からしたに向かって伸びており、そこから横に、枝のように、1軸と2軸からなる平面がぶら下がってる状態で、そのぶら下がってる平面それぞれに対して一つ、一つに行列 mx の平面と引き算を行うという感じです。何行ってるかわかんないですね。
28
r = dx**2 + dy**2
3次元の r 配列の要素である行列はデータポイントとの距離を表しています。
29
z = r.argmin(axis=0)
3次元の r 配列の0軸方向で一番値の小さいインデックスを求めています。これは、xy座標上のある領域において、その領域では、どのデータポイントが一番近いか求めているということです。
33
ax.plot(df.x, df.y, lw=0, marker='o')
plot 関数により、データポイントの描画を行います。 scatter でも良いでしょう。
34
ax.pcolormesh(mx, my, z)
pcolormesh 関数により、領域に色を付けます

3 出力結果

出力結果は以下になります。青い点がある領域内における点に一番近い青い点は、その領域内にある青い点になっていることがわかります。

f:id:keimina:20191215203434p:plain

図1. 出力結果(ボロノイ図)

4 マウスクリックしたときに一番近いデータポイントを表示するプログラムを実装する

Matplotlib で何らかのデータの点の分布を描画したときに、一番近いデータポイントを表示したい時があります。これは、ボロノイ図アルゴリズムを応用し簡単に実現できます。以下にコードだけ貼り付けておきます。

4.1 コード

# -*-coding:utf-8 -*-
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

np.random.seed(0)

# グリッドに関する設定
GRID_M = 100
GRID_N = 100

STEP_x = 1.0/GRID_N
STEP_y = 1.0/GRID_M
mx, my = np.meshgrid(np.linspace(0, 1, GRID_N), np.linspace(0, 1, GRID_M))

# データに関する設定
DATA_N = 10
df = pd.DataFrame(np.random.random((DATA_N, 2)), columns=['x', 'y'])

print(df)

# データをブロードキャスト(データを3階のテンソルへ変換)
x_broad = np.broadcast_to(df.x[:, np.newaxis, np.newaxis], (DATA_N, GRID_M, GRID_N))
y_broad = np.broadcast_to(df.y[:, np.newaxis, np.newaxis], (DATA_N, GRID_M, GRID_N))

# 2乗距離が最小の平面のインデックスを求める
dx = x_broad - mx
dy = y_broad - my

r = dx**2 + dy**2
z = r.argmin(axis=0)

# 描画
fig, ax = plt.subplots()
ax.plot(df.x, df.y, lw=0, marker='o')
ax.pcolormesh(mx, my, z)

def onclick(event):
        if (event.xdata is None) or (event.ydata is None):
                return
        gx = int(event.xdata*GRID_N)
        gy = int(event.ydata*GRID_M)
        index = int(z[gy, gx])
        print(df.loc[[index]])

cid = fig.canvas.mpl_connect('button_press_event', onclick)

# アスペクト比を 1:1 にする
ax.set_aspect('equal')
ax.set_xlim(0,1)
ax.set_ylim(0,1)
plt.show(True)

f:id:keimina:20191215203255g:plain

図2. マウスクリックした座標に一番近いデータポイントの座標を表示するプログラム

5 まとめ

ボロノイ図を書くプログラムを実装しました。またその応用としてマウスクリックした座標に一番近いデータポイントの座標を表示するプログラムを作りました。

以上です。

Author: MINAGAWA KEI

Created: 2019-12-15 日 20:24

Validate