CN101336831A

CN101336831A - 实时三维医学超声图像的重建方法

Info

Publication number: CN101336831A
Application number: CNA2008100302330A
Authority: CN
Inventors: 郭境峰; 王海潮; 蔡泽杭
Original assignee: Shantou Institute of Ultrasonic Instruments Co Ltd
Current assignee: Shantou Institute of Ultrasonic Instruments Co Ltd
Priority date: 2008-08-13
Filing date: 2008-08-13
Publication date: 2009-01-07
Anticipated expiration: 2028-08-13
Also published as: CN101336831B

Abstract

一种实时三维医学超声图像的重建方法，依次包括下述步骤：(1)对来自超声三维传感器的体积图像数据进行预处理，获得重建图像所需的体积图像数据；(2)对步骤(1)获得的体积图像数据的数据点进行分类，将反映同一物质的数据点归入同一类；(3)对数据点赋值，包括对数据点赋颜色值和不透明度值；(4)在图形处理器中对体积图像数据进行三维坐标变换，将立方体数据变换成双曲面数据；(5)在图形处理器中进行重采样，将数据点的原始坐标转换成数据点在双曲面数据中的坐标；(6)图像合成。本发明既不丢失原始图像信息，又能实现图像的实时显示，实现三维医学超声图像的实时重建，并使最后形成的图像能够真实反映被测物体的实际形态。

Description

实时三维医学超声图像的重建方法

技术领域

本发明涉及图像数据处理，具体地说，涉及一种实时三维医学超声图像的重建方法。

背景技术

传统的医学二维图像的可视化和解释已被广为接受，但体积图像技术发明后，由于每个体积图像数据(即三维图像数据)都包含一系列的二维切片(即二维图像数据)，必须观看的二维切片的数量呈爆炸式的增加，因而迫切需要更直观地观察体积图像数据。

三维可视化的主要目的是将医学影像设备得到的一系列的二维切片(这一系列的二维切片构成体积图像数据)，使用计算机图形学的技术，构造出一个器官的三维模型，并且非常逼真的显示给观看者，使观看者能够直观地观察体积图像数据。随着医疗影像设备的发展，医学影像的数据量越来越大，这给传统的三维重建与绘制技术带来了非常大的挑战。

目前三维可视化可分为二大类：面绘制(Surface Rendering)和体绘制(Volumn Rendering)。面绘制包括经典的Marching Cubes算法，Cuberille算法等；体绘制包括经典的Ray Casting算法、Splatting算法、Shear Warp算法等。

面绘制方法首先由体积图像数据构造出中间几何图元(如曲面、平面等)，然后再由传统计算机图形学实现画面绘制，产生显示在屏幕上的二维图像。这时，可以将中间几何图元的生成过程看作是映射，但是这时的映射只是将原始数据中的部分属性映射成平面或曲面，在绘制平面或者曲面的过程中丢失了部分原始信息，因而这种方法构造出的可视化图形虽然可以产生比较清晰的等值面图像，但是不能反映整个原始数据场的全貌及细节。

体绘制方法与面绘制方法完全不同，它不构造中间几何图元，而是直接由体积图像数据产生显示在屏幕上的二维图像，是近年来得到迅速发展的一种体积图像数据可视化方法，它直接对精炼及处理后的数据进行绘制，以生成二维图像。体绘制方法不会丢失原始数据，相对于面绘制方法具有图像质量高的优点，但其计算量很大，计算时间较长，难以达到实时显示的效果。

另外，在医学超声图像领域，典型的超声波成像数据不能从笛卡尔网格获得，而必须由声学网格来获得，如图1所示。超声三维传感器一般由凸形传感器芯部与电动机构成，制作时需在凸形传感器芯部外围安装一个外壳进行保护，通过该外壳直接接触被测物体(如人体)，同时外壳内壁与凸形传感器芯部必须存在一定空隙，保证凸形传感器芯部可以在外壳内自由转动。扫描时，电动机带动凸形传感器芯部绕固定轴进行摆动，通过控制摆动的总角度与总时间、以及每次摆动的角度与时间等各种参数来控制凸形传感器芯部的运动范围与运动时间；凸形传感器芯部每摆动一个角度，对被测物体(如人体)发出一个超声波脉冲，并接收回波，从而获得一个二维切片；凸形传感器芯部在摆动的过程中，每次发射超声波脉冲的角度并不一致，因而得到的各个二维切片之间是不平行的；由于在凸形传感器芯部扫描方向与摆动方向形成的面都是扇形面，因此形成的体积图像数据是双曲面数据。但是，图像处理装置从超声三维传感器得到的体积图像数据实际上是立方体数据，即从超声三维传感器得到的体积图像数据由一系列相互平行的二维切片所组成，这些二维切片在内存中按时间次序进行摆放，根据该立方体数据形成的图像与被测物体不相符，会导致形状失真与变形，不能真实反映被测物体的实际形态。

