張欣，徐永瀟，王兵

(1.河北大學(xué) 電子信息工程學(xué)院，河北 保定 071002；2. 河北大學(xué) 數(shù)學(xué)與信息科學(xué)學(xué)院，河北 保定 071002)

肺癌目前已經(jīng)是致死率非常高的惡性腫瘤，而且肺癌的發(fā)病率也在不斷上升[1].計算機(jī)斷層掃描技術(shù)(computed tomography,CT)和計算機(jī)輔助診斷(computer aided diagnosis,CAD)為肺癌的診斷提供了更多的方法[2].為了提高醫(yī)生對腫瘤診斷的精確性和可靠性，需要對肺及肺部病變部位進(jìn)行準(zhǔn)確的分析.有些CT圖像中腫瘤粘連著肺邊界，其密度與胸廓或心臟等部位的密度相差很小,并且腫瘤的形狀、大小和粘連的位置都存在很大的差異[3]，導(dǎo)致對肺實質(zhì)與腫瘤的分割沒有一個全面有效的方法.

針對肺實質(zhì)的分割與邊界修復(fù)問題，多采用閾值法[4]、聚類法[5]、邊緣檢測法和區(qū)域生長[6]等方法對CT肺圖像進(jìn)行粗分割，然后對肺邊界粘連型腫瘤造成的肺邊界缺陷部分進(jìn)行修復(fù).目前的修復(fù)方法有滾球法、鏈碼法[7]、凸包算法[8]等.文獻(xiàn)[5]提出了一種基于灰度積分投影與模糊C均值聚類并結(jié)合滾球法修復(fù)缺陷部分的分割方法，但滾球半徑的選擇對肺邊界修復(fù)效果影響較大.文獻(xiàn)[7]提出了一種基于改進(jìn)鏈碼和Bresenham算法結(jié)合的缺陷肺實質(zhì)修復(fù)方法，但該算法容易出現(xiàn)過分割現(xiàn)象.凸包算法的修復(fù)效果受到缺陷處形狀的影響較大，當(dāng)需要修復(fù)的肺實質(zhì)邊界是曲率變化較小的平滑曲線時，凸包算法可以作為首要的修復(fù)方法.代雙鳳等[8]提出一種基于3D區(qū)域增長和凸包算法相結(jié)合的分割方法，但是針對缺陷位于二維肺實質(zhì)邊界曲率變化比較大處的情況，例如當(dāng)缺陷位于肺葉尖端時，二維肺部CT圖像中缺少指導(dǎo)肺邊界缺陷修復(fù)的信息，效果不理想.因此，本文提出一種基于三維曲面重建[9]的方法，實驗證明該方法用于修復(fù)肺實質(zhì)邊界曲率變化較大處的缺陷效果良好，修復(fù)后的肺實質(zhì)有助于準(zhǔn)確分割CT圖像中的肺腫瘤.

1 算法描述與分析

基于三維曲面重建算法的肺邊界缺陷修復(fù)方法主要由3部分組成：CT圖像提取肺實質(zhì)、肺實質(zhì)的三維重建、修復(fù)具有缺陷的三維肺實質(zhì)，具體流程如圖1所示.

圖1 肺實質(zhì)修復(fù)流程Fig.1 Flow chart of pulmonary parenchyma restore

1.1 CT圖像提取肺實質(zhì)

從CT圖像中提取出肺實質(zhì)區(qū)域：首先將肺部CT圖像進(jìn)行格式和灰度轉(zhuǎn)換，再進(jìn)行高斯濾波去除噪聲，得到較為平滑的人體肺部CT圖像；然后為了更加快速地提取序列CT圖像中的肺實質(zhì)區(qū)域，使用質(zhì)心灰度法改進(jìn)區(qū)域生長算法；最后通過形態(tài)學(xué)操作填補肺結(jié)節(jié)和氣管等造成的空洞，得到肺實質(zhì)區(qū)域.

由于傳統(tǒng)的區(qū)域生長算法對每一張二維圖像都需要人為選取種子點，操作起來較為繁瑣.本文使用相鄰層生長區(qū)域的質(zhì)心位置結(jié)合鄰域灰度閾值篩選的方法自動確定當(dāng)前層切片圖像的生長種子點，適用于從肺部序列CT圖像中提取肺實質(zhì).具體步驟如下.

步驟1：在某層肺實質(zhì)區(qū)域內(nèi)選取1個或多個種子點，標(biāo)記為已生長區(qū)域R[t]，令t為迭代次數(shù).

步驟2：設(shè)平均灰度值即為生長區(qū)域的相似性質(zhì)，計算已生長區(qū)域R[t]的平均灰度值

(1)

其中amount(R[t])表示R[t]內(nèi)像素點的數(shù)量，I(y)表示像素y點的灰度值.

步驟5：在相鄰切片上選擇qc+1的8鄰域像素點中滿足min|I(Qc+1)-μc|條件的像素點qo,其中Qc+1為qc+1的8鄰域像素點集，I(Qc+1)表示點集Qc+1中像素點的灰度值.將qo作為該層切片的初始種子點，重復(fù)步驟2～5，直至到達(dá)最底層切片中肺部消失.此時該切片上沒有肺實質(zhì)，中心點的灰度值與相鄰切片的平均灰度值相差較大，設(shè)置合理的終止閾值φ，當(dāng)無法滿足條件I(qo)<φ時停止生長，即得到完整序列的二維肺實質(zhì)圖像.

1.2 肺實質(zhì)的三維重建

因為后續(xù)缺陷修復(fù)需要三維曲面信息，故先采用面繪制中最具有代表性的MC(marching cubes)算法[10]進(jìn)行肺實質(zhì)的三維重建，即通過閾值判斷提取三維模型表面的灰度等值面實現(xiàn)重建.再使用拉普拉斯平滑算法對生成的三維模型進(jìn)行平滑處理，最終得到表面光滑的三維肺實質(zhì).具體步驟如下.

步驟1：讀取經(jīng)過預(yù)處理的二維肺實質(zhì)圖像.在三維數(shù)據(jù)場中掃描所有空間像素點并構(gòu)建體素(體素是由相鄰的2個切片垂直對應(yīng)的8個像素點組成的立方體).由于灰度等值面與不同體素相交后，體素內(nèi)部會形成不同的剖分面，將立方體中所有可能產(chǎn)生的剖分面進(jìn)行分類，創(chuàng)建立方體構(gòu)型的索引表用于響應(yīng)體素.

步驟2：提取表面體素.與等值面相交的體素為表面體素，即體素棱邊與等值面的交點大于等于1.通過線性插值法求出上述交點即為等值面片的頂點，計算出等值面片頂點的法向量用以提供三維模型的表面光照變化信息，使得視覺效果更好.

步驟3：將等值面片進(jìn)行連接，選定光照模型并進(jìn)行平滑處理，完成三維肺實質(zhì)的重建.

1.3 修復(fù)肺實質(zhì)缺陷

根據(jù)三維曲面信息修復(fù)肺實質(zhì)邊界缺陷的方法分為：點云數(shù)據(jù)集構(gòu)建和Possion三維曲面重建.

1.3.1 構(gòu)建點云數(shù)據(jù)集

在肺實質(zhì)缺陷處周圍的表面進(jìn)行采樣，將所得點集記為P.具體步驟如下.

步驟1：采用凸包理論中常用的Graham掃描法和閾值判斷法相結(jié)合的方法[8]找出序列CT圖像中肺邊界的缺陷部位.然后將缺陷部位擴(kuò)大τ個像素的距離并采集二維缺陷的邊界信息，將所得信息映射到重建的三維肺實質(zhì)上，即可得到肺實質(zhì)缺陷處的三維輪廓點集.對三維點集進(jìn)行插值和連接，最終得到肺實質(zhì)缺陷處的三維輪廓曲線.

步驟2：計算三維輪廓點集中質(zhì)心與最遠(yuǎn)點的距離r.構(gòu)建一個以質(zhì)心為球心、2r為半徑的球體，對球體區(qū)域與缺陷三維肺實質(zhì)區(qū)域進(jìn)行布爾運算中的“交運算”，獲得感興趣區(qū)域ROI(region of interest).針對ROI區(qū)域進(jìn)行操作可以提高運算效率.

