前回の続きです。今回は、前回に引き続き、ポリゴンに収まる点の集合を作成します。アルゴリズムは簡単で、ランダムに点を打ち、それがポリゴンの内側であれば点を残す、ポリゴンの外側であれば点を消すを繰り返すだけです。「ポリゴンの内側か？」を判定する関数を作成し、その関数に点の位置を次々と問い合わせればよいでしょう。
※ここでは「ポリゴン」は多角形の意味で使っています。
前回の記事:
keimina.hatenablog.jp
前々回の記事：
keimina.hatenablog.jp
N角形(N>=3)の領域内の点の集合を得ることができますので、前回や前々回の記事で作成した点の集合も、この記事の方法だけで作成することができるようになります。
目次
1. 点がポリゴンの領域の内側にある点か判定する方法
点がポリゴンの領域の内側にあるか判定する方法は以下に色々書かれていました。
上記の回答の一つを見たところ、点がポリゴンの内側にあるとき、その点とポリゴンを構成する点から算出される角度の総和は必ず360度になるらしいです。図を書いてイメージをつかみます。以下が、点A,点B,点C,,,点Jを結んだポリゴンの内部に点Pがあった時の図になります。点Bから点Cあたりでガタガタと移動する箇所があります。点Bから点Cに移動する時のように角度θiがマイナスになるときもあるのですが、その後のガタガタで相殺されます。θiの総和が360°のとき左回りに１回転できたということになります。※かなりざっくりした説明で厳密ではありません。
一方、点がポリゴンの外側にあるとき、360度にはならないと書かれています。これも図を書いてイメージをつかみます。θiの箇所についてみて見ると前半は右回りで、後半は左回りとなっています。右回りと、左回りで相殺されますので、θiの総和は0°になることがわかります。
ポリゴンのような多角形の図形は、円の一つだと考えられます。そのため多角形を円とみなして考えるとわかりやすいです。円形の壁の中に、自分がいたとして、ぐるっと見回して壁に囲まれていたら、円の中だとわかります。一方、ぐるっと見回して左には壁があったけど右には壁がなかったなどの場合は壁の中でないことがわかります。算出される角度の総和が360°だったというのは、壁をたどると左周りに１回転できたということと同じです。１回転できなかった場合は、壁を辿っても１回転できなかったということですので、壁がつながっていない、すなわち壁の内側ではないということがわかります。
今回は本アルゴリズムを使用します。(アルゴリズムがちゃんと機能するのかも確認します)
2. 角度の求め方
角度を求めるには以下の数学公式を用います。
cosθ = (PA・PB)/(|PA|*|PB|)
θ = arccos(cosθ)
3. 角度を求めるための関数を作成する
numpyで角度を求める関数を実装すると以下になります。
コード：
import numpy as np from numpy import inner from numpy.linalg import norm def angle(P, A, B): ''' args ---- P: np.array(x, y)) A: np.array(x, y)) B: np.array(x, y)) returns ------- theta: x (degree) ''' # 全体をシフトさせる PA = A - P PB = B - P # cosθを求める costheta = inner(PA, PB)/norm(PA)/norm(PB) # 角度を求める(度) theta = np.arccos(costheta)/np.pi * 180 return theta # テスト P = np.array((0,0)) A = np.array((1,1)) B = np.array((1,0)) print(angle(P, A, B))
出力：
45.00000000000001
4. 時計回りか、反時計回りかを判定する
3. によって角度の絶対値はもとめることができました。しかし、反時計回りか時計回りか区別できません。時計回りか、反時計回りか区別できないと、角度の総和が常に360°になってしまうため、区別できるようにする必要があります。そこで、点A(a, b)と点B(c, d)と点P(0, 0)によってできるベクトルPAとPBのなす角θについて、時計回りと反時計回りを以下のように定義します。
(i) d > (b/a) * c
ならば、反時計回り方向の角度
(ii) d < (b/a) * c
ならば時計回り方向の角度
※説明を簡単にするため点P(0, 0)としています。点Pが(0, 0)以外のときでも、全体をシフトさせて点Pが(0, 0)になるように点Aと点Bを移動させても、角度に関する情報は失われないので、P(0, 0)でも問題ありません。
また(i)(ii)の不等号の式の両辺にaをかけて、数式を変形すると以下のように書き直せます。
(i)' a*d - b*c > 0
(ii)' a*d - b*c < 0
図で表すと以下になります。点Bが緑色の領域にあればPABの角度は反時計回り方向。緑色の領域にあれば、時計回り方向というふうに定義しています。
なぜ、この定義をするのかは、説明しづらいですが、図1、図2を実際に手でなぞって見ると、まず回転方向が必要なことがわかります。また0°<θ<180°のときは足し算で-180°<θ<0°のときは引き算のように計算方法を変えなければいけないこともわかると思います。そのため、このような定義が必要であるということがわかります。
5. 角度を求めるための関数2(反時計回り、時計回りの情報を追加)
3.で作成したangle関数を、4.の定義を使用して、反時計周りならプラスを、時計回りならマイナスの角度を返すように修正します。
コード：
import numpy as np from numpy import inner from numpy.linalg import norm def angle2(P, A, B): ''' args ---- P: np.array(x, y)) A: np.array(x, y)) B: np.array(x, y)) returns ------- theta: x (degree) ''' # 全体をシフトさせる PA = A - P PB = B - P # cosθを求める costheta = inner(PA, PB)/norm(PA)/norm(PB) # 角度を求める(ラジアン) theta = np.arccos(costheta) # 角度を求める(度) theta = np.arccos(costheta)/np.pi * 180 # if PA[0]*PB[1] - PA[1]*PB[0] > 0: # 条件(i)'に該当するときはなにもしない(ここではthetaはプラス) pass else: # 条件(ii)'に該当するときはthetaをマイナスにする theta *= -1.0 return theta # テスト P = np.array((0,0)) A = np.array((1,1)) B = np.array((1,0)) print(angle2(P, A, B))
出力：
-45.00000000000001
6. 点がポリゴンの領域の内側にある点かどうかを判定する関数を作成しテストする
準備が整いましたので、点がポリゴンの領域の内側にある点かどうかを判定する関数を作成し、機能するかテストします。
点(x, y)が1.のアルゴリズムを使用しポリゴンの領域内の点であるかどうかを判定するis_inside_polygon関数を実装します。以下のようになります。
この関数は[(x1, y1), (x2, y2),,,]のような点(x, y)の集合pointsと点p(x, y)を受け取り、pointsからなるポリゴンの領域内に点pが含まれるか判定します。
コード：
import random import matplotlib.pyplot as plt import numpy as np from numpy import inner from numpy.linalg import norm def angle2(P, A, B): ''' args ---- P: np.array((x, y)) A: np.array((x, y)) B: np.array((x, y)) returns ------- theta: float (dgree) ''' # 全体をシフトさせる PA = A - P PB = B - P # cosθを求める costheta = inner(PA, PB)/norm(PA)/norm(PB) # 角度を求める(度) theta = np.arccos(costheta)/np.pi * 180 # 時計回りか、反時計回りか判定し、thetaをプラスかマイナスにして返す if PA[0]*PB[1] - PA[1]*PB[0] > 0: # 条件(i)'に該当するときはなにもしない(ここではthetaはプラス) pass else: # 条件(ii)'に該当するときはthetaをマイナスにする theta *= -1.0 return theta def is_inside_polygon(points, P): ''' args ---- points: list of np.array((x, y)) P: np.array((x, y)) returns ------- theta: x (degree) ''' points = [np.array(p) for p in points] P = np.array(P) angle_total = 0.0 for A, B in zip(points, points[1:] + points[0:1]): # Notice! "points[1:] + points[0:1]" means something like "points[1:].extend(points[0:1])" angle_total += angle2(P, A, B) if np.isclose(abs(angle_total), 360): # 360°の場合はポリゴンの内側なのでTrueを返す return True else: # ポリゴンの外側 return False # 三角形 #points = ((0,0), (1,0), (0, 1)) # リボン? #points = ((1.0, 0.0), (0.5, 0.5), (1.0, 1.0), (-1.0, 1.0), (-0.5, 0.5), (-1.0, 0.0)) # 十字 points = ((1.0, -3.0), (1.0, -1.0), (3.0, -1.0), (3.0, 1.0), (1.0, 1.0), (1.0, 3.0), (-1.0, 3.0), (-1.0, 1.0), (-3.0, 1.0), (-3.0, -1.0), (-1.0, -1.0), (-1.0, -3.0)) points = [(i[0]/3, i[1]/3) for i in points] N = 5000 for _ in range(N): # 点を打つ場所をランダムに決める(ただし、-1<x<1, -1<y<1とする) x = (random.random() - 0.5) * 2 y = (random.random() - 0.5) * 2 p = (x, y) # 決めた点(x, y)が三角形の内側にあればグラフに点を打つ if is_inside_polygon(points, p): plt.plot(x, y, color='k', marker='.') # グラフを表示する plt.gca().set_xlim(-2,2) plt.gca().set_ylim(-2,2) plt.gca().set_aspect('equal') plt.show()
7. 出力結果
上から順に三角形、リボン、十字形になります。ランダムに点を打っているのですが、pointsで定義された多角形の内側にある点だけ残ることがわかります。
※形を変えるには、多角形の定義(#points... の箇所)を変えてください。
8. まとめ
ポリゴンに収まる点の集合を作成できました。
以上です。
