0 引言

纹理合成和纹理映射是计算机虚拟现实的主要技术，能够模拟物体表面的颜色细节或几何细节，从而大幅度增强场景的真实感.基于样图的纹理合成技术^[1-6]，合成速度较快，内存占用小，但其效率的提升是以降低绘制精度为代价的，且图像纹理不支持遮挡、阴影、轮廓等重要的效果.近几年，研究者们提出了采用几何纹理代替图像纹理来表示物体的表面细节，从而进一步提高了绘制的精度.几何纹理具有更丰富、细致的表现能力，随着显卡硬件的提升，几何纹理在虚拟现实技术中得到了广泛应用.几何纹理和图像纹理一样都具有自相似性，但几何纹理是由不规则拓扑连接的网格构成的，因此其数据结构比图像纹理更为复杂，所以在合成过程中往往需要大量人工交互才能得到较理想的合成结果.同时，基于全局邻域搜索的几何纹理合成方式，使得几何纹理合成过程中普遍存在高计算量、高存储占用等问题.因而，自动、快速和高效的几何纹理合成成为几何纹理是否能够得到广泛应用的关键.

近年来，随着计算机硬件技术的飞速发展，图形处理单元GPU的计算能力呈几何级数增长，促使通用工程计算软件由CPU中央处理向CPU/GPU协同处理的方向发展.GPU具有出色的图元运算能力、浮点计算能力、低能耗高带宽、可靠的并行架构以及灵活的可编程性，为解决几何纹理合成中的计算瓶颈提供了可靠的现实基础.本文提出了一种新的基于GPU加速的几何纹理合成方法，该方法基于块的几何纹理合成方法，引入GPU多线程并发技术，通过将最佳几何子块匹配和几何子块拼接等移植到GPU上执行，从而加快几何纹理合成速度，高效合成大尺寸的新的几何纹理.本文方法不仅内存占用更少，而且在保证合成质量的同时，可显著提高几何纹理的合成速度.

1 相关工作

2007年QIN等^[7]提出了一种欧拉3D纹理的合成方法.通过构建输入样本的欧拉矩阵，将输入样本按照矩阵中元素的顺序排列，以生成体纹理，亦作欧拉3D纹理.该方法在进行纹理合成时无须人工干预，是一种全自动的体纹理合成方法，但由于无法很好地解决欧拉矩阵中体素的色彩更新问题，容易产生比较明显的缝隙.2008年，严志程等^[8]提出了一种基于方向场的体纹理合成方法.在输入的体素样本中搜索邻域，优化邻域匹配算法逐个合成体素；在匹配体素样本像素邻域和体素邻域时，通过考虑各体素上的方向能够生成各向异性体纹理.2010年，PIETRONI等^[9]提出了在三维模型表面上合成体数据的方法.同年，WANG等^[10]提出了新的体纹理合成方法，即采用一个随机向量表示体素与体素之间具有平滑颜色过渡的混合重叠区域，并对此向量按照待合成体纹理的体积进行分区，通过定义分区之间的颜色距离公式，以合成新的体纹理.2011年，江巨浪等^[11]提出了一种基于二维样图的体纹理快速生成算法，根据纹理的空间分布特征，设计二维样图在三维空间中的运动路径，然后使样图中的像素通过其运动轨迹对三维空间着色生成体纹理空间，最终合成体纹理.该方法仅对大理石和木材纹理有着较好的合成效果，无法通用.上述方法^[7-11]的合成速度较快，但由于大多体素结构过于简单且内容单一，合成结果缺乏丰富和细致的表现力，因此研究者们提出了采用几何体代替体素合成几何纹理.

