楊澤富，戰(zhàn)蔭偉，楊榮騫

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

0 引 言

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

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

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

1 算法設(shè)計(jì)基礎(chǔ)

1.1 體數(shù)據(jù)空間結(jié)構(gòu)表示

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

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

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

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

其中，Nx和Ny表示單張切片的長(zhǎng)和寬，Nz表示體數(shù)據(jù)中切片的張數(shù)。

圖1 CT圖像序列和體數(shù)據(jù)

1.2 FM方法介紹

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

(1)

其中

因此將式(1)展開可得

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

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

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

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

(2)從已知點(diǎn)出發(fā)向外訪問其四鄰域點(diǎn)，如果鄰域點(diǎn)是未知點(diǎn)則將其狀態(tài)改為實(shí)驗(yàn)點(diǎn)；

(3)根據(jù)式(1)計(jì)算并更新所有實(shí)驗(yàn)點(diǎn)的值；

(4)在所有實(shí)驗(yàn)點(diǎn)中選擇值最小的一點(diǎn)，將其狀態(tài)改為已知點(diǎn)；

(5)返回步驟(2)繼續(xù)循環(huán)操作。

2 改進(jìn)的FM方法

圖3為改進(jìn)的FM方法流程圖，本節(jié)將根據(jù)該流程圖具體描述改進(jìn)的FM方法。

圖3 算法流程

2.1 自動(dòng)設(shè)置初始種子點(diǎn)集

FM方法需要手動(dòng)設(shè)置種子點(diǎn)，屬于半自動(dòng)分割模式，本文設(shè)計(jì)一種自動(dòng)設(shè)置初始種子點(diǎn)集的方式，將FM方法改進(jìn)為全自動(dòng)分割模式。因?yàn)镃T圖像中腦血腫區(qū)域與相鄰正常組織之間存在差異，所以可以使用閾值處理提取種子點(diǎn)，但考慮到頭皮組織對(duì)閾值處理的影響，因此在閾值處理前應(yīng)該先提取腦實(shí)質(zhì)，圖4描述了具體的流程。

圖4 自動(dòng)設(shè)置初始種子點(diǎn)集流程

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

(2)

隨后根據(jù)提取的顱骨計(jì)算腦實(shí)質(zhì)區(qū)域的二值模板，即統(tǒng)計(jì)顱骨內(nèi)邊界中包含的所有體素，并將它們作為二值模板的前景，剩下的作為模板的背景。

體數(shù)據(jù)乘法濾波器對(duì)輸入的兩份大小相等的體數(shù)據(jù)，進(jìn)行體素級(jí)別的乘法操作，即將兩份體數(shù)據(jù)中坐標(biāo)相同的體素對(duì)應(yīng)相乘，并將結(jié)果作為一份新的體數(shù)據(jù)輸出。由于上述二值模板的前景對(duì)應(yīng)腦實(shí)質(zhì)區(qū)域且值為1，背景值為0，因此體數(shù)據(jù)乘法濾波器的輸出結(jié)果為腦實(shí)質(zhì)區(qū)域。下一步對(duì)腦實(shí)質(zhì)區(qū)域進(jìn)行閾值處理，將[40，100]設(shè)為閾值區(qū)間[12]，則有閾值函數(shù)

(3)

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

2.2 確定感興趣區(qū)域包圍盒

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

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

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

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

圖5 二維包圍盒計(jì)算過程

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

(4)

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

(5)

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

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

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

體數(shù)據(jù)下采樣前后同一相對(duì)位置坐標(biāo)之間存在比例關(guān)系。假設(shè)空間坐標(biāo)系x軸和y軸方向采樣因子等于n，z軸方向采樣因子等于1，則包圍盒R2的起點(diǎn)坐標(biāo) (x0,y0,z0) 在保持相對(duì)位置不變的情況下，經(jīng)下采樣后將變?yōu)?x0/n,y0/n,z0)，長(zhǎng)寬高也分別變?yōu)閍/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}

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

圖6 根據(jù)采樣因子逆推包圍盒

圖7 感興趣區(qū)域金字塔迭代分割

2.3 ROI金字塔迭代分割

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

ROI金字塔迭代分割算法：

輸入：腦血腫CT圖像序列，即原始體數(shù)據(jù)；初始種子點(diǎn)集。

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

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

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

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

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

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

