血管的三维重建

血管的三维重建
1 摘要
序列图像的三维重建在各学科中都起到至关重要的作用，本次讨论的是血管的三维重建。

首先，假设该管道是由球心沿着某一曲面的球滚动包络而成，故本次的主要目的是求出中轴线坐标及半径。

现有100张平行切片图像，本次建立的模型可分为四步;第一步，采集图形边界点数据。

由于每张图片都是512*512的矩阵，故此数据很大，采用imread()函数将其读入矩阵A中。

第二步，最大内切圆寻找及半径的确定。

提出两种方案,分别是切线法和最大覆盖法; 从上述两种方法分析及考虑到我们所使用的工具和材料, 可以得出方法二更加直观, 计算机实现更容易, 计算复杂度更低, 所以我们采用后者。

根据以上算法, 我们抽取了所有的切片图进行半径的提取, 然后再求其平均值, 求其均值得到球的半径为29.6345。

第三步，轨迹的搜索。

在第二步中求出了血管的半径，轨迹的搜索就可以建立在半径确定的基础上, 当然我们也可以求出每一个切面图形的最大内切圆, 然后得到每个圆心的坐标, 即中轴线坐标, 但这样做计算机的运算量会很大, 同时由于最大内切圆搜索法的稳定性不高, 从而会造成搜索的不精确, 所以采用定半径搜索。

本文提出了三种方法, 分别为网格法、蒙特卡罗法和非线性规划法;本次采用非线性规划来实现。

第四步，绘制中轴线空间曲线图和在XOY、YOZ、XOZ 三个平面的投影图。

由定理1: 切片上血管截面图的头部顶点在XOY 平面上的投影点一定会落在
中
轴线在XOY 平面上的投影曲线上（在论文中以证明），并得出推论：切片上血管截面中中位线与中轴线在XOY 面上的投影重合。

最后可由中轴线和血管半径在作图软件中达到血管的三维重建，本次的模型还存在一定的不足，其假设为管道中轴线与每个切面有且只有一个交点，事实上还存在有多个交点的情况，但为了简化模型在此做了一定的假设，故会存在一定的误差。

关键词：三维重建内切圆半径轨迹（中轴线）
注：求边界时采用了老师的思想和程序。

2 问题重述
假设某些血管可视为一类特殊的管道，该管道的表面是由球心沿着某一曲线(称为中轴线)的球滚动包络而成。

例如圆柱就是这样一种管道，其中轴线为直线，由半径固定的球滚动包络形成。

现有某血管管道的相继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、YZ、ZX平面的投影图。

全部图象请从网上下载。

关于BMP图象格式可参考：
1．《VisualC+ +数字图象处理》第12页2.3.1节。

何斌等编著，人民邮电出版社，2001年4月。

2．http：／／／home／mxr／gfx／2d／BMP.txt
3 模型的假设
1.医学上, 血管不存在严重扭曲, 没有折皱。

2.血管可视为等径管道。

3.管道中轴线与每张切片有且只有一个交点。

4.切片间距以及图象象素的尺寸均为1。

5.对切片拍照的过程中不存在误差，数据误差仅与切片数字图象的分辨率有
关。

6．血管的表面是由半径固定的球的球心沿着某一曲线( 称为中轴线) 滚动包络而成。

7.切片足够薄, 其厚度对计算的影响可以忽略不计。

8.切片包含的其它圆的半径一定小于r。

9.不能被切片包含的圆的半径一定大于r。

4 符号说明
X 血管的最大内切圆的圆心的x的坐标值
Y 血管的最大内切圆的圆心的y的坐标值
Z 血管的最大内切圆的圆心的z的坐标值
R 血管最大内切圆的的半径
R血管最大内切圆整个圆周的离散坐标
ij
A 图像数据储存的三维矩阵
f R血管最大内切圆心匹配点数的函数
()
ij
5 分析与建立模型
5.1 分析
1）假定管道中轴线与每张切片有且只有一个交点, 球的半径固定, 切片的间距以及图像象素的尺寸均为1。

