劉明威， 顧力栩

(上海交通大學(xué) 生物醫(yī)學(xué)工程學(xué)院，上海 200240)

基于LUT的快速3D氣道樹骨架線提取

劉明威， 顧力栩

(上海交通大學(xué) 生物醫(yī)學(xué)工程學(xué)院，上海 200240)

為提高肺部支氣管骨架線的提取效率，提出并使用一種基于look-up-table(LUT) 的腐蝕細(xì)化算法. 分析建立腐蝕模型，并根據(jù)該模型優(yōu)化建立了LUT，以該LUT為依據(jù)通過索引查找對(duì)原始數(shù)據(jù)進(jìn)行快速腐蝕細(xì)化，對(duì)得到的腐蝕結(jié)果進(jìn)行剪枝處理以得到最終的骨架線. 實(shí)驗(yàn)結(jié)果表明：提取過程中LUT的應(yīng)用從根本上降低了腐蝕細(xì)化中判斷的復(fù)雜度，將復(fù)雜的簡(jiǎn)單點(diǎn)判斷問題轉(zhuǎn)化為L(zhǎng)UT中的查詢問題，從而極大地優(yōu)化了腐蝕細(xì)化中關(guān)鍵的腐蝕過程. 相比傳統(tǒng)方法，基于LUT的腐蝕細(xì)化算法顯著提高了骨架線的提取速度，較傳統(tǒng)細(xì)化法提速近22.95倍.

計(jì)算機(jī)輔助診斷； 氣管樹； 虛擬支氣管鏡；骨架線；查找表

肺癌是威脅人類生命健康的惡性腫瘤疾病，中國(guó)目前總體肺癌發(fā)病率增長(zhǎng)了19.48%，這很可能是加劇的環(huán)境污染問題如大規(guī)模霧霾引起的[1]. 早期診斷能大幅度提高病患存活率，但目前早期診斷率僅15%[2]. 常見肺癌早期診斷手段包括CT成像與光學(xué)內(nèi)窺鏡. CT成像擁有密度分辨率高、操作簡(jiǎn)單等優(yōu)勢(shì)，但傳統(tǒng)CT只能提供單層的斷面信息，不直觀，且在診斷過程中病人需承受大量放射輻射；光學(xué)內(nèi)窺鏡可在零輻射下觀察到管腔內(nèi)部真實(shí)紋理色彩，但操作上要比CT成像復(fù)雜得多.

隨著計(jì)算機(jī)輔助診斷技術(shù)的發(fā)展，虛擬支氣管鏡VB[3](virtual bronchoscopy)技術(shù)在肺癌診斷領(lǐng)域正逐漸得以普遍應(yīng)用. 該技術(shù)利用從CT數(shù)據(jù)中重建的氣管模型，對(duì)肺部疾病診斷有重大意義. 臨床表明VB不僅可以顯示有無氣道異常如支氣管狹窄等，還可以配合導(dǎo)航系統(tǒng)引導(dǎo)經(jīng)支氣管鏡肺活檢，明顯提高支氣管鏡在肺部外周病灶診斷效率. 相對(duì)于傳統(tǒng)內(nèi)窺鏡而言，VB引導(dǎo)下的導(dǎo)航技術(shù)具有以下優(yōu)點(diǎn)：1)定位準(zhǔn)確，診斷效能高；2)可快速引導(dǎo)支氣管鏡，縮短檢查時(shí)間，減少病患痛苦；3)可減少支氣管鏡檢查次數(shù)；4)可提供自由的虛擬觀測(cè)角度，降低支氣管異常的辨別難度； 5)可部分代替X線引導(dǎo)的支氣管鏡活檢，減少輻射；6)可避免盲目穿刺，降低并發(fā)癥發(fā)生率. 大量文獻(xiàn)證實(shí)虛擬內(nèi)窺鏡能在不失診斷準(zhǔn)確率的基礎(chǔ)上，保障診斷的微創(chuàng)性[4-5].

在虛擬內(nèi)窺鏡實(shí)現(xiàn)過程中，快速的骨架線提取一直是其中的難點(diǎn)和要點(diǎn). 而在VB中，一個(gè)良好提取的骨架線不僅可以提供直觀的支氣管拓?fù)浣Y(jié)構(gòu)，還可以為管腔內(nèi)瀏覽起到引導(dǎo)作用. 與人體其他管腔如結(jié)腸、血管相比，支氣管結(jié)構(gòu)更為復(fù)雜[6]，體現(xiàn)在：1)復(fù)雜的拓?fù)浣Y(jié)構(gòu)，擁有多級(jí)分支，多通路；2)各級(jí)分支粗細(xì)不一；3)個(gè)體差異性，不同人支氣管雖在級(jí)數(shù)上有一定相似性，但大小和分支位置卻存在差異. 目前三維骨架提取方式主要有距離變換法、中軸變換法和腐蝕細(xì)化算法. 然而由于支氣管的復(fù)雜性，這些方式在支氣管骨架線的提取上都略顯性能不足. 隨著當(dāng)今4DCT[7]數(shù)據(jù)的逐步普及，傳統(tǒng)方法的運(yùn)算速度更顯得不盡人意.

本文提出一種基于腐蝕細(xì)化的改進(jìn)三維骨架線提取方法，以克服傳統(tǒng)腐蝕細(xì)化方法中復(fù)雜度高、耗時(shí)長(zhǎng)的問題. 本方法充分利用了當(dāng)今計(jì)算機(jī)多核以及高內(nèi)存的優(yōu)勢(shì)，在傳統(tǒng)三維腐蝕過程中引入LUT. 相關(guān)實(shí)驗(yàn)證實(shí)，本方法能夠顯著提高骨架線的提取速度.

