血管切片的三維重建

摘要：為了利用血管切片圖象重建血管的三維形態(tài)，我們首先編程對(duì)切片圖象進(jìn)行由BMP數(shù)據(jù)格式向文本格式的轉(zhuǎn)換，而這樣(zhèyàng)獲得的數(shù)據(jù)文件將較大，共約50M。顯然在利用計(jì)算機(jī)做進(jìn)一步讀寫與處理前，要著重面對(duì)的問題是對(duì)大量數(shù)據(jù)的處理，而其中卻有大量冗余信息存在，于是第一頁(yè)，共三十一頁(yè)。我們?cè)趯ふ颐繌堁芮衅行妮S點(diǎn)的過程中，利用多種優(yōu)化算法以簡(jiǎn)化問題，并確定出100個(gè)中心軸點(diǎn)。以此100個(gè)中心軸點(diǎn)為樣點(diǎn)，依靠樣條插值，利用Matlab軟件對(duì)分別對(duì)平面及空間曲線進(jìn)行插值，先后建立了模型一、模型二，從而擬合得到XY，YZ，ZX平面的投影(tóuyǐng)曲線以及中第二頁(yè)，共三十一頁(yè)。心軸線，得到血管半徑為29個(gè)像素單位，并模擬給出血管的三維空間形態(tài)(如以下圖)，進(jìn)而對(duì)結(jié)果進(jìn)行了廣泛的分析與評(píng)價(jià)。同通過矩陣投影到XY，YZ，ZX平面(píngmiàn)而獲得的圖象進(jìn)行比較，可以驗(yàn)證模型的正確性，以及模型的先進(jìn)性。第三頁(yè)，共三十一頁(yè)。

血管(xuèguǎn)切片的三維重建圖第四頁(yè)，共三十一頁(yè)。一、問題重述

斷面可用于了解生物組織、器官等的形態(tài)。例如，將樣本染色后切成厚約1mm的切片，在顯微鏡下觀察該橫斷面的組織形態(tài)結(jié)構(gòu)(jiégòu)。如果用切片機(jī)連續(xù)不斷地將樣本切成數(shù)十、成百的平行切片，可依次逐片觀察。根據(jù)拍照并采樣得到的平行切片數(shù)字圖象，運(yùn)用計(jì)算機(jī)可重建組織、器官等準(zhǔn)確的三維形態(tài)。第五頁(yè)，共三十一頁(yè)。假設(shè)某些血管可視為一類(yīlèi)特殊的管道，該管道的外表是由球心沿著某一曲線〔稱為中軸線〕的球滾動(dòng)包絡(luò)而成。例如圓柱就是這樣一種管道，其中軸線為直線，由半徑固定的球滾動(dòng)包絡(luò)形成。

現(xiàn)有某管道的相繼100張平行切片圖象，記錄了管道與切片的交。圖象文第六頁(yè)，共三十一頁(yè)。件名依次為0.bmp、1.bmp、…、99.bmp，格式均為BMP，寬、高均為512個(gè)象素〔pixel〕。為簡(jiǎn)化起見，假設(shè)：管道中軸線與每張切片有且只有一個(gè)交點(diǎn)；球半徑固定；切片間距以及(yǐjí)圖象象素的尺寸均為1。

取坐標(biāo)系的Z軸垂直于切片，第1張切片為平面Z=0，第100張切片為平面第七頁(yè)，共三十一頁(yè)。=99。Z=z切片圖象中象素的坐標(biāo)依它們?cè)谖募谐霈F(xiàn)的前后次序?yàn)?/p>

〔-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ù)以上所給信息(xìnxī),如何重建血管的三維形態(tài),是一個(gè)重要而且實(shí)用的問題。解決方案如下:第八頁(yè)，共三十一頁(yè)。二、問題分析

