鮑 潮，姜 杉，張國彬，楊志永

(天津大學(xué)機(jī)械工程學(xué)院，天津 300350)

在我國，肺癌一直是發(fā)病率和死亡率最高的癌癥[1]．作為放射治療重要分支之一，近距離粒子放射治療由于其治療時間短、劑量控制率高等優(yōu)點被廣泛用于眾多部位癌癥的治療，例如肺部、肝部、前列腺等部位[2-4]．近距離放射治療是在醫(yī)學(xué)圖像(MRI、CT、US等)引導(dǎo)下，將所有放射性粒子通過空心穿刺針植入到病人病變組織內(nèi)，借助粒子的放射性以達(dá)到殺死癌細(xì)胞的目的．手術(shù)前要求物理師根據(jù)病人影像在治療計劃系統(tǒng)(treatment planning system，TPS)上制定出最優(yōu)的治療計劃，在術(shù)后根據(jù)植入粒子后的影像判斷治療效果[5]．

將所有粒子放置在術(shù)前計劃位置以保證最優(yōu)劑量分布是近距離手術(shù)成功的關(guān)鍵．然而由于針刺入過程中多種因素的存在，如腫瘤腫脹、組織變形、針變形等導(dǎo)致術(shù)前計劃并沒有很好地在術(shù)中執(zhí)行．這種亞優(yōu)化的劑量分布不僅會影響治療效果還會對腫瘤周圍重要器官造成威脅，因此，為了解決術(shù)前計劃失效問題，本文提出一種基于CT圖像的自動針道識別技術(shù)．有了術(shù)中實際的針道位置，術(shù)前計劃中的理想粒子便可以根據(jù)其術(shù)前所屬的針全部導(dǎo)入到術(shù)中對應(yīng)的實際針路徑上，通過調(diào)整粒子位置和數(shù)量便可以獲得一個最優(yōu)的術(shù)中劑量計劃以指導(dǎo)粒子植入．

由于術(shù)中實際針道位置對醫(yī)生來說非常重要，很多學(xué)者都對此進(jìn)行了研究．Barva等[6]利用最大化并行積分投影算法并且聯(lián)合分層網(wǎng)格加速技術(shù)實現(xiàn)了電極軸的提?。M管該方法實現(xiàn)了較高的精度，但是由于其較長的處理時間(MATLAB上運行10min左右)導(dǎo)致在臨床上并不是一個實際的手段．為了改進(jìn)上述問題，Uhercik等[7]提出一種多分辨率的并行積分投影算法，實現(xiàn)電極軸快速定位的同時且保持相同的精度．Qiu等[8]使用三維霍夫變換算法在實際病例數(shù)據(jù)中實現(xiàn)了1mm的針尖定位誤差且運行時間是2s．而后通過對霍夫變換進(jìn)行改進(jìn)，Qiu等[9]又提出在粗分辨率下進(jìn)行隨機(jī)霍夫變換，將算法運行時間減少至1s．Okazawa等[10]將霍夫變換推廣到超聲圖像中曲線針識別問題．Ding等[11]通過兩幅正交投影圖像獲取三維超聲圖像中針的位置信息．Yan等[12]還提出利用形狀信息和水平集技術(shù)也能夠成功分割出針道.

上述算法各有優(yōu)勢，但是這些算法也存在一些缺陷，例如對干擾非常敏感、計算復(fù)雜且只關(guān)注于單針的分離．為了解決現(xiàn)有算法中所存在的這些問題，本文提出一種更加穩(wěn)定、高效的針道檢測算法．該算法通過結(jié)合圖像預(yù)處理技術(shù)與改進(jìn)的隨機(jī)樣本一致性(random sample consensus，RANSAC)算法實現(xiàn)了針形狀的近似．在 RANSAC算法中，利用預(yù)檢驗技術(shù)加速了原始算法的迭代過程從而提高了針分割的效率．與此同時，利用局部 RANSAC[13]和主成分分析(principal component analysis，PCA)[14]相結(jié)合的措施保證了分割的精度．此外，通過循環(huán)淘汰最優(yōu)模型的局內(nèi)點實現(xiàn)了一張切片上多針的分離．最后，該算法通過模擬數(shù)據(jù)及病例數(shù)據(jù)得到了全面的驗證．

