基于MATLAB的血管三維重建源代碼

1、附錄1:圖像二值矩陣的0-1互換的matlab程序代碼(zhuanhua.m)functionb0=zhuanhua(b0)%圖像二值矩陣的0-1互換fori=1:512forj=1:512ifb0(i,j)=1b0(i,j)=0;elseb0(i,j)=1;endendend附錄2：求各切片的最大內(nèi)切圓的半徑及圓心坐標(biāo)matlab程序代碼(ff.m)functionr,zhongxindian=ff%輸出各切片最大內(nèi)切圓半徑及圓心坐標(biāo)a=zeros(512,512);b=zeros(512,512);fori=1:512forj=1:512a(i,j)=i-257;%橫坐標(biāo)的對應(yīng)b(i,j)

2、=j-257;%縱坐標(biāo)的對應(yīng)endend%圖像在xyz面上的x軸、y軸坐標(biāo)zhongxindian=zeros(100,2);r=zeros(100,1);fork=0:99t=strcat('f:/',int2str(i),'.bmp');b=imread(t);b=zhuanhua(b);%將01互換blunkuo=edge(b,'sobel');%提取輪廓bgujia=bwmorph(b,'skel',inf);%提取骨架%尋找內(nèi)切圓x0,y0,v0=find(b0lunkuo);a0,b0,c0=find(b0gujia)