发明内容

本发明所要解决的技术问题是提供一种实时三维医学超声图像的重建方法，采用这种方法重建图像，既不丢失原始图像信息，又能实现图像的实时显示，并且形成的图像能够真实反映被测物体的实际形态。采用的技术方案如下：

一种实时三维医学超声图像的重建方法，依次包括下述步骤：

(1)对来自超声三维传感器的体积图像数据进行预处理，获得重建图像所需的体积图像数据；

在获得来自超声三维传感器的体积图像数据后，首先对该体积图像数据进行预处理，包括原始体积图像数据的格式转换、剔除冗余数据，并导出所需的体积图像数据。

超声三维传感器得到的数据一般为16位或者16位以上，但屏幕只能显示8位灰阶图像，这就需要我们把原始体积图像数据进行格式转换，使其符合屏幕的显示要求。

冗余数据一般出现在超声三维传感器扫描的开始位置与结束位置，可采取一些方法(如去除法和清零法)进行剔除，这样既能保证图像质量，又可减少了数据量。

(2)对步骤(1)获得的体积图像数据的数据点进行分类，将反映同一物质的数据点归入同一类；

体积图像数据中包含了来自被测物体(如人体或动物体)中多种物质(如人体或动物体中的皮肤、血液、骨骼、心脏等组织或器官)的数据，由于这些物质的密度各不相同，对超声波的反射也不相同，因此超声三维传感器从不同物质探测到的数据也不相同。按被探测的物质种类对体积图像数据的数据点进行分类，其目的是找出数据点与各种物质之间的对应关系，使反映同一物质的数据点归入同一类，而反映不同物质的数据点则不在同一类。

对体积图像数据的数据点进行分类可采用下述方法：根据体积图像数据所在应用领域的背景知识，或者对全部数据点进行统计，设定若干阈值d_m(m＝1，2...，n，即m为正整数)；如果各数据点以f(x_i，y_j，z_k)表示，则将满足条件d_m≤f(x_i，y_j，z_k)＜d_m+1的数据点归入同一类。

(3)对数据点赋值，包括对数据点赋颜色值和不透明度值；

由于数据点本身具备灰度信息而不具备颜色信息，并且需要显示体积图像数据内部的信息，因此需要给数据点赋颜色值与不透明度值。

对数据点赋值可通过转换函数来进行，从而给数据点赋上颜色值与不透明度值。从数据点的灰度值转换到颜色值，其实质是对一个预先设定好的颜色表(颜色表可根据对当前的被探测的物质进行统计分析计算后得出)进行查找，例如一个最简单的线性查找表Color[i]＝i(i＝0，1...，255)，则体积图像数据中的最小值对应查找表中的Color[0]＝0，最大值对应查找表中的Color[255]＝255，其余的值则按线性插值进行查找；赋不透明度值的方法也是对一个预先设定好的不透明度表进行查找。

通常，对同一类数据点赋上相同的颜色，即同一物质在图像中采用一种颜色表示，颜色值则对应颜色的饱和度。由于是人为进行赋值，得到的最终图像的颜色将是伪彩色。

不透明度表示数据点对光的吸收度，值为1时表示完全不透明，值为0时表示完全透明。

(4)在图形处理器中对体积图像数据进行三维坐标变换；

由于典型的超声三维传感器都是凸形传感器，但得到的体积图像数据是立方体数据，因此需要对体积图像数据进行三维坐标变换，将立方体数据变换成双曲面数据，经过三维坐标变换，对每个二维切片都进行了顶点坐标变换，并重新组织了二维切片之间的位置关系。三维坐标变换可分两步进行，首先立方体数据经过第一次变换先形成锥形数据，然后再经过第二次变换后形成双曲面数据，这样形成的数据才会使最终在屏幕上形成的图像与被测物体相符合，不会导致形状失真与变形，能够真实反映被测物体的实际形态。

在一具体方案中，进行三维坐标变换，具体可利用图形处理器(GPU)中的可编程处理器，依照下述步骤对体积图像数据进行纹理坐标变换(纹理坐标变换是三维坐标变换中的一种)：

①对可编程顶点处理器中的顶点信息不行进处理，直接传递到可编程片段处理器中，并保存原始的纹理坐标；

②选取XY轴形成的平面，依据超声三维传感器的参数信息(如探头半径、探头角度、扫描线数、显示角度、显示点数等)，对Y方向数据进行缩放，在XY平面上形成梯形；利用步骤①保存的纹理坐标来对XZ方向进行坐标变换，首先对Z方向数据进行缩放，在XZ平面上形成梯形；经上述坐标变换后形成锥形数据；

③求出对当前Y对应X变换后形成的角度，对X方向进行缩放，并进行平移，产生弧形，就完成了XY方向上的坐标变换；求出对当前Z对应X变换后形成的角度，对X方向进行缩放，并进行平移，产生弧形，就完成了XZ方向上的坐标变换；经上述坐标变换后形成双曲面数据；

④根据超声三维传感器的参数信息重新计算体积图像数据中心。由于经过纹理坐标变换后，体积图像数据中心的位置发生变化，因此应当重新计算体积图像数据中心，以确保体积图像数据能绕体积图像数据中心进行旋转、平移、缩放等动作。

在另一具体方案中，进行三维坐标变换，具体可利用图形处理器(GPU)中的可编程处理器，依照下述步骤对体积图像数据进行纹理坐标变换(纹理坐标变换是三维坐标变换中的一种)：

②选取XY轴形成的平面，依据超声三维传感器的参数信息(如探头半径、探头角度、扫描线数、显示角度、显示点数等)，对Y方向数据进行缩放，在XY平面上形成梯形；

③求出对当前Y对应X变换后形成的角度，对X方向进行缩放，并进行平移，产生弧形，就完成了XY方向上的坐标变换；

④利用步骤①保存的纹理坐标来对XZ方向进行坐标变换，首先对Z方向数据进行缩放，在XZ平面上形成梯形；

⑤求出对当前Z对应X变换后形成的角度，对X方向进行缩放，并进行平移，产生弧形，就完成了XZ方向上的坐标变换；

⑥根据超声三维传感器的参数信息重新计算体积图像数据中心。由于经过纹理坐标变换后，体积图像数据中心的位置发生变化，因此应当重新计算体积图像数据中心，以确保体积图像数据能绕体积图像数据中心进行旋转、平移、缩放等动作。

(5)在图形处理器中进行重采样；

经过三维坐标变换后，数据点的坐标通常会发生变化，也就是说，数据点在双曲面数据中的坐标与其在立方体数据中的原始坐标不同，根据数据点的原始坐标不能在双曲面数据中找到相应数据点的数据值(包括颜色值和不透明度值)，因此必须进行重采样，将数据点的原始坐标转换成数据点在双曲面数据中的坐标。优选重采样是在图形处理器(GPU)的可编程片段处理器中进行的，可编程片段处理器根据数据点的原始坐标计算得到数据点在双曲面数据中的坐标，并一一对应地用数据点在双曲面数据中的坐标替换原始坐标。每个数据点对片段处理器来说相当于一个运算核，片段处理器对各运算核进行并行计算，从而提高运算速度以达到实时的目的。

(6)图像合成。

图像合成可在图形处理器中进行。图像合成一般采用光线投射法的原理进行，其目的是沿某一像素点发出一条射线，计算该射线上各数据点的颜色值及不透明度值，并进行合成，形成该像素点的最终颜色与最终透明度。合成的原理是：体积图像数据可以看成是多个二维切片组成的，这些二维切片经过坐标变换后投影到显示屏幕上，则多个切片上的点会映射到显示屏幕上的同一个点，显示屏幕上的这个点的最终显示颜色与最终透明度由多个二维切片上的点来混合决定。

经过图像合成后，可在显示屏幕上显示图像，该图像能够反映被测物体的实际形态。

上述实时三维医学超声图像的重建方法中，其关键在于在图形处理器中对体积图像数据进行三维坐标变换和重采样，而对体积图像数据进行预处理、对体积图像数据的数据点进行分类、对数据点赋值和图像合成等步骤可采用常规技术手段实现。

为了使屏幕上显示的图像能够较逼真地模拟出自然界的光照效果，本发明还通过图形处理器(GPU)的可编程处理器，构建光照模型(可用CG(C for graphics)显卡编程语言表达光照模型)，在重采样前对各数据点进行光照量的计算，即在步骤(3)对数据点赋值后，对各数据点进行光照量的计算，并将光照量与数据点的颜色值进行混合，从而模拟出较逼真的光照效果。

用CG语言表达的光照模型如下：

surfaceColor＝emi ssive+ambient+diffuse+specular；

其中surfaceColor表示光照量，光照量surfaceColor是放射光分量emissive、环境光分量ambient、漫反射光分量diffuse和镜面反射光分量specular等光照分量之和；各分量具体描述如下：

1)放射光分量emissive表示了由物体表面所发出的光，它的作用是独立于所有光源的，是一种在计算完其它所有光照分量后添加的颜色，它与其它分量没有联系，用K_e表示，即emissive＝K_e；

2)环境光分量ambient代表了入射光在一个场景里经过多次折射后、看起来像是来自四面八方一样的光，并不依赖于光源的位置，只依赖于物体材质的反射能力与照射在物体上的环境光的颜色，表示为ambient＝K_a*globalAmbient，其中K_a为环境光分量的系数，globalAmbient为入射环境光的颜色；

