欢迎关注更多精彩
关注我,学习常用算法与数据结构,一题多解,降维打击。
本文实现的算法来自下面的论文
参考文献Bilateral Normal Filtering for Mesh Denoising ,https://ieeexplore.ieee.org/document/5674028
算法原理及过程
算法原理
先利用双边滤波对面法向进行调整,然后根据法向去调整顶点,使得面的法向尽量接近调整后的法向。
法向调整
n k + 1 ( f i ) = 1 K p ∑ j ∈ Ω ( i ) A j W s ( ∣ ∣ c i − c j ∣ ∣ ) W r ( ∣ ∣ n k ( f i ) − n k ( f j ) ∣ ∣ ) ⋅ n k ( f j ) n^{k+1}(f_i) = \frac{1}{K_p} \displaystyle \sum_{j\in Ω(i)}A_{j}W_s(||c_i-c_j||)W_r(||n^k(f_i)-n^k(f_j)||)\cdot n^k(f_j) nk+1(fi)=Kp1j∈Ω(i)∑AjWs(∣∣ci−cj∣∣)Wr(∣∣nk(fi)−nk(fj)∣∣)⋅nk(fj)
上述式子中
k代表迭代次数
n k ( f i ) 表 示 经 过 k 次 迭 代 后 第 i 个 面 的 法 向 n^k(f_i) 表示经过k次迭代后第i个面的法向 nk(fi)表示经过k次迭代后第i个面的法向
c i 代 表 第 i 个 面 的 中 心 c_i代表第i个面的中心 ci代表第i个面的中心
Ω ( i ) 表 示 第 i 个 面 的 所 有 邻 面 Ω(i)表示第i个面的所有邻面 Ω(i)表示第i个面的所有邻面
A j 表 示 第 j 个 面 的 面 积 A_j表示第j个面的面积 Aj表示第j个面的面积
W s ( x ) 是 一 个 高 斯 函 数 = e − x 2 2 σ s 2 , 其 中 x = c i − c j 的 长 度 , σ s 是 可 调 参 数 , 本 次 实 现 中 取 值 为 所 有 邻 面 中 心 差 值 的 平 均 值 W_s(x)是一个高斯函数 = e^{-\frac{x^2}{2\sigma^2_s}}, \\其中x=c_i-c_j的长度,\sigma_s是可调参数,\\本次实现中取值为所有邻面中心差值的平均值 Ws(x)是一个高斯函数=e−2σs2x2,其中x=ci−cj的长度,σs是可调参数,本次实现中取值为所有邻面中心差值的平均值
W r ( x ) 是 一 个 高 斯 函 数 = e − x 2 2 σ r 2 , 其 中 x = n k ( f i ) − n k ( f j ) 的 长 度 , σ r 是 可 调 参 数 , 本 次 实 现 中 取 值 为 1 W_r(x)是一个高斯函数 = e^{-\frac{x^2}{2\sigma^2_r}}, \\其中x=n^k(f_i)-n^k(f_j)的长度,\sigma_r是可调参数,\\本次实现中取值为1 Wr(x)是一个高斯函数=e−2σr2x2,其中x=nk(fi)−nk(fj)的长度,σr是可调参数,本次实现中取值为1
K p = ∑ j ∈ Ω ( i ) A j W s ( ∣ ∣ c i − c j ∣ ∣ ) W r ( ∣ ∣ n k ( f i ) − n k ( f j ) ∣ ∣ ) K_p = \displaystyle \sum_{j\in Ω(i)}A_{j}W_s(||c_i-c_j||)W_r(||n^k(f_i)-n^k(f_j)||) Kp=j∈Ω(i)∑AjWs(∣∣ci−cj∣∣)Wr(∣∣nk(fi)−nk(fj)∣∣)
法向调整可以运用多次,具体参数可调k。
每次计算完法向后,要对法向进行归一化。
点顶恢复
想要根据法向计算出所有顶点是非常困难的。这里采用单步高斯迭代的方式,每次只调整一个点顶,使得法向尽可能接近计算出来的法向。
根据定义我们有如下公式
{ n f ⋅ ( x j − x i ) = 0 n f ⋅ ( x k − x j ) = 0 n f ⋅ ( x i − x k ) = 0 \left\{\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}\right. ⎩⎨⎧nf⋅(xj−xi)=0nf⋅(xk−xj)=0nf⋅(xi−xk)=0
x i , x j , x k 是 某 个 面 的 三 点 , n f 是 调 整 后 的 法 向 , 上 述 公 式 表 示 每 条 边 都 要 垂 直 于 法 线 x_i, x_j, x_k 是某个面的三点,n_f是调整后的法向,上述公式表示每条边都要垂直于法线 xi,xj,xk是某个面的三点,nf是调整后的法向,上述公式表示每条边都要垂直于法线
我 们 的 目 标 是 使 上 述 公 式 越 接 近 于 0 越 好 我们的目标是使上述公式越接近于0越好 我们的目标是使上述公式越接近于0越好
可以得到一个能量方程
E = ∑ f k ∑ i , j ∈ f k ( n k ⋅ ( x j − x i ) ) 2 E=\displaystyle \sum_{f_k}\displaystyle \sum_{i,j \in f_k}(n_k\cdot (x_j-x_i))^2 E=fk∑i,j∈fk∑(nk⋅(xj−xi))2
现在要做的就是求一组顶点,使得E最小。
这里每次只移动一个点,使得E最小。
对上述方程进行求导并使其等于0,可以推导出某一点的计算公式
x i n e w = x i + 1 N i ∑ f j ∈ Ω ( i ) n j ⋅ ( n j ⋅ ( c j − x i ) ) x^{new}_i=x_i+\frac {1}{N_i} \displaystyle \sum_{f_j \in Ω(i)}n_j \cdot (n_j \cdot (c_j-x_i)) xinew=xi+Ni1fj∈Ω(i)∑nj⋅(nj⋅(cj−xi))
N i 表 示 顶 点 的 度 N_i 表示顶点的度 Ni表示顶点的度
Ω ( i ) 是 第 i 个 点 所 有 邻 面 Ω(i)是第i个点所有邻面 Ω(i)是第i个点所有邻面
x i 是 原 来 要 变 化 的 点 x_i是原来要变化的点 xi是原来要变化的点
c j 是 f i 中 心 点 的 坐 标 c_j 是f_i中心点的坐标 cj是fi中心点的坐标
n j 是 第 j 个 面 的 目 标 法 向 n_j 是第j个面的目标法向 nj是第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
原模型
平滑后模型
从上图中可以看到有很明显平滑效果。
本人码农,希望通过自己的分享,让大家更容易学懂计算机知识。