针对并行体绘制实现剪切波数据三维可视化的访存方法

摘要

本发明提供一种针对并行体绘制实现剪切波数据三维可视化的访存方法，提出了4基访存优化策略，包括纹理存储器策略、全局存储器策略、全局+共享存储器策略、针对频繁调用剪切波数据的视角矩阵而设计的常量存储器策略、绘制质量的评价方法与准则。通过本发明，可以提高基于CUDA架构下采用体绘制方法实现的剪切波数据三维可视化的绘制效率。从而更好地帮助科研人员准确了解板块构造、动力学过程以及地震的发震机理等。进而有助于地质勘探、大陆板块、地震解释等领域的发展。

说明书

技术领域

本发明涉及可视化技术领域，尤其是涉及一种针对并行体绘制实现剪切波数据三维可视 化的访存方法。

背景技术

地震波是目前我们所知道的唯一能够穿透地球内部的波。它是在断层破裂时，使断层两 侧的岩石体发生相对移动，并释放出大量的能量。大部分的能量在克服摩擦力而做功，从而 转换为热能；而另一部分则以弹性波的形式以震源为中心向四周的地壳传播。其中，伴随着 断层两侧岩体快速移动产生的推力而形成的压缩波，被称为P波；沿着断层面相对位移而形 成的弹性波被称为剪切波（shear-wave），简称为S波。地震波由震源发出后会通过介质向各 个方向传播，在传播过程中还会发生反射、折射和散射。因各种不同的地下矿产资源在构造 上都会具有不同的特征属性，因此通过分析在地表上接收到的地震波到达的时间及形态，就 能对地下岩层的结构、深度、形态等做出精确推断，从而实现矿产资源的勘察。

由于剪切波只受控于剪切模量和密度结构，且剪切模量的大小与物质的温度的呈指数关 系，使得通过对剪切波的研究能帮助科研人员准确了解板块构造、动力学过程以及地震的发 震机理等，因此对剪切波的三维可视化能提供更多的地下结构和地球动力学信息，有助于地 质勘探、大陆板块、地震解释等领域的发展。

当前对于剪切波的研究主要集中在数据的处理和反演上，对于地震剪切波的解释和阐述 主要还是基于定性的结果。然而，目前主要的图像显示都是以二维静态为主，主要应用的软 件为GMT（Generic Mapping Tools,Wessel,2012），通过给定层面上的速度大小，在平面上 展示每层或者每个剖面上的结果，不能给出整体的三维构架状状态。这对三维的空间展示和 认识造成了一定的困难。如果能够定量的给出三维的可视化图像，特别是动态的三维剪切波 演化图像，将能在很大程度上增加对大陆形变、动力学演化以及地震孕震机理的认识和理解。

二十世纪八十年代以来，地震数据规模迅速膨胀，由过去的几万道数据，发展到今天的 百万道数据。同时研究中使用的地震属性也越来越多。随着可视化技术在地震勘探领域的广 泛应用，出现了大量相关可视化研究工作。

Gerald D.Kidd第一次对地震数据进行三维建模通过直接集成的构造方法，结合地层和深 度的数据，更加细致的展现了小区域的具体细节，在短时间内提供了更高质量的地震数据诠 释。Castanie,L等人针对石油和天然气勘探设计一种体数据漫游系统，通过基于沿体数据和 高效体积分的方法对大几千兆的体数据进行积分，并使用通用的多模态体积渲染框架呈现高 质量的绘制，最后使用交互式的漫游方法对绘制结果进行显示。Huw James等人则通过加入 光照模型对地震数据的显示质量进行了深入研究，在体绘制模型中引入光照模型能将体数据 以高光和阴影的方式区分开来，从而吸纳详细的内部结构和数据耦合的空间关系。

