1、血管的三维重建摘要本文探讨血管的三维重建，由血管的相继 100 张平行切片图像计算血管的中轴线与半径，并绘制血管在三个坐标平面上的投影。由于血管的表面是由球心沿着某一曲线（称为中轴线）的球滚动包络而成，由此我们得出结论：每个切片一定包含滚动球的大圆，并且它一定为切片的最大内切圆，而最大圆所对应的半径即为血管半径，所以求血管半径就转化为求每一个切片内部的点到切片外部轮廓的最大半径。首先，读取 100 张血管切面图，把它们转换成 Logical 矩阵，从中提取切片截面轮廓点构成一个新的矩阵。然后找到原图片矩阵中像素点的内点（切片图片中轮廓线中的点），从而得到内点到切片轮廓点的最小距离矩阵和最小距离
2、中的最大值矩阵，最大值即为血管半径。最后计算所有切片的血管半径，并对这些半径求平均值，得到平均血管半径为：29.6799m 。由 100 张切片的最大内切圆圆心坐标拟合得出中轴线方程以及其在三个坐标平面内的投影曲线方程。由中轴线得到血管的三维立体重建图，用平面 去截血)74,920(,iZ管的三维立体重建图，得到新的 4 张截面图。把它们分别与题设中的对应截面进行内点个数对比。我们定义两张切片所共同拥有的内点个数与原切片内点个数的比值为重合度。计算得到平均重合度为：98.19% 。关键词：血管半径 中轴线 切片 重建（来自作者：欢迎各界人士批评指正， 学术交流邮箱 。文章作于 2011 年 8
3、 月 10 日，陕西科技大学理学院实验室）1 问题的重述断面可用于了解生物组织、器官等的形态。例如，将样本染色后切成厚约1m 的切片，在显微镜下观察该横断面的组织形态结构。如果用切片机连续不断地将样本切成数十、成百的平行切片， 可依次逐片观察。根据拍照并采样得到的平行切片数字图象，运用计算机可重建组织、器官等准确的三维形态。假设某些血管可视为一类特殊的管道，该管道的表面是由球心沿着某一曲线（称为中轴线）的球滚动包络而成。例如圆柱就是这样一种管道，其中轴线为直线，由半径固定的球滚动包络形成。现有某管道的相继 100 张平行切片图像，记录了管道与切片的交。图像文件名依次为 0.bmp、1.bmp、
4、 99.bmp，格式均为 BMP，宽、高均为 512个象素（pixel）。为简化起见，假设：管道中轴线与每张切片有且只有一个交点；球半径固定；切片间距以及图像象素的尺寸均为 1。计算管道的中轴线与半径，给出具体的算法，并绘制中轴线在 XY、YZ、ZX平面的投影图。 2 模型假设 (1) 假设血管中轴线与每张切片有且只有一个交点。(2) 假设血管的半径固定不变。(3) 假设在对切片拍照的过程中不存在误差。(4) 假设血管不存在严重扭曲。3 符号说明第 张切片中轮廓线的最大内切圆圆心kCk第 张切片中轮廓线的最大内切圆半径R第 张切片中第 个内点与第 个轮廓点的距离ijkdkij100 张切片图片
5、转换以后的三维 0-1 矩阵1052PA100 张切片图片的轮廓线生成的矩阵B100 张切片轮廓的最大内切圆圆心坐标310cen平均血管半径R第 张切片图片中所有内点的集合kk第 张切片图片中轮廓点的集合kk原切片图片的上内点及边界点的集合S重新切片得到的内点及边界点的集合TR-square 多项式拟合的指标数4 问题的分析根据题目整个管道是由球心沿着某一曲线（称为中轴线）的球滚动包络而成。基于几何原理可以得出以下两条定理。（1）球的任意截面都是圆。（2）经过球心的球截面是所有的截面圆当中半径最大的圆。基于上述两个定理可以得出：每张切片的最大内切圆的圆心位于血管的中轴线上，该圆的半径等于血管半
6、径。分别求出 100 张切片轮廓线的最大内切圆圆心 与半径 。kCkR对所有最大内切圆的圆心 做拟合，即可求出血管的中轴线。kC为减少误差，我们用求均值的方法 来得到平均血管半径 。1010)(kRR根据拟合出的中轴线方程，即可求出中轴线在 XY、YZ、ZX 平面上的投影。5 模型的建立与求解5.1 模型建立5.1.1 求血管的半径（1）导入数据，转换存储方式为了计算方便，利用 MATLAB1软件将 100 张 BMP 文件转换存储为一个三维 0-1 矩阵 （ 0 代表黑色像素点，1 代表白色像素点）。512PA（2）求切片轮廓线上各点坐标利用 MATLAB 软件的内部函数 edge，求得所有
7、切片图片的轮廓线生成的矩阵 1052PB（3）求切片轮廓线的最大内切圆的半径首先求出每个内点距轮廓线的距离，取其中的最小值，即为以这些内点为圆心的轮廓线的内切圆半径；其次找出所有内点确定的最小内切圆中半径最大的一个，即为轮廓线的最大内切圆半径。具体步骤：Step1 在三维 0-1 矩阵 中遍历搜索内点,将其存储于集合 中 ； 1052PA)(kStep2 在轮廓线矩阵 中遍历搜索轮廓点，将其存储于集合 中； B )(Step3 计算集合 中第 个内点 到集合 中第 个轮廓点 的距离 )(kiyx,)(kjqpB,，22qypxdijk首先取每个内点到轮廓点的最小距离，即为以这些内点为圆心的内切
8、圆半径；其次找出所有内点确定的最小内切圆中半径最大的一个，即为轮 廓线的最大内切圆半径 minaxijkjikdRStep4 计算平均血管半径 1010)(kR（4）求切片轮廓线的最大内切圆心坐标将 100 张轮廓线的最大内切圆半径所对应的圆心坐标 记录于矩阵),(wzCk中。)3,10(cen)10,.2(,)3,(2,)1,(kcenwzkk5.1.2 由圆心坐标求中轴线方程及曲线投影（1）利用 MATLAB 画图工具由 100 组圆心坐标画出中轴线以及中轴线在XY、YZ、ZX 面的投影，由于图像文件以象素为单位存储数据，这使得图像的原始数据是离散化的。对于离散的位图来说，将变得不可分辨，
9、可能切片轮廓线的最大圆不止一个。为了减少误差，将中轴线在 XY、YZ、ZX 面的曲线方程进行多项式拟合，建立拟合方程： )(1XFY2Z)(3（2）由建立的投影曲面得到中轴线的方程 );(;213TFZYX5.2 模型求解根据所建模型求出 100 张切片的最大内切圆的圆心坐标、半径以及血管半径表 1 最大内切圆半径及圆心坐标最大内切圆圆心坐标切片序号最大内切圆半径（ R/m）横 坐标(X/m)纵 坐标(Y/m)竖 坐标（Z/m）0 29.069 -160 1 01 29.069 -160 1 12 29 -160 2 23 29.069 -160 2 34 29.069 -160 2 45 2
10、9.069 -160 2 56 29.155 -160 2 67 29.155 -160 2 78 29.275 -160 2 89 29.275 -160 3 910 29.428 -160 3 1011 29.428 -160 4 1112 29.614 -160 4 1213 29.614 -160 5 1314 29.833 -160 6 1415 29.833 -160 7 1516 29.833 -160 8 1617 30 -160 9 1718 30 -160 10 1819 30 -160 11 1920 30 -160 12 2021 30 -160 13 2122 30
11、-160 14 2223 30 -160 15 2324 30 -160 15 2425 29.833 -160 16 2526 29.614 -159 27 2627 29.614 -159 27 2728 29.614 -159 28 2829 29.614 -158 34 2930 29.614 -158 34 3031 29.614 -157 40 3132 29.614 -155 48 3233 29.614 -155 48 3334 29.614 -155 48 3435 29.732 -153 55 3536 29.732 -152 58 3637 29.732 -152 58
12、3738 29.732 -152 58 3839 29.614 -152 58 3940 29.547 -148 68 4041 29.547 -149 66 4142 29.53 -140 84 4243 29.53 -140 84 4344 29.682 -135 92 4445 29.682 -135 92 4546 29.682 -135 92 4647 29.698 -115 116 4748 29.698 -113 118 4849 29.698 -112 119 4950 29.698 -111 120 5051 29.698 -111 120 5152 29.698 -111
13、120 5253 29.698 -111 120 5354 29.682 -111 120 5455 29.53 -66 150 5556 29.53 -59 153 5657 29.732 -70 148 5758 29.967 -70 148 5859 29.732 -54 155 5960 29.732 -54 155 6061 29.614 -31 162 6162 29.614 -31 162 6263 29.614 -6 166 6364 29.833 -6 166 6465 29.833 8 167 6566 30 10 167 6667 30 11 167 6768 30 12
14、 167 6869 30 12 167 6970 30 12 167 7071 29.833 26 166 7172 29.614 26 166 7273 29.614 46 163 7374 29.732 65 158 7475 29.732 65 158 7576 29.732 71 156 7677 29.682 96 145 7778 29.682 96 145 7879 29.698 127 124 7980 29.698 127 124 8081 29.732 122 128 8182 29.732 122 128 8283 29.732 122 128 8384 29.698 1
15、27 124 8485 29.732 150 101 8586 29.732 150 101 8687 29.732 154 96 8788 29.698 137 115 8889 29.698 139 113 8990 29.682 160 88 9091 29.682 162 85 9192 29.547 176 59 9293 29.614 183 40 9394 29.732 180 49 9495 29.614 183 40 9596 29.614 183 40 9697 29.614 185 33 9798 30 190 0 9899 30 190 1 99（1）血管半径为 =29
16、.6799m。R（2）中轴线在 XY 平面投影曲线的拟合方程(R-square:0.9939)256-104.7X8.59- 0.17693X.8-X04.26.09-X1.748632 45-87819- Y（3）中轴线在 YZ 平面投影曲线的拟合方程(R-square:1) 8625 3456-57-19-103.Y.70. 81.Y-3.6290.17Y2.16 Z（4）中轴线在 ZX 平面投影曲线的拟合方程(R-square:1)92.8.460Z0.52- 0.5368Z.29-Z107.561.23-Z.61-19 45-78-3 X（5）中轴线方程： ，92.8.460Z0.52
17、- 0.5368Z.29-Z107.5610.23-Z1.65-19,03Y.7. 1.Y-3Y-32,-04.X8. 0.7693X.8-X104.2610.9-X1.74158162 456-78-625 5-89- 45-87-4-XZY0 100 200300 400 500250300350400450020406080100x中中中中中yz0 10 2030 40 50250303504045002040608010x中中中中中中中yz图 1 中轴线曲线 图 2 中轴线拟合曲线50 100 150 200 250 300 350 400 4502402602803003203403
18、60380400420440xy中中中中XY中中中中中中中50 100 150 200 250 300 350 400 450240260280300320340360380400420440xy中中中中XY中中中中中中中中中图 3 中轴线在 XY 平面的投影 图 4 中轴在 XY 平面的投影拟合曲线240 260 280 300 320 340 360 380 400 420 4400102030405060708090100yz中中中中YZ中中中中中中中240 260 280 300 320 340 360 380 400 420 4400102030405060708090100yz中中
19、中中YZ中中中中中中中中中图 5 中轴线在 YZ 平面的投影 图 6 中轴线在 YZ 平面的投影拟合曲线0 10 20 30 40 50 60 70 80 90 10050100150200250300350400450zx中中中中ZX中中中中中中中0 10 20 30 40 50 60 70 80 90 105010150202503035040450zx中中中中ZX中中中中中中中中中图 7 中轴线在 ZX 平面的投影 图 8 中轴线在 ZX 平面的投影拟合曲线图 9 曲线模拟的血管6 模型检验与分析用多项式拟合得到重建后血管的三维空间曲线，再用平面 去截)75,021(,iZ空间曲线得到新
20、的 4 层切片平面分别与题目中对应切片进行对比，在这里我们称 的TS面积与 的面积比值为切片重合度。TS第 01、25、50、75 张图片分别与对应的新截取图片的对比图. . .图 1 第 01 张切片 图 11 第 25 张切片 图 12 第 50 张切片 图 13 第 75 张切片图. . .图 14 截后第 01 张切片图 15 截后第 25 张切片图 16 截后第 50 张切片图 17 截后第 75 张切片根据上面曲线拟合得到的血管重新切片后，我们随机取出四张重切的图片，与原始的切片比较，得到两切片的重合度。计算得：第 01 张切片的重合度为 99.20% ，第 25 张切片的重合度为
21、 98.65% ，第 50 张切片的重合度为 97.30% ，第 75 张切片的重合度为 97.61% ，计算平均重合度得到 98.19%，说明本模型较为可靠。7 模型的优点与缺点优点在于本模型的算法采用遍历搜索的方法得到最大内切圆的半径的精度较高，但是由于遍历搜索的运算次数异常庞大，导致模型的缺点是实践难度较大。从模型检验的数据可以看出计算的精度，但是整个运算在 Pentium(R) Dual-Core 处理器上耗费的时间已经超过 120 分钟。参考文献1刘会灯，朱 飞.MATLAB 编程基础与典型应用M 2008 年 3 月附录% program 1 （找出切片的边界线的图）pic=zer
22、os(512,512,100);pic2=zeros(512,512,100);for i=0:99s=sprintf(d:Apics%d.bmp,i);pic(:,:,i+1)=imread(s);pic2(:,:,i+1)=edge(pic(:,:,i+1);imshow(pic2(:,:,i+1);Endm,n,i+1=find(pic2(:,:,i+1);% program 2 （求半径，找中轴线与切片交点）ss=100;lens=1000*ones(512,512,ss);centres=zeros(ss,2);r=zeros(ss,1);for k=1:ssfor i=1:512f
23、or j=1:512if pic(i,j,k)=0for m=1:512for n=1:512if pic2(m,n,k)=1t1=sqrt(i-m)*(i-m)+(j-n)*(j-n);if lens(i,j,k)t1lens(i,j,k)=t1;endendendendendendendfor i=1:512for j=1:512if lens(i,j,k)=1000lens(i,j,k)=-1;endendendr(k,1)=max(max(lens(:,:,k);for j=1:512for i=1:512if r(k,1)-lens(i,j,k)s1(i,j)s2(i,j)=s1(i,j);endendif s2(i,j)-29.67990.00000001plot(i,j);hold onaxis equalpc(i,j,z)=1;endendendend% program 6 （计算原切片与新切片的重合度 p）common=0;inner=0;p=0;for i=1:512for j=1:512if a1(i,j)=0 endif a1(i,j)=0inner=inner+1;endendendp=common/inner;% All programs is over.%
