李明豪, 尹 皓, 趙海娜, 劉東權(quán)

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

2(四川大學(xué)華西醫(yī)院 超聲醫(yī)學(xué)科, 成都 610041)

頸部淋巴結(jié)會受甲狀腺病變細(xì)胞轉(zhuǎn)移的影響而發(fā)生病變. 發(fā)生腫瘤轉(zhuǎn)移的淋巴結(jié), 腫瘤細(xì)胞會通過淋巴管首先聚集于邊緣竇, 繼而逐漸侵入淋巴結(jié)實質(zhì); 其可破壞淋巴結(jié)正常的動靜脈結(jié)構(gòu), 誘導(dǎo)新生血管的形成. 腫瘤細(xì)胞的轉(zhuǎn)移會影響頸部淋巴結(jié)的血流情況. 基于超聲造影(contrast-enhanced ultrasound, CEUS)可以通過造影劑判斷血流灌注的信息, 輔助醫(yī)生進(jìn)行病情診斷[1].

傳統(tǒng)的灌注曲線(時間強(qiáng)度曲線, TIC)是基于病灶區(qū)域每幀取強(qiáng)度值的平均值而沿時間軸成時間序列而得. 然而, 由于感興趣區(qū)域(region of interest, ROI)的灌注強(qiáng)度可能并不均勻, 這樣將使提取的特征不準(zhǔn)確[2].因此, 本文對感興趣區(qū)域進(jìn)行像素級分析, 即視頻的一個坐標(biāo)將沿時間軸擬合一條時間強(qiáng)度曲線. 所有感興趣區(qū)域內(nèi)的曲線對一個灌注特征(達(dá)峰時間、峰值等)進(jìn)行提取可以形成一個二維矩陣. 二維矩陣可以進(jìn)行編碼顯示成圖, 也可以進(jìn)一步提取更深層次的信息來進(jìn)行可視化.

本文主要的創(chuàng)新實驗是, 基于達(dá)峰時間(time to peak,TTP)二維矩陣進(jìn)行研究, 猜測感興趣區(qū)域中相距不遠(yuǎn)的達(dá)峰時間接近的兩坐標(biāo)可能存在流向先后的關(guān)聯(lián)性并因此尋找血流流向的信息, 來進(jìn)行灌注流向的成像可視化.

1 相關(guān)技術(shù)與研究現(xiàn)狀

1.1 研究現(xiàn)狀

探頭固定時, 頸部淋巴結(jié)采集的造影視頻, 目標(biāo)跟蹤過程是不考慮尺度變化的. 比如, Zhang等人在2011年認(rèn)為頸部淋巴結(jié)在探頭不動的情況主要會受呼吸作用影響而做類似周期性的運(yùn)動, 而直接利用基于塊匹配的方法來進(jìn)行目標(biāo)跟蹤[3]. 其利用了最小化絕對差和(sum of absolute difference, SAD). 基于光流法也能進(jìn)行目標(biāo)跟蹤, 但是對光照要求太過嚴(yán)格[4]. Ta等人在2012年利用了基于互信息的方法來進(jìn)行目標(biāo)跟蹤[5].然而互信息始終計算量太大, 并且也不能很好地處理異常幀. 而后, 林細(xì)林等人2018年利用了感知壓縮跟蹤[6], 該方法通過改進(jìn)可以處理異常幀的情況[7], 但壓縮感知跟蹤的跟蹤結(jié)果位置不夠精確, 甚至位置還會跑得很遠(yuǎn)[8]. 由此本文實驗了KCF (核相關(guān)濾波)算法來進(jìn)行目標(biāo)跟蹤, 其利用了相關(guān)濾波理論與核函數(shù)技術(shù), 使得目標(biāo)跟蹤很魯邦而且計算量不會很大[9].

