螺旋CT的三維空間面皮參數平面化方法

摘 要:將三維空間面皮轉換為二維平面可以使獲取面皮特征、顱骨與面皮的關聯信息更方便。通過改進的輪廓線提取方法剔除原始切片冗余數據，對重構得到的面皮數據進行空間坐標系的規范，并以最近鄰區域(NNR)劃分方法將空間面皮的點云數據點集化，最后結合新劃分的平面點數據集的特點提出了線性插補的參數平面化方法和基于鏈表的參數平面化方法。實驗證明兩種方法的效果良好，驗證了其正確性和可行性。

關鍵詞:螺旋CT;輪廓線;最近鄰區域劃分;參數平面

中圖分類號:TP391.41 文獻標志碼:A

文章編號:1001-3695(2010)03-1195-03

doi:10.3969/j.issn.1001-3695.2010.03.110

Methods of spiral CT three-dimensional face parametric plane

MAO Xiao-lin1， GENG Guo-hua1， ZHOU Ming-quan2

(1. College of Information Science Technology， Northwest University， Xi’an 710127， China; 2.College of Information Science Technology， Beijing Normal University， Beijing 100875， China)

Abstract:Three-dimensional space face into two-dimensional plane to obtaining face features，information between skull and face is good efforts.This paper improved extraction methods to contour the raw data to remove redundant data，reconstruct and adjust face spaces coordinate system，coordinating of the point cloud data after unification select a set of nearest neighbor region divided methods.According to characteristics of data plane，gave a linear interpolation parametric plane method and a based list parametric plane method.The experiment proves that two methods have good effects，and proves its accuracy and the feasibility.

Key words:spiral CT; contour; nearest neighbor region divided; parametric plane

0 引言

人臉面皮特征是顱面形態學的一個研究熱點，對人類學、醫學整形、法醫學等都有很大的貢獻。大多數研究者在三維空間中對面皮數據進行分析，這樣常會出現以下問題:在空間中對于面皮點的選擇依靠主觀性和研究者視覺判斷，造成標志不一致問題;三維空間中信息的提取比較困難，同時導致準確率降低等。將三維空間面皮轉換到二維平面進行研究可以避免這些問題。文獻[1]中提到了將人顱骨面皮進行參數化平面表示，但該文獻中要求CT數據是序列掃描，即圖像彼此所在的平面是平行的。然而高精度的CT數據則是采用螺旋掃描方法獲得，其特點是各層圖像數據之間有重復和發生重疊現象，因此每個切片的數據是不平行的。這種數據對醫學研究的作用很大，但對于顱面形態學的研究是不方便的。本文工作就是將螺旋CT得到的三維空間面皮數據參數平面化。

1 參數平面化思想來源

頭部面皮可以看成是空間三維的一個簡單曲面。從數學和幾何角度講，此類三維空間簡單曲面都是二維流形的三維嵌入問題，即該三維曲面每個點都與某參數對之間存在一個對應關系，因此，三維空間曲面轉換為二維平面的出發點就是設計并計算出這樣一種映射，同時不丟失該曲面的三維幾何信息。但是，在笛卡爾坐標系下，對于封閉曲面無論怎么設計也不有一個坐標變量對與第三個坐標變量存在一個對應關系。然而，對此類封閉三維空間曲面采用球面極坐標系問題則有所突破。在極坐標系下球面上任意一點極徑的長度r與兩個夾角數據對(φ，α)之間存在一個映射關系，其中φ∈[0，2π]，α∈[0，π]。由此可知，雖然不能將球面展開為在一個空間平面上，但可根據這種對應關系將球面展開在其參數所變化的仿射區域內。三維笛卡兒坐標系下部分球面，如圖1(a)所示。根據式(1)，定義V1、V2兩個向量V1=(0，0，1)T，V2=(x-0，y-0，z-0)T，則α=arccos((v1，v2)|v1||v2|)，可計算函數r=f(φ，α)。

x=r sin α cos φy=r sin α sin φ(1)

z=r cos α

于是得到參數平面如圖1(b)所示。其中(φ，α)點顏色值與r大小對應。圖1(a)是球面一部分，r是定值，因此顏色相同。參數化平面展開的原則是:

a)三維曲面上的點經過適當的映射，使該曲面上的坐標點在其二維參數空間中的參數對與第三個坐標變量之間有一種函數關系。

b)函數關系是一對一的，可以經逆變換將二維參數平面中的點轉換為三維空間坐標系中的點。