步驟3：采集ROI區(qū)域的點云數(shù)據(jù)，刪除位于ROI區(qū)域的缺陷輪廓曲線內(nèi)點云，得到點云數(shù)據(jù)集P.

1.3.2 Possion三維曲面重建算法修復(fù)肺實質(zhì)

Possion曲面重建算法是一種隱式曲面重建算法[11]，具有局部和全局?jǐn)M合的特點，魯棒性強(qiáng)，可以使重建出來的肺實質(zhì)三維輪廓更精準(zhǔn)平滑.該算法的核心就是構(gòu)建空間元素的向量場來近似模擬三維肺實質(zhì)表面的指示函數(shù)[9]，將向量場的求解轉(zhuǎn)化為求解Possion等式的問題，通過求解的Possion等式得到曲面重建所需要的等值面.具體步驟如下.

步驟1：節(jié)點混合函數(shù)的定義.首先對P建立一個節(jié)點為o的八叉樹?.然后令基函數(shù)F∶R3→R表示函數(shù)空間，該基函數(shù)具有節(jié)點積分和尺度縮放性質(zhì).所得函數(shù)空間具有多尺度結(jié)構(gòu)[9].設(shè)q為函數(shù)空間中的變量，為函數(shù)空間中的?的每個節(jié)點o添加一個對應(yīng)于單位積分的混合函數(shù)Fo，將混合函數(shù)以節(jié)點o的中心和寬度進(jìn)行展開可得

(2)

其中o.c和o.w為節(jié)點的中心和寬度.所有節(jié)點的混合函數(shù)Fo線性求和可以表示向量場V.

步驟2：基函數(shù)的定義.采用盒式濾波器的n維卷積近似基函數(shù)F，該方法能夠使節(jié)點函數(shù)Fo(q)的線性求和更加精確的表示向量場V，基函數(shù)

(3)

步驟3：向量場的構(gòu)造.先針對當(dāng)前節(jié)點臨近的8個節(jié)點采用三次線性插值的方法進(jìn)行分配，然后再構(gòu)造梯度域函數(shù)的向量場

(4)

對點云數(shù)據(jù)P進(jìn)行平均采樣，所得樣本數(shù)據(jù)s∈P，其中，樣本數(shù)據(jù)中的點s.p∈s，s.N為s的表面內(nèi)法向量，NgbrD(s)表示深度為D的八叉樹中最靠近s.p的節(jié)點區(qū)域，αo,s為權(quán)重系數(shù).因為三維模型表面都存在一個幾乎不會發(fā)生變化的指示函數(shù)，可以看做一個常量.因此指示函數(shù)的梯度在任何除了表面的節(jié)點上都為0，即可以用上面構(gòu)造的向量場近似估計三維肺實質(zhì)ROI區(qū)域表面.

(5)

步驟5：等值面的提取.通過均值法計算等值，再從指示函數(shù)中提取等值面，使用MC算法對等值面進(jìn)行曲面重建即可將三維肺實質(zhì)ROI區(qū)域的表面修復(fù)完整.對修復(fù)面與帶有缺陷的肺實質(zhì)進(jìn)行“并運算”即可將肺實質(zhì)修復(fù)完好.最后將切片輪廓從修復(fù)完好的肺實質(zhì)中提取出來，進(jìn)行填充處理后制成掩模，該掩膜用于從原始CT圖像中分割肺實質(zhì).

2 實驗分析

實驗數(shù)據(jù)共選取了10組肺部帶有粘連型腫瘤的三維CT圖像數(shù)據(jù).一部分?jǐn)?shù)據(jù)來源于河北大學(xué)附屬醫(yī)院CT室，數(shù)據(jù)格式為DICOM，CT圖像掃描厚度為3 mm，分辨率為512×512像素. 另一部分來源于美國國家癌癥研究中心公開的DSB2017肺癌預(yù)測數(shù)據(jù)集(data science bowl)，CT圖像的格式為DICOM，掃描厚度為2 mm，分辨率為512×512像素.實驗硬件平臺基于Windows10操作系統(tǒng)，Intel(R) Core(TM) i5處理器，8 G內(nèi)存， 2.50 GHz主頻，軟件平臺基于Matlab R2016a，VTK，Visual Studio 2017.

