1、樊知轩,阙介民,魏存峰,等.基于 CT 数据的二次曲面拟合算法研究J.CT 理论与应用研究,2023,32(1):35-42.DOI:10.15953/j.ctta.2022.040.FAN Z X,QUE J M,WEI C F,et al.Research on Quadric Surface Fitting Algorithm Based on CT DataJ.CT Theory andApplications,2023,32(1):35-42.DOI:10.15953/j.ctta.2022.040.(in Chinese).基于 CT 数据的二次曲面拟合算法研究樊知轩1a,2,阙介
2、民2,魏存峰2,3,刘宝东2,3,魏彪1b1.重庆大学 a）自动化学院；b）光电工程学院，重庆 4000442.中国科学院高能物理研究所（北京市射线成像技术与装备工程技术研究中心）,北京 1000493.中国科学院大学核科学与技术学院,北京 100049摘要：二次曲面工件在工业中比较常见，为测量物体内部的二次曲面，本文采用工业计算机断层成像（CT）技术获取物体的断层图像序列，利用 U-net 图像分割网络获得断层图像上的目标区域，对分割结果的边缘进行检测和曲线拟合并堆叠成三维点集，通过曲面拟合获取物体曲面的三维空间坐标信息。研究结果表明，本文的方法能够有效地实现边界提取和界面参数的拟合工作，拟
3、合误差在 1 以内，相比传统方法有较大改进。关键词：工业 CT；深度学习；二次曲面；曲面拟合DOI:10.15953/j.ctta.2022.040中图分类号：0 242;TP 391.41文献标识码：A随着工业技术与数字技术的不断发展，在工业 4.0 建设要求的不断推进下，传统的生产模式已经发生了根本性的变化，在生产过程中对产品进行实时数字化观察显得十分重要，工业用 X 射线计算机断层成像（X-ray computed tomography，X-CT），在晶体生长、工件装配等领域都有着重要的应用。在晶体生长领域，单晶体的制备过程主要采用 Bridgman 晶体生长方法1，由于锥部生长过程中晶
4、界的出现和冷却过程中裂纹的出现，使大尺寸、高纯单晶沿选择取向生长极为困难。对于这些问题，Tandjaoui 等2开发了一种使用X 射线成像技术观察多晶硅生成过程中的固液界面的装置，在ESRF（欧洲同步辐射装置）的 BM05 光束线上用同步辐射 X 射线照相术研究了固/液界面的形态，图像空间的分辨率为 15 m；Dong 等3利用X 光成像技术对铝锌合金定向凝固过程中固液界面形态转变过程进行了二维和三维原位观察。X-CT 技术是一种集自动控制、机械、光电、计算机图像处理等于一体的无损检测技术，它以图像方式，直观、清晰地呈现被测物质（体）的内部结构等信息，用于检测物体内部的缺陷，如空洞、裂隙等，在
5、材料科学、航空航天、核工业、机械工业、地质及考古等领域具有重要的应用价值。在测量方面，工业 CT 也可以无损地获取物体的内部形态，有助于获取内部的几何参数。二次曲面是在三维坐标系(x,y,z)下三元二次代数方程对应的所有图形的统称。在欧氏三维空间里坐标 x,y,z 之间的二次方程（系数为实数，且二次项系数不全为零）所表示的曲面，最常见的二次曲面是球面、直圆柱面和直圆锥面，此外，二次曲面还包括椭球面、双曲面和抛物面。二次曲面也是工程中常见的一种曲面，例如在晶体生长中，因为晶体的类型、生长工艺以及相关因素影响，界面呈现微凸4或者微凹5的曲面形式。在二次曲面拟合方面，曹冰等6利用基于曲率小波变换的多
6、尺度分析方法，消除大量冗余信息，提高轮廓的提取效率，为三维重构提供了条件。传统曲面拟合算法一般采用 Lukcs 等7提出的针对给定类型的二次曲面拟合方法，使用迭代法最小化目标函数进行求解，迭代速度较慢。邹斌等8提出了工业CT 图像管道拟合，利用 Facet 模型提取管道边缘的基础上，将三维圆柱面、椭圆柱面管道收稿日期：20220309。基金项目：中国科学院科研仪器设备研制项目（固液界面可视可控的光学晶体生长设备的研制（E12821V1）。第 32 卷第 1 期CT 理论与应用研究Vol.32,No.12023 年 1 月（3542）CT Theory and ApplicationsJan.,
7、2023的拟合转化为中心轴截面上的二维轮廓曲线圆、椭圆的拟合，得到了较好的结果。在用序列断层图像重构工业 CT 三维图像的过程中，一般要把目标物体从断层图像中分离出来，然后重构物体的三维模型9，因此，能否准确分割目标物体会直接影响到三维重构图像的质量。目前人们已经提出了各种类型的分割方法，比较常用的方法有阈值比较法10、区域生长法11、像素迭代聚类法12等。近年来，卷积神经网络在计算机视觉方面取得了重大进展，Shelhamer 等13提出了全卷积神经网络（fully convolutional network，FCN）之后，基于此网络的各种变形，U-net14、PSPNet15、Deeplab
8、16-18等网络解决了FCN 存在的细节不清晰、缺乏空间一致性等问题，其中 U-net在医学图像分割领域表现出了强大的实力。本文提出一种用于工业 CT 图像的二次曲面拟合算法，并通过实验证明所提出的方法的有效性。1方法1.1模体设计及实验方法为了研究用于工业 CT 图像的二次曲面拟合算法，实验中使用 X 射线 CT 系统对二次曲面进行观察。要想得到较好的观察结果，需要射线具有足够高的能量保证能够穿透样品。经过计算，采用最高电压为 450 kV 的 X 射线源工业 CT。本文通过仿真进行前期验证。设计感兴趣区域为直径 150 mm，高 14.6 mm，曲面曲率为 1/200 的球弧面，感兴趣区域
9、上方为线吸收率与感兴趣区域相差较小的背景物质，二者由厚度为 15 mm 的外围材料紧密包裹，最外层为 5 mm 厚的高吸收率物质。创建三维模体模拟这个系统，模体中每个体素的值为该位置元素的吸收率，具体数值与每个体素对应的实际大小有关，在本次仿真中，为获得较好的效果，设置每个体素对应 0.1 mm。模体如图 1 所示。1.2二次曲面提取与拟合在 CT 重建出包含物体内部信息的三维图像序列后，感兴趣区域与背景之间的分界面为需要拟合的二次曲面。为了在三维空间获得完整的二次曲面，需要对 CT 重建出的图像序列进行分割。受 CT 成像分辨率和部分容积效应的影响，难以精确定位二次曲面的最高点，而曲面中远离
10、最高点的部位重建相对准确，因此这里避开二次曲面的最顶端，重点采集下部的数据，使用基于最小二乘法的曲面拟合算法将曲面拟合成一个多项式，这样可以利用除顶部之外部分的信息来推断整个曲面的形态信息。如图 2 所示，整个拟合流程分为以下 4 个步骤。（1）图像分割：将三维图像序列中每一层的感兴趣区域分割出来，得到二值化分割图。（2）椭圆拟合：使用边界跟踪算法提取每层分割图中感兴趣区域的边缘，对提取到的边缘使用基于直接最小二乘椭圆拟合算法，得到每层拟合的椭圆边界点。（3）数据点堆叠：精简椭圆边界数据，剔背景高吸收外壳感兴趣区域外围材料图 1模体示意图Fig.1Schematic of the model椭
11、圆拟合：对感兴趣区域边缘进行提取井将其拟合成相圆数据点堆叠：将多余的数据点去除，井将剩余的数据点堆叠到一起曲面拟合：使用最小二乘法进行曲面拟合图像分割：使用 U-net 网络对图像进行分割，获取感兴趣区域图 2拟合流程图Fig.2Fitting process36CT 理论与应用研究32 卷除误差过大的数据点。保留下来的椭圆数据是一组二维序列图，需要根据 1 mm 的层间距给每层赋以高度值，最终得到二次曲面的三维点集。（4）曲面拟合：使用基于最小二乘法的曲面拟合算法获得二次曲面。在工业应用中，需要检测的材料往往尺寸较大且密度较高，X 射线难以穿透，导致 CT 重建图像噪声较高，且工件装配中结合
12、紧密，传统算法（如阈值分割法、区域生长法等）难以很好地将外壳与材料分开。当前各种基于深度学习的图像处理方法成为研究热点，在图像分割领域，卷积神经网络（convolutional neural networks，CNNs）具有许多传统算法不具备的优势，其中 U-net 是一种效果优异的轻量化卷积网络，可以实现对二次曲面快速准确分割。考虑到不同层的 CT 图像可能会受到不同程度噪声的干扰，因此在拟合出每层的椭圆数据点之后需要对数据进行精简。针对本文任务，这里设计了一种剔除重复与病态数据的精简算法，在提高拟合精度的同时大大降低了计算复杂度，精简算法在 2.1 节中详细介绍。1.3U-net 图像分割
13、网络U-net 网络在 FCN 基础上发展来，本文研究工作使用的 U-net 网络的结构如图 3 所示。U-net 由下采样和上采样两部分组成，其中 VGG16 作为下采样的骨干网络，每步下采样由一组卷积（由两个不进行填充的 33 卷积层组成）、一个激活函数和一个步长为 2 的 22 最大池化层组成，其中激活函数通常使用 ReLU（rectified linear unit，ReLU）。这样，每对特征图进行一步下采样操作，特征图的分辨率都会下降一倍，通道数提升一倍。在进行 4 次下采样后，采用 11 卷积代替 VGG16 中的全连接层，经过两次 11 卷积层后进行上采样。上采样可以看作是下采样
14、的逆运算，使用反卷积代替下采样过程中的最大池化层，实现特征图分辨率的倍增。与此同时，U-net 网络将上下采样过程中尺寸相同的特征图进行跨越连接，能够有效避免训练过程中的梯度消失。对于语义分割网络来说，网络的结构越复杂，特征提取能力越强，能够获得更大的感受野。一般来说，网络的浅层部分更容易关注输入图像的纹理信息，而深层部分会获取更加高维的特征信息。在本文任务中，纹理信息更为重要，利用跨越连接的结构可以让网络在上采样过程中依然能够关注到纹理信息。众所周知，当卷积网络进行下采样时，一定会损失某些纹理信息，丢失的信息无法在上采样过程中获得，但是 U-net 网络的跨越连接能够使特征图在上采样过程中不
15、断地从前序网络中获得丢失的纹理信息。本文所涉及到的图像语义内容较为简单，不同图像中的物体组成结构比较固定，不需要对信息进行过多的过滤，所有细节信息的重要程度大致相同，使用跨越连接结构可以将高级和低级的特征都利用上。同时网络结构较为简单，在面对较小的数据集时基本不会出现过拟合。1 024 1 024Bottleneck Conv512 512512 512 512 5121/41/81/8256 256256 256 256 256128 128128 128128128646464 64 64Softmax图 3U-net 网络结构图Fig.3U-net network structure d
16、iagram1 期樊知轩等：基于 CT 数据的二次曲面拟合算法研究371.4最小二乘法拟合最小二乘法（least square method，LSM）通过最小化误差的平方和寻找数据的最优函数匹配。本文使用最小二乘法进行椭圆拟合与曲面拟合19。椭圆方程可以用以下二阶多项式来表示：a1x2+a2xy+a3y2+a4x+a5y+a6=0。（1）a=(a1,a2,a3,a4,a5,a6)x=(x2,xy,y2,x,y,1)T a x=0令，方程可以表示为，所以拟合椭圆的最优化问题可以表示为：min?a x?2，（2）x a a=0其中表示数据样本集合，矩阵维度为 6n，每列表示 1 个样本。表示椭圆方程的参数。为了避免出现这种无意义的解，需要对参数添加约束条件。z=a0+a1x+a2y+a3x2+a4xy+a5y2在拟合曲面时，选择二次方程作为待拟合曲面的标准方程。对样本数据集中观测点计算误差平方并求和为：S=mi=1(f(x,y)z)2。（3）aj=(j=0,1,2,5)Xv为了寻找最优化的多项式参数估计值，对于给定的观测点求解误差平方和取最小值的参数。通过构建对应的范德蒙矩阵，利用矩阵运算
