血管切片的三维重建.doc

血管切片的三维重建摘 要本文利用血管100张切片图，通过分析其几何特性,给出了确定其管道中轴线和半径的数学模型，并进行了血管的三维重建。对血管的三维重建问题，本文假定血管为等径管道，管道的表面是由球心沿着某一曲线（称为中轴线）的球滚动包络而成。首先读取100张血管切面图,把它们转换成二值矩阵,提取出边界矩阵、骨架矩阵，通过搜索每个切片截面最大内切圆,该内切圆圆心即为切片截面与管道中轴线的交点,内切圆半径即为管道半径,再通过对各个交点进行曲线拟合求出中轴线方程。利用多项式拟合进行了中轴线在三平面上投影的精确定位。本文较好的进行了三维血管的重建,得出所求平均半径为29.49,100张切片最大内切圆的圆心坐标见表1,通过建立空间三维曲线的参数方程，将Z作为参数，再用MATLAB中的polyfit函数，选取适当的拟合次数，使偏差平方和尽量小，但拟合多项式的最高次数不能太高，分别进行X，Z和Y，Z的多项式拟合，从而得到中轴线的参数方程，再对X,Y进行多项式拟合，最终得到中轴线及其在,和平面上的投影图、散点图。本文最大的亮点，对所建模型进行了很好的检验，在所求中轴线方程的基础上，求得血管空间曲面方程，令Z=0:1:99，对其进行切割，得到新的截痕，再对截痕内部进行填充，得到100张拟合三维血管管道新的平行切片图像。本文定义了重合度，通过计算新切片与原对应切片坐标相同点的个数所占百分比，计算重合度，最终所求最高重合度为80.25%，验证了模型的正确性。关键词：MATLAB图像处理；图像骨架；最小二乘曲线拟合；三维重建；重合度1.问题重述断面可用于了解生物组织、器官等的形态。例如将样本染色后切成厚约1的切片,在显微镜下观察该横断面的组织形态结构。如果用切片机连续不断地将样本切成数十、成百的平行切片,可依次逐片观察。根据拍照并采样得到的平行切片数字图象,运用计算机可重建组织、器官等准确的三维形态。假设某些血管可视为一类特殊的管道,该管道的表面是由球心沿着某一曲线(称为中轴线)的球滚动包络而成。例如圆柱就是这样一种管道,其中轴线为直线,由半径固定的球滚动包络形成。现有某管道的相继100张平行切片图象,记录了管道与切片的交。图象文件名依次为0.bmp、1.bmp、99.bmp,宽、高均为512个象素(pixel)。为简化起见,假设:管道中轴线与每张切片有且只有一个交点;球半径固定;切片间距以及图象象素的尺寸均为1。取坐标系的Z轴垂直于切片,第1张切片为平面Z=0,第100张切片为平面Z=99。试计算管道的中轴线与半径,给出具体的算法,并绘制中轴线在XY、YZ、ZX平面的投影图。2.模型假设1）假设样本血管可视为一类特殊的管道，该管道的表面是由球心沿着某一曲线的球滚动包络而成，球半径固定。2）医学上,血管不存在严重扭曲。3）假设管道中轴线与每张切片有且只有一个交点，即将切片视为无厚度的切平面。4）假设切片间距以及图像像素的尺寸均为1。5）假设管道中轴线处处连续，且充分光滑。6）假设两点间距离舍入服从“四舍五入”原则。7）中轴线上任两点处的法截面圆不相交。3.符号说明第i张切片的位图信息矩阵，为方便起见，有时也称其为第i个切平面；：管道中轴线参数方程，其中；：管半径，即沿中轴线滚动球的半径；：中心点，即管道中轴线与第i个切平面的交点，根据上述方程可知，其中；：的候选点的集合，为定义在某切平面上所有被判定为最大圆面的圆心；：，所有候选点的集合；4.模型建立及求解4.1.问题分析问题第一部分需要求出管道的中轴线方程和半径,第二部分需要绘制中轴线在各个平面上的投影。解决问题第一部分的关键在于发现以下定理：定理：在一条粗细均匀血管的任何横截面的图象内，其包含的最大内切圆的圆心位于中轴线上，该圆的半径等于滚动球的半径。基于：1）球的任意截面都是圆2）经过球心的球截面是所有截圆当中半径最大的圆【证明】：假设中轴线上存在另一点，以其为球心，以R为半径作球，则该球与此切平面相交成的圆面的半径不大于R，若为R，则可知与此切平面的距离为0，换句话说，在此切平面上。这与题设中轴线与每个切平面有且只有一个交点矛盾。因此可证明在该包络区域中可容纳的最大圆面是以为球心，以R为半径的圆。则可通过滚球法求解，即将管道近似看作是一个半径固定的球体滚动而成的，中轴线是球心滑过的曲线，是连续的。等距平行切割血管，中轴线与每张切片有且仅有一个交点，也就是每张切片上有且仅有一个球心，那么在每张切片上总可以找到且只能找到一个以球心为圆心，球半径为半径的圆，而且是此切片的最大内切圆，反过来也是成立的。因此，我们只需找到每张切片中的球心坐标就可以用多项式曲线拟合得到中轴线方程，通过寻找100个切平面的最大内切圆得到100个半径，再利用平均法求滚动球半径r，而中轴线在XY，YZ，ZX平面的投影图只需令Z=0，X=0，Y=0就可以得到。考虑到实际图象边界上的点是连续的,只有位置而没有大小,且点的位置以坐标来表示。在转换成BMP格式图象时,象素表示的图象边界是离散的,一般成锯齿状,图象范围与实际图象有误差,包括图象转换的系统误差,即点取舍引起的误差和计算的舍入误差。最大圆面不唯一是由于数字图像离散化的结果，【命题一】假设坐标及半径均为连续精确的，即使两个相近的数字也是可以分辨的，但实际情况并非如此理想。由于题设切片间距以及图像像素的尺寸均为1，而图像像素是图中可分辨的最小尺度。例如以为球心，为半径的球与第i+1个切平面相交的圆面半径为：，在数字图像中将被舍入为29，其结果是该圆面也可能被误判为最大圆面，这样所得实际为在上的投影。从而造成存在多个最大圆面。管半径的不唯一是由于题设管道中轴线与每张切片有且只有一个交点，而实际中切片总是有厚度的，再加上数字图像离散化的影响：对于管半径为29的血管，其切片图像产生1个单位的偏差是符合实际情况的。而且从计算结果也可以看出，出现最大半径为30的切平面的位置处在中轴线沿Z轴变化较平缓部分，这部分也是最容易产生偏差的。图像的骨架即图像中央的骨骼部分，是描述图像几何及拓扑性质的重要特征之一。求一图像骨架的过程通常称为对图像细化的过程，对被处理的图像进行细化有助于突出形状特点和减少冗余信息量，骨架是用1个点与1个点集的最小距离来定义的，可写成：其中：距离度量可以是欧氏的、城区的、或棋盘的。因最近距离取决于所用的距离量度，所以图像重建的结果也和所用的距离量度有关，由于中轴定义中保留有点到边界的最小距离，因此可以利用许多中心位于中轴线上的圆盘区域的并集覆盖与该中轴线对应的区域，该圆盘的中心为距离边界最远的点，这些圆盘称为最大圆，推广到三维则称为最大球，中轴即为最大球的球心，细化方法：即逐次去掉物体的边界点，同时必须保证区域的连通性不被破坏，最后保留的点组成物体的中轴，这种方法多采用模板方式来实现，抽取的中轴可以保留物体的拓扑结构。提取截面的骨架,对骨架上每一个点分别计算它到所有边界点的距离,取距离中最小的一个,这个距离就是以该点为圆心的、与轮廓线内切的若干个圆中最小的圆的半径。在这些由骨架上的点确定的最小内切圆中半径最大的一个就是轮廓线的最大内切圆。4.2.模型的建立4.2.1.问题第一部分：计算管道的中轴线与半径由于所给原始数据是512512单色BMP图像格式，便于观察，但不便于模型处理，因此首先将图像转化为512512的二维0-1矩阵形式，其0-1分布与图像对应，1表示黑色像素点，左上角为，为记录第i张切片的图象信息的0-1矩阵，为方便起见，本文中也用来表示第i张切平面。为的集合，记为根据题中建立的坐标系，某点的0-1值与矩阵M中元素的对应关系为：由假设，中轴线与每张切片有且只有一个交点，设该交点为，根据题中建立的坐标系可以设中轴线参数方称为：，其中由题设也可写为：，其中因此中轴线与第i张切片的交点。4.2.1.1.管道半径的算法步骤Step1：将中0-1进行互换；Step2：由“管道中轴线与每张切片有且只有一个交点”,知截面为单连通的区域,轮廓线为一闭合曲线。提取第i张图片轮廓（调用matlab中的edge函数，用二值矩阵来存储轮廓线）；Step3：提取第i张图片的骨架（调用matlab中bwmorph函数）；Step4：第i张图片的骨架上任一点到轮廓上每一点的最短距离，求这些最短距离中的最大值，即为第i张图片的对应最大内切圆的半径Ri，再找出这个最大值对应的骨架上的点的坐标即为圆心Ci；Step5：用for语句循环引用图片文件(strcat(f:/,int2str(i),.bmp),将上述step1-4循环可得到100张图片的圆心坐标和对应最大内切圆的半径；Step5:由确定管道半径R；Step6:将各圆心坐标用polyfit进行多项式拟合得到的曲线即为管道的中轴线。程序代码见附录2。4.2.1.2.中轴线的拟合的具体步骤Step1：由空间曲线的参数方程可设中轴线参数方程为：，其中Step2：已知圆心坐标，分别对x和z、y和z进行多项式拟合，调用polyfit函数返回了不同次数多项式的降幂排列的系数向量和偏差平方和的平方根，根据偏差平方和的平方根的大小确定拟合多项式的最高次数。由上述步骤得到了中轴线的参数方程。程序代码见附录3。4.2.2.问题第二部分：绘制中轴线在XY、YZ、ZX平面的投影图4.2.2.1.绘制中轴线步骤Step1：在得到中轴线参数方程的基础上，调用polyval函数可得到z=0：0.1：99对应的x、y坐标，即中轴线上点的空间坐标；Step2：使用MATLAB中的plot3函数，plot3（x，y，z）绘制中轴线；Step3：使用MATLAB中的hold on命令，绘制问题第一部分求得100个最大内切圆的圆心的散点图，与绘制的中轴线进行比较。程序代码见附录5。4.2.2.2.绘制中轴线在XY、YZ、ZX平面的投影图步骤Step1：调用MATLAB中的p,S = polyfit(x,y,n)函数，x，y分别为100个最大内切圆圆心的X,Y坐标，选取偏差平方和较小的最高拟合次数，绘制中轴线在XOY平面上的投影图，使用MATLAB中hold on命令，绘制散点图，与拟合投影图进行比较；Step2：同上，调用MATLAB中的p,S = polyfit(x,z,n)函数，x，z分别为100个最大内切圆圆心的X,Z坐标，选取偏差平方和较小的最高拟合次数，绘制中轴线在XOZ平面上的投影图，使用MATLAB中hold on命令，绘制散点图，与拟合投影图进行比较；Step3：调用MATLAB中的p,S = polyfit(z,y,n)函数，z，y分别为100个最大内切圆圆心的Z,Y坐标，选取偏差平方和较小的最高拟合次数，绘制中轴线在ZOY平面上的投影图，使用MATLAB中hold on命令，绘制散点图，与拟合投影图进行比较。程序代码见附录5。4.3.模型的求解4.3.1.管道半径及中轴线方程的求解对100个半径取平均值,得到管道半径:R=29.49。100张血管切片图对应的最大内切圆的园心坐标及半径见表1：表1 血管切片图对应球心坐标及半径表切片序号对应球半径对应球坐标横坐标(X)纵坐标(Y)竖坐标（Z）029.07-16100128.28-161-11229-16102329.07-16113429.07-16114529.07-16115629-16106729.02-16137829-16108928.86-161091028.86-1616101128.86-1617111228.86-1618121329.02-1619131429.02-16111141529.02-16112151629.02-16113161729.02-16115171829.02-16116181929.02-16117192029.02-16118202129.02-16119212229.02-16120222329.02-16121232429.07-16120242529.07-16120252629.07-16120262729.15-16029272829.27-16029282929.27-1602829302961616161737373736155685368-120111434429.7-119112444529.7-118113454629.7-118113464730.41-117114474830.41-116115484930.41-116115495030.41-115116505130.41-115116515230.41-114117525330.41-113118535430-112119545529.73-112119555629.7-112119565729.7-112119575829.53-82141585929.55-52155596029.55-52155606129.61-32161616229.61-32161626329.61-32161636429.61-36160646529.61-36160656629.43-27162666729.41-36160676829.27-27162686929.4345162697029.6145162707129.6145162717229.6145162727329.6164157737429.7367156747529.7364157757629.5580151767729.5380151777829.5380151787929.7134116798029.7135115808129.7135115818229.7136114828329.7137113838429.7137113848529.7138112858629.7138112868729.7138112878829.7139111888929.7139111899029.5317166909129.5317166919229.5317166929329.5317166939429.7318142949529.6118623959629.6118623969729.6118623979829.6118623989929.431871799平均半径29.49由上述计算得到100个圆心的空间坐标，调用polyfit函数返回了不同次数多项式的降幂排列的系数向量和偏差平方和的平方根（程序见附录五）,本文分别进行了x、z和y、z的1到10次的多项式拟合，返回偏差平方和的的平方根数值见下表：表2 x、z和y、z的1到10次的多项式拟合偏差平方和的的平方根数值最高次数x，z拟合偏差平方和的平方根y，z拟合偏差平方和的平方根1493.19447.942215.27299.393195.87105.74128.37105.44512491.9046119.5991.9047112.791.718112.4191.0719108.5390.96910105.4989.866通过计算观察比较，本文选取最高拟合次数遵循以下两个原则：1） 偏差平方和尽量小；2） 拟合最高次数不能太高。基于以上原则，在表二的数据中，分别选取x，z的7次拟合多项式，y，z的5次拟合多项式，得到中轴线的参数方程如下：向量px、py见下表：表3 向量px、py表pxpx(1)px(2)px(3)px(4)px(5)px(6)px(7)px(8)数值5.35E-10-1.75E-072.19E-05-0.001330.042474-0.65473.9331-165.83pypy(1)py(2)py(3)py(4)py(5)py(6)数值4.33E-07-0.0001090.008233-0.19062.2354-4.20754.3.2.绘制中轴线及其在XY、YZ、ZX平面的投影图4.3.2.1.绘制中轴线及其散点图图1 血管中轴线图4.3.2.2.绘制中轴线在X0Z面的投影图及其散点图图2 中轴线在X0Z面的投影图及其散点图4.3.2.3.绘制中轴线在Y0Z面的投影图及其散点图图3 中轴线在Y0Z面的投影图及其散点图4.3.2.4.绘制中轴线在X0Y面的投影图及其散点图在得到100个圆心的空间坐标，调用polyfit函数返回了不同次数多项式的降幂排列的系数向量和偏差平方和的平方根（程序见附录五）,进行了x、y的1到10次的多项式拟合，返回偏差平方和的的平方根数值见下表：表4 x、y的1到10次的多项式拟合偏差平方和的的平方根数值最高次数x，y拟合偏差平方和的平方根1525.262145.373134.25480.314577.453662.594757.034851.561949.761047.941同样遵循本问前面所述选取最高拟合多项式的2个原则，本文选取x，y的8次拟合多项式，绘制中轴线在XOY平面上的投影图如下：图4 中轴线在X0Y面的投影图及其散点图 4.3.2.5.利用100张血管切片图对血管进行三维重建用Matlab中的plot3命令，利用100张血管切片图对血管进行三维重建，得到如图5所示的血管三维图：图5 血管三维重建图5.误差分析及模型检验5.1.误差分析1）切片数字图像产生偏差实际图象边界上的点是连续的，在转换成bmp图象时，象素表示的图象边界是离散的，成锯齿状，与实际图象有误差（舍入误差）。其次，由显微镜的工作原理可知，一定厚度切片的显微图像是由其中所有组织的成像叠加而成的，而切片厚度是不可避免的，这将造成切片内不同层面图像的干扰，也就是说，每张切片图中的图像实际上是多个层面图像的叠加。然而，题中假设管道中轴线与每张切片有且只有一个交点，即不考虑切片厚度，因此以下叙述也将切片称为切平面，忽略其厚度。2）同一张切片上的最大内切圆不唯一。3）求各血管切片图对应的最大内切圆半径时存在舍入误差，以及由100张切片图对应圆的半径取平均得到球的半径时存在误差。5.2.模型检验5.2.1.检验思想在所求中轴线方程的基础上，求得血管空间曲面方程，令Z=0:1:99，对其进行切割，得到新的截痕，再对截痕内部进行填充，得到100张拟合三维血管管道新的平行切片图像。本文通过计算所得拟合图像与原图数据的对比，定义重合度(n:新切片与原对应切片坐标相同点的个数，即对应二值矩阵中0的重合个数；m：原图二值矩阵中0的个数)。5.2.2.检验具体步骤Step1:在得到中轴线参数方程的基础上，求得中轴线上任意一点的切线的方向向量；Step2:与该切线垂直的平面的方程（本文假定球半径固定为平均值29.49）为；Step3:而血管柱面上任意一点需满足Step4：从而得到血管空间曲面方程为Step5：再得到血管空间曲面方程，从进行切割，得到100个新的截痕；Step6：再对得到的100个截痕内部进行填充，得到100张新的切片图；Step7：计算新切片与原始切片的重合度。程序代码见附录6，附录7。5.2.3.检验所得结果本文随机选取第39号，第48号，第60号，第65号图，分别提取其原图轮廓和新切片图的轮廓进行对比见下图：图7 39号拟合图图6 39号原图图9 48号拟合图图8 48号原图图11 60号拟合图图10 60号原图图13 65号拟合图图12 65号原图图15 48号叠合图图14 39号叠合图图17 65号叠合图图16 60号叠合图新切片与原始切片的重合度见下表：表5 新切片与原始切片的重合度表拟合图序号10394856606586与原图重合度(%)55.0964.8774.3979.3180.2580.0563.62由上图可知,用中轴线方程重建的截面与原截面相比，重建图有一部分与原图不重合，分析这种现象产生的原因：在求最大内切圆半径时取的是平均值,显然偏离了中轴线,理想的算法应该是将几个同时为最大的点拟合,得到样条曲线,在样条曲线上取中点,才是比较合理的结果。6模型优缺点6.1.模型的优点该模型融合多个切片的信息来确定血管的中轴线，可以较精确地从血管的切片数字图像中得到血管中轴线和管半径。本模型所采用的算法简明且易于实现，效率较高，可以在题目所允许的范围内保证较高的精度，并且采取了对比的方法对模型进行了检验，取到了较好的效果。模型用到的数学知识是简单易懂的。多项式曲线拟合内容现在已广泛为人所知且在实际中用得很广。模型结果直观形象，具有非常高的实用和推广价值6.2模型的缺点该模型在求时花费时间较大。综上所述，本模型基本上较好地解决本题假设下基于切片数字图像的血管的三维重建问题，结果稳定性和精度令人满意。参考文献1陈凌钧,骆岩林.等径管道的三维重建J. 高校应用数学学报.Vol.13 Ser A Suppl.1998,87-90附录附录1：图像二值矩阵的0-1互换的matlab程序代码（zhuanhua.m）function b0=zhuanhua(b0) %图像二值矩阵的0-1互换for i=1:512 for j=1:512 if b0(i,j)=1 b0(i,j)=0; else b0(i,j)=1; end endend附录2：求各切片的最大内切圆的半径及圆心坐标matlab程序代码（ff.m）function r, zhongxindian=ff %输出各切片最大内切圆半径及圆心坐标a=zeros(512,512);b=zeros(512,512);for i=1:512 for j=1:512 a(i,j)=i-257; %横坐标的对应 b(i,j)=j-257; %纵坐标的对应 endend %图像在xyz面上的x轴、y轴坐标zhongxindian=zeros(100,2);r=zeros(100,1);for k=0:99 t=strcat(f:/,int2str(i),.bmp);b=imread(t);b=zhuanhua(b);%将01互换blunkuo=edge(b,sobel);%提取轮廓bgujia=bwmorph(b,skel,inf);%提取骨架%寻找内切圆x0,y0,v0=find(b0lunkuo);a0,b0,c0=find(b0gujia);m=length(a0);n=length(x0);juli=zeros(m,n);cunfang=zeros(m,2);for i=1:m for j=1:n p1=a0(i);q1=b0(i); p2=x0(j);q2=y0(j);juli(i,j)=sqrt(a(p1,q1)-a(p2,q2)2+(b(p1,q1)-b(p2,q2)2);%骨架上的各个点到轮廓的距离 end zx,zxxh=min(juli(i,:);%骨架上一点到轮廓的最短距离即以骨架上各个点为圆心的内切园的半径 cunfang(i,1)=zx; cunfang(i,2)=zxxh;endzd,zdxh=max(cunfang(:,1);%寻找半径中最大的半径和其对应的圆心坐标g=a0(zdxh);h=b0(zdxh);zhongxindian(k+1,1)=a(g,h);zhongxindian(k+1,2)=b(g,h);r(k+1)=zd;end附录3：通过计算不同次数多项式拟合的偏差平方和确定拟和次数的matlab程序代码（pczx.m）function j=pczx(z,t) %根据不同次数的多项式拟合与原图数据偏差平方和的大小来确定多项式拟和的次数delta=zeros(10,1);for k=1:10p,s=polyfit(z,t,k);delta(k)=s.normrendi,j=min(delta);附录4：根据轮廓画出血管的三维图像的matlab程序代码for b=0:99 %提取原图的轮廓，根据轮廓画出血管的三维图像 m1=imread(int2str(b),.bmp); m(:,:,b+1)=edge(m1,sobel);endfor k=0:99 for i=1:512 for j=1:512 if (m(i,j,k+1)=1) plot3(i,j,k+1,r-.); hold on end end endendgrid ontitle(血管三维图)rotate3dhold off附录5：绘制中轴线及在各平面的投影图matlab程序代码format longpx=polyfit(z,x,7);%x,z的7次多项式拟合x1=polyval(px,z);py=polyfit(z,y,5);%y,z的5次多项式拟合y1=polyval(py,z);figure(1); %画中心轴线图plot3(x1,y1,z)grid onxlabel(X轴);ylabel(Y轴);zlabel(Z轴);title(血管中轴线图);figure(2); %画中心轴线在xoz平面上的投影plot(z,x1,-r)ylabel(Z轴);xlabel(X轴)title(血管中轴线XOZ平面投影图);grid onfigure(3);%画中心轴线在yoz平面上的投影 plot(z,y1,-b)xlabel(Z轴);ylabel(Y轴);title(血管中轴线YOZ平面投影图);grid onfigure(4);%画中心轴线在xoy平面上的投影plot(x1,y1,-g)xlabel(X轴);ylabel(Y轴);title(血管中轴线XOY平面投影图);grid on附录6：求第pn张拟合图的轮廓的二值矩阵的matlab程序代码（dian.m）function pnjj=dian(px,py,pn) %输出pnjj，为第pn张拟合图片的轮廓二值矩阵a=zeros(991);b=zeros(991);q=zeros(991);w=zeros(991);r=zeros(1,2);s=zeros(1,2);k=1;for c=0:0.1:99a(k)=7*px(1)*c.6+6*px(2)*c.5+5*px(3)*c.4+4*px(4)*c.3+3*px(5)*c.2+2*px(6)*c+px(7);b(k)=4*py(1)*c.3+3*py(2)*c.2+2*py(3)*c+py(4); %中心轴线方程关于z的导数即a(k),b(k),1为z在k处的切线的方向向量q(k)=px(1)*c.7+px(2)*c.6+px(3)*c.5+px(4)*c.4+px(5)*c.3+px(6)*c.2+px(7)*c+px(8);w(k)=py(1)*c.5+py(2)*c.4+py(3)*c.3+py(4)*c.2+py(5)*c; %中心轴线方程在z=k处的x,y值k=k+1;end%提取新的截痕u=;v=;syms x yk=1;for i=0:0.1:99 m=a(k)*(x-q(k)+b(k)*(y-w(k)+(pn-i); n=(x-q(k)2+(y-w(k)2+(pn-i)2-29.492; g,h=solve(m,n); r=double(g); s=double(h);if (abs(imag(r(1)0.01) %去除复数根 u=u;real(r(1)+256 real(r(2)+256; v=v;real(