基于三角面法向的双边滤波网格平滑算法实现

欢迎关注更多精彩
关注我，学习常用算法与数据结构，一题多解，降维打击。

本文实现的算法来自下面的论文
参考文献Bilateral Normal Filtering for Mesh Denoising ，https://ieeexplore.ieee.org/document/5674028

算法原理及过程

算法原理

先利用双边滤波对面法向进行调整，然后根据法向去调整顶点，使得面的法向尽量接近调整后的法向。

法向调整

$n^{k+1}(f_i)=\frac{1}{K_p}\sum _{j\in \Omega (i)}{A}_j{W}_s(||c_i-c_j||)W_r(||n^k(f_i)-n^k(f_j)||)\cdot n^k(f_j)$

上述式子中

k代表迭代次数

$n^k(f_i)表示经过k次迭代后第i个面的法向$

$c_i代表第i个面的中心$

$\Omega (i)表示第i个面的所有邻面$

$A_j表示第j个面的面积$

$$\begin{gathered} W_s(x)是一个高斯函数=e^{-\frac{x^2}{2{\sigma }_s^2}}, \\ 其中x=c_i-c_j的长度，{\sigma }_s是可调参数， \\ 本次实现中取值为所有邻面中心差值的平均值 \end{gathered}$$

$$\begin{gathered} W_r(x)是一个高斯函数=e^{-\frac{x^2}{2{\sigma }_r^2}}, \\ 其中x=n^k(f_i)-n^k(f_j)的长度，{\sigma }_r是可调参数， \\ 本次实现中取值为1 \end{gathered}$$

$K_p=\sum _{j\in \Omega (i)}{A}_j{W}_s(||c_i-c_j||)W_r(||n^k(f_i)-n^k(f_j)||)$

法向调整可以运用多次，具体参数可调k。
每次计算完法向后，要对法向进行归一化。

点顶恢复

想要根据法向计算出所有顶点是非常困难的。这里采用单步高斯迭代的方式，每次只调整一个点顶，使得法向尽可能接近计算出来的法向。

根据定义我们有如下公式

