三维网格处理(补洞，去噪，测地线，细化)

参考https://www.cnblogs.com/shushen/p/5759679.html

　　在逆向工程中，由于设备或模型的原因，我们获取得到的三维模型数据往往并不完整，从而使得生成的网格模型存在孔洞，这对后续的模型分析会造成影响。下面介绍一种基于径向基函数（RBF：Radial Basis Function）的三角网格补洞方法。

Step 1：检测孔洞边界

　　三角网格是由一系列顶点（V）以及由这些顶点所构成的三角面片（F）所组成，由三角面片可以得到网格的边（E）。通常一条边连接两个三角面片，这种边称为网格内部边，而如果某条边仅连接一个三角面片，那么称这条边为网格边界边，所有的边界边按顺序连接之后就形成了网格的孔洞。

Step 2：初始化网格

　　为了使孔洞填充简单、健壮，可以采用最小角度法进行网格修补，具体步骤如下：

　　（1）得到孔洞边界点信息，计算边界边长度的平均值l。

　　（2）计算每个边界点的两条相邻边的夹角大小。

　　（3）找出具有最小夹角的边界点，计算它的两个相邻边界点的距离s，判断s < 2×l是否成立：若成立，则按下图所示增加一个三角形，若不成立则增加两个三角形。

　　（4）更新边界点信息。

　　（5）判断孔洞是否修补完整，若未完整，转（2），否则结束。

Step 3：最小二乘网格

　　初始化补洞得到的网格质量不是很好，我们需要优化网格顶点的位置，优化的条件是：

其中di为顶点vi的1环邻域顶点数。

　　上式可以用一个线性方程组来描述：LV = 0，其中L是Laplace矩阵，具体形式为

，矩阵L的秩等于n – k，n为网格顶点的数目，k为网格连通区域的个数，如果网格是全连通的，那么矩阵L的秩是n – 1。因此如果我们要使方程有解，需要加入m（m ≥ k）个控制顶点坐标vs作为方程的边界条件，在实际中我们将初始化网格的边界顶点作为控制顶点。

　　其实上述线性方程组的求解等价于如下能量函数最小化的求解：

　　能量函数的第一部分是使得网格顶点尽量光滑，即每个顶点位于其1环邻域顶点的中心，第二部分是为了控制顶点的位置满足要求。

　　最小二乘网格的优势是能够生成高质量的光滑网格，生成过程仅需要网格的拓扑连接关系和少数控制点的坐标信息。

关于控制点的思想，感觉就是b样条的思路。参考https://zhuanlan.zhihu.com/p/50626506

b样条是贝塞尔曲线的一般化。

先在平面上取两点，A和B。问题：如何在A、B之间做插值线段？如何把它们连起来？

嗯……当然是直接连起来就好了。

在将AB连起来以后，我们需要搞清线段AB上每一点的含义。换言之，如何用已知的A、B位置，来表明线段上任意一点？

比如，怎么表示AB正中间的点？

显而易见，答案是 1/2A+1/2B。

然后加点难度，如何表示偏向A一侧七分之一处的点？

6/7A+1/7B。 。因为P更靠近A一侧，所以A的值更大一些。

很简单，对吧。我们可以按照这个规律，将这种表述方式，推广到上的任意一点：从A向B慢慢移动，最一开始由于从A出发，A占完全的比重1，而B是0。随着点P缓慢向右推移，B的占比逐渐增大，A的占比逐渐减小。最终到达B点的时候，A成为了0，而B占完全的比重1。我们将这个比重值抽象成一个值t，可以获得下列公式：

这样，我们在已经知道AB位置的情况下，就可以通过修改t值，来简单的表述线段AB上的任意点了。

两个点的插值再怎么变化也只能做出一条直线。所以接下来，我们要引入一个新点C。

C的出现，使得线段的状态变得不确定了起来。为什么这么说呢？我们先考虑比较简单的情况，A和B的比值保持不变，因此还是t和1-t，这个比值，更多的时候我们将其称之为权重。同样的，C也需要一个权重。然后C的值——我们设为y——想要表示成这个样子：一开始y=0，随着t值增高，y先缓慢升高，在t=0.5的时候达到1，然后逐渐降低，在t=1的时候达到0。最终得到的是一个二次函数图像的形状。

