http://www.ma-xy.com
目录
第一章 血管重建 1
1.1 题目要求 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 血管重建 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2.1 模型假设 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2.2 符号说明 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2.3 问题的分析 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2.4
模型的建立与求解
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2.5 程序 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.2.6 结果 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1
http://www.ma-xy.com
第一章 血管重建
1.1 题目要求
A 题：血管的三维重建
断面可用于了解生物组织、器官等的形态。例如，将样本染色后切成厚约 1µm 的切片，在
显微镜下观察该横断面的组织形态结构。如果用切片机连续不断地将样本切成数十、成百的平行
切片，可依次逐片观察。根据拍照并采样得到的平行切片数字图像，运用计算机可重建组织、器
官等准确的三维形态。
假设某些血管可视为一类特殊的管道，该管道的表面是由球心沿着某一曲线 (称为中轴线)
的球滚动包络而成。例如圆柱就是这样一种管道，其中轴线为直线，由半径固定的球滚动包络形
成。
现有某管道的相继 100 张平行切片图像，记录了管道与切片的交。图像文件名依次为 0.bmp、
1.bmp、…、99.bmp，格式均为 BMP，宽、高均为 512 个象素 (pixel)。为简化起见，假设：管道
中轴线与每张切片有且只有一个交点；球半径固定；切片间距以及图像象素的尺寸均为 1。
取坐标系的 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)
试计算管道的中轴线与半径，给出具体的算法，并绘制中轴线在 XY 、Y Z、ZX 平面的投
影图。图 (1.1) 是 100 张平行切片图像中的 5 张，全部图像请从网上 (http://mcm.edu.cn) 下载。
关于 BMP 图像格式可参考
¬
¬
《Visual C++ 数字图像处理》第 12 页 2.3.1 节。何斌等编著，人民邮电出版社，2001 年 4 月。
http://www.dcs.ed.ac.uk/home/mxr/gfx/2d/BMP.txt
1
http://www.ma-xy.com
1.2 血管重建 第一章 血管重建
图 1.1: 6 张切片示意图
1.2 血管重建
1.2.1 模型假设
1. 假设血管无严重扭曲；
2. 假设管道中轴线与每张切片有且只有一个交点；
3. 假设球半径固定；
4. 假设切片间距以及图像象素的尺寸均为 1；
5. 假设切片拍摄不存在误差，数据误差仅与切片数字图像的分辨率有关。
1.2.2 符号说明
1.2.3 问题的分析
计算管道的中轴线与半径，给出具体的算法，并绘制中轴线在 XY 、Y Z、ZX 平面的投影
图。
1.2.4 模型的建立与求解
取坐标系的 z 轴垂直于切片，第 1 张切片为平面 z = 0，第 100 张切片为平面 z = 99，并
以切片左顶点为坐标原点 o，z = 0 切片面为 xoy 面建立坐标系 oxyz，如图 (1.2) 所示
图 1.2: 血管重建坐标系示意图
如果我们想重构三维空间中的血管，我们就必须要知道滚动球的大小以及中轴线。就单一切
面而言，我们需要求解球心坐标以及球半径，我们对单一切面进行分析，如图 (1.3) 所示
http://www.ma-xy.com 2 http://www.ma-xy.com
http://www.ma-xy.com
第一章 血管重建 1.2 血管重建
图 1.3: 单一切面坐标系
我们共有 100 张切片，第 k 个切片的球中心为 (x
k
, y
k
)，球半径为 r
k
。下面，我们的主要工
作就是求解 x
k
, y
k
, r
k
，即切片的最大内切圆。对切片内部任意一个点 (x
i
, y
i
)，求出它到轮廓线
上所有点 (m
j
, n
j
) 的距离，并取其最小值，由于所有内点都对应一个最小值，在这些最小值中取
最大值，即为：最大内切圆的半径。假设切片内部有 p 个点，第 i 个点的坐标为 (x
i
, y
i
)，切片轮
廓线上有 q 个点，第 j 个点的坐标为 (m
j
, n
j
)，则切片内部第 i 点到轮廓线上第 j 点的距离为
d
i,j
=
√
(x
i
− m
j
)
2
+ (y
i
− n
j
)
2
第 i 个点对应的最小距离 r
i
为
r
i
= min
j
d
ij
其中：j = 1, 2, . . . , q。最大内切圆半径为
r
k
= max
i
r
i
其中：k = 1, 2, . . . , 100 为切片序号。将 100 个切片所对应的最大内切圆半径求取平均值，即为：
血管的半径
r =
1
100
100
∑
k=1
r
k
在求解出中轴线及球的大小之后，我们就可以重建血管了。我们用函数拟合的方法来求解中
轴线函数，并拟合中轴线在 3 个平面上的投影。这里我们不介绍函数拟合的方法。
1.2.5 程序
1 %%%%%%%%%%%%%%%%%%%%%血 管 重 建 主 程 序%%%%%%%%%%%%%%%%%%%%%%
2 %% %%%%%%%%%%%%%%%%%%求 解 最 大 内 切 圆%%%%%%%%%%%%%%%%%%
3 %[m, medg] = read_qiepian ; %读 取 切 片 数 据
4 cl c , c le ar
5 load ( ’ 100 pic . mat ’ )
6 %计 算 100 张 切 片 的 圆 心 坐 标 ， 将 其 记 录 在 j i e g u o 中 的 前 两 列 ， 半 径 记 录 在 j i e g u o 中 的 第 三 列
7 j i e g u o = ze r o s (100 , 3) ;
8 f o r k = 1:10 0
9 s e c t i o n = m( : , : , k) ;
http://www.ma-xy.com 3 http://www.ma-xy.com
http://www.ma-xy.com
1.2 血管重建 第一章 血管重建
10 [ x , y , r ] = MaxInCircle ( se c t io n ) ;% 求 解 每 个 切 片 的 内 切 圆 半 径 坐 标x^k , y^k ， 及 半 径 r^k
11 ji e g u o (k , 1) = x ;
12 ji e g u o (k , 2) = y ;
13 ji e g u o (k , 3) = r ;
14 end
15 save ( ’ ji e g u o . mat ’ , ’ j i e g u o ’ )
16 %% %%%%%%%%%%%%%%%%%%求 中 轴 线 的 拟 合 曲 线 及 其 在 三 个 坐 标 平 面 的 投 影%%%%%%%%%%%%%%%%%%
17 % 1. 绘 制 中 轴 线
18 f i g u r e
19 z = ( 1 :10 0 ) ;
20 plo t3 ( j i e g u o ( : , 1) ’ , ji e g u o ( : , 2) ’ , z , ’ ro ’ ) ;
21 hold on
22 % 拟 合 中 轴 线
23 p1 = p o l y f i t ( ( 1 : 100) ’ , j i e g u o ( : , 2) , 6) ;
24 p2 = p o l y f i t ( ( 1 : 100) ’ , j i e g u o ( : , 1) , 8) ;
25 x = p1 (1 ) * z.^6+p1 (2 ) * z.^5+p1 (3 ) * z.^4+p1 (4 ) *z.^3+p1 (5 ) *z.^2+p1 (6 ) *z+p1 ( 7 ) ;
26 y = p2 (1 ) * z.^8+p2 (2 ) * z.^7+p2 (3 ) * z.^6+p2 (4 ) *z.^5+p2 (5 ) *z.^4+p2 (6 ) *z.^3+p2 (7 ) * z.^2+p2 (8 ) * z+
p2 (9 ) ;
27 plo t3 ( y , x , z , ’b ’ )
28 gr i d on
29 x l a b e l x , y l a b e l y , zl a b e l z
30 legend ( ’ 中 轴 线 ’ , ’ 拟 合 中 轴 线 ’ )
31 t i t l e ( ’ 拟 合 中 轴 线 ’ )
32 hold o f f
33
34 % 2. 中 轴 线 在 xoz 面 上 的 投 影
35 f i g u r e
36 pl o t ( j i e g u o ( : , 2) ’ , z , ’ ro ’ ) ;
37 hold on
38 % 拟 合 中 轴 线 在 xoz 面 上 的 投 影
39 p = p o l y f i t ( ( 1 : 1 00 ) ’ , ji e guo ( : , 2) , 6) ; %以 6 次 多 项 式 拟 合 ， 横 坐 标 是 1 到 100
40 x = p(1 ) *z.^6+p(2) *z.^5+p( 3) *z.^4+p ( 4) *z.^3+p (5 ) *z.^2+p ( 6) *z+p ( 7) ;
41 pl o t (x , z , ’b ’ )
42 gr i d on
43 x l a b e l x , y l a b e l z
44 legend ( ’ xoz 上 的 投 影 ’ , ’ 拟 合 xoz 上 投 影 ’ )
45 t i t l e ( ’ 拟 合 中 轴 线 在 xoz 上 投 影 ’ )
46 hold o f f
47
48 % 3. 中 轴 线 在 yoz 面 上 的 投 影
49 f i g u r e
50 pl o t ( j i e g u o ( : , 1) ’ , z , ’ ro ’ ) ;
51 hold on
52 % 拟 合 中 轴 线 在 yoz 面 上 的 投 影
53 p = p o l y f i t ( ( 1 : 1 00 ) ’ , ji e guo ( : , 1) , 8) ; % 以 8 次 多 项 式 拟 合 ， 横 坐 标 是 1 到 100
54 y = p(1 ) *z.^8+p(2) *z.^7+p( 3) *z.^6+p ( 4) *z.^5+p (5 ) *z.^4+p ( 6) *z .^3+p ( 7) *z.^2+p ( 8) *z+p ( 9) ;
55 pl o t (y , z , ’b ’ )
56 x l a b e l y , y l a b e l z
57 legend ( ’ yoz 面 上 的 投 影 ’ , ’ 拟 合 yoz 面 上 的 投 影 ’ )
58 gr i d on
59 t i t l e ( ’ 拟 合 中 轴 线 在 yoz 上 投 影 ’ )
60 hold o f f
61
http://www.ma-xy.com 4 http://www.ma-xy.com
http://www.ma-xy.com
第一章 血管重建 1.2 血管重建
62 % 4. 中 轴 线 在 xoy 面 上 的 投 影
63 f i g u r e
64 pl o t ( j i e g u o ( : , 1) ’ , j i e g u o ( : , 2) ’ , ’ ro ’ ) ;
65 hold on
66 % 拟 合 中 轴 线 在 xoy 面 上 的 投 影
67 p1 = p o l y f i t ( ( 1 : 10 0) ’ , ji e g u o ( : , 2 ) , 6 ) ;
68 p2 = p o l y f i t ( ( 1 : 10 0) ’ , ji e g u o ( : , 1 ) , 8 ) ;
69 x = p1 (1 ) * z.^6+p1 (2 ) * z.^5+p1 (3 ) * z.^4+p1 (4 ) *z.^3+p1 (5 ) *z.^2+p1 (6 ) *z+p1 ( 7 ) ;
70 y = p2 (1 ) * z.^8+p2 (2 ) * z.^7+p2 (3 ) * z.^6+p2 (4 ) *z.^5+p2 (5 ) *z.^4+p2 (6 ) *z.^3+p2 (7 ) * z.^2+p2 (8 ) * z+
p2 (9 ) ;
71 pl o t (y , x , ’b ’ )
72 gr i d on
73 x l a b e l x , y l a b e l y
74 legend ( ’ xoy 面 上 的 投 影 ’ , ’ 拟 合 xoy 面 上 的 投 影 ’ )
75 t i t l e ( ’ 拟 合 中 轴 线 在 xoy 上 投 影 ’ )
76 %%%%%%%%%%%%%%%%%%带 有 拟 合 中 轴 线 的 轮 廓 切 面 图 空 间 效 果 图%%%%%%%%%%%%%%%%%%
77 c l c ; c l e ar ; c l o s e a l l ;
78 load ( ’ 100 pic . mat ’ ) ; load ( ’ ji e g u o . mat ’ ) ;
79 n=100;m1=s i z e (m, 1 ) ;m2=s i z e (m, 2 ) ;
80 n=zer o s (m1,m2, n) ;
81 f o r k=0:99
82 n ( : , : , k+1)=edge (m( : , : , k+1) ) ;
83 end
84 f o r k=0:5:99
85 f o r i =1:2:512
86 f o r j =1:2:512
87 i f (n( i , j , k+1)==1)
88 pl ot 3 ( i −257,j − 257,k+1, ’b . ’ ) ; hold on
89 end , end , end , end
90 hold on
91
92 p1=p o l y f i t (( 1 : 1 0 0 ) ’ , j i e g u o ( : , 2 ) , 6) ;
93 p2=p o l y f i t (( 1 : 1 0 0 ) ’ , j i e g u o ( : , 1 ) , 8) ;
94 z =(1:100) ;
95 x=p1 ( 1) *z.^6+p1 ( 2) *z.^5+p1 (3 ) * z.^4+p1 (4 ) * z.^3+p1 (5 ) * z.^2+p1 (6 ) * z+p1 ( 7 ) ;
96 y=p2 ( 1) *z.^8+p2 ( 2) *z.^7+p2 (3 ) * z.^6+p2 (4 ) * z.^5+p2 (5 ) * z.^4+p2 (6 ) * z.^3+p2 (7 ) *z.^2+p2 (8 ) *z+p2
(9 ) ;
97 plo t3 ( y , x , z , ’ r ’ )
98 gr i d on
99 t i t l e ( ’ 拟 合 中 轴 线 的 轮 廓 切 面 图 ’ )
100 hold o f f
101 %%%%%%%%%%%%%%%%%%%计 算 平 均 半 径%%%%%%%%%%%%%%%%%%
102 load ( ’ ji e g u o . mat ’ ) ;
103 r = mean( ji e g u o ( : , 3 ) ) ;
104 %%%%%%%%%%%%%%%%%%函 数 对 管 道 表 面 进 行 三 维 重 建 ( 透 明 化 处 理 )%%%%%%%%%%%%%%%%%%
105 c l c ; c l e ar ; c l o s e a l l ;
106 load ( ’ 100 pi c . mat ’ ) ; load ( ’ jieg u o .mat ’ ) ;
107 n=100;m1=s i z e (m,1 ) ;m2=s i z e (m,2 ) ;
108 p1=p o l y f i t (( 1 : 1 0 0 ) ’ , j i e g u o ( : , 2 ) , 6) ;
109 p2=p o l y f i t (( 1 : 1 0 0 ) ’ , j i e g u o ( : , 1 ) , 8) ;
110 z =(1:100) ;
111 x=p1 ( 1) *z.^6+p1 ( 2) *z.^5+p1 (3 ) * z.^4+p1 (4 ) * z.^3+p1 (5 ) * z.^2+p1 (6 ) * z+p1 ( 7 ) ;
http://www.ma-xy.com 5 http://www.ma-xy.com
http://www.ma-xy.com
1.2 血管重建 第一章 血管重建
112 y=p2 ( 1) *z.^8+p2 ( 2) *z.^7+p2 (3 ) * z.^6+p2 (4 ) * z.^5+p2 (5 ) * z.^4+p2 (6 ) * z.^3+p2 (7 ) *z.^2+p2 (8 ) *z+p2
(9 ) ;
113 r=mean( ji e g u o ( : ,3 ) ) ;
114 X= [ ] ;Y= [ ] ;Z = [ ] ;
115 f o r t = 0: 0.5:4 * pi ;
116 xx=cos ( t ) ;
117 yy=s i n ( t ) ;
118 x1=r*xx+x ;
119 y1=r*yy+y ;
120 z1=z ;
121 X=[X; x1 ] ;
122 Y=[Y; y1 ] ;
123 Z=[Z ; z1 ] ;
124 end
125 f i g u r e
126 s u r f (X,Y, Z)
127 shading f l a t
128 aa =0.5;
129 alpha ( aa) ;
130 hold on
131 p lo t3 (x , y , z , ’ r ’ )
132 hold o f f
133 gr i d on
134 t i t l e ( ’ 使 用 s ur f 曲 面 函 数 按 照 拟 合 中 轴 线 和 半 径 模 拟 血 管 ’ )
135
求解切片圆心坐标及圆半径的程序 MaxInCircle 如下
1 f u nc t io n [ x , y , r ] = MaxInCircle ( s e c t i o n )
2 % 此 函 数 用 与 求 解 单 一 切 片 的 最 大 内 切 圆 圆 心 坐 标 x^k , y^k 及 半 径 r^k
3 %%%%%%输 入%%%%%
4 % se c t i o n ： 单 一 切 片；
5 %%%%%%输 出%%%%%
6 % x ： 圆 心 横 坐 标 ；
7 % y ： 圆 心 纵 坐 标 ；
8 % r ： 圆 半 径；
9 s e c t i o n = 1 − s e c t i o n ;
10 BW = edge ( se ct i on , ’ s obe l ’ ) ;
11 BW2 = bwmorph( s e ct i on , ’ s k e l ’ , i n f ) ;
12 [ x , y , ~ ] = f i n d (BW2) ; %轮 廓 内 点
13 [m, n , ~ ] = f i n d (BW) ; %轮 廓 边 界 点
14 d = ze r o s ( lengt h (x) , le ngt h (m) ) ;
15 r = z e ros ( l ength (m) , 1) ;
16 f o r i = 1 : len gth ( x )
17 f o r j = 1 : lengt h (m)
18 d( i , j )=s q rt ( ( x ( i )−m( j ) ) ^2+(y ( i )−n( j ) ) ^2) ;
19 end
20 r ( i ) = min(d( i , : ) ) ; % 最 小 半 径
21 end
22 [ r , ind ] = max( r ) ; % 内 切 圆 半 径
23 x = x( ind ) − 256; % 内 切 圆 横 坐 标
24 y = y( ind ) − 256; % 内 切 圆 纵 坐 标
25
http://www.ma-xy.com 6 http://www.ma-xy.com
