一種新的三維表面重構(gòu)方法

楊昆朋，張榮國，劉 霞，劉小君

(1.太原科技大學(xué)計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，太原 030024;2.邢臺(tái)學(xué)院數(shù)學(xué)與信息技術(shù)學(xué)院，河北邢臺(tái)054001; 3.合肥工業(yè)大學(xué)機(jī)械與汽車工程學(xué)院，合肥 230009)

一種新的三維表面重構(gòu)方法

楊昆朋1，張榮國1，劉 霞2，劉小君3

(1.太原科技大學(xué)計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，太原 030024;2.邢臺(tái)學(xué)院數(shù)學(xué)與信息技術(shù)學(xué)院，河北邢臺(tái)054001; 3.合肥工業(yè)大學(xué)機(jī)械與汽車工程學(xué)院，合肥 230009)

針對CT圖像，給出了一種三維表面重構(gòu)的方法。首先用中值濾波的方法進(jìn)行預(yù)處理后，根據(jù)曲率信息提取輪廓線上的特征點(diǎn)，然后按照最短對角線法進(jìn)行連接，其中在選取上下輪廓線兩點(diǎn)間最短距離時(shí)，采用歐式距離進(jìn)行實(shí)現(xiàn)，同時(shí)為了使構(gòu)造的三角形視覺效果較好，分別計(jì)算當(dāng)前點(diǎn)角的余弦值，從中選擇較小的。最后將該方法與原方法分別進(jìn)行實(shí)驗(yàn)，驗(yàn)證了使用部分輪廓特征點(diǎn)也能實(shí)現(xiàn)三維表面重構(gòu)，重構(gòu)效果較為滿意。

CT;輪廓線;曲率;最短對角線

三維表面重構(gòu)技術(shù)中主要兩類不同的方法實(shí)現(xiàn):體繪制和面繪制。由于體繪制的計(jì)算量偏大，而且對圖像計(jì)算和顯示的硬件要求比較高。因此都采用面繪制的方法進(jìn)行三維表面重構(gòu)。主要方法有Lorensen[1]在1987年提出了移動(dòng)立方體(marching cubes，MC)算法。MC采用等值面的思想，用三角形擬合一定灰度閾值的等值面。但MC算法在生成曲面的時(shí)候很容易產(chǎn)生二義性問題，因此針對這個(gè)問題，很多人對此做了大量的改進(jìn)。張榮國等[2]提出基于共軛梯度的B樣條主動(dòng)輪廓邊緣提取方法。王慶霞等[3]根據(jù)四叉樹結(jié)構(gòu)的層次化模型簡化生成三維地形網(wǎng)格，錢蘇斌[4]提出用行掃描線法提取帶孔的輪廓線，再結(jié)合最小對角線完成三角面片的連接，姚民倉等[5]為解決輪廓線的分支問題，提出了改進(jìn)的最短對角線法，李運(yùn)鋒等[6]針對輪廓線拼接的復(fù)雜性，提出方向包圍盒的方法進(jìn)行輪廓點(diǎn)的拼接，最后用最短對角線法實(shí)現(xiàn)表面的重構(gòu)。然而這些方法都采用了輪廓線上的所有像素點(diǎn)，然而在重構(gòu)表面時(shí)，有些卻不是必須的。

因此，在獲得物體的輪廓線之后，根據(jù)曲率的大小選擇合適的像素點(diǎn)，然后采用改進(jìn)的最短對角線法進(jìn)行三維表面的重構(gòu)。

1 輪廓線的提取及三維重構(gòu)

1.1 CT預(yù)處理

在獲取CT切片的過程中，由于各種因素的影響，使得圖像中存在一定的噪聲，圖像中目標(biāo)物體部分邊緣也可能局部不清晰。因此在獲得圖像的二值圖像后，采用水平集中Chan-Vase圖像分割模型的方法(式1)對感興趣的區(qū)域進(jìn)行提取:

