基于MATLAB的血管三维重建及模型检验

三维血管重建优秀论文含源代码Prepared on 21 November 2021血管的三维重建摘要本文探讨血管的三维重建，由血管的相继100张平行切片图像计算血管的中轴线与半径，并绘制血管在三个坐标平面上的投影。

由于血管的表面是由球心沿着某一曲线（称为中轴线）的球滚动包络而成，由此我们得出结论：每个切片一定包含滚动球的大圆，并且它一定为切片的最大内切圆，而最大圆所对应的半径即为血管半径，所以求血管半径就转化为求每一个切片内部的点到切片外部轮廓的最大半径。

首先，读取100张血管切面图，把它们转换成Logical矩阵，从中提取切片截面轮廓点构成一个新的矩阵。

然后找到原图片矩阵中像素点的内点（切片图片中轮廓线中的点），从而得到内点到切片轮廓点的最小距离矩阵和最小距离中的最大值矩阵，最大值即为血管半径。

最后计算所有切片的血管半径，并对这些半径求平均值，得到平均血管半径为：μm。

由100张切片的最大内切圆圆心坐标拟合得出中轴线方程以及其在三个坐标平面内的投影曲线方程。

由中轴线得到血管的三维立体重建图，用平面)74=iiZ去截血,0(,=2449,,管的三维立体重建图，得到新的4张截面图。

把它们分别与题设中的对应截面进行内点个数对比。

我们定义两张切片所共同拥有的内点个数与原切片内点个数的比值为重合度。

计算得到平均重合度为：% 。

关键词：血管半径中轴线切片重建（来自作者：欢迎各界人士批评指正，。

文章作于2011年8月10日，陕西科技大学理学院实验室）1问题的重述断面可用于了解生物组织、器官等的形态。

例如，将样本染色后切成厚约1μm的切片，在显微镜下观察该横断面的组织形态结构。

如果用切片机连续不断地将样本切成数十、成百的平行切片，可依次逐片观察。

根据拍照并采样得到的平行切片数字图象，运用计算机可重建组织、器官等准确的三维形态。

假设某些血管可视为一类特殊的管道，该管道的表面是由球心沿着某一曲线（称为中轴线）的球滚动包络而成。

血管的三维重建摘要本文对于血管的三维重建问题，通过分析给定图片的像素数据，通过以下方法完成了血管的三维重建。

首先我们假定血管为等径管道，将血管看作半径不变的球沿着一条中轴线滚动形成。

根据题目中提供的100张BMP图片的像素点坐标数据，通过分析其几何特性，给出寻找其中轴线和半径的方法，由此找到每张图片中最大内切圆的圆心和半径，利用平均法，得到半径的平均值约为29.25 m.，即为管道半径。

接下来利用100个圆心的坐标利用最小二乘法拟合得到中轴线的方程：x t=7.9021t7−1.6321t6+0.1378t5−0.0056t4+0.0001t3y t=−2.9429t7+0.6330t6−0.0546t5+0.0023t4−0.0001t3z t=t通过对圆心坐标点绘制出的曲线与拟合后的曲线进行对比，发现拟合的曲线与得到的圆心坐标点非常吻合，从而更好的反映了中轴线上各点。

然后再用枚举法以得到的圆心坐标点为球心，画出100个球面上的点组成的还原图与100张切片的边界叠加成的还原图进行比对，形状相符，成功的还原了血管的三维图像。

关键词：血管；三维重建；像素；还原；最小二乘拟合；滚动法；枚举法1.问题的重述断面可用于了解生物组织、器官等的形态。

例如，将样本染色后切成厚约1 m的切片，在显微镜下可以观察该横断面的组织形态结构。

如果用切片机连续不断地将样本切成数十、成百的平行切片，可依次逐片观察。

根据拍照并采样得到的平行切片数字图象，运用计算机可重建组织、器官等准确的三维形态。

假设某些血管可视为一类特殊的管道，该管道的表面是由球心沿着某一曲线（称为中轴线）的球滚动包络而成。

例如圆柱就是这样一种管道，其中轴线为直线，由半径固定的球滚动包络形成。

现有某管道的相继100张平行切片图像，记录了管道与切片的交。

图像文件名依次为0.bmp、1.bmp、…、 99.bmp，格式均为BMP，宽、高均为512个像素（pixel）。

血管三维重建的数学模型刘俊健、吴友玲、徐勇哲[摘要]:本文探讨了血管的三维重建问题.我们首先利用数字图象处理的相关技术,从100张图中读取出相应数据,然后以定位所需的相关参数(包括中心轴线和半径)为目标,主要运用了包络面相关理论及几何拓扑网络中的相关知识得出了如下结论：切片一定包含一个滚动球的大圆,并且它恰好是切片内的最大圆.运用这些结论设计出有效的算法,利用数学软件(Matlab 及Maple)和计算机技术,对100个切片编程求解,找出切片所包含的最大圆,即得到中心轴线与100个切片交点的坐标及大圆半径.经过误差分析确定出管道半径为69849.29=r ,并给出了中心轴线上点的坐标,最后我们还绘制出中心轴线在zx yz xy ,,平面上的投影图.关键词：三维重建;最大圆;中轴线;有障碍物的最大空隙问题1 问题的提出生物组织连续切片的计算机三维重建技术是把一系列切片的图象，通过计算机进行处理,从而得到该生物组织立体结构的一种方法.在切片图象内部寻找定位结构,是此核技术的关键,同时也是目前国内外三维重建技术的共同难题.解决此难题的一个重要方向是通过提取切片中二维图象的数据信息,建立模型,以达到定位效果.假设某类血管表面可视为由球心沿着一曲线（称为中轴线）的球滚动包络而成.现有一此类血管相继100张平行切片图象,记录了管道与切片的交.图象格式为BMP,宽、高均为512个象素（pixel ）.假设管道中轴线与每张切片有且只有一个交点,球半径固定,切片间距以及图象象素的尺寸均为1.如果取坐标系的Z 轴垂直于切片,第1张切片为平面0=Z ,第100张切片为平面99=Z .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) 假设血管为实心体,其表面为一包络面；2) 假设管道中轴线与每张切片有且只有一个交点,滚动球的半径固定； 3) 假设切片间距以及图象的象素尺寸均为1； 4) m 表示中轴线,是一空间曲线； 5) r 表示滚动球的半径；6) S 表示球心沿中轴线滚动形成的包络面；7) V 表示由曲面S 与平面0=Z ,99=Z 所围成的管道体； 8) i F 表示管道体V 与平面i Z =的切片； 9) i e 表示i F 的边界.没有说明的符号在文中初次引用时均作出了说明.3 问题的分析与模型的建立首先分析题目中所给出的100张图片.对这100张黑白图片通过计算机图象处理,我们可以观察到其切片呈现一个旋转渐变的过程,如图（一）因此我们估计该血管的形状类似一条螺旋管.只要将这一百张切片在空间中定位,由于切片间距非常小（为一个象素）,所以就可以近似地还原出血管的三维结构.现先对每张图作出一些数据上的分析.针对BMP 图象的存储格式,我们将所给的图片转为二值图象处理.所谓二值图象,就是只有黑白两个灰度级,即象素灰度级非1即0.我们可以从每张图中读取出一个512阶的0-1矩阵.其中0表示象素为黑色.1表示象素为白色.我们可以给出矩阵元素的位置与题目中所定坐标系之间的坐标变换.即矩阵的第i 行,第j 列的元素坐标相应转化为)256,257(i j --.现在我们把问题转化为数学语言.在笛卡尔坐标系内,有一半径为r 的球,沿空间的某一曲线m （称为中轴线）,滚动包络形成空间曲面S ,现有100张平面图,99,,1,0===Z Z Z ,其中0Z ,99Z 与曲面S 围成一实心管道体,现题目就是要求我们求出该管道体的中轴线m 及滚动球的半径,并要给出中轴线m 在XY,YZ,ZX 平面上的内投影图.现设管道V 中轴线m 的方程为：⎪⎩⎪⎨⎧===)()()(t z z t y y t x x 其中99)(,0)(],,[990990==∈t Z t Z t t t 由于管道体中心轴线m 与每张切片有且只有一个交点,则不妨设中心轴m 与切片)990( =i F i 有且只有点)),(),((i t Y t X O i ,i Z 与1+i Z 的间距为1个象素,在如此小的间距下,我们可以用100个交点i O 近似地刻画中轴线.结论1：管道体V 与平面i Z =的切片i F 一定包含一个过滚动球球心,半径为r 的圆A . 现证明如下：以中心轴m 与平面i F 的唯一交点)),(),((i i Y t X O i 为球心（亦即包含一个滚动球的大圆）,r 为半径所作出的球A 与平面i Z =交于圆A ,其中圆A 的半径为r ,由于球A 在管道体V 内,所以球A 与平面i Z =的切片也在管道体与i Z =的切片内,亦即圆在i F 内,结论1得证.结论2：管道体V 与平面i Z =的切片i F 内的最大圆恰好是滚动球中过球心半径为r 的圆A ,并且最大圆唯一.现证明如下： 由中轴线m 的方程：⎪⎩⎪⎨⎧===)()()(t Z z t Y y t X x ],[990t t t = )1( 其中99)(,0)(990==t Z t Z . 有曲面S （即包络面）的方程来：⎩⎨⎧=-+-+-=-+-+-0)('))(()('))(()('))(())(())(())((2222t Z t Z z t Y t Y y t X t X x r t Z z t Y y t X x )2(现考虑i Z =,取't t =使得i t Z =)'(则⎩⎨⎧=-+-=-+-0)('))(()('))(())'(())'((0000222t Y t Y y t X t X x r t Y y t X x )3( （3）式的几何意义是,在切片i F 的边界曲线i S 上的点与中轴线和切片的交点的最近距离为r ,即可作以))'(),'((t Y t X 为圆心,半径为r 的圆.边界曲线i S 和中轴线m 与切片i F 的交点)),(),((i t Y t X O i 的距离达到最大值r .又由于每一个i 属于集合{}99,,1,0 ,i z =与中轴线⎪⎩⎪⎨⎧===)()()(t Z z t Y y t X x 只有唯一交点,故只有唯一点t 满足i t Z =)(,所以切片包含的最大圆唯一.结论2得证.至此,我们只要在每个切片上找出一个能被切片i F 包含的最大圆的圆心,即为该切片与中轴线的交点.就一般包络面而言,切片的最大圆应是唯一的,但由于数据存储误差,我们可能会找到多于一个的最大圆.所以我们需要一个能找出多人最大圆的可行算法.利用以上的分析和结论,我们试图求出管道体半径r 及每个切片与中轴线的100个交点,由于BMP 存储方式的限制,每个切片实际是由有限多个点组成.我们利用几何拓扑网络中的有障碍物的最大空隙问题解法,把障碍物视为半径为0的圆,因此我们可以在切片上作如下约束优化：222)()(R s.t.Rmax i i y y x x ---=其中,i i i i F y x S y x ∈∈),( ),(i F 为管道体切片,i S 为i F 的边界,且i F 、i S 均为限点集.上述模型可简化为如下：))()((22i i y y x x Min Max r -+-=其中i i i i F y x S y x ∈∈),( ),( i F 为管道体切片,i S 为i F 的边界,且i F ,i S 均为限点集,关于r 的求解和圆心坐标,算法及求解结果将于下部分给出.4 模型的求解与算法的实现我们根据上面的在切片上进行约束优化思想,利用计算机编程来寻找最优解.由于切片图象是用BMP格式存储的,因此它都有是由一些有限的点构成,我们用Matlab读入每张图,用一个512⨯512的0-1矩阵表示,其中0表示象素为黑色的点.去掉所有非零元素,并记下所有零元素在矩阵中出现的位置,再转化为题中的笛卡儿坐标系的坐标,这样,就可以用描点方法作出切片的用点集表示的图,如果我们对100个切片上的点都在空间中描出来,那么就可以大体反映出管道体的三维结构.受图片存储的精度所限,只能在表示切片的点集上寻找最大圆心.我们根据上节中所列出来的约束优化思想,再利用计算机编程的方法来寻找最大圆.首先我们要作出切片的边界.由于图上的点都是呈阵列排布的,只要依次取切片上的点,选出与1相邻的0点即为切片的边界.现在按如下算法寻找最大圆的圆心和半径:1．取出切片点集中的一个点,再算出该点与边界上各点的距离,找出最小值,记为h,这就是以该点为圆心所能作出的包含于切片的最大圆的半径.2．继续取点集中的下一个点,再算出该点与边界上各点的距离,再找出最小值与 h进行比较,若比h大则把值赋给h,若与h相等则把值赋给h并记下坐标,同时还把这个点都记起来,这样就保证若最大圆不唯一,也可把所有最大圆找出来.3．取遍切片点集上的所有点运行第二步,最后得出的h即为该切片中最大圆的半径,相应该点的坐标即为切片与中轴线交点的坐标.()．源程序见附录我们记每张切片上点的数目为n1,切片边界上点的数目为n2,显然n2远小于n1,我们要进行n1重循环,每一次循环作n2次比较,其复杂性不超过O()21n,所以,此算法为有效算法.我们用Matlab编程进行搜索,利用Matlab强大的矩阵运算功能,在编程时进行向量化,大大提高了程序的效率,对100张图分别算出如下半径和交点的坐标：注：上表中（X,Y,Z）表示与中轴线交点的笛卡儿坐标.R表示在切片上最大圆半径.由于图象精度问题及算法的误差.每张切片上算出的半径有微小差异,我们根据误差产生的原因进行分析（详细分析过程下节给出）,取出所有半径中最大值r=29.69849,即为题目所要求的半径,把所得的100个中轴线上的点在空间中描出来,得到结果如图（二）所示：图（二）为了不影响最终结果的客观性和全面性,我们对数据可以不作任何拟合,以图表的方式给出结果,以供有需要人士随时查阅.我们已经求出中轴线上100个紧密的点,因此中轴线在XY,YZ,ZX平面上的投影应该是一些紧密的点,为直观起见,我们把这些点用线段连接起来,得到的投影图如下：在XY平面的投影在XZ平面的投影在YZ平面的投影５ 结果分析由于BMP 图象是用点阵来表示的,用它读取出来的数据只能近似还原出原来的结构,因此得到的中轴线上的点必然有一定的误差,不过由于中轴线上的点非常密,作成的图表的效果也很不错.对于得到的100个半径的数据,我们取其中的最大值作为管道的半径,是基于如下的误差分析：由于程序算法可知,我们得到的最大圆的半径是切片点集中的所有点到边界点集的最短路距离中的最大值,所以由该点作出的圆不会超出切片的范围内,但由于我们找最大圆时,圆心都定在象素点上,而实际切片包含的最大圆的圆心可能不在象素点上,如图（三）所示：在正方形点集图上,每一个交点代表一个象素.我们取象素点为圆心找最大圆,只能找到半径为1的圆,其中中间的四个交点都可能是圆心,它实际包含的最大圆的半径应为1.5,圆心落在大正方形的中心点（非象素点）上.所以,当所得的大圆与切片只有一个切点时,大圆半径一定仍可以扩大,但所扩大的长度不会超过一个象素.而且有两段边界出现平行时,会出现不止一个最大圆圆心.而我们在程序中也考虑了进去,但并没有发现有出现这种情况.因此,只要提高图象的分辨率,便可减少误差.但由于点与点之间总有间隔.分辨率不可能无穷大,所以误差总是存在的.因此对于我们得出的100个半径值,我们取出其中最大值便是误差最小的,而我们这批数据找到的最大值为29.69849,最小值为28.86174,相差为0.83675,误差没有超过一个象素.CD图（三）６模型的改进与推广考虑到生物结构的完整性和连续性,我们当然希望得到的100个点能连成较为光滑的空间曲线,但实验数据在坐标平面上往往表现为一些离散的点.直接用这些观测点依序连结的折线来反映其变化情势常常过于粗糙.因此需要给出光滑的实验曲线,通常要求这条实验曲线与观测数据的残差平方和最小.对于本模型中的数据点,我们可作进一步的调整,从误差分析,我们可以知道许多切片上取得“最大圆”实质并非真正最大圆,但它必与已切片相切于一点,而与该点切线垂直的方向上的圆点,与切片的点相距应小于１个象素,所以其圆心可以适当向这个方向微调.经过调整修复后,我们还原出近似的三维血管图象如下：在误差分析时,我们发现误差很大程度上是由图象存储的精度决定的,只要数据足够精确,我们的模型将会得到更准确的结果,误差将会大大的减少.随着计算机技术的飞速发展和分辨电子显微学的发展,我们必能得出更为完善的结果.参考文献1.周培德、卢开澄. 计算几何. 北京：清华大学出版社.19992.夏良正.数学图象处理.南京：东南出版社. 1998.83.张宜华、史惠康.精通MA TLAB5. 北京：清华大学出版社.1999,６4.李世奇、杜梦琴.计算机代数系统应用程序设计.重庆：重庆大学生出版社.1999.45.数学手册编写组.数学手册.北京：高等教育出版社.19796.何昆等.生物组织细胞连续切三维重建的研究进展. 军事医学科学院刊. 2000年第24卷北京.2000.3 编辑:许智杰方创彬。

A题血管的三维重建问题摘要：本论文讨论基于切片的血管三维重建问题。

其背景是：采取存储二维切片信息，使用时再利用切片信息重建原物体三维形态的方法，可以有效地保存和利用三维信息。

此技术在实际中有很大的用途，在医学和其他领域有广泛的应用。

如要将人体全部三维信息，包含内部错综复杂的结构，完整地存储在计算机中，以现在的技术也是有一定难度的，但若改用存储人体切片信息，使用时重建再现的方法，则是利用现有技术可以解决的。

本论文基于题中对血管形态的假设，建立管道中轴线参数方程，并综合考虑实际情况中由于切片厚度及数字图像离散化带来的偏差，通过在每张切片图像中搜索其中阴影区域所能包含的最大圆面，确定管半径为R=29，在此基础上，将每张切片图像中阴影区域所能包含的半径大于等于R的圆面圆心作为中轴线与各切片交点（即中心点）的候选点集合。

本模型使用了三种改进算法对该候选点集进行筛选以确定实际交点。

最终迭代算法简述如下：1．对每个切片，建立中心点的候选点集，并取点集的中位点为中心点初值2．利用得到的中心点建立中轴线方程3．利用中轴线方程推导导数信息，根据导数信息比例选取中心点的候选点集的某点作为中心点的新值4．重复步骤2、3，直至结果达到较稳定状态为止5．输出中心点及中轴线方程在模型建立中，对选取侯选点集、求中位点、利用导数信息进行比例选取均给出完整的算法，并且对半径确定、候选点选取、采用导数作为比例选取依据等问题给出详尽的证明。

考虑到实际血管的中轴线应充分光滑，计算最终中轴线参数表达式时采取了六阶多项式拟合。

最后用还原的血管形态模拟切片过程可以得到一系列数字图像，与原切片图像进行比较，可以检验模型的合理性及精度。

该模型最终计算结果如下。

血管中轴线示意图从模型结果中看出，中心点分布均匀稳定，模拟检验的切片数字图像与原切片的数字图像吻合较好，模型结果精度及稳定性符合要求。

本模型算法简明，理论严密，比例选取算法使结果中心点尽可能收敛于真实中心点，迭代算法保证了结果的精度和稳定性，符合题目要求。

血管的三维重建摘要对于血管的三维重建，本文研究了血管这一类特殊管道的中轴线及其半径的算法，绘制中轴线在XY 、YZ 、ZX 平面的投影图这些问题，问题分为三部分。

针对第一部分，先将100张切片图片在MATLAB 中导出生成0-1矩阵数据，在计算100张切片的最大内切圆半径及对应圆心坐标，为减小误差求100张切片最大内切圆的平均半径41666.29 d 。

中轴线的曲线方程可在MATLAB 中拟合得到。

针对第二部分，得到中轴线曲线方程在MATLAB 中绘制出中轴线方程的空间曲线，之后将其投影在XY 、YZ 、ZX 平面上。

针对第三部分，对100张切片进行叠加重合，得到血管的三维立体图，再通过MATLAB 对血管的三维立体图进行优化完成血管的三维重建。

关键词：MATLAB 软件管道半径中轴线曲线方程一、问题重述1.1基本情况断面可用于了解生物组织、器官等的形态。

如果用切片机连续不断地将样本切成数十、成百的平行切片，可依次逐片观察。

根据拍照并采样得到的平行切片数字图象，运用计算机可重建组织、器官等准确的三维形态。

1.2相关信息假设某些血管可视为一类特殊的管道，该管道的表面是由球心沿着某一曲线（称为中轴线）的球滚动包络而成。

现有某管道的相继100张平行切片图象，记录了管道与切片的交。

图象文件名依次为0.bmp、1.bmp、…、99.bmp，格式均为BMP，宽、高均为512个象素（pixel）。

取坐标系的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）。

1.3提出的问题问题一：计算出管道的中轴线与半径，给出具体的算法。

血管的三维重建模型摘要：本文对血管三维重建中，中轴线及球的半径确定问题进行了讨论。

首先，根据问题及图象处理提取有效数据，给出两种可行算法，利用上述数据建立了最大最小方法和二次规划方法。

搜索中心点，并给出全局和局部搜索，得到各切片中心点坐标(见表1)，并通过插值方式得到中轴线图象及其各投影。

最后对模型给出检验方式。

一 、问题的重述假设某些血管可视为一类特殊的管道,该管道的表面是由球心沿着某一曲线 (称为中轴线)的球（命名为包络球）滚动包络而成。

现有某管道的相继100张平行切片图象,记录了管道与切片的交。

假设：管道中轴线与每张图片有且只有一个交点；球半径固定；切片间距以及图象象素的尺寸均为1.取坐标的Z 轴垂直于切片，第1张切片为平面0=Z ，第100张切片为平面99=Z . 计算管道的中轴线与半径，给出具体的算法，并绘制中轴线在XY 、YZ 、ZX 平面的投影图。

二、模型假设与符号说明1、 基本假设：（1） 该管道的表面为一定长半径的球沿一固定的曲线运动所得曲面族包络的光滑表面。

（2） 该管道的中轴线连续而且光滑。

（3） 该管道的中轴线与每个切面有且只有一个交点。

（4） 图象象素的尺寸为1. （5） 切片的间距尺寸为1.2、 符号说明：L 中轴线R 包络球的半径()z y x O i ,, 中轴线与第i 个切片的交点（定为此切片的中心）i S 第i 个切片切得的图形 i D 第i 个切片的图象数据矩阵三、问题分析及建模准备 问题分析：通常血管的表面可认为是连续且光滑的曲面，断面可用于了解其形态等特性。

本问题给出的是一些离散的切面，要求重建出原图中轴线和求出包络球半径。

因为每一个切面与中轴线L 有且只有一个交点i O ，如果找出所有i O ,就可以用插值或拟合的方式作出L 的近似图象，其在坐标平面上的投影就很容易画出。

问题的关健转变为求每个平面上的i O . 建模准备：1、 图象的读取由于切片图象中只有黑、白两种颜色的象素,而且所给的BMP 格式图象文 件是512×512象素的.因此,把图象读取为一个512×512的数字矩阵；用数字1表示黑色的象素，用数字0表示白色的象素。

血管的三维重建1 摘要序列图像的三维重建在各学科中都起到至关重要的作用，本次讨论的是血管的三维重建。

首先，假设该管道是由球心沿着某一曲面的球滚动包络而成，故本次的主要目的是求出中轴线坐标及半径。

现有100张平行切片图像，本次建立的模型可分为四步;第一步，采集图形边界点数据。

由于每张图片都是512*512的矩阵，故此数据很大，采用imread()函数将其读入矩阵A中。

第二步，最大内切圆寻找及半径的确定。

提出两种方案,分别是切线法和最大覆盖法; 从上述两种方法分析及考虑到我们所使用的工具和材料, 可以得出方法二更加直观, 计算机实现更容易, 计算复杂度更低, 所以我们采用后者。

根据以上算法, 我们抽取了所有的切片图进行半径的提取, 然后再求其平均值, 求其均值得到球的半径为29.6345。

第三步，轨迹的搜索。

在第二步中求出了血管的半径，轨迹的搜索就可以建立在半径确定的基础上, 当然我们也可以求出每一个切面图形的最大内切圆, 然后得到每个圆心的坐标, 即中轴线坐标, 但这样做计算机的运算量会很大, 同时由于最大内切圆搜索法的稳定性不高, 从而会造成搜索的不精确, 所以采用定半径搜索。

本文提出了三种方法, 分别为网格法、蒙特卡罗法和非线性规划法;本次采用非线性规划来实现。

第四步，绘制中轴线空间曲线图和在XOY、YOZ、XOZ 三个平面的投影图。

由定理1: 切片上血管截面图的头部顶点在XOY 平面上的投影点一定会落在中轴线在XOY 平面上的投影曲线上（在论文中以证明），并得出推论：切片上血管截面中中位线与中轴线在XOY 面上的投影重合。

最后可由中轴线和血管半径在作图软件中达到血管的三维重建，本次的模型还存在一定的不足，其假设为管道中轴线与每个切面有且只有一个交点，事实上还存在有多个交点的情况，但为了简化模型在此做了一定的假设，故会存在一定的误差。

关键词：三维重建内切圆半径轨迹（中轴线）注：求边界时采用了老师的思想和程序。

血管的三维重建数学建模
首先，血管的三维重建通常是通过医学影像学来实现的。

医学
影像学包括CT、MRI等技术，这些技术可以提供血管的断层扫描图像。

在这些图像的基础上，可以利用图像处理的方法，如边缘检测、分割等技术，来提取血管的形状和结构信息。

其次，几何建模是血管三维重建的关键环节。

在图像处理的基
础上，需要进行几何建模，将提取到的血管形状转化为数学模型。

这涉及到曲面重建、体素网格生成等技术，以及对血管内部结构的
建模。

另外，数学算法在血管三维重建中也起着重要作用。

例如，曲
面重建可以利用曲面拟合算法，体素网格生成可以利用体细胞自动
机等算法。

此外，对血管的分支、扭曲等特征的识别和建模也需要
借助数学算法来实现。

除此之外，血管的三维重建数学建模还涉及到计算机图形学、
计算几何学等领域的知识。

这些知识和技术的综合运用，可以实现
对血管形状、结构和特征的全面建模和重建。

总的来说，血管的三维重建数学建模是一个复杂而多样化的过程，涉及到多个学科和领域的知识。

通过综合运用图像处理、几何建模、数学算法等技术，可以实现对血管的全面、准确的三维重建和建模。

血管的三维重建摘要本题要求我们通过100张血管的平行切片图像，运用matlab进行数字图像处理,计算血管的中轴线和半径，绘制中轴线在XY、YZ、ZX平面的投影图，以达到将血管的三位形态表现出来的目的。

血管可视为一种特殊的管道，而管道的表面可视为由球心沿着某一曲线(中轴线)的球滚动包络而成。

血管中轴线与每张切片有且只有一个交点，且等径。

首先，我们确定每个切片中截面的轮廓线所包含区域内最大内切圆的圆心即为切片和血管中轴线的交点，半径即为血管半径；然后，我们用imread函数读取100张血管的切面图,位图文件就转换成了二值矩阵（像素点取值只限于0，1）,用函数edge可以提取出边界矩阵，用函数bwmorph可以提取出骨架矩阵。

我们抽取其中部分平行切片图，matlab编程计算，得出20组截面的最大内切圆半径和坐标（表1）的数据，以这些数据为基础，运用高次多项式拟合，绘制出中轴线在XY、YZ、ZX平面的投影图，对投影图进行3D旋转，即可以看到血管中轴线的空间大致形态。

关键词：MATLAB数字图像处理最大内切圆三维重建 matlab一 问题重述断面可用于了解生物组织、器官等的形态。

例如，将样本染色后切成厚约如m 的切片，在显微镜下观察该横断面的组织形态结构。

如果用切片机连续不断地将样本切成数十、成百的平行切片，可依次逐片观察。

根据拍照并采样得到的平行切片数字图象，运用计算机可重建组织、器官等准确的三维形态。

假设某些血管可视为一类特殊的管道，该管道的表面是由球心沿着某一曲线(称为中轴线)的球滚动包络而成。

例如圆柱就是这样一种管道，其中轴线为直线，由半径固定的球滚动包络形成。

现有某管道的相继100张平行切片图象，记录了管道与切片的交。

图象文件名依次为0.bmp 、1.bmp 、…、99.bmp ，格式均为BMP ，宽、高均为512个象素(pixel)。

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

基于模型的血管重建及应用丘文峰①鲍苏苏②冯天亮①吴应江①潘家辉②①广东医学院信息工程学院，523808，广东湛江霞山文明东路2号②华南师范大学计算机学院，520631，广州市天河区石牌摘要血管可视化对于临床诊断、手术计划起着至关重要的作用。

本文通过分析现有的血管可视化技术，选择基于模型的血管重建技术，利用vmtk工具包设计了一个针对血管树的中心线提取方案。

该方案由用户交互选择计算每一分支中心线并最终合并得到完整的中心线。

通过对64排CT腹部数据中的肝脏血管可视化，表明该方案可以得到较好的血管模型中心线。

关键词血管重建中心线骨架化图像分割1 引言血管重建和分析对于临床应用是十分重要的，最典型的例子是用于心血管疾病和涉及血管的外科手术方案制定，例如肝脏外科手术。

心脑血管病普通认为在西方国家和日本发病率高，然而2010年9月上旬，国际研究人员在北京提出最新证据，通过对50万亚洲人口的研究表明：心脑血管疾病正开始对中国及其他亚洲国家的健康和财富构成严重影响，且有年青化趋势。

血管动力学研究对于理解心脑血管疾病病理机制，并由此制定治疗方案十分重要。

血管动力学研究需要采用计算流体力学（CFD）方法进行模拟仿真。

而这依赖于通过影像设备，如CT、MR 采集的数据，以血管提取技术和三维可视化技术得到精确的血管几何模型。

在与血管相关的手术规划中，血管直径信息往往与规划相关。

血管三维可视化大致可以分成如下几步：通过图像分割，提取血管（部分可视化技术不需要）；对血管数据进行三维可视化；根据具体临床问题，采用建模分析技术，对血管结构进行定量分析。

由于部分可视化技术，可以直接对原始体数据操作，为此进一步划分两种不同的血管三维可视化策略（见图1）。

血管提取阈值法区域生长法水平集法等血管可视化MC算法曲面重建算法基于模型的算法血管定量分析伪彩血管直径度量血管与周边组织，病灶的关系分割结果（二值体数据）血管三维模型(a) 含血管提取的血管可视化策略血管可视化直接体绘制最大密度投影血管交互浏览体数据（b) 不含血管提取的血管可视化策略图1 两种不同的血管可视化策略图像分割作为图像分析图像理解的前期步骤，是计算机图像处理领域中最基本、也是最困难的问题之一。

α文章编号:100127445(2003)0320206203血管三维图像重建的数学方法黄新民(广西大学数学与信息科学学院,广西南宁530004)摘要:使用M A TLAB 读出血管切片图的数据,算出其半径及中轴线方程,建立了血管曲面的参数方程,实现了图形的三维重建.关键词:血管;参数方程;三维重建中图分类号:O 24211 文献标识码:A2001年全国大学生数学建模竞赛中的A 题:血管的三维重建,要求学生从给定的100张平行切片图中,编程读出图像中的数据,并以此计算管道的中轴线与半径,给出具体算法,并绘制中轴线在X Y ,Y Z ,ZX 平面的投影图.该题的具体计算已经在文[1～6]中有许多详细而细致的讨论.但是所有这些文章都没有说明使用的计算机语言及相应的程序.按题目所附参考材料看,出题人建议的是使用C ++程序.事实上许多软件都可以处理图形图像,而且还比使用C ++编程更容易.本文使用M A TLAB 比较简单地完成题目所要求的计算.原竞赛题目中只要求给出图像的中轴线及其在三个坐标平面上的投影,实际上并没有真正实现血管的三维图像重建,本文将讨论实现图形三维重建的一个方法.要将图形数据读入,在M A TLAB 中只要使用命令i m read 即可实现.程序虽然简单,但数据量较大,因此要通过几次计算才能将所有数据读完,再用命令save 存储到硬盘上,以便下一步计算时使用.对于读出的每一张图,我们需要读出血管截面数据及其边界数据,以便下一步从中读出每张内切圆的圆心坐标与内切圆半径.M A TLAB 记录的图形数据是这样的:如果记录图形数据的是一个二维数组A (i ,j ),则第i 行第j 列上的数值为一个0到255的整数(占用一个字节的存储空间),这个值对应于一个有三列元素的长为256的颜色向量表中的一行,该行的三分量就分别对应于红、绿、兰颜色的颜色值.本题中由于只有黑白两色,所以数组中每个元素只有值为零及为1两个.在程序中我们用值为1记血管截面上的点,因此就用周围四个点的值和为4(因此每一个点的值都是1)来找出血管截面的内点,去掉内点后我们就得到了血管截面的边界点.不过我们得到的还不能算是真正的血管截面的边界.因为图形扫描成数据点时,血管截面中的每一个点实际代表了一小片面积的中心.我们不妨称这个边界为血管截面的内边界.而元素值为0的点(也就是血管截面外的点)的边界称为外边界.将内外边界求平均表示真正的边界才更准确.读出边界及所围的区域数据后,我们从区域中每取出一个点,计算其与边界的最小距离(这就是在这个点所作的与边界相切的圆半径),再从所有内切圆的半径中找出最大时的圆心坐标与半径,就得到小球在这个血管截面处时的球心坐标与球半径.求出所有各张图的球心坐标与球半径,就得到了中轴线坐标及球半径(严格地说球半径应该是一个常数,但是由于扫描与计算误差,因此只能取所有半径的平均值为球心半径).此程序的编写也不复杂.但实际使用时不难发现由于数据量大,需要计算很长时间(当然我们将程序运行后就可以静候结果),为加快计算速度可以使用各种办法.在本文中我们发现圆半径不小于29,第28卷第3期2003年9月广西大学学报(自然科学版)Journal of Guangxi U niversity (N at Sci Ed )V o l .28,N o.3　Sep t .,2003　α收稿日期:20021208;修订日期:20030810作者简介:黄新民(1945),男,广西贺州人,广西大学教授.因此将内点与边界点相同x 坐标的边界点比较其y 坐标,当两个点的坐标之差小于29时就将这个点从内点中排除.对坐标y 相同的点亦作同样处理.由于减法比计算距离的平方运算要快,因此这样做可以加快计算过程.算出中轴线坐标及球半径后,只要使用简单的二维作图命令p lo t 就可以画出中轴线向三个坐标平面的投影图,完成建模题的要求.但是这样做出的图形并没有实现图形的真正三维重建.血管的三维图像重建有多种方法,最简单的是将每张图的轮廓线按其高度在三维空间中复原.本文将介绍实现图像三维重建的另外一种方法.这种方法与学生在大学中所学的高等数学及线性代数的知识密切相关,而且重建的图像可以使用光照、颜色等加工,使图形显得十分美丽,值得介绍.设已求出中轴线方程:x =x (z ),y =y (z ),0≤z ≤99,下面我们求半径为r 的球沿此中轴线运动而得的包络面方程.设(X ,Y ,Z )是包络面上任意一点,则存在z 使得(X -x (z ))2+(Y -y (z ))2+(Z -z )2=r 2.(1)上式的z 为一个参数,当它从0运动到99时其轨迹就形成了管道.将(1)式两边对参数z 求导,则得第二个方程x ′(z )(X -x (z ))+y ′(z )(Y -y (z ))+(Z -z )=0.(2)将(2)式代入(1),化简后得(1+x ′(z )2)(X -x (z ))2+2x ′(z )y ′(z )(X -x (z ))(Y -y (z ))+(1+y ′(z )2)(Y -y (z ))2=r 2,(3)(3)式左方是一个二次型,容易求出其特征值分别是Κ1=1,Κ2=1+x ′(z )2+y ′(z )2,其相应的特征向量是p 1=-y ′(z )x ′(z ), p 2=x ′(z )y ′(z ).因此作变换X -x (z )=-y ′(z )u +x ′(z )v , Y -y (z )=x ′(z )u +y ′(z )v ,(4)之后,(3)式将化简为(x ′(z )2+y ′(z )2)(u 2+(1+x ′(z )2+y ′(z )2)v 2)=r 2,(5)因此,(u ,v )可建立参数式u =r co s (t )x ′(z )2+y ′(z )2, v =r sin (t )(x ′(z )2+y ′(z )2)(1+x ′(z )2+y ′(z )2).将上式代入(4)式,立即得到管道曲面的参数方程X =x (z )+-ry ′(z )co s (t )x ′(z )2+y ′(z )2+rx ′(z )sin (t )(x ′(z )2+y ′(z )2)(1+x ′(z )2+y ′(z )2),Y =y (z )+rx ′(z )co s (t )x ′(z )2+y ′(z )2+ry ′(z )sin (t )(x ′(z )2+y ′(z )2)(1+x ′(z )2+y ′(z )2),Z =z -x ′(z )2+y ′(z )21+x ′(z )2+y ′(z )2r sin (t ).其中0≤z ≤99,0≤t ≤2Π,在中轴线坐标及球半径r 求出来后,可以通过数据拟合得到中轴线的参数方程(由于数据是近似的,因此作数据拟合是合适的)x =x (z ), y =y (z ),　0≤z ≤99.然后就可以画出血管的三维图形了.通过选择观察点及颜色、光照等可以作出十分精致的三维图像.下面的程序即为作图程序(其中x (z ),y (z )是分别拟合成六次多项式).图1是此程序运行后得到的M A TLAB 作出的血管三维图.load cendata %调入中轴线坐标数据,假设存储在变量B 中,共100行,分别表示x ,y ,r 坐标x =B (:,1);%中轴线的x 坐标向量y =B (:,2);%中轴线的y 坐标向量r =sum (B (:,3)) 100;%内切圆半径平均值702第3期黄新民:血管三维图像重建的数学方法z =[0:99]’;%z 坐标向量p 1=po lyfit (z ,x ,6);%以z 为自变量将x 拟合成六次多项式p 2=po lyfit (z ,y ,6);%以z 为自变量将y 拟合成六次多项式zz =[0:0.4:99]’;%重新定义步长更小的高度坐标向量zz xx =po lyval (p 1,zz );%用拟合多项求对应于zz 的坐标xx yy =po lyval (p 2,zz );%用拟合多项求对应于zz 的坐标yy dx 1=po lyder (p 1);%对x 的多项式微分多项式dy 1=po lyder (p 2);%对y 的多项式微分多项式dx 2=po lydval (dx 1,zz );%求x’(z )dy 2=po lydval (dy 1,zz );%求y’(z)图1　用M A TLAB 作出 的血管三维图R 2=dx 2,^2+dy 2.^2;r 1=sqrt (R 2);r 2=r 1.3sqrt (1+R 2);t =linspace (0,23p i ,251);%定义角度s =size (t );X =xx 3ones (s )+r 3((-dy 2. r 1)3co s (t )+(dx 2. r 2)3sin (t ));Y =yy 3ones (s )+r 3((dx 2. r 1)3co s (t )+(dy 2. r 2)3sin (t ));Z =zz 3ones (s )-r 3(r 1. r 2)3sin (t );m esh (X ,Y ,Z ) %三维网线作图surface (X ,Y ,Z )%三维表面作图axis equal%按实际比例值作图axis off%不画坐标轴,使图形更逼真h idden off%消隐,去掉曲面上的网线,使图形更逼真view ([35,88])%选择观察角shading interp%使作图时片与片之平滑过渡co lo r m ap ho t%色彩处理(不一定恰当,请读者自己尝试)caxis ([-70,450])%重设颜色,使颜色接近红色cam ligh t (200,180)%指定光源位置,这些值为尝试值,不一定合适ligh ting gouraud %设置照明方式参考文献:[1]　廖武斌,邓俊晔,王　丹.管道切片的三维重建[J ].工程数学学报,2002,19(5):22228..[2]　胡亦斌,向　杰,程　翔.利用切片的二维空间相关操作实现血管的三维重建[J ].工程数学学报,2002,19(5):29234.[3]　徐　晋,刘雪峰,柏容刚.血管的三维重建[J ].工程数学学报,2002,19(5):35240.[4]　柳海东,陈　璐,江　浩.管道切片的三维重建[J ].工程数学学报,2002,19(5):41246.[5]　丁峰平,周立丰,李孝朋.血管管道的三维重建[J ].工程数学学报,2002,19(5):47253.[6]　汪国昭,陈凌钧.血管三维重建的问题[J ].工程数学学报,2002,19(5):54258.3D recon struction of blood vesselHU AN G X in 2m in(Co llege of M athem atics and Info r m ati on Science ,Guangxi U niversity ,N anning 530004,Ch ina )Abstract :M A TLAB is u sed to ob tain data of b lood vessel slices .T he radiu s of the b lood vessel is deter m ined and the axes cu rve has calcu lated .T he p aram etric equati on is ob tained and the 3D i m age of b lood vessel su rface is recon structed .Key words :b lood vessel ;param etric equati on ;3D recon structi on(责任编辑　刘海涛)802广西大学学报(自然科学版)第28卷　。

血管的三维重建模型摘要：本文对血管三维重建中，中轴线及球的半径确定问题进行了讨论。

首先，根据问题及图象处理提取有效数据，给出两种可行算法，利用上述数据建立了最大最小方法和二次规划方法。

搜索中心点，并给出全局和局部搜索，得到各切片中心点坐标(见表1)，并通过插值方式得到中轴线图象及其各投影。

最后对模型给出检验方式。

一 、问题的重述假设某些血管可视为一类特殊的管道,该管道的表面是由球心沿着某一曲线 (称为中轴线)的球（命名为包络球）滚动包络而成。

现有某管道的相继100张平行切片图象,记录了管道与切片的交。

假设：管道中轴线与每张图片有且只有一个交点；球半径固定；切片间距以及图象象素的尺寸均为1.取坐标的Z 轴垂直于切片，第1张切片为平面0=Z ，第100张切片为平面99=Z . 计算管道的中轴线与半径，给出具体的算法，并绘制中轴线在XY 、YZ 、ZX 平面的投影图。

二、模型假设与符号说明1、 基本假设：（1） 该管道的表面为一定长半径的球沿一固定的曲线运动所得曲面族包络的光滑表面。

（2） 该管道的中轴线连续而且光滑。

（3） 该管道的中轴线与每个切面有且只有一个交点。

（4） 图象象素的尺寸为1. （5） 切片的间距尺寸为1.2、 符号说明：L 中轴线R 包络球的半径()z y x O i ,, 中轴线与第i 个切片的交点（定为此切片的中心）i S 第i 个切片切得的图形 i D 第i 个切片的图象数据矩阵三、问题分析及建模准备 问题分析：通常血管的表面可认为是连续且光滑的曲面，断面可用于了解其形态等特性。

本问题给出的是一些离散的切面，要求重建出原图中轴线和求出包络球半径。

因为每一个切面与中轴线L 有且只有一个交点i O ，如果找出所有i O ,就可以用插值或拟合的方式作出L 的近似图象，其在坐标平面上的投影就很容易画出。

问题的关健转变为求每个平面上的i O . 建模准备：1、 图象的读取由于切片图象中只有黑、白两种颜色的象素,而且所给的BMP 格式图象文 件是512×512象素的.因此,把图象读取为一个512×512的数字矩阵；用数字1表示黑色的象素，用数字0表示白色的象素。

血管的三维重建摘要本文以血管的三维重建为研究对象，对100张平行切片图像进行分析，利用这些宽、高均为512象素的切片，计算管道的半径和确定中轴线方程，并在此基础上画出重建后的血管三维图像，主要内容如下：对于问题一，计算管道的半径，由于血管表面是由球心沿着某一曲线（称为中轴线）的球滚动包络而成，可以得出结论：切片中包含的最大圆的半径即血管半径，所以问题转化为求每一切片上的最大内切圆的半径。

为了便于计算，运用Matlab imread 函数，将BMP 格式文件转化为0-1矩阵，然后运用edge bwmorph 、函数确定轮廓和骨架的位置，并求解骨架上每一点到边缘的最短距离。

这些最短距离中的最大值即为最大内切圆半径也就是血管半径。

最后对所有的半径取平均值，得出结果：100()1=29.41666100k k RR ==∑对于问题二，根据问题一中求出的100个圆心坐标及半径求解中轴线方程，运用Matlab 软件对圆心所形成的曲线进行n 阶多项式拟合。

为使中轴线较为光滑，在Matlab 拟合工具箱多次试验后，取最高阶次=7n 。

由于z 轴值是逐层单调递增的,为简化方程的计算,取t 为参变量，分别对其投影在YZ 、ZX 平面上进行多项式拟合，最后得到中轴线在平面投影上拟合的曲线方程如下：()()()-107-76-55432-107-86-55-3432-3.2310 1.16910-1.628100.00108-0.035260.5706-3.105+5.243=3.06110-9.62310+1.3610-0.640610+0.01912-0.298+1.89-1.63.3=y t t t t t t t t f x t t t t t t t t z t t ⎧=⨯+⨯⨯+⎪+⎪⎪=⨯⨯⨯⨯⎨⎪⎪⎪⎩最后根据方程画出中轴线图形，YZ YX ZX 、、平面的投影在拟合工具箱中可以直接得到。

对于问题三，根据问题一、二求出的中轴线的参数方程和100张切片的最大内切圆的半径，运用Matlab 软件画出血管的三维立体图。

基于MATLAB的血管三维重建
刘波
【期刊名称】《湖北民族学院学报（自然科学版）》
【年(卷),期】2008(026)004
【摘要】根据血管序列切片的二值图像特征,利用MATLAB丰富的矩阵运算和图像处理命令,将血管三维重建过程分成半径搜索、交点定位和轴线拟合三个主要步骤,并编制通用M程序包实现从数据采集到模型渲染的全程自动计算.最后应用给出的方法完成了对100张序列切片图像的计算机三维重建.
【总页数】3页(P398-400)
【作者】刘波
【作者单位】湖北民族学院数学系,湖北恩施445000
【正文语种】中文
【中图分类】O242.1;TP391.41;R322.1+2
【相关文献】
1.基于MATLAB的CT图像三维重建的研究与实现 [J], 张振东
2.一种基于MATLAB的快速MRI图像三维重建方法 [J], 王军;王昕
3.基于MATLAB的脑CT图像三维重建研究 [J], 裴智军;高关心;谢伟
4.基于MATLAB的粗粒土组构的三维重建 [J], 袁强; 汪向兰; 石文杰; 杨曼; 陈琰
5.基于MATLAB的运动恢复结构三维重建算法的研究与实现 [J], 李美燕;黄世玲因版权原因，仅展示原文概要，查看原文内容请购买。

第19卷　建模专辑2002年02月工　程　数　学　学　报J OU RNAL OF EN GIN EERIN G MA THEMA TICSVol.19Supp.Feb.2002文章编号:100523085(2002)0520035206血管的三维重建徐　晋(电子工程与信息科学系),　刘雪峰(数学系),　柏容刚(电子工程与信息科学系)指导老师:窦　斗(数学系)(中国科学技术大学,合肥230026)编者按:本文分析了每张切片与管道曲面的交线是一族圆的包络线,且这族圆中半径最大者即为最大内切圆这一几何特性,建立了相应的算法,将中轴线拟合成Bézier曲线,并用两种算法对模型进行了检验,方法有一定特色。

摘　要:对血管的三维重建问题,我们假定血管为等径管道,通过分析其几何特性,给出了确定其管道中轴线和半径的数学模型———搜索每个切片截面,求最大内切圆,该内切圆圆心即为切片截面与管道中轴线的交点,该内切圆半径即为管道半径,再通过拟合各个交点求出轴心线。

本模型中,我们确立了两种有效的误差分析方法;并由此发现由于中轴线与切片交角过小会使结果产生较大偏差。

为解决此问题,我们从其它方向重新对血管进行切割,再进行处理求解,得到更加精确的结果。

关键词:血管;等径管道;旋转切面分类号:AMS(2000)65D17 中图分类号:O242.1 文献标识码:A1　问题重述(略)条件假设1).血管的表面是由半径固定、圆心连续变化的一族球滚动形成的包络面。

2).医学上,血管不存在严重扭曲。

3).管道中轴线与每张切片有且只有一个交点。

2　问题分析与模型建立根据假设,血管可视为表面是由球心沿着某一曲线(称为中轴线)的球滚动包络而成的管道。

根据所查文献(参考文献[6]),可知这种管道有如下几何特性:定理　等径管道每个切片的轮廓线是一族半径、圆心连续变化的圆的包络线,而这族圆中半径最大的圆的圆心即为管道的中轴线与切片的交点,半径即为管道半径。