用gpu實(shí)現(xiàn)壓縮感知超聲成像的方法
【專利摘要】本發(fā)明公開了一種用GPU實(shí)現(xiàn)壓縮感知超聲成像的方法，主要解決壓縮感知理論框架下成像重建時(shí)間較慢的問題。其實(shí)現(xiàn)步驟為：1.根據(jù)設(shè)定的分辨率，對(duì)探測(cè)區(qū)域進(jìn)行離散化并對(duì)該區(qū)域進(jìn)行寬帶脈沖掃描，得到回波向量和觀測(cè)矩陣，進(jìn)而建立超聲成像數(shù)學(xué)模型；2.將回波向量與觀測(cè)矩陣分塊并復(fù)制到GPU顯存中；3.在GPU中計(jì)算迭代步長；4.將迭代步長帶入快速迭代收縮閾值算法求解出重建觀測(cè)場(chǎng)景散射強(qiáng)度；5.將該散射強(qiáng)度復(fù)制到主存中取模值并排列成一個(gè)二維矩陣，得到重建的超聲圖像。本發(fā)明相對(duì)傳統(tǒng)的快速迭代收縮閾值算法，重建時(shí)間從分鐘級(jí)別降低到毫秒級(jí)別，極大提高了超聲成像的實(shí)時(shí)性，可用于超聲實(shí)時(shí)處理領(lǐng)域。
【專利說明】用GPU實(shí)現(xiàn)壓縮感知超聲成像的方法
【技術(shù)領(lǐng)域】
[0001] 本發(fā)明屬于圖像處理【技術(shù)領(lǐng)域】，特別涉及一種用GPU實(shí)現(xiàn)快速成像的方法，可用 于B超成像。
【背景技術(shù)】
[0002] 醫(yī)學(xué)超聲成像經(jīng)過60余年的發(fā)展，它具有相對(duì)安全、實(shí)時(shí)性好、無創(chuàng)、便攜、價(jià)格 低廉等優(yōu)點(diǎn)，其與X射線診斷技術(shù)、電子計(jì)算機(jī)斷層掃描CT成像技術(shù)、核磁共振成像技術(shù)一 起稱為現(xiàn)代醫(yī)學(xué)四大影像技術(shù)，已令億萬患者受益。
[0003] 但是超聲成像仍然存在一些不足，如分辨率不高，多為毫米級(jí)；受噪聲干擾嚴(yán)重， 圖像質(zhì)量較差；實(shí)時(shí)性一般。
[0004] 近年來，在信號(hào)處理領(lǐng)域興起的壓縮感知理論吸引了諸多學(xué)者的關(guān)注，該理論指 出，只要信號(hào)在某一個(gè)空間W上具有稀疏性，就可以利用觀測(cè)矩陣對(duì)其以遠(yuǎn)低于奈奎斯特 采樣速率進(jìn)行觀測(cè)，進(jìn)而利用優(yōu)化手段高概率地從混疊觀測(cè)中重構(gòu)原信號(hào)，這使得傳感器 的采樣成本大大降低。而且通過恰當(dāng)選擇空間W，信號(hào)的稀疏性越強(qiáng)，精確恢復(fù)原信號(hào)的可 能性就越大，這樣就會(huì)在提高圖像分辨率、抑制噪聲方面有出色的表現(xiàn)。從近幾年國內(nèi)外發(fā) 表的文獻(xiàn)來看，對(duì)壓縮感知理論的研究已經(jīng)涉及眾多領(lǐng)域如壓縮感知CS雷達(dá)成像、醫(yī)學(xué)圖 像處理、光譜分析、遙感圖像處理等，具有非常廣闊的應(yīng)用前景。
[0005] 由于病灶區(qū)域與正常組織的密度特征有明顯差別，可認(rèn)為超聲圖像在空間域是稀 疏的，將壓縮感知理論應(yīng)用到超聲成像可以較好的解決超聲成像分辨率不高和噪聲干擾嚴(yán) 重的問題，但是壓縮感知理論的問題在于重建過程中觀測(cè)矩陣維度巨大、導(dǎo)致計(jì)算復(fù)雜度 高，圖像的重建時(shí)間過長。
[0006] 針對(duì)這個(gè)問題，以色列學(xué)者A.Beck等人在論文"Afastiterative shrinkage-thresholdingalgorithmforlinearinverseproblems，'SIAMJ.IMAGESCI ENCES，Vol. 2.No. 1，pp. 183-202, 2009中提出了快速迭代收縮閾值算法，簡稱FISTA。利用 該算法，超聲成像的基本框架可以表示為：
[0007]
【權(quán)利要求】
1. 一種用GPU實(shí)現(xiàn)壓縮感知超聲成像的方法，其特征在于，包括以下步驟： (1) 將超聲探測(cè)區(qū)域二維離散化，得到N個(gè)離散化的像素點(diǎn)，其中N=TXS，T表示軸 向像素的個(gè)數(shù)，S表示側(cè)向像素的個(gè)數(shù)； (2) 將超聲寬帶脈沖發(fā)射信號(hào)在頻域均勻采樣得到W個(gè)頻點(diǎn)，按頻點(diǎn)順序依次對(duì)已 離散化的二維探測(cè)區(qū)域進(jìn)行一次平面波掃描，每次掃描得到一個(gè)長度為A的局部觀測(cè) 回波向量bt，并將這W個(gè)局部觀測(cè)回波向量按從上到下順序合成長度為M=AXW的觀 測(cè)回波向量b= 同時(shí)保存由這W個(gè)頻點(diǎn)產(chǎn)生的回波聲場(chǎng)強(qiáng)度矩陣 W1,..., Ψ"，其中，A表示超聲線陣的陣元個(gè)數(shù)，矩陣Wt的寬度為A，長度為N， I^t^W； (3) 將回波聲場(chǎng)強(qiáng)度矩陣Ψρ. . .，Ψ,，. . .，Ψ"，按照從上到下順序排列成一個(gè)大小為 MXN的觀測(cè)矩陣Ψ;將離散化的二維探測(cè)區(qū)域按照行優(yōu)先的順序排列成一個(gè)目標(biāo)向量X; (4) 根據(jù)回波向量b和觀測(cè)矩陣Ψ定義基于壓縮感知的超聲成像數(shù)學(xué)模型： 尤* = &^η?η{||ψχ-碌+IlxlJ 1) X 其中礦為重建觀測(cè)場(chǎng)景散射強(qiáng)度，X為目標(biāo)向量，λ為正則化參數(shù)，|Ψλ·-Ζ|表示向量 ψχ-b二范數(shù)的平方，I|x|I1表示目標(biāo)向量X的一范數(shù)； (5) 在GPU中對(duì)上述數(shù)學(xué)模型進(jìn)行求解，得到重建觀測(cè)場(chǎng)景散射強(qiáng)度X% (5a)初始化：n= 0，ε=KT3,其中η表示第η次迭代，ε表示迭代終止條件； (5b)將觀測(cè)矩陣Ψ按行均勻等分成寬度為Μ/Κ，長度為N的K個(gè)子矩陣Ψ= {Ψ31，. . .，Wsi，. . .，Ψ3Κ}，要保證每個(gè)子矩陣數(shù)據(jù)量小于顯存的容量上限，1彡i彡K; (5c)將每個(gè)子矩陣Wsi以列優(yōu)先順序向量化為一維行向量CLMi并采用流操作拷貝 到顯存中，此時(shí)觀測(cè)矩陣Ψ轉(zhuǎn)換成一個(gè)寬度為K，長度為M/KXN的向量化矩陣d_MA= {d_ M1,. . .，CLMi,. . .，d_MAK}，其中觀測(cè)矩陣Ψ與向量化矩陣它們存在--對(duì)應(yīng)的關(guān)系；將 回波向量b等分成長度為Μ/Κ的K塊，b=Ib1,. ..，h. ..，bK}，并拷貝到顯存中； (5d)根據(jù)向量化矩陣d_MA計(jì)算基于梯度下降算法的迭代步長μ: (5dl)在GPU中建立一個(gè)重復(fù)執(zhí)行K次的循環(huán)，每次循環(huán)創(chuàng)建一個(gè)k_st印內(nèi)核函數(shù)， 求出向量(LMAi對(duì)應(yīng)觀測(cè)矩陣Ψ每一列元素模值的和； (5d2)循環(huán)結(jié)束后，求出向量d_sum的二范數(shù)，取倒數(shù)即得到迭代步長μ; (5e)將回波向量b、向量化矩陣d_MA和迭代步長μ帶入快速迭代收縮閾值算法中，經(jīng) 過多次梯度下降和快速閾值收縮過程，直到目標(biāo)向量滿足迭代終止條件，得到重建觀測(cè)場(chǎng) 景散射強(qiáng)度Χ# ; (6) 將重建場(chǎng)景散射強(qiáng)度礦拷貝到主機(jī)端內(nèi)存取模值，并按照先行后列的順序排列成 一個(gè)二維矩陣，得到重建的超聲圖像。
2. 如權(quán)利要求1所述的用GPU實(shí)現(xiàn)壓縮感知超聲成像方法，其特征在于，所述步驟（2) 中，頻點(diǎn)取kt時(shí)的回波聲場(chǎng)強(qiáng)度矩陣Ψ,，通過下式計(jì)算：
其中，Ain(kt)表示超聲寬帶脈沖發(fā)射信號(hào)在頻點(diǎn)取值為kt時(shí)的幅度；表示超聲 寬帶脈沖發(fā)射信號(hào)在頻點(diǎn)取值為kt時(shí)離散二維探測(cè)區(qū)域中各個(gè)像素點(diǎn)返回的相位，A表示 超聲寬帶脈沖發(fā)射信號(hào)方位角的單位向量，指定為軸向方向表示超聲寬帶脈沖發(fā)射信 號(hào)從超聲線陣軸心位置到離散二維探測(cè)區(qū)域各個(gè)像素點(diǎn)距離的向量； 表示格林函數(shù)，表示超聲線陣軸心位置到到超聲線陣各個(gè)陣元距離的向 量，1彡t彡W，1彡j彡N，1彡m彡A。
3. 如權(quán)利要求1所述的用GPU實(shí)現(xiàn)壓縮感知超聲成像方法，其特征在于，所述步驟 (5dl)中，在GPU中建立一個(gè)重復(fù)執(zhí)行K次的循環(huán)，每次循環(huán)執(zhí)行的步驟如下： (5dl.a)在GPU中開辟N/V個(gè)線程塊和長度為N的向量d_sum，每個(gè)線程塊中設(shè)有M/K個(gè)線程，其中向量《^_.、'"/?=[(1_!;111丨11丨，...，(1_51011, /，".，(1_511丨11\]￡；]1\服表示實(shí)數(shù)域，1彡9彡1 I^V^N； (5dl.b)用第z個(gè)線程塊讀取M/KXV個(gè)元素，即向量CLMi中第M/KXVX(z-l)+l到 第M/KXVXz的元素；在該線程塊中分V次取出這些元素，每次按順序取M/K個(gè)元素，依次 分配給線程1?線程M/K，每個(gè)線程計(jì)算各自元素的模值；將這M/K個(gè)模值累加給向量d_ sum的第d_sumvx(z_lHp個(gè)元素,其中1彡z彡N/V，1彡p彡V。
4. 如權(quán)利要求1或3所述的用GPU實(shí)現(xiàn)壓縮感知超聲成像方法，其特征在于，所述步驟 (5dl)中，常量V的取值要滿足V對(duì)N取余值為0,優(yōu)選的，V取值為2?4可以達(dá)到理想的 性能。
5. 如權(quán)利要求1所述的用GPU實(shí)現(xiàn)壓縮感知超聲成像方法，其特征在于，所述步驟 (5e)，按如下步驟進(jìn)行： (5el)將回波向量b、向量化矩陣d_MA和迭代步長μ帶入下式，更新梯度下降序列un:un=Yn-IiX cublas (d_MAH · (cublas(d_MA · yn)~b)) 3) 其中yn是快速收縮向量，初始值為〇,長度為N ;d_MAH表示向量化矩陣d_MA的共軛轉(zhuǎn) 置；μ為迭代步長；cublas是CUDA線性代數(shù)庫CUBLAS Library中所支持的庫函數(shù)，它可以 執(zhí)行向量化矩陣與向量的乘法運(yùn)算； (5e2)將梯度下降序列Un帶入下式，更新目標(biāo)場(chǎng)景散射強(qiáng)度Xn: xn =Sr(Un) 4) 其中Sr為閾值函數(shù)：
其中Γ為閾值，Γ=λμ，λ取值在2e4?4e4之間，e表示科學(xué)計(jì)數(shù)，取值為10 ; signO表不取符號(hào)函數(shù)； (5e3)判斷收斂條件I|xn-xn_」12〈ε是否成立： 若成立，則停止計(jì)算，重建觀測(cè)場(chǎng)景散射強(qiáng)度礦=Xn ; 若不成立，令η=η+1，更新快速收縮向量yn為： Yn = Xn-l+(Xn-2_Xn-i)X (U1Vt2 5) 其中 |x" =^X(x"[j]-x" ,[)'])'，當(dāng)η= 1 時(shí)，xQ = 0;j表示系數(shù)向量xn 和Xlri 中 的第j個(gè)元素，xn[j]表示向量Xn的第j個(gè)元素的值，XlriU]表示向量Xlri的第j個(gè)元素 值；tpt2是兩個(gè)數(shù)值不同的加速因子，L的初始值為1，6=0.5 + 0.5x0+4x6x00 新為h=t2，返回步驟（5el)。
6.如權(quán)利要求5所述的用GPU實(shí)現(xiàn)壓縮感知超聲成像方法，其特征在于，所述步驟 (5e2)中，更新目標(biāo)向量xn的計(jì)算是在內(nèi)核函數(shù)kernel中完成的，即xn =Sr^(Un)通過使用 多線程進(jìn)行并行計(jì)算得到。
【文檔編號(hào)】A61B8/00GK104306022SQ201410578176
【公開日】2015年1月28日 申請(qǐng)日期:2014年10月24日 優(yōu)先權(quán)日:2014年10月24日
【發(fā)明者】林杰, 韓亭玉, 賀玉高, 石光明 申請(qǐng)人:西安電子科技大學(xué)