毕业论文三维血管重建建模论文含源代码.doc

血管的三维重建
摘要
本文探讨血管的三维重建，由血管的相继100张平行切片图像计算血管的中轴线与半径，并绘制血管在三个坐标平面上的投影。
由于血管的表面是由球心沿着某一曲线（称为中轴线）的球滚动包络而成，由此我们得出结论：每个切片一定包含滚动球的大圆，并且它一定为切片的最大内切圆，而最大圆所对应的半径即为血管半径，所以求血管半径就转化为求每一个切片内部的点到切片外部轮廓的最大半径。
首先，读取100张血管切面图，把它们转换成Logical矩阵，从中提取切片截面轮廓点构成一个新的矩阵。然后找到原图片矩阵中像素点的内点（切片图片中轮廓线中的点），从而得到内点到切片轮廓点的最小距离矩阵和最小距离中的最大值矩阵，最大值即为血管半径。最后计算所有切片的血管半径，并对这些半径求平均值，得到平均血管半径为：29.6799μm。由100张切片的最大内切圆圆心坐标拟合得出中轴线方程以及其在三个坐标平面内的投影曲线方程。
由中轴线得到血管的三维立体重建图，用平面去截血管的三维立体重建图，得到新的4张截面图。把它们分别与题设中的对应截面进行内点个数对比。我们定义两张切片所共同拥有的内点个数与原切片内点个数的比值为重合度。计算得到平均重合度为：98.19% 。

关键词：血管半径 中轴线 切片 重建

（来自作者：欢迎各界人士批评指正，学术交流邮箱nibz@。文章作于2011年8月10日，陕西科技大学理学院实验室）
1问题的重述
断面可用于了解生物组织、器官等的形态。例如，将样本染色后切成厚约1μm的切片，在显微镜下观察该横断面的组织形态结构。如果用切片机连续不断地将样本切成数十、成百的平行切片， 可依次逐片观察。根据拍照并采样得到的平行切片数字图象，运用计算机可重建组织、器官等准确的三维形态。
假设某些血管可视为一类特殊的管道，该管道的表面是由球心沿着某一曲线（称为中轴线）的球滚动包络而成。例如圆柱就是这样一种管道，其中轴线为直线，由半径固定的球滚动包络形成。
现有某管道的相继100张平行切片图像，记录了管道与切片的交。图像文件名依次为0.bmp、1.bmp、…、 99.bmp，格式均为BMP，宽、高均为512个象素（pixel）。为简化起见，假设：管道中轴线与每张切片有且只有一个交点；球半径固定；切片间距以及图像象素的尺寸均为1。
计算管道的中轴线与半径，给出具体的算法，并绘制中轴线在XY、YZ、ZX平面的投影图。 
2模型假设 
(1) 假设血管中轴线与每张切片有且只有一个交点。
(2) 假设血管的半径固定不变。
(3) 假设在对切片拍照的过程中不存在误差。
(4) 假设血管不存在严重扭曲。
3符号说明
第张切片中轮廓线的最大内切圆圆心
第张切片中轮廓线的最大内切圆半径
第张切片中第个内点与第个轮廓点的距离
100张切片图片转换以后的三维0-1矩阵
100张切片图片的轮廓线生成的矩阵
100张切片轮廓的最大内切圆圆心坐标
平均血管半径
第张切片图片中所有内点的集合
第张切片图片中轮廓点的集合
原切片图片的上内点及边界点的集合
重新切片得到的内点及边界点的集合
R-square 多项式拟合的指标数

