對血管切片做三維重建的一種方法

1、<p>　　對血管切片做三維重建的一種方法 </p><p>　　作者：王斯剛　馮有前　趙學(xué)軍　王錦江 </p><p>　　【關(guān)鍵詞】 血管切片 </p><p>　　關(guān)鍵詞: 血管切片;中軸線;搜索算法;內(nèi)切圓;擬合 </p><p>　　Abstract:Firstly，taking use of searching a

2、lgorithm to calcu-late coordinate of point of intersection for kinds of blood vessel’s axis demarcation line and radius of incircle.And then，averaging the radius of incircle which was calculated by a piece of dicing.It c

3、aught pipline radius of blood vessel by simulating all points of intersection. </p><p><b>　　0 引言 </b></p><p>　　在醫(yī)學(xué)、天文觀測、工業(yè)非破壞性試驗(yàn)等一些工業(yè)問題中，需要確定某個(gè)空間物體的形狀，但由于技術(shù)上的限制，無法將物體分離出來并直觀的顯現(xiàn).解決此

4、類問題的一個(gè)可行的方法是用等間隔的平行平面去截取這個(gè)物體，得到一組平行截面，通過采樣得到平行截面的數(shù)字圖像，進(jìn)行數(shù)據(jù)處理后可運(yùn)用計(jì)算機(jī)重建物體的三維形態(tài).在2001年全國大學(xué)生數(shù)學(xué)建模競賽題中，A題給出了某血管管道的相繼N張平行切片圖像，圖像格式、尺寸及坐標(biāo)均已給定，這些圖片記錄了血管管道與切片的交點(diǎn)，并且假定該血管管道是由一個(gè)小球沿某路徑移動(dòng)而形成的.取z軸垂直于切片，規(guī)定第一張切片平面為z=0，第N張切片平面為z=N-1.要求解決的

5、主要問題是:求解血管管道的中軸線與半徑，同時(shí)給出具體算法;以及繪制中軸線在XY，YZ，ZX平面的投影圖.現(xiàn)在我們試圖利用這組平行截面所提供的信息來重建被截血管管道. </p><p><b>　　1 一些假設(shè) </b></p><p>　　題中已給出的假設(shè):（1）管道的中軸線與每張切片有且只有一個(gè)交點(diǎn);（2）形成管道的球半徑固定;（3）切片間距以及圖像像素的尺寸均

6、為1. </p><p>　　求解所需的假設(shè):（1）假設(shè)球是沿光滑曲線運(yùn)動(dòng)的;（2）將像素抽象成為一個(gè)點(diǎn)，由于像素寬度與切片間距相等，所以由原數(shù)據(jù)得到的是一個(gè)均勻的三維點(diǎn)陣.（3）像素點(diǎn)的坐標(biāo)位置在像素點(diǎn)的中心. </p><p><b>　　2 問題分析 </b></p><p>　　此類問題的關(guān)鍵是:如何利用所已知的一組等間距的平行截面

7、，根據(jù)其平行截面的形狀確定球心軌跡（此軌由于采用本文所用的方法，參與“2001年全國大學(xué)生數(shù)學(xué)建模競賽”的建模小組獲得了“國家一等獎(jiǎng)”. 跡是一個(gè)工程上可接受其精度的近似值）.由于問題中假設(shè)管道中軸線與每張切片有且只有一個(gè)交點(diǎn)，則切片圖像在XZ，YZ平面的投影是以z為自變量的函數(shù).可以通過求出單個(gè)截面與中軸線的交點(diǎn)，再將所有截面與中軸線的交點(diǎn)進(jìn)行曲線擬合得到所要求的結(jié)果.具體分析思路如下. </p><p>　

