血管三维重建模型

全国一等奖血管三维重建模型指导教师：覃思义学生：朱天飞丁理杰孙红芳血管三维重建模型摘要本文对重建血管三维图象的问题建立了数学模型。

通过对问题的分析，我们认为首先求出半径是解决该问题的关键，在此基础上应该用尽量少的搜索次数找出每个切片的轴心。

我们用“选优搜索”法首先求出了血管半径在28～30个象素点之间。

每幅切片图上的轴心一定在“有效象素点”（即血管切片边缘内的点）上，我们对图象的象素点进行了分析，为了减少搜索范围，得到了对“有效象素点”的筛选方法，得到了“筛选点集”。

由此确定了一个很小的搜索区域。

利用半径的平均值，在此范围内通过较少的搜索次数求得各个切片的与中轴线的交点。

我们对所得结果进行了分析。

对明显的偏差点进行了纠偏，对纠偏后结果的进行了有效性检验，所有轴心坐标的平均有效率为99.8％。

为了得到较精确的中轴线在X-Y、Y-Z、X-Z面的投影，对纠偏后的所有坐标进行分段线性插值。

由此得到了较精确的投影和三维散点图。

最后我们给出了该算法的评价，提出了更有效的搜索算法。

一．问题的重述与分析在对生物组织、器官等的形态的研究中，需要知道组织、器官的三维形态。

其中的一种途径是通过样本切片来重建其三维图象。

其大体的步骤为：用切片机对研究的样本进行平行切片；将切片拍照并采样得到平行切片的数字图象；运用计算机建立三维重建的算法，根据二维数字图象的信息重现原样本的三维形态。

我们现在研究的对象是一种血管样本。

这种血管的特点是：血管管道的表面是由球心沿着某一曲线的球滚动包络而成。

现在已经得到了该血管的100张平行切片的图象，每张图象的象素点为512×512个，格式为bmp。

为了简化，已有下面三个重要的假设：1．管道中轴线与每张切片有且只有一个交点。

2．滚动包络成管道的球的半径是固定。

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

需要建立计算机算法，通过图片所包含的信息，对样本进行三维重建，再现血管的三维形态。

最后需要得到血管管道的半径和中轴线三维空间的图象以及它分别在XY、YZ 、XZ平面投影的图象。

如何建立有效的算法，是我们所要解决的问题。

问题的一个重要要求是求得中轴线的三维和它投影的二维图象。

又根据上面的假设1，可知100张图片中包含且仅包含有100个中轴线上的点，它们是图象包含的最重要的信息，我们称之为“中轴点”。

所以，我们将建立算法的目的是“提取每幅图象中的中轴点”。

根据假设1、2，通过分析，我们得到一个重要的结论：结论一：每幅图上有且只有一个最大内切园，该圆的圆心即是这幅图上的中轴点。

我们将每幅图上的这个圆称为“关键圆”。

根据这个结论，问题转化为在每幅图上查找其切面围线内的最大内切圆。

在上述的分析基础上，我们另外发现，图象中的对查找中轴点有着决定性作用的象素点集中在关键圆的周围。

可先将这些最有效点提取出来，再在点集所在的范围内查找关键圆点。

求出每个面上的中轴点，对其进行处理，可以得到所要求的中轴线图象。

二．定义与符号说明中轴点：切片与血管中轴线相交的点。

切面：切片图象上血管边缘所包围的面。

切面围线：图上血管切片的边缘所围成的封闭曲线。

关键圆：切面围线内的最大内切圆。

原图：所提供的100幅切片图。

复合图：将100幅原图重叠，实质是管道在xy平面上的投影图。

R：血管管道的半径。

图阵：原图经过转化的矩阵，大小为512×512，包含着原图象信息。

三．基本假设⑴管道中轴线与每张切片有且只有一个交点。

⑵滚动包络成管道的球的半径是固定。

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

四．模型的建立与求解1）准备工作(1)获取初始数据：对所给出的100幅bmp格式的图，通过数学软件matlab的imread命令，转化成matlab能处理的图阵。

该图阵是一个元素值只有0 和1 的矩阵，每幅图转换后的大小为512×512。

矩阵中为0的点对应图上的灰色度高的点（即切面上的点），元素为1的点对应灰色度低的点（即切面外的点）。

（2）我们将100幅原图重叠，得到一个复合图，如下图该图是管道在X-Y平面上的投影，同时也将它转换成图阵。

我们认为首先应确定出血管管道半径，然后在此基础上进行中轴点的搜索。

（3）对关键圆进行搜索时，必定使用到原图的象素点的信息。

每幅原图有效点（即黑色点）点的数量非常大，如果全使用到的话，其工作量将很大。

所以，我们考虑先对图阵上的有效点进行筛选，再作进一步的处理。

通过对图象的分析，结合结论一，我们发现切面上分布在最大内切圆附近的象素点，（插图进行说明）尤其是其中分布在切面围线上的点，对这种最大内切圆的查找有着重要的作用。

所以，我们建立“筛选算法”对有效点进行筛选。

通过将每幅图的切面围线与复合图的切面围线进行比较，其重合的点即为要找的切点。

算法如下：1．将复合图的边缘图求出（可以通过matlab的edge命令得到）。

2．将第n幅图的边缘图求出（n＝1，2…100）。

3．比较原图与复合图的边缘图，取其重合的点作为要找的筛选点。

4．重复2到4，直到找完所有的图的筛选点。

于是对于第i幅图，就有筛选点集S i与其对应。

2）确定血管半径：针对求解血管半径的问题，我们专门给出了一个算法。