虽然三维可视化技术在地震、勘探等领域内已得到广泛应用，但由于传统体绘制算法的 复杂及计算量大，且需要研究人员有一定的可视化专业基础，以及需要掌握GPU的硬件特性 和渲染绘制编程语言，因此给研究人员利用三维可视化的方法对地震数据进行观测造成了一 定的难度，现在的地震领域科研人员任然是依靠多个层面及切面的方式对地震体数据进行分 析。早期的基于Cg语言的GPU架构虽然能给体绘制算法带来显著的效率提升，但是GPU 的编程模型还是针对纹理数据进行处理，因此在使用该架构进行地震数据可视化时，需要将 地震数据先通过特殊的数据结构，根据GPU的硬件架构的限制，开发出只适合该模型的可视 化工具。由于可移植性差，且对开发人员要求较高，因此对地震数据的三维可视化的难度上 并没有带来多大的改观。

CUDA（统一计算设备）是NVIDIA推出的GPU通用计算产品，它的编程更简单，功 能更强大，架构体系也更加适用于GPU通用计算，而基于CUDA实现的三维可视化算法 较好地解决了传统GPU加速存在的缺陷，达到了很高的实时交互速度以及不错的绘制质量。 但是难点是绘制时基于CUDA架构的访存效率的提升，进而更进一步提高程序的运行效率。

发明内容

本发明针对体数据及光线投射算法的特点，提供一种针对并行体绘制实现剪切波数据三 维可视化的访存方法。

本发明的技术方案为一种针对并行体绘制实现剪切波数据三维可视化的访存方法，基于 CUDA架构实现剪切波数据三维可视化，所述CUDA架构包括纹理存储器、全局存储器、共 享存储器、常量存储器，访存的数据对象包括体数据的三维纹理、体数据和视角矩阵，

当需要对体数据的三维纹理进行访存的实现方式时，采用纹理存储器策略，包括利用纹 理寄存器建立体数据的三维纹理，在对纹理坐标的拾取中通过硬件插值滤波过程提供额外的 浮点处理能力；

当需要对体数据进行访存的实现方式时，采取全局存储器策略，包括利用全局存储器建 立体数据；在遍历体数据的过程中，通过为以采样点位置为中心周围8个体素的三线性插值 得到采样点的值；并在满足合并访问条件时对全局存储器进行合并访问，当对访问全局存储 器发生非合并访问时，通过共享存储器充当缓冲区域，调整Half-warp中线程的访问顺序，将 对全局存储器的非合并访问优化为合并访问；

当需要对视角矩阵进行访存的实现方式时，采用常量存储器策略，包括从OpenGL内建 模型矩阵中获得视角矩阵，再将视角矩阵拷贝到常量存储器；在射线需要与视角矩阵进行计 算时直接调用视角矩阵。

而且，用一阶导数的方差值来度量边缘曲线的光滑度，作为剪切波数据三维可视化的图 片质量的评判准则。

本发明设计了纹理存储器策略、全局存储器策略、全局+共享存储策略，以及常量存储器 策略，共同实现并行体绘制中对剪切波数据三维可视化的访存优化。

附图说明

图1是本发明实施例的纹理建立流程图。

图2是本发明实施例的纹理拾取流程图。

图3是本发明实施例的体数据在全局存储器中的建立流程图。

图4是本发明实施例的体数据在局存储器中的采样流程图。

图5是本发明实施例的float3型数据合并访问示意图。

图6是本发明实施例的全局+共享策略的采样流程图。

具体实施方式

本发明技术方案可由本领域技术人员结合计算机软件技术实现自动运行。以下结合附图 和实施例详细说明本发明技术方案。

本发明提供一种针对并行体绘制实现剪切波数据三维可视化的访存方法，基于CUDA架 构实现剪切波数据三维可视化，所述CUDA架构包括纹理存储器、全局存储器、共享存储器、 常量存储器，访存的数据对象包括体数据的三维纹理、体数据和视角矩阵。针对剪切波数据 绘制过程中体数据的三维纹理的建立、体数据建立和遍历、访问时段不对齐和间隔访问问题 以及视角矩阵的存取实现访存优化。

本发明中，对体数据的三维纹理进行访存的实现方式，包括利用纹理寄存器建立体数据 的三维纹理，在对纹理坐标的拾取中通过硬件插值滤波过程提供额外的浮点处理能力。