2004年，BHAT等^[12]首次提出了一种基于几何体的几何纹理合成方法.该方法主要借鉴三维模型表面上二维纹理合成的方法，采用几何体取代二维纹理中的像素单位，通过相邻区域局部体素的三维空间坐标值构建相邻关系，结合网格编辑技术生成新的几何纹理.该方法需要大量的人工交互，且耗时较长.2006年，ZHOU等^[13]提出了基于网格缝合(Mesh Quilting)在曲面上合成几何纹理的方法.该方法是Image Quilting在三维上的拓展，通过距离能量公式标定几何纹理块之间的相似度，作为邻域匹配的度量单位，以此找到最佳匹配输出；最后采用Graph Cut方法在几何纹理块之间找到一条最佳分割路径，并结合网格编辑技术合成新的几何纹理.该方法具有很好的适用性，但是在模型表面进行合成时耗时较长.2009年，韩建伟等^[14]提出了三维曲面上基于Wang Tiles的几何纹理合成方法.该方法根据给定的几何纹理预计算一组Wang Tiles，然后在不同的Polycube化的目标物体上用其即时生成新的几何纹理.该方法生成的几何纹理中，Wang Tiles可以重用到不同的目标物体上，同时占用的存储空间及计算量更小，速度更快.但是模型参数化和几何纹理映射到模型表面上亦需要大量的人工交互，而且，该方法不适用于结构不规整的几何纹理.2011年，MA等^[15]提出一种基于离散元素的几何纹理合成方法，该方法是由几何纹理数据驱动的.通过对几何纹理样本中每个元素的位置、大小、形状、方向等信息进行编码，对于给定的几何纹理合成空间，该方法采用基于样本的邻域相似性度量和能量优化求解器来合成所需的输出.2012年，MA等^[16]提出了另一种动态元素纹理的几何纹理合成方法，用于合成具有运动效果的几何纹理.由于该类方法是由几何纹理数据驱动的，所以只能针对特定几何纹理合成到特定的输出空间上，具有一定的局限性.2013年，ALMERAJ等^[17]提出了基于块的几何纹理合成方法，该方法仅适用于样本中元素是离散排列的几何纹理合成，并且合成结果容易出现较大裂缝.

上述纹理合成方法均无法避免合成过程中的高计算量、高内存占用和高耗时等问题，因此，有学者提出使用GPU加速几何纹理合成过程，以解决合成过程中产生的海量计算问题，文献[18]提出了一种采用GPU加速网格节点着色合成肝脏体纹理的合成方法.然而，现有几何纹理合成算法都比较复杂，无法很好地将算法并行化，使得目前基于GPU的相关研究仍较少.本文结合文献[13]和^[17]的基本思想，提出了一种基于GPU的几何纹理合成算法.GPU并行多线程技术可加速几何纹理的合成，同时该方法适用于处理多种类型的样本几何纹理，具有较好的通用性.

2 基于GPU的几何纹理合成

现有的基于块的几何纹理合成方法将整个样本几何纹理作为输入，通过求得相似度量能量公式最优解获得纹理块之间的最佳放置位置，进而采用网格编辑技术将纹理块进行融合拼接，以合成大尺寸的几何纹理.其核心是对样本块之间重叠区域的几何相似进行全局搜索，得到最优匹配，但是全局搜索范围受限于样本几何纹理的尺寸.因此，合成结果会存在较大的空隙(离散型几何纹理)，或者出现局部无法对齐(连通型几何纹理)等问题.

本文提出了一种基于GPU的几何纹理合成算法，见图 1，该算法以光栅扫描的方式进行几何纹理合成.通过网格化样本几何纹理进行子块划分，此时样本子块作为合成的基本单位，扩充了搜索匹配的全局空间，进而能够得到更好的合成结果.其子块划分虽然能够得到较好的合成结果，但邻域搜索匹配过程中的计算量将急剧增加，使得几何纹理合成过程耗时更长.因此，引入GPU并行处理技术，将计算量最大的最佳几何子块匹配和几何子块拼接进行并行处理，通过采用多线程并发技术处理这两部分产生的海量计算，从而大幅提高几何纹理的合成速率.

2.1 样本几何纹理预处理

给定一个包围盒大小为l_in×w_in×h_in的样本几何纹理M_in={V_in, F_in}，其中V_in是顶点集合；F_in是面片集合；l_in，w_in和h_in分别表示M_in的长度、宽度和高度(厚度).输出的是一个大小为l_out×w_out×h_out的新的几何纹理M_out，其中l_out，w_out和h_out分别为合成几何纹理的长度、宽度和高度.按输入样本的l_in和w_in将样本几何纹理M_in等分为m×n的网格Grid_in={X_in, Y_in, Cell}，在几何纹理合成过程中就可以不考虑样本几何纹理的高度(厚度，h_in)，其中X_in，Y_in为Cell在Grid_in的位置；Cell={V, F}为网格单元内顶点数据和面片的集合，即包含所有样本几何子块.类似地，合成结果M_out划分为Grid_out={X_out, Y_out, Cell}，并要求M_in与M_out是具有相同大小的Cell.需要指出的是，本文算法要求事先指定M_out的尺寸.

