醫(yī)學專題血管切片的三維重建

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

2、像素單位，并模擬給出血管的三維空間形態(tài)(如下圖)，進而對結(jié)果進行了廣泛的分析與評價。同通過矩陣投影到XY，YZ，ZX平面而獲得的圖象進行比較，可以驗證模型的正確性，以及模型的先進性。2021/7/20 星期二3血管切片的三維重建圖2021/7/20 星期二4一、問題重述斷面可用于了解生物組織、器官等的形態(tài)。例如，將樣本染色后切成厚約1m m的切片，在顯微鏡下觀察該橫斷面的組織形態(tài)結(jié)構。如果用切片機連續(xù)不斷地將樣本切成數(shù)十、成百的平行切片，可依次逐片觀察。根據(jù)拍照并采樣得到的平行切片數(shù)字圖象，運用計算機可重建組織、器官等準確的三維形態(tài)。2021/7/20 星期二5假設某些血管可視為一類特殊的管道

3、，該管道的表面是由球心沿著某一曲線（稱為中軸線）的球滾動包絡而成。例如圓柱就是這樣一種管道，其中軸線為直線，由半徑固定的球滾動包絡形成?，F(xiàn)有某管道的相繼100張平行切片圖象，記錄了管道與切片的交。圖象文2021/7/20 星期二6件名依次為0.bmp、1.bmp、 99.bmp，格式均為BMP，寬、高均為512個象素（pixel）。為簡化起見，假設：管道中軸線與每張切片有且只有一個交點；球半徑固定；切片間距以及圖象象素的尺寸均為1。取坐標系的Z軸垂直于切片，第1張切片為平面Z=0，第100張切片為平面2021/7/20 星期二7=99。Z=z切片圖象中象素的坐標依它們在文件中出現(xiàn)的前后次序為（

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

5、找到且只能找2021/7/20 星期二9到一個以球心為圓心，球半徑為半徑的圓，而且是此切片的最大內(nèi)切圓，反過來也是成立的。因此，我們只需找到每張切片中的球心坐標就可以用樣條插值得到中軸線，通過尋找最大內(nèi)切圓得到半徑，而中軸線在XY，YZ，ZX平面的投影圖只需令Z=0，X=0，Y=0就可以得到。2021/7/20 星期二10三、模型假設1、假設血管管道的表面是由球心沿著某一曲線（中軸線）的球滾動而成的，也就是血管管道半徑（即球體半徑）固定且中軸線是連續(xù)的。2、我們很合理的假設所有數(shù)據(jù)均是準確的。用大量象素能夠非常近似地描繪一個圖形。2021/7/20 星期二113、因血管管道可以看作是一個半徑固

6、定的球體滾動包絡而成的，因此我們認為中軸線，中軸線的一階導數(shù)及二階導數(shù)都是連續(xù)的。4、數(shù)據(jù)精確到單位像素。因切片厚度為1uM，因此像素單位亦為uM，對于的信息距陣來說，此精確程度已足以解釋問題。5、切片與血管中軸線的交點存在且唯一。2021/7/20 星期二12四、符號說明r：表示血管管道半徑。ri：表示第i張切片中計算得到的管道半徑。ro：表示r的初始估計值（ror）。Sjk：表示第 i張切片中第j個內(nèi)點到第k個邊界點的距離。Sj ：表示第i張切片中第j個內(nèi)點到邊界的最短距離。(xi,yi,zi)：表示第 i張切片得到的最大內(nèi)切圓的圓心坐標。2021/7/20 星期二13五、模型構成 首先我

7、們用第一個C程序把BMP圖象文件格式轉(zhuǎn)換為保存了圖象點陣信息的文本文件。每個BMP文件對應一個文本文件（從00.txt到99.txt），每個文本文件均為512行，512列，按圖象的視覺直觀順序，以0表示白色象素點， 1表示黑色象素點( 如右圖)，以此方便后續(xù)程序的讀寫及人為地糾錯。2021/7/20 星期二14接著，第二個程序?qū)崿F(xiàn)了各個球心坐標求取。確定每張切片與中軸線的交點的坐標及半徑：按理論來說，血管截面邊界是圓滑曲線，對一個球體過球心切割，無論怎樣切，都必得到一個大圓，由假設5，每個切片包含且僅包含一個這樣的大圓，即它的最大內(nèi)切圓，圓心即球心，圓半徑即球半徑。下面我們確定每張切片中最大內(nèi)

8、切2021/7/20 星期二15圓的圓心。對第i 張切片上的第j個內(nèi)點，求到邊界點 k的距離Sjk,從中選取一個最小距離Sj,=minSjk再從中選取一個最大的,記為ri=maxSj.不難理解，這就是第i張切片的最大內(nèi)切圓的半徑。相應的內(nèi)點即為圓心（中軸線與截面的交點）。對100張切片搜索后就得到100個球心點。但我們要面臨的問題是切片是由大量的象素點近似描繪出來的，會2021/7/20 星期二16給我們的計算帶來一定的誤差，但是由于象素很小，誤差就不會很大，我們依然可以應用上述理論，人為的進行誤差分析和修改，更準確地得到球心(x y z)及由假設4可知半徑r=29個像素單位（29=ri=29