8、　2.1 對單個(gè)截面的分析 文中所給的單個(gè)截面，我們可以認(rèn)為是形成管道的球體沿某一空間曲線穿過一截平面，球體被截平面連續(xù)截取所得的圓在該截平面上形成的包絡(luò)即為截面.由截面形成過程可知，截平面所截得的包含球心的圓就是截面的最大內(nèi)切圓，其圓心坐標(biāo)為截面與中軸線的交點(diǎn)坐標(biāo).由單個(gè)截面所給的信息，我們便可確定每個(gè)切片與管道中軸線的交點(diǎn)，并通過求解其最大內(nèi)切圓的平均半徑得出球體半徑.形成管道的平移球體穿過截面時(shí)，所形成的相交圖形是一系列半徑可變

9、的圓，可設(shè)其半徑函數(shù)為r=r（t），t1 </p><p>　　2.2 對所有截面的分析 單個(gè)截面所提供的是局部坐標(biāo)所包含的內(nèi)容，按上述同樣思路求得所有切片與中軸線的交點(diǎn).所有截面的信息則包含了連接各切片的中軸線的內(nèi)容，將所有的交點(diǎn)擬合后便可得到中軸線的近似軌跡. </p><p><b>　　3 問題求解 </b></p><p

10、>　　基于上述分析，對于一個(gè)給定的截面形狀，就可以得到這個(gè)截面與管道中軸線的唯一交點(diǎn)，若求出每個(gè)截面與中軸線的交點(diǎn)，則通過空間曲線的擬合就可以得到原管道的中軸線. </p><p>　　3.1 求解流程 （1）將截面的BMP圖像中的數(shù)據(jù)提取出來并保存在一個(gè)二維數(shù)組dot［x］［y］中;（2）對dot［x］［y］進(jìn)行邊界點(diǎn)識(shí)別，除去原點(diǎn)陣dot［x］［y］中的非邊界點(diǎn)，形成一個(gè)只包含邊界點(diǎn)的二維數(shù)組Bord

11、er［x］［y］;（3）搜索邊界內(nèi)的點(diǎn)，找出最大內(nèi)切圓圓心，此點(diǎn)即為管道中軸線與截面的交點(diǎn)Ci;（4）對所有截面進(jìn)行同樣處理，求出所有的Ci （i=0，1，2……N-1）;（5）將這N個(gè)空間點(diǎn)分別投影到XZ，YZ，XY平面上并進(jìn)行擬合;（6）將擬合后的曲線方程兩兩聯(lián)立，即得到三個(gè)所求中軸線的空間曲線方程;（7）再將這三個(gè)曲線平均，得到最終的計(jì)算結(jié)果 三維的中軸線. </p><p>　　3.2 搜索最大內(nèi)切圓圓心

12、Ci 由切片及截面形成的過程可知，在截面上必定存在一個(gè)最大內(nèi)切圓，中軸線與截面的交點(diǎn)即為該內(nèi)切圓的圓心. </p><p>　　算法描述:（1）將邊界點(diǎn)的坐標(biāo)按順時(shí)針（逆時(shí)針也可）的方向依次寫入一個(gè)坐標(biāo)序列（Xi ，Yi ）（i=0，1，2……N-1，N為邊界點(diǎn)總數(shù)）中;（2）對任意一個(gè)截面內(nèi)的點(diǎn)（即黑色的點(diǎn)），記錄其到邊界的最短距離dmin p ;（3）掃描所有截面內(nèi)的點(diǎn)，對于每一個(gè)點(diǎn)按上述方法求出到邊界的

13、最短距離，找出這些最短距離中最大的一個(gè)所對應(yīng)的點(diǎn)，則此點(diǎn)即為最大的內(nèi)切圓圓心，亦即管道中軸線與此截面的交點(diǎn)Ci . </p><p>　　為了進(jìn)一步精確求得圓心位置，我們進(jìn)行了以下修正:將所取的參考點(diǎn)在其周圍的8個(gè)方向上都減小1/2步長，再以上述同樣的方法作圓，這樣所得的半徑最大的圓就與最大內(nèi)切圓更加接近，所確定的該參考點(diǎn)即為圓心. </p><p>　　3.3 擬合中軸線 通過MATL

