一种基于渲染的快速体数据骨架提取方法

一种基于渲染的快速体数据骨架提取方法
【专利摘要】本发明涉及一种基于渲染的快速体数据骨架提取方法，包括以下步骤：在上下前后左右6个视点下，用RayCasting直接体绘制算法对三维体数据进行渲染，得到6个2D图像；对所有直接体绘制获得的2D图像提取其轮廓及部分关键特征线；将获得的2D轮廓线及关键特征线反投影回3D空间，获得3D候选曲线；对空间采样，计算空间采样点的曲线密度，构建曲线密度场；获得曲线密度场后，用三维Sobel算子采样点的梯度；根据采样点的梯度值，用最小二乘方法计算采样点的曲线方向；将采样点的曲线方向赋值给其所对应的三维体数据的体素；根据体素的方向，构建包围盒，计算包围盒体素平均值获得表征节点；连接表征节点并平滑曲线最终得到三维体数据骨架。
【专利说明】
【技术领域】
[0001] 本发明属于骨架提取及体绘制【技术领域】，具体涉及基于渲染的快速体数据骨架提 取方法。 一种基于渲染的快速体数据骨架提取方法

【背景技术】
[0002] 骨架作为形状表示的一种有效形式在三维模型的各个研究领域中被广泛采用。 Bluml967年给出了骨架的最初定义：骨架（中轴）是模型内部各个最大内切球中心的集 合。它还有一个grassfire的模拟定义，从模型表面开始点火，各个方向上的火的相遇点所 构成的集合。因为模型的骨架很好的保留了模型的拓扑连接性及其几何形状信息，所以经 常被用于碰撞检测、三维动画、模型渲染、模型表面重建、模型检索、三维形状识别与检索、 模型编辑、简化以及三维分割等应用中，也有研究人员采用骨架为模型的分解做矫正。不同 的应用，对于骨架应该保存的信息要求不完全相同，故而抽取思路也不完全相同。
[0003] 对三维模型骨架的研究由来已久，出现过很多方法，有的是源于对二维图像的扩 展，有的是针对三维模型提出的，大体上来说，面向三维网格模型或者三维体素模型的骨架 提取方法主要包括：细化、距离变换、分解等。
[0004] 提取三维体素模型骨架的拓扑细化方法的基本思想是：由外向内，层层剥离不改 变模型拓扑特性的体素来实现骨架抽取，它检测物体数据的拓扑相关性，而不是形状的几 何特征。这种方法在保持物体拓扑性质的同时，在离散域中反复的逐层剥离去物体的表面， 其中判断一个体素是否需要剥离的检测过程比较耗时，这种发方法需要完全分割好的、紧 凑的并且相互连接的物体模型。Bertrand提出一种简单的并行细化方法抽取骨架。Ju采 用边细化边去除冗余体素的方法抽取骨架。Melada提出结合欧式距离测度与细化过程的骨 架抽取方法。但这些方法产生的骨架包含杂枝，甚至破坏了骨架的连续性。
[0005] 基于距离变换的骨架提取方法：它一般的计算对象要求也必须是体素表示的模 型，通过计算每个体素的距离来求取模型的脊点。Zhou提出了边界距离变换和种子点距离 变换，采用体素编码的方式抽取并优化骨架结构，但骨架在分支处的中心性较差。Wan的方 法采用精确的边界距离变换，建立极小扩张树提取骨架，从而使得到的骨架点尽量不偏离 中心，但会出现中心线的"扰动"现象。基于Level Set的骨架提取方法稳定性较好，也具有 较高的拓扑无关性，能够有效消除尖点、骨架断裂等问题，但是即便采用了 Fast Marching Method策略，其计算复杂度依然很高。Fan引入snake模型调整骨架点的位置，改善了骨架 点的中心性，但复杂度也同时增大。
[0006] 相对于基于细化和距离变换两种方法而言，基于形状分解或者分割的三维骨架提 取方法将骨架提取和形状分割结合起来，骨架提取与形状分解二者相辅相成具有显著优 势。基于拓扑连接信息，先行对三维模型进行分解；基于分解结果可以降低骨架提取的复杂 度、提高骨架的精确度。按照拓扑结构，正确分解三维模型能够保证快速准确地获取骨架， 反之，骨架作为拓扑结构的参考能够保证模型分解的准确性。基于分解的方法源于对形体 的认识，Hoffman和Richards把形体自然分解成基本部件，把分解用于分类。Rom等提出利 用中轴逐步分解二维形体的方法，先分解形体，在提取分解区域的中轴，而中轴是二维形体 的自然表示，但是这些方法难以提取中心良好的三维立体形体的分级骨架。Minimal Rule Lien等提出基于模型近似图面体分解，通过计算三维模型顶点的凹凸性展开分解，并通过 迭代生成骨架，但得到的骨架中心性较差，也难以建立层次关系。
[0007] 基于Reeb图思想。该类算法首先在模型上定义一个连续函数，计算每个模型顶点 的函数值，将具有相同函数值的顶点聚合成一个顶点，得到模型的骨架。Hilaga等人提出了 基于Reeb图思想的多分辨率Reeb图（MRG)，将Reeb图用于骨架提取，并将映射函数定义为 顶点到整个模型表面所有顶点的最短距离与面积之乘积的和。Julien Tierny等人使用了 更加合理的标量函数计算方法。
[0008] 亦有学者将二维图像领域中的Voronoi图技术引入到三维骨架抽取中，Dey等人 提出过一种利用Voronoi图直接近似中轴的算法。因为对三维模型而言只有部分Voronoi 顶点能够汇聚成骨架，所以他通过定义角度和比例这两个与大小、比重都无关的筛选标准 来实现它的中轴近似算法，并证明该方法能够保证收敛。Voronoi图方法对于复杂的三维模 型，理论上可以实现，但计算量庞大。Voronoi图方法多用于生成多尺度骨架，但对于边界噪 声非常敏感，导致Voronoi图过于稠密，需要额外进行剪枝处理，而且规则复杂，导致该方 法实用性不光。
[0009] 体绘制技术是一种重要的科学可视化技术，采用该技术，不经过复杂的特征算法， 就能够交互地浏览数据内部结构特征。该项技术广泛应用于医学、气象和地质等领域。目前 直接体绘制算法主要有Ray Casting算法、Shear-Warp算法、足迹算法、Texture Mapping 算法等。LevoyM等将体数据递归均勻子分成一个八叉树结构，在光线投射算法中遍历八叉 树来跳过包含当前采样点的最大空区域，以此达到加速的效果。
[0010] 综上所述，现有的技术存在着计算复杂度大，计算时间长，骨架连续性、中心性差、 有扰动，骨架对噪声敏感等缺点，本发明提供一种快速、精确、噪声扰动不敏感的骨架提取 方法。


