平面多边形Wachspress坐标、中值坐标的计算及其等高线
1. Wachspress坐标等高线绘制
1.1 程序
# 计算有向面积
def dir_acr(point1,point2,point3):
result1 = np.linalg.det([[1, 1, 1], point1[0],point1[1]])
result2 = np.linalg.det([[1, 1, 1], point2[0],point2[1]])
result3 = np.linalg.det([[1, 1, 1], point3[0],point3[1]])
return result1,result2,result3
# 射线法判断点是否在给定区域内
def is_in_poly(p, poly):
"""
:param p: [x, y]
:param poly: [[], [], [], [], ...]
:return:
"""
px, py = p
is_in = False
for i, corner in enumerate(poly):
next_i = i + 1 if i + 1 < len(poly) else 0
x1, y1 = corner
x2, y2 = poly[next_i]
if (x1 == px and y1 == py) or (x2 == px and y2 == py): # if point is on vertex
is_in = True
break
if min(y1, y2) < py <= max(y1, y2): # find horizontal edges of polygon
x = x1 + (py - y1) * (x2 - x1) / (y2 - y1)
if x == px: # if point is on edge
is_in = True
break
elif x > px: # if point is on left-side of line
is_in = not is_in
return is_in
points = [[1.0,3.0,4.0,0.0], [0.0,0.0,3.0,2.0]]
n = len(points[0])
x_min = min(points[0])
x_max = max(points[0])
y_min = min(points[1])
y_max = max(points[1])
x = np.arange(x_min,x_max,0.01)
y = np.arange(y_min,y_max,0.01)
# X,Y = np.meshgrid(x,y)
Z=[]
# 记v为内点
for k in y:
for j in x:
v = [j,k]
w = []
if not is_in_poly(v, np.array(points).T.tolist()): # 判断点是否为内点
Z.append(np.nan)
continue
for i in range(n):
if i==0:
p_x = [points[0][-1]]+points[0][0:2]
p_y = [points[1][-1]]+points[1][0:2]
elif i==3:
p_x = points[0][i-1:i+1]+[points[0][0]]
p_y = points[1][i-1:i+1]+[points[1][0]]
else:
p_x = points[0][i-1:i+2]
p_y = points[1][i-1:i+2]
p = [p_x]+[p_y]
px_v1 = p_x[0:2]+[v[0]]
py_v1 = p_y[0:2]+[v[1]]
p_v1 = [px_v1]+[py_v1]
px_v2 = p_x[1:3]+[v[0]]
py_v2 = p_y[1:3]+[v[1]]
p_v2 = [px_v2]+[py_v2]
a1,a2,a3 = dir_acr(p,p_v1,p_v2)
wi = a1/(a2*a3)
w.append(wi)
if i==3:
w2 = wi
lambda2 = w2/sum(w)
Z.append(lambda2)
X,Y = np.meshgrid(x,y)
# Z =
#设置打开画布大小,长10,宽6
plt.figure(figsize=(10,6))
#填充颜色,f即filled
plt.plot(points[0]+[1],points[1]+[0])
# plt.contourf(X,Y,Z)
#画等高线
z = np.mat(Z)
z = np.array(z)
z.shape = (300,400)
contour = plt.contour(X,Y,z,[0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9])
# plt.clabel(contour,fontsize=10,colors=('k'))
plt.show()
2.中值坐标等高线绘制
2.1 步骤及原理
- 知识储备:numpy的三角函数,np.cos(),np.sin(),np.tan(),反三角函数,np.arccos(),np.arcsin(),np.arctan()
-
λ
i
=
ω
i
∑
j
=
1
k
ω
j
,
ω
j
=
tan
α
i
−
1
2
+
tan
α
i
2
∣
∣
v
i
−
v
∣
∣
\lambda_i = \frac{\omega_i}{\sum_{j=1}^k\omega_j},\omega_j=\frac{\tan{\frac{\alpha_{i-1}}{2}}+\tan{\frac{\alpha_{i}}{2}}}{||v_i-v||}
λi=∑j=1kωjωi,ωj=∣∣vi−v∣∣tan2αi−1+tan2αi
2.2 程序
def distance(point1,point2): # 计算点与点之间的距离
# 如果传入的是列表,则需要转化为np.array()数组,因为列表无法进行加减法
d_xy = np.array(point1)-np.array(point2)
d_xy = np.sum(np.power(d_xy,2))
d = np.sqrt(d_xy)
return d
def alpha(a,b,c): # 角度alpha的计算
cos_alpha = (np.power(b,2)+np.power(c,2)-np.power(a,2))/(2*b*c)
alpha = np.arccos(cos_alpha)
return alpha
# 给定一个多边形顶点points,记长度为n,随机选一个内点v
# 首先计算内点v到p_{i-1},pi,p_{i+1}的长度,然后再计算 p_{i-1}pi,pip_{i+1}的长度a1,a2,另外三条线记为l1,l2,l3
# 根据求得的长度计算角度,然后求正切值
points = [[1.0,3.0,4.0,0.0], [0.0,0.0,3.0,2.0]]
n = len(points[0])
x_min = min(points[0])
x_max = max(points[0])
y_min = min(points[1])
y_max = max(points[1])
x = np.arange(x_min,x_max,0.01)
y = np.arange(y_min,y_max,0.01)
# X,Y = np.meshgrid(x,y)
Z_2=[]
for k in y:
for j in x:
v = [j,k] # v为内点
w = []
if not is_in_poly(v, np.array(points).T.tolist()): # 判断点是否为内点
Z_2.append(np.nan)
continue
for i in range(n):
if i==0:
p_x = [points[0][-1]]+points[0][0:2]
p_y = [points[1][-1]]+points[1][0:2]
elif i==3:
p_x = points[0][i-1:i+1]+[points[0][0]]
p_y = points[1][i-1:i+1]+[points[1][0]]
else:
p_x = points[0][i-1:i+2]
p_y = points[1][i-1:i+2]
p = [p_x]+[p_y]
# print(p)
l1_point = [p_x[0],p_y[0]]
l2_point = [p_x[1],p_y[1]]
l3_point = [p_x[2],p_y[2]]
l1 = distance(l1_point,v)
l2 = distance(l2_point,v)
l3 = distance(l3_point,v)
a1 = distance(l1_point,l2_point)
a2= distance(l3_point,l2_point)
alpha1 = alpha(a1,l1,l2)
alpha2 = alpha(a2,l2,l3)
wi = (np.tan(alpha1/2)+np.tan(alpha2/2))/l2
w.append(wi)
if i == 3: # 修改需要进行计算的坐标顶点
w2 = wi
lambda_ = w2/sum(w)
Z_2.append(lambda_)
# 第1个坐标顶点的等高线
X,Y = np.meshgrid(x,y)
#设置打开画布大小,长10,宽6
plt.figure(figsize=(10,6))
#填充颜色,f即filled
plt.plot(points[0]+[1],points[1]+[0])
# plt.contourf(X,Y,Z)
#画等高线
z2 = np.mat(Z_2)
z2 = np.array(z2)
z2.shape = (300,400)
contour = plt.contour(X,Y,z2,[0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9])
# plt.clabel(contour,fontsize=10,colors=('k'))
plt.show()