本发明公开了一种用GPU实现压缩感知超声成像的方法，主要解决压缩感知理论框架下成像重建时间较慢的问题。其实现步骤为：1.根据设定的分辨率，对探测区域进行离散化并对该区域进行宽带脉冲扫描，得到回波向量和观测矩阵，进而建立超声成像数学模型；2.将回波向量与观测矩阵分块并复制到GPU显存中；3.在GPU中计算迭代步长；4.将迭代步长带入快速迭代收缩阈值算法求解出重建观测场景散射强度；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线性代数库CUBLAS Library中所支持的库函数，它可以执行向量化矩阵与向量的乘法运算；
(5e2)将梯度下降序列un带入下式，更新目标场景散射强度xn：xn＝SΓ(un) 4)
其中Γ为阈值，Γ＝λμ，λ取值在2e4～4e4之间，e表示科学计数，取值为10；sign()表示取符号函数；
(5e3)判断收敛条件||xn-xn-1||2<ε是否成立：
若成立，则停止计算，重建观测场景散射强度X*＝xn；
yn＝xn-1+(xn-2-xn-1)×(1-t1)/t2 5)其中 当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)通过使用多线程进行并行计算得到。
用GPU实现压缩感知超声成像的方法
技术领域
[0001] 本发明属于图像处理技术领域，特别涉及一种用GPU实现快速成像的方法，可用于B超成像。
背景技术
[0002] 医学超声成像经过60余年的发展，它具有相对安全、实时性好、无创、便携、价格低廉等优点，其与X射线诊断技术、电子计算机断层扫描CT成像技术、核磁共振成像技术一起称为现代医学四大影像技术，已令亿万患者受益。
[0003] 但是超声成像仍然存在一些不足，如分辨率不高，多为毫米级；受噪声干扰严重，图像质量较差；实时性一般。
[0004] 近年来，在信号处理领域兴起的压缩感知理论吸引了诸多学者的关注，该理论指出，只要信号在某一个空间Ψ上具有稀疏性，就可以利用观测矩阵对其以远低于奈奎斯特采样速率进行观测，进而利用优化手段高概率地从混叠观测中重构原信号，这使得传感器的采样成本大大降低。而且通过恰当选择空间Ψ，信号的稀疏性越强，精确恢复原信号的可能性就越大，这样就会在提高图像分辨率、抑制噪声方面有出色的表现。从近几年国内外发表的文献来看，对压缩感知理论的研究已经涉及众多领域如压缩感知CS雷达成像、医学图像处理、光谱分析、遥感图像处理等，具有非常广阔的应用前景。
[0005] 由于病灶区域与正常组织的密度特征有明显差别，可认为超声图像在空间域是稀疏的，将压缩感知理论应用到超声成像可以较好的解决超声成像分辨率不高和噪声干扰严重的问题，但是压缩感知理论的问题在于重建过程中观测矩阵维度巨大、导致计算复杂度高，图像的重建时间过长。
[0006] 针对这个问题，以色列学者A.Beck等人在论文“A fast iterative shrinkage-thresholding algorithm for linear inverse problems”SIAM J.IMAGE SCIENCES,Vol.2.No.1,pp.183-202,2009中提出了快速迭代收缩阈值算法，简称FISTA。利用该算法，超声成像的基本框架可以表示为：
[0007] 其中，X*为重建观测场景散射强度，x为目标向量，b是采样后超声阵元接收的回波数据，λ为正则化参数， 表示向量ΨX-b二范数的平方，||x||1表示向量x的一范数。
[0008] 上述FISTA算法虽然复杂度低，适合求解大规模矩阵的重建问题，同时具备全局收敛性，收敛速度快，使得迭代时间大大缩短。但是，在面对数量级为千兆GB的数据时，仍然需要数分钟到半小时才能恢复一幅图像，这并不能满足对超声成像的实时性要求。
发明内容
[0009] 本发明的目的在于针对上述已有技术的不足，提出一种用GPU实现压缩感知超声成像的方法，以缩短压缩感知理论应用于超声成像时的重建时间，满足超声成像的实时性要求。
[0010] 本发明的技术方案是这样实现的：
[0011] 一.技术原理
[0012] 2007年，英伟达 公司提出了统一计算设备架构CUDA，使图形处理器GPU具备了容易上手的快速并行计算能力，它采用单指令多线程SIMT的编程模式，具备数百倍乃至上千倍CPU核心数的CUDA核心数，如英伟达精视TM系列显卡GTX770拥有1536个CUDA核心，单精度的计算峰值达到3.6TFLOPs。本发明首先建立起压缩感知超声成像的数学模型；再将CUDA技术应用到FISTA算法的重建，在重建之前改进了原FISTA算法中迭代步长的计算方式并实现了并行处理。通过采用矩阵分块思想和CUDA流技术减少数据交换时的延迟，采用指令级并行操作集ILP实现一个线程处理多组数据，调用线性代数函数库实现矩阵乘法运算的性能最优，最终实现了快速成像。
[0013] 二、技术方案
[0014] 根据上述原理，本发明的实现步骤如下：
[0015] (1)将超声探测区域二维离散化，得到N个离散化的像素点，其中N＝T×S，T表示轴向像素的个数，S表示侧向像素的个数；
[0016] (2)将超声宽带脉冲发射信号在频域均匀采样得到W个频点，按频点顺序依次对已离散化的二维探测区域进行一次平面波扫描，每次扫描得到一个长度为A的局部观测回波向量bt，并将这W个局部观测回波向量按从上到下顺序合成长度为M＝A×W的观测回波向量b＝{b1,...,bt,...,bW}，同时保存由这W个频点产生的回波声场强度矩阵Ψ1,...,Ψt,...,ΨW，其中，A表示超声线阵的阵元个数，矩阵Ψt的宽度为A，长度为N，1≤t≤W；
[0017] (3)将回波声场强度矩阵Ψ1,...,Ψt,...,ΨW，按照从上到下顺序排列成一个大小为M×N的观测矩阵Ψ；将离散化的二维探测区域按照行优先的顺序排列成一个目标向量x；
[0018] (4)根据回波向量b和观测矩阵Ψ定义基于压缩感知的超声成像数学模型：
[0019]
[0020] 其中X*为重建观测场景散射强度，x为目标向量，λ为正则化参数， 表示向量Ψx-b二范数的平方，||x||1表示目标向量x的一范数；
[0021] (5)在GPU中对上述数学模型进行求解，得到重建观测场景散射强度X*：
[0022] (5a)初始化：n＝0，ε＝10-3,其中n表示第n次迭代，ε表示迭代终止条件；
[0023] (5b)将观测矩阵Ψ按行均匀等分成宽度为M/K，长度为N的K个子矩阵Ψ＝{Ψs1,...,Ψsi,...,ΨsK}，要保证每个子矩阵数据量小于显存的容量上限，1≤i≤K；
[0024] (5c)将每个子矩阵Ψsi以列优先顺序向量化为一维行向量d_MAi并采用流操作拷贝到显存中，此时观测矩阵Ψ转换成一个宽度为K，长度为M/K×N的向量化矩阵d_MA＝{d_MA1,...,d_MAi,...,d_MAK}，其中观测矩阵Ψ与向量化矩阵它们存在一一对应的关系；将回波向量b等分成长度为M/K的K块，b＝{b1,...,bi,...,bK}，并拷贝到显存中；
[0025] (5d)根据向量化矩阵d_MA计算基于梯度下降算法的迭代步长μ：
[0026] (5d1)在GPU中建立一个重复执行K次的循环，每次循环创建一个k_step内核函数,求出向量d_MAi对应观测矩阵Ψ每一列元素模值的和；
[0027] (5d2)循环结束后，求出向量d_sum的二范数，取倒数即得到迭代步长μ；
[0028] (5e)将回波向量b、向量化矩阵d_MA和迭代步长μ带入快速迭代收缩阈值算法中，经过多次梯度下降和快速阈值收缩过程，直到目标向量满足迭代终止条件，得到重建观测场景散射强度X*；
[0029] (6)将重建场景散射强度X*拷贝到主机端内存取模值，并按照先行后列的顺序排列成一个二维矩阵，得到重建的超声图像。
[0030] 本发明与现有技术相比有如下优点：
[0031] 本发明是基于CUDA架构下的GPU并行计算超声成像过程，通过优化FISTA算法、对观测矩阵分块划分和充分利用GPU中线程块的性能，将成像速度提升了180倍，重建时间由现在的分钟级减少至毫秒级，达到了准实时的超声成像。
附图说明
[0032] 图1是本发明的实现流程图；
[0033] 图2是本发明对探测区域离散化的示意图；
[0034] 图3是本发明使用的超声阵列与探测区域位置示意图；
[0035] 图4是本发明将回波数据排列成列向量的示意图；
[0036] 图5是本发明中的矩阵分块示意图；
[0037] 图6是本发明中进行快速迭代收缩阈值算法的子流程图；
[0038] 图7是本发明中计算核函数kernel的子流程图；
[0039] 图8是用本发明仿真得到的重建图像与原始图像的对比图。
具体实施方式
[0040] 本发明采用CUDA语言，可以在 任何一款支持CUDA架构的GPU设备上实现，建议选择计算能力3.0以上的显卡，计算能力越强，对显卡优化越好。为满足大规模矩阵向量乘法的要求，建议显存容量在2GB以上，可以搭载多块GPU实现显存容量高于观测矩阵的大小，这样可以减少数据交换带来的延迟。
[0041] 以下结合附图对本发明作进一步详细描述。
[0042] 参照图1，本发明的实现步骤如下：
[0043] 步骤一：对探测区域进行离散化。
[0044] 本发明对探测区域离散化采用等间隔采样，得到离散化的二维探测区域N＝T×S，如图2所示，其中N表示离散化总像素点数，T表示轴向像素点数，S表示侧向像素点数。
[0045] 步骤二：对二维离散化的探测区域进行宽带脉冲平面波扫描，得到回波向量b和观测矩阵Ψ。
[0046] 具体实现如下：
[0047] (2a)根据图3所示的超声阵列与探测区域的相对位置建立直角坐标系(x,y)，其中x表示侧向，y表示轴向；
[0048] (2b)在直角坐标系中，将超声阵列固定在轴向坐标为零的位置，即y＝0，且将阵列中心与探测区域中心对齐，超声阵列的长度为(A-1)×d，则第l个阵元的横坐标xl为：
[0049]
[0050] 其中A为超声线阵的阵元个数，d为相邻两个阵元之间的间隔，1≤l≤A；
[0051] (2c)将超声宽带脉冲发射信号在频域均匀采样得到W个频点，按频点顺序依次对已离散化的二维探测区域进行一次平面波扫描，每次扫描得到一个长度为A的局部观测回波向量bt，1≤t≤W，并将这W个局部观测回波向量按图4所示从上到下顺序合成长度为M＝A×W的观测回波向量b＝{b1,...,bt,...,bW}；
[0052] (2d)计算探测区域中各个像素点在第m个超声线阵阵元处产生的声场强度[0053]
[0054] 其中Ain(kt)表示超声宽带脉冲发射信号在频点取值为kt时的幅度； 表示超声宽带脉冲发射信号在频点取值为kt时离散二维探测区域中各个像素点返回的相位， 表示超声宽带脉冲发射信号方位角的单位向量，指定为轴向方向；表示超声宽带脉冲发射信号从超声线阵轴心位置到离散二维探测区域各个像素点距离的向量； 表示格林函数， 表示超声线阵轴心位置到到超声线阵各个阵元距离的向量，1≤m≤A；
[0055] (2e)由所述声场强度 计算出当频点为kt时，宽度为A，长度为N回波声场强度矩阵Ψt：
[0056]
[0057] 其中1≤t≤W,1≤j≤N,1≤m≤A；
[0058] (2f)将回波声场强度矩阵Ψ1,...,Ψt,...,ΨW按照从上到下顺序排列成一个大小为M×N的观测矩阵Ψ；
[0059] (2g)将离散化的二维探测区域按照行优先顺序排列成一个长度为N的目标向量x。
[0060] 步骤三：构建基于压缩感知的超声成像数学模型。
[0061]
[0062] 其中X*为重建观测场景散射强度，x为目标向量，λ为正则化参数， 表示向量Ψx-b二范数的平方，||x||1表示目标向量x的一范数。
[0063] 步骤四：在GPU中对上述数学模型进行求解，得到重建观测场景散射强度X*。
[0064] 具体实现如下：
[0065] (4a)初始化：n＝0，ε＝10-3,其中n表示第n次迭代，ε表示迭代终止条件；
[0066] (4b)如图5所示，将观测矩阵Ψ按行均匀等分成宽度为M/K，长度为N的K个子矩阵，Ψ＝{Ψs1,...,Ψsi,...,ΨsK}，要保证每个子矩阵的数据量小于显存的容量上限，1≤i≤K；
[0067] (4c)将子矩阵Ψsi以列优先顺序向量化为一维行向量d_MAi并拷贝到显存中，此时观测矩阵Ψ转换成一个宽度为K，长度为M/K×N的向量化矩阵d_MA＝{d_MA1,...,d_MAi,...,d_MAK}，它们存在一一对应的关系；将回波向量b等分成长度为M/K的K块，b＝{b1,...,bi,...,bK}，拷贝到显存中；
[0068] (4d)根据向量化矩阵d_MA计算基于梯度下降算法的迭代步长μ：
[0069] (4d1)在GPU中建立一个重复执行K次的循环，每次循环创建一个k_step内核函数求出向量d_MAi对应观测矩阵Ψ每一列元素绝对值的和，具体步骤如下：
[0070] (4d1.a)在GPU中开辟N/V个线程块和长度为N的向量d_sum，每个线程块中设有M/K个线程，其中向量 V要满足V对N取余值为0，V取值为2～4可以达到理想的性能， 表示实数域，1≤q≤N，1≤V≤N；
[0071] (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；
[0072] (4d2)循环结束后，求出向量d_sum的二范数，取倒数即得到迭代步长μ。
[0073] (4e)将回波b、观测矩阵Ψ和迭代步长μ带入快速迭代收缩阈值算法中，进行梯度下降和快速阈值收缩，得到重建观测场景散射强度X*。
[0074] 参见图6，本发明步骤具体步骤如下：
[0075] (4e1)将回波向量b、向量化矩阵d_MA和迭代步长μ带入下式，更新梯度下降序列un：
[0076] un＝yn-μ×cublas(d_MAH·(cublas(d_MA·yn)-b))，
[0077] 其中 表示复数域；yn是快速收缩变量，初始值为0,长度为N；d_MAH表示向量化矩阵d_MA的共轭转置；μ为迭代步长；cublas是CUDA线性代数库CUBLASLibrary中所支持的库函数，它可以执行矩阵与向量乘法运算，1≤j≤N；
[0078] (4e2)参见图7，在kernel内核函数中开辟N个线程，将梯度下降序列un带入kernel内核函数，用第j个线程对梯度下降序列un的元素uj进行阈值操作，得到目标元素xj：xj＝SΓ(uj)，其中SΓ为阈值函数：
[0079]
[0080] 其中Γ为阈值，Γ＝λμ，λ取值在2e4～4e4之间，e表示科学计数，取值为10；sign()表示取符号函数，1≤j≤N；
[0081] 将目标元素x1,...,xj,...,xN组合为目标向量xn，
[0082] (4e3)判断收敛条件||xn-xn-1||2<ε是否成立：
[0083] 若成立，则停止计算，重建观测场景散射强度X*＝xn；
[0084] 若不成立，令n＝n+1，更新快速收缩变量yn为：
[0085] yn＝xn-1+(xn-2-xn-1)×(1-t1)/t2，
[0086] 其中 当n＝1时，x0＝0；j表示系数向量xn和xn-1中的第j个元素，xn[j]表示向量xn的第j个元素的值，xn-1[j]表示向量xn-1的第j个元素值；
[0087] 其中t1、t2是两个数值不同的加速因子，t1初始值为1， 将t1更新为t1＝t2，返回步骤(4e1)。
[0088] 步骤五：将重建场景散射强度X*拷贝到主机端内存取模值，并按照先行后列的顺序排列成一个二维矩阵，即得到重建的超声图像。
[0089] 本发明的效果可以通过以下仿真实验说明：
[0090] 1.仿真内容
[0091] 为了验证本发明在CUDA架构下并行计算的优越性，通过一组仿真实验对本发明方法所需的时间与CPU串行所需的时间定量分析，仿真测试平台为：
[0092] CPU：Intel Core i7 4770，内存：16GB
[0093] 显卡：NVIDIA GeForce GTX770，显卡显存：4GB，计算能力3.0
[0094] 软件平台：windows7 64位操作系统，VS2010 x64编译器，CUDA版本：5.5[0095] 2.仿真结果及分析
[0096] 仿真实验用到两组数据，第一组数据观测矩阵大小为1GB，图像大小为128×128，目标点个数为20个；第二组数据观测矩阵大小为8GB，图像大小为256×512，目标点个数为100个。图像的分辨率为100um，要高于传统超声成像分辨率，迭代次数设定为30次，第一组仿真得到的重建图像与原始图像的对比如图8所示。
[0097] 图8表明，CUDA架构下GPU的重建图像与原始图像在主观上评价基本一致，点目标的位置与原始图像完全吻合，从而达到了较高的重建效果。
[0098] 将CUDA架构下GPU并行计算时间和CPU串行时间进行对比，如表1和表2所示：
[0099] 表1 观测矩阵为1GB重建时间
[0100]
[0101] 从表1看出，本发明方法对超声成像的重构速度提高明显，达到了每秒2.5帧图像，基本上解决了压缩感知框架下处理大规模矩阵耗时较长的问题，实现了准实时性。
[0102] 表2 观测矩阵为8GB重建时间
[0103]
[0104] 表2中，GPU的加速效果比表1中的加速比差，其主要原因在于显存不足导致每次迭代都需要进行数据交换，将数据从主机端内存交换到显存端。尽管如此，GPU加速处理后仍然比CPU串行快13倍，而且通过分析数据可以发现在CPU+GPU的方案中数据交换的时间占到总时间的92％，因此，若可以利用提升硬件解决显存不足的问题，加速比仍然可以达到180倍左右。