其中，μ、v是非負(fù)常數(shù)代表各能量項(xiàng)的權(quán)重系統(tǒng)，μ0是初始圖像，u是逼近初始圖像μ0的一個(gè)分段圖像，Ω是圖像I(x，y)定義的區(qū)域，C是分割前景和背景區(qū)域的邊界曲線。

1.2 輪廓點(diǎn)提取與精簡

通過上述方法在獲得CT圖像的輪廓線之后，然后再采用逐行掃描的方法對輪廓線進(jìn)行掃描，獲

1.3 三維重構(gòu)

輪廓線連接表面重構(gòu)自70年代提出以來，有許多人從不同的角度提出了各種算法，如最小表面積法、最大體積法、同步前進(jìn)法、最短對角線法等。由于最短對角線法具有簡單、容易實(shí)現(xiàn)，而且當(dāng)上下兩條輪廓線的大小和形狀相近，相互對中情況效果是比較好的。

在相鄰兩條輪廓線中，設(shè)上層輪廓線上的點(diǎn)和下層輪廓線上的點(diǎn)分別表示為:Ci={p(i)，i=1，2，…m} Ci+1={q(j)，j=1，2，…n}，如圖1所示，設(shè)Pi是上層輪廓線Ci上的任意一點(diǎn)，在下層輪廓線Ci+1上找到距離Pi最小的像素點(diǎn)，設(shè)為Qj，則以跨距PiQj為基礎(chǔ)，如果PiQj+1＜Pi+1Qj，則連接Pi，Qj+1形成三角面片。得其像素點(diǎn)。由于每一條輪廓線線上都由許多的像素點(diǎn)構(gòu)成，然而在進(jìn)行三維表面重構(gòu)的時(shí)候并不是需要所有的輪廓線上的像素點(diǎn)，因此采用文獻(xiàn)[7]中的方法對輪廓線上的特征點(diǎn)進(jìn)行精簡，進(jìn)而獲得該像素點(diǎn)在直角坐標(biāo)系下的二維坐標(biāo)值。

采用文獻(xiàn)[7]中曲率的計(jì)算方法，具體如下:

圖1 最短對角線法示意圖Fig.1 The shortest diagonal method

在計(jì)算得到各個(gè)像素點(diǎn)的曲率值之后，人工進(jìn)行選取適當(dāng)?shù)姆秶?，使符合條件的像素點(diǎn)都落在選取的區(qū)間中。

采用最短對角線法連接三角面片時(shí)，在選取上下兩條輪廓線兩點(diǎn)間的最短距離時(shí)，采用公式(3)歐式距離進(jìn)行選擇。

其中Q.x、Q.y分別表示當(dāng)前像素點(diǎn)的坐標(biāo)值，P.x、P.y分別表示離當(dāng)前點(diǎn)最近的坐標(biāo)值。同時(shí)為了在三角化的過程中，取得較好的視覺效果，使構(gòu)造的三角形盡量為銳角三角形，避免出現(xiàn)如圖2中的情況:

圖2 最短對角線中一種特殊情況Fig.2 A special case of the shortest diagonal

因此分別計(jì)算∠piqjpi+1和∠qjpipi+1的余弦值，如果前者小于后者，則連接pi、qj+1，否則連接pi+1、qj.同時(shí)根據(jù)其歐式距離的大小來選擇上下輪廓線中距離最短的像素點(diǎn)，在連接像素點(diǎn)構(gòu)造三角面片時(shí)，采用Delaunay三角網(wǎng)逐點(diǎn)插入的方法構(gòu)造，采用如下所示的數(shù)據(jù)結(jié)構(gòu):

(1)點(diǎn)數(shù)據(jù)結(jié)構(gòu)，用于存儲(chǔ)像素點(diǎn)的坐標(biāo)，其定義如下:

typedef struct{double x，y，z;像素點(diǎn) x，y，z坐標(biāo)}POINT3D

