三維血管重建論文含源代碼

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

2、陣，最大值即為血管半徑。最后計算所有切片的血管半徑，并對這些半徑求平均值，得到平均血管半徑為：29.6799m。由100張切片的最大內(nèi)切圓圓心坐標(biāo)擬合得出中軸線方程以及其在三個坐標(biāo)平面內(nèi)的投影曲線方程。由中軸線得到血管的三維立體重建圖，用平面去截血管的三維立體重建圖，得到新的4張截面圖。把它們分別與題設(shè)中的對應(yīng)截面進(jìn)行內(nèi)點個數(shù)對比。我們定義兩張切片所共同擁有的內(nèi)點個數(shù)與原切片內(nèi)點個數(shù)的比值為重合度。計算得到平均重合度為：98.19% 。關(guān)鍵詞：血管半徑 中軸線 切片 重建（來自作者：歡迎各界人士批評指正，。文章作于2011年8月10日，陜西科技大學(xué)理學(xué)院實驗室）1問題的重述斷面可用于了解生物組

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

4、設(shè)：管道中軸線與每張切片有且只有一個交點；球半徑固定；切片間距以及圖像象素的尺寸均為1。計算管道的中軸線與半徑，給出具體的算法，并繪制中軸線在XY、YZ、ZX平面的投影圖。 2模型假設(shè) (1) 假設(shè)血管中軸線與每張切片有且只有一個交點。(2) 假設(shè)血管的半徑固定不變。(3) 假設(shè)在對切片拍照的過程中不存在誤差。(4) 假設(shè)血管不存在嚴(yán)重扭曲。3符號說明第張切片中輪廓線的最大內(nèi)切圓圓心第張切片中輪廓線的最大內(nèi)切圓半徑第張切片中第個內(nèi)點與第個輪廓點的距離100張切片圖片轉(zhuǎn)換以后的三維0-1矩陣100張切片圖片的輪廓線生成的矩陣100張切片輪廓的最大內(nèi)切圓圓心坐標(biāo)平均血管半徑第張

5、切片圖片中所有內(nèi)點的集合第張切片圖片中輪廓點的集合原切片圖片的上內(nèi)點及邊界點的集合重新切片得到的內(nèi)點及邊界點的集合R-square多項式擬合的指標(biāo)數(shù)4問題的分析根據(jù)題目整個管道是由球心沿著某一曲線（稱為中軸線）的球滾動包絡(luò)而成?；趲缀卧砜梢缘贸鲆韵聝蓷l定理。（1）球的任意截面都是圓。（2）經(jīng)過球心的球截面是所有的截面圓當(dāng)中半徑最大的圓?；谏鲜鰞蓚€定理可以得出：每張切片的最大內(nèi)切圓的圓心位于血管的中軸線上，該圓的半徑等于血管半徑。分別求出100張切片輪廓線的最大內(nèi)切圓圓心與半徑。對所有最大內(nèi)切圓的圓心做擬合，即可求出血管的中軸線。為減少誤差，我們用求均值的方法來得到平均血管半徑。根據(jù)擬合出

6、的中軸線方程，即可求出中軸線在XY、YZ、ZX平面上的投影。5模型的建立與求解5.1模型建立求血管的半徑（1）導(dǎo)入數(shù)據(jù)，轉(zhuǎn)換存儲方式為了計算方便，利用MATLAB1軟件將100張BMP文件轉(zhuǎn)換存儲為一個三維0-1矩陣（0代表黑色像素點，1代表白色像素點）。（2）求切片輪廓線上各點坐標(biāo)利用MATLAB軟件的內(nèi)部函數(shù)edge，求得所有切片圖片的輪廓線生成的矩陣（3）求切片輪廓線的最大內(nèi)切圓的半徑首先求出每個內(nèi)點距輪廓線的距離，取其中的最小值，即為以這些內(nèi)點為圓心的輪廓線的內(nèi)切圓半徑；其次找出所有內(nèi)點確定的最小內(nèi)切圓中半徑最大的一個，即為輪廓線的最大內(nèi)切圓半徑。具體步驟：Step1在三維0-1矩陣中

7、遍歷搜索內(nèi)點,將其存儲于集合中 ； Step2 在輪廓線矩陣中遍歷搜索輪廓點，將其存儲于集合中； Step3計算集合中第個內(nèi)點到集合中第個輪廓點的距離 ，首先取每個內(nèi)點到輪廓點的最小距離，即為以這些內(nèi)點為圓心的內(nèi)切圓半徑；其次找出所有內(nèi)點確定的最小內(nèi)切圓中半徑最大的一個，即為輪 廓線的最大內(nèi)切圓半徑Step4計算平均血管半徑（4）求切片輪廓線的最大內(nèi)切圓心坐標(biāo)將100張輪廓線的最大內(nèi)切圓半徑所對應(yīng)的圓心坐標(biāo)記錄于矩陣中。 由圓心坐標(biāo)求中軸線方程及曲線投影（1）利用MATLAB畫圖工具由100組圓心坐標(biāo)畫出中軸線以及中軸線在XY、YZ、ZX面的投影，由于圖像文件以象素為單位存儲數(shù)據(jù)，這使得圖像的