然而空間坐標系中的面皮數據比較復雜，特別是耳輪處的數據，造成同一對(φ，α)對應不同的點，因此直接使用球面極坐標將無法滿足一對一的關系。本文將設計其他參數化方法解決此問題。

2 目的和流程

2.1 參數平面化的目的

將三維空間面皮數據參數平面化的目的是:參數平面上可以迅速獲得三維空間兩點之間的參數距離，空間任意點到三維幾何中心的參數距離; 將三維面貌特征點直觀地展現到二維平面上，便于拾取;對不同的頭部數據，歸一化后直接進行參數平面的對照得到異同信息。

2.2 螺旋CT參數平面化流程

螺旋CT數據參數平面化的方法步驟為:

a)輪廓線提取，剔除CT冗余;

b)重建單層面皮，獲取點云數據;

c)規范空間坐標系;

d)平面點集劃分，面皮點云數據合理劃分到系列平行的平面上，每個平面形成一個像素點的集合;

e)參數平面化。

本文具體處理流程如圖2所示。

3 螺旋CT預處理與點集劃分

3.1 改進的輪廓線提取算法剔除冗余影像

獲取單層面皮需對螺旋CT進行邊緣檢測和輪廓線提取。文獻[2]中算法健壯性差，容易造成邊緣點數據丟失，本文在此基礎上提出改進的輪廓線提取算法。

首先對CT圖像進行快速模糊C均值聚類(FCM)分割[3]。FCM是一種典型的無監督模糊聚類方法，此方法是將圖像分為c(c≥1)類并求每類的聚類中心，使得用隸屬度函數定義的目標函數達到最小。FCM法在聚類過程中不考慮圖像的灰度分布，因此具有很好的魯棒性。它的目標函數為

JFCM=∑j∈N∑Ck=1μqjk‖yj-vk‖2，1

其中:N代表圖像中像素點的集合;y1代表圖像中像素點j的灰度;μjk代表像素點j對于第k類的隸屬度;vk代表第k類的聚類中心;‖yj-vk‖代表像素j與第k類的距離。FCM方法可以有效地防止邊緣點丟失。對 FCM分割后的二值化圖像輪廓線進行提取。每個點是邊界像素點需要滿足兩個條件:

a)邊界點的自身灰度值為1;

b)邊界點的四鄰域內至少有一個像素點的灰度值為0。

設一個棧，棧底是輪廓線的初始點，獲得初始點的方法就是:過CT圖像幾何中心作一條平行于冠狀面的直線，與坐標軸的交點記做o，由o向右前進，碰到第一個灰度值為1的像素點定為初始點放入棧底。以初始點為中心進行如下步驟:

a)初始點開始逆時針進行掃描(圖3)，中心圓點為初始點，那么點7是灰度值為0的像素點，面皮數據是閉合的，所以八個鄰域內必有灰度為1的點。從點0、1、2、3、4、5、6的順序開始搜索灰度值為1的點。

b)找到灰度為1的點后進棧，設置為新的中心點。此時將掃描線由最近一次的掃描方向順時針轉90°后作為新中心點的掃描方向，繼續逆時針掃描，這樣不會漏掉相鄰且灰度為1的點。

c)重復b)，假設當一個點為中心點掃描時，八鄰域內僅有一個灰度為1的點并且這個點之前已經進棧，則說明此點為一個離散點，不考慮其存在;然后取出棧頂像素點作為中心點，將掃描線逆時針旋轉45°后繼續掃描。當碰到棧頂與棧底點為同一個點時掃描結束。

從棧底到棧頂的連續點得到了CT圖像輪廓線。根據棧內點的坐標值，保留對應在其原始CT中點的灰度信息，剔除其他冗余影像。CT圖像處理過程如圖4所示。

3.2 重構及規范坐標系

螺旋CT數據在采集過程中會造成圖像的位置和角度等發生變化，切片和冠狀面有個角度。利用CT切片提取等值面后進行表面三維重建，不受坐標數值以及之間切片平面相互關系的影響，使得利用螺旋掃描與序列掃描在此類處理上沒有差別。

面繪制算法Marching Cubes[4]重構后對面皮點數據相對坐標和拓撲序列沒有影響。去掉重構后的三角面片信息，保留面皮的點云數據，如圖5所示。不失一般性，空間點云數據使用三維笛卡爾坐標系，坐標系Z軸與人腦的冠狀面垂直。對在三維空間中的面皮點云數據采用包圍盒[5]的方法。