我們主要求解的是血管管道的中軸線與半徑。管道可以近似地看作(kànzuò)是一個(gè)半徑固定的球體滾動(dòng)而成的，中軸線是球心滑過的曲線，是連續(xù)的。我們等距平行切割血管，中軸線與每張切片有且僅有一個(gè)交點(diǎn)，也就是每張切片上有且僅有一個(gè)球心，那么在每張切片上總可以找到且只能找第九頁(yè)，共三十一頁(yè)。到一個(gè)以球心為圓心，球半徑為半徑的圓，而且是此切片(qiēpiàn)的最大內(nèi)切圓，反過來也是成立的。因此，我們只需找到每張切片(qiēpiàn)中的球心坐標(biāo)就可以用樣條插值得到中軸線，通過尋找最大內(nèi)切圓得到半徑，而中軸線在XY，YZ，ZX平面的投影圖只需令Z=0，X=0，Y=0就可以得到。第十頁(yè)，共三十一頁(yè)。三、模型假設(shè)

1、假設(shè)血管管道的外表是由球心沿著某一曲線〔中軸線〕的球滾動(dòng)而成的，也就是血管管道半徑〔即球體半徑〕固定且中軸線是連續(xù)的。

2、我們很合理的假設(shè)所有數(shù)據(jù)均是準(zhǔn)確(zhǔnquè)的。用大量象素能夠非常近似地描繪一個(gè)圖形。第十一頁(yè)，共三十一頁(yè)。3、因血管管道可以看作是一個(gè)半徑固定的球體滾動(dòng)包絡(luò)而成的，因此我們認(rèn)為中軸線，中軸線的一階導(dǎo)數(shù)及二階導(dǎo)數(shù)都是連續(xù)的。

4、數(shù)據(jù)精確(jīngquè)到單位像素。因切片厚度為1uM，因此像素單位亦為１uM，對(duì)于５１２＊５１２的信息距陣來說，此精確(jīngquè)程度已足以解釋問題。

5、切片與血管中軸線的交點(diǎn)存在且唯一。第十二頁(yè)，共三十一頁(yè)。四、符號(hào)說明(shuōmíng)

r：表示血管管道半徑。

ri：表示第i張切片中計(jì)算得到的管道半徑。

ro：表示r的初始估計(jì)值〔ro<r〕。

Sjk：表示第i張切片中第j個(gè)內(nèi)點(diǎn)到第k個(gè)邊界點(diǎn)的距離。

Sj：表示第i張切片中第j個(gè)內(nèi)點(diǎn)到邊界的最短距離。

(xi,yi,zi)：表示第i張切片得到的最大內(nèi)切圓的圓心坐標(biāo)。第十三頁(yè)，共三十一頁(yè)。五、模型構(gòu)成

首先我們用第一個(gè)C程序把BMP圖象文件格式轉(zhuǎn)換為保存了圖象點(diǎn)陣信息的文本文件。每個(gè)BMP文件對(duì)應(yīng)一個(gè)文本文件〔從00.txt到99.txt〕，每個(gè)文本文件均為512行，512列，按圖象的視覺直觀順序，以0表示白色象素點(diǎn)，1表示黑色象素點(diǎn)(如右圖)，以此方便后續(xù)程序的讀寫及人為(rénwéi)地糾錯(cuò)。第十四頁(yè)，共三十一頁(yè)。接著，第二個(gè)程序?qū)崿F(xiàn)了各個(gè)球心坐標(biāo)求取。確定每張切片與中軸線的交點(diǎn)的坐標(biāo)及半徑：按理論來說，血管截面邊界是圓滑曲線(qūxiàn)，對(duì)一個(gè)球體過球心切割，無論怎樣切，都必得到一個(gè)大圓，由假設(shè)5，每個(gè)切片包含且僅包含一個(gè)這樣的大圓，即它的最大內(nèi)切圓，圓心即球心，圓半徑即球半徑。下面我們確定每張切片中最大內(nèi)切第十五頁(yè)，共三十一頁(yè)。圓的圓心。對(duì)第i張切片上的第j個(gè)內(nèi)點(diǎn)，求到邊界點(diǎn)k的距離Sjk,從中選取一個(gè)最小距離Sj,=min{Sjk}再?gòu)闹羞x取一個(gè)最大的,記為ri=max{Sj}.不難理解，這就是(jiùshì)第i張切片的最大內(nèi)切圓的半徑。相應(yīng)的內(nèi)點(diǎn)即為圓心〔中軸線與截面的交點(diǎn)〕。對(duì)100張切片搜索后就得到100個(gè)球心點(diǎn)。但我們要面臨的問題是切片是由大量的象素點(diǎn)近似描繪出來的，會(huì)第十六頁(yè)，共三十一頁(yè)。給我們的計(jì)算帶來一定的誤差，但是由于象素很小，誤差就不會(huì)很大，我們依然可以應(yīng)用上述理論，人為的進(jìn)行誤差分析和修改，更準(zhǔn)確地得到球心(xyz)及由假設(shè)4可知半徑r=29個(gè)像素(xiànɡsù)單位〔29=<ri<=29.69〕。第十七頁(yè)，共三十一頁(yè)。XYZXYZXYZXYZXYZ(-16100)(-16101)(-16102)(-16103)(-16104)

(-16105)(-16106)(-16117)(-16118)(-16119)

(-161210)(-161211)(-161212)(-161413)(-161514)

(-161615)(-161816)(-1611017)(-1611318)(-1611719)

(-1611820)(-1611921)(-1612022)(-1612023)(-1612024)

(-1612025)(-1612026)(-1602927)(-1603028)(-1593529)

(-1593530)(-1593531)(-1584032)(-1574433)(-1564834)

(-1555135)(-1564836)(-1564837)(-1526038)(-1506539)

(-1506540)(-1388841)(-1369142)(-1369143)(-1369144)

(-1369145)(-1369146)(-11911247)(-11811348)(-11711449)

(-11611550)(-11511651)(-11411752)(-11311853)(-11211954)(-10412655)(-9613256)(-7114757)(-7114758)(-6015259)

(-6015260)(-4615761)(-2016362)(-2016363)(-2016364)

(-1316465)(-1316466)(-1316467)(3816368)(4316269)

(4816170)(5316071)(6015872)(6015873)(6715674)

(7515375)(8015176)(8714877)(8714878)(11613179)

(11912980)(13111981)(13111982)(13211883)(14410684)

(14410685)(14510586)(14510587)(1519888)(1519889)

(1638190)(1667691)(1755892)(1755893)(1775394)

(1804595)(1804596)(1814297)(1833598)(1843199)

在此之后我們找到兩種求中軸線的方案(fāngàn)，并對(duì)其進(jìn)行分析第十八頁(yè)，共三十一頁(yè)。模型1：

把從切片中得到的100個(gè)球心分別投影(tóuyǐng)到XY,XZ,YZ平面，直接應(yīng)用三次樣條插值中Ｍ連續(xù)性方程如下：第十九頁(yè)，共三十一頁(yè)。由假設(shè)3知一階導(dǎo)數(shù)連續(xù)，又可得到n-1個(gè)等式，那么可求解。

相應(yīng)地求出三條三次樣條插值曲線Y(z),Y(x),Z(x)沿垂直于它們所在平面的方向擴(kuò)展到空間，將得到三個(gè)曲面，它們兩兩相交得到三條曲線。如果三條曲線能夠重合，那無疑就是我們要求的中軸線，但是由于曲線Y(z),Y(x),Z(x)都是三次樣條插值擬合出的曲線。在插值過程(guòchéng)中，由于分別利用空間點(diǎn)的投影數(shù)據(jù)，因此插值結(jié)果將分別喪失一維信息，這樣由第二十頁(yè)，共三十一頁(yè)。投影點(diǎn)擬合的曲線，將極有可能不能復(fù)原成三維結(jié)果，即使復(fù)原也不會(huì)完全通過100個(gè)先前得到樣本點(diǎn)，因此誤差將會(huì)較大。如果能找到一個(gè)準(zhǔn)那么將三條曲線合并成一條，問題就解決了，但是三條曲線中的任意一條都滿足過樣本點(diǎn)且處處連續(xù)，我們已有的數(shù)據(jù)只有樣本點(diǎn)，所以找不到這樣的一個(gè)準(zhǔn)那么。那么我們只能任意選取兩條曲線沿它們所在平面的垂直方向(fāngxiàng)擴(kuò)展得到的曲面相交得到的曲線作為中軸線。這樣我們可以解決重建血管的三維形態(tài)的問題，但是它存在著很大弊端，由上面的分析我第二十一頁(yè)，共三十一頁(yè)。們知道，因?yàn)?yīnwèi)中軸線是由兩條三次樣條插值曲線分別沿著它們所在平面的垂直方向擴(kuò)展的兩個(gè)曲面相交而得到的，而不能由那三個(gè)曲面的兩兩相交得到的三條曲線重合得到，也就不可防止地存在較大誤差，這種方案也就存在著一定的局限性。因此，在此根底上我們提出第二種方案。

