一种快速高质量最大密度投影实现方法

摘要

本发明提供一种快速高质量最大密度投影实现方法，其结合现有改进方案，根据光线进入数据块的24种情况，利用八叉树结构对数据块进行递归处理，并利用八叉树最大最小值结构对透明数据块或子节点进行有效地跳过，同时采用并行计算提高CPU利用效率，有效减少运算时间，此外利用单元边界采样，减少数据访问量和计算量。从而大大减少数据的计算，实现实时高质量的最大密度投影。

说明书

技术领域

本发明提供一种投影实现方法，尤其涉及一种快速高质量最大密度投影实现方法。

背景技术

MIP(maximum intensity projection，最大密度投影)是常用的三维可视化技术之一，其优点是简单可行，能有效地从NMRI(Nuclear Magnetic Resonance Imaging，核磁共振成像)和CT(computed tomography，电子计算机X射线断层扫描技术)数据中获取血管、骨骼和软组织等结构的可视化图像，在血管和肿瘤组织等疾病的诊断工作中具有十分重要的辅助作用。

MIP的基本原理是：从投影平面上每一点发出的平行投影光线，光线从某个角度穿过体数据，每条光线上的最大值即是投影面上对应该光线的像素点的值。从原理得知，该算法要求每一条光线上所有的像素值都被求出进行比较，每个像素的值需要由体数据插值得到，因而该算法的计算量比较大；加之投影角度变化时，对于不同图像层而言，其数据寻址顺序与其存储顺序不一样，实际的计算量更大，可见速度和图像质量是MIP算法应解决的核心问题。

尽管现在计算机的计算速度、内存及其他加速设备的发展日新月异，但是对于庞大的体数据而言，MIP速度始终难以尽如人意。为解决这种状况，学者们基于现有硬件环境，根据不同的应用需求，在生成的图像质量和速度之间寻求可能的最佳方案，以尽可能地满足用户需求，同时保证一定的图像质量。由此，形成了许多不同的MIP算法，常用的有：光线投射法(Ray Casting)、错切-变形法(Shear-Warp)、抛雪球法(Splatting)和基于硬件的3D纹理映射方法(Haedware-Assisted 3Dtexture-mapping)。错切-变形法是速度最快的软件MIP算法，但是该算法存在一定的不足，如当视角为45°时会出现阶梯状走样，图像放大后变模糊等问题。抛雪球法是“以物体空间为序”的算法，其优点是能按照提数据存储顺序来存取对象，同时只有与图像相关的体素才被投射和现实，可大大减少体数据的访问，但是由于象平面上的投影面随着视点改变而随意缩放和旋转，因而计算对周围像素影像的范围和对其每一点所影像的大小是十分费时的，同时设计既能保证图像质量又高效的重构函数往往比较困难。

光线投射法是以图像空间为序的体绘制方法，它是从图像空间的每一个像素出发，按视线方向发射一条射线，射线穿过体数据进行采样，采样点中灰度最大值即为图像空间的像素值。该算法要求遍历每个体素，而且当视角改变时，数据采样点前后关系随之变化，需要重新采样，计算庞大，为此学者们提出了不少优化方法。例如，与抛雪球算法结合，采用并行计算；采用空间数据结构(金字塔结构、八叉树结构和k-d树结构等)研究数据内部的相关性来提高计算效率；将体数据分块存储，符合计算机的存储结构，减少数据访问时间；采用非均匀采样减少数据访问和计算量。

传统的最大密度投影存在以下问题：速度和图像质量不能并取，速度快的方法需要牺牲图像质量来换取，而图像质量高的速度慢，无法满足实际需要。

发明内容

本发明提供一种快速高质量最大密度投影实现方法，其结合现有改进方案，根据光线进入数据块的24种情况，利用八叉树结构对数据块进行递归处理，并利用八叉树最大最小值结构对透明数据块或子节点进行有效地跳过，同时采用并行计算提高CPU利用效率，有效减少运算时间，此外利用单元边界采样，减少数据访问量和预算量。从而大大减少数据的计算，实现实时高质量的最大密度投影。