時間強(qiáng)度曲線可以反映血液中的造影劑濃度變化.Zhang等人早期就對感興趣區(qū)域提取亮度灰度值的平均值, 隨幀數(shù)生成時間強(qiáng)度曲線(time intensity curves,TIC), 提取出達(dá)峰時間、峰值(peak value, PI)等灌注參數(shù)來分析視頻[3]. Gaia分析類風(fēng)濕性關(guān)節(jié)炎時, 指出可以利用像素級的時間強(qiáng)度曲線, 即每個坐標(biāo)點(diǎn)擬合一條曲線. 然后對每一個灌注參數(shù)可以生成一個二維矩陣, 可以編碼成像來幫助醫(yī)生診斷病情. 其擬合曲線時利用的是Gamma公式[10]. 王本剛等人分析肝臟病例也進(jìn)行了基于像素的時間強(qiáng)度曲線分析, 不過只是利用了savgol濾波器進(jìn)行曲線擬合[2]. 這些文獻(xiàn)沒有研究怎么樣將沒有很好灌注趨勢的曲線篩除, 本文將對此進(jìn)行研究.

針對灌注流向問題, 本文嘗試?yán)昧骶€與向量場來進(jìn)行實驗. 自Brian 等提出線積分卷積 (line integral convolution, LIC)的方法后, 對向量場進(jìn)行可視化大部分是利用了LIC方法及其變體[11]. 雖然Wegenkittl 等人提出 Oriented LIC 方法來使結(jié)果紋理含有方向指向信息[12]. 但本文病例感興趣區(qū)域太小, 因此不能用該方法表示方向. 本文使用了Brain提出的基本的LIC方法, 同時利用顏色編碼來表示坐標(biāo)的先后關(guān)系.

1.2 相關(guān)技術(shù)

第1步. 對淋巴結(jié)進(jìn)行目標(biāo)跟蹤.

實驗中所采集的超聲造影視頻都為雙模式并排顯示, 如圖1所示, 左為超聲灰階B模式, 右為超聲造影模式. 目標(biāo)跟蹤在灰階B模式上面進(jìn)行, 然后再將位置信息應(yīng)用到彩色造影模式中以分析灌注情況.

圖1 雙模式視頻幀

由于淋巴結(jié)主要是因為呼吸會產(chǎn)生運(yùn)動, 但沒有明顯尺度變化. 所以可以由醫(yī)生進(jìn)行感興趣區(qū)域的標(biāo)記, 然后提取出掩膜mask以便在視頻中進(jìn)行目標(biāo)跟蹤, 如圖2所示.

圖2 提取mask掩膜示意圖

這里手工標(biāo)定的幀由醫(yī)生選取, 被用來做目標(biāo)跟蹤的起始幀, 被稱為“參考幀”. 參考幀不一定是原視頻的第1幀.

本文選擇核相關(guān)濾波(kernelized correlation filters,KCF)[9]和Kalman濾波結(jié)合的方式進(jìn)行目標(biāo)跟蹤[13].

第2步. 進(jìn)行像素級的時間強(qiáng)度曲線擬合.

基本思路是: 在B模式對淋巴結(jié)進(jìn)行目標(biāo)跟蹤后,將位置信息用到造影模式. 造影模式ROI每一個坐標(biāo)沿著時間軸都能擬合一條時間強(qiáng)度曲線. 本文利用Gamma公式模型, 同時也研究怎么篩選符合灌注特征的曲線.

第3步. 提取灌注特征參數(shù). 因為是基于像素的分析, 所以對符合標(biāo)準(zhǔn)的曲線提取灌注參數(shù)(達(dá)峰時間、峰值等), 形成二維矩陣. 對二維矩陣可以編碼成像顯示, 輔助醫(yī)生診斷病情.

第4步. 利用TTP參數(shù)二維矩陣可視化灌注流向情況. 設(shè)種子點(diǎn)為TTP值低于設(shè)定閾值的坐標(biāo)點(diǎn). 先找種子坐標(biāo)點(diǎn), 每個種子點(diǎn)會找一條流線, 這些流線可以進(jìn)行可視化. 實驗也嘗試了利用流線信息生成向量場,并利用線積分卷積的方式可視化.

2 病灶跟蹤

2.1 淋巴結(jié)關(guān)于目標(biāo)跟蹤方面的特點(diǎn)

探頭固定的情況下, 視頻灰階B模式相鄰兩幀之間的運(yùn)動變化很小可以忽略, 整體運(yùn)動以呼吸運(yùn)動的影響為主, 而呼吸運(yùn)動可以看成周期運(yùn)動[3]. 在視頻中淋巴結(jié)很少發(fā)生巨大的尺度變化(即縮放變化), 也不會有大的形變(指形狀改變成另外的樣子).

灰階B模式采用的是基波, 造影模式采用的是諧波. 造影劑的濃度變化會使彩色造影部分的亮度跟著改變. 這也會對灰階B模式產(chǎn)生影響, 但亮度不會發(fā)生很大變化.

視頻中會出現(xiàn)短暫的異常幀, 持續(xù)時間很短, 并且之后目標(biāo)形態(tài)會恢復(fù)正常. 異常幀常是由于信號噪聲、病人的較強(qiáng)的呼吸引來的其他平面信息等因素而導(dǎo)致的, 主要表現(xiàn)是目標(biāo)細(xì)節(jié)暫時丟失.

總結(jié)下來, 沒有異常幀時可以利用簡單目標(biāo)跟蹤方法(如基于灰度的模板匹配)來處理. 有異常幀時就要考慮如何檢測與處理異常幀.

2.2 目標(biāo)跟蹤模型

傳統(tǒng)的基于灰度的模板匹配(如最小化絕對差和SAD、最小化平均絕對差MAD、最大化歸一化互相關(guān)NCC等)進(jìn)行目標(biāo)跟蹤, 一幀一幀更換參考模板的可以漸變地跟蹤目標(biāo). 但是, 其并沒有提出好的檢測異常幀的方法.

基于要處理異常幀的思路, 本文選擇了相關(guān)濾波的方法進(jìn)行目標(biāo)跟蹤.

相關(guān)濾波的思想是利用信號的相關(guān)性, 設(shè)計一個濾波模板, 利用模板與目標(biāo)候選區(qū)做相關(guān)運(yùn)算生成響應(yīng)圖, 響應(yīng)圖上值越大的位置越可能是待搜索幀的目標(biāo)位置. 相關(guān)濾波跟蹤的優(yōu)點(diǎn)是可以利用卷積計算在傅里葉域的性質(zhì)減少計算量.

用KCF (核相關(guān)濾波)進(jìn)行目標(biāo)跟蹤, 利用HOG特征會對光照變化有魯邦性, 異常幀檢測可以利用響應(yīng)圖的峰值旁瓣比(peak to sidelobe ratio,PSR)[9].

峰值旁瓣比公式為[14]:

其中,gmax為相關(guān)響應(yīng)圖最大峰值, 峰值周圍半徑為11像素大小的圓形區(qū)域之外為旁瓣, μs與δs分別為旁瓣的均值與方差.

2.3 卡爾曼濾波

卡爾曼濾波器是一種高效的遞歸濾波器, 能夠從一系列不完全包含噪聲的測量中, 估計出動態(tài)系統(tǒng)的狀態(tài)[13]. 其常用在軌跡預(yù)測與目標(biāo)定位中. 對于一般的線性隨機(jī)系統(tǒng), 可表示為:

其中,xk∈Rn表示狀態(tài)向量,zk∈Rm表示測量向量,與Gk表示狀態(tài)轉(zhuǎn)移矩陣,Hk表示測量矩陣,wk是過程噪聲,vk是 測量噪聲.wk的協(xié)方差矩陣為Qk,vk的協(xié)方差矩陣為Rk. 實驗中假定Qk與Rk均為零均值高斯白噪聲, 且兩者互不相關(guān).

上述模型是線性的Kalman濾波模型, 常有5個基本公式. 用Kalman濾波進(jìn)行數(shù)據(jù)的濾波, 可以從后驗結(jié)果得到最終的數(shù)據(jù)濾波結(jié)果, 比如濾波后的曲線數(shù)據(jù); 用Kalman濾波進(jìn)行預(yù)測, 可以從先驗結(jié)果里得到預(yù)測結(jié)果, 比如預(yù)測的位置信息.

2.4 總的目標(biāo)跟蹤思路

視頻中大部分時間淋巴結(jié)的運(yùn)動速度都是較慢的,但突發(fā)情況時會發(fā)生較快的運(yùn)動. 較快的運(yùn)動就會有一些幀中目標(biāo)變得模糊, 或產(chǎn)生偽影. 因此建立線性的卡爾曼模型來檢測目標(biāo)是否運(yùn)動過快, 一旦速度超過閾值, 就不對KCF濾波模板進(jìn)行更新[7].

