一、概述
在做完了 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 catmull−clarksubdivision 之后，紧接着继续做 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′=P⋃Eps
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′=(1−nβ)⋅O+βi=0∑n−1Viwhereβ=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=0∑n−1Vi
这两者是一致的，这样
α
=
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=2n∑i=1nfpi−n⋅O
这样，就与第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)
2898
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