本发明为解决上述技术问题所采用的技术方案为：

一种快速高质量最大密度投影实现方法，包括以下步骤：

A.数据预处理，将数据分块整理存储并建立最大值最小值八叉树；

B.光线初始化，计算出光线的进入点并根据进入点位置将所有光线分配给其进入的数据块；

C.根据光线方向将数据块分组，采用并行计算分别处理每组数据块；

D.处理组内每个数据块的所有光线，当光线穿出数据块时，判断其是否已经穿出数据体；

若没有穿出数据体，根据采样位置将光线分配给下一个数据块后进入E步骤；

E.继续下一个数据块处理；

F.所有分组处理完成。

所述B步骤中的光线初始化，先利用移动模板的方法得到每条光线的进入深度后，再通过叠加的方法计算得到每条光线的进入点。

所述A步骤中的的数据块大小为2ⁿ*2ⁿ*2ⁿ，n大于1，数据块小于CPU的缓存容量。

所述B步骤的所述移动模板大小为(2ⁿ+m)*(2ⁿ+m)*(2ⁿ+m)，其中m为大等于1的自然数，所述移动模板中心与每个待投影的数据块的中心重叠，所述移动模板记录有模板中光线的数目和模板数据，所述模板数据包括光线相对于模板中心的2D坐标偏移值和光线进入深度。

所述光线进入深度和进入点的计算方法为：

三维坐标表示为vector(x，y，z)，vector.x、vector.y、vector.z分别表示其X、Y、Z分量，

图像表示大小size(cx，cy)，size.cx、size.cy分别表示图像的宽和高，

体数据坐标系约定为Voxel坐标系，投影平面空间坐标系约定为View坐标系；体数据到投影平面空间的三维坐标变换为VoxeltoView，投影平面空间到体数据的三维坐标变换为ViewtoVoxel；

Voxel坐标vecVoxel与对应的View坐标vecView之间的转换关系为：

vecVoxel＝ViewtoVoxel×vecView

vecView＝VoxeltoView×vecVoxel

投影大小为imagesize，根据VoxeltoView计算数据块间在X、Y、Z三个方向的3D View坐标偏移DeltaBX、DeltaBY，DeltaBZ：

$\left(\begin{array}{c}\mathrm{DeltaBX}=\mathrm{VoxeltoView}\times \mathrm{Voxel}\left(\mathrm{32,0,0}\right)-\mathrm{VoxeltoView}\times \mathrm{Voxel}\left(\mathrm{0,0,0}\right)\\ \mathrm{DeltaBY}=\mathrm{VoxeltoView}\times \mathrm{Voxel}\left(\mathrm{0,32,0}\right)-\mathrm{VoxeltoView}\times \mathrm{Voxel}\left(\mathrm{0,0,0}\right)\\ \mathrm{DeltaBZ}=\mathrm{VoxeltoView}\times \mathrm{Voxel}\left(\mathrm{0,0,32}\right)-\mathrm{VoxeltoView}\times \mathrm{Voxel}\left(\mathrm{0,0,0}\right)\end{array}\right)$

当前数据块中心在投影空间中的坐标为CurPos，初始值为(0，0，0)，对于顺序为(bx，by，bz)的数据块而言，其对应的位置为：

Pos＝CurPos+DeltaBX×bx+DeltaBY×by+DeltaBZ×bz

对于该数据块内的第i个光线，其进入深度为：

Depth[i]＝Pos.z+EntryDepth[i]

其在投影面上对应的2D坐标索引值：

Index[i]＝Pos.x+Pos.y×imagesize.cx+PixelOffset[i]

以数据块为单位进行循环，计算出数据块内所有光线在的进入深度Depth[i]，根据其在投影面上的索引值Index[i]得到投影面上对应的进入深度Depth[index[i]]，如果Depth[i]＜Depth[index[i]]，则将更新Depth[index[i]]的值为Depth[i]；

得到投影面上所有光线的进入深度后，根据ViewtoVoxel求出光线在X、Y、Z三个方向3D Voxel坐标偏移DeltaX、DeltaY、DeltaZ：