由于对样本几何纹理进行了网格化处理，因此，本文算法只需通过访问Cell的位置坐标即可快速读取样本几何纹理中任意位置的几何子块数据，加快数据存取的速度.另外，根据几何子块分别在Grid_in和Grid_out中的相对空间关系，算得将Grid_in中的Cell平移到Grid_out中指定位置的平移向量T，平移后的Cell表示为Cell′=Cell{V_in·T, F_in}.虽然Cell与Cell′的主要区别在于空间位置，但是Cell′可以和Cell共享顶点、面片和材质等数据，因而不需要给Cell′分配额外的空间保存顶点、面片和材质等数据.同样，可对Grid_out中所有已填充的Cell做相同的处理.通过建立Grid_in和Grid_out之间的顶点、面片和材质等数据的共享机制，将样本数据拷贝至GPU显存后进行数据共享，因此不再需要额外的CPU和GPU的数据交换，从而有效地减少了CPU与GPU之间的数据交互量.

2.2 最佳匹配几何子块查找

以光栅扫描的方式进行几何纹理合成，在样本纹理中首先随机选取一个子块填充到合成输出空间的初始位置，而与其相邻的未填充Cell (即当前待合成区域)定义为种子区域Seed.然后取与Seed相邻的已合成Cell的边界区域作为后续几何结构相似性度量的匹配区域P_out，P_out的宽度在本文所有实验中设置为Cell宽度的10%.将P_out沿水平方向或垂直方向，以P_out的宽度为单位在M_in中进行平移，在样本几何纹理中进行最佳匹配几何子块查找，P_out每次平移后与M_in重叠的区域标记为P_in(t).此时，P_in(t)是与P_out具有相同大小的区域.由于本文仅将P_out的移动范围限定在M_in内，因此需要在M_in中找到具有最优相似度匹配的平移量t，此时的P_in(t)即为最佳匹配的几何子块.在本文方法中，度量P_out与P_in(t)之间几何结构相似程度的公式定义如下：

$ E\left( t \right)=\sum\limits_{v_{\text{out}}^{i}\in {{P}_{\text{out}}}}{E\left( v_{\text{out}}^{i}, {{P}_{\text{in}}}\left( t \right) \right)}, $

(1)

其中：v_outⁱ表示P_out的顶点；P_in(t)表示P_out每次平移后与M_in的重叠区域；E(v_outⁱ, P_in(t))为P_out中各顶点v_outⁱ到P_in(t)中所有面片的最小相似误差，定义如下：

$ E\left( {v_{{\rm{out}}}^i,{P_{{\rm{in}}}}\left( t \right)} \right) = \mathop {\min }\limits_{f_{{\rm{in}}}^j \in {P_{{\rm{in}}}}\left( t \right)} D\left( {v_{{\rm{out}}}^i,f_{{\rm{in}}}^j} \right), $

(2)

其中：f_in^j为P_in(t)的面片；D(v_outⁱ, f_in^j)用于计算P_out中各顶点v_outⁱ到P_in(t)中各面片f_in^j之间的距离(相似误差)，定义如下：

$ D\left( {v_{{\rm{out}}}^i,f_{{\rm{in}}}^j} \right) = \left( {1 + \lambda {\rm{Dist}}\left( {v_{{\rm{out}}}^i,f_{{\rm{in}}}^j} \right)} \right) \times \left( {1 + \left\| {\mathit{\boldsymbol{n}}\left( {v_{{\rm{out}}}^i} \right) - \mathit{\boldsymbol{n}}} \right.\left. {\left( {f_{{\rm{in}}}^j} \right)} \right\|} \right), $

(3)

其中：n(v_outⁱ)是顶点v_outⁱ的法向量；n(f_in^j)是面片f_in^j的法向量；λ是顶点到面片最短距离Dist (·, ·)在相似度误差计算中所占的权重，本文中λ=1.

当P_out与P_in(t)之间的相似度误差E(t)为最小时，认为该重叠区域中的几何结构是最相似的.然后在样本几何纹理中取包含P_in(t)的几何子块m_out(除重叠区域且与当前所选Cell相同的部分)放到Seed所在的位置，合成区域顺延至下一个Cell.

2.3 几何子块拼接

