【Python】Catmull-Clark 细分算法的python实现

最新推荐文章于 2023-03-11 20:14:25 发布

小周ipython

最新推荐文章于 2023-03-11 20:14:25 发布

阅读量776

点赞数 3

分类专栏： CAGD 文章标签：算法 python 几何学

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

版权

CAGD 专栏收录该内容

13 篇文章 1 订阅

订阅专栏

一、资料参考来源

二、Catmull-Clark细分规则

2.1 几何规则

新面点 —— 一个多边形面上旧顶点的均值；
新边点 —— 一条边两端点的均值以及该边所相邻两个面的面点的均值的平均；
新顶点 —— 由一个旧顶点以及其周围的新面点的均值、新边点计算得到，通式为 $\frac{Q}{n}+\frac{2R}{n}+\frac{S(n-3)}{n}$ ，其中，n为与旧顶点相关联的多边形面的数量，Q为旧顶点相关联的多边形面的面点的均值，R为旧顶点所在边的边点的均值，S为旧顶点。

2.2 拓扑规则

将每个新面点连接到对应旧面的边的新边点；
将每个新顶点连接到所有与旧顶点关联的旧边的新边点。

三、编程思路

Step1. 计算新面点；
Step2. 获取边的参数以及边中点；
Step3. 计算内部新边点： $edge\_point = \frac{（边中点+边的邻接面的两个面点的均值）}{2},$ 同时，对于边界处的边的边点计算为 $edge\_point = 边中点$ 即可；
Step4. 根据已知面点求平均面点Q；
Step5. 计算平均边中点R；
Step6. 计算每个点涉及几个平面；
Step7. 首先对原始坐标顶点是边界顶点还是内部顶点做出判断，然后再根据
$Q ， R ，原始坐标点$ 或 $\frac{3}{4}v_0+\frac{1}{8}(v_1+v_2)$ 计算新的坐标顶点；
Step8. 将新边点与新面点添加到新顶点之后，并且给出相应的点的下标列表；
Step9. 根据新面点与新边点在更新点中的下标，以及旧顶点，得到新的面列表；
Step10. 根据新顶点和新面，依次绘制图像，图像绘制方法为：使用循环，将面依次绘制出来。

四、代码

本代码在Catmull-Clark subdivision的Python代码基础上添加了边界处新顶点生成的规则，并且修改了边界处新边点生成的规则，使得代码能够在闭合曲面以及开曲面上均能实现细分。（部分代码已添加我的个人理解）

# 输入的坐标点的下标进行前后大小调整
def adjust_position(pointnum1,pointnum2):
    
    if pointnum1>pointnum2:
        v = pointnum1
        pointnum1 = pointnum2
        pointnum2 = v

    return pointnum1,pointnum2

def center_point(p1, p2):
    """
    returns a point in the center of the
    segment ended by points p1 and p2
    """
    cp = []
    for i in range(3):
        cp.append((p1[i] + p2[i]) / 2)

    return cp


def sum_point(p1, p2):
    """
    adds points p1 and p2
    """
    sp = []
    for i in range(3):
        sp.append(p1[i] + p2[i])

    return sp


def div_point(p, d):
    """
    divide point p by d
    """
    sp = []
    for i in range(3):
        sp.append(p[i] / d)

    return sp


def mul_point(p, m):
    """
    multiply point p by m
    """
    sp = []
    for i in range(3):
        sp.append(p[i] * m)

    return sp


def get_face_points(input_points, input_faces):   # 计算新面点并返回
    """
    From http://rosettacode.org/wiki/Catmull%E2%80%93Clark_subdivision_surface

    1. for each face, a face point is created which is the average of all the points of the face.
    """

    NUM_DIMENSIONS = 3  # 维度


    # face_points will have one point for each face

    face_points = []

    for curr_face in input_faces:  # 对于每一个输入的面的顶点坐标
        face_point = [0.0, 0.0, 0.0]   # 提前预置计算新面点的存储空间
        for curr_point_index in curr_face:
            curr_point = input_points[curr_point_index]
            # add curr_point to face_point
            # will divide later
            for i in range(NUM_DIMENSIONS):
                face_point[i] += curr_point[i]
        # divide by number of points for average
        num_points = len(curr_face)   # 一个面由多少个顶点构成
        for i in range(NUM_DIMENSIONS):
            face_point[i] /= num_points
        face_points.append(face_point)  #　计算出每个面的面点并进行存储

    return face_points