取坐标的z轴垂直于切片, 第一张切片记为平面z = 0, 第二张切片记为z = 1, 第100张切片的平面记为z = 99。

要求根据以上所给条件计算管道的中轴线与半径, 并绘制中轴线在X Y、YZ、ZX 平面的投影图。

5.2 几个重要的结论
对于某个切痕,可以作出它在空间中的结构图。

如图1所示,设某一切片的切痕与中轴线的交点为O,切下这一切痕的同时也截下中心在点O 的小球的一个大圆K,过点O 的血管的轴截面S也是同一个小球的一个大圆。

这两个大圆的交点即为小球的直径的两个端点,且这两个端点为切痕边界与大圆K 的两个切点,因此可得出以下结论:
1）每一个切痕中存在一个最大内切圆且最大内切圆的圆心在管道的中轴线上。

2）每个切痕的最大内切圆的半径均相等, 且等于管道的半径。

根据上面的结论,要求血管的中轴线,只需找出每一个切痕的最大内切圆的圆心。

图1 某个切痕的空间结构图
5.3 模型的建立
通过对血管的三维重现的分析可知, 我们的模型应由三部分组成:
5.3.1 采集图形边界点数据
边界点数据的采集方法有人工采集法和计算机采集法, 对于数据量较小的,
又具有某些特殊意义的点可以由人工采集, 若数量比较大且采集边界又非常复杂时, 则采取机器采集。

5.3.2 最大内切圆寻找及半径的确定
方法1 切线法
此方法可以从切面的外围轮廓线分析着手。

由题意分析知道, 所给的图片切面是由无数个球切面组成的。

而且外围轮廓线与最大内切圆有且仅有两个交点, 所以经过这两点的外围轮廓线的两条切线平行且间距最大。

基于上述分析, 我们可以通过找到这两条切线来找到最大内切圆的圆心及半径。

在实际操作中, 由于对图片的象素提取的离散性, 我们在计算导数时是用差分来代替。

方法2 最大覆盖法
最大覆盖法就是在切面中找到最佳的圆心位置和半径长度, 从而使得由这个圆心和半径所决定的圆面, 能最大面积地覆盖管道切面的图形, 这样搜索到的圆一定是最大内切圆, 这个圆的圆心就是我们所要找的球的球心, 这个圆的半径就是我们所要找的球的半径。

从上述两种方法分析及考虑到我们所使用的工具和材料, 可以得出方法二更加直观, 计算机实现更容易, 计算复杂度更低, 所以我们采用后者。

具体实现中, 我们先得到任意一张图片的象素矩阵, 然后将用于匹配的圆根据其圆心和半径将其圆周离散( 即以象素表示) , 并映射到
512⨯的图中, 其中圆周上的点为0, 其余的点为1, 即形成另一个象素矩512
512⨯两个矩阵在相同位置点上的值进行逻辑或运算, 如果其值为0, 阵。

这512
则为匹配点, 即此点在管道切面图形内, 否则其在切面图形外。

这样搜索到的匹配点最多同时半径最大的就是所要找的最大内切圆。

图2和图3 就是分别从图1. bmp 和图89. bmp 搜索到的最大匹配圆( 内部白色部位为管道切面图形, 深色圈为最大匹配圆) 。

根据以上算法, 我们抽取了所有的切片图进行半径的提取, 然后再求其平均值, 求其均值得到球的半径为29.6354。

图2 1. bmp 的最大匹配图图3 89. bmp 的最大匹配图
5.3.3 轨迹的搜索
5.3.3.1 目标函数的确立
在求出半径以后, 轨迹的搜索就可以建立在半径确定的基础上, 当然我们也可以求出每一个切面图形的最大内切圆, 然后得到每个圆心的坐标, 即中轴线坐标, 但这样做计算机的运算量会很大, 同时由于最大内切圆搜索法的稳
定性不高, 从而会造成搜索的不精确, 所以采用定半径搜索。