结合A的比值t，B的比值1-t，曲线整体的公式为：

最终的整体曲线效果类似下图，最初A的比值为1，后面随着A减少和B、C的逐渐增大，A、B、C各占0.5，再之后B的比值逐渐增大，A、C开始减少，最终曲线轨迹完全贴合B。

好了，现在我们有了一条插值曲线，只需要通过修改t值，我们就可以表述曲线上的任意一点。

不过，真正的B-样条曲线的规则没有这么简单，它的插值规则是不确定的。不过不管怎么变，目标是一样的：最终，只需要通过修改t值，就可以表述任意曲线上的任意点。

做曲线也要按照数学的基本规律啊，对不对？当然我们提供的基本参数也是很重要的……

到目前为止，我们讲的思路是：首先做两点直线，然后做三点曲线，按这个规律接下来我们应该逐步推理多点曲线的情况，再从中获取B样条的表达方式，但那就太复杂了。所以接下来我们的思路会反过来：首先介绍B样条曲线的一些基本参数，通过这些参数就可以生成B样条曲线。不论从哪种角度出发，最终的目标都是一样的，只需要修改t值，就可以表述曲线上的任意点。

好了，接下来开始说B样条曲线的基本参数。

B样条曲线的基本参数中其实就几样，t，阶数，控制点列表，节点表，基本函数表。

t值的作用前面已经说了很多了，在此不提。控制点列表代表一系列需要用户提供的顶点。实际上上面直线中的AB，曲线中的ABC，都可以理解为控制点。

我们还没接触的，实际上只有阶数，节点值和基本函数表了。

剩下的参数中最基本的当然是阶数。每个t值都是通过控制点和权重相乘计算得出的结果。而阶数越高，生成每个t值所需要的控制点数越多。

所以，阶数到底是什么？阶数=所有权重中t值的最高次幂。

觉得费劲？有个更简单的：阶数=生成t值所需要的控制点数-1。

像前面所提到的曲线，生成t值，需要A、B、C三个控制点。因此它是一个3-1=2阶的曲线（不过它不是B-样条）。而用比较绕的最高幂次方式理解，同样是该曲线，A的权重是 t，B的权重是 

 ，C的权重是 

 ，最高幂次是 

 中的2，因此同样可以得出是2阶曲线。

在理解了阶数的概念以后，就可以说节点表了。

节点表是生成基本函数表的关键参数，大小严格等于控制点数量+阶数+1。节点表的参数是人为设置的，如果不是今天在讲B-样条曲线，你想怎么设置都可以。但B-样条还是要遵守一些规律的。

一般设置的方法有两种：顺序方法和Clamped方法。前者用于制作标准的B-样条开曲线和闭曲线，后者用于制作一种比较实用的B-样条曲线。

顺序列表只需要从0-1线性递增设置即可，Clamped列表则需要将前后各阶数+1个节点设置成0。

举个例子就能讲的很清楚：假设曲线有6个控制点，阶数是3阶，那么节点表大小=6+3+1=10。

如果是顺序列表，只需要按顺序设置：

如果是Clamped列表，由于是3阶，前面3+1个参数均设置为0，后面3+1个参数均设置为1，然后剩余参数均匀递增：

基本函数表本质上是一个递归方程，但同时它也是一个中间参数。我想首先将基本函数表的公式列出来，看了十有八九会头大，不过没关系，后面我会一点一点解剖这个公式：

这个大家会很熟悉，线性插值嘛。

t的话最容易理解，就是前面反复提到的t值。

knot代表节点，也就是我们前面提到的节点表。

而 knot_i 代表节点表中的第 i个元素。

而 

 就是基本函数表的参数了。没错，基本函数表是一个二维数组。参数i和deg分别表示第几个元素和阶数。 所以

的意思就是用户输入值为t时，基本函数表在第deg阶的第i个元素的值。

这时候我们再纵观整个长公式，发现了吗？整个公式实际上是两个部分相加。而加号两侧的公式格式一致：一个通过节点表 knot和i经过一系列计算得出的权重值和一个比当前更低一阶的基本函数表值的乘积。