此时，m_out和其周围已合成区域的几何误差最小.然而，重叠区域中仍然存在因几何结构无法对齐导致的错位问题(见图 2(a))，需要对处于重叠区域中的P_out和P_in(t)∈m_out进行几何结构对齐.如果能够在P_in(t)中找到与P_out中的各顶点v_outⁱ具有最小相似误差的面片f_in^j，则能够以此顶点v_outⁱ和对应的面片f_in^j建立几何子块之间的几何结构对齐关系；否则，重叠区域中无法对齐的几何结构部分则需要去除.通过几何结构对齐之后，在P_out和P_in(t)∈m_out之间建立几何结构对齐关系集S={(v_outⁱ, f_in^j)}.然而在实际的纹理空间中，由于P_out和P_in(t)∈m_out的一部分区域相互重叠而其他区域可能是分离的，S中的几何对齐关系不一定都是合理的.因此，本文通过以下3个规则对S进行优化：1)去除v_outⁱ和f_in^j两者法向量反向的对齐关系；2)去除Dist (v_outⁱ, f_in^j)远大于局部面片边长平均值的对齐关系；3)将v_outⁱ所在的面片和f_in^j投影到同一平面上，去除投影不重叠的对齐关系.

在进行几何结构对齐之后，本文采用拉普拉斯网格编辑技术^[20]对几何结构对齐后的P_in(t)∈m_out和P_out进行变形处理，这使得P_in(t)∈m_out和P_out在保持局部几何细节不变的同时，能够解决m_out与其周围已合成区域之间存在的几何结构错位问题，从而产生准确的放置位置(见图 2(b)).在进行变形处理之前，需要1)分别将P_in(t)∈m_out和P_out中所有顶点的空间坐标转换为拉普拉斯坐标；2)计算P_out中各顶点v_outⁱ到与其建立几何对齐关系的面片f_in^j上距离最近的点h_inⁱ，并将它们的中点c_in=(v_outⁱ+h_inⁱ)/2作为形变的基准点.根据式(4)和(5)，通过对其进行二次极小化问题求解，分别计算P_in(t)∈m_out和P_out中所有顶点变形后的顶点坐标：

$ \begin{matrix} {{E}_{\text{out}}}\left( \left\{ {{w}^{i}} \right\} \right)=\sum\limits_{i=1}^{{{N}_{\text{out}}}}{{{\left\| \mathscr{L}\left( {{w}^{i}} \right)-\mathscr{L}\left( v_{\text{out}}^{i} \right) \right\|}^{2}}+} \\ \mu \sum\limits_{i=1}^{m}{{{\left\| {{w}^{i}}-{{c}^{i}} \right\|}^{2}}, } \\ \end{matrix} $

(4)

$ \begin{matrix} {{E}_{\text{in}}}\left( \left\{ {{w}^{i}} \right\} \right)=\sum\limits_{i=1}^{{{N}_{\text{in}}}}{{{\left\| \mathscr{L}\left( {{w}^{i}} \right)-\mathscr{L}\left( v_{\text{in}}^{i} \right) \right\|}^{2}}+} \\ \mu \sum\limits_{i=1}^{m}{{{\left\| {{\alpha }^{i}}v_{\text{in}}^{i, 1}+{{\beta }^{i}}v_{\text{in}}^{i, 2}+{{\gamma }^{i}}v_{\text{in}}^{i, 3}-{{c}^{i}} \right\|}^{2}}, } \\ \end{matrix} $

(5)

其中：$\mathscr{L}$(·)用于计算各顶点的拉普拉斯坐标，详见文献[13]；μ=1为平衡系数；α，β和γ分别为面片f_in^j各顶点对应的权重系数.

实际上，计算P_out和P_in(t)∈m_out的顶点vⁱ形变坐标的主要目是确定顶点在vⁱ和h_inⁱ之间所有可能位置上的形变误差，使其产生最佳形变.此外，通过进一步拆分式(4)和(5)，可用于表示将点vⁱ变形到vⁱ和h_inⁱ之间的某一点p_jⁱ的变形误差，进而将最佳形变坐标的求解转换为求每个顶点的最小变形误差：

$ {{E}_{\text{out}}}\left( {{v}^{i}} \right)={{\left\| \mathscr{L}\left( p_{j}^{i} \right)-\mathscr{L}\left( {{v}^{i}} \right) \right\|}^{2}}+\mu {{\left\| p_{j}^{i}-{{c}^{i}} \right\|}^{2}}, $

