一种基于OpenCL的三维体数据可视化并行渲染方法与流程

本发明属于图像渲染技术领域，具体涉及一种基于opencl的三维体数据可视化并行渲染方法。

背景技术：

作为科学计算可视化的重要组成部分，体绘制在超声医学成像中发挥着重要的作用。体绘制以计算机图形学、图像处理、计算机视觉等学科为基础，将采集到的三维离散数据以二维图像的方式显示在屏幕上。操作者通过操作鼠标可以达到对成像平面进行旋转和缩放的效果，屏幕上的图像随着鼠标的运动做相应的变换，来完成对三维空间信息的二维图像渲染。借助于opencl并行编程技术，可以实现对三维数据的实时可视化。

图形处理器gpu(graphicsprocessingunit)是显卡的计算核心。与cpu(centralprocessingunit)所不同的是gpu主要用于图形渲染和通用计算。随着gpu并行处理能力的提升、存储能力的加强，基于gpu的光线投射技术己成为当前流行的直接体绘制技术，广泛应用于体数据可视化领域，大大提高了算法的执行效率，改善了用户的交互体验。

为了增加成像的真实感，可以引入光照效果，这就涉及到如何计算体数据中采样点的法向量的问题。由于人体组织对超声激励脉冲的随机漫反射，医学超声设备接收到的超声回波信号中会包含斑点噪声。如果直接采用传统的有限差分方法或简单的sobel算子计算体数据中采样点位置处的梯度，很大程度上会受到斑点噪声的影响，梯度的计算值不准确，在光照的作用下导致成像目标表面呈现出凹凸不平的效果。

技术实现要素：

针对现有技术中存在的以上问题，本发明提供了一种能够高效利用计算资源，提升渲染效果的基于opencl的三维体数据可视化并行渲染方法。

为了达到上述发明目的，本发明采用的技术方案为：

一种基于opencl的三维体数据可视化并行渲染方法，包括以下步骤：

s1、读取用于图像渲染的三维体数据；

s2、对三维体数据进行坐标归一化处理；

s3、对视点进行坐标变换；

s4、基于光线投射算法，对成像平面的每一个像素点引出一条射线，采用opencl并行计算射线与体数据的交点；

s5、选取衰减函数作为自适应采样方法的采样密度函数，在射线上非均匀采样得到采样点；

s6、采用三维b样条函数对采样点进行平滑插值和梯度计算，得到对应三维体数据的采样点位置处的取值和梯度值；

s7、采用光照模型计算采样点的颜色值和阻光度并进行累加求和，得到成像平面上对应像素点的取值进行三维可视化渲染。

进一步地，所述步骤s2具体为：

设定三维体数据的长、宽、高分别为h，w，d，其中最大值为max，将三维体数据的每一个维度坐标限定为(-h/max，h/max)，(-w/max，w/max)，(-d/max，d/max)。

进一步地，所述步骤s3具体为：

对视点进行逆向的旋转和平移操作，将包含平移和旋转操作的变换矩阵表示为：

其中，平移变换矩阵为

绕x轴旋转矩阵为

绕y轴旋转矩阵为

绕z轴旋转矩阵为

得到逆变换矩阵表示为：

其中，tx，ty，tz分别表示x，y，z轴的平移距离，分别表示绕x，y，z轴的旋转角度，m0至m11表示变换矩阵系数。

进一步地，所述步骤s4具体为：

设定视点坐标为引出的视线方向为将视线方程表示为三维体数据平面的法向量为平面上的点为将平面方程表示为得到射线与体数据的交线段表示为：

进一步地，所述步骤s5中选取衰减函数y＝e^-αx作为自适应采样方法的采样密度函数，其中α表示衰减系数，x表示距离，在射线上非均匀采样的间距为：

其中，a为衰减常数，δ为采样点密度，i为采样点序号。

进一步地，所述步骤s6中采用三维b样条函数表示为：

其中，ni，p表示b样条基函数，p，q，r分别表示每个维度的b样条基函数的幂次，pi，j，k表示三个维度的控制点，u，v，ω分别表示x，y，z轴的偏移量。

进一步地，所述三维b样条一阶、二阶、三阶函数控制点的矩阵分别表示为：

梯度矩阵表示为

本发明具有以下有益效果：

(1)本发明通过采用衰减函数作为自适应采样方法的采样密度函数进行非均匀采样，解决了传统光线投射算法采用的均匀采样造成的计算量过大、性能需求过高的问题，进而提高了渲染执行效率，能够更高效的运用计算资源；