当需要建立体数据的三维纹理而对文理存储器进行访存时，此时采用纹理存储器策略。 在体绘制算法需要大量精细计算获得较高的结果图像质量时，采用纹理存储器策略支持对体 数据进行频繁的读取。

实施例设计的纹理存储器策略如下：

实施例中利用纹理存储器对体数据的三维纹理进行存储的具体方案如图1所示，利用纹 理寄存器建立体数据的三维纹理，由体数据类型得到通道类型描述符，再由通道类型描述符 和体数据大小得到设备内存3D Array指针，最后由体数据主机内存指针、通道类型描述符和 设备内存3D Array指针得到体数据3D Texture纹理。

实施例中利用纹理存储器对体数据的三维纹理进行读取的具体方案如图2所示的纹理坐 标的拾取，是利用CUDA提供的tex3D()函数来根据体数据3D Texture纹理和采样点位置对浮 点型纹理坐标进行拾取，可记为Tex3D(tex,pos)，其中tex为体数据3D Texture纹理，pos为浮 点型纹理坐标；然后进行硬件插值滤波，最终的到采样点sample，其中的硬件插值滤波过程 即为纹理存储器的特性之一，它不会占用可编程单元，提供了额外的浮点处理能力。

本发明中，对体数据进行访存的实现方式，包括利用全局存储器建立体数据；在遍历体 数据的过程中，通过为以采样点位置为中心周围8个体素的三线性插值得到采样点的值；并 在访问全局存储器时，通过共享存储器充当缓冲区域，调整Half-warp中线程的访问顺序，将 对全局存储器的非合并访问优化为合并访问。

当需要建立体数据及对全局存储器中的三维体数据进行遍历而对全局存储器进行访存 时，采取全局存储器策略。实施例设计的全局存储器策略如下：

全局存储器占据了显示设备显存的绝大部分，拥有较高的存储空间，且整个grid中的任 意thread都能对全局存储器进行任意位置的读写，同时能提供很高带宽，但是没有缓存机制， 具有较高的访存延迟，因此在对全局存储器进行访存是采取全局存储器策略。

实施例中利用全局存储器对体数据进行存储的具体方案如图3所示，其中包括利用全局 存储器建立体数据的过程，体数据设备端数组指针d_PitchedPtr为CUDA架构中现有 cudaPitchedPtr结构体的一个实例，包含了体数据在全局存储器中的指针ptr和体数据对齐后 的pitch值。设体数据大小为volumeSzie，为体数据设备端数组指针d_PitchedPtr分配 volumeSzie大小空间，然后根据体数据主机内存指针进行HosttoDevice（主机到设备）方向 数据拷贝。

实施例中利用全局存储器对体数据进行读取的具体方案如图4所示，其中包括利用全局存 储器中遍历体数据的过程，体数据的三维包括depth维（深度）、height维（高度）和width维 （宽度）。设体素位置为pos(x,y,z)，根据体数据设备端数组指针d_PitchedPtr和体素位置 pos(x,y,z)依次得到体数据depth维、体数据height维和体数据width维，从而得到体素voxel。

采样点的值为以采样点位置为中心，周围8个体素的三线性插值得到的，因此需要以图 3所示的方法得到8个体素，再插值得到最终的采样点。

一个Half-warp是16个线程的集合，并且在CUDA架构中是以一个Half-warp为单位进 行指令发射的。通过能够满足一定条件，使得只需要一次传输就能同时处理这些线程的请求 的话，来提升全局存储器的访存效率。这个条件即为合并访问（coalesced access）条件。不 同计算能力的CUDA设备在合并访问条件上也有不同要求。

在计算能力为1.0和1.1的设备上，一个half-warp对段的访问是线性且连续的，即第i 个线程必须访问段里的第i个字，并且half-warp访问的段地址需要对齐到每个线程访问的字 长的16倍，而且只支持对字长为32bit、64bit和128bit的数据进行合并访问。如果满足合并 访问条件，但其中有线程不访存，那么看作一次合并访问；如果不满足合并访问条件，则该 half-warp被分解为16次的串行访问。

