三維血管重建優(yōu)秀論文_含源代碼

1、血管的三維重建摘要木文探討血管的三維重建，由血管的相繼100張平行切片圖像計算血管的 中軸線與半徑，并繪制血管在三個坐標平面上的投影。由于血管的表面是由球心沿著某一曲線（稱為屮軸線）的球滾動包絡(luò)而成, 由此我們得出結(jié)論：每個切片一定包含滾動球的大圓，并且它一定為切片的最大 內(nèi)切圓，而最人圓所對應(yīng)的半徑即為血管半徑，所以求血管半徑就轉(zhuǎn)化為求每一 個切片內(nèi)部的點到切片外部輪廓的最大半徑。首先，讀取100張血管切面圖，把它們轉(zhuǎn)換成logical矩陣，從中提取切片 截面輪廓點構(gòu)成一個新的矩陣。然后找到原圖片矩陣中像素點的內(nèi)點（切片圖片 屮輪廓線屮的點），從而得到內(nèi)點到切片輪廓點的最小距離矩陣和最小距離

2、屮的 最大值矩陣，最大值即為血管半徑。最后計算所有切片的血管半徑，并對這些半 徑求平均值，得到平均血管半徑為：29.6799u mo由100張切片的最大內(nèi)切圓 圓心坐標擬合得出中軸線方程以及其在三個坐標平面內(nèi)的投影曲線方程。由中軸線得到血管的三維立體重建圖,用平面z = i,（i = 0,24,49,74）去截血管 的三維立體重建圖，得到新的4張截面圖。把它們分別與題設(shè)中的對應(yīng)截面進行 內(nèi)點個數(shù)對比。我們泄義兩張切片所共同擁有的內(nèi)點個數(shù)與原切片內(nèi)點個數(shù)的比 值為重合度。計算得到平均重合度為：9&19%。關(guān)鍵詞：血管半徑屮軸線切片重建（來自作者：歡迎各界人士批評指正，學術(shù)交流郵箱nibz

3、。文章 作于2011年8月10日，陜西科技大學理學院實驗室）1問題的重述斷面可用于了解生物組織、器官等的形態(tài)。例如，將樣本染色后切成厚約 的切片，在顯微鏡下觀察該橫斷面的組織形態(tài)結(jié)構(gòu)。如果用切片機連續(xù)不 斷地將樣本切成數(shù)十、成百的平行切片，可依次逐片觀察。根據(jù)拍照并采樣得 到的平行切片數(shù)字圖象，運用計算機可重建組織、器官等準確的三維形態(tài)。假設(shè)某些血管可視為一類特殊的管道，該管道的表面是曲球心沿著某一曲 線(稱為中軸線)的球滾動包絡(luò)而成。例如圓柱就是這樣一種管道，其中軸線為 直線，由半徑固定的球滾動包絡(luò)形成?，F(xiàn)有某管道的相繼100張平行切片圖像，記錄了管道與切片的交。圖像文件 名依次為0.bmp

4、、1.bmp、99.bmp,格式均為bmp,寬、高均為512個 象素(pixel) o為簡化騷見，假設(shè)：管道中軸線與每張切片有且只有一個交點; 球半徑固定；切片間距以及圖像象素的尺寸均為1。計算管道的中軸線與半徑，給出具體的算法，并繪制中軸線在xy、yz、zx 平面的投影圖。2模型假設(shè)(1) 假設(shè)血管中軸線與每張切片有且只有一個交點。(2) 假設(shè)血管的半徑固定不變。(3) 假設(shè)在對切片拍照的過程中不存在誤差。(4) 假設(shè)血管不存在嚴重扭曲。3符號說明第£張切片中輪廓線的最大內(nèi)切圓圓心第r張切片屮輪廓線的最大內(nèi)切圓半徑d.jk第£張切片中第，個內(nèi)點與第丿個輪廓點的距離“512

5、x512x100100張切片圖片轉(zhuǎn)換以后的三維01矩陣刖5陀站00100張切片圖片的輪廓線生成的矩陣c"®00x3100張切片輪廓的最大內(nèi)切圓圓心坐標rf,平均血管半徑第£張切片圖片中所有內(nèi)點的集合第k張切片圖片屮輪廓點的集合str-square原切片圖片的上內(nèi)點及邊界點的集合 重新切片得到的內(nèi)點及邊界點的集合 多項式擬合的指標數(shù)4問題的分析根據(jù)題忖整個管道是由球心沿著某一曲線（稱為中軸線）的球滾動包絡(luò)而成。 基于幾何原理可以得出以下兩條定理。（1）球的任意截面都是圓。（2）經(jīng)過球心的球截面是所有的截面圓當中半徑最大的圓。基于上述兩個定理可以得出：每張切片的最大內(nèi)