最終, KCF的模板大小為長寬96 像素. 如果圖像太小, 模板邊長要適當(dāng)變小. PSR的最小閾值為40,Kalman預(yù)測位置與相關(guān)濾波檢測位置的歐式距離最大閾值為60 像素.

如果上述兩閾值條件都滿足, 才更新KCF的濾波模板. 對于Kalman模型的參數(shù)更新, 如果某幀不滿足PSR最小閾值就用Kalman預(yù)測位置更新; 如果滿足PSR最小閾值, 就用KCF檢測位置更新. 目標(biāo)跟蹤中的更新操作流程如圖3所示.

圖3 目標(biāo)跟蹤中的更新操作流程

3 進(jìn)行像素級的時間強(qiáng)度曲線擬合

3.1 二維平滑窗口的選擇

時間強(qiáng)度曲線是每個坐標(biāo)都沿時間軸提取一條曲線. 但是會出現(xiàn)大量坐標(biāo)的數(shù)據(jù)不能曲線擬合成功.

針對該問題, 猜測淋巴結(jié)內(nèi)部各結(jié)構(gòu)會隨著呼吸相對于淋巴結(jié)整體發(fā)生很細(xì)微的運(yùn)動. 考慮利用一個二維平滑窗口, 將呼吸引起的細(xì)微運(yùn)動的軌跡籠罩住,窗口內(nèi)提取平均值. 平均值的變化就能符合灌注趨勢了. 這里平滑窗口高寬是“29×29”.

3.2 TIC曲線數(shù)據(jù)的提取

TIC (time-intensity curve)曲線即時間強(qiáng)度曲線. 將造影圖片從RGB顏色空間轉(zhuǎn)換為“HSV”顏色空間[15].生成H分量的直方圖之后發(fā)現(xiàn)數(shù)據(jù)分布集中在很窄的一段, 如圖4, 即造影部分色調(diào)都比較接近. 因此, 圖像的色調(diào)不存在太多信息, 飽和度又不能說明亮度變化.最終選擇了“HSV”的“V”分量(縮放到[0, 255])來表示“強(qiáng)度”.

圖4 一幀造影與其HSV空間的H分量的直方圖

每一坐標(biāo)點(diǎn)可以沿時間軸得到一條TIC曲線的原始數(shù)據(jù). 當(dāng)TIC的原始數(shù)據(jù)得到之后, 會首先采用S-G平滑濾波(Savitsky-Golay濾波器)[2], 多項式階數(shù)是2,窗口大小是81 (原視頻幀率是10 fps, 所以這里代表8.1 s). 濾波目的是減弱噪聲.

3.3 基于Gamma公式的曲線擬合

Gamma模型公式[10]為:

其中,A≥0,C＞0,k＞0 ,B是初始值即t=0時刻的y值.B為固定值, 不屬于可變參數(shù),AT(arrival time)為造影劑到達(dá)時間.

利用Gamma公式進(jìn)行曲線擬合, 可以描述造影劑在血管中的洗入(wash-in, 表現(xiàn)為曲線上升)與洗出(wash-out, 表現(xiàn)為曲線下降)的時間段情況, 便于提取相關(guān)的灌注特征參數(shù), 如圖5所示.

圖5 原數(shù)據(jù)和Gamma擬合圖

3.4 曲線數(shù)據(jù)的篩選之粗篩選與細(xì)篩選

并非ROI內(nèi)所有的點(diǎn)都能很好地擬合曲線, 考慮原因是并非所有區(qū)域都是液態(tài)的有血流經(jīng)過的區(qū)域.因此, 需要將不符合的曲線數(shù)據(jù)篩除.

這里考慮用“粗篩選-細(xì)篩選”的方式, 只有兩種篩選都通過的曲線才能用于以后的實驗.

第1步, 粗篩選.

粗篩選大致分兩個思路, 兩個子篩選條件都滿足了才算通過粗篩選. 一個是真實值的判斷, 一個是縮放值的判斷.

(1)真實值判斷, 就是對曲線原數(shù)據(jù)的峰值、起始值判斷大小.

比如在實驗中, 峰值不能小于30, 曲線起始值不大于100. 這里注意由于利用HSV的方法, 所有曲線中強(qiáng)度值范圍可統(tǒng)一為[0, 255], 如圖6所示.