3)漫反射光分量diffuse代表从一个表面相等地向所有方向反射出去的方向光，通常物体表面微观尺寸上是非常粗糙的，当入射光线到达物体表面时，光线会向各个方向反射，表示为diffuse＝K_d*lightColor*max(N·L，0)，其中K_d是物体的漫反射颜色，lightColor是入射漫反射光的颜色，N是规范化的表面法向量(通常情况下我们并不能直接获得每个顶点的表面法向量，这时需要用纹理点的表面法向量来代替)，L是规范化的指向光源的向量；

4)镜面反射光分量specular代表了从一个表面主要的反射方向附近被反射的光，它非常依赖于观察者的位置，不仅受光源和物体的反射颜色性质的影响，也受物体表面的光泽度的影响，表示为

specular＝K_s*lightColor*facing*(max(N·H)，0)^shininess，

其中K_s是物体材质的镜面反射颜色，lightColor是入射镜面反射光的颜色，N是规范化的表面法向量，H是视点向量与光源向量的点积，shininess表示物体的光泽度。

在GPU中的可编程处理器中对各数据点的各个光照分量进行计算后，对各个光照分量进行累加，得到各数据点的光照量。

如果需要使最后在屏幕上显示的图像具有聚光灯和/或半透明效果，还需计算距离衰减，距离衰减量attenuationFactor表示为attenuationFactor＝1/(K_c+K_l*d+k_q*d*d)，其中d是物体到光源的距离，K_c、K_l、k_q为控制衰减量的常量；然后将距离衰减量与上述各个光照分量累加，得到光照量。

进行重采样之后，将各数据点的光照量与颜色值进行加权求和(光照量与颜色值各自的权重可根据需要设定，当要求图像较亮时光照量的权重较大，当要求图像较暗时光照量的权重较小)，可得出各数据点的光亮度值，从而使最后在屏幕上显示的图像具有较理想的光照效果。在这种情况下，步骤(6)的图像合成一般采用光线投射法的原理进行，其目的是沿某一像素点发出一条射线，计算该射线上各数据点的光亮度值及不透明度值，并进行合成，形成该像素点的最终光亮度与最终透明度。

为了使屏幕上显示的图像具有凹凸感，增强图像的真实感，优选在上述对各数据点光照量的计算中还具有明暗计算，利用明暗计算突出显示不同物质之间的边界面。明暗计算的步骤如下：

(a)计算各数据点的梯度值(梯度是指灰阶的变化程度)；

梯度值的计算可在中央处理器中进行，由如下中心差分公式求出，其中f(x_i，y_i，z_i)表示源数据点的灰度值，G_f(x_i，y_i，z_i)表示目标数据点的梯度值：

G_{f} (x_{i}, y_{i}, z_{i}) = [\begin{matrix} \frac{1}{2} (f (x_{i - 1}, y_{i}, z_{i}) - f (x_{i + 1}, y_{i}, z_{i})) \\ \frac{1}{2} (f (x_{i}, y_{i - 1}, z_{i}) - f (x_{i}, y_{i + 1}, z_{i})) \\ \frac{1}{2} (f (x_{i}, y_{i}, z_{i - 1}) - f (x_{i}, y_{i}, z_{i + 1})) \end{matrix}]

(b)将步骤(a)得到的梯度值作为纹理保存进图形处理器(GPU)中；

(c)图形处理器的可编程顶点处理器接收到纹理后，将纹理传递给可编程片段处理器；可编程片段处理器接收纹理后，用纹理点的表面法向量来代替顶点的表面法向量，计算光照量。这样可减小中央处理器(CPU)向图形处理器(GPU)传输的数据量，确保实时的性能。

优选在步骤(5)进行重采样的同时还进行图像后处理，图像后处理包括平滑、增强、灰阶变换等，可在图形处理器中与重采样一起实现，以节省硬件的资源。

平滑、增强的实现有多种方法，可以在空间域或频率域上进行，可以针对不同的组织采用不同的方法，比较常用的三种方法是：中值法，最小方差法，低通滤波法。在重采样的时候，不仅对当前点进采样，还需对其周边相关的点进行采样，这样得到一系列与当前点相关的点，然后通过各种平滑方法对其进行运算，达到平滑或增强的效果；这些平滑方法本身并不复杂但计算量很大，因此将平滑方法整合进图形处理器(GPU)进行计算，达到实时平滑或增强的效果。例如，在频域上实现平滑算法，一般需作快速傅立叶(FFT)转换，对于体积图像数据来说计算量非常大，因此可用图形处理器(GPU)来实现快速傅立叶(FFT)转换，再根据需要对其进行高通、低通或者其它滤波，运算速度可提高8～10倍；又如，采用空间域的中值法实现平滑算法，对任一个当前点来说(边缘上的点需特殊处理)，先计算出它在三个坐标轴上的相邻点的偏移量，以512×256×128的体积图像数据来说，它在三个方向上的偏移量分别为1/512、1/256、1/128，根据这个偏移量计算得到当前点周围相邻26个点的纹理坐标偏移量，然后对这些纹理坐标分别进行三维纹理采样，得到各点的采样值，再进行中值的查找来代替当前点的采样值。

