第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层所有切片的图形,方法如下:
建模专辑 管道切片的三维重建 20 -150100 图7轴线在YZ平面上的投影图 140 X80 图8轴线在ZX平面上的投影图 欲产生第Z层切片,由Z-29,2-28,…,z,z+1,…,Z+29相应轴心坐标计算出对应 的所有圆的象索点,并将它们迭加便可形成第Z层切片。进一步,将所得图形与原切片作比较 计算相异点所占比例,经计算,第40层切片的相对误差为3%。这说明我们的模型精度是较高 但是建立第一个模型是比较粗糙的。通过控制 半径的步长大小,可以得到一系列半径非常接近的 最大圆,以此我们能够给出对于一张切片上最大圆 半径的一个范围,发现多数是靠近30这个理想半径 的。出现这种问题的原因是由于切片上的图象是由 离散的象素点表示的。不同方向的半径是不相同 的,见图9。由于问题中需要求解中轴线的方程,这图9aBMP格式存储的圆图9b理论圆 条曲线在数学上是没有粗细且唯一确定的。因此对于中轴线与切片的交点应该是唯一的。但 是用来表示交点的象素点经筛选后仍不唯一。可选这些象素点的几何中心的整数点为轴心坐 标 由于实际的图形转化为BMP格式图象文件时会造成图象的失真,正如前面提到过计算 机是根据一定的算法去存储图象的,并且算法的对象是离散的点,这对于实际连续的实物来说 误差不可避免,所以是一种系统误差。而模型二所采用的算法恰是对离散BMP象素点的处 理,可以根本消除系统误差
工程数学学报 第19卷 6模型的推广 如果考虑中轴线与切片的交点不止一个,比如有两个交点,仍然可以利用我们的离散模型 求解,因为不含轴心的切片中一定不含有最大圆。 我们对于离散问题讨论较多,但对于轴心坐标及半径是实数的情况,没有过多的涉及。问 题本身是连续的,用离散方法去模拟,不可避免地会出现误差。如果知道象素圆的生成算法 我们就可以运用本模型的思想,将模型推广到实数领域,从而使模型的精度提高。 参考文献 [1]管伟光体视化技术及其应用[M],北京:电子工业出版社,1998 [2]孙家广,扬长贵计算机图形学[M].北京:清华大学出版社,1995 [3] Kenneth.R. castleman,数字图象处理[M]北京:电子工业出版社,1981 3D Reconstruction of Pipeline Slice LIAO Wu-peng, DENG Jun-ye, WANG Dan Advisor: Mathematical Modeling Tutor Group Hohai University, Nanjing 210098 Abstract: This article proves that there exists the biggest cirele on each slice, whose center is on the pipeline center curve. By con- tructing the continuous model and the discrete one, we are able to determine the biggest circle on 0. bmp slice, of which radius is 30 and center coordinate is(0, 160). Using the predetermined condition that the biggest cirele should be contained in each slice, we find out all possible coordinates of the pipeline center curve on corresponding slices. For cach slice, basing on the principle that slice which lies on over 29-under 29 layers of the slice must contain a corresponding section of the sphere of radius 30, we filtrate the coordinates of the pipeline center curve and determine the more precise one. However, there are still some coordinates of the pipeline center curve left to be filtrated, such as slices from 71 to 99, due to the information is not enough. Further, we use the pointed end property to determine the coordinates of the pipeline center curve on other slices, The projecting figures on each coordi ate plane are smooth and fluent. In the end, by using the coordinates to reconstruct all the slices and comparing the slices with the original ones, we found that the error of different pixels is less than 3%. The result is satisfied Key words: continuous model: discrete model: pointed end property
第19卷建模专辑 工程数学学报 Vol 19 Sup 2002年02月 JOURNAL OF ENGINEERING MATHEMATICS Feb.2002 文章编号:10053085(2002)05-0029-06 利用切片的二维空间相关操作实现 血管的三维重建 胡亦斌,向杰,程翔 指导老师:王若鹏 (北京大学物理系,北京100871) 编者按:本文利用相关作为两个图形相似程度的度量,通过快速傅立叶变换及反变换对切片进行空间相关操作,具有一定特 色。这种方法能够快速确定出中轴线与切片的交点,所得到的计算结果较为准确。请注意文中X轴的方向与题意 要求是相反的。 摘要:相关( Correlation)作为两个图形相似程度的度量,被广泛的用于图形图像自动识别中。为对血管的二维切片图像进 行分析并重构出血管以及血管中轴线的三维空间形貌,我们利用快速傅立叶变换(FFT)及反变换对切片进行空间 相关操作,几乎一步即可确定出中轴线与切片的交点,从而给出中轴线的空间坐标。我们求出了血管的半径,利用 这些结果,绘出了血管中轴线的三维曲线及其投影线,并且利用计算机软件画出了血管的三维造型,在该造型中作 血管切片,结果与初始的切片数据一致。文中分析了相关法进行图像处理的优点与局限性,对利用近代光学信息处 理的手段进行切片三维重建的思路进行了讨论。 关键字:傅立叶变换,相关操作,三维重建 分类号:AMS(2000)65D17 中图分类号:0242.1 文献标识码:A 1问题分析与半径的获得 血管可以看成无穷多个等径并且圆心相距无穷小的球包络面组成。因此,切片上的二维 图形就应该是由无穷多个球被截的圆叠加而成。这些圆都是被截球的大圆或者小圆,其半径 有一极大值R,R同时也是球的半径。这样一个半径R的圆是球心在切片平面内的球被截而 成的,其圆心为中轴线与切片平面的交点。假设,中轴线与每张切片有且只有一个交点,所以 每一张切片图上包含且只包含一个半径为R的圆。我们只要找到这个圆,就可以定出中轴线 与切片平面交点的坐标,用这些交点坐标我们可以建立起中轴线的空间形态。 现在我们把问题分解为两部分:第一步,找到半径R;第二步,在每一张切片图中找到半径 为R的圆的圆心。第二个问题我们使用了空间相关操作方法来解决,这将在下一部分中详细 描述。我们先来解决第一个问题 图1是一张叠加图,是由0.bmp到99.bmp的相应像素进行或运算所得。这张图代表了 待求血管在X-Y平面上的投影的形状。我们不难看出,除去这一“弯月”的两端之外,这一图 形的宽度是保持不变的,这一宽度就是我们所要求的球的直径2R。我们找到直径的方法是这 样的:在图形的中段也就是宽度很均匀部分的中段(y=260-280)左侧固定一个在轮廓上的点
第19卷 A,在A点对应的另一侧轮廓线上从与A同一Y坐标的 点向附近寻找一点B,使得|AB在这些点中最短,可 以证明|AB|就是我们要求的直径。这样,我们就可以求 出2R=60,R=30(像素)。 2血管中轴线的确定 在已获得球半径的基础上,我们对切片图样进行相 关操作来确定半径为R的圆的圆心的位置 斤谓相关操作是指利用两函数f(x,y)、g(x,y) 作相关积分时的特性,在样品图样中定位某种模版图样 的位置方法我们以一简单情形为例简单的说明相关操 作。假设样品具有如下图形: 图1所有血管切片的叠加图 f(x)=2a(x-x),其中g(x)为一个单位方势垒函数:g(x)=11x15a2是若 干个在x轴上的离散分布的固定点。为了求出在样品f(x)中,和g(x)相同的图样的位置。我 们可以对这两个函数进行相关运算 g(x)f(x)=∫g(x’-x)f(x)dx ∫g(x’-x)g(x'-x;)dx,(1) 作变量替换:x”=x’-x1,则 g(x)④f(x)=∑g(x”+x1-x)g(x”)d /mg(x”-(x-x1))g(x)"dx” (2) 其中c(x)=g(x)g(x),为角形函数。凡原来在f(x)中有方垒的地方,g(x)④f(x)中 就有一个尖峰尖峰的位置正是各离散点的中心,也就是在f(x)图中g(x)所在位置。 粗略的讲,相关函数是进行相关运算的两个函数重叠程度的描述,在完全重合处,相关有 极大值。一般说来,由于相同的函数曲线才能完全重叠,故它们的自相关往往比不同函数间 的关联强得多由上面的简单例子可以看到,当f(x)代表一些相同的信息g(x)的离散排列 时,相关函数g(x)f(x)的意义类似于用g(x)在f(x)中去搜索,凡遇到与自身相同的信 息之处记录下来一个信号峰。从这一简单例子中得到的结论可以很容易的推广到高维复杂情 形 由于g(x,y)f(x,y)台G‘(a,v)F(n,v),其中公式中台两边表示一对傅立叶变换 对,G、F分别表示g(x,y)和f(x,y)的傅立叶变换函数。因此,在具体计算上,我们可以先将 样品图样进行傅立叶变换求出它的频谱,再将其与模版图样的频谱的复共轭相乘,并将积作反 傅立叶变换,即可得到所求的两图样的相关 我们将模版g(x,y)取为大圆,用此模版在切片图样f(x,y)中搜索相同的图样。由大圆 存在的唯一性可知,这样找到的唯一一个最大亮点就应是大圆的圆心,也即中轴线在此平面上 的坐标,这样问题也就得到了解决。我们利用 MATLAB提供的傅立叶变换和反傅立叶变换的
建模专辑 利用切片的二维空间相关操作实现血管的三维重建 函数,就可以方便的求解此问题,整个计算在处理上是非常简单快捷的。 在具体选择模版时还有以下两点是值得注意的: 首先,我们选择半径为R=30的空心圆而非实心圆作为模版(见图2)。之所以这样做是 因为会使所得相关图样的对比度较大,从而更加突出满足条件的高亮区。出现这种情况的原因 是因为,如前面所述,两个图样的相关函数的大小是由两个图样的交叠程度决定的,当我们选 取实心圆作为模版时,尽管模版和大圆的交叠程度很大,但是样品其它区域和模版圆的交叠程 度也很大,尤其是和大圆相近的区域,它们交叠程度与大圆交叠程度的差别仅在于模版圆的边 界,而这与交叠部分相比只占很小一部分,因而在相干图样上产生的亮度差别是很小的。而如 果采用空心圆作为模版,一方面,在整个图形中只有大圆能够和它整个的相叠,产生最大交叠 另一方面,样品图中其它区域和模版图样的未交叠部分占整个模版图形的比重就会大大的增 加,从而使得不同区域交叠程度的很小差别就会起到很大作用,体现在结果中就是使得图样的 对比度明显的加大,相关信号的尖峰更加明显 图2上面两个图为两种R=30的模板,左侧为空心圆模板、右侧为实心圆模板 下面两个图为使用两种不同的模板对叠加“弯月”图作相关操作的结果 左侧为空心圆模板结果,右侧为实心模板的结果 第二,如图2所示,我们将模版图象分为四份,分别放置在图样的四角。这种作法是为了保 证所获得的相干图样的坐标位置与原样品图样的坐标位置一致,使得在最终处理数据时更方
