腦血腫CT圖像分割的改進快速行進方法

楊澤富，戰蔭偉，楊榮騫

(1.廣東工業大學 計算機學院，廣東 廣州 510006；2.華南理工大學 材料科學與工程學院，廣東 廣州 510640)

0 引 言

腦血腫CT圖像的分割是規劃腦出血手術的基礎，由于手動分割及半自動分割均非常耗時，所以應采用效率更高的全自動分割方法。目前醫學圖像分割領域已有大量的研究工作，其中多數算法用于分割二維切片。如文獻[1]利用極限樹設計分類器識別病灶，解決了手動分割重復性差的問題且穩定性較強。文獻[2]使用隨機森林識別血腫區域，但對小面積區域識別較差。文獻[3]利用離散曲率波變換對圖像進行多尺度幾何分析，進而區分正常腦組織和病變組織。文獻[4]雖然改進并減少水平集的演化時間，但仍需初始化零水平集，分割效率較低。

臨床上對病灶進行三維重建可以輔助放射科醫師更加直觀地了解病癥，減小誤診率。因為體數據是三維重建的基礎，所以針對體數據分割技術的研究逐漸成為熱點。如文獻[5]采用有監督學習方式分析病灶的鄰域數據從而實現自動分割，分割精度較高。文獻[6]結合社會群體優化算法和模糊熵理論提出一種半自動圖像處理方法，在ISLES 2015競賽中有較好的表現。文獻[7,8]使用卷積神經網絡識別病灶，分割結果較為準確，但算法對數據集的依賴性較強，穩定性較差。以上算法雖然都能取得較好的分割結果，但分割速度較慢，為此文獻[9]提出一種特殊的三維水平集算法，通過增加多個零水平集來加快分割速度。文獻[10]在FM方法的基礎上提出一種Eikonal方程快速解法，以此來加快FM方法的分割速度。

文獻[9,10]雖然分割速度較快，但距離臨床時間要求仍有一定的距離，由于手術計劃的制定過程是串行的，下一步的執行必須使用上一步的結果，因此如果分割速度較慢，將增加整個手術計劃的制定時間，降低醫生的工作效率，為此本文提出一種分割精度較高且速度更快的分割方法，輔助醫生更加快捷地制定手術計劃。

1 算法設計基礎

1.1 體數據空間結構表示

圖1左側展示了一組連續的CT圖像切片序列，右側為該序列按掃描順序排列后得到的掃描部位體數據。體數據由體素組成，體素是基本體積元素，可以理解為三維空間內排列的點或者一小塊區域。設Nx、Ny和Nz均為大于1的正整數，體數據中每一張CT圖像切片可視作一個大小為Nx×Ny的二維網格。若將切片記為Ik，k表示切片的序號，則可定義

Ik={V(i,j,k);1≤i≤Nx,1≤j≤Ny}

其中，V(i,j,k) 表示空間坐標為 (i,j,k) 的體素，(i,j,k)∈Z3,Z表示整數集。同理，體數據可視為一個規模大小等于Nx×Ny×Nz的三維空間區域，若用V表示體數據，參照二維切片的集合定義，則有

V={V(i,j,k);1≤i≤Nx,1≤j≤Ny,1≤k≤Nz}

其中，Nx和Ny表示單張切片的長和寬，Nz表示體數據中切片的張數。

圖1 CT圖像序列和體數據

1.2 FM方法介紹

FM方法的基本思想是根據輸入的速度圖像進行輪廓曲線演化，并在時間圖像中提取演化后的輪廓曲線，其輸出結果為一幅二值掩模圖像。其中速度圖像上每點表示輪廓曲線在該點的演化速度，時間圖像上每點表示輪廓曲線到達該點所需的時間。如果用F(x,y) 表示速度圖像上點 (x,y) 的值，用T(x,y) 表示時間圖像上點 (x,y) 的值，則可滿足Eikonal方程

(1)

其中

因此將式(1)展開可得

