針對肺部病變及支氣管干擾等因素導致的肺實質分割困難的問題，本文提出一種融合表面波（surfacelet）變換與脈沖耦合神經網絡（PCNN）的肺實質分割算法。首先，通過 surfacelet 變換對三維肺部計算機斷層掃描數據進行多尺度多方向分解，利用局部修正拉普拉斯算子選擇處理后的子帶系數增強圖像的邊緣特征；然后，經 surfacelet 逆變換得到增強后的圖像作為 PCNN 的反饋輸入；最后，通過循環迭代完成肺實質的分割。所提算法對公開數據集中的樣本進行了測試。結果表明，本文算法的分割性能優于 surfacelet 變換邊緣提取算法、三維區域生長算法和三維 U 形網絡（U-NET）算法，能夠有效抑制肺部病變及支氣管的干擾，得到更完整的肺實質圖像。
引用本文： 張華海, 白培瑞, 郭子楊, 杜令豪, 李昶, 任延德, 楊凱, 劉慶一. 一種融合表面波變換與脈沖耦合神經網絡的三維肺實質分割算法. 生物醫學工程學雜志, 2020, 37(4): 630-640. doi: 10.7507/1001-5515.201908060 復制
引言
從計算機斷層掃描圖像（computed tomography，CT）中準確地分割肺實質是計算機輔助診斷（computer aided diagnosis，CAD）技術的關鍵步驟，對于早期診斷肺癌等肺部疾病具有重要意義[1-2]。另外，肺的邊沿分割精度直接影響肺部體積測量的結果，而肺部體積的測量對研究氣胸、肺氣腫等疾病具有重要意義。肺實質分割方法主要分成三類：第一類是基于傳統方法的肺實質分割，如閾值法[3]、區域生長[4]、邊緣檢測[5]、活動輪廓[6]等，各有一定缺陷。基于閾值的方法對噪聲和成像偽影比較敏感[7]，基于區域生長的方法易出現欠分割或過分割現象，基于邊緣檢測的方法易產生偽邊界或邊界丟失現象[8]。第二類是不同方法相結合的肺實質分割。比如，Nizami 等[9]提出使用小波框架與 K-均值（K-means）聚類相結合分割肺實質的算法；Birkbeck 等[10]提出將統計學習與來自鄰近解剖結構（心臟、肝臟、脾臟和肋骨）的解剖學約束相結合的算法；Zhou 等[11]提出利用增強的閾值算法自動分割肺實質，然后通過紋理感知的活動輪廓模型平滑血管結構和肺部邊界的算法。第三類是基于深度學習的方法。Kumar 等[12]提出基于 U 形網絡（U-NET）的深度卷積神經網絡，可以實現對肺實質的自動分割，該模型采用數據增加（data augmentation）技術，使具有跳接結構的深度網絡能夠根據少量可用的醫學圖像進行像素分類。Kamal 等[13]提出一種融合卷積和遞歸神經網絡的三維混合 U-NET（three dimensional dense UNET，3D-denseUNET）框架，用于分割 CT 圖像中的肺部腫瘤。Harrison 等[14]引入一種基于全卷積神經網絡（fully convolutional neural network，FCN）的漸進式整體嵌套網絡（progressive holistically-nested networks，P-HNNs），用于分割肺組織病變。此外，近幾年出現了將傳統特征提取方法與深度學習方法融合的思路。比如，Ngo 等[15]提出將深度置信網絡與距離正則化水平集相結合，利用深度置信網絡的學習結果驅動零水平集曲線演化進行肺實質的分割。但是，深度學習分割方法的性能依賴于可用的醫學圖像數量，對數據樣本質量和計算能力的要求較高[16]。
本文提出一種融合表面波（surfacelet）變換與脈沖耦合神經網絡（pulse coupled neural network，PCNN）的三維肺實質分割方法，利用 surfacelet 變換和 PCNN 優勢互補以提高肺實質分割的準確性。Surfacelet 變換是一種多分辨率變換，能夠有效提取相鄰幀間的多尺度多方向結構信息，可以通過對變換域各子帶信息進行濾波以增強信號細節[17]。脈沖耦合神經網絡是一種特殊的神經網絡，能夠考慮像素之間灰度的幾何分布選擇合適的歸類，有效分離對象和背景分布之間的重疊，準確提取感興趣區域。所以，利用 surfacelet 變換增強 CT 圖像的邊緣與細節信息，可以彌補 PCNN 只考慮圖像灰度信息的不足，從而能夠減少 CT 圖像中肺部病變及支氣管的干擾，實現肺實質的精確分割。
1 方法
1.1 基于 surfacelet 變換的圖像增強
Surfacelet 變換能夠利用不同尺度、不同頻率的方向子帶準確地捕獲 N 維信號中的面狀奇異，具有平移不變性、完全重構、低冗余、多尺度多方向等性質[18]。Surfacelet 變換的主要思想是將 N 維方向濾波器組（N-dimensional directional filter bank，NDFB）與多尺度金字塔結構結合，利用方向濾波器組實現金字塔形狀的 N 維頻域方向子帶劃分，利用多尺度金字塔結構實現多尺度分析與重構。本文主要考慮三維圖像分割，所以下面介紹 surfacelet 變換的三維形式。
1.1.1 三維方向濾波器組
三維方向濾波器組（three dimensional directional filter bank，3D-DFB）能夠將三維信號的頻譜分解為不同的方向子帶。假設 n1、n2 和 n3 分別表示三維空間域沿垂直方向、水平方向和深度方向相互正交的坐標軸，ω1、ω2 和 ω3 是三維頻域對應的坐標軸，則 3D-DFB 沿 ω1 軸的理想頻率分解示意圖如圖 1 所示[19]，分別表示沿（n1，n2）平面的二維理想楔形頻率子帶劃分、沿（n1，n3）平面的二維理想楔形頻率子帶劃分和金字塔形方向子帶，將兩個二維楔形頻率子帶劃分進行重疊交叉，就會形成頂點對稱的包含多個金字塔形方向子帶（或通道）的沙漏濾波器（hourglass-shaped filter），深灰色部分表示其中的一個金字塔形子帶（pyramid-shaped subband）。類似地，沿 ω2 軸和 ω3 軸的沙漏濾波器可以通過旋轉沙漏濾波器得到。
圖1
沿 ω1 軸的金字塔形頻率子帶分解示意圖
Figure1.
Schematic diagrams of pyramid-shaped frequency sub-band decomposition along ω1 axis
假設 a 表示如圖 1 所示沿（n1，n2）平面的某一層二維濾波器的序號，l2 表示沿（n1，n2）方向的分解層數，b 表示如圖 1 所示沿（n1，n3）平面的某一層二維濾波器的序號，l3 表示沿（n1，n3）方向的分解層數，則如圖 1 所示沿（n1，n2）平面的楔形濾波器可以表示為 
（ω1，ω2），0 ≤ a<
，沿（n1，n3）平面的楔形濾波器可以表示為 
（ω1，ω3），0 ≤ b<
。兩個方向的分解層數決定沙漏濾波器包含的子帶數（或通道數），沙漏濾波器可以根據分解層數和層序號表示為 
。比如，l2 = 2，l3 = 2，則如圖 1 所示的金字塔形方向子帶共有 
=24 = 16 個，其中用深灰色顯示的子帶下標號為 a = 1，b = 1，所以該子帶可以用 
表示。如圖 1 所示的重疊交叉操作可以表示為兩個楔形濾波器做乘積運算，如式（1）所示：
 |
