技术特征：
1.一种基于cuda的波数域三维超声全矩阵成像方法，其特征在于，包括如下步骤：步骤一：设置面阵超声换能器，以直接接触的方式与各向同性均匀介质的试件耦合，使用全矩阵模式采集原始三维全矩阵数据，确定所用试件的声速、面阵超声换能器的参数以及所成三维图像的大小范围，以面阵超声换能器左下角的阵元的中心为坐标原点，以面阵超声换能器的两边为x轴和y轴，以垂直换能器平面的声波传播方向为z轴，建立原始数据坐标系；步骤二：确定五维傅里叶变换每个维度的计算长度，基于weyl’s identity编写cuda核函数，计算原始数据坐标系和图像坐标系之间的波数转换关系；步骤三：拷贝cpu内存中的原始三维全矩阵数据至gpu显存中，于cuda核函数中将三维全矩阵数据重排列为五维全矩阵数据，数据索引按维度组成；使用行-列算法分解五维傅里叶变换，按照矩阵转置的坐标关系改变数据在存储空间中的索引，将五维傅里叶变换转换为计算多次二维和三维傅里叶变换；使用快速五维傅里叶变换，将原始数据转换到波数域，并且对频率维度做矩形窗截断；步骤四：对于波数域五维矩阵中发射维度的每组三维数据，按照步骤二中原始数据坐标系和图像坐标系之间的波数转换关系为插值索引，选取空间维度上相邻的若干个点为参照，使用插值算法分别计算波数域五维矩阵在三维图像下的stolt映射；步骤五：对于步骤四中计算得出的波数域stolt映射，提取其中每组三维数据分别对齐叠加，获得三维波数域图像矩阵，做三维快速反傅里叶变换，将图像矩阵由波数域重新变换到空间域，生成三维超声全矩阵图像，截断感兴趣区域以外的部分；步骤六：对于步骤五所生成三维超声全矩阵图像，将其从gpu显存中拷贝到cpu内存中，压入vtk模块的数据处理流水线，使用最大密度投影渲染为二维图像并显示在窗口，并且能够实时更新数据，交互改变视角位置和焦点远近。2.根据权利要求1所述的基于cuda的波数域三维超声全矩阵成像方法，其特征在于，步骤一中所述的试件为均匀各向同性介质，其声速设为c，使用的面阵超声换能器的阵元数分布为其中离散下标w1,w2分别代表面阵超声换能器在x轴和y轴方向上的阵元数量，而x轴和y轴的阵元中心间距pitch分别为试件及成像区域摆放在z轴正半轴，成像区域的像素点数在x，y和z轴方向分别设为n
x
,n
y
,n
z
，其像素点的间距分别为d
x
,d
y
,d
z
；成像区域的原点与原始数据坐标系原点重合；获得的原始三维全矩阵超声数据记为e
raw
(t0,v,u)。3.根据权利要求1所述的基于cuda的波数域三维超声全矩阵成像方法，其特征在于，所述的步骤二具体为：由深度、宽度和像素间隔求出五维傅里叶变换每个维度的所需要的点数值，分别为n
t
,其中离散下标t代表时间，v1,v2代表接收阵元在x和y轴方向的分量，u1,u2代表发射阵元在x和y轴方向的分量；使用全矩阵采集进行成像时，有通过weyl’s identity和包含格林函数的前向传播模型，所述的weyl’s identity在三维空间中的具体表达式为：
其中，g(ω,x,y,z)为三维空间中平面波分解的格林函数，k为空间波数，k
x
,k
y
分别为空间波数k在x和y轴上的分量，ω为角频率，x,y,z分别为x，y和z轴方向上的空间位移；根据成像系统简化的前向传播模型，发射-接收信号对，也即接收数据即为e(t,v1,v2,u1,u2)，推导出波数域中原始坐标系和图像坐标系的坐标轴波数转换关系：推导出波数域中原始坐标系和图像坐标系的坐标轴波数转换关系：推导出波数域中原始坐标系和图像坐标系的坐标轴波数转换关系：将上式改写为用来表示k,k
v1
,k
v2
，为了与原始数据的空间波数做区分，这里以k
g
,k
v1g
,k
v2g
代替：k
v1g
＝k
x-k
u1g
,k
v2g
＝k
y-k
v1g
,其中，k
x
，k
y
，k
z
分别为空间波数k在x、y、z轴上的分量，其中离散下标x,y,z分别代表x，y和z轴，分别为接收阵元在x和y轴方向的波数分量以及发射阵元在x和y轴方向的波数分量。4.根据权利要求1所述的基于cuda的波数域三维超声全矩阵成像方法，其特征在于，所述的步骤三具体为：三维原始数据e
raw
(t0,v,u)转换为五维数据e(t,v1,v2,u1,u2)，其中索引t0代表原始数据中的时间，索引v,u分别代表全矩阵采集的接收和发射阵元。表示为离散矩阵形式为e
raw
[t0][v][u]和e[t][v1][v2][u1][u2]，其数据索引的关系式如下：t＝t
0000
其中，函数floor代表向下取整的函数；使用行-列算法分解五维傅里叶变换，转换为分别执行多次三维和二维傅里叶变换，将傅里叶变换的结果及中间步骤记为e
n
(t,v1,v2,u1,u2)其计算流程如下：(1)次二维点快速傅里叶变换：
(2)转置数据矩阵(2)转置数据矩阵(3)次三维n
t
,点快速傅里叶变换：截断数据，保留ω∈[2πf
min
,2πf
max
]的部分，设其点数为5.根据权利要求1所述的基于cuda的波数域三维超声全矩阵成像方法，其特征在于，所述的步骤四中的插值算法为八线性插值。6.根据权利要求1所述的基于cuda的波数域三维超声全矩阵成像方法，其特征在于，所述的步骤五中波数域中沿所有发射阵元叠加，并做三维快速反傅里叶变换，其表达式如下所示：式如下所示：其中，f(z,x,y)即为所计算的三维超声图像，截断感兴趣区域之外的部分，得到感兴趣区域之内的图像i(z,x,y)。7.根据权利要求1所述的基于cuda的波数域三维超声全矩阵成像方法，其特征在于，步骤六中所述的压入vtk模块的数据处理流水线，使用最大密度投影渲染为二维图像，具体包括如下步骤：a：初始化vtk模块参数；b：三维图像数据做滤波操作；c：使用最大密度投影将三维图像数据渲染为二维图像；d：设置vtk交互窗口，在窗口中实时显示和刷新二维图像。
技术总结
本发明公开了一种基于CUDA的波数域三维超声全矩阵成像方法。本发明将原始全矩阵数据表示为五维矩阵的形式，通过一次五维快速傅里叶变换将得到的原始全矩阵数据转换到波数域，借助于Weyl’s Identity实现原始数据坐标系和成像坐标系之间的波数转换，通过八线性插值得到图像坐标系的映射，再通过一次三维快速反傅里叶变换得到三维图像。本发明采用VTK模块，使用MIP(最大密度投影)的光线投影方法得到渲染后的二维图像，并且可以实时改变三维图像的观察角度和距离远近。本发明具有对可变感兴趣区域的实时高效率、高质量三维成像和二维渲染能力，能够灵活移植到具备CUDA运算平台的计算机设备，可以应对大量三维全矩阵数据的快速处理。理。理。
技术研发人员：陈沐 金浩然 杨克己
受保护的技术使用者：浙江大学
技术研发日：2023.02.02
技术公布日：2023/3/16