結合上述方程式，圖2描述了FM方法中輪廓曲線的演化過程。

圖2 FM方法中的輪廓曲線演化

在FM方法中輪廓曲線從初始種子點開始進行演化，圖2中使用黑色網格點表示已知點，灰色網格點表示實驗點，白色網格點表示未知點。在輪廓曲線演化過程中，首先設置一點為種子點，并將該種子點作為已知點，同時將其余點設為未知點，然后通過下述步驟進行曲線演化：

(1)設置已知點的初始值，一般將初始值設為0，表示曲線演化到該點所需的時間等于0；

(2)從已知點出發向外訪問其四鄰域點，如果鄰域點是未知點則將其狀態改為實驗點；

(3)根據式(1)計算并更新所有實驗點的值；

(4)在所有實驗點中選擇值最小的一點，將其狀態改為已知點；

(5)返回步驟(2)繼續循環操作。

2 改進的FM方法

圖3為改進的FM方法流程圖，本節將根據該流程圖具體描述改進的FM方法。

圖3 算法流程

2.1 自動設置初始種子點集

FM方法需要手動設置種子點，屬于半自動分割模式，本文設計一種自動設置初始種子點集的方式，將FM方法改進為全自動分割模式。因為CT圖像中腦血腫區域與相鄰正常組織之間存在差異，所以可以使用閾值處理提取種子點，但考慮到頭皮組織對閾值處理的影響，因此在閾值處理前應該先提取腦實質，圖4描述了具體的流程。

圖4 自動設置初始種子點集流程

在自動設置種子點的過程中，體數據的每張切片大小均為512×512，設采樣因子為4對體數據進行下采樣操作，即在切片的橫向和縱向上每隔4個體素采樣一次，得到大小為128×128的較小切片。然后對下采樣結果進行閾值處理提取顱骨，由于顱骨的HU最低約為150[11]，與其它組織差異較為明顯，因此可將閾值函數設置為

(2)

隨后根據提取的顱骨計算腦實質區域的二值模板，即統計顱骨內邊界中包含的所有體素，并將它們作為二值模板的前景，剩下的作為模板的背景。

體數據乘法濾波器對輸入的兩份大小相等的體數據，進行體素級別的乘法操作，即將兩份體數據中坐標相同的體素對應相乘，并將結果作為一份新的體數據輸出。由于上述二值模板的前景對應腦實質區域且值為1，背景值為0，因此體數據乘法濾波器的輸出結果為腦實質區域。下一步對腦實質區域進行閾值處理，將[40，100]設為閾值區間[12]，則有閾值函數

(3)

然后對閾值處理的結果進行數學形態學二值腐蝕操作，消除離散的噪聲點同時收縮種子點區域，得到待篩選的種子點集合，即圖4(e)白色區域。在確認種子點時，逐一判斷每張切片上的種子點個數，如果數量大于設定值20，說明該切片上包含血腫塊，保留種子點；反之則舍去該切片上的所有種子點。

2.2 確定感興趣區域包圍盒

分割過程中應該先提取體數據內的腦血腫區域作為感興趣區域(region of interest，ROI)來減少計算量。ROI一般通過包圍盒的形式表示。為了便于理解三維包圍盒的計算過程，先介紹二維包圍盒的計算過程。

二維包圍盒的計算即在圖像上對ROI加框。若用R1表示二維包圍盒，則可定義

R1={p(x,y);x0≤x≤x0+l,y0≤y≤y0+w}

其中，p(x,y) 表示圖像坐標等于 (x,y) 的像素，(x0,y0) 為包圍盒矩形的左上角起點坐標，l和w分別為包圍盒的長和寬，如圖5所示。

圖5 二維包圍盒計算過程

二維包圍盒的計算過程為：首先對原始圖像進行中值濾波去除脈沖噪聲，然后計算圖像的梯度得到梯度圖像，再通過Sigmoid函數對梯度圖像進行非線性映射求得速度圖像。其中Sigmoid函數公式為