分析FM方法可以發(fā)現(xiàn)，其總是從種子點(diǎn)出發(fā)向外尋找邊界，且如果出發(fā)位置距離邊界越近，演化計(jì)算的次數(shù)越少，反之則越多。因此如果通過某種方式使種子點(diǎn)盡量貼近邊界內(nèi)側(cè)，就可大量減少計(jì)算量，縮短FM方法的分割時(shí)間。本文根據(jù)該點(diǎn)提出上述ROI金字塔迭代分割算法，即除初始種子點(diǎn)集之外，ROI金字塔每一層的種子點(diǎn)集，都是通過對(duì)上一層分割結(jié)果先上采樣后腐蝕獲得的，這樣處理得到的種子點(diǎn)都將包含于分割目標(biāo)內(nèi)部，且種子點(diǎn)集合的最外側(cè)點(diǎn)將貼近于目標(biāo)邊界，又因?yàn)榉N子點(diǎn)集合的內(nèi)部區(qū)域均為已知點(diǎn)，F(xiàn)M方法會(huì)忽略內(nèi)部區(qū)域直接從最外側(cè)種子點(diǎn)開始計(jì)算，所以只需要少量的計(jì)算即可完成分割?？偠灾?，在ROI金字塔迭代分割中，通過使用上一層的分割結(jié)果來(lái)估計(jì)下一層的邊界，即上一層的分割可視為粗分割，當(dāng)前層在粗分割的基礎(chǔ)上使用FM方法再做一次細(xì)分割，如此層層迭代分割下去，從而避免大量多余的計(jì)算，達(dá)到提高分割速度的目的。

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

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

3.1 評(píng)估指標(biāo)

本文對(duì)分割結(jié)果的評(píng)估指標(biāo)主要包括：Dice相似系數(shù)(Dice coefficient)，靈敏度(Sensitivity)，特異性(Specifi-city)和準(zhǔn)確率(Accuracy)。Dice相似系數(shù)用來(lái)衡量真實(shí)結(jié)果與預(yù)測(cè)結(jié)果之間的相似性，預(yù)測(cè)結(jié)果越接近真實(shí)結(jié)果，Dice值越大。若用U1表示Dice的大小，則

(6)

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

(7)

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

(8)

準(zhǔn)確率表示分割正確的像素個(gè)數(shù)占總像素個(gè)數(shù)的比例，若用U4表示準(zhǔn)確率，則

(9)

3.2 實(shí)驗(yàn)對(duì)比

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

圖8 腦血腫圖像分割結(jié)果對(duì)比

圖8(a)為原始CT圖像；圖8(b)為專家手工分割結(jié)果，并將其作為實(shí)驗(yàn)分割的金標(biāo)準(zhǔn)；圖8(c)為文獻(xiàn)[10]的分割結(jié)果，文獻(xiàn)[10]是目前關(guān)于FM方法最新的改進(jìn)方法；圖8(d)為本文方法的分割結(jié)果。通過對(duì)比可以發(fā)現(xiàn)，本文方法的結(jié)果與專家手工分割結(jié)果較為接近，分割較為精確。

為了更加客觀地反映本文分割算法的有效性，將本文方法分割結(jié)果的各項(xiàng)評(píng)估指標(biāo)與原FM方法，文獻(xiàn)[9]和文獻(xiàn)[10]的各項(xiàng)評(píng)估指標(biāo)進(jìn)行對(duì)比，其中文獻(xiàn)[9]為最新的水平集分割腦血腫方法，對(duì)比結(jié)果見表1。

表1 本文方法與其它方法的評(píng)估指標(biāo)對(duì)比

由表1可知，在本文方法分割結(jié)果的各項(xiàng)評(píng)估指標(biāo)中，除TNR一項(xiàng)之外，其它評(píng)估指標(biāo)均為最優(yōu)，說明本文方法是有效的。

在臨床應(yīng)用中，分割的時(shí)間效率也是需要考慮的重要因素。因?yàn)槊刻證T數(shù)據(jù)集的張數(shù)不同，所以本文選用單張CT圖像的平均處理時(shí)間來(lái)評(píng)估分割的時(shí)間效率。為了進(jìn)一步說明本文方法在時(shí)間效率上的優(yōu)勢(shì)，將其與其它方法進(jìn)行對(duì)比，對(duì)比結(jié)果見表2。從表2可以看出，本文方法分割單張CT圖像平均用時(shí)0.043 s，分割效率明顯優(yōu)于其它3種方法，特別是在處理的圖像數(shù)量較大時(shí)，分割效率的優(yōu)勢(shì)將更加明顯。

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

4 結(jié)束語(yǔ)

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