基于改進(jìn)的Live-Wire算法在ARPlanner中的應(yīng)用

汪欣，康世功，郎錦義

1. 北京全域醫(yī)療技術(shù)集團(tuán)有限公司，北京 100013；2. 美中腫瘤診療技術(shù)創(chuàng)新研究院，北京 100013 3. 四川省腫瘤醫(yī)院，四川 成都 610041

引言

患者進(jìn)行放療之前，需要對患者制定放療計劃，而精準(zhǔn)的靶區(qū)及危及器官的勾畫是精準(zhǔn)放療的前提[1]。在精準(zhǔn)勾畫前提下，靶區(qū)附近危及器官受照劑量才會更精確，從而提高治療的精度，降低副反應(yīng)發(fā)生率。醫(yī)生在勾畫靶區(qū)以及危及器官時需要反復(fù)修改，因而花費(fèi)了大量時間，降低了勾畫靶區(qū)的工作效率。一些危及器官雖然經(jīng)過醫(yī)生的反復(fù)修改，但是勾畫的精度仍遠(yuǎn)遠(yuǎn)不夠。

目前，靶區(qū)的勾畫主要分為三種[2-4]：自動勾畫、人機(jī)交互式半自動勾畫、手動勾畫。由于醫(yī)學(xué)圖像的復(fù)雜性，自動勾畫難以得到準(zhǔn)確的器官勾畫結(jié)果，容易造成輪廓的偏離；手動勾畫需要消耗醫(yī)生大量的時間與精力，隨著CT設(shè)備的不斷提升，層數(shù)數(shù)量也在不斷的增加，這給醫(yī)生手動勾畫帶來了大量的工作，勾畫一個患者往往需要數(shù)個小時甚至幾天，需要不斷的修改，由于耗時長，精準(zhǔn)度又受主觀的影響較大，極大地降低了醫(yī)生的工作效率；人機(jī)交互式勾畫，可以通過醫(yī)生的經(jīng)驗配合圖像處理的技術(shù)，因此本文采用Live-Wire算法對患者的靶區(qū)以及危及器官進(jìn)行交互式勾畫，可以極大地提高醫(yī)生的工作效率與勾畫的精準(zhǔn)度。

1 方法

Live-Wire[5-7]首先是快速地定位到該圖像區(qū)域的邊緣上，通過與用戶交互，更精確地完成整個感興趣區(qū)域的勾畫、分割。Live-Wire算法的實現(xiàn)通過兩步，第一步是對構(gòu)造成本函數(shù)的計算，第二步是對最短路徑的搜索。首先將待分割的圖像映射成一個加權(quán)有向圖，圖像中的像素點對應(yīng)有向圖中的結(jié)點。連接圖像中的相鄰像素的邊對應(yīng)于連接圖中相鄰結(jié)點的邊，并且在邊上定義構(gòu)造函數(shù)，使得圖像中的器官邊緣具有較小的權(quán)重。給予非邊緣區(qū)域更大的權(quán)重，然后根據(jù)圖像的特征計算每個像素的Laplace過零點特征函數(shù)、梯度特征函數(shù)、梯度方向函數(shù)局部代價，然后，基于加權(quán)有向圖尋路算法Dijkstra執(zhí)行最短路徑搜索，并根據(jù)初始種子點和目標(biāo)點搜索最短路徑，獲得感興趣區(qū)域的邊緣。傳統(tǒng)的Live-Wire算法的梯度幅值由水平方向和垂直方向來計算。本文對Live-Wire算法進(jìn)行了改進(jìn)，在梯度幅值[8-10]的計算由原算法中的水平方向和垂直方向改進(jìn)為由水平方向、45°方向、垂直方向和135°方向來計算。

1.1 改進(jìn)的構(gòu)造成本函數(shù)

每個像素的構(gòu)造成本函數(shù)[11]計算公式如(1)所示。

其中，公式(1)中函數(shù)l(p,q)表示像素p與相鄰像素q之間的成本函數(shù)；函數(shù)fZ(q)表示像素q處的Laplace零交叉點特征函數(shù)；函數(shù)fG(q)表示像素q的梯度特征函數(shù)；函數(shù)fD(p,q)表示點(p,q)的梯度方向。經(jīng)過多次實驗，權(quán)值wZ=0.43,wG=0.43,wD=0.14。

1.1.1 Laplace過零點特征函數(shù)fZ(q)的計算

