基于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纹理坐标顺序地(比如沿顺时针方向或逆时针方向)组成每张多边形切片的顶点数组。其中，计算交点的具体方法是设体数据在局部坐标系下数据区间为, ，卜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的夹角，"为反射指数，其本质实质上是反映本文档来自技高网...

【技术保护点】
基于ＧＰＵ加速的三维医学图像显示方法，包括以下步骤：　步骤１：读入医学ＤＩＣＯＭ图像序列文件并以体数据的方式保存到系统内存；　步骤２：利用ＯｐｅｎＧＬ或者ＤｉｒｅｃｔＸ的三维图形库编程扩展接口函数ＡＰＩ，将步骤１保存在系统内存中 的体数据加载入ＧＰＵ显存，成为ＧＰＵ可以访问的３Ｄ纹理；　步骤３：计算生成代理几何体，具体分为以下几个步骤：　步骤３－１：计算体数据包络盒的顶点在视点坐标系下的坐标；　首先将体数据包络盒八个顶点的局部坐标转换成世界坐标，然 后通过视点世界坐标参数矩阵，将体数据包络盒八个顶点的世界坐标变换到视点坐标系中；　步骤３－２：计算代理几何体所包含的多边形切片数目；　在视点坐标系下，首先计算体数据包络盒的八个顶点在视线方向的最大值ｚ↓［ｍａｘ］和最小值ｚ↓［ｍ ｉｎ］；然后确定采样率Δｈ，即相邻两层多边形切片间的距离；最后计算代理几何体所包含的多边形切片数目Ｎ↓［Ｓｈｃｅｓ］为：　Ｎ↓［Ｓｈｃｅｓ］＝（ｚ↓［ｍａｘ］－ｚ↓［ｍｉｎ］／Δｈ　步骤３－３：计算代理几何体每张多边形切片的顶点 的世界坐标；　对于代理几何体中的每张多边形切片，通过计算每张多边形切片所在平面与体数据包络盒在世界坐标系下的交点，然后将这些交点的世界坐标和３Ｄ纹理坐标顺序地组成每张多边形切片的顶点数组；　步骤４：对步骤３所得的每张多边形切片， 逐像素进行光照计算和颜色计算，具体分为以下几个步骤：　步骤４－１：计算代理几何中每张多边形切片的每一像素点的体数据梯度；　一个三维数据ｆ（ｘ，ｙ，ｚ）的梯度由一个偏微分梯度算子表示为▽ｆ（ｘ，ｙ，ｚ）＝（＊ｆ／＊ｘ，＊ｆ／＊ｙ， ＊ｆ／＊ｚ）；具体采用基于６邻域采样差分的梯度计算方法，计算位于（ｘ，ｙ，ｚ）处的梯度▽ｆ（ｘ↓［ｉ］，ｙ↓［ｊ］，ｚ↓［ｋ］）：▽ｆ（ｘ↓［ｉ］，ｙ↓［ｊ］，ｚ↓［ｋ］）＝１／２（ｆ（ｘ↓［ｉ＋１］，ｙ↓［ｊ］，ｚ↓［ｋ］）－ｆ（ｘ↓［ｉ－１］，ｙ↓［ｊ］，ｚ↓［ｋ］），ｆ（ｘ↓［ｉ］，ｙ↓［ｊ＋１］，ｚ↓［ｋ］）－ｆ（ｘ↓［ｉ］，ｙ↓［ｊ－１］，ｚ↓［ｋ］），ｆ（ｘ↓［ｉ］，ｙ↓［ｊ］，ｚ↓［ｋ＋１］）－ｆ（ｘ↓［ｉ］，ｙ↓［ｊ］，ｚ↓［ｋ－１］））　步骤４－２：采 用颜色传递函数计算代理几何中每张多边形切片的每一像素点的颜色；　其中，所述颜色传递函数是一个类似ｃ＝ｆ（ｖ）的一维函数，ｃ为输...

【技术特征摘要】

【专利技术属性】
技术研发人员：解梅，张帆，
申请(专利权)人：电子科技大学，
类型：发明
国别省市：90[中国|成都]