MATLAB的血管三维重建源代码Word格式.docx
《MATLAB的血管三维重建源代码Word格式.docx》由会员分享，可在线阅读，更多相关《MATLAB的血管三维重建源代码Word格式.docx（8页珍藏版）》请在冰豆网上搜索。
zhongxindian=zeros（100,2）;
r=zeros（100,1）;
fork=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）;
m
n
p1=a0（i）;
q1=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（g,h）;
r（k+1）=zd;
附录3：
通过计算不同次数多项式拟合的偏差平方和确定拟和次数的matlab程序代码（pczx.m）
functionj=pczx（z,t）%根据不同次数的多项式拟合与原图数据偏差平方和的大小来确定多项式拟和的次数
delta=zeros（10,1）;
fork=1:
10
[p,s]=polyfit（z,t,k）;
delta（k）=s.normr
[i,j]=min（delta）;
附录4：
根据轮廓画出血管的三维图像的matlab程序代码
forb=0:
99%提取原图的轮廓，根据轮廓画出血管的三维图像
m1=imread（[int2str（b）,'
]）;
m（:
:
b+1）=edge（m1,'
fori=1:
if（m（i,j,k+1）==1）
plot3（i,j,k+1,'
r-.'
holdon
gridon
title（'
血管三维图'
）
rotate3d
holdoff
附录5：
绘制中轴线及在各平面的投影图matlab程序代码
formatlong
px=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）;
%画中心轴线在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）
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:
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.^4+py（3）*c.^3+py（4）*c.^2+py（5）*c;
%中心轴线方程在z=k处的x,y值
k=k+1;
%提取新的截痕
u=[];
v=[];
symsxy
fori=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.49^2;
[g,h]=solve（m,n）;
r=double（g）;
s=double（h）;
if（abs（imag（r
（1）））<
0.01）%去除复数根
u=[u;
[real（r
（1））+256real（r
（2））+256]];
v=[v;
[real（s
（1））+256real（s
（2））+256]];
%根据新的切平面的轮廓坐标得到新轮廓的图像矩阵
plot（v（:
1）,u（:
1）,'
r.'
v（:
2）,u（:
2）,'
axis（[05120512]）;
pnj=imread（strcat（int2str（pn）,'
lk=edge（pnj,'
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;
%画原图轮廓
imshow（lk）
%画新图轮廓
imshow（pnjj）
%画原图与新图的轮廓图对比
imshow（pnjj|lk）
pnjj=zhuanhua（pnjj）;
附录7：
求拟合图与原切片图的重合度的matlab程序代码（baifenbi1.m）
functionbaifenbi=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;
%再填充右边界的左半部分
forj=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;
升级会员 