数学建模血管的三维重建问题

1、A 题血管的三维重建问题摘要：本论文讨论基于切片的血管三维重建问题。其背景是：采取存储二维切片信息，使用时再利用切片信息重建原物体三维形态的方法，可以有效地保存和利用三维信息。此技术在实际中有很大的用途，在医学和其他领域有广泛的应用。如要将人体全部三维信息，包含内部错综复杂的结构，完整地存储在计算机中，以现在的技术也是有一定难度的，但若改用存储人体切片信息，使用时重建再现的方法，则是利用现有技术可以解决的。本论文基于题中对血管形态的假设，建立管道中轴线参数方程，并综合考虑实际情况中由于切片厚度及数字图像离散化带来的偏差，通过在每张切片图像中搜索其中阴影区域所能包含的最大圆面，确定管半径为R=2

2、9，在此基础上，将每张切片图像中阴影区域所能包含的半径大于等于R 的圆面圆心作为中轴线与各切片交点（即中心点）的候选点集合。本模型使用了三种改进算法对该候选点集进行筛选以确定实际交点。最终迭代算法简述如下：1对每个切片，建立中心点的候选点集，并取点集的中位点为中心点初值2利用得到的中心点建立中轴线方程3利用中轴线方程推导导数信息，根据导数信息比例选取中心点的候选点集的某点作为中心点的新值4重复步骤2、3，直至结果达到较稳定状态为止5输出中心点及中轴线方程在模型建立中，对选取侯选点集、求中位点、利用导数信息进行比例选取均给出完整的算法，并且对半径确定、候选点选取、采用导数作为比例选取依据等问题给

3、出详尽的证明。考虑到实际血管的中轴线应充分光滑，计算最终中轴线参数表达式时采取了六阶多项式拟合。最后用还原的血管形态模拟切片过程可以得到一系列数字图像，与原切片图像进行比较，可以检验模型的合理性及精度。该模型最终计算结果如下。 血管中轴线示意图从模型结果中看出，中心点分布均匀稳定，模拟检验的切片数字图像与原切片的数字图像吻合较好，模型结果精度及稳定性符合要求。本模型算法简明，理论严密，比例选取算法使结果中心点尽可能收敛于真实中心点，迭代算法保证了结果的精度和稳定性，符合题目要求。利用本模型可解决简单的切片三维重建问题，如应用于在医学、地质、地理等领域进行粗略的分析和三维重建。一、问题的提出及背

4、景断面可用于了解生物组织、器官等的形态。例如，将样本染色后切成厚约1m m 的切片，在显微镜下观察该横断面的组织形态结构。如果用切片机连续不断地将样本切成数十、成百的平行切片，可依次逐片观察。根据拍照并采样得到的平行切片数字图像，运用计算机可重建组织、器官等准确的三维形态。假设某些血管可视为一类特殊的管道，该管道的表面是由球心沿着某一曲线（称为中轴线）的球滚动包络而成。例如圆柱就是这样一种管道，其中轴线为直线，由半径固定的球滚动包络形成。现有某管道的相继100张平行切片图像，记录了管道与切片的交面。为简化起见，假设：管道中轴线与每张切片有且只有一个交点，球半径固定，切片间距以及图像像素的尺寸均

5、为1。计算管道的中轴线与半径，给出具体的算法，绘制中轴线在XY、YZ、ZX平面的投影图。二、问题的分析本问题属于基于切片信息的三维形态重建问题，根据所给的切片数字信息可以直接作图再现血管的三维形态，如下： 图1. 血管三维再现图从上图中基本可以了解血管的三维形态，但由于显微镜的工作原理以及数字图像离散化会造成一定的阴影和误差，其实际形态与再现图是有一定偏离的。根据题中假设，假设某些血管可视为一类特殊的管道，该管道的表面是由球心沿着某一曲线滚动的半径固定的球包络而成，该曲线称为中轴线，球半径也称为管半径。本模型将基于上述假设重建血管的三维形态，使其精确化，参数化。示意图如下： 图2. 血管三维模

