血管的三维重建.docx

1、血管的三维重建血管的三维重建1摘要序列图像的三维重建在各学科中都起到至关重要的作用，本次讨论的是血管 的三维重建。首先，假设该管道是由球心沿着某一曲面的球滚动包络而成，故本 次的主要目的是求岀中轴线坐标及半径。现有100平行切片图像，本次建立的模型可分为四步；第一步，采集图形边界 点数据。由于每图片都是512*512的矩阵，故此数据很大，采用imread（）函数将 其读入矩阵A中。第二步，最大切圆寻找及半径的确定。提出两种方案分别是切线法和最大 覆盖法；从上述两种方法分析及考虑到我们所使用的工具和材料可以得出方 法二更加直观计算机实现更容易计算复杂度更低.所以我们采用后者。 根据以上算法，我们

2、抽取了所有的切片图进行半径的提取.然后再求其平均值. 求其均值得到球的半径为29. 6345。第三步，轨迹的搜索。在第二步中求出了血管的半径，轨迹的搜索就可以建 立在半径确定的基础上.当然我们也可以求出每一个切面图形的最大切圆然 后得到每个圆心的坐标，即中轴线坐标，但这样做计算机的运算量会很大.同 时由于最大切圆搜索法的稳定性不髙.从而会造成搜索的不精确.所以采用定 半径搜索。本文提岀了三种方法.分别为网格法、蒙特卡罗法和非线性规划法； 本次采用非线性规划来实现。第四步，绘制中轴线空间曲线图和在XOY. YOZ. XOZ三个平面的投影图。 由定理1:切片上血管截面图的头部顶点在XOY平面上的投

3、影点一定会落在中 轴线在X0Y平面上的投影曲线上（在论文中以证明），并得出推论：切片上血管 截面中中位线与中轴线在XOY面上的投影重合。最后可由中轴线和血管半径在作图软件中达到血管的三维重建，本次的模 型还存在一定的不足，其假设为管道中轴线与每个切面有且只有一个交点，事实 上还存在有多个交点的情况，但为了简化模型在此做了一定的假设，故会存在一 定的误差。关键词：三维重建切圆半径轨迹（中轴线） 注：求边界时采用了老师的思想和程序。2问题重述假设某些血管可视为一类特殊的管道，该管道的表面是由球心沿着某一曲线 (称为中轴线)的球滚动包络而成。例如圆柱就是这样一种管道，其中轴线为直线, 由半径固定的球