包圍盒幾何中心x、y、z坐標分別平移tx、ty、tz(式(3))，使包圍盒的幾何中心與圖像顯示屏幕空間的中心重合，并設定其為空間坐標系的原點，使模型在屏幕空間的坐標可以進行旋轉、放大、縮小及其逆變化。對于此問題的研究都進行同樣的坐標系規范，有利于不同套數據的歸一化。

[x′，y′，z′，1]T=1 0 0 tx0 1 0 ty0 0 1 tz0 0 0 1[x，y，z，1]T(3)

3.3 最近鄰區域(NNR)劃分點集

參數平面化前提是將面皮點云數據劃分到一系列平行于人腦冠狀面的平面上，使得每個平面上的點數據形成一個集合。面皮點數據在正視場景時，包圍盒確立模型在坐標系X、Y、Z軸的范圍，高度在Z軸的坐標范圍記為[θmin，θmax]。X、Y軸的單位為原始切片中相鄰像素點之間的距離D。動態選取劃分集合的個數為N，于是Z軸單位就是d=(θmax-θmin)/N，等間距的平面表示為集合{θ1，θ2，θ3，…，θN}，其中，θmin=θ1，θmax=θN。本文根據文獻[6]提出了最近鄰區域劃分思想對點數據劃分集合。

NNR劃分思想:要將給定的點數據p1p2…pN劃分到已知的若干區域(集合)中。對于任意點pi計算其到每個區域的距離，將pi劃分到距離最近的區域(集合)中。根據NNR劃分思想對點數據劃分的具體算法如下:

Define p(z)為p點Z軸坐標值

for i=1:N

if 1≤i≤N-1

then 掃描在θi和θi+1面之間的點

if p(z)-θi<θi+1-p(z)

then p(z)=θi，p點劃分到θi面

else p(z)=θi+1，p點劃分到θi+1面

end

end

end

對于平面點集個數N較大時，劃分所造成點位置移動較小，對于后續研究產生的誤差不大。劃分后的平面點數據如圖6所示。

4 三維空間面皮數據參數平面化

針對劃分后平面點數據的分布情況，本文在參數化平面的原則下設計了基于線性插補的參數平面轉換方法和基于鏈表的參數平面轉換方法。

4.1 基于線性插補的參數平面轉換方法

θi可當做其中一個參數，由于有少量數據如圖6(b)和(c)都是不閉合的，平面角φ不能作為另一個參數。設平面上點數據灰度值為0，利用線性插補方法將不閉合的曲線段插補成一條連續的曲線。

對圖6(b)的情況，過斷開的曲線段的兩端點A、B畫一條直線，方程為yA-yB=kAB(xA-xB)，以端點A為中心點，在八鄰域內找一個插補點。選取插補的點(x，y)方法:灰度值為1;過(x，y)點斜率為-1/kAB的直線與線段AB有交點，且到AB垂直距離最短。此時令插補點灰度值為0.5，以新得到的插補點做中心點繼續上面的操作，直到當某一個中心點八鄰域內離直線最近的點是線段另一端點B時停止。經過線性插補后的平面數據成為閉合并包圍平面原點的曲線。

平面數據如圖6(c)有這樣的特征:曲線和過平面原點的任意一條直線的交點(不包括相切的情況)不超過一個。處理此情況的方法是將曲線點數據關于切片原點對稱復制，即點(x，y)經過原點對稱復制后得到點(x，y)，它們滿足

yx=YXx2+y2=X2+Y2 (4)

然后用線性插補法使閉合曲線包圍平面原點。于是對平面θi上的點按照平面極坐標的思想，X軸與閉合曲線交點為Q0，曲線上的點順時針記為{Q0，Q1，Q2，…，Qm，Q0}。則累加得到閉合的弧長Arc=∑mk=1|Qk-Qk-1|+|Qm-Q0|，設每個點對應的平面角為φk，且設Q0=Qm+1，φ0=0，φm+1=2π，有φk=φk-1+(|Qk-Qk-1|)/Arc。線性插補法參數化平面具體步驟如下:

a)對于平面數據中不閉合、離散數據段的情況，用直線插補的方法處理為閉合且包圍原點的曲線。

b)順時針計算相鄰兩點距離，累加求得弧長。根據公式求出每個點對應的φk。三維空間坐標系的中點對坐標單位可以得到比例系數K=dD，則r=x2+y2+(Kθi)2，設計函數f(Qi，θi)=r。

c)對于Q、θ都采取逼近的方式得到一條連續的曲線。

4.2 基于鏈表的參數平面轉換方法

