本发明属于三维网格处理技术领域，具体涉及一种基于法向量平均的常平均曲率曲面的构造方法。
背景技术：
平均曲率处处为常数的曲面被称为常平均曲率(cmc)曲面，极小曲面是cmc曲面的一种特殊情况，它的平均曲率处处为零。在许多建筑设计和工程创作中，cmc曲面因具有较高的美学价值和结构稳定性而备受欢迎，建筑中的张力膜结构可以被建模为极小曲面。充气建筑，包括充气圆顶和外壳，也被建模为cmc曲面。
cmc曲面构造是计算机图形学中的一个具有挑战性的问题。构建cmc的方法有许多种，通常使用的cmc表面建模方法是优化某些能量函数，例如willmore能，找到曲面的离散近似，通常是三角网格的形式。例如，一个cmc曲面可以倍计算为在体积约束下面积函数最小的曲面。基于这一性能，可以提出许多通过在网格上求解约束优化问题来计算cmc曲面的方法。
通常，一个网格模型由网格顶点和网格连通度来定义。当使用网格进行形状表示时，网格质量的两个主要标准是网格顶点的分布和边连通度。在许多几何任务中，如数值模拟和差分表面性质的估计，都需要高质量的网格，如网格光滑度。网格质量差是现有cmc曲面建模方法的一个主要问题。这些方法缺乏一种优化网格连通性的机制。通常，这种方法只更新网格顶点的位置，同时保持初始网格的相同的边缘连通性，从而将初始网格迭代到最终的形状。其结果是，由于元素表面形状不佳，而导致最终的网格是真实cmc表面的较差表示。另外，利用曲率流构建cmc曲面也是一种有效方法，例如共形曲率流，单位法向流等，利用曲率流构建cmc曲面具有鲁棒性以及高效性。
技术实现要素：
本发明的目的在于提供一种基于法向平均的cmc曲面构建方法。该方法可以解决现有的优化能量函数方法对初始网格质量要求比较高,适应性差，效率低等问题。该方法具有较强的鲁棒性和适应性，具有较快的收敛速度，并且对于一些复杂的三维网格模型，能收敛达到很好的cmc曲面效果，并且不会出现三角形翻转，网格撕裂等现象。
本发明所采用的技术方案是：一种基于法向平均的常平均曲率曲面构建方法，其特征在于，包括以下步骤：
步骤1：读取网格，计算网格每个面的单位法向量以及其邻居面平均单位法向量；
步骤2：根据每个面的单位法向量以及其邻居面平均单位法向量构造旋转矩阵；
步骤3：计算网格顶点的拉普拉斯矩阵，利用旋转矩阵计算网格顶点旋转后坐标矩阵，利用拉普拉斯矩阵和顶点旋转后的坐标求解线性方程组，得到顶点的实际新坐标；
步骤4：更新网格顶点坐标，输出显示网格，至此完成一次迭代；判断终止条件，若符合条件则停止迭代，得到最终的输出模型，否则将本次输出网格作为输入返回执行步骤1。
本发明的有益效果在于：本发明提供了一种基于法向平均的cmc曲面的构建方法，通过法向量平均求解线性方程组迭代更新网格定点位置来将网格变形为cmc曲面。其中可以通过参数控制法向量平均方法，使得方法对于不同的网格具有适应性以及高效性，提供了一个更加鲁棒和高效性的cmc曲面构建方法。
附图说明
图1为本发明实施例的流程图；
图2(a)为本发明实施例输入的原始封闭模型；
图2(b)为本发明实施例迭代60次输出的模型；
图2(c)为本发明实施例迭代120次输出的模型；
图2(d)为本发明实施例迭代180次输出的模型；
图2(e)为本发明实施例迭代240次输出的模型
图2(f)为本发明实施例最终的输出cmc模型；
图3为本发明实施例构建拉普拉斯矩阵时所使用的三角形内角示意图；
图4(a)为本发明实施例输入的原始带边界模型；
图4(b)为本发明实施例迭代60次输出的模型；
图4(c)为本发明实施例迭代120次输出的模型；
图4(d)为本发明实施例最终的cmc曲面模型。
具体实施方式
为了便于本领域普通技术人员理解和实施本发明，下面结合附图及实施例对本发明作进一步的详细描述，应当理解，此处所描述的实施示例仅用于说明和解释本发明，并不用于限定本发明。
请见图1，本发明提供的一种基于法向平均的常平均曲率曲面构建方法，包括以下步骤：
步骤1：读取网格，计算网格每个面的单位法向量以及其邻居面平均单位法向量；
步骤1：在本实施方式中，输入的三维网格模型如图2(a)所示。
步骤101：计算计算每个面的单位法向量currentnormal。每个面的单位法向量为从某一顶点出发的两条边的叉积。计算方法为：
其中pi,i＝1,2,3，为每个面的三个顶点。
步骤102：计算每个面的邻居平均单位法向量targetnormal。根据参数的不同，计算平均法向量的方法可以有多种，主要包括算术平均以及加权平均方法。
算术平均法向量计算方法为：
其中，nj为面j的法向量，neighbor(i)为面i的n圈邻居面；
加权平均法向量计算方法为：
其中，nj为面j的的法向量，areaj为其面积。
利用不同的方法对算术平均方法和加权平均方法改进，可以使得本方法的适用性更强。
对算数平均方法的改进包括所有面算术平均，相邻加速算数平均，动态圈数加速算数平均方法。计算方式如下：
所有面算术平均方法：
其中，face为网格所有面；
相邻加速算数平均方法：
其中，f表示网格所有面的个数；acclerator＝0.1，为加速因子；将结果单位化，即得到平均单位法向量。
动态圈数加速算数平均方法：
其中，k满足的条件为：
对加权平均方法的改进包括所有面加权平均，相邻加速加权平均，动态圈数加速加权平均方法。计算方法如下：
所有面加权平均方法：
相邻加速加权平均方法：
其中，f表示网格所有面的个数；acclerator＝0.1，为加速因子；将结果单位化，即得到平均单位法向量。
动态圈数加速加权平均方法：
其中，k满足的条件为：
计算每个面的邻居平均单位法向量targetnormal，本例采用的平均方法为算术平均方法。
步骤2：根据每个面的单位法向量以及其邻居面平均单位法向量构造旋转矩阵；
计算rot旋转矩阵的方法为：
其中(x,y,z)＝currentnormal×targetnormal，θ为向量currentnormal与targetnormal之间的夹角。
步骤3：构建网格顶点的拉普拉斯矩阵lefta，利用旋转矩阵rot计算旋转之后的顶点坐标矩阵rightb，求解线性方程组leftax＝rightb，即得到网格顶点的新坐标矩阵x。
步骤301：计算网格顶点的拉普拉斯矩阵lefta。lefta是大小为v×v的稀疏矩阵，其中v为网格顶点数。稀疏矩阵lefta的元素构造如下：
其中，xi～xj表示顶点xi与xj相邻接，αj与βj为边ij所在的两个三角形中与边ij相对的两个内角，如图3所示。
步骤302：根据每个网格顶点的实际坐标与旋转矩阵rot求得旋转后的坐标矩阵rightb。rightb是一个大小为v×3的矩阵，其构造方法为：
其中，n1(i)表示顶点xi的1-ring邻居顶点，xj表示与顶点xi相邻接的所有顶点，edgeij＝xi-xj，为边ij的向量表示，等式右边第一项中rotij表示从顶点xi出发的半边所代表的面的旋转矩阵。
weigthhalfedgeij＝cotγ/2；
其中，γ表示半边ij所代表的面存在时，三角形边ij所对的内角。
同理，等式右边第二项中rotji表示从指向顶点xi的半边所代表的面的旋转矩阵。
weigthhalfedgeji＝cotγ/2；
其中，γ表示半边ji所代表的面存在时，三角形边ji所对的内角。
步骤303：求解线性方程组leftax＝rightb，即得到网格顶点的新坐标矩阵x。
步骤4：根据顶点坐标矩阵x更新网格顶点，输出网格，至此完成一次迭代。然后计算计算顶点平均曲率，当顶点平均曲率足够平均时，表示已得到最终的常平均曲率曲面，迭代停止，否则将得到的网格作为输入网格返回步骤1，直到完成迭代，输出cmc曲面结果，迭代过程如图2(b)、2(c)、2(d)、2(e)所示，最终输出模型如图2(f)所示。
步骤401：计算每个顶点的voronoi面积，计算方法如下：
其中，xj表示与顶点xi相邻接的所有顶点,αj与βj分别为边ij所在的两个三角形中与边ij相对的两个内角。
步骤402：计算每个顶点的离散平均曲率，方法如下：
步骤403：判断网格每个顶点的平均曲率是否为常数，若是，表示已得到最终的常平均曲率曲面，迭代停止，否则将得到的网格作为输入网格返回步骤1，直到完成迭代，输出cmc曲面结果。
在本实施例中，通过法向平均方法，将输入网格变形为cmc曲面。本实施方式提供的cmc曲面构建方法，能够有效地解决当前cmc曲面构建方法存在的适应性差，效率低等问题。图4(a)、4(b)、4(c)、4(d)展示了将本方法应用在有边界模型得到的迭代过程，证明了本方法处理有边界模型的有效性。
应当理解的是，本说明书未详细阐述的部分均属于现有技术。
应当理解的是，上述针对较佳实施例的描述较为详细，并不能因此而认为是对本发明专利保护范围的限制，本领域的普通技术人员在本发明的启示下，在不脱离本发明权利要求所保护的范围情况下，还可以做出替换或变形，均落入本发明的保护范围之内，本发明的请求保护范围应以所附权利要求为准。