2.1 應(yīng)用本文方法修復(fù)肺實質(zhì)缺陷

基于Possion三維曲面重建算法的肺實質(zhì)修復(fù)方法實驗過程效果如圖2所示.首先選取帶有邊界粘連型腫瘤的肺部CT圖像，缺陷位于左肺上半部分的尖端處標(biāo)注的地方，如圖2a所示.使用改進(jìn)的區(qū)域生長法提取肺部序列圖像中的肺實質(zhì)圖像，經(jīng)過多次實驗對比后設(shè)置生長閾值θ=25,φ=180.再使用形態(tài)學(xué)操作填補肺實質(zhì)中由肺結(jié)節(jié)和血管等造成的孔洞，得到肺實質(zhì)的掩模圖像如圖2b所示.通過掩模得到肺實質(zhì)的初步分割結(jié)果圖如2c所示，從該結(jié)果可以明顯看到缺陷部位.然后使用MC算法對所有切片初分割后的肺實質(zhì)進(jìn)行三維重建，如圖2d所示.提取二維肺實質(zhì)缺陷邊界并映射到三維空間，經(jīng)過多次實驗對比后設(shè)置τ=3.構(gòu)建球體后提取ROI區(qū)域，同時標(biāo)記缺陷輪廓邊界，如圖2e所示.再對肺實質(zhì)的缺陷部位進(jìn)行重采樣獲得點云，如圖2f所示.使用本文的Possion曲面重建算法對缺陷部位進(jìn)行三維曲面重建，依據(jù)文獻(xiàn)[9]中的經(jīng)驗值設(shè)置八叉樹的最大深度D=8，如圖2g所示.將重建后的缺陷部位與原始缺陷肺實質(zhì)進(jìn)行結(jié)合，結(jié)果如圖2h所示.最后從修復(fù)完整的肺實質(zhì)中提取二維切片，如圖2i所示.對二維切面進(jìn)行填充獲得掩模后再對CT圖像中肺實質(zhì)進(jìn)行分割，得到的最終的肺實質(zhì)分割結(jié)果如圖2j所示.從分割結(jié)果可以看出，采用本文的修復(fù)方法可以修復(fù)標(biāo)注的肺實質(zhì)缺陷，且修復(fù)的邊緣更加平滑，同時將腫瘤完整的包含在了肺實質(zhì)中.

a.原始肺部CT圖像；b.區(qū)域生長得到掩摸；c.通過掩摸提取肺實質(zhì)；d.三維重建肺實質(zhì)；e.提取缺陷部位；f.重采樣獲得點云；g.三維重建缺陷部位；h.修復(fù)完整的肺實質(zhì)；i.提取切片處輪廓；j.最終分割結(jié)果.圖2 本文方法修復(fù)肺實質(zhì)的實驗效果Fig.2 Experimental effect of restoring lung parenchyma with the method used in this paper

2.2 本文的修復(fù)方法與經(jīng)典的“滾球法”和“凸包算法”相比較

在肺實質(zhì)的修復(fù)方法中，最經(jīng)典的2種方法分別是“滾球法”和“凸包算法”.它們具有計算速度快、原理簡單易操作的優(yōu)點，但是在修復(fù)效果上都存在著一定的缺陷.將本文所使用的修復(fù)方法與滾球法和凸包算法應(yīng)用到同一張切片上比較，修復(fù)效果如圖3所示.選用的切片原圖為圖3a，從圖3b和圖3c可以看出使用滾球法修復(fù)肺實質(zhì)的效果很依賴滾球半徑.當(dāng)滾球半徑選的較小時，如r=5，會出現(xiàn)肺實質(zhì)的“欠分割”，如圖3b中圈起來的部分所示，缺陷并沒有得到修復(fù).當(dāng)滾球半徑選的較大時，如r=15，從圖3c中上圓標(biāo)注的位置可以看出缺陷處得到一定程度的修復(fù)，但是在肺門處卻出現(xiàn)了“過分割”.從圖3d中可以看出使用凸包算法修復(fù)缺陷肺實質(zhì)時有一定的局限性.對肺實質(zhì)邊界曲率變化較大處的缺陷，由于二維圖像中可用于修復(fù)缺陷的參考特征信息較少，所以凸包算法的修復(fù)效果差.此外，從標(biāo)記的地方可以看出在其他正常邊界處容易出現(xiàn)“過分割”的現(xiàn)象.而本文采取的修復(fù)方法可以將粘連型腫瘤導(dǎo)致的肺實質(zhì)缺失部分補充出來，修復(fù)后的肺實質(zhì)可以完整的包含腫瘤，且邊緣連續(xù)性和平滑性都較好，如圖3e所示.從修復(fù)方法的復(fù)雜程度來看，滾球法和凸包算法的過程都較為簡單，而且僅在二維圖像上即可完成，算法的運行速度也較快.本文所使用的方法運算過程較為復(fù)雜，需要三維建模等較為耗時的運算，但可以直接從三維角度觀察肺實質(zhì)，使得修復(fù)后的效果更加直觀清晰.

