2015-01-02

[Triangulation」湖に島があってその中に池がある

python 幾何

昨日は、輪郭(Ａの字の外周〜湖があって、その中に「なかぬけ」として島(Ａの字のぐるり)があるときに、Triangulationを使ってプロットした。
それを扱うためのオブジェクトとしてdict(dictionaryタイプ)というのを使っているので、それの使い方に慣れつつ、トポロジー的に一段階複雑な、「湖に島があって、その中に池」というのを作ってみる

import numpy as np # ベクトル・行列・高次元アレイを使う
import scipy as sp # 線形代数・関数・数値計算
# matplotlibのpylabインタフェースは、MATLABの利用経験があるユーザがmatplotlibを簡単に習得できるように設計されている。というわけで以下２つを入れる
import matplotlib as mp 
import pylab as pl
import os # コンピュータのファイルの出し入れのため
os.chdir('C:\\Users\\ryamada\\Desktop\\python') # Windowsでは'\'が'\\'になることに注意
#from numpy import array

R = 1
n = 100
#t =  np.linspace(0,1,n) * 2 * np.pi
t = np.arange(n) * 2 * np.pi/n

x = R * sp.cos(t)
y = R * sp.sin(t)
mp.pyplot.plot(x,y)
mp.pyplot.xlim(x.min()-0.01, x.max()+0.01) # x,y軸の範囲を存在点座標の範囲より0.01広げる
mp.pyplot.ylim(y.min()-0.01, y.max()+0.01)
mp.pyplot.axis('off')
mp.pyplot.axes().set_aspect('equal') # x,y軸の長さを揃える
mp.pyplot.plot(x,y, 'b.')

R2 = 0.3
ctr_2x = 0.2
ctr_2y = 0.1
x2 = R2 * sp.cos(t) + ctr_2x
y2 = R2 * sp.sin(t) + ctr_2y
R3 = 0.1
ctr_3x = 0.22
ctr_3y = 0.11
x3 = R3 * sp.cos(t) + ctr_3x
y3 = R3 * sp.sin(t) + ctr_3y

output = {'vertices': None, 'holes': None, 'segments': None}

vertices = []
for i in range(n) :
	vertices.append([float(x[i]), float(y[i])])
for i in range(n) :
	vertices.append([float(x2[i]), float(y2[i])])
for i in range(n) :
	vertices.append([float(x3[i]), float(y3[i])])
output['vertices']=np.array(vertices)

segments = []
for i in range(n) :
	tmp = i-1
	if(i==0) :
		tmp = n-1
	segments.append([tmp,i])
for i in range(n) :
	tmp = i-1
	if(i==0):
		tmp = n-1
	segments.append([n+tmp, n+i])
for i in range(n) :
	tmp = i-1
	if(i==0):
		tmp = n-1
	segments.append([2*n+tmp, 2*n+i])
output['segments']=np.array(segments)

holes = []
N_holes = 1
for k in range(N_holes):
	holes.append([ctr_3x+R3*1.5, ctr_3y])
output['holes'] = np.array(holes)


from scipy.spatial import ConvexHull
import matplotlib.pyplot as plt
%matplotlib inline
lake_superior = output
vertices_ls = lake_superior['vertices']
vertices_ls.shape

%time hull = ConvexHull(vertices_ls)

plt.figure(figsize=(18, 18)) # 図のサイズ指定
plt.xlim(vertices_ls[:,0].min()-0.01, vertices_ls[:,0].max()+0.01) # x,y軸の範囲を存在点座標の範囲より0.01広げる
plt.ylim(vertices_ls[:,1].min()-0.01, vertices_ls[:,1].max()+0.01)
plt.axis('off')
plt.axes().set_aspect('equal') # x,y軸の長さを揃える
plt.plot(vertices_ls[:,0], vertices_ls[:,1], 'b.')
for simplex in hull.simplices: # hull.simplicesの要素(辺)を一つずつ描く
    plt.plot(vertices_ls[simplex, 0], vertices_ls[simplex, 1], 'r-')
    
plt.show()


from scipy.spatial import Voronoi, voronoi_plot_2d
vor = Voronoi(vertices_ls)

plt.figure(figsize=(10, 10))
ax = plt.subplot(111, aspect='equal')
voronoi_plot_2d(vor, ax=ax)
plt.xlim(-3,3) # 範囲は、適当に描かれた絵を見たあとで調整使用
plt.ylim(-3,3)
plt.show()

from scipy.spatial import Delaunay, delaunay_plot_2d

tri = Delaunay(vertices_ls)

plt.figure(figsize=(18, 18))
ax = plt.subplot(111, aspect='equal')
delaunay_plot_2d(tri, ax=ax)
plt.show()

from triangle import triangulate, plot as tplot
cndt = triangulate(lake_superior, 'p')

plt.figure(figsize=(18, 18))
ax = plt.subplot(111, aspect='equal')
tplot.plot(ax, **cndt)
plt.show()


cncfq20adt = triangulate(lake_superior, 'pq20a.001D')

plt.figure(figsize=(18, 18))
ax = plt.subplot(111, aspect='equal')
tplot.plot(ax, **cncfq20adt)
plt.show()

from scipy.spatial import minkowski_distance
X = cncfq20adt['triangles'][:,0]
Y = cncfq20adt['triangles'][:,1]
Z = cncfq20adt['triangles'][:,2]

Xvert = [cncfq20adt['vertices'][x] for x in X]
Yvert = [cncfq20adt['vertices'][y] for y in Y]
Zvert = [cncfq20adt['vertices'][z] for z in Z]

lengthsXY = minkowski_distance(Xvert, Yvert)
lengthsXZ = minkowski_distance(Xvert, Zvert)
lengthsYZ = minkowski_distance(Yvert, Zvert)

from scipy.sparse import lil_matrix
from scipy.sparse.csgraph import shortest_path
nvert = len(cncfq20adt['vertices'])
G = lil_matrix((nvert, nvert))
for k in range(len(X)):
    G[X[k], Y[k]] = G[Y[k], X[k]] = lengthsXY[k]
    G[X[k], Z[k]] = G[Z[k], X[k]] = lengthsXZ[k]
    G[Y[k], Z[k]] = G[Z[k], Y[k]] = lengthsYZ[k]
dist_matrix, pred = shortest_path(G, return_predecessors=True, directed=True, unweighted=False)
index = 2 # ウェブサイトで使っている指定番地では数が大きすぎるので適当な数を入れる
path = [2] # 同上

while index != 45: # ウェブサイトで使っている指定番地では数が大きすぎるので適当な数を入れる
    index = pred[45, index] # 同上
    path.append(index)
    
Xs = [cncfq20adt['vertices'][x][0] for x in path]
Ys = [cncfq20adt['vertices'][x][1] for x in path]
plt.figure(figsize=(18,18))
ax = plt.subplot(111, aspect='equal')
tplot.plot(ax, **cncfq20adt)
ax.plot(Xs, Ys, '-', linewidth=5, color='blue'); \
plt.show()

http://www.cs.cmu.edu/~quake/triangle.refine.html
http://docs.scipy.org/doc/scipy-0.14.0/reference/sparse.csgraph.html