血管的三维重建.doc

1血管的三维重建 [摘要] 本文分别利用轮廓提取[1]、轮廓跟踪[1]两种方法，把100 张切片图轮廓上的点的坐标提取出来，然后利用二维搜索求出所有切片图的最大内切圆半径以及圆心的坐标，此半径即是所求血管的半径，此圆心即是切片图与中轴线的交点根据这100 个交点的坐标利用Matlab 软件[2]分别画出了血管、中轴线的三维立体图以及中轴线在各个平面上的投影本文最后对两种模型的结果进行了误差分析，得出了比较精确的结果:血管半径: r ≈ 29[关键词]轮廓提取 轮廓跟踪 二维搜索 三维重建一 问题的重述断面可用于了解生物组织、器官等的形态例如，将样本染色体切成组织切片，在显微镜下观察该横断面的组织形态结构如果用切片机连续不断的将样本切成数十成百的平行切片，可以逐片观察根据拍照采样得到的平行切片数字图像，可运用计算机重建其组织器官的准确三维形态假设某些血管可视为一类特殊的管道，该管道的表面是由球心沿着某一曲线(称为中轴线)的球滚动包络而成现有某管道的相继100 张平行切片图像，记录了管道与切片的交图像格式为BMP,宽高均为512 个象素(pixel)，切片间距和图像象素间尺寸均为1试计算管道的中轴线与半径，给出具体的算法，并绘制中轴线在XY、YZ、ZX 平面的投影。

二 问题的分析本问题要求重建血管模型解决此问题的关键在于如何将100 张图片转化为便于处理的数据结构我们考虑先将对应图像文件中的特征数据读出，筛选出数量较少的关键点坐标针对本问题BMP 文件格式的特点，可以用Matlab 软件[2]编程扫描出图像中对应象素点的索引值鉴于所得数据十分庞大，应该只保留其中代表血管截面的数据在图像给出的所有象素中，与本问题的解决相关的仅仅是血管边缘的象素我们分别采用轮廓提取算法和轮廓跟踪算法[1]，通过编程得出边缘点的坐标观察最终所得数据，我们可以直接构建出管道三维立体图原问题指出，所要重建的血管由半径固定的球滚动包络形成，血管中轴线与每张切片有且只有一个交点,所以对各切片图血管截面最大内切圆求解得到的最大内切圆半径即为血管半径，圆心即为血管截面与中轴线交点求得圆心坐标后，中轴线也可得出三 问题的假设1. 所给图片边缘连续，即各象素都相邻，无孤立象素块；2．所给各切片图血管截面单连通；3．血管壁的厚度不计；4．图像中的每个象素均可视为质点四 符号约定i AA i = 0,1,2,...,99 第i张切片图的象素索引矩阵；i A i = 0,1,2,...,99 第i张切片图中包含血管截面的最小区域；i .bmp i = 0,1,2,...,99 第i张切片图的文件名；i PX i = 0,1,2,...,99 第i张切片图中血管截面边缘象素横坐标向量；i PY i = 0,1,2,...,99 第i张切片图中血管截面边缘象素纵坐标向量；p 切片图中血管截面边缘象素的数目；m n a , 象素点(m, n) 的索引值；；2pp pp = 1,2,3,..., p 血管截面边缘象素的标号；dist( pp) 任一内部点到第pp 象素的距离；i f i = 0,1,2,...,99 第i张切片图血管截面最大内切圆圆心坐标；i r i = 0,1,2,...,99 第i张切片图血管截面的最大内切圆半径；r 血管半径。

五 模型的建立及求解模型一第一步：读取 BMP 文件1．获取切片图的象素索引矩阵所给切片图的 BMP 文件为512× 512象素黑白色的二值图像，运用Matlab5.3 软件读取BMP文件，返回象素索引矩阵，记为i AA 象素索引矩阵为512× 512的二值矩阵，索引值1代表黑色，2代表白色2．简化象素索引矩阵.象素索引矩阵对切片中血管截面象素的索引值为1，对其它象素的索引值为2 为了减小运算量，我们只需要分离出血管截面的象素① 搜索出象素索引矩阵中索引值为1 的象素的最小、最大x 坐标和y 坐标， 分别记为min X ,max X;minY ,maxY ②将矩形区域min X → max X ,minY → maxY 定义为i A 第二步：轮廓的提取1．找出血管截面部分，对该部分进行轮廓提取按行列顺序搜索区域i A 中任意点周围相邻的8 个点的索引值，如果均为1，则将该点坐标记录下来例如：在区域i A 中，点(m, n)和周围的8 个点索引值均为1，表明该点为血管截面上的非边缘点，将m n a , 记为3，由所给问题可以得出：1 , = m n a (m, n)为血管截面的边缘点；2 , = m n a (m, n)为血管截面的外部点；3 , = m n a (m, n)为血管截面的内部点。

