第19卷建模专辑 工程数学学报 Vol 19 Supp 2002年02月 JOURNAL OF ENGINEERING MATHEMATICS Feb.2002 文章编号:1005-3085(2002)05-002207 管道切片的三维重建 廖武鹏,邓俊晔,王丹 指导教师:数模教练组 (河海大学,南京210098) 编者按:该论文根据问题以离散形式给出数据而所求轴心轨迹及切片轮廓实质是连续曲线的特点,并充分利用生成球的某 位置在上、下半径距离范围内的切片都有截点的规律,就确定边界点及求切片最大内切圆半径提出连续算法。在 后一种算法中还讨论了某切片最大圆心实际计算中出现的不唯一的情况下如何筛选的问题 摘要:文中证明了所有切片含有过轴心的大圆,该大圆直径一定与切片边界相交。通过构造连续型模型和离散型模型,从 0BMP中定出轴心为(0,160)和半径为30的最大圆,并相继在其它切片中运用最大圆必包含在切片中的先决条 件,找出相应切片中所有可能的轴心坐标,进一步对每一切片待选的轴心坐标,根据其球体必在上29~下29层切片 中存在相应半径的圆(在上下29层中存在半径为7.68在24层存在半径为18的该球体的相应的截面圆)的特征, 筛选上述待选轴心坐标,比较准确地定出了0到70层的轴心坐标。对于71至99层由于上29层的信息不全,还存 在不少待选点,再应用切片尖端特性(在70层左下角的点只能由半径较小的圆包络而成,由此定出99层的轴心坐 标)确定其余切片的轴心坐标。绘制出的三维图形和各坐标面的投影图是光滑流畅的。最后文中用所得轴心坐标 重新构造各切片,与原切片比较,相异象素点误差不足3%,结果令人满意 关键词:连续模型;离散模型;尖端特性 分类号:AMS(2000)65D1 中图分类号:O242.1 文献标识码:A 1问题的重述(略) 2模型的假设(略)。 3问题的分析 问题第一部分需要求出管道的中轴线方程和半径,第二部分需要绘制中轴线在各个平面 上的投影。 对于第一个问题,由于一张切片的厚度为1个单位,故切片间的中轴线可以做线性化处 理,即用一条线段代替。基于100张切片所提供的信息,利用计算机搜索球体的半径,并找到 每张切片上中轴线与切片交点的坐标,称交点坐标为轴心坐标。利用这些坐标,求出两张切片 间中轴线方程。将100个轴心连接起来,形成一条完整的空间曲线。 对于第二个问题,由得出的中轴线方程,求出每段中轴线段在XY、YZ、ZX平面上的投影 (一段线段),综合每一段投影即可得中轴线在三个坐标面上的投影(将每段投影连接起来)。 考虑到实际图象边界上的点是连续的,只有位置而没有大小,且点的位置以坐标(任何实
建模专辑 管道切片的三维重建 数对)来表示。在转换成BMP格式图象时,象素表示的图象边界是离散的,一般成锯齿状,图象 范围与实际图象有误差,包括图象转换的系统误差,即点取舍引起的误差和计算的舍入误差。 针对BMP图象,我们试图就圆这一特定的图形,反演图象变换,消除系统误差,找出实际 图形理论上的圆心和半径等。就一般平面图形可根据现有计算机图形生成技术,力求较好地再 现原始图形的特征指标。由于所给问题数据是离散的,我们可以将所求中轴线看成是象素坐标 上的点,但问题本身是连续的,进一步将连续点的处理与离散的象素图形综合考虑,求出相应 的轴心坐标,并检查所求结果与所给图形的误差大小。 模型的建立 从几何特征,先明确下面重要的结论: 定理1任一切片的边界必是滚动球在切片处相应的截面图所形成的包络图。 由定理1,我们可得下面定理2 定理2任意切片Z(Z表示切片的高度),一定包含球心在Z±1,Z±2,…,Z±29位置 上球体在切片处半径分别为√R2-2(i=1,2,…,29)的截面圆(这里R设为球体半径)。 题目中存在的不确定性分析: 注意到问题所给信息不能排除下述可能直径一定的球体可以产生比球体直径大得多的 切片。如管道是由球心沿竖直方向上升较慢的螺旋线的球滚动而成的。在此情形下,实际无 法确定中轴线坐标与球半径。针对所求问题,不妨设中轴线在水平方向的螺旋绕动变化相对 竖直方向变化较慢,使得任一切片最大内含圆的半径与滚动球半径相同。最大内含圆的圆心 即为切片与中轴线的交点。 根据题目中欲求的管道半径,设管道横截面是一个圆,其半径与球体半径相同,则有 定理3若切片与中轴线有交点,且管道的横断面是圆,则该切片必含有半径与球体相同 的最大圆,且圆心在交点处。 由上述假设可建立第一个模型。 模型一、连续模型 为了寻找球体的半径,需要对100张BMP格式的图象文件进行搜索,这里涉及到一个对 BMP格式文件的处理问题。由于BMP格式文件在计算机中是以二进制数进行存储的。图象 保存在一个二维的由0或1组成的矩阵中。0和1分别对应于图象中的黑白象素。一张BMP 格式的文件包含了512×512个象素,每一个象素都有自己的一个确定的坐标。一个16进制 的数代表4个二进制位,每个二进制位可以记录一个象素。为方便使用数据,我们将16进制 的文件编程转化为由0和1组成的文本文件,然后利用该文本文件计算每张图片上管道与切 片的交所形成的边界点坐标。 在说明计算边界点坐标方法之前需引进图象处理技术中的四邻域概念 四邻域:某个象素的左、右、上、下四个象素称为该象素的四邻域,如图1,象素D,B,F H称为象素E的四邻域。 具体寻找边界点坐标的算法:由于我们将图象信息用0和1两个不同的灰度值表示黑白 象素。对图象进行逐行搜索,当遇到灰度值为零的象素点时,再搜索其四邻域内的点。若在其 四邻域中有一个象素的灰度值为1,则该点为边界点。将每张切片的边界点坐标保存在一个 二维数组中,为求解半径所用
程数学学报 第19卷 寻找球体半径的算法:由于每张切片图象有一组确定 的边界坐标。每个被截的球体在切片上的投影均是圆,在 这些圆中过球心的圆的半径最大。图象的边界是由这些大 大小小的圆包络而形成的。不妨以50.BMP(具有一般性) 为例说明问题 第i张切片的边界点坐标保存在数组x[1],x[2],…, x[n];y[1],y[2],…,y[n]中。 取 minx=minx[1],x[2],…,x「n]}, 图1象素E的四邻域 x=maxx[1],x[2],…,x[n], y=miny[1],y[2],…,y[n], max y= maxy[1l, y[2],,y[n]l 以minx,maxx,miny,maxy为边界作矩形区域 D,在D的内部逐行搜索如果遇到一个值为0的点,再搜 索其四邻域的点,如果四邻域内所有象素点所对应的值 为0,则该点一定是边界区域所包围内部图形的点,称之 为内点。设内点为A,以其为圆心作一个半径为R的圆 (R可以从零开始递增),如果边界上所有的点都不落在 图2最大圆的搜索过程 该圆内,记录下圆心的坐标,并继续增加半径R直到有一部分边界的点落入圆内为止。最后寻 找半径最大的圆(该圆必与边缘相切),与边缘相切的最大圆即是过球心的最大圆。 用 Pascal语言编程计算出每张切片上的中轴线与切线交点的坐标以及该切面上最大圆的 半径的范围。此处给出前面几个切片的交点坐标和半径的范围 表1轴心的坐标与半径 切片 革径范图29.8(2.72994(090309928194(29.829)(29.8,29 坐标(0.160.5(0.16021(013981(0.1,159.9)(0.1,101)1(01,1602 模型二、离散方法 由前定理知,所有切片中均含有半径为R的最大圆。对0.BMP所对应的切面,沿X轴方 向在其内部取定坐标(这里选用题目所给的象素坐标,即几何坐标位置理解为象素坐标的中 心)。用通用 Bresenham算法计算出某一给定圆的边界坐标点,判断其是否含于切片中。根据电 脑画图的特点,这里选用的半径、坐标都为整数。用这种方法处理象素文件可以克服由原实际 图象转换到BMP数据图象的系统误差求出所有最大圆的圆心坐标,找出最大半径R。 利用计算机图形学中的 Bresenham算法作半径为R的圆的BMP格式图象。 具体算法如下 1)初始象素点取x0=0,y0=R,取f0=1,go=2×(1-R) 2)对=0到[+3]-1,有x+1=x1+1,且 若|f1≤g:1,则y+1=y,f+1=f1+2x1+3,g+1=f+1-2y1+1, 若1f1>1g:|,则y2+1=y-1,f+1=g1+2x:+3,g+1=f1+1-2y1+3, (x1+1,y+1)即为圆边界象素点坐标
建模专辑 管道切片的三维重建 利用 WINDOWS画图软件 Paintbrush生成半径为30的圆,并与由 Bresenham算法所作半 径为30的圆进行比较,发现两者边缘的象素点重合,见图3。我们有理由相信,切片上的图象 (一系列圆叠加而成)也是利用这种方法生成的。 图3 a Paintbrush生成的图 图3 b Bresenham算法作的图 寻找最大圆算法:对0.BMP图像运用 Bresenham算法计算出最大内含圆半径是30,圆心 为(0,160)。由于采用了离散型的模型,容许一个象素的误差,因此30是可以接受的。以此为基 础去寻找1.BMP到99.BMP中的最大圆。为了加快搜索的速度,对图象进行逐行扫描时,记 录扫描线与图象边缘的最左边与最右边交点。它们的横坐标分别为x[i],x[j]。这两个边缘 点之间的距离d=1x[j]-x[i1如果d<60则停止搜索,做下一行的扫描。如果d≥60, 则判断横坐标介于x[i]与x[j]之间的所有点是否为最大圆的圆心。判断方法是以这些点为 圆心,以30为半径作一圆周,检测边缘象素点是否在圆周以内,若是则此点不是最大圆的圆心 坐标,否则该点为最大圆的圆心坐标。这些圆心的坐标为待选的轴心坐标,需对其进行筛选,筛 选出的轴心坐标所对应的最大圆为最优圆 扫描线与边缘相交时可能会遇到如图4所示 的情况。扫描线1,2,3与边缘的交点有2个,3个 和4个的情况。如扫描线2与边缘相交三次,则我 们应该选择B2与C2之间的边缘点计算距离,扫2 描线3与边缘相交四次,应取A3与B3,C3与D33 之间的点计算距离 寻找最优圆算法:管道视为球体滚动包络而 成,对于某一个确定的球体被下29~上29张切 图4扫描线与切片边界相交情况 片将其均匀分割。形成一系列的同轴圆,假设球心经过第0张切片,球顶部分被第29张切片所 截,见图5。最大圆的圆心与所截的第29个小圆圆心位于同一个轴上。最大圆的圆心坐标就是 第29个小圆的圆心坐标。 如果找到第29个小圆在第29张切片上的位置,则最大圆圆心(即中轴线与切片的交点) 便确定了 下面阐述具体步骤 对0.BMP到99BMP均有多个待选轴心坐标,在第Z层中,以待选轴心坐标所作的球体 被上下i(i=1~29)张切片所截的小圆应在第Z±i张切片中,若有任意一张切片所截的小 圆不在第Z±i张切片中,则该待选轴心坐标不可能是最优圆的轴心坐标。 0.BMP到70.BMP图象用以上方法处理效果较好。对于0~28的几张切片一般只有三 四个待选轴心坐标需要进行筛选,虽然下方切片的信息不完整,但由于待选轴心坐标较少,完
工程数学学报 第19卷 全可以筛选出最优圆。但对于71.BMP到99.BMP 第29层切片 的切片经筛选后仍有不少待选点,这是因为71.BMP 到99.BMP的信息量不够,故被淘汰的最大圆不多。 这时需应用尖端特性进行二次筛选。以99.BMP为 第0层切片 例,注意到70.BMP切片左下角因不可能是轴心在第 98层和第42层切片处的球体在70层截面小圆的包 络,而是更小截面圆的包络,这个小圆只能来自第99 和第41层中球体,经计算知第41层不存在这样的球 体,故此小圆一定存在于轴心在第99层的球体中。 其他一些靠近第99层的切片也类似处理。表2 图5截面圆与轴心的关系图 为一部分切片的轴心坐标。 表2部分轴心坐标 016010 445 下面需要构造中轴线方程,在问题分析中我们提到过,由于每张切片的厚度很薄,考虑最 简单的方法,在每段切片中的中轴线用直线段线性化表示。得到直线的点法式方程将这99个 方程联立可得完整的中轴线的方程。 问题的第二部分要求绘出中轴线在XY,YZ,ZX三个坐标平面上投影,在每段中轴线的 点法式方程中消去相应的变量,则得到该段中轴线在这个平面上的投影方程。分别联立每个坐 标面上方程就构成了三个坐标面上的中轴线投影方程。下面是我们利用 EXCEL所绘制的管 道三维模型中轴线的投影图象 50 图6轴线在XY平面上的投影图 5模型的分析和检验 由所得中轴线坐标,可生成30-70层所有切片的图形,方法如下: