Loop 细分曲面(loop subdivision)详解,附Python完整代码

一、概述

在做完了 c a t m u l l − c l a r k s u b d i v i s i o n catmull-clark\quad subdivision catmullclarksubdivision 之后,紧接着继续做 l o o p s u b d i v i s i o n loop\quad subdivision loopsubdivision,但是跟之前一样,网上的帖子都是直接给出了各个控制点新增或更新时,所有相关点的权重,并没有给出具体过程是怎么做的,而且在做 l o o p loop loop细分 的时候并没有像之前一样有完整的“保真代码”(能够体现出具体操作步骤的代码),所以也没有比较直接的方法去了解每一步的过程,所以只能去看最原始的 l o o p s u b d i v i s i o n loop \quad subdivision loopsubdivision 论文了,原文链接为Smooth Subdivision Surfaces Based on Triangles,这篇论文是于1987年撰写的一篇硕士论文,而 catmull细分 的原始论文于 1978 年撰写,链接为Recursively generated B-spline surfaces on arbitrary topological meshes,那时候的论文是多么的朴实无华,看完之后真的受益良多,非常推荐大家看一下。同样,这里只记录怎么做,不讨论为什么,具体理论推导见原始论文吧。

二、Loop Subdivision 过程

loop细分曲面过程就是插入边点和更新原始点的过程,设每次细分之后的点集为 P ′ P' P, 新插入的点集为 E p s Eps Eps,原始点集更新位置后的点集为 P P P,那么:
P ′ = P ⋃ E p s P'=P\bigcup Eps P=PEps

1. 计算边点(edge_points)

在这里插入图片描述
在论文中提到,利用如上图所示的mask对网格中所有点进行遍历,对每一条边计算得到新的 e d g e _ p o i n t edge\_point edge_point,即:
e p 2 = 1 8 ( V 0 + V 2 ) + 3 8 ( V 1 + V 3 ) ep_2=\frac{1}{8}(V_0+V_2)+\frac{3}{8}(V_1+V_3) ep2=81(V0+V2)+83(V1+V3)
其中 e p 2 ep_2 ep2表示该边属于两个面。另外,原始论文中并没有提到某条边只属于一个面的情况,但是目前看到的帖子里,都说如果一个边只属于一个面,那么这条边新增的点为该边中点,(假设上图中,边 E ( V 1 , V 3 ) E(V_1,V_3) E(V1,V3)只属于一个面)则:
e p 1 = 1 2 ( V 1 + V 3 ) ep_1=\frac{1}{2}(V_1+V_3) ep1=21(V1+V3)
根据上述对边点的计算公式,我将此边点计算方法与catmull-clark方法进行了类比,理解得到了一个中间过程。

catmull-clark方法在计算边点时,先计算了每个四点面的中心点,再基于中心点与边的中点计算得到最终的边点,详细过程见Catmull 细分曲面 (Catmull-Clark subdivision) 详解 附Python 完整代码这篇文章。

根据这个思想,我先对每条边属于的两个三角面所组成的四点面计算面点 f a c e _ p o i n t ( f p ) face\_point(fp) face_point(fp),以及每条边的中点 e d g e _ c e n t e r ( e c ) edge\_center(ec) edge_center(ec),最后求这两点的平均点即为 e d g e _ p o i n t ( e p ) edge\_point(ep) edge_point(ep)。过程如下图所示。
在这里插入图片描述
这样来推导一下,最后每个点的权重是否符合前面所说的mask所示的分布。

(1)计算 f p fp fp的公式为:
f p = 1 4 ( V 0 + V 1 + V 2 + V 3 ) fp=\frac{1}{4}(V_0+V_1+V_2+V_3) fp=41(V0+V1+V2+V3)
(2)计算 e c ec ec的公式为:
e c = 1 2 ( V 1 + V 3 ) ec=\frac{1}{2}(V_1+V_3) ec=21(V1+V3)
(3)计算 e p ep ep的公式为:
e p = 1 2 ( f p + e c ) = 1 2 ( 1 4 ( V 0 + V 1 + V 2 + V 3 ) + 1 2 ( V 1 + V 3 ) ) = 1 8 ( V 0 + V 2 ) + 3 8 ( V 1 + V 3 ) \begin{aligned} ep&=\frac{1}{2}(fp+ec) \\&=\frac{1}{2}(\frac{1}{4}(V_0+V_1+V_2+V_3)+\frac{1}{2}(V_1+V_3)) \\&=\frac{1}{8}(V_0+V_2)+\frac{3}{8}(V_1+V_3) \end{aligned} ep=21(fp+ec)=21(41(V0+V1+V2+V3)+21(V1+V3))=81(V0+V2)+83(V1+V3)
正好,权重分布规则一致,这样一来,如果对于只属于一个面的边来说,中点就是边点,只需要把面点设置为 ( 0 , 0 , 0 ) (0, 0, 0) (0,0,0)即可,实现的时候就比较方便了,感觉更有条理了。

2. 更新原始点

(1)更新规则
在这里插入图片描述
对于每个点,利用如上图所示的mask进行计算,来更新原始点的位置坐标,设点 O ′ O' O O O O更新后的点,则更新公式如下:
O ′ = ( 1 − n β ) ⋅ O + β ∑ i = 0 n − 1 V i w h e r e β = 1 n [ 5 8 − ( 3 8 + 1 4 cos ⁡ 2 π n ) 2 ] O'=(1-n\beta)\cdot O+\beta \sum_{i=0}^{n-1}V_i \\where\quad \beta=\frac{1}{n}[\frac{5}{8}-(\frac{3}{8}+\frac{1}{4}\cos\frac{2\pi}{n})^2] O=(1nβ)O+βi=0n1Viwhereβ=n1[85(83+41cosn2π)2]
这里的 n n n是与上图所示的 k k k是一致的。

