全国大学生数学建模竞赛一等奖血管切片的三维重建

全国大学生数学建模竞赛一等奖----血管切片的三维重建血管切片的三维重建王月娇钱志刚刘磊指导教师花强（河北大学数学与计算机学院，保定，071002）摘要：为了利用血管切片图象重建血管的三维形态，我们首先编程对切片图象进行由BMP数据格式向文本格式的转换，而这样获得的数据文件将较大，共约50M。

显然在利用计算机做进一步读写与处理前，要着重面对的问题是对大量数据的处理，而其中却有大量冗余信息存在，于是我们在寻找每张血管切片中心轴点的过程中，利用多种优化算法以简化问题，并确定出100个中心轴点。

以此100个中心轴点为样点，依靠样条插值，利用Matlab 软件对分别对平面及空间曲线进行插值，先后建立了模型一、模型二，从而拟合得到XY，YZ，ZX平面的投影曲线以及中心轴线，得到血管半径为29个像素单位，并模拟给出血管的三维空间形态(如下图)，进而对结果进行了广泛的分析与评价。

同通过矩阵投影到XY，YZ，ZX平面而获得的图象进行比较，可以验证模型的正确性，以及模型的先进性。

血管切片的三维重建图一、问题重述断面可用于了解生物组织、器官等的形态。

例如，将样本染色后切成厚约1m 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平面的投影图只需令Z=0，X=0，Y=0就可以得到。

三、模型假设1、假设血管管道的表面是由球心沿着某一曲线（中轴线）的球滚动而成的，也就是血管管道半径（即球体半径）固定且中轴线是连续的。

2、我们很合理的假设所有数据均是准确的。

用大量象素能够非常近似地描绘一个图形。

3、因血管管道可以看作是一个半径固定的球体滚动包络而成的，因此我们认为中轴线，中轴线的一阶导数及二阶导数都是连续的。

4、数据精确到单位像素。

因切片厚度为1uM，因此像素单位亦为１uM，对于５１２＊５１２的信息距阵来说，此精确程度已足以解释问题。

5、切片与血管中轴线的交点存在且唯一。

四、符号说明r：表示血管管道半径。

ri：表示第i张切片中计算得到的管道半径。

ro：表示r的初始估计值（ro<r）。

Sjk：表示第 i张切片中第j个内点到第k个边界点的距离。

Sj ：表示第i张切片中第j个内点到边界的最短距离。

(xi,yi,zi)：表示第 i张切片得到的最大内切圆的圆心坐标。

五、模型构成首先我们用第一个C程序把BMP图象文件格式转换为保存了图象点阵信息的文本文件。

每个BMP文件对应一个文本文件（从00.txt到99.txt），每个文本文件均为512行，512列，按图象的视觉直观顺序，以0表示白色象素点， 1表示黑色象素点( 如右图)，以此方便后续程序的读写及人为地纠错。

接着，第二个程序实现了各个球心坐标求取。

确定每张切片与中轴线的交点的坐标及半径：按理论来说，血管截面边界是圆滑曲线，对一个球体过球心切割，无论怎样切，都必得到一个大圆，由假设5，每个切片包含且仅包含一个这样的大圆，即它的最大内切圆，圆心即球心，圆半径即球半径。

下面我们确定每张切片中最大内切圆的圆心。

对第i 张切片上的第j 个内点，求到边界点 k的距离Sjk,从中选取一个最小距离Sj,=min{Sjk}再从中选取一个最大的,记为ri=max{Sj}.不难理解，这就是第i张切片的最大内切圆的半径。

相应的内点即为圆心（中轴线与截面的交点）。

对100张切片搜索后就得到100个球心点。

但我们要面临的问题是切片是由大量的象素点近似描绘出来的，会给我们的计算带来一定的误差，但是由于象素很小，误差就不会很大，我们依然可以应用上述理论，人为的进行误差分析和修改，更准确地得到球心(x y z)及由假设4可知半径r=29个像素单位（29=<ri<=29.69）。