def get_edges_faces(input_points, input_faces):
    """

    Get list of edges and the one or two adjacent faces in a list.
    also get center point of edge

    Each edge would be [pointnum_1, pointnum_2, facenum_1, facenum_2, center]

    """

    # will have [pointnum_1, pointnum_2, facenum]

    edges = []

    # get edges from each face

    for facenum in range(len(input_faces)):
        face = input_faces[facenum]  # 获取第 facenum 个面
        num_points = len(face)
        # loop over index into face
        for pointindex in range(num_points):
            # 如果不是最后一个点，则边是当前点和下一个点的连线
            if pointindex < num_points - 1:
                pointnum_1 = face[pointindex]
                pointnum_2 = face[pointindex + 1]
            else:
                # 如果是最后一个点，则边是当前点和第一个点的连线
                pointnum_1 = face[pointindex]
                pointnum_2 = face[0]
            # order points in edge by lowest point number
            # 给点换位置，下标小的点在前，下标大的点在后，在最后的新面的获取中有用，防止数据混乱
            if pointnum_1 > pointnum_2:
                temp = pointnum_1
                pointnum_1 = pointnum_2
                pointnum_2 = temp
            edges.append([pointnum_1, pointnum_2, facenum])

    # sort edges by pointnum_1, pointnum_2, facenum
    # 对边的列表集合进行顺序排序，先按照 pointnum_1 排序，再按照 pointnum_2 排序，最后按照 facenum 排序
    edges = sorted(edges)

    # merge edges with 2 adjacent faces
    # [pointnum_1, pointnum_2, facenum_1, facenum_2] or
    # [pointnum_1, pointnum_2, facenum_1, None]

    num_edges = len(edges)
    eindex = 0
    merged_edges = []

    while eindex < num_edges:
        e1 = edges[eindex]
        # check if not last edge
        if eindex < num_edges - 1:
            e2 = edges[eindex + 1]
            if e1[0] == e2[0] and e1[1] == e2[1]:
                merged_edges.append([e1[0], e1[1], e1[2], e2[2]])
                eindex += 2
            else:
                merged_edges.append([e1[0], e1[1], e1[2], None])
                eindex += 1
        else:
            merged_edges.append([e1[0], e1[1], e1[2], None])
            eindex += 1

    # add edge centers

    edges_centers = []

    for me in merged_edges:
        p1 = input_points[me[0]]
        p2 = input_points[me[1]]
        cp = center_point(p1, p2)
        edges_centers.append(me + [cp])

    return edges_centers


def get_edge_points(input_points, edges_faces, face_points):
    """
    for each edge, an edge point is created which is the average
    between the center of the edge and the center of the segment made
    with the face points of the two adjacent faces.
    """

    edge_points = []

    for edge in edges_faces:
        # get center of edge
        cp = edge[4]
        # get center of two facepoints
        fp1 = face_points[edge[2]]
        # if not two faces just use one facepoint
        # should not happen for solid like a cube
#         if edge[3] == None:
#             fp2 = fp1
#         else:
#             fp2 = face_points[edge[3]]
#         cfp = center_point(fp1, fp2)
#         # get average between center of edge and
#         # center of facepoints
#         edge_point = center_point(cp, cfp)
#         edge_points.append(edge_point)
        if edge[3] == None:
            edge_points.append(cp)
        else:
            fp2 = face_points[edge[3]]
            cfp = center_point(fp1, fp2)
            # get average between center of edge and
            # center of facepoints
            edge_point = center_point(cp, cfp)
            edge_points.append(edge_point)

    return edge_points