(4)

X表示輸入像素大小，Y表示經變換后輸出的像素大小，P和Q分別表示輸出像素的最大值和最小值，α和β為大于零的常數。若令P=1，Q=0，式(4)可化簡為

(5)

式中：當X→-∞時，Y→1，當X→+∞時，Y→0，即對于定義域X∈(-∞,+∞)，值域Y∈(0,1)。 此處采用式(5)對梯度圖像做歸一化處理，經多次實驗驗證，當α=0.3，β=2.0時處理結果較為穩定。隨后根據速度圖像使用FM方法分割，遍歷輸出的二值圖像，從而確定包圍盒起點坐標、長度和寬度，求得包圍盒。

在二維平面上，包圍盒是矩形，而在三維空間中，包圍盒是長方體。若用R2表示三維包圍盒，則

R2={V(x,y,z);x0≤x≤x0+a,y0≤y≤y0+b;z0≤z≤z0+c}

體數據下采樣前后同一相對位置坐標之間存在比例關系。假設空間坐標系x軸和y軸方向采樣因子等于n，z軸方向采樣因子等于1，則包圍盒R2的起點坐標 (x0,y0,z0) 在保持相對位置不變的情況下，經下采樣后將變為(x0/n,y0/n,z0)，長寬高也分別變為a/n，b/n和c，從而得到下采樣后的包圍盒R′2，即

R′2={V(x,y,z);(x0/n)≤x≤(x0/n)+(a/n),(y0/n)≤y≤(y0/n)+(b/n),z0≤z≤z0+c}

因此可以先求體數據下采樣之后的包圍盒，再根據采樣因子逆推下采樣之前的包圍盒，從而節省計算時間。逆推包圍盒的求解過程如圖6所示。

圖6 根據采樣因子逆推包圍盒

圖7 感興趣區域金字塔迭代分割

2.3 ROI金字塔迭代分割

上文介紹了快速求取ROI的方法，在該方法的基礎上參考圖像金字塔原理可以建立一種特殊的金字塔，ROI金字塔。ROI金字塔是一系列以金字塔形狀排列的體數據集合，金字塔的每一層是三維的體而不是二維的面。ROI金字塔的底部為原始體數據ROI，往上每一層都通過對相鄰層進行下采樣獲得，層數越高，體數據尺寸越小，圖7中虛線框內為具體示例。

ROI金字塔迭代分割算法：

輸入：腦血腫CT圖像序列，即原始體數據；初始種子點集。

輸出：分割后的二值圖像序列。

(1)根據前文介紹的方法提取體數據ROI，并通過下采樣建立ROI金字塔；

(2)從金字塔頂層的體數據開始，輸入前文求得的初始種子點集，使用三維FM方法分割；

(3)采用最近鄰插值法對分割結果上采樣，上采樣結果的大小和金字塔下一層相等，然后使用半徑為一的探針做二值腐蝕操作，提取腐蝕結果中的白色體素作為金字塔下一層的種子點集；

(4)輸入新的種子點集，使用三維FM方法對金字塔的下一層進行分割；

(5)返回步驟(3)，迭代分割直至金字塔底層為止。

分析FM方法可以發現，其總是從種子點出發向外尋找邊界，且如果出發位置距離邊界越近，演化計算的次數越少，反之則越多。因此如果通過某種方式使種子點盡量貼近邊界內側，就可大量減少計算量，縮短FM方法的分割時間。本文根據該點提出上述ROI金字塔迭代分割算法，即除初始種子點集之外，ROI金字塔每一層的種子點集，都是通過對上一層分割結果先上采樣后腐蝕獲得的，這樣處理得到的種子點都將包含于分割目標內部，且種子點集合的最外側點將貼近于目標邊界，又因為種子點集合的內部區域均為已知點，FM方法會忽略內部區域直接從最外側種子點開始計算，所以只需要少量的計算即可完成分割?？偠灾赗OI金字塔迭代分割中，通過使用上一層的分割結果來估計下一層的邊界，即上一層的分割可視為粗分割，當前層在粗分割的基礎上使用FM方法再做一次細分割，如此層層迭代分割下去，從而避免大量多余的計算，達到提高分割速度的目的。