【发明内容】

[0011] 本发明要解决的技术问题是：克服了现有技术计算时间长、骨架中心性差、对噪声 敏感的不足，提供一种提取体数据骨架的方法，能够快速地提取出三维体数据模型的骨架， 并且能得到较好的效果。
[0012] 本发明解决上述技术问题的技术方案为：基于渲染的快速体数据骨架提取方法， 利用NVIDIA GPU硬件的计算能力，提高了计算效率，该方法主要利用具有视觉意义的2D轮 廓线及部分关键特征线，并根据能量最小化函数，计算对应的三维特征线，对空间采样后构 建曲线密度场进而求取采样点的曲线方向，并将采样点方向赋值给相应体素以构建包围盒 得到表征节点，连接这些表征节点并平滑曲线得到体数据骨架。
[0013] 如图1所示，具体实现如下步骤：
[0014] 1) 3D候选曲线的提取
[0015] 步骤（1)、在X、Y、Z3个坐标轴的正负方向（上下前后左右）上设置六个视点，用 RayCasting直接体绘制算法对三维体数据进行渲染，得到6幅2D图像；该方法是对图像的 每一个像素，沿固定方向（通常是视线方向）发射一条光线，光线穿越体数据的整个图像序 列，并在这个过程中，对图像序列进行采样获取颜色信息，同时依据光线吸收模型将颜色值 进行累加，直至光线穿越整个图像序列，最后得到的颜色值就是渲染图像的颜色。该过程可 以在GPU上实现，利用CUDA(计算统一设备架构）的硬件加速，大大提高渲染速度。
[0016] 步骤（2)、对步骤⑴得到的6幅2D图像，使用轮廓检测算法提取2D图像的轮廓， 直接提取的轮廓需要收缩以获得有效轮廓线，另外为表示更多细节需要使用Canny算子对 2D图像进行部分关键特征线的提取。Canny算子是一阶算子，其边缘检测可以分为三个步 骤：①高斯平滑函数。目的是为了平滑以消除噪声；②一阶差分卷积模板。目的是为了达到 边缘增强；③非极大值抑制（匪S)。目的是保留梯度方向上的最大值。其贡献在于给出了 一种寻找最优算子的思路，即确立了最优准则。
[0017] 步骤（3)、2D曲线在二维上是具有明显视觉意义的线，其表征了体数据的拓扑和 几何结构。为了更准确的表示体数据的结构信息将这些2D曲线反投回三维空间，找到每 条2D曲线对应3D候选曲线。然而2D曲线与3D曲线是一对多的关系，为了求解一条最有 3D对应曲线，定义了一个能量最小化函数，因此求解3D曲线问题转化为一个全局最小化问 题。3D曲线上点有位置信息，不透明度信息，颜色信息，都是在最小化问题需要考虑的因子。 由于求解的最优3D曲线是曲线上的相连点之间能量最小，因此能量最小化问题可以转化 为最短路径问题，利用动态规划算法可以求解该问题，并且借助GPU并行计算能力，大大提 高计算效率。在GPU上每个线程负责一条3D线的计算过程。
[0018] 2)曲线方向场生成
[0019] 步骤（4)、将2D曲线反投射回3D空间后，对3D空间进行采样，求得的3D曲线在采 样空间内会经过若干采样点，利用卷积曲面的方法，对经过采样点一定范围的曲线进行积 分求得采样点的曲线密度，以构建曲线密度场。曲线密度场中每个采样点的密度值是衡量 其空间范围内的曲线的多和少，由于很多采样点周围并没有曲线经过，因此其曲线密度值 为零，对这些点则可以通过预计算剔除，只需计算剩余采样点的曲线密度值即可。该过程是 在GPU上实现，在预计算阶段，每个线程控制一条曲线，将其经过范围的采样点标记，没被 标记的点则是被剔除的点。在采样点曲线密度求解阶段，每个线程控制一个剩余采样点的 计算，大大提高了计算速度。
[0020] 步骤（5)、构建完曲线密度场后，用Sobel算子求得采样点的梯度，由于步骤（4)已 经剔除了值为零的采样点，所以此过程只需计算非零采样点的梯度值。梯度表征了曲线密 度场的变化率最大的方向，因此采样点的梯度大致垂直曲线方向，然后在利用最小二乘方 法可以在梯度基础上求得采样点的曲线方向，该方向与3D曲线的走势基本一致。该过程是 在GPU上实现，每个线程控制一个非零采样点的计算过程，大大提高了计算速度。
[0021] 3)体数据模型骨架的生成
[0022] 步骤（6)、生成曲线方向场后，由于空间采样点与体素在空间位置上有对应关系， 将非零采样点的方向赋值到其对应的体素上，以具有方向的体素的位置和方向为依据构建 包围盒，包围盒的构建过程是以连通性为依据，在体素方向上上下一定距离L处构建两个 个垂直于体素方向的平面，两个平面之间所有联通的体素归为一组，然后将包围盒内所有 体素求一个平均值并作归中处理，得到表征节点。
[0023] 步骤（7)、以步骤（6)得到的表征节点为节点，以表征节点之间的距离为权，构建 无向完全图G，并用Prim算法求解图G的最小生成森林，检测出森林中树的叶节点和分叉 点，从而进一步生成Polyline曲线。对于得到的Polyline曲线，由于是把各个表征节点点 进行连接得到，这样的曲线是不够圆滑的，采用B样条拟合方法，Polyline上的点作为控制 点，计算出更多的插值点，求得B样条曲线，使得曲线更光滑，得到最终体数据模型骨架。
[0024] 本发明与现有技术相比的优点在于：本发明能够快速地提取出三维体数据模型的 骨架并且能得到较好的效果。本发明主要有两点贡献：第一，给出了一种基于渲染求解曲线 方向场的方法，该方法是从RayCasting直接体绘制算法渲染的2D图像出发，提取其轮廓和 部分关键特征线，进而投射回3D空间，并根据3D候选曲线利用卷积曲面方法构建曲线密度 场，在此基础上进而生成了曲线方向场，所有过程均在GPU上并行进行，具有快速高效并较 准确的优点。第二，将曲线方向赋值到体素上，并有体素位置和方向构建包围盒求得表征节 点，连接这些节点并利用B样条进行平滑拟合得到体数据模型的骨架，该过程降低了计算 量使得骨架提取过程快速而准确。