8、原始數(shù)據(jù)是離散化的。對于離散的位圖來說，將變得不可分辨，可能切片輪廓線的最大圓不止一個。為了減少誤差，將中軸線在XY、YZ、ZX面的曲線方程進(jìn)行多項式擬合，建立擬合方程：（2）由建立的投影曲面得到中軸線的方程5.2 模型求解根據(jù)所建模型求出100張切片的最大內(nèi)切圓的圓心坐標(biāo)、半徑以及血管半徑表1最大內(nèi)切圓半徑及圓心坐標(biāo)切片序號最大內(nèi)切圓半徑（R/m）最大內(nèi)切圓圓心坐標(biāo)橫 坐標(biāo)(X/m)縱 坐標(biāo)(Y/m)豎 坐標(biāo)（Z/m）029069-16011229-16022329069-16024529155-160267

9、29.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-16015242529.833-16016252629614614-15928

10、2829296146146146146147327327327326145475475353682682682698-115

11、116474829.698-113118484929.698-112119495029.698-111120505129.698-111120515229.698-111120525329.698-111120535429.682-111120545529.53-66150555629.53-59153565729.732-70148575829.967-70148585929.732-54155596029.732-54155606129.614-31162616229.614-31162626329.614-6166636429.833-6166646529.833816765663010

12、16766673011167676830121676869301216769703012167707129.83326166717229.61426166727329.61446163737429.73265158747529.73265158757629.73271156767729.68296145777829.68296145787929.698127124798029.698127124808129.732122128818229.732122128828329.732122128838429.698127124848529.732150101858629.73215010186872

13、9.73215496878829.698137115888929.698139113899029.68216088909129.68216285919229.54717659929329.61418340939429.73218049949529.61418340959629.61418340969729.614185339798301900989930190199（1）血管半徑為=29.6799m。（2）中軸線在XY平面投影曲線的擬合方程(R-square:0.9939)（3）中軸線在YZ平面投影曲線的擬合方程(R-square:1)（4）中軸線在ZX平面投影曲線的擬合方程(R-square

14、:1)（5）中軸線方程： 圖1中軸線曲線 圖2中軸線擬合曲線圖3中軸線在XY平面的投影 圖4中軸在XY平面的投影擬合曲線 圖5中軸線在YZ平面的投影 圖6中軸線在YZ平面的投影擬合曲線 圖7中軸線在ZX平面的投影 圖8中軸線在ZX平面的投影擬合曲線圖9曲線模擬的血管6模型檢驗與分析用多項式擬合得到重建后血管的三維空間曲線，再用平面去截空間曲線得到新的4層切片平面分別與題目中對應(yīng)切片進(jìn)行對比，在這里我們稱的面積與的面積比值為切片重合度。第01、25、50、75張圖片分別與對應(yīng)的新截取圖片的對比圖.圖1 第01張切片 圖11 第25張切片 圖12 第50張切片 圖13 第75張切片圖.圖14截后第

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

16、ntium(R) Dual-Core處理器上耗費的時間已經(jīng)超過120分鐘。參考文獻(xiàn)1劉會燈，朱 飛.MATLAB編程基礎(chǔ)與典型應(yīng)用M 2008年3月附錄% program1（找出切片的邊界線的圖）pic=zeros(512,512,100);pic2=zeros(512,512,100);for i=0:99 s=sprintf('d:Apics%d.bmp',i); pic(:,:,i+1)=imread(s); pic2(:,:,i+1)=edge(pic(:,:,i+1); imshow(pic2(:,:,i+1);end%program 2（求半徑，找中軸線與切片交點）

17、ss=100;lens=1000*ones(512,512,ss);centres=zeros(ss,2);r=zeros(ss,1);for k=1:ss for i=1:512 for j=1:512 if pic(i,j,k)=0 for m=1:512 for n=1:512 if pic2(m,n,k)=1 t1=sqrt(i-m)*(i-m)+(j-n)*(j-n); if lens(i,j,k)>t1 lens(i,j,k)=t1; end end end end end end end for i=1:512 for j=1:512 if lens(i,j,k)=1000

18、 lens(i,j,k)=-1; end end end r(k,1)=max(max(lens(:,:,k); for j=1:512 for i=1:512 if r(k,1)-lens(i,j,k)<0.0001 centres(k,1)=i; centres(k,2)=j; end end endenddisp('This program has been done!')%program 3（多項式擬合空間中軸線）format longcenn=zeros(100,3);x=cen(:,1);y=cen(:,2);z=cen(:,3);px=polyfit(z,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;cenn(:,3)=z;%program 4（繪制血管三維圖像）t=linspace(0,pi,25);p=linspace(0,2*pi,25);theta,phi=meshgrid(t,p);for i=1:100 x=29.6799*sin(theta).*sin(phi)+cenn(i,1); y=29.6799*sin(theta).*cos(phi)+cenn(