1 骨架線基本定義與現(xiàn)有方法綜述

1.1 骨架線基本定義

骨架(Skeleton)是指與原物體具有一致連通性和拓?fù)浣Y(jié)構(gòu)的細(xì)曲線的一種理想表示，就是位于物體內(nèi)部，且能體現(xiàn)其拓?fù)涮卣鞯暮?jiǎn)化圖形[8]. 在三維情況下提取的骨架線需要滿足如下要求[9]：1) 連通性，骨架線應(yīng)與目標(biāo)物體擁有相同的連通性；2) 居中性 ，骨架線盡量保持在管腔的中央；3) 單體素性，組成骨架線的體素要求為單元寬度. 對(duì)于簡(jiǎn)單的柱狀組織，中心線只有唯一的一條，而對(duì)于支氣管樹或血管樹這種復(fù)雜的管狀結(jié)構(gòu)，中心線會(huì)有多條，也被稱作為骨架線.

1.2 現(xiàn)有方法綜述

就目前而言，二維的中心線提取算法已經(jīng)相對(duì)成熟. 相對(duì)而言，三維中心線的提取方式只有距離變換法、中軸變換法和腐蝕細(xì)化算法.

距離變換法[10]是指根據(jù)不同目的，將物體的體素點(diǎn)進(jìn)行距離標(biāo)識(shí). 目前主流的方法是基于雙距離場(chǎng)的方法[11]：首先對(duì)物體進(jìn)行邊界距離變換，得到體素到達(dá)邊界的邊界距離場(chǎng)DFB (distance from boundary)；其次指定源點(diǎn)，進(jìn)行源點(diǎn)距離變換，得到源距離場(chǎng)DFS (distance from source)；最后通過邊界距離場(chǎng)中局部最大值的篩選，配合源距離場(chǎng)進(jìn)行骨架點(diǎn)的連接. 中軸變換法[12]在二維上可以用最大內(nèi)切圓盤的相關(guān)術(shù)語定義，而在三維上則表現(xiàn)為最大內(nèi)切球中心的集合. 其中圓盤中心到達(dá)邊界的距離稱為圓盤法線，一個(gè)目標(biāo)內(nèi)點(diǎn)的最大內(nèi)切圓在邊界上至少有兩個(gè)切點(diǎn)，而每個(gè)對(duì)稱點(diǎn)有兩個(gè)或多個(gè)圓盤法線. 腐蝕細(xì)化算法[13]的思想則更為樸素，根據(jù)物體歐拉特性[14]以及連通性的不變性，給出在刪除后不影響原物體拓?fù)浣Y(jié)構(gòu)的簡(jiǎn)單點(diǎn)，并以此為依據(jù)均勻、對(duì)稱地從不同方向?qū)δ繕?biāo)物體邊緣上的簡(jiǎn)單點(diǎn)(Simple Voxel, SV)[15]進(jìn)行剝離，直至僅剩下單像素寬的體素集合，去除該集合中的偽骨架點(diǎn)后得到最終的骨架線.

以上方法各自均有不同的局限性和一些共同不足. 距離變換法中，在構(gòu)建雙距離場(chǎng)時(shí)不僅因需要對(duì)各體素點(diǎn)以迭代的方式進(jìn)行距離標(biāo)記而導(dǎo)致高時(shí)耗，還因需要與目標(biāo)物體同樣大小的內(nèi)存空間儲(chǔ)存距離場(chǎng)而帶來高內(nèi)存消耗. 同時(shí)，依據(jù)不同的歐氏距離定義建立的距離場(chǎng)會(huì)產(chǎn)生不同的候選中心點(diǎn)，直接影響骨架線結(jié)果的一致性[16]，且候選中心點(diǎn)通常是不連續(xù)的，連接時(shí)會(huì)造成骨架線在曲率較大位置的離心偏移，難以保證居中性. 中軸變換法中，由于中軸變換是理想的數(shù)學(xué)模型描述，主要基于數(shù)學(xué)的幾何推理證明，從而導(dǎo)致使用機(jī)器語言進(jìn)行實(shí)現(xiàn)的難度大、算法復(fù)雜度高，常見的模型有地表火模型[17]，距離曲面脊線模型[18]. 以上兩種方法存在共同的缺陷就是易受邊緣噪聲干擾，且難以保證骨架的連續(xù)性. 腐蝕細(xì)化算法的不足在于需要不斷地對(duì)邊緣的簡(jiǎn)單點(diǎn)進(jìn)行判斷，現(xiàn)有的判斷模型復(fù)雜度均不夠理想，導(dǎo)致當(dāng)目標(biāo)物體體素點(diǎn)過多時(shí)，帶來巨大的時(shí)耗.

2 改進(jìn)的腐蝕細(xì)化方法