(2)三角形數(shù)據(jù)結(jié)構(gòu)，用于存儲(chǔ)三角形三個(gè)頂點(diǎn)，其定義如下:

typedef struct{POINT3D p[3];三角片三個(gè)頂點(diǎn)的坐標(biāo)}TRIANGLE

采用此方法把所有的CT輪廓線進(jìn)行連接之后，構(gòu)造出三角網(wǎng)格模型，同時(shí)借助于OpenGL開放式三維環(huán)境，對構(gòu)成的三維網(wǎng)格進(jìn)行光照渲染等處理。

2 其算法的具體實(shí)現(xiàn)步驟

根據(jù)最短對角線法進(jìn)行三維表面重構(gòu)的具體算法步驟如下:

(1)利用公式(1)獲得CT的輪廓線，并把獲得的像素點(diǎn)的坐標(biāo)值存入二維數(shù)組中;

(2)根據(jù)公式(2)計(jì)算各個(gè)像素點(diǎn)的曲率值，并把符合設(shè)置曲率區(qū)間的像素點(diǎn)的坐標(biāo)值保存到另一二維數(shù)組中;

(3)在上輪廓線Ci中選擇一點(diǎn)Pi，并在下輪廓線Ci+1中選擇一點(diǎn)Qj，在選擇點(diǎn)Pi+1還是選擇Qj+1時(shí)，根據(jù)公式(3)分別計(jì)算其歐式距離D(PiQj+1)、D(QjPi+1)，并人工設(shè)置閾值范圍 T，同時(shí)計(jì)算∠piqjpi+1和 ∠qjpipi+1，如果 D(PiQj+1) ＜ T 且∠qjpipi+1的余弦值小于∠piqjpi+1，則保留并連接三角形△PiQjQj+1，否則連接△PiPi+1Qj;

(4)重復(fù)(3)，當(dāng)某條輪廓線上只剩下最后一個(gè)點(diǎn)時(shí)，另一條輪廓線上的點(diǎn)都與這個(gè)點(diǎn)進(jìn)行連接，直到把所有點(diǎn)都連入三角形為止。

(5)直到把所有的輪廓線處理完，構(gòu)造出三角網(wǎng)格模型。

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

實(shí)驗(yàn)平臺(tái):Intel(R)Pentium(R)CPU，主頻為2.13GHZ，內(nèi)存為2GB，操作系統(tǒng)為Windows 7，軟件環(huán)境:VS2008結(jié)合OpenGL，搭建應(yīng)用平臺(tái)。選取曲率的范圍為[-0.1，0.1]之間。對大小為256×256的CT斷層圖像進(jìn)行三維表面重構(gòu)實(shí)驗(yàn)。采用原始圖像Head 183張CT切片，其中圖3為感興趣區(qū)域，圖5(a)、圖6(a)表面重構(gòu)，圖5(b)、圖6(b)側(cè)面效果圖，圖5(c)、圖6(c)為加入光照后的三維表面;采用原始圖像Shoulder 166張CT切片，其中圖7為感興趣區(qū)域，圖9(a)、圖10(a)表面重構(gòu)，圖9 (b)、圖10(b)側(cè)面效果圖，圖9(c)、圖10(c)加入光照后的三維表面。其中圖4(a)為采用Head圖像上的所有像素點(diǎn)，4(b)為精簡后的像素點(diǎn)，圖8 (a)為采用Head圖像上的所有像素點(diǎn)，8(b)為精簡后的像素點(diǎn)，不同方法的比較見表1.

其中PR=(原始點(diǎn)數(shù)-精簡后的點(diǎn)數(shù))/原始點(diǎn)數(shù)。

圖3 Head感興趣區(qū)域Fig.3 Head interested region

圖4 Head像素點(diǎn)精簡前后Fig.4 The pixels of Head pixels before and after streamline

圖5 Head三維表面重構(gòu)(本方法)Fig.5 Head 3d surface reconstruction(The method in this study)

