三维血管重建优秀论文_含源代码

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

2、屮的 最大值矩阵，最大值即为血管半径。最后计算所有切片的血管半径，并对这些半 径求平均值，得到平均血管半径为：29.6799u mo由100张切片的最大内切圆 圆心坐标拟合得出中轴线方程以及其在三个坐标平面内的投影曲线方程。由中轴线得到血管的三维立体重建图,用平面z = i,（i = 0,24,49,74）去截血管 的三维立体重建图，得到新的4张截面图。把它们分别与题设中的对应截面进行 内点个数对比。我们泄义两张切片所共同拥有的内点个数与原切片内点个数的比 值为重合度。计算得到平均重合度为：9&19%。关键词：血管半径屮轴线切片重建（來自作者：欢迎各界人士批评指正，学术交流邮箱nibz

3、。文章 作于2011年8月10日，陕西科技大学理学院实验室）1问题的重述断面可用于了解生物组织、器官等的形态。例如，将样本染色后切成厚约 的切片，在显微镜下观察该横断面的组织形态结构。如果用切片机连续不 断地将样本切成数十、成百的平行切片，可依次逐片观察。根据拍照并采样得 到的平行切片数字图象，运用计算机可重建组织、器官等准确的三维形态。假设某些血管可视为一类特殊的管道，该管道的表面是曲球心沿着某一曲 线(称为中轴线)的球滚动包络而成。例如圆柱就是这样一种管道，其中轴线为 直线，由半径固定的球滚动包络形成。现有某管道的相继100张平行切片图像，记录了管道与切片的交。图像文件 名依次为0.bmp

4、、1.bmp、99.bmp,格式均为bmp,宽、高均为512个 象素(pixel) o为简化骚见，假设：管道中轴线与每张切片有且只有一个交点; 球半径固定；切片间距以及图像象素的尺寸均为1。计算管道的中轴线与半径，给出具体的算法，并绘制中轴线在xy、yz、zx 平面的投影图。2模型假设(1) 假设血管中轴线与每张切片有且只有一个交点。(2) 假设血管的半径固定不变。(3) 假设在对切片拍照的过程中不存在误差。(4) 假设血管不存在严重扭曲。3符号说明第£张切片中轮廓线的最大内切圆圆心第r张切片屮轮廓线的最大内切圆半径d.jk第£张切片中第，个内点与第丿个轮廓点的距离“512

5、x512x100100张切片图片转换以后的三维01矩阵刖5陀站00100张切片图片的轮廓线生成的矩阵c"®00x3100张切片轮廓的最大内切圆圆心坐标rf,平均血管半径第£张切片图片中所有内点的集合第k张切片图片屮轮廓点的集合str-square原切片图片的上内点及边界点的集合 重新切片得到的内点及边界点的集合 多项式拟合的指标数4问题的分析根据题忖整个管道是由球心沿着某一曲线（称为中轴线）的球滚动包络而成。 基于几何原理可以得出以下两条定理。（1）球的任意截面都是圆。（2）经过球心的球截面是所有的截面圆当中半径最大的圆。基于上述两个定理可以得出：每张切片的最大内

6、切圆的圆心位于血管的屮轴 线上，该圆的半径等于血管半径。分别求出100张切片轮廓线的最人内切圆圆心q与半径他。 对所有最大内切圆的圆心做拟合，即可求出血管的中轴线。100为减少谋差，我们用求均值的方法r=耳來得到平均血管半径r。 100根据拟合出的中轴线方程，即可求出中轴线在xy、yz、zx平而上的投影。5模型的建立与求解5.1模型建立5.1.1求血管的半径（1）导入数据，转换存储方式为了计算方便，利用matlab软件将100张bmp文件转换存储为一个三维矩阵pa5i2x512x1oo （0代表黑色像素点，1代表白色像素点）。（2）求切片轮廓线上各点坐标利用matlab软件的内部函数edge,

7、求得所有切片图片的轮廓线生成的矩阵pb512x512x100（3）求切片轮廓线的最大内切圆的半径首先求出每个内点距轮廓线的距离，取其屮的最小值，即为以这些内点为圆 心的轮廓线的内切圆半径；其次找出所有内点确定的最小内切圆中半径最大的一 个，即为轮廓线的最人内切圆半径。具体步骤：stepl在三维01矩阵pa5i2x5i2xi00'i'遍历搜索内点，将其存储于集合5）屮;step2在轮廓线矩阵pb512x5i2x1（）（）中遍历搜索轮廓点，将其存储于集合屮中； step3计算集合几）中第i个内点a（x,y）到集合屮屮第丿个轮廓点b（p,q）的距离 dijk = 7（x-p）2+g-

8、?）2，首先取每个内点到轮廓点的最小距离，即为以这些内点为圆心的内切圆 半径；其次找出所有内点确定的最小内切圆中半径最大的一个，即为轮 廓线的最大内切圆半径rk = maxmind做i j jstep4计算平均血管半径100斤=100（4）求切片轮廓线的最大内切圆心坐标将100张轮廓线的最大内切圆半径所对应的圆心坐标c,（z,w）记录于矩阵100.3）中 °een（u）=乙< cen（k 2） = w，伙=1,2,100）c叫.3）二 k,5. 1. 2由圆心坐标求中轴线方程及曲线投影（1）利用matlab画图工具由100组圆心坐标画出中轴线以及中轴线在xy、yz、zx面的投影