考慮到基于距離場(chǎng)的算法對(duì)于邊緣噪聲的敏感性和內(nèi)存的大量消耗，以及傳統(tǒng)腐蝕細(xì)化算法中對(duì)邊緣點(diǎn)反復(fù)判斷的高復(fù)雜度，本文提出一種基于LUT的快速三維腐蝕細(xì)化算法. 該方法保留了傳統(tǒng)腐蝕細(xì)化算法的穩(wěn)定性，以及結(jié)果中對(duì)目標(biāo)物體拓?fù)浣Y(jié)構(gòu)和連通性的完好保留性，并通過使用LUT的思想簡(jiǎn)化細(xì)化過程，避免重復(fù)判斷，顯著提高了細(xì)化速度. 該方法首先從原始數(shù)據(jù)中用區(qū)域增長(zhǎng)的方式分割出支氣管二值圖像，建立腐蝕細(xì)化模型并建立LUT，然后查找LUT進(jìn)行腐蝕細(xì)化，最后對(duì)細(xì)化的結(jié)果進(jìn)行剪枝處理.

2.1 支氣管二值數(shù)據(jù)的獲取

采用有監(jiān)控的動(dòng)態(tài)閾值區(qū)域生長(zhǎng)配合形態(tài)學(xué)閉操作來獲取組成支氣管，流程及效果如圖1所示. 首先以氣管壁與空氣CT值的差異為依據(jù)，通過應(yīng)用DMQ (dynamic marking QUEUE)的數(shù)據(jù)結(jié)構(gòu)進(jìn)行增長(zhǎng)，提高了提取的速度. 在頂層CT中選擇氣管入口作為增長(zhǎng)原點(diǎn)壓入DMQ，目標(biāo)點(diǎn)離開DMQ時(shí)對(duì)其6鄰域進(jìn)行增長(zhǎng)判斷，若符合則壓入DMQ. 在判斷目標(biāo)點(diǎn)是否歸入增長(zhǎng)區(qū)域時(shí)，采用了18鄰域平均值比較的方式，最大可能增加分割目標(biāo)點(diǎn)數(shù). 對(duì)每一輪新增的點(diǎn)與上一輪增長(zhǎng)結(jié)果進(jìn)行比較，以進(jìn)行溢出監(jiān)控. 該方法得到一個(gè)26連通的初步結(jié)果，為平滑邊緣、填充管腔的空隙，在其基礎(chǔ)上使用公式IMG ·M = (IMG ? M) ⊕M進(jìn)行形態(tài)學(xué)閉操作，用同一模版M分別對(duì)IMG進(jìn)行膨脹腐蝕操作得到最終的二值分割結(jié)果.

2.2 細(xì)化方法

細(xì)化過程具體分為4步：1)腐蝕模型的建立； 2)基于腐蝕模型建立LUT(提高速度的關(guān)鍵步驟)；3)基于LUT進(jìn)行快速細(xì)化； 4)對(duì)腐蝕結(jié)果進(jìn)行剪枝處理.

圖1 支氣管分割流程

2.2.1 腐蝕模型的建立

在三維空間V3中，設(shè)中心點(diǎn)為c(x,y,z),對(duì)于任意點(diǎn)p(i,j,k)有如下的鄰域定義：

Nb6(p) = {p(i,j,k)∈V3| |i-x| + |j-y|+ |k-z| = 1},

Nb26(p) = {p(i,j,k)∈V3| |i-x|<2 ∩|j-y|<2∩|k-z| <2},

Nb18(p) = {p(i,j,k)∈V3| |i-x| + |j-y|+ |k-z| < 3}∩ Nb26(p).

將屬于目標(biāo)物體的體素定義為前景點(diǎn)Sfg，其余定義為背景點(diǎn)Sbg. 在前景點(diǎn)集合中存在子集Ssub-fg，若其中所有元素均符合以下四種判斷準(zhǔn)則[19]，則稱之為簡(jiǎn)單點(diǎn)SV(simplevoxel).

a) ?Sfgin Nb26and Sum(Sfgin Nb26)>1；

b) ?Sbgin Nb6；

c) ?Sfgin Nb26and is 26Connected(Sfg)；

d) ?Sbgin Nb6and is 6Connected By Nb18(Sbg).

a、b條件可以直接通過相應(yīng)鄰域中的前景點(diǎn)、背景點(diǎn)分布獲得，而c、d條件的判斷才是關(guān)鍵.

對(duì)于條件c，提出一種26鄰域區(qū)域增長(zhǎng)法，如圖2所示. 根據(jù)條件c的描述易證明：若Nb26中的前景點(diǎn)可以彼此26連通，那么在該Nb26中，一定有且只有一個(gè)26連通的前景點(diǎn)區(qū)域. 因此，取26鄰域中的任意Sfg，在該Nb26中進(jìn)行26鄰域的前景點(diǎn)區(qū)域增長(zhǎng)，得到增長(zhǎng)結(jié)果集合S′. 將其與26鄰域中的前景點(diǎn)構(gòu)成的集合S進(jìn)行比較，若兩者點(diǎn)數(shù)相同，則符合條件c.

在條件d的判斷上，并不能完全與條件c等同，因?yàn)樵贜b18中的Sbg中，6連通的數(shù)量并不唯一. 充分考慮到這一點(diǎn)，利用Nb18中Sbg的6鄰域局部區(qū)域增長(zhǎng)法，使用列表套列表的結(jié)構(gòu)來儲(chǔ)存連通的6鄰域數(shù)組. 在遍歷結(jié)束后，將Nb6中的背景點(diǎn)集合S與連通的6鄰域列表中的Si依次進(jìn)行比較，得到是否可以刪除.

圖2 條件c判斷

2.2.2 鄰域分布分析與LUT的建立

