第1页 / 共8页
第2页 / 共8页
第3页 / 共8页
第4页 / 共8页
第5页 / 共8页
第6页 / 共8页
第7页 / 共8页
第8页 / 共8页
血管的三维重建源程序，对于matlab初学者.doc
附录1：图像二值矩阵的0-1互换的matlab程序代码（zhuanhua.m）
附录2：求各切片的最大内切圆的半径及圆心坐标matlab程序代码（ff.m）
附录3：通过计算不同次数多项式拟合的偏差平方和确定拟和次数的matlab程序代码（pczx.m）
附录4：根据轮廓画出血管的三维图像的matlab程序代码
附录5：绘制中轴线及在各平面的投影图matlab程序代码
附录6：求第pn张拟合图的轮廓的二值矩阵的matlab程序代码（dian.m）
附录7：求拟合图与原切片图的重合度的matlab程序代码（baifenbi1.m）
附录 1：图像二值矩阵的 0-1 互换的 matlab 程序代码（zhuanhua.m）
function b0=zhuanhua(b0) %图像二值矩阵的 0-1 互换
for i=1:512
for j=1:512
if b0(i,j)==1
b0(i,j)=0;
else b0(i,j)=1;
end
end
end
附录 2：求各切片的最大内切圆的半径及圆心坐标 matlab 程序代码（ff.m）
function [r, zhongxindian]=ff %输出各切片最大内切圆半径及圆心坐标
a=zeros(512,512);
b=zeros(512,512);
for i=1:512
for j=1:512
a(i,j)=i-257; %横坐标的对应
b(i,j)=j-257; %纵坐标的对应
end
end
%图像在 xyz 面上的 x 轴、y 轴坐标
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);
for i=1:m
for j=1: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);%骨架上的各个点到轮廓的距离
end
[zx,zxxh]=min(juli(i,:));%骨架上一点到轮廓的最短距离即以骨架上各个点为圆心的内切
园的半径
cunfang(i,1)=zx;
cunfang(i,2)=zxxh;
end
[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;
end
附录 3：通过计算不同次数多项式拟合的偏差平方和确定拟和次数的 matlab 程
序代码（pczx.m）
function j=pczx(z,t) %根据不同次数的多项式拟合与原图数据偏差平方和的大小来确定多
项式拟和的次数
delta=zeros(10,1);
for k=1:10
[p,s]=polyfit(z,t,k);
delta(k)=s.normr
end
[i,j]=min(delta);
附录 4：根据轮廓画出血管的三维图像的 matlab 程序代码
for b=0:99 %提取原图的轮廓，根据轮廓画出血管的三维图像
m1=imread([int2str(b),'.bmp']);
m(:,:,b+1)=edge(m1,'sobel');
end
for k=0:99
for i=1:512
for j=1:512
if (m(i,j,k+1)==1)
plot3(i,j,k+1,'r-.');
hold on
end
end
end
end
grid on
title('血管三维图')
rotate3d
hold off
附录 5：绘制中轴线及在各平面的投影图 matlab 程序代码
format long
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)
grid on
xlabel('X 轴');
ylabel('Y 轴');
zlabel('Z 轴');
title('血管中轴线图');
figure(2); %画中心轴线在 xoz 平面上的投影
plot(z,x1,'-r')
ylabel('Z 轴');
xlabel('X 轴')
title('血管中轴线 XOZ 平面投影图');
grid on
figure(3);%画中心轴线在 yoz 平面上的投影
plot(z,y1,'-b')
xlabel('Z 轴');
ylabel('Y 轴');
title('血管中轴线 YOZ 平面投影图');
grid on
figure(4);%画中心轴线在 xoy 平面上的投影
plot(x1,y1,'-g')
xlabel('X 轴');
ylabel('Y 轴');
title('血管中轴线 XOY 平面投影图');
grid on
附录 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,2);
s=zeros(1,2);
k=1;
for c=0:0.1:99
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;
end
%提取新的截痕
u=[];
v=[];
syms x y
k=1;
for i=0:0.1:99
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))+256 real(r(2))+256]];
v=[v;[real(s(1))+256 real(s(2))+256]];
end
k=k+1;
end
%根据新的切平面的轮廓坐标得到新轮廓的图像矩阵
plot(v(:,1),u(:,1),'r.',v(:,2),u(:,2),'r.')
axis([0 512 0 512]);
pnj=imread(strcat(int2str(pn),'.bmp'));
lk=edge(pnj,'sobel');
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;
end
figure(1);%画原图轮廓
imshow(lk)
figure(2);%画新图轮廓
imshow(pnjj)
figure(3);%画原图与新图的轮廓图对比
imshow(pnjj|lk)
pnjj=zhuanhua(pnjj);
附录 7：求拟合图与原切片图的重合度的 matlab 程序代码（baifenbi1.m）
function baifenbi=baifenbi1(pnjj,pn) %输出拟合图与原切片图的重合度
%填充新图
juzheng=pnjj; %pnjj 为通过 dian.m 得到的轮廓边界二值矩阵
%先填充左边界的右半部分
for i=1:511
for j=1:511
if(pnjj(i,j)==0&pnjj(i,j+1)~=0)
pnjj(i,j+1)=pnjj(i,j);
end
end
end
you=pnjj;
%再填充右边界的左半部分
for i=1:512
for j=512:-1:2
if(juzheng(i,j)==0&juzheng(i,j-1)~=0)
juzheng(i,j-1)=juzheng(i,j);
end
end
end
zuo=juzheng;
shijijuzheng=you|zuo; %通过矩阵的或运算得到填充后的新图
imshow(you|zuo)
%原图的黑点的个数
biaozhunjuzheng=imread(strcat(int2str(pn),'.bmp'));
nbiao=0;
for i=1:512
for j=1:512
if(biaozhunjuzheng(i,j)==0)
nbiao=nbiao+1;
end
end
end
%新图与原图重合部分黑点的个数
chonghegs=0;
for i=1:512
for j=1:512
if(biaozhunjuzheng(i,j)==0&shijijuzheng(i,j)==0)
chonghegs=chonghegs+1;
end
end
end
%求百分比
baifenbi=chonghegs/nbiao;
相关推荐
2023年江西萍乡中考道德与法治真题及答案.doc
2012年重庆南川中考生物真题及答案.doc
2013年江西师范大学地理学综合及文艺理论基础考研真题.doc
2020年四川甘孜小升初语文真题及答案I卷.doc
2020年注册岩土工程师专业基础考试真题及答案.doc
2023-2024学年福建省厦门市九年级上学期数学月考试题及答案.doc
2021-2022学年辽宁省沈阳市大东区九年级上学期语文期末试题及答案.doc
2022-2023学年北京东城区初三第一学期物理期末试卷及答案.doc
2018上半年江西教师资格初中地理学科知识与教学能力真题及答案.doc
2012年河北国家公务员申论考试真题及答案-省级.doc
2020-2021学年江苏省扬州市江都区邵樊片九年级上学期数学第一次质量检测试题及答案.doc
2022下半年黑龙江教师资格证中学综合素质真题及答案.doc
