CN101606854A - 一种高精度实时超声图像扫描变换方法 - Google Patents

Abstract

本发明公开了一种高精度实时超声图像扫描变换方法，使用基于GPU的并行计算架构和SIMD矢量运算，通过绘制合适的代理几何模型，将扫描变换模块中的坐标变换和图像插值算法全部放在GPU中作32位浮点数高精度计算，同时使用优化的高精度的双三次插值算法，获得更平滑的图像质量。本发明可以在目前普通中低档显卡上，对512条线数据，每条线上1024个采样点这样大规模的数据，达到每秒钟60帧以上的扫描变换速度，可以充分挖掘利用当代GPU所提供的强大的并行计算能力，从而将传统的扫描变换模块升级为下一代高保真的扫描变换模块。

Description

一种高精度实时超声图像扫描变换方法

技术领域

本发明涉及一种高精度实时超声图像扫描变换方法，用于医用超声诊断设备。

背景技术

超声波回波成像技术目前已经被广泛应用于军事、医疗等领域，通过向目标区域发射超声波，然后使用接收装置接收反射回来的回波信号，并通过信号处理技术和图像处理技术，抑制回波信号中的无用部分，最终形成目标区域的图像。

在跟我们每个人的日常生活息息相关的医疗领域，超声波回波成像技术更是获得了长足的发展，目前各种医用超声诊断设备如B超等已经广泛应用于各个医院的临床诊断中，除了传统的黑白超可以观察病人的内部组织和器官的解剖结构外，彩超通过使用多普勒效应可以对血管内的血流成像，大大提高了超声诊断设备的临床应用范围。

在目前的超声诊断设备中，扫描变换是其中非常重要的一个模块。超声诊断设备上目前大部分使用的探头，一次是沿着某个方向发射和接收超声数据的，通过多次发射和接收，形成了一帧完整的图像。对于常用的凸阵探头，采集到的数据形成了一个扇形的几何形状，如附图1所示，探头第一次沿着线段AB的方向发射超声数据，然后接收的时候在线段AB上进行等间隔采样，图中的空心点即为采样点。探头第二次沿着CD方向发射、接收并采样，第三次沿着EF方向放射、接收并采样，第四次沿着GH方向发射、接收并采样。这四个方向的采样数据，形成了一帧图像，其几何形状为一个扇形ACEGHFDB，圆心在O点。但是对于目前的显示设备而言，都是光栅设备，只能显示一个一个离散的像素点，并且这些像素点只能排列在矩形网格上，如图中的实心点所示。因此要把超声诊断设备采集到的扇形数据，显示在PQRS这样一个矩形区域里面，必须通过一个变换过程，依据已知的空心点采样值，计算出每个实心点处的像素值，这个变换过程叫做扫描变换。

扫描变换分为两个部分，一个部分为极坐标和直角坐标之间的坐标变换，这一部分主要的运算量在于每次变换都需要计算反三角函数和开根号函数；另外一个部分为图像插值，因为空心点和实心点的坐标大部分情况下都不重合，因此需要通过插值算法来从多个空心点的采样值计算一个实心点的像素值。传统上来讲，扫描变换是超声诊断设备中非常耗费资源的计算密集型的模块，为了简化起见，图1中假设只扫描了四次，或者说只获取了四条线数据，每条线上只采样了四个点，对于实际的情况，往往是要采集512线以上，每条线上采样点的个数也在512个以上。对于这么大的数据量，计算坐标变换和图像插值需要耗费很大量的计算资源，而这些计算又必须实时地完成才能达到临床检查的要求。为了满足实时扫描变换的需求，大部分超声诊断设备上都专门用FPGA或者DSP芯片来完成扫描变换，但是受硬件资源的限制，这些计算都是采用定点数计算，反三角函数和开根号函数等都是用的近似算法或者查表算法，计算精度受到很大的影响，另外图像插值算法用的是最近邻插值算法或者双线性插值算法，更高阶的插值算法因为运算量太大，无法达到实时的计算速度。