(6)

$ {{E}_{\text{in}}}\left( {{v}^{i}} \right)={{\left\| \mathscr{L}\left( p_{j}^{i} \right)-\mathscr{L}\left( {{v}^{i}} \right) \right\|}^{2}}+\mu {{\left\| {{\alpha }^{i}}{{v}^{i, 1}}+{{\beta }^{i}}{{v}^{i, 2}}+{{\gamma }^{i}}{{v}^{i, 3}}-{{c}^{i}} \right\|}^{2}}, $

(7)

$ \left( 1+\left\| v_{\text{out}}^{i}-v_{\text{out}}^{j} \right\| \right)\left( 1+\text{Dist}\left( v_{\text{out}}^{i}, {{C}_{\text{in}}} \right)+\text{Dist}\left( v_{\text{out}}^{j}, {{C}_{\text{in}}} \right) \right). $

(8)

为了准确放置P_in(t)∈m_out并与P_out拼接成一个整体，须保持几何纹理合成结果结构的完整性和连续性，将Graph Cut算法^[19]扩展到网格数据进行处理，即P_in(t)∈m_out和P_out的面片为图的节点；而连接相邻节点之间边的权重为v_outⁱ到P_in(t)的最短距离，可由式(8)计算得到；若2个三角形fⁱ和f^j不共边，则fⁱ和f^j在图中所在结点之间的权重为无穷大.通过在P_out和m_out之间找到一条最优的拼接路径，从而将P_out和m_out进行无缝平滑的合并，如图 2(c)所示.最后，采用YU等^[20]的方法沿着最优路径融合P_out和P_in(t)∈m_out，以获得最终的几何纹理合成结果，如图 2(d)所示.

2.4 GPU并行设计与实现

在进行GPU并行程序设计之前，首先需要明确在串行算法中可被并行化的部分，即原始算法中具有并发性的步骤.判断是否具有并发性的主要依据：多个计算步骤之间互不干扰、没有共享数据操作.

根据上述准则对合成过程进行分析，发现最佳几何子块查找和几何子块拼接的内部步骤可以进行并行化处理，但是最佳子块查找和几何子块拼接存在顺序依赖，因此仅分别对上述2个步骤进行并行化.为方便程序实现，本文将算法中可并行的部分进行分解，分解后的操作简单地采用矩阵表示，其中矩阵中元素的个数即为线程并发量的数目，即为矩阵中每个元素都申请一个线程处理其中的操作.本文对计算量进行了预估，线程并发量以万为单位，下文将进行具体描述.

2.4.1 并行的最佳几何子块查找

在2.2节中，由于P_out中各顶点到P_in(t)所有面片的距离的计算是互不依赖的.因而，可采用GPU并行处理技术对式(3)进行并行化计算.分别构建P_out的顶点矩阵M_V=[v_out¹  v_out² … v_out^m-1  v_out^m]^-1和面片矩阵M_F=[f_in¹  f_in² … f_in^n-1  f_inⁿ]，其中m为P_out中顶点的个数；n为P_in(t)中面片的个数，将已构建的顶点矩阵和面片矩阵乘积后得到相似度误差并行计算矩阵M_E，描述如下：

$ {{\mathit{\boldsymbol{M}}}_{S}}=\left[\begin{matrix} D\left( v_{\rm{out}}^{1}, f_{\rm{in}}^{1} \right) & D\left( v_{\rm{out}}^{1}, f_{\rm{in}}^{2} \right) & \ \cdots & D\left( v_{\rm{out}}^{1}, f_{\rm{in}}^{n} \right) \\ D\left( v_{\rm{out}}^{2}, f_{\rm{in}}^{1} \right) & D\left( v_{\rm{out}}^{2}, f_{\rm{in}}^{2} \right) & \cdots & D\left( v_{\rm{out}}^{2}, f_{\rm{in}}^{n} \right) \\ \vdots & \vdots & {} & \vdots \\ D\left( v_{\rm{out}}^{m-1}, f_{\rm{in}}^{1} \right) & D\left( v_{\rm{out}}^{m-1}, f_{\rm{in}}^{2} \right) & \cdots & D\left( v_{\rm{out}}^{m-1}, f_{\rm{in}}^{n} \right) \\ D\left( v_{\rm{out}}^{m}, f_{\rm{in}}^{1} \right) & D\left( v_{\rm{out}}^{m}, f_{\rm{in}}^{2} \right) & \cdots & D\left( v_{\rm{out}}^{m}, f_{\rm{in}}^{n} \right) \\ \end{matrix} \right]. $