def get_avg_face_points(input_points, input_faces, face_points):
    """

    for each point calculate

    the average of the face points of the faces the point belongs to (avg_face_points)

    create a list of lists of two numbers [facepoint_sum, num_points] by going through the
    points in all the faces.

    then create the avg_face_points list of point by dividing point_sum (x, y, z) by num_points

    """

    # initialize list with [[0.0, 0.0, 0.0], 0]

    num_points = len(input_points)

    temp_points = []

    for pointnum in range(num_points):
        temp_points.append([[0.0, 0.0, 0.0], 0])

    # loop through faces updating temp_points

    for facenum in range(len(input_faces)):
        fp = face_points[facenum]
        for pointnum in input_faces[facenum]:
            tp = temp_points[pointnum][0]
            temp_points[pointnum][0] = sum_point(tp, fp)  # 计算一个点所连接的所有面的面点之和
            temp_points[pointnum][1] += 1  # 统计一个点与多少个面相接
            
    # divide to create avg_face_points
    avg_face_points = []

    for tp in temp_points:
        afp = div_point(tp[0], tp[1])
        avg_face_points.append(afp)

    return avg_face_points


def get_avg_mid_edges(input_points, edges_faces):
    """

    the average of the centers of edges the point belongs to (avg_mid_edges)

    create list with entry for each point
    each entry has two elements. one is a point that is the sum of the centers of the edges
    and the other is the number of edges. after going through all edges divide by
    number of edges.

    """

    # initialize list with [[0.0, 0.0, 0.0], 0]

    num_points = len(input_points)

    temp_points = []

    for pointnum in range(num_points):
        temp_points.append([[0.0, 0.0, 0.0], 0])

    # go through edges_faces using center updating each point

    for edge in edges_faces:
        cp = edge[4]
        for pointnum in [edge[0], edge[1]]:
            tp = temp_points[pointnum][0]
            temp_points[pointnum][0] = sum_point(tp, cp)
            temp_points[pointnum][1] += 1

    # divide out number of points to get average

    avg_mid_edges = []

    for tp in temp_points:
        ame = div_point(tp[0], tp[1])
        avg_mid_edges.append(ame)

    return avg_mid_edges


def get_points_faces(input_points, input_faces):
    # initialize list with 0

    num_points = len(input_points)

    points_faces = []

    for pointnum in range(num_points):
        points_faces.append(0)  # 提前预设存储空间，方便计算

    # loop through faces updating points_faces
    # 使用循环统计每个点属于多少个面
    for facenum in range(len(input_faces)):
        for pointnum in input_faces[facenum]:
            points_faces[pointnum] += 1
    
    # 返回每个点与多少个面相接的列表
    return points_faces


'''
添加顶点更新边界条件
'''
# 获取顶点的邻接点
def get_ad_point(input_points, edges_faces):
    n = len(input_points)
    ad_points = []  # 用于存储所有顶点的邻接点
    for pointnum in range(n):
        list1 = []   # 暂时存放某个顶点的邻接点
        for edge in edges_faces:
            if pointnum==edge[0]:
                list1.append(edge[1])
            elif pointnum==edge[1]:
                list1.append(edge[0])
        ad_points.append(list1)
        
    return ad_points


# 获取顶点是边界点还是内部点，内部点为False，边界点为True
def get_point_bool(input_points, edges_faces):
    # 获取邻接点
    ad_points = get_ad_point(input_points, edges_faces)
    
    n = len(input_points)
    point_bol = []   # 用于存放顶点是内部点还是边界点的布尔值
    for i in range(n):
        point_bol.append(False)
    
    v0_ad_point = []  # 用于存放每个v0的两个边界邻点
    for pointnum in range(n):
        v0_ad = ad_points[pointnum]
        v0 = pointnum
        v0_ad_list = []  # 用于存放v0的两个边界邻点
        for vi in v0_ad:
            pointnum1, pointnum2 = adjust_position(v0, vi)
            for edge in edges_faces:
                if pointnum1==edge[0] and pointnum2==edge[1]:
                    if None in edge:
                        point_bol[pointnum]=True
                        v0_ad_list.append(vi)
        v0_ad_point.append(v0_ad_list)
    
    return point_bol, v0_ad_point