分析腐蝕模型可以看出，判斷一個(gè)前景點(diǎn)是否為簡(jiǎn)單點(diǎn)是由a、b、c、d這4個(gè)條件的滿足情況決定的. 進(jìn)一步分析，條件a為非孤立點(diǎn)的判斷，取決于Nb26；條件b為邊界點(diǎn)的判斷，取決于Nb6；條件c為Nb26的背景點(diǎn)的26連通情況，取決于Nb26；條件d為Nb26前景點(diǎn)的18連通情況，同樣取決于Nb26. 因此可得出結(jié)論，對(duì)于任意體素點(diǎn)p，可根據(jù)其Nb26(p)分布情況VC(voxelcombination)得到其簡(jiǎn)單性.

在二值圖像前提下，任意體素點(diǎn)只有0或1兩種取值. 如圖3所示，若對(duì)26鄰域指定排列順序，則可等同于一個(gè)26bits的二進(jìn)制數(shù)，進(jìn)而任意一種VC將對(duì)應(yīng)唯一的26bits二進(jìn)制數(shù)，而相對(duì)的任意一個(gè)26bits二進(jìn)制數(shù)也能唯一映射為一種VC. 基于以上分析，可聯(lián)想到建立LUT輔助腐蝕細(xì)化. 在腐蝕細(xì)化過程中，將目標(biāo)點(diǎn)的26鄰域VC轉(zhuǎn)化為26bits的索引，在LUT中進(jìn)行查詢判斷. 文獻(xiàn)[20]顯示LUT在二維的細(xì)化中得到了很好的應(yīng)用，值得一提的是，由于三維計(jì)算量巨大，空間分布復(fù)雜性以及受拘束于過去計(jì)算機(jī)的內(nèi)存，三維LUT鮮有人嘗試.

為避免建立LUT過程中的重復(fù)判斷，對(duì)體素的鄰域的空間分布狀態(tài)進(jìn)行特征分析. 對(duì)于任意立方體，給出如下對(duì)于相似體素分布SVC(similarvoxelcombination)的定義：

1) 通過任意旋轉(zhuǎn)變換可以完全重合；

2) 通過左右手坐標(biāo)系對(duì)換可以完全重合，即

圖3 26鄰域轉(zhuǎn)換

(1)

(2)

(3)

條件1)可解釋為依次以立方體每一面朝向的坐標(biāo)系方向?yàn)橛^察方向，而后以該方向?yàn)檩S向，使用如式(1)～(3)中對(duì)應(yīng)方向的旋轉(zhuǎn)矩陣對(duì)立方體進(jìn)行圍繞該軸向的三維旋轉(zhuǎn)，為保證立方體空間位置不變，這里旋轉(zhuǎn)的角度應(yīng)該分別取0°、90°、180°、270°，這樣通過旋轉(zhuǎn)可得4*6 共24種SVC；條件2)可解釋為條件1)基礎(chǔ)上的鏡像變換，因此同樣存在24種SVC. 因此可得出結(jié)論，對(duì)于各向異性的立方體，存在24+24共48種SVC，應(yīng)用在二值化的Nb26情況下，粗略估算在全空間2^26中存在2^26/48約為1398101種SVC.

(4)

而在二值化的Nb26情況下，并非所有分布都是各向異性，存在著大量的自反對(duì)稱，經(jīng)過式(4)的計(jì)算，其中i表示Nb26中前景點(diǎn)的個(gè)數(shù)，R(i) 表示Nb26中重復(fù)的組合個(gè)數(shù)，由于二值互補(bǔ)性，因此只需計(jì)算13次(0～12個(gè)前景點(diǎn)的情況)，外加13個(gè)前景點(diǎn)的情況. 針對(duì)i取值不同時(shí)需要具體列出對(duì)稱情況進(jìn)行計(jì)算，但由于其復(fù)雜的旋轉(zhuǎn)對(duì)稱性，需要使用計(jì)算機(jī)輔助邏輯進(jìn)行計(jì)算，最終得到不同前景點(diǎn)個(gè)數(shù)下SVC如表1所示，SVC一共有1426 144種，并在實(shí)際運(yùn)算遍歷中加以檢驗(yàn).

表1 26鄰域中前景體素點(diǎn)與SVC關(guān)系表

基于以上分析，本文提出多線程并行配合雙布爾數(shù)組的算法建立LUT，數(shù)組大小為2^26，下標(biāo)代表該索引值所唯一對(duì)應(yīng)的VC. 其中一個(gè)數(shù)組為L(zhǎng)UT，是后續(xù)腐蝕細(xì)化中所依賴的參照，其取值表示該VC下的中心點(diǎn)是否為簡(jiǎn)單點(diǎn)；另一個(gè)數(shù)組為MT(marking table)，其取值用來表示該VC是否已被計(jì)算. 將所有2^26個(gè)索引等分為多個(gè)索引段，平均分配給多個(gè)子線程，配合多核CPU進(jìn)行計(jì)算，以達(dá)到最高的計(jì)算效率.