【专利附图】

【附图说明】
[0025] 图1为本发明整体流程图；
[0026] 图2为RayCasting直接体绘制原理图；
[0027] 图3为射线穿越体数据示意图；
[0028] 图4为体绘制渲染的2D图像效果图；
[0029] 图5为提取的2D曲线效果图；
[0030] 图6为2D曲线反投射回3D空间示意图；
[0031] 图7为3D候选曲线效果图；
[0032] 图8为计算采样点曲线密度场示意图；
[0033] 图9为采样空间曲线密度场效果图；
[0034] 图10为采样点梯度效果图；
[0035] 图11为采样点曲线方向效果图；
[0036] 图12为体素方向效果图；
[0037] 图13为表征节点效果图；
[0038] 图14为体数据模型骨架效果图。

【具体实施方式】
[0039] 下面结合附图以及本发明的【具体实施方式】进一步说明本发明。
[0040] 如图1所示，给定体数据由体素组成，本发明针对三维体数据的骨架提取包括以 下步骤：
[0041] 步骤（1)如图2(只设了 3个）在X、Y、Z3个坐标轴的正负方向（上下前后左右） 上设置六个相机，用RayCasting直接体绘制算法对三维体数据进行渲染，得到6幅2D图 像；基于图像空间的算法，对图像的每一个像素，沿固定方向（通常是视线方向）发射一条 光线，光线穿越体数据的整个图像序列，并在这个过程中，对图像序列进行采样获取颜色信 息，同时依据光线吸收模型将颜色值进行累加，直至光线穿越整个图像序列，最后得到的颜 色值就是渲染图像的颜色。采用了从前面到背面进行采样合成，其中颜色和不透明度的累 加公式如公式（1):
[0042] Cj = (l-α 卜广)。，。-/ (1)
[0043] α J = (1-α 卜广）α ηΛ
[0044] 其中，Q和α i分别是在体纹理上采样所得到的颜色值和不透明度，其实也就是体 素中蕴含的数据；C/和α /表示累加的颜色值和不透明度。i表示当前步骤，i-Ι表示前 一步骤。
[0045] 如图3所示，假定光线从f点投射到立方体中，并从1点投出，在立方体中穿越的 距离为m。当光线从f点投射到立方体中，穿越距离为n (n〈m)时进行采样，则存在公式（2):
[0046] t = tstart+d* Δ (2)
[0047] 其中tstart表示立方体表面被投射点f的体纹理坐标；d表示投射方向；Λ表示采 样间隔，随着η的增加而递增；t为求得的采样纹理坐标。通过求得的采样纹理坐标就可以 在体纹理上查询体素数据。直到n>m，或者透明度累加超过1，一条射线的采样过程才结束。
[0048] 该过程可以在GPU上实现，利用CUDA (计算统一设备架构）的硬件加速，大大提高 渲染速度。得到6个2D渲染结果，如图4。
[0049] 步骤（2)、对步骤⑴得到的6幅2D图像，使用轮廓检测算法提取2D图像的轮廓， 直接提取的轮廓需要收缩以获得有效轮廓线，另外为表示更多细节需要使用Canny算子对 2D图像进行部分关键特征线的提取。Canny算子是一阶算子，其边缘检测可以分为三个步 骤：①高斯平滑函数。目的是为了平滑以消除噪声；②一阶差分卷积模板。目的是为了达到 边缘增强；③非极大值抑制（匪S)。将提取的轮廓线与用Canny提取出特征线融合得到效 果图图5,图5中同一颜色的线代表为一 2D曲线。
[0050] 步骤（3)、2D曲线在二维上是具有明显视觉意义的线，其表征了体数据的拓扑和 几何结构。为了更准确的表示体数据的结构信息，将这些2D曲线反投回三维空间，如图6其 中由空心的点连接成的曲线是一条2D曲线，需要找到2D曲线对应3D候选曲线，在图6中由 实线连接黑色实心点而成的曲线是得到的3D曲线。由图6可知一个2D点反投射回3D空间 对应多个3D点，即图6中一个空心点对应多个实心点。定义了一个能量最小化函数（3):