(9)

其中，M_E为m×n的矩阵，M_E中的第i行第j列元素M_E(i, j)=D(v_outⁱ, f_in^j)，即为P_out的顶点v_outⁱ到P_in(t)中面片f_in^j的相似度误差，且M_E中每一行具有最小值的元素即为E(v_outⁱ, P_in(t)).经并行计算后，对M_E每行中具有最小值的元素进行累加求和，即可得到一次几何结构相似度匹配结果E(t).

实际上，每次平移后几何结构相似性度量的计算都是互不依赖的.因此，对几何纹理合成过程中穷尽的最佳匹配几何子块查找，同样可进行并行化处理，并行矩阵M描述为：

$ \mathit{\boldsymbol{M=}}\left( \begin{matrix} E\left( {{t}_{1, 1}} \right) & E\left( {{t}_{1, 2}} \right) & \cdots & E\left( {{t}_{1, k-1}} \right) & E\left( {{t}_{1, k}} \right) \\ E\left( {{t}_{2, 1}} \right) & E\left( {{t}_{2, 2}} \right) & \cdots & E\left( {{t}_{2, k-1}} \right) & E\left( {{t}_{2, k}} \right) \\ \vdots & \vdots & {} & \vdots & \vdots \\ E\left( {{t}_{h-1, 1}} \right) & E\left( {{t}_{h-1, 2}} \right) & \cdots & E\left( {{t}_{h-1, k-1}} \right) & E\left( {{t}_{h-1, k}} \right) \\ E\left( {{t}_{h, 1}} \right) & E\left( {{t}_{h, 2}} \right) & \cdots & E\left( {{t}_{h, k-1}} \right) & E\left( {{t}_{h, k}} \right) \\ \end{matrix} \right), $

(10)

其中：h为M_in网格化后Cell的行数；k为M_in网格化后Cell的列数；M中第i行第j列元素M(i, j)=E(t_{i, j})为P_out与M_in中第i行第j列P_in(t_{i, j})的相似度误差.经并行计算后，M中具有最小值的元素即为最小相似误差，其下标为M_in中对应P_in(t_{i, j})的下标，以此P_in(t_{i, j})为起始边界的样本几何子块，即从M_in中找到的最佳匹配几何子块m_out.

2.4.2 并行的几何子块拼接

在2.3节中，几何结构对齐关系集S的优化过程实质是进行对齐关系判别，这种关系判别相互之间互不影响，因此亦能采用GPU技术进行并行化处理.在相似度误差矩阵M_E的基础上构建几何结构对齐矩阵M_S：

$ {\mathit{\boldsymbol{M}}_S} = {\left[ {\begin{array}{*{20}{c}} {{\rm{ind}}_v^1}&{{\rm{ind}}_v^2}& \cdots &{{\rm{ind}}_v^{m - 1}}&{{\rm{ind}}_v^m}\\ {{\rm{ind}}_f^1}&{{\rm{ind}}_f^2}& \cdots &{{\rm{ind}}_f^{m - 1}}&{{\rm{ind}}_f^m} \end{array}} \right]^{ - 1}}, $

(11)

其中：M_S为m×2价矩阵；m为P_out中顶点的个数；ind_vⁱ是顶点索引，为矩阵M_E行的索引值；ind_fⁱ是面索引，为矩阵M_E每行中D(v_outⁱ, f_in^j)值最小的列的索引值.采用GPU多线程并发处理几何结构对齐关系集S中的对齐关系进行判别，将不符合对齐关系(ind_vⁱ, ind_fⁱ)的面索引ind_fⁱ置零，M_S即最终优化后的几何结构对齐关系集S.

2.3节中几何子块对齐区域变形过程中形变坐标计算时互不依赖，因此式(4)和(5)同样可采用GPU技术进行并行化处理，从而快速计算P_out和P_in(t)∈m_out各个顶点的形变坐标.需要指出的是，在进行变形之前需要先计算每个点的拉普拉斯坐标.另外，P_out和P_in(t)∈m_out的形变坐标并行计算过程完全一致，主要区别在于式(6)和(7)用于计算变形误差值.并行化的变形矩阵M_w描述为：