9、，由于图像文件以象素为单位存储数据，这使得图像的原始数据是离散化的。对于离散的位图來说，将变得不可分辨，可能切片轮廓线的最大 圆不止一个。为了减少误差，将中轴线在xy、yz、zx面的曲线方程进行多项式 拟合，建立拟合方程：y 話（x）x = “(z)（2）由建立的投影曲面得到中轴线的方程<r = f,(t)；z = f2(ry95.2模型求解根据所建模型求岀100张切片的最大内切圆的圆心坐标、半径以及血管半径表1最人内切圆半径及圆心处标切片序号最大内切圆半径（r/u m）最大内切圆圆心处标横坐标(x/ u m)纵坐标(y/ u m)竖坐标(z/u m)0290

10、69-16011229-16022329069-16024529155-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

11、-16015242529.833-1601625262961461461461461461461461473215355353629732152583738296145475475353

12、682682682698-115116474829.698-113118484929.698-112119495029.698111120505129.698-111120515229.698-111120525329.698-111120535429.682-111120545529.53-66150555629.53-59153565729.732-70148575829.967-70148585929.732-54155596029.732-54155606129.614-311626

13、16229.614-31162626329.6146166636429.833-6166646529.83381676566301016766673()11167676830121676869301216769703012167707129.83326166717229.61426166727329.61446163737429.73265158747529.73265158757629.73271156767729.68296145777829.68296145787929.698127124798029.698127124808129.732122128818229.73212212882

14、8329.732122128838429.698127124848529.732150101858629.732150101868729.73215496878829.698137115888929.698139113899029.68216088909129.68216285919229.54717659929329.61418340939429.73218049949529.61418340959629.61418340969729.61418533979830190098993()190199(1) 血管半径为斤=29.6799 um。(2) 中轴线在xy平面投影曲线的拟合方程(rsqu

15、are:0.9939)r = 6.23xlo',8x9 -1.58xloj4x8 +1.74xlo-,1x7 -1.09xl0 8x6 +4.26xlo'6x5 -0.00108x4 +0.17693x3 18.059x2 +1041.7x 25566(3) 屮轴线在yz平面投影曲线的拟合方程(r-square:1)z = -6.32x1 o'15 y9 +1.95x10-'1y8 -2.66x io 8 y7 + 2.11 x io 中轴线方程：x =6.23x10 ,8x9 -1.58x10',4x8 -i-lxlo x7 -1.09x108x6 +

16、4.26x10'6x5 -0.00108x4 +0.17693x3 -18.059x2 +1041.7x-25566,z = -6.32 x 1015 y9 +1.95 x 10 11y8 - 2.66 x 10 8 y7 + 2.11 x 1o*5 y6 - 0.01075y5 + 3.6293y4 - 813.88y3+ 17 x 105 y2-9.75 xlo6y +3.60x10%x =1.29x1013z9 -5.96x 1 o'11z8 +1.16x 108z7 -1.23x 1 o'6z6+7.56x1 o'5z5 -0.00269z4 + 0.0

17、53686z3-0.55402z2 + 2.4609z + 92.988, y6 -0.01075y5 + 3.6293y4 -813.88y3 + l17x105y2 -9.75x106y + 3.60x108(4) 屮轴线在zx平面投影曲线的拟合方程(r-square:1)x = 1.29xlo ,3z9 -5.96xlo-"z8 +1.16xlo'8z7 -1.23xlo'6z6 +7.56x105z5 -0.00269z4 + 0.053686z3-0.55402z2 + 2.4609z+ 92.988图1中轴线曲线图2中轴线拟合曲线图3屮轴线在xy平而的投影

18、图4屮轴在xy平面的投影拟合1111线图5中轴线在yz平面的投影图6中轴线在yz平面的投影拟合曲线x图7中轴线在zx平血的投影图8中轴线在zx平面的投影拟合曲线曲线模拟的血管250、20015010050、ox50、-100-150图9曲线模拟的血管6模型检验与分析用多项式拟合得到重建后血管的三维空i'可曲线，再用平而z = /,(/ = 1,25,50,75)去截空间曲线得到新的4层切片平而分别与题目中对应切片进行对比，在这里我们称sr>t的面积与宇的面积比值为切片重合度。第01、25、50、75张图片分别与对应的新截取图片的对比图图1第01张切片 图11第25张切片 图12第

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

20、型检验的数据可以看出计算的精度，但是整个运算在pentium(r) dual-core 处理器上耗费的时间己经超过120分钟。参考文献刘会灯，朱 飞.matlab编程基础与典型应用m 2008年3月附录% % program 1 (找出切片的边界线的图)pic=zeros(512,512,100); pic2=zeros(512,512,100);for i=0:99s=sprintf(,d:apics%d.bmp,i);pic(:,:,i+1)=imread(s);pic2(:,:,i+1)=edge(pic(:,:,i+1); imshow(pic2(:,:,i+1);endm,n,i+1

21、=find(pic2(:,:ji4-1);% program 2 (求半径，找中轴线与切片交点)ss=100;lens=1000*ones(512,512,ss);centres=zeros(ss,2);r=zeros(ss,1);for k=1:ssfor i=1:512for j=1:512if pic(i,j,k)=0for m=1:512for n=1:512if pic2(m,n,k)=1t1=sqrt(i-m)*(i-m)+(j-n)*(j-n);if iens(i,j,k)>t1iens(i,j,k)=t1;endendendendendendendfor i=1:512f

22、or j=1:512if lens(i,j,k)=1000lens(i,j,k)=-1;endendendr(k,1)=max(ma x(len s(:,:,k);forj=1:512for i=1:512if r(k,1)-lens(i,j,k)<0.0001cen tres(k,1)=i;cen tres(k,2)=j;endendendenddisp('this program has been done!*)% program 3 (多项式拟合空间屮轴线)format longcenn=zeros(100,3);x=cen(:,1);y=ce n(:,2);z=ce n(:,3);px=polyfit(乙 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 ;ce nn (:,3)=z;% % program 4(绘制血管三维图像)t=lin space(0,pi,25);p=li