Research on Quadric Surface Fitting Algorithm Based on CT Data
-
摘要: 二次曲面工件在工业中比较常见，为测量物体内部的二次曲面，本文采用工业计算机断层成像（CT）技术获取物体的断层图像序列，利用U-net图像分割网络获得断层图像上的目标区域，对分割结果的边缘进行检测和曲线拟合并堆叠成三维点集，通过曲面拟合获取物体曲面的三维空间坐标信息。研究结果表明，本文的方法能够有效地实现边界提取和界面参数的拟合工作，拟合误差在1％ 以内，相比传统方法有较大改进。Abstract: Quadric surface is a common type of workpiece shape in industry. It can be imaged by computed tomography (CT). A sample object was scanned and its slices were reconstructed to measure its internal quadric surface. We used the U-net image segmentation network to obtain the interested region and then detected the edge in the segmentation results and obtained curve fitting results. The curves were stacked into a three-dimensional point set. The three-dimensional spatial coordinate information for the internal quadric surface was computed through surface fitting. The results demonstrate that the proposed method can effectively extract the internal quadric surface parameters, and the fitting error can be controlled within 1%, which is superior to that achieved with a traditional algorithm.
-
Keywords:
- ICT /
- deep learning /
- quadric surface /
- curving fitting
-
伦琴于1895年发现X射线，而后英国工程师Hounsfield基于X射线发明出计算机断层成像（computed tomography，CT）机。自CT机问世几十年来，无论是在硬件性能（X光机、探测器）、扫描方式，还是重建算法上都获得了长足发展，在医疗、工业和安全检查等领域发挥着重要作用。更短的扫描和重建时间，以及更好的成像质量成为该领域研究人员追求的目标[1]。
双能CT是由Alvarez[2]在1976年首次提出，可以获取物体的等效原子序数和电子密度信息，从而更加准确地识别物质。多能谱CT是在双能CT的基础上，增加了参与扫描和重建的能区数量，引入了更多能谱信息，重建结果更加准确。
随着高计数率电子元器件的迅速发展，以光子计数探测器为核心的多能谱X射线成像技术成为射线检测领域的重要发展方向[3-4]。与传统能量积分探测器相比，光子计数探测器可以通过设置能量阈值来选择不同能区的信号，同时它将每个光子当作独立事件进行统计，一次扫描便能得到多段能谱，减小了辐射剂量[5]。但在实际的应用上也存在多能段光子数减少、信噪比下降，探测器性能不足造成的噪声等问题[6]，对影响光子计数X射线CT成像精度的因素及其规律进行基础性研究十分必要。
本文对系统能谱的获取、模体的扫描和重建过程进行仿真实验，以重建得到的等效原子序数和电子密度的标准差作为成像质量的评价标准，得出探测器能量分辨率、能谱阈值和能谱范围等因素对成像精度影响的有关结论，为后续相关技术的工程应用化提供参考。
1. 双能重建基本原理
本文主要研究探测器能量分辨率、能谱阈值和能谱范围等对重建图像精度的影响规律。以双能重建为例，根据光子计数探测原理[7]，进行双能成像的模拟。双能CT成像的重点研究内容之一是双能分解方法，常见的双能分解方法包含投影域前处理、图像域后处理和迭代方法[8]。
考虑到未来的工程化需求，本文选择的双能重建方法为基于基材料模型的双能投影预处理分解方法，基材料选取的是碳和铝。
材料模型公式：
$$ \mu \left(E\right)={{b}}_{1}{\mu }_{1}\left(E\right)+{{b}}_{2}{\mu }_{2}\left(E\right) \text{，} $$ (1) $ {\mu }_{1} $ (E)、$ {\mu }_{2} $ (E)分别为两种基材料的线性衰减系数；$ {b}_{1}$ 、$ {b}_{2}$ 对应两种基材料的分解系数，对于某一固定的物质，$ {b}_{1}$ 、$ {b}_{2}$ 是两个常数。双能分解即求解如下的双能投影积分方程组：
$$ \left\{\begin{aligned} {P}_{L}=\;& -{\ln}\bigg(\int \mathit{{S}_{L}} (E)\exp \Big({-{B}}_{1}{\mu }_{1}({E})-\\& {{B}}_{2}{\mu }_{2}({E})\Big){\rm d}{E}\bigg)+{\ln}\int \mathit{{S}_{L}}(E){\rm d}E\\ {P}_{H}=\;&-{\ln}\bigg(\int \mathit{{S}_{H}}(E) \exp \Big(-{{B}}_{1}{\mu }_{1}({E})-\\ & {{B}}_{2}{\mu }_{2}({E})\Big){\rm d}{E}\bigg)+{\ln}\int \mathit{{S}_{H}}(E){\rm d}E\end{aligned}\right. \text{，} $$ (2) 其中，SH、SL分别为高低能能谱，纵坐标记录的是光子数，PH、PL分别为高低能投影，
$ {{B}}_{1}=\displaystyle \int {{b}}_{1}\mathrm{d}{l}， {{B}}_{2}=\displaystyle \int {{b}}_{2}\mathrm{d}{l} $ 。求解
$ {B}_{1}$ 、$ {B}_{2}$ 的这个过程被称为投影分解过程。由于$ {B}_{1}$ 、$ {B}_{2}$ 为$ {b}_{1}$ 、$ {b}_{2}$ 的线积分投影值，求解出$ {B}_{1}$ 、$ {B}_{2}$ 后，根据CT重建的原理，利用滤波反投影重建算法[9]，便可计算出$ {b}_{1}$ 和$ {b}_{2}$ ，由此可以计算材质的等效原子序数和电子密度信息，以完成材料的探测识别。计算公式：$$ {Z}_{\rm{eff}}={\left(\frac{{b}_{1}{\rho }_{e1}{Z}_{1}^{n}+{b}_{2}{\rho }_{e2}{Z}_{2}^{n}}{{b}_{1}{\rho }_{e1}+{b}_{2}{\rho }_{e2}}\right)}^{\tfrac{1}{n}} \text{，} $$ (3) $$ {\rho }_{e}={{b}_{1}\rho }_{e1}+{{b}_{2}\rho }_{e2} ， $$ (4) 式（3）和式（4）中Z1、Z2分别为两种基材料的原子序数，ρe1、ρe2分别为两种基材料的电子密度，n一般取在3~4之间。
2. 光子计数CT成像系统能谱仿真
在Ubuntu系统下利用Geant 4软件，模拟高速电子撞击靶材射出X射线的过程。真空管内靶材为1 cm厚的钨，靶面倾角为25°，射线源射出能量为160 keV的电子，滤波片为2.5 mm厚的铝，射线穿过区域皆为真空，光子穿过虚拟探测器后被记录能量，得到射线源入射能谱。图1为能谱发射的几何示意，图2为得到的入射能谱。
通过对入射能谱与基于能量分辨率高斯函数的卷积来模拟探测器能量分辨率的影响，得到系统能谱[10]。后续实验所需的高低能能谱均基于系统能谱得到。探测器能量分辨率R通常用半高全宽（full width at half maximum，FWHM）W表示，定义：
$$ R=\frac{W}{E}\text{，} $$ (5) 其中，E是相对应的能量。W可以通过标准差σ与高斯分布的关系计算得到：
$$ W=2.355\;\sigma 。 $$ (6) 因此，标准差可以表示为：
$$ \sigma =\frac{R\times E}{2.355} 。 $$ (7) 根据探测器指定能量下的标准能量分辨率R0和能量E0，可由式（8）获得其他能量下的能量分辨率。由式（7）获得每个能量点下的标准差σ，并使用该标准差对原始入射能谱进行高斯模糊[11]。
$$ R={R}_{0}\sqrt{\frac{{E}_{0}}{E}} 。 $$ (8) 3. 投影数据加噪仿真
由于光子计数的特性遵循泊松分布[12-13]，所以本文为投影数据添加泊松噪声，用于模拟真实环境。泊松分布的概率函数见式（9），该函数表示在给定时间间隔或空间区域内发生k次事件的概率，其中，
$ \lambda $ 是泊松分布的均值和方差。$$ P\Big(X=k\Big)=\frac{{\lambda }^{k}{\exp}{(-\lambda )}}{k!} 。 $$ (9) 具体加噪过程。
（1）光子穿过模体后的衰减规律：
$$ {n}=\int S\left(E\right){\exp}{\left(\int -l\mu \left(E\right)\mathrm{d}l\right)}\mathrm{d}E \text{，} $$ (10) 其中，n为入射后探测器接收到的光子数，S(E)为系统能谱，exp为自然常数，l表示射线穿过模体的厚度，μ(E)为水模的线性衰减系数。
将各个单能（E1、E2、E3、···）下的光子数ni=
$ S\left({E}_{\mathrm{i}}\right){\exp}{\left(\int -l\mu \left({E}_{\mathrm{i}}\right)\mathrm{d}l\right)} $ 分别代入式（9）中的$ \lambda $ ，将对应泊松分布下生成的随机数作为该能量下加噪后的光子数ni_noise。（2）将能谱各单能段加噪后的光子数ni_noise累加，并与入射前探测器接收到的光子数n0按下式进行变换，便可得到加噪后的投影Pnoise：
$$ {P}_{\mathrm{n}\mathrm{o}\mathrm{i}\mathrm{s}\mathrm{e}}=\mathrm{l}\mathrm{n}\left(\frac{{n}_{0}}{{n}_{\mathrm{n}\mathrm{o}\mathrm{i}\mathrm{s}\mathrm{e}}}\right) \text{，} $$ (11) 其中，入射前探测器接收到的光子数
$ {n}_{0}= \displaystyle \int S\left(E\right)\mathrm{d}E $ 。4. 实验仿真和结果分析
4.1 实验模体和成像几何
设置的仿真模体为一圆柱状物体，材料为钙（原子序数20，密度1.55 g/cm3），半径为3 cm，模体截面如图3所示。
进行扫描仿真实验的探测系统几何结构如图4所示。射线源为扇形束，扇束张角为22.6°；探测器线性等距排列，共256个探测单元，总长512 mm；射线源、模体、探测器中心位于同一条直线；射线源到模体中心和探测器中心的距离分别为640 mm和
1280 mm；每隔1° 扫描一次模体，每次扫描发射的光子数为105。4.2 能量分辨率的影响
市场上性能相对较好的光子计数探测器材料为碲锌镉（CZT），其在60~100 keV下的能量分辨率大约是10%。因此参考碲锌镉材料的分辨率作为最优值，从10% @ 100 keV到40% @ 100 keV依次递减能量分辨率，间隔10%取4组能量分辨率进行实验，并配以理想能量分辨率（0%）作为对比，高低能能量阈值选为80 keV，各实验组参数见表1。各实验组高斯模糊化之后的高低能能谱见图5。
表 1 能量分辨率实验对照表Table 1. Experimental comparison table of energy resolution编号 能量分辨率 低能区光子数 高能区光子数 1 0% @ 100 keV 82978 17286 2 10% @ 100 keV 82964 17271 3 20% @ 100 keV 82964 17218 4 30% @ 100 keV 82963 17132 5 40% @ 100 keV 82932 17005 穿过模体的光子由所设能量阈值被划入不同的能量通道，从而得到高低能投影数据。为投影数据添加泊松噪声后，进行双能重建[14-15]。
上述实验组的成像结果见表2。不同能量分辨率下的重建结果标准差曲线如图6所示。
表 2 能量分辨率实验成像结果表Table 2. Imaging results of energy resolution experiment编
号能量
分辨率等效原子
序数等效原子
序数标准差电子密度 电子密
度标准差1 0% @
100 keV19.538647 0.268186 1.518764 0.039321 2 10% @
100 keV19.549584 0.291602 1.517544 0.042370 3 20% @
100 keV19.544331 0.295689 1.518037 0.042150 4 30% @
100 keV19.544904 0.313434 1.518503 0.044175 5 40% @
100 keV19.570808 0.350212 1.515818 0.048496 ![]()
图 6 不同能量分辨率下的成像结果对比Figure 6. Comparison of imaging results at different energy resolutions上述结果表明，随着探测器能量分辨率的提高，成像结果受噪声影响程度下降，等效原子序数相较电子密度而言，受噪声影响程度更大。不难分析，光子计数探测器能量分辨率越低，表示其对目标能区光子的统计区域越大，这会使得能谱重叠范围增加，不利于能量分解。但研究发现，能量分辨率在30%甚至40%的情况下，密度和原子序数的标准差变化并不十分明显。
图7是不同探测器能量分辨率的重建结果对比图，电子密度图的灰度显示窗口为[1.44 1.60]，等效原子序数图的灰度显示窗口为[18.3 20.7]。图8为同一排像素灰度曲线比较。
![]()
图 7 不同能量分辨率的重建结果对比Figure 7. Comparison of reconstruction results at different energy resolutions4.3 能量阈值及能量范围的影响
本部分考察的是高低能能谱重合程度对于成像的影响，在探测器能量分辨率相同的条件下，通过选取不同能量阈值，改变能量范围，可以得到重合度不一的高低能能谱，重合度的定义为：
$$ {\text{重合度}} = \frac{\text{能谱重叠部分面积}}{\text{全能谱范围面积}}。 $$ (12) 探测器能量分辨率选用30% @ 100 keV，具体对照实验见表3。各实验组能谱对比见图9。
表 3 能量范围实验对照表Table 3. Experimental comparison table of energy range编号 低能范围EL 高能范围EH 重合度Overlap_ratio 低能段光子数 高能段光子数 1 0~70 keV 110~160 keV 0.0025 74978 4419 2 0~80 keV 100~160 keV 0.0134 82323 7438 3 0~80 keV 80~160 keV 0.0648 82963 17132 4 0~100 keV 60~160 keV 0.1526 62561 37533 实验结果见表4，图10为重建结果的均值变化。图11为重建切片图，电子密度图的灰度显示窗口为[1.41 1.61]，等效原子序数图的灰度显示窗口为[18.6 20.5]。图12为同一排像素的灰度比较。
表 4 能量范围实验结果Table 4. Results of energy range experiment编号 重合度Overlap_ratio 等效原子序数 等效原子序数标准差 电子密度 电子密度标准差 1 0.0025 19.550129 0.316059 1.517544 0.049246 2 0.0134 19.551358 0.298083 1.517176 0.043675 3 0.0648 19.544904 0.313434 1.518503 0.044175 4 0.1526 19.554059 0.409670 1.517345 0.055390 ![]()
图 10 不同能谱重合度下的成像结果对比Figure 10. Comparison of imaging results at different energy overlap-ratios![]()
图 11 不同能量范围的重建结果对比Figure 11. Comparison of reconstruction results within different energy ranges从结果可以看出，设置能谱阈值改变能谱范围，使得高低能能谱重合度发生改变，会对成像质量造成影响。盲目减小高低能能谱重合度，并不一定使重建结果更好，因为当光子数减少时，噪声会增大。如何平衡能谱重合度和噪声是下一步的研究工作。
5. 小结
光子计数探测器作为多能谱成像的关键器件，成为研究能谱成像的一个重要内容。本文通过实验，模拟光子计数探测器在不同能量分辨率、不同能谱阈值和能量范围下的成像实验。通过比较重建的等效原子序数图和电子密度图的标准差，评价成像结果的好坏，得到上述因素对光子计数X射线CT成像精度的影响规律。
（1）针对安检CT中常用的160 kV射线源能量，光子计数探测器能量分辨率越高，成像效果越好，但是能量分辨率对于材料的分辨差异并不显著。
（2）在探测器能量分辨率确定的情况下，恰当选取能谱阈值对于成像很有帮助。通过设置合适的阈值，使高低能能谱重合度减小，往往可以得到更好的成像结果。但是随着能谱变窄、重合度下降，噪声也会加大，需要平衡这两种因素，否则更低的能谱重合度也有可能带来更差的成像结果。
本文的研究成果将有助于能谱CT系统的总体设计，如从双能成像材料分辨的角度，并不一定过高追求能量分辨率，这对于能量成像CT系统的硬件选型具有积极的指导意义。
-
表 1 拟合结果
Table 1 Fitting results
数据精简方式 曲面高度Δh/mm 曲面曲率δ/％ 模体值 拟合值 相对误差/％ 模体值 拟合值 相对误差/％ 进行异常数据剔除 14.6 14.67 0.40 0.0973 0.0978 0.40 不进行不变层去除 14.6 15.04 3.01 0.0973 0.1002 3.01 不进行稳定层去除 14.6 13.39 8.29 0.0973 0.0892 8.29 
下载: 导出CSV
-
[1] SENGUTTUVAN N, AOSHIMA M, SUMIYA K, et al. Oriented growth of large size calcium Fluoride single crystals for optical lithography[J]. Journal of Crystal Growth, 2005, 280(3/4): 462−466.
[2] TANDJAOUI A, MANGELINCK-NOEL N, REINHART G, et al. Investigation of grain boundary grooves at the solid-liquid interface during directional solidification of multi-crystalline silicon: In situ characterization by X-ray imaging[J]. Journal of Crystal Growth, 2013, 377(15): 203−211.
[3] DONG Y, SHUAI S, ZHENG T, et al. In-situ observation of solid-liquid interface transition during directional solidification of Al-Zn alloy via X-ray imaging[J]. Journal of Materials Science and Technology, 2020, 39: 113−123. doi: 10.1016/j.jmst.2019.06.026
[4] 陈福泉, 林樽达, 余炳和. 盐封下降法生长KRS-5晶体[J]. 激光与红外, 1981,(9): 50−56. CHEN F Q, LIN Z D, YU B H. Growth of KRS-5 crystal by salt sealing down method[J]. Laser & Infrared, 1981, (9): 50−56. (in Chinese).
[5] 李晓辉. 坩埚下降法CaF_2基晶体生长的数值模拟及缺陷研究[D]. 北京: 中国科学院大学, 2020. LI X H. Numerical simulation and defect analysis for the growth of alkaline-earth fluorides[D]. Beijing: University of Chinese Academy of Sciences, 2020. (in Chinese).
[6] 曹冰, 何朝明, 李柏林. 工业CT图像轮廓数据提取研究[J]. 机电工程技术, 2009,(9): 3. CAO B, HE C M, LI B L. Extraction of contour data based on industrial CT image[J]. Mechanical & Electrical Engineering Technology, 2009, (9): 3. (in Chinese).
[7] LUKÁCS G, MARTIN R R, MARSHALL A D. Faithful least-squares fitting of spheres, cylinders, cones and tori for reliable segmentation[C]//Computer Vision - ECCV'98, 5th European Conference on Computer Vision, Freiburg, Germany, June 2-6, 1998, Proceedings, Volume I. Springer Berlin Heidelberg, 1998.
[8] 邹斌, 曾理, 马睿. 工业CT图像管道拟合[J]. 计算机工程与应用, 2010,46(11): 208−210. doi: 10.3778/j.issn.1002-8331.2010.11.063 ZOU B, ZENG L, MA R. Fiting of industrial CT image pipes[J]. Computer Engineering and Applications, 2010, 46(11): 208−210. (in Chinese). doi: 10.3778/j.issn.1002-8331.2010.11.063
[9] 卢艳平, 王珏, 覃仁超. 一种剥皮算法在工业CT图像分割中的应用[J]. 无损检测, 2005,27(5): 235−237. doi: 10.3969/j.issn.1000-6656.2005.05.004 LU Y P, WANG J, QIN R C. Application of peeling method to segmentation of industrial computed tomographic image[J]. Nondestructive Testing, 2005, 27(5): 235−237. (in Chinese). doi: 10.3969/j.issn.1000-6656.2005.05.004
[10] KOHLER R. A segmentation system based on thresholding[J]. Computer Graphics and Image Processing, 1981, 15(4): 319−338. doi: 10.1016/S0146-664X(81)80015-9
[11] ADAMS R, BISHOF L. Seeded region growing[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1994, 16(6): 641-647.
[12] SHEN X, SPANN M. Segmentation of 2D and 3D images through a hierarchical clusteringbased on region modelling[C]//Image Processing, 1997. Proceedings. International Conference on. IEEE, 1997.
[13] SHELHAMER E, JONATHAN L, DARRELL T. Fully convolutional networks for semantic segmentation[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2017, 39(4): 640−651.
[14] RONNEBERGER O, FISCHER P, BROX T. U-net: Convolutional networks for biomedical image segmentation[J]. Springer International Publishing, 2015.
[15] ZHAO H, SHI J, QI X, et al. Pyramid scene parsing network[J]. 2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR). Honolulu: IEEE, 2017, 6230-6239.
[16] CHEN L C, PAPANDREOU G, KOKKINOS I, et al. Semantic image segmentation with deep convolutional nets and fully connected CRFs[C]//Computer Science. 2014: 357-361.
[17] CHEN L C, PAPANDREOU G, KOKKINOS I, et al. DeepLab: Semantic image segmentation with deep convolutional nets, atrous convolution, and fully connected CRFs[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2018, 40(4): 834−848. doi: 10.1109/TPAMI.2017.2699184
[18] CHEN L C, PAPANDREOU G, SCHROFF F, et al. Rethinking atrous convolution for semantic image segmentation[J]. arXiv preprint arXiv: 1706.05587, 2017.
[19] FITZGIBBON A, PILU M, FISHER R B. Direct least square fitting of ellipses[J]. IEEE Transactions on Pattern Analysis & Machine Intelligence, 1999.
[20] MICELI A, THIERRY R, BETTUZZI M, et al. Comparison of simulated and measured spectra of an industrial 450 kV X-ray tube[J]. Nuclear Instruments & Methods in Physics Research. Section A, Accelerators, Spectrometers, Detectors and Associated Equipment, 2007, 580(1): 123-126.
-
期刊类型引用(2)
1. 李梦雨，段诗苗，张雷，周咏春. 提高肺SBRT精度的多窗位动态组合诊疗手段. 中国CT和MRI杂志. 2024(04): 63-65 . 
百度学术
2. 李玉辉，刘龙进，徐乐意，刘远高. 胸部CT不同图像算法对人工智能辅助诊断软件肺结节检出效果的影响研究. 中国医学装备. 2023(12): 10-14 . 百度学术
其他类型引用(3)
计量
- 文章访问数: 386
- HTML全文浏览量: 161
- PDF下载量: 53
- 被引次数: 5
