MATLAB的血管三维重建源代码Word格式.docx

1、zhongxindian=zeros（100,2）;r=zeros（100,1）;for k=0:99 t=strcat（f:/,int2str（i）,.bmp）;b=imread（t）;b=zhuanhua（b）;%将01互换blunkuo=edge（b,sobel%提取轮廓bgujia=bwmorph（b,skel,inf）;%提取骨架%寻找内切圆x0,y0,v0=find（b0lunkuo）;a0,b0,c0=find（b0gujia）;m=length（a0）;n=length（x0）;juli=zeros（m,n）;cunfang=zeros（m,2）;mn p1=a0（i）;q1=

2、b0（i）; p2=x0（j）;q2=y0（j）;juli（i,j）=sqrt（a（p1,q1）-a（p2,q2）2+（b（p1,q1）-b（p2,q2）2）;%骨架上的各个点到轮廓的距离 zx,zxxh=min（juli（i,:）;%骨架上一点到轮廓的最短距离即以骨架上各个点为圆心的内切园的半径 cunfang（i,1）=zx; cunfang（i,2）=zxxh;zd,zdxh=max（cunfang（:,1）;%寻找半径中最大的半径和其对应的圆心坐标g=a0（zdxh）;h=b0（zdxh）;zhongxindian（k+1,1）=a（g,h）;zhongxindian（k+1,2）=b

3、（g,h）;r（k+1）=zd;附录3：通过计算不同次数多项式拟合的偏差平方和确定拟和次数的matlab程序代码（pczx.m）function j=pczx（z,t） %根据不同次数的多项式拟合与原图数据偏差平方和的大小来确定多项式拟和的次数delta=zeros（10,1）;for k=1:10p,s=polyfit（z,t,k）;delta（k）=s.normri,j=min（delta）;附录4：根据轮廓画出血管的三维图像的matlab程序代码for b=0:99 %提取原图的轮廓，根据轮廓画出血管的三维图像 m1=imread（int2str（b）,）; m（:,:,b+1）=edg

4、e（m1, for i=1: if （m（i,j,k+1）=1） plot3（i,j,k+1,r-. hold ongrid ontitle（血管三维图）rotate3dhold off附录5：绘制中轴线及在各平面的投影图matlab程序代码format longpx=polyfit（z,x,7）;%x,z的7次多项式拟合x1=polyval（px,z）;py=polyfit（z,y,5）;%y,z的5次多项式拟合y1=polyval（py,z）;figure（1）; %画中心轴线图plot3（x1,y1,z）xlabel（X轴ylabel（Y轴zlabel（Z轴血管中轴线图figure（2）

5、; %画中心轴线在xoz平面上的投影plot（z,x1,-r血管中轴线XOZ平面投影图figure（3）;%画中心轴线在yoz平面上的投影 plot（z,y1,-b血管中轴线YOZ平面投影图figure（4）;%画中心轴线在xoy平面上的投影plot（x1,y1,-g血管中轴线XOY平面投影图附录6：求第pn张拟合图的轮廓的二值矩阵的matlab程序代码（dian.m）function pnjj=dian（px,py,pn） %输出pnjj，为第pn张拟合图片的轮廓二值矩阵a=zeros（991）;b=zeros（991）;q=zeros（991）;w=zeros（991）;r=zeros（1

6、,2）;s=zeros（1,2）;k=1;for c=0:0.1:a（k）=7*px（1）*c.6+6*px（2）*c.5+5*px（3）*c.4+4*px（4）*c.3+3*px（5）*c.2+2*px（6）*c+px（7）;b（k）=4*py（1）*c.3+3*py（2）*c.2+2*py（3）*c+py（4）; %中心轴线方程关于z的导数即a（k）,b（k）,1为z在k处的切线的方向向量q（k）=px（1）*c.7+px（2）*c.6+px（3）*c.5+px（4）*c.4+px（5）*c.3+px（6）*c.2+px（7）*c+px（8）;w（k）=py（1）*c.5+py（2）*c.

7、4+py（3）*c.3+py（4）*c.2+py（5）*c; %中心轴线方程在z=k处的x,y值k=k+1;%提取新的截痕u=;v=;syms x yfor i=0: m=a（k）*（x-q（k）+b（k）*（y-w（k）+（pn-i）; n=（x-q（k）2+（y-w（k）2+（pn-i）2-29.492; g,h=solve（m,n）; r=double（g）; s=double（h）;if （abs（imag（r（1）0.01） %去除复数根 u=u;real（r（1）+256 real（r（2）+256; v=v;real（s（1）+256 real（s（2）+256;%根据新的切平面

8、的轮廓坐标得到新轮廓的图像矩阵plot（v（:,1）,u（:,1）,r.,v（:,2）,u（:,2）,axis（0 512 0 512）;pnj=imread（strcat（int2str（pn）,lk=edge（pnj,pnjj=zeros（512）;u=round（u）;v=round（v）;for t=1:length（u（:,1）pnjj（u（t,1）,v（t,1）=1;pnjj（u（t,2）,v（t,2）=1;%画原图轮廓imshow（lk）%画新图轮廓imshow（pnjj）%画原图与新图的轮廓图对比imshow（pnjj|lk）pnjj=zhuanhua（pnjj）;附录7：求拟

9、合图与原切片图的重合度的matlab程序代码（baifenbi1.m）function baifenbi=baifenbi1（pnjj,pn） %输出拟合图与原切片图的重合度%填充新图juzheng=pnjj; %pnjj为通过dian.m得到的轮廓边界二值矩阵%先填充左边界的右半部分511 if（pnjj（i,j）=0&pnjj（i,j+1）=0） pnjj（i,j+1）=pnjj（i,j）;you=pnjj;%再填充右边界的左半部分 for j=512:-1:2 if（juzheng（i,j）=0&juzheng（i,j-1）=0） juzheng（i,j-1）=juzheng（i,j）;zuo=juzheng;shijijuzheng=you|zuo; %通过矩阵的或运算得到填充后的新图imshow（you|zuo）%原图的黑点的个数biaozhunjuzheng=imread（strcat（int2str（pn）,nbiao=0; if（biaozhunjuzheng（i,j）=0） nbiao=nbiao+1;%新图与原图重合部分黑点的个数chonghegs=0; if（biaozhunjuzheng（i,j）=0&shijijuzheng（i,j）=0） chonghegs=chonghegs+1;%求百分比baifenbi=chonghegs/nbiao;