'''
添加顶点更新边界条件
'''


def get_new_points(input_points, edges_faces, points_faces, avg_face_points, avg_mid_edges):
    
    point_bol, v0_ad_point = get_point_bool(input_points, edges_faces)
#     print(point_bol)
    """

    m1 = (n - 3.0) / n
    m2 = 1.0 / n
    m3 = 2.0 / n
    new_coords = (m1 * old_coords)
               + (m2 * avg_face_points)
               + (m3 * avg_mid_edges)

    """

    new_points = []

    for pointnum in range(len(input_points)):
        old_coords = input_points[pointnum]
        if point_bol[pointnum]:
            v_x = 3/4*old_coords[0]
            v_y = 3/4*old_coords[1]
            v_z = 3/4*old_coords[2]
            for j in v0_ad_point[pointnum]:
                vi = input_points[j]
                v_x+=1/8*vi[0]
                v_y+=1/8*vi[1]
                v_z+=1/8*vi[2]
            new_coords = [v_x, v_y, v_z]
        else:
            n = points_faces[pointnum]
            m1 = (n - 3.0) / n
            m2 = 1.0 / n
            m3 = 2.0 / n
            p1 = mul_point(old_coords, m1)
            afp = avg_face_points[pointnum]
            p2 = mul_point(afp, m2)
            ame = avg_mid_edges[pointnum]
            p3 = mul_point(ame, m3)
            p4 = sum_point(p1, p2)
            new_coords = sum_point(p4, p3)

        new_points.append(new_coords)

    return new_points


def switch_nums(point_nums):
    """
    Returns tuple of point numbers
    sorted least to most
    """
    if point_nums[0] < point_nums[1]:
        return point_nums
    else:
        return (point_nums[1], point_nums[0])


