血管切片的三維重建

1、摘　要：為了利用血管切片圖象重建血管的三維形態(tài)，我們首先編程對切片圖象進行由BMP數(shù)據(jù)格式向文本格式的轉(zhuǎn)換，而這樣獲得的數(shù)據(jù)文件將較大，共約50M。顯然在利用計算機做進一步讀寫與處理前，要著重面對的問題是對大量數(shù)據(jù)的處理，而其中卻有大量冗余信息存在，于是,我們在尋找每張血管切片中心軸點的過程中，利用多種優(yōu)化算法以簡化問題，并確定出100個中心軸點。以此100個中心軸點為樣點，依靠樣條插值，利用Matlab 軟件對分別對平面及空間曲線進行

2、插值，先后建立了模型一、模型二，從而擬合得到XY，YZ，ZX平面的投影曲線以及中,心軸線，得到血管半徑為29個像素單位，并模擬給出血管的三維空間形態(tài)(如下圖)，進而對結(jié)果進行了廣泛的分析與評價。同通過矩陣投影到XY，YZ，ZX平面而獲得的圖象進行比較，可以驗證模型的正確性，以及模型的先進性。,,血管切片的三維重建圖,一、問題重述斷面可用于了解生物組織、器官等的形態(tài)。例如，將樣本染色后切成厚約1m m的切片，在顯微鏡下觀察該橫斷面的組織

3、形態(tài)結(jié)構(gòu)。如果用切片機連續(xù)不斷地將樣本切成數(shù)十、成百的平行切片，可依次逐片觀察。根據(jù)拍照并采樣得到的平行切片數(shù)字圖象，運用計算機可重建組織、器官等準確的三維形態(tài)。,假設(shè)某些血管可視為一類特殊的管道，該管道的表面是由球心沿著某一曲線（稱為中軸線）的球滾動包絡(luò)而成。例如圓柱就是這樣一種管道，其中軸線為直線，由半徑固定的球滾動包絡(luò)形成?，F(xiàn)有某管道的相繼100張平行切片圖象，記錄了管道與切片的交。圖象文,件名依次為0.bmp、1.bmp、…、

4、 99.bmp，格式均為BMP，寬、高均為512個象素（pixel）。為簡化起見，假設(shè)：管道中軸線與每張切片有且只有一個交點；球半徑固定；切片間距以及圖象象素的尺寸均為1。取坐標系的Z軸垂直于切片，第1張切片為平面Z=0，第100張切片為平面,,=99。Z=z切片圖象中象素的坐標依它們在文件中出現(xiàn)的前后次序為（-256，-256，z），（-256，-255，z），…（-256，255，z），（-255，-256，z），（-255，

5、-255，z），…（-255，255，z），……（255，-256，z）,（ 255，-255，z），…（255，255，z）。 根據(jù)以上所給信息,如何重建血管的三維形態(tài),是一個重要而且實用的問題。解決方案如下:,二、問題分析 我們主要求解的是血管管道的中軸線與半徑。管道可以近似地看作是一個半徑固定的球體滾動而成的，中軸線是球心滑過的曲線，是連續(xù)的。我們等距平行切割血管，中軸線與每張切片有且僅有一個交點，也就是每張

6、切片上有且僅有一個球心，那么在每張切片上總可以找到且只能找,到一個以球心為圓心，球半徑為半徑的圓，而且是此切片的最大內(nèi)切圓，反過來也是成立的。因此，我們只需找到每張切片中的球心坐標就可以用樣條插值得到中軸線，通過尋找最大內(nèi)切圓得到半徑，而中軸線在XY，YZ，ZX平面的投影圖只需令Z=0，X=0，Y=0就可以得到。,三、模型假設(shè)1、假設(shè)血管管道的表面是由球心沿著某一曲線（中軸線）的球滾動而成的，也就是血管管道半徑（即球體半徑）固定且中軸

7、線是連續(xù)的。2、我們很合理的假設(shè)所有數(shù)據(jù)均是準確的。用大量象素能夠非常近似地描繪一個圖形。,3、因血管管道可以看作是一個半徑固定的球體滾動包絡(luò)而成的，因此我們認為中軸線，中軸線的一階導(dǎo)數(shù)及二階導(dǎo)數(shù)都是連續(xù)的。4、數(shù)據(jù)精確到單位像素。因切片厚度為1uM，因此像素單位亦為１uM，對于５１２＊５１２的信息距陣來說，此精確程度已足以解釋問題。5、切片與血管中軸線的交點存在且唯一。,四、符號說明r：表示血管管道半徑。ri：表示第i張切片