9、.69）。2021/7/20 星期二17X 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)(-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

10、22)(-161 20 23)(-161 20 24)(-161 20 25)(-161 20 26)(-160 29 27)(-160 30 28)(-159 35 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 112

11、 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)(-13 164 65)(-13 164 66)(-13 164 67)(38 163 68)(43 162 69)(48 161 70)(53 160 71)(60

12、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 83)(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 9

13、8)(184 31 99)在此之后我們找到兩種求中軸線的方案，并對其進行分析2021/7/20 星期二18模型1： 把從切片中得到的100個球心分別投影到XY,XZ,YZ平面，直接應用三次樣條插值中連續(xù)性方程如下：2021/7/20 星期二19由假設3知一階導數(shù)連續(xù)，又可得到n-1 個等式，則可求解。相應地求出三條三次樣條插值曲線Y(z),Y(x),Z(x)沿垂直于它們所在平面的方向擴展到空間，將得到三個曲面，它們兩兩相交得到三條曲線。如果三條曲線能夠重合，那無疑就是我們要求的中軸線，但是由于曲線Y(z),Y(x),Z(x)都是三次樣條插值擬合出的曲線。在插值過程中，由于分別利用空間點的投影數(shù)

14、據(jù)，因此插值結(jié)果將分別喪失一維信息，這樣由2021/7/20 星期二20投影點擬合的曲線，將極有可能不能還原成三維結(jié)果，即使還原也不會完全通過100個先前得到樣本點，因此誤差將會較大。如果能找到一個準則將三條曲線合并成一條，問題就解決了，但是三條曲線中的任意一條都滿足過樣本點且處處連續(xù)，我們已有的數(shù)據(jù)只有樣本點，所以找不到這樣的一個準則。那么我們只能任意選取兩條曲線沿它們所在平面的垂直方向擴展得到的曲面相交得到的曲線作為中軸線。這樣我們可以解決重建血管的三維形態(tài)的問題，但是它存在著很大弊端，由上面的分析我2021/7/20 星期二21們知道，因為中軸線是由兩條三次樣條插值曲線分別沿著它們所在平

15、面的垂直方向擴展的兩個曲面相交而得到的，而不能由那三個曲面的兩兩相交得到的三條曲線重合得到，也就不可避免地存在較大誤差，這種方案也就存在著一定的局限性。因此，在此基礎上我們提出第二種方案。模型2：當100個中心軸樣點坐標已知后，考慮采用樣條插值將已知樣點光滑連接，但實際問題是一空間曲線的插值問題，因此與一般的一維插值有技術上的差別。我們將Z軸坐標z視為x及y坐標的參數(shù)，當z做單調(diào)變化時,則x(z),y(z)可分別看作相應此參量的橫、縱坐標方向的變化量，2021/7/20 星期二22這樣我們就可以利用一維樣條插值思想將空間曲線模擬出來。當中軸線確定后，再分別令x=0，y=0，z=0，便可以得到中

16、軸線在三個坐標平面上的投影，(見下圖)而且還可在中軸線上取充分多點，利用Matlab畫出半徑為29的球，由此而形成的球體包絡線便可展現(xiàn)三維血管的空間形態(tài)，且能得到分別在三個坐標面上的投影。（見后圖）2021/7/20 星期二232021/7/20 星期二24六、算法分析1、估計圓心 我們可以用計算機搜索第i張切片半徑ri=maxjminkSjk，記錄下該內(nèi)點就是此切片中最大內(nèi)切圓的圓心，即中軸線上的一點，最小距離中的最大值就是據(jù)此切片得到的管道半徑。但是這樣運算量非常大，所以我們進行如下優(yōu)化。2021/7/20 星期二25首先我們對第一個切片圖象進行分析，找到半徑一個初始值。第一個切片圖象可以

17、近似地看作一個圓，在截面的邊界點中選取4個不要挨的太近且不在兩條平行線上的點，每兩點相連，得到兩條線段，分別作它們的垂直平分線必相交于一點O0，計算四個所選取的四個邊界點到這交點的距離，取一個最小的作為r0的值。因為圓上任意一條弦的垂直平分線必過圓心，所以O0點可以近似的看作圓心，那么管道的真實半徑一定大于r0，但不會相差很多。于是我們就可以在計算機搜索第一張切片時加上約束條件：如果內(nèi)點到邊界點的距離出現(xiàn)一個小于r0就將這個內(nèi)點舍2021/7/20 星期二26去，搜索下一個；如果只有一個邊界點到內(nèi)點的距離達到最小值，也將這個內(nèi)點舍去。這樣我們就大大地降低了計算量，通過對第一個切片的搜索可以得到該切片中的球心（x1，y1，0）及半徑r1,由此我們可以適當選取新的更為接近管道真實半徑的r0值，對以后的99個切片做如上相同約束條件的搜索。記錄下每一個切片得到的圓心坐標和半徑，對100個圓心進行樣條插值，得到的插值函數(shù)就是我們所要求的中軸線。2021/7/20 星期二27七、結(jié)果模擬 利用距陣投影即把我們得到的100張圖象信息文本文件相應象素進行并位運算,合并出100張圖的信息,然后按照bmp格