我们通过定圆( 半径为R ) 来找其中轴线, 也就是用定圆覆盖到切面图形上去, 找到匹配点数最多的一个位置, 从而得到此定圆圆心的位置。

具体实现时只要用定圆的圆心位置进行变化, 设其为(,)A x y , 则由A 点可以得到整个圆周的离散坐标ij R , 令
()ij f R 为匹配点数的函数, 其计算方法与最大内切圆求法相同, 即将ij R 根据圆
心A 的坐标和半径R ( 固定) 离散化形成矩阵后和切面图形矩阵作逻辑与操作, 从达到匹配最大, 则优化目标函数就是:
max[()]ij f R (512512)A ∈⨯坐标面上的点
5.3.3.2 目标函数要求求出一点圆心A , 使得定圆覆盖切面图形最大, 即多元函数极值最优解, 这是一个以X 和Y 为自变量二元函数, 这样可以通过以下三种方法来求得:
1) 用直接搜索求最优解( 网格法)
搜索过程中, 对每一个圆心的坐标X 和Y , 在其取值的范围内均取100 个步长, 分为1002个网格, 这样, 在一定的精度范围内, 可以求的一个较好的最优解。

2) 蒙特卡罗法
蒙特卡罗法, 也就是随机实验试点法。

它的基本思想是: 在函数的可行域内随机地选取实验点, 由于随机取得的点在区域中分配比较均匀, 所以对函数的大致形态能较好的体现。

模型中, 随机点用以下方法产生的。

0(10)(1)x x x x rand =+-⨯ 0(10)(1)y y y y rand =+-⨯
其中, ( x0, x1) 为X 的取值范围, ( y0, y1) 为Y 的取值范围。

3) 非线性规划
非线性规划即无约束优化,以数值迭代为基本思想,基本步骤为选取初值A(X0,Y0),进行k 次迭代并求出迭代解,由迭代解得到搜索方向和步长,如果k + 1次迭代符合给定的迭代终止条件,则停止迭代,得出最优解;否则继续迭代。

非线性规划的关键是搜索方向、步长和初值。

我们用拟牛顿法来选定搜索方向,拟牛顿法是在牛顿法基础上,克服牛顿法中黑赛阵不仅计算复杂、而且会出现变态、不正定等情况,同时保持了较快收敛的优点,从而得到最好的下降方向。

搜索步长的确定使用线性搜索的方法,或更为有效的插值方法。

由于此搜索程序的初值对于程序正确有效的搜索影响很大,又因为相邻的切片只有一个象素单位的距离, 可以认为中轴线的变化很小,所以我们可以把前一次确定的球心作为下一次搜索的起点,从而大大提高搜索的效率和准确性。

但由于转角处变化较大,在这种情况下搜索起点会郴够逼近而导致优化搜索的失效。

我们采用回溯技术来避免这种现象。

即当球心间距前后相差较大时( 采用工程上的观点, 以6 倍为基线),回退一次搜索,并以当前球心作为回退搜索的起点,这样就可以把各个球心间距变得比较匀称,从而相应的消除了转角上的搜索失效。

从上述对三种方法的分析可以得到网格法和蒙特卡罗法实现思路简单, 程序容易实现,但网格法搜索的精度不高, 误差较大,搜索时间长,同样蒙特卡罗法的实现对采点的数目要求很高,计算量大; 非线性规划法的实现复杂,但搜索速度快,计算量少,而且通过MATLAB 的优化工
具箱的函数可以很方便地实现,所以非线性规划法较优。

轨迹搜索结束以后,我们可以得到中轴线的100 个点,用所求半径的球定位在这些点上和轮廓线相交就可以得到图4。

从上图可以看到把图形重组后出现了不平滑点,也就是说直接用这些中轴线上的点构成的管道还不平滑,所以我们采用曲线拟合的方法来精确定位中轴线。

