几何建模与处理-极小曲面全局方法、曲面参数化

sumsmile
1,555 阅读4分钟

欢迎关注公众号:sumsmile /专注图形渲染、Android开发

内容主要参考中科大刘利刚老师的课程《几何建模与处理基础》, 代码参考谢聪同学的作业。

极小曲面和曲面参数化从结果看，是两种不同的应用场景，原理有相同的部分，故放在一起。

基于C++ Clion开发, 代码:

github.com/summer-go/G…

工程依赖pmp(Polygon Mesh Processing Library)，pmp用法比较简单，不用刻意学，用到哪儿点进去看下就明白了

运行效果

all


模型


极小曲面


参数化后贴上纹理

极小曲面全局法

全局法是相对局部法而言，局部法是每次迭代一个点，逐步改变mesh的形状，全局法是一次生成极小曲面，相比局部法更快、效果更好(前提是mesh顶点数不大)

假设所有的点最后收敛，可得：

from-刘利刚《几何建模与处理》

感觉上图中公式viv_i项应该是vi∑j∈N(i)ωij{v_i\sum_{j \in N(i)} \omega_{ij}}


即L(vi)=vi∑j∈N(i)ωij−∑j∈N(i)ωijvj=0L(v_i) = {v_i\sum_{j \in N(i)} \omega_{ij} - \sum_{j \in N(i)} \omega_{ij}v_j} = 0

注意，极小曲面需要固定住边界点

即边界处的点不能移动，为已知数，将已知数移到右边，可联立方程：

（I−▽cot−weight）∗Pnew→=B边界→∗W（I - \bigtriangledown_{cot-weight}）* \overrightarrow{P_{new}} = \overrightarrow{B_{边界}} * W

这个算式不太正确，但是方便我自己理解

可结合代码和下面推导过程理解：

推导示意

代码处理分5步，参考下面的注释

// 方程组求解极小曲面 uniform_weight = false
void HarmonicSurface::harmonic_3d(bool uniform_weight) {
    std::vector<pmp::Vertex> inner_v;

    // 1. 整理出非边界的点
    for (auto v : mesh_.vertices())
    {
        if (!mesh_.is_boundary(v))
        {
            v_idx[v] = inner_v.size();
            inner_v.push_back(v);
        }
    }

    Eigen::MatrixX3d B = Eigen::MatrixX3d::Zero(inner_v.size(), 3);
    std::vector<Eigen::Triplet<double>> L_triplets;
    for (auto v : mesh_.vertices())
    {
        if (mesh_.is_boundary(v))
            continue;
        double w = 0;
        double ww = 0;
        for (auto h : mesh_.halfedges(v))
        {
            w = uniform_weight ? 1: e_cotan[mesh_.edge(h)];
            ww += w;
            auto vv = mesh_.to_vertex(h);
            if (v_idx[vv] == -1)
                // 2. 将边界点的权重累加到右边
                B.row(v_idx[v]) -= w * Eigen::Vector3d(mesh_.position(vv));
            else
                // 3. 非边界点的权重放到方程矩阵中(左边)
                L_triplets.push_back({v_idx[v], v_idx[vv], w});
        }
        L_triplets.push_back({v_idx[v], v_idx[v], -ww});
    }

    Eigen::SparseMatrix<double> L(inner_v.size(), inner_v.size());
    L.setFromTriplets(L_triplets.begin(), L_triplets.end());

    Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> solver(L);
    // 4. 解方程得到非边界点的新坐标
    Eigen::MatrixXd X = solver.solve(B);
    if (solver.info() != Eigen::Success)
    {
        std::cerr << "SurfaceSmoothing: Could not solve linear system\n";
    }
    else
    {
        // 5. 更新非边界点的坐标
        for (auto v : inner_v)
            mesh_.position(v) = X.row(v_idx[v]);
    }

    // 先不管下面的逻辑
    for(auto v: mesh_.vertices())
    {
        v_tex[v] = pmp::TexCoord(mesh_.position(v)[0], mesh_.position(v)[1]);
        v_tex[v] = (v_tex[v]/2)+pmp::TexCoord(0.5, 0.5);
    }

    mesh_.remove_vertex_property(v_idx);
}

