VEX in Houdini Laplacian and Taubin Smoothing

VEX in HoudiniLaplacian and Taubin Smoothing

Laplacian算子

Laplacian
算子是n维欧几里德空间中的一个二阶微分算子，定义为梯度grad（）的散度div（）。

定义

如果f是二阶可微的实函数，则f的拉普拉斯算子定义为：

(1) f的拉普拉斯算子也是笛卡儿坐标系xi中的所有非混合二阶偏导数求和：

(2) 作为一个二阶微分算子，拉普拉斯算子把C函数映射到C函数，对于k ≥ 2。表达式(1)（或(2)）定义了一个算子Δ : C(R) → C(R)，或更一般地，定义了一个算子Δ : C(Ω) → C(Ω)，对于任何开集Ω。

对于阶跃状边缘，导数在边缘点出现零交叉，即边缘点两旁二阶导数取异号。据此，对数字图像{f（i，j）}的每个像素，取它关于x轴方向和y轴方向的二阶差分之和，表示为

运算模板

函数的拉普拉斯算子也是该函数的黑塞矩阵的迹,可以证明，它具有各向同性，即与坐标轴方向无关，坐标轴旋转后梯度结果不变。如果邻域系统是4 邻域，Laplacian 算子的模板为：

0	1	0
1	-4	1
0	1	0

如果邻域系统是8 邻域，Laplacian 算子的模板为：

1	1	1
1	-8	1
1	1	1

前面提过，Laplacian 算子对噪声比较敏感，所以图像一般先经过平滑处理，因为平滑处理也是用模板进行的，所以，通常的分割算法都是把Laplacian 算子和平滑算子结合起来生成一个新的模板。

平滑的算法多数是基于伞状结构的操作，所谓伞状结构，是以三角网格中某一个顶点P为中心，取所有与其相邻的顶点P1...Pn-1和边组成的一阶邻域结构。如下图所示：



拉普拉斯平滑

在前面的博客中，介绍了最为基础的拉普拉斯平滑算法的实现，简单的拉普拉斯平滑算法的原理是将每个顶点都移动到相邻顶点的平均位置，即采用所谓伞状算子：



在伞状结构中表示这样的过程如下图：



其具体的实现逻辑表述如下：

初始化Mesh的邻接点结构集

新建临时点集，用来存储点平滑后的位置

对所有Mesh中的顶点P

初始化临时向量为零向量

获取P的邻域点集Adj(P)

对所有领域点T，将其位置加到临时向量里

临时向量/=领域点集数

将临时向量的位置存入临时点集

对所有Mesh中的顶点P ，将P的位置修改为临时点集中对应点的位置

拉普拉斯平滑算法有很多进一步的变形，首先在求取平均位置时，可以采用不同的加权策略，例如对不同的邻接点采用不同的权值。一般来说，距离中心点P较远的邻接点，我们可以让他对P平滑后的位置影响小一点。这样就可以采用一种距离的倒数为权值的拉普拉斯平滑。



有时为了控制平滑的速率，也会引入参数lambda来控制平滑的速率，即从原来所执行的：


转变成

同时，平滑算法往往可以反复对Mesh执行，使得Mesh越来越光顺，迭代次数T也是平滑算法中重要的参数。

int n = neighbourcount(0, @ptnum);

// weights should add to 1, so 1/neighbors is an easy choice
float weight = 1.0 / n;

// vector average toward neighboring points
vector avg = {0,0,0};

foreach(int j; neighbours(0, @ptnum)) {
avg += weight * (point(0, "P", j) - @P);
}

// displace current point by the avg neighbor vector
// which is scaled by the factor lambda
float lambda = 1.0;
@P += lambda * avg;

基于曲率的平滑

基于曲率的平滑与上面介绍的拉普拉斯平滑的区别在于，基于曲率的平滑是沿着顶点P的法向量所在直线去移动P的，从下图可以看出，一般的拉普拉斯平滑会讲点P移动到类似于重心的位置，而基于曲率的平滑的移动位置是与法向量方向相反的，这样做的目的是要更好的保持模型原来的大致形状。



通过适当的选择拉普拉斯平滑中的权值ω，我们可以让基于曲率的平滑与拉普拉斯平滑达到算法上的统一。经过推导可知ω选择如下值的时候，即等价于将P如上图2中进行移动。


　　其中α和β的是P与其第i个邻接点Pi所构成边所对的两个角，（这里假设网格为流形，即每条边都只有至多两个三角形共有），下图显示了这两个角的位置：



这样点P的位置仍然由如下公式来决定，只是修改了权值计算方式：



Taubin平滑

拉普拉斯平滑虽然能够让Mesh的表面光顺，但迭代次数一旦多了，就会使得模型整体发生收缩现象。



Taubin平滑的原理基于了一部分数字信号处理方面的知识，用简单的话来表述，就是使用一个负收缩因子μ将拉普拉斯平滑照成的收缩再放大回去，其算法的主体采用了2个与拉普拉斯平滑算法相似的过程，一个过程采用正因子λ（0~1），另一个过程采用负因子μ（-1~0），每一次迭代都必须重复这两个过程。即有：

int n = neighbourcount(0, @ptnum);

// weights should sum to 1, so 1/neighbors is an easy choice
float weight = 1.0 / n;

// vector average toward neighboring points
vector avg = {0,0,0};

foreach(int j; neighbours(0, @ptnum)) {
avg += weight * (point(0, "P", j) - @P);
}

// displace current point by the avg neighbor vector
// which is scaled by the factor lambda
float lambda = 1.0;
@P += lambda * avg;

int n = neighbourcount(0, @ptnum);

// weights should sum to 1, so 1/neighbors is an easy choice
float weight = 1.0 / n;

// vector average toward neighboring points
vector avg = {0,0,0};

foreach(int j; neighbours(0, @ptnum)) {
avg += weight * (point(0, "P", j) - @P);
}

// displace current point by the avg neighbor vector
// which is scaled by the factor lambda
float lambda = -1.0;
@P += lambda * avg;