1 材料與方法

1.1 材 料

本文所用到的肺癌近距離粒子植入手術(shù)病例均由天津醫(yī)科大學(xué)第二附屬醫(yī)院提供．所用的 18G植入針為日本 HAKKO有限公司制造．算法是在一臺配置為Intel(R)Core(TM)i3-2120 3.30GHz的計算機(jī)上運行，使用軟件為MATLAB．

1.2 方 法

1.2.1 感興趣區(qū)域自動提取

在進(jìn)行術(shù)中針道識別步驟之前，每張CT片子上腫瘤所在區(qū)域都已經(jīng)由物理師勾畫好了，即每張切片上的腫瘤區(qū)域信息是已知的．正是基于此，才有了本文的感興趣區(qū)域(region of interesting，ROI)自動提取方法：首先根據(jù)每張切片上腫瘤邊界信息求出包裹腫瘤最小邊界的矩形，而后將矩形整體擴(kuò)大 5mm以保證所有的針道信息都在 ROI內(nèi)，如果還有未包含的針道信息則忽略．最后將擴(kuò)大后矩形所包括的區(qū)域從原圖像中剪裁下來就得到了所述的ROI．

1.2.2 圖像預(yù)處理

為了提高植入針與周圍背景組織的對比度，本文在ROI引入了一個灰度映射函數(shù)．該函數(shù)為

式中：if是轉(zhuǎn)換后像素點的灰度；i是轉(zhuǎn)換前像素點的灰度；ilow_in和ihigh_in分別是輸入圖像的最小和最大灰度值；ilow_out和ihigh_out分別是輸出圖像的最小和最大灰度值；為轉(zhuǎn)化參數(shù)．if、i、ilow_in、ilow_out、ihigh_in、ihigh_out代表的都是歸一化的灰度值，范圍在 0～1之間．該函數(shù)的意義就是把原圖像灰度進(jìn)行映射得到新圖像的過程．原圖像中灰度值低于ilow_in的像素點轉(zhuǎn)換后被賦予了灰度值ilow_out，同理灰度值高于ihigh_in的像素點在轉(zhuǎn)換后圖像中灰度值為ihigh_out．而參數(shù)定義的則是中間灰度值的轉(zhuǎn)化形式．當(dāng)≤1時，靠近ilow_in附近低灰度變化范圍被拓寬，靠近ihigh_in附近高灰度變化范圍被壓縮；同理，當(dāng)時，低灰度變化范圍被壓縮，而高灰度變化范圍則被拓寬．在本文中，ilow_in和ihigh_in分別被設(shè)置為 0.15和 1.00．輸出圖像的灰度范圍設(shè)置為[0，1]．被賦值為1.7能夠拓寬高灰度值范圍，從而增大針和周圍組織的對比度．圖1為原ROI圖像和經(jīng)過對比度增強(qiáng)后的ROI圖像．

圖1 對比度增強(qiáng)后圖像Fig.1 Contrast enhanced image

經(jīng)過對比度增強(qiáng)后能夠得到感興趣區(qū)域灰度分布直方圖，見圖 2(a)．基于閾值，原圖像I(x)?[0,1]可被分為兩個相鄰部分：Xn(針像素點)和Xb(背景像素點)，即

圖2 閾值處理后圖像Fig.2 Threshold processed image

式中：X是ROI內(nèi)所有像素點；I(x)是某一個點的灰度值．由于不同病例對應(yīng)的灰度分布不一致，因此本文的值是根據(jù)經(jīng)驗來確定的．經(jīng)過閾值處理(式(2))后得到的結(jié)果如圖 2(b)所示，灰度值高于的像素點假定為針的候選點．盡管得到了所有的針候選點，但是Xn中仍然會包含許多偽點．如何利用RANSAC算法從這些點中提取出最準(zhǔn)確針道位置如下文所述．

1.2.3 針軸與針尖識別

