血管的三维重建.doc

血管的三维重建模型【摘要】本文对血管的三维重建问题进行了讨论。采用寻找最大内切圆的方法确定中轴线的离散坐标和半径。在建立三维图过程中运用切片叠加法和拟合曲线还原法，并进行了对比。首先将100张二值图象读入，通过编程求取中轴线的离散坐标和血管半径。其半径=29.416663pixel。将100张平行切片图象叠加在一起，从而得到血管的三维图。根据血管管道中轴线的离散点，采用多项式拟合数据曲线，得到中轴线的拟合曲线的参数方程用求出的半径还原出血管的三维图。将两张血管三维图进行对比和误差分析，得出最理想的三维图象。关键词：多项式拟合 血管三维图 切片叠加法 一、 问题重述 断面可用于了解生物组织、器官等的形态。例如，将样本染色后切成厚约1m的切片，在显微镜下观察该横断面的组织形态结构。如果用切片机连续不断地将样本切成数十、成百的平行切片， 可依次逐片观察。根据拍照并采样得到的平行切片数字图象，运用计算机可重建组织、器官等准确的三维形态。假设某些血管可视为一类特殊的管道，该管道的表面是由球心沿着某一曲线（称为中轴线）的球滚动包络而成。例如圆柱就是这样一种管道，其中轴线为直线，由半径固定的球滚动包络形成。现有某管道的相继100张平行切片图象，记录了管道与切片的交。图象文件名依次为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平面的投影图。二、 问题分析断面可用于了解生物组织、器官等的形状。假设某些血管可视为一类特殊的管道，该管道的表面是由球心沿着某一曲线（称为中轴线）的球滚动包络而成。现有某管道的相继100张平行切片图像，记录了管道与切片的交。首先运用matlab将每张二值图象转换成的矩阵。由于每张切片均过血管管道的中轴线，所以每张张图象中被染色部分的最大内切圆圆心必过血管管道的中轴线。求出全部最大内切圆圆心的坐标，运用拟合求出三维管道的中轴线方程。通过matlab求出每个最大内切圆的半径，求平均值得到血管管道半径。从而可以得到血管的三维图象。再将所得血管的三维图象切成厚1的100张切片。与题中所示切片对比，进行误差分析。三、 模型的假设1、切片机连续将样本切成平行切片，可依次逐片观察。2、血管视为特殊的管道，该管道的表面是由球心沿着某一曲线的球滚动包络而成。3、管道中轴线与每张切片有且只有一个交点，球半径固定，切片间距以及图象象素的尺寸为1。4、题中给出的100中平行切片图象真是有效。四、 符号说明符号 意义二值图象的0,1矩阵。图像中染色部分的最大内切圆的圆心坐标。图像中染色部分的最大内切圆的半径。血管管道的半径。血管管道的曲线方程。五、 模型的建立与求解5.1模型的建立：基本算法：Step1.读入数据,运用将100张二值图象读入，生成100个0,1矩阵 。Step2.求最大内切圆,在每个0,1矩阵中找到最大的内切圆，并确定其半径和圆心坐标分别为 ， 。Step3.血管管道半径,对所有半径求平均值，得到血管管道的半径。Step4.血管管道中轴线,对所有圆心坐标进行拟合得到三维曲线，即为血管管道中轴线的曲线方程。具体操作：1）我们用matlab7.11编程，用函数读入图象生成。2）最大内切圆的寻找及半径的确定，我们将用于匹配的圆根据其圆心和半径将其圆周离散(即以象素表示)，并映射512512的图中，其中圆周上的点为0，其余的点为1，即形成另一个象素矩阵。这两个矩阵在相同位置点上的值进行逻辑或运算，如果其值为0，则为匹配点，即此点在管道切面图形内否则其在切面图形外。这样搜索到的匹配点最多同时半径最大的就是所要找的最大内切圆。3)血管中轴线的多项式拟合：一般多项式拟合的目标是找出一组多项式系数，使得多项式能够较好地拟合原始数据。多项式拟合虽不能保证每个样本点都在拟合曲线上，但能使得整体的拟合误差较小。针对求出的圆心坐标与图象相对应的关系作具体的多项式拟合，这里我们取多项式的阶次=9。建立管道中轴线的参数式方程为：5.2模型的求解：5.2.1切片叠加法 由于样本被切成厚度为的切片，所以对100张图象进行如下图所示的叠加，可以模拟出大致的血管三维图。 矩阵m第100层m(:,:,100)矩阵m第50层m(:,:,50)矩阵m第1层m(:,:,1)第1幅图（0.bmp）第50幅图（49.bmp）第100幅图（99.bmp） 利用上述切片叠加法得到血管的粗略图1：图1 血管的三维图 5.2.2拟合曲线还原法导入二值图象后得到0,1矩阵，利用matlab编程得到每张图象最大内切圆圆心三维坐标和半径，求得平均半径=29.416663pixel具体见见附录中表1。由于管道长度很短，故得到的100个点彼此间的排列很紧密，因而我们可以用折线连接的方法得出中轴线的近似曲线，见以下图2。图2 中轴线的近似曲线图下面对所有圆心坐标进行数据曲线拟合，拟合程序见附录，得到管道中轴线的参数方程：其中、的拟合曲线图下表3。表3 参数方程的拟合曲线图血管管道中轴线的拟合曲线见以下图4。图4 中轴线的拟合曲线根据题中要求，画出中轴线在XY、YZ、XZ平面的投影图。六、模型检验和误差分析61 模型的检验：从管道的三维重建图沿z轴作切面，就可以得到切面。这样我们可以先从中取出第5张切片和第50张切片，再和第5张切面图片(即4bmp)和第50张切面图片(即49 bmp)的外轮廓线分别作在同一张图片上，验证模型的准确性。从两张可以得出结论：按照我们建立的模型求出来的解与实际的情况基本符合，切面得到的球的所有球的切面构成了整张的切面图形，外轮廓基本重合，这有力的证明了我们的模型的准确性和算法的合理性。62 误差的分析：当然整个模型存在一定的误差：(1)切面图本身存在象素上的误差(这是题目本身所带入的，非模型所致)。(2)数值离散化时带入误差，主要体现在搜索最大内切圆和中轴线轨迹时圆周的离散化运算中。(3)在最小内切圆搜索和轨迹搜索时，由于步长、初始点等因素的限制而造成误差。七、 模型的评价与改进7.1模型的优点： 1、模型的基本算法思路结构清晰，能够较好的拟合出血管管道中轴线的曲线方程。2、对血管的三维重现，我们不光运用拟合曲线还原血管的三维曲线的方法，另外还运用了切片叠加的方法从侧面佐证了拟合出的三维图的真实性。3、血管中轴线的多项式拟合，能够较好地拟合原始数据，保证了每个个样本点基本在拟合曲线上 ，能使得整体的拟合误差较小。4、本文写作方面大量使用的图形结合，能够在建模思想的叙述上更清晰。7.2模型的缺点：血管管道的半径固定，使得我们拟合出来的血管三维图与血管的实际图象有很大的偏差。当然，这也是题目中给出的假设，是无法改变的。八、参考文献姜起源，谢金星，叶俊.数学模型（第三版）.北京：高等教育出版社.2010张德丰.MALAB概率与数理统计分析.北京：机械工业出版社.2010徐士良.数值分析与算法.北京：机械工业出版社.2007韩中庚.数学建模方法及其应用.北京：高等教育出版社.2009何斌.Visual C+数字图像处理.北京：人民邮电出版社.2001张宏宇，金文仲，赵国栋，谢瑞明，刘岩，谢杰. 探讨CT血管分析技术结合仿真内镜在肾积水病因诊断中的应用价值.当代医学2011年3月第17卷第7期总第234九、附录附录1表1 最大内切圆圆心坐标XYRXYRXYRXYR19625729.068884269627729.0688845114237329.698497632141429.73213729625628.284271279627729.0688845214237329.698497733740829.54657339625829.017236289728629.1547595314337429.698497833740829.52964649625829.068884299728629.2745625414437529.698497933740829.52964659625829.068884309728529.274565514537629.681648039137429.41088269625829.068884319829129.427885614537629.206168139237329.69848579625729329929629.614195719340729.410888239237329.69848589626029.017236339929629.614195818140129.529658339337229.69848599625729349929629.614195917539829.5296468439437129.698485109625728.8617393510030029.614196020541229.5465738539437129.698485119626328.8617393610331129.732146120541229.5465738639537029.698485129626428.8617393710331129.732146222541829.6141868739537029.698485139626528.8617393810331129.732146322541829.6141868839537029.698485149626629.0172363910431429.732146422541829.6141868939636929.698485159626829.0172364010431429.614196522141729.6141869039636929.681644169626929.0172364110631929.546576622141729.6141869142832329.529646179627029.0172364210732229.546576723041929.4278789242832329.529646189627229.0172364310832429.529656822141729.4108829342832329.529646199627329.0172364410832429.529656923041929.2745629442832329.529646209627429.0172364511333429.529657030241929.4278789543829929.732137219627529.0172364611934429.410887130241929.6141869644328029.614186229627629.0172364711934429.410887230241929.6141869744328029.614186239627729.0172364814037129.698497330241929.6141869844328029.614186249627829.0172364914137229.698497432141429.6141869944328029.614186259627729.0688845014137229.698497532441329.73213710044427429.427878附录最大内切圆心和半径的程序：clcJ0=imread(C:Documents and SettingsAdministrator桌面第四题训练题3数据0.BMP);for i=1:1:512; for j=1:1:512; j0(i,j)=1-J0(i,j); endendlk=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);cf=zeros(m,2);for i=1:m for j=1:n p1=a0(i);q1=b0(i); p2=x0(j);q2=y0(j);jl(i,j)=sqrt(p1-p2)2+(q1-q2)2); end zx,zxxh=min(jl(i,:); cf(i,1)=zx; cf(i,2)=zxxh;endzd,zdxh=max(cf(:,1); g=a0(zdxh)h=b0(zdxh)r=zdL=g h r;save C:Documents and SettingsAdministrator桌面 c.txt L -ascii;曲线的拟合程序：z=29.068884 28.284271 29.017236 29.068884 29.068884 29.068884 29 29.017236 29 28.861739 28.861739 28.861739 28.861739 29.017236 29.017236 29.017236 29.017236 29.017236 29.017236 29.017236 29.017236 29.017236 29.017236 29.017236 29.068884 29.068884 29.068884 29.154759 29.274562 29.274562 29.427878 29.614186 29.614186 29.614186 29.614186 29.732137 29.732137 29.732137 29.732137 29.614186 29.546573 29.546573 29.529646 29.529646 29.529646 29.410882 29.410882 29.698485 29.698485 29.698485 29.698485 29.698485 29.698485 29.698485 29.681644 29.206164 29.410882 29.529646 29.529646 29.546573 29.546573 29.614186 29.614186 29.614186 29.614186 29.614186 29.427878 29.410882 29.274562 29.427878 29.614186 29.614186 29.614186 29.614186 29.732137 29.732137 29.546573 29.529646 29.529646 29.410882 29.698485 29.698485 29.698485 29.698485 29.698485 29.698485 29.698485 29.698485 29.698485 29.681644 29.529646 29.529646 29.529646 29.529646 29.732137 29.614186 29.614186 29.614186 29.614186 29.427878;x=96 96 96 96 96 96 96 96 96 96 96 96 96 96 96 96 96 96 96 96 96 96 96 96 96 96 96 97 97 97 98 99 99 99 100 103 103 103 104 104 106 107 108 108 113 119 119 140 141 141 142 142 143 144 145 145 193 181 175 205 205 225 225 225 221 221 230 221 230 302 302 302 302 321 324 321 337 337 337 391 392 392 393 394 394 395 395 395 396 396 428 428 428 428 438 443 443 443 44