所有子線程將共享LUT和MT，每個(gè)子線程將同時(shí)并行遍歷被分配的索引段. 在具體計(jì)算索引時(shí)，每一次簡(jiǎn)單點(diǎn)判斷的成本很高，但應(yīng)用雙數(shù)組算法可有效避免重復(fù)計(jì)算. 對(duì)于一個(gè)索引值，首先在MT中查找判斷該索引值是否已被計(jì)算，若尚未計(jì)算，則將該索引值轉(zhuǎn)換為其對(duì)應(yīng)的VC，并根據(jù)前文的判斷條件和方法得出其是否為簡(jiǎn)單點(diǎn)；而后通過正方體的旋轉(zhuǎn)對(duì)稱性，獲取該VC的所有其他47個(gè)SVC，并轉(zhuǎn)換為索引序列IA(index array)；最后在LUT中將IA中的索引值賦予與該VC相同的簡(jiǎn)單性，同時(shí)在MT中將IA中的這些索引值標(biāo)記為已判斷. 雙數(shù)組的優(yōu)勢(shì)也在此，通過標(biāo)記MT避免了遍歷到SVC所對(duì)應(yīng)索引值時(shí)潛在的重復(fù)簡(jiǎn)單性判斷，從而大大減少了判斷的次數(shù)，整體加速了LUT的建立過程. 為更為直觀，下面給出該流程框架的偽代碼.

函數(shù) 1: CreateLUT(創(chuàng)建LUT 的主函數(shù))

輸入: 需要的線程數(shù)目 threadNumber.

輸出: 計(jì)算好的lutArray.

for eachiinthreadNumberdo

subThreadsi=start_new_thread(MultiThreadCalculate(lutArray,mtArray,rangei) )

endfor

wait_sub_threads(subThreads)

write_to_bin_file(lutArray)

函數(shù)2:MultiThreadCalCulate(用于計(jì)算指定區(qū)段的子函數(shù))

輸入: 公共使用的lutArray，用于標(biāo)記是否計(jì)算過的mtArray，需要遍歷的區(qū)段range.

輸出: 計(jì)算指定區(qū)段range后的lutArray， 標(biāo)記指定區(qū)段range后的mtArray

foreachiinRange

ifnotmtArrayithen

isSimple=isSimplePoint(i)

similarIndexArray=getSimilarIndex(i)

foreachjinsimilarIndexArray

lutArrayj=isSimple

mtArrayj=true

endfor

endif

endfor

2.2.3 基于LUT進(jìn)行迭代腐蝕細(xì)化

在建立好LUT的基礎(chǔ)上，傳統(tǒng)的腐蝕細(xì)化過程便得以優(yōu)化，其中對(duì)于SV的判斷可以直接依據(jù)LUT的結(jié)果，而不是重新對(duì)前文4個(gè)定義條件的判斷. 在進(jìn)行細(xì)化前先將圖像中的前景點(diǎn)存入臨時(shí)數(shù)組PointListfg，以避免重復(fù)遍歷，提高效率. 而后進(jìn)行如下操作：

1)在PointListfg中根據(jù)Nb26體素分布情況尋找邊緣點(diǎn)，得到PointListborder.

2)遍歷PointListborder中的點(diǎn)，將Nb26體素分布轉(zhuǎn)換為26bits索引值，根據(jù)LUT的判斷值進(jìn)行簡(jiǎn)單點(diǎn)的判斷刪除.

3)重復(fù)步驟1)、2)直至PointListborder中再也無法找到簡(jiǎn)單點(diǎn)為止.

2.2.4 后續(xù)剪枝處理

實(shí)驗(yàn)證實(shí)，通過上述腐蝕細(xì)化算法所得結(jié)果會(huì)在一些較粗的主氣管上出現(xiàn)一些額外的毛刺，這無疑會(huì)對(duì)后續(xù)的內(nèi)窺鏡路徑的規(guī)劃帶來不必要的分枝，需要進(jìn)行剪枝處理.

根據(jù)毛刺分支長(zhǎng)度短的特點(diǎn)，設(shè)計(jì)了一種尺寸自適應(yīng)的剪枝算法. 記當(dāng)前的前景點(diǎn)集合為S，首先根據(jù)骨架線的體素?cái)?shù)設(shè)置適當(dāng)?shù)募糁﹂撝礣hreshold，設(shè)閾值百分比為th%，這里th取經(jīng)驗(yàn)值1.7，因此經(jīng)計(jì)算得到 Threshold=S×th%. 接著遍歷細(xì)化后的骨架線，得到所有端點(diǎn)PointListend，以PointListend中的點(diǎn)為起點(diǎn)進(jìn)行遍歷，直至分支點(diǎn)，若其長(zhǎng)度小于Threshold則將該支刪除.

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

3.1 實(shí)驗(yàn)環(huán)境與實(shí)驗(yàn)數(shù)據(jù)

實(shí)驗(yàn)數(shù)據(jù)采用從上海肺科醫(yī)院獲得的5組肺CT圖像序列，每組序列均包含有最少265層到最多468層CT斷層圖像，每層圖像的分辨率均是512*512 且灰度值在-2048到4196區(qū)間內(nèi). 用于實(shí)驗(yàn)的PCCPU為Intel(R)Core(TM)i7-3610QM，RAM為6GB，OS為Win7 64-bit. 算法的實(shí)現(xiàn)基于VTKv6.1.0開源庫(kù)，配合QT在VS2008IDE下進(jìn)行開發(fā)和實(shí)驗(yàn). 3.2 LUT的建立

為獲取最適合本實(shí)驗(yàn)環(huán)境的線程數(shù)量，以得到最快的LUT建立速度，以線程數(shù)作為輸入?yún)?shù)進(jìn)行實(shí)驗(yàn)，得到如圖4所示的結(jié)果.

圖4 LUT建立線程時(shí)間函數(shù)

可見當(dāng)線程數(shù)少于5個(gè)時(shí)，增加線程可顯著減少用時(shí)，5個(gè)線程的處理過程表現(xiàn)最為優(yōu)秀，當(dāng)線程多于5個(gè)時(shí)，用時(shí)無明顯減少.

