血管的三维重建
血管的三维重建问题
摘要
越来越多的在医学领域及生物领域需要大量的切片工作,来研究
切片表面的组织形态,因而如何能将这些切片在利用计算机进行复
原,显得很重要。
在对管道的半径的求解中,我们采用的是枚举法,即逐个将100
切片的最大内切圆的半径都求出来,最后取其平均值,即为我们所求
的管道的半径。在半径求出的同时把最大内切圆的圆心也求出来。要
求得切片的最大内切圆,首先对切片图像0-1矩阵转化,然后利用edge
函数找到图形的轮廓,利用bwmorph 函数求得图形的轮廓,然后利
用find 函数求出轮廓与骨架上的坐标。利用两点间的距离公式,我们
就能求出骨架上每个点到轮廓的距离,在所有点到轮廓的最短距离中
的那个最大值,就是我们要求的最大内切圆半径。此圆心即为中轴线
上的点。
在中轴线与中轴线在XY ,YZ ,ZX 平面的投影的计算中,我们主要
采用多项式拟合polyfit 的方法,从得出拟合曲线可以看出与实际图
形很接近。最后我们拟合出了血管的三维形态。
关键词:最大内切圆 半径 中轴线 多项式拟合
一.问题重述
如今在医学领域,生物学领域等都需要了解生物组织、器官等
的断面形态,以便于能够准确的研究该样品,得出具有说服力的研究
成果。例如,将样本染色后切成厚约1 m的切片,在显微镜下观察
该横断面的组织形态结构。如果用切片机连续不断地将样本切成数
十、成百的平行切片, 可依次逐片观察。这时候我们如果我们需要
知道血管的三维形态,根据拍照并采样得到的平行切片数字图象,运
用计算机就可重建组织、器官等的三维形态。假设某些血管可视为一
类特殊的管道,该管道的表面是由球心沿着某一曲线(称为中轴线)
的球滚动包络而成。
本文给出了某管道的相继100张平行切片图象,记录了管道与切
片的交。图象文件名依次为0.bmp 、1.bmp 、…、 99.bmp,格式均为
BMP ,宽、高均为512个象素(pixel )。我们要做的就是根据这些切
片计算管道的中轴线与半径,给出具体的算法,并绘制中轴线在XY 、
YZ 、ZX 平面的投影图。
二、问题分析
本题是一个求最优解的问题,通过具体的算法得出最接近实际血
管的中轴线与半径,并能较准确的绘制出中轴线在XY 、YZ 、ZX 平面
的投影图,得出直观的图形。需要处理大量的图片,得出大量的数据
信息,所以为处理该问题的程序中需要大量的循环,所以程序的完成
与运行皆很复杂。
本题管道为由球心沿着某一曲线的球滚动包络而成,且管道中轴
线与每张切片被假设有且只有一个交点,则每张切片一定包含球的最
大截面,即过球心的圆面,则每张切片的最大内切圆就是过球心的圆
面 。
由于本题中取坐标系的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 )。
则每个切片的最大内切圆的圆心的z 坐标已给,再加上x,y, 则中心
轴就是每个切片的最大内切圆的圆心(x,y,z)连成的光滑曲线,而管
道的半径就是这个球的半径,即最大内切圆的半径。进而再求出中轴
线在XY 、YZ 、ZX 平面的投影图。那么这道题的问题最后就转化为求
每张切片的最大内切圆的圆心与半径问题。
三、模型假设
1、管道的表面是由球心沿着某一曲线(称为中轴线)的球滚动包络
而成。
2、管道中轴线与每张切片有且只有一个交点。
3、球半径固定。
4、切片间距以及图象象素的尺寸均为1。
四、符号说明
五、模型的建立与求解
问题一:求管道的半径与每张切片的最大内切圆的圆心
1. 要求得管道的半径,首先需要求出每个切片的最大内切圆,
为了得出比较精确的结果,我们采用枚举法。对与100张切片,考虑
到求得的最大内切圆半径由于误差的存在肯定会有所不同,我们将所
有切片的最大内切圆半径求出后,然后取其平均值。
2. 为得到每张图片的最大内切圆半径,首先先将二值图片转换
成0-1矩阵,矩阵横纵坐标对应原图像的直角坐标系位置(1代表黑
像素,0代表白像素),为了后面find 函数寻找矩阵中为1的坐标并
记录。 通过相关函数找出轮廓(edge)和骨架(bwmorph )记录轮廓
q(i,j)和骨架p(i,j)的坐标,
利用两点间的公式去求骨架上的点到边界上所有点的最小距离即骨
架上此点的内切圆的半径。
(1)
骨架上所有点的最大内切圆即为该切片的最大内切圆,
(2)
(这里i 代表第i 张切片)
程序见附录1
其半径即为最大内切圆半径。通过循环找出最大半径并记录该半
径所对的坐标点。基于此我们能求出所有切片的最大内切圆的半径与
圆心,数据见附录5
通过求得的100张切片的最大内切圆的半径我们取其平均值得出更加精确的管道的半径。
问题二:求中心轴与中心轴在XY ,YZ ,ZX 平面的投影
(x,y,z )为切片的最大内切圆的圆心坐标
1. 利用问题一中求得数据(x,y,z ), 然后根据多项式拟合polyfit
得出x,y,z 的拟合方程为:
中轴线在Z-Y 平面的投影拟合多项式:
中轴线在X-Y 平面的投影拟合多项式:
中轴线在Z-X 平面的投影拟合多项式:
程序见附录2
得到的拟合曲线如下:
图1 为中轴线在XY 平面的投影,曲线代表多项式拟合的结果,折线代表真实数据的图形
图2 为中轴线在ZY 平面的投影,曲线代表多项式拟合的结果,折线代表真实数据的图形
图3为中轴线在ZX 平面的投影图,曲线代表多项式拟合的结果,折线代表真实数据的图形
从以上拟合图形中我们可以看出,通过多项式拟合的投影图与通
过坐标点得出的投影图基本相同,所以可以近似认为拟合的投影图即
为我们所要求得的中轴线在XY ,YZ ,ZX 平面的投影。
2. 同样利用多项式拟合polyfit 得出管道的中轴线在三维坐标中的
拟合曲线如下:(程序见附录3)
图4为中轴线在三维坐标的立体图
3. 利用球坐标以及以上求出的最大半径、以及拟合时在z 时的坐标,
x1,y1利用plot3,程序如附录四:
图5血管的三维重建
六、模型评价及总结
解决本题的首要前提是要理解像素的概念,在解此题的最初由
于对图片的像素没有明确的理解而浪费了很长时间。在我们的模型中
还有尚有不完美的地方,在求切片的最大内切圆半径与圆心的时候,
没有编写出能够一下就将所有切片的数据都求出的程序,而是仅仅编
写出一个图片的程序,然后逐一求出,在这个问题上花费的时间长,
且过程繁琐。但是在进行曲线拟合时,由于运用help知道了相关
函数的用法,运用polyfit 进行拟合,很大程度上节省了时间,拟合
出的相对来说比较满意。立体切片的三维重建在医学与生物学等领域
的作用很大,因为这些领域往往会研究一些物体切片的组织结构,所
以此问题的实现能够解决很多现实问题。为一些研究领域提供很大的
帮助。
七、参考文献
【1】董辰辉,MATLAB2008全程指南,电子工业出版社。
【2】张红云等,图像骨架的提取的应用,2010年第四期。【3】
http://wenku.baidu.com/view/04b1d31f964bcf84b8d57b06.htm
【4】http://wenku.baidu.com/view/04b1d31f964bcf84b8d57b06.html
附录
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,白色为0,为了后面find 函数寻找矩阵中 为1的坐标并记录
end
end
lk=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: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;
end
[zd,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,polyval(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 long
px=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')
4. t=linspace(0,pi,25);存储平均分pi 为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 end
5、