z
内边界重叠成的三维图像
x
y
图 4
5.3.4 绘制中轴线空间曲线图和在XOY 、YOZ 、XOZ 三个平面的投影图。

定理1: 切片上血管截面图的头部顶点在XOY 平面上的投影点一定会落在中 轴线在XOY 平面上的投影曲线上。

证明: 切片与XOY 平面平行, 点0, 点0c 分别为与切片相切的两个球的球心。

M, N 为两切点, 0 点投影在XOY 平面的投影点为P, 则OP 垂直于平面xoy, 又M 为切点, 所以OM 垂直于切面, 因为切面平行于XOY 平面, 所以OM 延长线垂直于XOY 平面于Pc, 又过平面外一点作平面的垂线有且只有一条。

故P 与Pc 重合。

即M 点与0 点在XOY 面上的投影点是重合的。

推广: 考虑切面与球面相交为一个圆的情况, Mc 为该截面的圆心, 如示意图2. 这时OM 仍垂直切面,Mc 与O 点在XOY 平面上的投影重合。

由此得推论: 切片上血管截面中中位线与中轴线在XOY 面上的投影重合。

6 模型的求解
6.1.1 导入血管片图像数据，存储为三维矩阵
利用Matlab 软件将100张切片图像数据储存为三维矩阵A 。

Matlab 程序如下： for i=0：1：99
Imname=sprintf （'%d.bmp',i ）;
A(:,:,i+1)=imread(imname);
End
由于ＢＭＰ格式文件在计算机中是以二进制数进行储存的，每张切片保存在一个二维的有０和１组成的矩阵中，其中０和１分别对应于图像中的黑象素和
512 个象白象素。

根据问题中的描述可知每一张ＢＭＰ格式的切片包含了512
素，每一个象素都有自己的一个确定的坐标。

在转换为矩阵存储后，Ａ（：，：，i+1）即代表了第ｉ张切片图像，此时象素坐标则对应地转换为矩阵的列与行。

因此，以下的求解过程中，均以象素在矩阵中的位置作为坐标，其中列对应于ｘ横坐标，行对应于ｙ坐标。

6.1.2 求出每张切片中血管图像的边界
首先，引进图像处理技术中的四邻域概念。

四邻域:某个像素的左、右、上、下四个象素称为该象素的四邻域。

如图5，象素E、S、W、N称为象素O的四邻域。

N
W E
O
S
图5 四邻域
60708090100110120130
图6 第0张切片中血管图像的内边界
求血管图像内边界的算法：逐点找出所有边界点坐标，即对图像进行逐行搜索，当遇到灰度值为0（黑）的象素点时，再搜索其四邻域的四个点，若在其四邻域中有一个象素的灰度值为1（白），则该点就是一个内边界点。

I=A（：，：，i）；
E = ones(size(I))；
for i = 1:512
p = find(I(i,:)==0)；
if ~isempty(p)
E(i,p(1)) = 0；
E(i,p(end)) = 0；
for j = 2:length(p)-1
if I(i-1,p(j))==1 || I(i,p(j)-1)==1 ||
I(i,p(j)+1)==1 || ... I(i+1,p(j))==1
E(i,p(j)) = 0；
end
end
end
end
imshow（E）
求血管图像外边界的算法：逐点找出所有边界点坐标，即对图像进行逐行搜索，当遇到灰度值为0（黑）的象素点时，再搜索其四邻域的四个点，若在其四邻域中任何一点象素的灰度值为1（白），则该四邻域中的这一点就是一个外边界点。

I=A（：，：，i）；
E = ones(size(I))；
for i = 1:512
p = find(I(i,:)==0);
if ~isempty(p)
for j = 1:length(p) if I(i-1,p(j))==1 E(i-1,p(j)) = 0; end
if I(i,p(j)-1)==1 E(i,p(j)-1) = 0; end
if I(i,p(j)+1)==1 E(i,p(j)+1) = 0; end
if I(i+1,p(j))==1 E(i+1,p(j)) = 0; end end end
end
最后将100张切片的血管图像边界叠加在一起，就得到了血管的三维叠加图，如图8。

.在上述的算法程序中，稍作修改便可以将每张切片的边界点保存在一个二维数组中，为下一步求解半径所用。

在程序中使用下列语句：
[R,C]=find(E==0); 其中R 表示行，Ｃ表示列，即矩阵E 的第R 行第Ｃ列是边界像素。

6.2.1球半径长度值的确定
球的半径通过对每一幅切片最大内切圆的搜索, 分别得到R1、R2、图7 管道中轴线拟合后的立体还原图,、R99、R100, 为减少搜索过程和离散化过程所引入的误差, 我们采用求均值的方法来得到球的半径, 可以得到半径R 为: 29.529。

6.2.2 血管管道的中轴线
管道的中轴线由轨道搜索而得，下图是由100点拟合而得的中轴线，表1中给出了直接搜索到的100 点中轴线在切面上的点的坐标(依次为x 、y 、z )
500
050
100
x
中轴线三维透视图
y
z
图 7
表一 中轴线在切面上的坐标 X
Y Z
X Y Z X Y Z X
Y Z
257 96 0 277 96 25 369 138 50 422 276 75 257 96 1 277 96 26 369 138 51 420 295 76 257 96 2 286 97 27 370 139 52 413 324 77 257 96 3 286 97 28 372 141 53 405 344 78 257 96 4 286 97 29 376 145 54 395 363 79 257 96 5 292 98 30 389 162 55 391 369 80 257 96 6 292 98 31 391 164 56 376 388 81 258 96 7 292 98 32 403 184 57 376 388 82 258 96 8 292 98 33 404 185 58 376 388 83 258 96 9 297 99 34 404 186 59 375 389 84 258 96 10 297 99 35 409 197 60 363 401 85 259 96 11 301 100 36 410 200 61 362 402 86 259 96 12 305 101 37 412 205 62 362 402 87 259 96 13 311 103 38 412 205 63 350 411 88 260 96 14 317 105 39 413 208 64 347 414 89 261 96 15 320 106 40 420 237 65 343 417 90 261 96 16 345 119 41 419 231 66 343 417 91 262 96 17 343 118 42 420 237 67 338 420 92 263 96 18 343 118 43 420 237 68 315 432 93 264 96 19 348 121 44 421 244 69 310 434 94 265 96 20 348 121 45 421 245 70 310 434 95 266 96 21 348 121 46
422 254 71 302 437 96 267 96 22 369 138 47 422 256 72 299 438 97 268 96 23 369 138 48 422 258 73 292 440 98 277 96
24
369
138
49
422
261
74
288
441
99
6.3绘制中轴线空间曲线图和在XOY 、YOZ 、XOZ 三个平面的投影图 中轴线在X Y 、YZ 、ZX 平面的投影图首先根据搜索到的中轴线在每一切面上的坐标投影到三个平面上, 然后根据投影点进行拟合，如下三图分别是投影点经过拟合得到的曲线。

250
300
350400
450
010*******
400
500x
y
中轴线在X Y 平面投影图
02040
6080100
250
300350400
450z
x
中轴线在ZX 平面投影图
20
40
60
80
100
010*******
400
500z
y
中轴线在ZY 平面投影图
6.4 血管的三维重建
这儿通过模型求解的中轴线的最佳拟合线得到管道的三维重建如图8，从图8可以发现此时的三维重建已经非常的好, 每个切面的外轮廓都能在重建图 中看到。

