利用Matlab实现原木CT断层图像的三维重建

利用Matlab实现原木CT断层图像的三维重建

张汝楠,孙丽萍.

（东北林业大学，黑龙江.哈尔滨.150040）

摘要：运用MATLAB7.0软件中的图象处理工具箱实现了原木CT 断层图像的三维表面重建及体重建(构)，原理简单，编程方便，重建速度快，显示效果良好。关键词：MATLAB ；CT 图像；表面重建；体重建（构）中图分类号：S781.1 TP391.41 文献标识码：A

文章编号：1001-036X(2008)04-0024-03

Three dimensions reconstruction of log CT image by MATLAB

ZHANG Ru-nan,SUN Li-ping

(Northeast Forestry Vnirersity,Harbin 150040,China)

Abstract: In order to solve how to reconstruct CT image using image processing toolbox of MATLAB7.0 is discussed. MABLAB gives us a convenient method to reconstruct and acquire a good result.Key words: MABLAB ；CT image ；surface reconstruction ；volume reconstruction

原木分等一直是林业生产和木材加工的重要环节，长期以来都是以人工检测为主，受人为因素影响较大，效率较低，难以形成自动化生产线。

在计算机科学技术的发展，以及计算机断层扫描（CT）、核磁共振扫描（MIT）等技术的基础上，本文研究了原木CT断层图像的三维重建，在应用计算机断层扫描（CT）技术对原木进行扫描、并获取原木的二维断层图像的基础上，解决了二维断层图像不能直观、准确表现原木整体情况的问题。三维重建的步骤如图1。

1 原木CT图像预处理

为了有利于从CT图像中准确地提取有用的信息，需要对原始图像进行预处理，以突出有效的图像信息，消除或减少噪声的干扰。1.1 CT 图像格式的转换与读写

原始的CT图像是采用DICOM3.0标准进行存储的，不能被MATLAB所识别，因此必须进行图像格式的转换，转换后的图像为256色BMP图像[1]。

MATLAB实现图像格式的转换程序如下：info=dicominfo(‘2’);dicomread(‘2’);j=dicomread(‘2’);

imwrite(j,‘222.bmp’,‘bmp’);1.2 图像增强

图像增强是将图像中感兴趣的特征有选择的突出，并衰减不需要的特征,从而改善图像的视觉效果，便于人与计算机的分析和处理。

图像增强主要包括直方图修改、图像平滑、图像边

收稿日期:2008-04-05

作者简介:张汝楠(1982-),女,硕士研究生。研究方向：信息检测与

处理技术。

图1 原木CT图像三维重建流程图

缘锐化和彩色增强等。

直方图灰度变换法：照片或电子方法得到的图像，常表现出低对比度即整个图像偏亮或偏暗，为此需要对图像中的每一个像素的灰度级进行标度变换，扩大图像灰度范围，以达到改善图像质量的目的。这一灰度调整过程可以用 imadjust()函数实现。

直方图均衡化是一种使输出图像直方图近似服从均匀分布的变换算法。均匀量化的自然图像的灰度直方图通常灰度的比例很大，经直方图调整后，各灰度等级的比例更加平衡，原图像灰度集中的区域拉开从而增大反差，使原图像的细节成分更加清楚。在MATLAB7.0中直方图均衡化通过histeq()函数实现。图2是原木CT 断层图象及均衡化后的图象。

种根据图像的局部方差来调整滤波器输出的线性降噪滤波器。当局部方差大时，滤波器的效果较弱，反之滤波器的效果较强。在MATLAB7.0中用wiener2()函数实现[2]，

如图5。

图2 原木CT断层图象及均衡化后的图象

图3 对比度自适应直方图均衡化后的图象

图4 原木CT断层图像及中值滤波后的图像

?? ??? 图5 原木CT断层图像及自适应滤波后的图像

对比度自适应直方图均衡化：对图像的某个部分均衡化可用adapthisteq()函数。histeq()函数作用整幅图像时，可用adapthisteq()函数对图像的一个小区域进行操作，使其对比度得到加强，再与指定的直方图相匹配，如图3。

