利用MATLAB实现CT断层图像的三维重建

1、文章编号:1004-4140(200402-0024-06利用MATLAB实现CT断层图像的三维重建曾筝1,董芳华1,陈晓2,周宏2,周建中1(1 华中科技大学水电与数字化工程学院,武汉430074; 2 总后军需装备研究所,北京100088摘要:医学图像三维重建为人体结构提供了真实、直观的反映,便于医学人员对病灶的观察及手术的进行。但图像三维重建编程实现困难,不易被非计算机专业人士所掌握。本文运用MA TLAB6.5软件中的图像处理工具箱实现了CT断层图像的三维表面重建及体重建,原理简单,编程实现方便。在对头部CT图片进行的三维表面重建及体重建实验中,重建速度快,显示效果良好,便于各类非计算

2、机专业人士推广应用。关键词:CT图像;表面重建;体重建;MA TLAB中图分类号:TP391 文献标识码:AThree Dimensions Reconstruction of CT Image by MATLAB ZENG Zheng1,DONG Fang-hua1,CHEN Xiao2,ZHOU Hong2,ZHOU Jian-zhong1(1. College of Hydroelectricity & Digital Engineering, HUST, Wuhan 430074, China;2. The Quartermaster Research Institute of

3、 the General Logistics Department of CPLA, Beijing 100088,ChinaAbstract: Three dimensions reconstruction of CT image gives a intuitionistic image of body. It is convenient to observation and operation. But it is difficult to those people not major in computer. In order to solve this problem, how to

4、reconstruct of CT image using image processing toolbox of MATLAB6.5 is discussed. MABLAB gives us a convenient method to reconstruct and acquires a good result. It can be mastered by anybody.Key words: CT image; surface reconstruction; volume reconstruction; MABLAB1引言随着计算机科学技术的发展及其与其他各学科间的相互交叉渗透,特别是

5、计算机断层扫描(CT、核磁共振扫描(MIT等技术逐渐应用于临床医学,使得获取各种二维医学图像的能力得到了空前的提高1。如何将这些二维数据转化为更容易被接受和处理的信息,正是医学数据三维可视化所要解决的问题。三维重建技术的实现方法包括两种2:一种是通过几何单元拼接拟合物体表面来描述物体的三维结构,称为表面重建;另一种是直接将体像素以一定的颜色和透明度投影到显示平面的方法,称为体重建。比较而言,表面重建运算量小,表面显示清晰,但对边缘检测的要求比较高;而体重建直接基于体数据进行显示,避免了重建过程中所造成的伪像痕迹,但运算量较大3。无论采用哪种重建方法,对于非计算机专业的人士,编程实现都是很困难的

6、,2期曾筝等:利用MA TLAB实现CT断层图像的三维重建 25 因此,本文介绍了基于MATLAB进行三维表面重建及体重建的方法,实现起来灵活、方便,便于掌握。三维重建步骤如图1所示。 图1 三维重建流程图为了有利于从CT图像中准确地提取出有用的信息,需要对原始图像进行预处理,以突出有效的图像信息,消除或减少噪声的干扰。(1 CT图像格式的转换与读写。原始的CT图像是采用DICOM3.0标准进行存储的,不能被MATLAB所识别,因此必须进行图像格式的转换。在这里,使用Visual C+6.0自行开发转换程序,在正确读取DICOM图像后,通过选择合适的窗宽、窗位,将窗宽范围内的值通过线性或非线性

7、变换转换为小于256的值,将CT图像转换为256色BMP图像4。使用MATLAB中的图像读入函数imread( ,可以读取BMP图像,同时可以使用图像写出函数imwrite( 及图像显示函数image( 、imshow( 对图像进行写出和显示5。(2 图像增强。图像增强就是根据某种应用的需要,人为地突出输入图像中的某些信息,从而抑制或消除另一些信息的处理过程6。使输入图像具有更好的图像质量,有利于分析及识别。图像增强技术主要包括直方图修改、图像平滑、图像边缘锐化及伪彩色增强等。分别采用以下函数实现:灰度直方图均衡化。均匀量化的自然图像的灰度直方图通常在低灰度区间上频率较大,使得图像中较暗区域中

8、的细节看不清楚,采用直方图修整可使原图像灰度集中的区域拉开或使灰度分布均匀,从而增大反差,使图像的细节清晰。直方图均衡化在MATLAB 中使用histeq( 函数实现。灰度变换法。照片或电子方法得到的图像,常表现出低对比度即整个图像偏亮或偏暗,为此需要对图像中的每一个像素的灰度级进行标度变换,扩大图像灰度范围,以达到改善图像质量的目的。这一灰度调整过程可以用imadjust( 函数实现。平滑与锐化滤波。平滑技术用于平滑图像中的噪声,基本采用在空间域上的求平均值或中值,或在频域上采取低通滤波。在MATLAB中,各种滤波方法都是在空间域中通过不同的卷积模板即滤波算子实现,可用fspecial( 函

9、数创建预定义的滤波算子,然后用filter2( 或26 CT 理论与应用研究 13卷conv2( 函数在实现卷积运算的基础上进行滤波。中值滤波是一种基于排序统计理论的抑制噪声的非线性信号处理技术,其在除去图像中的孤立点、线的噪声的同时,很好地保护了图像的边缘信息,适用于一些线性滤波器无法胜任地场合。在MATLAB 中使用medfilt( 函数实现中值滤波7。适当运用上述方法对原始图像进行处理,将使原始图像变得较清晰,能够较真实地反映图像的结构特征,便于三维重建的处理及显示。2CT 图像三维表面重建计算机三维表面重建是指首先运用图像技术从二维图像中分割出兴趣区的轮廓曲线,然后经图形处理,得到其三

10、维结构,从而再现原物体的空间结构。因此,对于三维表面重建而言,边界轮廓的提取尤为重要8。为了便于面部边界的提取,先对各CT 图片进行颜色处理,去掉非有效区,如头发、支架等部分,并使其色素尽量减少,图2a 为经过格式转换的头部断层图像。在提取边界时,首先采用逐行扫描图片的办法,通过比较相邻点的像素值,找到图片边界上的一个点,作为切片边界的起点。然后从边界起点开始,逐点判断与之相邻的八个点,如果某点为图片的边界点则记录下,并开始下一步判断,直到获得所有的边界点。图2b 通过上述方法得到的面部边界轮廓曲线。通过这种方法得到的轮廓曲线是以大量的点坐标形式存储,在保证拟合精度的情况下,希望精简数据点,以

11、减少存储空间和提高计算机处理速度。因此将图2b 中的面部边缘轮廓曲线以中心点为极点转换到极坐标系中。然后,将图像画入以极角0,2ð为横坐标,极径为纵坐标的直角坐标系中,如图3a 所示。采用12阶傅立叶级数对该图像进行拟合,得到拟合后的图像如图3b 所示。 图2面部边界提取(a原始CT 图像(b面部边界轮廓y0 50 100 150 200 250 300 350 400 75808590 95100 105 110 115(b 拟合后的边界曲线 (a 拟合前的边界曲线2期 曾 筝等:利用MA TLAB 实现CT 断层图像的三维重建 27在这里,面部模型可以采用以下数学表示:=+=12

12、1,0,sin cos (2(n n i n i i i M n b n a a 式中n i a ,n i b ,为各层边界曲线傅立叶级数拟合的系数:d n f a i n i =cos (1,d n f b i n i sin (1,=。这样,某一层断层面部边界曲线仅使用25个数据就可以表示,所占用的存储空间很小。 下面,用这些表征面部边界轮廓曲线的数据进行三维表面重建。 (1 重建数据的采集。运用上述傅立叶级数的系数,求出边界上若干个点x ,y 向坐标值,并为其加上适当的z 坐标值。xo=0: pi/180:2*pi; % x 的值在0,2中选取yo=yo+a(i*cos(i-1*xo+b

13、(i*sin(i-1*xo; % 通过傅立叶系数求y 值,其中yo 初始值为a 0 consx=consx;yo.*cos(xo; %将x ,y 值从极坐标系转换到直角坐标系 consy=consy;yo.*sin(xo;consz=consz;ones(1,length(xo*iLayer*(-4.0; %为每一切片层赋予z 坐标值,iLayer 为层数(2 边界轮廓曲线表面绘制。surf(consx,consy,consz; %利用surf( 函数进行三维表面绘制。(3 设置图像的颜色及阴影效果。colormap(gray; %利用colormap( 函数为图像定义颜色集 shading

14、flat; %利用shading 定义显示图像的颜色阴影 (4 设置图像光照效果。light('Position',-80,-262,-200,'style','infinite' %利用light( 函数为图像设置光照效果 light('Position',-500,-0,-4500,'style','infinite' light('Position',5000,100,-300,'style','infinite'(5 设置图像的显示效果view

15、(-144,20; %利用view( 函数定义观察者视角 lighting gouraud; %利用lighting 定义显示图像的光线阴影 axis equal; %利用axis 定义显示图像的轴图4 头部CT 图像三维表面重建(a view(-144,20(b view(3(c view(-70,3028 CT理论与应用研究 13卷当函数view(AZ,EL取不同的值时,可以得到图像不同的视角,其中AZ表示图像水平方向旋转的角度,EL表示图像垂直方向的高度,MATLAB6.5自行对view(2和view(3的方位进行了定义。图4a为view(-144,20得到的图像;图4b为view(3得

16、到的图像;图4c为view(-70,30得到的图像。采用view( 函数可以从不同视角对三维重构图像进行观察。根据上述表面重建步骤,对头部切片进行表面重建,结果如图4所示。由于采用上述的傅立叶级数拟合含有耳朵的CT图像层时,耳朵的轮廓信息无法很好地表达,所以对这些层进行处理时,首先去除耳朵。从图4中可以看出,运用MATLAB程序在进行CT图像边界轮廓提取的基础上得到三维表面重建图像。重建速度快、效果好;但是面绘制的缺点是信息的丢失比较大,运算量与景物和物体形状有关。3CT图像三维体重建体绘制通过计算所有体素对光线的作用得到二维投影图像,基于体绘制的三维体重建方法计算量不依赖于景物的复杂程度和物

17、体形状的复杂程度,也不需要对切片的边界轮廓进行提取,其计算过程不依赖于视点,处理三维采样信号方便,便于显示物体的内部结构。但是,三维体重建所需数据量大,运算速度较慢。下面,介绍利用MATLAB进行头部CT图像三维体重建的过程。(1 重建数据的采集。对现有的n幅头部CT图像数据进行三维数据集D的构造,得到的数据集D为一个x×y×n的矩阵。image1 = imread('01.bmp' %使用imread( 函数读入现有的n幅图像image2 = imread('02.bmp'imagen = imread('n.bmp'D=c

18、at(3,image1,image2,image3,imagen; %使用cat( 函数创建三维矩阵D(2 重建数据预处理采用上述方法构造的三维数据集D,数据量大,在体重建中速度慢,并且可能在计算中超出内存。因而,可以根据实际情况,对数据集D进行预处理,减少数据量。x y z D = reducevolume(D,a b c; %使用reducevolume( 函数减少数据量,其中a,b,c为x,y和z轴数据抽取的比例,根据数据情况自行定义。D = smooth3(D; %使用smooth( 函数对数据进行平滑处理(3 计算数据集在显示平面累计投影。fv = isosurface(x,y,z,

19、D,isovalue; %使用isosurface( 函数计算数据集在显示平面累计投影,isovalue根据实际情况自行定义(4 构造三维体重建碎片p = patch(fv, FaceColor', 'yellow', 'EdgeColor', 'none' %使用patch( 函数对碎片进行构造,并对图像的颜色,光线进行定义,其中fv是第(3步中得到的。(5 设置图像的颜色、阴影及显示效果。2期 曾 筝等：利用 MATLAB 实现 CT 断层图像的三维重建 colormap(gray; 29 %利用 colormap( 函数为图像定义颜

20、色集 view(3; %利用 view( 函数定义观察者视角 lighting gouraud; %利用 lighting 定义显示图像的光线阴影 axis equal; %利用 axis 定义显示图像的轴 daspect(x y z; %使用 daspect( 定义 x、y、z 轴的显示比例 根据上述体重建步骤， 对头部切片 进行体重建，结果如图 5 所示。 从图 5 中可以看出，运用 MATLAB 程序在构造体重建碎片的基础上，实现体三维重建。 重建速度较表面重建慢，但是体重建尽量保护了更多的表面信息，耳朵部分也不需要任何特 殊处理，就能得到很好的重建及显示效果。 图5 头部 CT 图像三维体重建 结论 本文介绍了运用 MATLAB6.5 图像处理工具箱进行 CT 断层图像的三维重建的原理及实 现方法。并结合实际应用，实现了头部 CT 图像的三维表面重建及体重建，编程实现简单， 显示效果理想，大大提高了实验效率，