4问题的分析
根据题目整个管道是由球心沿着某一曲线（称为中轴线）的球滚动包络而成。
基于几何原理可以得出以下两条定理。
（1）球的任意截面都是圆。
（2）经过球心的球截面是所有的截面圆当中半径最大的圆。
基于上述两个定理可以得出：每张切片的最大内切圆的圆心位于血管的中轴线上，该圆的半径等于血管半径。
分别求出100张切片轮廓线的最大内切圆圆心与半径。
对所有最大内切圆的圆心做拟合，即可求出血管的中轴线。
为减少误差，我们用求均值的方法来得到平均血管半径。
根据拟合出的中轴线方程，即可求出中轴线在XY、YZ、ZX平面上的投影。
5模型的建立与求解
5.1模型建立
5.1.1 求血管的半径
（1）导入数据，转换存储方式
为了计算方便，利用MATLAB[1]软件将100张BMP文件转换存储为一个三维0-1矩阵（0代表黑色像素点，1代表白色像素点）。
（2）求切片轮廓线上各点坐标
利用MATLAB软件的内部函数edge，求得所有切片图片的轮廓线生成的矩阵
（3）求切片轮廓线的最大内切圆的半径
首先求出每个内点距轮廓线的距离，取其中的最小值，即为以这些内点为圆心的轮廓线的内切圆半径；其次找出所有内点确定的最小内切圆中半径最大的一个，即为轮廓线的最大内切圆半径。
具体步骤：
Step1 在三维0-1矩阵中遍历搜索内点,将其存储于集合中 ；
Step2 在轮廓线矩阵中遍历搜索轮廓点，将其存储于集合中；
Step3 计算集合中第个内点到集合中第个轮廓点的距离
，
首先取每个内点到轮廓点的最小距离，即为以这些内点为圆心的内切圆半径；其次找出所有内点确定的最小内切圆中半径最大的一个，即为轮 廓线的最大内切圆半径

Step4 计算平均血管半径

（4）求切片轮廓线的最大内切圆心坐标
将100张轮廓线的最大内切圆半径所对应的圆心坐标记录于矩阵中。

5.1.2 由圆心坐标求中轴线方程及曲线投影
（1）利用MATLAB画图工具由100组圆心坐标画出中轴线以及中轴线在XY、YZ、ZX面的投影，由于图像文件以象素为单位存储数据，这使得图像的原始数据是离散化的。对于离散的位图来说，将变得不可分辨，可能切片轮廓线的最大圆不止一个。为了减少误差，将中轴线在XY、YZ、ZX面的曲线方程进行多项式拟合，建立拟合方程：

（2）由建立的投影曲面得到中轴线的方程