图像平滑用于去除图像噪声。基本采用在空间域上的求平均值或中值，或在频域上采取低通滤波。在MATLAB中，各种滤波方法都是在空间域中通过不同的卷积模板即滤波算子实现，可用fspecial()函数创建预定义的滤波算子，然后用filter2()或conv2()函数在实现卷积运算的基础上进行滤波。中值滤波是一种典型的低通滤波器，它在保护图像边缘的同时，能除去图像中的孤立点、线的噪声。在MATLAB7.0中用medfilt2()函数实现中值滤波(如图4)。自适应滤波是一

图像锐化是为了突出图像的边缘信息，加强图像的轮廓特征，便于人眼观察和计算机识别。可将锐化分为线性锐化和非线性退化锐化，包括微分算子、拉普拉斯

算子、Sobel算子和prewitt算子等[3]。

适当运用上述方法对原始图像进行处理，可使原始图像变得较清晰，能够较真实地反映图像的结构特征，便于三维重建的处理及显示。　

2 原木CT图像表面重建

2.1 重建数据的采集[5]

运用上述傅立叶级数的系数，求出边界上若干个点 x，y向的坐标值，并为其加上适当的z坐标值。

X 0=[0: pi/180:2*pi]; %x的值在[0,2π]中选取

Y 0=y 0+a(i)*cos((i-1)*x 0)+b(i)*sin((i-1)*x 0); %通过傅立叶系数求y值，其中y 0初始值为a 0consx=[consx;y 0*cos(x 0)];

%将x，y值从极坐标系转换到直角坐标系consy=[consy ;y 0*sin(x 0)];

consz=[consz;ones(1,length(x 0))*iLayer*(-4.0)]; %为每一切片层赋予z坐标值，iLayer为层数2.2 边界轮廓曲线表面绘制

Surf （consx,consy ,consz);

%利用surf()函数进行三维表面绘制

2.3 设置图像的颜色及阴影效果

colormap(gray);

%利用colormap()函数为图像定义颜色集shading flat;

%利用 shading 定义显示图像的颜色阴影

2.4 设置图像光照效果

light(‘Position’,[-80,-262,-200],’style’, ’infinite’);

%利用light()函数为图像设置光照效果

light(‘Position’,[-500,-0,-4500],’style’, ’infinite’);

light(‘Position’,[5000,100,-300],’style’, ’infinite’);

2.5 设置图像的显示效果

view(-144,20);

%利用view()函数定义观察者视角

lighting gouraud; %利用lighting定义显示图像的光线阴影axis equal;

%利用axis定义显示图像的轴

当函数view(AZ,EL)取不同的值时，可以得到图像不同的视角，其中 AZ 表示图像水平方向旋转的角度，EL 表示图像垂直方向的高度。

3 原木CT图像体重构

体绘制法是将三维空间的离散数据直接转换为二维图像而不必生成中间几何图元，通过计算所有体素对光线的作用得到二维投影图像。基于体绘制的三维体重建方法计算量不依赖于景物的复杂程度和物体形状的复杂程度，也不需要对切片的边界轮廓进行提取，其计算过程不依赖于视点，处理三维采样信号方便，便于显示物体的内部结构，但是，三维体重建所需数据量大，运算速度较慢。

3.1 体重构数据的采集

对现有的n幅原木CT图像数据进行三维数据集D的构造，得到的数据集D为一个x×y×n的矩阵。

image1 = imread(‘1.bmp’);

%使用imread()函数读入现有的n幅图像

image2 = imread(‘2.bmp’);

：

：

imagen = imread(‘n.bmp’);

D＝cat(3,image1,image2,image3,… imagen);

%使用cat()函数创建三维矩阵D

3.2 体重构数据预处理

采用上述方法构造的三维数据集D，数据量大，在体重建中速度慢，并且可能在计算中超出内存。因而,可以根据实际情况，对数据集D进行预处理，减少数据量。

[x y z D] = reducevolume(D,[a b c]);

%使用reducevolume()函数减少数据量，其中

a,b,c为x,y和z轴数据抽取的比例，根据数据

情况自行定义。

D = smooth3(D);

%使用smooth()函数对数据进行平滑处理

3.3 计算数据集在显示平面累计投影

fv = isosurface(x,y,z,D,isovalue);

%使用isosurface()函数计算数据集在显示平

面累计投影，isovalue 根据实际情况自行定义。 3.4 构造三维体重构碎片处理

p=patch(fv,‘FaceColor’,‘yellow’,‘EdgeColor’,‘none’);

%使用patch()函数对碎片进行构造，并对图像的 颜色、光线进行定义，其中fv是4.3中得到的。

