血管三维图像再现中的数学方法
作者：张志文 张秀梅 吴秀香
【关键词】 血管
摘要：序列图像的机三维重建是应用数学和计算机技术在医学与生物学领域的重要应用之一。依据一根血管的一组平行切片图，运用有效的数学方法，在计算机上成功再现了其三维图像。
关键词：位图；三维重建；中轴线；曲线拟合
血管三维重建问题来源于序列图像的计算机三维重建［1，2］。序列图像的计算机三维重建是应用数学和计算机技术在医学与生物学领域的重要应用之一。
为研究生物体的复杂结构，常将其本身或局部做成切片。切片图像序列把生物体内部的各种复杂结构和变化一层层地暴露出来了，人们通过依次对每张切片图像的观察、分析和比较，综合起来可以形成对生物体内部结构的立体认识。从几何角度看，这种综合就是由切片图像序列恢复生物体内部结构的几何形状，称此为序列图像的三维重建。这项工作，过去是在人脑中进行的，专业人员通过观察，凭经验在自己头脑里想象出生物体的内部结构和几何形状。在今天当然把这项繁杂的工作交由计算机完成，实行序列图像的三维重建的计算机化、自动化。序列图像的计算机三维重建是切片制作的逆过程，很复杂，需要综合运用图像处理、图形学、计算机辅助几何设计等多学科的方法，是当前研究的前沿和热点课题之一。
血管是血液流通的通路，其在生命活动中的重要性众所周知，诊断师在临床中经常需要了解血管的分布、走向等重要信息。理想的血管可以看成是粗细均匀的管道，如何建立其数学模型是图像三维重建的重要一环。
2001年全国大学生数学建模竞赛题目为“血管三维重建问题”，要较好解决该问题，远非建模竞赛要求的三天所能完成。因此，竞赛结束后产生的优秀中都存在不完善之处，本研究的目的正是基于此而产生的。
1 问题的提出
假设某些血管可视为一类特殊的管道，该管道的表面是由球心沿着某一曲线（称为中轴线）的球滚动包络而成。例如圆柱就是这样一种管道，其中轴线为直线，由半径固定的球滚动包络形成。
现有某管道的相继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）。
为了在计算机上再现血管的三维形态，需要计算管道的中轴线与半径。
2 管道半径的求解
求解管道半径有多种方法。但建模竞赛的优秀论文中以“平均法”和“抽样法”居多［1］。所谓平均法就是求出每张横断面图像内的最大内切圆半径，再取管道半径为它们的算术平均值。此方法的缺点是要以每张图像内的每个象素点作为圆心，令每张横断面图像内的内切圆半径值由小到大，动态地逼近最大内切圆半径的求解过程，其计算量相当庞大，计算机程序运行困难，优秀论文的作者中许多人提到这一点。
所谓抽样法就是利用滚动球半径是常数，取前几片横断面图像内的最大内切圆半径的平均值为管道半径的值。此方法正是意识到平均法计算量的庞大而提出的。其缺点是抽样样本数目选取的合理性较难确定，因样本数目少而存在计算误差。
本研究求解管道半径的方法是：首先将100张图片叠合，形成如图2(a)所示的一张图片。由于切片垂直Z轴的，此图片是血管在XOY面上的投影图像。因此，它也是滚动球在XOY面上的投影――滚动圆，沿中轴线的投影线滚动形成的二维包络图， 且滚动球的半径与滚动圆的半径相同，因此只需求出滚动圆半径。具体图像叠合的方法是首先运用Photoshop软件打开0.bmp 图片，再将1.bmp~99.bmp图片（共99张图片）通过叠加控件叠合。
为简化求半径的复杂程度，取图2(a)区域的一部分（图2(b)）。利用计算机编程搜索图像边缘点得到边缘曲线AB和AB，记录其坐标,即两条曲线上象素点的坐标(个数有限)。在凸弧AB上任取一点M,扫描计算凹弧AB上所有点到M点的距离，确定其最小值，此即所求滚动圆直径，本研究计算的半径结果是29.5个单位。
此方法的主要优点是选点M是任意的，即无限制，而且只需扫描AB一侧所有点，大大简化计算量。事实上在编程计算时，还可使弧AB更短，这样做显然并不改变计算的精确性。
查阅到的获奖论文中，大多是对一张图片（而非叠合图片）进行扫描，计算涉及两侧间所有点间距离，还要涉及最小和最大值的比较问题，然后对几张图片做上述类似工作，再取均值。相对而言，本研究方法无论是在计算量还是在精确度上都较为优化。
3 中轴线方程的求解
31 切片图像最大内切圆的存在性证明
依据基本假设：视血管表面为一类由球心沿着某一曲线的球滚动包络而成的特殊的管道，且管道中轴线与每一张切片有且只有一个交点，可知滚动球是沿中轴线严格单调上升的。
切片图像的产生是滚动球上升过程中与该切片所在平面相交而产生的一系列二维区域（圆域）的并集，从球与该平面接触到球离开该平面过程中，由滚动球严格单调上升性，球心与该平面相交且仅相交一次，此时形成最大的圆域，亦即该切片图像内存在唯一最大内切圆。
上述证明过程表明最大内切圆的圆心为中轴线与该切片的交点，最大内切圆的半径为滚动球的半径。
32 切片图像最大内切圆圆心的求解
求解每张切片图像最大内切圆圆心的方法有多种。但建模竞赛的优秀中以“枚举法”和“平行切线法”居多［1］。
所谓枚举法就是求每张横断面的图像内的最大内切圆的圆心时，以位于图像内每一个象素为圆心作圆．遍历所有象素点后再作确定。此种方法，由于每次做圆的过程半径是由小到大动态的，最大圆是经过两次循环获得的，量巨大。
所谓平行切线法就是横断面的图像边界上的两点的连线如果同时垂直边界在这两点处的切线，则这两点连线有可能是最大内切圆的直径。发现所有具有这样性质的点对，并检验之，以确定最大内切圆的圆心。此方法的缺点是缺乏合理性的理论证明，事实上未见有论文对此证明。另外在象素表示边界线这种离散状态下，判断“垂直”性方法并不明确。
本研究求每张切片最大内切圆的具体方法：首先扫描确定区域，如图3阴影部分所示。再以区域内每一象素点为圆心，29.5为半径按照计算机图形学中圆的Bresenham算法［3］做圆，尽管圆心变化，但因半径固定，此算法得到的圆周象素点个数及圆周的象素表示形状固定。程序中做一计数器，将区域内所有象素点的初始值赋值为0，每次做圆过程中，将在圆周上的区域内的点的计数器值加1，扫描区域内所有点后，区域内计数器值最大的点即为最大内切圆圆心。
上述做法的合理性是：如图3所示，圆O是最大内切圆，而圆O为某一个非最大圆。由于半径是固定的，圆心象素点的计数器值增加当且仅当它位于以它的圆周上区域点为圆心的圆周上。由此可知圆心点的计数器值等于其圆周上点位于区域内象素点的个数。从图3可看出，圆O的圆心计数器值明显大于圆O的圆心计数器的值。
该方法的优点是最大内切圆的圆周不必全部在区域内部，即最大内切圆的圆心象素点的计数器的值可能小于其圆周点象素的个数，而只需计数器的值最大即可。这样可忽略真实切片边缘进行数码转换（象素显示）时的误差。在bmp格式下，二维切片边缘并不是平滑的曲线，而是齿状形式，在对血管切片获取的过程中就有误差产生，不应该通过判断其是否全在黑色区域内来找到最大内切圆，如果某一圆只有一象素点不在黑色区域，可能该圆是实际血管切片（未转化为bmp格式之前血管切片）的最大内切圆，且实际只能存在一个这样的圆。
33 中轴线方程的拟合求解
利用所求得的最大内切圆圆心坐标，通过Matlab软件中的曲线拟合函数ployfit得到中轴线方程。由于方程的特殊性，本研究采取分段拟合的方式。
计算结果为：t=0~29X(t)=0.00090547846488*t.^3+0.0067784235874*t.^2+0.0784535895423*t.-0.295780088597090
Y(t)=-0.006008085052*t.^2+0.147566585759*t+149.264564456547Z(t)=t
t=30~99X(t)=0.000036785469*t^4-0.010418543628*t.^3+1.057784235874*t.^2-35.057946695578*t.+469.5578462318
Y(t)=0.0000458326*t^4+0.00098784455346*t^3-0.345864521232*t.^2+22.175546489658987*t-202.11102121284
Z(t)=t
依据上述方程，利用Matlab软件或Pro/engineer软件都可在计算机上再现该血管的中轴曲线及血管本身的三维图像。
1 汪国昭，陈凌均．血管三维重建的问题．工程数学学报，2002，19 supp（2）： 55～58
2 应荣超．人下颌下腺淋巴管、血管和腺导管的计算机三维重建技术．解剖学报，1991，22 （4）：342～346
3 潘云鹤． 计算机图形学―原理、方法及应用．第一版．高等出版社，2002，37～45
4 郎锐．数字图像处 visual C++实现． 北京希望出版社，2002，58～64
5 张志涌 ．精通Matlab6.5． 第一版． 北京航空航天大学出版社，2003，57～68
