三維血管重建建模含源代碼

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

2、新的矩陣。然后找到原圖片矩陣中像素點(diǎn)的內(nèi)點(diǎn)（切片圖片中輪廓線中的點(diǎn)），從而得到內(nèi)點(diǎn)到切片輪廓點(diǎn)的最小距離矩陣和最小距離中的最大值矩陣，最大值即為血管半徑。最后計(jì)算所有切片的血管半徑，并對這些半徑求平均值，得到平均血管半徑為：29.6799m。由100張切片的最大內(nèi)切圓圓心坐標(biāo)擬合得出中軸線方程以及其在三個(gè)坐標(biāo)平面內(nèi)的投影曲線方程。由中軸線得到血管的三維立體重建圖，用平面去截血管的三維立體重建圖，得到新的4張截面圖。把它們分別與題設(shè)中的對應(yīng)截面進(jìn)行內(nèi)點(diǎn)個(gè)數(shù)對比。我們定義兩張切片所共同擁有的內(nèi)點(diǎn)個(gè)數(shù)與原切片內(nèi)點(diǎn)個(gè)數(shù)的比值為重合度。計(jì)算得到平均重合度為：98.19% 。關(guān)鍵詞：血管半徑 中軸線 切片

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

4、像，記錄了管道與切片的交。圖像文件名依次為0.bmp、1.bmp、 99.bmp，格式均為BMP，寬、高均為512個(gè)象素（pixel）。為簡化起見，假設(shè)：管道中軸線與每張切片有且只有一個(gè)交點(diǎn)；球半徑固定；切片間距以及圖像象素的尺寸均為1。計(jì)算管道的中軸線與半徑，給出具體的算法，并繪制中軸線在XY、YZ、ZX平面的投影圖。 2模型假設(shè) (1) 假設(shè)血管中軸線與每張切片有且只有一個(gè)交點(diǎn)。(2) 假設(shè)血管的半徑固定不變。(3) 假設(shè)在對切片拍照的過程中不存在誤差。(4) 假設(shè)血管不存在嚴(yán)重扭曲。3符號(hào)說明第張切片中輪廓線的最大內(nèi)切圓圓心第張切片中輪廓線的最大內(nèi)切圓半徑第張切片中第

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

6、0張切片輪廓線的最大內(nèi)切圓圓心與半徑。對所有最大內(nèi)切圓的圓心做擬合，即可求出血管的中軸線。為減少誤差，我們用求均值的方法來得到平均血管半徑。根據(jù)擬合出的中軸線方程，即可求出中軸線在XY、YZ、ZX平面上的投影。5模型的建立與求解5.1模型建立5.1.1求血管的半徑（1）導(dǎo)入數(shù)據(jù)，轉(zhuǎn)換存儲(chǔ)方式為了計(jì)算方便，利用MATLAB1軟件將100張BMP文件轉(zhuǎn)換存儲(chǔ)為一個(gè)三維0-1矩陣（0代表黑色像素點(diǎn)，1代表白色像素點(diǎn)）。（2）求切片輪廓線上各點(diǎn)坐標(biāo)利用MATLAB軟件的內(nèi)部函數(shù)edge，求得所有切片圖片的輪廓線生成的矩陣（3）求切片輪廓線的最大內(nèi)切圓的半徑首先求出每個(gè)內(nèi)點(diǎn)距輪廓線的距離，取其中的最小值

7、，即為以這些內(nèi)點(diǎn)為圓心的輪廓線的內(nèi)切圓半徑；其次找出所有內(nèi)點(diǎn)確定的最小內(nèi)切圓中半徑最大的一個(gè)，即為輪廓線的最大內(nèi)切圓半徑。具體步驟：Step1在三維0-1矩陣中遍歷搜索內(nèi)點(diǎn),將其存儲(chǔ)于集合中 ； Step2 在輪廓線矩陣中遍歷搜索輪廓點(diǎn)，將其存儲(chǔ)于集合中； Step3計(jì)算集合中第個(gè)內(nèi)點(diǎn)到集合中第個(gè)輪廓點(diǎn)的距離 ，首先取每個(gè)內(nèi)點(diǎn)到輪廓點(diǎn)的最小距離，即為以這些內(nèi)點(diǎn)為圓心的內(nèi)切圓半徑；其次找出所有內(nèi)點(diǎn)確定的最小內(nèi)切圓中半徑最大的一個(gè)，即為輪 廓線的最大內(nèi)切圓半徑Step4計(jì)算平均血管半徑（4）求切片輪廓線的最大內(nèi)切圓心坐標(biāo)將100張輪廓線的最大內(nèi)切圓半徑所對應(yīng)的圓心坐標(biāo)記錄于矩陣中。 5.1.2由圓