8、中計算得到的管道半徑。ro：表示r的初始估計值（ro<r）。Sjk：表示第 i張切片中第j個內(nèi)點到第k個邊界點的距離。Sj ：表示第i張切片中第j個內(nèi)點到邊界的最短距離。(xi,yi,zi)：表示第 i張切片得到的最大內(nèi)切圓的圓心坐標。,五、模型構(gòu)成 首先我們用第一個C程序把BMP圖象文件格式轉(zhuǎn)換為保存了圖象點陣信息的文本文件。每個BMP文件對應(yīng)一個文本文件（從00.txt到99.txt），每個文本文件均為512行

9、，512列，按圖象的視覺直觀順序，以0表示白色象素點， 1表示黑色象素點( 如右圖)，以此方便后續(xù)程序的讀寫及人為地糾錯。,接著，第二個程序?qū)崿F(xiàn)了各個球心坐標求取。確定每張切片與中軸線的交點的坐標及半徑：按理論來說，血管截面邊界是圓滑曲線，對一個球體過球心切割，無論怎樣切，都必得到一個大圓，由假設(shè)5，每個切片包含且僅包含一個這樣的大圓，即它的最大內(nèi)切圓，圓心即球心，圓半徑即球半徑。下面我們確定每張切片中最大內(nèi)切,圓的圓心。對第i 張切片

10、上的第j個內(nèi)點，求到邊界點 k的距離Sjk,從中選取一個最小距離Sj,=min{Sjk}再從中選取一個最大的,記為ri=max{Sj}.不難理解，這就是第i張切片的最大內(nèi)切圓的半徑。相應(yīng)的內(nèi)點即為圓心（中軸線與截面的交點）。對100張切片搜索后就得到100個球心點。但我們要面臨的問題是切片是由大量的象素點近似描繪出來的，會,給我們的計算帶來一定的誤差，但是由于象素很小，誤差就不會很大，我們依然可以應(yīng)用上述理論，人為的進行誤差分析和修改，

11、更準確地得到球心(x y z)及由假設(shè)4可知半徑r=29個像素單位（29=<ri<=29.69）。,X Y Z X Y Z X Y Z X Y Z X Y Z (-161 0 0) (-161 0 1) (-161 0 2) (-161 0 3) (-161 0 4)(-161 0 5)(-161 0 6)(-161 1 7)(-161 1 8)(-161 1 9)(-161 2 10)

12、(-161 2 11)(-161 2 12)(-161 4 13)(-161 5 14)(-161 6 15)(-161 8 16)(-161 10 17)(-161 13 18)(-161 17 19)(-161 18 20)(-161 19 21)(-161 20 22)(-161 20 23)(-161 20 24)(-161 20 25)(-161 20 26)(-160 29 27)(-160 30 28)(-159 3

13、5 29)(-159 35 30)(-159 35 31)(-158 40 32)(-157 44 33)(-156 48 34)(-155 51 35)(-156 48 36)(-156 48 37)(-152 60 38)(-150 65 39)(-150 65 40)(-138 88 41)(-136 91 42)(-136 91 43)(-136 91 44)(-136 91 45)(-136 91 46)(-119 1

14、12 47)(-118 113 48)(-117 114 49)(-116 115 50)(-115 116 51)(-114 117 52)(-113 118 53)(-112 119 54)(-104 126 55)(-96 132 56)(-71 147 57)(-71 147 58)(-60 152 59)(-60 152 60)(-46 157 61)(-20 163 62)(-20 163 63)(-20 163 64

15、)(-13 164 65)(-13 164 66)(-13 164 67)(38 163 68)(43 162 69)(48 161 70)(53 160 71)(60 158 72)(60 158 73)(67 156 74)(75 153 75)(80 151 76)(87 148 77)(87 148 78)(116 131 79)(119 129 80)(131 119 81)(131 119 82)(132 118 8

16、3)(144 106 84)(144 106 85)(145 105 86)(145 105 87)(151 98 88)(151 98 89)(163 81 90)(166 76 91)(175 58 92)(175 58 93)(177 53 94)(180 45 95)(180 45 96)(181 42 97)(183 35 98)(184 31 99)在此之后我們找到兩種求中軸線的方案，并對其進行分析,模型1：

17、 把從切片中得到的100個球心分別投影到XY,XZ,YZ平面，直接應(yīng)用三次樣條插值中Ｍ連續(xù)性方程如下：,由假設(shè)3知一階導(dǎo)數(shù)連續(xù)，又可得到n-1 個等式，則可求解。相應(yīng)地求出三條三次樣條插值曲線Y(z),Y(x),Z(x)沿垂直于它們所在平面的方向擴展到空間，將得到三個曲面，它們兩兩相交得到三條曲線。如果三條曲線能夠重合，那無疑就是我們要求的中軸線，但是由于曲線Y(z),Y(x),Z(x)都是三次樣條插值擬合出的曲線。在插值過程中，由