灰阶变换相当于灰阶拉伸或者压缩，把需要变换的灰阶函数转换成查找表并通过纹理送入GPU，对重采样后的值通过该查找表来实现变换功能。

本发明提供一种基于图形处理器的体绘制方法，既不丢失原始图像信息，能够反映整个原始体积图像数据的全貌及细节，又可对可编程硬件图形处理器进行编程，充分利用图形处理器的高速运算能力和高度并行性，对体积图像数据进行快速运算，实现图像的实时显示，从而实现三维医学超声图像的实时重建，而且，通过三维坐标变换和重采样，能够使体积图像数据及其数据点与被测物体相符，从而使最后形成的图像能够真实反映被测物体的实际形态。另外，由于三维医学超声图像的重建过程交由图形处理器进行处理，中央处理器不再负担三维医学超声图像重建过程中繁重的运算任务，因此中央处理器的资源得到释放，可以更专心的处理控制、分支等比较擅长的任务，使硬件资源都得到合理的安排与利用。

附图说明

图1是超声三维传感器的成像网格示意图；

图2是本发明优选实施例实时三维医学超声图像的重建方法的流程图；

图3是本发明另一优选实施例三维坐标变换过程的示意图。

具体实施方式

如图2所示，这种实时三维医学超声图像的重建方法依次包括下述步骤：

(1)对来自超声三维传感器的体积图像数据进行预处理，获得重建图像所需的体积图像数据；本步骤可在中央处理器CPU中进行。

超声三维传感器得到的数据一般为16位或者16位以上，但屏幕只能显示8位灰阶图像，这就需要我们把原始体积图像数据进行格式转换(将16位或者16位以上的数据转换为8位数据)，使其符合屏幕的显示要求。

(2)对步骤(1)获得的体积图像数据的数据点进行分类，将反映同一物质的数据点归入同一类；本步骤可在中央处理器CPU中进行。

对体积图像数据的数据点进行分类可采用下述方法：根据体积图像数据所在应用领域的背景知识，或者对全部数据点进行统计，设定若干阈值d_m(m＝1，2...，n)；如果各数据点以f(x_i，y_j，z_k)表示，则将满足条件d_m≤f(x_i，y_j，z_k)＜d_m+1的数据点归入同一类。

(3)对数据点赋值；包括：

(3-1)赋颜色值(本步骤可在中央处理器CPU中进行)，对一个预先设定好的颜色表(颜色表可根据对当前的被探测的物质进行统计分析计算后得出)进行查找；对同一类数据点赋上相同的颜色，即同一物质在图像中采用一种颜色表示，颜色值则对应颜色的饱和度。由于是人为进行赋值，得到的最终图像的颜色将是伪彩色。

(3-2)赋不透明度值(本步骤可在中央处理器CPU中进行)，对一个预先设定好的不透明度表进行查找。不透明度表示数据点对光的吸收度，值为1时表示完全不透明，值为0时表示完全透明。

(3-3)通过图形处理器(GPU)的可编程处理器，构建光照模型(采用CG(C for graphics)显卡编程语言表达光照模型)，进行光照量的计算，并将光照量与数据点的颜色值进行混合，从而模拟出较逼真的光照效果。

光照量的计算中还具有明暗计算，明暗计算的步骤如下：

(a)计算各数据点的梯度值(梯度是指灰阶的变化程度)；

G_{f} (x_{i}, y_{i}, z_{i}) = [\begin{matrix} \frac{1}{2} (f (x_{i - 1}, y_{i}, z_{i}) - f (x_{i + 1}, y_{i}, z_{i})) \\ \frac{1}{2} (f (x_{i}, y_{i - 1}, z_{i}) - f (x_{i}, y_{i + 1}, z_{i})) \\ \frac{1}{2} (f (x_{i}, y_{i}, z_{i - 1}) - f (x_{i}, y_{i}, z_{i + 1})) \end{matrix}]

(b)将步骤(a)得到的梯度值作为纹理保存进图形处理器(GPU)中；该部分的计算机源程序(用OpenGL语言表达)如下：

glGenTextures(1，&GL_TiduTexture)；//申请创建纹理空间并设定纹理参数

glBindTexture(GL_TEXTURE_3D，GL_TiduTexture)；

glTexParameteri(GL_TEXTURE_3D，GL_TEXTURE_WRAP_S，GL_CLAMP)；

glTexParameteri(GL_TEXTURE_3D，GL_TEXTURE_WRAP_T，GL_CLAMP)；

glTexParameteri(GL_TEXTURE_3D，GL_TEXTURE_WRAP_R，GL_CLAMP)；

glTexParameteri(GL_TEXTURE_3D，GL_TEXTURE_MAG_FILTER，L_LINEAR)；

