第28卷第3期2003年9月广西大学学报(自然科学版)
Journal of Guangxi U niversity (N at Sci Ed ) V o l . 28,N o. 3 Sep t . , 2003
α文章编号:100127445(2003) 0320206203
血管三维图像重建的数学方法黄新民
(广西大学数学与信息科学学院, 广西南宁530004)
摘要:使用M A TLAB 读出血管切片图的数据, , , 实现了图形的三维重建.
关键词:血管; 参数方程; 三维重建
中图分类号:O 24211 :A
2001A 题:血管的三维重建, 要求学生从给定的100张平行切片图中, , 并以此计算管道的中轴线与半径, 给出具体算法, 并绘制中轴线在X Y ,
. Y Z , ZX 该题的具体计算已经在文[1～6]中有许多详细而细致的讨论. 但是所有这些文章都没有说明使用的计算机语言及相应的程序. 按题目所附参考材料看, 出题人建议的是使用C ++程序. 事实上许多软件都可以处理图形图像, 而且还比使用C ++编程更容易. 本文使用M A TLAB 比较简单地完成题目所要求的计算. 原竞赛题目中只要求给出图像的中轴线及其在三个坐标平面上的投影, 实际上并没有真正实现血管的三维图像重建, 本文将讨论实现图形三维重建的一个方法.
要将图形数据读入, 在M A TLAB 中只要使用命令i m read 即可实现. 程序虽然简单, 但数据量较大, 因此要通过几次计算才能将所有数据读完, 再用命令save 存储到硬盘上, 以便下一步计算时使用.
对于读出的每一张图, 我们需要读出血管截面数据及其边界数据, 以便下一步从中读出每张内切圆的圆心坐标与内切圆半径. M A TLAB 记录的图形数据是这样的:如果记录图形数据的是一个二维数组A (i , j ) , 则第i 行第j 列上的数值为一个0到255的整数(占用一个字节的存储空间) , 这个值对应于一个有三列元素的长为256的颜色向量表中的一行, 该行的三分量就分别对应于红、绿、兰颜色的颜色值. 本题中由于只有黑白两色, 所以数组中每个元素只有值为零及为1两个. 在程序中我们用值为1记血管截面上的点, 因此就用周围四个点的值和为4(因此每一个点的值都是1) 来找出血管截面的内点, 去掉内点后我们就得到了血管截面的边界点.
不过我们得到的还不能算是真正的血管截面的边界. 因为图形扫描成数据点时, 血管截面中的每一个点实际代表了一小片面积的中心. 我们不妨称这个边界为血管截面的内边界. 而元素值为0的点(也就是血管截面外的点) 的边界称为外边界. 将内外边界求平均表示真正的边界才更准确.
读出边界及所围的区域数据后, 我们从区域中每取出一个点, 计算其与边界的最小距离(这就是在这个点所作的与边界相切的圆半径) , 再从所有内切圆的半径中找出最大时的圆心坐标与半径, 就得到小球在这个血管截面处时的球心坐标与球半径. 求出所有各张图的球心坐标与球半径, 就得到了中轴线坐标及球半径(严格地说球半径应该是一个常数, 但是由于扫描与计算误差, 因此只能取所有半径的平均值为球心半径) .
此程序的编写也不复杂. 但实际使用时不难发现由于数据量大, 需要计算很长时间(当然我们将程序运行后就可以静候结果) , 为加快计算速度可以使用各种办法. 在本文中我们发现圆半径不小于29, α收稿日期:20021208; 修订日期:20030810
作者简介:黄新民(1945) , 男, 广西贺州人, 广西大学教授.
第3期黄新民:血管三维图像重建的数学方法207因此将内点与边界点相同x 坐标的边界点比较其y 坐标, 当两个点的坐标之差小于29时就将这个点从内点中排除. 对坐标y 相同的点亦作同样处理. 由于减法比计算距离的平方运算要快, 因此这样做可以加快计算过程.
算出中轴线坐标及球半径后, 只要使用简单的二维作图命令p lo t 就可以画出中轴线向三个坐标平面的投影图, 完成建模题的要求. 但是这样做出的图形并没有实现图形的真正三维重建. 血管的三维图像重建有多种方法, 最简单的是将每张图的轮廓线按其高度在三维空间中复原. 本文将介绍实现图像三维重建的另外一种方法. 这种方法与学生在大学中所学的高等数学及线性代数的知识密切相关, 而且重建的图像可以使用光照、颜色等加工, 使图形显得十分美丽, 值得介绍.
设已求出中轴线方程:x =x (z ) , y =y (z ) , 0≤z ≤99, 下面我们求半径为得的包络面方程. 设(X , Y , Z ) 是包络面上任意一点, z (X -x (z ) ) 2+(Y -y (z ) 2-z 22. (1)
上式的z 为一个参数, 当它从0运动到99. 1) 式两边对参数z 求导, 则得第二个方程
) (x ) ) y z ) (Y -y (z ) ) +(Z -z ) =0. (2)
将(2) 式代入, (1+x ′(z ) 2X x (z ) ) 2+2x ′(z ) y ′(z ) (X -x (z ) ) (Y -y (z ) ) +(1+y ′(z ) 2) (Y -y (z ) ) 2=r 2, (3) (3) 式左方是一个二次型, 容易求出其特征值分别是Κ(z ) 2+y ′(z ) 2, 其相应的特征向量1=1, Κ2=1+x ′
是
p 1=(z ) -y ′, p 2=(z ) x ′(z ) x ′(z ) y ′.
因此作变换
(z ) u +x ′(z ) v , Y -y (z ) =x ′(z ) u +y ′(z ) v , X -x (z ) =-y ′(4)
(5) 之后, (3) 式将化简为(z ) 2+y ′(z ) 2) (u 2+(1+x ′(z ) 2+y ′(z ) 2) v 2) =r 2, (x ′
因此, (u , v ) 可建立参数式
, v =u =222(x ′(z ) +y ′(z ) (z ) +y ′(z ) 2) (1+x ′(z ) 2+y ′(z ) 2) x ′
将上式代入(4) 式, 立即得到管道曲面的参数方程X =x (z ) ++, (z ) 2+y ′(z ) 2(x ′(z ) 2+y ′(z ) 2) (1+x ′(z ) 2+y ′(z ) 2) x ′
Y =y (z ) ++, 222(z ) +y ′(z ) (x ′(z ) +y ′(z ) 2) (1+x ′(z ) 2+y ′(z ) 2) x ′
Z =z -22
(z ) 2+y ′(z ) 21+x ′r sin (t ) .
其中0≤z ≤99, 0≤t ≤2Π, 在中轴线坐标及球半径r 求出来后, 可以通过数据拟合得到中轴线的参数方程(由于数据是近似的, 因此作数据拟合是合适的)
x =x (z ) , y =y (z ) , 0≤z ≤99.
然后就可以画出血管的三维图形了. 通过选择观察点及颜色、光照等可以作出十分精致的三维图像. 下面的程序即为作图程序(其中x (z ) , y (z ) 是分别拟合成六次多项式) . 图1是此程序运行后得到的
. M A TLAB 作出的血管三维图
load cendata %调入中轴线坐标数据, 假设存储在变量B 中, 共100行, 分别表示x , y , r 坐标x =B (:, 1) ;
y =B (:, 2) ;
100; r =sum (B (:, 3) ) %中轴线的x 坐标向量%中轴线的y 坐标向量%内切圆半径平均值
208
z =[0:99]’;
p 1=po lyfit (z , x , 6) ;
p 2=po lyfit (z , y , 6) ;
zz =[0:0. 4:99]’;
xx =po lyval (p 1, zz ) ;
yy =po lyval (p 2, zz ) ;
dx 1=po lyder (p 1) ;
dy 1=po lyder (p 2) ;
dx 2=po lydval (dx 1, zz ) ;
dy 2=po lydval (dy 1, zz ) ;
R 2=dx 2, ^2+dy 2. ^2;
r 1=sqrt (R 2) ;
r 2=r 1. 3sqrt (1+R 2) ;
t =linspace (0, 23p i , 251) ;
s =size (t ) ; 广西大学学报(自然科学版) %z 坐标向量%以z 为自变量将x 拟合成六次多项式%以z 为自变量将y 拟合成六次多项式%重新定义步长更小的高度坐标向量zz %用拟合多项求对应于zz 的坐标xx %用拟合多项求对应于zz 的坐标yy %对x 的多项式微分多项式%对y 的多项式微分多项式%求x’(z ) %求y’(z
) 第28卷 %X =xx 3ones (s ) +r . 1) ) +(2. r 2) 3sin (t ) ) ;
Y =yy 3() 3(() 3co s (t ) +(dy 2. r 2) 3sin (t ) ) ;
Z =zz 3ones (r (r 1. r 2) 3sin (t ) ;
m esh (X , Y , Z ) %三维网线作图
%三维表面作图surface (X , Y , Z )
axis equal
axis off
h idden off
view ([35, 88])
shading interp
co lo r m ap ho t
caxis ([-70, 450])
cam ligh t (200, 180)
ligh ting gouraud %按实际比例值作图%不画坐标轴, 使图形更逼真%消隐, 去掉曲面上的网线, 使图形更逼真%选择观察角%使作图时片与片之平滑过渡%色彩处理(不一定恰当, 请读者自己尝试) %重设颜色, 使颜色接近红色%指定光源位置, 这些值为尝试值, 不一定合适%设置照明方式图1 用M A TLAB 作出 的血管三维图
参考文献:
[1] 廖武斌, 邓俊晔, 王 丹. 管道切片的三维重建[J ]. 工程数学学报, 2002, 19(5) :22228. .
[2] 胡亦斌, 向 杰, 程 翔. 利用切片的二维空间相关操作实现血管的三维重建[J ]. 工程数学学报, 2002, 19(5) :292
34.
[3] 徐 晋, 刘雪峰, 柏容刚. 血管的三维重建[J ]. 工程数学学报, 2002, 19(5) :35240.
[4] 柳海东, 陈 璐, 江 浩. 管道切片的三维重建[J ]. 工程数学学报, 2002, 19(5) :41246.
[5] 丁峰平, 周立丰, 李孝朋. 血管管道的三维重建[J ]. 工程数学学报, 2002, 19(5) :47253.
[6] 汪国昭, 陈凌钧. 血管三维重建的问题[J ]. 工程数学学报, 2002, 19(5) :54258.
3D recon struction of blood vessel
HU AN G X in 2m in
(Co llege of M athem atics and Info r m ati on Science , Guangxi U niversity , N anning 530004, Ch ina )
Abstract :M A TLAB is u sed to ob tain data of b lood vessel slices . T he radiu s of the b lood vessel is deter m ined and the axes cu rve has calcu lated . T he p aram etric equati on is ob tained and the 3D i m age of b lood vessel su rface is recon structed .
(责任编辑 刘海涛) Key words :b lood vessel ; param etric equati on ; 3D recon structi on
第28卷第3期2003年9月广西大学学报(自然科学版)
Journal of Guangxi U niversity (N at Sci Ed ) V o l . 28,N o. 3 Sep t . , 2003
α文章编号:100127445(2003) 0320206203
血管三维图像重建的数学方法黄新民
(广西大学数学与信息科学学院, 广西南宁530004)
摘要:使用M A TLAB 读出血管切片图的数据, , , 实现了图形的三维重建.
关键词:血管; 参数方程; 三维重建
中图分类号:O 24211 :A
2001A 题:血管的三维重建, 要求学生从给定的100张平行切片图中, , 并以此计算管道的中轴线与半径, 给出具体算法, 并绘制中轴线在X Y ,
. Y Z , ZX 该题的具体计算已经在文[1～6]中有许多详细而细致的讨论. 但是所有这些文章都没有说明使用的计算机语言及相应的程序. 按题目所附参考材料看, 出题人建议的是使用C ++程序. 事实上许多软件都可以处理图形图像, 而且还比使用C ++编程更容易. 本文使用M A TLAB 比较简单地完成题目所要求的计算. 原竞赛题目中只要求给出图像的中轴线及其在三个坐标平面上的投影, 实际上并没有真正实现血管的三维图像重建, 本文将讨论实现图形三维重建的一个方法.
要将图形数据读入, 在M A TLAB 中只要使用命令i m read 即可实现. 程序虽然简单, 但数据量较大, 因此要通过几次计算才能将所有数据读完, 再用命令save 存储到硬盘上, 以便下一步计算时使用.
对于读出的每一张图, 我们需要读出血管截面数据及其边界数据, 以便下一步从中读出每张内切圆的圆心坐标与内切圆半径. M A TLAB 记录的图形数据是这样的:如果记录图形数据的是一个二维数组A (i , j ) , 则第i 行第j 列上的数值为一个0到255的整数(占用一个字节的存储空间) , 这个值对应于一个有三列元素的长为256的颜色向量表中的一行, 该行的三分量就分别对应于红、绿、兰颜色的颜色值. 本题中由于只有黑白两色, 所以数组中每个元素只有值为零及为1两个. 在程序中我们用值为1记血管截面上的点, 因此就用周围四个点的值和为4(因此每一个点的值都是1) 来找出血管截面的内点, 去掉内点后我们就得到了血管截面的边界点.
不过我们得到的还不能算是真正的血管截面的边界. 因为图形扫描成数据点时, 血管截面中的每一个点实际代表了一小片面积的中心. 我们不妨称这个边界为血管截面的内边界. 而元素值为0的点(也就是血管截面外的点) 的边界称为外边界. 将内外边界求平均表示真正的边界才更准确.
读出边界及所围的区域数据后, 我们从区域中每取出一个点, 计算其与边界的最小距离(这就是在这个点所作的与边界相切的圆半径) , 再从所有内切圆的半径中找出最大时的圆心坐标与半径, 就得到小球在这个血管截面处时的球心坐标与球半径. 求出所有各张图的球心坐标与球半径, 就得到了中轴线坐标及球半径(严格地说球半径应该是一个常数, 但是由于扫描与计算误差, 因此只能取所有半径的平均值为球心半径) .
此程序的编写也不复杂. 但实际使用时不难发现由于数据量大, 需要计算很长时间(当然我们将程序运行后就可以静候结果) , 为加快计算速度可以使用各种办法. 在本文中我们发现圆半径不小于29, α收稿日期:20021208; 修订日期:20030810
作者简介:黄新民(1945) , 男, 广西贺州人, 广西大学教授.
第3期黄新民:血管三维图像重建的数学方法207因此将内点与边界点相同x 坐标的边界点比较其y 坐标, 当两个点的坐标之差小于29时就将这个点从内点中排除. 对坐标y 相同的点亦作同样处理. 由于减法比计算距离的平方运算要快, 因此这样做可以加快计算过程.
算出中轴线坐标及球半径后, 只要使用简单的二维作图命令p lo t 就可以画出中轴线向三个坐标平面的投影图, 完成建模题的要求. 但是这样做出的图形并没有实现图形的真正三维重建. 血管的三维图像重建有多种方法, 最简单的是将每张图的轮廓线按其高度在三维空间中复原. 本文将介绍实现图像三维重建的另外一种方法. 这种方法与学生在大学中所学的高等数学及线性代数的知识密切相关, 而且重建的图像可以使用光照、颜色等加工, 使图形显得十分美丽, 值得介绍.
设已求出中轴线方程:x =x (z ) , y =y (z ) , 0≤z ≤99, 下面我们求半径为得的包络面方程. 设(X , Y , Z ) 是包络面上任意一点, z (X -x (z ) ) 2+(Y -y (z ) 2-z 22. (1)
上式的z 为一个参数, 当它从0运动到99. 1) 式两边对参数z 求导, 则得第二个方程
) (x ) ) y z ) (Y -y (z ) ) +(Z -z ) =0. (2)
将(2) 式代入, (1+x ′(z ) 2X x (z ) ) 2+2x ′(z ) y ′(z ) (X -x (z ) ) (Y -y (z ) ) +(1+y ′(z ) 2) (Y -y (z ) ) 2=r 2, (3) (3) 式左方是一个二次型, 容易求出其特征值分别是Κ(z ) 2+y ′(z ) 2, 其相应的特征向量1=1, Κ2=1+x ′
是
p 1=(z ) -y ′, p 2=(z ) x ′(z ) x ′(z ) y ′.
因此作变换
(z ) u +x ′(z ) v , Y -y (z ) =x ′(z ) u +y ′(z ) v , X -x (z ) =-y ′(4)
(5) 之后, (3) 式将化简为(z ) 2+y ′(z ) 2) (u 2+(1+x ′(z ) 2+y ′(z ) 2) v 2) =r 2, (x ′
因此, (u , v ) 可建立参数式
, v =u =222(x ′(z ) +y ′(z ) (z ) +y ′(z ) 2) (1+x ′(z ) 2+y ′(z ) 2) x ′
将上式代入(4) 式, 立即得到管道曲面的参数方程X =x (z ) ++, (z ) 2+y ′(z ) 2(x ′(z ) 2+y ′(z ) 2) (1+x ′(z ) 2+y ′(z ) 2) x ′
Y =y (z ) ++, 222(z ) +y ′(z ) (x ′(z ) +y ′(z ) 2) (1+x ′(z ) 2+y ′(z ) 2) x ′
Z =z -22
(z ) 2+y ′(z ) 21+x ′r sin (t ) .
其中0≤z ≤99, 0≤t ≤2Π, 在中轴线坐标及球半径r 求出来后, 可以通过数据拟合得到中轴线的参数方程(由于数据是近似的, 因此作数据拟合是合适的)
x =x (z ) , y =y (z ) , 0≤z ≤99.
然后就可以画出血管的三维图形了. 通过选择观察点及颜色、光照等可以作出十分精致的三维图像. 下面的程序即为作图程序(其中x (z ) , y (z ) 是分别拟合成六次多项式) . 图1是此程序运行后得到的
. M A TLAB 作出的血管三维图
load cendata %调入中轴线坐标数据, 假设存储在变量B 中, 共100行, 分别表示x , y , r 坐标x =B (:, 1) ;
y =B (:, 2) ;
100; r =sum (B (:, 3) ) %中轴线的x 坐标向量%中轴线的y 坐标向量%内切圆半径平均值
208
z =[0:99]’;
p 1=po lyfit (z , x , 6) ;
p 2=po lyfit (z , y , 6) ;
zz =[0:0. 4:99]’;
xx =po lyval (p 1, zz ) ;
yy =po lyval (p 2, zz ) ;
dx 1=po lyder (p 1) ;
dy 1=po lyder (p 2) ;
dx 2=po lydval (dx 1, zz ) ;
dy 2=po lydval (dy 1, zz ) ;
R 2=dx 2, ^2+dy 2. ^2;
r 1=sqrt (R 2) ;
r 2=r 1. 3sqrt (1+R 2) ;
t =linspace (0, 23p i , 251) ;
s =size (t ) ; 广西大学学报(自然科学版) %z 坐标向量%以z 为自变量将x 拟合成六次多项式%以z 为自变量将y 拟合成六次多项式%重新定义步长更小的高度坐标向量zz %用拟合多项求对应于zz 的坐标xx %用拟合多项求对应于zz 的坐标yy %对x 的多项式微分多项式%对y 的多项式微分多项式%求x’(z ) %求y’(z
) 第28卷 %X =xx 3ones (s ) +r . 1) ) +(2. r 2) 3sin (t ) ) ;
Y =yy 3() 3(() 3co s (t ) +(dy 2. r 2) 3sin (t ) ) ;
Z =zz 3ones (r (r 1. r 2) 3sin (t ) ;
m esh (X , Y , Z ) %三维网线作图
%三维表面作图surface (X , Y , Z )
axis equal
axis off
h idden off
view ([35, 88])
shading interp
co lo r m ap ho t
caxis ([-70, 450])
cam ligh t (200, 180)
ligh ting gouraud %按实际比例值作图%不画坐标轴, 使图形更逼真%消隐, 去掉曲面上的网线, 使图形更逼真%选择观察角%使作图时片与片之平滑过渡%色彩处理(不一定恰当, 请读者自己尝试) %重设颜色, 使颜色接近红色%指定光源位置, 这些值为尝试值, 不一定合适%设置照明方式图1 用M A TLAB 作出 的血管三维图
参考文献:
[1] 廖武斌, 邓俊晔, 王 丹. 管道切片的三维重建[J ]. 工程数学学报, 2002, 19(5) :22228. .
[2] 胡亦斌, 向 杰, 程 翔. 利用切片的二维空间相关操作实现血管的三维重建[J ]. 工程数学学报, 2002, 19(5) :292
34.
[3] 徐 晋, 刘雪峰, 柏容刚. 血管的三维重建[J ]. 工程数学学报, 2002, 19(5) :35240.
[4] 柳海东, 陈 璐, 江 浩. 管道切片的三维重建[J ]. 工程数学学报, 2002, 19(5) :41246.
[5] 丁峰平, 周立丰, 李孝朋. 血管管道的三维重建[J ]. 工程数学学报, 2002, 19(5) :47253.
[6] 汪国昭, 陈凌钧. 血管三维重建的问题[J ]. 工程数学学报, 2002, 19(5) :54258.
3D recon struction of blood vessel
HU AN G X in 2m in
(Co llege of M athem atics and Info r m ati on Science , Guangxi U niversity , N anning 530004, Ch ina )
Abstract :M A TLAB is u sed to ob tain data of b lood vessel slices . T he radiu s of the b lood vessel is deter m ined and the axes cu rve has calcu lated . T he p aram etric equati on is ob tained and the 3D i m age of b lood vessel su rface is recon structed .
(责任编辑 刘海涛) Key words :b lood vessel ; param etric equati on ; 3D recon structi on