def cmc_subdiv(input_points, input_faces):
    # 1. for each face, a face point is created which is the average of all the points of the face.
    # each entry in the returned list is a point (x, y, z).
    # 计算新面点，返回结果的顺序与给定的input_faces的顺序是一样的，防止求新边点时数据混乱
    face_points = get_face_points(input_points, input_faces) 

    # get list of edges with 1 or 2 adjacent faces，并且计算边的中点
    
    # facenum_1代表第一个邻接面，facenum_2代表第二个邻接面
    # pointnum_1代表边的一个端点，pointnum_2代表边的另一个端点
    # edges_faces 格式如下
    # [pointnum_1, pointnum_2, facenum_1, facenum_2, center] or
    # [pointnum_1, pointnum_2, facenum_1, None, center]
    edges_faces = get_edges_faces(input_points, input_faces)

    # get edge points, a list of points
    # 计算新边点
    edge_points = get_edge_points(input_points, edges_faces, face_points)
    
    # the average of the face points of the faces the point belongs to (avg_face_points)
    # 计算平均面点，求 Q
    avg_face_points = get_avg_face_points(input_points, input_faces, face_points)
    
    # the average of the centers of edges the point belongs to (avg_mid_edges)
    # 计算平均边中点，求R
    avg_mid_edges = get_avg_mid_edges(input_points, edges_faces)

    # how many faces a point belongs to
    # 一个点属于多少面
    points_faces = get_points_faces(input_points, input_faces)

    """

    m1 = (n - 3) / n
    m2 = 1 / n
    m3 = 2 / n
    new_coords = (m1 * old_coords)
               + (m2 * avg_face_points)
               + (m3 * avg_mid_edges)

    """

    new_points = get_new_points(input_points, edges_faces, points_faces, avg_face_points, avg_mid_edges)

    """

    Then each face is replaced by new faces made with the new points,

    for a triangle face (a,b,c):
       (a, edge_point ab, face_point abc, edge_point ca)
       (b, edge_point bc, face_point abc, edge_point ab)
       (c, edge_point ca, face_point abc, edge_point bc)

    for a quad face (a,b,c,d):
       (a, edge_point ab, face_point abcd, edge_point da)
       (b, edge_point bc, face_point abcd, edge_point ab)
       (c, edge_point cd, face_point abcd, edge_point bc)
       (d, edge_point da, face_point abcd, edge_point cd)

    face_points is a list indexed by face number so that is
    easy to get.

    edge_points is a list indexed by the edge number
    which is an index into edges_faces.

    need to add face_points and edge points to 
    new_points and get index into each.

    then create two new structures

    face_point_nums - list indexes by facenum
    whose value is the index into new_points

    edge_point num - dictionary with key (pointnum_1, pointnum_2)
    and value is index into new_points

    """

    # add face points to new_points
    # face_point_nums 存储面点在 new_points 中的下标
    face_point_nums = []

    # point num after next append to new_points
    next_pointnum = len(new_points)

    for face_point in face_points:
        new_points.append(face_point)
        face_point_nums.append(next_pointnum)
        next_pointnum += 1

    # add edge points to new_points
    # edge_point_nums 采用字典的key-value方式存储边点在 new_points 中的下标
    edge_point_nums = dict()

    for edgenum in range(len(edges_faces)):
        pointnum_1 = edges_faces[edgenum][0]
        pointnum_2 = edges_faces[edgenum][1]
        edge_point = edge_points[edgenum]
        new_points.append(edge_point)
        edge_point_nums[(pointnum_1, pointnum_2)] = next_pointnum
        next_pointnum += 1

    # new_points now has the points to output. Need new
    # faces

    """

    just doing this case for now:

    for a quad face (a,b,c,d):
       (a, edge_point ab, face_point abcd, edge_point da)
       (b, edge_point bc, face_point abcd, edge_point ab)
       (c, edge_point cd, face_point abcd, edge_point bc)
       (d, edge_point da, face_point abcd, edge_point cd)

    new_faces will be a list of lists where the elements are like this:

    [pointnum_1, pointnum_2, pointnum_3, pointnum_4]

    """

    new_faces = []

    for oldfacenum in range(len(input_faces)):
        oldface = input_faces[oldfacenum]
        # 4 point face
        if len(oldface) == 4:
            a = oldface[0]
            b = oldface[1]
            c = oldface[2]
            d = oldface[3]
            face_point_abcd = face_point_nums[oldfacenum]
            edge_point_ab = edge_point_nums[switch_nums((a, b))]
            edge_point_da = edge_point_nums[switch_nums((d, a))]
            edge_point_bc = edge_point_nums[switch_nums((b, c))]
            edge_point_cd = edge_point_nums[switch_nums((c, d))]
            new_faces.append((a, edge_point_ab, face_point_abcd, edge_point_da))
            new_faces.append((b, edge_point_bc, face_point_abcd, edge_point_ab))
            new_faces.append((c, edge_point_cd, face_point_abcd, edge_point_bc))
            new_faces.append((d, edge_point_da, face_point_abcd, edge_point_cd))
            
        elif len(oldface) == 3:
            a = oldface[0]
            b = oldface[1]
            c = oldface[2]
            face_point_abc = face_point_nums[oldfacenum]
            edge_point_ab = edge_point_nums[switch_nums((a, b))]
            edge_point_ca = edge_point_nums[switch_nums((c, a))]
            edge_point_bc = edge_point_nums[switch_nums((b, c))]
            new_faces.append((a, edge_point_ab, face_point_abc, edge_point_ca))
            new_faces.append((b, edge_point_bc, face_point_abc, edge_point_ab))
            new_faces.append((c, edge_point_ca, face_point_abc, edge_point_bc))

    return new_points, new_faces


