三维血管重建建模论文含源代码

1、 血管的三维重建摘要本文探讨血管的三维重建，由血管的相继100平行切片图像计算血管的中轴线与半径，并绘制血管在三个坐标平面上的投影。由于血管的表面是由球心沿着某一曲线（称为中轴线）的球滚动包络而成，由此我们得出结论：每个切片一定包含滚动球的大圆，并且它一定为切片的最大切圆，而最大圆所对应的半径即为血管半径，所以求血管半径就转化为求每一个切片部的点到切片外部轮廓的最大半径。首先，读取100血管切面图，把它们转换成Logical矩阵，从中提取切片截面轮廓点构成一个新的矩阵。然后找到原图片矩阵中像素点的点（切片图片中轮廓线中的点），从而得到点到切片轮廓点的最小距离矩阵和最小距离中的最大值矩阵，最大值

2、即为血管半径。最后计算所有切片的血管半径，并对这些半径求平均值，得到平均血管半径为：29.6799m。由100切片的最大切圆圆心坐标拟合得出中轴线方程以与其在三个坐标平面的投影曲线方程。由中轴线得到血管的三维立体重建图，用平面去截血管的三维立体重建图，得到新的4截面图。把它们分别与题设中的对应截面进行点个数对比。我们定义两切片所共同拥有的点个数与原切片点个数的比值为重合度。计算得到平均重合度为：98.19% 。关键词：血管半径 中轴线 切片 重建（来自 欢迎各界人士批评指正，学术交流nibz .。文章作于2011年8月10日，科技大学理学院实验室）1问题的重述断面可用于了解生物组织、器官等的形

3、态。例如，将样本染色后切成厚约1m的切片，在显微镜下观察该横断面的组织形态结构。如果用切片机连续不断地将样本切成数十、成百的平行切片， 可依次逐片观察。根据拍照并采样得到的平行切片数字图象，运用计算机可重建组织、器官等准确的三维形态。假设某些血管可视为一类特殊的管道，该管道的表面是由球心沿着某一曲线（称为中轴线）的球滚动包络而成。例如圆柱就是这样一种管道，其中轴线为直线，由半径固定的球滚动包络形成。现有某管道的相继100平行切片图像，记录了管道与切片的交。图像文件名依次为0.bmp、1.bmp、 99.bmp，格式均为BMP，宽、高均为512个象素（pixel）。为简化起见，假设：管道中轴线与

4、每切片有且只有一个交点；球半径固定；切片间距以与图像象素的尺寸均为1。计算管道的中轴线与半径，给出具体的算法，并绘制中轴线在XY、YZ、ZX平面的投影图。 2模型假设 (1) 假设血管中轴线与每切片有且只有一个交点。(2) 假设血管的半径固定不变。(3) 假设在对切片拍照的过程中不存在误差。(4) 假设血管不存在严重扭曲。3符号说明第切片中轮廓线的最大切圆圆心第切片中轮廓线的最大切圆半径第切片中第个点与第个轮廓点的距离100切片图片转换以后的三维0-1矩阵100切片图片的轮廓线生成的矩阵100切片轮廓的最大切圆圆心坐标平均血管半径第切片图片中所有点的集合第切片图片中轮廓点的

5、集合原切片图片的上点与边界点的集合重新切片得到的点与边界点的集合R-square多项式拟合的指标数4问题的分析根据题目整个管道是由球心沿着某一曲线（称为中轴线）的球滚动包络而成。基于几何原理可以得出以下两条定理。（1）球的任意截面都是圆。（2）经过球心的球截面是所有的截面圆当中半径最大的圆。基于上述两个定理可以得出：每切片的最大切圆的圆心位于血管的中轴线上，该圆的半径等于血管半径。分别求出100切片轮廓线的最大切圆圆心与半径。对所有最大切圆的圆心做拟合，即可求出血管的中轴线。为减少误差，我们用求均值的方法来得到平均血管半径。根据拟合出的中轴线方程，即可求出中轴线在XY、YZ、ZX平面上的投影。

6、5模型的建立与求解5.1模型建立5.1.1求血管的半径（1）导入数据，转换存储方式为了计算方便，利用MATLAB1软件将100BMP文件转换存储为一个三维0-1矩阵（0代表黑色像素点，1代表白色像素点）。（2）求切片轮廓线上各点坐标利用MATLAB软件的部函数edge，求得所有切片图片的轮廓线生成的矩阵（3）求切片轮廓线的最大切圆的半径首先求出每个点距轮廓线的距离，取其中的最小值，即为以这些点为圆心的轮廓线的切圆半径；其次找出所有点确定的最小切圆中半径最大的一个，即为轮廓线的最大切圆半径。具体步骤：Step1在三维0-1矩阵中遍历搜索点,将其存储于集合中 ； Step2 在轮廓线矩阵中遍历搜索

7、轮廓点，将其存储于集合中； Step3计算集合中第个点到集合中第个轮廓点的距离 ，首先取每个点到轮廓点的最小距离，即为以这些点为圆心的切圆半径；其次找出所有点确定的最小切圆中半径最大的一个，即为轮 廓线的最大切圆半径Step4计算平均血管半径（4）求切片轮廓线的最大切圆心坐标将100轮廓线的最大切圆半径所对应的圆心坐标记录于矩阵中。5.1.2由圆心坐标求中轴线方程与曲线投影（1）利用MATLAB画图工具由100组圆心坐标画出中轴线以与中轴线在XY、YZ、ZX面的投影，由于图像文件以象素为单位存储数据，这使得图像的原始数据是离散化的。对于离散的位图来说，将变得不可分辨，可能切片轮廓线的最大圆不止