4、滚动包络形成。现有某血管管道的相继100平行切片图象，记录了管道与切片的交点。图象 文件名依次为0.bmp、l.bmp、99. bmp,格式均为BMP,宽、高均为512个象 素(pixel) o先假设：管道中轴线与每切片有且只有一个交点；球半径固定；切 片间距及图象象素的尺寸均为lo取坐标系的Z轴垂直于切片，第1切片为平面 Z=0,第100切片为平面Z=99o Z=z切片图象中象素的坐标依它们在文件中出现 的前后次序为(256, 一256, z), (256, 一255, ,(一256, 255, z)(255, 一256, z), (255, 一255, ,(一255, 255, z)(25

5、5, 256, z), (255, 255, z), -(255, 255, z)。试计算管道的中轴线方程与半径，给出具体的算法，并绘制中轴线在XY、 YZ、ZX平面的投影图。全部图象请从网上下载。关于BMP图象格式可参考：1.VisualC+数字图象处理第12页2. 3. 1节。何斌等编著，人民邮电，2001 年4月。2.http： / / . des. ed. ac. uk / home / mxr / gfx / 2d / BMP.txt3模型的假设1医学上.血管不存在严重扭曲.没有折皱。2血管可视为等径管道。3.管道中轴线与每切片有且只有一个交点。4切片间距以及图象象素的尺寸均为lo5

6、.对切片拍照的过程中不存在误差，数据误差仅与切片数字图象的分辨率有关。6.血管的表面是由半径固定的球的球心沿着某一曲线(称为中轴线)滚动包 络而成。7 切片足够薄其厚度对计算的影响可以忽略不计。8.切片包含的其它圆的半径一定小于r0o9.不能被切片包含的圆的半径一定大于r。4符号说明XYZR血管的最大切圆的圆心的x的坐标值 血管的最大切圆的圆心的y的坐标值 血管的最大切圆的圆心的z的坐标值 血管最大切圆的的半径R.血管最大切圆整个圆周的离散坐标A图像数据储存的三维矩阵/（鸟）血管最大切圆心匹配点数的函数5分析与建立模型5. 1分析1）假定管道中轴线与每切片有且只有一个交点.球的半径固定，切片的

7、间 距以及图像象素的尺寸均为1。取坐标的z轴垂直于切片，第一切片记为平面z =0,第二切片记为z = 1,第100切片的平面记为z = 99O要求根据以上所给 条件计算管道的中轴线与半径.并绘制中轴线在X Y、YZ、ZX平面的投影图。5.2几个重要的结论对于某个切痕可以作出它在空间中的结构图。如图1所示设某一切片的切 痕与中轴线的交点为0,切下这一切痕的同时也截下中心在点0的小球的一个大 圆K过点0的血管的轴截面S也是同一个小球的一个大圆。这两个大圆的交点即 为小球的直径的两个端点且这两个端点为切痕边界与大圆K的两个切点因此 可得岀以下结论：1） 每一个切痕中存在一个最大切圆且最大切圆的圆心在

8、管道的中轴线上。2） 每个切痕的最大切圆的半径均相等.且等于管道的半径。根据上面的结论要求血管的中轴线只需找出每一个切痕的最大切圆的圆心。图1 某个切痕的空间结构图5.3模型的建立通过对血管的三维重现的分析可知，我们的模型应由三部分组成：5. 3.1采集图形边界点数据边界点数据的采集方法有人工采集法和计算机采集法.对于数据量较小的. 又具有某些特殊意义的点可以由人工采集.若数量比较大且采集边界又非常复 杂时，则采取机器采集。5. 3.2最大切圆寻找及半径的确定，左，左此方法可以从切面的外围轮廊线分析着手。由题意分析知道所给的图片切面是由无数个球切面组成的。而且外围轮廓线与最大切圆有且仅有两个交

9、点.所 以经过这两点的外围轮廓线的两条切线平行且间距最大。基于上述分析.我们可 以通过找到这两条切线来找到最大切圆的圆心及半径。在实际操作中.由于对图 片的象素提取的离散性.我们在计算导数时是用差分来代替。方法2最大覆盖法最大覆盖法就是在切面中找到最佳的圆心位置和半径长度.从而使得由这 个圆心和半径所决定的圆面能最大面积地覆盖管道切面的图形这样搜索到 的圆一定是最大切圆.这个圆的圆心就是我们所要找的球的球心这个圆的半 径就是我们所要找的球的半径。从上述两种方法分析及考虑到我们所使用的工具 和材料.可以得出方法二更加直观计算机实现更容易计算复杂度更低.所 以我们采用后者。具体实现中，我们先得到任

10、意一图片的象素矩阵然后将用于 匹配的圆根据其圆心和半径将其圆周离散（即以象素表示），并映射到 512x512的图中，其中圆周上的点为,其余的点为1.即形成另一个象素矩阵。这512x512两个矩阵在相同位置点上的值进行逻辑或运算.如果其值为0,则为匹配点.即此点在管道切面图形.否则其在切面图形外。这样搜索到的匹配 点最多同时半径最大的就是所要找的最大切圆。图2和图3就是分别从图1. bmp 和图89. bmp搜索到的最大匹配圆（部白色部位为管道切面图形深色圈为最 大匹配圆）。根据以上算法.我们抽取了所有的切片图进行半径的提取然后再 求其平均值.求其均值得到球的半径为29.6354。图2 1. b

11、mp的最大匹配图 图3 89. bmp的最大匹配图5.3.3轨迹的搜索5. 3. 3. 1目标函数的确立在求出半径以后，轨迹的搜索就可以建立在半径确定的基础上，当然我 们也可以求岀每一个切面图形的最大切圆.然后得到每个圆心的坐标.即中轴 线坐标但这样做计算机的运算量会很大，同时由于最大切圆搜索法的稳定性 不髙.从而会造成搜索的不精确，所以采用定半径搜索。我们通过定圆（半径为 R ）来找其中轴线.也就是用定圆覆盖到切面图形上去.找到匹配点数最多的 一个位置.从而得到此定圆圆心的位置具体实现时只要用定圆的圆心位置进行 变化，设其为Agy） 则由A点可以得到整个圆周的离散坐标坷，令/（倚） 为匹配点

12、数的函数，其计算方法与最大切圆求法相同，即将心 根据圆心A的 坐标和半径(固定)离散化形成矩阵后和切面图形矩阵作逻辑与操作，从达 到匹配最大，则优化目标函数就是：max/() Ae(512x512)坐标面上的点5. 3. 3. 2目标函数要求求出一点圆心A ,使得定圆覆盖切面图形最大，即多 元函数极值最优解，这是一个以X和Y为自变量二元函数，这样可以通过以下 二种肯泮非耒得_1)用直接搜索求最优解(网格法) 搜索过程中，对每一个圆心的坐标X和Y .在其取值的围均取100个步长.分 为1002个网格这样.在一定的精度围.可以求的一个较好的最优解。2)蒙特卡罗法蒙特卡罗法.也就是随机实验试点法。它

