基于三維曲面重建的CT圖像肺邊界缺陷修復(fù)

張欣，徐永瀟，王兵

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

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

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

1 算法描述與分析

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

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

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

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

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

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

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

(1)

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

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

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

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

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

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

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

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

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

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

在肺實(shí)質(zhì)缺陷處周?chē)谋砻孢M(jìn)行采樣，將所得點(diǎn)集記為P.具體步驟如下.

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

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

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

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

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

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

(2)

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

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

(3)

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

(4)

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

(5)

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

2 實(shí)驗(yàn)分析

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

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

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

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

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

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

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

.

(6)

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

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

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

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