(2)本发明通过采用三维b样条函数对采样点进行平滑插值和梯度计算，解决了传统的有限差分方法和sobel算子容易受斑点噪声的影响，不能保证一阶导数的连续性，影响光照成像效果的问题，进而在保证较低计算复杂度的前提下，提供了更大的自由度，能够得到更平滑的渲染结果。

附图说明

图1是本发明基于opencl的三维体数据可视化并行渲染方法流程图；

图2是本发明实施例中衰减函数的采样曲线示意图；

图3是本发明实施例中固定衰减常数的采样曲线示意图；

图4是本发明实施例中固定采样密度的采样曲线示意图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。

如图1所示，本发明实施例提供了一种基于opencl的三维体数据可视化并行渲染方法，包括以下步骤s1至s7：

s1、读取用于图像渲染的三维体数据；

在本实施例中，本发明首先获取用于图像渲染的三维体数据，不同应用场景获取的三维体数据的长、宽、高并不是固定的值，会随着深度、角度、采样率等参数的变化而变化，体数据在世界坐标系中的坐标也会随之发生变化。

s2、对三维体数据进行坐标归一化处理；

在本实施例中，由于体数据在世界坐标系中的坐标会发生变化，导致不能有效地确定视点的初始位置。

例如，体数据的中心置于世界坐标系的原点，x轴指向水平右手方向，y轴指向正上方，x-y平面位于计算机屏幕，z轴指向屏幕正前方。如果体数据的大小为256x256x256，那么其在每个坐标轴占据的空间为(-128～+128)的范围。如果视点位于屏幕正前方，那么视点的坐标可以设定为(x，y，z)＝(0，0，150)。如果体数据的大小为64x64x64，那么视点的坐标可以设定为(x，y，z)＝(0，0，50)。视点需要适应体数据的大小而改变，无法事先预设。当视点距离体数据过近时，类似焦距过小，只能对体数据的局部进行渲染，造成图像偏大；当视点距离体数据过远时，类似焦距过大，造成成像偏小。

基于上述原因，本发明对三维体数据进行坐标归一化处理，从而便于确定视点的初始位置。

对三维体数据进行坐标归一化处理具体为：

设定三维体数据的长、宽、高分别为h，w，d，其中最大值为max，将三维体数据的每一个维度坐标限定为(-h/max，h/max)，(-w/max，w/max)，(-d/max，d/max)，从而实现将三维体数据的每一个维度坐标限定在(-1，-1，-1)到(+1，+1，+1)的空间范围内。

s3、对视点进行坐标变换；

在本实施例中，步骤s2坐标归一化后的体数据被当作一个立方体位于世界坐标系的中心，即作为原点位置。根据步骤s2中的坐标系设置，视点的初始坐标为(x，y，z)＝(0，0，4)，相应的成像平面的初始位置(x，y，z)＝(0，0，-2)。

在与用户的交互过程中，为了方便计算，本发明并不是对体数据立方体进行旋转和平移操作，而是对视点和成像平面进行逆向的旋转和平移操作，将包含平移和旋转操作的变换矩阵表示为：

其中，平移变换矩阵为

绕x轴旋转矩阵为

绕y轴旋转矩阵为

绕z轴旋转矩阵为

视点坐标通过变换矩阵的处理，可以独立地绕x轴，y轴和z轴旋转。得到逆变换矩阵表示为：

其中，tx，ty，tz分别表示x，y，z轴的平移距离，分别表示绕x，y，z轴的旋转角度，m0至m11表示变换矩阵系数，即为平移变换矩阵、x轴旋转矩阵、y轴旋转矩阵、z轴旋转矩阵合并为总的变换矩阵的各项系数。

本发明通过步骤s2和s3可以确定三维体数据的空间归一化坐标，视点的坐标，以及成像平面的坐标。

s4、基于光线投射算法，对成像平面的每一个像素点引出一条射线，采用opencl并行计算射线与体数据的交点；

在本实施例中，本发明采用光线投射算法，从视点投射光线到成像平面的每一个像素点，每一条采样光线穿过三维体数据场，沿着这条射线选择多个采样点，通过插值方法计算出每个采样点的颜色值和阻光度，按照从前到后的顺序进行采样点的颜色值和阻光度的合成，直到阻光度累加至1，从而得到二维投影图像。

光线投射算法具有很高的数据并行度，从视点到成像平面的每一个像素点投射的光线的计算是相对独立的，利用基于opencl的并行编程技术可以实现三维体数据的实时渲染。