$ {{\mathit{\boldsymbol{M}}}_{w}}=\left( \begin{matrix} {{E}_{s}}\left( {{c}^{1}}-l \right) & \cdots & {{E}_{s}}\left( {{c}^{1}} \right) & \cdots & {{E}_{s}}\left( {{c}^{1}}+l \right) \\ {{E}_{s}}\left( {{c}^{2}}-l \right) & \cdots & {{E}_{s}}\left( {{c}^{2}} \right) & \cdots & {{E}_{s}}\left( {{c}^{2}}+l \right) \\ \vdots & {} & \vdots & {} & \vdots \\ {{E}_{s}}\left( {{c}^{m-1}}-l \right) & \cdots & {{E}_{s}}\left( {{c}^{m-1}} \right) & \cdots & {{E}_{s}}\left( {{c}^{m-1}}+l \right) \\ {{E}_{s}}\left( {{c}^{m}}-l \right) & \cdots & {{E}_{s}}\left( {{c}^{m}} \right) & \cdots & {{E}_{s}}\left( {{c}^{m}}+l \right) \\ \end{matrix} \right), $

(12)

其中：M_w是m×(2l+1)价矩阵；m表示P_out或P_in(t)∈m_out上顶点的个数；2l表示P_out或P_in(t)∈m_out上顶点中vⁱ到h_inⁱ距离的最大值；E_S表示计算变形误差值的方法，计算P_out变形误差值时E_S=E_out(vⁱ)，而计算P_in(t)∈m_out的变形误差值时E_S=E_in(vⁱ)；M_w中的第i行第j列元素M_w(i, j)=E_s(cⁱ+j-l)，其中(cⁱ+j-l)为拟变形的位置(0 < i≤m, 0≤j≤2l)，表示将vⁱ拟变形至(cⁱ+j-l)处的形变误差.经过并行计算之后，将(cⁱ+j-l)中的j替换为M_w中每行具有最小值元素的列索引值，即所求顶点的最佳形变位置，而对M_w每行中具有最小值的元素进行累计求和，即E_out({wⁱ})或E_in({wⁱ})的结果.由于P_out或P_in(t)∈m_out上的每个顶点vⁱ到与之对应的h_inⁱ的距离不一样，M_w每行的维数将不统一(l的值不唯一)，使得无法构建M_w.因此，需要求得P_out或P_in(t)∈m_out的顶点vⁱ到与之对应的h_inⁱ的距离的最大值，并将P_out或P_in(t)∈m_out所有顶点的变形区间都统一为这个最大值，以此保证M_w的顺利构建.需要指出的是，本文将变形的起点设在vⁱ和h_inⁱ之间的中点c_i，并以c_i为基准点，分别向vⁱ和h_inⁱ移动.另外，待变形区域顶点的拉普拉斯坐标和变形坐标基准点计算的并行化，可参照几何结构对齐关系集优化的并行化处理过程.

最后，通过并行的Graph Cut^[21]在P_out和m_out之间找到最优拼接路径，并采用YU等^[20]的方法将P_out和m_out拼接成一个整体，以此保持几何纹理合成结果结构的完整性和连续性.

3 实验结果与分析

在Windows8.1下用VS2013、OpenGL、CUDA6.5实现了本文的算法，运行环境是2.5 GHz的Intel (R) Core (TM) I5-3230M CPU、4G内存，显卡NVIDIA GT750M，384个流处理器，2 G显存，0.97 GHz GPU，计算能力3.0.分别从几何纹理的合成效果和合成耗时上验证算法的可行性.

图 3给出的是本文方法与现有几何纹理合成方法的比较.其中，第1列为样本几何纹理；第2列为所比较方法合成的几何纹理结果；最后1列为使用本文方法合成的几何纹理结果.由图 3(b1)，文献[17]的合成结果在去除重叠区之后，需要对纹理中的元素进行旋转，这将导致纹理内部元素分布不均匀，元素排列不紧凑，出现比较大的缝隙和空洞，并且局部破坏了样本的几何结构.本文方法元素排列紧凑有致，很好地保存了样本几何的结构特征，合成结果更饱满.文献[16]的合成结果排列紧凑有序，但破坏了样本的结构特征，与样本不具有自相似性；本算法元素的排列亦紧凑有序，并且保存了样本的几何特征.需要指出的是，文献[16-17]仅适用于内部元素离散排列的几何纹理，具有较大的局限性.相较于文献[16-17]，本算法对内部元素离散排列的几何纹理具有更好的合成效果.