$\left(\begin{array}{c}\mathrm{DeltaX}=\mathrm{ViewtoVoxel}\times \mathrm{View}\left(\mathrm{1,0,0}\right)-\mathrm{ViewtoVoxel}\times \mathrm{View}\left(\mathrm{0,0,0}\right)\\ \mathrm{DeltaY}=\mathrm{ViewtoVoxel}\times \mathrm{View}\left(\mathrm{0,1,0}\right)-\mathrm{ViewtoVoxel}\times \mathrm{View}\left(\mathrm{0,0,0}\right)\\ \mathrm{DeltaZ}=\mathrm{ViewtoVoxel}\times \mathrm{View}\left(\mathrm{0,0,1}\right)-\mathrm{ViewtoVoxel}\times \mathrm{View}\left(\mathrm{0,0,0}\right)\end{array}\right)$

同时求出Voxel(0，0，0)对应的View坐标Pos：

Pos＝ViewToVoxel×View(0，0，0)

则对于投影平面上点(x，y)，令其进入深度为depth，其进入点为：

voxel＝Pos+x×DeltaX+y×DeltaY+depth×DeltaZ

所述n为5，所述m为4，所述数据块大小为32*32*32，所述移动投影模板大小为36*36*36。

所述A步骤中的八叉树为对每个数据块建立的4级最大值最小值八叉树。

所述C步骤中处理每个数据块时透明体素跳过的方法为：

令第一个数据块及其子节点的起点为(x0，y0，z0)，棱长为2^(5-lev)的正方体，其中lev是其在八叉树中的深度，其内一点P(x，y，z)所在的子节点索引值为：

index＝((x-x0)＞＞(5-lev))+(((y-y0)＞＞(5-lev))＜＜1)

+(((z-z0)＞＞(5-lev))＜＜2)

其中lev＞＝1，结合面的索引值计算出该节点上所有入射光线所进入的子节点，计算好24种入射情况作为一个查找表，记为EntryList；

计算好光线在当前节点的子节点间的24种入射情况所依次穿过的子节点作为一个查找表，记为TraverseNodeList；

处理每个数据块时，数据块的最大值小于光线的当前值则跳过该数据块，否则计算光线进入的子节点索引及面索引，通过查找EntryList得到进入情况属于哪一种，根据进入情况查找TraverseNodeList得到光线依次穿过的1级子节点，依次处理完每个1级子节点即完成该数据块的处理；

处理1级子节点，若属于当前1级子节点的2级子节点的最大值小于当前光线的值则跳过，否则根据进入的2级子节点索引及面索引，，通过查找EntryList得到进入情况属于哪一种，根据进入情况查找TraverseNodeList得到光线依次穿过的2级子节点，采用递归的方式实现子节点的处理。

所述C步骤中的插值方法为单元边界采样方式的双线性插值。

该发明运行速度快，一般能够满足实时性要求。

本发明在初始化光线部分，在利用移动模板的方法得到每条光线的进入深度后如何快速有效地计算出每条光线的进入点。基本思想都是通过叠加的方法计算得到，而不是通过坐标变换或者其他方式，叠加的方法充分利用了点与点之间的变化关系，有效地减少了运算了，如果采用坐标的变换的关系的话，对于三维坐标变换，需要用一个4*4的浮点变换矩阵乘以一个三维浮点向量，这个运算量相当大，严重影响运算效率。本发明的voxel＝Pos+x×DeltaX+y×DeltaY+depth×DeltaZ，在循环处理中，前两项三维浮点向量乘法可以由三维浮点向量的一次加法代替，那么上式的计算就主要在最后一个三维浮点向量乘法了。

