文档介绍：该【血管的三维重建 】是由【才艺人生】上传分享，文档一共【21】页，该文档可以免费在线阅读，需要了解更多关于【血管的三维重建 】的内容，可以使用淘豆网的站内搜索功能，选择自己适合的文档，以下文字是截取该文章内的部分文字，如需要获得完整电子版，请下载此文档到您的设备，方便您编辑和打印。血管的三维重建
1
血管的三维重建
1摘要
序列图像的三维重建在各学科中都起到至关重要的作用,本次讨论的是血管的三维重建。首先,假设该管道是由球心沿着某一曲面的球滚动包络而成,故本次的主要目的是求出中轴线坐标及半径。
现有100张平行切片图像,本次建立的模型可分为四步;第一步,采集图形边界点数据。由于每张图片都是512*512的矩阵,故此数据很大,采用imread()函数将其读入矩阵A中。
第二步,最大内切圆寻找及半径的确定。提出两种方案,分别是切线法和最大覆盖法;从上述两种方法分析及考虑到我们所使用的工具和材料,可以得出方法二更加直观,计算机实现更容易,计算复杂度更低,所以我们采用后者。
根据以上算法,我们抽取了所有的切片图进行半径的提取,然后再求其平均值,。
第三步,轨迹的搜索。在第二步中求出了血管的半径,轨迹的搜索就可以建立在半径确定的基础上,当然我们也可以求出每一个切面图形的最大内切圆,然后得到每个圆心的坐标,即中轴线坐标,但这样做计算机的运算量会很大,同时由于最大内切圆搜索法的稳定性不高,从而会造成搜索的不精确,所以采用定半径搜索。本文提出了三种方法,分别为网格法、蒙特卡罗法和非线性规划法;本次采用非线性规划来实现。
第四步,绘制中轴线空间曲线图和在XOY、YOZ、XOZ三个平面的投影图。由定理1:切片上血管截面图的头部顶点在XOY平面上的投影点一定会落在中
轴线在XOY平面上的投影曲线上(在论文中以证明),并得出推论:切片上血管截面中中位线与中轴线在XOY面上的投影重合。
最后可由中轴线和血管半径在作图软件中达到血管的三维重建,本次的模型还存在一定的不足,其假设为管道中轴线与每个切面有且只有一个交点,事实上还存在有多个交点的情况,但为了简化模型在此做了一定的假设,故会存在一定的误差。
关键词:三维重建内切圆半径轨迹(中轴线)
注:求边界时采用了老师的思想和程序。
1
Z血管的最大内切圆的圆心的z的坐标值
R血管最大内切圆的的半径
血管最大内切圆整个圆周的离散坐标
A图像数据储存的三维矩阵
血管最大内切圆心匹配点数的函数
5分析与建立模型
1)假定管道中轴线与每张切片有且只有一个交点,球的半径固定,切片的间距以及图像象素的尺寸均为1。取坐标的z轴垂直于切片,第一张切片记为平面z=0,第二张切片记为z=1,第100张切片的平面记为z=99。要求根据以上所给条件计算管道的中轴线与半径,并绘制中轴线在XY、YZ、ZX平面的投影图。
对于某个切痕,可以作出它在空间中的结构图。如图1所示,设某一切片的切痕与中轴线的交点为O,切下这一切痕的同时也截下中心在点O的小球的一个大圆K,过点O的血管的轴截面S也是同一个小球的一个大圆。这两个大圆的交点即为小球的直径的两个端点,且这两个端点为切痕边界与大圆K的两个切点,因此可得出以下结论:
1)每一个切痕中存在一个最大内切圆且最大内切圆的圆心在管道的中轴线上。
2)每个切痕的最大内切圆的半径均相等,且等于管道的半径。
根据上面的结论,要求血管的中轴线,只需找出每一个切痕的最大内切圆的圆心。
图1某个切痕的空间结构图
通过对血管的三维重现的分析可知,我们的模型应由三部分组成:
边界点数据的采集方法有人工采集法和计算机采集法,对于数据量较小的,
1
又具有某些特殊意义的点可以由人工采集,若数量比较大且采集边界又非常复杂时,则采取机器采集。
方法1切线法
此方法可以从切面的外围轮廓线分析着手。由题意分析知道,所给的图片切面是由无数个球切面组成的。而且外围轮廓线与最大内切圆有且仅有两个交点,所以经过这两点的外围轮廓线的两条切线平行且间距最大。基于上述分析,我们可以通过找到这两条切线来找到最大内切圆的圆心及半径。在实际操作中,由于对图片的象素提取的离散性,我们在计算导数时是用差分来代替。
方法2最大覆盖法
最大覆盖法就是在切面中找到最佳的圆心位置和半径长度,从而使得由这个圆心和半径所决定的圆面,能最大面积地覆盖管道切面的图形,这样搜索到的圆一定是最大内切圆,这个圆的圆心就是我们所要找的球的球心,这个圆的半径就是我们所要找的球的半径。从上述两种方法分析及考虑到我们所使用的工具和材料,可以得出方法二更加直观,计算机实现更容易,计算复杂度更低,所以我们采用后者。具体实现中,我们先得到任意一张图片的象素矩阵,然后将用于匹配的圆根据其圆心和半径将其圆周离散(即以象素表示),并映射到的图中,其中圆周上的点为0,其余的点为1,即形成另一个象素矩阵。这两个矩阵在相同位置点上的值进行逻辑或运算,如果其值为0,则为匹配点,即此点在管道切面图形内,否则其在切面图形外。这样搜索到的匹配点最多同时半径最大的就是所要找的最大内切圆。(内部白色部位为管道切面图形,深色圈为最大匹配圆)。根据以上算法,我们抽取了所有的切片图进行半径的提取,然后再求其平均值,。
在求出半径以后,轨迹的搜索就可以建立在半径确定的基础上,当然我们也可以求出每一个切面图形的最大内切圆,然后得到每个圆心的坐标,即中轴线坐标,但这样做计算机的运算量会很大,同时由于最大内切圆搜索法的稳定
3
性不高,从而会造成搜索的不精确,所以采用定半径搜索。我们通过定圆(半径为)来找其中轴线,也就是用定圆覆盖到切面图形上去,找到匹配点数最多的一个位置,从而得到此定圆圆心的位置。具体实现时只要用定圆的圆心位置进行变化,设其为,则由A点可以得到整个圆周的离散坐标,令为匹配点数的函数,其计算方法与最大内切圆求法相同,即将根据圆心的坐标和半径(固定)离散化形成矩阵后和切面图形矩阵作逻辑与操作,从达到匹配最大,则优化目标函数就是:
坐标面上的点
,使得定圆覆盖切面图形最大,即多元函数极值最优解,这是一个以X和Y为自变量二元函数,这样可以通过以下三种方法来求得:
1)用直接搜索求最优解(网格法)
搜索过程中,对每一个圆心的坐标X和Y,在其取值的范围内均取100个步长,分为1002个网格,这样,在一定的精度范围内,可以求的一个较好的最优解。
2)蒙特卡罗法
蒙特卡罗法,也就是随机实验试点法。它的基本思想是:在函数的可行域内随机地选取实验点,由于随机取得的点在区域中分配比较均匀,所以对函数的大致形态能较好的体现。模型中,随机点用以下方法产生的。
其中,(x0,x1)为X的取值范围,(y0,y1)为Y的取值范围。
3)非线性规划
非线性规划即无约束优化,以数值迭代为基本思想,基本步骤为选取初值A(X0,Y0),进行k次迭代并求出迭代解,由迭代解得到搜索方向和步长,如果k+1次迭代符合给定的迭代终止条件,则停止迭代,得出最优解;否则继续迭代。非线性规划的关键是搜索方向、步长和初值。我们用拟牛顿法来选定搜索方向,拟牛顿法是在牛顿法基础上,克服牛顿法中黑赛阵不仅计算复杂、而且会出现变态、不正定等情况,同时保持了较快收敛的优点,从而得到最好的下降方向。搜索步长的确定使用线性搜索的方法,或更为有效的插值方法。由于此搜索程序的初值对于程序正确有效的搜索影响很大,又因为相邻的切片只有一个象素单位的距离,可以认为中轴线的变化很小,所以我们可以把前一次确定的球心作为下一次搜索的起点,从而大大提高搜索的效率和准确性。但由于转角处变化较大,在这种情况下搜索起点会郴够逼近而导致优化搜索的失效。我们采用回溯技术来避免这种现象。即当球心间距前后相差较大时(采用工程上的观点,以6倍为基线),回退一次搜索,并以当前球心作为回退搜索的起点,这样就可以把各个球心间距变得比较匀称,从而相应的消除了转角上的搜索失效。从上述对三种方法的分析可以得到网格法和蒙特卡罗法实现思路简单,程序容易实现,但网格法搜索的精度不高,误差较大,搜索时间长,同样蒙特卡罗法的实现对采点的数目要求很高,计算量大;非线性规划法的实现复杂,但搜索速度快,计算量少,而且通过MATLAB
3
的优化工具箱的函数可以很方便地实现,所以非线性规划法较优。轨迹搜索结束以后,我们可以得到中轴线的100个点,用所求半径的球定位在这些点上和轮廓线相交就可以得到图4。从上图可以看到把图形重组后出现了不平滑点,也就是说直接用这些中轴线上的点构成的管道还不平滑,所以我们采用曲线拟合的方法来精确定位中轴线。
图4
、YOZ、XOZ三个平面的投影图。
定理1:切片上血管截面图的头部顶点在XOY平面上的投影点一定会落在中
轴线在XOY平面上的投影曲线上。
证明:切片与XOY平面平行,点0,点0c分别为与切片相切的两个球的球心。M,N为两切点,0点投影在XOY平面的投影点为P,则OP垂直于平面xoy,又M为切点,所以OM垂直于切面,因为切面平行于XOY平面,所以OM延长线垂直于XOY平面于Pc,又过平面外一点作平面的垂线有且只有一条。故P与Pc重合。即M点与0点在XOY面上的投影点是重合的。推广:考虑切面与球面相交为一个圆的情况,Mc为该截面的圆心,,Mc与O点在XOY平面上的投影重合。由此得推论:切片上血管截面中中位线与中轴线在XOY面上的投影重合。
6模型的求解
,存储为三维矩阵
利用Matlab软件将100张切片图像数据储存为三维矩阵A。Matlab程序如下:
fori=0:1:99
Imname=sprintf('%',i);
4
A(:,:,i+1)=imread(imname);
End
由于BMP格式文件在计算机中是以二进制数进行储存的,每张切片保存在一个二维的有0和1组成的矩阵中,其中0和1分别对应于图像中的黑象素和白象素。根据问题中的描述可知每一张BMP格式的切片包含了个象素, 每一个象素都有自己的一个确定的坐标。在转换为矩阵存储后,A(:,:,i+1)即代表了第i张切片图像,此时象素坐标则对应地转换为矩阵的列与行。因此,以下的求解过程中,均以象素在矩阵中的位置作为坐标,其中列对应于x横坐标,行对应于y坐标。
首先,引进图像处理技术中的四邻域概念。
四邻域:某个像素的左、右、上、下四个象素称为该象素的四邻域。如图5,象素E、S、W、N称为象素O的四邻域。
N
WE
O
S
图5四邻域
5
60
70
80
90
100
110
120
130
220
230
240
250
260
270
280
290
300
图6第0张切片中血管图像的内边界
求血管图像内边界的算法:逐点找出所有边界点坐标,即对图像进行逐行搜索,当遇到灰度值为0(黑)的象素点时,再搜索其四邻域的四个点,若在其四邻域中有一个象素的灰度值为1(白),则该点就是一个内边界点。
I=A(:,:,i);
E=ones(size(I));
fori=1:512
p=find(I(i,:)==0);
if~isempty(p)
E(i,p(1))=0;
E(i,p(end))=0;
forj=2:length(p)-1
ifI(i-1,p(j))==1||I(i,p(j)-1)==1||I(i,p(j)+1)==1||...I(i+1,p(j))==1
E(i,p(j))=0;
end
end
end
end
imshow(E)
求血管图像外边界的算法:逐点找出所有边界点坐标,即对图像进行逐行搜索,当遇到灰度值为0(黑)的象素点时,再搜索其四邻域的四个点,若在其四邻域中任何一点象素的灰度值为1(白),则该四邻域中的这一点就是一个外边界点。
I=A(:,:,i);
E=ones(size(I));
fori=1:512
p=find(I(i,:)==0);
6
if~isempty(p)
forj=1:length(p)
ifI(i-1,p(j))==1
E(i-1,p(j))=0;
end
ifI(i,p(j)-1)==1
E(i,p(j)-1)=0;
end
ifI(i,p(j)+1)==1
E(i,p(j)+1)=0;
end
ifI(i+1,p(j))==1
E(i+1,p(j))=0;
end
end
end
end
最后将100张切片的血管图像边界叠加在一起,就得到了血管的三维叠加图,如图8。.在上述的算法程序中,稍作修改便可以将每张切片的边界点保存在一个二维数组中,为下一步求解半径所用。在程序中使用下列语句:
[R,C]=find(E==0);
其中表示行,表示列,即矩阵E的第行第列是边界像素。
球的半径通过对每一幅切片最大内切圆的搜索,分别得到R1、R2、图7管道中轴线拟合后的立体还原图,、R99、R100,为减少搜索过程和离散化过程所引入的误差,我们采用求均值的方法来得到球的半径,可以得到半径R为:。
管道的中轴线由轨道搜索而得,下图是由100点拟合而得的中轴线,表1中给出了直接搜索到的100点中轴线在切面上的点的坐标(依次为x、y、z)
图7
7