专利名称：基于gpu加速的三维医学图像显示方法
技术领域：
本发明属于医学图像处理技术领域，土要利用GPU (显示处理单元)的强大的并行流处 理能力，来对体绘制进行加速，相较亍传统CPU的体绘制技术，该方法使得体绘制能交互式 实时显示。
背景技术：
医学图像在医生诊断中的辅助作用越来越明显，但是仅可以从二维截面方向对人体进行 观察。为了提高医疗诊断和治疗规划的准确性和科学性，需要由二维断层图像序列转变为具 有直观立体效果的三维图像。然而目前的三维医疗辅助t会断系统都需要工作站级别的运行平
台，才能基本满足实时性的要求，由于价格问题使得这种系统难以推广。目能国际上计算机
图形图像学的科学家们提出将大规模的数据处理、运算等任务放在GPU运行作为甜沿的研 究方向，具有CPU不可比拟的速度优势，在消费级的GPU硬件上实现大规模数据的交互式 可视化。
当甜通用的体绘制技术， 一般是基于传统的CPU进行计算，对于传统消费级别大众PC，
很难实时进行交互式显示。如果要交互式实时显示运算， 一般要在并行图形工作站上进行运 行，成本相较大大提高。
发明内容
本发明针对现有基于CPU的医学图像显示方法存在运算速度慢、无法交互式执行的弊 端，提供一种基于GPU加速的三维医学图像显示方法，利用GPU强大的流运算能力进行加 速，从而实现交互式体绘制的效果。
本发明技术方案如下
基于GPU加速的三维医学图像显示方法，如图1所示，包括以下,:p驟 歩骤1:读入医学DICOM图像序列文件并以体数据的方式保存到系统内存。 体数据是由有顺序的2维医学DICOM图像序列构成。将这些医学图像的图像相关信息 读入系统内存(图像分辨率，层间距，图像像素信息)。为了更好的视觉效果，体数据应该首 先进行预处理，比如用图像滤波器进行除噪，并且通过插值层间距数据來得到更加细致的效 果。
歩骤2:利用OpenGL或者DirectX的三维图形库编程扩展接口函数API,将歩骤1保存在系统内存中的体数据加载入GPU显存，成为GPU可以访问的3D纹理。具体接口函数API 可以采用OpenGL库的扩展函数glTexImage3DEXT来将体数据加载为3D纹理。
歩骤3:计算生成代理几何体，具体分为以下几个歩骤
歩骤3-1:计算体数据包络盒的顶点在视点坐标系下的坐标。
这个操作涉及到一个顶点坐标变换操作，其实质就是将体数据包络盒的顶点坐标由体数 据的局部坐标系转换到视点坐标系。首先将体数据包络盒八个顶点的局部坐标转换成世界坐 标，然后通过视点世界坐标参数矩阵，将体数据包络盒八个顶点的世界坐标变换到视点坐标 系中。
歩骤3-2:计算出代理儿何体所包含的多边形切片数目。
在视点坐标系下，首先计算体数掂包络盒的八个顶点在视线方向的最大值zma,和最小值
—7m。；然后确定采样率A力，即相邻两层多边形切片W的距离；最后计算代理几何体所包含的 多边形切片数目A^"、为
步骤3-3:计算代理几何体每张多边形切片的顶点的世界坐标。
对于代理儿何体中的每张多边形切片，通过计算每张多边形切片所在平面与体数据包络
盒在世界坐标系下的交点，然后将这些交点的出:界坐标和3D纹理坐标顺序地(比如沿顺时
针方向或逆时针方向)组成每张多边形切片的顶点数组。其中，计算交点的具体方法是
设体数据在局部坐标系下数据区间为[-x,,,m, xm, ], [imin, }mll，]，卜z,，, 2mm ]，视点位置为
'视线方向向量为^ : fe,，>W，^),采样率为厶A 。对于每张多边形切
片，即垂直于视线方向、间隔为A/7的平行平面与立方体边界线段切割的交集区域。 一张多边
形切片所在平面可以表示成方程
d/厂p0) = 0
根据视线方向和所经过的体数据内点A,的位置，顺序计算体数据边界线段与平面
= 0的交点。由于我们将体数据局部坐标系和世界坐标系重合，极大方便我们的计 算，因为边界线段都是与坐标轴平行的。
^骤4:对歩骤3所得的每张多边形切片，逐像素进行光照计算和颜色计算，具体分为
以下几个歩骤步骤4-1:计算代理几何中每张多边形切片的每一像素点的体数据梯度。
—个三维数据/Oc,少,Z)的梯度由一个偏微分梯度算子表示为V/(X,y,Z)—i，！,！)；
&
具体采用基于6邻域采样差分的梯度计算方法，计算位于(X,)'，Z)处的梯度V/(X,，^，^):
V/(乂,,7,,zJ =1 /2(,(x,+1'少,，z》—,(x, i,-"'4)，、V''z*)-/化')',(x''X'zw)—/化'少"z",》
歩骤4-2:采用颜色传递凼数计算代理儿何中每张多边形切片的每一像素点的颜色。
由于体数据图像为灰度图像，在视觉上效果以及色度分辨率不利于对细节进行分辨，所
以在本方法中通过颜色传递函数进行体数据体素灰度到RGB颜色空间的映射。
颜色传递函数是一个类似c二/0)的一维函数，c为输出的颜色值，v为体数据的灰度值。
颜色传递函数的作用是用来将灰度值的体数据映射到RGB颜色空间的体数据。可以通过一个颜色映射表來进行实现，在GPU里面，可以通过l维纹理来存放这个颜色映射表。在医学图像序列毕.面，灰度数据是从O到N离散存放的。在这早，N是一个山颜色存储位数i来决定
的自然数# = 2'。通过预计算出函数映射表的1维纹理，通过对这个纹理进行采样，这样很容易得到输出纹理。这段代码在GPU eg语舎単面只是一个简单的指令texld(float x)。
步骤4-3:根据Phong光照模型计算代理几何体中每张多边形切片的每一像素点的光栅化像素着色值。
反射光分为三个部分，环境光，漫反射光，镜面反射光。
环境光本质是指光源经过环境等作用，然后间接对物体的影响，其本质是在物体和环境之间多次反射，最终达到平衡时的一种光强。本质上说，这种光强在物体上各点是不相同的，但是在简巾.光照模型中近似地认为同一环境下的环境光，其光强分布是均匀的，而且环境光在任何一个方向上的光强强度都相同。表现出來是一个常数。
漫反射光漫反射光的定义如下，假设当光源來自一个方向时，漫反射光将光均匀向各
方向反射并传播丌米，与反射光的强度与视点无关，因为漫反射光本质是由反射表面的微表
面凹凸不平所引起的，因而漫反射光的空间分布是均匀的。运算的式子如下；假设入射光的强度为/,,，物体表面微元上点P的法向为W，从点P指向入射光的光源的单位向量为Z，
两者间的夹角为6 ，由光强能量的余弦定律可以计算得到漫反射光强为乙=/,，>< & xcos(S), 0e(O,"/2)在上式子中，《,是与物体有关的光的漫反射系数，系数满
足0<《，<1 。当丄、V为单位向量时，也可用如下向量运算的形式进行表达
6乙=;x & x &' AO在有多个光源的情况下'可以有如下的表示/rf = &Z ; x (丄,.W)
镜面反射光最典型的例子就是对于理想镜面，镜面的反射光集中在一个方向，并遵守反射定律。将其推广，对一般的光滑表面，反射光集中在一个空间区间范围内，且由反射定律决定的反射方向光强最大。因此，对于同一点来说，从不同位置所观察到的镜面反射光强是不同的，这个与视线与镜面反射光方向有关。镜面反射光强的式子可表示为
=x《、x cos"W), P e (O,;r / 2)其中X、是与物体表面光学性质有关的镜面反射系数，P为
视线方向F与反射方向A的夹角，"为反射指数，其本质实质上是反映了物体表面的光泽程度， 一般为1 2000，通过式子可以看出，数目越大物体表面越光滑。镜面反射光将会在反射方向附近形成很亮的光斑，称为高光现象。同样，如果将F和i 都格式化为单位向量，那
么镜面反射光强可表示为/、 /p.《(i .r)"其中，A 口J由向量反射计算式子i^2W(7V.丄)-丄计算。与甜面漫反射光相向的，当对十多个光源的情形，镜面反射光强可以表示为
/、=《^[/,(7 ,J/)"]
歩骤5:通过Alpha混合将代理几何体中所有多边形切片合成二维医学图像。
在关闭深度测试的情况下，打丌Alpha混合，将将代理几何体中所有多边形切片合成三维医学图像。为了提高渲染速度，可以将代理几何体保存到显示列表，这样就能用单--指令渲染多个顶点，从而提高绘制效率。
Alpha混合是一个将代理几何体若十多边形层的颜色进行混合的方法。图形GPU硬件用的光照模型是针对多边形顶点的，而不能被用在体绘制上面，所以必须关闭光照。并且，由于打丌了 Alpha混合，就必须渲染代理几何体単面的每一个多边形层，所以必须关闭深度测试。下面进行简单解释。
Alpha混合OpenCH,中的绝大多数特效都与某些类型的(一般是色彩)混合有关。混色的定义为，将某个象素的颜色和已绘制在屏幕上与其对应的象素颜色相互结合。至于如何结合这两个颜色则依赖于颜色的alpha通道的分量值，以及/或者所使用的混色函数。Alpha通常是位于颜色值末尾的第4个颜色组成分量。大多数情况都认为Alpha分量代表材料的透明度。这就是说，alpha值为0.0时所代表的材料是完全透明的。alpha值为1.0时所代表的材料则是完全不透明的。其实混合的基本原理是就将要分色的图像各象素的颜色以及背景颜色均按照RGB规则各自分离之后，根据 一 图像的RGB颜色分量+alpha值+背景的RGB颜色分量^l-alpha值) 一 这样一个简单公式来混合之后，最后将混合得到的RGB分量重新合并，公式如下Co/or/ma, = Co/or/ra , xa + Co/o;ti x (1 -a) 。 OpenGL按照卜.面的公式计算这两个象素
的混色结果。这个公式一般会生成混合的透明/半透明的效果。
深度测试像素片断测试其实就是测试每一个像素，只有通过测试的像素才会被绘制，没有通过测试的像素则不进行绘制。OpenGL提供了多种测试操作，利用这些操作可以实现一些特殊的效果。"深度测试"的概念在绘制三维场景的时候特别有用。在不使用深度测试的时候，如果我们先绘制一个距离较近的物体，再绘制距离较远的物体，则距离远的物体因为后绘制，会把距离近的物体覆盖掉，这样的效果并不是我们所希望的。如果使用了深度测试，则情况就会有所不同每当一个像素被绘制，OpenGL就记录这个像素的"深度"(深度可以理解为该像素距离观察者的距离。深度值越大，表示距离越远)，如果有新的像素即将覆盖原來的像素时，深度测试会检查新的深度是否会比原来的深度值小。如果是，则覆盖像素，绘制成功；如果不是，则不会覆盖原来的像素，绘制被取消。这样一来，即使我们先绘制比较近的物体，再绘制比较远的物体，则远的物体也不会覆盖近的物体了。实际上，只要存在深度缓冲区，无论是否启用深度测试，OpenGL在像素被绘制时都会尝试将深度数据写入到缓冲区内，除非调用了 glDepthMask(GL—FALSE)來禁止写入。这些深度数据除了用于常规的测试外，还可以有一些有趣的用途，比如绘制阴影等等。除了深度测试，OpenGL还提供了剪战测试、Alpha测试和模板测试。
本发明提供了--种基于GPU加速的医学图像显示方法，相比于现有的基于CPU的医学图像显示方法，本发明具有很高的运算速度，可在普通消费级别大众PC上实现实时交互式显示；而无须使用图形工作站，使得成本大大降低。
图1为木发明流程示意图。
具体实施例方式
本发明采用上述技术方案，分别针对MRI和CT医学DICOM图像序列进行显示验证，均具有较好的显示效果。需要说明的是，由于MRI图像山于噪声较大，需要进行除噪预处理;而CT图像较为清晰，无需进行除噪预处理，可直接以体数据的方式进行显示。
权利要求
1、基于GPU加速的三维医学图像显示方法，包括以下步骤步骤1读入医学DICOM图像序列文件并以体数据的方式保存到系统内存；步骤2利用OpenGL或者DirectX的三维图形库编程扩展接口函数API，将步骤1保存在系统内存中的体数据加载入GPU显存，成为GPU可以访问的3D纹理；步骤3计算生成代理几何体，具体分为以下几个步骤步骤3-1计算体数据包络盒的顶点在视点坐标系下的坐标；首先将体数据包络盒八个顶点的局部坐标转换成世界坐标，然后通过视点世界坐标参数矩阵，将体数据包络盒八个顶点的世界坐标变换到视点坐标系中；步骤3-2计算代理几何体所包含的多边形切片数目；在视点坐标系下，首先计算体数据包络盒的八个顶点在视线方向的最大值zmax和最小值zmin；然后确定采样率Δh，即相邻两层多边形切片间的距离；最后计算代理几何体所包含的多边形切片数目NShces为<maths id="math0001" num="0001" ><math><![CDATA[ <mrow><msub> <mi>N</mi> <mi>Shces</mi></msub><mo>=</mo><mfrac> <mrow><mo>(</mo><msub> <mi>z</mi> <mi>max</mi></msub><mo>-</mo><msub> <mi>z</mi> <mi>min</mi></msub><mo>)</mo> </mrow> <mi>Δh</mi></mfrac> </mrow>]]></math></maths>步骤3-3计算代理几何体每张多边形切片的顶点的世界坐标；对于代理几何体中的每张多边形切片，通过计算每张多边形切片所在平面与体数据包络盒在世界坐标系下的交点，然后将这些交点的世界坐标和3D纹理坐标顺序地组成每张多边形切片的顶点数组；步骤4对步骤3所得的每张多边形切片，逐像素进行光照计算和颜色计算，具体分为以下几个步骤步骤4-1计算代理几何中每张多边形切片的每一像素点的体数据梯度；一个三维数据f(x，y，z)的梯度由一个偏微分梯度算子表示为<maths id="math0002" num="0002" ><math><![CDATA[ <mrow><mo>▿</mo><mi>f</mi><mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>,</mo> <mi>z</mi> <mo>)</mo></mrow><mo>=</mo><mrow> <mo>(</mo> <mfrac><mrow> <mo>∂</mo> <mi>f</mi></mrow><mrow> <mo>∂</mo> <mi>x</mi></mrow> </mfrac> <mo>,</mo> <mfrac><mrow> <mo>∂</mo> <mi>f</mi></mrow><mrow> <mo>∂</mo> <mi>y</mi></mrow> </mfrac> <mo>,</mo> <mfrac><mrow> <mo>∂</mo> <mi>f</mi></mrow><mrow> <mo>∂</mo> <mi>z</mi></mrow> </mfrac> <mo>)</mo></mrow><mo>;</mo> </mrow>]]></math> id="icf0002" file="A2009100598640002C2.tif" wi="45" he="9" top= "204" left = "139" img-content="drawing" img-format="tif" orientation="portrait" inline="yes"/></maths>具体采用基于6邻域采样差分的梯度计算方法，计算位于(x，y，z)处的梯度 id="icf0003" file="A2009100598640002C3.tif" wi="23" he="4" top= "220" left = "157" img-content="drawing" img-format="tif" orientation="portrait" inline="yes"/><maths id="math0003" num="0003" ><math><![CDATA[ <mrow><mo>▿</mo><mi>f</mi><mrow> <mo>(</mo> <msub><mi>x</mi><mi>i</mi> </msub> <mo>,</mo> <msub><mi>y</mi><mi>j</mi> </msub> <mo>,</mo> <msub><mi>z</mi><mi>k</mi> </msub> <mo>)</mo></mrow><mo>=</mo><mn>1</mn><mo>/</mo><mn>2</mn><mrow> <mo>(</mo> <mi>f</mi> <mrow><mo>(</mo><msub> <mi>x</mi> <mrow><mi>i</mi><mo>+</mo><mn>1</mn> </mrow></msub><mo>,</mo><msub> <mi>y</mi> <mi>j</mi></msub><mo>,</mo><msub> <mi>z</mi> <mi>k</mi></msub><mo>)</mo> </mrow> <mo>-</mo> <mi>f</mi> <mrow><mo>(</mo><msub> <mi>x</mi> <mrow><mi>i</mi><mo>-</mo><mn>1</mn> </mrow></msub><mo>,</mo><msub> <mi>y</mi> <mi>j</mi></msub><mo>,</mo><msub> <mi>z</mi> <mi>k</mi></msub><mo>)</mo> </mrow> <mo>,</mo> <mi>f</mi> <mrow><mo>(</mo><msub> <mi>x</mi> <mi>i</mi></msub><mo>,</mo><msub> <mi>y</mi> <mrow><mi>j</mi><mo>+</mo><mn>1</mn> </mrow></msub><mo>,</mo><msub> <mi>z</mi> <mi>k</mi></msub><mo>)</mo> </mrow> <mo>-</mo> <mi>f</mi> <mrow><mo>(</mo><msub> <mi>x</mi> <mi>i</mi></msub><mo>,</mo><msub> <mi>y</mi> <mrow><mi>j</mi><mo>-</mo><mn>1</mn> </mrow></msub><mo>,</mo><msub> <mi>z</mi> <mi>k</mi></msub><mo>)</mo> </mrow> <mo>,</mo> <mi>f</mi> <mrow><mo>(</mo><msub> <mi>x</mi> <mi>i</mi></msub><mo>,</mo><msub> <mi>y</mi> <mi>j</mi></msub><mo>,</mo><msub> <mi>z</mi> <mrow><mi>k</mi><mo>+</mo><mn>1</mn> </mrow></msub><mo>)</mo> </mrow> <mo>-</mo> <mi>f</mi> <mrow><mo>(</mo><msub> <mi>x</mi> <mi>i</mi></msub><mo>,</mo><msub> <mi>y</mi> <mi>j</mi></msub><mo>,</mo><msub> <mi>z</mi> <mrow><mi>k</mi><mo>-</mo><mn>1</mn> </mrow></msub><mo>)</mo> </mrow> <mo>)</mo></mrow> </mrow>]]></math></maths>步骤4-2采用颜色传递函数计算代理几何中每张多边形切片的每一像素点的颜色；其中，所述颜色传递函数是一个类似c＝f(v)的一维函数，c为输出的颜色值，v为体数据的灰度值；颜色传递函数的作用是用来将灰度值的体数据映射到RGB颜色空间的体数据；步骤4-3根据Phong光照模型计算代理几何体中每张多边形切片的每一像素点的光栅化像素着色值；所述Phong光照模型中反射光分为三个部分，环境光，漫反射光，镜面反射光；所述环境光采用同一环境下的环境光，其光强分布是均匀的，而且环境光在任何一个方向上的光强强度都相同；所述漫反射光的光强为Id＝Ip×Kd×cos(θ)，θ∈(0，π/2)；其中Ip为入射光的强度；Kd是与物体有关的光的漫反射系数，且满足0＜Kd＜1；θ为物体表面微元上点P的法向N与点P指向入射光的光源的单位向量L的夹角；所述镜面反射光的光强为Id＝Ip×Ks×cosn(θ′)，θ′∈(0，π/2)；其中Ip为入射光的强度；Ks是与物体表面光学性质有关的镜面反射系数，θ′为视线方向V与反射方向R的夹角，n为反射指数；步骤5通过Alpha混合将代理几何体中所有多边形切片合成三维医学图像；在关闭深度测试的情况下，打开Alpha混合，将将代理几何体中所有多边形切片合成三维医学图像。为了提高渲染速度，可以将代理几何体保存到显示列表，这样就能用单一指令渲染多个顶点，从而提高绘制效率。
2、 根据权利要求1所述的基于GPU加速的二维医学图像显示方法，其特征在于，歩骤 1读入医学D1C0M图像序列文件并以体数据的方式保存到系统内存时，为了更好的视觉效 果，先对体数据应改进行除U呆或插值层间数据的预处理。
3、 根据权利要求1所述的基于GPU加速的三维医学图像显示方法，其特征在于，歩骤 2中具体接口函数API采用OpenGL库的扩展函数glTexImage3DEXT來将体数据加载为3D 纹理。
4、 根据权利要求1所述的基于GPU加速的三维医学图像显示方法，其特征在于，歩骤 3-3中交点的具体计算方法是设体数据在局部坐标系下数据区间为[-xmil', x■ ]， [-ym,n. ym,n ]. [-—7,細,zmm ],视点位置为^^:",,<,,，_y_，^. .)，视线方向向量为^:(.、,，U,〃)，采样率为A/7; —张多边形切片所 在平面可以表示成方程 7r .(x — p0) = 0根据视线方向和所经过的体数据内点/;。的位置，顺序计算体数据边界线段与平面 (A'r (x — / 0) = 0白勺交点。
全文摘要
基于GPU加速的三维医学图像显示方法，属于医学图像处理技术领域。首先将医学DICOM图像序列文件以体数据的方式保存到系统内存；然后利用OpenGL或者DirectX的三维图形库编程扩展接口函数API，将体数据加载入GPU显存；再计算生成代理几何体，并代理几何体中的多边形切片逐像素进行光照计算和颜色计算；最后通过Alpha混合将代理几何体中所有多边形切片合成三维医学图像。本发明相比于现有的基于CPU的医学图像显示方法，本发明具有很高的运算速度，可在普通消费级别大众PC上实现实时交互式显示；而无须使用图形工作站，使得成本大大降低。
文档编号G06T15/00GK101593345SQ20091005986
公开日2009年12月2日 申请日期2009年7月1日 优先权日2009年7月1日
发明者帆 张, 梅 解 申请人:电子科技大学