假設(shè)在數(shù)據(jù)集X中有N個相互關(guān)聯(lián)的數(shù)據(jù)點．RANSAC算法從X中隨機(jī)抽取一個樣本容量為m的子集，m為計算模型參數(shù)所需最小的數(shù)據(jù)子集，其值與模型本身定義有關(guān)．通過這個最小子集可以確定出一個模型，然后利用誤差閾值檢驗初始數(shù)據(jù)集X并且確定該模型對應(yīng)的所有局內(nèi)點．此抽樣和檢驗的過程將會重復(fù)直到滿足結(jié)束的標(biāo)準(zhǔn)，也就是在一定的置信概率P下至少能夠篩選出一個不包含局外點的子集所需的最小的采樣次數(shù)．最后，包含局內(nèi)點最多的模型就是最終模型．假設(shè)數(shù)據(jù)集中局內(nèi)點的占有率為，那么選出一個正確模型概率為．此時，采樣k次，k個樣本數(shù)據(jù)中至少有一個局外點概率為．因此，為了保證此概率低于以下，采樣次數(shù)k必須滿足

通過上述步驟可以得出RANSAC算法的時間復(fù)雜度T為

式中：Tx為每次抽樣所需時間；Tc為每個數(shù)據(jù)點驗證模型所用的時間；N為數(shù)據(jù)集X中總數(shù)據(jù)點數(shù)．一般來說，Tx和Tc對一特定的問題來說都為定值，因此，T的大小主要取決于k和N．由式(3)可知k與及m呈指數(shù)關(guān)系．在相同的置信概率下，k和N的值會隨著數(shù)據(jù)集及模型的復(fù)雜而增大．為了抑制T的增長，提出了預(yù)檢驗技術(shù)．

假設(shè)有n個隨機(jī)樣本點用于預(yù)檢驗，且在一定置信概率tP下，正確模型的n次實驗中至少存在nf個局內(nèi)點．一旦發(fā)現(xiàn)局內(nèi)點數(shù)少于nf，則判斷此模型為錯誤模型，不再參與后續(xù)所有數(shù)據(jù)點的檢驗．因此，對正確模型而言，通過預(yù)檢驗的概率tP為

式中為從n中抽取i個局內(nèi)點的組合種類數(shù)．在實際計算時，最優(yōu)的nf和Pt值是在最低Pt限制下，通過測試不同的nf以找到最大nf及其所對應(yīng)的Pt值而得到．?dāng)?shù)據(jù)點是否屬于局內(nèi)點則是根據(jù)該點誤差的標(biāo)準(zhǔn)偏差來判斷[15]，即

當(dāng)n≥2m時，標(biāo)準(zhǔn)偏差s可以表示為

式中di為數(shù)據(jù)點的誤差項，如圖3所示．

圖3 RANSAC算法原理示意Fig.3 Schematic diagram of RANSAC algorithm principle

在預(yù)檢驗步驟中，正確模型有可能被判定是錯誤而被剔除．因此，選出一個正確模型的概率就改變?yōu)?，且此時式(3)也變?yōu)?/p>

由此可見，預(yù)檢驗步驟對算法效率的影響轉(zhuǎn)移到了采樣次數(shù)上．當(dāng)Pt＜1時，相應(yīng)的采樣次數(shù)k要大于一般的 RANSAC算法．盡管k值有所增大，但是相比于預(yù)檢驗所淘汰的模型個數(shù)而言，隨著k值增加而增加的模型個數(shù)是少之又少的，因為大部分的模型都被局外點所污染．這一點將會在后續(xù)實驗中得到進(jìn)一步的證實．

由于前面所述停止算法的條件過于理想，錯誤的模型也有可能獲得較多的局內(nèi)點，進(jìn)而導(dǎo)致錯誤的識別．于是本文引入了局部 RANSAC和 PCA相結(jié)合的方法，以解決精度問題．通過預(yù)檢驗步驟，能夠獲得一個最初解以及與之相對應(yīng)的內(nèi)點集．所謂的局部RANSAC就是在這個亞優(yōu)化內(nèi)點集上進(jìn)行采樣和驗證新模型的內(nèi)點數(shù)．不同于一般性 RANSAC算法，局部 RANSAC采樣是在內(nèi)點集上進(jìn)行，因此不需要最小化采樣次數(shù)．一般來說，為了節(jié)省時間，局部 RANSAC采樣數(shù)kL設(shè)置為 10～20之間．經(jīng)過局部 RANSAC算法，可以得到一個優(yōu)化后的模型和數(shù)據(jù)集．然而這并不是最優(yōu)的模型，因為該模型的模型參數(shù)是通過采樣所得的最小子集計算而來的．因此，為了獲得最優(yōu)的模型參數(shù)引入了 PCA算法，以最小化數(shù)據(jù)集內(nèi)各點的誤差．

