《【数学与应用数学】论文——血管三维重建的数学模型》由会员分享，可在线阅读，更多相关《【数学与应用数学】论文——血管三维重建的数学模型（9页珍藏版）》请在金锄头文库上搜索。
1、27 血管三维重建的数学模型血管三维重建的数学模型 摘要摘要:本文探讨了血管的三维重建问题.我们首先利用数字图象处理的相关技术,从 100 张图中读取出相应数据,然后以定位所需的相关参数(包括中心轴线和半径)为目标,主要运用了包络面相关理论及几何拓扑网络中的相关知识得出了如下结论：切片一定包含一个滚动球的大圆,并且它恰好是切片内的最大圆.运用这些结论设计出有效的算法,利用数学软件(Matlab 及 Maple)和计算机技术,对 100 个切片编程求解,找出切片所包含的最大圆,即得到中心轴线与 100 个切片交点的坐标及大圆半径.经过误差分析确定出管道半径为,并给出了中心轴线上点的坐标,最后我们
2、还绘制出中心轴线在平面上的69849.29rzxyzxy,投影图.关键词关键词：三维重建;最大圆;中轴线;有障碍物的最大空隙问题1 1 问题的提出问题的提出生物组织连续切片的计算机三维重建技术是把一系列切片的图象，通过计算机进行处理,从而得到该生物组织立体结构的一种方法.在切片图象内部寻找定位结构,是此核技术的关键,同时也是目前国内外三维重建技术的共同难题.解决此难题的一个重要方向是通过提取切片中二维图象的数据信息,建立模型,以达到定位效果.假设某类血管表面可视为由球心沿着一曲线（称为中轴线）的球滚动包络而成.现有一此类血管相继 100 张平行切片图象,记录了管道与切片的交.图象格式为 BMP
3、,宽、高均为512 个象素（pixel）.假设管道中轴线与每张切片有且只有一个交点,球半径固定,切片间距以及图象象素的尺寸均为 1.如果取坐标系的 Z 轴垂直于切片,第 1 张切片为平面,0Z第 100 张切片为平面.切片图象中象素在文件中出现的前后次序为99ZzZ （-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 平面的投影
4、图.2 2 基本假设和符号说明基本假设和符号说明28 1) 假设血管为实心体,其表面为一包络面；2) 假设管道中轴线与每张切片有且只有一个交点,滚动球的半径固定；3) 假设切片间距以及图象的象素尺寸均为 1；4) 表示中轴线,是一空间曲线；m5) 表示滚动球的半径；r6) 表示球心沿中轴线滚动形成的包络面；S7) 表示由曲面与平面,所围成的管道体；VS0Z99Z8) 表示管道体与平面的切片；iFViZ 9) 表示的边界.ieiF没有说明的符号在文中初次引用时均作出了说明.3 3问题的分析与模型的建立问题的分析与模型的建立首先分析题目中所给出的 100 张图片.对这 100 张黑白图片通过计算机
5、图象处理,我们可以观察到其切片呈现一个旋转渐变的过程,如图（一）因此我们估计该血管的形状类似一条螺旋管.只要将这一百张切片在空间中定位,由于切片间距非常小（为一个象素）,所以就可以近似地还原出血管的三维结构.现先对每张图作出一些数据上的分析.针对 BMP 图象的存储格式,我们将所给的图片转为二值图象处理.所谓二值图象,就是只有黑白两个灰度级,即象素灰度级非 1 即 0.我们可以从每张图中读取出一个 512 阶的 0-1矩阵.其中 0 表示象素为黑色.1 表示象素为白色.我们可以给出矩阵元素的位置与题目中所定坐标系之间的坐标变换.即矩阵的第 行,第列的元素坐标相应转化为.ij)256,257(i
6、j29 现在我们把问题转化为数学语言.在笛卡尔坐标系内,有一半径为的球,沿空间的某一r曲线（称为中轴线）,滚动包络形成空间曲面,现有 100 张平mS面图,其中,与曲面围成一实心管道体,现题目就是要求99, 1, 0ZZZ0Z99ZS我们求出该管道体的中轴线及滚动球的半径,并要给出中轴线在 XY,YZ,ZX 平面上的内mm投影图.现设管道中轴线的方程为：Vm其中 )()()(tzztyytxx99)(, 0)(,990990tZtZttt由于管道体中心轴线与每张切片有且只有一个交点,则不妨设中心轴与切片mm有且只有点,与的间距为 1 个象素,在如此小的间距)990(iFi),(),(itYtX
7、OiiZ1iZ下,我们可以用 100 个交点近似地刻画中轴线.iO结论 1：管道体与平面的切片一定包含一个过滚动球球心,半径为的圆.ViZ iFrA现证明如下：以中心轴与平面的唯一交点为球心（亦即包含一个滚动球的大miF),(),(iiYtXOi圆）,为半径所作出的球与平面交于圆,其中圆的半径为,由于球在管道rAiZ AArA体内,所以球与平面的切片也在管道体与的切片内,亦即圆在内,结论 1VAiZ iZ iF得证.结论 2：管道体与平面的切片内的最大圆恰好是滚动球中过球心半径为ViZ iF的圆,并且最大圆唯一.rA现证明如下：由中轴线的方程：m)()()(tZztYytXx,990ttt )
8、 1 (其中.99)(, 0)(990tZtZ有曲面（即包络面）的方程来：S 0)( )()( )()( )()()()(2222tZtZztYtYytXtXxrtZztYytXx)2(30 现考虑,取使得则iZ tt itZ) ( 0)( )()( )()()(0000222tYtYytXtXxrtYytXx)3(（3）式的几何意义是,在切片的边界曲线上的点与中轴线和切片的交点的最近iFiS距离为,即可作以为圆心,半径为的圆.r)(),(tYtXr边界曲线和中轴线与切片的交点的距离达到最大值.iSmiF),(),(itYtXOir又由于每一个 属于集合,与中轴线只有唯一交点,故只有i99,
9、1 , 0iz )()()(tZztYytXx唯一点 满足,所以切片包含的最大圆唯一.titZ)(结论 2 得证.至此,我们只要在每个切片上找出一个能被切片包含的最大圆的圆心,即为该切片与iF中轴线的交点.就一般包络面而言,切片的最大圆应是唯一的,但由于数据存储误差,我们可能会找到多于一个的最大圆.所以我们需要一个能找出多人最大圆的可行算法.利用以上的分析和结论,我们试图求出管道体半径及每个切片与中轴线的 100 个交点,由r于 BMP 存储方式的限制,每个切片实际是由有限多个点组成.我们利用几何拓扑网络中的有障碍物的最大空隙问题解法,把障碍物视为半径为 0 的圆,因此我们可以在切片上作如下约
10、束优化：222)()(R s.t.R maxiiyyxx其中,iiiiFyxSyx),( ),(为管道体切片,为的边界,且、均为限点集.iFiSiFiFiS上述模型可简化为如下：)()(22 iiyyxxMinMaxr其中 为管道体切片,为的边界,且,均为iiiiFyxSyx),( ),(iFiSiFiFiS限点集,关于的求解和圆心坐标,算法及求解结果将于下部分给出.r4 模型的求解与算法的实现模型的求解与算法的实现31 我们根据上面的在切片上进行约束优化思想,利用计算机编程来寻找最优解.由于切片图象是用 BMP 格式存储的,因此它都有是由一些有限的点构成,我们用 Matlab读入每张图,用一
11、个 512512 的 0-1 矩阵表示,其中 0 表示象素为黑色的点.去掉所有非零元素,并记下所有零元素在矩阵中出现的位置,再转化为题中的笛卡儿坐标系的坐标,这样,就可以用描点方法作出切片的用点集表示的图,如果我们对 100 个切片上的点都在空间中描出来,那么就可以大体反映出管道体的三维结构.受图片存储的精度所限,只能在表示切片的点集上寻找最大圆心.我们根据上节中所列出来的约束优化思想,再利用计算机编程的方法来寻找最大圆.首先我们要作出切片的边界.由于图上的点都是呈阵列排布的,只要依次取切片上的点,选出与 1 相邻的 0 点即为切片的边界.现在按如下算法寻找最大圆的圆心和半径:1 取出切片点集
12、中的一个点,再算出该点与边界上各点的距离,找出最小值,记为 h,这就是以该点为圆心所能作出的包含于切片的最大圆的半径.2 继续取点集中的下一个点,再算出该点与边界上各点的距离,再找出最小值与 h 进行比较,若比 h 大则把值赋给 h,若与 h 相等则把值赋给 h 并记下坐标,同时还把这个点都记起来,这样就保证若最大圆不唯一,也可把所有最大圆找出来.3 取遍切片点集上的所有点运行第二步,最后得出的 h 即为该切片中最大圆的半径,相应该点的坐标即为切片与中轴线交点的坐标. 源程序见附录我们记每张切片上点的数目为 n1,切片边界上点的数目为 n2,显然 n2 远小于 n1,我们要进行 n1 重循环,
13、每一次循环作 n2 次比较,其复杂性不超过 O,所以,此算法为有 21n效算法.我们用 Matlab 编程进行搜索,利用 Matlab 强大的矩阵运算功能,在编程时进行向量化,大大提高了程序的效率,对 100 张图分别算出如下半径和交点的坐标：32 注：上表中（X,Y,Z）表示与中轴线交点的笛卡儿坐标.R 表示在切片上最大圆半径.由于图象精度问题及算法的误差.每张切片上算出的半径有微小差异,我们根据误差产生的原因进行分析（详细分析过程下节给出）,取出所有半径中最大值 r=29.69849,即为题目所XYZ R0160029.000000160129.000000160229.000000160
14、329.000000160429.000000160529.000001160628.861741160729.000001160829.000001160929.0000021601029.0000021601129.0000031601229.0000041601329.0000051601429.0000061601529.0000081601629.00000XYZ R101601729.00000151601829.00000171601929.00000181602029.00000191602129.00000201602229.00000211602329.0000020160
15、2429.01724201602529.01724291592629.01724291592729.06888351582829.06888351582929.06888401573029.12044351583129.15476481553229.12044481553329.12044XYZ R511543429.12044571523529.12044601513629.12044601513729.12044651493829.15476651493929.15476651494029.15476881374129.15476911354229.15476941334329.15476911354429.41088911354529.41088911354629.4108811 41164729.6984911 61144829.6984911 71134929.6984911 81125029.69849XYZ R1181125129.698491181125229.698491181125329.698491191115429.41088132955529.20616132955629.20616146725729.41088147705829.41088152595929.15476152596029.15476161306129.12044162256229.12