在计算能力1.2及以上的设备上，合并访问条件的限制被大大放宽，首先，对字长的访 问被扩展为8bit、16bit、32bit、64bit和128bit，其次，对段的访问不需要顺序访问和对齐访 问，只要在一个段内，都符合合并访问条件。即使跨越了两个段，也只会被分解为两次满足 合并访问条件的传输，而不像1.0和1.1的设备被分解为16次，大大提高了全局存储器的访 问效率。

为解决对体数据的采样时访问时段不对齐和间隔访问带来的访问效率下降和全局存储器 带宽的浪费，本发明提出一种共享存储器与全局存储器相结合的策略。当对全局存储器发生 非合并访问时，可通过共享存储器充当缓冲区域，调整Half-warp中线程的访问顺序，从而优 化为合并访问。实施例设计了全局+共享策略如下：

如图5所示，利用共享存储器及CUDA提供的线程同步机制，可以将对float3型数据的 一次性读取的非合并访问转化为合并访问。

以256threads/block为例，转化伪代码为：

//步骤1

_share_float smem_data[256*3];

smem_data[threadIdx.x]=gmem_data[threadId];

smem_data[threadIdx.x+256]=gmem_data[threadId+256];

smem_data[threadIdx.x+512]=gmem_data[threadId+512]

_syncthreads();

//步骤2

Calc();

//步骤3

gmem_data[threadId]=smem_data[threadIdx.x];

gmem_data[threadId+256]=smem_data[threadIdx.x+256];

gmem_data[threadId+512]=smem_data[threadIdx.x+512];

其中，256threads/block表示大小为256的线程块；_share_float smem_data[256*3]定义 256×3的存储共享存储器数据的数组，threadId表示进程的Id，smem_data[]表示全局存储器 中数据的数组。

步骤1为取数据到共享存储器。将float3型数据以3个float的形式从全局存储器取到共 享存储器，并使用CUDA架构中现有的syncthreads()同步线程。在这一步中实现了非合并访 问到合并访问的转化。如图所示，每个线程读取以自己线程ID为起始偏移位置的0，256， 512三个float，使得线程t0、t1、t2、t3…t255能连续的读取256个float型数据，从而实现 了符合合并访问条件的访问。

步骤2为数据计算。与原始的calc()中的计算步骤一致，只是数据换成了从共享存储器 smem_data中根据threadId读取。函数calc()表示对数据做的处理，可以是任何处理。

步骤3为数据写回。与步骤1类似，每个线程以自己线程ID为起始偏移位置的0，256， 512从共享存储器中向全局存储器写回三个float。

该策略在全局存储器中建立体数据的流程与全局存储器策略一样，见图6：根据体数据 设备端数组指针d_PitchedPtr和体素位置pos(x,y,z)依次得到体数据depth维、体数据height 维和体数据width维后，(float)Offset=0、(float)Offset=256、(float)Offset=512，Offset 表示偏移量，每个线程以自己线程ID为起始偏移位置的0，256，512从共享存储器中向全局 存储器写回三个float，即Float3voxel。

结合上述对通过共享存储器转化合并访问的分析，全局+共享策略的采样将单个float3型 体素的读取分解为3个float数据并发读取。采样过程如图6所示。

与全局存储器策略不同的是，配合共享存储器后，在从全局存储器读取体素时，将float3 型数据分解成3个float，再通过256threads/block来进行符合合并访问条件的并发读取，数据 从全局存储器读取到共享存储器后，再以float3型格式读取数据。

本发明中，对视角矩阵进行访存的实现方式，包括从OpenGL内建模型矩阵中获得视角 矩阵，再将视角矩阵拷贝到常量存储器；在射线需要与视角矩阵进行计算时直接调用视角矩 阵。当需要运用到视角矩阵进行计算而对常量寄存器时，采用常量存储器策略。

实施例设计的常量存储器策略针对的是用于控制图像视角变换的视角矩阵，实施例中利 用常量存储器对视角矩阵进行存储的方式为：