獲得最理想的針軸模型之后，針尖提取就比較容易了．方法主要是在對比度增強(qiáng)后的圖像上沿著針軸方向?qū)γ總€點灰度值進(jìn)行遍歷，找到灰度突變點則為針尖點．至此，針軸和針尖信息提取完畢．

1.2.4 多針提取

RANSAC算法一次只能提取一根針，為了獲得一張切片上所有針道，筆者提出了連續(xù)內(nèi)點刪除方法．該算法首先列出所有經(jīng)過當(dāng)前預(yù)檢驗 RANSAC算法確定的最優(yōu)模型所對應(yīng)的局內(nèi)點，然后從整體數(shù)據(jù)中移除這些被標(biāo)記的點，最后使用更新后的數(shù)據(jù)進(jìn)行下一輪的針道檢測．正是由于數(shù)據(jù)的更新，對當(dāng)前要提取的針道而言局內(nèi)點的占有率也隨之改變?yōu)?/p>

式中：nneedles是當(dāng)前切片總針數(shù)；是圖像預(yù)處理后噪聲占比；nen是針對應(yīng)的最多局內(nèi)點數(shù)；N0是初始數(shù)據(jù)量；Nj是經(jīng)過j次提取后數(shù)據(jù)點數(shù)．該連續(xù)內(nèi)點刪除方法會一直運行，直到找不到滿足最優(yōu)模型條件的樣本子集為止．至此，整個識別步驟到此結(jié)束，流程如圖4所示．

圖4 針道識別流程Fig.4 Flow chart of needle detection

2 實驗與結(jié)果

2.1 模擬數(shù)據(jù)實驗

模擬數(shù)據(jù)集是由計算機(jī)生成的 1500個隨機(jī)點，且每套數(shù)據(jù)集的局內(nèi)點均來自于一條已知的直線．不同比例(0.1～0.8)的高斯噪聲均勻地分布在指定區(qū)域．在這里給出了算法的計算時間和拾取精度，拾取精度主要指的是算法識別出的針軸和理論針軸之間的夾角．

圖5(a)為本文改進(jìn)的RANSAC算法(P-RANSAC)和一般性RANSAC算法(G-RANSAC)在局內(nèi)點占比為 0.4時拾取結(jié)果．可以看出 P-RANSAC拾取結(jié)果更接近于理論直線．從圖 5(b)可以看出 G-RANSAC算法的拾取時間隨著局內(nèi)點下降呈指數(shù)升高，而PRANSAC算法運行時間卻只是發(fā)生了略微的增長．對于角度偏差，從圖 5(c)可以看出P-RANSAC不論在何種局內(nèi)點占比下都要小于 G-RANSAC，且精度一直控制在0.3°以內(nèi)．

圖5 模擬數(shù)據(jù)實驗結(jié)果Fig.5 Experimental results of simulation data

2.2 病例實驗

為了驗證針道識別算法在實際圖像環(huán)境中的性能，本文通過 12組肺癌病例加以驗證，圖 6為病例實驗結(jié)果．所有病例切片層厚均為 5mm，用到的植入針均為 18G的八光針．實驗中以手動分割作為參照的金標(biāo)準(zhǔn)，且識別算法在每套病例中均進(jìn)行 15次以獲得均值作為記錄．

圖6 病例數(shù)據(jù)實驗結(jié)果Fig.6 Experimental results of patient data

圖 6統(tǒng)計了所有植入針的信息識別情況．分析圖 6(a)數(shù)據(jù)可以看出，超過 95%的針尖拾取誤差小于 0.8mm，且最大針尖誤差不超過 0.85mm．分析圖6(b)數(shù)據(jù)可以看出，超過 97%的角度偏差小于 1.3°，且所有角度偏差都控制在1.4°以內(nèi)．由此可見本文提出的針道識別算法的拾取精度完全滿足臨床要求．