本发明在透明体素跳过中，通过建立的最大最小值八叉树，能够快速有效地跳过透明体素，有效减少计算量。通过建立由投影方向确定的24种入射情况，快速有效地将光线的当前采样点定位到Block的八叉树子节点中，并利用该子节点的最大值与光线的当前灰度比较，来确定是否跳过该节点。透明体素跳过的思想是体绘制计算光线进入点的有效方法，对于体绘制而言，当一个传递函数确定以后，在传递函数中阻光度为零的灰度认为是透明的，可以直接不参与计算。最大密度投影的透明体素跳过与体绘制有着异曲同工之妙，根据其成像原理，可以看出只有最大灰度才会对绘制结果产生影响，因此若光线的当前灰度已经大于将要处理的节点的最大灰度，那么则认为该节点对该光线没有作用，可以认为是透明的。关于最大密度中透明体素跳过的部分暂未见。

附图说明

图1为本发明实施例投影模板示意图；

图2为本发明实施例Block分组示意图；

图3为本发明实施例光线进入八叉树子节点的24种情况示意图；

图4为本发明实施例八叉树子节点排列顺序示意图；

图5为本发明实施例八叉树节点的面排列顺序

图6为本发明实施例等间距采样示意图；

图7为本发明实施例单元边界采样示意图。

具体实施方式

下面根据附图和实施例对本发明作进一步详细说明：

为了便于叙述和理解，给出下面的约定，在以后各个章节的叙述和说明中均遵照本小节的约定，

坐标约定：三维坐标表示为vector(x，y，z)，vector.x、vector.y、vector.z分别表示其X、Y、Z分量。

坐标约定：图像表示大小size(cx，cy)，size.cx、size.cy分别表示图像的宽和高。

坐标系及三维坐标变换约定：体数据坐标系约定为Voxel坐标系，投影平面空间坐标系约定为View坐标系；体数据到投影平面空间的三维坐标变换为VoxeltoView，投影平面空间到体数据的三维坐标变换为ViewtoVoxel。Voxel坐标vecVoxel与对应的View坐标vecView之间的转换关系为：

vecVoxel＝ViewtoVoxel×vecView                   (2-1)

vecView＝VoxeltoView×vecVoxel             (2-2)

数据块大小约定：数据块大小约定32*32*32。

数据块大小约定：数据块八叉树分级约定数据块本身为0级，每分一级级数加1。

数据预处理：为减少计算量和数据访问时间，需对体数据做如下预处理。

数据分块存储：当投影角度变化时，数据寻址顺序与其存储顺序不一样，导致实际的计算量更大。为降低数据访问时间，提高运算效率，学者们根据计算机硬件存储结构的特点，提出了体数据分块(Bricking)^[1]的方法，将体数据分成2ⁿ*2ⁿ*2ⁿ大小的Block(数据块)(例如，32*32*32)，使得每一个Block的大小在CPU的Cache容量范围内，这样利用Cache就可以完成每个数据块的处理，无需在内存与Cache之间反复传输数据。

Bricking除了可以降低数据访问时间、提高运算效率外，在服务器处理模式下，更有其独特优点：对于归档完毕的病人数据，需要做MIP的只需将数据整理一遍，保存在服务器上，就可以反复调用；并且对于VR(volume rendering，体绘制)而言，数据预处理方式都有前述优点。

最大值最小值八叉树建立：由MIP算法的基本原理可知，需要将光线上所有像素值计算出来进行比较求出最大值。而一般的医学提数据中与绘制结果无关的体素(即“透明”体素)占到总体素的70％-95，如何有效地跳过“透明体素”的处理，是影响运算效率的关键。当数据以Block方式存储时，处理每一个Block时，若已知Block的最大值小于光线的当前最大值，则该Block可以直接跳过，而无需处理。此外，空间数据结构的运用利于研究数据体内部的相关性，有效地跳过“透明”体素的处理，其中八叉树存储是一种较好的优化数据表达策略。在本文方法中，通过对每个Block建立4级最大值最小值8叉树，以利于在处理每个block时，能够对Block的每个子节点分别进行区别判断，减少计算，对于最大值最小值八叉树的运用在后面相应章节给出。

与Bricking优点类似，最大值最小值八叉树在服务器处理模式下，也具有独特的优点：在归档整理数据时，可以将该数据体以Block为单位，计算每个Block的最大值最小值八叉树并保存，之后即可反复调用；并且“透明”体素跳过的思想在VR中也适用，最大值最小值八叉树也可以有效减少VR的计算量。