18、于分別利用空間點的投影數(shù)據(jù)，因此插值結(jié)果將分別喪失一維信息，這樣由,投影點擬合的曲線，將極有可能不能還原成三維結(jié)果，即使還原也不會完全通過100個先前得到樣本點，因此誤差將會較大。如果能找到一個準則將三條曲線合并成一條，問題就解決了，但是三條曲線中的任意一條都滿足過樣本點且處處連續(xù)，我們已有的數(shù)據(jù)只有樣本點，所以找不到這樣的一個準則。那么我們只能任意選取兩條曲線沿它們所在平面的垂直方向擴展得到的曲面相交得到的曲線作為中軸線。這樣我們可以

19、解決重建血管的三維形態(tài)的問題，但是它存在著很大弊端，由上面的分析我,們知道，因為中軸線是由兩條三次樣條插值曲線分別沿著它們所在平面的垂直方向擴展的兩個曲面相交而得到的，而不能由那三個曲面的兩兩相交得到的三條曲線重合得到，也就不可避免地存在較大誤差，這種方案也就存在著一定的局限性。因此，在此基礎(chǔ)上我們提出第二種方案。模型2：當100個中心軸樣點坐標已知后，考慮采用樣條插值將已知樣點光滑連接，但實際問題是一空間曲線的插值問題，因此與一般的

20、一維插值有技術(shù)上的差別。我們將Z軸坐標z視為x及y坐標的參數(shù)，當z做單調(diào)變化時,則x(z),y(z)可分別看作相應(yīng)此參量的橫、縱坐標方向的變化量，,這樣我們就可以利用一維樣條插值思想將空間曲線模擬出來。當中軸線確定后，再分別令x=0，y=0，z=0，便可以得到中軸線在三個坐標平面上的投影，(見下圖)而且還可在中軸線上取充分多點，利用Matlab畫出半徑為29的球，由此而形成的球體包絡(luò)線便可展現(xiàn)三維血管的空間形態(tài)，且能得到分別在三個坐標面

21、上的投影。（見后圖）,六、算法分析1、估計圓心 我們可以用計算機搜索第i張切片半徑ri=maxj{mink{Sjk}}，記錄下該內(nèi)點就是此切片中最大內(nèi)切圓的圓心，即中軸線上的一點，最小距離中的最大值就是據(jù)此切片得到的管道半徑。但是這樣運算量非常大，所以我們進行如下優(yōu)化。,首先我們對第一個切片圖象進行分析，找到半徑一個初始值。第一個切片圖象可以近似地看作一個圓，在截面的邊界點中選取4個不要挨的太近且不在兩條平行線上的點，每兩點

22、相連，得到兩條線段，分別作它們的垂直平分線必相交于一點O0，計算四個所選取的四個邊界點到這交點的距離，取一個最小的作為r0的值。因為圓上任意一條弦的垂直平分線必過圓心，所以O(shè)0點可以近似的看作圓心，那么管道的真實半徑一定大于r0，但不會相差很多。于是我們就可以在計算機搜索第一張切片時加上約束條件：如果內(nèi)點到邊界點的距離出現(xiàn)一個小于r0就將這個內(nèi)點舍,去，搜索下一個；如果只有一個邊界點到內(nèi)點的距離達到最小值，也將這個內(nèi)點舍去。這樣我們就大

23、大地降低了計算量，通過對第一個切片的搜索可以得到該切片中的球心（x1，y1，0）及半徑r1,由此我們可以適當選取新的更為接近管道真實半徑的r0值，對以后的99個切片做如上相同約束條件的搜索。記錄下每一個切片得到的圓心坐標和半徑，對100個圓心進行樣條插值，得到的插值函數(shù)就是我們所要求的中軸線。,七、結(jié)果模擬 利用距陣投影即把我們得到的100張圖象信息文本文件相應(yīng)象素進行并位運算,合并出100張圖的信息,然后按照bmp格式重新

24、存入一個文件,得到一張新的位圖(左圖),與我們得到的結(jié)果相比(右圖),發(fā)現(xiàn)具有驚人的一致性,從而驗證了我們的重建結(jié)果是可信的。,左圖：距陣XY投影疊加圖 右圖：三維重建XY投影圖,八、模型特點（1）算法合理易行，便于接受。（2）模型用到的數(shù)學(xué)知識是簡單易懂的。一階導(dǎo)數(shù)的差分公式及三次樣條函數(shù)的內(nèi)容現(xiàn)在已廣泛為人所知且在實際中用得很廣。（3）只要有一臺計算機以及連續(xù)系列BMP格式切片圖象，即可重建血管的三維形態(tài)。（4）模型結(jié)果直觀