Paper 實現 - Implicit Fairing of Irregular Meshes using Diffusion and Curvature Flow

Desbrun, Mathieu & Meyer, Mark & Schröder, Peter & Barr, Alan. (2001). Implicit Fairing of Irregular Meshes Using Diffusion and Curvature Flow. SIGGRAPH. 99.

2. Implicit fairing

這一節介紹implicit fairing，一種用於網格平滑的擴散方程的隱式積分。證明幾點相對於explicit methods的優點。雖然本節僅限於使用擴散項的線性近似，但隱式光順將作為一種穩健而有效的數值方法貫穿全文，即使對於非線性運算元也是如此。我們從建立框架和定義符號開始。

2.1 標記以及定義

X表示網格；xi表示網格上的頂點，eij表示連線xi和xj的邊，N1(i)表示頂點的一鄰域。

在曲面光順（surface fairing）的文獻中，大多技術使用了受限能量最小。最常使用的公式為，表面的總曲率：

\[E(S) = \int_S\kappa_1^2 + \kappa_2^2dS \]

由於主曲率是非線性的。在實際運用中更加傾向於membrane或thin-plate函式。

\[E_{membrane}(X) = \frac{1}{2}\int_\Omega X_u^2 + X_v^2 dudv \\ E_{thin plate}(X) = \frac{1}{2}\int_\Omega X_{uu}^2 + 2X_{uv}^2 + X_{vv}^2 dudv \]

2.2 Diffusion equation for mesh fairing

正如我們指出的，在網格中衰減噪聲的一種常見方法是通過擴散過程：

\[\frac{\part X}{\part t} = \lambda L(X) \]

將上面的公式隨時間積分，一個小擾動將在其附近迅速分散，平滑高頻，而主要形狀只會輕微退化。拉普拉斯運算元可以用傘形運算元在每個頂點處線性逼近（被稱為umbrella operator，雨傘操作） ：

\[L(x_i) = \frac{1}{m}\sum_{j\in N_1(i)}x_j - x_i \]

其中m為鄰域頂點個數。

那麼網格的擴散方程可以表示為：

\[X^{n+1} = (I + \lambda dtL)X^n \]

對於傘運算元，穩定性準則要求λdt<1。如果時間步長不滿足這個條件，表面會出現漣漪，並最終在整個表面上產生越來越大的振盪。另一方面，如果滿足這個條件，隨著n的增加，初始網格會變得越來越平滑。

Laplacian用矩陣形式進行表示如下：

\[\begin{bmatrix} \Delta x_1 \\ \vdots \\ \Delta x_n \end{bmatrix} = L \times \begin{bmatrix} x_1 \\ \vdots \\ x_n \end{bmatrix} \]

其中：

\[L_{ij} = \left\{\begin{matrix} -\sum_{x_k \in N_1(x_i)} \frac{1}{m} = -1, i=j\\ \frac{1}{m}, x_j \in N_1(x_i) \\ 0, otherwise \end{matrix}\right. \]

對於顯示的Laplacian平滑，是根據當前的狀態預估出特定時間之後的狀態。

2.3 Time-shifted evaluation

2.2節中提供的方法是顯式的方法，也叫做前向euler方法。不足的地方是，當網格很大的時候，需要很多次迭代才能夠有個顯著的平滑效果。

implicit的方法，可以避免上述的情況，具體公式如下：

\[(I - \lambda dt L)X^{n+1} = X^n \]

然後，通過增加值λdt，可以安全地獲得後續平滑。但解決一個線性系統是要付出代價的。

2.4 Solving the sparse linear system

幸運的是矩陣\(A = (I-\lambda dt L)\)是稀疏矩陣，能夠高效的進行求解。我們可以用一個預處理的雙共軛梯度（PBCG）迭代求解該系統。可以使用eigen內建的迭代求解器進行求解，具體參見：https://blog.csdn.net/xuezhisdc/article/details/54634080

2.6 Filter improvement

可以使用：

\[(I + \lambda dt L^2)X^{n+1} = X^n \]

甚至可以是一階和二階的混搭（這個是根據理解得到的，不一定正確）：

\[(I - (\lambda + \mu) dt L + \lambda \mu dt L^2)X^{n+1} = X^n \]

3. Automatic anti-shrink fairing

純粹的diffusion方法，會引入收縮。此處提出了基於volume的方法控制收縮。

3.1 Volume Computation

具體實現可以參見：三角網格體積計算.

3.2 Exact volume preservation

在一次迭代步驟之後，網格將會得到一個新的體積\(V^n\)。我們想要將體積恢復到之前的狀態\(V^0\)，而減小平滑後縮放的問題。我們對每個頂點進行了很簡單的縮放的步驟，讓每個頂點都乘以一個係數\(\beta = (V^0/V^n)^{1/3}\)。也可以將不變數擴充套件為其他的變數，如表面積。

4. An accurate diffusion process

通過定義laplacian的新的近似方式，來實現擴充套件的過程。

4.1 Inadequacy of umbrella operator

雨傘操作基於的假設是：所有的網格的邊長是1，頂點周圍相鄰邊之間的角度是相等的。

對於不同頻率的表型，雨傘操作得到的結果可能是一致的。如上圖所示。

需要定義離散Laplacian更好的近似擴散。

4.2 Simulation of the 1D heat equation

非數學專業，理解不了公式的原理，直接給出公式吧：

\[(x_{uu})_i = {2 \over \delta + \Delta}({x_{i-1} - x_i \over \delta} + {x_{i+1} - x_i \over \Delta} ) \]

公式中括號內的值有點像二階導數求解。

4.3 Extension to 3D

雖然對於1D的形式並不很理解，但是將其擴充套件為3D的形式還是很容易給出來的，將邊的長度最為式子中的\(\delta\)值，就可以得到下式：

\[L(x_i) = {2 \over E}\sum_{j\in N_1(i)}{x_j - x_i \over |e_{ij}|}, \ with\ E = \sum_{j\in N_1(i)} |e_{ij}| \]

命個名叫：scale-dependent umbrella operator。

這種方式可以更好的對irregular mesh進行平滑。如下圖，圖c就是採用上面的方法實現的：

由於該實現是非線性的，需要額外的實現進行求解。參見原文。

5 Curvature flow for noise removal

使用曲率流替代diffusion。

5.1 Diffusion vs. curvature flow

曲率流通過以等於平均曲率κ的速度沿曲面法線n移動來平滑曲面：

\[{\partial x_i \over \partial t} = \bar{\kappa_i}\bold{n}_i \]

平均曲率計算可以參見：PMP SRC Algo - 微分幾何相關基礎.

5.2 Curvature normal calculation

一個好主意是檢查法向量的散度，因為它是平均曲率的定義（\(\bar{\kappa}=div\ \bold{n}\)）。如果一個頂點周圍的所有面法線都相同，則該頂點不應移動（零曲率）。考慮到這一點，我們選擇了下面曲率法線的微分幾何定義:

\[{\triangledown A \over 2 A} = \bar{\kappa}\bold{n} \]

式中，A是點P周圍需要曲率的小區域的面積，∇是相對於P的（x，y，z）座標的導數。根據這個定義，我們將得到平坦區域的零向量。如圖8所示，我們看到在平面上移動中心頂點席不會改變表面面積。另一方面，將其移動到平面上方或下方將始終增加區域性區域。因此，對於區域性平坦曲面，無論其價態、相鄰面的長寬比或頂點周圍的邊長，我們都具有所需的零面積梯度性質。

得到離散化的形式如下（具體怎麼來的真心看不懂呀，不過應該和cot Laplacian的推導比較類似，可以參見：DGP - 2. Discrete differential geometry）：

\[-\bar{\kappa}\bold n = \frac{1}{4A} \sum_{j\in N_1(i)}(\cot\alpha_j + \cot \beta_j)(x_j - x_i) \]

5.3 Boundaries

對於非閉合曲面或具有孔的曲面，我們可以為邊界上的頂點定義一種特殊處理方式。平均曲率的概念對這些頂點沒有意義。相反，我們希望平滑邊界，這樣隨著迭代的進行，孔本身的形狀會變得越來越圓。然後我們可以使用公式：

\[L(x_i) = {2 \over E}\sum_{j\in N_1(i)}{x_j - x_i \over |e_{ij}|}, \ with\ E = \sum_{j\in N_1(i)} |e_{ij}| \]

僅限於兩個相鄰的將平滑邊界曲線本身。

另一種可能的方法是建立一個虛擬頂點，儲存但不顯示，最初放置在閉合邊界上所有頂點的重心處。一組與該頂點相鄰並一個接一個地連線邊界頂點的面也是虛擬建立的。然後我們可以使用基本演算法，而不必對邊界進行任何特殊處理，因為現在每個頂點都有一個閉合區域。

具體如何實現？？？

5.4 Implementation

簡單而言上述的方法就是利用cot Laplacian替代之前的uniform Laplacian實現。

\[L(x_i) = \frac{1}{2A(x_i)}\cdot\frac{1}{\sum_{j\in \Omega(i)}(\cot \alpha_{ij} + \cot \beta_{ij})}\sum_{j\in \Omega(i)}(\cot \alpha_{ij} + \cot \beta_{ij})(f_j - f_i) \]

那麼對應的Laplacian矩陣中：

\[L_{ij} = \left\{\begin{matrix} -\frac{1}{2A(x_i)}, i=j\\ \frac{1}{2A(x_i)}\cdot\frac{(\cot \alpha_{ij} + \cot \beta_{ij})}{\sum_{j\in \Omega(i)}(\cot \alpha_{ij} + \cot \beta_{ij})}, x_j \in N_1(x_i) \\ 0, otherwise \end{matrix}\right. \]

原始碼實現

$ git clone --recursive https://github.com/pmp-library/pmp-template.git

然後替換MyViewer.h為：

#include <pmp/visualization/MeshViewer.h>

class SurfaceSmoothingTest;
class MyViewer : public pmp::MeshViewer
{
public:
    //! constructor
    MyViewer(const char* title, int width, int height);

protected:
    //! this function handles keyboard events
    void keyboard(int key, int code, int action, int mod) override;

    virtual void process_imgui();

private:
    std::shared_ptr<SurfaceSmoothingTest> spSmoothing_;
};

替換MyViewer.cpp為:

//=============================================================================
// Copyright (C) 2013-2019 The pmp-library developers
//
// This file is part of the Polygon Mesh Processing Library.
// Distributed under the terms of the MIT license, see LICENSE.txt for details.
//
// SPDX-License-Identifier: MIT
//=============================================================================

#include "MyViewer.h"

#include "pmp/SurfaceMesh.h"
#include "pmp/algorithms/SurfaceTriangulation.h"
#include "pmp/algorithms/SurfaceNormals.h"
#include "pmp/Types.h"

#include <Eigen/Dense>
#include <Eigen/Sparse>
#include "pmp/algorithms/DifferentialGeometry.h"

#include <imgui.h>

using namespace pmp;
class SurfaceSmoothingTest
{
public:
    SurfaceSmoothingTest(SurfaceMesh& mesh);

    ~SurfaceSmoothingTest();

    void explicit_smoothing(Scalar timestep, bool use_uniform_laplace, bool rescale=false);

    void implicit_smoothing(Scalar timestep, bool use_uniform_laplace, bool rescale=false);

private:
    void compute_weights();
private:
    SurfaceMesh& mesh_;

    unsigned int how_many_vertex_weights_;
    unsigned int how_many_edge_weights_;
    bool use_uniform_laplace_;
};

using SparseMatrix = Eigen::SparseMatrix<double>;
using Triplet = Eigen::Triplet<double>;

SurfaceSmoothingTest::SurfaceSmoothingTest(SurfaceMesh& mesh) : mesh_(mesh), use_uniform_laplace_(false)
{
    how_many_vertex_weights_ = 0;
    how_many_edge_weights_ = 0;
}

SurfaceSmoothingTest::~SurfaceSmoothingTest()
{
    auto eweight = mesh_.edge_property<Scalar>("e:laplace");
    auto vweight_laplace = mesh_.vertex_property<Scalar>("v:laplace");
    auto vweight_coeff = mesh_.vertex_property<Scalar>("v:coeff");
    if (eweight)
        mesh_.remove_edge_property(eweight);
    if (vweight_coeff)
        mesh_.remove_vertex_property(vweight_coeff);
    if (vweight_laplace)
        mesh_.remove_vertex_property(vweight_laplace);
}

void SurfaceSmoothingTest::compute_weights()
{
    auto eweight = mesh_.edge_property<Scalar>("e:laplace");
    auto vweight_laplace = mesh_.vertex_property<Scalar>("v:laplace");
    auto vweight_coeff = mesh_.vertex_property<Scalar>("v:coeff");

    if (use_uniform_laplace_)
    {
        for (auto e : mesh_.edges())
            eweight[e] = 1.0;
        for (auto v : mesh_.vertices())
        {
            vweight_laplace[v] = 1.0 / mesh_.valence(v);
            vweight_coeff[v] = 1.0;
        }
    }
    else
    {
        for (auto e : mesh_.edges())
            eweight[e] = std::max(0.0, cotan_weight(mesh_, e));

        Scalar value(0.0);
        for (auto v : mesh_.vertices())
        {
            value = Scalar(0.0);
            for (auto he : mesh_.halfedges(v))
            {
                value += eweight[mesh_.edge(he)];
            }
            vweight_laplace[v] = value == 0.0 ? 0.0 : Scalar(1.0) / value;
            vweight_coeff[v] = 0.5 / voronoi_area(mesh_, v);
        }
    }

    how_many_edge_weights_ = mesh_.n_edges();
    how_many_vertex_weights_ = mesh_.n_vertices();
}

void SurfaceSmoothingTest::implicit_smoothing(Scalar timestep, bool use_uniform_laplace, bool rescale)
{
    use_uniform_laplace_ = use_uniform_laplace;
    if (!mesh_.n_vertices())
        return;

    compute_weights();

    // properties
    auto points = mesh_.get_vertex_property<Point>("v:point");
    auto eweight = mesh_.edge_property<Scalar>("e:laplace");
    auto vweight_laplace = mesh_.vertex_property<Scalar>("v:laplace");
    auto vweight_coeff = mesh_.vertex_property<Scalar>("v:coeff");
    auto idx = mesh_.add_vertex_property<int>("v:idx", -1);

    // collect free (non-boundary) vertices in array free_vertices[]
    // assign indices such that idx[ free_vertices[i] ] == i
    unsigned i = 0;
    std::vector<Vertex> free_vertices;
    free_vertices.reserve(mesh_.n_vertices());
    for (auto v : mesh_.vertices())
    {
        if (!mesh_.is_boundary(v))
        {
            idx[v] = i++;
            free_vertices.push_back(v);
        }
    }
    const unsigned int n = free_vertices.size();

    // A * X = B
    // (I - \lambda dt L) * x(n+1) = X(n)
    SparseMatrix A(n, n);
    Eigen::MatrixXd B(n, 3);

    // nonzero elements of A as triplets: (row, column, value)
    std::vector<Triplet> triplets;

    // setup matrix A and rhs B
    dvec3 b;
    double ww;
    Vertex v, vv;
    Edge e;
    for (unsigned int i = 0; i < n; ++i)
    {
        v = free_vertices[i];

        // rhs row
        b = static_cast<dvec3>(points[v]);
        B.row(i) = (Eigen::Vector3d)b;

        // lhs row
        for (auto h : mesh_.halfedges(v))
        {
            vv = mesh_.to_vertex(h);
            e = mesh_.edge(h);

            // fixed boundary vertex -> right hand side
            if (!mesh_.is_boundary(vv))
            {
                triplets.emplace_back(i, idx[vv],  -timestep * vweight_coeff[v] * eweight[e] * vweight_laplace[v]);
            }
        }

        // center vertex -> matrix
        triplets.emplace_back(i, i, 1 + timestep * vweight_coeff[v]);
    }

    // build sparse matrix from triplets
    A.setFromTriplets(triplets.begin(), triplets.end());

    // solve A * X = B
    Eigen::LeastSquaresConjugateGradient<Eigen::SparseMatrix<double> > solver;
    solver.setTolerance(0.001f);
    solver.compute(A);

    Eigen::MatrixXd X = solver.solve(B);
    Scalar volume_before = mesh_volume(mesh_);

    if (solver.info() != Eigen::Success)
    {
        std::cerr << "SurfaceSmoothing: Could not solve linear system\n";
        return;
    }
    else
    {
        // copy solution
        for (unsigned int i = 0; i < n; ++i)
        {
            points[free_vertices[i]][0] = (float)X.row(i)[0];
            points[free_vertices[i]][1] = (float)X.row(i)[1];
            points[free_vertices[i]][2] = (float)X.row(i)[2];
        }
    }

    if (rescale)
    {
        // restore original surface area
        Scalar volume_after = mesh_volume(mesh_);
        Scalar scale = std::pow(volume_before / volume_after, 1.0 / 3.0);
        for (auto v : mesh_.vertices())
            mesh_.position(v) *= scale;
    }

    // clean-up
    mesh_.remove_vertex_property(idx);
}

void SurfaceSmoothingTest::explicit_smoothing(Scalar timestep, bool use_uniform_laplace, bool rescale) 
{
    use_uniform_laplace_ = use_uniform_laplace;
    timestep;
    if (!mesh_.n_vertices())
        return;

    compute_weights();

    // properties
    auto points = mesh_.get_vertex_property<Point>("v:point");
    auto eweight = mesh_.edge_property<Scalar>("e:laplace");
    auto vweight_laplace = mesh_.vertex_property<Scalar>("v:laplace");
    auto vweight_coeff = mesh_.vertex_property<Scalar>("v:coeff");
    auto idx = mesh_.add_vertex_property<int>("v:idx", -1);

    // collect free (non-boundary) vertices in array free_vertices[]
    // assign indices such that idx[ free_vertices[i] ] == i
    unsigned i = 0;
    std::vector<Vertex> free_vertices;
    free_vertices.reserve(mesh_.n_vertices());
    for (auto v : mesh_.vertices())
    {
        if (!mesh_.is_boundary(v))
        {
            idx[v] = i++;
            free_vertices.push_back(v);
        }
    }
    const unsigned int n = free_vertices.size();

    // X = A * B
    // x(n+1) = (I + \lambda dt L) * X(n)
    SparseMatrix A(n, n);
    Eigen::MatrixXd B(n, 3);

    // nonzero elements of A as triplets: (row, column, value)
    std::vector<Triplet> triplets;

    // setup matrix A and rhs B
    dvec3 b;
    Vertex v, vv;
    for (unsigned int i = 0; i < n; ++i)
    {
        v = free_vertices[i];

        // rhs row
        b = static_cast<dvec3>(points[v]);
        B.row(i) = (Eigen::Vector3d)b;

        // lhs row
        for (auto h : mesh_.halfedges(v))
        {
            vv = mesh_.to_vertex(h);

            // fixed boundary vertex -> right hand side
            if (!mesh_.is_boundary(vv))
            {
                triplets.emplace_back(i, idx[vv], timestep * vweight_laplace[v] * eweight[mesh_.edge(h)] * vweight_coeff[v]);
            }
        }

        // center vertex -> matrix
        triplets.emplace_back(i, i, 1 - timestep * vweight_coeff[v]);
    }

    // build sparse matrix from triplets
    A.setFromTriplets(triplets.begin(), triplets.end());

    // solve X = A * B
    Eigen::MatrixXd X = A * B;
    Scalar volume_before = mesh_volume(mesh_);
    // copy solution
    for (unsigned int i = 0; i < n; ++i)
    {
        points[free_vertices[i]] = X.row(i);
    }

    if (rescale)
    {
        // restore original surface area
        Scalar volume_after = mesh_volume(mesh_);
        Scalar scale = std::pow(volume_before / volume_after, 1.0/3.0);
        for (auto v : mesh_.vertices())
            mesh_.position(v) *= scale;
    }

    // clean-up
    mesh_.remove_vertex_property(idx);
}

//////////////////////////////////////////////////////////////////

MyViewer::MyViewer(const char* title, int width, int height)
    : MeshViewer(title, width, height)
{
    set_draw_mode("Smooth Shading");
    spSmoothing_ = std::make_shared<SurfaceSmoothingTest>(mesh_);
}

void MyViewer::keyboard(int key, int scancode, int action, int mods)
{
    if (action != GLFW_PRESS) // only react on key press events
        return;

    switch (key)
    {
        // add your own keyboard action here
        default:
        {
            MeshViewer::keyboard(key, scancode, action, mods);
            break;
        }
    }
}

void MyViewer::process_imgui()
{
    MeshViewer::process_imgui();

    ImGui::Spacing();
    ImGui::Spacing();

    if (ImGui::CollapsingHeader("Smoothing", ImGuiTreeNodeFlags_DefaultOpen))
    {
        static float timestep = 0.001;
        ImGui::PushItemWidth(100);
        ImGui::SliderFloat("TimeStep", &timestep, 0.0001, 10.0);
        ImGui::PopItemWidth();

        static bool brescale = false;
        ImGui::Checkbox("Enable Rescale", &brescale);

        static bool buse_uniform = false;
        ImGui::Checkbox("Use Unifrom", &buse_uniform);

        ImGui::Spacing();

        if (ImGui::Button("Explicit Smoothing"))
        {
            spSmoothing_->explicit_smoothing(timestep, buse_uniform, brescale);
            update_mesh();
        }

        ImGui::Spacing();

        if (ImGui::Button("Implicit Smoothing"))
        {
            spSmoothing_->implicit_smoothing(timestep, buse_uniform, brescale);
            update_mesh();
        }

        if (ImGui::Button("Reload mesh"))
        {
            load_mesh(filename_.c_str());
            update_mesh();
        }
    }
}

//=============================================================================

然後利用cmake進行編譯構建即可。

示例結果如下：