a.原圖;b.r=5，滾球法;c.r=15，滾球法;d.凸包算法;e.本文方法.圖3 本文的修復(fù)方法與經(jīng)典的“滾球法”和“凸包算法”相比較Fig.3 Patching algorithm in this paper is compared with the classical “ball rolling method” and “convex hull algorithm”

2.3 分割正確率

使用本文方法對選取的10組肺部CT數(shù)據(jù)進(jìn)行全肺分割，計算分割正確率

.

(6)

其中正確分割的圖像由放射科醫(yī)生進(jìn)行判斷.經(jīng)過實驗統(tǒng)計計算，表1為應(yīng)用本文方法和應(yīng)用滾球法與凸包算法對全肺進(jìn)行分割結(jié)果的正確率統(tǒng)計表.從統(tǒng)計結(jié)果可以看出，應(yīng)用本文的方法對肺實質(zhì)進(jìn)行分割的正確率優(yōu)于滾球法和凸包算法.對失敗的幾個案例進(jìn)行分析可知其主要原因是粘連的腫瘤過大，導(dǎo)致肺實質(zhì)的邊界缺陷過大，無法進(jìn)行三維重建.這種情況下繼續(xù)使用本文算法則會出現(xiàn)欠分割的現(xiàn)象，所以對于有大腫瘤粘連的肺實質(zhì)的修復(fù)方法需要進(jìn)一步研究.

表1 不同分割方法的正確率比較

3 結(jié)束語

在實際的計算機(jī)輔助臨床診斷中，準(zhǔn)確分割肺實質(zhì)對后續(xù)肺腫瘤的分割以及肺腫瘤的檢測和分析等都具有十分重要的意義.本文針對肺實質(zhì)邊界曲率變化較大處因腫瘤粘連導(dǎo)致的肺實質(zhì)邊界缺陷，由于該位置處的缺陷在二維圖像上缺少用于指導(dǎo)修復(fù)的特征信息，所以經(jīng)典的凸包算法或滾球法的修復(fù)效果較差，故提出一種基于三維曲面重建的修復(fù)方法.從實驗結(jié)果看出，該方法能夠通過提取肺實質(zhì)缺陷處的三維表面信息對修復(fù)曲面進(jìn)行指導(dǎo)，修復(fù)的效果更理想，只是相應(yīng)的計算時間較長.因此，若缺陷位于肺實質(zhì)邊界較平滑的位置，使用凸包算法進(jìn)行修復(fù)即可得到良好的修復(fù)效果，并且凸包算法需要的人工干預(yù)較少，計算速度更快.若缺陷位于肺實質(zhì)邊界曲率變化較大的位置，采用本文的曲面重建修復(fù)方法可以得到更好的修復(fù)效果.臨床醫(yī)生對使用本文方法進(jìn)行肺實質(zhì)修復(fù)的結(jié)果進(jìn)行分析后認(rèn)為，該方法能夠有效地修復(fù)此類肺實質(zhì)缺陷，對臨床病理學(xué)的分析研究有一定的輔助作用.通過本文的修復(fù)方法得到的肺實質(zhì)可以完整地包含腫瘤，對腫瘤的分割也具有參考價值，后續(xù)將從肺腫瘤的三維重建方向進(jìn)行深入研究，以實現(xiàn)肺腫瘤精確和有效地分割.