文章信息
- 罗瑞, 葛浙东, 陈龙现, 刘传泽, 周玉成, 徐文琦.
- Luo Rui, Ge Zhedong, Chen Longxian, Liu Chuanze, Zhou Yucheng, Xu Wenqi.
- 基于X射线木材断层扫描反投影算法中的最优滤波器选取
- Selection of the Optimum Filter in the Back Projection Algorithm Based on X-Ray Wood Tomography
- 林业科学, 2018, 54(11): 143-148.
- Scientia Silvae Sinicae, 2018, 54(11): 143-148.
- DOI: 10.11707/j.1001-7488.20181120
-
文章历史
- 收稿日期:2018-04-03
- 修回日期:2018-07-23
-
作者相关文章
CT(computed tomography),即计算机断层扫描,20世纪60年代首先应用于医学领域,具有穿透力强、分辨率高、检测速度快、检测结果直观且无需破坏被检测物等特点,是一种先进的非接触式检测技术。近年来,随着科学技术快速发展,CT在农业、林业、建筑等领域也显示出其独特的优越性(余晓锷等,2013),已有学者将CT技术应用到木材内部结构无损探伤研究中,通过重建断层图像,实现了对木材内部缺陷位置、大小和形状的检测(葛浙东等,2016a; 2016b;戚玉涵等,2016)。断层图像重建方法主要有解析法(Liu et al., 2006)和迭代法(Donath et al., 2006)两大类,而解析法中的滤波反投影算法应用最为广泛(闫镔等,2014),基于该算法可实现图像快速重建并获得高质量重建图像。在X射线计算机断层扫描图像重建中,滤波器是图像处理最重要也是最关键的环节,图像处理时由于存在图像数据重复、缺失等问题以及在图像二维傅里叶变换中的Gibbs现象(胡君杰等,2013),X射线断层成像过程通常携带噪声,并且存在图像不连续问题,因此滤波反投影算法的关键是对滤波函数的选取(黄亚等,2012)。骆岩红(2014)提出将指数滤波函数引入图像重建过程,研究不同α值对重建图像质量的影响,该方法虽然对星状伪影具有较好的滤波效果,但是却降低了重建图像的空间分辨率。张斌(2009)、翟静(2008)通过赋予R-L、S-L滤波函数不同权重对滤波函数进行折中处理,提出一种混合滤波器,在一定程度上提高了重建图像的空间分辨率和密度分辨率,改善了重建图像质量,但该方法运算过程复杂且效率较低。石本义等(2010)从加权平均思想出发,设计一种兼顾密度分辨率和空间分辨率的新滤波函数,保证了重建图像整体性能,但其重建图像的空间分辨率不及R-L滤波函数,密度分辨率不及S-L滤波函数。Pelt等(2014)提出最小化重建图像的投影与测量投影的平方差来计算新滤波器, 并通过求解最小二乘意义上的线性系统来寻找最优滤波器,该方法可最大限度减少投影误差并有效提高重建效率,但却很难保证重建图像的质量。目前,国内外滤波器共有60多种,在X射线计算机断层扫描图像重建算法中常用的有5种,即R-L、S-L、Cosine、Hamming和Hanning滤波器。对于不同扫描对象,不同滤波器的滤波效果是不同的,多数滤波函数设计通过Sheep-Logan模型仿真测试检验(张东平等,2007);但在木材断层扫描成像中,由于木材内部结构复杂,密度分布不均,其CT成像过程普遍存在大量噪声,重建图像存在严重伪影、伪迹。
鉴于此,本研究在木材X射线计算机断层扫描图像重建中,通过试验比较几种滤波器的重建图像效果,分析滤波函数重建图像质量,对反投影算法中的最优滤波器选取问题进行探讨,旨在获取更高质量的木材断层图像,为木材断层成像技术研究提供科学方法。
1 材料与方法 1.1 试验材料杉木(Cunninghamia lanceolata)木段取自浙江省杭州市,直径221 mm,高415 mm,去皮,外观笔直,表面有大小不一的裂缝沿树干径向分布,并存在节子及不同程度的腐朽。
1.2 试验设备试验平台为课题组自行开发的CT成像设备,如图 1所示,包括X射线发射器、平板接收器和支架车。X射线发射器(IXS160型),英国VJ TECHNOLOGIES公司生产,输入电压220 V,最高连续输出功率500 W,射线管工作电压在10~160 kV范围内连续可调,本研究选择120 kV。平板接收器(DS_LINX V3型,由硅光二极管阵列和闪烁晶体组成),英国SENS-TECH公司生产,长778 mm、宽201 mm,有效检测区域512 mm,像元尺寸0.4 mm× 0.4 mm,由此平板接收器线性等距分布1 280个探测点。支架车长2 800 mm、宽704 mm、高955 mm,设置转台,直径500 mm。木段置于转台中心随转台匀速旋转。平板接收器距转台中心365 mm,X射线发射器距转台中心410 mm。
CT成像设备围绕木段旋转1周,可得3 600行、1 280列原始投影数据,行坐标表示旋转角度,列坐标表示探测点编号。平板接收器存在损坏探测点,需对原始投影数据进行预处理。判断坏点位置采用最小值法,即认为原始投影数据中各行最小值所在列为损坏探测点位置。采集空扫数据,对空扫数据各列取平均值,填补于原始投影数据中损坏探测点所在列。对原始投影数据进行预处理后,选择R-L、S-L、Cosine、Hamming和Hanning滤波器分别与预处理后的数据进行卷积,对投影数据进行滤波运算。由于空间域中卷积运算复杂、耗时长,因此将卷积运算变换到频域作乘积运算,先将经傅里叶变换到频域的投影数据与滤波函数频域离散形式作乘积,再将频域运算结果作傅里叶逆变换得到空间域滤波后的投影数据,其结果与卷积运算等效,该算法有可效减少运算量,提高运算效率。卷积运算后,对滤波后的投影数据进行反投影重建,即得到木段断层扫描图像。
1.4 滤波器在木材断层扫描过程中,投影数据由于受环境、设备和试验材料的影响,存在噪声、振荡等干扰。为最大程度滤除干扰数据,呈现清晰、伪影少的木材断层重建图像,滤波器的选取至关重要。本研究选取R-L、S-L、Cosine、Hamming和Hanning滤波器,对5种滤波器的工作特性和滤波效果进行推理预测。进一步,分别采用5种滤波器对木段作反投影重建,获取不同滤波器下的木段断层重建图像,并对其滤波效果进行评价与分析。
R-L滤波器计算公式如下:
$ {H_{{\rm{R - L}}}} = |\rho |W\left(\rho \right) = |\rho |{\rm{rect}}\left({\rho /2B} \right)。$ | (1) |
该滤波器采用矩形窗rect(ρ/2B),函数结构简单、易于实现,矩形窗实时带宽较大,最大栅栏效应误差为3.92 dB;但由于矩形窗主瓣宽度大,过渡带相对较长,最大旁瓣相对主瓣大,所以平稳度相对较差,存在明显的Gibbs现象,即振荡响应。当重建图像过程中存在噪声时,重建图像质量相对较差。
S-L滤波器计算公式如下:
$ {H_{{\rm{S - L}}}}\left(\rho \right) = |\rho |{\rm{sinc}}\left({\frac{\rho }{{2B}}} \right){\rm{rect}}\left({\frac{\rho }{{2B}}} \right)。$ | (2) |
该滤波器在R-L滤波器的基础上引入了sinc窗函数,有效减轻了R-L滤波的振荡响应,降低了|ρ|=1/2d处的幅值,从而可更好补偿|ρ|=1/2d处的混迭现象;但由于HS-L(ρ)在低频段偏离了|ρ|,所以重建图像质量不及R-L滤波器。
Cosine滤波器计算公式如下:
$ {H_{{\rm{Cosine}}}} = {\rm{cos}}(\frac{{{\rm{ \mathsf{ π} }}\rho }}{{2B}}){\rm{rect}}\left({\frac{\rho }{{2B}}} \right)。$ | (3) |
该滤波器取理想滤波器的余弦形式并引入矩形窗,将高频段信号全部滤除,对中低频段也部分舍去,可以更精确地重建图像;但由于杉木断层成像过程中舍去了部分有效数据,所以导致重建图像存在瑕疵。
Hamming滤波器计算公式如下:
$ {H_{{\rm{Hamming}}}} = 0.54 + 0.46{\rm{cos}}(\frac{{{\rm{ \mathsf{ π} }}\rho }}{B}){\rm{rect}}\left({\frac{\rho }{{2B}}} \right)。$ | (4) |
该滤波器是加权系数为0.54的余弦窗函数,其加权系数使得旁瓣更小,从而可达到更高的幅值识别精度。
Hanning滤波器计算公式如下:
$ {H_{{\rm{Hamming}}}} = 0.5 + 0.5{\rm{cos}}(\frac{{{\rm{ \mathsf{ π} }}\rho }}{B}){\rm{rect}}\left({\frac{\rho }{{2B}}} \right)。$ | (5) |
该滤波器是加权系数为0.5的余弦窗函数,且是一个平滑的时间函数,旁瓣之间可以相互抵消,使得主瓣作用加强,从而可降低矩形窗的不连续性,减小泄露,提高频谱幅值精度。
2 结果与分析分别获取杉木木段投影数据和空扫数据,对投影数据和空扫数据进行预处理后,不采取滤波处理,直接反投影重建,成像结果如图 2a所示,该断层重建图像基本呈现木段生长轮、裂纹等信息,但图像存在大量噪声,生长轮界限不明显,裂缝缝隙特征被噪声严重干扰。
分别采用R-L、S-L、Cosine、Hamming和Hanning滤波器对木段作反投影重建。采用R-L滤波器获得的木段反投影重建图像(图 2b),滤波后可有效去除伪影、伪迹,明显减少噪声干扰,木段断层生长轮清晰,4处裂缝基本表征,但数据预处理过程中均值计算造成的无效数据点仍然存在,裂缝结构特征不够清晰。采用S-L滤波器获得的木段反投影重建图像(图 2c),与R-L滤波后的重建图像相比,生长轮特征不清晰。采用Hanning滤波器获得的木段反投影重建图像(图 2d),呈现木段断层生长轮和4处裂缝,有效减少了木段断层图像周围的干扰点,相比图 2b,生长轮边界模糊。采用Hamming滤波器获得的木段反投影重建图像(图 2e),木段断层图像周围的干扰点有效减少,但与图 2b、c、d、f相比,裂缝结构和生长轮清晰度最差,图像质量最差。
采用滤波器进行反投影重建可以不同程度提高重建图像质量,为进一步分析滤波后的重建图像质量,引入归一化均方距离d、归一化平均绝对距离r比较图 2中5幅滤波后断层图像之间差异:
$ d = \frac{{\sum\limits_{u = 1}^M {} \sum\limits_{v = 1}^N {} {{({t_{u, v}} - {r_{u, v}})}^2}}}{{2\sum\limits_{u = 1}^M {} \sum\limits_{v = 1}^N {} {{({t_{u, v}} - \bar t)}^2}}}; $ | (6) |
$ r = \frac{{\sum\limits_{u = 1}^M {} \sum\limits_{v = 1}^N {} |{t_{u, v}} - {r_{u, v}}|}}{{\sum\limits_{u = 1}^M {} \sum\limits_{v = 1}^N {} |{t_{u, v}}|}}。$ | (7) |
式中:tu, v为检测模型图像第u行、第v列的像素密度;ru, v为重建断层图像第u行、第v列的像素密度;
归一化均方距离d越大,表示检测模型图像与重建断层图像的误差越大,当d=0时,重建图像能够真实反映检测模型图像,当出现个别较大偏差时,d偏大,因此d可准确判断少数点出现较大误差的情况。归一化平均绝对距离r越大,表示误差越大,当r=0时,表示无误差,因此r可反映小误差普遍存在的情况。
选择高质量重建图像(图 3)作为检测模型,与图 2b、c、d、e、f进行对比,结果如表 1所示。
由表 1可知,采用不同滤波器对同一组投影数据作滤波处理,重建断层图像质量差距很大。经R-L滤波器滤波后的数据,断层重建图像的d和r分别为0.075 7和0.032 6,将图 2b与检测模型图像进行对比,2幅图像相似度最高,图像质量最好。经Hamming滤波器滤波后的数据,断层重建图像的d和r分别为0.105 1和0.043 6,将图 2e与检测模型图像进行对比,2幅图像相似度最低,图像质量最差。
4 结论对比R-L、S-L、Cosine、Hamming和Hanning 5种滤波器的木材断层重建图像,分析滤波函数重建图像质量,可以发现,采用R-L滤波器得到的重建图像误差最小,图像最清晰,空间分辨率高,能有效去除伪影,滤波后的木材断层图像可准确判断缺陷形状、位置和大小,满足各科研院所、高等院校、检测机构等对木材内部缺陷检测和分析的要求。滤波函数的选择对重建图像质量影响非常大,因此多种滤波器的比较和分析将是下一步研究的重点内容。
葛浙东, 侯晓鹏, 鲁守银, 等. 2016a. 基于反投影坐标快速算法的木材CT检测系统研究. 农业机械学报, 47(3): 335-341, 327. (Ge Z D, Hou X P, Lu S Y, et al. 2016a. Wood CT detection system based on quick algorithm of back projection coordinates. Transactions of the Chinese Society for Agricultural Machinery, 47(3): 335-341, 327. [in Chinese]) |
葛浙东, 侯晓鹏, 李早芳, 等. 2016b. 断层扫描技术在木材无损检测中的应用. 木材工业, 30(3): 49-52. (Ge Z D, Hou X P, Li Z F, et al. 2016b. Application of tomography in non-destructive testing of wood. China Wood Industry, 30(3): 49-52. [in Chinese]) |
胡君杰, 马晨欣, 闫镔. 2013. CT图像重建中滤波函数的优化. CT理论与应用研究, 22(1): 85-92. (Hu J J, Ma C X, Yan B. 2013. Optimization of filtering function in CT image reconstruction. CT Theory and Applications, 22(1): 85-92. [in Chinese]) |
黄亚, 张祥志, 祝江威, 等. 2012. CT断层图像重建的新滤波函数. 核电子学与探测技术, 32(12): 1388-1393. (Huang Y, Zhang X Z, Zhu J W, et al. 2012. New filter function for reconstruction of CT tomographic images. Nuclear Electronics & Detection Technology, 32(12): 1388-1393. DOI:10.3969/j.issn.0258-0934.2012.12.011 [in Chinese]) |
骆岩红. 2014. CT图像重建滤波反投影算法中指数滤波器的研究. 计算机科学, 41(C1): 220-223. (Luo Y H. 2014. Research on exponential filter in CT image reconstruction filtering back projection algorithm. Computer Science, 41(C1): 220-223. [in Chinese]) |
戚玉涵, 徐佳鹤, 张星梅, 等. 2016. 基于扇形X射线束的立木CT成像系统. 林业科学, 52(7): 121-128. (Qi Y H, Xu J H, Zhang X M, et al. 2016. Tiling CT imaging system based on fan X-ray beam. Scientia Silvae Sinicae, 52(7): 121-128. [in Chinese]) |
石本义. 2010.基于GPU的计算机断层成像技术研究.武汉: 华中科技大学硕士学位论文. (Shi B Y. 2010. Research on GPU-based computed tomography. Wuhan: MS thesis of Huazhong University of Science and Technology.[in Chinese]) http://cdmd.cnki.com.cn/Article/CDMD-10487-1012013280.htm |
闫镔, 李磊. 2014. CT图像重建算法. 北京: 科学出版社, 14-18. (Yan B, Li L. 2014. CT Image reconstruction algorithm. Beijing: Science Press, 14-18. [in Chinese]) |
余晓锷, 龚剑, 马建华, 等. 2013. CT原理与技术. 北京: 科学出版社, 95-97. (Yu X E, Gong J, Ma J H, et al. 2013. CT principles and techniques. Beijing: Science Press, 95-97. [in Chinese]) |
张斌. 2009.滤波反投影图像重建算法中插值和滤波器的研究.太原: 中北大学硕士学位论文. (Zhang B. 2009. Interpolation and filter research in filter back-projection image reconstruction algorithm. Taiyuan: MS thesis of North University of China.[in Chinese]) http://cdmd.cnki.com.cn/Article/CDMD-10110-2009170835.htm |
张东平, 张定华, 张丰收, 等. 2007. 滤波器对滤波反投影重建图像质量的影响. 机械设计与制造, (9): 66-68. (Zhang D P, Zhang D H, Zhang F S, et al. 2007. Influence of filter on filtered back projection reconstruction image quality. Mechanical Design & Manufacture, (9): 66-68. [in Chinese]) |
翟静. 2008.三维锥束CT中滤波反投影算法的研究.太原: 中北大学硕士学位论文. (Zhai J. 2008. Research on filtered back projection algorithm in 3D cone-beam CT. Taiyuan: MS thesis of North University of China.[in Chinese]) http://cdmd.cnki.com.cn/Article/CDMD-10110-2008142435.htm |
Donath T, Beckmann F, Schreyer A. 2006. Automated determination of the center of rotation in tomography data. Journal of the Optical Society of America A:Optics Image Science and Vision, 23(5): 1048-1057. DOI:10.1364/JOSAA.23.001048 |
Liu T, Malcolm A A. 2006. Comparison between four methods for central ray determination with wire phantom in micro-computed-tomography systems. Optical Engineering, 45(6): 711-725. |
Pelt D M, Batenburg K J. 2014. Improving filtered back projection reconstruction by data-dependent filtering. IEEE Transactions on Image Processing, 23(11): 4750-4762. DOI:10.1109/TIP.2014.2341971 |
Stephen G A, Daniel J S, et al. 1990. Calculation of the rotational centers in computed tomography sinogram. IEEE Trans on Nuclear Science, 37(4): 1525-1540. DOI:10.1109/23.55866 |