血管三维重建01年A题 数学建模.doc

A题 血管的三维重建摘要对于血管的三维重建问题，关键是找出理论假设下血管的中轴线和血管的直径。通过这两个参数的确定就可以基本上绘制出血管图来。首先，对给出的100张血管断面的二进制图像进行取反操作，应用软件获得图片中血管图像内部点和边缘点的坐标，根据程序求得最大内切圆的半径及其圆心坐标。具体是（1）求出内部任意一点与边缘点的距离，取距离最小的值作为以该点为圆心的内切圆半径；（2）在这些半径中找到数值最大的值即为该图像的最大内切圆半径，对应的点坐标即为最大内切圆的圆心坐标。对所有求得的半径取平均值得到。然后，根据100个圆心坐标进行多项式曲线拟合，得到中轴线方程为： 剩余标准差为1.5518，说明此回归模型的显著性好。 绘制出曲线图，并投影到X-Y，X-Z，Y-Z坐标面上。关键词：边界提取最大内切圆法多项式拟合1 问题重述断面可用于了解生物组织、器官等的形态。例如，将样本染色后切成厚约的切片，在显微镜下观察该横断面的组织形态结构。如果用切片机连续不断地将样本切成数十、成百的平行切片， 可依次逐片观察。根据拍照并采样得到的平行切片数字图象，运用计算机可重建组织、器官等准确的三维形态。现有某管道的相继100张平行切片图像，记录了管道与切片的交。图像文件名依次为0.bmp、1.bmp、 99.bmp，格式均为BMP，宽、高均为512个象素（pixel）。试计算管道的中轴线与半径，给出具体的算法，并绘制中轴线在X-Z、Y-Z、X-Y平面的投影图。2 模型假设（1）血管可视为一类特殊的管道，该管道的表面是由球心沿着某一曲线（称为中轴线）的球滚动包络而成；（2）球半径固定；（3）管道中轴线与每张切片有且只有一个交点；（4）图片间距以及图像像素尺寸均为1。3 符号说明：切片数，从0到99。：切片图像最大内切圆圆心范围坐标。：切片图像边缘点的坐标。：切片图像的最大内切圆半径。：切片图像最大内切圆圆心的横坐标。：切片图像最大内切圆圆心的纵坐标。：从（-9，-9）点到（10,10）点矩形范围内所有横、纵坐标为整数的点。：切片图像中任意一点到的最短距离4 问题分析本题给出了某管道的相继100张平行切片图像，记录了管道与切片的交。由于该血管是由球心沿着某一曲线的球滚动包络而成，且球的半径固定。图像格式均为BMP，宽、高均为512个象素（pixel），图像象素的尺寸为1，对图像进行取反操作，一是函数求解需要，二是得到最大内切圆半径为图像内部球体半径更符合要求。所以我们得到以下两个结论：（1）每张图片图形的最大内切圆的半径即为滚动球体的半径；（2）圆心坐标即为滚动球体的球心坐标，由于切片间距也为1，故可以得到滚动球体的三维坐标。根据每张图片最大内切圆的圆心拟合出管道的中轴线，并绘制出来，将绘制出的中轴线分别投影到X-Y、Y-Z、Z-X平面即可得到中轴线在X-Y、Y-Z、Z-X平面的投影图。5模型建立与求解首先用画图板对100张图片取反色，然后导入到matlab中。（1） 求得切片图像边缘点坐标及第一张切片图像的内部点坐标；（2） 求中任意一点到的最短距离；（3） 在中取最大值即为第一张切片图的最大内切圆半径，对应的内的一点的坐标即为最大内切圆的圆心坐标（，）；（4） 求下一张图像的最大内切圆圆心范围；（5） 求中任意一点到的最短距离；（6） 在中取最大值即为第Z张切片图的最大内切圆半径，对应的内点的坐标即为最大内切圆的圆心坐标（，）；（7） 检验所求圆心是否在原来预计的圆心范围内。若不符合与则 令 返回到（5）（8） 利用上面得到的中心点坐标通过Matlab拟合出中轴线方程，并绘制出。求得的最大内切圆半径和圆心坐标见附录表一。通过对100个最大内切圆半径求平均值，得到。利用matlab多项式拟合函数polyfit程序五进行三维曲线的拟合，得到中轴线方程为：剩余标准差为1.5518，数值小，说明此回归模型的显著性好。中轴线三维图为：图一 中轴线三维图中轴线X-Y面上的投影：图二 中轴线X-Y面上的投影中轴线X-Z面上的投影：图三 中轴线X-Z面上的投影中轴线Y-Z面上的投影：图四 中轴线Y-Z面上的投影6 模型检验画出由图片信息绘制三维点图程序，根据求得中轴线和半径绘制三维血管图。实际血管切片图像边缘包络出的血管图在各个坐标面上的投影：图五 实际血管X-Y平面投影图六 实际血管X-Z平面投影图七 实际血管Y-Z平面投影通过球体沿着中轴线滚出的血管图在各个坐标面上的投影：图八 理论血管X-Y平面投影图九 理论血管X-Z平面投影图十 理论血管Y-Z平面投影将模拟出来的实际血管投影图和理论血管投影图进行重叠，重合的部分不显示，进行投影比较：实际血管和理论血管X-Y平面投影图比较：图十一 X-Y平面投影图比较实际血管和理论血管X-Z平面投影图比较：图十二 X-Z平面投影图比较实际血管和理论血管Y-Z平面投影图比较：图十三 Y-Z平面投影图比较选取Y-Z平面投影，用matlab求出图形面积。 图十四 实际Y-Z平面投影 图十五 X-Z平面投影图比较（经处理）算的实际血管投影面积为：实际与理论面积偏差：偏差比误差较小，说明计算得到的理论血管比较符合实际血管。7 模型评价此模型优点在于对原始图片进行取反操作能较为准确的得到血管的三维图，多使用matlab函数进行求解，方便快捷。缺点在于模型检验是通过整体面积比较。没有做出算出的血管的切片和原始图片进行比较，准确性不高。8 模型推广此模型可以由实际的图片切片得到三维图的关键数据信息已达到重建三维模型的目的。可广泛应用于有数学特征的模型的重建上。比如，可以去点等半径球体假设，球体半径按照规律可变，或者血管的法向切面为其他形体的实体的重建。参考文献1姜启源，谢金星，叶俊，数学模型（第三版），北京：高等教育出版社，2003.82阳明盛、熊西文、林建华MATLAB基础及数学软件，大连理工大学出版社，2003年8月第一版参考网站：1中国数学建模网，附录：表一 最大内切圆半径和圆心坐标序号（Z坐标）半径X坐标Y坐标029257416129257416229257416329257416429257416529257416628.8617257416729258416829258416929258416102925841611292594161229259416132925941614292604161529261416162926141617292624161829263416192926441620292654162129266416222926741623292684162429.01722774162529.01722764162629.01722864152729.06892864152829.06892924142929.06892924143029.12042974133129.15482924143229.12043054113329.12043054113429.12043084103529.12043144083629.12043174073729.12043174073829.15483224053929.15483224054029.15483224054129.15483453934229.15483483914329.15483513894429.41093483914529.41093483914629.41093483914729.69853713724829.69853733704929.69853743695029.69853753685129.69853753685229.69853753685329.69853753685429.41093763675529.20623893515629.20623893515729.41094033285829.41094043265929.15484093156029.15484093156129.12044182866229.12044192816329.15484202756429.15484202756529.15484202756629.12044192816729.06894212686829.01724212686929.01724222367029.12044192117129.12044161987229.12044131887329.12044101807429.12044101807529.15484081757629.15484051687729.15483881397829.15483881397929.20623821318029.20623761248129.41093761248229.68163761248329.68163751238429.68163631118529.68163621108629.68163621108729.41093621108829.2062347988929.2062347989029.1548333899129.1548310789229.1548310789329.1548310789429.1548310789529.1204283709629.1204278699729.1204278699829.1204278699929.120428871Matlab:程序一：求中心点范围程序：function a=g(xx,yy)for j=1:20for i=1:20a(j*20+i-20,1)=xx-10+i;a(j*20+i-20,2)=yy-10+j;endEnd程序二：求最大内切圆半径和圆心坐标程序：BW10=imread(0.bmp);BW20=bwmorph(BW10,remove);i0,j0=find(BW10=1);b(:,1)=j0;b(:,2)=512-i0;r xx yy=ff(a,b)程序三：计算最大半径程序：function r xx yy=ff(a,b)m=size(a,1);n=size(b,1);w=zeros(m,1)+40;h=0;for i=1:m for j=1:n D(i,j)=sqrt(a(i,1)-b(j,1)2+(a(i,2)-b(j,2)2); if D(i,j)w(i,1) w(i,1)=D(i,j); x(i,1)=a(i,1); y(i,2)=a(i,2); end endendr index=max(w);xx=a(index,1);yy=a(index,2);程序四：求每张图片中个点的坐标程序：BW10=imread(0.bmp);BW20=bwmorph(BW10,remove);p0,q0=find(BW20=1);i0,j0=find(BW10=1)程序五：曲线拟合程序：p=polyfit(z,x,6);q=polyfit(z,y,6);xx=polyval(p,z);yy=polyval(q,z)plot3(x,y,z,.r,xx,yy,z)程序六：连续读图程序：I=zeros(99,0); for i=0:99url=strcat(num2str(i); imageurl=strcat(url,.bmp); I=imread(imageurl); I=imresize(I,512 512,bilinear);figure;imshow(I);Imgstore(:,:,:,i)=I;end 程序七：轮廓提取程序：BW1=imread(0.bmp);imshow(BW1);BW2=bwmorph(BW1,remove);BW3=bwmorph(BW1,skel,inf);imshow(BW2)figure,imshow(BW3)p,q=find(BW2=1);plot(q,-p,.r)程序八：根据图片信息绘制三维点图程序：clear;for b=0:99 m1=imread(int2str(b),.bmp); m(:,:,b+1)=edge(m1,sobel);endfor k=0:99 for i=1:512 fo