glTexParameteri(GL_TEXTURE_3D，GL_TEXTURE_MIN_FILTER，GL_LINEAR)；

glTexImage3DEXT(GL_TEXTURE_3D，0，GL_RGBA，iWidth，iHeight，iDepth，0，

GL_RGBA，GL_UNSIGNED_BYTE，NULL)；

glBindTexture(GL_TEXTURE_3D，GL_TiduTexture)；//绑定梯度3D纹理

glTexSubImage3DEXT(GL_TEXTURE_3D，0，0，0，0，iWidth，iHeight，iDepth，

GL_LUMINANCE，GL_UNSIGNED_BYTE，TiduPtr)；

(c)图形处理器的可编程顶点处理器接收到纹理后，将纹理传递给可编程片段处理器；可编程片段处理器接收纹理后，用纹理点的表面法向量来代替顶点的表面法向量，计算光照量。

用CG语言表达的光照模型如下：

surfaceColor＝emissive+ambient+diffuse+specular

+attenuationFactor；

其中surfaceColor表示光照量，光照量surfaceColor是放射光分量emissive、环境光分量ambient、漫反射光分量diffuse、镜面反射光分量specular和距离衰减量attenuationFactor等光照分量之和；各分量具体描述如下：

3)漫反射光分量diffuse代表从一个表面相等地向所有方向反射出去的方向光，通常物体表面微观尺寸上是非常粗糙的，当入射光线到达物体表面时，光线会向各个方向反射，表示为diffuse＝K_d*lightColor*max(N·L，0)，其中K_d是物体的漫反射颜色，lightColor是入射漫反射光的颜色，N是规范化的表面法向量(用纹理点的表面法向量来代替顶点的表面法向量，进行计算)，L是规范化的指向光源的向量；

specular＝K_s*lightColor*facing*(max(N·H)，0)^shininess，

其中K_s是物体材质的镜面反射颜色，lightColor是入射镜面反射光的颜色，N是规范化的表面法向量(用纹理点的表面法向量来代替顶点的表面法向量，进行计算)，H是视点向量与光源向量的点积，shininess表示物体的光泽度；

5)距离衰减量attenuationFactor表示为attenuationFactor＝1/(K_c+K_l*d+k_q*d*d)，其中d是物体到光源的距离，K_c、K_l、k_q为控制衰减量的常量；

在图形处理器(GPU)中的可编程处理器中对各数据点的各个光照分量进行计算后，对各个光照分量进行累加，得到各数据点的光照量。

(4)中央处理器将体积图像数据传输给图形处理器，在图形处理器中对体积图像数据进行三维坐标变换；

三维坐标变换具体可利用图形处理器(GPU)中的可编程处理器，依照下述步骤对体积图像数据进行纹理坐标变换(纹理坐标变换是三维坐标变换中的一种)：

在可编程片段处理器中，在[0，1]区间外的数据属于冗余数据或者重复数据，经过坐标变换后，可以利用剪切面对纹理坐标[0，1]外的部分进行剪切，剪切后的数据不参加重采样的计算，从而提高运算性能。

(5)在图形处理器中进行重采样；

经过三维坐标变换后，数据点的坐标通常会发生变化，也就是说，数据点在双曲面数据中的坐标与其在立方体数据中的原始坐标不同，根据数据点的原始坐标不能在双曲面数据中找到相应数据点的数据值(包括颜色值、不透明度值和光照量)，因此必须进行重采样，将数据点的原始坐标转换成数据点在双曲面数据中的坐标。重采样是在图形处理器(GPU)的可编程片段处理器中进行的，可编程片段处理器根据数据点的原始坐标计算得到数据点在双曲面数据中的坐标，并一一对应地用数据点在双曲面数据中的坐标替换原始坐标。每个数据点对片段处理器来说相当于一个运算核，片段处理器对各运算核进行并行计算，从而提高运算速度以达到实时的目的。

进行重采样之后，将各数据点的光照量与颜色值进行加权求和(光照量与颜色值各自的权重可根据需要设定，当要求图像较亮时光照量的权重较大，当要求图像较暗时光照量的权重较小)，可得出各数据点的光亮度值，从而使最后在屏幕上显示的图像具有较理想的光照效果。

本实施例用于获得光亮度值(包括光照量的计算及其重采样、光照量与颜色值的混合)的部分计算机源程序(用CG语言表达)如下：

float3 ambient＝Ka*globalAmbient； //环境光分量

float3 emissive＝Ke； //放射光分量

//得到重采样的梯度值，进行明暗度计算

float4 normal4＝tex3Dproj(TiduTexture，inTex.yzxw)；

float3 N＝normal4.xyz； //使用纹理坐标代替顶点坐标

N＝N*2-1.0； //放大梯度

N＝normalize(N)； //规范化的表面法向量

float3 L＝normalize(lightPosition-N)；//片段位置到光源的规范化向量

float diffuse＝dot(oNormal，L)；//漫反射光分量

float3 V＝normalize(eyePosition-N)；//视点到片段位置的规范化向量

float3 H＝normalize(L+V)； //规范化V与L的半角向量

float specular＝pow(max(dot(N，H)，0)，shininess)；//镜面反射光分量

float4 light＝lit(diffuse，specular，shininess)；//综合计算

float4 lightcolor＝float4((emissive+ambient+Kd*light.y+

Ks*light.z)，1.0)；

inTex.yz＝inTex.yz*inTex.w； //预先把坐标乘以缩放因子

float val＝tex3Dproj(USTexture，inTex.yzxw)；//三维重建采样

sColor0＝tex1D(ColorMap，val)；//颜色映射1D纹理采样

//光亮度值输出

sColor0＝float4(sColor0.x*lightcolor.x，sColor0.y*lightcolor.y，

sColor0.z*lightcolor.z，sColor0.w)；

(6)图像后处理；

图像后处理包括平滑、增强、灰阶变换等，可在图形处理器中与重采样一起实现，以节省硬件的资源。

平滑、增强的实现有多种方法，可以在空间域或频率域上进行，可以针对不同的组织采用不同的方法，比较常用的三种方法是：中值法，最小方差法，低通滤波法。在重采样的时候，不仅对当前点进采样，还需对其周边相关的点进行采样，这样得到一系列与当前点相关的点，然后通过各种平滑方法对其进行运算，达到平滑或增强的效果；这些平滑方法本身并不复杂但计算量很大，因此将平滑方法整合进图形处理器(GPU)进行计算，达到实时平滑或增强的效果。例如，在频域上实现平滑算法，一般需作快速傅立叶(FFT)转换，对于体积图像数据来说计算量非常大，因此可用图形处理器(GPU)来实现快速傅立叶(FFT)转换，再根据需要对其进行高通、低通或者其它滤波，运算速度可提高8～10倍；又如，采用空间域的中值法实现平滑算法，对任一个当前点来说(边缘上的点需特殊处理)，先计算出它在三个坐标轴上的相邻点的偏移量，以512×256×128的体积图像数据来说，它在三个方向上的偏移量分别为1/512、1/256、1/128，根据这个偏移量通过三角函数公式计算计算得到当前点周围相邻26点的纹理坐标偏移量，然后对这些纹理坐标分别进行三维纹理采样，得到各点的采样值，再进行中值的查找来代替当前点的采样值。

(7)图像合成；图像合成在图像处理器中进行。

图像合成一般采用光线投射法的原理进行，其目的是沿某一像素点发出一条射线，计算该射线上各数据点的光亮度值及不透明度值，并进行合成，形成该像素点的最终光亮度与最终透明度。合成的原理是：体积图像数据可以看成是多个二维切片组成的，这些二维切片经过坐标变换后投影到显示屏幕上，则多个切片上的点会映射到显示屏幕上的同一个点，显示屏幕上的这个点的最终显示颜色与最终透明度由多个二维切片上的点来混合决定。

本实施例采用显卡生产厂商NVIDIA提供的图形处理器(GPU)，并利用NVIDIA的CG(C for graphics)显卡编程语言与OPENGL三维图形接口对图形处理器进行编程。

在另一具体方案中，如图3所示，三维坐标变换可分两步进行，首先立方体数据A经过第一次变换先形成锥形数据B，然后再经过第二次变换后形成双曲面数据C。具体是利用图形处理器(GPU)中的可编程处理器，依照下述步骤，对体积图像数据进行纹理坐标变换(纹理坐标变换是三维坐标变换中的一种)：

②选取XY轴形成的平面，依据超声三维传感器的参数信息(如探头半径、探头角度、扫描线数、显示角度、显示点数等)，对Y方向数据进行缩放，在XY平面上形成梯形；利用步骤①保存的纹理坐标来对XZ方向进行坐标变换，首先对Z方向数据进行缩放，在XZ平面上形成梯形；经上述坐标变换后形成锥形数据B；

③求出对当前Y对应X变换后形成的角度，对X方向进行缩放，并进行平移，产生弧形，就完成了XY方向上的坐标变换；求出对当前Z对应X变换后形成的角度，对X方向进行缩放，并进行平移，产生弧形，就完成了XZ方向上的坐标变换；经上述坐标变换后形成双曲面数据C；

Claims

1、一种实时三维医学超声图像的重建方法，其特征在于依次包括下述步骤：

(3)对数据点赋值，包括对数据点赋颜色值和不透明度值；

(4)在图形处理器中对体积图像数据进行三维坐标变换，将立方体数据变换成双曲面数据；

(5)在图形处理器中进行重采样，将数据点的原始坐标转换成数据点在双曲面数据中的坐标；

(6)图像合成。