3、;m=length(a0);n=length(x0);juli=zeros(m,n);cunfang=zeros(m,2);fori=1:mforj=1:np1=a0(i);q1=b0(i);p2=x0(j);q2=y0(j);juli(i,j)=sqrt(a(p1,q1)-a(p2,q2)A2+(b(p1,q1)-b(p2,q2)A2);%骨架上的各個(gè)點(diǎn)到輪廓的距離endzx,zxxh=min(juli(i,:);%骨架上一點(diǎn)到輪廓的最短距離即以骨架上各個(gè)點(diǎn)為圓心的內(nèi)切園的半徑cunfang(i,1)=zx;cunfang(i,2)=zxxh;endzd,zdxh=max(cunfang(:

4、,1);%尋找半徑中最大的半徑和其對應(yīng)的圓心坐標(biāo)g=a0(zdxh);h=b0(zdxh);zhongxindian(k+1,1)=a(g,h);zhongxindian(k+1,2)=b(g,h);r(k+1)=zd;end附錄3：通過計(jì)算不同次數(shù)多項(xiàng)式擬合的偏差平方和確定擬和次數(shù)的matlab程序代碼(pczx.m)functionj=pczx(z,t)%根據(jù)不同次數(shù)的多項(xiàng)式擬合與原圖數(shù)據(jù)偏差平方和的大小來確定多項(xiàng)式擬和的次數(shù)delta=zeros(10,1);fork=1:10p,s=polyfit(z,t,k);delta(k)=s.normrendi,j=min(delta);附錄4

5、：根據(jù)輪廓畫出血管的三維圖像的matlab程序代碼forb=0:99%提取原圖的輪廓，根據(jù)輪廓畫出血管的三維圖像m1=imread(int2str(b),'.bmp');m(:,:,b+1)=edge(m1,'sobel');endfork=0:99fori=1:512forj=1:512if(m(i,j,k+1)=1)plot3(i,j,k+1,'r-.');holdonendendendendgridontitle('血管三維圖')rotate3dholdoff附錄5:繪制中軸線及在各平面的投影圖matlab程序代碼forma

6、tlongpx=polyfit(z,x,7);%x,z的7次多項(xiàng)式擬合x1=polyval(px,z);py=polyfit(z,y,5);%y,z的5次多項(xiàng)式擬合y1=polyval（py,z）;figure（1）;%畫中心軸線圖plot3（x1,y1,z）gridonxlabel（'X軸'）；ylabel（'Y軸'）;zlabel（'Z軸'）;title（'血管中軸線圖'）;figure（2）;%畫中心軸線在xoz平面上的投影plot（z,x1,'-r'）ylabel（'Z軸'）；xlabel（

7、'X軸'）title（'血管中軸線XOZ平面投影圖'）;gridonfigure（3）;%畫中心軸線在yoz平面上的投影plot（z,y1,'-b'）xlabel（'Z軸'）；ylabel（'Y軸'）；title（'血管中軸線YOZ平面投影圖'）;gridonfigure（4）;%畫中心軸線在xoy平面上的投影plot（x1,y1,'-g'）xlabel（'X軸'）；ylabel（'Y軸'）;title（'血管中軸線XOY平面投影圖'）

8、;gridon附錄6:求第pn張擬合圖的輪廓的二值矩陣的matlab程序代碼（dian.m）functionpnjj=dian（px,py,pn）%輸出pnjj,為第pn張擬合圖片的輪廓二值矩陣a=zeros（991）;b=zeros(991);q=zeros(991);w=zeros(991);r=zeros(1,2);s=zeros(1,2);k=1;forc=0:0,1:99a(k)=7*px(1)*c.A6+6*px(2)*c.A5+5*px(3)*c.A4+4*px(4)*c.A3+3*px(5)*c.A2+2*px(6)*c+px(7);b(k)=4*py(1)*c.A3+3*py

9、(2)*c.A2+2*py(3)*c+py(4);%中心軸線方程關(guān)于z的導(dǎo)數(shù)即a(k),b(k),1為z在k處的切線的方向向量q(k)=px(1)*c.A7+px(2)*c.A6+px(3)*c.A5+px(4)*c.A4+px(5)*c.A3+px(6)*c.A2+px(7)*c+px(8);w(k)=py(1)*c.A5+py(2)*c.A4+py(3)*c.A3+py(4)*c.A2+py(5)*c;%中心軸線方程在z=k處的x,y值k=k+1;end%提取新的截痕u=;v=;symsxyk=1;fori=0:0.1:99m=a(k)*(x-q(k)+b(k)*(y-w(k)+(pn-i

10、);n=(x-q(k)A2+(y-w(k)A2+(pn-i)A2-29.49A2;g,h=solve(m,n);r=double(g);s=double(h);if(abs(imag(r(1)<0.01)%去除復(fù)數(shù)根u=u;real(r(1)+256real(r(2)+256;v=v;real(s(1)+256real(s(2)+256;endk=k+1;end%根據(jù)新的切平面的輪廓坐標(biāo)得到新輪廓的圖像矩陣plot(v(:,1),u(:,1),'r.',v(:,2),u(:,2),'r.')axis(05120512);pnj=imread(strcat(

11、int2str(pn),'.bmp');lk=edge(pnj,'sobel');pnjj=zeros(512);u=round(u);v=round(v);fort=1:length(u(:,1)pnjj(u(t,1),v(t,1)=1;pnjj(u(t,2),v(t,2)=1;endfigure(1);%畫原圖輪廓imshow(lk)figure(2);%畫新圖輪廓imshow(pnjj)figure(3);%畫原圖與新圖的輪廓圖對比imshow(pnjj|lk)pnjj=zhuanhua(pnjj);附錄7：求擬合圖與原切片圖的重合度的matlab程序代碼

12、(baifenbi1.m)functionbaifenbi=baifenbi1(pnjj,pn)%輸出擬合圖與原切片圖的重合度%填充新圖juzheng=pnjj;%pnjj為通過dian.m得到的輪廓邊界二值矩陣%先填充左邊界的右半部分fori=1:511forj=1:511if(pnjj(i,j)=0&pnjj(i,j+1)=0)endendendyou=pnjj;%再填充右邊界的左半部分fori=1:512forj=512:-1:2if(juzheng(i,j)=0&juzheng(i,j-1)=0)juzheng(i,j-1)=juzheng(i,j);endendendzuo=juzheng;shijijuzheng=you|zuo;%通過矩陣的或運(yùn)算得到填充后的新圖imshow(you|zuo)%原圖的黑點(diǎn)的個(gè)數(shù)biaozhunjuzheng=imread(strcat(int2str(pn),'.bmp');nbiao=0;fori=1:512forj=1:512