视角矩阵是3x4大小的float型矩阵，在光线投射算法中，每一条射线在计算采样点时都 需要根据当前的视角矩阵计算出眼睛所在的位置以及射线的归一化法向量。

CUDA架构中提供了64KB大小的常量存储器用于存储需要被频繁访问的只读数据。一 个half-warp的线程只需要一个周期即可访问常量寄存器并得到所需数据，效率大大高于从显 存中读取数据。

OpenGL（现有的图形程序接口）中内建的模型矩阵为4x3大小的float型数组，通过现 有OpenGL函数glMatrixMode(GL_MODELVIEW)将当前激活矩阵设定为模型矩阵，再通过 现有OpenGL函数glGetFloatv(GL_MODELVIEW_MATRIX,modelView)即可将模型矩阵拷贝 到事先声明的float[16]数组modelView当中，再对modelView求逆即可得到视角矩阵。

再通过CUDA架构提供的cudaMemcpyToSymbol()方法来将主机中的数据拷贝到常量存 储器中。

在射线需要与视角矩阵进行计算时从常量存储器直接调用即可。

为了比较绘制效率，实施例还设计了以图像边缘的一阶导数的方差模型，作为衡量不同 优化策略下的体绘制的效率。由于方差可以用来度量随机变量和其数学期望（即均值）之间 的偏离程度，而图像边缘的一阶导数为图像边缘曲线各处斜率的体现，如果图像质量比较高， 则边缘曲线会比较光滑，较少出现大的起伏和断点，即较少有斜率大小的特别大的改变，因 此这种情况下变圆曲线一阶导数的方差值也会比较低。因此，本发明用一阶导数的方差值来 度量边缘曲线的光滑度，以此来评价绘制质量。

实施例的评价准则模型设计如下：

①使用canny算子提取图像边缘。确定临界值k，读取每一列像素值，从左到右， 若像素值开始小于k，则认为进入图像区域，记录边界点的坐标值，得到边界坐标值数组。

②根据①中确定的边界坐标值数组，用d＝diff(s_i,s_j)，计算每两个像素s_i、s_j之间 的一阶导数值（边缘线的斜率），得到斜率数组。

$d=\mathrm{diff}(s_i,s_j)=\frac{\mathrm{\Delta y}}{\mathrm{\Delta x}}=(y_i-y_j)/(x_i-x_j)$

其中，（x_i，y_i）为像素点s_i的坐标值，（x_j，y_j）为像素点s_j的坐标值。统计的图 像为从整体图像中取得很小一部分，大小在100*100像素左右。将100*100的区域内的 图像边缘的一阶导进行计算。

③根据②中得到的斜率数组（即存储边缘曲线上像素点的一阶导数的数组，设数组 的长度为n），用计算斜率的方差值。其中，d_i为边缘曲线上任一像 素点处的一阶导数，即边缘曲线上像素点的个数，为斜率数组的平均值，

利用上述的评价方法来纵向比较针对体数据的前三种策略在相同条件下的绘制效率。比 较读取128x128x128的体数据时每秒渲染的帧数（fps），得出，全局存储器策略绘制效率太 低，导致算法性能不理想，纹理存储器策略和全局存储器策略绘制效率不错，在对体数据读 取不频繁的情况下，纹理存储器策略绘制效率比全局+共享策略略高。而随着对体数据采样间 距的缩小，即对存储器的读取次数逐渐增加的情况下，纹理策略的绘制效率下降非常快，而 全局+共享策略对实验条件的变化没有产生较大变化。从此实验可知，当算法需要对体数据频 繁读取时，全局+共享策略的要优于纹理策略。

同样，比较使用和未使用常量存储器策略时读取大小为128x128x128的体数据的fps，实 验结果表明，使用常量存储器优化视角矩阵的存储后，绘制效率几乎提升了一倍，因此由时 间结果可知，常量存储器策略在提升光线投射算法的性能上是非常有效的。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技 术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不 会偏离本发明的精神或者超越所附权利要求书所定义的范围。