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

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个象素（pixe

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

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

6、轴线方程，即可求出中轴线在xy、yz、zx平面上的投影。5模型的建立与求解5.1模型建立5.1.1求血管的半径（1）导入数据，转换存储方式为了计算方便，利用matlab1软件将100张bmp文件转换存储为一个三维0-1矩阵（0代表黑色像素点，1代表白色像素点）。（2）求切片轮廓线上各点坐标利用matlab软件的内部函数edge，求得所有切片图片的轮廓线生成的矩阵（3）求切片轮廓线的最大内切圆的半径首先求出每个内点距轮廓线的距离，取其中的最小值，即为以这些内点为圆心的轮廓线的内切圆半径；其次找出所有内点确定的最小内切圆中半径最大的一个，即为轮廓线的最大内切圆半径。具体步骤：step1在三维0-1

7、矩阵中遍历搜索内点,将其存储于集合中 ； step2 在轮廓线矩阵中遍历搜索轮廓点，将其存储于集合中； 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-1602452915

9、5-16026729.155-16027829.275-16028929.275-160391029.428-1603101129.428-1604111229.614-1604121329.614-1605131429.833-1606141529.833-1607151629.833-1608161730-1609171830-16010181930-16011192030-16012202130-16013212230-16014222330-16015232430-16015242529.833-160162526296146

106146146146146146147327327327326145475475353682682682/p>

11、.698-115116474829.698-113118484929.698-112119495029.698-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.8338167

12、6566301016766673011167676830121676869301216769703012167707129.83326166717229.61426166727329.61446163737429.73265158747529.73265158757629.73271156767729.68296145777829.68296145787929.698127124798029.698127124808129.732122128818229.732122128828329.732122128838429.698127124848529.732150101858629.732150

13、101868729.73215496878829.698137115888929.698139113899029.68216088909129.68216285919229.54717659929329.61418340939429.73218049949529.61418340959629.61418340969729.614185339798301900989930190199（1）血管半径为=29.6799m。（2）中轴线在xy平面投影曲线的拟合方程(r-square:0.9939)（3）中轴线在yz平面投影曲线的拟合方程(r-square:1)（4）中轴线在zx平面投影曲线的拟合方程(

14、r-square:1)（5）中轴线方程： 图1中轴线曲线 图2中轴线拟合曲线图3中轴线在xy平面的投影 图4中轴在xy平面的投影拟合曲线 图5中轴线在yz平面的投影 图6中轴线在yz平面的投影拟合曲线 图7中轴线在zx平面的投影 图8中轴线在zx平面的投影拟合曲线图9曲线模拟的血管6模型检验与分析用多项式拟合得到重建后血管的三维空间曲线，再用平面去截空间曲线得到新的4层切片平面分别与题目中对应切片进行对比，在这里我们称的面积与的面积比值为切片重合度。第01、25、50、75张图片分别与对应的新截取图片的对比图.图1 第01张切片 图11 第25张切片 图12 第50张切片 图13 第75张切片

15、图.图14截后第01张切片图15截后第25张切片图16截后第50张切片图17截后第75张切片根据上面曲线拟合得到的血管重新切片后，我们随机取出四张重切的图片，与原始的切片比较，得到两切片的重合度。计算得：第01张切片的重合度为99.20% ，第25张切片的重合度为98.65% ，第50张切片的重合度为97.30% ，第75张切片的重合度为97.61% ，计算平均重合度得到98.19%，说明本模型较为可靠。7模型的优点与缺点优点在于本模型的算法采用遍历搜索的方法得到最大内切圆的半径的精度较高，但是由于遍历搜索的运算次数异常庞大，导致模型的缺点是实践难度较大。从模型检验的数据可以看出计算的精度，但

16、是整个运算在pentium(r) dual-core处理器上耗费的时间已经超过120分钟。参考文献1刘会灯，朱 飞.matlab编程基础与典型应用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

17、=100;lens=1000*ones(512,512,ss);centres=zeros(ss,2);r=zeros(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=1:512 for i=1:512 if r(k,1)-lens(i,j,k)s1(i,j) s2(i,j)=s1(i,j); end end if s