光线初始化：为使用并行计算提高运算速度，需要计算出投影平面上每一条光线进入体数据的进入点，同时将光线进行排序分配给每个Block。

模板建立：为快速计算出投影平面上每条光线的进入点，采用模板平移^[2]的方法。文献2为计算出VR中光线的进入点，对每个Block四级子节点分别建立模板。对于MIP而言，仅需知道Block级的进入点即可，无须深入到其子节点。

计算投影模板采用的Block大小为36*36*36，采用36而不用32为了避免移动投影模板，离散化到投影面坐标出现漏点的问题，并假设投影模板中心与每个待投影的blcok的中心是重叠的。如图1所示，模板是包含数据块投影的最小矩形，它需要记录模板中光线的数目和模板数据，模板数据主要分为两部分PixelOffset和EntryDepth，其中PixelOffset为相对于模板中心的2D坐标偏移值，EntryDepth指光线进入深度，即进入点到投影面的距离。

创建光线，根据投影模板，通过循环Block，可以计算出整个图像每个像素进入深度。具体做法如下：

假设投影大小imagesize，根据VoxeltoView计算Block间在X、Y、Z三个方向的3D View坐标偏移DeltaBX、DeltaBY，DeltaBZ：

$\left(\begin{array}{c}\mathrm{DeltaBX}=\mathrm{VoxeltoView}\times \mathrm{Voxel}\left(\mathrm{32,0,0}\right)-\mathrm{VoxeltoView}\times \mathrm{Voxel}\left(\mathrm{0,0,0}\right)\\ \mathrm{DeltaBY}=\mathrm{VoxeltoView}\times \mathrm{Voxel}\left(\mathrm{0,32,0}\right)-\mathrm{VoxeltoView}\times \mathrm{Voxel}\left(\mathrm{0,0,0}\right)\\ \mathrm{DeltaBZ}=\mathrm{VoxeltoView}\times \mathrm{Voxel}\left(\mathrm{0,0,32}\right)-\mathrm{VoxeltoView}\times \mathrm{Voxel}\left(\mathrm{0,0,0}\right)\end{array}\right)---(2-3)$

当前Block中心在投影空间中的坐标为CurPos，初始值为(0，0，0)。那么对于顺序为(bx，by，bz)的Block而言，其对应的位置则为：

Pos＝CurPos+DeltaBX×bx+DeltaBY×by+DeltaBZ×bz            (2-4)

对于该Blcok内的第i个光线而言，其进入深度为：

Depth[i]＝Pos.z+EntryDepth[i]                              (2-5)

其在投影面上对应的2D坐标索引值：

Index[i]＝Pos.x+Pos.y×imagesize.cx+PixelOffset[i]         (2-6)

以Block为单位进行循环，计算出Block为内所有光线在的进入深度Depth[i]，根据其在投影面上的索引值Index[i]得到投影面上对应的进入深度Depth[index[i]]，如果Depth[i]＜Depth[index[i]]，则将更新Depth[index[i]]的值为Depth[i]。

得到投影面上所有光线的进入深度后，根据ViewtoVoxel求出光线在X、Y、Z三个方向3D Voxel坐标偏移DeltaX、DeltaY、DeltaZ：

$\left(\begin{array}{c}\mathrm{DeltaX}=\mathrm{ViewtoVoxel}\times \mathrm{View}\left(\mathrm{1,0,0}\right)-\mathrm{ViewtoVoxel}\times \mathrm{View}\left(\mathrm{0,0,0}\right)\\ \mathrm{DeltaY}=\mathrm{ViewtoVoxel}\times \mathrm{View}\left(\mathrm{0,1,0}\right)-\mathrm{ViewtoVoxel}\times \mathrm{View}\left(\mathrm{0,0,0}\right)\\ \mathrm{DeltaZ}=\mathrm{ViewtoVoxel}\times \mathrm{View}\left(\mathrm{0,0,1}\right)-\mathrm{ViewtoVoxel}\times \mathrm{View}\left(\mathrm{0,0,0}\right)\end{array}\right)---(2-7)$