def graph_output(output_points, output_faces, fig, m):

    ax = fig.add_subplot(111, projection='3d')


    """

    Plot each face

    """
    ax.set_title('{0}次细分'.format(str(m)),fontsize=18)
    for facenum in range(len(output_faces)):
        curr_face = output_faces[facenum]
        xcurr = []
        ycurr = []
        zcurr = []
        for pointnum in range(len(curr_face)):
            xcurr.append(output_points[curr_face[pointnum]][0])
            ycurr.append(output_points[curr_face[pointnum]][1])
            zcurr.append(output_points[curr_face[pointnum]][2])
        xcurr.append(output_points[curr_face[0]][0])
        ycurr.append(output_points[curr_face[0]][1])
        zcurr.append(output_points[curr_face[0]][2])
        ax.plot(xcurr, ycurr, zcurr , color='b')
#         break

# #　正方体
# # input_points 代表输入的三维坐标点
# input_points = [
#     [-1.0, 1.0, 1.0],
#     [-1.0, -1.0, 1.0],
#     [1.0, -1.0, 1.0],
#     [1.0, 1.0, 1.0],
#     [1.0, -1.0, -1.0],
#     [1.0, 1.0, -1.0],
#     [-1.0, -1.0, -1.0],
#     [-1.0, 1.0, -1.0]
# ]

# # input_faces 代表输入的面的各个顶点的下标，通过面顶点的下标来找面
# input_faces = [
#     [0, 1, 2, 3],
#     [3, 2, 4, 5],
#     [5, 4, 6, 7],
#     [7, 0, 3, 5],
#     [7, 6, 1, 0],
#     [6, 1, 2, 4],
# ]

# # 梯形
# # input_points 代表输入的三维坐标点
# input_points = [
#     [-1.0, 1.0, 3.0],
#     [-1.0, -1.0, 3.0],
#     [1.0, -1.0, 1.0],
#     [1.0, 1.0, 1.0],
#     [1.0, -1.0, -1.0],
#     [1.0, 1.0, -1.0],
#     [-1.0, -1.0, -1.0],
#     [-1.0, 1.0, -1.0]
# ]

# # input_faces 代表输入的面的各个顶点的下标，通过面顶点的下标来找面
# input_faces = [
#     [0, 1, 2, 3],
#     [3, 2, 4, 5],
#     [5, 4, 6, 7],
#     [7, 0, 3, 5],
#     [7, 6, 1, 0],
#     [6, 1, 2, 4],
# ]

# # 锥体
# # input_points 代表输入的三维坐标点
# input_points = [
#     [1,-1,0],
#     [1,1,0],
#     [-1,0,0],
#     [0,0,1],
# ]

# # input_faces 代表输入的面的各个顶点的下标，通过面顶点的下标来找面
# input_faces = [
#     [0,1,2],
#     [0,1,3],
#     [0,2,3],
#     [1,2,3],
# ]

# # 平面
# # input_points 代表输入的三维坐标点
# input_points = [
#     [0.0,0.0,0.0],
#     [1.0,0.0,0.0],
#     [2.0,0.0,0.0],
#     [2.0,1.0,0.0],
#     [1.0,1.0,0.0],
#     [0.0,1.0,0.0],

# ]

# # input_faces 代表输入的面的各个顶点的下标，通过面顶点的下标来找面
# input_faces = [
#     [0,1,4,5],
#     [1,2,3,4],
# ]

#　曲面
input_points = [
    [1.0,0.0,0.0],
    [2.0,0.0,0.0],
    [2.0,1.0,0.0],
    [1.0,1.0,0.0],
    [0.0,1.0,1.0],
    [0.0,0.0,1.0]
]

input_faces = [
    [0,1,2,3],
    [0,3,4,5]
]


iterations = 4

plt.ion()  # 维持多张图展示，不因为plt.show()而停止展示
output_points, output_faces = input_points, input_faces
fig = plt.figure(1,figsize=(10,8))
plt.clf()
graph_output(output_points, output_faces, fig, 0)
plt.pause(1)

for i in range(iterations):
    output_points, output_faces = cmc_subdiv(output_points, output_faces)

    fig = plt.figure(1,figsize=(10,8))
    plt.clf()
    graph_output(output_points, output_faces, fig, i+1)
    plt.pause(1)
#     break