3.3 骨架線提取結(jié)果與比較

為比較基于LUT的方法與傳統(tǒng)腐蝕細(xì)化算法的效率，采用如下的對(duì)比試驗(yàn)，使用相同的數(shù)據(jù)集進(jìn)行重復(fù)試驗(yàn)100次，表2列出了實(shí)驗(yàn)結(jié)果.

在腐蝕細(xì)化階段，基于LUT的平均時(shí)間為482.74 ms，而基于傳統(tǒng)腐蝕細(xì)化算法的平均時(shí)間為11 077.7 ms，可見本文方法提高了22.95倍速度.

3.4 骨架線的提取效果

圖5為本文算法提取骨架線的結(jié)果與基于距離變換算法[21]的結(jié)果，其中圖5(a)、(b)、(c)三組結(jié)果中，上方為本文算法的結(jié)果，下方為距離變換算法的結(jié)果. 通過對(duì)比可以看出，基于距離變換的方法會(huì)在末端發(fā)生斷裂，不連續(xù)，在同一支氣管段上可能產(chǎn)生多分支，不方便路徑規(guī)劃，且操作中需要手工選定端點(diǎn)，既不方便也無法保證連續(xù)性；而本文方法在支氣管骨架線的提取上既可以保證骨架線在細(xì)末端的連續(xù)性，不發(fā)生斷裂，也無需手動(dòng)，具備可重復(fù)性. 圖5中的三組數(shù)據(jù)顯示，該算法很好地保持了支氣管拓?fù)浣Y(jié)果，同時(shí)也保障骨架線的中心性、連續(xù)性、一致性和可重復(fù)性，為內(nèi)窺鏡的漫游提供可能.

表2 實(shí)驗(yàn)結(jié)果對(duì)比表

(a) 512*512*265 (b) 512*512*429 (c) 512*512*438

4 實(shí)驗(yàn)討論

4.1 基于LUT算法的優(yōu)越性

本文同時(shí)還與基于距離場(chǎng)[21]算法得出的速度結(jié)果進(jìn)行比較，由于使用不同的數(shù)據(jù)集，因此定義指標(biāo)ts(thining speed)來表示算法的效率作為比較依據(jù)，ts表示毫秒內(nèi)處理的立方點(diǎn)數(shù), 單位為103voxels·ms-1. 經(jīng)計(jì)算，文獻(xiàn)[21]的速度為220.16 ts，而本文的速度為45 253.39 ts，提高了205倍. 結(jié)合本文與以上兩種優(yōu)化過的傳統(tǒng)骨架線提取方式[21-22]的比較，可以充分看出基于LUT的算法能在保證骨架線拓?fù)浣Y(jié)構(gòu)完整的基礎(chǔ)上，顯著提高運(yùn)行速度.

總的來說，基于LUT的算法在本質(zhì)上優(yōu)化了傳統(tǒng)細(xì)化過程中的腐蝕細(xì)化環(huán)節(jié)，從而在整體上提高了骨架線提取的速度. 不僅如此，LUT的另一大優(yōu)勢(shì)是將LUT的建立與LUT的參照分離，這樣在簡(jiǎn)單點(diǎn)的定義發(fā)生變化時(shí)，可以直接更新LUT的數(shù)據(jù)，而不需要更改后續(xù)整個(gè)程序.

4.2 參數(shù)與復(fù)雜度分析

基于LUT的腐蝕細(xì)化算法中，主要的輸入?yún)?shù)體現(xiàn)于兩點(diǎn)：一是LUT建立過程中線程數(shù)的選擇，二是應(yīng)用LUT進(jìn)行腐蝕細(xì)化時(shí)，不同輸入物體的體素?cái)?shù)量.

在LUT建立過程中，線程數(shù)是重要的輸入?yún)?shù)，并影響最終LUT建成的時(shí)耗. 實(shí)驗(yàn)結(jié)果表明5個(gè)線程為最優(yōu)，分析其原因?yàn)椋寒?dāng)線程池大小為CPU數(shù)+1時(shí)，可達(dá)到最優(yōu)利用率，即當(dāng)執(zhí)行密集型任務(wù)時(shí)既使偶爾因頁缺失故障或者其他原因?qū)е聲和r(shí)，此額外線程也可確保CPU的時(shí)鐘周期不被浪費(fèi). 實(shí)驗(yàn)PC 為四核CPU，因此5個(gè)線程可達(dá)到最優(yōu)化. 算法的空間復(fù)雜度為2*2^26*1 byte = 128 MB，用于緩存雙數(shù)組中的LUT和MT，在LUT建成結(jié)束后，得到的結(jié)果為一個(gè)64MB的二進(jìn)制文件，可以直接作為查找腐蝕的依據(jù). 在時(shí)間上，對(duì)腐蝕細(xì)化的判斷條件c的判斷由于應(yīng)用了26鄰域內(nèi)區(qū)域增長(zhǎng)法，只用一重循環(huán)，是一個(gè)線性過程，將該步驟的時(shí)間復(fù)雜度降低為O(n)；而對(duì)于條件d的判斷，在區(qū)域增長(zhǎng)上的平均時(shí)間復(fù)雜度依舊為O(n)，而在列表比較上的平均時(shí)間復(fù)雜度為O(n2)，兩步之間為串行，因而整體的復(fù)雜度依舊為O(n2)，但較傳統(tǒng)方法相比依舊有所提升.