在具體設計包含多個金字塔形狀方向頻域子帶的沙漏濾波器時，需要兩個步驟：
第一步，設計沙漏濾波器，首先采用平滑凱塞爾窗對二維辛格（sinc）函數進行截取，做二維傅里葉變換分別得到扇形結構的二維濾波器 E（ω1，ω2）和 E（ω1，ω3），然后將二者相乘得到沿 ω1 軸近似理想形狀的沙漏濾波器 K（ω1，ω2，ω3），如式（2）所示：
 |
另外，沿 ω2 軸和 ω3 軸的沙漏濾波器 K（ω2，ω3，ω1）和 K（ω3，ω1，ω2），可以由該濾波器旋轉得到。將該沙漏濾波器代入如式（3）所示計算，對分解沙漏濾波器進行頻率響應校正，如式（3）所示：
 |
用類似方法，可以獲得沿 ω2 軸和 ω3 軸的沙漏分解濾波器 H2（ω1，ω2，ω3）和 H3（ω1，ω2，ω3）。以上三個濾波器可以彼此旋轉得到，即 H2（ω1，ω2，ω3）=H1（ω2，ω3，ω1），H3（ω1，ω2，ω3）=H1（ω3，ω1，ω2）。如果要實現精確重建，約束條件如式（4）所示：
 |
重建（或合成）沙漏濾波器可以通過對分解濾波器進行坐標翻轉獲得。由于沙漏濾波器關于頂點對稱，所以重建濾波器與分解濾波器相同，如式（5）所示：
 |
第二步，利用二維迭代重采樣棋盤格（iteratively resampled checkerboard，IRC）（符號記為：IRC）濾波器組對三維沙漏濾波器進行頻帶分解，獲得一定數量的金字塔形子帶通道。如果用 
表示待劃分子帶的三維沙漏濾波器，它可以表示為兩個二維濾波器的乘積，如式（6）所示：
 |
其通道數為 
=20 = 1。該沙漏濾波器的子帶分解過程如圖 2 所示。如圖 2 沿 ω1 軸三維沙漏濾波器分解示意圖所示，輸入信號 x（n）經三維沙漏濾波器 
做帶通濾波，將輸出信號 y（n）的頻譜限制在沙漏形的帶通范圍內。然后，采用 IRC 濾波器以級聯的方式進行子帶切分。首先，利用 
沿（n1，n2）平面切分，子帶通道數量為 
個，每個通道對應的輸出信號用za（n）表示（0≤a<
）。然后，對這 
個子帶信號沿（n1，n3）平面切分。如果切分層數為 l3，則輸出的金字塔形子帶通道數量為 
個，每個通道對應的輸出信號用 
表示（0≤b<
）。IRC 由一定層數的二進制樹級聯組成，如圖 2 中 
濾波器結構圖所示。沿（n1，n2）平面切分的 
采用 l2 層的二進制樹級聯分解，由雙通道棋盤格濾波器組（checkerboard filter，CBF）（符號記為：CBF）實現分解，R 為重采樣矩陣，負責在每個通道末尾對信號重采樣。為得到如圖 1 所示的理想金字塔形子帶分解，IRC 濾波器需要滿足以下兩個條件：
圖2
沿 ω1 軸三維沙漏濾波器金字塔形多通道子帶分解過程示意圖
Figure2.
Schematic flow chart of decomposing the three-dimensional hourglass filter along ω1 axis into multi-channel pyramid-shaped sub- bands
（1）等價采樣矩陣。
第 k（0≤k<
）個通道的采樣矩陣 M 需如式（7）所示：
 |
其中，k 為通道數，l2 表示子帶通道的數量為 
個。
（2）楔形優化。設 
為 
第 k 個通道的等價濾波器，則如式（8）所示：
 |
的結構與 
一樣，負責完成沿（n1，n3）平面的子帶切分，切分層數為 l3。以上的分解過程與重建過程是完全對稱的。
1.1.2 Surfacelet 變換
Lu 等[19]將 NDFB 與多尺度分解相結合，提出了 surfacelet 變換。三維 surfacelet 變換及其逆變換如圖 3 所示，其中Li（ω）為低通濾波器，Hi（ω）為高通濾波器，分別用來實現多尺度分解得到對應信號的低頻和高頻部分。↑I 表示上采樣，↓I 表示下采樣。S（ω）為抗混疊濾波器，用來抑制上采樣操作造成的混頻干擾。信號在通過第一層的高通濾波器 H0（ω）后經 3D-DFB 實現多方向分解，通過第一層的低通濾波器 L0（ω）后，先進行 1.5 倍的下采樣操作（即 2 倍上采樣后進行 3 倍下采樣操作，能夠大幅減少頻域混疊現象），后續分解采用 2 倍下采樣，通過重復如圖 3 所示虛線框中的分解過程，對信號實現更多級尺度的分解。此過程中，重建過程與分解過程是完全對稱的。
圖3
三維 surfacelet 變換的分解與重建過程示意圖
Figure3.
Schematic flow chart of the three-dimensional surfacelet transformation and its inverse transformation
多維信號的采樣操作通常定義為一個晶格，一個 N 維信號的晶格是一個大小為 N × N 的非奇異整數矩陣 I ［符號記為：LAT（I）］，如式（9）所示：
 |
其中，n =（n1，n2，
，nN）T 為整數向量，輸入信號 x ［n］經過 I 倍的下采樣（或上采樣）得到的下采樣輸出 yd ［n］（或上采樣輸出 yu ［n］）的表達式如式（10）所示：
 |
低通濾波器 Li（ω）的定義如式（11）所示：
 |
其中，d0 = 6N/2，d1 = 2N/2，N 代表信號的維數，
是一個沿著 ωn 軸的 1 維低通濾波器，如式（12）所示：
 |
其中，
為通帶頻率，當 i = 0，1 時，分別取值 π/3、π/4；
為阻帶頻率，當 i = 0，1 時，分別取值 2π/3、π/2。同樣，抗混疊濾波器 S（ω）的表達式為參數化形式與如式（13）所示相同的一維低通濾波器的可分離乘積，如式（14）所示：
 |
 |
其中，通帶頻率為 
，阻帶頻率為 
。借助理想的多尺度重構條件，高通濾波器 Hi（ω）可以由低通濾波器 Li（ω）獲得[19]，如式（15）所示：
 |
