CN115808469A - 基于cuda的波数域三维超声全矩阵成像方法 - Google Patents
基于cuda的波数域三维超声全矩阵成像方法 Download PDFInfo
- Publication number
- CN115808469A CN115808469A CN202310052077.2A CN202310052077A CN115808469A CN 115808469 A CN115808469 A CN 115808469A CN 202310052077 A CN202310052077 A CN 202310052077A CN 115808469 A CN115808469 A CN 115808469A
- Authority
- CN
- China
- Prior art keywords
- dimensional
- data
- wave number
- matrix
- image
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 239000011159 matrix material Substances 0.000 title claims abstract description 77
- 238000003384 imaging method Methods 0.000 title claims abstract description 62
- HPTJABJPZMULFH-UHFFFAOYSA-N 12-[(Cyclohexylcarbamoyl)amino]dodecanoic acid Chemical compound OC(=O)CCCCCCCCCCCNC(=O)NC1CCCCC1 HPTJABJPZMULFH-UHFFFAOYSA-N 0.000 title claims abstract 12
- 238000000034 method Methods 0.000 claims abstract description 22
- 238000009877 rendering Methods 0.000 claims abstract description 18
- 238000012545 processing Methods 0.000 claims abstract description 17
- 238000013507 mapping Methods 0.000 claims abstract description 10
- 238000006243 chemical reaction Methods 0.000 claims abstract description 9
- 238000004364 calculation method Methods 0.000 claims description 25
- 238000012360 testing method Methods 0.000 claims description 23
- 230000006870 function Effects 0.000 claims description 14
- 230000008569 process Effects 0.000 claims description 7
- 238000001914 filtration Methods 0.000 claims description 5
- 238000002604 ultrasonography Methods 0.000 claims description 5
- 230000002452 interceptive effect Effects 0.000 claims description 4
- 238000004613 tight binding model Methods 0.000 claims description 4
- 230000008878 coupling Effects 0.000 claims description 3
- 238000010168 coupling process Methods 0.000 claims description 3
- 238000005859 coupling reaction Methods 0.000 claims description 3
- 238000003860 storage Methods 0.000 claims description 3
- 238000005520 cutting process Methods 0.000 claims description 2
- 238000000354 decomposition reaction Methods 0.000 claims description 2
- 238000006073 displacement reaction Methods 0.000 claims description 2
- 238000009826 distribution Methods 0.000 claims description 2
- 238000003825 pressing Methods 0.000 claims description 2
- 230000001131 transforming effect Effects 0.000 claims description 2
- 230000017105 transposition Effects 0.000 claims description 2
- 230000008859 change Effects 0.000 abstract description 3
- 230000007547 defect Effects 0.000 description 18
- 238000001514 detection method Methods 0.000 description 13
- 239000000523 sample Substances 0.000 description 8
- 238000005516 engineering process Methods 0.000 description 7
- 238000009659 non-destructive testing Methods 0.000 description 4
- 238000005070 sampling Methods 0.000 description 4
- 238000011156 evaluation Methods 0.000 description 3
- 238000001228 spectrum Methods 0.000 description 3
- 230000000007 visual effect Effects 0.000 description 3
- 230000009286 beneficial effect Effects 0.000 description 2
- 230000005540 biological transmission Effects 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 238000007689 inspection Methods 0.000 description 2
- 230000003993 interaction Effects 0.000 description 2
- 229910052751 metal Inorganic materials 0.000 description 2
- 239000002184 metal Substances 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000012805 post-processing Methods 0.000 description 2
- 229910052782 aluminium Inorganic materials 0.000 description 1
- XAGFODPZIPBFFR-UHFFFAOYSA-N aluminium Chemical compound [Al] XAGFODPZIPBFFR-UHFFFAOYSA-N 0.000 description 1
- 230000015572 biosynthetic process Effects 0.000 description 1
- 238000007796 conventional method Methods 0.000 description 1
- 238000012937 correction Methods 0.000 description 1
- 238000013506 data mapping Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000009434 installation Methods 0.000 description 1
- 239000003208 petroleum Substances 0.000 description 1
- 238000004451 qualitative analysis Methods 0.000 description 1
- 238000011158 quantitative evaluation Methods 0.000 description 1
- 230000035945 sensitivity Effects 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 238000003786 synthesis reaction Methods 0.000 description 1
- 238000012285 ultrasound imaging Methods 0.000 description 1
Images
Landscapes
- Ultra Sonic Daignosis Equipment (AREA)
Abstract
本发明公开了一种基于CUDA的波数域三维超声全矩阵成像方法。本发明将原始全矩阵数据表示为五维矩阵的形式,通过一次五维快速傅里叶变换将得到的原始全矩阵数据转换到波数域,借助于Weyl’s Identity实现原始数据坐标系和成像坐标系之间的波数转换,通过八线性插值得到图像坐标系的映射,再通过一次三维快速反傅里叶变换得到三维图像。本发明采用VTK模块,使用MIP(最大密度投影)的光线投影方法得到渲染后的二维图像,并且可以实时改变三维图像的观察角度和距离远近。本发明具有对可变感兴趣区域的实时高效率、高质量三维成像和二维渲染能力,能够灵活移植到具备CUDA运算平台的计算机设备,可以应对大量三维全矩阵数据的快速处理。
Description
技术领域
本发明属于工业超声无损检测技术领域,涉及一种基于CUDA的波数域三维超声全矩阵成像的实现方法。
背景技术
作为重要的无损检测技术,超声检测从一开始的无损探伤、无损检测的定性分析,已经转向定量无损评价的方向发展,这也对超声检测技术的定量检测特性提出了更高的要求。随着现代工业结构的复杂化以及工业设备运行条件的极端化,对工件与设备质量的要求也越来越高,同时也要求超声检测技术能够发现和检测出更加细微的缺陷,并且对检出缺陷的位置、形态和类型都有了更高的要求。在这种情况下,传统的常规超声已经无法满足要求。相控阵换能器是由多组小型的晶片阵元按照一定的几何方式排列而成,通常有线阵、面阵以及环形阵列等结构。相控阵超声检测技术凭借其波束的柔性合成与自适应控制能力,在检测效率、检测范围、声束可达性和灵敏度方面具有明显的技术优势,也使其在近二十年来得到了快速发展和广泛应用。
在相控阵超声无损检测中,全矩阵数据采集通过相控阵换能器依次激励单个换能器阵元发射超声波,然后用所有阵元接收回波信号并存储,充分利用了相控阵超声换能器中全部阵元发射接收的信号数据,包含了换能器阵元从不同角度入射到被测工件内的声波信息,拥有更加全面的被测工件内部信息,对于缺陷位置、尺寸和形态的描述更全面、更准确。尤其是当缺陷位置不可知时,相控阵超声全矩阵成像无需反复调整声束偏转角度与聚焦位置,只需按照原有的数据采集方式即可获得未知缺陷的全面信息。这一特性大大提高了相控阵超声成像的检测精度。并且,全矩阵成像技术不需要预先计算延时聚焦法则,也不需要根据不同的缺陷类型和缺陷位置调整波束合成方案,只需按照全矩阵采集方式对所有被测工件进行检测与信号收集,然后通过上位机软件进行全区域聚焦成像,即可一次性完成一定区域内所有缺陷的检测。
相比于传统的二维超声成像检测,三维成像以其对于空间缺陷更为全面准确的表征,更为广泛的成像区域获得越来越多的关注。二维图像往往只能获取试件的单个截面信息,不能获得试件缺陷整体的评估,而在石油管道、流水线检测和螺纹紧固件等检测方面,三维成像逐渐取代二维,成为不可或缺的无损检测手段。三维图像可以有效定位缺陷在试件中的空间位置,能够渲染出缺陷的形状和特征,便于定量无损评估试件的破坏部位和损坏程度。
全矩阵三维成像具备诸多便利,但是其更大的数据量和计算消耗给后处理能力提出了更高的要求,严重影响了成像的效率和精度。相较于传统的时域算法,波数域算法借助于快速傅里叶变换较低的计算复杂度,大幅提高全矩阵三维成像的计算效率,而采用较为精确的八线性插值,可以有效弥补插值所带来的相位误差,提高感兴趣区域内缺陷成像的精度。为了在保证成像精度的基础上进一步提高计算效率,满足实时在线检测的要求,充分利用快速傅里叶变换并行计算和灵活的移植能力,波数域算法可以进一步利用具备CUDA平台的GPU来加速运算过程,并且借助于VTK模块实现GPU的二维渲染,达成原始数据到图像显示与键鼠交互的完整流程。
发明内容
本发明的目的在于提供一种基于CUDA的波数域三维超声全矩阵成像方法,采用本发明的方法可以提高三维全矩阵成像的效率和精度,渲染二维最大密度投影,实现实时在线三维超声无损检测的目标。根据所成三维图像的大小和范围设置五维快速傅里叶变换的基本参数,根据面阵超声换能器的中心频率和带宽设置波数域的截断滤波,基于Weyl’sIdentity获得原始数据坐标系和图像坐标系之间的波数转换关系,通过八线性插值将波数域原始数据映射到波数域图像,采用三维快速反傅里叶变换得到所需的三维图像。通过VTK模块渲染三维图像的最大密度投影,并通过交互窗口显示,三维图像的观察视角和距离可以实时改变,数据可以实时刷新,有利于实时三维超声无损定量评价。
本发明采用的技术方案是:
本发明首先提供了一种基于CUDA的波数域三维超声全矩阵成像方法,其包括如下步骤:
步骤一:设置面阵超声换能器,以直接接触的方式与各向同性均匀介质的试件耦合,使用全矩阵模式采集原始三维全矩阵数据,确定所用试件的声速、面阵超声换能器的参数以及所成三维图像的大小范围,以面阵超声换能器左下角的阵元的中心为坐标原点,以面阵超声换能器的两边为X轴和Y轴,以垂直换能器平面的声波传播方向为Z轴,建立原始数据坐标系;
步骤二:确定五维傅里叶变换每个维度的计算长度,基于Weyl’s Identity编写CUDA核函数,计算原始数据坐标系和图像坐标系之间的波数转换关系;
步骤三:拷贝CPU内存中的原始三维全矩阵数据至GPU显存中,于CUDA核函数中将三维全矩阵数据重排列为五维全矩阵数据,数据索引按维度组成;使用行-列算法分解五维傅里叶变换,按照矩阵转置的坐标关系改变数据在存储空间中的索引,将五维傅里叶变换转换为计算多次二维和三维傅里叶变换;使用快速五维傅里叶变换,将原始数据转换到波数域,并且对频率维度做矩形窗截断;
步骤四:对于波数域五维矩阵中发射维度的每组三维数据,按照步骤二中原始数据坐标系和图像坐标系之间的波数转换关系为插值索引,选取空间维度上相邻的若干个点为参照,使用插值算法分别计算波数域五维矩阵在三维图像下的Stolt映射;
步骤五:对于步骤四中计算得出的波数域Stolt映射,提取其中每组三维数据分别对齐叠加,获得三维波数域图像矩阵,做三维快速反傅里叶变换,将图像矩阵由波数域重新变换到空间域,生成三维超声全矩阵图像,截断感兴趣区域以外的部分;
步骤六:对于步骤五所生成三维超声全矩阵图像,将其从GPU显存中拷贝到CPU内存中,压入VTK模块的数据处理流水线,使用最大密度投影渲染为二维图像并显示在窗口,并且能够实时更新数据,交互改变视角位置和焦点远近。
作为本发明的优选实施方案,步骤一中所述的试件为均匀各向同性介质,其声速设为c,使用的面阵超声换能器的阵元数分布为其中离散下标w1,w2分别代表面阵超声换能器在X轴和Y轴方向上的阵元数量,而X轴和Y轴的阵元中心间距pitch分别为试件及成像区域摆放在Z轴正半轴,成像区域的像素点数在X,Y和Z轴方向分别设为Nx,Ny,Nz,其像素点的间距分别为dx,dy,dz;成像区域的原点与原始数据坐标系原点重合;获得的原始三维全矩阵超声数据记为eraw(t0,v,u)。
作为本发明的优选实施方案,所述的步骤二具体为:
由深度、宽度和像素间隔求出五维傅里叶变换每个维度的所需要的点数值,分别为Nt,其中离散下标t代表时间,v1,v2代表接收阵元在X和Y轴方向的分量,u1,u2代表发射阵元在X和Y轴方向的分量;使用全矩阵采集进行成像时,有
通过Weyl’s Identity和包含格林函数的前向传播模型,所述的Weyl’s Identity在三维空间中的具体表达式为:
其中,G(ω,x,y,z)为三维空间中平面波分解的格林函数,k为空间波数,kx,ky分别为空间波数k在X和Y轴上的分量,ω为角频率,x,y,z分别为X,Y和Z轴方向上的空间位移;根据成像系统简化的前向传播模型,发射-接收信号对,也即接收数据即为e(t,v1,v2,u1,u2),
其关于时间的频谱E(ω,v1,v2,u1,u2)可以用如下关系式计算:
其中,函数G为格林函数,f(x,y,z)为前向传播模型中的缺陷散射点分布,也即所需要获得的三维图像。将Weyl’s Identity代入前向传播模型公式,以v1,v2,u1,u2为变量改写关系式:
其中,kv1,分别为接收阵元在X和Y轴方向的波数分量以及发射阵元在X和Y轴方向的波数分量,F(kx,ky,kz)为三维图像f(x,y,z)的波数域形式。因此,e(t,v1,v2,u1,u2)的五维傅里叶变换形式可以表示为如下关系式:
其中又有如下替换关系式:
kv1g=kx-ku1g,
kv2g=ky-kv1g,
对于给定的面阵探头参数和初始化成像区域,其空间波数和各个方向的波数分量已知,如下关系式表示:
fs为接收全矩阵信号的采样频率。根据步骤三做傅里叶变换之后根据探头带宽做矩形截断,其截断最小最大频率分别设为fmin,fmax,相对应的角频率范围ω∈[2πfmin,2πfmax]。因此截断之后的角频率为:
根据步骤三和步骤四,为了充分利用插值映射的原始数据,三维图像在Z轴方向的空间波数kz设为如下关系式:
作为本发明的优选实施方案,所述的步骤三具体为:三维原始数据eraw(t0,v,u)转换为五维数据e(t,v1,v2,u1,u2),其中索引t0代表原始数据中的时间,索引v,u分别代表全矩阵采集的接收和发射阵元。表示为离散矩阵形式为eraw[t0][v][u]和e[t][v1][v2][u1][u2],其数据索引的关系式如下:
t=t0
其中,函数floor代表向下取整的函数;使用行-列算法分解五维傅里叶变换,转换为分别执行多次三维和二维傅里叶变换,将傅里叶变换的结果及中间步骤记为En(t,v1,v2,u1,u2)其计算流程如下:
由于缺陷成像检测更关注相对幅值,忽略正比系数。代入步骤二中预先计算的波数参量,由于空间波数在原始数据和图像数据的转换非线性,并且锁定ku1,ku2之后仅剩ω,kv1,kv2三个参量,本发明采用八线性插值为插值算法,通过相应CUDA核函数完成计算。
作为本发明的优选实施方案,所述的步骤五中波数域中沿所有发射阵元叠加,并做三维快速反傅里叶变换,其表达式如下所示:
其中,f(z,x,y)即为所计算的三维超声图像,截断感兴趣区域之外的部分,得到感兴趣区域之内的图像I(z,x,y)。
作为本发明的优选实施方案,步骤六中所述的压入VTK模块的数据处理流水线,使用最大密度投影渲染为二维图像,具体包括如下步骤:
a:初始化VTK模块参数;
b:三维图像数据做滤波操作;
c:使用最大密度投影将三维图像数据渲染为二维图像;
d:设置VTK交互窗口,在窗口中实时显示和刷新二维图像。
本发明的有益效果主要表现在:
(1)本发明利用快速傅里叶变换将常规的时域延时叠加转换为波数域处理,大幅降低以往时域算法的计算复杂度,从而提高三维超声成像的效率和精度,实现实时高帧率三维超声全矩阵成像。
(2)本发明采用了高效简洁的波数域计算公式,能够整体计算图像中所有像素点的值,公式稳定性高,方法简便易于理解,避免常规时域算法分别延时叠加计算每个像素点的低效率,从而从原理上提高三维超声全矩阵成像的速度。
(3)本发明采用了波数域的频谱处理技术,降低了滤波等提升图像质量的手段的计算时间,从而有效降低计算时间,以较小的计算代价获取较高的成像效果。
(4)本发明采用了自主编写和修改的插值方法,解决难以在软件中实现复杂插值方法的问题,从而有效解决三维波数域超声成像的软件实现问题,并且其计算速度较快,同样降低了计算消耗。
(5)本发明利用CUDA的GPU并行处理能力进一步提高三维成像的效率,由于所提出的利用快速傅里叶变换转换为波数域处理、频谱处理技术和插值方法均在本发明中以CUDA平台的方法实现,主要计算在GPU中进行,解决常规CPU方法中串行计算的低效率,大幅节省成像时间。
(6)本发明的算法具备极强的便捷性和可移植能力,采用具备CUDA平台的计算设备,较之普通FPGA(现场可编程门阵列)避免烧写代码和代码更新升级周期相对较长的问题,便于纠错和个人定制,使得能够在复杂灵活的现场环境实现三维实时超声成像。
(7)本发明的算法不仅仅提供三维超声图像,更能够在获得三维超声图像之后,进一步渲染出二维映射图像以便显示在屏幕,同时具备实时更新图像源和相机视角的能力,渲染图窗可以实时交互,这避免了常规方法只计算三维图像无法给用户提供直观实时的成像检测能力,大大有利于用户在现场及后处理中的读图与判别能力。
附图说明
图1是本发明方法整体流程图。
图2是三维检测模型的安装位置示意图。
图3是使用VTK模块实时渲染显示二维图像流程图。
图4是使用三维实时超声成像及实时渲染二维图像结果图。
具体实施方式
以下结合附图对本发明作进一步说明。
本发明的实施例如下:
如图1所示,本发明利用一种基于CUDA的波数域三维超声全矩阵成像方法的操作流程可以分为以下几个步骤:
步骤一:设置面阵相控阵探头,以直接接触的方式与各向同性均匀介质的试件耦合,使用全矩阵模式采集原始数据。各部件及相关参数如图2所示,超声换能器为8×8的正方形面阵探头,每个阵元的中心间距在X和Y轴方向的分量分别为dx=1.5mm,dy=1.5mm,每个阵元的形状为正方形,阵元与阵元间的缝隙在X轴和Y轴方向分别为s=0.1mm。试件选取为金属铝,其声速为c=6350m/s,三维成像区域为面阵探头下方,其原点与面阵探头的原点重合,深度为hz,宽度在X和Y轴方向上分别为hx和hy,每个方向的像素间隔分别为δx,δy,δz,由此计算出每个方向的像素数Nx=ceil(hx/δx),Ny=ceil(hy/δy),Nz=ceil(hz/δz)。
步骤二:由深度、宽度和像素间隔求出五维傅里叶变换每个维度的所需要的点数值,分别为Nt,通过Weyl’s Identity和包含格林函数的前向传播模型,推导出波数域中原始坐标系和图像坐标系的坐标轴波数转换关系:
具体实施需要改变变量的组成形式,将上式改写为用kx,ky,kz来表示k,为了与原始数据的空间波数做区分,这里以kg,kv1g,kv2g代替,其中下标的g,v1g,v2g代表在成像坐标系下的空间波数,X轴方向的波数和Y轴方向的波数:
kv1g=kx-ku1g,
kv2g=ky-kv1g,
由于超声信号的带限特性,根据探头带宽做矩形截断,其截断最小最大频率分别设为fmin,fmax,相对应的角频率范围ω∈[2πfmin,2πfmax]。因此截断之后的角频率为:
步骤三:将原始全矩阵三维数据变换为五维矩阵,通过多次二维和三维傅里叶变换代替直接五维傅里叶变换,将数据转换至波数域。
所述的三维原始数据eraw(t0,v,u)转换为五维数据e(t,v1,v2,u1,u2),表示为离散矩阵形式为eraw[t0][v][u]和e[t][v1][v2][u1][u2],其数据索引的关系式如下:
t=t0
使用行-列算法分解五维傅里叶变换,转换为分别执行多次三维和二维傅里叶变换,将傅里叶变换的结果及中间步骤记为En(t,v1,v2,u1,u2)其计算流程如下:
步骤四:计算得到的波数域数据通过kg,kv1g,kv2g的转换关系,对于每一次发射所对应的三维数据矩阵,分别进行八线性插值计算,转换至对应的图像坐标系,X、Y和Z轴方向的插值索引关系式如下:
(1)所有数据点依次并行读取,首先获取插值所需的立方体在三个方向X,Y和Z轴上坐标最小的点:
indz=floor((kz-min(k))/(k(1)-k(0)))
borda=(1-indx)/Nx
bordb=(1-indy)/Ny
a=kv1g-indx
b=kv2g-indy
c=kg-indz
步骤五:用步骤四得到的每组发射对应的插值后的图像坐标系下的波数域数据,对齐叠加,和为一组三维数据矩阵,该矩阵做三维快速反傅里叶变换即得到三维超声全矩阵图像,根据每个方向的感兴趣区域大小Nx_roi,Ny_roi,Nz_roi截断图像。
步骤五中矩阵做三维快速反傅里叶变换,其表达式如下所示:
其中,f(z,x,y)即为所计算的三维超声图像,截断感兴趣区域之外的部分,得到I(z,x,y)。
步骤六:将数据由GPU显存中拷贝至内存中,使用VTK模块的最大密度投影渲染二维图像,在数据映射过程中缩小网格间距,提升图像质量,设置相机视角和图像窗口,增加窗口交互功能。渲染过程中为了提高渲染效果,使用10倍升采样因子。本实施例使用VTK模块实时渲染显示二维图像的流程如图3所示,具体包括:
a:初始化模块参数,设置VTK数据原点和网格间距,设置VTK体积映射的内存空间和采样间距,设置VTK体积的属性,包括颜色和不透明度等。
b:三维图像数据中间处理,进行升降采样和滤波等操作,得到处理后的三维图像。
c:处理后的三维图像数据使用最大密度投影渲染为二维图像数据,需要将体积映射和属性添加至VTK体绘制变量。并且设置VTK渲染变量及相机视角,将三维图像数据导入渲染变量,获得二维图像数据。
d:设置VTK交互窗口变量,导入二维图像数据,根据原始全矩阵数据的变化,实时显示和刷新窗口图像,相机视角同样可以实时改变。
设计一个斜向排列具有5个孔洞型缺陷的方形金属试件,5个孔洞大小一致,在倾斜方向上等距排列。以Nt=1024,为例,选择三维成像区域的像素数分别为Nx=64,Ny=64,Nz=512,五维傅里叶变换之后频率维度截断最小最大频率设为fmin=3MHz,fmax=7MHz。选取缺陷所处的大致范围为感兴趣区域,其大小为Nx_roi=64,Ny_roi=64,Nz_roi=300,通过波数域的计算处理,使用VTK渲染,成像结果如图4所示,可以完整有层次地利用像素点的强度在图像中展现设计的5个孔洞的缺陷,并且试件的方形底面也能完整的计算出来,图像质量较高,没有条纹感。
由于快速傅里叶变换的低算法复杂度和GPU的并行计算能力,可以实现大范围实时高帧率三维超声全矩阵成像、渲染和显示。
以上所述实施例仅表达了本发明的几种实施方式,其描述较为具体和详细,但并不能因此而理解为对本发明专利范围的限制。应当指出的是,对于本领域的普通技术人员来说,在不脱离本发明构思的前提下,还可以做出若干变形和改进,这些都属于本发明的保护范围。因此,本发明的保护范围应以所附权利要求为准。
Claims (7)
1.一种基于CUDA的波数域三维超声全矩阵成像方法,其特征在于,包括如下步骤:
步骤一:设置面阵超声换能器,以直接接触的方式与各向同性均匀介质的试件耦合,使用全矩阵模式采集原始三维全矩阵数据,确定所用试件的声速、面阵超声换能器的参数以及所成三维图像的大小范围,以面阵超声换能器左下角的阵元的中心为坐标原点,以面阵超声换能器的两边为X轴和Y轴,以垂直换能器平面的声波传播方向为Z轴,建立原始数据坐标系;
步骤二:确定五维傅里叶变换每个维度的计算长度,基于Weyl’s Identity编写CUDA核函数,计算原始数据坐标系和图像坐标系之间的波数转换关系;
步骤三:拷贝CPU内存中的原始三维全矩阵数据至GPU显存中,于CUDA核函数中将三维全矩阵数据重排列为五维全矩阵数据,数据索引按维度组成;使用行-列算法分解五维傅里叶变换,按照矩阵转置的坐标关系改变数据在存储空间中的索引,将五维傅里叶变换转换为计算多次二维和三维傅里叶变换;使用快速五维傅里叶变换,将原始数据转换到波数域,并且对频率维度做矩形窗截断;
步骤四:对于波数域五维矩阵中发射维度的每组三维数据,按照步骤二中原始数据坐标系和图像坐标系之间的波数转换关系为插值索引,选取空间维度上相邻的若干个点为参照,使用插值算法分别计算波数域五维矩阵在三维图像下的Stolt映射;
步骤五:对于步骤四中计算得出的波数域Stolt映射,提取其中每组三维数据分别对齐叠加,获得三维波数域图像矩阵,做三维快速反傅里叶变换,将图像矩阵由波数域重新变换到空间域,生成三维超声全矩阵图像,截断感兴趣区域以外的部分;
步骤六:对于步骤五所生成三维超声全矩阵图像,将其从GPU显存中拷贝到CPU内存中,压入VTK模块的数据处理流水线,使用最大密度投影渲染为二维图像并显示在窗口,并且能够实时更新数据,交互改变视角位置和焦点远近。
3.根据权利要求1所述的基于CUDA的波数域三维超声全矩阵成像方法,其特征在于,所述的步骤二具体为:
由深度、宽度和像素间隔求出五维傅里叶变换每个维度的所需要的点数值,分别为Nt,其中离散下标t代表时间,v1,v2代表接收阵元在X和Y轴方向的分量,u1,u2代表发射阵元在X和Y轴方向的分量;使用全矩阵采集进行成像时,有
通过Weyl’s Identity和包含格林函数的前向传播模型,所述的Weyl’s Identity在三维空间中的具体表达式为:
其中,G(ω,x,y,z)为三维空间中平面波分解的格林函数,k为空间波数,kx,ky分别为空间波数k在X和Y轴上的分量,ω为角频率,x,y,z分别为X,Y和Z轴方向上的空间位移;根据成像系统简化的前向传播模型,发射-接收信号对,也即接收数据即为e(t,v1,v2,u1,u2),
推导出波数域中原始坐标系和图像坐标系的坐标轴波数转换关系:
kv1g=kx-ku1g,
kv2g=ky-kv1g,
4.根据权利要求1所述的基于CUDA的波数域三维超声全矩阵成像方法,其特征在于,所述的步骤三具体为:三维原始数据eraw(t0,v,u)转换为五维数据e(t,v1,v2,u1,u2),其中索引t0代表原始数据中的时间,索引v,u分别代表全矩阵采集的接收和发射阵元。表示为离散矩阵形式为eraw[t0][v][u]和e[t][v1][v2][u1][u2],其数据索引的关系式如下:
t=t0
其中,函数floor代表向下取整的函数;使用行-列算法分解五维傅里叶变换,转换为分别执行多次三维和二维傅里叶变换,将傅里叶变换的结果及中间步骤记为En(t,v1,v2,u1,u2)其计算流程如下:
5.根据权利要求1所述的基于CUDA的波数域三维超声全矩阵成像方法,其特征在于,所述的步骤四中的插值算法为八线性插值。
7.根据权利要求1所述的基于CUDA的波数域三维超声全矩阵成像方法,其特征在于,步骤六中所述的压入VTK模块的数据处理流水线,使用最大密度投影渲染为二维图像,具体包括如下步骤:
a:初始化VTK模块参数;
b:三维图像数据做滤波操作;
c:使用最大密度投影将三维图像数据渲染为二维图像;
d:设置VTK交互窗口,在窗口中实时显示和刷新二维图像。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202310052077.2A CN115808469A (zh) | 2023-02-02 | 2023-02-02 | 基于cuda的波数域三维超声全矩阵成像方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202310052077.2A CN115808469A (zh) | 2023-02-02 | 2023-02-02 | 基于cuda的波数域三维超声全矩阵成像方法 |
Publications (1)
Publication Number | Publication Date |
---|---|
CN115808469A true CN115808469A (zh) | 2023-03-17 |
Family
ID=85487349
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202310052077.2A Pending CN115808469A (zh) | 2023-02-02 | 2023-02-02 | 基于cuda的波数域三维超声全矩阵成像方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN115808469A (zh) |
Cited By (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN117710184A (zh) * | 2024-01-04 | 2024-03-15 | 浙江大学 | 基于fpga和gpu的超声全聚焦波数域成像方法、系统及装置 |
CN117825512A (zh) * | 2024-01-04 | 2024-04-05 | 浙江大学 | 面向钢轨轨头缺陷的高速超声动态成像检测方法及装置 |
Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN112067698A (zh) * | 2020-09-14 | 2020-12-11 | 南昌航空大学 | 一种时频结合快速全聚焦超声成像方法 |
CN112684005A (zh) * | 2020-12-10 | 2021-04-20 | 苏州热工研究院有限公司 | 基于二维矩阵换能器的全聚焦检测方法 |
CN114487117A (zh) * | 2022-02-18 | 2022-05-13 | 浙江大学 | 一种超声相控阵全矩阵数据非递归高效成像方法 |
CN114544775A (zh) * | 2022-02-14 | 2022-05-27 | 浙江大学 | 一种多层结构孔洞缺陷的超声相控阵高效相位偏移成像方法 |
CN114858926A (zh) * | 2022-05-10 | 2022-08-05 | 浙江大学 | 一种基于频域逆时偏移的全矩阵数据高质量成像方法 |
-
2023
- 2023-02-02 CN CN202310052077.2A patent/CN115808469A/zh active Pending
Patent Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN112067698A (zh) * | 2020-09-14 | 2020-12-11 | 南昌航空大学 | 一种时频结合快速全聚焦超声成像方法 |
CN112684005A (zh) * | 2020-12-10 | 2021-04-20 | 苏州热工研究院有限公司 | 基于二维矩阵换能器的全聚焦检测方法 |
CN114544775A (zh) * | 2022-02-14 | 2022-05-27 | 浙江大学 | 一种多层结构孔洞缺陷的超声相控阵高效相位偏移成像方法 |
CN114487117A (zh) * | 2022-02-18 | 2022-05-13 | 浙江大学 | 一种超声相控阵全矩阵数据非递归高效成像方法 |
CN114858926A (zh) * | 2022-05-10 | 2022-08-05 | 浙江大学 | 一种基于频域逆时偏移的全矩阵数据高质量成像方法 |
Non-Patent Citations (6)
Title |
---|
吴海腾: "基于相控阵超声成像的圆柱类部件自动化无损检测理论与实践的研究", 《中国博士学位论文全文数据库 工程科技Ⅰ辑》, no. 6, 15 June 2017 (2017-06-15) * |
张海燕等: "频率-波数域的薄铝板缺陷检测研究", 《声学技术》, vol. 39, no. 4, 31 August 2020 (2020-08-31) * |
曾邱毓晨等: "渐进式频率-波数域全聚焦超声成像", 《声学学报》, vol. 47, no. 5, 30 September 2022 (2022-09-30) * |
李洋等: "基于全矩阵数据的阵列超声检测与评价技术研究和应用概述", 《航空制造技术》, vol. 62, no. 14, 15 July 2019 (2019-07-15) * |
金浩然: "圆柱类部件在线相控阵超声成像理论与技术的研究", 《中国博士学位论文全文数据库 工程科技Ⅰ辑》, no. 12, 15 December 2017 (2017-12-15) * |
钟芳桃等: "航空金属薄板分层缺陷的频率-波数域超声相控阵全聚焦法成像检测", 《无损检测》, vol. 43, no. 7, 10 July 2021 (2021-07-10) * |
Cited By (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN117710184A (zh) * | 2024-01-04 | 2024-03-15 | 浙江大学 | 基于fpga和gpu的超声全聚焦波数域成像方法、系统及装置 |
CN117825512A (zh) * | 2024-01-04 | 2024-04-05 | 浙江大学 | 面向钢轨轨头缺陷的高速超声动态成像检测方法及装置 |
CN117825512B (zh) * | 2024-01-04 | 2024-06-11 | 浙江大学 | 面向钢轨轨头缺陷的高速超声动态成像检测方法及装置 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN115808469A (zh) | 基于cuda的波数域三维超声全矩阵成像方法 | |
US5546807A (en) | High speed volumetric ultrasound imaging system | |
Aykin et al. | Three-dimensional target reconstruction from multiple 2-d forward-scan sonar views by space carving | |
US20100074496A1 (en) | Multi-dimensional empirical mode decomposition (emd) method for image texture analysis | |
US20120013710A1 (en) | System and method for geometric modeling using multiple data acquisition means | |
CN102901444B (zh) | 一种基于mp小波滤波的零件尺寸检测方法及其检测系统 | |
CN101292883B (zh) | 超声三维快速成像方法及其装置 | |
CN101209210B (zh) | 超声波图像取得装置 | |
US20190380685A1 (en) | Methods and system for shading a two-dimensional ultrasound image | |
CN103077559B (zh) | 基于序列图像的果穗三维重建方法 | |
Nocerino et al. | Surface reconstruction assessment in photogrammetric applications | |
CN104090033A (zh) | 基于ebsd图谱的粗晶材料fdtd超声检测仿真模型建立方法 | |
CN101448461B (zh) | 超声波诊断装置及边界提取方法 | |
US20170032702A1 (en) | Method and Apparatus For Generating an Ultrasound Scatterer Representation | |
CN114820952A (zh) | 肺部超声可视化三维重建方法和系统 | |
Gavryushkina | The potential and problems of volumetric 3D modeling in archaeological stratigraphic analysis: A case study from Chlorakas-Palloures, Cyprus | |
Castanié et al. | VolumeExplorer: Roaming large volumes to couple visualization and data processing for oil and gas exploration | |
CN108078590B (zh) | 基于超声频谱多普勒的血流动力学可视化方法与系统 | |
CN114067058A (zh) | 多角度sar立体成像方法 | |
US11484286B2 (en) | Ultrasound evaluation of anatomical features | |
KR101131994B1 (ko) | 원자력 발전소 자동 초음파 신호를 평가하기 위한 실시간 비쥬얼 시스템 | |
Gaits et al. | Ultrasound volume reconstruction from 2D Freehand acquisitions using neural implicit representations | |
CN101635046B (zh) | 基于统一计算设备架构技术的图像处理方法和装置 | |
Yeh et al. | Imaging of internal cracks in concrete structures using the surface rendering technique | |
Frisken et al. | Efficient estimation of 3D euclidean distance fields from 2D range images |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination |