CN104306022B - 用gpu实现压缩感知超声成像的方法 - Google Patents
用gpu实现压缩感知超声成像的方法 Download PDFInfo
- Publication number
- CN104306022B CN104306022B CN201410578176.5A CN201410578176A CN104306022B CN 104306022 B CN104306022 B CN 104306022B CN 201410578176 A CN201410578176 A CN 201410578176A CN 104306022 B CN104306022 B CN 104306022B
- Authority
- CN
- China
- Prior art keywords
- rightarrow
- vector
- centerdot
- matrix
- value
- 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.)
- Expired - Fee Related
Links
Classifications
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B8/00—Diagnosis using ultrasonic, sonic or infrasonic waves
- A61B8/13—Tomography
- A61B8/14—Echo-tomography
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B8/00—Diagnosis using ultrasonic, sonic or infrasonic waves
- A61B8/44—Constructional features of the ultrasonic, sonic or infrasonic diagnostic device
- A61B8/4483—Constructional features of the ultrasonic, sonic or infrasonic diagnostic device characterised by features of the ultrasound transducer
- A61B8/4488—Constructional features of the ultrasonic, sonic or infrasonic diagnostic device characterised by features of the ultrasound transducer the transducer being a phased array
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B8/00—Diagnosis using ultrasonic, sonic or infrasonic waves
- A61B8/52—Devices using data or image processing specially adapted for diagnosis using ultrasonic, sonic or infrasonic waves
- A61B8/5207—Devices using data or image processing specially adapted for diagnosis using ultrasonic, sonic or infrasonic waves involving processing of raw data to produce diagnostic data, e.g. for generating an image
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B8/00—Diagnosis using ultrasonic, sonic or infrasonic waves
- A61B8/52—Devices using data or image processing specially adapted for diagnosis using ultrasonic, sonic or infrasonic waves
- A61B8/5215—Devices using data or image processing specially adapted for diagnosis using ultrasonic, sonic or infrasonic waves involving processing of medical diagnostic data
Landscapes
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Engineering & Computer Science (AREA)
- Medical Informatics (AREA)
- Surgery (AREA)
- Pathology (AREA)
- Radiology & Medical Imaging (AREA)
- Biophysics (AREA)
- Biomedical Technology (AREA)
- Heart & Thoracic Surgery (AREA)
- Physics & Mathematics (AREA)
- Molecular Biology (AREA)
- Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
- Animal Behavior & Ethology (AREA)
- General Health & Medical Sciences (AREA)
- Public Health (AREA)
- Veterinary Medicine (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Gynecology & Obstetrics (AREA)
- Ultra Sonic Daignosis Equipment (AREA)
Abstract
本发明公开了一种用GPU实现压缩感知超声成像的方法,主要解决压缩感知理论框架下成像重建时间较慢的问题。其实现步骤为:1.根据设定的分辨率,对探测区域进行离散化并对该区域进行宽带脉冲扫描,得到回波向量和观测矩阵,进而建立超声成像数学模型;2.将回波向量与观测矩阵分块并复制到GPU显存中;3.在GPU中计算迭代步长;4.将迭代步长带入快速迭代收缩阈值算法求解出重建观测场景散射强度;5.将该散射强度复制到主存中取模值并排列成一个二维矩阵,得到重建的超声图像。本发明相对传统的快速迭代收缩阈值算法,重建时间从分钟级别降低到毫秒级别,极大提高了超声成像的实时性,可用于超声实时处理领域。
Description
技术领域
本发明属于图像处理技术领域,特别涉及一种用GPU实现快速成像的方法,可用于B超成像。
背景技术
医学超声成像经过60余年的发展,它具有相对安全、实时性好、无创、便携、价格低廉等优点,其与X射线诊断技术、电子计算机断层扫描CT成像技术、核磁共振成像技术一起称为现代医学四大影像技术,已令亿万患者受益。
但是超声成像仍然存在一些不足,如分辨率不高,多为毫米级;受噪声干扰严重,图像质量较差;实时性一般。
近年来,在信号处理领域兴起的压缩感知理论吸引了诸多学者的关注,该理论指出,只要信号在某一个空间Ψ上具有稀疏性,就可以利用观测矩阵对其以远低于奈奎斯特采样速率进行观测,进而利用优化手段高概率地从混叠观测中重构原信号,这使得传感器的采样成本大大降低。而且通过恰当选择空间Ψ,信号的稀疏性越强,精确恢复原信号的可能性就越大,这样就会在提高图像分辨率、抑制噪声方面有出色的表现。从近几年国内外发表的文献来看,对压缩感知理论的研究已经涉及众多领域如压缩感知CS雷达成像、医学图像处理、光谱分析、遥感图像处理等,具有非常广阔的应用前景。
由于病灶区域与正常组织的密度特征有明显差别,可认为超声图像在空间域是稀疏的,将压缩感知理论应用到超声成像可以较好的解决超声成像分辨率不高和噪声干扰严重的问题,但是压缩感知理论的问题在于重建过程中观测矩阵维度巨大、导致计算复杂度高,图像的重建时间过长。
针对这个问题,以色列学者A.Beck等人在论文“Afastiterativeshrinkage-thresholdingalgorithmforlinearinverseproblems”SIAMJ.IMAGESCIENCES,Vol.2.No.1,pp.183-202,2009中提出了快速迭代收缩阈值算法,简称FISTA。利用该算法,超声成像的基本框架可以表示为:
其中,X*为重建观测场景散射强度,x为目标向量,b是采样后超声阵元接收的回波数据,λ为正则化参数,表示向量ΨX-b二范数的平方,||x||1表示向量x的一范数。
上述FISTA算法虽然复杂度低,适合求解大规模矩阵的重建问题,同时具备全局收敛性,收敛速度快,使得迭代时间大大缩短。但是,在面对数量级为千兆GB的数据时,仍然需要数分钟到半小时才能恢复一幅图像,这并不能满足对超声成像的实时性要求。
发明内容
本发明的目的在于针对上述已有技术的不足,提出一种用GPU实现压缩感知超声成像的方法,以缩短压缩感知理论应用于超声成像时的重建时间,满足超声成像的实时性要求。
本发明的技术方案是这样实现的:
一.技术原理
2007年,英伟达公司提出了统一计算设备架构CUDA,使图形处理器GPU具备了容易上手的快速并行计算能力,它采用单指令多线程SIMT的编程模式,具备数百倍乃至上千倍CPU核心数的CUDA核心数,如英伟达精视TM系列显卡GTX770拥有1536个CUDA核心,单精度的计算峰值达到3.6TFLOPs。本发明首先建立起压缩感知超声成像的数学模型;再将CUDA技术应用到FISTA算法的重建,在重建之前改进了原FISTA算法中迭代步长的计算方式并实现了并行处理。通过采用矩阵分块思想和CUDA流技术减少数据交换时的延迟,采用指令级并行操作集ILP实现一个线程处理多组数据,调用线性代数函数库实现矩阵乘法运算的性能最优,最终实现了快速成像。
二、技术方案
根据上述原理,本发明的实现步骤如下:
(1)将超声探测区域二维离散化,得到N个离散化的像素点,其中N=T×S,T表示轴向像素的个数,S表示侧向像素的个数;
(2)将超声宽带脉冲发射信号在频域均匀采样得到W个频点,按频点顺序依次对已离散化的二维探测区域进行一次平面波扫描,每次扫描得到一个长度为A的局部观测回波向量bt,并将这W个局部观测回波向量按从上到下顺序合成长度为M=A×W的观测回波向量b={b1,...,bt,...,bW},同时保存由这W个频点产生的回波声场强度矩阵Ψ1,...,Ψt,...,ΨW,其中,A表示超声线阵的阵元个数,矩阵Ψt的宽度为A,长度为N,1≤t≤W;
(3)将回波声场强度矩阵Ψ1,...,Ψt,...,ΨW,按照从上到下顺序排列成一个大小为M×N的观测矩阵Ψ;将离散化的二维探测区域按照行优先的顺序排列成一个目标向量x;
(4)根据回波向量b和观测矩阵Ψ定义基于压缩感知的超声成像数学模型:
其中X*为重建观测场景散射强度,x为目标向量,λ为正则化参数,表示向量Ψx-b二范数的平方,||x||1表示目标向量x的一范数;
(5)在GPU中对上述数学模型进行求解,得到重建观测场景散射强度X*:
(5a)初始化:n=0,ε=10-3,其中n表示第n次迭代,ε表示迭代终止条件;
(5b)将观测矩阵Ψ按行均匀等分成宽度为M/K,长度为N的K个子矩阵Ψ={Ψs1,...,Ψsi,...,ΨsK},要保证每个子矩阵数据量小于显存的容量上限,1≤i≤K;
(5c)将每个子矩阵Ψsi以列优先顺序向量化为一维行向量d_MAi并采用流操作拷贝到显存中,此时观测矩阵Ψ转换成一个宽度为K,长度为M/K×N的向量化矩阵d_MA={d_MA1,...,d_MAi,...,d_MAK},其中观测矩阵Ψ与向量化矩阵它们存在一一对应的关系;将回波向量b等分成长度为M/K的K块,b={b1,...,bi,...,bK},并拷贝到显存中;
(5d)根据向量化矩阵d_MA计算基于梯度下降算法的迭代步长μ:
(5d1)在GPU中建立一个重复执行K次的循环,每次循环创建一个k_step内核函数,求出向量d_MAi对应观测矩阵Ψ每一列元素模值的和;
(5d2)循环结束后,求出向量d_sum的二范数,取倒数即得到迭代步长μ;
(5e)将回波向量b、向量化矩阵d_MA和迭代步长μ带入快速迭代收缩阈值算法中,经过多次梯度下降和快速阈值收缩过程,直到目标向量满足迭代终止条件,得到重建观测场景散射强度X*;
(6)将重建场景散射强度X*拷贝到主机端内存取模值,并按照先行后列的顺序排列成一个二维矩阵,得到重建的超声图像。
本发明与现有技术相比有如下优点:
本发明是基于CUDA架构下的GPU并行计算超声成像过程,通过优化FISTA算法、对观测矩阵分块划分和充分利用GPU中线程块的性能,将成像速度提升了180倍,重建时间由现在的分钟级减少至毫秒级,达到了准实时的超声成像。
附图说明
图1是本发明的实现流程图;
图2是本发明对探测区域离散化的示意图;
图3是本发明使用的超声阵列与探测区域位置示意图;
图4是本发明将回波数据排列成列向量的示意图;
图5是本发明中的矩阵分块示意图;
图6是本发明中进行快速迭代收缩阈值算法的子流程图;
图7是本发明中计算核函数kernel的子流程图;
图8是用本发明仿真得到的重建图像与原始图像的对比图。
具体实施方式
本发明采用CUDA语言,可以在任何一款支持CUDA架构的GPU设备上实现,建议选择计算能力3.0以上的显卡,计算能力越强,对显卡优化越好。为满足大规模矩阵向量乘法的要求,建议显存容量在2GB以上,可以搭载多块GPU实现显存容量高于观测矩阵的大小,这样可以减少数据交换带来的延迟。
以下结合附图对本发明作进一步详细描述。
参照图1,本发明的实现步骤如下:
步骤一:对探测区域进行离散化。
本发明对探测区域离散化采用等间隔采样,得到离散化的二维探测区域N=T×S,如图2所示,其中N表示离散化总像素点数,T表示轴向像素点数,S表示侧向像素点数。
步骤二:对二维离散化的探测区域进行宽带脉冲平面波扫描,得到回波向量b和观测矩阵Ψ。
具体实现如下:
(2a)根据图3所示的超声阵列与探测区域的相对位置建立直角坐标系(x,y),其中x表示侧向,y表示轴向;
(2b)在直角坐标系中,将超声阵列固定在轴向坐标为零的位置,即y=0,且将阵列中心与探测区域中心对齐,超声阵列的长度为(A-1)×d,则第l个阵元的横坐标xl为:
其中A为超声线阵的阵元个数,d为相邻两个阵元之间的间隔,1≤l≤A;
(2c)将超声宽带脉冲发射信号在频域均匀采样得到W个频点,按频点顺序依次对已离散化的二维探测区域进行一次平面波扫描,每次扫描得到一个长度为A的局部观测回波向量bt,1≤t≤W,并将这W个局部观测回波向量按图4所示从上到下顺序合成长度为M=A×W的观测回波向量b={b1,...,bt,...,bW};
(2d)计算探测区域中各个像素点在第m个超声线阵阵元处产生的声场强度
其中Ain(kt)表示超声宽带脉冲发射信号在频点取值为kt时的幅度;表示超声宽带脉冲发射信号在频点取值为kt时离散二维探测区域中各个像素点返回的相位,表示超声宽带脉冲发射信号方位角的单位向量,指定为轴向方向;表示超声宽带脉冲发射信号从超声线阵轴心位置到离散二维探测区域各个像素点距离的向量;表示格林函数,表示超声线阵轴心位置到到超声线阵各个阵元距离的向量,1≤m≤A;
(2e)由所述声场强度计算出当频点为kt时,宽度为A,长度为N回波声场强度矩阵Ψt:
其中1≤t≤W,1≤j≤N,1≤m≤A;
(2f)将回波声场强度矩阵Ψ1,...,Ψt,...,ΨW按照从上到下顺序排列成一个大小为M×N的观测矩阵Ψ;
(2g)将离散化的二维探测区域按照行优先顺序排列成一个长度为N的目标向量x。
步骤三:构建基于压缩感知的超声成像数学模型。
其中X*为重建观测场景散射强度,x为目标向量,λ为正则化参数,表示向量Ψx-b二范数的平方,||x||1表示目标向量x的一范数。
步骤四:在GPU中对上述数学模型进行求解,得到重建观测场景散射强度X*。
具体实现如下:
(4a)初始化:n=0,ε=10-3,其中n表示第n次迭代,ε表示迭代终止条件;
(4b)如图5所示,将观测矩阵Ψ按行均匀等分成宽度为M/K,长度为N的K个子矩阵,Ψ={Ψs1,...,Ψsi,...,ΨsK},要保证每个子矩阵的数据量小于显存的容量上限,1≤i≤K;
(4c)将子矩阵Ψsi以列优先顺序向量化为一维行向量d_MAi并拷贝到显存中,此时观测矩阵Ψ转换成一个宽度为K,长度为M/K×N的向量化矩阵d_MA={d_MA1,...,d_MAi,...,d_MAK},它们存在一一对应的关系;将回波向量b等分成长度为M/K的K块,b={b1,...,bi,...,bK},拷贝到显存中;
(4d)根据向量化矩阵d_MA计算基于梯度下降算法的迭代步长μ:
(4d1)在GPU中建立一个重复执行K次的循环,每次循环创建一个k_step内核函数求出向量d_MAi对应观测矩阵Ψ每一列元素绝对值的和,具体步骤如下:
(4d1.a)在GPU中开辟N/V个线程块和长度为N的向量d_sum,每个线程块中设有M/K个线程,其中向量V要满足V对N取余值为0,V取值为2~4可以达到理想的性能,表示实数域,1≤q≤N,1≤V≤N;
(4d1.b)用第z个线程块读取M/K×V个元素,即向量d_MAi中第M/K×V×(z-1)+1到第M/K×V×z的元素;在该线程块中分V次取出这些元素,每次按顺序取M/K个元素,依次分配给线程1~线程M/K,每个线程计算各自元素的模值;将这M/K个模值累加给向量d_sum的第d_sumV×(z-1)+p个元素,其中1≤z≤N/V,1≤p≤V;
(4d2)循环结束后,求出向量d_sum的二范数,取倒数即得到迭代步长μ。
(4e)将回波b、观测矩阵Ψ和迭代步长μ带入快速迭代收缩阈值算法中,进行梯度下降和快速阈值收缩,得到重建观测场景散射强度X*。
参见图6,本发明步骤具体步骤如下:
(4e1)将回波向量b、向量化矩阵d_MA和迭代步长μ带入下式,更新梯度下降序列un:
un=yn-μ×cublas(d_MAH·(cublas(d_MA·yn)-b)),
其中表示复数域;yn是快速收缩变量,初始值为0,长度为N;d_MAH表示向量化矩阵d_MA的共轭转置;μ为迭代步长;cublas是CUDA线性代数库CUBLASLibrary中所支持的库函数,它可以执行矩阵与向量乘法运算,1≤j≤N;
(4e2)参见图7,在kernel内核函数中开辟N个线程,将梯度下降序列un带入kernel内核函数,用第j个线程对梯度下降序列un的元素uj进行阈值操作,得到目标元素xj:xj=SΓ(uj),其中SΓ为阈值函数:
其中Γ为阈值,Γ=λμ,λ取值在2e4~4e4之间,e表示科学计数,取值为10;sign()表示取符号函数,1≤j≤N;
将目标元素x1,...,xj,...,xN组合为目标向量xn,
(4e3)判断收敛条件||xn-xn-1||2<ε是否成立:
若成立,则停止计算,重建观测场景散射强度X*=xn;
若不成立,令n=n+1,更新快速收缩变量yn为:
yn=xn-1+(xn-2-xn-1)×(1-t1)/t2,
其中当n=1时,x0=0;j表示系数向量xn和xn-1中的第j个元素,xn[j]表示向量xn的第j个元素的值,xn-1[j]表示向量xn-1的第j个元素值;
其中t1、t2是两个数值不同的加速因子,t1初始值为1,将t1更新为t1=t2,返回步骤(4e1)。
步骤五:将重建场景散射强度X*拷贝到主机端内存取模值,并按照先行后列的顺序排列成一个二维矩阵,即得到重建的超声图像。
本发明的效果可以通过以下仿真实验说明:
1.仿真内容
为了验证本发明在CUDA架构下并行计算的优越性,通过一组仿真实验对本发明方法所需的时间与CPU串行所需的时间定量分析,仿真测试平台为:
CPU:IntelCorei74770,内存:16GB
显卡:NVIDIAGeForceGTX770,显卡显存:4GB,计算能力3.0
软件平台:windows764位操作系统,VS2010x64编译器,CUDA版本:5.5
2.仿真结果及分析
仿真实验用到两组数据,第一组数据观测矩阵大小为1GB,图像大小为128×128,目标点个数为20个;第二组数据观测矩阵大小为8GB,图像大小为256×512,目标点个数为100个。图像的分辨率为100um,要高于传统超声成像分辨率,迭代次数设定为30次,第一组仿真得到的重建图像与原始图像的对比如图8所示。
图8表明,CUDA架构下GPU的重建图像与原始图像在主观上评价基本一致,点目标的位置与原始图像完全吻合,从而达到了较高的重建效果。
将CUDA架构下GPU并行计算时间和CPU串行时间进行对比,如表1和表2所示:
表1观测矩阵为1GB重建时间
从表1看出,本发明方法对超声成像的重构速度提高明显,达到了每秒2.5帧图像,基本上解决了压缩感知框架下处理大规模矩阵耗时较长的问题,实现了准实时性。
表2观测矩阵为8GB重建时间
表2中,GPU的加速效果比表1中的加速比差,其主要原因在于显存不足导致每次迭代都需要进行数据交换,将数据从主机端内存交换到显存端。尽管如此,GPU加速处理后仍然比CPU串行快13倍,而且通过分析数据可以发现在CPU+GPU的方案中数据交换的时间占到总时间的92%,因此,若可以利用提升硬件解决显存不足的问题,加速比仍然可以达到180倍左右。
Claims (5)
1.一种用GPU实现压缩感知超声成像的方法,其特征在于,包括以下步骤:
(1)将超声探测区域二维离散化,得到N个离散化的像素点,其中N=T×S,T表示轴向像素的个数,S表示侧向像素的个数;
(2)将超声宽带脉冲发射信号在频域均匀采样得到W个频点,按频点顺序依次对已离散化的二维探测区域进行一次平面波扫描,每次扫描得到一个长度为A的局部观测回波向量bt,并将这W个局部观测回波向量按从上到下顺序合成长度为M=A×W的观测回波向量b={b1,...,bt,...,bW},同时保存由这W个频点产生的回波声场强度矩阵Ψ1,...,Ψt,...,ΨW,其中,A表示超声线阵的阵元个数,矩阵Ψt的宽度为A,长度为N,1≤t≤W;
(3)将回波声场强度矩阵Ψ1,...,Ψt,...,ΨW,按照从上到下顺序排列成一个大小为M×N的观测矩阵Ψ;将离散化的二维探测区域按照行优先的顺序排列成一个目标向量x;
(4)根据回波向量b和观测矩阵Ψ定义基于压缩感知的超声成像数学模型:
其中X*为重建观测场景散射强度,x为目标向量,λ为正则化参数,表示向量Ψx-b二范数的平方,||x||1表示目标向量x的一范数;
(5)在GPU中对上述数学模型进行求解,得到重建观测场景散射强度X*:
(5a)初始化:n=0,ε=10-3,其中n表示第n次迭代,ε表示迭代终止条件;
(5b)将观测矩阵Ψ按行均匀等分成宽度为M/K,长度为N的K个子矩阵Ψ={Ψs1,...,Ψsi,...,ΨsK},要保证每个子矩阵数据量小于显存的容量上限,1≤i≤K;
(5c)将每个子矩阵Ψsi以列优先顺序向量化为一维行向量d_MAi并采用流操作拷贝到显存中,此时观测矩阵Ψ转换成一个宽度为K,长度为M/K×N的向量化矩阵d_MA={d_MA1,...,d_MAi,...,d_MAK},其中观测矩阵Ψ与向量化矩阵它们存在一一对应的关系;将回波向量b等分成长度为M/K的K块,b={b1,...,bi,...,bK},并拷贝到显存中;
(5d)根据向量化矩阵d_MA计算基于梯度下降算法的迭代步长μ:
(5d1)在GPU中建立一个重复执行K次的循环,每次循环创建一个k_step内核函数求出向量d_MAi对应观测矩阵Ψ每一列元素绝对值的和,具体步骤如下:
(5d1.a)在GPU中开辟N/V个线程块和长度为N的向量d_sum,每个线程块中设有M/K个线程,其中向量V要满足V对N取余值为0,V取值为2~4可以达到理想的性能,表示实数域,1≤q≤N,1≤V≤N;
(5d1.b)用第z个线程块读取M/K×V个元素,即向量d_MAi中第M/K×V×(z-1)+1到第M/K×V×z的元素;在该线程块中分V次取出这些元素,每次按顺序取M/K个元素,依次分配给线程1~线程M/K,每个线程计算各自元素的模值;将这M/K个模值累加给向量d_sum的第d_sumV×(z-1)+p个元素,其中1≤z≤N/V,1≤p≤V;
(5d2)循环结束后,求出向量d_sum的二范数,取倒数即得到迭代步长μ;
(5e)将回波向量b、向量化矩阵d_MA和迭代步长μ带入快速迭代收缩阈值算法中,经过多次梯度下降和快速阈值收缩过程,直到目标向量满足迭代终止条件,得到重建观测场景散射强度X*;
(6)将重建场景散射强度X*拷贝到主机端内存取模值,并按照先行后列的顺序排列成一个二维矩阵,得到重建的超声图像。
2.如权利要求1所述的用GPU实现压缩感知超声成像方法,其特征在于,所述步骤(2)中,频点取kt时的回波声场强度矩阵Ψt,通过下式计算:
其中,Ain(kt)表示超声宽带脉冲发射信号在频点取值为kt时的幅度;表示超声宽带脉冲发射信号在频点取值为kt时离散二维探测区域中各个像素点返回的相位,表示超声宽带脉冲发射信号方位角的单位向量,指定为轴向方向;表示超声宽带脉冲发射信号从超声线阵轴心位置到离散二维探测区域各个像素点距离的向量;
表示格林函数,表示超声线阵轴心位置到到超声线阵各个阵元距离的向量,1≤t≤W,1≤j≤N,1≤m≤A。
3.如权利要求1所述的用GPU实现压缩感知超声成像方法,其特征在于,所述步骤(5d1)中,常量V的取值要满足V对N取余值为0,V取值为2~4。
4.如权利要求1所述的用GPU实现压缩感知超声成像方法,其特征在于,所述步骤(5e),按如下步骤进行:
(5e1)将回波向量b、向量化矩阵d_MA和迭代步长μ带入下式,更新梯度下降序列un:
un=yn-μ×cublas(d_MAH·(cublas(d_MA·yn)-b))3)
其中yn是快速收缩向量,初始值为0,长度为N;d_MAH表示向量化矩阵d_MA的共轭转置;μ为迭代步长;cublas是CUDA线性代数库CUBLASLibrary中所支持的库函数,它可以执行向量化矩阵与向量的乘法运算;
(5e2)将梯度下降序列un带入下式,更新目标场景散射强度xn:
xn=SΓ(un)4)
其中SΓ为阈值函数:
其中Γ为阈值,Γ=λμ,λ取值在2e4~4e4之间,e表示科学计数,取值为10;sign()表示取符号函数;
(5e3)判断收敛条件||xn-xn-1||2<ε是否成立:
若成立,则停止计算,重建观测场景散射强度X*=xn;
若不成立,令n=n+1,更新快速收缩向量yn为:
yn=xn-1+(xn-2-xn-1)×(1-t1)/t25)
其中当n=1时,x0=0;j表示系数向量xn和xn-1中的第j个元素,xn[j]表示向量xn的第j个元素的值,xn-1[j]表示向量xn-1的第j个元素值;t1、t2是两个数值不同的加速因子,t1的初始值为1,将t1更新为t1=t2,返回步骤(5e1)。
5.如权利要求4所述的用GPU实现压缩感知超声成像方法,其特征在于,所述步骤(5e2)中,更新目标向量xn的计算是在内核函数kernel中完成的,即xn=SΓ(un)通过使用多线程进行并行计算得到。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201410578176.5A CN104306022B (zh) | 2014-10-24 | 2014-10-24 | 用gpu实现压缩感知超声成像的方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201410578176.5A CN104306022B (zh) | 2014-10-24 | 2014-10-24 | 用gpu实现压缩感知超声成像的方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN104306022A CN104306022A (zh) | 2015-01-28 |
CN104306022B true CN104306022B (zh) | 2016-05-25 |
Family
ID=52361446
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201410578176.5A Expired - Fee Related CN104306022B (zh) | 2014-10-24 | 2014-10-24 | 用gpu实现压缩感知超声成像的方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN104306022B (zh) |
Families Citing this family (11)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN104777484B (zh) * | 2015-02-13 | 2016-04-27 | 西安交通大学 | 压缩自适应波束合成的平面波超声成像和微泡成像的方法与系统 |
CN104688271B (zh) * | 2015-03-27 | 2017-04-26 | 清华大学 | 合成聚焦超声成像方法和装置 |
CN105891808B (zh) * | 2016-05-27 | 2019-02-26 | 复旦大学 | 一种多目标定位的声波发射装置 |
CN106597388B (zh) * | 2016-11-24 | 2019-05-03 | 北京华航无线电测量研究所 | 一种两侧滑窗取平均一维检测的fpga实现方法 |
CN106802418B (zh) * | 2017-01-19 | 2019-05-03 | 重庆大学 | 一种合成孔径压缩感知超声成像中的高效能稀疏字典的设计方法 |
CN108670304B (zh) * | 2018-06-06 | 2021-03-02 | 东北大学 | 一种基于改进dmas算法的超声平面波成像方法 |
CN109276276B (zh) * | 2018-08-24 | 2021-06-08 | 广东省医疗器械质量监督检验所 | 基于Labview平台的超声内窥成像系统及方法 |
CN109709547A (zh) * | 2019-01-21 | 2019-05-03 | 电子科技大学 | 一种实波束扫描雷达加速超分辨成像方法 |
CN110780274B (zh) * | 2019-11-04 | 2022-04-01 | 电子科技大学 | 一种用于扫描雷达的改进l1正则化方位超分辨成像方法 |
CN114216582B (zh) * | 2021-12-08 | 2022-09-06 | 电子科技大学长三角研究院(湖州) | 一种高精度快速温度场重建方法、系统、设备及终端 |
CN116244078B (zh) * | 2023-02-27 | 2023-12-01 | 青岛中海潮科技有限公司 | 基于多线程和simd的水下声场快速计算方法 |
Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102058416A (zh) * | 2010-12-14 | 2011-05-18 | 哈尔滨工业大学 | 基于压缩感知的微波热声成像装置及方法 |
CN103235298A (zh) * | 2013-05-08 | 2013-08-07 | 西安电子科技大学 | 基于稀疏阵列的微波关联成像系统与成像方法 |
CN103983968A (zh) * | 2014-03-20 | 2014-08-13 | 西安电子科技大学 | 基于分布式压缩感知的全极化sar超分辨成像方法 |
Family Cites Families (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
EP1839264A2 (en) * | 2005-01-18 | 2007-10-03 | Trestle Corporation | System and method for creating variable quality images of a slide |
US8435180B2 (en) * | 2007-09-17 | 2013-05-07 | Siemens Medical Solutions Usa, Inc. | Gain optimization of volume images for medical diagnostic ultrasonic imaging |
ITGE20120048A1 (it) * | 2012-05-04 | 2013-11-05 | Esaote Spa | Metodo di ricostruzione di immagini biomediche |
-
2014
- 2014-10-24 CN CN201410578176.5A patent/CN104306022B/zh not_active Expired - Fee Related
Patent Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102058416A (zh) * | 2010-12-14 | 2011-05-18 | 哈尔滨工业大学 | 基于压缩感知的微波热声成像装置及方法 |
CN103235298A (zh) * | 2013-05-08 | 2013-08-07 | 西安电子科技大学 | 基于稀疏阵列的微波关联成像系统与成像方法 |
CN103983968A (zh) * | 2014-03-20 | 2014-08-13 | 西安电子科技大学 | 基于分布式压缩感知的全极化sar超分辨成像方法 |
Non-Patent Citations (4)
Title |
---|
A iterative thresholding algorithm for linear inverse problems with multi-constraints and its applications;Sergry voronin,etc;《ELSEVIER》;20121231;第109-130页 * |
A new iterative firm-thresholding algorithm for inverse problems with sparsity constraints;Sergey Voronin,etc;《ElSEVIER》;20131231;第151-164页 * |
Iterative thresholding algorithms;Massimo Fornasier,etc;《ELSEVIER》;20081231;第187-208页 * |
压缩感知理论及其研究进展;石光明,刘丹华,等;《电子学报》;20090531;第37卷(第5期);第1070-1081页 * |
Also Published As
Publication number | Publication date |
---|---|
CN104306022A (zh) | 2015-01-28 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN104306022B (zh) | 用gpu实现压缩感知超声成像的方法 | |
CN104306023B (zh) | 基于压缩感知的超声成像快速实现方法 | |
Wang et al. | Seismic trace interpolation for irregularly spatial sampled data using convolutional autoencoder | |
CN104933683B (zh) | 一种用于磁共振快速成像的非凸低秩重建方法 | |
CN109683161A (zh) | 一种基于深度admm网络的逆合成孔径雷达成像的方法 | |
CN106934778B (zh) | 一种基于小波域结构和非局部分组稀疏的mr图像重建方法 | |
Wang et al. | Three-dimensional reconstruction from a multiview sequence of sparse ISAR imaging of a space target | |
CN107798697A (zh) | 一种基于卷积神经网络的医学图像配准方法、系统及电子设备 | |
CN107657217A (zh) | 基于运动目标检测的红外与可见光视频的融合方法 | |
CN106663316A (zh) | 一种基于块稀疏压缩感知的红外图像重构方法及其系统 | |
CN109377520A (zh) | 基于半监督循环gan的心脏图像配准系统及方法 | |
CN103455989B (zh) | 一种结合超声图像提高有限角度ct成像质量的方法 | |
CN104181528B (zh) | 基于bp优化的压缩感知多层isar成像方法 | |
CN109557540A (zh) | 基于目标散射系数非负约束的全变差正则化关联成像方法 | |
Shi et al. | A fast and accurate basis pursuit denoising algorithm with application to super-resolving tomographic SAR | |
CN106054189A (zh) | 基于dpKMMDP模型的雷达目标识别方法 | |
Janiszewski | Asymptotically hyperbolic black holes in Horava gravity | |
CN105662357A (zh) | 磁共振温度成像方法及系统 | |
CN106037795A (zh) | 一种优化全波反演法重建声速图像的方法 | |
Xie et al. | Super-resolution of Pneumocystis carinii pneumonia CT via self-attention GAN | |
CN102183761B (zh) | 星载干涉合成孔径雷达数字高程模型重建方法 | |
JP5451414B2 (ja) | 被検体情報処理装置および被検体情報処理方法 | |
CN104424625A (zh) | 一种gpu加速cbct图像重建方法和装置 | |
CN106769734A (zh) | 一种超声波聚焦式河流泥沙浓度在线测量方法 | |
CN103985111B (zh) | 一种基于双字典学习的4d‑mri超分辨率重构方法 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
C14 | Grant of patent or utility model | ||
GR01 | Patent grant | ||
CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20160525 Termination date: 20211024 |
|
CF01 | Termination of patent right due to non-payment of annual fee |