2、根据权利要求1所述的实时三维医学超声图像的重建方法，其特征是：步骤(1)中，对来自超声三维传感器的体积图像数据的预处理包括原始体积图像数据的格式转换、剔除冗余数据，并导出所需的体积图像数据。

3、根据权利要求1所述的实时三维医学超声图像的重建方法，其特征是：步骤(2)中，对体积图像数据的数据点进行分类采用下述方法：根据体积图像数据所在应用领域的背景知识，或者对全部数据点进行统计，设定若干阈值d_m，其中m为正整数；如果各数据点以f(x_i，y_j，z_k)表示，则将满足条件d_m≤f(x_i，y_j，z_k)＜d_m+1的数据点归入同一类。

4、根据权利要求1所述的实时三维医学超声图像的重建方法，其特征是：步骤(4)中，三维坐标变换分两步进行，首先立方体数据经过第一次变换先形成锥形数据，然后再经过第二次变换后形成双曲面数据。

5、根据权利要求4所述的实时三维医学超声图像的重建方法，其特征是：步骤(4)中进行三维坐标变换，具体是，利用图形处理器中的可编程处理器，依照下述步骤对体积图像数据进行纹理坐标变换：

②选取XY轴形成的平面，依据超声三维传感器的参数信息，对Y方向数据进行缩放，在XY平面上形成梯形；利用步骤①保存的纹理坐标来对XZ方向进行坐标变换，首先对Z方向数据进行缩放，在XZ平面上形成梯形；经上述坐标变换后形成锥形数据；

④根据超声三维传感器的参数信息重新计算体积图像数据中心。

6、根据权利要求1所述的实时三维医学超声图像的重建方法，其特征是：步骤(4)中进行三维坐标变换，具体是，利用图形处理器中的可编程处理器，依照下述步骤对体积图像数据进行纹理坐标变换：

②选取XY轴形成的平面，依据超声三维传感器的参数信息，对Y方向数据进行缩放，在XY平面上形成梯形；

⑥根据超声三维传感器的参数信息重新计算体积图像数据中心。

7、根据权利要求1所述的实时三维医学超声图像的重建方法，其特征是：步骤(5)中，重采样是在图形处理器的可编程片段处理器中进行的，可编程片段处理器根据数据点的原始坐标计算得到数据点在双曲面数据中的坐标，并一一对应地用数据点在双曲面数据中的坐标替换原始坐标。

8、根据权利要求1～7任一项所述的实时三维医学超声图像的重建方法，其特征是：通过图形处理器的可编程处理器，构建光照模型，在步骤(3)对数据点赋值后，对各数据点进行光照量的计算；光照模型用CG语言表达如下：

surfaceColor＝emissive+ambient+diffuse+specular；

其中surfaceColor表示光照量，光照量surfaceColor是放射光分量emissive、环境光分量ambient、漫反射光分量diffuse和镜面反射光分量specular之和；各分量具体描述如下：

1)放射光分量emissive用K_e表示，即emissive＝K_e；

2)环境光分量ambient表示为ambient＝K_a*globalAmbient，其中K_a为环境光分量的系数，globalAmbient为入射环境光的颜色；

3)漫反射光分量diffuse表示为

diffuse＝K_d*lightColor*max(N·L，0)，

其中K_d是物体的漫反射颜色，lightColor是入射漫反射光的颜色，N是规范化的表面法向量，L是规范化的指向光源的向量；

4)镜面反射光分量specular表示为

specular＝K_s*lightColor*facing*(max(N·H)，0)^shininess，

其中K_s是物体材质的镜面反射颜色，lightColor是入射镜面反射光的颜色，N是规范化的表面法向量，H是视点向量与光源向量的点积，shininess表示物体的光泽度；

在步骤(5)进行重采样之后，将各数据点的光照量与颜色值进行加权求和，得出各数据点的光亮度值。

9、根据权利要求8所述的实时三维医学超声图像的重建方法，其特征是：计算各数据点的光照量时，还计算距离衰减，距离衰减量attenuationFactor表示为attenuationFactor＝1/(K_c+K_l*d+k_q*d*d)，其中d是物体到光源的距离，K_c、K_l、k_q为控制衰减量的常量；光照量surfaceColor是放射光分量emissive、环境光分量ambient、漫反射光分量diffuse、镜面反射光分量specular和距离衰减量attenuationFactor之和。

10、根据权利要求8所述的实时三维医学超声图像的重建方法，其特征是：所述对各数据点光照量的计算中还具有明暗计算，明暗计算的步骤如下：

(a)计算各数据点的梯度值，梯度值在中央处理器中进行；

(b)将步骤(a)得到的梯度值作为纹理保存进图形处理器中；

(c)图形处理器的可编程顶点处理器接收到纹理后，将纹理传递给可编程片段处理器；可编程片段处理器接收纹理后，用纹理点的表面法向量来代替顶点的表面法向量，作为规范化的表面法向量N，计算光照量。