作者：王斯刚 冯有前 赵学军 王锦江
【关键词】 血管切片
关键词: 血管切片;中轴线;搜索算法;内切圆;拟合
Abstract:Firstly，taking use of searching algorithm to calcu-late coordinate of point of intersection for kinds of blood vessel’s axis demarcation line and radius of incircle.And then，averaging the radius of incircle which was calculated by a piece of dicing.It caught pipline radius of blood vessel by simulating all points of intersection.
0 引言
在医学、天文观测、工业非破坏性试验等一些工业问题中，需要确定某个空间物体的形状，但由于技术上的限制，无法将物体分离出来并直观的显现.解决此类问题的一个可行的方法是用等间隔的平行平面去截取这个物体，得到一组平行截面，通过采样得到平行截面的数字图像，进行数据处理后可运用计算机重建物体的三维形态.在2001年全国大学生数学建模竞赛题中，A题给出了某血管管道的相继N张平行切片图像，图像格式、尺寸及坐标均已给定，这些图片记录了血管管道与切片的交点，并且假定该血管管道是由一个小球沿某路径移动而形成的.取z轴垂直于切片，规定第一张切片平面为z=0，第N张切片平面为z=N-1.要求解决的主要问题是:求解血管管道的中轴线与半径，同时给出具体算法;以及绘制中轴线在XY，YZ，ZX平面的投影图.现在我们试图利用这组平行截面所提供的信息来重建被截血管管道.
1 一些假设
题中已给出的假设:（1）管道的中轴线与每张切片有且只有一个交点;（2）形成管道的球半径固定;（3）切片间距以及图像像素的尺寸均为1.
求解所需的假设:（1）假设球是沿光滑曲线运动的;（2）将像素抽象成为一个点，由于像素宽度与切片间距相等，所以由原数据得到的是一个均匀的三维点阵.（3）像素点的坐标位置在像素点的中心.
2 问题分析
此类问题的关键是:如何利用所已知的一组等间距的平行截面，根据其平行截面的形状确定球心轨迹（此轨由于采用本文所用的方法，参与“2001年全国大学生数学建模竞赛”的建模小组获得了“国家一等奖”. 迹是一个工程上可接受其精度的近似值）.由于问题中假设管道中轴线与每张切片有且只有一个交点，则切片图像在XZ，YZ平面的投影是以z为自变量的函数.可以通过求出单个截面与中轴线的交点，再将所有截面与中轴线的交点进行曲线拟合得到所要求的结果.具体分析思路如下.
2.1 对单个截面的分析 文中所给的单个截面，我们可以认为是形成管道的球体沿某一空间曲线穿过一截平面，球体被截平面连续截取所得的圆在该截平面上形成的包络即为截面.由截面形成过程可知，截平面所截得的包含球心的圆就是截面的最大内切圆，其圆心坐标为截面与中轴线的交点坐标.由单个截面所给的信息，我们便可确定每个切片与管道中轴线的交点，并通过求解其最大内切圆的平均半径得出球体半径.形成管道的平移球体穿过截面时，所形成的相交图形是一系列半径可变的圆，可设其半径函数为r=r（t），t1
2.2 对所有截面的分析 单个截面所提供的是局部坐标所包含的内容，按上述同样思路求得所有切片与中轴线的交点.所有截面的信息则包含了连接各切片的中轴线的内容，将所有的交点拟合后便可得到中轴线的近似轨迹.
3 问题求解
基于上述分析，对于一个给定的截面形状，就可以得到这个截面与管道中轴线的唯一交点，若求出每个截面与中轴线的交点，则通过空间曲线的拟合就可以得到原管道的中轴线.
3.2 搜索最大内切圆圆心Ci 由切片及截面形成的过程可知，在截面上必定存在一个最大内切圆，中轴线与截面的交点即为该内切圆的圆心.
算法描述:（1）将边界点的坐标按顺时针（逆时针也可）的方向依次写入一个坐标序列（Xi ，Yi ）（i=0，1，2……N-1，N为边界点总数）中;（2）对任意一个截面内的点（即黑色的点），记录其到边界的最短距离dmin p ;（3）扫描所有截面内的点，对于每一个点按上述方法求出到边界的最短距离，找出这些最短距离中最大的一个所对应的点，则此点即为最大的内切圆圆心，亦即管道中轴线与此截面的交点Ci .
为了进一步精确求得圆心位置，我们进行了以下修正:将所取的参考点在其周围的8个方向上都减小1/2步长，再以上述同样的方法作圆，这样所得的半径最大的圆就与最大内切圆更加接近，所确定的该参考点即为圆心.
3.3 拟合中轴线 通过MATLAB软件中相应功能，应用最小二乘拟合将所有切片与中轴线交点在空间进行拟合.具体拟合方法如下:对于上述方法求出的C（xi ，yi ，zi ），i=0，1，2，…N-1，将其投影到各个坐标平面上，即得到XY平面上的（xi ，yi ），YZ平面上的（yi ，zi ），XZ平面上的（xi ，zi ），应用最小二乘法将二维点用一个多项式拟合，应用MATLAB的PULYFIT函数选取适当的拟合次数进行求解，得到中轴线在三个坐标面上的投影曲线的方程［2］ .
使用MATLAB软件画出两两联立的空间曲线，选择最优的一条，得到最终的空间曲线为:XZ平面与YZ平面的投影拟合方程联立.y=0.4371-0.7134z+0.1308z2 -0.0082z3 +0.0003z4x=-160.4458-2.1263z+0.12019z2 -0.2649z3 +0.0301z4 -0.002z5 +0.0001z6具体图像见图1.
图1 血管管道中轴线 略
4 算法改进及推广
4.1 算法改进 在上述问题中寻找最大内切圆时，程序采用的是对边界内的所有点进行逐个搜索求解，这样处理程序效率不高，为了减少在寻找最大内切圆时所搜索的点数，设各个边界点均有按序的标号，我们将具体的算法改进如下:首先，找出内切圆圆心.在所求切片截面中选取一内部点作为基准点，由搜索法求得边界上所有点到此点的最短距离r1 .若有 r2 -r1 &<ε成立（r2 为其他边界点到该基准点的距离，ε为任意小的正数），再判断满足上述条件的点中，是否有标号相差较大（大于5）的两点.若存在至少两点满足上述条件，则可确定该基准点为内切圆圆心.其次，选定r1 在大于一个规定的值a的范围内搜索.由前边分析可知，内切圆半径随x，y变化，会在中间某处出现最大值，所以，若r
1 a的方向移动，可确定作为参考的内切圆圆心.最后，在此参考圆心的周围取一矩形区域，在区域内搜索半径更大的内切圆圆心，并逐步迭代［3］ .这样就将所需搜索的范围大大压缩.经改进后的算法在处理更复杂的数据群时，更能体现出优越性.
4.2 算法的推广 考虑到实际中生物组织、器官等的形态是一类特殊的管道，该管道的表面与我们上述讨论有所不同，是由半径在连续变化的球体其球心沿着某一曲线滚动包络而成.由于管道半径是连续变化的，所以在较小的变化范围内，可近似视为其球体半径是固定的，按问题所提供的算法，再进行相应的处理.
参考文献:
［1］黄友谦，李岳生.数值逼近［M］.第2版.北京:高等教育出版社，1999:97-131.
［2］王沫然.MATLAB5.X与科学计算［M］.北京:清华大学出版社，2000:35-57.
［3］寿纪麟.数学建模 方法与范例［M］.西安:西安交通大学出版社，1999:73-127.1