得到结果：

Nefertiti_face 极小曲面

Balls 极小曲面

曲面参数化及纹理映射

曲面参数化使用的一个场景是，将3D的曲面按一定的规则展开成二维平面，比如纹理贴图。

from-刘利刚《几何建模与处理》

from-刘利刚《几何建模与处理》

这一段就来看看这个规则是什么。把这段弄清楚可加深理解"渲染里纹理贴图的逻辑"。

from-刘利刚《几何建模与处理》

曲面参数化

曲面参数化的计算和极小曲面计算的原理是一致的，都是假设曲面平滑收敛了。

即 Laplace-Beltrami算子的结果为0.

参数化考虑规整性，将原曲面参数化到矩形或者圆，更适合后续的处理。

  • 单位圆参数化

有些曲面接近圆，可考虑将纹理做成圆贴图来采样

单位圆参数化计算，将边界点按极坐标逆时针(或顺时针,结合坐标系考虑)排序，映射函数如下:

void HarmonicSurface::mapping_tex_boundary_to_circle(){
    std::vector<std::tuple<pmp::Vertex, double>> b_v;

    // 1. 将边界按连接顺序，添加到vector中，按边长计算权重w
    sort_boundary(b_v);
    for(auto& [v, w]: b_v)
        // 2. 按权重w取占2π的比例
        v_tex[v] = pmp::TexCoord (std::cos(w*2*M_PI), std::sin(w*2*M_PI));
}

void HarmonicSurface::sort_boundary(std::vector<std::tuple<pmp::Vertex, double>>& b_v) {
    pmp::Halfedge b_h;
    for(auto h: mesh_.halfedges())
        if(mesh_.is_boundary(h)) {
            b_h = h;
            break;
        }
    double boundary_len = mesh_.edge_length(mesh_.edge(b_h));
    b_v.push_back({mesh_.from_vertex(b_h), boundary_len});
    // 注意，按照mesh的生成规则，边界的half_edge按照next可依次遍历成一个环，但非边界的half_edge按next递归只能遍历单个三角形

    for(auto h = mesh_.next_halfedge(b_h); h != b_h; h = mesh_.next_halfedge(h)) {
        boundary_len += mesh_.edge_length(mesh_.edge(h));
        b_v.push_back({mesh_.from_vertex(h), boundary_len});
    }
    for(auto& [v, w]: b_v)
        w /= boundary_len;
}

注意，按照mesh的生成规则，边界的half_edge按照next可依次遍历成一个环，但非边界的half_edge按next递归只能遍历单个三角形，参考：Half-Edge Data Structures

  • 矩形参数化

实现方案不固定，原理类似，也是参照边长权重，规则的分布到矩形的四条边上。

比较好理解的方案 from 王宸

实际的代码有点差异，没有将点均匀分配到4条边上

结合图片和代码理解，注释写的比较详细了 

// ratio 用来处理正方形、长方形的逻辑
void HarmonicSurface::mapping_tex_boundary_to_rectangle(double ratio){
    std::vector<std::tuple<pmp::Vertex, double>> b_v;
    sort_boundary(b_v);
    // 1. 第一个点放到(1, 0) 处
    v_tex[std::get<0>(b_v[0])] = pmp::TexCoord(1, 0);
    for(size_t i = 1; i < b_v.size(); i++)
    {
        auto& v = std::get<0>(b_v[i]);
        auto& w = std::get<1>(b_v[i]);
        auto& w1 = std::get<1>(b_v[i-1]);

        // 2. 映射到矩形的四个顶点，先固定好四个顶点，对比前一个点和当前点处于的位置
        if(w1 < 1.0/8 && w > 1.0/8)
            v_tex[v] = pmp::TexCoord(1, ratio*1);
        else if(w1 < 3.0/8 && w > 3.0/8)
            v_tex[v] = pmp::TexCoord(-1, ratio*1);
        else if(w1 < 5.0/8 && w > 5.0/8)
            v_tex[v] = pmp::TexCoord(-1, ratio*-1);
        else if(w1 < 7.0/8 && w > 7.0/8)
            v_tex[v] = pmp::TexCoord(1, ratio*-1);

        // 3. 映射到矩形的四条边上
        // 右边
        else if(w < 1.0/8 || w > 7.0/8)
            v_tex[v] = pmp::TexCoord(1, ratio*std::sin(w*2*M_PI));

        // 左边
        else if (w >= 3.0/8 && w < 5.0/8)
            v_tex[v] = pmp::TexCoord(-1, ratio*std::sin(w*2*M_PI));

        // 上边
        else if (w >= 1.0/8 && w < 3.0/8)
            v_tex[v] = pmp::TexCoord(std::cos(w*2*M_PI), ratio*1);

        // 下边
        else if (w >= 5.0/8 && w < 7.0/8)
            v_tex[v] = pmp::TexCoord(std::cos(w*2*M_PI), -ratio*1);
    }
}