模型2：當(dāng)100個(gè)中心軸樣點(diǎn)坐標(biāo)后，考慮采用樣條插值將樣點(diǎn)光滑連接，但實(shí)際問題是一空間曲線的插值問題，因此與一般的一維插值有技術(shù)上的差異。我們將Z軸坐標(biāo)z視為x及y坐標(biāo)的參數(shù)，當(dāng)z做單調(diào)變化時(shí),那么x(z),y(z)可分別看作相應(yīng)此參量的橫、縱坐標(biāo)方向的變化量，第二十二頁(yè)，共三十一頁(yè)。這樣我們就可以利用一維樣條插值思想將空間曲線模擬出來。當(dāng)中軸線確定后，再分別令x=0，y=0，z=0，便可以得到中軸線在三個(gè)坐標(biāo)平面上的投影，(見以下圖)而且還可在中軸線上取充分多點(diǎn)，利用Matlab畫出半徑為29的球，由此而形成的球體包絡(luò)線便可展現(xiàn)三維血管(xuèguǎn)的空間形態(tài)，且能得到分別在三個(gè)坐標(biāo)面上的投影?！惨姾髨D〕第二十三頁(yè)，共三十一頁(yè)。第二十四頁(yè)，共三十一頁(yè)。六、算法分析

1、估計(jì)圓心

我們可以用計(jì)算機(jī)搜索第i張切片(qiēpiàn)半徑ri=maxj{mink{Sjk}}，記錄下該內(nèi)點(diǎn)就是此切片(qiēpiàn)中最大內(nèi)切圓的圓心，即中軸線上的一點(diǎn)，最小距離中的最大值就是據(jù)此切片(qiēpiàn)得到的管道半徑。但是這樣運(yùn)算量非常大，所以我們進(jìn)行如下優(yōu)化。第二十五頁(yè)，共三十一頁(yè)。首先我們對(duì)第一個(gè)切片圖象進(jìn)行分析，找到半徑一個(gè)初始值。第一個(gè)切片圖象可以近似地看作一個(gè)圓，在截面的邊界點(diǎn)中選取4個(gè)不要挨的太近且不在兩條平行線上的點(diǎn)，每?jī)牲c(diǎn)相連，得到兩條線段，分別作它們的垂直平分線必相交于一點(diǎn)O0，計(jì)算四個(gè)所選取的四個(gè)邊界點(diǎn)到這交點(diǎn)的距離，取一個(gè)最小的作為r0的值。因?yàn)閳A上任意一條弦的垂直平分線必過圓心，所以O(shè)0點(diǎn)可以近似的看作圓心，那么(nàme)管道的真實(shí)半徑一定大于r0，但不會(huì)相差很多。于是我們就可以在計(jì)算機(jī)搜索第一張切片時(shí)加上約束條件：如果內(nèi)點(diǎn)到邊界點(diǎn)的距離出現(xiàn)一個(gè)小于r0就將這個(gè)內(nèi)點(diǎn)舍第二十六頁(yè)，共三十一頁(yè)。去，搜索下一個(gè)；如果只有一個(gè)邊界點(diǎn)到內(nèi)點(diǎn)的距離到達(dá)最小值，也將這個(gè)內(nèi)點(diǎn)舍去。這樣我們(wǒmen)就大大地降低了計(jì)算量，通過對(duì)第一個(gè)切片的搜索可以得到該切片中的球心〔x1，y1，0〕及半徑r1,由此我們(wǒmen)可以適中選取新的更為接近管道真實(shí)半徑的r0值，對(duì)以后的99個(gè)切片做如上相同約束條件的搜索。記錄下每一個(gè)切片得到的圓心坐標(biāo)和半徑，對(duì)100個(gè)圓心進(jìn)行樣條插值，得到的插值函數(shù)就是我們(wǒmen)所要求的中軸線。第二十七頁(yè)，共三十一頁(yè)。七、結(jié)