醫(yī)學(xué)超聲實(shí)時(shí)成像中周期性偽影的抑制①

程?hào)|耀,尹 皓,劉東權(quán)

(四川大學(xué) 計(jì)算機(jī)學(xué)院,成都 610065)

隨著醫(yī)療科技不斷的發(fā)展和進(jìn)步,超聲圖像診斷是一項(xiàng)重要的醫(yī)學(xué)圖像診斷手段,由于超聲成像的圖像質(zhì)量較差,需要采用相對(duì)應(yīng)的方法來(lái)對(duì)超聲圖像進(jìn)行增強(qiáng)[1],當(dāng)前,超聲成像技術(shù)在臨床醫(yī)學(xué)診斷方面應(yīng)用的越來(lái)越廣泛,超聲成像技術(shù)作為一種非侵入性,對(duì)人體無(wú)害的病理檢測(cè)技術(shù)[2],一直受到國(guó)內(nèi)外相關(guān)科學(xué)研究人員的關(guān)注.超聲診斷通過(guò)超聲儀向人體發(fā)射超聲波,并通過(guò)不斷的掃描來(lái)對(duì)不同的人體組織進(jìn)行探測(cè),當(dāng)探頭接收到由人體組織和器官反射回來(lái)的超聲信號(hào)后,再通過(guò)放大信號(hào)以及一系列信息的處理,最終形成人體組織器官的超聲實(shí)時(shí)成像,三維成像可以提供直觀的立體圖像信息,在心肌損傷定位,胸腹部腫瘤的檢測(cè)等方面有著重要的作用[3].由于超聲探頭的材質(zhì)例如超聲陶瓷換能器制造的局限性以及為了降低成本和提高幀速率,而采用的低通道掃描方案,往往會(huì)使原始的超聲成像包含束狀條紋偽影,這些偽影往往會(huì)掩蓋人體的組織信息,從而引起醫(yī)生錯(cuò)誤的診斷.常規(guī)的偽影抑制方法都是在空間域來(lái)進(jìn)行的,對(duì)于周期性偽影難以實(shí)現(xiàn)最佳的去除效果.目前,隨著醫(yī)學(xué)超聲圖像處理技術(shù)不斷地成熟,具體的研究方向包括降噪,分割,壓縮等圖像處理技術(shù),還包括超聲圖像采集技術(shù)[2].本研究提出了一種基于圖像的頻率域的偽影抑制算法,根據(jù)大量的實(shí)驗(yàn)結(jié)果進(jìn)行定性的分析及定量的MSE和SNR 來(lái)驗(yàn)證該算法的準(zhǔn)確性,其結(jié)果對(duì)超聲成像的周期性偽影有著顯著的抑制效果,進(jìn)而幫助醫(yī)生得到更清晰的超聲成像,以提高醫(yī)生的診斷效率和效果.

1 超聲圖像通用偽影抑制介紹

在傳統(tǒng)的超聲圖像處理領(lǐng)域,圖像噪聲抑制是數(shù)字圖像處理中的重要環(huán)節(jié)和步驟,去噪效果的好壞直接影響后續(xù)圖像處理的相關(guān)工作.通常,抑制偽影的算法分為兩大類,第一類是基于空間域?yàn)V波的圖像噪聲抑制算法,第二類是基于頻率域?yàn)V波的圖像噪聲抑制算法,本研究將列出一些常用的空間域?yàn)V波去噪算法以及接下來(lái)將使用的頻率域偽影抑制算法的設(shè)計(jì)思想,并通過(guò)對(duì)比空間域偽影抑制算法和頻率域偽影抑制算法的優(yōu)缺點(diǎn)及采用空間域?yàn)V波算法的局限性,詳細(xì)闡述本研究使用基于頻率域的超聲實(shí)時(shí)成像的周期性偽影抑制算法的原因.

1.1 空間域偽影抑制

圖像空間域偽影抑制常常使用的是平滑處理,平滑主要利用卷積運(yùn)算對(duì)圖像鄰域的像素灰度進(jìn)行平均化,從而減少圖像中的噪點(diǎn)[4].圖像的平滑處理分為線性平滑和非線性平滑,線性平滑包括均值濾波和高斯濾波,非線性平滑包括中值濾波和雙邊濾波.目前隨著空間域?yàn)V波技術(shù)不斷地發(fā)展,越來(lái)越多的改進(jìn)濾波技術(shù)不斷的應(yīng)用于圖像平滑之中,接下來(lái)介紹并分析目前改進(jìn)的典型濾波算法.

以迭代均值濾波算法[5]為例,它改進(jìn)了傳統(tǒng)的均值濾波算法,它利用鄰域內(nèi)已經(jīng)估算的重構(gòu)灰度級(jí)參與當(dāng)前像素灰度級(jí)的運(yùn)算,提高了圖像重構(gòu)的質(zhì)量和效率.雖然迭代均值濾波算法能夠有效地抑制噪聲,但是也會(huì)破壞圖像的細(xì)節(jié)部分,因此對(duì)于具有強(qiáng)邊界信息的圖像,此算法會(huì)破壞圖像的邊緣信息,這就導(dǎo)致圖像整體變得模糊.