8、心坐標(biāo)求中軸線方程及曲線投影（1）利用MATLAB畫圖工具由100組圓心坐標(biāo)畫出中軸線以及中軸線在XY、YZ、ZX面的投影，由于圖像文件以象素為單位存儲(chǔ)數(shù)據(jù)，這使得圖像的原始數(shù)據(jù)是離散化的。對于離散的位圖來說，將變得不可分辨，可能切片輪廓線的最大圓不止一個(gè)。為了減少誤差，將中軸線在XY、YZ、ZX面的曲線方程進(jìn)行多項(xiàng)式擬合，建立擬合方程：（2）由建立的投影曲面得到中軸線的方程5.2 模型求解根據(jù)所建模型求出100張切片的最大內(nèi)切圓的圓心坐標(biāo)、半徑以及血管半徑表1最大內(nèi)切圓半徑及圓心坐標(biāo)切片序號(hào)最大內(nèi)切圓半徑（R/m）最大內(nèi)切圓圓心坐標(biāo)橫 坐標(biāo)(X/m)縱 坐標(biāo)(Y/m)豎 坐標(biāo)（Z/m）029

9、069-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-16014222

10、330-16015232430-16015242529.833-16016252629614614614614614614614614732732732732614547547

1153682682682698-115116474829.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-5415

12、5606129.614-31162616229.614-31162626329.614-6166636429.833-6166646529.83381676566301016766673011167676830121676869301216769703012167707129.83326166717229.61426166727329.61446163737429.73265158747529.73265158757629.73271156767729.68296145777829.68296145787929.698127124798029.698127124808129.732122128

13、818229.732122128828329.732122128838429.698127124848529.732150101858629.732150101868729.73215496878829.698137115888929.698139113899029.68216088909129.68216285919229.54717659929329.61418340939429.73218049949529.61418340959629.61418340969729.614185339798301900989930190199（1）血管半徑為=29.6799m。（2）中軸線在XY平面投影

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

15、重合度。第01、25、50、75張圖片分別與對應(yīng)的新截取圖片的對比圖.圖1 第01張切片 圖11 第25張切片 圖12 第50張切片 圖13 第75張切片圖.圖14截后第01張切片圖15截后第25張切片圖16截后第50張切片圖17截后第75張切片根據(jù)上面曲線擬合得到的血管重新切片后，我們隨機(jī)取出四張重切的圖片，與原始的切片比較，得到兩切片的重合度。計(jì)算得：第01張切片的重合度為99.20% ，第25張切片的重合度為98.65% ，第50張切片的重合度為97.30% ，第75張切片的重合度為97.61% ，計(jì)算平均重合度得到98.19%，說明本模型較為可靠。7模型的優(yōu)點(diǎn)與缺點(diǎn)優(yōu)點(diǎn)在于本模型的算法

16、采用遍歷搜索的方法得到最大內(nèi)切圓的半徑的精度較高，但是由于遍歷搜索的運(yùn)算次數(shù)異常龐大，導(dǎo)致模型的缺點(diǎn)是實(shí)踐難度較大。從模型檢驗(yàn)的數(shù)據(jù)可以看出計(jì)算的精度，但是整個(gè)運(yùn)算在Pentium(R) Dual-Core處理器上耗費(fèi)的時(shí)間已經(jīng)超過120分鐘。參考文獻(xiàn)1劉會(huì)燈，朱 飛.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

17、); pic2(:,:,i+1)=edge(pic(:,:,i+1); imshow(pic2(:,:,i+1);end%program 2（求半徑，找中軸線與切片交點(diǎn)）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 le

18、ns(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 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（多項(xiàng)

19、式擬合空間中軸線）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)