Surfacelet 變換采用非整數的 1.5 倍下采樣雖然能夠減少頻域混疊現象，但是冗余度較高，所以本文直接采用 2 倍下采樣，兩者冗余比約為 4.0∶3.4。本文使用像素為 512×512 的 CT 數據，假設其中一個數據有 64 幀，進行 2 層 surfacelet 分解，高低頻子帶大小分別為 256×256×32。高頻子帶在粗尺度得到 3 個主要方向子帶，每個方向子帶在細尺度內有 16 個不同的子帶，這些子帶分為 128×128×8、128×64×16、64×128×16 三種不同尺寸。
1.2 PCNN 模型
PCNN 模型是在哺乳動物視覺皮層模型的基礎上修改得到的一種人工神經網絡結構。與傳統的神經元相比，PCNN 模型的每個神經元與圖像的像素相對應，不需要任何訓練。前饋神經網絡采用一種單向多層結構，各神經元從輸入層開始，接收前一級輸入，并輸出到下一級，直至輸出層，整個網絡中無反饋；反饋神經網絡每個神經元同時將自身的輸出信號作為輸入反饋給其他神經元，需要工作一段時間才能達到穩定；相比于傳統的前饋神經網絡、反饋神經網絡，PCNN 模型具有同步脈沖發放特性、線性相加和非線性相乘調制耦合等特性，對輸入信息處理能力強，適合于圖像處理[20]。傳統的 PCNN 模型結構比較復雜，控制參數多。本文采用一種基于最小交叉熵的簡化 PCNN 模型，該模型的優點在于能夠利用交叉熵原理自動計算模型的最佳迭代次數，提高了算法的自動化程度[21]。該 PCNN 模型表達式如式（16）～（20）所示：
 |
 |
 |
 |
 |
其中，n 為迭代次數。ij 表示神經元在網絡中位置的下標，kl 表示相鄰神經元的下標。Iij 為外部刺激，即與神經元相關聯的像素的灰度值。ζij（n）為神經元的連接輸入，主要接收來自局部相鄰神經元的輸出信息 Ykl （n ? 1）。VL 表示連接輸入的振幅，W 為 3×3 的內部連接矩陣，表示神經元之間的連接權值。Uij 為內部活動項，β 是連接強度系數，需要手動設置參數。Tij 為動態閾值，α = C/μ，C 為常數，μ 為圖像平均灰度。α 在手動設置參數后，可以根據不同圖像的平均灰度自適應調整。Yij 為神經元的脈沖輸出。可以看出，PCNN 模型的點火像素控制主要由如式（20）所示的動態二值處理過程完成。與傳統的根據圖像灰度分布統計特性選擇全局閾值或局部閾值進行二值化處理不同，PCNN 模型的二值化處理是迭代尋優過程中的一個步驟，閾值是自適應變化的，而且還充分考慮了鄰域像素的灰度值。
PCNN 模型中單個神經元的結構如圖 4 所示。每個神經元均由接受域、調制部分與脈沖發生器組成。其中，接受域接收來自周圍其他神經元的輸入即連接輸入 ζij（n）的信息，調制部分對連接輸入 ζij（n）與饋送輸入 Fij（n）進行調制，得到神經元的內部活動項 Uij，脈沖發生器用于判斷該像素是否點火。從功能上看，PCNN 模型神經元的脈沖發生器與傳統神經網絡中的激勵函數類似，區別在于傳統神經網絡的激勵函數存在于每層結構和輸出層，一般采用非線性函數來表示上層節點的輸出與下層節點的輸入之間的映射關系，而脈沖發生器是根據當前像素與相鄰像素之間的灰度值相近程度判斷該像素對應的神經元是否點火。
圖4
簡化 PCNN 的神經元模型
Figure4.
A neuron model of the simplified PCNN
當分割開始時，需要設置初始閾值（對于 CT 圖像，該閾值一般設置為 255 的灰度值）和迭代次數 n。在第一次迭代過程中，每個神經元接收到與之對應像素灰度值的刺激，此時 ζij 為 0，Uij 不會大于 T0，所以在第一次迭代中沒有神經元脈沖。從第二次迭代開始，Tij（n）以指數形式衰減，隨著迭代次數的增加，使得 Uij（n）大于 Tij（n），因此會有神經元脈沖，脈沖的神經元會產生連接輸入 ζij（n）到相鄰神經元，如果使得相鄰神經元的內部活動項大于閾值，即說明它們具有相似的灰度值，從而產生同步脈沖，并且每個神經元在一個脈沖周期內只能點火一次。但是，在 PCNN 模型的迭代過程中，最優分割結果可能出現在預定迭代次數之前。所以，本文應用最小交叉熵準則來自動選擇最優分割結果[21-22]，在產生最優分割結果時自動停止迭代，可以節省計算時間，提高分割性能。最小交叉熵如式（21）所示：
 |
 |
 |
其中，f 表示圖像灰度值，h（f）表示灰度值為 f 的像素所占總像素的比值，Z 為灰度值上界（255），T = Tij（n），μ1（T）和 μ2（T）為類內均值，分別代表分割結果中目標和背景的平均灰度，其計算公式如式（22）～（23）所示。當 E（f）為最小值時的迭代次數及閾值即為最佳迭代次數 nbest 及最佳閾值 Tbest，而與其對應的輸出 Yij（nbest）即為最佳結果。
1.3 本文算法
本文利用 surfacelet 變換的多尺度、多方向分解特性及局部修正拉普拉斯算子（local sum-modified-laplacian，LSML）增強三維 CT 圖像中的邊緣及細節信息，然后使用 PCNN 算法分割肺實質，算法流程如圖 5 所示，具體步驟如下：
圖5
本文算法流程圖
Figure5.
Flow chart of the proposed algorithm
（1）使用 surfacelet 變換分解三維 CT 肺部數據，得到一系列高頻與低頻子帶，分解的角度級別越高，產生的方向子帶越多。實驗表明，當 l = 2 即產生 3 × 22l = 48 個方向子帶時，性能最優。
（2）對高頻子帶系數采取自適應閾值濾波處理，去除噪聲干擾，設 
為高頻子帶系數，自適應閾值 
的表達式如式（24）所示：
 |
 |
 |
 |
其中，mean 為求平均值運算，c 為常數，經實驗驗證，c = 
時去噪效果最好，
為調整因子，
是由極大似然估計得到的信號方差，M、N、K 分別為三維子帶系數的長寬高，
為噪聲方差，
為第一層分解系數的噪聲方差，通過中值估計法得到，如式（28）所示：
 |
為了降低稀疏系數幅值的波動，對低頻子帶系數 
做歸一化處理，如式（29）所示：
 |
其中，lmax 為低頻子帶系數中的最大值。
（3）LSML（以符號 LSML 表示）能夠增強子帶圖像的邊緣細節信息。利用 LSML 分別對得到的 surfacelet 變換域的低頻和高頻子帶逐幀進行處理，在二維圖像上處理的優點是能夠在增強圖像邊緣特征的同時降低內存限制，LSML 的表達式如式（30）所示：
 |
 |
 |
其中，ML（.）為能量函數，
是三維子帶其中一幀在
處的分解系數，
是一個加權模板，用來突出中心像素，
。
（4）根據如式（15）所示的完全重構條件，經過 surfacelet 逆變換重建過程可以重建得到增強后的圖像 P，并將其作為 PCNN 的輸入。運行 PCNN 簡化模型，計算每次迭代后的二值圖像 Q 與圖像 P 的交叉熵，將得到最小交叉熵時的迭代次數作為 PCNN 的最佳迭代次數，分割肺實質。
（5）基于形態學操作的方法填充肺實質孔洞，如式（33）所示：
 |