文獻(xiàn)[6]提出了一種自適應(yīng)各項(xiàng)異性高斯濾波算法,該算法首先計(jì)算圖像中各個(gè)像素位置的梯度值以及梯度方向,根據(jù)梯度值的大小來(lái)決定采用的濾波類型,采用各項(xiàng)異性高斯濾波,通過(guò)該點(diǎn)的梯度方向角得到各項(xiàng)異性高斯算子的長(zhǎng)軸方向,此算法在邊緣信息的維護(hù)方面強(qiáng)于傳統(tǒng)的高斯濾波算法.此算法的局限性在于只能夠?qū)Ω咚乖肼曔M(jìn)行平滑,對(duì)于在空間域上周期性分布的噪聲,該算法就失去了作用效果.

文獻(xiàn)[7]提出了一種可以保留圖像細(xì)節(jié)信息的中值濾波改進(jìn)算法,通過(guò)檢查待處理像素點(diǎn)是否存在于角點(diǎn)以及線段中,并根據(jù)結(jié)果來(lái)決定是否進(jìn)行保留或者采用中值濾波算法進(jìn)行處理,該算法可以保留圖像中的細(xì)節(jié)信息同時(shí)可以避免圖像的模糊.該算法的局限性在于當(dāng)噪聲像素個(gè)數(shù)大于窗口內(nèi)像素個(gè)數(shù)的一半時(shí),會(huì)導(dǎo)致排序后的中間值仍為噪聲,因此去噪效果非常差.

文獻(xiàn)[8]提出了一種改進(jìn)的雙邊濾波和閾值函數(shù)的圖像增強(qiáng)算法,該算法對(duì)帶有噪聲的圖像進(jìn)行小波分解,從而獲取圖像的低頻和高頻系數(shù),采用改進(jìn)的雙邊濾波的Retinex 算法對(duì)圖像的低頻系數(shù)進(jìn)行處理,采用閾值函數(shù)對(duì)高頻系數(shù)進(jìn)行處理,最后進(jìn)行小波反變換重構(gòu)來(lái)實(shí)現(xiàn)圖像中噪聲的消除以及圖像細(xì)節(jié)部分的增強(qiáng).但對(duì)于周期性偽影,由于其在小波域上均存在于高頻分量以及低頻分量中,因此難以實(shí)現(xiàn)對(duì)其進(jìn)行分離從而達(dá)到抑制的目的.

1.2 頻率域偽影抑制

上面主要介紹了4 種基于空間域偽影抑制的濾波算法,但它們?cè)趯?duì)特定分布的偽影在進(jìn)行抑制的同時(shí)往往會(huì)具有局限性,例如周期性偽影,在空間域上采用平滑操作,雖然會(huì)在一定程度上降低偽影的灰度級(jí),但是通過(guò)平滑后的圖像會(huì)變得模糊,圖像偽影附近的有效信息會(huì)產(chǎn)生丟失,為了解決這一問(wèn)題,在CT 領(lǐng)域,常常需要獲取紋理在頻率域上的相位分布和幅值分布,來(lái)達(dá)到偽影抑制的目的[9].對(duì)于帶有周期性偽影的超聲圖像需要通過(guò)變換域處理的方式來(lái)對(duì)其進(jìn)行抑制,變換域處理主要指的是圖像的頻率域變換及小波變換,通常情況下,頻率域具有和空間域的互補(bǔ)性,例如周期性偽影的抑制在空間域處理往往需要對(duì)圖像進(jìn)行分割以及融合,而這些偽影在頻率域就顯得具有一定規(guī)律,這些特定在頻率域表現(xiàn)得更為直觀,它可以解釋空間域?yàn)V波的一些屬性,因此只需要找出周期性偽影在頻率域上的分布并對(duì)其進(jìn)行衰減[10]就可以達(dá)到去除的目的,本課題將對(duì)于帶有周期性偽影的超聲圖像來(lái)進(jìn)行基于頻率域上[11]算法研究.

2 正弦波模擬偽影抑制

這一節(jié)主要使用正弦波來(lái)模擬周期性偽影,并通過(guò)正弦波的疊加實(shí)驗(yàn)來(lái)判斷周期性偽影在頻域的表現(xiàn)特征.

2.1 二維傅立葉變換

本研究對(duì)周期性偽影的抑制在超聲圖像的頻率域來(lái)進(jìn)行處理主要是通過(guò)圖像的二維傅立葉變換來(lái)完成的,它是由20世紀(jì)60年代的Cooley 和Tukey 提出的一種算法[12],能夠明顯提高運(yùn)算效率.二維快速傅立葉變換(FFT)是圖像處理領(lǐng)域中常用的分析和計(jì)算工具,常常應(yīng)用于頻譜分析,頻域?yàn)V波和卷積計(jì)算等相關(guān)圖像分析當(dāng)中[13].二維傅立葉變換可以看出兩個(gè)一維傅立葉變換的結(jié)合,一維傅立葉變換的作用是將一個(gè)一維信號(hào)分解成若干個(gè)復(fù)指數(shù)波,即有:

再根據(jù)歐拉方程:

對(duì)于一個(gè)正弦波,想要確定它就要定義3 個(gè)參數(shù),分別是頻率w,幅值A(chǔ),相位φ,因此在頻率域中,一維坐標(biāo)代表頻率,每個(gè)坐標(biāo)對(duì)應(yīng)的函數(shù)值為F(u)的一個(gè)復(fù)數(shù),它的幅度值|F(u)|.一維傅立葉變換是一種基變換,在空間域中,基是一族脈沖信號(hào)δ(x-n),在頻率域中,基是e-j2πu,這組基是正交基.一維信號(hào)可以看成一個(gè)序列,傅立葉變換能將其分解成若干個(gè)簡(jiǎn)單函數(shù)之和.對(duì)于圖像來(lái)說(shuō),可以將它看成一個(gè)二維信號(hào),一個(gè)二維信號(hào)可以分解為若干個(gè)復(fù)平面波之和.對(duì)于一個(gè)長(zhǎng)度為M,寬度為N的圖像f(x,y),其二維離散傅立葉變換由式(3)決定:

類似于一維中的情景,要完成表達(dá)式的計(jì)算,需要對(duì)u,v進(jìn)行計(jì)算,其中對(duì)應(yīng)的二維傅立葉反變換表達(dá)式有:

通過(guò)式(4),可以計(jì)算出每一個(gè)平面波在圖像中的成分.二維傅立葉變換實(shí)質(zhì)是將圖像與每個(gè)不同頻率的不同方向的復(fù)平面波做內(nèi)積運(yùn)算,也就是求在基向量e-j2π(ux+vy)上投影的過(guò)程.根據(jù)式(4)可以得到傅立葉變換的一些特性: 平移性,旋轉(zhuǎn)不變性,卷積性[14].根據(jù)式(4)可以推斷出在頻域中原點(diǎn)平移到(u0,v0)時(shí),其空間域要乘以一個(gè)正的指數(shù)項(xiàng)ej2π(u0x/M+v0y/N),空間域中原圖像原點(diǎn)移動(dòng)到(x0,y0)時(shí),其頻率域F(u,v)要乘以一個(gè)負(fù)的指數(shù)項(xiàng)e-j2π(u0x/M+v0y/N),在進(jìn)行圖像的傅立葉變換時(shí),往往需要將F(u,v)原點(diǎn)平移到頻域中心,u0=M/2,v0=N/2,代入公式中得到將圖像的頻域原點(diǎn)平移到中心就需要在f(x,y)上乘以(-1)(x+y)即可實(shí)現(xiàn),體現(xiàn)出二維FFT 變換結(jié)果均具有中心共軛對(duì)稱性[14].二維傅立葉變換的旋轉(zhuǎn)不變性指的是如果f(x,y)旋轉(zhuǎn)了一個(gè)角度,那么f(x,y)旋轉(zhuǎn)后的圖像的傅立葉變換也會(huì)旋轉(zhuǎn)相同的角度,用極坐標(biāo)代替平面直角坐標(biāo)系即有:

替換法則有f(x,y)→f(r,θ)→F(w,φ),如果f(x,y)被旋轉(zhuǎn)同一角度,那么有f(r,θ+θ0)→F(w,w+θ0)[14]得出結(jié)論一幅圖像的空間域上的旋轉(zhuǎn)角度和頻率域上的旋轉(zhuǎn)角度是一致的.同時(shí)二維傅立葉變換還有卷積特性,即空間域上的卷積運(yùn)算等于頻率域上的乘積運(yùn)算,該性質(zhì)的好處在于將需要進(jìn)行復(fù)雜的卷積運(yùn)算轉(zhuǎn)化為乘積運(yùn)算,極大地提高了計(jì)算效率,這也是本研究采用頻率域抑制偽影的原因.

2.2 正弦波的模擬實(shí)驗(yàn)

根據(jù)上述的結(jié)論,為了更好地達(dá)到頻率域上周期性偽影抑制的效果,本研究模擬了兩個(gè)正弦波的疊加后然后通過(guò)二維傅立葉變換后的頻譜圖,再通過(guò)上述的特性來(lái)抑制特定頻率的正弦波.設(shè)計(jì)的步驟如實(shí)驗(yàn)1.

實(shí)驗(yàn)1.正弦波模擬偽影實(shí)驗(yàn)1)編程分別實(shí)現(xiàn)一幅256×256,頻率為1/18 的水平正弦波和垂直正弦波;2)將水平正弦波和垂直正弦波進(jìn)行疊加,得到混合正弦波;3)分別對(duì)水平,垂直正弦波進(jìn)行二維傅立葉變換,觀察各自得到的頻率域譜.4)將混合正弦波做二維傅立葉變換,觀察得到的頻譜,并比較之前獲取的傅立葉頻譜;5)找出頻譜上水平方向的亮點(diǎn)并置為0,這些亮點(diǎn)根據(jù)傅立葉變換的特性可以判斷出屬于垂直正弦波;6)對(duì)頻譜進(jìn)行逆傅立葉變換到空間域,觀察到最終結(jié)果.

上述設(shè)計(jì)步驟通過(guò)豎直正弦波模擬周期性偽影并進(jìn)行傅立葉變換和逆傅立葉變換,來(lái)驗(yàn)證在頻率域抑制周期性偽影的效果,如圖1、圖2.上述兩個(gè)正弦波分別代表頻率為1/18 的垂直于水平正弦波,按照步驟進(jìn)行疊加得到混合正弦波,這里采用OpenCV 提供的add 接口來(lái)實(shí)現(xiàn)正弦波的疊加,如圖3.分別對(duì)垂直和水平正弦波及混合正弦波進(jìn)行二維傅立葉變換,得到垂直和水平正弦波頻譜,如圖4、圖5.通過(guò)以上觀察發(fā)現(xiàn),當(dāng)豎直正弦波旋轉(zhuǎn)90 度到水平位置后,其頻域也會(huì)旋轉(zhuǎn)90 度,這驗(yàn)證了空間域和頻域變換之間對(duì)應(yīng)旋轉(zhuǎn)不變性的特征,圖6 是混合正弦波的頻譜.下一步通過(guò)將混合正弦波水平中線上的數(shù)值置為0,再通過(guò)逆向傅立葉變換觀察復(fù)原后的正弦波,如圖7、圖8.

圖1 頻率為1/18 的垂直正弦波

圖2 頻率為1/18 的水平正弦波

圖3 混合正弦波

圖4 垂直正弦波頻譜

圖5 水平正弦波頻譜

圖6 混合正弦波頻譜

圖7 去除水平方向上脈沖點(diǎn)后的頻譜

圖8 復(fù)原后的空間域正弦波

3 超聲實(shí)時(shí)成像中周期性偽影抑制實(shí)驗(yàn)

本研究將對(duì)超聲實(shí)時(shí)成像進(jìn)行分幀處理,對(duì)每一幀的超聲圖像進(jìn)行2D-FFT 到其頻率域上,并通過(guò)設(shè)計(jì)算法達(dá)到抑制周期性偽影抑制的效果,再對(duì)所有的超聲圖像進(jìn)行合幀處理,完成超聲視頻的偽影抑制.

3.1 超聲實(shí)時(shí)成像的幀率及OpenCV 分幀合幀

超聲實(shí)時(shí)就是一系列捕獲的超聲圖像以給定的頻率顯示,通過(guò)在一序列的特定幀處停止可獲得單個(gè)視頻幀,也就是圖像,對(duì)于超聲實(shí)時(shí)成像中帶有的周期性偽影,本研究采用的方法為使用OpenCV 這一計(jì)算機(jī)視覺(jué)開(kāi)源庫(kù)[15]來(lái)對(duì)成像進(jìn)行離線分幀處理,OpenCV提供了視頻幀的讀取,顯示函數(shù)以及獲得視頻幀率的函數(shù)[16],幀率代表以幀為單位的圖像連續(xù)出現(xiàn)在顯示器上的頻率以及單位時(shí)間內(nèi)通過(guò)圖像的數(shù)量,對(duì)超聲實(shí)時(shí)成像中的偽影進(jìn)行處理首先需要對(duì)單個(gè)成像幀中的圖像進(jìn)行采樣,采樣完成后對(duì)單個(gè)圖像進(jìn)行偽影抑制,才能達(dá)到理想的效果.

OpenCV 分幀的主要步驟為首先使用OpenCV 提供的接口VideoCapture 來(lái)打開(kāi)一個(gè)視頻,獲取視頻的相關(guān)信息,輸出視頻的高度和寬度以及幀率,視頻的每一幀寬度和高度都是一致的,根據(jù)上述信息來(lái)獲取單個(gè)超聲圖像.合幀的主要步驟為創(chuàng)建VideoWriter 對(duì)象,用來(lái)讀取圖像并輸出到指定文件中,包括設(shè)置參數(shù)例如幀數(shù)以及圖像尺寸大小,從而獲取到處理后的超聲成像.圖9 是使用單個(gè)幀內(nèi)的部分超聲圖像作為樣本.

圖9 是人體某組織通過(guò)超聲探頭掃描的通道數(shù)據(jù)的視頻通過(guò)分幀后的圖像呈現(xiàn),可以發(fā)現(xiàn)這些圖像普遍具有周期性的掃描線,這些是由超聲探頭采用的低通道掃描導(dǎo)致的,同時(shí),這些偽影會(huì)遮擋該人體組織中的有用信息,會(huì)對(duì)醫(yī)生的診斷造成干擾.本研究將采取在頻率域?qū)@些周期性掃描線進(jìn)行抑制,這樣也會(huì)提高偽影抑制的效率.

圖9 帶有周期性偽影的人體組織(低通道掃描)

3.2 偽影脈沖檢測(cè)算法及結(jié)果的定量分析

在大部分圖像中,鄰近像素高度相關(guān),距離較遠(yuǎn)的像素相關(guān)性較弱,在頻域上也存在相關(guān)特性[17],根據(jù)之前做的正弦波模擬偽影實(shí)驗(yàn),可以發(fā)現(xiàn)具有周期性的正弦波在頻率域上的方向和空間域的方向垂直,對(duì)于豎直周期性偽影,其在頻率域上的方向是在水平方向上,通常表現(xiàn)為具有對(duì)稱性的峰值[18],為了更加準(zhǔn)確的定位周期性偽影在頻率域上的位置,需要設(shè)計(jì)一個(gè)算法模型,該算法模型能準(zhǔn)確地找出偽影在頻率域上的坐標(biāo).可以采用滑動(dòng)窗口掃描結(jié)合合適的閾值大小這一算法模型用來(lái)選取疑似偽影脈沖頻率域上的點(diǎn)并進(jìn)行幅值降低處理.

該算法模型的實(shí)現(xiàn)過(guò)程如算法1.

作為典型的東北城市，長(zhǎng)春應(yīng)如何依據(jù)入室盜竊案件發(fā)案特點(diǎn)和規(guī)律做好防范工作，將是其能否躋身全國(guó)安全城市的重要指標(biāo)。具體來(lái)說(shuō)公安機(jī)關(guān)應(yīng)從以下幾個(gè)方面做出努力。

算法1.滑動(dòng)窗口閾值掃描算法模型1)定義一個(gè)掃描路徑,被掃描點(diǎn)到頻域中心的距離由大到小來(lái)進(jìn)行搜索;2)定義一個(gè)大小為3×3 的滑動(dòng)窗口,依次從左往右,從上往下開(kāi)始滑動(dòng);3)找出滑動(dòng)窗口內(nèi)中心點(diǎn)的幅值大于當(dāng)前點(diǎn)的上下左右各點(diǎn)幅度值的所有點(diǎn).4)定義一個(gè)閾值T,閾值大小由超聲圖像的動(dòng)態(tài)范圍決定,例如動(dòng)態(tài)范圍為60 db,閾值可以設(shè)置為圖像動(dòng)態(tài)范圍的k 倍(0≤k ＜ 1),k 由超聲圖像偽影抑制結(jié)果評(píng)估模型決定,對(duì)系數(shù)-預(yù)期結(jié)果做采樣,找出最佳閾值系數(shù)大致范圍,然后設(shè)置截止半徑R,截止半徑由圖像頻域幅度峰值到左右兩側(cè)下一個(gè)峰值的距離確定,用來(lái)保證頻域中心點(diǎn)的數(shù)據(jù)不被破壞,R 的范圍由采樣決定.如果該點(diǎn)所在的滑動(dòng)窗口內(nèi)的動(dòng)態(tài)范圍大于T,那么將該點(diǎn)保存至集合中.5)將滿足步驟4 的點(diǎn)存入集合P,P 中的點(diǎn)就是疑似周期性偽影在頻譜上的坐標(biāo)點(diǎn),分別是P1,P2,P3,…, Pn,并且r(P1)＞ r(P2)＞ r(P3)＞…＞ r(Pn),r(P)為P 集合中各個(gè)點(diǎn)到頻譜中心點(diǎn)的距離.6)對(duì)P 中的點(diǎn)進(jìn)行降低幅值處理,壓低后的幅值大小為小于等于通過(guò)閾值計(jì)算后的大小.

通過(guò)該算法模型,可以找出頻域上疑似周期性偽影的點(diǎn)集,然后就可以通過(guò)對(duì)找到的點(diǎn)集進(jìn)行幅值降低的處理,再通過(guò)逆向傅立葉變換得到抑制周期性偽影的超聲圖像.然后根據(jù)客觀準(zhǔn)則來(lái)定量的計(jì)算超聲圖像的偽影抑制效果,定義兩個(gè)集合,分別用來(lái)標(biāo)記偽影區(qū)域A 和非偽影區(qū)域B,假設(shè)O為無(wú)偽影超聲圖像,G為原始含偽影圖像,其大小為M×N,Gn為偽影抑制后的圖像,需要統(tǒng)計(jì)G中偽影區(qū)域A及非偽影區(qū)域B和Gn中原偽影區(qū)域An及非偽影區(qū)域Bn然后通過(guò)計(jì)算區(qū)域之間偏離程度的方法來(lái)分析圖像的偽影抑制質(zhì)量,通常采用的方法是均方誤差估計(jì)算法,它通過(guò)計(jì)算輸入輸出圖像的均方值MSE(mean square error)來(lái)判定[19].對(duì)于該圖像中的任意一點(diǎn)(x,y),它的誤差值表達(dá)式有:

均方誤差可以表示為:

均方誤差值越小,其偽影抑制效果越好.同時(shí)也可以采用信噪比SNR(signal-to-noise ratio)[19]作為另一種對(duì)偽影進(jìn)行抑制所產(chǎn)生效果的評(píng)價(jià)準(zhǔn)則,同時(shí)也是使用最廣泛的客觀測(cè)量法之一,即有:

根據(jù)式(8),對(duì)于偽影區(qū)域A,假設(shè)其區(qū)域大小為S1,則非偽影區(qū)域B大小為S2,則其均方誤差為:

其信噪比為:

對(duì)于偽影區(qū)域,其處理前后的均方誤差應(yīng)盡可能大,對(duì)于非偽影區(qū),其均方誤差應(yīng)盡可能小.接下來(lái),分別通過(guò)以上兩種評(píng)估模型MSE和SNR可以來(lái)驗(yàn)證這個(gè)算法模型的準(zhǔn)確性.

3.3 算法模型的驗(yàn)證和實(shí)現(xiàn)

根據(jù)上述算法模型,該算法實(shí)現(xiàn)的功能在于輸入一幅帶有周期性偽影的(M×N)超聲圖像,通過(guò)2DFFT 后變換到頻率域,采用第3.2 節(jié)中設(shè)計(jì)的算法,再使用2D-IFFT 反變換到空間域,得出去除周期性偽影的超聲圖像,經(jīng)過(guò)本研究提供的算法模型測(cè)試,通過(guò)處理后的超聲圖像偽影得到去除,算法作用效果明顯.整體周期性偽影抑制模型如圖10.

圖10 整體周期性偽影抑制模型

根據(jù)算法實(shí)現(xiàn)功能模塊來(lái)對(duì)該算法模型進(jìn)行驗(yàn)證.使用OpenCV 編程實(shí)現(xiàn)該算法模型思想,首先在一幅超聲圖像中(大小為M×N),分別對(duì)其偽影區(qū)域A和非偽影區(qū)域B進(jìn)行標(biāo)記,作為后續(xù)定量分析的參考,如圖11、圖12.

圖11 采集到的低通道超聲圖像

在圖12 中左側(cè)方框?yàn)閭斡皡^(qū),右側(cè)方框?yàn)榉莻斡皡^(qū).獲得該頻譜后,再由算法模型確定包含這些點(diǎn)的坐標(biāo),為了驗(yàn)證這些點(diǎn)坐標(biāo)的準(zhǔn)確性,本研究需要獲取到位于水平中心線上所有點(diǎn)的幅值曲線函數(shù),即當(dāng)x=M/2 時(shí),來(lái)判斷該點(diǎn)的大致范圍是否坐落于該曲線上,如圖13.