在腐蝕細(xì)化的過程中，由于LUT已經(jīng)確立，因此主要的輸入?yún)?shù)為物體的體素?cái)?shù). 對(duì)3.3中的實(shí)驗(yàn)結(jié)果進(jìn)行線性擬合，得到t=0.007 2v-152.91，擬合優(yōu)度R2為0.988 5，其中t表示腐蝕細(xì)化的時(shí)耗，單位為ms，v表示物體的體素?cái)?shù). 因此可得，腐蝕細(xì)化的時(shí)耗正比于輸入物體的體素?cái)?shù). 算法的空間復(fù)雜度取決于物體的體素?cái)?shù)，會(huì)隨著CT采集時(shí)分辨率的不同而有所變化，另外需考慮LUT占用的64MB；在時(shí)間復(fù)雜度上，查詢LUT接近于O(1)，對(duì)體素點(diǎn)進(jìn)行擦除的復(fù)雜度也為O(1)，因此整體復(fù)雜度仍為O(1).

4.3 基于LUT算法的不足

與基于其他理論的骨架線提取方式如距離變換法相比，基于LUT的腐蝕細(xì)化法不足在于最終得到的骨架線是組成骨架線的離散點(diǎn)，缺少點(diǎn)與點(diǎn)之間拓?fù)潢P(guān)系的描述，為后續(xù)骨架線拓?fù)浣Y(jié)構(gòu)分析帶來了一定的困難. 同時(shí)剪枝過程也有一定的局限性，無法精確區(qū)分細(xì)末較短的支氣管與偽分支. 而在LUT的使用中依舊可以作適當(dāng)優(yōu)化，如在查詢表的尺寸上，若是對(duì)內(nèi)存有苛刻的要求，可將LUT壓縮為8MB，以bit的形式來儲(chǔ)存每一個(gè)索引的索引值. 由于對(duì)bit的檢索需要額外的比較，其效率可能會(huì)有一定下降，這是一個(gè)棘手的問題.

受實(shí)驗(yàn)數(shù)據(jù)源限制，未進(jìn)行大規(guī)模驗(yàn)證. 同時(shí)受實(shí)驗(yàn)室硬件條件限制，缺少在其他配置相當(dāng)?shù)挠?jì)算機(jī)上進(jìn)行重復(fù)驗(yàn)證，也沒有獲取該算法在更差運(yùn)行環(huán)境中的表現(xiàn). 更加充足的數(shù)據(jù)源和豐富的實(shí)驗(yàn)平臺(tái)將能夠?qū)υ撍惴ㄓ幸粋€(gè)更全面的評(píng)估，這也是未來工作的目標(biāo).

5 結(jié) 語

基于LUT的腐蝕細(xì)化算法在傳統(tǒng)的三維腐蝕細(xì)化算法中引入LUT的思想，通過鄰域區(qū)域增長(zhǎng)法優(yōu)化了簡(jiǎn)單點(diǎn)腐蝕模型的判斷，并通過二值化26鄰域旋轉(zhuǎn)對(duì)稱性的分析進(jìn)行快速的LUT建立，在腐蝕細(xì)化中使用LUT取代傳統(tǒng)算法中對(duì)簡(jiǎn)單點(diǎn)的復(fù)雜重復(fù)判斷過程，從而顯著提高了骨架線的提取速度. 實(shí)驗(yàn)結(jié)果表明：本文方法在速度上是傳統(tǒng)細(xì)化法方法[22]的22.95倍，傳統(tǒng)距離場(chǎng)方法[21]的近220倍，可以完好地保留肺部氣管樹的拓?fù)浣Y(jié)構(gòu)，具備易操作性、連續(xù)性和可重復(fù)性，為內(nèi)窺鏡瀏覽導(dǎo)航提供可能.

[1] CHEN W, ZHENG R, ZENG H, et al. Annual report on status of cancer in China, 2011[J]. Chinese Journal of Cancer Research, 2015, 27(1): 2.

[2] 韋春暉. 肺癌早期診斷進(jìn)展[J]. 臨床肺科雜志, 2010, 15(8):1136-1138.

WEI Chunhui. Progress in early diagnosis of lung cancer[J]. Journal of Clinical Pulmonary Medicine,2010,15(8):1136-1138.

[3] JUNG S Y, PAE S Y, CHUNG S M, et al. Three-dimensional CT with virtual bronchoscopy: a useful modality for bronchial foreign bodies in pediatric patients[J]. European Archives of Oto-Rhino-Laryngology, 2012, 269(1): 223-228.

[4] WERNER H, SANTOS J R L D, FONTES R, et al. Virtual bronchoscopy for evaluating cervical tumors of the fetus[J]. Ultrasound in Obstetrics & Gynecology, 2013, 41(1): 90-94.

[5] DE WEVER W, VANDECAVEYE V, LANCIOTTI S, et al. Multidetector CT-generated virtual bronchoscopy: an illustrated review of the potential clinical indications[J]. European Respiratory Journal, 2004, 23(5): 776-782.

[6] DENIZ A, HOFFMAN E A, GEOFFREY M L, et al. Segmentation and analysis of the human airway tree from three-dimensional X-ray CT images.[J]. IEEE Transactions on Medical Imaging, 2003, 22(8):940-50.

[7] COLE A J, O'Hare J M, McMahon S J, et al. Investigating the potential impact of four-dimensional computed tomography (4DCT) on toxicity, outcomes and dose escalation for radical lung cancer radiotherapy[J]. Clinical Oncology, 2014, 26(3): 142-150.

