三维血管重建郑志文 陈波 朱玲

三维血管重建

摘要

本文通过对血管样本切片图像的分析，具体采用枚举法，并利用Matlab模拟实现血管的构建，得到了该血管的三维形态。

对于计算管道的半径，通过读取每张图片，将其变换为0-1矩阵。运用Matlab 找到每张切片所对应的矩阵部分，找出该部分的圆心。将得到的100张切片的圆心取平均，即为该管道的半径为29.675个象素。

对于计算管道的中轴线，首先，将找到的每张切片的圆心运用Matlab相关函数命令，得到了中轴线的空间曲线（见图一）。然后，通过Matlab编程将每张图片叠加，得到了三维血管重建图（见图二）。最后，运用Matlab相关命令，得到了中轴线在XX、YZ、ZX平面的投影图（见附录中图三、图四、图五）。

对于本模型的误差，划定一个较小的范围，即来自求切片圆心的确定，在此刻画一个反应误差大小的量，对本模型进行检验。通过编程统计出100张图片中不符合总体圆心规律点的个数所占的比例大小为5%以下，即在误差许可范围内，本模型合理。

本模型的优点是考虑到了提取图像的边缘点，减少了计算量；程序所涉及的语句简单且计算100张图片时间不到2分钟;算法的时间复杂度为3n,空间复杂度为(2)

O,相对较小。

关键词：轨迹内点边界点

一、问题重述

断面可用于了解生物组织、器官等的形态。例如，将样本染色后切成厚约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、YZ、ZX平面的投影图。

二、问题分析

2.1名词解释

像素（Pixel）：像素是指基本原色素及其灰度的基本编码。是由Picture(图像) 和 Element（元素）这两个单词的字母所组成的，是用来计算数码影像的一种单位，如同摄影的相片一样，数码影像也具有连续性的浓淡阶调，我们若把影像放大数倍，会发现这些连续色调其实是由许多色彩相近的小方点所组成，这些小方点就是构成影像的最小单位“像素”（Pixel）。这种最小的图形的单元能在屏幕上显示通常是单个的染色点。越高位的像素，其拥有的色板也就越丰富，越能表达颜色的真实感[1]。

边界点：灰度图像中黑白部分交替地方的所有黑色点的集合。

内点：灰度图像中除去黑色部分中边界点的所有点的集合。

2.2问题分析

根据题意中假设血管可以视为一类特殊的管道，该管道的表面是由球心沿着某以曲线（称为中轴线）的球滚动包络而成和假设1这两条线索，可以分析得到以下几点：

1.管道的中轴线是由滚动的球的球心轨迹连线而成，而且每张切片必过

中轴线上的点，即切片都是切过球心的点。

2.每张切片只存在一个最大的内切圆，且过球心的点（即与中轴线的交

点）为该最大内切圆的圆心。

通过细致的分析得出，该问题的关键之处在于求解每张图片最大内切圆的圆心。由于给出的数据是一些BMP格式的图像，为了从图像中获取相应的信息，利用将图像处理成数据的专业数学软件Matlab，把每张图片读取成一个0—1二维

矩阵。在矩阵中，1代表的灰度图像中的白色部分，0代表灰度图像中的黑色部分。通过对图像的处理，问题简化为对我们熟知的数字矩阵进行操纵。回归到本问题的重心是在于求解每张切片的圆心以及半径，对每张切片对应的0—1二维矩阵来说，即要在0的部分中求圆心以及不超过0的部分中求一个最大圆的半径。

针对如何计算管道的半径，首先，根据一个圆中只有圆心到各个圆边点的距离最大为出发点，找出每张图片的内点和边界点；然后，通过Matlab 编程找到每个内点到所有边界点的最短距离；最后，取最短距离集合中的最大距离，即为该张图片的圆心。

对于管道的中轴线，首先，将每张图片的圆心x,y 坐标分别存储在一个数组里；然后，调用Matlab 中的一个画图函数即可得到该血管的中轴线所对应的示意图；最后，继续调用函数便可简单地绘制出该中轴线在XX 、YZ 、ZX 平面所对应的投影图。

对于血管的重构，采用将所有切片图一起叠加，即得到一个相似的血管重建图。

三、 模型假设

1、假设管道中轴线与每张切片有且只有一个交点，球半径固定，切片间距以及图象象素的尺寸均为1；

2、假设每张切片都是等间距、等像素、平行的切片；

3、假设每张切片的厚度可以忽略不计；

4、假设血管内外壁的厚度可以忽略不计，即中轴线到血管内壁的距离等于球半径；

5、假设每张切片厚薄均匀。

四、 模型的建立与求解

5.1模型的准备