本发明从视点发出一条射线与体数据立方体相交，到达成像平面上的一个像素点，交点连接的线段为交线段。

设定视点坐标为引出的视线方向为将视线方程表示为三维体数据平面的法向量为平面上的点为将平面方程表示为l与立方体在x，y，z三个坐标轴的取值范围有关，对于归一化的坐标，l∈[-1，1]。

视线与立方体的交点位于表平面上，得到射线与体数据的交线段表示为：

通过求解参数t的取值范围即可得到交线段，进而求得视点与成像平面像素点之间的连线穿过三维体数据的前后交点。

s5、选取衰减函数作为自适应采样方法的采样密度函数，在射线上非均匀采样得到采样点；

在本实施例中，传统的均匀化采样方法采用了恒定的步长，即每两个相邻的采样点之间的间隔是固定不变的。设视点到第i个采样点的距离向量为pi，视点到目标像素的方向向量为d，那么在传统的均匀采样中，相邻两个采样点满足：pi＝pi-1+λd，其中λ是两个采样点之间的间隔距离，这个值在均匀采样中是固定常数。

然而在实际的应用场景下，较远距离的采样点并不需要和较近距离的采样点采用相同的采样间隔。较远的距离，采样点之间的间隔应当比较近采样点之间的间隔更大，这样才能保证近距离的采样分辨率高，保证距离视点较近的部分渲染完全，然而相对较远距离的，相对不重要的采样点之间采样分辨率相对较低。

因此，本发明采用自适应采样方法，并选取衰减函数y＝e^-αx作为自适应采样方法的采样密度函数，其中α表示衰减系数，x表示距离，在步骤s4得到的交线段上进行非均匀采样，可以从距离视点较近的交点开始采样，到距离视点较远的交点结束，也可以反向进行采样。

如图2所示，当采样密度δy均匀分布时，采样点密度δx的分布满足非均匀采样的要求，即靠近原点的位置采样间距小，远离原点的位置采样间距大。

衰减函数的逆函数可以表示为：

共n+1个采样点，相应的x的取值为采样点的间距λi，共有n个采样点间的间距。

非均匀采样的间距λi为：

其中，a为衰减常数，δ为采样点密度，i为采样点序号。

对体数据坐标进行归一化处理后，上式中的采样点间距λ满足约束条件：0＜λ，衰减常数a也需要满足约束条件0＜a＜1，得到采样密度δ的取值范围为：

通过设置合适的采样密度δ和衰减常数a，可以得到非均匀采样中的采样间距λ随采样点变化的曲线，其随着与视点之间的距离逐渐变大，满足非均匀采样的要求。在在非均匀采样中，本发明设定最大采样点的个数为100个点，则采样密度需满足δ＜0.01。这里选择δ＝0.005，设定初始采样间距λ0＝0.02，从而得到a＝0.25。

固定衰减常数a＝0.35，不同的采样密度δ条件下采样间距λ随采样点的变化曲线如图3所示。

固定采样密度δ＝8×10^-3，不同衰减常数a的条件下采样间距λ随采样点的变化曲线如图4所示。

根据图3和图4可以看出，选择参数a＝0.35和δ＝8×10^-3可以达到最优的采样效果。

s6、采用三维b样条函数对采样点进行平滑插值和梯度计算，得到对应三维体数据的采样点位置处的取值和梯度值；

在本实施例中，通过步骤s5确定了当前采样点的三维坐标它的取值范围为[-1，1]的归一化坐标，因此本发明将其还原到对应的三维体数据的下标取值(x，y，z)。三维体数据的每一个点对应的下标取值为整数，而采样点的坐标为浮点数，因此需要做插值计算。计算用于光照模型的该点位置的梯度，也就是在x，y，z三个方向求一阶导数，也需要做类似的处理。为了获得平滑的插值和梯度计算结果，本发明采用b样条函数。

一维b样条函数表示为：

其中，ni，p表示b样条基函数，p表示b样条基函数的幂次，pi表示第i个控制点。常用的线性b样条为一次函数，有2个控制点，二次b样条函数有3个控制点，三次b样条函数有4个控制点。

对于3d体数据，将一维b样条函数扩展到三维b样条函数表示为：