从2000年以来，随着三维游戏产业的迅猛发展，普通的民用显卡的计算能力越来越强，并且允许开发人员使用类似C语言的语法直接对显卡进行编程，如美国NVidia公司的Cg、CUDA编程语言，Microsoft公司的HLSL编程语言，OpenGL架构委员会的GLSL编程语言等。跟CPU相类似，显卡的计算核心被叫做GPU(Graphics Processing Unit)，目前GPU除了应用于三维游戏当中的计算之外，也越来越多地被应用于通用计算。西门子公司的专利“Volume Renderingin the Acoustic Grid Methods and Systems for Ultrasound Diagnostic Imaging”(美国专利号：6852081)提出了一种基于GPU的扫描变换方法，首先在CPU中生成一个几何模型，把AC、CE、EG、HF、FD、DB用直线连接起来，用多边形ACEGHFDB来近似地逼近扇形ACEGHFDB，如图1所示，多边形ACEGHFDB叫做代理几何模型。AB、CD、EF、GH四条线上的采样数据被作为一个二维纹理，传输给GPU的纹理显存。代理几何模型上的每个顶点都被赋予一个纹理坐标，然后通过GPU内部硬件实现的纹理映射和双线性插值，来完成扫描变换过程。虽然这个方法可以在GPU上实现实时的扫描变换，但是仍然存在一些不足：第一，其代理几何模型只是对原始扇形的一个逼近，会造成一定的计算误差；第二，只是代理几何模型的顶点(A、C、E、G、H、F、D、B)作了精确的坐标变换计算，代理几何模型内部的点均是通过GPU内部的纹理坐标插值机制来计算坐标变换的，是一个近似算法；第三，代理几何模型的顶点坐标和纹理坐标仍然需要CPU来计算，其中会牵涉到反三角函数和开根号函数，当采样的线数和帧频非常高时，会对CPU造成一定的计算负担；第四，图像插值是利用的GPU内置的双线性插值算法，插值所得到的函数曲面只是分段连续的，精度不够高。

考虑到上述问题，在超声诊断设备上提供一种高精度的，能够实时完成计算的扫描变换方法，是非常有意义的。

发明内容

本发明的目的是解决目前超声诊断设备中所使用的扫描变换模块精度不够高的问题，提供一种高精度实时超声图像扫描变换方法，使用基于GPU的计算架构，通过绘制合适的代理几何模型，将扫描变换模块中的坐标变换和图像插值算法全部放在GPU中作32位浮点数高精度计算，同时使用高精度的双三次插值算法，获得优化的图像质量。

按照本发明提供的技术方案，一种高精度实时超声图像扫描变换方法步骤如下：

步骤1：获取超声探头所采集的每帧图像，由LN条线数据组成，形成一个扇形，内径为r，外径为r+l，第1条线数据和第LN条线数据所夹的圆心角为2θ，LN条线数据为等角度采样，每条采样线上有PN个采样点，为等距离采样。将所述每帧图像作为二维纹理传输至GPU显存，并保存为线数据纹理；

步骤2：生成代理几何模型并绘制，包括：

使用六边形ABQRHG作为代理几何模型，其顶点坐标分别为：

(A_x，A_y)＝(-rsinθ，rcosθ)

(B_x，B_y)＝(-(r+l)sinθ，(r+l)cosθ)

(Q_x，Q_y)＝(-(r+l)sinθ，r+l)

(R_x，R_y)＝((r+l)sinθ，r+l)

(H_x，H_y)＝((r+l)sinθ，(r+l)cosθ)

(G_x，G_y)＝(rsinθ，rcosθ)

其中A和G、B和H、Q和R的x坐标相反，y坐标相同，因此只需计算A、B、Q三个点的坐标，其它三个点的坐标根据上述关系求得；

将六边形ABQRHG使用三角带形式分解为四个三角形，顶点顺序为G，H，A，R，B，Q，三角形分别为：GHA、HAR、ARB和RBQ；

使用OpenGL或者DirectX的API函数绘制三角带；

步骤3：在GPU中对代理模型绘制后的每个光栅点(Fragment)计算坐标变换，包括：

将代理模型绘制后的每个光栅点从显示器的直角坐标(x，y)转换为所述线数据纹理所在的极坐标(s，t)，计算公式为： $s=(\sqrt{x^2+y^2}-r)(\mathrm{PN}-1)/l,$ t＝(arctan(x/y)+θ)(LN-1)/2θ，其中t代表在第t条线数据上，s代表在线数据的第s个采样点上；

判断(s，t)是否满足0≤t≤(LN-1)且0≤s≤(PN-1)，如果不满足则终止计算；

步骤4：在GPU中计算双三次图像插值，包括：

