建模专辑 管道切片的三维重建 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所示,我们将模版图象分为四份,分别放置在图样的四角。这种作法是为了保 证所获得的相干图样的坐标位置与原样品图样的坐标位置一致,使得在最终处理数据时更方