圖6 真實值判斷示意圖

(2)縮放值判斷, 將曲線原數(shù)據(jù)值縮放到[0,n],n為正整數(shù), 這里為8. 對數(shù)據(jù)中的第2峰值、后面(1/2時間之后)最大值、起伏過程中的上升幅度、起伏過程中的異常上升總次數(shù)等進(jìn)行判斷.

比如實驗中, 第2峰值在數(shù)據(jù)后面不能大于5, 在前面不能相比之前的波谷上升超過3. 后端1/5時間段最大值不能是6. 上升幅度大于3為異常上升, 異常上升次數(shù)不能多于3個, 如圖7所示.

圖7 縮放值判斷示意圖

第2步, 細(xì)篩選.

滿足粗篩選的曲線再次用細(xì)篩選進(jìn)行處理. 細(xì)篩選是利用Gamma公式去擬合數(shù)據(jù), 根據(jù)擬合情況對數(shù)據(jù)進(jìn)行最后篩選.

因為“R-squared”一般應(yīng)用于線性擬合中的擬合度評價準(zhǔn)則, 所以沒有使用在這里[16]. 本文的擬合度的評價準(zhǔn)則如下:

(1)NRMSE, 即歸一化的RMSE(均方根誤差).

這里RMSE公式為:

其中,ri是 原始數(shù)據(jù),fi是擬合后的數(shù)據(jù).

而NRMSE的公式為:

NRMSE的最大閾值為0.15.

(2) Adjusted Cosine即修正的余弦相似性[17]. 公式為:

其中,u即原數(shù)據(jù)r與擬合數(shù)據(jù)f合在一起共同的平均值. 使用修正余弦相似性有一個重要條件, 需要兩個時間序列一樣長.

考慮到擬合的曲線可能會直接平滑掉起伏頻率高的數(shù)據(jù)段, 可能出現(xiàn)NRMSE小但是跟原數(shù)據(jù)根本就不像的情況. 因此, 利用余弦相似度進(jìn)行衡量, 這里設(shè)置的最小閾值為0.75.

4 提取灌注特征參數(shù)

4.1 提取2-D灌注參數(shù)特征

對于一條TIC曲線可以提取灌注參數(shù), 而參數(shù)大致分兩類[5]: (1)基于時間的, 如TTP (達(dá)峰時間)、MTT(平均渡越時間); (2)基于數(shù)值的, 如AUC (曲線下面積)、PI (峰值).

在視頻中, 像素級提取時間強(qiáng)度曲線, 則對于一個灌注特征參數(shù)可得到一個二維矩陣, 可稱為“灌注參數(shù)二維矩陣”.

另外, 曲線擬合成功與失敗的坐標(biāo)需要記錄下來,下文需要使用.

4.2 灌注參數(shù)二維矩陣的可視化

曲線擬合成功的點(diǎn)才能進(jìn)行正常的顏色編碼來可視化, 失敗的坐標(biāo)點(diǎn)將會設(shè)置為黑色.

首先會將參數(shù)二維矩陣中成功的坐標(biāo)的值統(tǒng)一縮放到[0, 240]之間. 因為在HSV顏色空間中, 色調(diào)分量H從“0到240”會對應(yīng)“紅到藍(lán)”. 將S和V統(tǒng)一設(shè)置為1 (或滿值), 不同坐標(biāo)點(diǎn)只需設(shè)置各自的H分量, 便可得到一張彩色圖.

然后, 對縮放后的二維矩陣進(jìn)行編碼可顯示成圖片. 不過, 顏色編碼方式需結(jié)合實際情況指定.

比如, 對于TTP (達(dá)峰時間)來說, 達(dá)峰時間越小就說明達(dá)峰時間早, 越早越是醫(yī)生感興趣的越偏向暖色(紅色). 此這一類數(shù)據(jù)采取“藍(lán)到紅, 大到小”的編碼方法, 如圖8.

圖8 TTP二維矩陣成像圖

而對于AUC (曲線下面積), 值越大說明灌注總量越大, 越表示醫(yī)生感興趣的而偏暖色(紅色). 因此這一類數(shù)據(jù)采取“藍(lán)到紅, 小到大”的編碼方法, 如圖9.

圖9 AUC二維矩陣成像圖

