血管的三维重建.docx

1、血管的三维重建 血管的三维重建问题 摘要越来越多的在医学领域及生物领域需要大量的切片工作，来研究切片表面的组织形态，因而如何能将这些切片在利用计算机进行复原，显得很重要。在对管道的半径的求解中，我们采用的是枚举法，即逐个将100切片的最大内切圆的半径都求出来，最后取其平均值，即为我们所求的管道的半径。在半径求出的同时把最大内切圆的圆心也求出来。要求得切片的最大内切圆，首先对切片图像0-1矩阵转化，然后利用edge函数找到图形的轮廓，利用bwmorph函数求得图形的轮廓，然后利用find函数求出轮廓与骨架上的坐标。利用两点间的距离公式，我们就能求出骨架上每个点到轮廓的距离，在所有点到轮廓的最短距

2、离中的那个最大值，就是我们要求的最大内切圆半径。此圆心即为中轴线上的点。在中轴线与中轴线在XY，YZ，ZX平面的投影的计算中，我们主要采用多项式拟合polyfit的方法，从得出拟合曲线可以看出与实际图形很接近。最后我们拟合出了血管的三维形态。 关键词：最大内切圆 半径 中轴线 多项式拟合 一问题重述 如今在医学领域，生物学领域等都需要了解生物组织、器官等的断面形态，以便于能够准确的研究该样品，得出具有说服力的研究成果。例如，将样本染色后切成厚约1 m的切片，在显微镜下观察该横断面的组织形态结构。如果用切片机连续不断地将样本切成数十、成百的平行切片， 可依次逐片观察。这时候我们如果我们需要知道血

3、管的三维形态，根据拍照并采样得到的平行切片数字图象，运用计算机就可重建组织、器官等的三维形态。假设某些血管可视为一类特殊的管道，该管道的表面是由球心沿着某一曲线（称为中轴线）的球滚动包络而成。 本文给出了某管道的相继100张平行切片图象，记录了管道与切片的交。图象文件名依次为0.bmp、1.bmp、 99.bmp，格式均为BMP，宽、高均为512个象素（pixel）。我们要做的就是根据这些切片计算管道的中轴线与半径，给出具体的算法，并绘制中轴线在XY、YZ、ZX平面的投影图。 二、问题分析 本题是一个求最优解的问题，通过具体的算法得出最接近实际血管的中轴线与半径，并能较准确的绘制出中轴线在XY

4、、YZ、ZX平面的投影图，得出直观的图形。需要处理大量的图片，得出大量的数据信息，所以为处理该问题的程序中需要大量的循环，所以程序的完成与运行皆很复杂。本题管道为由球心沿着某一曲线的球滚动包络而成，且管道中轴线与每张切片被假设有且只有一个交点，则每张切片一定包含球的最大截面，即过球心的圆面，则每张切片的最大内切圆就是过球心的圆面 。由于本题中取坐标系的Z轴垂直于切片，第1张切片为平面Z=0，第100张切片为平面Z=99。Z=z切片图象中象素的坐标依它们在文件中出现的前后次序为（-256，-256，z），（-256，-255，z），（-256，255，z），（-255，-256，z），（-255

5、，-255，z），（-255，255，z），（ 255，-256，z），（ 255，-255，z），（255，255，z）。则每个切片的最大内切圆的圆心的z坐标已给，再加上x,y,则中心轴就是每个切片的最大内切圆的圆心(x,y,z)连成的光滑曲线，而管道的半径就是这个球的半径，即最大内切圆的半径。进而再求出中轴线在XY、YZ、ZX平面的投影图。那么这道题的问题最后就转化为求每张切片的最大内切圆的圆心与半径问题。 三、模型假设 1、管道的表面是由球心沿着某一曲线（称为中轴线）的球滚动包络而成。2、管道中轴线与每张切片有且只有一个交点。3、球半径固定。4、切片间距以及图象象素的尺寸均为1。 四、符