3 實驗結果與分析

實驗系統配置為Intel(R) Core(TM) i5-8400 CPU @ 2.80 GHz，16 GB RAM，操作系統為Windows 10專業版，軟件環境為Microsoft Visual Studio 2015。實驗數據來源于中山大學附屬第一醫院神經外科的8套臨床CT數據集，共計2536張，數據集層厚均為0.5 mm，單張CT圖像大小為512×512。實驗中使用專家的手工分割結果作為分割的金標準，同時將本文方法與原FM方法，以及文獻[9,10]進行對比，以此驗證本文方法的有效性。

3.1 評估指標

本文對分割結果的評估指標主要包括：Dice相似系數(Dice coefficient)，靈敏度(Sensitivity)，特異性(Specifi-city)和準確率(Accuracy)。Dice相似系數用來衡量真實結果與預測結果之間的相似性，預測結果越接近真實結果，Dice值越大。若用U1表示Dice的大小，則

(6)

靈敏度又稱召回率(Recall)，或稱真陽性率(true positive rate，TPR)，用來衡量算法對正樣本的識別能力，若用U2表示TPR的值，則

(7)

特異性又稱為真陰性率(true negative rate，TNR)，用來衡量算法對負樣本的識別能力，若用U3表示TNR的值，則

(8)

準確率表示分割正確的像素個數占總像素個數的比例，若用U4表示準確率，則

(9)

3.2 實驗對比

為了直觀地觀察算法分割效果，圖8給出了包括本文方法在內的3種不同方法，分別提取腦血腫輪廓的結果對比。

圖8 腦血腫圖像分割結果對比

圖8(a)為原始CT圖像；圖8(b)為專家手工分割結果，并將其作為實驗分割的金標準；圖8(c)為文獻[10]的分割結果，文獻[10]是目前關于FM方法最新的改進方法；圖8(d)為本文方法的分割結果。通過對比可以發現，本文方法的結果與專家手工分割結果較為接近，分割較為精確。

為了更加客觀地反映本文分割算法的有效性，將本文方法分割結果的各項評估指標與原FM方法，文獻[9]和文獻[10]的各項評估指標進行對比，其中文獻[9]為最新的水平集分割腦血腫方法，對比結果見表1。

表1 本文方法與其它方法的評估指標對比

由表1可知，在本文方法分割結果的各項評估指標中，除TNR一項之外，其它評估指標均為最優，說明本文方法是有效的。

在臨床應用中，分割的時間效率也是需要考慮的重要因素。因為每套CT數據集的張數不同，所以本文選用單張CT圖像的平均處理時間來評估分割的時間效率。為了進一步說明本文方法在時間效率上的優勢，將其與其它方法進行對比，對比結果見表2。從表2可以看出，本文方法分割單張CT圖像平均用時0.043 s，分割效率明顯優于其它3種方法，特別是在處理的圖像數量較大時，分割效率的優勢將更加明顯。

表2 單張CT圖像的平均分割時間/s

4 結束語

本文對FM方法進行改進，提出一種三維全自動腦血腫CT圖像序列分割方法，實現了全自動、迅速、準確的腦血腫區域輪廓提取。通過簡單的閾值處理提取種子點，解決了FM方法需要手動設置種子點的問題，使算法得以全自動執行。建立一種特殊的感興趣區域金字塔進行迭代分割，并采用先做粗分割，再做細分割的設計思想，從而在保證分割精度的同時，極大地減少了算法的計算量。實驗結果表明，本文方法對比其它分割方法，不僅在分割精度上占優，而且在分割效率上也有著明顯的優勢，說明其具有高效性，實用性和有效性，能夠為臨床手術計劃的制定提供有利條件。