5 基于達(dá)峰時間的灌注流向估計與可視化

5.1 灌注方向的求法

灌注大致方向?qū)τ诓∏榈呐袛嗍怯袔椭? 所以,為了可視化灌注方向, 本課題利用了達(dá)峰時間二維矩陣進(jìn)行估計. 總體思路是, 造影會從TTP小的地方流向TTP大一點(diǎn)的地方.

5.2 流線種子點(diǎn)的選取

針對達(dá)峰時間二維矩陣, 設(shè)置一個TTP最大閾值,將低于最大閾值的區(qū)域提取出來作為種子點(diǎn).

對于TTP二維矩陣查找最大值maxValue和最小值minValue, 然后閾值Th為:

其中,frac為“0”到“1”之間的小數(shù), 這里frac設(shè)置為0.3.

5.3 每條流線的生成

流線是一條其上任意一點(diǎn)的矢量方向都與流線相切的曲線. 流線方程定義如下[18]:

其中, μ(τ) 代 表空間中點(diǎn)的位置,v(μ(τ))為流線上該點(diǎn)的切線方向,v是τ 的函數(shù), τ可以是時間或弧長等參量.

本文一條流線可以利用一組具有先后順序的采樣點(diǎn)來確定. 這里利用鏈表的形式給出, 為了理解, 這里借用C++的容器表示.

class Point{

int x;

int y;

};

std::vector＜Point＞ streamlineList; //即一條流線的點(diǎn)組成的list.

而這里尋找流線的思想: 每個點(diǎn)以自己為起點(diǎn)找下一個符合的點(diǎn), 并以找到的點(diǎn)為起點(diǎn)再找下一個點(diǎn).如圖10(a), 一條流線的例子, 從A點(diǎn)開始, 直到F點(diǎn)結(jié)束就形成了一條流線的采樣點(diǎn).

圖10 尋找流線示意圖

每次一個點(diǎn)為起點(diǎn)尋找下一個點(diǎn)時, 找大于該點(diǎn)且時間差為“0.2”或“0.3”秒的最小值(兩個點(diǎn)不能同時發(fā)生, 所以時間差要有最小閾值), 搜索范圍并不是采用固定的8鄰域, 而是以當(dāng)前點(diǎn)為中心, 一層一層(1、2到3的順序)地搜尋外圍的鄰域點(diǎn), 直到找到一個TTP大于該點(diǎn)且超過一定數(shù)值的最小值.

如圖10(b), 以當(dāng)前點(diǎn)P搜索下一個候選點(diǎn)的過程如下:

(1)設(shè)當(dāng)前參考點(diǎn)為P, 并設(shè)置預(yù)定標(biāo)準(zhǔn)為下一點(diǎn)是該點(diǎn)值是比當(dāng)前參考點(diǎn)大0.2的點(diǎn)的最小值.

(2)搜索P點(diǎn)的8鄰域點(diǎn), 如圖中標(biāo)記為“1”的那些點(diǎn). 如果找到預(yù)定標(biāo)準(zhǔn)的候選點(diǎn)就結(jié)束; 否則將當(dāng)前點(diǎn)和其8鄰域點(diǎn)加入一個點(diǎn)群pointGroup, 然后進(jìn)行第(3)步;

(3)繼續(xù)搜索pointGroup點(diǎn)群的4鄰域點(diǎn)(比如上次搜索標(biāo)記為“1”, 這次搜索標(biāo)記為“2”的那些點(diǎn)). 當(dāng)前參考點(diǎn)依舊為中心點(diǎn)P (即第(1)步的點(diǎn)). 如果找到預(yù)定標(biāo)準(zhǔn)的點(diǎn)就結(jié)束; 否則就將其周圍點(diǎn)加入pointGroup點(diǎn)群, 重復(fù)本步驟. 另外, 此步驟中pointGroup點(diǎn)群的長度大于限制值或ROI內(nèi)所有點(diǎn)都已經(jīng)遍歷也會終止.

5.4 流線顯示

找到每條流線采樣點(diǎn)中, 一條流線上每相鄰的兩個點(diǎn)可能在圖上空間不相鄰, 可以利用Bresenham直線算法用一條近似直線的像素點(diǎn)連接兩點(diǎn)[19].

