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

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

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

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

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

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

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

1.1 研究現(xiàn)狀

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

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

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

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

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

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

圖1 雙模式視頻幀

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

圖2 提取mask掩膜示意圖

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

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

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

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

第3步. 提取灌注特征參數(shù). 因?yàn)槭腔谙袼氐姆治? 所以對(duì)符合標(biāo)準(zhǔn)的曲線(xiàn)提取灌注參數(shù)(達(dá)峰時(shí)間、峰值等), 形成二維矩陣. 對(duì)二維矩陣可以編碼成像顯示, 輔助醫(yī)生診斷病情.

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

2 病灶跟蹤

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

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

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

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

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

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

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

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

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

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

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

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

2.3 卡爾曼濾波

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

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

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

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

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

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

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

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

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

3.1 二維平滑窗口的選擇

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

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

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

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

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

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

3.3 基于Gamma公式的曲線(xiàn)擬合

Gamma模型公式[10]為:

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

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

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

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

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

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

第1步, 粗篩選.

粗篩選大致分兩個(gè)思路, 兩個(gè)子篩選條件都滿(mǎn)足了才算通過(guò)粗篩選. 一個(gè)是真實(shí)值的判斷, 一個(gè)是縮放值的判斷.

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

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

圖6 真實(shí)值判斷示意圖

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

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

圖7 縮放值判斷示意圖

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

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

因?yàn)椤癛-squared”一般應(yīng)用于線(xiàn)性擬合中的擬合度評(píng)價(jià)準(zhǔn)則, 所以沒(méi)有使用在這里[16]. 本文的擬合度的評(píng)價(jià)準(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合在一起共同的平均值. 使用修正余弦相似性有一個(gè)重要條件, 需要兩個(gè)時(shí)間序列一樣長(zhǎng).

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

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

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

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

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

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

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

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

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

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

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

圖8 TTP二維矩陣成像圖

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

圖9 AUC二維矩陣成像圖

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

5.1 灌注方向的求法

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

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

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

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

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

5.3 每條流線(xiàn)的生成

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

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

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

class Point{

int x;

int y;

};

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

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

圖10 尋找流線(xiàn)示意圖

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

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

(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)加入一個(gè)點(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é)束; 否則就將其周?chē)c(diǎn)加入pointGroup點(diǎn)群, 重復(fù)本步驟. 另外, 此步驟中pointGroup點(diǎn)群的長(zhǎng)度大于限制值或ROI內(nèi)所有點(diǎn)都已經(jīng)遍歷也會(huì)終止.

5.4 流線(xiàn)顯示

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

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

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

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

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

6 生成向量場(chǎng)并可視化

6.1 生成向量場(chǎng)

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

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

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

6.2 線(xiàn)積分卷積

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

其生成過(guò)程需要輸入向量場(chǎng)和一張?jiān)肼晥D(通常是白噪聲)[11], 如圖12所示.

圖12 LIC圖生成示意圖

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

6.3 灌注流向的向量場(chǎng)成像

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

如圖13, 是對(duì)上面流線(xiàn)圖的例子進(jìn)行計(jì)算的結(jié)果.其紋理顯示了向量場(chǎng)的大致?tīng)顟B(tài).

圖13 LIC結(jié)果圖

7 實(shí)驗(yàn)結(jié)果分析

7.1 數(shù)據(jù)

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

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

如表1所示, 18例視頻都進(jìn)行了實(shí)驗(yàn). 有效坐標(biāo)點(diǎn)是上文曲線(xià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é)果

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

因?yàn)楸疚牡闹攸c(diǎn)是探尋灌注流向信息可視方法,所以實(shí)驗(yà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é)果. 其中, 每一行圖片為一個(gè)視頻的結(jié)果, 從上而下有效坐標(biāo)占比依次為“0–20%”“20%–50%”“50%–80%”“80%–100%”. 每一行從左到右分別為T(mén)TP成像圖、流線(xiàn)圖、流線(xiàn)技術(shù)圖(log增強(qiáng))、LIC結(jié)果圖.

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

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

7.4 總的結(jié)果分析

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

8 結(jié)論與展望

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

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

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

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