Laplace過零點特征函數(shù)[11-12]是一種二值化的邊緣檢測方法[13-16]。因為圖像邊緣的灰度變化較大，所以圖像的一階偏導(dǎo)數(shù)在邊緣處有局部最大值或最小值，二階偏導(dǎo)數(shù)在邊緣處會通過“0”點。過“0”點指兩個相鄰的像素二階偏導(dǎo)數(shù)由正變?yōu)樨?fù)，或者由負(fù)變?yōu)檎狞c。這一正一負(fù)的兩個像素點，值與“0”接近的點稱為過“0”點，取值0。Laplace過零點特征函數(shù)計算得到的結(jié)果是一個像素寬的窄帶，表示圖像內(nèi)目標(biāo)的邊緣。Laplace過零點特征函數(shù)中的過“0”點表示具有較強(qiáng)的邊緣特征的點。假設(shè)IL(q)是圖像I在q點處的函數(shù)值，則其公式如下：

其中，函數(shù)IL(q)是像素點q經(jīng)過Laplace變換后的值，函數(shù)fZ(q)的作用是突出醫(yī)學(xué)圖像中器官的邊緣。

1.1.2 改進(jìn)的梯度特征函數(shù)fG(q)的計算

圖像使用Laplace過零點特征函數(shù)計算得到的是一個二值矩陣，計算的結(jié)果只有極少的點值是0。由于計算得到的是一個二值矩陣，因此它僅僅能表示像素點是否屬于目標(biāo)邊緣而不能說明像素點屬于圖像邊緣概率的大小、梯度大小及邊緣特征的強(qiáng)弱。因此在公式(1)中引入了一個表示像素點圖像邊緣特性強(qiáng)度的梯度特征函數(shù)[11,17]。有很多種方法的算子可以計算圖像的梯度值，在本論文中采用的是高斯函數(shù)。由于在一副圖像中每個像素都具有八個相鄰的像素點，當(dāng)邊緣與檢測方向垂直才能很好地檢測到邊緣點，因此在上述方法中只用到水平方向和垂直方向來計算梯度的幅值在一定條件下存在漏檢情況，為了減少漏檢，本文在計算梯度幅值方法上做了改進(jìn)，將用水平方向Ix、45°方向I45°、垂直方向Iy、135°方向I135°作為梯度幅值的計算方向，以下為四個方向的算子。

梯度G的計算公式(4)如下：

梯度值越大說明圖像的邊緣特征越強(qiáng)，那么像素點的梯度特征函數(shù)得到的權(quán)值就應(yīng)該越小，

為使高梯度產(chǎn)生低能量，令

其中G′=G-min(G)，為了保證梯度特征最大值計算方式的統(tǒng)一，梯度特征函數(shù)值需要使用歐幾里得距離進(jìn)行處理。如果q是p的斜向的相鄰點，fG(q)需要除以；如果q是p的水平或垂直方向的相鄰點，fG(q)需要除以1，如公式(6)所示。

1.1.3 梯度方向fD(p,q)的計算

在邊緣提取中，保證圖像邊緣的平滑是必要的。在公式(1)中，梯度方向函數(shù)[17]就是用來控制目標(biāo)邊緣平滑的平滑約束因子。梯度方向函數(shù)給邊緣方向變化較大的點賦予一個較高的權(quán)值，給邊緣變化平滑的點賦予較小的權(quán)值。梯度方向值是一個由Ix和Iy決定的單位向量。假設(shè)D(p)是同p點處梯度方向垂直的單位向量，梯度方向函數(shù)的公式是：

其中

從公式(7)可以看出，如果與相鄰點的梯度方向具有相同或者相似梯度方向的點具有較小的權(quán)值，而與相鄰點梯度方向垂直的點具有較大的權(quán)值。其中公式(9)中的Ix(p)，Iy(p)，Ix(q)，Iy(q)分別表示像素p、像素q沿方向x、方向y的方向梯度。

由上述三個計算成本函數(shù)的公式可以方便的找到圖像的邊緣點。

1.2 計算最小成本路徑

當(dāng)成本函數(shù)確定之后，下一步就是使用加權(quán)有向圖的路徑尋找算法來尋找最佳路徑。算法的思想：只需每次設(shè)置種子點，然后重新計算從種子點到其他結(jié)點的最短路徑。也就是說，所有路徑的集合構(gòu)成從種子點擴(kuò)散的最短路徑樹。在本文中，Dijkstra算法[18-21]用于找到最短路徑，并且可以找到從單個源點到其他頂點的最短路徑。