其中，ni，p表示b样条基函数，p，q，r分别表示每个维度的b样条函数的幂次，pi，j，k表示三个维度的控制点，u，v，ω分别表示x，y，z轴的偏移量。本发明可以根据体数据的具体情况，选取一次线性样条，二次样条、三次样条曲线来拟合，不同的维度可以选取不同幂次的b样条。例如：如果x，y，z三个维度都选择三次样条，那么一共有4x4x4＝64个点参与计算，如果x维度选择一次样条，y维度选择二次样条，z维度选择三次样条，一共有2x3x4＝24个点参与计算。

三维b样条一阶、二阶、三阶函数控制点的矩阵分别表示为：

其中，参数u位于控制点i和i+1之间。三次b样条函数对参数u求导，得到梯度矩阵表示为

例如：计算得到三维体数据的下标坐标px，y，z，其中x，y，z的取值为正的浮点数；pi，j，k表示与之邻近的取值小于(x，y，z)的体素的最大整数坐标，其中u＝x-i，v＝y-j，ω＝z-k。以x方向为例，计算x方向的梯度包括以下步骤1至3：：

步骤1、如果x，y，z三个方向均采用三次b样条函数，一共需要64个点进行计算。pi-1，j-1，k-1，pi，j-1，k-1，pi+1，j-1，k-1，pi+2，j-1，k-1这四个点根据公式c3(u)和c′3(u)分别计算对应的点pi′，j-1，k-1位置处的插值和x方向梯度；

类似的，pi-1，j，k-1，pi，j，k-1，pi+1，j，k-1，pi+2，j，k-1四个点应用公式c3(u)和c′3(u)得到点pi′，j，k-1位置处的插值和x方向梯度；pi-1，j+1，k-1，pi，j+1，k-1，pi+1，j+1，k-1，pi+2，j+1，k-1四个点应用公式c3(u)和c′3(u)得到点pi′，j+1，k-1位置处的插值和x方向梯度；pi-1，j+2，k-1，pi，j+2，k-1，pi+1，j+2，k-1，pi+2，j+2，k-1四个点应用公式c3(u)和c′3(u)得到点pi′，j+2k-1位置处的插值和x方向梯度；以此类推，得到pi′，j-1，k，pi′，j，k，pi′，j+1，k，pi′，j+2，k；pi′，j-1，k+1，pi′，j，k+1，pi′，j+1，k+1，pi′，j+2，k+1；pi′，j-1，k+2，pi′，j，k+2，pi′j+1，k+2，pi′，j+2，k+2位置处的插值和x方向梯度，共计16个点；

步骤2、将pi′，j-1，k-1，pi′，j，k-1，pi′，j+1，k-1，pi′，j+2，k-1四个点应用公式c3(v)得到pi′，j′，k-1；将pi′，j-1，k，pi′，j，k，pi′，j+1，k，pi′，j+2，k四个点应用公式c3(v)得到pi′，j′，k；以此类推得到pi′，j′，k+1，pi′，j′，k+2共计4个点；

步骤3、将pi′，j′，k-1，pi′，j′，k，pi′，j′，k+1，pi′，j′，k+2四个点应用公式c3(ω)得到最终的采样点pi′，j′，k′＝px，y，z位置处的取值和x方向梯度值。

与上述步骤类似的，计算y方向梯度时，在步骤1中忽略计算c′3(u)，而在步骤2中计算c′3(v)；同理，计算z方向梯度时，在步骤1，2中忽略计算c′3(u)和c′3(v)，而在步骤3中计算c′3(ω)。

s7、采用光照模型计算采样点的颜色值和阻光度并进行累加求和，得到成像平面上对应像素点的取值进行三维可视化渲染。

在本实施例中，本发明根据步骤s6得到当前采样点位置处的取值，通过传递函数可以获得其对应的阻光度的值ai；再根据颜色编码可以获得其对应的rgb的值。

本发明根据坐标系中视点的坐标、光源的坐标、采样点的坐标，以及计算得到的采样点位置处梯度值，根据光照模型可以计算出当前采样点的颜色取值ci。光照模型可以选用phong模型、改进的blinn-phong模型、基于真实感的pbr模型等，将环境光、漫反射光以及镜面反射光进行融合，得到最终的颜色值。

对采样点的颜色值ci和阻光度ai进行累加求和，表示为

ci＝ci-1+(1-ai-1)ci

ai＝ai-1+(1-ai-1)ai

其中，ci和ai分别表示第i个采样点的颜色值和阻光度，ci-1和ai-1分别表示第i个采样点之前的颜色值和阻光度的累计值，它们的初始值均设置为0。

通过上述计算可以得到最终的成像平面上对应像素点的取值，进而进行三维可视化渲染。

本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。