实际上,原论文中的计算公式为:
O ′ = ( 1 − α ) ⋅ O + α ⋅ Q O'=(1-\alpha)\cdot O+\alpha\cdot Q O=(1α)O+αQ
其中
Q = 1 n ∑ i = 0 n − 1 V i Q=\frac{1}{n}\sum_{i=0}^{n-1}V_i Q=n1i=0n1Vi
这两者是一致的,这样
α = n β = 5 8 − ( 3 8 + 1 4 cos ⁡ 2 π n ) 2 \alpha=n\beta=\frac{5}{8}-(\frac{3}{8}+\frac{1}{4}\cos\frac{2\pi}{n})^2 α=nβ=85(83+41cosn2π)2

(2)实现方法

因为在第1步中计算了每相邻两个三角面构成的四点面的面点,那么在更新原始点时可利用下图所示方法来进行计算:
在这里插入图片描述
计算公式则变为:
O ′ = ( 1 − α ) ⋅ O + α ⋅ Q O'=(1-\alpha)\cdot O+\alpha\cdot Q O=(1α)O+αQ
则利用面点 f p fp fp来计算 Q Q Q
Q = ∑ i = 1 n f p i − n ⋅ O 2 n Q=\frac{\sum_{i=1}^nfp_i-n\cdot O}{2n} Q=2ni=1nfpinO
这样,就与第1步串联起来了,实现起来也更加方便。

三、实验结果

输入数据为一个正四面体:
在这里插入图片描述
细分1次:
在这里插入图片描述
细分2次:
在这里插入图片描述
细分3次:
在这里插入图片描述
细分4次:
在这里插入图片描述

四、完整代码


from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt
import numpy as np
import sys
import time
import math


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 sub_point(p1, 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_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]
        num_points = len(face)
        # loop over index into face
        for pointindex in range(num_points):
            # if not last point then edge is curr point and next point
            if pointindex < num_points - 1:
                pointnum_1 = face[pointindex]
                pointnum_2 = face[pointindex + 1]
            else:
                # for last point edge is curr point and first point
                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

    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_faces_points(input_points, input_faces, edges_faces):
    num_edges = len(edges_faces)
    faces_points = []

    for edge in edges_faces:
        if edge[3] != None:
            faces_point = [0.0, 0.0, 0.0]
            for face_num in [edge[2],edge[3]]:
                for point_num in input_faces[face_num]:
                    if point_num != edge[0] and point_num != edge[1]:
                        faces_point = sum_point(faces_point,mul_point(input_points[point_num], 0.25))
                    else:
                        faces_point = sum_point(faces_point,mul_point(input_points[point_num], 0.125))
            faces_points.append(faces_point)
            print(faces_point)
        else:
            faces_point = [0.0, 0.0, 0.0]
            faces_points.append(faces_point)

    return faces_points


def get_edge_points(edges_faces, faces_points):
    edge_points = []
    for i in range(len(edges_faces)):
        if edges_faces[i][3] != None:
            edge_points.append(center_point(edges_faces[i][4], faces_points[i]))
        else:
            edge_points.append(edges_faces[i][4])

    return edge_points


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][0] += 1
            points_faces[pointnum][1].append(facenum)

    return points_faces
#统计每个点属于多少个面

def get_face_sum(input_points, input_faces):
    face_sum = []
    for face in input_faces:
        p = sum_point(input_points[face[2]],sum_point(input_points[face[0]],input_points[face[1]]))
        face_sum.append(p)
    return face_sum

def get_new_points(input_points, points_faces, face_sum):
    new_points = []

    for i in range(len(input_points)):
        q = [0.0, 0.0, 0.0]
        beta = 5.0/8-pow(3.0/8+1.0/4*math.cos(2*math.pi/points_faces[i][0]),2)
        for point_face_num in points_faces[i][1]:
            q = sum_point(q, face_sum[point_face_num])
        q = sub_point(q, mul_point(input_points[i], points_faces[i][0]))
        q = div_point(q, 2*points_faces[i][0])
        v = mul_point(input_points[i],1-1*beta)
        new_points.append(sum_point(v, mul_point(q, 1*beta)))

    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 loop_subdiv(input_points, input_faces):
    edges_faces = get_edges_faces(input_points, input_faces)

    faces_points = get_faces_points(input_points, input_faces, edges_faces)
    #print(len(faces_points))
    edge_points = get_edge_points(edges_faces, faces_points)

    print(edge_points)
    points_faces = get_points_faces(input_points, input_faces)

    face_sum = get_face_sum(input_points, input_faces)

    new_points = get_new_points(input_points, points_faces, face_sum)

    next_pointnum = len(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_faces = []

    for oldfacenum in range(len(input_faces)):
        oldface = input_faces[oldfacenum]
        # 3 point face
        if len(oldface) == 3:
            a = oldface[0]
            b = oldface[1]
            c = oldface[2]

            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, edge_point_ca))
            new_faces.append((b, edge_point_bc, edge_point_ab))
            new_faces.append((c, edge_point_ca, edge_point_bc))
            new_faces.append((edge_point_ca, edge_point_ab, edge_point_bc))


    return new_points, new_faces



def graph_output(output_points, output_faces, fig):

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


    """

    Plot each face

    """

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


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 = [
    [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_faces = [
    [0, 3, 1],
    [0, 3, 2],
    [0, 2, 1],
    [3, 2, 1],
]


iterations = 4

plt.ion()
output_points, output_faces = input_points, input_faces


for i in range(iterations):
    output_points, output_faces = loop_subdiv(output_points, output_faces)
    fig = plt.figure(1)
    plt.clf()
    graph_output(output_points, output_faces, fig)
    plt.pause(1)
评论 5
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值