6、型图切片数字图像产生偏差的原因有两方面：首先，数字图像离散化造成图像的分辨率有限，如题中给出图像像素尺寸为1。这将造成局部细节被忽略，距离相近的两个点被视为一个点，点之间距离被舍入等等偏差。本模型中根据实际情况假设距离舍入服从四舍五入原则。其次，由显微镜的工作原理可知，一定厚度切片的显微图像是由其中所有组织的成像叠加而成的，而切片厚度是不可避免的，这将造成切片内不同层面图像的干扰，也就是说，每张切片图中的图像实际上是多个层面图像的叠加。然而，题中假设管道中轴线与每张切片有且只有一个交点，即不考虑切片厚度，因此以下叙述也将切片称为切平面，忽略其厚度。由于所给原始数据是512×512单色

7、BMP 图像格式，便于观察，但不便于模型处理，因此首先将图像转化为512×512的二维0-1矩阵i M 形式，其0-1分布与图像对应，1表示黑色像素点，左上角为i M 1, 1(，i M 为记录第i 张切片的图象信息的0-1矩阵，为方便起见，本文中也用i M 来表示第i 张切平面。M 为i M 的集合，记为990, =i M M i 跟据题中建立的坐标系，某点 , , (z y x 的0-1值与矩阵M 中元素的对应关系为：257257( , , (x y M z y x z ，由假设，中轴线与每张切片有且只有一个交点，设该交点为 990(i P i ，根据题中建立的坐标系可以设中轴线

8、参数方称为= ( ( (t Z z t Y y t X x ，其中990 (=t t t Z 由题设也可写为= ( (z Y y z X x ，其中990z 因此中轴线与第i 张切片的交点= (, (i Y i X P i 。三、模型的假设及符号说明模型假设：一.假设样本血管可视为一类特殊的管道，该管道的表面是由球心沿着某一曲线的球滚动包络而成，球半径固定。二.假设管道中轴线与每张切片有且只有一个交点，即将切片视为无厚度的切平面。三.假设切片间距以及图像像素的尺寸均为1。四.假设管道中轴线处处连续，且充分光滑。五.假设两点间距离舍入服从“四舍五入”原则。符号说明：:i M 第i 张切片的位图信

9、息矩阵，为方便起见，有时也称其为第i 个切平面:M 位图信息矩阵集合990, =i M M i (=z Y y z X x ：管道中轴线参数方程，其中990z R ：管半径，即沿中轴线滚动球的半径i P ：中心点，即管道中轴线与第i 个切平面的交点，根据上述方程可知( ( (i Y i X P i ，=，其中990i i SP ：i P 的候选点的集合，为定义在某切平面i M 上所有被判定为最大圆面 (R r 的圆心SP ：i SP SP =，所有候选点的集合四、模型的建立及求解模型的建立：由以上分析，中轴线可以通过确定iP (990i 的坐标来求得。因此问题即转化为如何确定i P (990i

10、 的坐标。由于中轴线为球心滚动轨迹，因此i P 为球心。而理想无厚度 切片的图像为切平面与沿中轴线滚动的球体相交成的圆面叠加而成，即右图中各圆周的包络区域：下面证明的命题对本模型有基础意义。【命题一【命题一】：在该包络区域中可容纳的最大圆面以i P 为圆心，R 为半径。【证明】：假设中轴线上存在另一点P ，以其为球心，以R为半径作球，则该球与此切平面相交成的圆面的半径不大于R ，若为R ，则可知P 与此切平面的距离为0，换句话说，P 在此切平面上。这与题设中轴线与每个切平面有且只有一个交点矛盾。因此可证明在该包络区域中可容纳的最大圆面是以i P 为球心，以R 为半径的圆。图3. 某切面圆面叠加

11、图因此理论上可以通过寻找某切平面中包络区域可容纳的最大圆面来确定中轴线与该切平面的交点位置以及管半径。具体算法如下： , (P M GetMaxR i 用来确定在某位图信息矩阵i M 中以某点P 为圆心的最大圆面半径。流程图见附图五步骤一：令r 0步骤二：若以P 为圆心，R 为半径的圆面包含在位图i M 中，继续下一步骤，否则转步骤四步骤三：r r+1，转步骤二步骤四：返回r-1以下将对i P 的选取和确定作进一步讨论：【算法一【算法一】：将某一切平面i M 所有的点代入 , (P M GetMaxR i ，将求得的半径排序，选取半径最大的圆面圆心作为i P ，循环即得所有中心点坐标及对应最大

12、圆面半径。流程图见附图六。A 题血管的三维重建问题使用该算法求得中心点坐标及对应最大圆面半径见附表一。所求出的i P 构成中轴线的 散点图及拟合曲线见图4：图4. 【算法一】散点图及拟合曲线由假设，血管中轴线充分光滑，由此可用三次多项式拟合散点数据，中轴线参数解析表达式为：9901097. 107. 41077. 11032. 11048. 190. 11038. 31022. 2213322234×+××+×=×××+×=z z z z y z z z x 从上图可以看出，由【算法一】得到的中轴线形态与血管三维再

13、现图基本吻合，但光滑度较差，可以明显看出插值曲线局部存在折点，因此有必要对算法进行改进。在模型计算过程中，注意到在每个切平面上的包络区域中可容纳的最大圆面并不惟一；管半径也不唯一，为29或30。这表面上与【命题一】矛盾，但实质上并不矛盾。最大圆面不唯一是由于数字图像离散化的结果，【命题一】假设坐标及半径均为连续精确的，即使两个相近的数字也是可以分辨的，但实际情况并非如此理想。由于题设切片间距以及图像像素的尺寸均为1，而图像像素是图中可分辨的最小尺度。例如以i P 为球心，29=R 为半径的球与第i+1个切平面相交的圆面半径为：9828. 28129 , (tan 2222=i i M P ce

14、 Dis R ，在数字图像中将被舍入为29，其结果是该圆面也可能被误判为最大圆面，这样所得1+i P 实际为i P 在1+i M 上的投影。从而造成存在多个最大圆面。管半径的不唯一是由于题设管道中轴线与每张切片有且只有一个交点，而实际中切片总是有厚度的，再加上数字图像离散化的影响：对于管半径为29的血管，其切片图像产生1个单位的偏差是符合实际情况的。而且从计算结果也可以看出，出现最大半径为30的切平面的位置处在中轴线沿Z 轴变化较平缓部分，这部分也是最容易产生偏差的。基于以上分析及【算法一】的结果，可以确定管半径R=29。若定义在某切平面i M 上所有被判定为最大圆面圆心的集合为i P 的候选

15、点集，用i SP 表示。调用 , (P M GetMaxR i 可以得到i M 上所有半径R r 即半径为29和30的圆面圆心的集合i SP 。【算法二】：综合某切片内所有i P 的侯选点集i SP 的信息，确定i P 的坐标，可以提高结果的精度。经计算，SP 中有为13596个点。 图5. 从切平面50M 上50P 的候选集50SP 的XY 投影图例如由上图可以看出50P 的候选集50SP 呈带状分布，而【算法一】选取的50P 点处于其候选点集边缘，其效果必然不理想。从算法来看，【算法一】是在50SP 中选取最先遇到的点，弃掉其后的所有点，这种选取原则造成【算法一】结果图中的中轴线产生折点。

16、如果可以将i P 所有的候选点集以某种次序排序，选取其中位点，该点应较接近其真实点i P，其结果相对于【算法一】的取法应有较大改善。具体算法如下： 所求出的i P 坐标见附表，中轴线的散点图及拟合曲线如右图 图6. 【算法二】散点图及拟合曲线A 题血管的三维重建问题由散点数据拟合生成的六次多项式拟合曲线表达式为：990z +×××+×+××=××+×××+××=776 . 0367. 01065. 51032. 31046. 51015. 11076. 51059.

17、 153. 11012. 21012. 11058. 21036. 21040. 7233445566922132445669z z z z z z y z z z z z z x 图7. 中轴线三视图从上图中可以看出，【算法二】结果基本满足题目要求，在模型的检验和结果分析中将对该算法进行具体结果分析及检验。对P 点在其候选区域的相对位置作进一步分析，有以下结论：【结论一】：当且仅当中轴线上某点与某切平面距离小于等于5，其在该切平面上的投影点有可能被识别为中轴线与该切平面的交点，产生误判。若点P 到某切平面i M 的距离5D ，则以点P 为球心，R=29为半径的球与切平面iM 的相交所得圆面的

18、半径5657. 285292222=D R r ，这在离散化的数字图像中将被视为29，该点被误判。另一方面，若3725. 28629, 52222=>D R r D r 将被舍入为28，因此可证明结论。由以上结论可对【算法二】作进一步改进。A 题血管的三维重建问题【算法三】：综合相邻多个切片信息，确定i P 的坐标，进一步提高结果的精度。由【结论一】可知处于某切平面i M 上下五层之间的中轴线段在该切平面i M 上的投影线段L 上的所有点均有可能被判为中轴线与该切平面的交点，这也可从上面图5中P 的候选点集呈带状分布得到初步验证。理想情况下图5中P 的候选点集应呈线状分布，但由于实际中偏

19、差不可避免，而且由图中可以看出实际带状分布与线状分布相当接近。以下讨论将其视为线状分布，设为直线L 。如右图所示，由于中轴线各段的斜率不同，将造成实际中轴线与切平面的交点并不位于L 的中点，斜率大的中轴线部分投影较长，实际交点应靠近斜率小的一端。由于实际离散化的结果，接近的点不被分辨，只有两点距离大于1才会被分辨。因此一段线段将被视为一串离散点，点的个数将取决于线段的长度。如图中，实际L 将被识别为5个离散点，而中轴线与切平面图8. 中位点误差说明图的交点位于第4个交点上。实际交点在点串中位置与其中轴线两端斜率关系叙述如下：如右图所示：若将切平面间的中轴线视为直线段，则若该段中轴线斜率为K ，

20、其在切平面上的投影长度为1/K，因此某切平面M 上下五层之间的中轴线段分别在该切平面M 上的投影线段d u L L , 的长度可以累加求得，d u L L L +=。而实际中轴线与切平面的交点即位于直线L 上的d u L L 比分点上。由此可得图9. 切线分界图【结论二】：实际中轴线与切平面的交点近似位于其候选点集（设其呈线状分布）的d u L L 比分点上。下面给出算法三的具体算法。首先利用【算法二】得到的中轴线来获取各点的导数信息：该曲线的参数方程可表示为：=t z t y y t x x ( (K=Tan( 11/K由前面假设，该曲线处处光滑可导，则在曲线上A 点（坐标为(000, ,

21、z y x 处）切线的方向为00,z y x dtdz dtdy dt dx 。因为1=dt dz ，所以有dzdx dt dz dz dx dt dx =dzdy dt dz dz dy dt dy =则切线的方向为1, ,00y x dtdy dt dx 。设切线与横切面夹角为，如图有：221+=dz dy dz dx tg tg 就是所要求得的斜率。图10. 切线投影图设tg k 1=，有22+=dz dy dz dx k 由于切平面i M 上直线L 理论上是由i M 上下相邻五层内中轴线上的点投影而成的，这样基于【结论一】、【结论二】就得到：+=51i i j ju kL ，=15i

22、i j jd kL ，=+=1551i i j ji i j jd u i kkL L t 我们基于【算法二】得到的中轴线求得每一层交点的k 值，进一步求得分配系数i t ，将位于d u L L 比分点上的分界点作为中轴线与切平面i M 的交点。其效果理论上应较【算法二】有所改善。为保证计算结果精度，采取迭代的方法，将第一步得到的中轴线结果作为输入再次进行运算。由于两次迭代结果十分相近，因此将第二次迭代结果作为最终结果。【算法三】求得的中心点坐标及导数信息见附表3、4，中轴线的散点图及拟合曲线如下图： 图11. 【算法三】散点图及拟合曲线由假设，血管的中轴线应充分光滑，由散点数据拟合生成的六次

23、多项式拟合曲线表达式为：990z +×××+×+××=××+×××+××=723. 0547. 01026. 31033. 11087. 31003. 11039. 51060. 108. 11065. 11041. 91027. 21011. 21064. 6223345566922133445669z z z z z z y z z z z z z x 五、模型的检验及结果分析由于本模型是在题目假设前提下利用切片信息对血管形态进行三维重建，因此模型结果的优劣可以

24、用还原的血管形态模拟得到切片数字图像，将其与原切片的数字图像相比，即可知模型重建结果的合理性和准确性。对于使用【算法二】和【算法三】求出的中轴线与切平面的两组交点 990(21=i P P i i 和，采用直观的模拟法进行检验以证明其合理性并比较两组结果的精度：具体步骤如下：让半径为R=29的球沿i P 所构成的中轴线连续滚动得到血管的三维实体，再用切平面99. 0=i z 截取该实体得到截面i M ，通过对两组截面i M 与i M 的对比，即可判断模型结果正确性和精度。以下是生成的截面i M 与原截面i M 的重合比较图样，由于篇幅限制，只给出较有代表性的几幅。图样说明：原图样为空白，而模拟

25、图样为阴影的部分。记为A：原图样为阴影，而模拟图样为空白的部分。记为 B ：原图样和模拟图样都为阴影的部分。记为C 则表中“不重合率”的计算公式为CB A B A =不重合率【算法二】【算法三】 图样（第5幅）不重合率4.80%4.51% 图样（第35幅）不重合率6.17%5.10% 图样（第45幅）不重合率5.05%4.35% 图样（第55幅）不重合率6.68%6.86% 图样（第65幅）不重合率5.74%6.18% 图样（第75幅）不重合率7.72%7.81% 图样（第99幅）不重合率40.42%40.57%注：出于排版的原因，上表中行间图样比例严格一致，各列中图样比例有放缩。由直观对比以

26、及不重合率分析可以看出：1由于模拟某一切片图象需要中轴线上与其上下相邻距离29之间的各点位置信息，但对于前29个和后29个切片，上述条件是不满足的。因此首尾两端的切片只能够重现出部分图象，这从第99张图样中可以明显看出。因此对模型的精度评价将主要考虑中部的再现，首尾仅作参考。2两种算法在生成中心点时都存在一定的误差，基本上误差都可以被控制在5%左右。3从【算法二】和【算法三】的对比中可以看出，【算法三】较【算法二】有所改善，但相差不大。但从第1幅图样可以看出，【算法三】生成的模拟切片图样都被限定在原图样的阴影部分内，而【算法二】则有一部分超出原图样阴影范围。这说明【算法二】生成的中心点较实际中

27、心点可能会有较大的偏移，而【算法三】使用中轴线的导数信息，对中心点的选取也更为合理和稳定。所以【算法三】较【算法二】更优。六、模型的进一步讨论由于原始切片数字图像精度只能达到1，理论上应有5. 0±的误差允许范围，因此模型结果中半径应处于5. 029±范围内。此外模型结果中求得中轴线与各切平面交点均为整数坐标，这将造成舍入误差。事实上模型中假设中轴线是过切面的某个象素点的。实际情况中，中轴线可能不过任何象素点，而是从相邻的象素点之间穿过，但理论上模型结果给出中轴线应充分接近实际中轴线，因此在不超过原始数据精度条件下，模型精度总是可以达到要求的。根据血管形态特点，其中轴线应该

28、是连续并充分光滑的，但求得的离散点插值结果不是很光滑，因此有必要进行多项式拟合，以得到充分光滑的曲线，这样处理之后拟和中轴线应更加接近实际中轴线的形态。七、模型的评价优点：该模型融合多个切片的信息来确定血管的中轴线，可以较精确地从血管的切片数字图像中得到血管中轴线和管半径。有较高的稳定性和精度。本模型所采用的算法简明且易于实现，效率较高，可以在题目所允许的范围内保证较高的精度，并且采取了仿真模拟对比的方法对模型进行了检验，取到了较好的效果。缺点：该模型在求i SP 时花费时间较大。综上所述，本模型基本上较好地解决本题假设下基于切片数字图像的血管的三维重建问题，结果稳定性和精度令人满意。参考文献

29、：1VisualC+数字图象处理第12页2.3.1节。何斌等编著，人民邮电出版社，2001年4月。2计算机图形学及其在工程中的应用王莉等编著，人民交通出版社，1988年10月。3应用Matlab建模与仿真陈桂明等编著，电子工业出版社，2000年5月。完成时间：2001年9月24日星期一附录附表1：算法一：管道中轴线与各切平面交点zP (990z X Y Z R X Y Z R X Y Z R XYZ R -1620029-161192530-145765029121657529-1620129-161192630-141835129161657629-1620229-160282730-138

30、885229271647729-1620329-160282830-136915329291647829-1620429-160272930-131985429371637929-1620529-160283030-1281025529481618029-1620629-160283130-1221095629571598129-1620729-159353230-1181135729671568229-1620829-159353330-1141175829751538329-1621929-159353430-1061245929801518429-16211029-158403530-1

31、051256029871488529-16211129-157443630-1001296129931458629-16211229-156483730-9613262291041398729-16211329-156483830-8713863291131338829-16221429-153573930-8414064291231268929-16221529-152604030-7714465291301209029-16231629-158404129-7314666291411099129-16231729-157444229-6714967291441069229-16241829

32、-156484329-5515468291491009329-16241929-154544429-491566929154949429-16252029-153574529-391597029157909529-16262129-152604629-311617129161849629-16262229-150654729-251627229163819729-16272329-149674829-141647329167749829-161202430-147724929-31657429172659929-16229-161192530-145765029121657529算法二：管道中

33、轴线与各切平面交点zP (990z X Y Z X Y Z X Y Z X Y Z -16200-1612225-113117507415375-16201-1602426-106123518015076-16102-1612827-102128528614877-16103-1593028-96133539414678-16214-1603329-911365410014079-16215-1583530-861405510813780-16216-1593831-791435611513381-16117-1584132-731455712112882-16118-1574533-6814

34、85812912283-16129-1554934-631525913611484-162310-1555335-561536014710485-161311-1545636-50155611519886-162412-1516137-42159621579187-162513-1516438-35161631618488-161614-1506839-26163641657889-162715-1477140-16164651687190-162816-1467541-7165661726691-162917-14377423166671736192-1611018-142814311165

35、681765693-1611119-139864421164691795094-1621220-136894531165701804495-1621321-133944640162711823796-1621522-1291004750160721853197-1611723-1271044860159731862498-1611924-1211094966157741861599-16200-1612225-113117507415375-1621-1602426-106123518015076算法三：第一次迭代管道中轴线与各切平面交点zP ，K(990z X Y Z K X Y Z K X

36、 Y ZKX Y Z K -162005.769178-16119253.248884-14576506.09655512165757.046504-162015.34368-16119263.438636-14183516.13204316165767.166938-162024.930488-16028273.623732-13888526.16335427164777.300573-162034.530246-16028283.80349-13691536.19086929164787.447924-162044.143798-16027293.977403-13198546.21500

37、537163797.609433-162053.772266-16028304.145095-128102556.23620848161807.785464-162063.417156-16028314.30629-122109566.25496257159817.97631-162073.080516-15935324.460788-118113576.27178667156828.182189-162082.765149-15935334.608448-114117586.28723975153838.403255-162192.474901-15935344.74918-10612459

38、6.30191580151848.639595-1621102.214986-15840354.882932-105125606.31645287148858.891247-1621111.992247-15744365.009686-100129616.33152393145869.158193-1621121.815012-15648375.129454-96132626.347843104139879.440378-1621131.692004-15648385.242273-87138636.366161113133889.737709-1622141.629878-15357395.

39、348201-84140646.3872631231268910.05006-1622151.630133-15260405.447322-77144656.4119651301209010.37729-1623161.687647-15840415.539737-73146666.4411121411099110.71924-1623171.792288-15744425.625569-67149676.4755661441069211.07572-1624181.932273-15648435.704958-55154686.5162061491009311.44656-1624192.0

40、96853-15454445.778068-49156696.563917154949411.83156-1625202.277431-15357455.845079-39159706.619578157909512.23052-1626212.467594-15260465.906195-31161716.684058161849612.64327-1626222.6627-15065475.961641-25162726.758202163819713859429-14967486.011665-14164736.842821167749813.50932-16

41、120243.055397-14772496.056537-3165746.938683172659913.96227-1625.769178-16119253.248884-14576506.09655512165757.046504算法三：第二次迭代管道中轴线与各切平面交点zP ，K (990z X Y Z KX Y Z KX Y Z KXYZ K2579505.89918327697253.221984373142506.113611411327757.0669972579515.46127727897263.417567378148516.147776409334767.1933252579625.03558428396273.608024383153526.177642407341777.3335492579634.62267828696283.792697387160536.203611404349787.4881752589544.22331228797293.97110439