此外，分析了算法的分割效率．從圖 6(c)數(shù)據(jù)可以看出算法平均拾取每根針時間是 0.236s，其他各組病例的平均針道識別時間在這一值下來回波動．經(jīng)統(tǒng)計，12組病例共植入了八光針 153根，算法成功拾取了145根，拾取率為94.78%．此外，本文對幾組病例中針道識別失敗的原因也進(jìn)行了仔細(xì)分析．主要分為兩個方面：一是由于部分病例背景像素(如肋骨等)干擾嚴(yán)重以及存在大量針道偽影；二是由于有些針只有很少一部分刺入了腫瘤，導(dǎo)致只有有限的針候選點可用于針道識別算法，這種情況經(jīng)常發(fā)生在腫瘤邊緣部位的針道．圖 7為某病例最后拾取出的結(jié)果．

圖7 某病例針道拾取效果Fig.7 Example of needle segmentation result in a patient’s case

3 討 論

本文針對實際臨床需求及相關(guān)研究現(xiàn)狀，提出了一種基于 CT圖像的、自動高效的針道識別算法，并通過模擬和實際病例實驗驗證了算法的性能．對于實際的病例實驗，所有的針尖定位誤差和針軸方向誤差均被控制在0.85mm和1.4°，這種亞毫米級的精度完全滿足臨床近距離手術(shù)要求．而且相比于文獻(xiàn)[8]實際病例針尖拾取精度 1mm，本文算法拾取結(jié)果要更優(yōu)．此外，就針軸方向誤差而言，本文實驗結(jié)果要小于 2.3°(文獻(xiàn)[16])．與此同時，算法識別整體失敗率為5.22%，這對于一個隨機(jī)迭代算法(RANSAC)而言是完全合理的．分割效率方面，本文提出的算法平均每根針分割時間為 0.236s，這相比于文獻(xiàn)[7]中7.3s及文獻(xiàn)[9]中 1.0s都要短很多，這足以說明本文算法效率是非常高的．

此外，對于肺部近距離粒子植入手術(shù)而言，在粒子植入前實現(xiàn)術(shù)前計劃的再計劃以達(dá)到最優(yōu)的術(shù)中劑量分布是非常關(guān)鍵的．因此，本文所提出的基于CT圖像的、快速精確的針道識別算法為劑量的修正提供了新思路．因為一旦所有針的實際位置都得知以后，物理師就能夠?qū)⑿g(shù)前計劃導(dǎo)入到術(shù)中，而后通過調(diào)整最初粒子位置或添加針以及粒子就可以改善術(shù)中劑量分布．另一方面，自動識別出的針軸方向和針尖位置可以作為一種圖像測量工具被醫(yī)生所使用．一旦發(fā)現(xiàn)針的植入位置出現(xiàn)較大錯誤，醫(yī)生便可立即采取相應(yīng)的措施進(jìn)行補(bǔ)救．最后，本文所提出的算法還具有很高的靈活性，只要重新定義模型函數(shù)及其相應(yīng)的誤差項便可用于其他參數(shù)化目標(biāo)的識別．

4 結(jié) 語

本文提出了一種改進(jìn)的 RANSAC算法用于 CT圖像的針道識別，該算法首先利用預(yù)檢驗技術(shù)加速了迭代過程，而后結(jié)合局部RANSAC和PCA算法保證了較高的精度，最后通過模擬和病例數(shù)據(jù)實驗驗證了該算法的性能．與此同時，所得的實驗結(jié)果(較短的分割時間和亞毫米級的精度)滿足臨床手術(shù)的要求.在未來工作中，筆者希望一個基于自動針道識別的可靠的術(shù)中劑量優(yōu)化模塊能夠整合到近距離放射治療計劃系統(tǒng)中去．

[1] 2017年中國最新癌癥數(shù)據(jù)[J].中國腫瘤臨床與康復(fù)，2017，24(6)：760.

Latest cancer statistics in China，2017[J].Chinese Journal of Clinical Oncology and Rehabilitation，2017，24(6)：760(in Chinese).