13、的基本思想是：在函数的可行域随机地 选取实验点.由于随机取得的点在区域中分配比较均匀所以对函数的大致形 态能较好的体现。模型中，随机点用以下方法产生的。x = xO + (xl-.vO) x rand()y = yO + (yl yO) x rand(其中，(xO, xl)为X的取值围，(yO, yl)为Y的取值围。3)非线性规划非线性规划即无约束优化，以数值迭代为基本思想基本步骤为选取初值 A(XO,YO).进行k次迭代并求岀迭代解由迭代解得到搜索方向和步长如果k+ 1 次迭代符合给定的迭代终止条件则停止迭代得岀最优解；否则继续迭代。非线 性规划的关键是搜索方向、步长和初值。我们用拟牛顿法来

14、选定搜索方向拟牛 顿法是在牛顿法基础上克服牛顿法中黑赛阵不仅计算复杂、而且会岀现变态、 不正定等情况，同时保持了较快收敛的优点从而得到最好的下降方向o搜索步长 的确定使用线性搜索的方法或更为有效的插值方法。由于此搜索程序的初值对 于程序正确有效的搜索影响很大，又因为相邻的切片只有一个象素单位的距离. 可以认为中轴线的变化很小所以我们可以把前一次确定的球心作为下一次搜索 的起点从而大大提髙搜索的效率和准确性。但由于转角处变化较大，在这种情况 下搜索起点会郴够逼近而导致优化搜索的失效。我们采用回溯技术来避免这种现 象。即当球心间距前后相差较大时(采用工程上的观点.以6倍为基线)回退一 次搜索并以当

15、前球心作为回退搜索的起点这样就可以把各个球心间距变得比 较匀称从而相应的消除了转角上的搜索失效。从上述对三种方法的分析可以得 到网格法和蒙特卡罗法实现思路简单，程序容易实现但网格法搜索的精度不高, 误差较大搜索时间长，同样蒙特卡罗法的实现对采点的数目要求很高计算量大; 非线性规划法的实现复杂但搜索速度快计算量少，而且通过MATLAB的优化工 具箱的函数可以很方便地实现所以非线性规划法较优。轨迹搜索结束以后我们 可以得到中轴线的100个点用所求半径的球定位在这些点上和轮廊线相交就可 以得到图4。从上图可以看到把图形重组后岀现了不平滑点也就是说直接用这些 中轴线上的点构成的管道还不平滑所以我们采用曲线拟合的方法来精确定位中 轴线。内边界重叠成的三维图像图45. 3.4绘制中轴线空间曲线图和在XOY、YOZ、X0Z三个平面的投影图。 定理1:切片上血管截面图的头部顶点在X0Y平面上的投影点一定会落在中 轴线在X0Y平面上的投影曲线上。证明：切片与X0Y平面平行，点0,点0c分别为与切片相切的两个球的球 心。M, N为两切点，0点投影在XOY平面的投影点为P,则0P垂直于平面xoy, 又M为切点所以()M垂直于切面.因为切面平行于XOY平面所以0M延长线垂 直于X