图 3中(b3)、(c3)为本文算法和文献[13]合成结果对比图，文献[13]的合成结果中缝合拼接处产生了比较明显的变形，使得合成结果不自然，这主要是因为文献[13]的方法以样本块作为拼接单位，使得全局搜索范围相对较小，同时块与块之间在较大区域范围内进行对齐，在进行样本块对齐和变形时无法取得较好的变形效果，令合成结果不理想.另外，该方法对内部元素排列不规则的样本(枫叶)无法取得较好的合成结果.而本文方法通过对样本几何纹理进行几何子块细分，并以几何子块作为合成的基本单元，扩大了最佳几何子块匹配全局搜索范围，缩小了几何结构对齐的区域范围，从而取得了较好的变形结果，使得合成结果更加平滑自然，合成效果更好.图 4给出了更多的几何纹理合成结果.

本算法在保证合成质量的同时，极大地提高了几何纹理的合成速度.为了对基于GPU加速算法的加速效果有更直观的认识，对最佳几何子块匹配查找加速前后算法的计算复杂度进行了比较分析.假设M_in划分为h行和k列，P_out中有m个顶点，P_in(t)中有n个面片，那么，串行的最佳匹配几何子块查找的计算复杂度为$c=\sum\limits_{1}^{h}{\sum\limits_{1}^{k}{m\times n}} $.而并行化处理后的最佳匹配几何子块查找的计算复杂度将变为1.这主要是因为，通过采用GPU多线程并发技术处理查找过程中的海量运算，能够大大降低最佳匹配几何子块查找的计算复杂度，在极大地提高查找速度的同时，有效缩短了几何纹理合成的耗时.

表 1对本文方法和文献[13, 16-17]中的方法在合成耗时上进行了统计.由表 1可知，上述方法合成一个2倍大小的几何纹理所花的时间基本相当，约为5 min，而本算法合成同样大小的几何纹理只需约20 s，合成速度提高了近20倍.如果采用更高性能的硬件，合成速度将会进一步提高.同时，实验发现，几何元素离散分布且排列规则有序的样本几何纹理，其合成速度要明显快于结构复杂、元素排列无规则的编织网样本几何纹理.表 1给出了枫叶样本的合成时间，其面片数和顶点数明显多于编织网样本，但花费的时间却要少于编织网样本.这是因为在进行样本块拼接时，花费了大量时间对结构复杂、元素排列无规则的编织网样本进行几何块或几何子块变形.为了验证上述结论，对几何纹理合成过程各个阶段的耗时进行了统计(见表 2).

从表 2中给出的统计数据可以看出，对于顶点和面片数据较多的样本(枫叶)，在几何子块匹配阶段耗时较多；内部元素具有连通性的几何纹理样本(编织环、编织网)在几何子块拼接阶段的耗时要大于内部元素离散分布的几何纹理样本(枫叶)；内部元素连通且排序不规则的几何样本纹理(编织网)的几何拼接耗时要大于内部元素连通且规则排列的样本几何纹理(编织环).由此可见，在合成过程中，合成结果和样本几何纹理共享样本的顶点数据，因不需要增加保存顶点数据的额外内存，不仅降低了内存占用，而且降低了CPU并加快了GPU的数据交互量.另外，GPU并行处理可有效加快几何纹理合成过程，提高合成速度.因此，相较于现有的基于块的几何纹理合成方法，本算法在保证合成质量的同时，大幅度提高了几何纹理的合成速度.

4 总结

现有方法在几何纹理合成过程中存在高计算量、高存储占用和高耗时等问题.本文提出了一种基于GPU加速的几何纹理合成方法.该算法可重用样本顶点数据的数据结构，优化存储以降低内存的占用率.通过对样本几何纹理进行子块划分，增大了邻域全局搜索匹配的范围；通过采用GPU多线程并发技术，设计了GPU并行加速算法，以加快几何纹理合成速度，最终能够快速高效生成新的且任意大尺寸的几何纹理.本文主要研究由小尺寸几何纹理样本生成大尺寸的几何纹理，将生成的大尺寸几何纹理快速高效地映射到具体的模型表面，有待进一步讨论.