图 8 血管的三维重建图
7 模型的检验
将根据中轴线方程重构的血管对应位置的100张切片与原题给出的100张切片进行比较,对它们每对重合的像素坐标进行比较分析, 以此作为模型取舍的判别依据. 令 i U 为原血管第i 层切片像素坐标集 i M 与重组血管第i 层切片像素
坐标集i N 的交集:
i i
i U M N =
10()1100%()
n i i i P N a n P M -==⨯∑
101
0()1100%
()()()1
100%()n i i i n i i i i
P U b n P M P N P U c n P M -=-==⨯-=
⨯∑∑
其中100n =, ()()()i i i P M P N P U 、、分别为i i i M N U 、、 中的像素个数. a 为各重组切片像素个数与其对应原切片像素个数的百分比的平均值, b 为各重组
切片像素集与其对应原切片像素集的交集中像素个数与原切片像素个数的百分比的平均值, c 为各重组切片像素集不在其对应的原切片像素集内的像素个数与原切片像素个数的百分比的平均值。

基于 a 与b 在本模型中所起的作用的不对称, 因此我们给出一个赋权的判别函数- 拟合度函数:
(,,)(1)F a b t t a a t =-⨯+⨯
该函数反映了拟合的精确程度, 数值越大拟合得就越好。

后吻合程度不一, 前半段i a 大于100%,
为了防止重组切片与原切片在前
后半段i a 小于100% , 而又不能由拟合度函数显示出来, 定义了第2 个判别函数- 拟合一致性函数:
1()11n
i i i H a a n ==--∑
该函数反映了拟合的前后一致性, 数值越大拟合一致性越高。

经过上机编程, 计算出半径精确度为1 的最大内接圆的各重组切片像素个数与其对应原切片像素个数的百分比的平均值、各重组切片像素集与其对应原切片像素集交集的像素个数与原切片像素个数的百分比的平均值、各重组切片像素集不在其对应的原切片像素集内的像素个数与原切片像素个数的百分比的平均值分别81.165%,76.137%,51.28%a b c ===. 取权0.16t =, 则得拟合度函数:
(,,)(1)0.481.65%0.676.3%78.482%F a b t t a t b ==-⨯+⨯=⨯+⨯=
当把半径的输出精确度提高到0.01以后,
93.77%,91.22%,(,,)92.24%a b F a b t ===
由于球半径为29.7, 所以各切片是由球心在该切片内的球以及该切片上下方各30 个切片所对应的球与该切片所在平面相交而成. 当根据100张重组切片
的有序叠加来还原血管三维形态时, 第71 张以上的重组切片因为部分欠缺了第101张以上切片所对应的球, 使得第71张以上的切片被重组后像素上有欠缺, 有可能不能完全还原. 因此, 第71 张以上的重组切片作判别检验时可以不计, 以免影响判别结果. 但对于第31张以下的重组切片, 由于中轴线在该段几 乎与Z 轴平行, 因而重组切片受像素欠缺的影响不明显, 可忽略。

基于上面的说明, 另外再提高半径精度到0.01 这样可得到各指标如下( 具体计算过程略) :
95.26%,93.23%, 2.03%,(,,)94.042%a b c F a b t ====
(,,)F a b t 从78.482% 变成了94.042% 显示似合精度有了很大的提高。

以不同数量级作为搜索精度进行搜索, 并计算出如表一的指标值.随着半径精确度的不断提高, 拟合度函数的值不断增大, 表明所重组的管道与原管道拟合得越来越好. 随着半径精确度的不断提高,拟合一致性函数的值不断增加, 表明所重组由表1 可看出当精度为0.01 时, 似合度为94.04%, 拟合一致性指标为95.60% , 达到很高的精度, 由此说明模型达到了很好的效果.根据不同精度进行搜索, 对所有搜索到的最大内接圆的半径计算算术平均值可得精确度为
10.1.0、和01 时的平均半径分别为29.229.729.7、
和. 半径精确度为0.10.01和1时半径均近似为29.7. 因此, 取半径r =29.7。

8 模型的评价与推广
1 模型实际上是在综合考虑了编程的可实现性和运算量的大小后, 采用改
进的覆盖法来求出中轴线与每个截面的交点位置以及管道半径, 如果在程序中增大搜索范围则所得的结果会更为精确, 但会增加运算量。