所以为什么说基本函数表是递归方程，这就是原因。当前deg阶的元素需要通过两个deg-1阶的元素获得，deg-1阶的元素则需要通过deg-2阶的元素获得……以此类推，直到递归deg次以后，回退到0阶为止。

但是，回退到0阶的时候怎么算呢？

B-样条算法规定，回退到0阶时使用以下公式：

即，0阶的时候，t处于第i个knot值和第i+1个knot值之间的时候，才会等于1，其他时候都等于0.

怎么样，是不是突然发现了knot节点表的意义？实际上节点表的作用就是确认一下t在元素表的哪个位置，然后我们只需要找出基本函数表在deg中能够递归到该位置的值，就可以在后续的最终计算中，得到t值在B-样条曲线的对应位置了。

函数表应该是我们实现准备好的。

Step 4：径向基函数（RBF）隐式曲面

　　径向基函数是一个仅依赖于离控制点c距离的函数，即h(x,c) = h(||x - c||)，距离通常是欧式距离，也可以是其他形式距离。径向基函数网络是一个三层BP网络，其可以表示为多个基函数的线性组合：

　　下面列出2个最常用的基函数形式：

　　（1）Gaussian：h(r) = e-(εr)2

　　（2）Polyharmonic spline：h(r) = rk，k = 1、3、5、… 或h(r) = rkln(r)，k = 2、4、6、…

　　利用径向基函数网络并通过给定控制点xi和对应的值fi之后，可以求解得到网络中的系数λi，因此径向基函数网络能够解决空间散乱数据点的平滑插值问题，函数的零等值面就是我们要求的曲面。

　　在实际求解时函数f(x)表达式中通常会增加一个一次多项式P(x)，如下：

，其中：

　　将控制点位置xi和值fi代入函数f(x)可以得到如下方程，求解之后即可确定隐式曲面f(x)。

　　控制点xi可分为三类：

　　（1）边界控制点：曲面通过的点，f(xi) = 0

　　（2）外部控制点：将边界控制点沿法向正方向移动一小段距离而得到的控制点，取f(xi) = -1

　　（3）内部控制点：将边界控制点沿法向负方向移动一小段距离而得到的控制点，取f(xi) = 1

　　计算时我们选择的基函数为h(r) = r^3，其实此时隐式曲面函数满足Δ^3f = 0，将基函数代入后得到隐式曲面的表达式如下：