参数化完成了，即曲面被展开到二维平面上，下面看如何计算每个点对应的纹理坐标

注意，原曲面点不做更改，生成的展开面需要另外保存，用来计算和纹理的映射关系

纹理映射

纹理映射

上面这张图很清晰的说明了纹理映射的逻辑

  1. 曲面展开到平面上，和极小曲面的计算类似。区别在于最后算出来的mesh结果用来映射纹理坐标。
  2. 计算平面上点的坐标，即对应纹理的采样坐标(注意坐标[-1, 1] 转换到纹理域[0, 1])
  3. opengl渲染(熟悉opengl的同学，很容易理解该逻辑)
// 1. 计算纹理坐标
Eigen::MatrixXd X = solver.solve(B);
if (solver.info() != Eigen::Success)
{
    std::cerr << "SurfaceSmoothing: Could not solve linear system\n";
}
else
{
    // copy solution
    for (auto v : inner_v)
        // 纹理映射和极小曲面处理不同
        // 实现纹理映射并不会真的改变顶点的坐标，只绑定每个点对应的纹理的坐标
        v_tex[v] = X.row(v_idx[v]);
}
for(auto v: mesh_.vertices())
    // 注意 点的坐标取值范围是[-1, 1]，约定俗成方便计算
    // 纹理坐标取值范围是[0, 1]
    // 将[-1, 1] 映射到 [0, 1]
    v_tex[v] = (v_tex[v]/2)+pmp::TexCoord(0.5, 0.5);

mesh_.remove_vertex_property(v_idx);

调度流程

if (ImGui::CollapsingHeader("Harmonic", ImGuiTreeNodeFlags_DefaultOpen))
{
    static int boundary_shape = 0;

    ImGui::RadioButton("none", &boundary_shape, 0);
    ImGui::SameLine();
    ImGui::RadioButton("circle", &boundary_shape, 1);
    ImGui::SameLine();
    ImGui::RadioButton("rect(2:1)", &boundary_shape, 2);
    ImGui::SameLine();
    ImGui::RadioButton("square", &boundary_shape, 3);

    if (ImGui::Button("Harmonic Tex"))
    {
        pmp_pupa::HarmonicSurface harmonicSurface(mesh_);

        // 映射到不同的平面形状
        if (boundary_shape == 1) // circle
            harmonicSurface.mapping_tex_boundary_to_circle();
        else if (boundary_shape == 2) // rect(2:1)
            harmonicSurface.mapping_tex_boundary_to_rectangle(0.5);
        else if (boundary_shape == 3) // square 正方形
            harmonicSurface.mapping_tex_boundary_to_rectangle(1.0);

        // 求参数化(极小曲面算法一样)
        harmonicSurface.harmonic_tex();

        // 更新mesh(mesh没有变更，这一步无效)
        update_mesh();
        // 将模型视图移到window正中间
        view_all();
        // 纹理使用pmp自带的checkerboard 纹理,即黑白格子
        mesh_.use_checkerboard_texture();
        // 设置绘制模式为Texture
        set_draw_mode("Texture");
    }

其他参考： Polygon Mesh Processing Library

参考资料

[1] Half-Edge Data Structures : 

jerryyin.info/geometry-pr…

[2] Polygon Mesh Processing Library: 

www.pmp-library.org/installatio…

欢迎关注公众号:sumsmile /专注图形渲染、Android开发

avatar
文章
阅读
粉丝