工程数学学报 全可以筛选出最优圆。但对于71.BMP到99BMP 第29层切片 的切片经筛选后仍有不少待选点,这是因为71.BMP 到99BMP的信息量不够,故被淘汰的最大圆不多。 这时需应用尖端特性进行二次筛选。以99.BMP为 第0层切片 例,注意到70.BMP切片左下角因不可能是轴心在第 98层和第42层切片处的球体在70层截面小圆的包 络,而是更小截面圆的包络,这个小圆只能来自第99 和第41层中球体,经计算知第41层不存在这样的球 体,故此小圆一定存在于轴心在第99层的球体中 其他一些靠近第99层的切片也类似处理。表2图5截面圆与轴心的关系图 为一部分切片的轴心坐标。 表2部分轴心坐标 Y 下面需要构造中轴线方程,在问题分析中我们提到过,由于每张切片的厚度很薄,考虑最 简单的方法,在每段切片中的中轴线用直线段线性化表示。得到直线的点法式方程将这99个 方程联立可得完整的中轴线的方程 问题的第二部分要求绘出中轴线在XY,YZ,ZX三个坐标平面上投影,在每段中轴线的 点法式方程中消去相应的变量,则得到该段中轴线在这个平面上的投影方程分别联立每个坐 标面上方程就构成了三个坐标面上的中轴线投影方程。下面是我们利用EXCE所绘制的管 道三维模型中轴线的投影图象。 图6轴线在Xy平面上的投影图 5模型的分析和检验 由所得中轴线坐标,可生成30-70层所有切片的图形,方法如下 o1995-2004 Tsinghua Tongfang Optical Disc Co, LId. All rights reserved
图 5 截面圆与轴心的关系图 全可以筛选出最优圆。但对于 71. BMP 到 99. BMP 的切片经筛选后仍有不少待选点 ,这是因为 71. BMP 到 99. BMP 的信息量不够 ,故被淘汰的最大圆不多。 这时需应用尖端特性进行二次筛选。以 99. BMP 为 例 ,注意到 70. BMP 切片左下角因不可能是轴心在第 98 层和第 42 层切片处的球体在 70 层截面小圆的包 络 ,而是更小截面圆的包络 ,这个小圆只能来自第 99 和第 41 层中球体 ,经计算知第 41 层不存在这样的球 体 ,故此小圆一定存在于轴心在第 99 层的球体中。 其他一些靠近第 99 层的切片也类似处理。表 2 为一部分切片的轴心坐标。 表 2 部分轴心坐标 X Y Z 0 160 0 0 160 1 0 160 2 0 160 3 0 160 4 0 160 5 114 116 50 119 111 51 123 106 52 X Y Z 128 101 53 132 95 54 136 89 55 140 83 56 144 76 57 147 70 58 151 62 59 153 55 60 71 - 170 90 X Y Z 65 - 173 91 57 - 176 92 51 - 179 93 44 - 181 94 37 - 183 95 31 - 185 96 24 - 187 97 17 - 188 98 12 - 189 99 下面需要构造中轴线方程 ,在问题分析中我们提到过 ,由于每张切片的厚度很薄 ,考虑最 简单的方法 ,在每段切片中的中轴线用直线段线性化表示。得到直线的点法式方程将这 99 个 方程联立可得完整的中轴线的方程。 问题的第二部分要求绘出中轴线在 X Y , Y Z , ZX 三个坐标平面上投影 ,在每段中轴线的 点法式方程中消去相应的变量 ,则得到该段中轴线在这个平面上的投影方程。分别联立每个坐 标面上方程就构成了三个坐标面上的中轴线投影方程。下面是我们利用 EXCEL 所绘制的管 道三维模型中轴线的投影图象。 图 6 轴线在 X Y 平面上的投影图 5 模型的分析和检验 由所得中轴线坐标 ,可生成 30 - 70 层所有切片的图形 ,方法如下 : 62 工 程 数 学 学 报 第 19 卷 © 1995-2004 Tsinghua Tongfang Optical Disc Co., Ltd. All rights reserved
建模专辑 管道切片的三维重建 图7轴线在YZ平面上的投影图 X80 0 图8轴线在ZX平面上的投影图 欲产生第z层切片,由Z-29,z-28,…z,Z+1,…z+29相应轴心坐标计算出对应 的所有圆的象素点,并将它们迭加便可形成第z层切片进一步,将所得图形与原切片作比较 计算相异点所占比例,经计算,第40层切片的相对误差为3%这说明我们的模型精度是较高 但是建立第一个模型是比较粗糙的。通过控制 半径的步长大小,可以得到一系列半径非常接近的 最大圆,以此我们能够给出对于一张切片上最大圆 半径的一个范围,发现多数是靠近30这个理想半径 的。出现这种问题的原因是由于切片上的图象是由 离散的象素点表示的。不同方向的半径是不相同 的,见图9。由于问题中需要求解中轴线的方程,这图9aBMP格式存储的圆图9b理论圆 条曲线在数学上是没有粗细且唯一确定的。因此对于中轴线与切片的交点应该是唯一的。但 是用来表示交点的象素点经筛选后仍不唯一。可选这些象素点的几何中心的整数点为轴心坐 由于实际的图形转化为BMP格式图象文件时会造成图象的失真,正如前面提到过计算 机是根据一定的算法去存储图象的,并且算法的对象是离散的点,这对于实际连续的实物来说 误差不可避免,所以是一种系统误差。而模型二所采用的算法恰是对离散BMP象素点的处 理,可以根本消除系统误差。 o1995-2004 Tsinghua Tongfang Optical Disc Co, LId. All rights reserved
图 7 轴线在 Y Z 平面上的投影图 图 8 轴线在 ZX 平面上的投影图 欲产生第 Z 层切片 ,由 Z - 29 , Z - 28 , …, Z , Z + 1 , …, Z + 29 相应轴心坐标计算出对应 的所有圆的象素点 ,并将它们迭加便可形成第 Z层切片。进一步 ,将所得图形与原切片作比较 , 计算相异点所占比例 ,经计算 ,第 40 层切片的相对误差为 3 %。这说明我们的模型精度是较高 的。 图 9a B M P 格式存储的圆 图 9b 理论圆 但是建立第一个模型是比较粗糙的。通过控制 半径的步长大小 ,可以得到一系列半径非常接 近的 最大圆 ,以此我们能够给出对于一张切片上最大圆 半径的一个范围 ,发现多数是靠近 30 这个理想半径 的。出现这种问题的原因是由于切片上的图象是由 离散的象素点表示的。不同方向的半径是不相同 的 ,见图 9。由于问题中需要求解中轴线的方程 ,这 条曲线在数学上是没有粗细且唯一确定的。因此对于中轴线与切片的交点应该是唯一的。但 是用来表示交点的象素点经筛选后仍不唯一。可选这些象素点的几何中心的整数点为轴心坐 标。 由于实际的图形转化为 BMP 格式图象文件时会造成图象的失真 ,正如前面提到过计算 机是根据一定的算法去存储图象的 ,并且算法的对象是离散的点 ,这对于实际连续的实物来说 误差不可避免 ,所以是一种系统误差。而模型二所采用的算法恰是对离散 BMP 象素点的处 理 ,可以根本消除系统误差。 建模专辑 管道切片的三维重建 72 © 1995-2004 Tsinghua Tongfang Optical Disc Co., Ltd. All rights reserved
工程数学学 第19卷 6模型的推广 如果考虑中轴线与切片的交点不止一个,比如有两个交点,仍然可以利用我们的离散模型 去求解,因为不含轴心的切片中一定不含有最大圆 我们对于离散问题讨论较多,但对于轴心坐标及半径是实数的情况,没有过多的涉及。问 题本身是连续的,用离散方法去模拟,不可避免地会出现误差。如果知道象素圆的生成算法, 我们就可以运用本模型的思想,将模型推广到实数领域,从而使模型的精度提高。 参考文献 []管伟光体视化技术及其应用[M]北京:电子工业出版社,1998 [2]孙家广扬长贵计算机图形学M]北京:清华大学出版社,1995 [3] Kenneth.R. Castleman.数字图象处理[M]北京:电子工业出版社,1981 3D Reconstruction of Pipeline Slice L IAO Wurpeng, DENGJurrye, WANG Dan Advisor: Mathematical Modeling Tutor Grou Hohai University, Nanjing 210098) Abstract: This article proves that there exists the biggest circle on each slice, whose center is on the pipeline center curve. By structing 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 circle should be contained in each slice, we find out all possible coordinates of the pipeline center curve on corresponding slices. For each slice, basing on the principle that a 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 he 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 sices. The projecting figures on each coordi- nate 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 o1995-2004 Tsinghua Tongfang Optical Disc Co, LId. All rights reserved
6 模型的推广 如果考虑中轴线与切片的交点不止一个 ,比如有两个交点 ,仍然可以利用我们的离散模型 去求解 ,因为不含轴心的切片中一定不含有最大圆。 我们对于离散问题讨论较多 ,但对于轴心坐标及半径是实数的情况 ,没有过多的涉及。问 题本身是连续的 ,用离散方法去模拟 ,不可避免地会出现误差。如果知道象素圆的生成算法 , 我们就可以运用本模型的思想 ,将模型推广到实数领域 ,从而使模型的精度提高。 参考文献 : [ 1 ] 管伟光. 体视化技术及其应用[ M ] . 北京 :电子工业出版社 ,1998 [ 2 ] 孙家广 ,扬长贵. 计算机图形学[ M ] . 北京 :清华大学出版社 ,1995 [ 3 ] Kenneth. R. Castleman. 数字图象处理[ M ] . 北京 :电子工业出版社 ,1981 3D Reconstruction of Pipeline Slice L IAO Wu2peng , DEN G J un2ye , WAN G Dan Advisor : Mathematical Modeling Tutor Group ( Hohai University , Nanjing 210098 ) Abstract : This article proves that there exists the biggest circle on each slice , whose center is on the pipeline center curve. By con2 structing 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 circle should be contained in each slice , we find out all possible coordinates of the pipeline center curve on corresponding slices. For each slice , basing on the principle that a 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 coordi2 nate 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 82 工 程 数 学 学 报 第 19 卷 © 1995-2004 Tsinghua Tongfang Optical Disc Co., Ltd. All rights reserved
第19卷建模专辑 工程数学学报 Vol 19 002年02月 JOURNAL OF EN GIN EERIN G MA THEMA TICS 文章编号:1005-3085(2002)05-0029-06 利用切片的二维空间相关操作实现 血管的三维重建 胡亦斌,向杰,程翔 指导老师:王若鹏 (北京大学物理系,北京100871) 编者按本文利用相关作为两个图形相似程度的度量,通过快速傅立叶变换及反变换对切片进行空间相关操作,具有一定特 色。这种方法能够快速确定出中轴线与切片的交点,所得到的计算结果较为准确。请注意文中X轴的方向与题意 要求是相反的 摘要:相关 Correlation)作为两个图形相似程度的度量,被广泛的用于图形图像自动识别中。为对血管的二维切片图像进 行分析并重构出血管以及血管中轴线的三维空间形貌,我们利用快速傅立叶变换(FFT)及反变换对切片进行空间 相关操作,几乎一步即可确定出中轴线与切片的交点,从而给出中轴线的空间坐标。我们求出了血管的半径,利用 这些结果,绘出了血管中轴线的三维曲线及其投影线,并且利用计算机软件画出了血管的三维造型,在该造型中作 管切片,结果与初始的切片数据一致。文中分析了相关法进行图像处理的优点与局限性,对利用近代光学信息处 理的手段进行切片三维重建的思路进行了讨论 关键字:傅立叶变换,相关操作,三维重建 分类号:AMS(2000)65D17 中图分类号:O242.1 文献标识码:A 1问题分析与半径的获得 血管可以看成无穷多个等径并且圆心相距无穷小的球包络面组成。因此,切片上的二维 图形就应该是由无穷多个球被截的圆叠加而成。这些圆都是被截球的大圆或者小圆,其半径 有一极大值R,R同时也是球的半径。这样一个半径R的圆是球心在切片平面内的球被截而 成的,其圆心为中轴线与切片平面的交点。假设,中轴线与每张切片有且只有一个交点,所以 每一张切片图上包含且只包含一个半径为R的圆。我们只要找到这个圆,就可以定出中轴线 与切片平面交点的坐标,用这些交点坐标我们可以建立起中轴线的空间形态。 现在我们把问题分解为两部分:第一步,找到半径R;第二步,在每一张切片图中找到半径 为R的圆的圆心。第二个问题我们使用了空间相关操作方法来解决,这将在下一部分中详细 描述。我们先来解决第一个问题。 图1是一张叠加图,是由0bmp到99.bmp的相应像素进行或运算所得。这张图代表了 待求血管在ⅹ-Y平面上的投影的形状。我们不难看出,除去这一“弯月”的两端之外,这一图 形的宽度是保持不变的,这一宽度就是我们所要求的球的直径2R。我们找到直径的方法是这 样的在图形的中段也就是宽度很均匀部分的中段(y=260~280)左侧固定一个在轮廓上的点 o1995-2004 Tsinghua Tongfang Optical Disc Co, LId. All rights reserved
第 19 卷 建模专辑 2002 年 02 月 工 程 数 学 学 报 JOURNAL OF EN GIN EERIN G MA THEMA TICS Vol. 19 Supp. Feb. 2002 文章编号 :100523085 (2002) 0520029206 利用切片的二维空间相关操作实现 血管的三维重建 胡亦斌 , 向 杰 , 程 翔 指导老师 : 王若鹏 (北京大学物理系 ,北京 100871) 编者按 :本文利用相关作为两个图形相似程度的度量 ,通过快速傅立叶变换及反变换对切片进行空间相关操作 ,具有一定特 色。这种方法能够快速确定出中轴线与切片的交点 ,所得到的计算结果较为准确。请注意文中 X 轴的方向与题意 要求是相反的。 摘 要 :相关(Correlation) 作为两个图形相似程度的度量 ,被广泛的用于图形图像自动识别中。为对血管的二维切片图像进 行分析并重构出血管以及血管中轴线的三维空间形貌 ,我们利用快速傅立叶变换( FFT) 及反变换对切片进行空间 相关操作 ,几乎一步即可确定出中轴线与切片的交点 ,从而给出中轴线的空间坐标。我们求出了血管的半径 ,利用 这些结果 ,绘出了血管中轴线的三维曲线及其投影线 ,并且利用计算机软件画出了血管的三维造型 ,在该造型中作 血管切片 ,结果与初始的切片数据一致。文中分析了相关法进行图像处理的优点与局限性 ,对利用近代光学信息处 理的手段进行切片三维重建的思路进行了讨论。 关键字 :傅立叶变换 ,相关操作 ,三维重建 分类号 : AMS(2000) 65D17 中图分类号 : O24211 文献标识码 : A 1 问题分析与半径的获得 血管可以看成无穷多个等径并且圆心相距无穷小的球包络面组成。因此 ,切片上的二维 图形就应该是由无穷多个球被截的圆叠加而成。这些圆都是被截球的大圆或者小圆 ,其半径 有一极大值 R ,R 同时也是球的半径。这样一个半径 R 的圆是球心在切片平面内的球被截而 成的 ,其圆心为中轴线与切片平面的交点。假设 ,中轴线与每张切片有且只有一个交点 ,所以 每一张切片图上包含且只包含一个半径为 R 的圆。我们只要找到这个圆 ,就可以定出中轴线 与切片平面交点的坐标 ,用这些交点坐标我们可以建立起中轴线的空间形态。 现在我们把问题分解为两部分 :第一步 ,找到半径 R ;第二步 ,在每一张切片图中找到半径 为 R 的圆的圆心。第二个问题我们使用了空间相关操作方法来解决 ,这将在下一部分中详细 描述。我们先来解决第一个问题。 图 1 是一张叠加图 ,是由 0. bmp 到 99. bmp 的相应像素进行或运算所得。这张图代表了 待求血管在 X - Y平面上的投影的形状。我们不难看出 ,除去这一“弯月”的两端之外 ,这一图 形的宽度是保持不变的 ,这一宽度就是我们所要求的球的直径 2R。我们找到直径的方法是这 样的 :在图形的中段也就是宽度很均匀部分的中段(y = 260～280) 左侧固定一个在轮廓上的点 © 1995-2004 Tsinghua Tongfang Optical Disc Co., Ltd. All rights reserved
工程数学学报 第19卷 A,在A点对应的另一侧轮廓线上从与A同一Y坐标的 点向附近寻找一点B,使得|AB在这些点中最短,可 以证明AB|就是我们要求的直径。这样我们就可以求 出2R=60,R=30(像素) 2血管中轴线的确定 在已获得球半径的基础上,我们对切片图样进行相 关操作来确定半径为R的圆的圆心的位置 所谓相关操作是指利用两函数f(x,y)、g(x,y 作相关积分时的特性,在样品图样中定位某种模版图样 的位置方法。我们以一简单情形为例简单的说明相关操 作。假设样品具有如下图形 图1所有血管切片的叠加图 f(x)=(x-xb,其中8(为一个单位方势垒函数:g(x)01x1<u2+是若 干个在x轴上的离散分布的固定点为了求出在样品f(x)中,和g(x)相同的图样的位置我 们可以对这两个函数进行相关运算 (x)⑨f(x)=∫:∞g(x-x)fx 作变量替换:x”=x’-x1,则 8(x)Of(x) (x”+x1-x)g(x”)d ∞g(x”-(x-x)g(x) 其中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,yG'(u,)F(u,v),其中公式中两边表示一对傅立叶变换 对,G、F分别表示g(x,y和f(x,y的傅立叶变换函数因此,在具体计算上我们可以先将 样品图样进行傅立叶变换求出它的频谱,再将其与模版图样的频谱的复共轭相乘,并将积作反 傅立叶变换,即可得到所求的两图样的相关 我们将模版g(x,y取为大圆,用此模版在切片图样f(x,y中搜索相同的图样。由大圆 存在的唯一性可知,这样找到的唯一一个最大亮点就应是大圆的圆心,也即中轴线在此平面上 的坐标,这样问题也就得到了解决我们利用 MATLAB提供的傅立叶变换和反傅立叶变换的 o1995-2004 Tsinghua Tongfang Optical Disc Co, LId. All rights reserved
图 1 所有血管切片的叠加图 A ,在 A 点对应的另一侧轮廓线上从与 A 同一 Y坐标的 一点向附近寻找一点 B ,使得| AB| 在这些点中最短 ,可 以证明| AB| 就是我们要求的直径。这样 ,我们就可以求 出 2R = 60 ,R = 30 (像素) 。 2 血管中轴线的确定 在已获得球半径的基础上 ,我们对切片图样进行相 关操作来确定半径为 R 的圆的圆心的位置。 所谓相关操作是指利用两函数 f ( x , y) 、g ( x , y) 作相关积分时的特性 ,在样品图样中定位某种模版图样 的位置方法。我们以一简单情形为例简单的说明相关操 作。假设样品具有如下图形 : f ( x) = 6 i =1 g ( x - x i ) ,其中 g ( x ) 为一个单位方势垒函数 : g ( x) = 1 | x | < a/ 2 0 | x | < a/ 2 x i 是若 干个在 x 轴上的离散分布的固定点。为了求出在样品 f ( x) 中 ,和 g ( x) 相同的图样的位置。我 们可以对这两个函数进行相关运算 : g ( x) Ý f ( x) = + ∞ - ∞g ( x’- x ) f ( x’) dx’= 6 N i =1 + ∞ - ∞g ( x’- x ) g ( x’- x i ) dx’ (1) 作变量替换 : x”= x’- x i ,则 g ( x) Ý f ( x) = 6 N i = 1 + ∞ - ∞g ( x”+ x i - x ) g ( x”) dx” = 6 N i = 1 + ∞ - ∞g ( x”- ( x - x i ) ) g ( x) ”dx” = 6 N i = 1 c ( x - x i ) (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 3 ( u , v) F( u , v) ,其中公式中 Ζ 两边表示一对傅立叶变换 对 , G、F 分别表示 g ( x , y) 和 f ( x , y) 的傅立叶变换函数。因此 ,在具体计算上 ,我们可以先将 样品图样进行傅立叶变换求出它的频谱 ,再将其与模版图样的频谱的复共轭相乘 ,并将积作反 傅立叶变换 ,即可得到所求的两图样的相关。 我们将模版 g ( x , y) 取为大圆 ,用此模版在切片图样 f ( x , y) 中搜索相同的图样。由大圆 存在的唯一性可知 ,这样找到的唯一一个最大亮点就应是大圆的圆心 ,也即中轴线在此平面上 的坐标 ,这样问题也就得到了解决。我们利用 MA TLAB 提供的傅立叶变换和反傅立叶变换的 03 工 程 数 学 学 报 第 19 卷 © 1995-2004 Tsinghua Tongfang Optical Disc Co., Ltd. All rights reserved