Dijkstra算法的輸入為：有向圖G和G中的源頂點S，用V來表示G中所有頂點的集合。用E來表示G中連接結(jié)點邊的集合。(u,v)表示從頂點u到v的路徑是連通的，邊緣的權(quán)重由權(quán)重函數(shù)w:E→[0,∞]定義，因此，w(u,v)是頂點u到頂點v的非負(fù)成本（cost），兩個頂點之間的距離可用于計算邊的成本。尋路算法還可以找到從一個頂點到圖中任何其他頂點的最短路徑。如圖1所示，結(jié)點A、B、C、D、E分別為平面上的5個結(jié)點，現(xiàn)在需要分別求出結(jié)點A到其他4個結(jié)點B、C、D、E的最短路徑，優(yōu)先隊列Q={A、B、C、D、E}。

圖1 有向圖G及各邊權(quán)重

尋路過程如下：

（1）設(shè)起點的距離為0，即鍵值為0，其他的點為∞，因為此時不知道如何到達(dá)，即鍵值為∞，從隊列中取出最小值，取出隊列中的最小值A(chǔ)，然后觀察A，把A加到集合S中，A在隊列里移出來，就不會再放回去了。

（2）接下來更新隊列中其他結(jié)點的鍵值，觀察跟A有邊連接的點，從A到B有一條，它的權(quán)值是10，沿AB邊走，它的花費(fèi)是10，總權(quán)值是0+10=10，更新隊列里B的權(quán)值為10，從A到C是另外一條邊，0+3=3，更新C的鍵值3，然后在隊列里找出最小值3即C點，說明這里的最短路徑是：從A到C，從隊列里刪除C。

（3）然后看所有從C出來的邊，有一條是到B的，權(quán)值為4+3=7，比之前的10小，所以從這條路徑到B會更好，更新B的鍵值為7，這里還有一條從C到D的，費(fèi)用為8的邊加上3為11，更新D的鍵值為11，然后還有一條從C到E的，費(fèi)用為3+2=5，更新E的鍵值為5，到現(xiàn)在所有頂點都有一條有限的路徑相連了。從這個隊列里找出最小值E，此時這是最短路徑，到達(dá)E的方案是從A->C->E，稱這是最短路徑。

（4）然后看從E點出發(fā)的邊，這里只有一條，從E到D，費(fèi)用是5+9=14，比之前的11要大，所以不走這條路徑。

（5）現(xiàn)在剩下B和D，A、C、E都被刪除了，取出最小值B的鍵值為7，這里有一條路徑B到C，費(fèi)用是7+1=8，要比3大，接著B到D，費(fèi)用為7+2=9，比11要小，此時我們找到了一條更短的路徑，所以當(dāng)前D的最短路徑的權(quán)值為9，即這條路徑是經(jīng)過A->C->B->D，移除B點，現(xiàn)在只剩D點了。

（6）觀察從D點出發(fā)的邊，這里只有一條，從D到E，費(fèi)用為9+7=16，比之前的5要大，所以E的最短路徑已經(jīng)確定了，現(xiàn)在刪除D點，隊列是空的。

其尋路結(jié)果，見表1。

表1 A點到B、C、D、E各點的最短路徑

在本文中使用的算法是Dijkstra算法的一種變體，輸入包括： 種子點，有向圖。輸出為： 輸入圖中的最小路徑樹，每一個結(jié)點指向它的前一個結(jié)點，沿著種子點到該結(jié)點的最小代價路徑。還將為每個結(jié)點分配總成本，對應(yīng)于從結(jié)點到種子的最小成本路徑。在有向圖中，每一個結(jié)點都有以下三種狀態(tài)：初始狀態(tài)、活動狀態(tài)、擴(kuò)展?fàn)顟B(tài)。當(dāng)全部結(jié)點都被擴(kuò)展時，本算法終止。當(dāng)算法開始時，圖中的全部結(jié)點進(jìn)行初始化。當(dāng)算法運(yùn)行時，所有活動結(jié)點都保持在優(yōu)先隊列Q中，該優(yōu)先隊列根據(jù)結(jié)點到種子的當(dāng)前總成本排序。本文的變體的Dijkstra算法流程，見圖2。

圖2 本文變體的Dijkstra算法流程圖

2 結(jié)果

本文在驗證算法過程中在200例不同癌種患者中隨機(jī)抽取了患者的CT圖像序列，實驗效果表明改進(jìn)后的勾畫效果比原算法勾畫效果更能滿足臨床要求，本文采取了4幅典型圖像進(jìn)行對比顯示，見圖3～6。

