1、36 血管的三维形态重建模型摘要:本文从几何的角度证明了管道的切面与中轴线的交点就是管道切面所包含的唯一的最大圆的圆心，并且这个最大圆的半径与动球大圆的半径相等，用 MATLAB 软件把图象分别读成 0-1 图象数字距阵，从中易得出管道切面象素的坐标距阵，为减少计算机的运算量，我们先用人工的方法分析切轴点的可能范围，再用 MATLAB 在可能范围内搜索，找出最大圆及其圆心（切轴点）坐标（各圆心的坐标参见附录二）.最后计算得到的管道半径结果为： （单位：象素）.根据计算机搜索所得的 100 个切轴点的坐标，用30rMAPLE 软件进行描点，描出了中轴线的空间形状图形，及中轴线分别在 xy、yz、
2、xz 平面上的投影（图形请参看 6.2 模型的结果）. 最后还模拟出了血管管道的空间图形（见附录三）.关键词：包络；切轴点；图象数字距阵；象素坐标距阵1 问题的提出断面可用于了解生物组织、器官的形态， 将样本染色后切成平行切片依次在显微镜下观察，然后根据拍照采样得到平行切片数字图象，运用计算机可重建组织、器官的三维形态.现有某血管样本的相继 100 张平行切片的图象（均为 BMP 格式）.要求计算出管道的中轴线及管道半径，给出具体算法，并绘制中轴线在 xy、yz、zx 平面的投影图.2 问题的分析该管道的表面是由球心沿着中轴线滚动包络而成 ， 给出 100 张平行切片图象，每张与中轴线有且只有
3、一个交点（切轴点）.因此，如果知道每一个切面的切轴点坐标，则可用这 100 个切轴点坐标描出中轴线的空间形状，进而求出半径，因此首要解决的问题就是如何求出切轴点坐标.3 模型的假设与名词说明1） 所给血管样本经过染色等处理后认为是实心的.2） 血管的表面是有球心沿着中轴线的球滚动包络而成，且球半径固定，即血管的大小是均匀的.3） 中轴线与每张切片有且只有一个交点.4） 平行切片间的间距和图象象素的尺寸均为 1.5） 切轴点：即切面与管道中轴线的交点.第一期（2002 年 10 月） 韶关学院学生数学建模论文集 No.137 5 模型的建立血管管道的表面是由半径固定的、球心沿着中轴线的球滚动而形
4、成的包络面.每个切平面与中轴线有唯一的交点，即切轴点（动球球心）在切片中是唯一的.若能证明切片的最大圆即为动球的大圆，且最大圆唯一，则只要在切片中找出最大圆，那么最大圆的半径和圆心坐标即为动球半径和切轴点坐标.5.1 切轴点位置特性的分析及有关证明为了便于计算切轴点，我们先证明下面两个命题：命题 1：每张切片图象必能包含动球的一个大圆.证明 ：设中轴线与平面 的唯一交点为 , 则点)9,.210(,kZ ),(KYXOkk为动球球心的一个位置， 此时动球与平面 的交在管道所围实体与 的交kOkZz的内部 （ 由于动球在管道内 ） ， 所以每张切片图象必能包含动球的一个大圆.（证毕）命题 2：每
5、张切片图象所包含的圆中，最大的圆恰是动球的大圆，且最大的圆是唯一的.证明：令中轴线 的参数方程为： , )(90ttzytx其中 , 分别满足 .0t9 )(,90t则动球球面方程为：-(1)222 )()()( RtzxyxR 为动球半径，从而动球的包络面 S 的方程为：-(2) 0)()()()(222 tzttyttxt我们用平面 Z=k,k=0,1,2,.99 去截管道，得切片图象的边界曲线：-（3）kztkttyttxt R)()()()( 222 将此曲线投影到 xy 平面，即是切片中管道的边界曲线.38 有（3）式易知，当且仅当 时，边界曲线上的点 到中轴线与切kt)(kyx,面
6、交点 的距离达到最大值为 R，所以切点边界曲线所包含的圆的最大半径是kt),(R，这说明切点所含圆中最大的圆恰是动球的大圆.又由于切片 与中轴线 只有唯一的一个交点，故只有唯一,Z)9,.210(的 t 满足 因此切片包含的最大圆只有一个.（证毕）)(kt5.2 切轴点位置的算法分析根据命题 1、命题 2 的证明，下面给出各个切面上最大圆的圆心和半径的算法分析.题目没有给出切片图象数字，只给出 100 张切片的 BMP 图象，且切片图象不太规则，难以用一般的区域内最大圆的计算方法求解，故考虑用计算机搜索，先用 MALAB 软件的imread 图象读入命令分别把 100 张图象读成阶为 5125
7、12 的 0-1 距阵 （0 区代表切片kT中血管的图象区，由于距阵很大，本文不予列出） ，接着利用 MATLAB 的距阵生成函数，由 中所有 0 元素的坐标生成管道切面图象的象素坐标距阵 ：为减少计算量，先分析kT kQ切轴点的可能范围 ，然后用计算机在 范围内进行搜索，寻找最大圆的圆心及半径.kAkA6 模型的求解6.1 算法实现1） 读入图片数据，生成 512512 的 0-1 距阵 ；kT2） 产生管道切面图象象素的坐标距阵 ，把 中元素为 0 的下标加入 ；kQkQ3） 分析可能的切轴点坐标范围 ；A4） 计算机在 范围内搜寻最大圆的圆心和半径；kA记 为切面 上含有的象点的集合SZ
8、a. 选定一象点 作为估计的圆心，选取一正整数 L,把以 为中心,边长为 2L 的0P0P正方形内的象点作为搜索范围 S.b. 取一正整数 r（例如 r=29）作为开始检验的半径值.c. 从 开始,由内向外在 S 中寻找一点 P，使以 P 为圆心，r 为半径的圆所包含的0P象素点都在 中.kd. 当在 c 中找到一点 P 时， ，返回 c.1r第一期（2002 年 10 月） 韶关学院学生数学建模论文集 No.139 e. 如在 c 中找不到 P 时， ，返回 c 时.1rf. 如当 时在 c 中能找到一点 ，但 时不能找出 P 点，则取iri1ir, 作为 上的圆心， 作为中轴线上对应iR)
9、,(iyxkZkyxQik,的球心，对切面 的检验停止.g. 取 的平均值的取整值作为 R.9210,.,5） 重复步骤 4）得出 100 张切片的切轴点坐标.6） 用 MAPLE 的作图函数描出中轴线的空间图形及其在 XY、YZ、XZ 片面上的投影.具体算法的 MATLAB 程序请看附录一.6.2 模型的结果利用上述的模型和方法计算得出如下结果：管道半径 r=30，各切轴点的坐标请参看附录二，血管管道的模拟图参见附录三，中轴线的空间形态及其分别在 YZ、XZ、XY 平面上的投影图如下：中轴线的空间形态中轴线在 XY 平面的投影 中轴线在 XZ 平面的投影中轴线在 YZ 平面的投影40 7 模
10、型的评价与推广本模型从切面与中轴线的交点的特性出发， 证明了交点即是切面图象的最大圆的圆心，并把图象读成数字距阵，经人为缩小范围后，利用计算机搜索求解.这样既节省了计算机的运行时间，又减少了计算误差.这种空间图形的重建方法是比较好的，在许多领域内也是可用或可参考的，如其它类似的组织、器官的三维重建，CT 扫描、医用诊断机的扫描图片分析等.参考文献：1程卫国 冯华等编著. MALAB5.3 应用指南. 北京：人民邮电出版社. 2000.32陈维桓编著. 微分几何初步. 北京：北京大学出版社 1997 10附录一（求切轴点坐标和最大圆半径的主程序）：clear;clc;file=input(输入文
11、件名:);tmp=imread(file,bmp);%tmp=imread(c:temp50.bmp,bmp);%b50=sparse(ones(size(tmp)-double(tmp);tmpl=ones(size(tmp)-double(tmp);ix,iy=find(b50);ind50=ix,iy;tmp_cen0=input(输入估计的圆心坐标),%tmp_cen0=124,350,%tmp_r=input(输入估计的半径,r=29),%tmp_r=29,%r=29%clear k;for k=0:15clear k1;for k1=tmp_cen0(1)-k:tmp_cen0(1
12、)+ktmp_cen1=k1,tmp_cen0(2)-k;ok1=1;clear i;clear j;for i=-tmp_r:tmp_rfor j=-tmp_r:tmp_rif i2+j2=0.8disp(center=),tmp_cen1,disp(半径 =),tmp_r,disp(转换成题目的坐标),disp(x=),tmp_cen1(2)-256,disp(y=),255-tmp_cen1(1),break;endendif ok1=1break;endclear k2;for k2=tmp_cen0(1)-k:tmp_cen0(1)+ktmp_cen2=k2,tmp_cen0(2)+
13、k;ok2=1;clear i,clear j;for i=-tmp_r :tmp_rfor j=-tmp_r :tmp_rif i2+j2=0.8disp(center=),tmp_cen2,disp(半径=),tmp_r,disp(转换成题目的坐标),disp(x=),tmp_cen2(2)-256,disp(y=),255-tmp_cen2(1),break;endendif ok2=1break;endclear k3;for k3=tmp_cen0(2)-k:tmp_cen0(2)+k42 emp_cen3=tmp_cen0(1)+k,k3;ok3=1;clear i,clear j
14、,for i=-tmp_r:tmprfor j=-tmp_r:tmp_rif i2+j2=0.8disp(center=),tmp_cen3,disp(半径=),tmp_r,disp(转换成题目的坐标),disp(x=),tmp_cen3(2)-256,disp(y=),255-tmp_cen3(1),break;endendif ok3=1break;endclear k4;for k4=tmp_cen0(2)-k:tmp_cen0(2)+kemp_cen4=tmp_cen0(1)+k,k4;ok4=1;clear i,clear j,for i=-tmp_r:tmprfor j=-tmp_
15、r:tmp_rif i2+j2=0.8disp(center=),tmp_cen4,disp(半径=),tmp_r,disp(转换成题目的坐标),第一期（2002 年 10 月） 韶关学院学生数学建模论文集 No.143 disp(x=),tmp_cen4(2)-256,disp(y=),255-tmp_cen4(1),break;endendif ok4=1break;endendspy(b50)附录二（各切面与中轴线交点的坐标）：图片 切轴点 图片 切轴点 图片 切轴点Z=0 (1,160) Z=34 (36,157) Z=68 (166,-21)Z=1 (1,159) Z=35 (41,
16、156) Z=69 (165,-21)Z=2 (1,159) Z=36 (52,153) Z=70 (166,-21)Z=3 (2,159) Z=37 (61,150) Z=71 (165,-31)Z=4 (2,159) Z=38 (66,148) Z=72 (163,-45)Z=5 (2,159) Z=39 (66,148) Z=73 (161,-55)Z=6 (2,159) Z=40 (66,148) Z=74 (157,-68)Z=7 (2,159) Z=41 (84,148) Z=75 (154,-77)Z=8 (2,159) Z=42 (84,139) Z=76 (154,-77)Z
17、=9 (2,159) Z=43 (84,139) Z=77 (149,-89)Z=10 (2,159) Z=44 (92,134) Z=78 (141,-104)Z=11 (3,159) Z=45 (92,134) Z=79 (138,-109)Z=12 (3,159) Z=46 (92,134) Z=80 (138,-109)Z=13 (3,159) Z=47 (117,113) Z=81 (120,-133)Z=14 (4,159) Z=48 (118,111) Z=82 (119,-134)Z=15 (8,159) Z=49 (120,110) Z=86 (107,-146)Z=16 (
18、8,159) Z=50 (115,115) Z=84 (106,-147)Z=17 (9,159) Z=51 (118,112) Z=85 (105,-148)Z=18 (10,160) Z=52 (112,118) Z=86 (105,-148)Z=19 (9,159) Z=53 (114,116) Z=87 (82,-165)Z=20 (8,159) Z=54 (120,110) Z=88 (77,-168)Z=21 (10,159) Z=55 (139,85) Z=89 (77,-168)Z=22 (14,159) Z=56 (139,85) Z=90 (77,-168)44 Z=23
19、(15,159) Z=57 (139,85) Z=91 (77,-168)Z=24 (21,159) Z=58 (141,82) Z=92 (54,-179)Z=25 (22,159) Z=59 (148,69) Z=93 (49,-181)Z=26 (21,159) Z=60 (153,58) Z=94 (45,-182)Z=27 (29,158) Z=61 (159,41) Z=95 (39,-184)Z=28 (29,158) Z=62 (156,50) Z=96 (27,-187)Z=29 (28,158) Z=63 (160,37) Z=97 (22,-188)Z=30 (29,158) Z=64 (161,33) Z=98 (7,-190)Z=31 (29,158) Z=65 (165,10) Z=99 (7,-190)Z=32 (30,157) Z=66 (166,-1)Z=33 (36,157) Z=67 (166,-1)附录三：（血管管道空间模拟图）
