模拟2---A题--血管的三维重建论文

A 三维血管的重建摘要本文根据三维重建理论的知识，运用计算机就血管的二维平行切片重建成血管的三维形态结构问题给予研究。问题一：计算管道的中轴线与半径，运用MATLAB将附件A中的图像全部传化成0-1的二维矩阵，再通过MATLAB编程，得到每一张图片的边缘轮廓线，找到其边缘轮廓线的最大内切圆，找出圆心（即半径固定的球的球心），将每一张图片得到的圆心拟合起来，就是管道的中轴线(见图 4 )；将每张图片的最大内切圆的半径相比取其较大者，所得即为管道的半径R，得R=29.2476问题二：绘制中轴线在XY、YZ、ZX平面的投影图，由第一问得到函数，运用MATLAB绘出图像（见正文）（）关键词：三维重建理论，二维矩阵，四邻域 一 问题的重述 断面可用于了解生物组织、器官等的形态。例如，将样本染色后切成厚约1um的切片，在显微镜下观察该横断面的组织形态结构。如果用切片机连续不断地将样本切成数十、成百的平行切片， 可依次逐片观察。根据拍照并采样得到的平行切片数字图象，运用计算机可重建组织、器官等准确的三维形态。假设某些血管可视为一类特殊的管道，该管道的表面是由球心沿着某一曲线（称为中轴线）的球滚动包络而成。例如圆柱就是这样一种管道，其中轴线为直线，由半径固定的球滚动包络形成。现有某管道的相继100张平行切片图像（模拟2附件A题bmp.zip），记录了管道与切片的交。图像文件名依次为0.bmp、1.bmp、 99.bmp，格式均为BMP，宽、高均为512个象素（pixel）。为简化起见，假设：管道中轴线与每张切片有且只有一个交点；球半径固定；切片间距以及图像象素的尺寸均为1。 取坐标系的Z轴垂直于切片，第1张切片为平面Z=0，第100张切片为平面Z=99。Z=z切片图像中象素的坐标依它们在文件中出现的前后次序为（-256，-256，z），（-256，-255，z），（-256，255，z），（-255，-256，z），（-255，-255，z），（-255，255，z），（ 255，-256，z），（ 255，-255，z），（255，255，z）。 试计算管道的中轴线与半径，给出具体的算法，并绘制中轴线在XY、YZ、ZX平面的投影图。二 模型的假设1.管道中轴线与每张切片有且只有一个交点；2.球半径固定；3.切片间距以及图像象素的尺寸均为1。4.管道的表面是由球心沿着某一曲线（称为中轴线）的球滚动包络而成；三 符号说明1.R表示管道的半径；2. a（i）表示每个图片转换成矩阵后所对应二维矩阵；3.b(i)表示边界矩阵；四 问题的分析1. 计算管道的中轴线与半径，运用MATLAB将附件A中的图像全部传化成0-1的二维矩阵，再通过MATLAB编程，得到每一张图片的边缘轮廓线，找到其边缘轮廓线的最大内切圆，找出圆心（即半径固定的球的球心），将每一张图片得到的圆心拟合起来，就是管道的中轴线(见图 )；将每张图片的最大内切圆的半径相比取其较大者，所得即为管道的半径R；2. 绘制中轴线在XY、YZ、ZX平面的投影图，由第一问得到的结果，运用MATLAB绘出图像（见正文），并对三个图像进行拟合，得到下面的图形方程五 模型的建立与求解 5.1 模型的建立首先求出滚动球的半径和中轴线的方程：1.将切片的图像数据转化为0-1的二维矩阵，并存储为三维矩阵；2.求出每张切片的边缘轮廓线：首先，引进图像处理技术中的四邻域概念：【1】四邻域:某个像素的上、下、左、右四个像素称为该像素的四邻域，如下图： 图1然后，求切片的边缘轮廓线：3.求出边缘轮廓线的最大内切圆的半径和圆心；定理：在一条粗细均匀血管的任何横截面的图像内，其包含的最大内切圆的圆心位于内切圆上，该圆的半径等于内切圆的半径。基于：（1） 球的任意截面都是圆 ；（2） 经过球心的截面是最大的；【2】4.运用MATLAB绘制中轴线在XY、YZ、ZX平面的投影图。5.2 模型的求解步骤一：将切片图象数据储存为三维矩阵 在MATLAB中运用imread函数将图片导成矩阵数组，并储存起来。求解程序为：for i=1:100 ai=imread(strcat(num2str(i-1),.bmp);end因为bmp格式的图片在存储时是以二进制格式储存的，每个图片保留在一个二进制的矩阵中，0代表黑点1代表白点，由题可得每个bmp格式的图片都有512 512个像素点每一个点上都对应了一个0或1。在转换成矩阵后每个图片都有一个对应矩阵a（i），每个像素点就对应一个矩阵坐标 。所以在下面的求解过程中我们以矩阵的列表示坐标的x轴矩阵的行表示y轴。步骤二：图形边界曲线的生成 任取一张切片图的矩阵a（i），在对应坐标系中任找一个0点Pij如果将其四周的四个点相加不为0则此点保留，反之不保留，找完所有的0点即可得到此图片的边界图矩阵b(i)。将100张图片都做同样的处理就可以得到全部的边界图形。下图为其中的第1、第98张切片的轮廓线图2 图3具体编程如下：【3】for i=1:100 ai=imread(strcat(num2str(i-1),.bmp);endb=a;d=;for i=1:100 ddd=; for j=2:(length(ai)-1) for k=2:(length(ai)-1) if ai(j-1,k)+ai(j+1,k)+ai(j,k-1)+ai(j,k+1)=0 bi(j,k)=1; ddd=ddd;j,k; end end end di=ddd;endc=;for i=1:100 ddd=; for j=1:length(ai) for k=1:length(ai) if bi(j,k)=0; ddd=ddd;j,k; end end end ci=ddd; endddd=;save(bianjieshuju.mat,a,b,c,d););步骤三：中轴线的生成取出第一张图片并在区域d1中任意取一点dij1找出它与区域边界的最小距离h1，然后再在dij附近找点（如：d(i+1)j1,d(i-1)j1.d(i+1)(j+1)1 ）确定其与边界的最小距离hn直到hn比其附近所有点对应的最小距离都大时，此时的点dij1即为要找的焦点；然后取出第二张图片找出dij2到边界的最短距离，然后再在dij2附近找点直到h”n比其附近所有点对应的最小距离都大时，并取出焦点d”ij2；再取出第三张图片在点2(d”ij2)-dij1附近找点直到h(3)n比其附近所有点对应的最小距离都大时，并取出焦点d”ij3；以此类推即可得到所有的焦点。图片中点到边界最短距离的最大值即为图片内接圆的半径，全部图片的半径的均值可看做是球的半径。将每张图片的焦点用matlb画出散点图，然后拟合为函数，此函数即为管道的中轴线函数。血管的中轴线：0100200300400500250300350400450020406080100xyz图4具体编程如下（见附录1）绘制中轴线在XY、YZ、ZX平面的投影图在上面我们已经找出了每个焦点的坐标,当我们将每个焦点z取0时焦点将回落在xy平面上，在平面直角坐标系中将每一个改进后的点在标注出来就可以拟合出中轴线在xy平面上的投影，同理当我们依次将x,y取0后，就可以拟合出中轴线在yz平面xz平面上的投影图。其拟合方程如下：所得图像： 六 模型的改进与评价优点：该模型结合多张切片的信息来确定血管中轴线，可以较精确的从血管切片数字图像中得到血管中轴线和管半径，有较高稳定性和精度。缺点：本模型的求解过程需要借助大量编程，且未知量多，数字处理量大，在编程过程中程序的调试需要花费大量精力，且易出错。参考文献【1】 /view/4b0ae621bcd126fff7050b8d.html【2】 /p-4949231.html【3】 Visual C+数字图象处理第12页2.3.1节。何斌等编著，人民邮电出版社，2001年4月；附录1：function zuobiao,zuida=qiuxin(chushizuobiao,lunkuo)%UNTITLED1 Summary of this function goes here% Detailed explanation goes heredianshu=length(lunkuo); zuobiao=chushizuobiao; kong=1; while kong zuida=0; for j=(zuobiao(1)-2):(zuobiao(1)+2) for k=j=(zuobiao(2)-2):(zuobiao(2)+2) juli=100; for i=1:dianshu dx=j-lunkuo(i,1); dy=k-lunkuo(i,2); dz=zuobiao(3)-lunkuo(i,3); juli=dx*dx+dy*dy+dz*dz; if juli0zuida zuida=juli; x0=j; y0=k; end end end if x0=zuobiao(1)&y0=zuobiao(2) kong=0; else zuobiao(1)=x0;zuobiao(2)=y0; end end cengshu=length(c);e=;m=10;for i=1:cengshu; % if i=1 zuobiao=di(1,:),i; lunkuo=; for ii=max(i-30,1):min(i+30,cengshu); lunkuo=lunkuo;cii,ii*ones(length(cii),1); end % dianshu=length(lunkuo); kong=1; while kong zuida=0; for j=(zuobiao(1)-m):(zuobiao(1)+m) for k=(zuobiao(2)-m):(zuobiao(2)+m) juli=10000; for i0=1:dianshu dx=j-lunkuo(i0,1); dy=k-lunkuo(i0,2); dz=zuobiao(3)-lunkuo(i0,3); juli0=dx*dx+dy*dy+dz*dz; if juli0zuida zuida=juli; x0=j; y0=k; end end end if x0=zuobiao(1)&y0=zuobiao(2) kong=0; else zuobiao(1)=x0;zuobiao(2)=y0; end end % e=e;zuobiao,zuida; % elseif i=2 zuobiao=e(1,1:2),i; lunkuo=; for ii=max(i-30,1):min(i+30,cengshu); lunkuo=lunkuo;cii,ii*ones(length(cii),1); end % dianshu=length(lunkuo); kong=1; while kong zuida=0; for j=(zuobiao(1)-m):(zuobiao(1)+m) for k=(zuobiao(2)-m):(zuobiao(2)+m) juli=10000; for i0=1:dianshu dx=j-lunkuo(i0,1); dy=k-lunkuo(i0,2); dz=zuobiao(3)-lunkuo(i0,3); juli0=dx*dx+dy*dy+dz*dz; if juli0zuida zuida=juli; x0=j; y0=k; end end end if x0=zuobiao(1)&y0=zuobiao(2) kong=0; else zuobiao(1)=x0;zuobiao(2)=y0; end end % e=e;zuobiao,zuida; % else zuobiao=e(i-1,1:3)*2-e(i-2,1:3); lunkuo=; for ii=max(i-30,1):min(i+30,cengshu); lunkuo=lunkuo;cii,ii*ones(length(cii),1); end % dianshu=length(lunkuo); kong=1; while kong zuida=0; for j=(zuobiao(1)-m):(zuobiao(1)+m) for k=(zuobiao(2)-m):(zuobiao(2)+m) juli=10000; for i0=1:dian