X Y Z X Y Z X Y Z X Y Z X Y Z(-161 0 0)(-161 0 1)(-161 0 2)(-161 0 3)(-161 0 4)(-161 0 5)(-161 0 6)(-161 1 7)(-161 1 8)(-161 1 9)(-161 2 10)(-161 2 11)(-161 2 12)(-161 4 13)(-161 5 14)(-161 6 15)(-161 8 16)(-161 10 17)(-161 13 18)(-161 17 19)(-161 18 20)(-161 19 21)(-161 20 22)(-161 20 23)(-161 20 24)(-161 20 25)(-161 20 26)(-160 29 27)(-160 30 28)(-159 35 29)(-159 35 30)(-159 35 31)(-158 40 32)(-157 44 33)(-156 48 34)(-155 51 35)(-156 48 36)(-156 48 37)(-152 60 38)(-150 65 39)(-150 65 40)(-138 88 41)(-136 91 42)(-136 91 43)(-136 91 44)(-136 91 45)(-136 91 46)(-119 112 47)(-118 113 48)(-117 114 49)(-116 115 50)(-115 116 51)(-114 117 52)(-113 118 53)(-112 119 54) (-104 126 55)(-96 132 56)(-71 147 57)(-71 147 58)(-60 152 59)(-60 152 60)(-46 157 61)(-20 163 62)(-20 163 63)(-20 163 64)(-13 164 65)(-13 164 66)(-13 164 67)(38 163 68)(43 162 69)(48 161 70)(53 160 71)(60 158 72)(60 158 73)(67 156 74)(75 153 75)(80 151 76)(87 148 77)(87 148 78)(116 131 79)(119 129 80)(131 119 81)(131 119 82)(132 118 83)(144 106 84)(144 106 85)(145 105 86)(145 105 87)(151 98 88)(151 98 89)(163 81 90)(166 76 91)(175 58 92)(175 58 93)(177 53 94)(180 45 95)(180 45 96)(181 42 97)(183 35 98)(184 31 99)在此之后我们找到两种求中轴线的方案，并对其进行分析。

模型1：把从切片中得到的100个球心分别投影到XY,XZ,YZ平面，直接应用三次样条插值中Ｍ连续性方程如下：由假设3知一阶导数连续，又可得到n-1 个等式，则可求解。

相应地求出三条三次样条插值曲线Y(z),Y(x),Z(x)沿垂直于它们所在平面的方向扩展到空间，将得到三个曲面，它们两两相交得到三条曲线。

如果三条曲线能够重合，那无疑就是我们要求的中轴线，但是由于曲线Y(z),Y(x),Z(x)都是三次样条插值拟合出的曲线。

在插值过程中，由于分别利用空间点的投影数据，因此插值结果将分别丧失一维信息，这样由投影点拟合的曲线，将极有可能不能还原成三维结果，即使还原也不会完全通过100个先前得到样本点，因此误差将会较大。

如果能找到一个准则将三条曲线合并成一条，问题就解决了，但是三条曲线中的任意一条都满足过样本点且处处连续，我们已有的数据只有样本点，所以找不到这样的一个准则。

那么我们只能任意选取两条曲线沿它们所在平面的垂直方向扩展得到的曲面相交得到的曲线作为中轴线。

这样我们可以解决重建血管的三维形态的问题，但是它存在着很大弊端，由上面的分析我们知道，因为中轴线是由两条三次样条插值曲线分别沿着它们所在平面的垂直方向扩展的两个曲面相交而得到的，而不能由那三个曲面的两两相交得到的三条曲线重合得到，也就不可避免地存在较大误差，这种方案也就存在着一定的局限性。

因此，在此基础上我们提出第二种方案。

模型2：当100个中心轴样点坐标已知后，考虑采用样条插值将已知样点光滑连接，但实际问题是一空间曲线的插值问题，因此与一般的一维插值有技术上的差别。

我们将Z轴坐标z视为x及y坐标的参数，当z做单调变化时,则x(z),y(z)可分别看作相应此参量的横、纵坐标方向的变化量，这样我们就可以利用一维样条插值思想将空间曲线模拟出来。

当中轴线确定后，再分别令x=0，y=0，z=0，便可以得到中轴线在三个坐标平面上的投影，(见下图)而且还可在中轴线上取充分多点，利用Matlab画出半径为29的球，由此而形成的球体包络线便可展现三维血管的空间形态，且能得到分别在三个坐标面上的投影。

（见后图）六、算法分析1、估计圆心我们可以用计算机搜索第i张切片半径ri=maxj{mink{Sjk}}，记录下该内点就是此切片中最大内切圆的圆心，即中轴线上的一点，最小距离中的最大值就是据此切片得到的管道半径。