1、1 断面可用于了解生物组织、器官等的形态。例如，将断面可用于了解生物组织、器官等的形态。例如，将样本染色后切成厚约样本染色后切成厚约1 1m mm m的切片，在显微镜下观察该横断的切片，在显微镜下观察该横断面的组织形态结构。如果用切片机连续不断地将样本切成面的组织形态结构。如果用切片机连续不断地将样本切成数十、成百的平行切片，数十、成百的平行切片，可依次逐片观察。根据拍照并可依次逐片观察。根据拍照并采样得到的平行切片数字图象，运用计算机可重建组织、采样得到的平行切片数字图象，运用计算机可重建组织、器官等准确的三维形态。器官等准确的三维形态。假设某些血管可视为一类特殊的管道，该管道的表面是假设某
2、些血管可视为一类特殊的管道，该管道的表面是由球心沿着某一曲线（称为中轴线）的球滚动包络而成。由球心沿着某一曲线（称为中轴线）的球滚动包络而成。例如圆柱就是这样一种管道，其中轴线为直线，由半径固例如圆柱就是这样一种管道，其中轴线为直线，由半径固定的球滚动包络形成。定的球滚动包络形成。现有某管道的相继现有某管道的相继100100张平行切片图象，记录了管道与张平行切片图象，记录了管道与切片的交。切片的交。2 为简化，假设：管道中轴线与每张切片有且只有一个交点；为简化，假设：管道中轴线与每张切片有且只有一个交点；球半径固定；切片间距以及图象象素的尺寸均为球半径固定；切片间距以及图象象素的尺寸均为1 1
3、。取坐标系的取坐标系的Z Z轴垂直于切片，第轴垂直于切片，第1 1张切片为平面张切片为平面Z=0Z=0，第第100100张切片为平面张切片为平面Z=99Z=99。Z=zZ=z切片图象中象素的坐标依它们在切片图象中象素的坐标依它们在文件中出现的前后次序为文件中出现的前后次序为（-256,-256,z）,（-256,-255,z）,（-256，255，z），），（-255,-256,z-255,-256,z）,（-255,-255,z-255,-255,z）,（-255-255，255255，z z）（255,-256,z 255,-256,z）,（255,-255,z 255,-255,z）,（
4、255255，255255，z z）。）。试计算管道的中轴线与半径，给出具体的算法，并绘制中轴试计算管道的中轴线与半径，给出具体的算法，并绘制中轴线在线在XYXY、YZYZ、ZXZX平面的投影图。平面的投影图。3Z=1Z=0Z=49Z=50Z=99Z=9841 1）血管的表面是由半径固定、球心沿着某一曲血管的表面是由半径固定、球心沿着某一曲线（称为中轴线）的球滚动而形成的包络面。线（称为中轴线）的球滚动而形成的包络面。2 2）中轴线上任两点处的法截面圆不相交。）中轴线上任两点处的法截面圆不相交。3 3）管道中轴线与每张切片平面有且只有一个交）管道中轴线与每张切片平面有且只有一个交点。点。5 二
5、值图象的象素矩阵为二值图象的象素矩阵为0-10-1矩阵，矩阵横矩阵，矩阵横纵坐标对应原图象的直角坐标系位置，纵坐标对应原图象的直角坐标系位置，0 0和和1 1分分别代表该位置上的黑或白象素。别代表该位置上的黑或白象素。00000110011100106 将管道看成是小球运动包络形成的，则管将管道看成是小球运动包络形成的，则管道的切片便是滚动的小球在平面道的切片便是滚动的小球在平面z z上的截圆形上的截圆形成的包络。成的包络。-30-25-20-15-10-5051015-20-1001020小 球 沿 螺 旋小 球 沿 螺 旋线 滚 动 时，线 滚 动 时，在 某 平 面 上在 某 平 面 上
6、的 截 圆 形 成的 截 圆 形 成包络。包络。78解决问题的关键在于发现解决问题的关键在于发现定理：在一条粗细均匀血管的任何横截面的图象内，其包定理：在一条粗细均匀血管的任何横截面的图象内，其包含的最大内切圆的圆心位于中轴线上，该圆的半径等于滚含的最大内切圆的圆心位于中轴线上，该圆的半径等于滚动球的半径。动球的半径。基于：基于：1 1）球的任意截面都是圆）球的任意截面都是圆2 2）经过球心的球截面是所有）经过球心的球截面是所有截圆当中半径最大的圆截圆当中半径最大的圆9求滚动球半径求滚动球半径r r的方法的方法1 1）平均法）平均法2 2）抽样法）抽样法3 3）极大似然法）极大似然法4 4）叠
7、加法）叠加法求出的半径在求出的半径在28.528.5 30.530.5之间，为什么？之间，为什么？10需考虑的细节：需考虑的细节：1 1）如何读入一个数字图象，即如何由）如何读入一个数字图象，即如何由bmpbmp图象文件图象文件（或其它格式的图象文件）得到象素矩阵？或其它格式的图象文件）得到象素矩阵？MATLABMATLAB函数函数imreadimread()()：m=imread(0.bmp)m=imread(0.bmp)或把图象序列转换为三维数组或把图象序列转换为三维数组for b=0:99for b=0:99 m(:,:,b+1)=imread(int2str(b),.bmp);m(:,
8、:,b+1)=imread(int2str(b),.bmp);endend11需考虑的细节：需考虑的细节：2 2）何谓边界点？）何谓边界点？四邻域的概念四邻域的概念 找边界点坐标的算法找边界点坐标的算法也可用也可用MATLABMATLAB函数函数:edge()edge()如：如：m=imread(0.bmp)m=imread(0.bmp)；bw=edge(m,sobel)bw=edge(m,sobel)3 3）最大内切圆的圆心和半径须取小数和整数有区别）最大内切圆的圆心和半径须取小数和整数有区别吗？吗？4 4）取整方法：四舍五入（）取整方法：四舍五入（round)round)，向上取整向上取整
9、(ceil)ceil)，向下取整向下取整(floor)floor)125 5）如何由切片图叠合作出血管的立体图？）如何由切片图叠合作出血管的立体图？方法方法1 1：利用利用plot3()plot3()clearclearfor b=0:99for b=0:99 m1=imread(int2str(b),.bmp);m1=imread(int2str(b),.bmp);m(:,:,b+1)=edge(m1,sobel);m(:,:,b+1)=edge(m1,sobel);endendfor k=0:99for k=0:99 for i=1:512 for i=1:512 for j=1:512
10、for j=1:512 if(m(i,j,k+1)=1)if(m(i,j,k+1)=1)plot3(i,j,k+1,b-.);hold on plot3(i,j,k+1,b-.);hold onend,end,end,endend,end,end,endrotate3drotate3d，hold offhold off135 5）如何由切片图叠合作出血管的立体图？）如何由切片图叠合作出血管的立体图？方法方法1 1：利用利用plot3()plot3()14方法方法2 2：利用利用patch()patch()for b=0:99m(:,:,b+1)=imread(int2str(b),.bmp,b
11、mp);endm1=m(1:4:512,1:4:512,:);for i=1:size(m1,1)for j=1:size(m1,2)for k=1:size(m1,3)if m1(i,j,k)=0 m1(i,j,k)=88;else m1(i,j,k)=0;end,end,end,endms=smooth3(m1);hiso=patch(isosurface(ms,5),FaceColor,1,.75,.65,EdgeColor,none);hcap=patch(isocaps(m1,5),FaceColor,interp,EdgeColor,none);colormap(default)v
12、iew(45,30),axis tightgriddaspect(1,1,.4)lightangle(45,30);lighting phongisonormals(ms,hiso)set(hcap,AmbientStrength,.6)set(hiso,SpecularColorReflectance,0,SpecularExponent,50)rotate3d155 5）如何由切片图叠合作出血管的立体图？）如何由切片图叠合作出血管的立体图？方法方法2 2：利用利用patch()patch()16 中轴线的建模：求中轴线中轴线的建模：求中轴线 与各横断面的交点与各横断面的交点和曲线拟合、逼近
13、。参赛者使用的方法主要有：和曲线拟合、逼近。参赛者使用的方法主要有：1 1）枚举法）枚举法2 2）平行切线法）平行切线法3 3）外推法）外推法4 4）滚球法）滚球法5 5）投影法）投影法6 6）变换法）变换法7 7）细化法）细化法17方法一：以黑色点为圆心，方法一：以黑色点为圆心，以以r0r0为半径画圆，判断此圆为半径画圆，判断此圆周 上 的 象 素 点 是 否 为 黑周 上 的 象 素 点 是 否 为 黑色若全为黑色，则增大半色若全为黑色，则增大半径再判断，直到圆周上有白径再判断，直到圆周上有白色象素点，记录该圆的圆心色象素点，记录该圆的圆心坐标和半径，然后，取下一坐标和半径，然后，取下一个
14、黑色象素点重复以上步骤，个黑色象素点重复以上步骤，直至黑色区域中所有象素点直至黑色区域中所有象素点都搜索完为止；所有记录点都搜索完为止；所有记录点中对应于最大半径的圆心坐中对应于最大半径的圆心坐标，就是该切片的最大内切标，就是该切片的最大内切圆圆心坐标，此半径即是血圆圆心坐标，此半径即是血管的半径管的半径 -200-1000100200-1000100200020406080100 xyz18方法一得出的中轴线投影：方法一得出的中轴线投影：0204060801001201401601800102030405060708090100y z 0102030405060708090100-200-1
15、50-100-50050100150200z x 中轴线在中轴线在ZXZX平面投影图平面投影图 中轴线在中轴线在YZYZ平面投影图平面投影图什么原因造成内切圆圆心的误差如此之大？什么原因造成内切圆圆心的误差如此之大？19 模型的检验模型的检验 不检验只能说完成问题的一半。由于所给图象只不检验只能说完成问题的一半。由于所给图象只是离散数据，以象素为单位。无论以何种方式建模，是离散数据，以象素为单位。无论以何种方式建模，都是近似计算，效果如何，检验很必要。还可通过检都是近似计算，效果如何，检验很必要。还可通过检验，发现模型的误差，修正模型，提高模型的正确性。验，发现模型的误差，修正模型，提高模型的
16、正确性。阅卷中发现的检验方法主要有：阅卷中发现的检验方法主要有：1 1）逐片比较）逐片比较2 2）法平面法）法平面法3 3）滚动法）滚动法20 模型检验时也发现内切圆圆心误差较大，特别是模型检验时也发现内切圆圆心误差较大，特别是后面的切片。后面的切片。Z=30Z=021（-255,-256,z）,（-255,-255,z）,（-255，255，z）在第z层，求出其所有最大内切圆的圆心，若是轴心坐标，则以它为中心的球体被上下i(i=129)张切片所截的小圆应在第Zi张切片中，否则，该圆心就不是轴心。m(:,:,b+1)=imread(int2str(b),.daspect(1,1,.断面可用于了
17、解生物组织、器官等的形态。rotate3d，hold off方法1：利用plot3()也可用MATLAB函数:edge()m1(i,j,k)=0;方法一：以黑色点为圆心，以r0为半径画圆，判断此圆周上的象素点是否为黑色若全为黑色，则增大半径再判断，直到圆周上有白色象素点，记录该圆的圆心坐标和半径，然后，取下一个黑色象素点重复以上步骤，直至黑色区域中所有象素点都搜索完为止；模型检验时也发现内切圆圆心误差较大，特别是后面的切片。误差大的原因和改进途径（-256,-256,z）,（-256,-255,z）,（-256，255，z），view(45,30),axis tight求出同一张切片上的所有最
18、大内切圆的圆心，然后求平均值。MATLAB函数imread()：m=imread(0.为简化，假设：管道中轴线与每张切片有且只有一个交点；例如，将样本染色后切成厚约1m m的切片，在显微镜下观察该横断面的组织形态结构。模型检验时也发现内切圆圆心误差较大，特别是模型检验时也发现内切圆圆心误差较大，特别是后面的切片。后面的切片。Z=50Z=4022 模型检验时也发现内切圆圆心误差较大，特别是模型检验时也发现内切圆圆心误差较大，特别是后面的切片。后面的切片。Z=60Z=57 23实际图象边界上的点是连续的，在转换成实际图象边界上的点是连续的，在转换成bmpbmp图象时，象素表示的图象边界是离散的，成
19、图象时，象素表示的图象边界是离散的，成锯齿状，与实际图象有误差（舍入误差）。锯齿状，与实际图象有误差（舍入误差）。2 2）同一张切片上的最大内切圆不唯一）同一张切片上的最大内切圆不唯一 1 1）方法一：取平均）方法一：取平均求出同一张切片上的所有最大内切圆的圆心，然求出同一张切片上的所有最大内切圆的圆心，然后求平均值。后求平均值。241 1）方法一：取平均）方法一：取平均求出同一张切片上的所有最大内切圆的圆心，然求出同一张切片上的所有最大内切圆的圆心，然后求平均值。后求平均值。XY平面上的投影平面上的投影YZ平面上的投影平面上的投影ZX平面上的投影平面上的投影25检验检验 Z=30Z=40Z=
20、50Z=60261021021021022),(),(1),(1NyNxNyNxyxfyxgNyxeNe27切片Z值32435765均方误差0.00190.00240.00480.0040 由计算出的均方误差情况可以看出，均方由计算出的均方误差情况可以看出，均方误差的平均值低于误差的平均值低于0.5%0.5%。还可以用信噪比，相异象素点所占的比例还可以用信噪比，相异象素点所占的比例等其它性能指标来度量重构图象的质量。等其它性能指标来度量重构图象的质量。282 2）方法二：筛选）方法二：筛选+尖点特性尖点特性 在第在第z z层，求出其所有最大内切圆的圆心，若是层，求出其所有最大内切圆的圆心，若是
21、轴心坐标，则以它为中心的球体被上下轴心坐标，则以它为中心的球体被上下i(i=1i(i=1 29)29)张张切片所截的小圆应在第切片所截的小圆应在第Z Z i i张切片中，否则，该圆心张切片中，否则，该圆心就不是轴心。就不是轴心。第第0层切片层切片第第29层切片层切片292 2）方法二：筛选）方法二：筛选+尖端特性尖端特性 0.0.bmp bmp 7070.bmpbmp用上述方法处理效果较好，但用上述方法处理效果较好，但71.71.bmp bmp 9999.bmpbmp经筛选后还有不少候选点，可应用经筛选后还有不少候选点，可应用尖端特性进行二次筛选。尖端特性进行二次筛选。如如9999层，由于层，由于7070层切片左下角不可能是第层切片左下角不可能是第9898或或4242层球体在层球体在7070层截圆的包络，而是更小截圆的包络层截圆的包络，而是更小截圆的包络，该小圆只能来自，该小圆只能来自9999或或4141层，经计算层，经计算4141层不可能，层不可能，故此小圆定在第故此小圆定在第9999层的球体中。层的球体中。30
