三角网格曲面参数化

发表于 2020-10-31 分类于图形学阅读次数： Changyan：

本文介绍文章 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
3

clear, close all
model_name = 'data/bunny.obj';
[v,f]=read_mesh(model_name);   % 导入模型

计算边界点

1
2
3

boundary=compute_boundary(f);   % 计算边界点，注意此时的边界是有序的
assert(~isempty(boundary),'仅适用于亏格数为0的开网格！');
MVC = sparse(length(v),length(v));

对非边界点循环计算均值重心坐标

for 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)$参数

%% 指定边界点的(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]';

构造最后一步中的稀疏线性系统

%% (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;

可视化结果

%% 可视化结果
uv(:,3) = 0;
plot_mesh(v,f)
figure
plot_mesh(uv,f)

将绿白相见棋盘格贴图到模型表面可以直观感受参数化的好坏