本文算法在原有算法的基礎(chǔ)上對梯度幅值的計算進(jìn)行了改進(jìn)，從水平方向和垂直方向來計算梯度的幅值改進(jìn)為從水平方向、45°方向、垂直方向、135°方向計算梯度幅值，添加醫(yī)學(xué)圖像在多個梯度方向都很重要，而原方法只能提取2個方向梯度，對對角線方向梯度不能有效提取，這將導(dǎo)致在邊緣提取的時候無法實現(xiàn)精準(zhǔn)提取，而本文中方法明顯對對角線方向?qū)崿F(xiàn)梯度提取，兼顧了多個方向，邊緣提取的精度得到有效保障。改進(jìn)之后的算法能夠使邊緣更精確的被檢測出來，本文對改進(jìn)的Live-Wire算法和原Live-Wire算法的效果進(jìn)行了對比，通過對比圖可以發(fā)現(xiàn)，改進(jìn)的Live-Wire算法能夠更精確的器官的邊緣進(jìn)行檢測。本文的算法能夠使靶區(qū)附近危及器官受照劑量更精確，提高治療的精度，降低副反應(yīng)發(fā)生率，提高了醫(yī)生勾靶的工作效率，為臨床勾畫工作降低了時間成本。

圖3 肺部勾畫效果示例

圖4 肝臟勾畫效果示例

圖5 腎臟勾畫效果示例

圖6 腦部勾畫效果示例

將本文改進(jìn)的交互式自動分割算法Live-Wire應(yīng)用到ARPlanner軟件（ARPlanner軟件為北京全域醫(yī)療技術(shù)集團(tuán)有限公司開發(fā)的放射治療計劃輪廓勾畫軟件，與醫(yī)用加速器配套，適用于腫瘤患者的放射治療。支持患者CT/MRI圖像輸入，并作三維圖像重建與顯示；支持任意切面解剖圖像的顯示；支持PET/CT、CT/MR、增強(qiáng)/平掃圖像之間的融合配準(zhǔn)；支持手動勾畫、半自動勾畫、自動勾畫靶區(qū)和器官輪廓；具有DICOM網(wǎng)關(guān)功能，支持設(shè)備間數(shù)據(jù)通訊）中，可以看到通過本文算法能夠更精確的勾畫出器官（表2）。通過與Monaco軟件（Monaco是由Elekta開發(fā)生產(chǎn)的放射治療計劃軟件，主要用于常規(guī)方法的輪廓繪制、影像處理、模擬、影像融合、計劃優(yōu)化等）半自動勾畫功能進(jìn)行對比，發(fā)現(xiàn)本文的算法勾畫精確度要高于Monaco中勾畫的精度。從而說明本文算法具有一定的臨床應(yīng)用價值。

表2 手動勾畫與本文算法勾畫以及Monaco中邊緣檢測勾畫對比表

3 結(jié)論

本論文將改進(jìn)的Live-Wire算法運(yùn)用到實際的放射治療計劃輪廓勾畫軟件中，極大地提高了腫瘤醫(yī)學(xué)勾畫的精準(zhǔn)度與速度，為醫(yī)生勾畫靶區(qū)以及危及器官提高了工作效率以及精準(zhǔn)度，為后續(xù)做計劃提供了精準(zhǔn)的前提，從而提高了放療的精準(zhǔn)度。本文的主要貢獻(xiàn)是：一是引入水平方向、45°方向、垂直方向、135°方向計算梯度幅值，添加醫(yī)學(xué)圖像在多個梯度方向都很重要，而原方法只能提取2個方向梯度，對對角線方向梯度不能有效提取，這將導(dǎo)致在邊緣提取的時候無法實現(xiàn)精準(zhǔn)提取，而本文中方法明顯對對角線方向?qū)崿F(xiàn)梯度提取，兼顧了多個方向，邊緣提取的精度得到有效保障。二是根據(jù)Dijkstra原始算法進(jìn)行變體，使得該算法在確定初始種子點以后，可以靈活自由的進(jìn)行路徑選擇。與當(dāng)前應(yīng)用比較廣泛的放射治療軟件Monaco中的半自動勾畫功能進(jìn)行了對比，發(fā)現(xiàn)本文算法的勾畫效果要優(yōu)于Monaco中的半自動勾畫算法勾畫的效果，有一定的臨床應(yīng)用價值。