三角网格曲面参数化
本文介绍文章 Parametrization and smooth approximation of surface triangulations 中的保形参数化算法,并给出MATLAB程序 曲面参数化问题:对于给定的三角网格曲面$\mathbf{S}$,找到一个映射$\phi(u,v): \mathbb{R}^2\longmapsto\mathbb{R}^3$,使得平面参数域中的点与曲面网格的点一一对应,求解参数$(u,v)$,即$\phi^{-1}$的过程称为参数化。
算法描述
对于三角网格曲面$\mathbf{S}$的每个内部点$x_i$,寻找1领域内的邻接点,将其记为$x_{j_1},x_{j_2},\dots,x_{j_m}$,定义$x_i$的子图的顶点集$X_i={x_i,x_{j_1},x_{j_2},\dots,x_{j_m}}$表示,如上图左侧所示。
将步骤1中的$X_i$通过某种方式投影到平面域$\mathbb{R}^2$中,得到平面域中的图$P_i={p,p_1,p_2,\dots,p_m}$。上图中间
这里的投影方式有多种选择,文章提及到最小二乘面和平均法向面,但是当曲面中三角面片之间的夹角过大时,这两种方法不具有稳定性。于是介绍了一种W-W提出的测地极坐标映射的方法。
$$ r_k=\Vert p_k-p\Vert =\Vert x_{j_k}-x_i\Vert $$
$$ \theta_k = ang(p_k,p,p_{k+1})=2\pi \frac{ang(x_{j_k},x_i,x_{j_{k+1}})}{\sum_kang(x_{j_k},x_i,x_{k_{k+1}})} $$ 这里$ang(a,b,c)$表示向量$\vec {ba}$与向量$\vec {bc}$之间的角度。实际上上面第二个等式就是将空间角度等比例放缩到平面中,因为平面周角为$2\pi$,所以前面需要乘$2\pi$。定义$p=(0,0)$,$p_1=(\Vert x_{j_1}-x_i\Vert,0)$,这里$p$相当于极点,$pp_1$相当于极轴,根据极坐标$(r_k,\theta_k)$则可计算出每个$p_k$。则得到平面域中的多边形${p,p_1,p_2,\dots,p_m}$及其内部一点$p$
计算$p$关于该平面多边形的广义重心坐标。这里采用均值重心坐标[1] 均值重心坐标的计算方式如下:
c++实现均值重心坐标可参考我的另一篇文章:均值重心坐标Mean Value Coordinates In 2D代码实现
在第3步可计算出投影到平面域中的每个内部点关于1-ring点的广义重心坐标。 记:$\lambda_i=[\dots,0,\dots,\lambda_{j_a},\dots,\lambda_{j_b},\dots,0,\dots]$,即只是$x_i$的邻接点对应序号中的元素为第三步计算出的重心坐标值,其余位置为0
固定边界点对应的参数(如参数化到矩形域,则将网格曲面上边界点参数有序定义为矩形域的边界) 假定后$N$个点为边界点,初始化${(u_i,v_i)}_{i=n+1}^N$
记$\Lambda=[\lambda_1,\lambda_2,\dots,\lambda_n]^{\rm T}$,求稀疏解线性方程组:
MATLAB代码
下面代码仅使用与开网格 建议在桌面端浏览器或大屏移动端查看代码
首先导入模型文件:
1
2
3clear, close all
model_name = 'data/bunny.obj';
[v,f]=read_mesh(model_name); % 导入模型计算边界点
1
2
3boundary=compute_boundary(f); % 计算边界点,注意此时的边界是有序的
assert(~isempty(boundary),'仅适用于亏格数为0的开网格!');
MVC = sparse(length(v),length(v));对非边界点循环 计算均值重心坐标
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37for i = setdiff(1:length(v),boundary) % 对非边界点循环 i
[~,indexes_of_adjacent_faces] = find(f==i); % 找到x_i的邻接面
adjacent_vertices = setdiff(unique(f(:,indexes_of_adjacent_faces)),i); % x_i的邻接点
%% 对空间中的x_i的无序邻接点排序
djacent_vertices_ordered = zeros(1,length(adjacent_vertices));
djacent_vertices_ordered(1) = adjacent_vertices(1); % 指定一个初始的
adjacent_faces = f(:,indexes_of_adjacent_faces); % 取出邻接面
for ii = 2:length(adjacent_vertices)
[~,y] = find(adjacent_faces == djacent_vertices_ordered(ii-1));
p = setdiff(unique(adjacent_faces(:,y)),[djacent_vertices_ordered i]);
djacent_vertices_ordered(ii) = p(1);
end
%% 投影到平面
adjacent_edges = v(:,djacent_vertices_ordered) - repmat(v(:,i),1,length(adjacent_vertices)); % X_J-X_i
adjacent_edges_1 = [adjacent_edges(:,2:end) adjacent_edges(:,1)]; % X_{J+1}-X_I
length_of_adjacent_edges = vecnorm(adjacent_edges,2); % 空间中邻接边的长度
cos_value = dot(adjacent_edges,adjacent_edges_1)./...
length_of_adjacent_edges./[length_of_adjacent_edges(2:end) length_of_adjacent_edges(1)];
beta = acos(cos_value); % 空间中的角度beta
alpha = beta/sum(beta)*(2*pi); % 计算对应的平面中的角度alpha
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 计算平面多边形,实际上用MVC的话,就不再需要多边形的具体值了,
% 因为只需要角度和长度,在上步骤已经得出
theta = alpha*triu(ones(length(alpha))-eye(length(alpha)));
polygon_in_plane = length_of_adjacent_edges.*[cos(theta);sin(theta)];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 广义重心坐标:均值重心坐标
tan_value = tan(alpha/2);
MVC_i = ([tan_value(end) tan_value(1:end-1)]+tan_value)./...
length_of_adjacent_edges;
MVC(i,adjacent_vertices) = MVC_i/(sum(MVC_i));
end筛选出有用的MVC
1
2% 只保留内部点的MVC
MVC = MVC(any(MVC, 2),:);初始化边界点的$(u,v)$参数
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25%% 指定边界点的(u,v)
% 假定考虑参数化到[0,1]X[0,1]的正方形中
chord = vecnorm(v(:,[boundary(2:end) boundary(1)])-v(:,boundary),2);
chord = chord/sum(chord)*4;
uv_ = chord*triu(ones(length(chord))-eye(length(chord)));
uv_1(1,:) = uv_(uv_ < 1 & uv_ >= 0); uv_1(2,:) = 0;
uv_2(2,:) = uv_(uv_ < 2 & uv_ >= 1)-1; uv_2(1,:) = 1;
if uv_2(2,1) < 1-uv_1(1,end)
uv_2(2,1) = 0;
else
uv_1(1,end) = 1;
end
uv_3(1,:) = 3-uv_(uv_ < 3 & uv_ >= 2); uv_3(2,:) = 1;
if 1-uv_3(1,1) < 1-uv_2(2,end)
uv_3(1,1) = 1;
else
uv_2(2,end) = 1;
end
uv_4(2,:) = 4-uv_(uv_ <= 4 & uv_ >= 3); uv_4(1,:) = 0;
if 1-uv_4(2,1) < uv_3(1,end)
uv_4(2,1) = 1;
else
uv_3(1,end) = 0;
end
uv_boundary = [uv_1 uv_2 uv_3 uv_4]';构造最后一步中的稀疏线性系统
1
2
3
4
5
6
7
8%% (I-Lambda_1)*uv_int = Lambda_2*uv_boundary
Lambda_1 = MVC(:,setdiff(1:length(v),boundary)); % 构造Lambda_1
Lambda_2 = MVC(:,boundary); % 构造Lambda_2
A = eye(size(Lambda_1)) - Lambda_1; % I-Lambda_1
uv_Int = A\Lambda_2*uv_boundary; % 解线性方程组
uv(setdiff(1:length(v),boundary),:) = uv_Int;
uv(boundary,:) = uv_boundary;可视化结果
1
2
3
4
5%% 可视化结果
uv(:,3) = 0;
plot_mesh(v,f)
figure
plot_mesh(uv,f)
将绿白相见棋盘格贴图到模型表面可以直观感受参数化的好坏