無二義性從輪廓線重建三維醫學表面的新方法
2009-01-01 00:00:00何金國
計算機應用研究 2009年5期
（1. 中央民族大學 理學院 北京100081； 2. 北京大學 信息科學中心 北京100871)
摘 要：設計了一種將輪廓線三維建模與多尺度二維圖像處理技術相結合的新方法。該方法的核心是邊緣距離圖算法。該方法對于復雜環境下的重建結果優于傳統方法，對于傳統方法較難解決的復雜環境下輪廓線對應和分支問題，能較好地解決。并且，當原始切片縱向分辨率較低時，該方法有較高的抵抗插值噪聲的能力。
關鍵詞：三維表面; 輪廓線; 圖像切片; 分支與對應問題
中圖分類號：TP391.9;TP302文獻標志碼：A
文章編號：1001-3695(2009)05-1997-04
New method for uniquely reconstructing 3D medical surfaces from contours
HE Jin-guo1，2
（1. College of Science Central University for Nationalities Beijing100081 China; 2. Center for Information Science Peking University Beijing 100871 China)
Abstract:This paper proposed a new contour reconstruction method combining multi-scale strategy and the core component of new method is edge distance graph algorithm.New methods could get better results under complicated environments since correspondence and branching problems were solved better. Also new methods were robust while running upon noise slices.
Key words：3D surfaces; contours; slices; correspondence and branching problem
0 引言
從一個序列的平行輪廓線出發，通過三角面片將它們連接而成為一個三維網格表面，就是輪廓線重建問題，有的文獻也稱之為輪廓線連接(contours connecting)問題。基于醫學圖像切片重建三維表面，以及從二維等高線恢復地形數據，都要用到輪廓線重建技術。
從輪廓線重建表面的研究早在20世紀70年代就開始了，近幾十年來，不斷有新的研究成果發表，從簡單到復雜，不斷深入。為了簡化研究，一般假設參與重建的輪廓線都是位于相互平行的平面序列上，這是符合實際應用背景的。另外，大量文獻把問題簡化為兩層平行輪廓線之間的重構問題。多層重構問題看做多個相鄰層重建結果的組合[1~3]。當相鄰層輪廓線的連接關系比較簡單，如每層上有一條輪廓線時，現有文獻提出許多有效的重建算法。這些算法能夠在輪廓線間存在的大量連接方式中選擇合理的、可接受的連接形式作為重建網格模型。這種選擇過程可以使用最優化方法來實現，優化目標可以是生成表面的面積最小、包圍的體積最大等。搜索全局最優解效率不高，一些啟發式的算法被提出來以求得局部最優解，如最短對角線法、最大體積法、相鄰輪廓線同步前進法等。這些算法雖然求得的不是全局最優解，但是效率高 重建結果與最優解差別不大，也是可接受的。
當相鄰層有多條輪廓線參與重建，連接關系比較復雜時，輪廓線重建技術存在分支和對應問題圍繞這兩個問題有許多的文獻發表[1，4~6]，但至今沒有十分理想的解決方案。文獻[7，8]對輪廓線重建技術領域的研究工作進行了較好的綜述性陳述。
輪廓線重建問題簡化為兩層平行輪廓線連接問題后，如前所述，依照參與重建的兩層相鄰輪廓線的數量，可以把輪廓線重建問題分為一對一、一對多或多對多。現有的文獻對一對一的輪廓線重建問題提出了許多很有成效的解決方案，有的文獻在此基礎上假設一對一的輪廓線重建問題已經解決[1，9]。在一對一的輪廓線重建算法的基礎上，研究者開始了對一對多或多對多重建問題的探索[8~13]。相對一對一的重建問題而言，一對多或多對多重建問題的困難之處在于輪廓線對應和分支問題[1，3，14]。可以認為，輪廓線重建技術中所需要進行研究的核心課題就是對應問題和分支問題。盡管這兩個問題至今尚未完全解決，但是許多文獻從各個不同角度提出了在不同情況下適用的解決方案[1，2，6，13]。
值得指出的是，將解決二維曲線演化問題的Level Set 方法推廣為解決三維曲面的演化方法（基于PDE模型）[8，15，16]，從而使得重建表面更為光滑，拓撲結構更穩定等，是輪廓線重建課題的一種全新方法。然而本質上，Level Set也不可能完全繞開對應和分支問題。具體地說，給出一個輪廓線集合（從二維醫學切片中提取），基于通常的Level Set方法重建的表面，由于對應問題的存在，重建的表面可能與真實表面在拓撲結構上不一致，這在許多應用如醫學應用上是不能接受的。
對應和分支問題中，對拓撲結構影響最大的是對應問題。
事實上，對于給定的輪廓線重建問題，可能有多個非奇異解，即不能惟一定義表面，脫離原始體數據單獨考慮輪廓線重建問題，在這種情況下非常困難。本文通過引入原始切片體數據信息，幫助解決輪廓線對應問題。實驗表明這種方法是有效的。本文的重建方法是這樣進行的：首先用人機交互的方式提取醫學圖像切片上參與重建的輪廓線；然后用插值法求出垂直方向切片，通過垂直切片上的邊緣信息，引入邊緣距離圖這一輪廓線相似程度度量機制，來解決對應問題；最后將參與重建的兩層輪廓線按對應結果分組，再用改進的BPLI算法重建表面。
1 對應問題和分支問題及新的解決方法:分組技術
對應問題是指異層輪廓線之間是否存在連接關系 即是否對應。 如圖1所示 共有三層輪廓線參與重建表面 相鄰兩層之間的輪廓線都存在連接關系 即相鄰兩層之間的輪廓線兩兩互相對應。分支問題是研究輪廓線出現一對多或多對多的對應關系時如何進行具體的分支連接，如圖1所示 頂層的兩條輪廓線與中間層的一條輪廓線對應 需要解決如何進行分支連接的問題。 圖1展示了一種合理的分支連接方案。
考慮參與重建的兩層輪廓線分別位于S1和S2層上，顯然，可以通過對應關系將參與重建的輪廓線分組。如圖2所示，組1中的S1層上的一條輪廓線與S2層上的兩條輪廓線都存在對應關系，而組2中的S1和S2層上各有一條輪廓線，它們互相對應，而組1與組2之間的輪廓線互相不對應。經過分組后，再采用一些直接從輪廓線重建三維表面的方法 如能較好地解決分支問題的改進后的BPLI方法，可以很好地解決重建表面與實際表面不符的問題[14]。
2 將輪廓線三維建模與多尺度二維圖像處理技術
相結合的新建模方法
2.1 多尺度初始輪廓線提取
初始輪廓線的提取工作需要使用交互式的方法來實現。最初，通過Canny算子提取原切片圖像的輪廓線。不過，這種自動提取的輪廓線是不閉合的，而且有許多是重建中不需要的輪廓線。要把感興趣的輪廓線區分出來需要很好的物質分類技術，目前還很難實現。本文用以下交互式方法來解決這一問題：設計一個能夠手工提取輪廓線的程序，由操作人員逐點進行邊緣提取，除了必要的一些編輯操作，如撤銷、存盤之外，以下兩種方法的使用很好地提高了手工提取的效率：
a) 自動校準。為了重建結果準確，希望提取的輪廓線與真正的邊緣盡量吻合。這一工作完全由手工做比較耗時，本文設計的輪廓線提取程序在每次邊緣點被手工選取之后，提取程序能夠將它與Canny算子提取的邊緣比對，然后選擇最為鄰近的Canny邊緣點。這一功能使得邊緣提取的手工操作可以更加簡便。
b) 邊緣跟蹤。邊緣跟蹤能夠提高輪廓線提取的自動化程度。對于連續性較好的輪廓線，使用邊緣跟蹤算法進行自動提取，既提高了速度，又增加了準確性。許多文獻對邊緣跟蹤算法進行了探討，本文選擇一種基于動態規劃機制的邊緣跟蹤算法。給定一段連續性較好的邊緣上的兩點xA與xB，該算法能自動提取出xA和xB間的邊緣線。
2.2 插值法求垂直方向切片
為了簡便起見，本文采用線性插值方法求垂直方向切片。不妨設原圖像切片平行于xy平面，一般而言，在醫學圖像中，原序列切片圖像在z方向的分辨率要大大低于xy平面的切片分辨率。為了得到垂直于xy平面的切片圖像，需要用線性插值的方法在z方向上求插值切片。通過這種方法能夠得到在z方向上與xy平面上具有相同分辨率的體數據之后，可以很方便地求任意方向的插值切片。本文只需要求垂直于xy平面且平行于x軸的切片序列。圖3(a)是原始切片示意圖，它們平行于xy平面；(b)是沿x軸的插值切片序列。
垂直方向切片的信息將用來幫助解決輪廓線對應問題。確切地說，使用垂直方向切片上的邊緣信息來解決輪廓線對應問題。這一方法的思路是:圖4(a)中有C0和C1兩層輪廓線。其中C0層上只有一條輪廓線a，而C1層上有兩條輪廓線b和c。這是一對多的輪廓線重建問題，因而需要解決對應和分支問題。根據對對應和分支問題的定義，首先需要通過輪廓線分組解決輪廓線對應問題，也就是要回答輪廓線a與b、c的對應情況。如圖4所示，本文利用垂直切片來解決該問題：如(b)所示，當垂直切片S上的邊緣連接a與b時，且沒有邊緣連接a與c時，說明輪廓線a只與b對應；(c)的情況則說明輪廓線a只與c對應；而(d)則說明a與b、c都對應，因為a、b與a、c之間都有S上的邊緣連接，此時出現一對多的對應問題，因而存在分支結構。
由于邊緣定位不可能非常精確，同時插值得到的垂直切片含有許多噪聲，本文采用多個縱向切片，結合下面將要闡述的一種閾值機制來進行輪廓線分組；同時，由于噪聲的存在，邊緣的連接性被破壞，垂直切片的邊緣也不可能總是與平行輪廓線重合。為了判斷兩條輪廓線之間是否有垂直切片上的邊緣連接，借助一種邊緣距離圖來進行。
2.3 構造邊緣距離圖
通過Canny算子，自動提取垂直切片的邊緣，然后通過邊緣跟蹤算法提取出所有長度大于一定閾值的邊緣。在本文中，這一閾值取5。由于這一閾值在小范圍內調整對結果影響不大，本文不對它進行討論。這樣，得到了一個邊緣序列C1，…，Ck。
定義1 以邊緣線C1，…，Ck為頂點，以C1，…，Ck之間最小歐氏距離為相應頂點之間邊的權值的無向完全圖，稱為邊緣距離圖G。
通過G可以計算出G中任何兩個頂點Ci與Cj之間的最短路徑，即任何兩個邊緣線之間的最短路徑，路徑長度記為phCiCj。
定義2 設S是一垂直切片，C是S上的邊緣線，A是S上任意一點，記disAC是在S上點A到線C上的最小歐氏距離。
定義3 設G是垂直切片M的邊緣距離圖，M的頂點是G上的邊緣線C1，…，Ck，A、B是M上任意兩點。其中：CA是離A最近的邊緣；CB是離B最近的邊緣。定義A、B之間的對應距離為
CoAB=disACA+disBCB+phCACB
對應距離度量的性質如圖5所示。其中：(a)是垂直切片的輪廓線圖；(b)的兩條白色橫線表示相鄰兩水平切片的位置。定義的對應距離必須滿足如下性質，即由于a和c在同一邊緣線上，則a與c的對應距離為0，小于a與b的對應距離，這一對應距離度量可以用來判斷a應該與b、c中的哪一條對應。具體的算法還要結合一種閾值機制進行，2.4節是詳細的算法步驟。
2.4 輪廓線分組步驟
2.3節闡述了輪廓線分組的基本思路，以下則是算法的詳細步驟。
假設要重建的輪廓線是閉合的、不相交的。在S1層有m條輪廓線，在S2層有n條輪廓線O1，…，On。
首先要計算S1層上輪廓線每一頂點的對應號。不妨設頂點A是待求對應號的點，方法如下：
a)找到垂直切片M，滿足A∈M∩S1。
b)生成M的邊緣距離圖，并計算M的每對邊緣之間的最小路徑。
c)找到點P1，…，Pj，滿足Pj∈∪nl=1Ol∩M(1≤j＜i)。其中：i是滿足條件的交點數。
d)計算對應距離CoAPj（1≤j＜i）。
e)找到Pj，滿足CoAPj=minil=1 CoAPl。
設Pj∈Oy，則點A的對應號為y。
當S1層上輪廓線C上所有點的對應號計算出來之后。定義S1層的輪廓線C和S2層上輪廓線Ox(1≤x≤n)的相似度為RCO=NCO/NC。其中：NCO是輪廓線C上所有對應號為x的點的個數；NC是C上所有點的個數。
如果RCO＞θ，θ為一個閾值，則認為C與Ox對應。這一閾值機制由于綜合了C上所有點的信息，提高了抗噪聲的能力。
如果在S2層上有多條輪廓線與C對應，由于相似度被多條輪廓線分享，這一閾值機制可能會失效。然而這一問題可以用如下方法簡單地克服：
對S2層上的每一條輪廓線，用以上方法檢查它在S1上的對應輪廓線，這樣丟失的對應關系被重新檢索到。
需要指出的是，這種對應問題檢索機制要求假設S1或S2層上的任何一條輪廓線至少能與另一層上的一條輪廓線對應。
對應關系確定之后，相鄰層的輪廓線可以根據對應結果分組。這樣，相似層的多對多輪廓線被分成幾個輪廓線組。在一個組內不同層的輪廓線一定互相對應，而不同組內的不同層輪廓線一定不對應。
2.5 輪廓線的三角網格重建
每一組輪廓線都獨立地使用文獻[14]的新方法進行三角網格重建，每個組的重建結果的并集為最終重建結果，重建方法在文獻[14]有詳細的討論，在此不再描述。需要指出的是，通過實驗發現，盡管對大部分數據直接用文獻[14]的方法能夠得到可接受的表面，但是存在一些數據，特別是拓撲結構比較復雜的數據，直接用文獻[14]的方法求解會因為對應問題沒有正確求解而導致拓撲結構錯誤。筆者將在實際結果中舉例說明這一情況，并對比結合輪廓線分組之后的結果來說明本文方法在解決對應問題方面的作用。
3 實驗結果
本章所使用的實驗數據是Visible Human Project中的CT切片圖像。為了測試本文的對應問題判別方法，首先設計了前文所述的輪廓線提取軟件，交互式地提取切片上感興趣區域的輪廓線；然后用線性插值的方法沿z軸求出垂直切片。生成的垂直切片用于計算邊緣距離圖，從而計算不同層的每對輪廓線之間的相似度。所有的算法都用Visual C++在Windows平臺上實現。
在實驗中取域值θ=0.6。實驗結果表明，對于直接用文獻[14]基于BPLI方法不能有效處理的對應問題，本文方法能給出正確的重建結果。
圖6是輪廓線對應問題判別算法的不同步驟。其中(a)和(b)是初始水平切片S1和S2。S1中提取出三條輪廓線，即(c)中的a1、a2、a3 而S2中提取出兩條輪廓線，即(c)中的b1、b2。(e)是某一垂直切片與S1和S2分別相交于直線a和b；(d)則是求出對應問題后的分組結果。S1和S2層的五條輪廓線被分成兩組V1和V2。
輪廓線分組對于重建結果準確性是很有幫助的。例如，圖7是某兩層數據的實際結果。其中(a)中輪廓線a1，a2∈S1，b1，b2∈S2，分組之后的重建結果如(c)所示，而不分組直接使用BPLI方法重建結果如(b)所示。圖7(b)由于對應錯誤而造成了嚴重的拓撲結構錯誤。
圖8是本文方法重建的一個比較復雜的三維表面，數據來源是Visible Human Project，整個重建結果包含19層連續的輪廓線。
4 結束語
本文提出結合原切片圖像信息的對應問題解決方法，并以此提出結合輪廓線分組的輪廓線重建方法。實驗表明這一方法解決對應問題是有效的。同時，由于這一方法引入了垂直切片的邊緣信息求解對應問題，理論上非相鄰層的對應問題也可以得到合理的解決。
垂直的插值切片不可避免地引入噪聲，對于某些能夠提供垂直方向的原始切片的醫學圖像來說，如MRI，這一方法應該更加適用。因此，本文的工作只是一個開始，進一步的工作包括MRI圖像的應用研究和不相鄰層輪廓線對應問題的應用研究等。
參考文獻：
［1］BAJAJ C L COYLE E J LIN K N. Surface and 3D triangular meshes from planar cross sections[C]//Proc of the 5th International Meshing Roundfable.[S.l.]: Sandia National Laboratories 1996:169-178.
[2]BAREQUET G SHAPIRO D TAL A. Multilevel sensitive reconstruction of polyhedral surfaces from parallel slices[J].The Visual Computer 200016(2): 116-133.
[3]HE Jin-guo ZHA Hong-bin SHI Qing-yun. Reconstruction of surfaces from medical slices using a multi-scale strategy[C]//Proc ofIEEE International Conference on Systems Man and Cybernetics.2001: 1993-1998.
[4]GEIGER B. Three-dimensional modeling of human organs and its application to diagnosis and surgical planning Technical Report 2105[R].1993.
[5]SHINAGAWA Y KUNII T L. The homotopy model: a generalized model for smooth surface generation from cross sectional data[J].The Visual Computer 19917(2-3): 72-86.
[6]BAREQUET G VAXMAN A. Geometric computing for volumetric data I: nonlinear interpolation between slices[C]// Proc ofACM Symposium on Solid and Physical Modeling，SPM’07. 2007: 97-107.
[7]MüLLER H KLINGERT A.Surface interpolation from cross sections[M]//Focus on Scienlific Visualization. London:Springer-Verlag 1993: 139-189.
[8]NILSSON O BREEN D MUSETH K. Surface reconstruction via contour metamorphosis: an eulerian approach with lagrangian particle tracking[C]// Proc of IEEE VIS. 2005:407-414.
[9]MEYERS D SKINNER S. Surfaces from contours[J].ACM Trans on Graphics 199211(3): 228-258.
[10]BOISSONNAT J. D GEIGER B. Three dimensional reconstruction of complex shapes based on the delaunay triangulation Technical Report 1697[R]. 1992.
[11]EKOULE A. B PEYRIN F C ODET C L.A triangulation algorithm from arbitrary shaped multiple planar contours[J].ACM Trans on Graphics 1991 10(2): 182-199.
[12]BAREQUET G GOODRICH M T LEVI-STEINER A et al. Contour interpolation by straight skeletons[J].Graphical Models 200466(4): 245-260.
[13]ZHANG Hong JOHN N YIM P J. Surface reconstruction from orthogonal contours[C]// Proc of SPIE. 2006: 98-104.
[14]何金國，查紅彬，基于BPLI從二維平行輪廓線重建三維表面的新算法[J]. 北京大學學報，200339(3): 408-411.
[15]BRAUDE I. Smooth 3D surface reconstruction from contours of biological data with MPU implicits[D]//[S.l.]:University Crossings，2005.
[16]WANG De-sheng OUBAY H O KENNETH M K et al. Efficient surface reconstruction from contours based ontwo-dimensional Delaunay triangulation[J]. International Journal for Numerical Me-thods in Engineering 2006 65(5):734-751.
[17]PITAS I. Digital image processing algorithm[M].[S.l.]: Prentice Hall，1993.