其中，X0 是像素值全為 0 的圖像，A 為待填充圖像，AC 為A 的補集，B 為四連通結構元。為了去除氣管區域，對分割結果采用基于八鄰域的區域統計方法，然后去除像素個數小于 1 500 的區域。
（6）基于體繪制方法實現三維重建。
2 實驗與結果分析
2.1 實驗數據
實驗環境：筆記本電腦（Inspiron，Dell Inc.，美國）。
使用的軟件系統：矩陣實驗室軟件 MATLAB R2017b（MathWorks Inc.，美國），微軟 Visual Studio 2013（Microsoft Inc.，美國）。
實驗數據：Kaggle 肺癌檢測挑戰賽上的公開數據集（網址：https://www.kaggle.com/c/data-science-bowl-2017/data）和肺部圖像數據庫聯盟（lung image database consortium，LIDC）公開數據集（網址：https://wiki.cancerimagingarchive.net）。從每個數據集中隨機選取各 10 組胸腔 CT 數據進行實驗，圖像大小為 512×512 像素。
2.2 實驗結果
為了有效評估本文所提算法性能，將其與三維 surfacelet 邊緣提取算法（簡稱：文獻[5]方法）、三維區域生長算法[23]和三維 U-NET 算法進行比較[24]。如圖 6 所示，為 4 幅 CT 肺部切片（來自 Kaggle 數據集，圖片尺寸為 512×512 像素）的分割結果。以手動分割結果為參照，當圖像較為簡單時，上述 4 種算法都能將完整肺實質分割出來，如圖 6 切片 1 所示。當肺實質與氣管存在粘連時，文獻[5]方法與三維區域生長算法誤將氣管分割為肺實質，而三維 U-NET 算法與本文算法能較好地將肺實質分割出來，如圖 6 切片 2 所示。當存在粘連且圖像灰度不均勻，文獻[5]方法誤將氣管分割為肺實質，而其他三種算法分割效果較好，如圖 6 切片 3 所示。當存在肺實質粘連區域，且圖像較模糊時，如圖 6 切片 4 所示，文獻[5]方法與三維區域生長算法在粘連區域分割失敗，無法將左右肺區分開來，三維 U-NET 算法在切片 4 出現了欠分割現象，部分肺實質沒有分割出來，而本文算法能夠分割得到完整肺實質。
圖6
四幅 CT 切片的分割結果比較
Figure6.
Segmentation results comparison of four CT slices
為了進行定量評估，本文從公開數據集 Kaggle 與 LIDC 中選取 20 組三維數據進行實驗。從分割結果中隨機選取 500 張二維 CT 肺部切片進行測試，并采用 3 個定量指標進行性能評估[25]，即 Dice 相似度系數（dice similarity coefficient，DSC）（以符號 DSC 表示）、過分割率（over-segmentation，OS）（以符號 OS 表示）和欠分割率（under-segmentation，US）（以符號 US 表示）。這些指標的定義如式（34）～（36）所示：
 |
其中，R1 表示算法分割結果，R2 表示手工分割結果。
 |
其中，FR 表示手工分割目標的面積或像素數，FS 表示誤分割為目標的面積或像素數。
 |