圖12 標(biāo)記示意圖(左側(cè)方框?yàn)閭斡皡^(qū),右側(cè)方框?yàn)榉莻斡皡^(qū))

圖13 水平方向上頻域的幅度分布函數(shù)

通過(guò)獲取到的幅值曲線函數(shù),可以發(fā)現(xiàn),位于中心點(diǎn)的幅值最高,其左右兩側(cè)也有部分峰值,可以確定該峰值是周期性偽影在頻域上的表現(xiàn)形式,將獲取的點(diǎn)集的坐標(biāo)來(lái)和峰值進(jìn)行對(duì)比,能夠確定通過(guò)該算法模型獲取頻域上疑似的偽影坐標(biāo)是正確的,這些點(diǎn)存入集合P中,d為集合P中距離頻域中心的最短距離,同時(shí)需要保留頻譜中心點(diǎn)上的有用信息,避免由于該算法誤判導(dǎo)致圖像的非偽影區(qū)紋理細(xì)節(jié)被破壞,假設(shè)頻域某點(diǎn)坐標(biāo)為(x0,y0)則其到中心點(diǎn)(M/2,N/2)的歐氏距離為:

對(duì)于一幅超聲圖像,其頻率域中心點(diǎn)范圍代表的是圖像的低頻信號(hào),其四周范圍代表的是圖像的高頻信息,根據(jù)R的范圍,來(lái)對(duì)R進(jìn)行采樣,分別計(jì)算原始超聲圖像和逆傅立葉變換之后的偽影區(qū)和非偽影區(qū)的MSE和SNR值來(lái)計(jì)算出最優(yōu)截止半徑,并觀察它們之間的相關(guān)性,偽影抑制效果以及圖像中有用信息是否完整.確定R值后,下一步就是對(duì)P中的點(diǎn)根據(jù)圖像的動(dòng)態(tài)范圍進(jìn)行幅度值上的降低.動(dòng)態(tài)范圍是一個(gè)用于定義在一定的范圍內(nèi)能夠捕捉圖像細(xì)節(jié)的術(shù)語(yǔ),通常用一幅圖像的最低值到最高溢出值之間的范圍來(lái)描述,動(dòng)態(tài)范圍越大,就能盡可能地保留圖像高光區(qū)和陰影區(qū)的信息.所以通過(guò)改變集合P中點(diǎn)的幅度值來(lái)達(dá)到符合預(yù)期的動(dòng)態(tài)范圍.假設(shè)一幅圖像G,它的頻率域幅值的最大值為max(G),最小值為min(G),它的全局動(dòng)態(tài)范圍計(jì)算公式為:

計(jì)算出動(dòng)態(tài)范圍,根據(jù)給定的閾值系數(shù)K(0≤K＜1),根據(jù)P集合確定滑動(dòng)窗口內(nèi)所有點(diǎn)的幅度值,重新計(jì)算其中心點(diǎn)動(dòng)態(tài)范圍,用來(lái)確定相應(yīng)的動(dòng)態(tài)范圍大小為:

再根據(jù)式(14)得出新的幅度值為:

將P集合中的點(diǎn)幅度值設(shè)置成符合當(dāng)前滑動(dòng)窗口內(nèi)的數(shù)值A(chǔ),通過(guò)修改幅度值來(lái)實(shí)現(xiàn)了周期性偽影的抑制效果,可以通過(guò)逆向傅立葉變換觀察算法的使用效果.為了確定閾值系數(shù)K的大小,需要對(duì)其范圍進(jìn)行采樣,在本次實(shí)驗(yàn)驗(yàn)證中,采樣周期為0.1,并從視頻中獲取1 張超聲圖像,依次對(duì)不同采樣點(diǎn)區(qū)域A(偽影區(qū)域)區(qū)域B(非偽影區(qū)域)的MSE和SNR作為偽影抑制標(biāo)準(zhǔn),找出最佳閾值系數(shù).確定閾值系數(shù)后,使用上面的單幀部分超聲圖像來(lái)進(jìn)行驗(yàn)證偽影抑制算法模型的準(zhǔn)確性.通過(guò)該算法的程序并選取4 個(gè)在閾值系數(shù)為0 的情況下不同截止半徑以及3 個(gè)在截止半徑為d時(shí)不同閾值系數(shù)的采樣點(diǎn)并進(jìn)行對(duì)比從而確定最佳截止半徑及閾值系數(shù),如圖14-圖20.左側(cè)為原始超聲圖像,中間為通過(guò)該算法改進(jìn)后的頻率域圖像,右側(cè)為通過(guò)頻率域圖像進(jìn)行逆傅立葉變換后的圖像.

圖14 截止半徑為0 的超聲低通道圖像效果

圖15 截止半徑為2 的超聲低通道圖像效果

圖16 截止半徑為6 的超聲低通道圖像效果

圖17 截止半徑為20 的超聲低通道圖像效果

圖18 閾值系數(shù)為0.9 的偽影抑制結(jié)果

圖19 閾值系數(shù)為0.6 的偽影抑制結(jié)果

圖20 閾值系數(shù)為0.1 的偽影抑制結(jié)果