[2] 李任飛，王月東，閆 龑，等.125I粒子植入治療非小細(xì)胞肺癌近期療效評估[J].介入放射學(xué)雜志，2014，23(1)：65-68.

Li Renfei，Wang Yuedong，Yan Yan，et al.Implantation of125I seeds for the treatment of non-small cell lung cancer：Evaluation of short-term[J].Journal of Interventional Radiology，2014，23(1)：65-68(in Chinese).

[3] 黃文藪，蔡明岳，曾昭吝，等.Tace聯(lián)合125I放射性粒子植入治療肝細(xì)胞癌門靜脈癌栓[J].介入放射學(xué)雜志，2015，24(6)：488-493.

Huang Wensou，Cai Mingyue，Zeng Zhaolin，et al.Transarterial chemoembolization combined with CT-guided125I seed implantation for the treatment of hepatocellular carcinoma associated with portal vein tumor thrombus[J].Journal of Interventional Radiology，2015，24(6)：488-493(in Chinese).

[4] 張麗霞，陳金燕，徐彩云，等.125I粒子植入治療高風(fēng)險局限性前列腺癌臨床研究[J].浙江醫(yī)學(xué)，2017，39(24)：2267-2269.

Zhang Lixia，Chen Jinyan，Xu Caiyun，et al.Efficacy of 125-iodine seeds branchytherapy for high-risk local prostate cancer[J].Zhejiang Medical Journal，2017，39(24)：2267-2269(in Chinese).

[5] Stewart A，Parashar B，Patel M，et al.American brachytherapy society consensus guidelines for thoracic brachytherapy for lung cancer[J].Brachytherapy，2016，15(1)：1-11.

[6] Barva M，Uhercik M，Mari J M，et al.Parallel integral projection transform for straight electrode localization in 3-D ultrasound images[J].IEEE Transaction on Ultrasonics，F(xiàn)erroelectrics，and Frequency Control，2008，55(7)：1559-1569.

[7] Uhercik M，Kybic J，Liebgott H，et al.Multiresolution parallel integral projection for fast localization of a straight electrode in 3D ultrasound images[C]//5th IEEE International Symposium on Biomedical Imaging.USA：IEEE，2008：33-36.

[8] Qiu W，Yuchi M，Ding M，et al.Needle segmentation using 3D Hough transform in 3D TRUS guided prostate transperineal therapy[J].Med Phys，2013，40(4)：305.

[9] Qiu W，Ding M，Yuchi M.Needle segmentation using 3D quick randomized Hough transform[C]// Proc SPIE Intell Netw Intell Syst(ICNIS).USA：IEEE，2008：449-452.

[10] Okazawa S H，Ebrahimi R，Chuang J，et al.Methods for segmenting curved needles in ultrasound images[J].Med Image Anal，2006，10(3)：330-342.

[11] Ding M，Cardinal H N，F(xiàn)enster A.Automatic needle segmentation in three-dimensional ultrasound images using two orthogonal two-dimensional image projections[J].Med Phys，2003，30(2)：222-234.

[12] Yan P，Cheeseborough J C，Chao K S C.Automatic shape-based level set segmentation for needle tracking in 3-D TRUS-guided prostate brachytherapy[J].Ultrasound Med Biol，2012，38(9)：1626-1636.

[13] Chum O，Matas J，Kittler J.Locally optimized ransac[C]// Lecture Notes in Computer Science.London，UK，2003：236-243.

[14] 許淑娜，李長坡.對主成分分析法三個問題的剖析[J].數(shù)學(xué)理論與應(yīng)用，2011，31(4)：116-121.

Xu Shuna，Li Changpo.Dissection to three typical issues of principal component analysis[J].Mathematical Theory and Applications，2011，31(4)：116-121(in Chinese).

[15] Torr P H S，Zisserman A.MLESAC：A new robust estimator with application to estimating image geometry[J].Comput Vis Image Und，2000，78(1)：138-156.

[16] Cool D W，Gardi L，Romagnoli C，et al.Temporalbased needle segmentation algorithm for transrectal ultrasound prostate biopsy procedures[J].Med Phys，2010，37(4)：1660-1673.