肺部四維計算機斷層掃描(4D-CT)能準確定位腫瘤靶區, 對于實現肺癌的精確放療有重要意義。但4D-CT圖像數量巨大, 單純靠人工分割腫瘤勾畫靶區很不現實, 而且肺4D-CT圖像復雜, 目前的自動分割技術難以保證分割精度。針對這些問題, 本文提出了一種基于星形先驗和圖割的肺4D-CT腫瘤自動分割技術。首先, 在初始相位圖像上選取目標種子點, 以此點為中心, 形成初始目標塊。其次, 采用運動估計中的完全搜索塊匹配算法, 獲得下一相位圖像中與初始目標塊最相似的目標塊, 同時計算出它們之間的運動位移, 以此類推, 得到所有相位的目標塊以及對應塊之間的運動位移。隨后, 利用這些運動位移和初始相位目標種子點的位置, 計算出其余各相位目標種子點, 將它們作為星形先驗的中心點。最后, 在各相位的目標塊上使用結合星形先驗的圖割算法, 得到腫瘤分割結果。視覺和量化兩方面評價結果均表明, 相較于傳統圖割算法, 本文提出的方法不僅能夠提高分割的準確性, 而且大幅提升了自動化程度。
引用本文： 申正文, 高圓圓, 張煜. 基于星形先驗和圖割的肺部四維計算機斷層掃描腫瘤自動分割. 生物醫學工程學雜志, 2016, 33(2): 295-302. doi: 10.7507/1001-5515.20160050 復制
引言
肺癌在20世紀末已經成為發病率和死亡率增長最快、對人群健康和生命威脅最大的惡性腫瘤之一[1]。雖然目前對于腫瘤的治療手段越來越多，但手術、放療和化療依然是三種最有效的主要治療手段，放射治療適應癥寬、療效較好，在腫瘤的治療中有著無可置疑的重要地位[2]。因為獨特的空間分辨率和密度分辨率優勢，計算機斷層掃描(computed tomography, CT)在放射治療中扮演著重要的角色。在肺部腫瘤的放射治療中，指導放療計劃和劑量投射的信息通常都來自于肺部CT圖像。近年來，隨著影像技術的進步，將時間因素納入三維計算機斷層掃描(three dimensional computed tomography, 3D-CT)的掃描和重建過程，就出現了四維計算機斷層掃描(four dimensional computed tomography, 4D-CT)。4D-CT圖像與放療技術相結合，對運動靶區的精確放療產生了深遠的影響。利用4D-CT圖像分析靶區及正常器官在呼吸過程中位置和體積的變化規律，讓放射線的劑量最大限度地集中到靶區內，可減少周圍正常組織和器官受到的照射，達到保護正常組織與器官的目的[3]。
放射治療中靶區的精確定位至關重要，靶區的精度會直接影響最終的放療結果，確定靶區可以采用圖像分割技術。目前肺4D-CT腫瘤分割存在兩大主要問題：①4D-CT生成的圖像數量巨大，通常有1 000~2 000張，甚至更多，如果僅靠放療科醫生手工分割每幅圖像中的腫瘤來獲得靶區，顯然費時費力，也不現實[4]；②肺部腫瘤經常與周圍肺壁、膈肌、血管等正常器官組織發生黏連，若采用自動分割技術，誤分割的概率很大，直接影響靶區精度。因此研究出一種自動化程度高、分割準確、魯棒性強的肺腫瘤分割算法十分必要。
自從Boykov和Jolly于2000年提出圖割算法以來[5]，圖割算法憑借其原理簡單、分割精確、全局最優以及能融入先驗知識并且可用于高維分割等優點，一直受到各方重視，成為圖像分割領域中的研究熱點。圖割算法是一種交互式的分割算法，用戶先手動選擇目標種子點和背景種子點，而后交由計算機自動分割出目標物體。相比于其他全自動分割算法，圖割算法因其具有交互性，在分割的準確性上具有明顯優勢，但與此同時也帶來了不可避免的缺陷，分割結果對選取的種子點十分敏感，與用戶輸入息息相關。如果直接將圖割算法用于放療計劃中分割腫瘤，放療科醫生必須精確選取種子點才能確保分割結果良好，這將意味著增加本就十分繁重的醫生工作。
提升圖割算法自動化程度的關鍵在于自動獲取種子點，考慮到腫瘤形態和位置會隨呼吸運動發生改變，并且4D-CT相鄰相位圖像之間存在一定的相似性，本文采用運動估計來自動獲取種子點。同時為了保證分割精度，在傳統圖割算法中加入星形先驗，獲得的種子點就作為星形中心點。具體過程如下：①在4D-CT初始相位的三維圖像上，由醫生選取目標種子點，以此種子點為中心，形成一個初始目標塊，腫瘤及其周邊皆被包含在其中[6]。②采用運動估計中常用的完全搜索塊匹配算法[7]，獲得下一相位圖像中與初始目標塊最相似的目標塊，同時估計出它們之間的運動位移，以此類推，可以得到所有相位的目標塊以及對應塊之間的運動位移。③利用這些運動位移和初始相位目標種子點的位置，計算出其余各相位目標種子點，作為星形先驗的中心點。④在各相位的目標塊上使用結合星形先驗的圖割算法，即可得到腫瘤分割結果。
本文提出的方法只需醫生在初始相位上選取目標種子點，即可獲得同一呼吸周期內所有相位的腫瘤分割結果。不僅大幅提高了圖割算法的自動化程度，降低了放療科醫生的工作量，而且從量化評價指標來看，能夠確保分割結果良好，有效解決了前述目前放療中肺4D-CT腫瘤分割存在的兩大主要問題。
1 理論與方法
1.1 圖割理論
1.1.1 傳統圖割理論
圖割是一種基于圖論的組合優化方法，它將一幅圖像映射成一個網絡圖，其中像素點對應圖中的節點，相鄰像素點之間的關系對應圖中節點之間的邊，在此基礎上建立能量函數，運用最大流/最小割算法對網絡圖進行切割，得到網絡圖的最小割，即原圖的分割結果。圖割理論不僅適用于二維圖像分割，還適用于高維分割，并且基本原理不變。
圖割理論詳細論述參見文獻[8-9]，現概括如下，設G=(V, E)是一個帶有權值的無向圖。其中V為點集，對應圖像中各像素點p的集合P。V還包含兩個額外添加的頂點，稱其為終端，一個為源點，另一個為匯點。E為邊集，對應圖像像素點之間的關系，E包含兩類邊，一類是相鄰像素{p, q}之間的邊，另一類是像素點p與終端之間的邊。給圖中每一條邊賦予權值，設定權值的原則為：像素間的差異越小則權值越大，像素間差異越大則權值越小。用N代表像素點集合P中的4鄰域系統，用1和0分別作為目標和背景的標記，fp∈{1, 0}表示分配給像素點p的標記，由此可定義能量函數：
| $ E\left(f \right)=\sum\limits_{p \in P} {{D_p}\left({{f_p}} \right) + \lambda \sum\limits_{\left({p, q} \right) \in N} {{V_{pq}}\left({{f_p}, {f_q}} \right)} } $ |
能量函數式(1)中右邊第一項為區域項，Dp(fp)是賦予像素點p標記fp的代價測度；右邊第二項為邊界項，Vpq(fp, fq)是賦予相鄰像素點p、q標記fp、fq的代價測度。λ為非負系數，決定了區域項和邊界項所占的比重，本文實驗中λ=300[10]。求E(f)的最小值利用最大流/最小割算法完成[11-13]，最終得到了網絡圖的最小割，也就是原始圖像的最終分割結果。
1.1.2 結合星形先驗的圖割理論
融入形狀先驗能夠獲得更精確更魯棒的分割結果，所以越來越多的研究在分割過程中融入了形狀先驗。Veksler[14]首先提出了星形形狀先驗。星形先驗來自星形凸集這一數學概念，在集合Y中，c是Y的中心點，若對于集合Y內的任意一點p，連接中心點c和p的直線都在集合Y內，則稱Y是星形凸集。文獻[14]給出了星形形狀具體定義與示例：只要能找到一個中心點c滿足條件：對任意一點p，點c和點p連線上的點q仍在目標內，則可認為目標Y是星形形狀。仍用1和0分別作為目標和背景的標記，若p標記為1，則點c和點p連線上的點q標記也為1。星形先驗約束項的能量函數可寫為：
| $ {S_{pq}}\left({{f_p}, {f_q}} \right)=\left\{ \begin{array}{l} \infty {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{if}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {f_{p=1}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{and}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {f_q}=0\ 0{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{oth}}erwise \end{array} \right. $ |
引入星形先驗后，圖割的能量函數變為如下形式：
| $ \begin{array}{l} E\left(f \right)=\sum\limits_{p \in P} {{D_p}\left({{f_p}} \right) + \lambda \sum\limits_{\left({p, q} \right) \in N} {{V_{pq}}\left({{f_p}, {f_q}} \right)} } + \ \sum\limits_{\left({p, q} \right) \in N} {{S_{pq}}\left({{f_p}, {f_q}} \right)} \end{array} $ |
圖像的最優分割對應能量函數式(3)的最小值。
結合星形先驗的圖割理論需要人為交互標記出星形形狀的中心，即圖像的目標種子點。下面將具體說明如何自動獲取目標種子點。
1.2 肺4D-CT各相位目標塊與種子點的獲取
1.2.1 獲取初始相位目標塊與種子點
實驗采用的肺4D-CT數據來源于德克薩斯安德森腫瘤中心DIR實驗室[15]。肺4D-CT圖像包含了肺自由呼吸時不同相位(即不同時刻)的3D-CT圖像，各相位圖像中腫瘤的形態位置和大小會隨著呼吸運動發生變化。因此放療科醫生首先需要根據各相位各層面圖像，對該腫瘤的形態位置和大小作出一個整體的判斷，然后再選取初始相位圖像的目標種子點，即腫瘤中心，還要確定目標塊N的大小，確保以選取的目標種子點為中心形成的N×N×N大小的目標塊，可以在各相位都能完全包含腫瘤。以冠狀面為例，如圖 1所示，目標種子點用紅色標出，以此點為中心建立的初始目標塊用藍框標出。
圖1
初始目標塊與種子點(紅點:種子點；藍框：初始目標塊)
Figure1.
Initial target block and object seed (Red point: selected seed; Blue block: target block)
1.2.2 其余各相位目標塊與種子點的自動獲取
本文采用運動估計中的完全搜索塊匹配算法來實現其余各相位目標塊與種子點的自動獲取。塊匹配算法是目前國際上使用廣泛的，針對序列圖像的運動估計算法[16]。而完全搜索塊匹配算法簡單可靠，并且能夠得到全局最優的結果[17]。完全搜索塊匹配算法有兩種常用的匹配準則：平均平方誤差(mean square error, MSE)和平均絕對誤差(mean absolute difference，MAD)。本文采用平均絕對誤差匹配準則，定義如下：
| $ \begin{array}{l} MAD\left({i, j, k} \right)=\frac{1}{{N \times N \times N}} \times \sum\limits_{n=1}^N {\sum\limits_{n=1}^N {\sum\limits_{n=1}^N {\left| {{f_p}\left({x, y, z} \right)} \right.} } } \ \left. {-{f_{p-1}}\left({x + i, y + i, z + k} \right)} \right| \end{array} $ |
式中，(i, j, k)為位移矢量，fp(x, y, z)表示坐標為(x, y, z)的體素在當前相位的灰度值，fp-1(x+i, y+j, z+k)表示坐標為(x+i, y+j, z+k)的體素在上一相位的灰度值，N×N×N為塊的大小。此方法在整幅圖像的范圍內，逐個體素計算MAD值，從所有點中找出MAD最小值的點，該點對應的塊即為要找的匹配塊，該點處的運動矢量即對應最佳運動矢量。
由此獲得肺4D-CT其余各相位的目標塊以及對應目標塊之間的運動位移。根據獲得的運動位移和初始相位目標種子點的位置，就可計算得出其余各相位的目標種子點。得到的目標種子點作為星形先驗中心，背景種子點默認設置為目標塊的頂點，即實現了結合星形先驗的自動化圖割算法。而且分割腫瘤是在各相位目標塊上，而非整幅三維圖像，此舉縮小了范圍，有利于提高分割精度。
1.3 分割結果量化評價
請經驗豐富的放療科醫生逐幅純手工勾畫腫瘤邊界，由此得到的分割結果作為評價本文實驗結果的金標準。我們分別采用Hausdorff距離(hausdorff distance，HD)和Dice相似系數(dice similarity coefficient，DSC)去量化評價分割結果[18-19]。
HD定義如下：給定兩個有限集合A={a1, a2, …, ap}和B={b1, b2, …, bq}，則A，B之間的HD定義為：
| $ H\left({A, B} \right)=\max \left({h\left({A, B} \right), h\left({B, A} \right)} \right) $ |
其中,
| $ h\left({A, B} \right)=\max {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \min \left\| {{a_i}-{b_j}} \right\| $ |
| $ h\left({B, A} \right)=\max {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \min \left\| {{b_i}-{a_j}} \right\| $ |
式中‖‖表示點集A，B之間的距離范數。函數h(A，B)和h(B，A)分別稱為前向和后向HD。如果h(A，B)=d，則表示A中所有點到B中點的距離不超過d，也就是說A中點都在B中點的距離為d的范圍之內。HD表征了兩個點集之間的不相似程度，因此HD越小則代表分割越精確。
DSC用來度量兩個集合之間的相似性，定義如下：
| $ Dice\left({A, B} \right)=\frac{{2\left({A \cap B} \right)}}{{A + B}} $ |
由上式可知，DSC介于0和1之間，值越大，說明分割結果與金標準越接近。
2 實驗結果
本文實驗采用了兩組肺4D-CT數據，均來自德克薩斯安德森腫瘤中心DIR實驗室[15]，每組數據由10個相位的3D-CT數據組成，包括了極端吸氣和呼氣相位。第一組數據各相位3D圖像大小為320×224×92像素，第二組數據大小為320×224×86像素。第一組數據中，腫瘤位于左肺中部，與周圍組織無黏連；第二組數據中，腫瘤位于右肺下部，與周圍組織多有黏連。基于醫生對腫瘤大小的估計，本實驗中第一組數據目標塊設定為45×45×45像素，第二組數據目標塊設定為75×75×75像素。針對這兩組數據，分別采用傳統圖割算法和本文方法進行分割，并比較兩種算法的分割結果。
2.1 目標塊與目標種子點自動獲取結果
醫生在初始相位的三維圖像上選擇目標種子點，以此點為中心，形成初始目標塊。利用完全搜索塊匹配算法依次得到其他相位的目標塊，同時計算出運動位移。利用運動位移和初始相位目標種子點位置，計算出其余各相位目標種子點。以冠狀面為例，如圖 2所示，分別顯示了兩組數據得到的各相位目標塊以及目標種子點，紅點代表目標種子點，藍框表示初始目標塊。
圖2
獲取的目標塊與目標種子點(紅點：目標種子點，藍框：初始目標塊)
Figure2.
Target blocks and object seeds (Red point: object seeds; Blue block: initial target blocks)
2.2 分割實驗結果
2.2.1 視覺評價
第一組數據相位4的傳統圖割算法和本文方法分割結果的對比圖，如圖 3所示。如本文之前所述，傳統圖割算法分割結果與用戶輸入(選取的種子點集和背景點集)息息相關，對比不難發現，傳統圖割算法若沒有適當選取種子點，會造成分割失敗。在精確選取種子點的情況下，傳統圖割算法雖然能成功分割出腫瘤，但無法去掉圖像中與腫瘤灰度值接近的部分。反觀本文方法，不僅有效克服了傳統圖割算法對于種子點敏感這一缺陷，而且能“干凈”地分割出腫瘤。如圖 4所示，以冠狀面為例，給出了第一組數據相位6的三維圖像在不同層面的典型分割結果，以及其中一層面在不同相位的典型分割結果。
圖3
第一組數據分割結果對比圖
Figure3.
Segmentation result views of the 1st group
圖4
第一組數據典型分割結果
Figure4.
Representative segmentation results of the 1st group
第一組數據中，腫瘤是孤立的，分割難度小，第二組數據中，腫瘤與周圍組織黏連，分割難度大。如圖 5所示為第二組數據相位7傳統圖割算法和本文方法分割結果的對比圖。從中可以看出，針對結構復雜的肺4D-CT腫瘤圖像，在精確選取種子點的情況下，傳統圖割算法也不能完全將腫瘤與黏連部分分割開(紅色箭頭所指處)。而本文方法是在目標塊上分割，并非是在整幅圖像上，所以排除了與腫瘤灰度值接近組織的干擾，同時因為加入了星形先驗，使得腫瘤與黏連部分成功分離。如圖 6所示，以冠狀面為例，給出了第二組數據相位2的三維圖像在不同層面的典型分割結果，以及其中一層面在不同相位的典型分割結果。
圖5
第二組數據分割結果對比圖(紅色箭頭所指即為腫瘤與周圍組織黏連處)
Figure5.
Segmentation result views of the 2nd group (The red arrows point to the adhesion of tumor and surrounding tissue)
圖6
第二組數據典型分割結果
Figure6.
Representative segmentation results of the 2nd group
2.2.2 量化評價
根據本文1.3節所述的HD和DSC兩個指標，分別量化評價兩組實驗結果。HD值越小，DSC值越接近于1，代表分割結果越接近金標準。
在第一組數據的10個相位中，分別計算傳統圖割算法成功分割結果和本文方法分割結果與金標準的HD和DSC。量化評價傳統圖割算法分割結果時，只取腫瘤部分。與傳統圖割算法成功分割結果(HD=22.05±0.99，DSC=0.853±0.008)相比，本文方法分割結果HD值更小(HD=21.79±0.92)，DSC值更接近于1(DSC=0.868±0.008)，說明本文方法分割結果更接近金標準，具體如表 1所示。
表1
第一組數據分割結果量化比較
Table1.
Quantitative evaluation results of 1st group
第二組數據中，腫瘤隨著呼吸運動與周圍組織相連，傳統圖割算法不能成功分割出腫瘤，無法與金標準量化比較，所以如表 2所示，只給出了本文方法與金標準的HD值(24.35±1.04)和DSC值(0.830±0.008)。雖無法與傳統圖割算法比較，卻也說明了在分割復雜肺4D-CT腫瘤圖像時，本文方法仍能獲得較好的結果。
表2
第二組數據本文方法分割結果量化評價
Table2.
Quantitative evaluation results of 2nd group
3 討論
針對肺4D-CT圖像腫瘤分割問題，本文提出了一種基于星形先驗和圖割的自動分割技術。其基本思想是先利用運動估計算法獲得各相位包含腫瘤的目標塊以及目標種子點，再使用結合星形先驗的圖割算法進行分割，達到提升傳統圖割算法自動化程度和分割精度的目的。
本文方法中，首先需要注意的是初始目標塊N的大小。初始目標塊過小，則無法完全包含腫瘤，后續運動估計和分割都失去了意義；初始目標塊過大，雖然可以確保完全包含腫瘤，但是卻增加了整個算法的運算時間，同時也會降低分割精度。因此，本文方法應用于臨床實踐時，需要醫生預估腫瘤大小以設定合適的初始目標塊大小。當前的圖像歸檔與傳輸系統(picture archiving and communication systems，PACS)及其他醫學圖像處理系統一般均包含目標測量功能，因此腫瘤大致尺寸可較易獲得。
另外，不同相位圖像間的運動估計很關鍵，只有確保了運動估計的準確性，才能保證得到合適的目標塊與種子點，提高分割精度。所以本文采用了運動估計中簡單可靠的，并能得到全局最優結果的完全搜索塊匹配算法。
即使是在應對腫瘤與周圍組織黏連這類復雜情況時，本文方法也能夠“干凈”地分割出腫瘤，因為：第一，分割是在目標塊上進行，而非整幅三維圖像，從而去掉了盡可能多的干擾信息，相當于圖像處理中的感興趣區域(region of interest，ROI)技術；第二，背景種子點設置為目標塊頂點和加入的星形先驗，給分割提供了特定的約束條件，從而加強了在目標塊內部腫瘤與黏連部分區分開的能力。
4D-CT能夠真實準確地反映肺部隨呼吸運動的變化規律，對于準確定位靶區、實現精確放療具有重要的應用價值。然而，4D-CT生成的圖像數量巨大，腫瘤靶區不能全靠醫生手工分割，并且還要保證分割的準確性。為此，我們提出了一種基于星形先驗和圖割的自動化分割方法，該方法不僅能大幅提高算法的自動化程度，而且能夠提升分割的準確性。實驗結果表明，本文提出的方法在視覺評價和定量評估方面均優于傳統圖割算法。
引言
肺癌在20世紀末已經成為發病率和死亡率增長最快、對人群健康和生命威脅最大的惡性腫瘤之一[1]。雖然目前對于腫瘤的治療手段越來越多，但手術、放療和化療依然是三種最有效的主要治療手段，放射治療適應癥寬、療效較好，在腫瘤的治療中有著無可置疑的重要地位[2]。因為獨特的空間分辨率和密度分辨率優勢，計算機斷層掃描(computed tomography, CT)在放射治療中扮演著重要的角色。在肺部腫瘤的放射治療中，指導放療計劃和劑量投射的信息通常都來自于肺部CT圖像。近年來，隨著影像技術的進步，將時間因素納入三維計算機斷層掃描(three dimensional computed tomography, 3D-CT)的掃描和重建過程，就出現了四維計算機斷層掃描(four dimensional computed tomography, 4D-CT)。4D-CT圖像與放療技術相結合，對運動靶區的精確放療產生了深遠的影響。利用4D-CT圖像分析靶區及正常器官在呼吸過程中位置和體積的變化規律，讓放射線的劑量最大限度地集中到靶區內，可減少周圍正常組織和器官受到的照射，達到保護正常組織與器官的目的[3]。
放射治療中靶區的精確定位至關重要，靶區的精度會直接影響最終的放療結果，確定靶區可以采用圖像分割技術。目前肺4D-CT腫瘤分割存在兩大主要問題：①4D-CT生成的圖像數量巨大，通常有1 000~2 000張，甚至更多，如果僅靠放療科醫生手工分割每幅圖像中的腫瘤來獲得靶區，顯然費時費力，也不現實[4]；②肺部腫瘤經常與周圍肺壁、膈肌、血管等正常器官組織發生黏連，若采用自動分割技術，誤分割的概率很大，直接影響靶區精度。因此研究出一種自動化程度高、分割準確、魯棒性強的肺腫瘤分割算法十分必要。
自從Boykov和Jolly于2000年提出圖割算法以來[5]，圖割算法憑借其原理簡單、分割精確、全局最優以及能融入先驗知識并且可用于高維分割等優點，一直受到各方重視，成為圖像分割領域中的研究熱點。圖割算法是一種交互式的分割算法，用戶先手動選擇目標種子點和背景種子點，而后交由計算機自動分割出目標物體。相比于其他全自動分割算法，圖割算法因其具有交互性，在分割的準確性上具有明顯優勢，但與此同時也帶來了不可避免的缺陷，分割結果對選取的種子點十分敏感，與用戶輸入息息相關。如果直接將圖割算法用于放療計劃中分割腫瘤，放療科醫生必須精確選取種子點才能確保分割結果良好，這將意味著增加本就十分繁重的醫生工作。
提升圖割算法自動化程度的關鍵在于自動獲取種子點，考慮到腫瘤形態和位置會隨呼吸運動發生改變，并且4D-CT相鄰相位圖像之間存在一定的相似性，本文采用運動估計來自動獲取種子點。同時為了保證分割精度，在傳統圖割算法中加入星形先驗，獲得的種子點就作為星形中心點。具體過程如下：①在4D-CT初始相位的三維圖像上，由醫生選取目標種子點，以此種子點為中心，形成一個初始目標塊，腫瘤及其周邊皆被包含在其中[6]。②采用運動估計中常用的完全搜索塊匹配算法[7]，獲得下一相位圖像中與初始目標塊最相似的目標塊，同時估計出它們之間的運動位移，以此類推，可以得到所有相位的目標塊以及對應塊之間的運動位移。③利用這些運動位移和初始相位目標種子點的位置，計算出其余各相位目標種子點，作為星形先驗的中心點。④在各相位的目標塊上使用結合星形先驗的圖割算法，即可得到腫瘤分割結果。
本文提出的方法只需醫生在初始相位上選取目標種子點，即可獲得同一呼吸周期內所有相位的腫瘤分割結果。不僅大幅提高了圖割算法的自動化程度，降低了放療科醫生的工作量，而且從量化評價指標來看，能夠確保分割結果良好，有效解決了前述目前放療中肺4D-CT腫瘤分割存在的兩大主要問題。
1 理論與方法
1.1 圖割理論
1.1.1 傳統圖割理論
圖割是一種基于圖論的組合優化方法，它將一幅圖像映射成一個網絡圖，其中像素點對應圖中的節點，相鄰像素點之間的關系對應圖中節點之間的邊，在此基礎上建立能量函數，運用最大流/最小割算法對網絡圖進行切割，得到網絡圖的最小割，即原圖的分割結果。圖割理論不僅適用于二維圖像分割，還適用于高維分割，并且基本原理不變。
圖割理論詳細論述參見文獻[8-9]，現概括如下，設G=(V, E)是一個帶有權值的無向圖。其中V為點集，對應圖像中各像素點p的集合P。V還包含兩個額外添加的頂點，稱其為終端，一個為源點，另一個為匯點。E為邊集，對應圖像像素點之間的關系，E包含兩類邊，一類是相鄰像素{p, q}之間的邊，另一類是像素點p與終端之間的邊。給圖中每一條邊賦予權值，設定權值的原則為：像素間的差異越小則權值越大，像素間差異越大則權值越小。用N代表像素點集合P中的4鄰域系統，用1和0分別作為目標和背景的標記，fp∈{1, 0}表示分配給像素點p的標記，由此可定義能量函數：
| $ E\left(f \right)=\sum\limits_{p \in P} {{D_p}\left({{f_p}} \right) + \lambda \sum\limits_{\left({p, q} \right) \in N} {{V_{pq}}\left({{f_p}, {f_q}} \right)} } $ |
能量函數式(1)中右邊第一項為區域項，Dp(fp)是賦予像素點p標記fp的代價測度；右邊第二項為邊界項，Vpq(fp, fq)是賦予相鄰像素點p、q標記fp、fq的代價測度。λ為非負系數，決定了區域項和邊界項所占的比重，本文實驗中λ=300[10]。求E(f)的最小值利用最大流/最小割算法完成[11-13]，最終得到了網絡圖的最小割，也就是原始圖像的最終分割結果。
1.1.2 結合星形先驗的圖割理論
融入形狀先驗能夠獲得更精確更魯棒的分割結果，所以越來越多的研究在分割過程中融入了形狀先驗。Veksler[14]首先提出了星形形狀先驗。星形先驗來自星形凸集這一數學概念，在集合Y中，c是Y的中心點，若對于集合Y內的任意一點p，連接中心點c和p的直線都在集合Y內，則稱Y是星形凸集。文獻[14]給出了星形形狀具體定義與示例：只要能找到一個中心點c滿足條件：對任意一點p，點c和點p連線上的點q仍在目標內，則可認為目標Y是星形形狀。仍用1和0分別作為目標和背景的標記，若p標記為1，則點c和點p連線上的點q標記也為1。星形先驗約束項的能量函數可寫為：
| $ {S_{pq}}\left({{f_p}, {f_q}} \right)=\left\{ \begin{array}{l} \infty {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{if}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {f_{p=1}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{and}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {f_q}=0\ 0{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{oth}}erwise \end{array} \right. $ |
引入星形先驗后，圖割的能量函數變為如下形式：
| $ \begin{array}{l} E\left(f \right)=\sum\limits_{p \in P} {{D_p}\left({{f_p}} \right) + \lambda \sum\limits_{\left({p, q} \right) \in N} {{V_{pq}}\left({{f_p}, {f_q}} \right)} } + \ \sum\limits_{\left({p, q} \right) \in N} {{S_{pq}}\left({{f_p}, {f_q}} \right)} \end{array} $ |
圖像的最優分割對應能量函數式(3)的最小值。
結合星形先驗的圖割理論需要人為交互標記出星形形狀的中心，即圖像的目標種子點。下面將具體說明如何自動獲取目標種子點。
1.2 肺4D-CT各相位目標塊與種子點的獲取
1.2.1 獲取初始相位目標塊與種子點
實驗采用的肺4D-CT數據來源于德克薩斯安德森腫瘤中心DIR實驗室[15]。肺4D-CT圖像包含了肺自由呼吸時不同相位(即不同時刻)的3D-CT圖像，各相位圖像中腫瘤的形態位置和大小會隨著呼吸運動發生變化。因此放療科醫生首先需要根據各相位各層面圖像，對該腫瘤的形態位置和大小作出一個整體的判斷，然后再選取初始相位圖像的目標種子點，即腫瘤中心，還要確定目標塊N的大小，確保以選取的目標種子點為中心形成的N×N×N大小的目標塊，可以在各相位都能完全包含腫瘤。以冠狀面為例，如圖 1所示，目標種子點用紅色標出，以此點為中心建立的初始目標塊用藍框標出。
圖1
初始目標塊與種子點(紅點:種子點；藍框：初始目標塊)
Figure1.
Initial target block and object seed (Red point: selected seed; Blue block: target block)
1.2.2 其余各相位目標塊與種子點的自動獲取
本文采用運動估計中的完全搜索塊匹配算法來實現其余各相位目標塊與種子點的自動獲取。塊匹配算法是目前國際上使用廣泛的，針對序列圖像的運動估計算法[16]。而完全搜索塊匹配算法簡單可靠，并且能夠得到全局最優的結果[17]。完全搜索塊匹配算法有兩種常用的匹配準則：平均平方誤差(mean square error, MSE)和平均絕對誤差(mean absolute difference，MAD)。本文采用平均絕對誤差匹配準則，定義如下：
| $ \begin{array}{l} MAD\left({i, j, k} \right)=\frac{1}{{N \times N \times N}} \times \sum\limits_{n=1}^N {\sum\limits_{n=1}^N {\sum\limits_{n=1}^N {\left| {{f_p}\left({x, y, z} \right)} \right.} } } \ \left. {-{f_{p-1}}\left({x + i, y + i, z + k} \right)} \right| \end{array} $ |
式中，(i, j, k)為位移矢量，fp(x, y, z)表示坐標為(x, y, z)的體素在當前相位的灰度值，fp-1(x+i, y+j, z+k)表示坐標為(x+i, y+j, z+k)的體素在上一相位的灰度值，N×N×N為塊的大小。此方法在整幅圖像的范圍內，逐個體素計算MAD值，從所有點中找出MAD最小值的點，該點對應的塊即為要找的匹配塊，該點處的運動矢量即對應最佳運動矢量。
由此獲得肺4D-CT其余各相位的目標塊以及對應目標塊之間的運動位移。根據獲得的運動位移和初始相位目標種子點的位置，就可計算得出其余各相位的目標種子點。得到的目標種子點作為星形先驗中心，背景種子點默認設置為目標塊的頂點，即實現了結合星形先驗的自動化圖割算法。而且分割腫瘤是在各相位目標塊上，而非整幅三維圖像，此舉縮小了范圍，有利于提高分割精度。
1.3 分割結果量化評價
請經驗豐富的放療科醫生逐幅純手工勾畫腫瘤邊界，由此得到的分割結果作為評價本文實驗結果的金標準。我們分別采用Hausdorff距離(hausdorff distance，HD)和Dice相似系數(dice similarity coefficient，DSC)去量化評價分割結果[18-19]。
HD定義如下：給定兩個有限集合A={a1, a2, …, ap}和B={b1, b2, …, bq}，則A，B之間的HD定義為：
| $ H\left({A, B} \right)=\max \left({h\left({A, B} \right), h\left({B, A} \right)} \right) $ |
其中,
| $ h\left({A, B} \right)=\max {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \min \left\| {{a_i}-{b_j}} \right\| $ |
| $ h\left({B, A} \right)=\max {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \min \left\| {{b_i}-{a_j}} \right\| $ |
式中‖‖表示點集A，B之間的距離范數。函數h(A，B)和h(B，A)分別稱為前向和后向HD。如果h(A，B)=d，則表示A中所有點到B中點的距離不超過d，也就是說A中點都在B中點的距離為d的范圍之內。HD表征了兩個點集之間的不相似程度，因此HD越小則代表分割越精確。
DSC用來度量兩個集合之間的相似性，定義如下：
| $ Dice\left({A, B} \right)=\frac{{2\left({A \cap B} \right)}}{{A + B}} $ |
由上式可知，DSC介于0和1之間，值越大，說明分割結果與金標準越接近。
2 實驗結果
本文實驗采用了兩組肺4D-CT數據，均來自德克薩斯安德森腫瘤中心DIR實驗室[15]，每組數據由10個相位的3D-CT數據組成，包括了極端吸氣和呼氣相位。第一組數據各相位3D圖像大小為320×224×92像素，第二組數據大小為320×224×86像素。第一組數據中，腫瘤位于左肺中部，與周圍組織無黏連；第二組數據中，腫瘤位于右肺下部，與周圍組織多有黏連。基于醫生對腫瘤大小的估計，本實驗中第一組數據目標塊設定為45×45×45像素，第二組數據目標塊設定為75×75×75像素。針對這兩組數據，分別采用傳統圖割算法和本文方法進行分割，并比較兩種算法的分割結果。
2.1 目標塊與目標種子點自動獲取結果
醫生在初始相位的三維圖像上選擇目標種子點，以此點為中心，形成初始目標塊。利用完全搜索塊匹配算法依次得到其他相位的目標塊，同時計算出運動位移。利用運動位移和初始相位目標種子點位置，計算出其余各相位目標種子點。以冠狀面為例，如圖 2所示，分別顯示了兩組數據得到的各相位目標塊以及目標種子點，紅點代表目標種子點，藍框表示初始目標塊。
圖2
獲取的目標塊與目標種子點(紅點：目標種子點，藍框：初始目標塊)
Figure2.
Target blocks and object seeds (Red point: object seeds; Blue block: initial target blocks)
2.2 分割實驗結果
2.2.1 視覺評價
第一組數據相位4的傳統圖割算法和本文方法分割結果的對比圖，如圖 3所示。如本文之前所述，傳統圖割算法分割結果與用戶輸入(選取的種子點集和背景點集)息息相關，對比不難發現，傳統圖割算法若沒有適當選取種子點，會造成分割失敗。在精確選取種子點的情況下，傳統圖割算法雖然能成功分割出腫瘤，但無法去掉圖像中與腫瘤灰度值接近的部分。反觀本文方法，不僅有效克服了傳統圖割算法對于種子點敏感這一缺陷，而且能“干凈”地分割出腫瘤。如圖 4所示，以冠狀面為例，給出了第一組數據相位6的三維圖像在不同層面的典型分割結果，以及其中一層面在不同相位的典型分割結果。
圖3
第一組數據分割結果對比圖
Figure3.
Segmentation result views of the 1st group
圖4
第一組數據典型分割結果
Figure4.
Representative segmentation results of the 1st group
第一組數據中，腫瘤是孤立的，分割難度小，第二組數據中，腫瘤與周圍組織黏連，分割難度大。如圖 5所示為第二組數據相位7傳統圖割算法和本文方法分割結果的對比圖。從中可以看出，針對結構復雜的肺4D-CT腫瘤圖像，在精確選取種子點的情況下，傳統圖割算法也不能完全將腫瘤與黏連部分分割開(紅色箭頭所指處)。而本文方法是在目標塊上分割，并非是在整幅圖像上，所以排除了與腫瘤灰度值接近組織的干擾，同時因為加入了星形先驗，使得腫瘤與黏連部分成功分離。如圖 6所示，以冠狀面為例，給出了第二組數據相位2的三維圖像在不同層面的典型分割結果，以及其中一層面在不同相位的典型分割結果。
圖5
第二組數據分割結果對比圖(紅色箭頭所指即為腫瘤與周圍組織黏連處)
Figure5.
Segmentation result views of the 2nd group (The red arrows point to the adhesion of tumor and surrounding tissue)
圖6
第二組數據典型分割結果
Figure6.
Representative segmentation results of the 2nd group
2.2.2 量化評價
根據本文1.3節所述的HD和DSC兩個指標，分別量化評價兩組實驗結果。HD值越小，DSC值越接近于1，代表分割結果越接近金標準。
在第一組數據的10個相位中，分別計算傳統圖割算法成功分割結果和本文方法分割結果與金標準的HD和DSC。量化評價傳統圖割算法分割結果時，只取腫瘤部分。與傳統圖割算法成功分割結果(HD=22.05±0.99，DSC=0.853±0.008)相比，本文方法分割結果HD值更小(HD=21.79±0.92)，DSC值更接近于1(DSC=0.868±0.008)，說明本文方法分割結果更接近金標準，具體如表 1所示。
表1
第一組數據分割結果量化比較
Table1.
Quantitative evaluation results of 1st group
第二組數據中，腫瘤隨著呼吸運動與周圍組織相連，傳統圖割算法不能成功分割出腫瘤，無法與金標準量化比較，所以如表 2所示，只給出了本文方法與金標準的HD值(24.35±1.04)和DSC值(0.830±0.008)。雖無法與傳統圖割算法比較，卻也說明了在分割復雜肺4D-CT腫瘤圖像時，本文方法仍能獲得較好的結果。
表2
第二組數據本文方法分割結果量化評價
Table2.
Quantitative evaluation results of 2nd group
3 討論
針對肺4D-CT圖像腫瘤分割問題，本文提出了一種基于星形先驗和圖割的自動分割技術。其基本思想是先利用運動估計算法獲得各相位包含腫瘤的目標塊以及目標種子點，再使用結合星形先驗的圖割算法進行分割，達到提升傳統圖割算法自動化程度和分割精度的目的。
本文方法中，首先需要注意的是初始目標塊N的大小。初始目標塊過小，則無法完全包含腫瘤，后續運動估計和分割都失去了意義；初始目標塊過大，雖然可以確保完全包含腫瘤，但是卻增加了整個算法的運算時間，同時也會降低分割精度。因此，本文方法應用于臨床實踐時，需要醫生預估腫瘤大小以設定合適的初始目標塊大小。當前的圖像歸檔與傳輸系統(picture archiving and communication systems，PACS)及其他醫學圖像處理系統一般均包含目標測量功能，因此腫瘤大致尺寸可較易獲得。
另外，不同相位圖像間的運動估計很關鍵，只有確保了運動估計的準確性，才能保證得到合適的目標塊與種子點，提高分割精度。所以本文采用了運動估計中簡單可靠的，并能得到全局最優結果的完全搜索塊匹配算法。
即使是在應對腫瘤與周圍組織黏連這類復雜情況時，本文方法也能夠“干凈”地分割出腫瘤，因為：第一，分割是在目標塊上進行，而非整幅三維圖像，從而去掉了盡可能多的干擾信息，相當于圖像處理中的感興趣區域(region of interest，ROI)技術；第二，背景種子點設置為目標塊頂點和加入的星形先驗，給分割提供了特定的約束條件，從而加強了在目標塊內部腫瘤與黏連部分區分開的能力。
4D-CT能夠真實準確地反映肺部隨呼吸運動的變化規律，對于準確定位靶區、實現精確放療具有重要的應用價值。然而，4D-CT生成的圖像數量巨大，腫瘤靶區不能全靠醫生手工分割，并且還要保證分割的準確性。為此，我們提出了一種基于星形先驗和圖割的自動化分割方法，該方法不僅能大幅提高算法的自動化程度，而且能夠提升分割的準確性。實驗結果表明，本文提出的方法在視覺評價和定量評估方面均優于傳統圖割算法。