6、切圓的圓心位于血管的屮軸 線上，該圓的半徑等于血管半徑。分別求出100張切片輪廓線的最人內(nèi)切圓圓心q與半徑他。 對所有最大內(nèi)切圓的圓心做擬合，即可求出血管的中軸線。100為減少謀差，我們用求均值的方法r=耳來得到平均血管半徑r。 100根據(jù)擬合出的中軸線方程，即可求出中軸線在xy、yz、zx平而上的投影。5模型的建立與求解5.1模型建立5.1.1求血管的半徑（1）導入數(shù)據(jù)，轉(zhuǎn)換存儲方式為了計算方便，利用matlab軟件將100張bmp文件轉(zhuǎn)換存儲為一個三維矩陣pa5i2x512x1oo （0代表黑色像素點，1代表白色像素點）。（2）求切片輪廓線上各點坐標利用matlab軟件的內(nèi)部函數(shù)edge,

7、求得所有切片圖片的輪廓線生成的矩陣pb512x512x100（3）求切片輪廓線的最大內(nèi)切圓的半徑首先求出每個內(nèi)點距輪廓線的距離，取其屮的最小值，即為以這些內(nèi)點為圓 心的輪廓線的內(nèi)切圓半徑；其次找出所有內(nèi)點確定的最小內(nèi)切圓中半徑最大的一 個，即為輪廓線的最人內(nèi)切圓半徑。具體步驟：stepl在三維01矩陣pa5i2x5i2xi00'i'遍歷搜索內(nèi)點，將其存儲于集合5）屮;step2在輪廓線矩陣pb512x5i2x1（）（）中遍歷搜索輪廓點，將其存儲于集合屮中； step3計算集合幾）中第i個內(nèi)點a（x,y）到集合屮屮第丿個輪廓點b（p,q）的距離 dijk = 7（x-p）2+g-

8、?）2，首先取每個內(nèi)點到輪廓點的最小距離，即為以這些內(nèi)點為圓心的內(nèi)切圓 半徑；其次找出所有內(nèi)點確定的最小內(nèi)切圓中半徑最大的一個，即為輪 廓線的最大內(nèi)切圓半徑rk = maxmind做i j jstep4計算平均血管半徑100斤=100（4）求切片輪廓線的最大內(nèi)切圓心坐標將100張輪廓線的最大內(nèi)切圓半徑所對應(yīng)的圓心坐標c,（z,w）記錄于矩陣100.3）中 °een（u）=乙< cen（k 2） = w，伙=1,2,100）c叫.3）二 k,5. 1. 2由圓心坐標求中軸線方程及曲線投影（1）利用matlab畫圖工具由100組圓心坐標畫出中軸線以及中軸線在xy、yz、zx面的投影