6、号说明一张切片内骨架点i到边缘点的距离第i张图的最大内切圆半径 半径的平均值 五、模型的建立与求解问题一：求管道的半径与每张切片的最大内切圆的圆心.要求得管道的半径，首先需要求出每个切片的最大内切圆，为了得出比较精确的结果，我们采用枚举法。对与100张切片，考虑到求得的最大内切圆半径由于误差的存在肯定会有所不同，我们将所有切片的最大内切圆半径求出后，然后取其平均值。.为得到每张图片的最大内切圆半径，首先先将二值图片转换成0-1矩阵，矩阵横纵坐标对应原图像的直角坐标系位置（1代表黑像素，0代表白像素），为了后面find函数寻找矩阵中为1的坐标并记录。 通过相关函数找出轮廓(edge)和骨架（bw

7、morph）记录轮廓q(i,j)和骨架p(i,j)的坐标，利用两点间的公式去求骨架上的点到边界上所有点的最小距离即骨架上此点的内切圆的半径。 （1） 骨架上所有点的最大内切圆即为该切片的最大内切圆， (2) (这里i代表第i张切片)程序见附录1其半径即为最大内切圆半径。通过循环找出最大半径并记录该半径所对的坐标点。基于此我们能求出所有切片的最大内切圆的半径与圆心，数据见附录5通过求得的100张切片的最大内切圆的半径我们取其平均值得出更加精确的管道的半径。问题二：求中心轴与中心轴在XY，YZ，ZX平面的投影 （x,y,z）为切片的最大内切圆的圆心坐标1.利用问题一中求得数据（x,y,z）,然后根

8、据多项式拟合polyfit得出x,y,z的拟合方程为：中轴线在Z-Y平面的投影拟合多项式：中轴线在X-Y平面的投影拟合多项式：中轴线在Z-X平面的投影拟合多项式：程序见附录2得到的拟合曲线如下：图1 为中轴线在XY平面的投影，曲线代表多项式拟合的结果，折线代表真实数据的图形图2 为中轴线在ZY平面的投影，曲线代表多项式拟合的结果，折线代表真实数据的图形图为中轴线在ZX平面的投影图，曲线代表多项式拟合的结果，折线代表真实数据的图形 从以上拟合图形中我们可以看出，通过多项式拟合的投影图与通过坐标点得出的投影图基本相同，所以可以近似认为拟合的投影图即为我们所要求得的中轴线在XY，YZ，ZX平面的投影

9、。2.同样利用多项式拟合polyfit得出管道的中轴线在三维坐标中的拟合曲线如下：（程序见附录3）图为中轴线在三维坐标的立体图3.利用球坐标以及以上求出的最大半径、以及拟合时在z时的坐标，x1,y1利用plot3，程序如附录四：图血管的三维重建 六、模型评价及总结 解决本题的首要前提是要理解像素的概念，在解此题的最初由于对图片的像素没有明确的理解而浪费了很长时间。在我们的模型中还有尚有不完美的地方，在求切片的最大内切圆半径与圆心的时候，没有编写出能够一下就将所有切片的数据都求出的程序，而是仅仅编写出一个图片的程序，然后逐一求出，在这个问题上花费的时间长，且过程繁琐。但是在进行曲线拟合时，由于运

10、用知道了相关函数的用法，运用polyfit进行拟合，很大程度上节省了时间，拟合出的相对来说比较满意。立体切片的三维重建在医学与生物学等领域的作用很大，因为这些领域往往会研究一些物体切片的组织结构，所以此问题的实现能够解决很多现实问题。为一些研究领域提供很大的帮助。 七、参考文献【1】董辰辉，MATLAB2008全程指南，电子工业出版社。【2】张红云等，图像骨架的提取的应用，2010年第四期。【3】附录1.jieguo=zeros(100,4); J0=imread(E:99.BMP); for i=1:512; for j=1:512; j0(i,j)=1-J0(i,j); 转化为黑色为1，白