【权利要求】
1. 一种基于渲染的快速体数据骨架提取方法，其特征在于包括如下步骤： (1) 在X、Y、Z，3个坐标轴的正负方向即上下前后左右上设置六个视点，用RayCasting 直接体绘制算法对三维体数据进行渲染，得到6个2D图像； (2) 利用轮廓线检测方法对得到的6个2D图像进行轮廓线的提取时，直接提取的轮廓 需要收缩以获得有效轮廓线，并用Canny算子提取所有2D图像的部分关键特征线，将轮廓 线与部分关键特征线融合以获得有效2D曲线，这些2D曲线体现了体数据的几何及拓扑结 构；所述部分关键特征线即关节部位的特征线； (3) 对步骤（2)提取的2D曲线反投射回三维，由于每条2D曲线与3D曲线是一对多的 关系，故定义一个能量最小化函数，利用动态规划算法求解最优解，得到2D曲线所对应的 最优3D曲线；该过程是在GPU上实现，每个线程控制一条曲线的求解过程，大大提高了求解 速度；该过程是在GPU上实现，每个线程控制一条曲线的求解过程，大大提高了求解速度； (4) 对空间进行采样，求得的3D曲线在采样空间内会经过若干采样点，采用卷积曲面 的方法，对经过采样点一定范围的曲线进行积分求得采样点的曲线密度，以构建曲线密度 场；由于很多采样点周围并没有曲线经过，因此其曲线密度值为零，对这些点则通过预计算 剔除，只需计算剩余采样点的曲线密度值即可；该过程是在GPU上实现，在预计算阶段，每 个线程控制一条曲线，将其经过范围的采样点标记，没被标记的点则是被剔除的点；在采样 点曲线密度求解阶段，每个线程控制一个剩余采样点的计算，大大提高了计算速度； (5) 获得曲线密度场后，用Sobel算子求得采样点的梯度，由于步骤（4)已经剔除了值 为零的采样点，所以此过程只需计算非零采样点的梯度值；梯度表征了场的变化率最大的 方向，因此采样点的梯度大致垂直曲线方向，利用最小二乘方法可以在梯度基础上求得采 样点的曲线方向，该过程是在GPU上实现，每个线程控制一个非零采样点的计算过程，大大 提高了计算速度； (6) 由于3D曲线表征了体数据的几何和拓扑结构，求得采样点的曲线方向后，将采样 点的曲线方向赋值给其所对应的体素，以体素的位置和方向为依据构建包围盒，包围盒的 构建过程是以连通性为依据，将包围盒内所有体素求平均值，得到表征节点； (7) 用Prim算法将所用表征节点用最小生成树连接成Polyline曲线，最后用B样条将 求得的Polyline曲线平滑拟合，得到最终体数据骨架。
2. 根据权利要求1所述的基于渲染的快速体数据骨架提取方法，其特征在于：所述步 骤（1)包括：从X，Y，Z正负轴方向设置6个相机，用RayCasting直接体绘制算法对三维体 数据进行渲染，得到6个2D图像，该过程是在GPU上实现，利用计算统一设备架构CUDA的 硬件加速，大大提高渲染速度，并利用PBO, FB0技术渲染到纹理，并保存到2D渲染结果。
3. 根据权利要求1所述的基于渲染的快速体数据骨架提取方法，其特征在于：所述步 骤⑵具体实现为：以这些表征节点为节点，以表征节点之间的距离为权，构建无向完全图 G，并用Prim算法求解图G的最小生成森林，检测出森林中树的叶节点和分叉点，从而进一 步生成Polyline曲线。采用B样条拟合方法，Polyline上的点作为控制点，计算出更多的 插值点，求得B样条曲线，得到更加的光滑体数据骨架。
【文档编号】G06T15/00GK104156997SQ201410363994
【公开日】2014年11月19日 申请日期:2014年7月28日 优先权日:2014年7月28日
【发明者】王莉莉, 戚庆林, 侯飞, 李帅, 郝爱民 申请人:北京航空航天大学