在-2＜u＜2的区间上取NX个离散的三次多项式f(u)的非零值，NX为2的正整数次方，用一维纹理传输至GPU显存，并保存为三次多项式纹理，其中f(u)的定义如下：

$f\left(u\right)=\left\{\begin{array}{cc}\frac{7}{6}{\left|u\right|}^3-2{\left|u\right|}^2+\frac{8}{9},& \left|u\right|<1\\ -\frac{7}{18}{\left|u\right|}^3+2{\left|u\right|}^2-\frac{10}{3}\left|u\right|+\frac{16}{9},& 1\&le;\left|u\right|<2\\ 0,& \left|u\right|\&GreaterEqual;2\end{array}\right.;$

分别使用纹理坐标(floor(s)-1，floor(t)-1)，(floor(s)，floor(t)-1)，(ceil(s)，floor(t)-1)，(ceil(s)+1，floor(t)-1)访问所述线数据纹理，得到4个长度为4的矢量v1，v2，v3，v4，其中，函数floor(v)代表不大于v的最大的整数，函数ceil(v)代表不小于v的最小的整数；

分别使用纹理坐标floor(t)-t-1和floor(s)-s-1访问所述三次多项式纹理，得到2个长度为4的矢量f1，f2；

利用GPU的SIMD矢量计算指令，计算得到一个长度为4的矢量temp＝(<v1·f1>，<v2·f1>，<v3·f1>，<v4·f1>)，<·>表示矢量点积计算；

利用GPU的SIMD矢量计算指令，计算待求的采样值I_val＝<temp·f2>；

最后，将采样值I_val输出至显示器的(x，y)坐标处显示。

所述线数据纹理采用RGBA四个通道，每个通道为32位浮点数的格式；对于第i条线上第j个采样点，其R通道存放(i，j)处的采样值，G通道存放(i+1，j)处的采样值，B通道存放(i+2，j)处的采样值，A通道存放(i+3，j)处的采样值。

所述三次多项式纹理采用RGBA四个通道，每个通道为32位浮点数的格式；对于在区间-2＜u＜2的u，R通道存储f(u)的值，G通道存储f(u+1)的值，B通道存储f(u+2)的值，A通道存储f(u+3)的值。

本发明的优点是：本发明使用目前普通民用显卡提供的GPU计算架构，可以对512条线数据，每条线上1024个采样点这样大规模的数据，达到每秒钟60帧以上的扫描变换速度。相比于已有的GPU方法，本发明使用了更简化的代理几何模型，并且无需指定每个顶点的纹理坐标，因此大大降低了CPU的负担，而将所有的计算全部放在GPU中以最高精度完成，没有使用任何的近似算法。同时，图像插值算法并没有利用GPU内置的双线性插值算法，而是使用更高精度的双三次插值算法，保证了最终图像的质量。本发明所提供的方法可以充分挖掘利用当代GPU所提供的强大的并行计算能力，从而将传统的扫描变换模块升级为下一代高保真的扫描变换模块。

附图说明

图1是传统扫描变换的示意图。

图2是本发明的实施流程图。

图3是代理几何模型的示意图。

图4是代理几何模型的坐标系及参数示意图。

图5是双线性插值的示意图。

图6是双三次插值的示意图。

具体实施方式

下面结合附图和实施例详细说明本发明技术方案中所涉及的各个细节问题。应指出的是，所描述的实施例仅旨在便于对本发明的理解，而对其不起任何限定作用。

如图2所示，本发明的实施流程分为4个步骤。在步骤1中，首先将探头采集到的线数据传输至GPU中的纹理显存。探头一次沿着某个方向发射和接收超声数据，在一个方向上接收到的数据叫做一条线数据，比如图1中的线段AB。一条线数据可以用一个一维数组在计算机中表示，假设一条线上采样点的个数是PN，那么这个一维数组的长度就是PN。假设形成一帧图像需要采集的线数为LN，则总共得到LN个一维的数组，或者说一个二维的数组，这个二维数组的大小为PN×LN。假设我们用lineArray来代表这个二维数组，则lineArray[i][j]表示第i条线上第j个采样点的数值，这里使用基于0的索引，即0≤i≤(LN-1)，0≤j≤(PN-1)。以图1为例，因为只有四条采样线，所以LN＝4，每条采样线上有四个采样点，所以PN＝4，lineArray[1][2]代表CD采样线上T采样点的数值。为了进行扫描变换，必须首先把lineArray作为一个二维纹理传输到GPU的纹理显存中，然后GPU的计算单元才能访问到它。目前对GPU进行编程的工具主要有OpenGL和DirectX，可以调用它们所提供的API函数来完成上述的纹理传输工作。

需要注意的是，本发明中使用RGBA四通道的纹理格式，每通道是32位的浮点数，这个一方面是为了保证高精度的计算，另外一方面也是为了加快图像插值算法的速度，具体会在步骤4中给出详细的说明。

在步骤2中，生成代理几何模型，然后通过OpenGL或者DirectX等图形API的绘制命令将代理几何模型传送至GPU，由GPU对其进行光栅化操作，生成一个一个的光栅片段(Fragment)，再由GPU的Fragment Shader对每一个Fragment进行处理。代理几何模型的选择非常重要，不能过于复杂，否则会对CPU造成很大的计算代理几何模型顶点坐标的负担，也会导致通过AGP总线或者PCIE总线给GPU传输数据太多，从而形成传输瓶颈。与之前的技术不同，本发明中使用非常简单的代理几何模型，如图3所示，六边形ABQRHG完全包含了ABHG扇形成像区域，可以选择它作为代理几何模型。直线QR是BH弧的切线，直线BQ和HR垂直于直线QR。生成六边形ABQRHG只需计算六个顶点的坐标，并且顶点A和G、B和H、Q和R都是关于六边形的中心线对称的，更进一步简化了CPU所需的计算。因为GPU对由三角形组成的几何模型的绘制效率最高，因此还需把六边形分解成一系列的三角形。为了减少传输给GPU的几何数据，这里使用三角带(Triangle Strip)来对六边形进行分解，所谓三角带，就是指的一系列的三角形，相邻两个三角形必须共享一条边，或者说两个顶点。三角带对顶点顺序是有规定的，假设给定1，2，3，......，n+2总共n+2个顶点，则1，2，3三个顶点组成一个三角形，2，3，4三个顶点组成一个三角形，依次类推，直到n，n+1，n+2三个顶点组成一个三角形，因此可以用n+2个顶点来表示n个三角形(n≥1)，而如果用普通的三角形表达方式，则需要3n个顶点才能表示n个三角形。分别连接AH、AR、RB三条直线段，将六边形分解为四个三角形，用三角带的形式表达，顶点顺序为G，H，A，R，B，Q，分别代表三角形GHA、HAR、ARB和RBQ。

用三角形来对代理几何模型进行分解后，还需要计算出每个顶点的坐标，如图4所示，首先建立坐标系，坐标原点定在扇形ABHG的圆心O处，x轴水平从左到右，y轴竖直从上到下。直线段OA的长度为r，AB的长度为l，OB与y轴正向的夹角为θ，这些参数都是已知值。从上面的这些条件，可以计算出顶点A的坐标为：(A_x，A_y)＝(-rsinθ，rcosθ)，顶点G的坐标为：(G_x，G_y)＝(rsinθ，rcosθ)，顶点B的坐标为：(B_x，B_y)＝(-(r+l)sinθ，(r+l)cosθ)，顶点H的坐标为：(H_x，H_y)＝((r+l)sinθ，(r+l)cosθ)，顶点Q的坐标为：(Q_x，Q_y)＝(-(r+l)sinθ，r+l)，顶点R的坐标为：(R_x，R_y)＝((r+l)sinθ，r+l)。从上述坐标计算公式可以看出，实际上只需计算A、B、Q三个顶点的坐标，G、H、R分别和A、B、Q关于y轴对称，因此无需计算，另外B和Q、H和R的x坐标相等，又可以节省昂贵的三角函数计算。计算得到代理模型的顶点坐标之后，就可以使用OpenGL或者DirectX的API函数将代理几何模型绘制并交给GPU进行后续处理。

与之前已有的技术相比，本发明所提出方法除了使用的代理模型更为简单，计算量更小之外，坐标变换的精度也更高。传统的技术在生成代理几何模型时，对于每个顶点都要计算并分配一个纹理坐标，这一方面增加了CPU的计算负担，另外一方面依赖纹理坐标的双线性插值去计算坐标变换，是一个粗糙的逼近。本发明中在生成代理几何模型时，并没有给顶点分配纹理坐标，而是在步骤3中，直接在GPU的Fragment Shader中用32位浮点数进行高精度的坐标变换计算。需要注意的是，步骤1和步骤2都是在CPU中完成的，从步骤3开始，所有的计算都是在GPU中完成。通过步骤2，六边形ABQRHG被绘制成为一系列离散的像素点，这些像素点的集合组成了六边形ABQRHG的内部和边界，不过在屏幕的分辨率下肉眼看是连续的。其中每个像素点叫做Fragment，记录有自己的坐标、颜色等信息，会送给GPU的Fragment Shader作进一步的处理。GPU的Fragment Shader是一个可编程的器件，可以使用GLSL、Cg、HLSL或者CUDA对其进行编程，以达到自己特定的计算目的。每个Fragment所记录的坐标是在屏幕坐标系的(x，y)坐标，是直角坐标系，而探头所采集的线数据用极坐标系更容易表达，如图5和图1，任取一个Fragment，假设为I，其直角坐标为(x，y)，线段OI的长度为ρ，OI与y轴正向的夹角为φ，则I可以用极坐标(ρ，φ)来表示，其中 $\mathrm{\&rho;}=\sqrt{x^2+y^2},$ φ＝arctan(x/y)，r≤ρ≤(r+l)，-θ≤φ≤θ。另外，还必须计算出I是处于第几条采样线上，是在采样线的第几个采样点上。假设采集一帧图像的总采样线的条数为LN，一条采样线上采样点的个数是PN，I处于第t条采样线的第s个采样点，0≤t≤(LN-1)，0≤s≤(PN-1)，则从极坐标(ρ，φ)可以很容易地推出来(s，t)的取值：s＝(ρ-r)(PN-1)/l，t＝(φ+θ)(LN-1)/2θ，最终得到的(s，t)和(x，y)之间的坐标变换关系为： $s=(\sqrt{x^2+y^2}-r)(\mathrm{PN}-1)/l,$ t＝(arctan(x/y)+θ)(LN-1)/2θ。

通过上述坐标变换计算得到的s和t不一定是整数，但是探头实际采样得到的数据，s和t必须是整数，因此需要使用图像插值算法得到任意(s，t)坐标处的采样值。传统的技术为了计算速度，往往采用最近邻插值或者双线性插值算法。所谓最近邻插值算法，是对s和t分别四舍五入，取最近的整数坐标处的采样值，这种方法速度最快，但是图像质量也最差。用的最为广泛的是双线性插值算法，在I的周围取最近的四个邻域采样点J、K、L、M，其(s，t)坐标分别为：(floor(s)，floor(t))、(ceil(s)，floor(t))、(ceil(s)，ceil(t))、(floor(s)，ceil(t))，其中floor(x)函数代表不大于x的最大的整数，ceil(x)函数代表不小于x的最小的整数。I点处的采样值由J、K、L、M四个整数坐标点的采样值加权平均计算得到，具体计算公式为：I_val＝J_val(1-ds)(1-dt)+K_valds(1-dt)+L_valdsdt+M_val(1-ds)dt，其中ds＝s-floor(s)，dt＝t-floor(t)，I_val，J_val，K_val，L_val，M_val分别是I、J、K、L、M点处的采样值。由于双线性插值算法的计算复杂度不高，在GPU硬件内部也内置了对双线性插值算法的支持，能达到实时的要求，并且图像质量也可以接受，因此目前所有主流的超声诊断设备上都使用其作为扫描变换模块的插值算法。

本发明使用的简单的代理几何模型，大大简化了前面几个步骤的计算负担，因此为使用更为复杂、精度更高的图像插值算法提供了可能。在步骤4中，不同于传统技术使用双线性插值算法，本发明在GPU中计算双三次图像插值。如图5所示，在I的周围取最近的16个邻域采样点A、C、E、G、W、J、M、N、V、K、L、U、B、D、F、H，其(s，t)坐标分别为：

(A_s，A_t)＝(floor(s)-1，floor(t)-1)

(C_s，C_t)＝(floor(s)-1，floor(t))

(E_s，E_t)＝(floor(s)-1，ceil(t))

(G_s，G_t)＝(floor(s)-1，ceil(t)+1)

(W_s，W_t)＝(floor(s)，floor(t)-1)

(J_s，J_t)＝(floor(s)，floor(t))

(M_s，M_t)＝(floor(s)，ceil(t))

(N_s，N_t)＝(floor(s)，ceil(t)+1)

(V_s，V_t)＝(ceil(s)，floor(t)-1)

(K_s，K_t)＝(ceil(s)，floor(t))

(L_s，L_t)＝(ceil(s)，ceil(t))

(U_s，U_t)＝(ceil(s)，ceil(t)+1)

(B_s，B_t)＝(ceil(s)+1，floor(t)-1)

(D_s，D_t)＝(ceil(s)+1，floor(t))

(F_s，F_t)＝(ceil(s)+1，ceil(t))

(H_s，H_t)＝(ceil(s)+1，ceil(t)+1)

I点处的采样值由上述16个整数坐标点的采样值加权平均计算得到，具体计算方法为：先由A、C、E、G四个采样点处的采样值A_val，C_val，E_val，G_val插值计算得到一个临时的采样值temp1，计算公式为

temp1＝A_valf(A_t-t)+C_valf(C_t-t)+E_valf(E_t-t)+G_valf(G_t-t)，

其中f(x)为一个三次多项式，其定义为：

$f\left(x\right)=\left\{\begin{array}{cc}\frac{7}{6}{\left|x\right|}^3-2{\left|x\right|}^2+\frac{8}{9},& \left|x\right|<1\\ -\frac{7}{18}{\left|x\right|}^3+2{\left|x\right|}^2-\frac{10}{3}\left|x\right|+\frac{16}{9},& 1\&le;\left|x\right|<2\\ 0,& \left|x\right|\&GreaterEqual;2\end{array}\right.$

接着分别由W、J、M、N四个采样点处的采样值W_val，J_val，M_val，N_val插值计算得到第二个临时的采样值temp2，计算公式为

temp2＝W_valf(W_t-t)+J_valf(J_t-t)+M_valf(M_t-t)+N_valf(N_t-t)；

由V、K、L、U四个采样点处的采样值V_val，K_val，L_val，U_val插值计算得到第三个临时的采样值temp3，计算公式为

temp3＝V_valf(V_t-t)+K_valf(K_t-t)+L_valf(L_t-t)+U_valf(U_t-t)；

由B、D、F、H四个采样点处的采样值B_val，D_val，F_val，H_val插值计算得到第四个临时的采样值temp4，计算公式为

temp4＝B_valf(B_t-t)+D_valf(D_t-t)+F_valf(F_t-t)+H_valf(H_t-t)；

最后由四个临时采样值temp1～temp4再次插值计算得到I点处的采样值，计算公式为I_val＝temp1·f(A_s-s)+temp2·f(W_s-s)+temp3·f(V_s-s)+temp4·f(B_s-s)。

从上面的计算过程可以看出双三次插值需要大量的运算，对于每个Fragment，光是计算三次多项式就需要20次，每个三次多项式至少需要7次乘法，加上加权平均需要的20次乘法，每个Fragment至少需要160次乘法，另外还需要16次内存访问。假设扫描变换后的图像为800×600，为了达到实时超声检查的需求，成像速度至少要高于30帧/秒，因此光是插值这一个步骤就至少需要800×600×30×160＝2.304×10⁹次乘法运算，再加上800×600×30×16＝2.304×10⁸次内存访问和一些加法运算，这在目前的CPU上是不可能的，如果用专用的FPGA或者ASIC硬件电路来实现，如此规模的电路需要很高的成本，因此目前高精度的双三次图像插值算法尚没有广泛应用。即使对于目前拥有高度并行化架构的GPU，其一次能够并行处理多个Fragment，但是其对显存的访问速度仍然远远落后于其计算速度，对于如此大的运算量和显存访问次数，仍然很难达到实时的要求。

要使得双三次插值在目前中低档的显卡上能够达到超声实时成像的需求，必须进行优化。本发明中对最耗时的两个部分结合GPU的特点进行了优化：第一个部分是三次多项式的计算，第二个部分是纹理显存的访问优化。对于三次多项式f(x)来说，其只在-2＜x＜2的范围内取非零值，在其它区间都是零，因此可以将-2＜x＜2的区间离散化成NX个点，预先算好这些离散点上的f(x)取值，然后将其作为一维纹理传输并储存在GPU的纹理显存中。GPU的纹理单元一般都是对长度为2的整数次方的纹理作了高度优化，因此在具体实施中，NX可以取值为512、1024、2048等，取值越大，三次多项式的计算精度就越高，但是所占用的纹理显存也越大，因此要根据所用显卡的实际情况来决定NX的取值。另外，GPU所采用的都是SIMD指令集，对于矢量运算提供了最优化的支持，仔细观察双三次插值中的5次插值计算，发现它们都是两个长度为4的矢量的点积计算，可以统写为：

<(val1，val2，val3，val4)·(f(x₁)，f(x₂)，f(x₃)，f(x₄)>，

其中val1～val4分别代表四个采样值，x₁～x₄代表对应的自变量取值，另外它们满足：x₄＝x₃+1＝x₂+2＝x₁+3，对于前四次插值计算，自变量只牵涉到t坐标，并且这四次所用到的x₁～x₄是对应相等的，第5次插值计算，自变量是用的s坐标，其所用x₁～x₄跟前面四次是不同的，因此只需两套x₁～x₄即可。对于目前支持OpenGL 2.0或者DirectX 9.0的比较新的GPU来讲，其纹理可以使用RGBA四个通道，每个通道32位浮点，一个纹理元素总共128位，因此可以同时放四个采样值或者f(x)的值。在步骤1生成线数据纹理时，第i条线上第j个采样点的R通道放(i，j)处的采样值，G通道放(i+1，j)处的采样值，B通道放(i+2，j)处的采样值，A通道放(i+3，j)处的采样值；同样，在生成三次多项式纹理时，对于任意一个x，R通道存储f(x)的值，G通道存储f(x+1)的值，B通道存储f(x+2)的值，A通道存储f(x+3)的值。

在具体实现优化的双三次插值算法时，要首先通过四次纹理显存访问，得到16个采样点A、C、E、G、W、J、M、N、V、K、L、U、B、D、F、H的采样值，这个可以通过使用纹理坐标

(floor(s)-1，floor(t)-1)，

(floor(s)，floor(t)-1)，

(ceil(s)，foor(t)-1)，

(ceil(s)+1，floor(t)-1)

访问线数据纹理，得到4个长度为4的RGBA矢量v1，v2，v3，v4，其中

v1＝(A_val，C_val，E_val，G_val)，v2＝(W_val，J_val，M_val，N_val)，

v3＝(V_val，K_val，L_val，U_val)，v4＝(B_val，D_val，F_val，H_val)；

再通过使用纹理坐标floor(t)-t-1和floor(s)-s-1对三次多项式纹理进行两次访问，得到2个长度为4的RGBA矢量f1，f2，其中

f1＝(f(A_t-t)，f(C_t-t)，f(E_t-t)，f(G_t-t))＝(f(W_t-t)，f(J_t-t)，f(M_t-t)，f(N_t-t))＝(f(V_t-t)，f(K_t-t)，f(L_t-t)，f(U_t-t))＝(f(B_t-t)，f(D_t-t)，f(F_t-t)，f(H_t-t))，f2＝(f(A_s-s)，f(W_s-s)，f(V_s-s)，f(B_s-s))。然后利用GPU内部提供的SIMD矢量点积指令，计算temp＝(<v1·f1>，<v2·f1>，<v3·f1>，<v4·f1>)，最后再计算出待求的采样值I_val＝<temp·f2>，并将I点处的采样值I_val输出至显示器。

通过上述的优化算法，一方面大大减少了访问纹理显存的次数，从原来的需要16次减少到现在的16/4+2＝6次，其中包括4次访问线数据纹理和2次访问三次多项式纹理；另外一方面可以充分利用GPU提供的SIMD矢量运算，将原来的160次乘法减少到5次矢量点积运算。通过这两个方面的优化，可以充分挖掘GPU的潜力，实现实时的双三次插值运算。

除了上述的优化之外，还可以通过减少要处理的Fragment的数量来进一步加快算法的速度。如图3所示，在绘制代理几何模型时，图中的阴影区域实际上并不在实际成像区域中，因此需要忽略掉阴影区域光栅化得到的那些Fragment。对于每个Fragment，在步骤3通过坐标变换计算得到它的(s，t)坐标，如果不满足0≤t≤(LN-1)，0≤s≤(PN-1)这两个条件，就可以说明此Fragment处于阴影区域中，不需要进行步骤4的插值运算，可以提前终止此Fragment的处理。

Claims (3)

1、一种高精度实时超声图像扫描变换方法，其特征在于，步骤如下：

步骤1：获取超声探头所采集的每帧图像，由LN条线数据组成，形成一个扇形，内径为r，外径为r+l，第1条线数据和第LN条线数据所夹的圆心角为2θ，LN条线数据为等角度采样，每条采样线上有PN个采样点，为等距离采样。将所述每帧图像作为二维纹理传输至GPU显存，并保存为线数据纹理；

步骤2：生成代理几何模型并绘制，包括：

使用六边形ABQRHG作为代理几何模型，其顶点坐标分别为：

(A_x，A_y)＝(-rsinθ，rcosθ)

(B_x，B_y)＝(-(r+l)sinθ，(r+l)cosθ)

(Q_x，Q_y)＝(-(r+l)sinθ，r+l)

(R_x，R_y)＝((r+l)sinθ，r+l)

(H_x，H_y)＝((r+l)sinθ，(r+l)cosθ)

(G_x，G_y)＝(r sinθ，rcosθ)

其中A和G、B和H、Q和R的x坐标相反，y坐标相同，因此只需计算A、B、Q三个点的坐标，其它三个点的坐标根据上述关系求得；

将六边形ABQRHG使用三角带形式分解为四个三角形，顶点顺序为G，H，A，R，B，Q，三角形分别为：GHA、HAR、ARB和RBQ；

使用OpenGL或者DirectX的API函数绘制三角带；

步骤3：在GPU中对代理模型绘制后的每个光栅点(Fragment)计算坐标变换，包括：

将代理模型绘制后的每个光栅点从显示器的直角坐标(x，y)转换为所述线数据纹理所在的极坐标(s，t)，计算公式为： $s=(\sqrt{x^2+y^2}-r)(\mathrm{PN}-1)/l,$ t＝(arctan(x/y)+θ)(LN-1)/2θ，其中t代表在第t条线数据上，s代表在线数据的第s个采样点上；

判断(s，t)是否满足0≤t≤(LN-1)且0≤s≤(PN-1)，如果不满足则终止计算；

步骤4：在GPU中计算双三次图像插值，包括：

在-2＜u＜2的区间上取NX个离散的三次多项式f(u)的非零值，NX为2的正整数次方，用一维纹理传输至GPU显存，并保存为三次多项式纹理，其中f(u)的定义如下：

$f\left(u\right)=\left\{\begin{array}{cc}\frac{7}{6}{\left|u\right|}^3-2{\left|u\right|}^2+\frac{8}{9},& \left|u\right|<1\\ -\frac{7}{18}{\left|u\right|}^3+2{\left|u\right|}^2-\frac{10}{3}\left|u\right|+\frac{16}{9},& 1\&le;\left|u\right|<2;\\ 0,& \left|u\right|\&GreaterEqual;2\end{array}\right.$

分别使用纹理坐标(floor(s)-1，floor(t)-1)，(floor(s)，floor(t)-1)，(ceil(s)，floor(t)-1)，(ceil(s)+1，floor(t)-1)访问所述线数据纹理，得到4个长度为4的矢量v1，v2，v3，v4，其中，函数floor(v)代表不大于v的最大的整数，函数ceil(v)代表不小于v的最小的整数；

分别使用纹理坐标floor(t)-t-1和floor(s)-s-1访问所述三次多项式纹理，得到2个长度为4的矢量f1，f2；

利用GPU的SIMD矢量计算指令，计算得到一个长度为4的矢量temp＝(<v1·f1>，<v2·f1>，<v3·f1>，<v4·f1>)，<·>表示矢量点积计算；

利用GPU的SIMD矢量计算指令，计算待求的采样值I_val＝<temp·f2>；

最后，将采样值I_val输出至显示器的(x，y)坐标处显示。

2、根据权利要求1所述的高精度实时超声图像扫描变换方法，其特征在于，所述线数据纹理采用RGBA四个通道，每个通道为32位浮点数的格式；对于第i条线上第j个采样点，其R通道存放(i，j)处的采样值，G通道存放(i+1，j)处的采样值，B通道存放(i+2，j)处的采样值，A通道存放(i+3，j)处的采样值。

3、根据权利要求1所述的高精度实时超声图像扫描变换方法，其特征在于，所述三次多项式纹理采用RGBA四个通道，每个通道为32位浮点数的格式；对于在区间-2＜u＜2的u，R通道存储f(u)的值，G通道存储f(u+1)的值，B通道存储f(u+2)的值，A通道存储f(u+3)的值。

Images