11、色为0，为了后面find函数寻找矩阵中 为1的坐标并记录 end endlk=edge(j0,sobel);gj=bwmorph(j0,skel,inf);x0,y0,v0=find(lk);找到轮廓灰色区域a0,b0,c0=find(gj);找到骨架灰色区域m=length(a0);骨架灰色区域个数n=length(x0);轮廓灰色区域个数jl=zeros(m,n);建立0矩阵为求内切圆半径cf=zeros(m,2);建立两列0矩阵为存放中心点坐标for i=1:mfor j=1:np1=a0(i);q1=b0(i);p2=x0(j);q2=y0(j);骨架、轮廓坐标赋值jl(i,j)=sq

12、rt(p1-p2)2+(q1-q2)2);endzx,zxxh=min(jl(i,:); 骨架上一点到轮廓的最短距离，即骨架上各个点的内切圆的半径cf(i,1)=zx;cf(i,2)=zxxh;endzd,zdxh=max(cf(:,1); 找到其中最大的半径，并把半径和圆心坐标存储x=a0(zdxh)-256; 与题目所给坐标轴对应y=b0(zdxh)-256; 与题目所给坐标轴对应jieguo(k+1,1)=x; x轴坐标jieguo(k+1,2)=y;jieguo(k+1,3)=k+1;jieguo(k+1,4)=zd;2.polyfit(Z,X,7)plot(Z,X,0:10:99,p

13、olyval(ans,0:10:99)p,s=polyfit(Z,X,7)polyfit(Z,Y,7)plot(Z,Y,0:10:99,polyval(ans,0:10:99)p,s=polyfit(Z,Y,7)plot(X,Y,-160:10:188,polyval(ans,-160:10:188)3.format longpx=polyfit(z,x,7);x1=polyval(px,z); 取拟合后再z点时的点py=polyfit(z,y,7);y1=polyval(py,z);figure(1);plot3(x1,y1,z,r).t=linspace(0,pi,25);存储平均分pi为

14、25份存储p=linspace(0,2*pi,25);存储平均分2pi为25份存储theta,phi=meshgrid(t,p);存储角度for i=1:100 ceen(:,1)=x1;ceen(:,2)=y1;ceen(:,3)=z; x=29.49288*sin(theta).*sin(phi)+ceen(i,1); y=29.49288*sin(theta).*cos(phi)+ceen(i,2); z1=29.49288*cos(theta)+ceen(i,3); hold on 使图像不被覆盖 plot3(x,y,z1); axis equal;hold off end5、 切片号

15、半径x坐标y坐标z坐标切片号半径x坐标y坐标z坐标129.069-160005130.414-11411750228.284-160215230.414-11411751329-160125330.414-11311852429.069-160135430.414-11211953529.069-160245530-11112054629.069-160255629.732-11112055729-160165729.698-11112056829.017-160475829.698-11112057929-160185929.53-81142581028.862-160196029.547-

16、51156591128.862-1607106129.547-51156601228.862-1608116229.614-31162611328.862-1609126329.614-31162621429.017-16010136429.614-31162631529.017-16012146529.614-35161641629.017-16013156629.614-35161651729.017-16014166729.428-26163661829.017-16016176829.411-35161671929.017-16017186929.275-26163682029.017

17、-16018197029.42846163692129.017-16019207129.61446163702229.017-16020217229.61446163712329.017-16021227329.61446163722429.017-16022237429.61465158732529.069-16021247529.73268157742629.069-16021257629.73265158752729.069-16021267729.54781152762829.155-15930277829.5381152772929.275-15930287929.538115278

18、3029.275-15929298029.698135117793129.428-15835308129.698136116803229.614-15740318229.698136116813329.614-15740328329.698137115823429.614-15740338429.698138114833529.614-15644348529.698138114843629.732-15355358629.698139113853729.732-15355368729.698139113863829.732-15355378829.698139113873929.732-152

19、58388929.698140112884029.614-15258399029.698140112894129.547-15063409129.5317267904229.682-13592419229.5317267914329.53-14868429329.5317267924429.682-119112439429.5317267934529.698-118113449529.73218243944629.698-117114459629.61418724954729.698-117114469729.61418724964830.414-116115479829.61418724974930.414-115116489929.61418724985030.414-1151164910029.4281881899