通過(guò)上述人體組織的偽影抑制仿真實(shí)驗(yàn),發(fā)現(xiàn)閾值系數(shù)K以及截止半徑R對(duì)超聲圖像的偽影區(qū)域A及非偽影區(qū)域B的MSE和SNR具有相關(guān)性,所以在該算法中,需要選取合適的截止半徑R以及閾值系數(shù)K來(lái)實(shí)現(xiàn)超聲圖像的偽影抑制效果.

在超聲實(shí)時(shí)成像領(lǐng)域,對(duì)于超聲圖像處理效果上的判斷通過(guò)定性分析以及定量分析兩者結(jié)合來(lái)完成的,本研究在表1 以及表2 中對(duì)于偽影的抑制效果和有效信息的完整程度分別給出了3 個(gè)評(píng)價(jià)標(biāo)準(zhǔn).偽影抑制效果評(píng)價(jià)分為3 個(gè)等級(jí),分別是很差,一般,良好.非偽影區(qū)域有效信息完整程度評(píng)價(jià)分為3 個(gè)等級(jí),分別是無(wú)效,一般,良好.定性的評(píng)價(jià)標(biāo)準(zhǔn)往往由超聲醫(yī)師的經(jīng)驗(yàn)來(lái)決定.本研究在前文中提到了將超聲圖像分為偽影區(qū)域以及非偽影區(qū)域,這是在超聲醫(yī)師的幫助下完成的.本文同時(shí)也給出了定量的分析,對(duì)于超聲醫(yī)師勾勒出的偽影區(qū)域以及非偽影區(qū)域,本文分別計(jì)算了在采用該算法之后兩個(gè)區(qū)域內(nèi)的MSE數(shù)值和SNR數(shù)值,通過(guò)對(duì)這些數(shù)據(jù)進(jìn)行分析結(jié)合超聲醫(yī)師的經(jīng)驗(yàn)完成超聲圖像質(zhì)量上的評(píng)估.

截止半徑R對(duì)該算法在超聲圖像上的作用起著重要的作用,當(dāng)截止半徑過(guò)小時(shí),頻域中心的數(shù)據(jù)會(huì)被破壞,不僅不能夠?qū)Τ晥D像進(jìn)行偽影抑制,而且還會(huì)破壞圖像中有效的醫(yī)學(xué)信息,結(jié)合超聲醫(yī)師對(duì)超聲圖像判斷的經(jīng)驗(yàn),有效信息的完整程度的評(píng)價(jià)等級(jí)為無(wú)效.當(dāng)隨著截止半徑不斷擴(kuò)大,部分頻域中心數(shù)據(jù)得以保留,但仍有部分?jǐn)?shù)據(jù)被破壞,此時(shí)即使能夠抑制部分偽影,有效的醫(yī)學(xué)信息仍然會(huì)被部分破壞,因此結(jié)合超聲醫(yī)師的經(jīng)驗(yàn),有效信息完整程度的評(píng)價(jià)等級(jí)為一般.隨著截止半徑不斷地?cái)U(kuò)大,會(huì)有部分偽影在頻域上的脈沖信號(hào)點(diǎn)被囊括在截止半徑之外,因此部分偽影無(wú)法得到抑制,此時(shí)偽影抑制的評(píng)價(jià)等級(jí)逐漸從一般到很差.

閾值系數(shù)K對(duì)于偽影抑制的效果同樣具有非常重要的作用,當(dāng)閾值系數(shù)過(guò)大時(shí),疑似偽影的脈沖點(diǎn)的幅度值沒(méi)有超過(guò)閾值,因此得以保留,此時(shí)偽影抑制的評(píng)價(jià)等級(jí)為很差.隨著閾值系數(shù)的降低,結(jié)合醫(yī)師的判斷,偽影抑制效果逐漸從一般轉(zhuǎn)變?yōu)榱己?

表1 為不同截止頻率在超聲圖像偽影抑制前后區(qū)域A(偽影區(qū))和區(qū)域B(非偽影區(qū))的均方誤差,信噪比結(jié)果以及偽影抑制效果和有效信息的完整程度的對(duì)比,表2 為不同閾值系數(shù)在超聲圖像偽影抑制前后區(qū)域A(偽影區(qū))和區(qū)域B(非偽影區(qū))的均方誤差和信噪比結(jié)果以及偽影抑制效果和有效信息的完整程度的對(duì)比(偽影抑制效果評(píng)估標(biāo)準(zhǔn))

表1 超聲圖像不同截止半徑處理前后區(qū)域A (偽影區(qū))和區(qū)域B (非偽影區(qū))的均方誤差和信噪比結(jié)果對(duì)比

通過(guò)表1 的結(jié)果可以判斷出在R＜ 3 時(shí),頻域中心有用信息被抑制,導(dǎo)致通過(guò)IFFT 后的圖像部分有用信息及偽影的損失,以在區(qū)域A及區(qū)域B的MSE值偏大,當(dāng)R＞10 時(shí),由于部分頻域上偽影亮點(diǎn)在截止半徑之內(nèi)沒(méi)有進(jìn)行任何壓低處理,導(dǎo)致通過(guò)IFFT 后仍然會(huì)有部分周期性偽影存在(如圖17),在此截止半徑下圖像區(qū)域A及區(qū)域B的MSE值會(huì)非常小.當(dāng)R位于3 和6 之間時(shí),由于部分周期性偽影分布在該區(qū)域內(nèi),對(duì)該區(qū)域內(nèi)的亮點(diǎn)幅值進(jìn)行適當(dāng)壓低,可以達(dá)到去除周期性偽影的效果.