3.5 设置图像的颜色、阴影及显示效果

colormap(gray);

%利用colormap()函数为图像定义颜色集

view(3);

%利用view()函数定义观察者视角

lighting gouraud;

%利用lighting 定义显示图像的光线阴影

axis equal;

%利用axis 定义显示图像的轴

daspect([x y z]);

%使用daspect()定义 x、y、z 轴的显示比例

4 结论

本文运用 MATLAB7.0软件中的图像处理工具箱，结合林业应用实际，实现了原木CT图像的三维表面重建及体重构，编程简单，显示效果理想,大大提高了实验

(下转第29页)

可变角探头作出入射角和反射回波振幅的关系曲线，各曲线上极大值所对应的入射角即为在该板中可能激励出的兰姆波入射角；从中又选出板端面和人工缺陷反射回波前沿陡削、无杂波，板边盲区和探头前沿盲区较小的入射角和频率作为可供选择的检测入射角和频率。

当入射角和频率选定以后，可以调整超声换能器发出符合要求的超声波，此时激励出来的兰姆波的群速度等于相速度，可以用于单板顺纹弹性模量的计算。

图1所示为其检测结构示意图：频率进行限定，才能使兰姆波的群速度等于相速度；另一方面，入射角是前提，没有合适的入射角就没有兰姆波的形成。没有兰姆波的形成，我们通过测试传播时间测出的声速值，必然与通过公式Cp=C

L

/sina计算出的声速值不吻合。

(2）入射角和频率的确定，我们是在有机玻璃斜楔块上用不同频率的可变角探头作出入射角和反射回波振幅的关系曲线图，之后在各曲线上找出极大值所对应的入射角，即为在该板中可能激励出的兰姆波入射角。

(3）如何发射并接收兰姆波、入射角度、楔形材料的选择及其能否在单板顺纹方向顺利传播是决定于兰姆波能否在单板顺纹弹性模量检测中应用的关键问题。

(4）采用兰姆波检测单板顺纹弹性模量，目前看来还是比较前沿的，之前由于多用于金属薄板的检测而没有人做过单板方面的研究，因此还面临很多困难和问题，诸如耦合剂的选择、单板的缺陷和含水率的影响等，需要我们在实验中进一步探讨与解决。

[参考文献]

[1]Forest Products Laboratory.Wood handbook wood-as an engineering material [M].USDA Forest Service,1999.

[2] 卢晓宁.速生杨木单板顺纹弹性模量预测模型[J].南京林业大学学报,2002,26(3):9-13.

[3] 林莉.超声无损表征薄层结构研究进展[J].无损探伤,2006,30(5): 1-4.

[4] 荆峰.铝薄板的蓝姆波检测试验[J].航天制造技术,2003,(5):42-46.

[5] 周乐文.金属薄板的超声波检测[J].云南大学学报,2002,24(1A): 249-250.

图1 单板超声检测弹性模量结构图

此外我们也可以利用此图所示的结构，通过超声发生/接收器，以及数字示波器，测试出兰姆波在单板中的传播时间，从而求出兰姆波在单板中的传播速度（单板长度/传播时间）；之后与公式Cp=C

L

/sina计算得出的兰姆波相速度进行比对，从而对结果进行验证。如果相差太大，则需要进一步选择入射角和频率，以及斜楔块的材料，直到两种方法测出的兰姆波的传播速度相近为止。

3 初步结论

(1）兰姆波的激励问题，归根结底就是入射角和频率的选择，只有选择了合适的频率，从而对超声发生器

效率，在原木无损检测中具有一定的使用价值和实际指导意义。在图像三维重建的其他领域内，对于类似的图像，都可以运用上述方法来实现三维重建，以获得重建对象的三维表面信息及体信息。　

[参考文献]

[1]王成波,陈伟,谢兵.DICOM图像与BMP图像的转换研究[J].医疗

卫生装备,2004,(1):13-17.

[2]罗军辉,冯平,哈力旦·A.MATLAB7.0在图像处理中的应用[M].北京:机械工业出版社,2005.

[3]贺兴华,周媛媛,王继阳,周晖.M A T L A B7.x图像处理[M].北京:人民邮电出版社,2006.

[4]曾筝,董芳华,陈晓,周宏,周建中.利用M A T L A B实现C T断层图像的三维重建[J].CT理论与应用研究,2004,(5):24-29.

(上接第26页)