可問題是, 方向怎么表示. 圖像小, 不考慮用箭頭表示流向, 而利用TTP參數(shù)成像圖里的顏色(色調(diào))不同進(jìn)行標(biāo)記.

另外, 對區(qū)域中每一個坐標(biāo)點(diǎn)統(tǒng)計有多少條流線經(jīng)過(遍歷所有曲線得出), 可以生成一個流線計數(shù)二維矩陣, 流線經(jīng)過比較多的區(qū)域有很大可能是有血管的區(qū)域. 由于流線計數(shù)圖亮度低的區(qū)域的細(xì)節(jié)看不清.可以利用log增強(qiáng)技術(shù)來提高低灰度區(qū)域的對比度.

流線圖、流線計數(shù)圖以及流線計數(shù)圖log增強(qiáng)圖都顯示在圖11中.

圖11 流線結(jié)果顯示

6 生成向量場并可視化

6.1 生成向量場

不同的種子點(diǎn)會有不同的流線, 可生成大小與原參數(shù)成像圖相同的向量場, 每個坐標(biāo)點(diǎn)有一個向量, 為了表示每一個坐標(biāo)點(diǎn)的造影劑在下一時刻的可能去向[18].

一條流線可以利用Bresenham直線算法填充空隙,然后從前到后依次指向, 每個非終止點(diǎn)可以指向流線上下一點(diǎn)而成向量. 而每一條流線的終點(diǎn)會設(shè)置成零向量.

可以用一張mask來表明某坐標(biāo)是否已經(jīng)因為一條流線經(jīng)過而被賦值過. 遍歷每一條流線, 然后在向量場相應(yīng)位置修改向量值. 只有流線經(jīng)過的點(diǎn)有非零的向量, 其他點(diǎn)向量都為零向量. 當(dāng)出現(xiàn)某個點(diǎn)有多條流線經(jīng)過的時候, 選取流線計數(shù)大的下一個點(diǎn)來作為下一個最優(yōu)指向.

6.2 線積分卷積

LIC技術(shù)根據(jù)向量場, 用一維濾波器沿流線卷積白噪聲圖像, 合成紋理. 該紋理與向量場的向量方向高度相關(guān).

其生成過程需要輸入向量場和一張噪聲圖(通常是白噪聲)[11], 如圖12所示.

圖12 LIC圖生成示意圖

其基本思想是, LIC會選擇噪聲圖像作為輸入背景紋理, 然后對每一個坐標(biāo)點(diǎn)正反向去各追蹤一條流線, 流線上的每一個像素點(diǎn)根據(jù)卷積核算出的權(quán)重來進(jìn)行卷積計算, 得出紋理值. 而采用不同的核函數(shù), 加之不同的背景噪聲, 最終得到的效果也不一樣[20], 本實驗選擇的是盒型核與白噪聲.

6.3 灌注流向的向量場成像

生成LIC圖時要對零向量區(qū)域進(jìn)行一步抑制, 縮小其亮度范圍, 利用了TTP成像圖的顏色來著色表示坐標(biāo)點(diǎn)的先后關(guān)系(越接近暖色即紅色, 則達(dá)峰時間越小).

如圖13, 是對上面流線圖的例子進(jìn)行計算的結(jié)果.其紋理顯示了向量場的大致狀態(tài).

圖13 LIC結(jié)果圖

7 實驗結(jié)果分析

7.1 數(shù)據(jù)

本課題所選擇的病例視頻由四川成都華西醫(yī)院提供, 造影劑采用意大利博萊科公司生產(chǎn)的第二代造影劑聲諾維. 所有病例來自于同一個邁瑞的機(jī)器. 采集期間, 探頭固定, 囑患者固定體位, 平穩(wěn)呼吸, 勿做吞咽動作. 目前實驗數(shù)據(jù)共18例病例, 12例良、6例惡.

7.2 曲線擬合篩選后的有效點(diǎn)面積占比

如表1所示, 18例視頻都進(jìn)行了實驗. 有效坐標(biāo)點(diǎn)是上文曲線擬合成功坐標(biāo)點(diǎn). 有效坐標(biāo)點(diǎn)占比是有效坐標(biāo)點(diǎn)數(shù)占ROI內(nèi)總點(diǎn)數(shù)的比重.