8、一个。为了减少误差，将中轴线在XY、YZ、ZX面的曲线方程进行多项式拟合，建立拟合方程：（2）由建立的投影曲面得到中轴线的方程5.2 模型求解根据所建模型求出100切片的最大切圆的圆心坐标、半径以与血管半径表1最大切圆半径与圆心坐标切片序号最大切圆半径（R/m）最大切圆圆心坐标横 坐标(X/m)纵 坐标(Y/m)竖 坐标（Z/m）029069-16011229-16022329069-16024529155-16026729.155-16027829.275-16028929.275-160391029.428

9、-1603101129.428-1604111229.614-1604121329.614-1605131429.833-1606141529.833-1607151629.833-1608161730-1609171830-16010181930-16011192030-16012202130-16013212230-16014222330-16015232430-16015242529.833-16016252629614614614614614-157

10、403132296146147327327327326145475475353682682682698-115116474829.698-113118484929.698-112119495029.69

11、8-111120505129.698-111120515229.698-111120525329.698-111120535429.682-111120545529.53-66150555629.53-59153565729.732-70148575829.967-70148585929.732-54155596029.732-54155606129.614-31162616229.614-31162626329.614-6166636429.833-6166646529.8338167656630101676667301116767683012167686930121676970301216

12、7707129.83326166717229.61426166727329.61446163737429.73265158747529.73265158757629.73271156767729.68296145777829.68296145787929.698127124798029.698127124808129.732122128818229.732122128828329.732122128838429.698127124848529.732150101858629.732150101868729.73215496878829.698137115888929.6981391138990

13、29.68216088909129.68216285919229.54717659929329.61418340939429.73218049949529.61418340959629.61418340969729.614185339798301900989930190199（1）血管半径为=29.6799m。（2）中轴线在XY平面投影曲线的拟合方程(R-square:0.9939)（3）中轴线在YZ平面投影曲线的拟合方程(R-square:1)（4）中轴线在ZX平面投影曲线的拟合方程(R-square:1)（5）中轴线方程： 图1中轴线曲线 图2中轴线拟合曲线图3中轴线在XY平面的投影 图4

14、中轴在XY平面的投影拟合曲线 图5中轴线在YZ平面的投影 图6中轴线在YZ平面的投影拟合曲线 图7中轴线在ZX平面的投影 图8中轴线在ZX平面的投影拟合曲线图9曲线模拟的血管6模型检验与分析用多项式拟合得到重建后血管的三维空间曲线，再用平面去截空间曲线得到新的4层切片平面分别与题目中对应切片进行对比，在这里我们称的面积与的面积比值为切片重合度。第01、25、50、75图片分别与对应的新截取图片的对比图.图1 第01切片 图11 第25切片 图12 第50切片 图13 第75切片图.图14截后第01切片图15截后第25切片图16截后第50切片图17截后第75切片根据上面曲线拟合得到的血管重新切片

15、后，我们随机取出四重切的图片，与原始的切片比较，得到两切片的重合度。计算得：第01切片的重合度为99.20% ，第25切片的重合度为98.65% ，第50切片的重合度为97.30% ，第75切片的重合度为97.61% ，计算平均重合度得到98.19%，说明本模型较为可靠。7模型的优点与缺点优点在于本模型的算法采用遍历搜索的方法得到最大切圆的半径的精度较高，但是由于遍历搜索的运算次数异常庞大，导致模型的缺点是实践难度较大。从模型检验的数据可以看出计算的精度，但是整个运算在Pentium(R) Dual-Core处理器上耗费的时间已经超过120分钟。参考文献1会灯，朱 飞.MATLAB编程基础与典

16、型应用M 2008年3月附录% program1（找出切片的边界线的图）pic=zeros(512,512,100);pic2=zeros(512,512,100);for i=0:99 s=sprintf('d:Apics%d.bmp',i); pic(:,:,i+1)=imread(s); pic2(:,:,i+1)=edge(pic(:,:,i+1); imshow(pic2(:,:,i+1);end%program 2（求半径，找中轴线与切片交点）ss=100;lens=1000*ones(512,512,ss);centres=zeros(ss,2);r=zeros(

17、ss,1);for k=1:ss for i=1:512 for j=1:512 if pic(i,j,k)=0 for m=1:512 for n=1:512 if pic2(m,n,k)=1 t1=sqrt(i-m)*(i-m)+(j-n)*(j-n); if lens(i,j,k)>t1 lens(i,j,k)=t1; end end end end end end end for i=1:512 for j=1:512 if lens(i,j,k)=1000 lens(i,j,k)=-1; end end end r(k,1)=max(max(lens(:,:,k); for j

18、=1:512 for i=1:512 if r(k,1)-lens(i,j,k)<0.0001 centres(k,1)=i; centres(k,2)=j; end end endenddisp('This program has been done!')%program 3（多项式拟合空间中轴线）format longcenn=zeros(100,3);x=cen(:,1);y=cen(:,2);z=cen(:,3);px=polyfit(z,x,9); x1=polyval(px,z);py=polyfit(z,y,9); y1=polyval(py,z);figure(1); plot3(x1,y1,z)cenn(:,1)=x1;cenn(:,2)=y1;cenn(:,3)=z;%program 4（绘制血管三维图像）t=linspace(0,pi,25);p=linspace(0,2*pi,25);theta,phi=meshgrid(t,p);for i=1:100 x=29.6799*sin(theta).*sin(phi)+cenn(i,1); y=29.6799*sin(theta).*cos(phi)+cenn(i,2); z=29.67