2．作出血管截面的轮廓在i A 中搜索所有1 , = m n a 的点的坐标(m, n)，并绘制出来，得到图1将图1 中矩形内区域放大得到图2由图2 可见折线右边的点均为血管截面内部点，而左边的点为边缘点将内部点，即3 , = m n a 对应的点剔除后，剩下的点便是边缘点然后依次找出m n a , 为1 的所有点的坐标，记为(PX ( pp),PY ( pp)) i i 通过采用轮廓提取算法，我们得出第i张切片层血管截面边缘坐标为(PX ( pp),PY ( pp)) i i ，边缘点数量为p ，( pp = 1,2,3,...p)仍以第77 张切片图为例，利用Matlab 软件编程绘出边缘图形（图3），边缘点总数为1226第三步:血管的三维重建用上述方法处理得到100 张切片的边缘轮廓图和边缘点的X,Y 坐标，由题意第i张切片层的Z 坐标为i通过Matlab 软件编程（见程序）依次读取文件i.bmp (i = 0,1,2...99)，提取出所有边缘点的三维坐标由于所得数据量充足，可以直接绘出三维图，见图4第四步：血管半径及中轴线的求解由于图 4 的三维模型与第77 张图的边缘图吻合良好，这里仍以之为例：31．求截面内部点到边缘的最短距离① 在Ai区域内沿行列方向逐点判定搜索内点(m,n)与边缘上第pp 点的距离dist( pp)如果1 , = m n a ，表示当前点在边缘上，不需计算；如果2 , = m n a ，表示当前点在切片外，不需计算；如果3 , = m n a ，表示当前点位于边缘内部,参与计算。

图1 图 2图 3 图4200 250 300 350 400 450 5000100200300400dis( pp) (m PX ( pp))2 (n PY ( pp))2 i i = − + − ② 得出dis( pp)中的最小值，即d = min(dis( pp))，并记录下对应点的坐标③ 逐次比较所有内部点返回的d 值,可以确定其中的最大值是当前切片层血管截面的最大内切圆半径2．求解中轴线因为血管截面的最大内切圆半径即为血管半径，即i r 为所有内部点返回d 值中的最大值① 先令i r 为0；② 对任一内部点求它与边缘上各点的距离dist( pp)；③ 求出d = min(dist)比较当前内部点返回的d 值与i r 的大小若r d i < ,则令r d i = ④ 循环至②，直到所有内部点搜索完毕在搜寻完成后， i r 表示当前切片层的最大内切圆半径， i f 表示该内切圆的圆心坐标在将所有切片层数据处理完成后，我们可以根据i f 构建出中轴线已知管道半径一定，所以各切片最大内切圆半径即为4管道半径所有数据见附录）因为本题目所给数据在边缘处离散，对各切片层求出的最大内切圆半径间存在较小误差。

为了减小误差，可对ri求平均值，求得r = 28.8879，近似得到r = 29模型二观察图2 可以发现，模型一所得边缘点由接近平行的两组点组成，边缘点数较多，直接影响最大内切圆半径的求解速度和求解精度为提高模型的求解效率和质量，我们采取轮廓跟踪算法来提取边界点的坐标具体算法是：只要点(m,n)上下左右四个相邻点的值均为1，就认为点(m,n)为内点，对内点做异于边缘点标志的记号图5通过图1 与图4、图2 与图5 对应比较，发现模型二处理的图像边缘更精确用模型二对该图进行计算，得出边缘点总数为865，只有模型一的2 / 3，大大降低了数据量，提高了100 张图片的处理速度下表分别列出了模型一、二对所有图像处理的情况模型 边缘点的总个数 程序运行时间（s）模型一 轮廓搜索 85804 263.1500模型二 轮廓跟踪 60704 133.5700由上表可以看出，采用模型二搜索边缘，效率提高了1/2根据模型一、二得出的结果分别绘出血管的立体图，如图6，图7 所示经模型二计算得到血管半径r = 29.1884，近似得到r = 29图 8、图9、图10 分别给出了血管中轴线在xy 平面， zx 平面， yz 平面上的投影。

六. 模型的评价本问题的特点是所给数据在边缘处离散，各切片层最大内切圆的圆心位置相互独立，且从初步模拟结果可以看出，部分相邻切片层间圆心坐标变化较大为了求得相对精确的结论，我们需要对每张图片进行计算，求出所有切片层最大内切圆因为计算量十分庞大，所以，合理地简化计算数据十分重要模型一提出的算法仅对所有数据做出初步选择，模型二在此基础上进一步简化了模型一得出的边缘点，取得了较好的效果七．模型结果的分析由理论分析，在每个模型中求到的半径i r 应该完全相同但是在结果中发现有差异，造成误差的原因是假设每张切片的边缘是连续的，而读取BMP 文件所得到的数据是离散的，对i r 求算术平均值后得出了比较精确的解在结果分析时发现，两个模型求得的半径也有相对较大的差异这是由两个模型的轮廓点的坐标的提取方法的差异（轮廓提取算法与轮廓跟踪算法的差异）造成的由于在轮廓跟踪算法中确定的轮5廓能更好地反映血管切片的边缘，故模型二求得的半径精确度更高，最后得出血管半径r = 29.1884，近似得到r = 29参考文献[1] 吕凤军，数字图像处理编程入门. 清华大学出版社，1999 年[2] 蒙以正，柳成茂. MATLAB 5.X 应用与技术. 科学出版社，1999 年图6 图7 图8 xy 平面上的投影0 20 40 60 80 100 120 140 160 180-200-150-100-50050100150200x图9 zx 平面上的投影图 10 yz 平面上的投影0 20 40 60 80 100020406080100120140160180200 150 100 50 0 -50 -100 -150 -2000102030405060708090100x6附录中轴线的三维立体图xy 平面上的投影0501001502000100200300400500050100x0 20 40 60 80 100 120 140 160 180-200-150-100-50050100150200xxz 平面上的投影 yz 平面上的投影0 20 40 60 80 100 120 140。