另外, 我们还可以用收缩或扩大管道半径的方法来建立模型, 即对于一个初始给定的半径R0 对每一个截面图象Z=i , 判断这个半径为R0的圆是否可以完全含于截面图象中。

如果每张截面图象均能完全覆盖此圆, 则把球半径R0 增大1/ 2pix el, 否则将R0减小1/ 2pixel 。

直到对某一个R0,当R0增大1/ 2pixel 时至少有一个圆不能完全含于所有的截面图象中, 则将这个R0作为管道半径。

用这种方法建模也是比较合理的, 但计算量较大, 留待以后进行具体计算。

2 本模型是建立在条件假设l之上的, 即假设血管是一个光滑均匀的管道, 但实际情况并非如此, 血管的外壁一般有很多不规则的突起纹理的分支和其它附着物, 这就使得我们的模型对解决实际问题有一定的局限性。

另外当实际情况不满足假设3时，即管道中轴线与每张切片的交点不止一个(即存在多个圆心)，
n≥个交点，可取此时不可以用该模型。

对该种情况作简单的处理：假如有n()1
这n个交点的重心或选取与前一张所得圆心距离最近的圆心作为当前切片最大内切圆的圆心来进行模拟计算。

9 参考文献
[1] 汪国昭，陈凌钧.血管三维重建的问题[J].工程数学学报（建模专辑），
Vol.19，No.5,2002:54-58.
[2]中国大学生数学建模竞赛，http：///mcm01/problems.htm.
[3]徐晋，刘雪峰，柏荣刚.血管的三维重建[J]. 工程数学学报. Vol.19.No.5,2002:35-40.
[4] 顶峰平，周立峰，李孝朋.血管管道的三维重建[J]. 工程数学学报，Vol.19,No.5,2002. 47-35.
[5] 赵小健，陈立璋，吴小波，张传林. 血管的三维重建[J].暨南大学学报（自然科学版），Vol.24,No.5,2003:43-46.
10 附录
1求解出来的圆心坐标（X、Y、Z）及半径值R
X Y Z R
25796029.5
25796129.5
25796229.5
25796329.483
25796429.483
25796529.483
25796629.347
25896729.5
25896829.5
25896929.5
258961029.5
259961129.5
259961229.5
259961329.5
260961429.5
261961529.5
261961629.5
262961729.5
263961829.5 264961929.5 265962029.5 266962129.5 267962229.5 268962329.5 277962429.517 277962529.517 277962629.517 286972729.568 286972829.568 286972929.568 292983029.602 292983129.652 292983229.602 292983329.602 297993429.602 297993529.602 3011003629.602 3051013729.602 3111033829.619 3171053929.619
319.51064029.619 3451194129.577 342.5117.54229.586 342.5117.54329.586 3481214429.822 3481214529.822 3481214629.822 368.5137.54730.056 3691384830.056 3691384930.056 3691385030.056 3691385130.056 369.5138.55230.056 371.5140.55330.056 3761455429.822 389161.55529.611 390.5163.55629.611 4031845729.822 403.51855829.888 4041865929.619 4091976029.619 4102006129.602
4122056229.602 4122056329.602 4132086429.602 4202376529.652 4192316629.602 4202376729.568 4202376829.517 4212446929.5 4212457029.5 4222547129.5 4222567229.5 4222587329.5 4222617429.5 4222767529.5 4202957629.347 4133247729.602 4053447829.619 394.5362.57929.611 391368.58029.611 3763888129.822 3763888230.048 3763888330.048
3753898430.048 3634018530.048 3624028630.048 3624028729.822 3504118829.611 346.54148929.611 342.54179029.586 342.54179129.586 3384209229.518 3154329329.619 3104349429.619 3104349529.602 3024379629.602 2994389729.602 2924409829.602 2884419929.602
平均值29.63454。