5.2 模型求解
根据所建模型求出100张切片的最大内切圆的圆心坐标、半径以及血管半径
表1最大内切圆半径及圆心坐标
切片序号
最大内切圆半径
（R/μm）
最大内切圆圆心坐标
横 坐标
(X/μm)
纵 坐标
(Y/μm)
竖 坐标
（Z/μm）
0
29.069
-160
1
0
1
29.069
-160
1
1
2
29
-160
2
2
3
29.069
-160
2
3
4
29.069
-160
2
4
5
29.069
-160
2
5
6
29.155
-160
2
6
7
29.155
-160
2
7
8
29.275
-160
2
8
9
29.275
-160
3
9
10
29.428
-160
3
10
11
29.428
-160
4
11
12
29.614
-160
4
12
13
29.614
-160
5
13
14
29.833
-160
6
14
15
29.833
-160
7
15
16
29.833
-160
8
16
17
30
-160
9
17
18
30
-160
10
18
19
30
-160
11
19
20
30
-160
12
20
21
30
-160
13
21
22
30
-160
14
22
23
30
-160
15
23
24
30
-160
15
24
25
29.833
-160
16
25
26
29.614
-159
27
26
27
29.614
-159
27
27
28
29.614
-159
28
28
29
29.614
-158
34
29
30
29.614
-158
34
30
31
29.614
-157
40
31
32
29.614
-155
48
32
33
29.614
-155
48
33
34
29.614
-155
48
34
35
29.732
-153
55
35
36
29.732
-152
58
36
37
29.732
-152
58
37
38
29.732
-152
58
38
39
29.614
-152
58
39
40
29.547
-148
68
40
41
29.547
-149
66
41
42
29.53
-140
84
42
43
29.53
-140
84
43
44
29.682
-135
92
44
45
29.682
-135
92
45
46
29.682
-135
92
46
47
29.698
-115
116
47
48
29.698
-113
118
48
49
29.698
-112
119
49
50
29.698
-111
120
50
51
29.698
-111
120
51
52
29.698
-111
120
52
53
29.698
-111
120
53
54
29.682
-111
120
54
55
29.53
-66
150
55
56
29.53
-59
153
56
57
29.732
-70
148
57
58
29.967
-70
148
58
59
29.732
-54
155
59
60
29.732
-54
155
60
61
29.614
-31
162
61
62
29.614
-31
162
62
63
29.614
-6
166
63
64
29.833
-6
166
64
65
29.833
8
167
65
66
30
10
167
66
67
30
11
167
67
68
30
12
167
68
69
30
12
167
69
70
30
12
167
70
71
29.833
26
166
71
72
29.614
26
166
72
73
29.614
46
163
73
74
29.732
65
158
74
75
29.732
65
158
75
76
29.732
71
156
76
77
29.682
96
145
77
78
29.682
96
145
78
79
29.698
127
124
79
80
29.698
127
124
80
81
29.732
122
128
81
82
29.732
122
128
82
83
29.732
122
128
83
84
29.698
127
124
84
85
29.732
150
101
85
86
29.732
150
101
86
87
29.732
154
96
87
88
29.698
137
115
88
89
29.698
139
113
89
90
29.682
160
88
90
91
29.682
162
85
91
92
29.547
176
59
92
93
29.614
183
40
93
94
29.732
180
49
94
95
29.614
183
40
95
96
29.614
183
40
96
97
29.614
185
33
97
98
30
190
0
98
99
30
190
1
99
（1）血管半径为=29.6799μm。
（2）中轴线在XY平面投影曲线的拟合方程(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模型检验与分析
用多项式拟合得到重建后血管的三维空间曲线，再用平面去截空间曲线得到新的4层切片平面分别与题目中对应切片进行对比，在这里我们称的面积与的面积比值为切片重合度。

第01、25、50、75张图片分别与对应的新截取图片的对比图
...
图1 第01张切片 图11 第25张切片 图12 第50张切片 图13 第75张切片图
...
图14截后第01张切片图15截后第25张切片图16截后第50张切片图17截后第75张切片

根据上面曲线拟合得到的血管重新切片后，我们随机取出四张重切的图片，与原始的切片比较，得到两切片的重合度。
计算得：
第01张切片的重合度为99.20% ，
第25张切片的重合度为98.65% ，
第50张切片的重合度为97.30% ，
第75张切片的重合度为97.61% ，
计算平均重合度得到98.19%，说明本模型较为可靠。
7模型的优点与缺点
优点在于本模型的算法采用遍历搜索的方法得到最大内切圆的半径的精度较高，但是由于遍历搜索的运算次数异常庞大，导致模型的缺点是实践难度较大。从模型检验的数据可以看出计算的精度，但是整个运算在Pentium(R) Dual-Core处理器上耗费的时间已经超过120分钟。

参考文献
[1]刘会灯，朱 飞.MATLAB编程基础与典型应用[M] 2008年3月

附录
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% program 1 （找出切片的边界线的图）
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);
pic2(:,:,i+1)=edge(pic(:,:,i+1));
imshow(pic2(:,:,i+1));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% program 2 （求半径，找中轴线与切片交点）
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
lens(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
end
end
disp('This program has been done!!!')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% program 3 （多项式拟合空间中轴线）
format long
cenn=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)+cenn(i,2);
z=29.6799*cos(theta)+cenn(i,3);
hold on
surf(x,y,z);
axis equal;
end
shading flat
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% program 5 （重新切片）
s1=zeros(512);s2=1000*ones(512);pc=zeros(512,512,100);
for z=1:100
for i=1:512
for j=1:512
for k=1:100
s1(i,j)=sqrt((i-cenn(k,1))*(i-cenn(k,1))+(j-cenn(k,2))*(j-cenn(k,2))+(z-cenn(k,3))*(z-cenn(k,3)));
if s2(i,j)>s1(i,j)
s2(i,j)=s1(i,j);
end
end
if s2(i,j)-29.6799<0.00000001
plot(i,j);
hold on
axis equal
pc(i,j,z)=1;
end
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% program 6 （计算原切片与新切片的重合度p）
common=0;inner=0;p=0;
for i=1:512
for j=1:512
if a1(i,j)==0 && r1(i,j)==1
common=common+1;
end
if a1(i,j)==0
inner=inner+1;
end
end
end
p=common/inner;
% All programs is over.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%