血管的三维重建.docx
《血管的三维重建.docx》由会员分享，可在线阅读，更多相关《血管的三维重建.docx（11页珍藏版）》请在冰豆网上搜索。
血管的三维重建
血管的三维重建问题
摘要
越来越多的在医学领域及生物领域需要大量的切片工作，来研究切片表面的组织形态，因而如何能将这些切片在利用计算机进行复原，显得很重要。
在对管道的半径的求解中，我们采用的是枚举法，即逐个将100切片的最大内切圆的半径都求出来，最后取其平均值，即为我们所求的管道的半径。
在半径求出的同时把最大内切圆的圆心也求出来。
要求得切片的最大内切圆，首先对切片图像0-1矩阵转化，然后利用edge函数找到图形的轮廓，利用bwmorph函数求得图形的轮廓，然后利用find函数求出轮廓与骨架上的坐标。
利用两点间的距离公式，我们就能求出骨架上每个点到轮廓的距离，在所有点到轮廓的最短距离中的那个最大值，就是我们要求的最大内切圆半径。
此圆心即为中轴线上的点。
在中轴线与中轴线在XY，YZ，ZX平面的投影的计算中，我们主要采用多项式拟合polyfit的方法，从得出拟合曲线可以看出与实际图形很接近。
最后我们拟合出了血管的三维形态。
关键词：
最大内切圆半径中轴线多项式拟合
一．问题重述
如今在医学领域，生物学领域等都需要了解生物组织、器官等的断面形态，以便于能够准确的研究该样品，得出具有说服力的研究成果。
例如，将样本染色后切成厚约1μm的切片，在显微镜下观察该横断面的组织形态结构。
如果用切片机连续不断地将样本切成数十、成百的平行切片，可依次逐片观察。
这时候我们如果我们需要知道血管的三维形态，根据拍照并采样得到的平行切片数字图象，运用计算机就可重建组织、器官等的三维形态。
假设某些血管可视为一类特殊的管道，该管道的表面是由球心沿着某一曲线（称为中轴线）的球滚动包络而成。
本文给出了某管道的相继100张平行切片图象，记录了管道与切片的交。
图象文件名依次为0.bmp、1.bmp、…、99.bmp，格式均为BMP，宽、高均为512个象素（pixel）。
我们要做的就是根据这些切片计算管道的中轴线与半径，给出具体的算法，并绘制中轴线在XY、YZ、ZX平面的投影图。
二、问题分析
本题是一个求最优解的问题，通过具体的算法得出最接近实际血管的中轴线与半径，并能较准确的绘制出中轴线在XY、YZ、ZX平面的投影图，得出直观的图形。
需要处理大量的图片，得出大量的数据信息，所以为处理该问题的程序中需要大量的循环，所以程序的完成与运行皆很复杂。
本题管道为由球心沿着某一曲线的球滚动包络而成，且管道中轴线与每张切片被假设有且只有一个交点，则每张切片一定包含球的最大截面，即过球心的圆面，则每张切片的最大内切圆就是过球心的圆面。
由于本题中取坐标系的Z轴垂直于切片，第1张切片为平面Z=0，第100张切片为平面Z=99。
Z=z切片图象中象素的坐标依它们在文件中出现的前后次序为
（-256，-256，z），（-256，-255，z），…（-256，255，z），
（-255，-256，z），（-255，-255，z），…（-255，255，z），
……
（255，-256，z），（255，-255，z），…（255，255，z）。
则每个切片的最大内切圆的圆心的z坐标已给，再加上x,y,则中心轴就是每个切片的最大内切圆的圆心（x,y,z）连成的光滑曲线，而管道的半径就是这个球的半径，即最大内切圆的半径。
进而再求出中轴线在XY、YZ、ZX平面的投影图。
那么这道题的问题最后就转化为求每张切片的最大内切圆的圆心与半径问题。
三、模型假设
1、管道的表面是由球心沿着某一曲线（称为中轴线）的球滚动包络而成。
2、管道中轴线与每张切片有且只有一个交点。
3、球半径固定。
4、切片间距以及图象象素的尺寸均为1。
四、符号说明
一张切片内骨架点i到边缘点ｊ的距离
第i张图的最大内切圆半径
半径的平均值
五、模型的建立与求解
问题一：
求管道的半径与每张切片的最大内切圆的圆心
１.要求得管道的半径，首先需要求出每个切片的最大内切圆，为了得出比较精确的结果，我们采用枚举法。
对与100张切片，考虑到求得的最大内切圆半径由于误差的存在肯定会有所不同，我们将所有切片的最大内切圆半径求出后，然后取其平均值。
２.为得到每张图片的最大内切圆半径，首先先将二值图片转换成0-1矩阵，矩阵横纵坐标对应原图像的直角坐标系位置（1代表黑像素，0代表白像素），为了后面find函数寻找矩阵中为1的坐标并记录。
通过相关函数找出轮廓（edge）和骨架（bwmorph）记录轮廓q（i,j）和骨架p（i,j）的坐标，
利用两点间的公式去求骨架上的点到边界上所有点的最小距离即骨架上此点的内切圆的半径。
（1）
骨架上所有点的最大内切圆即为该切片的最大内切圆，
（2）
（这里i代表第i张切片）
程序见附录1
其半径即为最大内切圆半径。
通过循环找出最大半径并记录该半径所对的坐标点。
基于此我们能求出所有切片的最大内切圆的半径与圆心，数据见附录5
通过求得的100张切片的最大内切圆的半径我们取其平均值得出更加精确的管道的半径
。
问题二：
求中心轴与中心轴在XY，YZ，ZX平面的投影
（x,y,z）为切片的最大内切圆的圆心坐标
1.利用问题一中求得数据（x,y,z）,然后根据多项式拟合polyfit得出x,y,z的拟合方程为：
中轴线在Z-Y平面的投影拟合多项式：
中轴线在X-Y平面的投影拟合多项式：
中轴线在Z-X平面的投影拟合多项式：
程序见附录2
得到的拟合曲线如下：
图1为中轴线在XY平面的投影，曲线代表多项式拟合的结果，折线代表真实数据的图形
图2为中轴线在ZY平面的投影，曲线代表多项式拟合的结果，折线代表真实数据的图形
图３为中轴线在ZX平面的投影图，曲线代表多项式拟合的结果，折线代表真实数据的图形
从以上拟合图形中我们可以看出，通过多项式拟合的投影图与通过坐标点得出的投影图基本相同，所以可以近似认为拟合的投影图即为我们所要求得的中轴线在XY，YZ，ZX平面的投影。
2.同样利用多项式拟合polyfit得出管道的中轴线在三维坐标中的拟合曲线如下：
（程序见附录3）
图４为中轴线在三维坐标的立体图
3.利用球坐标以及以上求出的最大半径、以及拟合时在z时的坐标，x1,y1利用plot3，程序如附录四：
图５血管的三维重建
六、模型评价及总结
解决本题的首要前提是要理解像素的概念，在解此题的最初由于对图片的像素没有明确的理解而浪费了很长时间。
在我们的模型中还有尚有不完美的地方，在求切片的最大内切圆半径与圆心的时候，没有编写出能够一下就将所有切片的数据都求出的程序，而是仅仅编写出一个图片的程序，然后逐一求出，在这个问题上花费的时间长，且过程繁琐。
但是在进行曲线拟合时，由于运用ｈｅｌｐ知道了相关函数的用法，运用polyfit进行拟合，很大程度上节省了时间，拟合出的相对来说比较满意。
立体切片的三维重建在医学与生物学等领域的作用很大，因为这些领域往往会研究一些物体切片的组织结构，所以此问题的实现能够解决很多现实问题。
为一些研究领域提供很大的帮助。
七、参考文献
【1】董辰辉，MATLAB2008全程指南，电子工业出版社。
【2】张红云等，图像骨架的提取的应用，2010年第四期。
【3】
附录
1.jieguo=zeros（100,4）;
J0=imread（'E:
\99.BMP'）;
fori=1:
512;
forj=1:
512;
j0（i,j）=1-J0（i,j）;转化为黑色为1，白色为0，为了后面find函数寻找矩阵中为1的坐标并记录
end
end
lk=edge（j0,'sobel'）;
gj=bwmorph（j0,'skel',inf）;
[x0,y0,v0]=find（lk）;找到轮廓灰色区域
[a0,b0,c0]=find（gj）;找到骨架灰色区域
m=length（a0）;骨架灰色区域个数
n=length（x0）;轮廓灰色区域个数
jl=zeros（m,n）;建立0矩阵为求内切圆半径
cf=zeros（m,2）;建立两列0矩阵为存放中心点坐标
fori=1:
m
forj=1:
n
p1=a0（i）;
q1=b0（i）;
p2=x0（j）;
q2=y0（j）;骨架、轮廓坐标赋值
jl（i,j）=sqrt（（p1-p2）^2+（q1-q2）^2）;
end
[zx,zxxh]=min（jl（i,:
））;骨架上一点到轮廓的最短距离，即骨架上各个点的内切圆的半径
cf（i,1）=zx;
cf（i,2）=zxxh;
end
[zd,zdxh]=max（cf（:
1））;找到其中最大的半径，并把半径和圆心坐标存储
x=a0（zdxh）-256;与题目所给坐标轴对应
y=b0（zdxh）-256;与题目所给坐标轴对应
jieguo（k+1,1）=x;x轴坐标
jieguo（k+1,2）=y;
jieguo（k+1,3）=k+1;
jieguo（k+1,4）=zd;
2.polyfit（Z,X,7）
plot（Z,X,0:
10:
99,polyval（ans,0:
10:
99））
[p,s]=polyfit（Z,X,7）
polyfit（Z,Y,7）
plot（Z,Y,0:
10:
99,polyval（ans,0:
10:
99））
[p,s]=polyfit（Z,Y,7）
plot（X,Y,-160:
10:
188,polyval（ans,-160:
10:
188））
3.formatlong
px=polyfit（z,x,7）;
x1=polyval（px,z）;取拟合后再z点时的点
py=polyfit（z,y,7）;
y1=polyval（py,z）;
figure
（1）;
plot3（x1,y1,z,'r'）
４.t=linspace（0,pi,25）;存储平均分pi为25份存储
p=linspace（0,2*pi,25）;存储平均分2pi为25份存储
[theta,phi]=meshgrid（t,p）;存储角度
fori=1:
100
ceen（:
1）=x1;ceen（:
2）=y1;ceen（:
3）=z;
x=29.49288*sin（theta）.*sin（phi）+ceen（i,1）;
y=29.49288*sin（theta）.*cos（phi）+ceen（i,2）;
z1=29.49288*cos（theta）+ceen（i,3）;
holdon使图像不被覆盖
plot3（x,y,z1）;
axisequal;
holdoff
end
5、
切片号
半径
x坐标
y坐标
z坐标
切片号
半径
x坐标
y坐标
z坐标
1
29.069
-160
0
0
51
30.414
-114
117
50
2
28.284
-160
2
1
52
30.414
-114
117
51
3
29
-160
1
2
53
30.414
-113
118
52
4
29.069
-160
1
3
54
30.414
-112
119
53
5
29.069
-160
2
4
55
30
-111
120
54
6
29.069
-160
2
5
56
29.732
-111
120
55
7
29
-160
1
6
57
29.698
-111
120
56
8
29.017
-160
4
7
58
29.698
-111
120
57
9
29
-160
1
8
59
29.53
-81
142
58
10
28.862
-160
1
9
60
29.547
-51
156
59
11
28.862
-160
7
10
61
29.547
-51
156
60
12
28.862
-160
8
11
62
29.614
-31
162
61
13
28.862
-160
9
12
63
29.614
-31
162
62
14
29.017
-160
10
13
64
29.614
-31
162
63
15
29.017
-160
12
14
65
29.614
-35
161
64
16
29.017
-160
13
15
66
29.614
-35
161
65
17
29.017
-160
14
16
67
29.428
-26
163
66
18
29.017
-160
16
17
68
29.411
-35
161
67
19
29.017
-160
17
18
69
29.275
-26
163
68
20
29.017
-160
18
19
70
29.428
46
163
69
21
29.017
-160
19
20
71
29.614
46
163
70
22
29.017
-160
20
21
72
29.614
46
163
71
23
29.017
-160
21
22
73
29.614
46
163
72
24
29.017
-160
22
23
74
29.614
65
158
73
25
29.069
-160
21
24
75
29.732
68
157
74
26
29.069
-160
21
25
76
29.732
65
158
75
27
29.069
-160
21
26
77
29.547
81
152
76
28
29.155
-159
30
27
78
29.53
81
152
77
29
29.275
-159
30
28
79
29.53
81
152
78
30
29.275
-159
29
29
80
29.698
135
117
79
31
29.428
-158
35
30
81
29.698
136
116
80
32
29.614
-157
40
31
82
29.698
136
116
81
33
29.614
-157
40
32
83
29.698
137
115
82
34
29.614
-157
40
33
84
29.698
138
114
83
35
29.614
-156
44
34
85
29.698
138
114
84
36
29.732
-153
55
35
86
29.698
139
113
85
37
29.732
-153
55
36
87
29.698
139
113
86
38
29.732
-153
55
37
88
29.698
139
113
87
39
29.732
-152
58
38
89
29.698
140
112
88
40
29.614
-152
58
39
90
29.698
140
112
89
41
29.547
-150
63
40
91
29.53
172
67
90
42
29.682
-135
92
41
92
29.53
172
67
91
43
29.53
-148
68
42
93
29.53
172
67
92
44
29.682
-119
112
43
94
29.53
172
67
93
45
29.698
-118
113
44
95
29.732
182
43
94
46
29.698
-117
114
45
96
29.614
187
24
95
47
29.698
-117
114
46
97
29.614
187
24
96
48
30.414
-116
115
47
98
29.614
187
24
97
49
30.414
-115
116
48
99
29.614
187
24
98
50
30.414
-115
116
49
100
29.428
188
18
99
升级会员 