$\{\begin{array}{l}{n}_f\cdot (x_j-x_i)=0\\ n_f\cdot (x_k-x_j)=0\\ n_f\cdot (x_i-x_k)=0\end{array}$

$x_i,x_j,x_k是某个面的三点，n_f是调整后的法向，上述公式表示每条边都要垂直于法线$

$我们的目标是使上述公式越接近于0越好$

可以得到一个能量方程

$E=\sum _{f_k}\sum _{i,j\in f_k}(n_k\cdot (x_j-x_i){)}^2$

现在要做的就是求一组顶点，使得E最小。

这里每次只移动一个点，使得E最小。

对上述方程进行求导并使其等于0，可以推导出某一点的计算公式

$x_i^{new}=x_i+\frac{1}{N_i}\sum _{f_j\in \Omega (i)}{n}_j\cdot (n_j\cdot (c_j-x_i))$

$N_i表示顶点的度$

$\Omega (i)是第i个点所有邻面$

$x_i是原来要变化的点$

$c_j是f_i中心点的坐标$

$n_j是第j个面的目标法向$

上述过程每次对所有顶点进行恢复，可以进行多次参数可调，本次实现采用8次。

算法实现

代码github https://github.com/LightningBilly/ACMAlgorithms/tree/master/%E5%9B%BE%E5%BD%A2%E5%AD%A6%E7%AE%97%E6%B3%95/%E4%B8%89%E8%A7%92%E7%BD%91%E6%A0%BC%E7%AE%97%E6%B3%95/MeshWork/src/hw3

#include <iostream>
#include <fstream>
#include <sstream>
#include "../../PolyMesh/IOManager.h"
#include <string>
#include <unistd.h>

using namespace std;
using namespace acamcad;
using namespace polymesh;

void hw3()
{
    puts("hw3");

    char buffer[500];
    getcwd(buffer, 500);
    printf("The current directory is: %s/../\n", buffer);
    string mesh_path = buffer;
    mesh_path += "/../src/hw3/bunny_random.obj";
    double SigmaCenter, SigmaNormal=1;
    double NormalIterNum=8, VertexIterNum=8;

    PolyMesh* mesh = new PolyMesh();
    loadMesh(mesh_path, mesh);

    mesh->updateMeshNormal();
    std::vector<MVector3> NewNormal(mesh->numPolygons());
    std::vector<double> FaceArea(mesh->numPolygons());
    std::vector<MPoint3> FaceCenter(mesh->numPolygons());

    for (MPolyFace* fh : mesh->polyfaces())
    {
        int f_id = (*fh).index();
        NewNormal[f_id] = (*fh).normal();
        std::vector<MVert*> P;
        for (FaceVertexIter vv_it = mesh->fv_iter(fh); vv_it.isValid(); ++vv_it)
        {
            P.push_back(*vv_it);
        }
        auto e12 = P[1]->position() - P[0]->position();
        auto e13 = P[2]->position() - P[0]->position();
        double area = cross(e12, e13).norm() * 0.5;
        FaceArea[f_id] = area;
        FaceCenter[f_id] = mesh->calculatFaceCenter(fh);
    }

    /*
     * 计算点标准差参数
     * 把所有邻面的中点作差
     * 最后除以面的3倍
     */
    SigmaCenter = 0;
    for (MPolyFace* fh : mesh->polyfaces())
    {
        int f_id = (*fh).index();
        for (FaceFaceIter nei_fh = mesh->ff_iter(fh); nei_fh.isValid(); ++nei_fh)
        {
            int ff_id = (*nei_fh)->index();
            SigmaCenter += (FaceCenter[f_id] - FaceCenter[ff_id]).norm();
        }
    }
    SigmaCenter /= mesh->numPolygons() * 3;

    /*
     * 利用法向迭代公式计算法向
     */
    for (int i = 0; i < NormalIterNum; i++)
    {
        for (MPolyFace* fh : mesh->polyfaces())
        {
            double Kp = 0;
            MVector3 NewN(0, 0, 0);
            int fh_id = (*fh).index();
            for (FaceFaceIter nei_fh = mesh->ff_iter(fh); nei_fh.isValid(); ++nei_fh)
            {
                int nei_fh_id = (*nei_fh)->index();
                double delta_center = (FaceCenter[fh_id] - FaceCenter[nei_fh_id]).norm();
                double delta_normal = (NewNormal[fh_id] - NewNormal[nei_fh_id]).norm();
                double Aj = FaceArea[nei_fh_id];
                double Ws = exp(-delta_center * delta_center / (2 * SigmaCenter * SigmaCenter));
                double Wr = exp(-delta_normal * delta_normal / (2 * SigmaNormal * SigmaNormal));
                NewN += Aj * Ws * Wr * NewNormal[nei_fh_id];
                Kp += Aj * Ws * Wr;
            }
            NewNormal[fh_id] = NewN / Kp;
            NewNormal[fh_id] /= NewNormal[fh_id].norm();
        }
    }

    // 调整点位置
    for (int i = 0; i < VertexIterNum; i++)
    {
        for (MVert* vh : mesh->vertices())
        {
            MPoint3 x_i = (*vh).position();
            MPoint3 delta_xi(0, 0, 0);
            int Nei_count = 0;
            for (VertexFaceIter fh = mesh->vf_iter(vh); fh.isValid(); ++fh)
            {
                Nei_count++;
                MPoint3 cj = mesh->calculatFaceCenter(*fh);
                MVector3 nj = NewNormal[(*fh)->index()];
                delta_xi = delta_xi + nj * (nj.data()[0] * (cj - x_i).data()[0] + nj.data()[1] * (cj - x_i).data()[1] + nj.data()[2] * (cj - x_i).data()[2]);
            }
            x_i = x_i + delta_xi / Nei_count;
            (*vh).setPosition(x_i);
        }
    }

    // 输出网格数据
    writeMesh("result.obj", mesh);
}

算法效果

原型网格 https://github.com/LightningBilly/ACMAlgorithms/blob/master/%E5%9B%BE%E5%BD%A2%E5%AD%A6%E7%AE%97%E6%B3%95/%E4%B8%89%E8%A7%92%E7%BD%91%E6%A0%BC%E7%AE%97%E6%B3%95/MeshWork/src/hw3/bunny_random.obj

平滑后的网格
https://github.com/LightningBilly/ACMAlgorithms/blob/master/%E5%9B%BE%E5%BD%A2%E5%AD%A6%E7%AE%97%E6%B3%95/%E4%B8%89%E8%A7%92%E7%BD%91%E6%A0%BC%E7%AE%97%E6%B3%95/MeshWork/src/hw3/result.obj

原模型

平滑后模型

从上图中可以看到有很明显平滑效果。

---

本人码农，希望通过自己的分享，让大家更容易学懂计算机知识。