平面θi上任意一點Pj，尋求函數關系F(p，θ)=r。對于面θi設置一個鏈表p，中心點是平面坐標系原點，步驟如下:

a)原點從正左方向開始順時針掃描周圍第一層的8個點，然后掃描第二層的16個點，同理繼續掃描其他層。同時設置一個層變化計數器M和每層像素點數的計數器N，第一層的(M，N)對應的數據就是(1，8)，第M層是(M，8×M)。

b)當碰到灰度值為0的數據點時鏈表p動態增加，此時下標為j，根據(M，N)計算出該點到所在切面原點的距離Lj，在M層上，N個像素與原點距離有M+1種情況:

DM2+W2，W={0，1，2，…，M}

設置臨時變量T從1到N累加:

A:0→M LT=DM2+W2

W={0，1，2，…，M}

B:M→0 LT=DM2+W2

W={M，M-1，…，0}

交替執行A、B直到T=N，停止本層的搜索。

c)重復進行b)，直到鏈表節點數等于該平面上數據點總數時停止。

綜上，可以得到這樣一個函數關系:

F(pj，θi)=(pj)2+(Kθi)2(5)

面θi所對應的點數[1，ni]歸一化到[0，1]。得到的數據鏈和θ用一種逼近的方法得到新的曲線。

5 實驗結果

實驗數據是西門子16排螺旋CT螺旋掃描采集到的306張512×512 CT切片。利用工具VC++6.0 與動態鏈接庫VTK、OpenGL重構面皮點數據為349 508個， 劃分1 600張平面點集合，給出了兩種方法的參數平面展開結果，如圖7所示。

為了形象化，面皮展開時采取數據中點對齊。參數平面上顏色信息代表空間面皮上的點與三維空間坐標原點的距離。亮度越高，表示距離越大，黑色說明距離為零，即沒有數據。根據實驗結果，對兩種方法的優缺點、適用性進行對比:

1)基于線性插補的參數平面轉換方法

a)優點:研究僅關心三維空間數據與參數平面上所對應點的關系，因此構造參數平面時的插補點信息不造成影響;參數化平面的橫坐標從原點向正方向的顏色信息反映的是在三維空間面皮數據中過X軸平行于矢狀面的平面順時針旋轉時面皮展開情況。

b)缺點:產生中間信息，計算量大。

c)適用性:用于測量特征點之間的歐式距離;在對比不同人面皮的整體或局部區域的幾何信息時，將平面圖像歸一化后直接進行疊加以灰度差異反映異同;有利于后續研究顱骨與面皮之間的關聯關系和配準。

2)基于鏈表的參數平面轉換方法

a)優點:對幾種類型的平面數據均適用;在計算每個點與三維空間坐標原點的距離時用(M，N)可以快速得到，沒有中間信息，計算量少。

b)缺點:對平面點數據進行了重排列，在進一步研究時需要進行重新調整。

c)適用性:明顯反映面皮特征分布，獲取特征點時更加方便。

以上兩種方法都保持與三維空間面皮點數據的一一對應，參數平面上點的灰度值代表了在空間中的位置信息。

6 結束語

本文將面皮數據參數化展開成二維圖像，充分利用螺旋掃描的CT數據高精度和適用性的特點，同時根據不同的給定集合平面的個數來動態實現。在不同數據進行對比時的完整性和歸一性，方便地實現對比整個面皮或者部分面皮數據，從而研究各自特征，能更加迅速地確定面皮上某些特征點之間的差異，這些都對研究顱面形態學有著重要的意義。

參考文獻:

[1]董建民，周明全.空間三維面貌的二維參數平面轉化算法[C]//北京圖像圖形技術與應用會議.2008:380-385.

[2]唐振軍，張顯全.一種二值圖像邊界提取算法[J].圖象處理， 2006，22(10-3):281-283.

[3]丁震，胡鐘山，楊靜宇，等.FCM算法用于灰度圖像分割的研究[J].電子學報，1997，25(5):39-43.

[4]LORENSEN W E， CLINE H.E. Marching Cubes:a high resolution 3D surface construction algorithm[C]//Proc ofSIGGRAPH’87.1987:163-169.

[5]GREENE N.Detecting intersection of a rectangular solid and a convex polyhedron[M]//Graphics Gems IV.San Diego:Academic Press Professional Inc，1994:74-82.

[6]TAN S L，BELATON B，RAJION Z A，et al.Filtering for fast craniofacial surface[J].Reconstruction Archives of Orofacial Science，2006，1:29-35.