圖6 Head三維表面重構(gòu)(原方法)Fig.6 Head 3d surface reconstruction(The original method)

圖7 Shoulder感興趣區(qū)域Fig.7 Shoulder interested region

圖8 Shoulder像素點(diǎn)精簡前后Fig.8 The pixels of Shoulder before and after streamline

圖9 Shoulder三維表面重構(gòu)(本方法)Fig.9 Shoulder 3D surface reconstruction (The method in this study)

表1 不同方法的比較Tab.1 The comparison of different methods

圖10 Shoulder三維表面重構(gòu)(原方法)Fig.10 Shoulder 3D surface reconstruction (The original method)

4 結(jié)束語

通過對輪廓線上的像素點(diǎn)采用曲率進(jìn)行特征點(diǎn)的提取，然后設(shè)置曲率區(qū)間從中選擇適當(dāng)?shù)南袼攸c(diǎn)，最后用最短對角線法進(jìn)行三維表面重構(gòu)。實(shí)驗(yàn)證明，使用部分像素點(diǎn)用該方法也能重構(gòu)三維表面，而且通過使用這兩種方法分別進(jìn)行三維表面重構(gòu)來看，其效果也較為滿意。

[1]CLINE H E，LORENSEN W E，LUDKE S，et al.Two algorithms for the three-dimensional reconstruction of tomograms[J].Medical physics，1988，15(3):320-327.

[2]張榮國，劉小君，劉焜，等.基于共軛梯度的B樣條主動(dòng)輪廓邊緣提取[J].中國圖象圖學(xué)學(xué)報(bào)，2010，15(1):103-108.

[3]王慶霞，張榮國，武妍，等.四叉樹結(jié)構(gòu)地形網(wǎng)格簡化算法的研究[J].太原科技大學(xué)學(xué)報(bào)，2013，34(1):1-5.

[4]錢蘇斌.基于輪廓線的任意形體三維重建[J].成都大學(xué)學(xué)報(bào):自然科學(xué)版，2013，32(3):262-266.

[5]姚民倉，秦新強(qiáng)，胡鋼.由斷層圖像輪廓重構(gòu)三維表面的新算法[J].微計(jì)算機(jī)信息，2009(25):214-216.

[6]李運(yùn)鋒，劉修國.基于方向包圍盒投影轉(zhuǎn)換的輪廓線拼接算法[J].計(jì)算機(jī)應(yīng)用，2011，31(12):3353-3356.

[7]HE X C，YUNG N H C.Corner detector based on global and local curvature properties[J].Optical Engineering，2008，47(5): 8-12.

A New Method of 3d Surface Reconstruction

YANG Kun-peng1，ZHANG Rong-guo1，LIU Xia2，LIU Xiao-jun3
(1.School of Computer Science and Technology，Taiyuan University of Science and Technology，Taiyuan 030024，China;2.School of Mathematics and Information Technology，Xingtai University，Heibei Xingtai 054001，China; 3.School of Mechanical and Automotive Engineering，Hefei University of Technology，Hefei 230009，China)

a method of 3d surface reconstruction for CT images was proposed.The feature points on the contour are extracted based on the curvature information.Then feature points are connected based on the shortest diagonal method.When selecting the shortest distance between the up and down contour points，euclidean distance was used to achieve.At the same time，the smaller current point angle was choosen to make the structure of the triangular visual effected better.Finally，the experiments verified that using partial contour feature points also can realize 3d surface reconstruction，and the result is satisfactory.

CT，contour，curvature，the shortest diagonal

TP391

A

10.3969/j.issn.1673-2057.2015.04.009

1673-2057(2015)04-0284-04

2014-12-29

國家自然基金(51375132);山西省自然基金(2013011017);山西省研究生創(chuàng)新項(xiàng)目(20143096)

楊昆朋(1988-)，男，碩士研究生，主要研究方向?yàn)橛?jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)。