9、，由于圖像文件以象素為單位存儲數(shù)據(jù)，這使得圖像的原始數(shù)據(jù)是離散化的。對于離散的位圖來說，將變得不可分辨，可能切片輪廓線的最大 圓不止一個。為了減少誤差，將中軸線在xy、yz、zx面的曲線方程進行多項式 擬合，建立擬合方程：y 話（x）x = “(z)（2）由建立的投影曲面得到中軸線的方程<r = f,(t)；z = f2(ry95.2模型求解根據(jù)所建模型求岀100張切片的最大內(nèi)切圓的圓心坐標、半徑以及血管半徑表1最人內(nèi)切圓半徑及圓心處標切片序號最大內(nèi)切圓半徑（r/u m）最大內(nèi)切圓圓心處標橫坐標(x/ u m)縱坐標(y/ u m)豎坐標(z/u m)0290

10、69-16011229-16022329069-16024529155-16026729.155-16027829.275-16028929.275-160391029.428-1603101129.428-1604111229.614-1604121329.614-1605131429.833-1606141529.833-1607151629.833-1608161730-1609171830-16010181930-16011192030-16012202130-16013212230-16014222330-16015232430

11、-16015242529.833-1601625262961461461461461461461461473215355353629732152583738296145475475353

12、682682682698-115116474829.698-113118484929.698-112119495029.698111120505129.698-111120515229.698-111120525329.698-111120535429.682-111120545529.53-66150555629.53-59153565729.732-70148575829.967-70148585929.732-54155596029.732-54155606129.614-311626

13、16229.614-31162626329.6146166636429.833-6166646529.83381676566301016766673()11167676830121676869301216769703012167707129.83326166717229.61426166727329.61446163737429.73265158747529.73265158757629.73271156767729.68296145777829.68296145787929.698127124798029.698127124808129.732122128818229.73212212882

14、8329.732122128838429.698127124848529.732150101858629.732150101868729.73215496878829.698137115888929.698139113899029.68216088909129.68216285919229.54717659929329.61418340939429.73218049949529.61418340959629.61418340969729.61418533979830190098993()190199(1) 血管半徑為斤=29.6799 um。(2) 中軸線在xy平面投影曲線的擬合方程(rsqu

15、are:0.9939)r = 6.23xlo',8x9 -1.58xloj4x8 +1.74xlo-,1x7 -1.09xl0 8x6 +4.26xlo'6x5 -0.00108x4 +0.17693x3 18.059x2 +1041.7x 25566(3) 屮軸線在yz平面投影曲線的擬合方程(r-square:1)z = -6.32x1 o'15 y9 +1.95x10-'1y8 -2.66x io 8 y7 + 2.11 x io 中軸線方程：x =6.23x10 ,8x9 -1.58x10',4x8 -i-lxlo x7 -1.09x108x6 +

16、4.26x10'6x5 -0.00108x4 +0.17693x3 -18.059x2 +1041.7x-25566,z = -6.32 x 1015 y9 +1.95 x 10 11y8 - 2.66 x 10 8 y7 + 2.11 x 1o*5 y6 - 0.01075y5 + 3.6293y4 - 813.88y3+ 17 x 105 y2-9.75 xlo6y +3.60x10%x =1.29x1013z9 -5.96x 1 o'11z8 +1.16x 108z7 -1.23x 1 o'6z6+7.56x1 o'5z5 -0.00269z4 + 0.0

17、53686z3-0.55402z2 + 2.4609z + 92.988, y6 -0.01075y5 + 3.6293y4 -813.88y3 + l17x105y2 -9.75x106y + 3.60x108(4) 屮軸線在zx平面投影曲線的擬合方程(r-square:1)x = 1.29xlo ,3z9 -5.96xlo-"z8 +1.16xlo'8z7 -1.23xlo'6z6 +7.56x105z5 -0.00269z4 + 0.053686z3-0.55402z2 + 2.4609z+ 92.988圖1中軸線曲線圖2中軸線擬合曲線圖3屮軸線在xy平而的投影

18、圖4屮軸在xy平面的投影擬合1111線圖5中軸線在yz平面的投影圖6中軸線在yz平面的投影擬合曲線x圖7中軸線在zx平血的投影圖8中軸線在zx平面的投影擬合曲線曲線模擬的血管250、20015010050、ox50、-100-150圖9曲線模擬的血管6模型檢驗與分析用多項式擬合得到重建后血管的三維空i'可曲線，再用平而z = /,(/ = 1,25,50,75)去截空間曲線得到新的4層切片平而分別與題目中對應(yīng)切片進行對比，在這里我們稱sr>t的面積與宇的面積比值為切片重合度。第01、25、50、75張圖片分別與對應(yīng)的新截取圖片的對比圖圖1第01張切片 圖11第25張切片 圖12第

19、50張切片 圖13第75張切片圖圖14截后第01張切片圖15截后第25張切片圖16截后第50張切片圖17截后第75張切片根據(jù)上面曲線擬合得到的血管重新切片后，我們隨機取出四張重切的圖片, 與原始的切片比校，得到兩切片的重合度。計算得：第01張切片的重合度為99.20% ,第25張切片的重合度為98.65% ,第50張切片的重合度為97.30% ,第75張切片的重合度為97.61% ,計算平均重合度得到98.19%,說明木模型較為可靠。7模型的優(yōu)點與缺點優(yōu)點在于本模型的算法采用遍歷搜索的方法得到最人內(nèi)切圓的半徑的精度 較高，但是由于遍歷搜索的運算次數(shù)界常龐大，導致模型的缺點是實踐難度較人。 從模

20、型檢驗的數(shù)據(jù)可以看出計算的精度，但是整個運算在pentium(r) dual-core 處理器上耗費的時間己經(jīng)超過120分鐘。參考文獻劉會燈，朱 飛.matlab編程基礎(chǔ)與典型應(yīng)用m 2008年3月附錄% % program 1 (找出切片的邊界線的圖)pic=zeros(512,512,100); pic2=zeros(512,512,100);for i=0:99s=sprintf(,d:apics%d.bmp,i);pic(:,:,i+1)=imread(s);pic2(:,:,i+1)=edge(pic(:,:,i+1); imshow(pic2(:,:,i+1);endm,n,i+1

21、=find(pic2(:,:ji4-1);% program 2 (求半徑，找中軸線與切片交點)ss=100;lens=1000*ones(512,512,ss);centres=zeros(ss,2);r=zeros(ss,1);for k=1:ssfor i=1:512for j=1:512if pic(i,j,k)=0for m=1:512for n=1:512if pic2(m,n,k)=1t1=sqrt(i-m)*(i-m)+(j-n)*(j-n);if iens(i,j,k)>t1iens(i,j,k)=t1;endendendendendendendfor i=1:512f

22、or j=1:512if lens(i,j,k)=1000lens(i,j,k)=-1;endendendr(k,1)=max(ma x(len s(:,:,k);forj=1:512for i=1:512if r(k,1)-lens(i,j,k)<0.0001cen tres(k,1)=i;cen tres(k,2)=j;endendendenddisp('this program has been done!*)% program 3 (多項式擬合空間屮軸線)format longcenn=zeros(100,3);x=cen(:,1);y=ce n(:,2);z=ce n(:,3);px=polyfit(乙 x,9);x1=polyval(px,z);py=polyfit(z,y,9);y1=polyval(py,z);figure(1);plot3(x1 ,y1,z)cenn(:,1 )=x1 ;cenn (:,2)=y1 ;ce nn (:,3)=z;% % program 4(繪制血管三維圖像)t=lin space(0,pi,25);p=li