function options = RBF_Create(x, y, type)

    % x: input data; y: input data value

    options = [];

    options.nodes = x;

    switch type

        case 'linear'

            options.phi = @rbfphi_linear;

        case 'cubic'

            options.phi = @rbfphi_cubic;

    end    

    [n, dim] = size(x);

    A = zeros(n,n);

    for i = 1:n

        r = normrow(bsxfun(@minus, x(i,:), x));

        temp = feval(options.phi, r);

        A(i,:) = temp';

        A(:,i) = temp;

    end

    % Polynomial part

    P = [ones(n,1) x];

    A = [ A      P

          P' zeros(dim+1, dim+1)];

    B = [y; zeros(dim+1, 1)];

    coeff = A\B;

    options.coeff = coeff;

end

% Radial Base Functions

function u = rbfphi_linear(r)

    u = r;

end

function u = rbfphi_cubic(r)

    u = r.*r.*r;

end

Step 5：牛顿插值

　　为了得到插值隐式曲面的网格，我们需要把最小二乘网格的顶点投射到隐式曲面上，这里我们采用牛顿迭代法：

　　最小二乘网格的顶点位置作为迭代初始条件，当||xk+1 – xk|| < ε时，迭代停止，ε为设定的误差。上式中梯度表达式如下：

效果：

受图像双边滤波算法的启发，[Fleishman et al. 2003]和[Jones et al. 2003]分别提出了利用双边滤波算法对噪声网格进行光顺去噪的算法，两篇文章都被收录于当年的SIGGRAPH，至今引用超500余次。虽然从今天看两篇文章的去噪效果还不算非常好，但是其中的思想是值得学习的。图像双边滤波算法可以参考http://blog.csdn.net/abcjennifer/article/details/7616663，图像双边滤波器由空间域核与值域核组成，在图像的特征区域，自身像素值与周围像素值差别较大，这时起主导的值域核决定周围像素值的权重系数大幅降低，因此双边滤波能够在去噪的同时保持细节特征。

　　[Fleishman et al. 2003]提出的网格去噪算法为一种迭代方法：首先定义每个顶点u的计算域，空间域定义为顶点u到计算域内其他顶点p的距离，值域定义为计算域内其他顶点p到顶点u切平面的带符号距离，然后利用双边滤波方法计算顶点u在法向上的移动距离，更新完网格所有顶点位置后即完成一次迭代。具体公式如下：

其中空间域核Wc(x) = e^(-x2/2σc2)，值域核Ws(x) = e^(-x2/2σs2)，N(u)为满足||u - pi|| < ρ = 2σc的顶点集{pi}，I(p)为顶点p到顶点u切平面的带符号距离，最后计算得到的I(u)即为顶点u需要在其法向上移动的距离。

效果：

　　[Jones et al. 2003] 提出的网格去噪算法不需要迭代，具体公式如下： 

其中f和g分别为空间域核以及值域核，核函数都为高斯函数，aq为三角片的q的面积，cq为三角片的q的中心，S为满足||cq - p|| < 2σf的三角片集，最后直接计算得到新顶点的位置p’。

　　那么如何得到∏q(p)，∏q(p)为顶点p在三角片q切平面上的投影点，但该投影点并不是在原始网格上计算。文章提出忽略双边滤波公式中的值域核函数g，并令∏q(p) = cq，即利用高斯滤波得到一个虚拟网格，而∏q(p)即为顶点p到虚拟网格中三角片q切平面上的投影点。

效果：

多次提到了核方法，参考https://blog.csdn.net/u014662865/article/details/80970470

https://blog.csdn.net/fjssharpsword/article/details/79092082

首先通俗理解一下核函数存在的意义，按照其他一些官方上的解释，核函数就是为了将低维空间上的点映射到高维空间上，是为了方便将不能用线性分割的数据转化成可以线性分割的数据（白话解释），如图所示：

如图所示，左面的图为为原空间，右面的图为映射后的空间，从图中也可以看出来，左面图要用一个椭圆才能将两个类别分割开来，而右面的图要用一个超平面就可以分割开，也如图上的共识所示，原空间点左边为（x1,x2）,经过某个函数或者某种计算方法，转化为特征空间上点坐标为(z1,z2,z3)，所以我们之前的说法是正确的，将低维空间转化到高维空间大概率可以对其中的点进行线性分割。

这样，我们第一步就理解了，就是在低维空间上的点通过某一函数转化为高维空间上，更有助于线性分类。那么我们开始理解下一步。

在核函数中，会涉及到两个点的内积的计算，本节主要论述计算内积的意义。

两个点之间的内积是有一定意义的，可以通过两个点的内积计算距离和两个向量之间的角度，计算方法如下所示：

其中的

字符代表着一个点。而根据两个点的距离，夹角等信息可以得到某些模型的超平面计算，讨论完内积的意义，下节将着重讨论核函数的重点内容。

由于我们原特征空间是线性不可分的，所以要将原特征空间中的点x经过

函数，转化到高维特征空间上，理论来讲，将低维空间转化到高维空间上，计算复杂度会变高很多，但是由于核函数的引用，巧妙的解决了这个问题。

我们需要优化的函数如下：

其中c为正负两类样本点连线的中点。如图所示：

w = (c+) -(c-),为红色线的向量

本函数计算绿色的线和红色线向量夹角，如果大于90度属于负类，否则属于正类。

计算过程如下：

其中的k(x,x1)，代表核函数，也就是高维空间两个点的内积，不同的内积定义会有不同的核函数。

所谓核函数，则是在原空间上两点内积的一个函数得到的，核函数样例如下所示：

将样本点从二维空间转化到三维空间：

从上述公式可以看出，核函数是一个可以用原特征空间上点内积的方式经过运算转化成高维空间点内积，而不必完全由高维空间上的点进行计算，从而达到降低运算复杂度的作用。

也可以参考大神的视频：

核函数怎么定义：

对一个二维空间做映射，选择的新空间是原始空间的所有一阶和二阶的组合，得到了五个维度；如果原始空间是三维，那么我们会得到 19 维的新空间，随着维度的增加，非线性组合的数目是呈爆炸性增长，维度无穷无尽。

核函数避免了求内积。

核函数有如此的好处，但构造核函数本身就是一件不容易的事。

巧妙的核技巧令人垂涎：通过核函数，用低维的计算量得到了高维的结果，没有增加计算复杂度的同时，得到了性质更好的高维投影。

只要涉及到内积运算，都可运用核函数替代来得到高维投影的内积，于是虽然寻找核函数困难，根据Mercer定理还是有规律可循。

核矩阵就是格拉姆矩阵。

被证明可用的核函数有：

总结：核方法在线性与非线性间架起一座桥梁，通过巧妙地引进核函数，避免了维数灾难，没有增加计算复杂度却实现了高维点积运算。核函数的定义需满足Mercer定理，被证明可用的主要有三类线性核、高斯核、多项核。

　基于两步法的网格去噪算法顾名思义包含两个步骤：首先对网格表面的法向进行滤波，得到调整后的网格法向信息，然后根据调整后的法向更新顶点坐标位置，下面介绍三篇该类型的文章。

　　[Sun et al. 2007]文章首先介绍了当前法向滤波方法以及顶点坐标更新方法，然后提出自己的法向滤波方法和顶点坐标更新方法。

，其中N(i)为三角片fi的邻域三角片集，有两种选择，见下图，f(x) = x^2，0 ≤ T ≤ 1，T的取值决定了邻域三角片N(i)中有多少三角片法向会参与计算，当T = 1时，只有与ni相等的邻域三角片法向会参与计算，而当T = 0时，所有邻域三角片法向都会参与计算。当网格有明显棱边时，T可以取相对较小值，而当网格比较光滑时，T可以取相对较大值。

Ⅰ基于顶点（vertex based）的三角片邻域；Ⅱ表示基于边（edge based）的三角片邻域

　　顶点坐标更新方法：

　　1.求解线性方程组：

　　2.梯度下降法：

其中Nv(i)为顶点i的1环邻域顶点集。

　　3. [Sun et al. 2007]文章顶点更新方法：

其中Fv(i)为顶点i的1环邻域三角片集。

效果：

　　[Zheng et al. 2011]文章提出了两种三角片法向更新方法，而其顶点坐标更新方法采用和[Sun et al. 2007]文章一样的方法。

效果：

参考https://blog.csdn.net/program_developer/article/details/80177487

当P=0时，也就是L0范数，L0范数并不是一个真正的范数，它主要被用来度量向量中非零元素的个数。用L-P定义可以得到的L-0的定义为： 

这里就有点问题了，我们知道非零元素的零次方为1，但零的零次方为0，非零数开零次方都是什么鬼，很不好说明L0的意义，所以在通常情况下，大家都用的是： 

表示向量x中非零元素的个数。

　但是当有角度接近0时，微分算子的权重会变成inf，因此文章又提出了一种优化后的微分算子表达式：

效果：

感觉就很像蒙特卡洛方法。

  那么通过求解该线性方程组就可以得到各个点到达目的地N点的概率，以上就是一维随机游走算法原理。

[Grady et al. 2006]提出了利用随机游走思想来分割二维图像，文章将图像考虑成一张图（Graph），每个像素对应图中一个节点，根据亮度差值定义节点间的权重（相当于一维随机游走中向左和向右的概率），然后用户指定前景（foreground）和背景（background）标签（相当于一维随机游走中N点和0点），通过求解线性方程组就可以得到各个像素点属于前景或背景的概率，如果将阈值概率设置为0.5，那么就可以分割得到期望的图像区域。

       [Lai et al. 2008]将这种思想扩展到三维网格分割，文章将网格中每个三角片对应图中一个节点，利用相邻三角片之间的二面角来定义节点之间的权重，具体如下：

同样通过求解线性方程组可以得到网格分割效果。

　　[Zhang et al. 2010]对[Lai et al. 2008]的网格分割算法做了部分改进，文章将网格中每个顶点对应图中一个节点，由于一个网格的三角片数量通常是顶点数量的2倍左右，这样求解的方程变量数就会减少一半左右，计算速度就会得到提高。

效果：

参考https://www.cnblogs.com/shushen/p/5174421.html

测地线又称为大地线，可以定义为空间曲面上两点的局部最短路径。测地线具有广泛的应用，例如在工业上测地线最短的性质就意味着最优最省，在航海和航空中，轮船和飞机的运行路线就是测地线。[Crane et al. 2013]提出了利用热运动方程来计算网格测地线的方法，可以想象一下，当一根烫的针尖接触到曲面上的一点时，热量会随着时间的推移而扩散，测地距离因此可以和热运动相联系。具体算法过程如下图所示：

　第一步：热运动方程用来描述热的传播状态：

　　将热运动方程离散化并整理后得到：

其中：id为单位矩阵，t为时间间隔，Δ为离散Laplacian算子，ut为t时刻的热状态，u0为初始时刻的热状态。

　　第二步：第一步计算得到的热梯度方向与测地距离的梯度方向相同，由Eikonal方程知道测地距离的梯度为单位向量，于是通过归一化热梯度我们得到测地距离的梯度：

function [D] = geodesics_in_heat(V, F, src)   

    % choose time step

    c = 5;

    t = c * mean(doublearea(V, F))/2;

    %% Step 1: Integrate the heat flow for some fixed time t

    L = cotmatrix(V, F);

    M = massmatrix(V, F, 'barycentric');

    nV = size(V, 1);

    u0 = zeros(nV, 1);

    u0(src) = 1;

    A = M - t*L;

    B = M*u0;

    nsrc = length(src);

    % 1.1 dirichlet condition

    hole = Cal_Boundary(F);

    if isempty(hole)

        boundary = [];

    else

        boundary = hole.boundary.edge(:,1)';

    end

    b = setdiff(boundary, src);

    nb = [b, src];

    Acons = sparse([1:length(nb)], nb, ones(1,length(nb)), length(nb), nV);

    Bcons = [zeros(length(b), 1); ones(nsrc, 1)];

    % 硬约束        

    r = setdiff([1:nV], nb);

    uD = [A(r,:);Acons]\[B(r,:);Bcons];

    % 1.2 neumann condition

    Acons = sparse([1:nsrc], src, ones(1,nsrc), nsrc, nV);

    Bcons = ones(nsrc, 1);

    % 硬约束  

    r = setdiff([1:nV], src);

    uN = [A(r,:);Acons]\[B(r,:);Bcons];

    % averaged boundary condition

    u = 0.5*(uN + uD);

    %% Step 2: Evaluate the vector field X

    G = grad(V, F); % nF*3 by nV matrix(梯度算子 - 所有三角片顶点基函数)

    grad_u = reshape(G*u, size(F,1), 3); % nF by nV matrix(所有三角片中u的梯度)

    grad_u_norm = sqrt(sum(grad_u.^2, 2));

    normalized_grad_u = bsxfun(@rdivide, grad_u, grad_u_norm+eps);

    X = -normalized_grad_u;

    %% Step 3: Solve the Poisson equation

    Div = div(V, F); % 散度算子

    div_X = Div*X(:); % #nV by #nF*3

    Lcons = sparse([1:nsrc], src, ones(1,nsrc), nsrc, nV);

    div_Xcons = zeros(nsrc, 1);

    % 硬约束

    r = setdiff([1:nV], src);

    D = [L(r,:);Lcons]\[div_X(r,:);div_Xcons];

    D = D - min(D);

end

参考https://www.cnblogs.com/shushen/p/5251070.html

　　下图描述了细分的基本思想，每次细分都是在每条边上插入一个新的顶点，可以看到随着细分次数的增加，折线逐渐变成一条光滑的曲线。曲面细分需要有几何规则和拓扑规则，几何规则用于计算新顶点的位置，拓扑规则用于确定新顶点的连接关系。下面介绍两种网格细分方法：Catmull-Clark细分和Loop细分。

Catmull-Clark subdivision：

　　Catmull-Clark细分是一种四边形网格的细分法则，每个面计算生成一个新的顶点，每条边计算生成一个新的顶点，同时每个原始顶点更新位置。下图为Catmull-Clark细分格式的细分掩膜，对于新增加的顶点位置以及原始顶点位置更新规则如下：

5.网格边界E-顶点位置：

　　设边界边的两个端点为v0、v1，则新增加的顶点位置为v = 1/2*(v0 + v1)。

效果：

function [VV, FF, S] = CC_subdivision(V, F, iter)

    % Catmull_Clark subdivision

    if ~exist('iter','var')

        iter = 1;

    end

    VV = V;

    FF = F;

    for i = 1:iter

        nv = size(VV,1);

        nf = size(FF,1);

        O = outline(FF);

        original = 1:nv;

        boundary = O(:,1)';

        interior = original(~ismember(original, boundary));

        no = length(original);

        nb = length(boundary);

        ni = length(interior);

        %% Sv

        Etmp = sort([FF(:,1) FF(:,2);FF(:,2) FF(:,3);FF(:,3) FF(:,4);FF(:,4) FF(:,1)],2);

        [E, ~, idx] = unique(Etmp, 'rows');

        Aeven = sparse([E(:,1) E(:,2)], [E(:,2) E(:,1)], 1, no, no);

        Aodd = sparse([FF(:,1) FF(:,2)], [FF(:,3) FF(:,4)], 1, no, no);

        Aodd = Aodd + Aodd';

        val_even = sum(Aeven,2);

        beta = 3./(2*val_even);

        val_odd = sum(Aodd,2);

        gamma = 1./(4*val_odd);

        alpha = 1 - beta - gamma;

        Sv = sparse(no,no);

        Sv(interior,:) = ...

            sparse(1:ni, interior, alpha(interior), ni, no) + ...

            bsxfun(@times, Aeven(interior,:), beta(interior)./val_even(interior)) + ...

            bsxfun(@times, Aodd(interior,:), gamma(interior)./val_odd(interior));

        Sboundary = ...

            sparse([O(:,1);O(:,2)],[O(:,2);O(:,1)],1/8,no,no) + ...

            sparse([O(:,1);O(:,2)],[O(:,1);O(:,2)],3/8,no,no);

        Sv(boundary,:) = Sboundary(boundary,:);

        %% Sf

        Sf = 1/4 .* sparse(repmat((1:nf)',1 ,4), FF, 1);

        i0 = no + (1:nf)';

        %% Se

        flaps = sparse([idx;idx], ...

                       [FF(:,3) FF(:,4);FF(:,4) FF(:,1);FF(:,1) FF(:,2);FF(:,2) FF(:,3)], ...

                       1);

        onboundary = (sum(flaps,2) == 2);

        flaps(onboundary,:) = 0;

        ne = size(E,1);

        Se = sparse( ...

                [1:ne 1:ne]', ...

                [E(:,1); E(:,2)], ...

                [onboundary;onboundary].*1/2 + ~[onboundary;onboundary].*3/8, ...

                ne, ...

                no) + ...

                flaps*1/16;

        %% new faces & new vertices

        i1 = no +   nf + (1:nf)';

        i2 = no + 2*nf + (1:nf)';

        i3 = no + 3*nf + (1:nf)';

        i4 = no + 4*nf + (1:nf)';

        FFtmp = [i0 i4 FF(:,1) i1; ...

                 i0 i1 FF(:,2) i2; ...

                 i0 i2 FF(:,3) i3; ...

                 i0 i3 FF(:,4) i4];

        reidx = [(1:no)'; no+(1:nf)'; no+nf+idx];

        FF = reidx(FFtmp);

        S = [Sv; Sf; Se];

        VV = S*VV;

    end

 end

Loop subdivision：

　　Loop细分是一种三角形网格的细分法则，它按照1-4三角形分裂，每条边计算生成一个新的顶点，同时每个原始顶点更新位置。下图为Loop细分格式的细分掩膜，对于新增加的顶点位置以及原始顶点位置更新规则如下：

4.网格边界E-顶点位置：

　　设边界边的两个端点为v0、v1，则新增加的顶点位置为v = 1/2*(v0 + v1)。

to be continued