首先， 把题中给的数据图片导入到Matlab 中，运用Matlab 的读图函数，运用for 循环语句连续的读出每张图片。由于图片为黑白，即可视为灰度图像，即读出的图像是以0-1二维矩阵存储的。在矩阵中，1代表的灰度图像中的白色部分，0代表灰度图像中的黑色部分。然后，分析矩阵中的数据并试着找到该问题的关键点，即找每一张图像的最大内切圆的圆心，接着将这些圆心点存储，运用matlab 相应的函数生成一条曲线，即为所求的管道的中轴线。最后，通过找每个内点到各个边界点的最短距离，然后取这些最短距离的最大值即为半径为动机，设计相应的模型算法对其编程求解。 5.2模型的建立

以任意一张切片的所有内点到边际点的最短距离为目标函数：

min ij S =(,)i i j j x y ??内点，(x ,y )边界点

以内点到边界点的最短距离的最大值为目标函数：

max {}ij S

5.3模型的算法

通过详细分析，把所要求解的图像细化为一个个像素点，把除去背景点的部

分分为边界点和内点。通过模型的准备，已经确定了每个像素点在矩阵中对应的

像素值。具体算法步骤为：

Step1：运用matlab中的imread()函数将100张图片的信息依次读入到各

自相应的矩阵;

Step2：再运用matlab中的edge()函数，将每张图片对应的边界点提取并

存入矩阵；

Step3:运用find()函数找到边界点和内点所对应的索引值即坐标，其中内

点是通过对原图片矩阵用符号~取反并减去边界点得到的；

Step4:搜索每个内点到所有边界点所对应的最短距离，并取所有最短距离集

合的最大距离，并记录其坐标。

5.4模型的求解

由于Matlab专业图像处理能力以及对矩阵的灵活操作，因此使用Matlab进

行编程求解，批量读取切片，下表给出了部分切片的半径以及坐标。

表一、部分切片的半径及坐标

Pi P0P1P2P3P4P5 R 29.0688837 29.06888371 29.017236 29.0688837 29.068884 29.0688837 (X,Y) (-160,1) (-160,1) (-160,2) (-160,2) (-160,2) (-160,2) Pi P23P24P25P26P27P28 R 30 30 29.832868 29.6141858 29.614186 29.6141858 (X,Y) (-160,14) (-160,15) (-160,16) (-160,17) (-159,26) (-159,27) Pi P51P52P53P54P55P56 R 29.6984848 29.69848481 29.698485 29.6816442 29.681644 29.6816442 (X,Y) (-119,112) (-118,113) (-116,115) (-111,120) (-98,131) (-91,136) Pi P93P94P95P96P97P98 R 29.6984848 29.69848481 29.698485 29.6816442 29.681644 29.6816442 (X,Y) (185,33) (183,40) (187,24) (187,24) (189,12) (190,0)

通过计算编程算出了一百张图片的半径平均值等于29.675个像素。又通过

程序确定出来的每张图片的圆心坐标，用Matlab画出了中轴线在空间的图像如

表二：

图一、中轴线的空间曲线

由中轴线画出了三维血管重建图，如图二：

图二、三维血管重建图

同样用Matlab画出了中轴线在XY、YZ、ZX的投影，分别见附录表三、表四、表五。

5.4误差分析和模型检验

现实管道的切面边缘应是连续的，而BMP图像是由离散的像素点组成，所以图像的最大内切圆可能有多个，而本模型只是取其中的一个，这就造成了一定的误差。

针对可能出现的误差，本模型定义了一个反应误差大小的量，对模型作了一定的检验。

分析中轴线在ZX、YZ上的投影，可以发现其在ZX的投影总体近似于一条单调递增的曲线，那么后一个图像的最大内切圆圆心的横坐标值应该大于前一个图像；同样中轴线在YZ上的投影具有一定的单调性，通过计算可以得到前69张图片的内切圆圆心纵坐标值递增，后面的图片纵坐标值递减。由于误差的存在，使得一些图片的内切圆坐标值不满足上述特征。通过对相邻两图像的最大内切圆圆心坐标的比较，可以定义两个个量分别纪录不满足上述特征的横、纵坐标组的个数，用以描述误差的大小。记这两个量为anx和any（0

五、模型的评价与改进

7.1模型的评价

(1)模型的优点

1、灵活运用了Matlab中工具箱的函数，快速简单的找到了边界点和内点，提高了计算的速度，使得本模型计算每一张切片的半径以及坐标不超过2秒。

2、本模型的算法简单，虽说是用枚举法，但巧妙地运用的相关函数对其缩小了搜索的范围。相比其他复杂高效的算法，有一定的优势。

(2)模型的缺点

1、本模型的缺点是内点到边界点的最短距离中的最大值可能同时有多个内点，

而本模型的算法中并没有对这些内点的坐标进行取平均值。

2、本模型算法中使用的内点过多，没有合理的剔除一些不可能是圆心点的点，以减少计算速度。

六、参考文献

[1] https://www.360docs.net/doc/2416865128.html,/view/575.htm

[2] 薛定宇，陈阳泉.高等应用数学问题的MATLAB求解（第二版）[M],北京:清

华大学出版社. 2008.10

七、附录

图三、中轴线在XZ的投影图四、中轴线在YZ投影

50100150200250300350400450

图五、中轴线在XY投影