同时求出Voxel(0，0，0)对应的View坐标Pos：

Pos＝ViewToVoxel×View(0，0，0)                 (2-8)

则对于投影平面上点(x，y)，假设其进入深度为depth，那么其进入点为：

voxel＝Pos+x×DeltaX+y×DeltaY+depth×DeltaZ    (2-9)

根据进入点坐标可以计算光线进入Block的索引值，达到将光线分配给Block的目的。

重建优化，并行计算：基于Block的光线投影法的最大的特点是：重建的Block之间相互独立，利于并行计算。采用并行计算可以有效提高计算机的CPU使用效率，减少运算时间。

当光线分配给每个Block之后，可以根据投影方向将Block分组由前到后排列建立BlockList，组内的Block之间是相互独立的，如图2所示。每个Block内光线的当前采样点在该Block内，初始值有光线初始化部分得到。通过处理由前到后排列的各组Block完成整个数据体的处理。在组内，Block间是相互独立的，可以按任意顺序处理，每个Block都可以用一个并行计算线程来独立完成。当一条光线离开当前处理Block组将进入下一个Block组，同时也被分配给下一个进入的Block，直至所有光线都穿出数据体完成MIP绘制。

透明体素跳过：在处理每个Block时，利用其4级最大值最小值八叉树，可有效跳过透明体素减少计算量。

最大值最小值八叉树的运用的关键是快速有效地将当前采样点定位到Block的子节点中。对于Block的节点而言，有六个入射面，每个面上有四个子节点面，那么就有24种入射情况，如图3所示。

八叉树的子节点的排列顺序如图4所示，对于第一个Block及其子节点而言，可以看成是一个起点为(x0，y0，z0)，棱长为2^(5-lev)的正方体，其中lev是其在八叉树中的深度，那么其内一点P(x，y，z)所在的子节点索引值为：

index＝((x-x0)＞＞(5-lev))+(((y-y0)＞＞(5-lev))＜＜1)  (2-10)

+(((z-z0)＞＞(5-lev))＜＜2)

其中lev＞＝1，结合面的索引值(如图5所示)可以计算出该节点上所有入射光线所进入的子节点，24种入射情况也就确定，而且不随入射角度变化，可以作为一个查找表而预先计算好，记为EntryList。

光线在当前节点的子节点间最多穿行三次，若投影角度确定，24种入射情况所依次穿过的子节点也是确定的，因此也可事先将其计算出作为一个查找表，记为TraverseNodeList。

处理每个Block时，若Block的最大值小于光线的当前值则跳过该Block，否则计算光线进入的子节点索引及面索引，通过查找EntryList得到进入情况属于哪一种，根据进入情况查找TraverseNodeList得到光线依次穿过的1级子节点，依次处理完每个1级子节点即完成该Block的处理。处理1级子节点与处理Block类似，若属于当前1级子节点的2级子节点的最大值小于当前光线的值则跳过，否则根据进入的2级子节点索引及面索引，，通过查找EntryList得到进入情况属于哪一种，根据进入情况查找TraverseNodeList得到光线依次穿过的2级子节点，以此类推。由于子节点的处理和父节点的处理方法类似，可以采用递归的方式实现。

插值是MIP重建中最关键的一步，直接影响重建图像的质量。传统的插值方法是采用等间距采样(如图6(a)所示)，利用三线性插值求出采样点的灰度值。本文采用单元边界采样方式，与等间距采样相比，该采样方式可以有效介绍数据访问量和计算。采样点在立方体单元格的边界时，根据三线性插值的权重关系可知，有一半的权重是0，实际相当于一个双线性插值，因此可以直接用双线性插值完成。双线性插值与三线性插值相比，其数据访问量和计算量都只有三线性插值的一半，有效减少地数据访问量和计算量。

本领域技术人员不脱离本发明的实质和精神，可以有多种变形方案实现本发明，以上所述仅为本发明较佳可行的实施例而已，并非因此局限本发明的权利范围，凡运用本发明说明书及附图内容所作的等效结构变化，均包含于本发明的权利范围之内。