【Python】平面多边形Wachspress坐标、中值坐标的计算及其等高线

小周ipython

已于 2022-07-08 12:10:54 修改

阅读量761

点赞数 1

分类专栏： CAGD 文章标签： python 平面几何学

于 2022-06-13 16:08:40 首次发布

本文链接：https://blog.csdn.net/weixin_45718987/article/details/125262073

版权

CAGD 专栏收录该内容

13 篇文章 1 订阅

订阅专栏

平面多边形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()
$\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||}$

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()