取R=4,來(lái)對(duì)P集合中[0,M/2-R]和[M/2+R,M]范圍內(nèi)的亮點(diǎn)通過(guò)閾值系數(shù)K來(lái)進(jìn)行幅值壓低處理獲取表2 的數(shù)據(jù).

表2 超聲圖像不同閾值系數(shù)處理前后區(qū)域A (偽影區(qū))和區(qū)域B (非偽影區(qū))的均方誤差和信噪比結(jié)果對(duì)比

由于當(dāng)閾值系數(shù)K不斷增大,集合P中滿足當(dāng)前閾值的亮點(diǎn)逐漸減少,當(dāng)K等于0.8 時(shí),所有亮點(diǎn)都不滿足超過(guò)閾值的條件,因此沒(méi)有對(duì)集合P中的任何亮點(diǎn)進(jìn)行壓低,所以處理前后的區(qū)域A及區(qū)域B的MSE值近似于0,證明K值與偽影抑制的效果呈負(fù)相關(guān),即K越大偽影抑制效果越差.

根據(jù)采樣結(jié)果,在當(dāng)前超聲圖像中能夠確定最優(yōu)的截止半徑R和閾值系數(shù)K后通過(guò)該算法進(jìn)一步實(shí)現(xiàn)的部分連續(xù)超聲圖像偽影抑制效果,如圖21、圖22.

圖21 未處理的超聲低通道圖像

圖22 偽影抑制后的超聲低通道圖像

實(shí)驗(yàn)結(jié)果表明,對(duì)于原始超聲低通道掃描圖像包含的周期性偽影,該算法基本能在保留圖像有用信息的基礎(chǔ)上將其去除,并且該算法效果明顯,實(shí)驗(yàn)結(jié)果證明了該算法模型的正確性.

4 結(jié)論及下一步研究

醫(yī)學(xué)超聲實(shí)時(shí)成像的研究是一項(xiàng)比較新型的研究方向,近幾年的發(fā)展非常迅速,從早期的A 型,M 型,B 型超聲診斷技術(shù)到目前的超聲多普勒成像技術(shù),三維成像技術(shù)的發(fā)展歷程來(lái)看,超聲成像技術(shù)的發(fā)展在于提高圖像質(zhì)量,已經(jīng)取得較大的突破[20],并且廣泛應(yīng)用于醫(yī)生的輔助診斷等領(lǐng)域,在超聲實(shí)時(shí)成像領(lǐng)域也有著廣泛的研究空間,伴隨著現(xiàn)代醫(yī)學(xué)的進(jìn)步以及發(fā)展,醫(yī)學(xué)超聲診斷儀的重要性不言而喻.超聲實(shí)時(shí)成像的質(zhì)量確面臨的巨大的挑戰(zhàn),由于超聲聲波固有的特點(diǎn)和技術(shù)缺陷,超聲實(shí)時(shí)成像上往往都帶有斑點(diǎn),周期性等類型的偽影,這些偽影往往會(huì)嚴(yán)重影響到醫(yī)生的診斷效果,所以如何抑制這些類型的偽影成為了一個(gè)新的課題.

本研究對(duì)于超聲實(shí)時(shí)成像的周期性偽影抑制提出了一種新的算法,用來(lái)解決目前遇到的難題,通過(guò)正弦波的模擬實(shí)驗(yàn)以及人體組織的算法仿真實(shí)驗(yàn),該算法能夠較好地解決人體組織上周期性偽影影響醫(yī)生診斷的難題,視頻中的周期性偽影得到了明顯的抑制,能夠提高醫(yī)生診斷的準(zhǔn)確度.偽影的來(lái)源是由超聲探頭制造上的局限性以及采用低通道掃描而產(chǎn)生的,可以發(fā)現(xiàn)該偽影垂直于人體的組織器官,因此有別于超聲波束作用于不同紋理的器官組織的表面產(chǎn)生回波導(dǎo)致的偽影.本研究當(dāng)中的偽影,與器官組織的紋理無(wú)關(guān),對(duì)于任何人體組織器官,在超聲領(lǐng)域內(nèi)采用該類型探頭以及掃描方式均會(huì)產(chǎn)生該類型的偽影,因此針對(duì)不同的器官組織不會(huì)有任何區(qū)別.

盡管如此,在本研究中,還有一些地方需要改進(jìn),例如對(duì)于截止頻率的判定以及閾值系數(shù)的選取由于采樣周期較大,只能確定一個(gè)大致的范圍,并不能準(zhǔn)確的定位出具體數(shù)值,在極端情況下可能會(huì)出現(xiàn)“誤判”,將非偽影區(qū)域當(dāng)成偽影給抑制掉,需要通過(guò)下一步研究來(lái)對(duì)算法模型進(jìn)行改進(jìn).

當(dāng)前,超聲圖像的修復(fù)技術(shù)尚未成熟,圖像修復(fù)不是一個(gè)單一學(xué)科的問(wèn)題,它涉及很多相關(guān)領(lǐng)域的研究[21],隨著科學(xué)技術(shù)的不斷發(fā)展,醫(yī)學(xué)超聲領(lǐng)域?qū)?huì)不斷地進(jìn)行深耕,新的醫(yī)學(xué)超聲技術(shù)也將會(huì)出現(xiàn),這一行業(yè)會(huì)隨著技術(shù)的不斷成熟和完善不斷蓬勃發(fā)展,在保證人體健康的過(guò)程中發(fā)揮越來(lái)越大的作用.