[8] 陳剛, 呂煊, 王志成, 等. 肺 CT 圖像的血管骨架化方法[J]. 計(jì)算機(jī)科學(xué), 2013, 40(5): 274-278.

CHEN Gang, LV Xuan, WANG Zhicheng, et al. Vessel Skeletonization Method for Lung CT Images[J].Computer Science, 2013, 40(5): 274-278.

[9] ZHOU Y, TOGA A W. Efficient skeletonization of volumetric objects[J]. IEEE Transactions on Visualization & Computer Graphics, 1999, 5(3):196-209.

[10]WAN M, DACHILLE F, KAUFMAN A. Distance-field based skeletons for virtual navigation[C]//Proceedings of the Conference on Visualization'01. IEEE Computer Society, 2001: 239-246.

[11]WANG S, WU J, WEI M, et al. Robust curve skeleton extraction for vascular structures[J]. Graphical Models, 2012, 74(4): 109-120.

[12]潘鵬, 賀三維, 吳艷蘭, 等. 曲邊多邊形中軸提取的新方法[J]. 測(cè)繪學(xué)報(bào), 2012, 41(2): 278-283.

PAN Peng, HE Sanwei, WU Yanlan, et al. A new method for extracting curved-polygon medial axis[J]. Acta Geodaetica et Cartographica Sinica, 2012, 41(2): 278-283.

[13]CHEN W, SUI L, XU Z, et al. Improved Zhang-Suen thinning algorithm in binary line drawing applications[C]// International Conference on Systems and Informatic. Yantai:IEEE, 2012: 1947-1950.

[14]CHOI W P, LAM K M, SIU W C. Extraction of the Euclidean skeleton based on a connectivity criterion[J]. Pattern Recognition, 2003, 36(3): 721-729.

[15]MICHEL C, GILLES B. New Characterizations of Simple Points in 2D, 3D, and 4D Discrete Spaces[J]. IEEE Transactions on Pattern Analysis & Machine Intelligence, 2008, 31(4):637-648.

[16]GENG H, YANG J, TAN W, et al. Fast 3D skeleton extraction of airways and applications to virtual bronchoscopy[C]// The 26th Chinese Control and Decision Conference (2014 CCDC). Changsha:IEEE, 2014: 3879-3884.

[17]LIU L, CHAMBERS E W, LETSCHER D, et al. Extended grassfire transform on medial axes of 2D shapes[J]. Computer-Aided Design, 2011, 43(11): 1496-1505.

[18]GOLDBERGERA L, AMARAL L A N, GLASS L, et al. Physiobank, physiotoolkit, and physionet components of a new research resource for complex physiologic signals[J]. Circulation, 2000, 101(23): e215-e220.

[19]PALAGYI K. A 3D fully parallel surface-thinning algorithm[J]. Theoretical Computer Science, 2008, 406(1-2):119-135.

[20]楊威, 郭科, 魏義坤. 一種有效的基于八鄰域查表的指紋圖像細(xì)化算法[J]. 四川理工學(xué)院學(xué)報(bào)(自然科學(xué)版), 2008, 21(2): 61-63.

YANG Wei, GUO Ke, WEI Yikun. An effective thinning algorithm for finger print based on 8-neighbors LUT[J]. Journal of Sichuan University of Science and Engineering(Natural Science Edition), 2008, 21(2): 61-63.

[21]JIANG G, GU L. An automatic and fast centerline extraction algorithm for virtual colonoscopy[C]// 27th Annual International Conference of the Engineering in Medicine and Biology Society. Shanghai:IEEE, 2005: 5149-5152.

[22]陳磊, 王勝軍, 鄭全錄, 等. 基于 CT 圖像的三維拓?fù)浼?xì)化算法及其在心臟 CAD 中的應(yīng)用[J]. 計(jì)算機(jī)應(yīng)用, 2007, 27(B06): 406-410.

CHEN Lei, WANG Shengjun, ZHENG Quanlu, et al. A 3D topological thinning algorithm based on CT image and its application in cardiac CAD, Computer Application, 2007, 27(B06): 406-410.

(編輯 王小唯， 苗秀芝)

A fast LUT-based airway skeleton extraction algorithm for virtual bronchoscopy

LIU Mingwei, GU Lixu

(School of Biomedical Engineering, Shanghai Jiao Tong University, Shanghai 200240, China)

In order to improve the efficiency of lung airway skeleton extraction, this paper introduces a new look-up-table (LUT) based thinning algorithm. This new approach consists of three major steps: the analysis and creation of the thinning model and the establishment of LUT based on above thinning model result. Thinning process by index-searching uses the LUT. Branch cutting based on the result of thinning obtains the final result. The usage of LUT index-searching during the thinning step transforms simple point judgments into LUT index search and has significantly improved the performance of the whole algorithm. Experimental results demonstrate that the new skeleton algorithm is 22.95 times faster than the existing thinning algorithm.

computer-aided diagnosis; lung airway; virtual bronchoscopy; skeleton; LUT

10.11918/j.issn.0367-6234.201601026

2016-01-01

國(guó)家自然科學(xué)基金(61271318)

劉明威(1990—)，男，碩士; 顧力栩(1964—)，男，教授，博士生導(dǎo)師

顧力翔，gulixu@sjtu.edu.cn.

文獻(xiàn)標(biāo)志碼: A 文章編號(hào): 0367-6234(2017)05-0134-07