表1 有效坐標(biāo)點(diǎn)占比結(jié)果

對于造影視頻, 符合灌注趨勢的坐標(biāo)應(yīng)是有血流經(jīng)過的坐標(biāo). 有效坐標(biāo)點(diǎn)占比太小則合格血管的覆蓋面積太少.

因為本文的重點(diǎn)是探尋灌注流向信息可視方法,所以實驗的數(shù)據(jù)不是很多. 但在已有的數(shù)據(jù)中, 大部分良性視頻的有效成功點(diǎn)數(shù)占比為“50%–100%”, 而有一半的惡性視頻的有效成功點(diǎn)數(shù)占比為“0–50%”. 如果發(fā)生了有效成功點(diǎn)占比低于50%的事情, 樣本極有可能是惡性的, 也可能是采集異常, 可標(biāo)記為“可疑的”.

7.3 4種不同面積占比的視頻的灌注可視化結(jié)果

如圖14, 觀察不同的擬合成功點(diǎn)占比下灌注情況可視化結(jié)果. 其中, 每一行圖片為一個視頻的結(jié)果, 從上而下有效坐標(biāo)占比依次為“0–20%”“20%–50%”“50%–80%”“80%–100%”. 每一行從左到右分別為TTP成像圖、流線圖、流線技術(shù)圖(log增強(qiáng))、LIC結(jié)果圖.

圖14 4種不同有效點(diǎn)面積占比的結(jié)果

當(dāng)有效坐標(biāo)點(diǎn)占比過于小時(低于50%時), 有關(guān)灌注情況的效果會很差. 因此, 建議觀察流向信息時,有效坐標(biāo)占比要超過50%.

7.4 總的結(jié)果分析

當(dāng)確定ROI內(nèi)基本是血管時, 可以將有效坐標(biāo)點(diǎn)占比小(比如低于50%)的視頻標(biāo)記為“可疑”視頻供醫(yī)生分析. 有效坐標(biāo)點(diǎn)(擬合曲線成功的坐標(biāo))占比太小,不建議進(jìn)行灌注可視化. 利用上述方法尋找到的流線,流線圖與流線計數(shù)圖非常類似于血管的分布, 應(yīng)具有一定的研究價值.

8 結(jié)論與展望

本文研究了對造影視頻的灌注特征的提取與可視化. 提取灌注參數(shù)(達(dá)峰時間、峰值等)來編碼成像可以協(xié)助醫(yī)生分析病情, 使用流線等方法可視化造影劑在ROI內(nèi)的流向信息. 結(jié)果中流線圖與流線計數(shù)圖中的線條類似血管, 應(yīng)具有一定的研究意義.

經(jīng)歷了目標(biāo)跟蹤、時間強(qiáng)度曲線擬合、灌注參數(shù)提取、根據(jù)達(dá)峰時間二維矩陣進(jìn)行灌注情況可視化共4個步驟. 目標(biāo)跟蹤根據(jù)視頻的特點(diǎn)而選擇了KCF與Kalman濾波結(jié)合的方式. 時間強(qiáng)度曲線擬合是像素級的, 經(jīng)歷了粗細(xì)篩選的方式, 將符合灌注特征的曲線留下來. 灌注參數(shù)提取主要是提取達(dá)峰時間、峰值等. 像素級的曲線擬合與提取參數(shù)可以形成一個二維矩陣,可顏色編碼成像顯示, 也可做一些處理提取特征. 根據(jù)達(dá)峰時間二維矩陣, 做了流向情況可視化, 最后可得到流向可視化.

做流向情況可視化時, 有效坐標(biāo)占比要大于50%才能有較好的可視化效果.

未來更進(jìn)一步的實驗將是: 更好的曲線擬合公式與篩選方式, 因為目前Gamma公式可能并非對灌注情況描述最佳的公式, 雖然達(dá)峰時間誤差小但峰值處并不能很好地貼合原數(shù)據(jù); 根據(jù)參數(shù)二維矩陣更好的尋找流線的方式, 目前尋找流線的方式還需要進(jìn)行更進(jìn)一步的標(biāo)準(zhǔn); 使用更多的數(shù)據(jù), 只有數(shù)據(jù)量足夠才能找到更精確的標(biāo)準(zhǔn).