三角网格曲面参数化

本文介绍文章 Parametrization and smooth approximation of surface triangulations 中的保形参数化算法,并给出MATLAB程序 曲面参数化问题:对于给定的三角网格曲面\(\mathbf{S}\),找到一个映射\(\phi(u,v): \mathbb{R}^2\longmapsto\mathbb{R}^3\),使得平面参数域中的点与曲面网格的点一一对应,求解参数\((u,v)\),即\(\phi^{-1}\)的过程称为参数化。


算法描述

  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}}\)表示,如上图左侧所示。

  2. 将步骤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\)

  1. 计算\(p\)关于该平面多边形的广义重心坐标。这里采用均值重心坐标[1] 均值重心坐标的计算方式如下:

    c++实现均值重心坐标可参考我的另一篇文章:均值重心坐标Mean Value Coordinates In 2D代码实现

  2. 在第3步可计算出投影到平面域中的每个内部点关于1-ring点的广义重心坐标。 记:\(\lambda_i=[\dots,0,\dots,\lambda_{j_a},\dots,\lambda_{j_b},\dots,0,\dots]\),即只是\(x_i\)的邻接点对应序号中的元素为第三步计算出的重心坐标值,其余位置为0

  3. 固定边界点对应的参数(如参数化到矩形域,则将网格曲面上边界点参数有序定义为矩形域的边界) 假定后\(N\)个点为边界点,初始化\({(u_i,v_i)}_{i=n+1}^N\)

  4. \(\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));
  • 对非边界点循环 计算均值重心坐标

    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
    37
    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)\)参数

    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)

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