其中，OR 表示手工分割目標的面積或像素數，OS 表示誤分割為背景的面積或像素數。
4 種算法對 500 張測試圖像的平均 DSC、平均 OS 和平均 US 如表 1 所示。可以看出，本文所提算法的 3 個定量指標均優于其他算法。與文獻[5]方法、三維區域生長算法、三維 U-NET 算法相比，平均 DSC 分別提升 2.7%、1.5%、0.2%，平均 OS 分別減少 1.53%、0.51%、0.17%，平均 US 分別減少 1.1%、1.67%、0.79%。所以，本文所提算法的分割精確性有較明顯的提升。由于三維 U-NET 分割算法需要大量帶標簽數據進行訓練，且訓練時間較長，所以本文所提算法在分割性能相近的情況下，具有計算成本小的優勢。
表1
4 種分割算法的定量評估
Table1.
Quantitative evaluation of four segmentation algorithms
如圖 7 所示，分別為 LIDC 公開數據集中三個肺部體數據（LIDC-IDEI-12，LIDC-IDEI-25，LIDC-IDEI-13）的三維重建結果，每個二維切片的尺寸為 512×512 像素，如圖 7 所示。可以看到，文獻[5]方法和三維區域生長算法受氣管、支氣管的影響，并受到圖像中偽影的干擾，分割效果存在部分失真（三角箭頭所指為誤分割的氣管部分）。三維 U-NET 算法可以較好地避免干擾，但是存在欠分割（圓點示意所指為缺失部分）。而通過與其他算法分割結果比較可以看出，本文算法能夠解決上述問題，得到的分割效果更好。
圖7
LIDC 公開數據集中三個肺部體數據的三維分割與重建結果
Figure7.
Three-dimensional segmentation and reconstruction results of three pulmonary volume data in LIDC
3 討論
由于 PCNN 模型在循環迭代時以圖像的灰度信息為依據，對灰度不均勻的圖像難以區別模糊的目標邊緣，分割質量較差。如圖 8 第一列所示的雙肺 CT 圖像，灰度信息相近的左右肺在結構上非常靠近，二者之間的邊緣比較模糊，如果單純采用 PCNN 模型進行分割會出現大概率的分割失敗。如圖 8 第二列所示，4 幅切片圖像中僅有第 4 行切片 8 的左右肺分割出來，其余三個切片圖像均分割失敗。本文利用改進的 surfacelet 變換算法可以增強圖像的邊緣與細節信息，將增強后的圖像輸入 PCNN 模型進行迭代，對模糊邊緣的目標可以得到更為魯棒有效的分割結果。此外，PCNN 模型分割不同的圖像時都需要設置特定的參數，而一個三維醫學影像數據包含幾十甚至上百張二維圖像，手工設置或調整參數耗時費力。本文所提算法可以根據圖像特征自適應調整 PCNN 模型的動態閾值控制參數［如式（17）所示的 α 參數］，而且能夠在迭代過程中根據最小交叉熵自動選擇最佳分割結果，提高了 PCNN 模型處理三維數據的靈活性。
圖8
本文方法與 PCNN 方法的比較
Figure8.
Comparison of the proposed method with the PCNN
如表 2 所示，為 4 種算法對 500 張切片的平均分割時間。這里，三維 U-NET 算法的計算時間指模型訓練好之后的分割時間，沒有包括數據標注和需要時間（約 51 h）。與文獻[5]方法、三維區域生長算法、三維 U-NET 算法相比，本文算法的計算時間分別增加 1.17、7.45、7.93 s。這是由于 surfacelet 變換對三維數據分解與重建帶來的計算成本增加，但是從分割算法的性能提升角度看，這種適量的計算成本增加是可以接受的。
表2
四種分割算法的分割時間
Table2.
Computational time of the four segmentation algorithms
4 結論
本文提出了一種融合 surfacelet 變換與 PCNN 模型的三維肺實質分割算法。該算法利用 surfacelet 變換的多尺度多方向分解特性增強 CT 圖像的邊緣與細節信息，有效地提升了 PCNN 模型分割肺實質的準確性。而且，PCNN 模型動態閾值控制參數的自動選擇策略，以及利用最小交叉熵自動選擇最佳分割結果可以有效減少人工干預，提高分割的魯棒性和靈活性。通過與 surfacelet 變換邊緣提取算法、三維區域生長分割算法和三維 U-NET 分割算法的實驗比較，表明本文所提分割算法對三維 CT 肺實質分割具有更高的分割性能。
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
引言
從計算機斷層掃描圖像（computed tomography，CT）中準確地分割肺實質是計算機輔助診斷（computer aided diagnosis，CAD）技術的關鍵步驟，對于早期診斷肺癌等肺部疾病具有重要意義[1-2]。另外，肺的邊沿分割精度直接影響肺部體積測量的結果，而肺部體積的測量對研究氣胸、肺氣腫等疾病具有重要意義。肺實質分割方法主要分成三類：第一類是基于傳統方法的肺實質分割，如閾值法[3]、區域生長[4]、邊緣檢測[5]、活動輪廓[6]等，各有一定缺陷。基于閾值的方法對噪聲和成像偽影比較敏感[7]，基于區域生長的方法易出現欠分割或過分割現象，基于邊緣檢測的方法易產生偽邊界或邊界丟失現象[8]。第二類是不同方法相結合的肺實質分割。比如，Nizami 等[9]提出使用小波框架與 K-均值（K-means）聚類相結合分割肺實質的算法；Birkbeck 等[10]提出將統計學習與來自鄰近解剖結構（心臟、肝臟、脾臟和肋骨）的解剖學約束相結合的算法；Zhou 等[11]提出利用增強的閾值算法自動分割肺實質，然后通過紋理感知的活動輪廓模型平滑血管結構和肺部邊界的算法。第三類是基于深度學習的方法。Kumar 等[12]提出基于 U 形網絡（U-NET）的深度卷積神經網絡，可以實現對肺實質的自動分割，該模型采用數據增加（data augmentation）技術，使具有跳接結構的深度網絡能夠根據少量可用的醫學圖像進行像素分類。Kamal 等[13]提出一種融合卷積和遞歸神經網絡的三維混合 U-NET（three dimensional dense UNET，3D-denseUNET）框架，用于分割 CT 圖像中的肺部腫瘤。Harrison 等[14]引入一種基于全卷積神經網絡（fully convolutional neural network，FCN）的漸進式整體嵌套網絡（progressive holistically-nested networks，P-HNNs），用于分割肺組織病變。此外，近幾年出現了將傳統特征提取方法與深度學習方法融合的思路。比如，Ngo 等[15]提出將深度置信網絡與距離正則化水平集相結合，利用深度置信網絡的學習結果驅動零水平集曲線演化進行肺實質的分割。但是，深度學習分割方法的性能依賴于可用的醫學圖像數量，對數據樣本質量和計算能力的要求較高[16]。
本文提出一種融合表面波（surfacelet）變換與脈沖耦合神經網絡（pulse coupled neural network，PCNN）的三維肺實質分割方法，利用 surfacelet 變換和 PCNN 優勢互補以提高肺實質分割的準確性。Surfacelet 變換是一種多分辨率變換，能夠有效提取相鄰幀間的多尺度多方向結構信息，可以通過對變換域各子帶信息進行濾波以增強信號細節[17]。脈沖耦合神經網絡是一種特殊的神經網絡，能夠考慮像素之間灰度的幾何分布選擇合適的歸類，有效分離對象和背景分布之間的重疊，準確提取感興趣區域。所以，利用 surfacelet 變換增強 CT 圖像的邊緣與細節信息，可以彌補 PCNN 只考慮圖像灰度信息的不足，從而能夠減少 CT 圖像中肺部病變及支氣管的干擾，實現肺實質的精確分割。
1 方法
1.1 基于 surfacelet 變換的圖像增強
Surfacelet 變換能夠利用不同尺度、不同頻率的方向子帶準確地捕獲 N 維信號中的面狀奇異，具有平移不變性、完全重構、低冗余、多尺度多方向等性質[18]。Surfacelet 變換的主要思想是將 N 維方向濾波器組（N-dimensional directional filter bank，NDFB）與多尺度金字塔結構結合，利用方向濾波器組實現金字塔形狀的 N 維頻域方向子帶劃分，利用多尺度金字塔結構實現多尺度分析與重構。本文主要考慮三維圖像分割，所以下面介紹 surfacelet 變換的三維形式。
1.1.1 三維方向濾波器組
三維方向濾波器組（three dimensional directional filter bank，3D-DFB）能夠將三維信號的頻譜分解為不同的方向子帶。假設 n1、n2 和 n3 分別表示三維空間域沿垂直方向、水平方向和深度方向相互正交的坐標軸，ω1、ω2 和 ω3 是三維頻域對應的坐標軸，則 3D-DFB 沿 ω1 軸的理想頻率分解示意圖如圖 1 所示[19]，分別表示沿（n1，n2）平面的二維理想楔形頻率子帶劃分、沿（n1，n3）平面的二維理想楔形頻率子帶劃分和金字塔形方向子帶，將兩個二維楔形頻率子帶劃分進行重疊交叉，就會形成頂點對稱的包含多個金字塔形方向子帶（或通道）的沙漏濾波器（hourglass-shaped filter），深灰色部分表示其中的一個金字塔形子帶（pyramid-shaped subband）。類似地，沿 ω2 軸和 ω3 軸的沙漏濾波器可以通過旋轉沙漏濾波器得到。
圖1
沿 ω1 軸的金字塔形頻率子帶分解示意圖
Figure1.
Schematic diagrams of pyramid-shaped frequency sub-band decomposition along ω1 axis
假設 a 表示如圖 1 所示沿（n1，n2）平面的某一層二維濾波器的序號，l2 表示沿（n1，n2）方向的分解層數，b 表示如圖 1 所示沿（n1，n3）平面的某一層二維濾波器的序號，l3 表示沿（n1，n3）方向的分解層數，則如圖 1 所示沿（n1，n2）平面的楔形濾波器可以表示為 （ω1，ω2），0 ≤ a<，沿（n1，n3）平面的楔形濾波器可以表示為 （ω1，ω3），0 ≤ b<。兩個方向的分解層數決定沙漏濾波器包含的子帶數（或通道數），沙漏濾波器可以根據分解層數和層序號表示為 。比如，l2 = 2，l3 = 2，則如圖 1 所示的金字塔形方向子帶共有 =24 = 16 個，其中用深灰色顯示的子帶下標號為 a = 1，b = 1，所以該子帶可以用 表示。如圖 1 所示的重疊交叉操作可以表示為兩個楔形濾波器做乘積運算，如式（1）所示：
|
在具體設計包含多個金字塔形狀方向頻域子帶的沙漏濾波器時，需要兩個步驟：
第一步，設計沙漏濾波器，首先采用平滑凱塞爾窗對二維辛格（sinc）函數進行截取，做二維傅里葉變換分別得到扇形結構的二維濾波器 E（ω1，ω2）和 E（ω1，ω3），然后將二者相乘得到沿 ω1 軸近似理想形狀的沙漏濾波器 K（ω1，ω2，ω3），如式（2）所示：
|
另外，沿 ω2 軸和 ω3 軸的沙漏濾波器 K（ω2，ω3，ω1）和 K（ω3，ω1，ω2），可以由該濾波器旋轉得到。將該沙漏濾波器代入如式（3）所示計算，對分解沙漏濾波器進行頻率響應校正，如式（3）所示：
|
用類似方法，可以獲得沿 ω2 軸和 ω3 軸的沙漏分解濾波器 H2（ω1，ω2，ω3）和 H3（ω1，ω2，ω3）。以上三個濾波器可以彼此旋轉得到，即 H2（ω1，ω2，ω3）=H1（ω2，ω3，ω1），H3（ω1，ω2，ω3）=H1（ω3，ω1，ω2）。如果要實現精確重建，約束條件如式（4）所示：
|
重建（或合成）沙漏濾波器可以通過對分解濾波器進行坐標翻轉獲得。由于沙漏濾波器關于頂點對稱，所以重建濾波器與分解濾波器相同，如式（5）所示：
|
第二步，利用二維迭代重采樣棋盤格（iteratively resampled checkerboard，IRC）（符號記為：IRC）濾波器組對三維沙漏濾波器進行頻帶分解，獲得一定數量的金字塔形子帶通道。如果用 表示待劃分子帶的三維沙漏濾波器，它可以表示為兩個二維濾波器的乘積，如式（6）所示：
|
其通道數為 =20 = 1。該沙漏濾波器的子帶分解過程如圖 2 所示。如圖 2 沿 ω1 軸三維沙漏濾波器分解示意圖所示，輸入信號 x（n）經三維沙漏濾波器 做帶通濾波，將輸出信號 y（n）的頻譜限制在沙漏形的帶通范圍內。然后，采用 IRC 濾波器以級聯的方式進行子帶切分。首先，利用 沿（n1，n2）平面切分，子帶通道數量為 個，每個通道對應的輸出信號用za（n）表示（0≤a<）。然后，對這 個子帶信號沿（n1，n3）平面切分。如果切分層數為 l3，則輸出的金字塔形子帶通道數量為 個，每個通道對應的輸出信號用 表示（0≤b<）。IRC 由一定層數的二進制樹級聯組成，如圖 2 中 濾波器結構圖所示。沿（n1，n2）平面切分的 采用 l2 層的二進制樹級聯分解，由雙通道棋盤格濾波器組（checkerboard filter，CBF）（符號記為：CBF）實現分解，R 為重采樣矩陣，負責在每個通道末尾對信號重采樣。為得到如圖 1 所示的理想金字塔形子帶分解，IRC 濾波器需要滿足以下兩個條件：
圖2
沿 ω1 軸三維沙漏濾波器金字塔形多通道子帶分解過程示意圖
Figure2.
Schematic flow chart of decomposing the three-dimensional hourglass filter along ω1 axis into multi-channel pyramid-shaped sub- bands
（1）等價采樣矩陣。第 k（0≤k<）個通道的采樣矩陣 M 需如式（7）所示：
|
其中，k 為通道數，l2 表示子帶通道的數量為 個。
（2）楔形優化。設 為 第 k 個通道的等價濾波器，則如式（8）所示：
|
的結構與 一樣，負責完成沿（n1，n3）平面的子帶切分，切分層數為 l3。以上的分解過程與重建過程是完全對稱的。
1.1.2 Surfacelet 變換
Lu 等[19]將 NDFB 與多尺度分解相結合，提出了 surfacelet 變換。三維 surfacelet 變換及其逆變換如圖 3 所示，其中Li（ω）為低通濾波器，Hi（ω）為高通濾波器，分別用來實現多尺度分解得到對應信號的低頻和高頻部分。↑I 表示上采樣，↓I 表示下采樣。S（ω）為抗混疊濾波器，用來抑制上采樣操作造成的混頻干擾。信號在通過第一層的高通濾波器 H0（ω）后經 3D-DFB 實現多方向分解，通過第一層的低通濾波器 L0（ω）后，先進行 1.5 倍的下采樣操作（即 2 倍上采樣后進行 3 倍下采樣操作，能夠大幅減少頻域混疊現象），后續分解采用 2 倍下采樣，通過重復如圖 3 所示虛線框中的分解過程，對信號實現更多級尺度的分解。此過程中，重建過程與分解過程是完全對稱的。
圖3
三維 surfacelet 變換的分解與重建過程示意圖
Figure3.
Schematic flow chart of the three-dimensional surfacelet transformation and its inverse transformation
多維信號的采樣操作通常定義為一個晶格，一個 N 維信號的晶格是一個大小為 N × N 的非奇異整數矩陣 I ［符號記為：LAT（I）］，如式（9）所示：
|
其中，n =（n1，n2，，nN）T 為整數向量，輸入信號 x ［n］經過 I 倍的下采樣（或上采樣）得到的下采樣輸出 yd ［n］（或上采樣輸出 yu ［n］）的表達式如式（10）所示：
|
低通濾波器 Li（ω）的定義如式（11）所示：
|
其中，d0 = 6N/2，d1 = 2N/2，N 代表信號的維數， 是一個沿著 ωn 軸的 1 維低通濾波器，如式（12）所示：
|
其中， 為通帶頻率，當 i = 0，1 時，分別取值 π/3、π/4； 為阻帶頻率，當 i = 0，1 時，分別取值 2π/3、π/2。同樣，抗混疊濾波器 S（ω）的表達式為參數化形式與如式（13）所示相同的一維低通濾波器的可分離乘積，如式（14）所示：
|
|
其中，通帶頻率為 ，阻帶頻率為 。借助理想的多尺度重構條件，高通濾波器 Hi（ω）可以由低通濾波器 Li（ω）獲得[19]，如式（15）所示：
|
Surfacelet 變換采用非整數的 1.5 倍下采樣雖然能夠減少頻域混疊現象，但是冗余度較高，所以本文直接采用 2 倍下采樣，兩者冗余比約為 4.0∶3.4。本文使用像素為 512×512 的 CT 數據，假設其中一個數據有 64 幀，進行 2 層 surfacelet 分解，高低頻子帶大小分別為 256×256×32。高頻子帶在粗尺度得到 3 個主要方向子帶，每個方向子帶在細尺度內有 16 個不同的子帶，這些子帶分為 128×128×8、128×64×16、64×128×16 三種不同尺寸。
1.2 PCNN 模型
PCNN 模型是在哺乳動物視覺皮層模型的基礎上修改得到的一種人工神經網絡結構。與傳統的神經元相比，PCNN 模型的每個神經元與圖像的像素相對應，不需要任何訓練。前饋神經網絡采用一種單向多層結構，各神經元從輸入層開始，接收前一級輸入，并輸出到下一級，直至輸出層，整個網絡中無反饋；反饋神經網絡每個神經元同時將自身的輸出信號作為輸入反饋給其他神經元，需要工作一段時間才能達到穩定；相比于傳統的前饋神經網絡、反饋神經網絡，PCNN 模型具有同步脈沖發放特性、線性相加和非線性相乘調制耦合等特性，對輸入信息處理能力強，適合于圖像處理[20]。傳統的 PCNN 模型結構比較復雜，控制參數多。本文采用一種基于最小交叉熵的簡化 PCNN 模型，該模型的優點在于能夠利用交叉熵原理自動計算模型的最佳迭代次數，提高了算法的自動化程度[21]。該 PCNN 模型表達式如式（16）～（20）所示：
|
|
|
|
|
其中，n 為迭代次數。ij 表示神經元在網絡中位置的下標，kl 表示相鄰神經元的下標。Iij 為外部刺激，即與神經元相關聯的像素的灰度值。ζij（n）為神經元的連接輸入，主要接收來自局部相鄰神經元的輸出信息 Ykl （n ? 1）。VL 表示連接輸入的振幅，W 為 3×3 的內部連接矩陣，表示神經元之間的連接權值。Uij 為內部活動項，β 是連接強度系數，需要手動設置參數。Tij 為動態閾值，α = C/μ，C 為常數，μ 為圖像平均灰度。α 在手動設置參數后，可以根據不同圖像的平均灰度自適應調整。Yij 為神經元的脈沖輸出。可以看出，PCNN 模型的點火像素控制主要由如式（20）所示的動態二值處理過程完成。與傳統的根據圖像灰度分布統計特性選擇全局閾值或局部閾值進行二值化處理不同，PCNN 模型的二值化處理是迭代尋優過程中的一個步驟，閾值是自適應變化的，而且還充分考慮了鄰域像素的灰度值。
PCNN 模型中單個神經元的結構如圖 4 所示。每個神經元均由接受域、調制部分與脈沖發生器組成。其中，接受域接收來自周圍其他神經元的輸入即連接輸入 ζij（n）的信息，調制部分對連接輸入 ζij（n）與饋送輸入 Fij（n）進行調制，得到神經元的內部活動項 Uij，脈沖發生器用于判斷該像素是否點火。從功能上看，PCNN 模型神經元的脈沖發生器與傳統神經網絡中的激勵函數類似，區別在于傳統神經網絡的激勵函數存在于每層結構和輸出層，一般采用非線性函數來表示上層節點的輸出與下層節點的輸入之間的映射關系，而脈沖發生器是根據當前像素與相鄰像素之間的灰度值相近程度判斷該像素對應的神經元是否點火。
圖4
簡化 PCNN 的神經元模型
Figure4.
A neuron model of the simplified PCNN
當分割開始時，需要設置初始閾值（對于 CT 圖像，該閾值一般設置為 255 的灰度值）和迭代次數 n。在第一次迭代過程中，每個神經元接收到與之對應像素灰度值的刺激，此時 ζij 為 0，Uij 不會大于 T0，所以在第一次迭代中沒有神經元脈沖。從第二次迭代開始，Tij（n）以指數形式衰減，隨著迭代次數的增加，使得 Uij（n）大于 Tij（n），因此會有神經元脈沖，脈沖的神經元會產生連接輸入 ζij（n）到相鄰神經元，如果使得相鄰神經元的內部活動項大于閾值，即說明它們具有相似的灰度值，從而產生同步脈沖，并且每個神經元在一個脈沖周期內只能點火一次。但是，在 PCNN 模型的迭代過程中，最優分割結果可能出現在預定迭代次數之前。所以，本文應用最小交叉熵準則來自動選擇最優分割結果[21-22]，在產生最優分割結果時自動停止迭代，可以節省計算時間，提高分割性能。最小交叉熵如式（21）所示：
|
|
|
其中，f 表示圖像灰度值，h（f）表示灰度值為 f 的像素所占總像素的比值，Z 為灰度值上界（255），T = Tij（n），μ1（T）和 μ2（T）為類內均值，分別代表分割結果中目標和背景的平均灰度，其計算公式如式（22）～（23）所示。當 E（f）為最小值時的迭代次數及閾值即為最佳迭代次數 nbest 及最佳閾值 Tbest，而與其對應的輸出 Yij（nbest）即為最佳結果。
1.3 本文算法
本文利用 surfacelet 變換的多尺度、多方向分解特性及局部修正拉普拉斯算子（local sum-modified-laplacian，LSML）增強三維 CT 圖像中的邊緣及細節信息，然后使用 PCNN 算法分割肺實質，算法流程如圖 5 所示，具體步驟如下：
圖5
本文算法流程圖
Figure5.
Flow chart of the proposed algorithm
（1）使用 surfacelet 變換分解三維 CT 肺部數據，得到一系列高頻與低頻子帶，分解的角度級別越高，產生的方向子帶越多。實驗表明，當 l = 2 即產生 3 × 22l = 48 個方向子帶時，性能最優。
（2）對高頻子帶系數采取自適應閾值濾波處理，去除噪聲干擾，設 為高頻子帶系數，自適應閾值 的表達式如式（24）所示：
|
|
|
|
其中，mean 為求平均值運算，c 為常數，經實驗驗證，c = 時去噪效果最好， 為調整因子， 是由極大似然估計得到的信號方差，M、N、K 分別為三維子帶系數的長寬高， 為噪聲方差， 為第一層分解系數的噪聲方差，通過中值估計法得到，如式（28）所示：
|
為了降低稀疏系數幅值的波動，對低頻子帶系數 做歸一化處理，如式（29）所示：
|
其中，lmax 為低頻子帶系數中的最大值。
（3）LSML（以符號 LSML 表示）能夠增強子帶圖像的邊緣細節信息。利用 LSML 分別對得到的 surfacelet 變換域的低頻和高頻子帶逐幀進行處理，在二維圖像上處理的優點是能夠在增強圖像邊緣特征的同時降低內存限制，LSML 的表達式如式（30）所示：
|
|
|
其中，ML（.）為能量函數，是三維子帶其中一幀在處的分解系數，是一個加權模板，用來突出中心像素，。
（4）根據如式（15）所示的完全重構條件，經過 surfacelet 逆變換重建過程可以重建得到增強后的圖像 P，并將其作為 PCNN 的輸入。運行 PCNN 簡化模型，計算每次迭代后的二值圖像 Q 與圖像 P 的交叉熵，將得到最小交叉熵時的迭代次數作為 PCNN 的最佳迭代次數，分割肺實質。
（5）基于形態學操作的方法填充肺實質孔洞，如式（33）所示：
|
其中，X0 是像素值全為 0 的圖像，A 為待填充圖像，AC 為A 的補集，B 為四連通結構元。為了去除氣管區域，對分割結果采用基于八鄰域的區域統計方法，然后去除像素個數小于 1 500 的區域。
（6）基于體繪制方法實現三維重建。
2 實驗與結果分析
2.1 實驗數據
實驗環境：筆記本電腦（Inspiron，Dell Inc.，美國）。
使用的軟件系統：矩陣實驗室軟件 MATLAB R2017b（MathWorks Inc.，美國），微軟 Visual Studio 2013（Microsoft Inc.，美國）。
實驗數據：Kaggle 肺癌檢測挑戰賽上的公開數據集（網址：https://www.kaggle.com/c/data-science-bowl-2017/data）和肺部圖像數據庫聯盟（lung image database consortium，LIDC）公開數據集（網址：https://wiki.cancerimagingarchive.net）。從每個數據集中隨機選取各 10 組胸腔 CT 數據進行實驗，圖像大小為 512×512 像素。
2.2 實驗結果
為了有效評估本文所提算法性能，將其與三維 surfacelet 邊緣提取算法（簡稱：文獻[5]方法）、三維區域生長算法[23]和三維 U-NET 算法進行比較[24]。如圖 6 所示，為 4 幅 CT 肺部切片（來自 Kaggle 數據集，圖片尺寸為 512×512 像素）的分割結果。以手動分割結果為參照，當圖像較為簡單時，上述 4 種算法都能將完整肺實質分割出來，如圖 6 切片 1 所示。當肺實質與氣管存在粘連時，文獻[5]方法與三維區域生長算法誤將氣管分割為肺實質，而三維 U-NET 算法與本文算法能較好地將肺實質分割出來，如圖 6 切片 2 所示。當存在粘連且圖像灰度不均勻，文獻[5]方法誤將氣管分割為肺實質，而其他三種算法分割效果較好，如圖 6 切片 3 所示。當存在肺實質粘連區域，且圖像較模糊時，如圖 6 切片 4 所示，文獻[5]方法與三維區域生長算法在粘連區域分割失敗，無法將左右肺區分開來，三維 U-NET 算法在切片 4 出現了欠分割現象，部分肺實質沒有分割出來，而本文算法能夠分割得到完整肺實質。
圖6
四幅 CT 切片的分割結果比較
Figure6.
Segmentation results comparison of four CT slices
為了進行定量評估，本文從公開數據集 Kaggle 與 LIDC 中選取 20 組三維數據進行實驗。從分割結果中隨機選取 500 張二維 CT 肺部切片進行測試，并采用 3 個定量指標進行性能評估[25]，即 Dice 相似度系數（dice similarity coefficient，DSC）（以符號 DSC 表示）、過分割率（over-segmentation，OS）（以符號 OS 表示）和欠分割率（under-segmentation，US）（以符號 US 表示）。這些指標的定義如式（34）～（36）所示：
|
其中，R1 表示算法分割結果，R2 表示手工分割結果。
|
其中，FR 表示手工分割目標的面積或像素數，FS 表示誤分割為目標的面積或像素數。
|
其中，OR 表示手工分割目標的面積或像素數，OS 表示誤分割為背景的面積或像素數。
4 種算法對 500 張測試圖像的平均 DSC、平均 OS 和平均 US 如表 1 所示。可以看出，本文所提算法的 3 個定量指標均優于其他算法。與文獻[5]方法、三維區域生長算法、三維 U-NET 算法相比，平均 DSC 分別提升 2.7%、1.5%、0.2%，平均 OS 分別減少 1.53%、0.51%、0.17%，平均 US 分別減少 1.1%、1.67%、0.79%。所以，本文所提算法的分割精確性有較明顯的提升。由于三維 U-NET 分割算法需要大量帶標簽數據進行訓練，且訓練時間較長，所以本文所提算法在分割性能相近的情況下，具有計算成本小的優勢。
表1
4 種分割算法的定量評估
Table1.
Quantitative evaluation of four segmentation algorithms
如圖 7 所示，分別為 LIDC 公開數據集中三個肺部體數據（LIDC-IDEI-12，LIDC-IDEI-25，LIDC-IDEI-13）的三維重建結果，每個二維切片的尺寸為 512×512 像素，如圖 7 所示。可以看到，文獻[5]方法和三維區域生長算法受氣管、支氣管的影響，并受到圖像中偽影的干擾，分割效果存在部分失真（三角箭頭所指為誤分割的氣管部分）。三維 U-NET 算法可以較好地避免干擾，但是存在欠分割（圓點示意所指為缺失部分）。而通過與其他算法分割結果比較可以看出，本文算法能夠解決上述問題，得到的分割效果更好。
圖7
LIDC 公開數據集中三個肺部體數據的三維分割與重建結果
Figure7.
Three-dimensional segmentation and reconstruction results of three pulmonary volume data in LIDC
3 討論
由于 PCNN 模型在循環迭代時以圖像的灰度信息為依據，對灰度不均勻的圖像難以區別模糊的目標邊緣，分割質量較差。如圖 8 第一列所示的雙肺 CT 圖像，灰度信息相近的左右肺在結構上非常靠近，二者之間的邊緣比較模糊，如果單純采用 PCNN 模型進行分割會出現大概率的分割失敗。如圖 8 第二列所示，4 幅切片圖像中僅有第 4 行切片 8 的左右肺分割出來，其余三個切片圖像均分割失敗。本文利用改進的 surfacelet 變換算法可以增強圖像的邊緣與細節信息，將增強后的圖像輸入 PCNN 模型進行迭代，對模糊邊緣的目標可以得到更為魯棒有效的分割結果。此外，PCNN 模型分割不同的圖像時都需要設置特定的參數，而一個三維醫學影像數據包含幾十甚至上百張二維圖像，手工設置或調整參數耗時費力。本文所提算法可以根據圖像特征自適應調整 PCNN 模型的動態閾值控制參數［如式（17）所示的 α 參數］，而且能夠在迭代過程中根據最小交叉熵自動選擇最佳分割結果，提高了 PCNN 模型處理三維數據的靈活性。
圖8
本文方法與 PCNN 方法的比較
Figure8.
Comparison of the proposed method with the PCNN
如表 2 所示，為 4 種算法對 500 張切片的平均分割時間。這里，三維 U-NET 算法的計算時間指模型訓練好之后的分割時間，沒有包括數據標注和需要時間（約 51 h）。與文獻[5]方法、三維區域生長算法、三維 U-NET 算法相比，本文算法的計算時間分別增加 1.17、7.45、7.93 s。這是由于 surfacelet 變換對三維數據分解與重建帶來的計算成本增加，但是從分割算法的性能提升角度看，這種適量的計算成本增加是可以接受的。
表2
四種分割算法的分割時間
Table2.
Computational time of the four segmentation algorithms
4 結論
本文提出了一種融合 surfacelet 變換與 PCNN 模型的三維肺實質分割算法。該算法利用 surfacelet 變換的多尺度多方向分解特性增強 CT 圖像的邊緣與細節信息，有效地提升了 PCNN 模型分割肺實質的準確性。而且，PCNN 模型動態閾值控制參數的自動選擇策略，以及利用最小交叉熵自動選擇最佳分割結果可以有效減少人工干預，提高分割的魯棒性和靈活性。通過與 surfacelet 變換邊緣提取算法、三維區域生長分割算法和三維 U-NET 分割算法的實驗比較，表明本文所提分割算法對三維 CT 肺實質分割具有更高的分割性能。
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