从假设一、二和结论一可以得出，理论上讲，从各个不同的切片图求出的血管半径是一样的但是该图是用离散代替连续，所以所搜索出的半径来说是一个范围。

从复合图上容易得知其直径在60个象素点左右。

由此，我们得出以下算法：1．任取一幅切片图，取其“筛选点集”。

2．在刚能包容“筛选点集”的矩形区域内，采用逐点取点的方法选取测试圆的圆心。

即从矩形的左上端开始，按照“从左到右”，再“从上到下”的方式选取圆心。

3．在25～35个象素点范围内，以从小到大的次序，依次取一个值作为测试圆的半径R i。

4．记录“筛选点集”中落在测试圆上的点数，如为2，继续往下执行；否则，跳至3。

5．判断落在圆上的两点与圆心是否共线，若是，将此时的其存入数组B中；否则，跳至3。

6．若矩形区域上的点未取完，跳至2继续执行；否则继续到7。

7．取数组B中最大的R i作为血管的半径。

算法的程序见附件三。

根据程序求出的半径为的范围是28～30个象素点，我们取其平均，R取30个象素点。

3）搜索中轴点的算法由于已得到点数较少的“筛选点集”，而且最大内切圆与切面围线的一对切点必在其中，所以可以根据最大内切圆与这对切点的有关几何关系进行关键圆圆心的搜索。

我们利用的几何关系有两个：1． 对于上面提到的最大内切圆，它与切面围线有且只有两个切点，切点到圆心的距离必为半径R 。

2． 圆的圆心为O(x 0,y 0)点。

过O 点与圆交于两点M 1（x 1,y 1）,M 2(x 2, y 2)则有(y 1－y 2)×(x 1-x 0)＝(y 1-y 0)×(x 1-x 2)具体算法：1． 确定第j （j 的初值取1）幅图的可能中轴点的搜索范围。

找出筛选点集S j 中的最小横、纵坐标min min ,y x 以及最大横、纵坐标max max y x .以A(min x ,max y ),B(min x ,min y ),C(max x x min y ,),D(max x ,min y )为顶点的矩形区域E 即为所求。

2． E 中取出点p i （按照先固定x ，以步长1使y 增加，当y 取遍区域内所有可能值时，再将x 以步长1增加，直到将该E 搜索完毕）。

3． 用一变量v 记录特殊筛选点(特殊筛选点：与p i 的距离等于最大内切圆半径R的筛选点)的个数，再用一变量a 记录该筛选点的坐标。

4． 如果v ＝2，则p i 可能是中轴点，故用利用上述第二个几何关系进行判定。

5． 即判定这三点是否共线。

如果共线则（说明p i 为中轴点）将p i 的坐标记录到变量z(j).ct 中；否则转2。

6． 如果j 小于100，则j ＝j ＋1，转1，否则转6。

7. 程序结束。

用上述算法得到了每幅切片上的可能中轴点坐标。

在理想情况下，每幅切片上的中轴点只有一个但由于是用离散点代替连续的情况，故每幅切片上的可能有多个点满足被程序判定为中轴点坐标。

故还要对求得的每幅切片的所有可能中轴点在进行选取。

算法如下：1． 取第j(j 的初值为1)幅切片的所有可能中轴点坐标z1(x 1,y 1),z2(x 2,y 2)……z k （有界且N N k ∈），求出所有这些点的横坐标的平均值以纵坐标的平均值（X,Y ）.2． 计算z 1，z 2，……z k 到（X ，Y ）的距离d i （],1[k i ∈）,从中取d i 最小的点作为第j 个切片的中轴点.j=j+1.3． 如果j 不大于100，则转1，否则，程序结束。

由以上算法可得到这100幅切片的中轴点坐标，以及中轴线在XY,YZ,ZX 平面上投影的散点图。

（见附录一）4）结果优化对结果进行分段线性插值：为了能给出较精确的中轴线投影，我们充分利用得到的坐标，对其进行插值。

由于不知到x ，y ，z 之间的具体函数关系，所以我们选择了分段线性插值。

设]100,1[∈k ， k<i t <k+1 2,1=i ……, T(k)=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡)()()(k z k y k x 则T(i t )= (k+1-i t )⨯T(k+1)+( i t -k)⨯T(k)利用上式我们可以得到99个插值点5．中轴线在xy ，yz ，zx 平面上5863764-2793.16701803x 85.2583412 x 0340.96081920-796x 0.00583119 311x 0.00002036-071x 0.000000040004x -0.000000023456++++=y用四次多项式进行拟合，曲线方程为389259683.4100400326908x 3.746630945227120x -0.2504109892569x 0.005489404992775x -0.0000300)(234++++=x z（3）在yz 平面上的投影曲线如下图，红色为散点图该曲线是分段拟合得出的，其中三段是多项式拟合，一段用直线拟合,其表达式如下z(y)=0.00249235673y -2.07086537382y + 573.8394245513y -53004.2749150970(]290,257[∈y ,∈z ]35,1[)z(y)=-0.000025080833y -0.027*********y +9.72082596635y +1154.71951234549 ]421,307[∈y ,∈z ]61,36[)421=y ,∈z ]73,62[z(y)=0.00000022098333y -0.00117224399862y +0.552071633403752y 30.0608345123813]242,421[∈y ,∈z ]100,73[)1）处理明显误差点通过以上的程序算出的中轴点坐标有小部分相邻切面的中轴点坐标是相等的（理论上，这种情况只有在两幅切片图完全相同才会出现），这与实际不符合，所以我们对上述情况进行了处理。