14、AB軟件中相應(yīng)功能，應(yīng)用最小二乘擬合將所有切片與中軸線交點(diǎn)在空間進(jìn)行擬合.具體擬合方法如下:對于上述方法求出的C（xi ，yi ，zi ），i=0，1，2，…N-1，將其投影到各個(gè)坐標(biāo)平面上，即得到XY平面上的（xi ，yi ），YZ平面上的（yi ，zi ），XZ平面上的（xi ，zi ），應(yīng)用最小二乘法將二維點(diǎn)用一個(gè)多項(xiàng)式擬合，應(yīng)用MATLAB的PULYFIT函數(shù)選取適當(dāng)?shù)臄M合次數(shù)進(jìn)行求解，得到中軸線在三個(gè)坐標(biāo)面上的投影曲線的方程［

15、2］ . </p><p>　　使用MATLAB軟件畫出兩兩聯(lián)立的空間曲線，選擇最優(yōu)的一條，得到最終的空間曲線為:XZ平面與YZ平面的投影擬合方程聯(lián)立.y=0.4371-0.7134z+0.1308z2 -0.0082z3 +0.0003z4x=-160.4458-2.1263z+0.12019z2 -0.2649z3 +0.0301z4 -0.002z5 +0.0001z6具體圖像見圖1. </p>

16、<p>　　圖1 血管管道中軸線　略 </p><p>　　4 算法改進(jìn)及推廣 </p><p>　　4.1 算法改進(jìn) 在上述問題中尋找最大內(nèi)切圓時(shí)，程序采用的是對邊界內(nèi)的所有點(diǎn)進(jìn)行逐個(gè)搜索求解，這樣處理程序效率不高，為了減少在尋找最大內(nèi)切圓時(shí)所搜索的點(diǎn)數(shù)，設(shè)各個(gè)邊界點(diǎn)均有按序的標(biāo)號，我們將具體的算法改進(jìn)如下:首先，找出內(nèi)切圓圓心.在所求切片截面中選取一內(nèi)部點(diǎn)作為基準(zhǔn)點(diǎn)，由

17、搜索法求得邊界上所有點(diǎn)到此點(diǎn)的最短距離r1 .若有 r2 -r1 &lt;ε成立（r2 為其他邊界點(diǎn)到該基準(zhǔn)點(diǎn)的距離，ε為任意小的正數(shù)），再判斷滿足上述條件的點(diǎn)中，是否有標(biāo)號相差較大（大于5）的兩點(diǎn).若存在至少兩點(diǎn)滿足上述條件，則可確定該基準(zhǔn)點(diǎn)為內(nèi)切圓圓心.其次，選定r1 在大于一個(gè)規(guī)定的值a的范圍內(nèi)搜索.由前邊分析可知，內(nèi)切圓半徑隨x，y變化，會(huì)在中間某處出現(xiàn)最大值，所以，若r </p><p>　　

18、1 a的方向移動(dòng)，可確定作為參考的內(nèi)切圓圓心.最后，在此參考圓心的周圍取一矩形區(qū)域，在區(qū)域內(nèi)搜索半徑更大的內(nèi)切圓圓心，并逐步迭代［3］ .這樣就將所需搜索的范圍大大壓縮.經(jīng)改進(jìn)后的算法在處理更復(fù)雜的數(shù)據(jù)群時(shí)，更能體現(xiàn)出優(yōu)越性. </p><p>　　4.2 算法的推廣 考慮到實(shí)際中生物組織、器官等的形態(tài)是一類特殊的管道，該管道的表面與我們上述討論有所不同，是由半徑在連續(xù)變化的球體其球心沿著某一曲線滾動(dòng)包絡(luò)而成.

19、由于管道半徑是連續(xù)變化的，所以在較小的變化范圍內(nèi)，可近似視為其球體半徑是固定的，按問題所提供的算法，再進(jìn)行相應(yīng)的處理. </p><p><b>　　參考文獻(xiàn): </b></p><p>　　［1］黃友謙，李岳生.數(shù)值逼近［M］.第2版.北京:高等教育出版社，1999:97-131. </p><p>　　［2］王沫然.MATLAB5.X與科