Differentiated Backprojection Reconstruction Method for Square Cross-section FOV Rotational CL
-
摘要:
随着现代工业的发展,计算机断层成像(CT)在无损检测中发挥了重大的作用。由于成像视野或X射线穿透能力的限制,CT难以实现对板状物体的扫描和成像,针对板状物体的检测,通过更改扫描几何和重建算法,发展出了计算机层析成像技术(CL)。CL重建借鉴了CT重建方法,有解析重建方法和迭代重建方法,滤波反投影类算法快速但是重建与真值差距很大;迭代方法通常比解析重建准确,但是耗时很长。本文借鉴圆轨道扇束CT的微分反投影(DBP)重建方法,推导了方形截面视野旋转CL的DBP重建方法,以仿真数据和实际印刷电路板扫描数据进行重建,DBP方法的重建图像质量与滤波反投影算法相似,且在处理投影截断问题上有一定优势。
Abstract:With advancements in modern industry, computed tomography (CT) has played a significant role in nondestructive testing. However, owing to the limitations in the imaging field of view and X-ray penetration capability, CT faces challenges in the scanning and imaging of plate-like objects. To inspect plate-like objects, computer lithography (CL) has been developed by modifying the scanning geometry and reconstruction algorithms. CL reconstruction methods are inspired by CT reconstruction techniques and include analytical and iterative reconstruction methods. Filtered backprojection algorithms are fast, but often result in significant discrepancies from the true values, and iterative methods, although more accurate than analytical reconstruction, are usually time-consuming. This study draws on the differentiated backprojection (DBP) reconstruction method of circular trajectory fan-beam CT to derive a DBP reconstruction method for a square cross-sectional field-of-view (FOV) rotational CL. Using both simulated data and actual printed circuit board scanning data for reconstruction, the image quality of the DBP method was found to be similar to that of the filtered back projection algorithm, with certain advantages in addressing projection truncation challenges.
-
自发现X射线后,医学上就开始用它来检查人体疾病。随着现代工业的发展,计算机断层成像(CT)在无损检测和逆向工程中发挥了重大作用。工业CT无损检测技术能够很好地检测出气孔、夹杂、裂纹等各种常见缺陷,与其他常规无损检测技术相比,工业CT无损检测技术的空间和密度分辨率高,不受工件材料种类和几何形状的限制,可生成材料缺陷的三维图像,极具研究价值,在化工、石油、机械、宇航、电子、核能等工业中得到了非常广泛的应用。
对于如印刷电路板(printed circuit board,PCB)等板壳状结构的物体,在平行和垂直于板面的两个方向上,X射线透射强度的差异极大,从而很难兼顾两个方向上检测信号的有效性,无法达到满意的检测能力。领域内针对特殊薄板构件的检测,发展出了计算机层析成像(computer laminography,CL)技术,主要用于大尺寸、低厚度板状物体的无损检测[1]。针对CL的图像重建方法从很多方面借鉴了CT的重建算法,有解析重建方法也有迭代重建方法。2010年,Yang等[2]实现了锥束圆轨迹扫描的一种滤波反投影重建算法,并将之应用在电路板检测中。随后在2013年,Myagotin等[3]就给出了平行束CL的FBP重建方法。Voropaev等[4]2016年提出的直接傅里叶逆变换法是在傅里叶域进行反投影,再使用三维逆傅里叶变换将反投影结果转换为空间域。该算法适用于平行束CL,在重建质量上与传统的滤波反投影方法相当,但重建时间大幅减少。2021年Sun等[5]提出了一种基于投影变换的锥束层析成像重建方法,通过将锥束层析成像数据转换为满足FDK算法要求的投影数据,再利用FDK算法进行重建,有效降低了图像伪影,提高了图像质量。这些解析算法的效率很高,但是重建图像都因扫描几何的限制和数据不完备而存在有限角伪影和混叠伪影。迭代重建算法可以引入先验知识到重建过程,能够适应各种成像模式,重建质量相对解析重建有所提高。迭代类的算法自20世纪80年代后期就开始被应用到CL成像中。常用的迭代重建算法有代数重建方法(algebraic reconstruction technique,ART)[6],联合迭代重建方法(simultaneous iterative reconstruction technique,SIRT)[7],联合代数重建方法(simultaneous ART,SART)[8],最大后验估计(maximum a posteriori,MAP)等算法[9]。
对于CL成像,层间混叠伪影极大地影响了重建图像质量。 Dobbins等人提出的矩阵反演融合法(matrix inversion tomosynthesis,MITS)[10],通过引入一组脉冲响应函数来模拟反投影时层间模糊的过程,进而结合线性系统理论和成像几何知识,将平面内的真实细节从层间模糊中分离出来,从而提高了层内结构的对比度。Lauritsch等[11]从滤波反投影的思想出发,通过对扫描几何的全面考虑,设计了三维滤波函数,显著提升了重建图像质量。为了抑制CL重建图像的中高衰减材料产生的混叠伪影,Levakhina等[12]提出了加权SART算法,通过定位可能在体素中产生严重伪影的投影数据,减少这些投影数据对更新项的贡献,来抑制层间混叠,提高细节对比度。Zhao等提出了一种基于边缘信息扩散的重建方法[13]。该方法利用了切片内的边缘信息,通过扩散过程修正切片间伪影,并与SART算法结合实现迭代重建,从而有效地抑制了CL重建中的切片间伪影,恢复了垂直边缘信息,提高了重建图像的质量。
由于CL成像特殊的扫描几何,投影数据的不完备会使得重建图像不可避免地存在不足,并且在CL成像的过程中由于实际情况的限制,可能得到的投影数据还会出现截断问题。例如在实现对PCB高分辨率成像时,射线源需要非常靠近被成像物体,而探测器的面积有限,这意味着CL不能扫描整块电路板,只能扫描感兴趣区域,这种扫描获得的投影数据就是有截断的。迭代重建算法通常需要耗费较长时间,不能满足工业检测的时效性要求,并且迭代算法很难只重建感兴趣区域,无法处理截断问题,所以在工业检测CL成像中迭代重建算法很难得到应用。DBP重建算法是一种解析算法,它的优点是重建速度快,并且在某些特定情况下[14],可以方便地处理投影数据的截断问题,并且重建图像质量与标准滤波反投影方法非常相似。潘晓川等[15]在2004年提出了螺旋锥束CT的DBP重建算法,提供了一种新的CT精确重建方法。同年,Noo等[14]推导了平行束几何和扇束圆轨道线探测器的DBP重建算法。2005年,Noo等[16]又推导了三维锥束的DBP重建算法。2013年,Kudo等[17]阐述了DBP在内部CT图像重建中的应用,通过将感兴趣区分解为一系列希尔伯特线,并利用希尔伯特变换的逆变换,实现了从截断投影数据中精确重建图像。同年,Tang等[18]提出的径向微分内部成像方法使用微分反投影和凸集投影算法进行图像重建。该方法能够在不降低图像质量的情况下减少CT扫描的辐射剂量,并有效地突出感兴趣区域的边缘信息。Yoseob等[19]于2019年使用截断的微分反投影数据作为神经网络的输入,用网络学习希尔伯特变换的卷积核和偏移量,获得了新的重建方法。该方法泛化能力强,并且可以有效地去除杯状伪影,并保留图像中的细节结构和纹理。2021年,Ernst等[20]提出了一种结合卷积神经网络和微分反投影的方法,这种方法先对稀疏投影的微分反投影结果进行Hilbert逆变换,然后使用卷积神经网络对其进行优化。这种方法可以有效地减少条带伪影和锥束伪影。综上可知,DBP重建算法在处理截断投影方面具有优势,在内部CT中应用广泛。
接下来的内容中,我们将首先介绍方形截面视野旋转CL成像的系统几何,然后推导相应的DBP重建方法。在第二节我们用计算机模拟的投影数据和扫描实际电路板得到的投影数据对我们提出的DBP重建方法进行验证,并与滤波反投影算法[21]进行对比。
1. CL成像几何及DBP重建方法
1.1 CL成像几何
方形截面视野旋转CL成像系统如图1所示,其主要特点是探测器水平排列并在旋转中保持方向不变,因此其成像视野的横截面呈现为正方形。我们将射线源
$ S $ 和探测器中心$ {O}_{\mathrm{d}} $ 的连线$ S{O}_{\mathrm{d}} $ 定义为中心射线,中心射线与旋转轴(黑色虚线所示)相交于点$ O $ 。$ u $ ,$ v $ 为探测器坐标轴,$ u $ 轴、$ v $ 轴的原点为探测器中心$ {O}_{\mathrm{d}} $ 。以$ O $ 为原点建立坐标系$ Oxyz $ ,$ x $ 轴,$ y $ 轴与探测器$ u $ 轴,$ v $ 轴平行,红色虚线轨迹为射线源运动轨道,与$ xOy $ 平面平行,蓝色虚线轨迹为探测器中心运动轨道,与$ xOy $ 平面平行,当射线源$ S $ 沿着射线源轨道逆时针方向旋转扫描时,探测器中心$ {O}_{\mathrm{d}} $ 相应地沿着探测器轨道进行逆时针方向旋转,为精确定位射线源和探测器的位置,引入旋转角$ \xi $ ,$ \xi $ 为橙色中心射线在射线源轨道平面上的投影与$ x $ 轴之间的夹角,倾斜角$ \alpha $ 为中心射线与$ xOy $ 平面的夹角,量化了射线的倾斜程度。对感兴趣区域进行成像是工业CL成像的常见任务,由于感兴趣区域处于物体之内,图1所示的方形截面视野旋转CL成像系统扫描所得的投影数据会不可避免地包含感兴趣区域之外的信息,这种扫描获得的投影数据就是截断的。领域内能有效解决截断问题的一种重建算法是DBP重建算法。
根据Noo F等人的工作[14,16],DBP重建算法的原理是通过微分反投影操作从投影获得待重建物体的Hilbert变换图像,再对Hilbert图像做逆Hilbert变换即可得到重建图像。所以DBP算法可以分为两步,在第一步中,先对投影数据沿垂直中心射线的方向求导,再进行反投影,所得结果即为物体的Hilbert变换。第二步是在微分反投影图像中在特定方向的线上应用Hilbert逆变换,就可以得到重建图像,这些线被称为PI线。下面我们将从扇束圆轨道线探测器的DBP重建公式出发,推导方形截面视野旋转CL成像系统的DBP重建方法。
1.2 圆轨道扇束DBP重建方法简介
对于如图2所示的圆轨道扇束等距线探测器CT成像,设待重建物体函数为
$ f=f\left(\mathbf{x}\right)=f\left(x,y\right) $ ,$ \mathbf{x}=(x,y) $ 为重建区域某一点的空间坐标,其投影数据表示为:$$ p\left(\xi ,l\right)=\int _{0}^{\mathrm{\infty }}f\left({\bf{\Lambda}}\left(\xi \right)+t\boldsymbol{\kappa }\left(\xi ,l\right)\right)dt , $$ (1) 其中
$ \xi $ 表示旋转角度,$ l $ 表示探测器坐标轴,$ {R}_{0} $ 为射线源到旋转中心的距离,$ {D}_{0} $ 为射线源到探测器中心的距离,则$ \mathit{\Lambda }\left(\xi \right)={R}_{0}(-\mathrm{s}\mathrm{i}\mathrm{n}\left(\xi \right),\mathrm{c}\mathrm{o}\mathrm{s}\left(\xi \right)) $ 为射线源坐标,$ \mathit{\kappa }\left(\xi ,l\right) = \displaystyle\frac{({D}_{0}\mathrm{s}\mathrm{i}\mathrm{n}\left(\xi \right) + l\mathrm{c}\mathrm{o}\mathrm{s}\left(\xi \right),-{D}_{0}\mathrm{c}\mathrm{o}\mathrm{s}\left(\xi \right) + l\mathrm{s}\mathrm{i}\mathrm{n}\left(\xi \right))}{\sqrt{{{D}_{0}}^{2} + {l}^{2}}} $ 为射线的单位方向向量。取β = (cosθ, sinθ)为PI线方向,θ为PI线与y轴的夹角,t = xβT代表重建区域中的点所在的PI线位置索引,则DBP重建公式为:
$$ \begin{split}&{f}\left(\mathbf{x}\right)=f\left(t{\boldsymbol{\beta} }+s{\boldsymbol{\beta }}^{\perp }\right)=\frac{-1}{\sqrt{\left(s-{L}_{t}\right)\left({U}_{t}-s\right)}}\times\\& \left(\int _{{L}_{t}}^{{U}_{t}}\sqrt{\left({s}{'}-{L}_{t}\right)\left({U}_{t}-{s}{'}\right)}\frac{{b}_{\theta }\left(t\boldsymbol{\beta }+s{\boldsymbol{\beta }}^{\perp }\right)}{-2{\pi }^{2}\left(s-{s}{'}\right)}d{s}{'}+{C}_{t}\right) \text{,} \end{split}$$ (2) 其中,
$ {\boldsymbol{\beta}} ^{\perp }=(-\mathrm{s}\mathrm{i}\mathrm{n}\theta ,\mathrm{cos}\theta ),s=\mathbf{x}{{\boldsymbol{\beta }}^{\perp }}^{\mathrm{T}} $ ,$ {U}_{t} $ 和$ {L}_{t} $ 为以$ t $ 为索引的PI线的边界,表示重建区域的大小,$ {b}_{\theta }\left(\mathbf{x}\right)={b}_{\theta }\left(t\boldsymbol{\beta }+s{\boldsymbol{\beta }}^{\perp }\right) $ 为投影数据微分后在PI线上的反投影结果,在区间$ [{L}_{t},{U}_{t}] $ 内,$$\begin{split} {b}_{\theta }\left(\mathbf{x}\right)=&\frac{1}{2}\int _{0}^{2\pi }\frac{{R}_{0}{D}_{0}}{{{T}_{0}}^{2}sgn}\left(\mathrm{sin}\left(\xi +\mathrm{arctan}\frac{{l}^{\mathrm{*}}}{{D}_{0}}-\theta \right)\right)\\& \frac{d}{dl}\left\{\frac{{D}_{0}}{\sqrt{{{D}_{0}}^{2}+{l}^{2}}}p\left(\xi ,l\right)\right\}{|}_{l={l}^{\mathrm{*}}}d\xi +\\& \frac{p\left({\xi }_{1},{l}_{1}^{\mathrm{*}}\right)}{\left \| \mathbf{x}-{\mathit{\Lambda }}_{1}\right \| }-\frac{p\left({\xi }_{2},{l}_{2}^{\mathrm{*}}\right)}{\left \|\mathbf{x}-{\mathit{\Lambda }}_{2}\right \| } \text{,} \end{split}$$ (3) 其中
$ {T}_{0}={R}_{0}+x\mathrm{s}\mathrm{i}\mathrm{n}\left(\xi \right)-y\mathrm{c}\mathrm{o}\mathrm{s}\left(\xi \right) $ 为点$ \mathbf{x} $ 在中心射线上的投影与射线源之间的距离,$ {l}^{*}=\displaystyle\frac{{D}_{0}}{{T}_{0}} (x\mathrm{c}\mathrm{o}\mathrm{s}\left(\xi \right)+ y\mathrm{s}\mathrm{i}\mathrm{n}\left(\xi \right)) $ 为点$ \mathbf{x} $ 对应的探测器坐标。公式(3)中最后两项的$ {\boldsymbol{\Lambda }}_{1} $ 和$ {\boldsymbol{\Lambda }}_{2} $ 对应PI线与射线源轨道的两个交点,$ {\xi }_{1}=\theta -\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{s}\mathrm{i}\mathrm{n}\left(\right(\mathbf{x}{\boldsymbol{\beta }}^{\mathrm{T}})/{R}_{0}) $ ,$ {\xi }_{2}= \pi +2\theta -{\xi }_{1} $ ,$ {\boldsymbol{\Lambda }}_{\mathit{i}} = {R}_{0}\left(-\mathrm{sin}({\xi }_{i}),\mathrm{cos}({\xi }_{i})\right),i=1,2 $ 。公式(2)中
$ {C}_{t} $ 为一常数,通过这一常数来保证重建结果沿PI线方向的投影与实际投影一致,设$ {\epsilon }_{t} $ 为一小量,在区间$ [{L}_{t}+{\varepsilon }_{t},{U}_{t}-{\varepsilon }_{t}] $ 外,已知$ f\left(t\boldsymbol{\beta }+s{\boldsymbol{\beta }}^{\perp }\right) $ 等于0,则$ {C}_{t} $ 的一种计算方式如下,$$ {C}_{t}=\displaystyle\frac{-p\left({\xi }_{0},{l}_{0}^{\mathrm{*}}\right)-\displaystyle\int _{{L}_{t}-{\varepsilon }_{t}}^{{U}_{t}-{\varepsilon }_{t}}\displaystyle\frac{1}{\sqrt{\left(s-{L}_{t}\right)\left({U}_{t}-s\right)}}\int _{{L}_{t}}^{{U}_{t}}\sqrt{\left({s}{'}-{L}_{t}\right)\left({U}_{t}-{s}{'}\right)}\displaystyle\frac{{b}_{\theta }\left(t\boldsymbol{\beta }+{s}{'}{\boldsymbol{\beta }}^{\perp }\right)/(-2\pi )}{\pi \left(s-{s}{'}\right)}d{s}{'}ds}{\displaystyle\int _{{L}_{t}-{\varepsilon }_{t}}^{{U}_{t}-{\varepsilon }_{t}}\displaystyle\frac{1}{\sqrt{\left(s-{L}_{t}\right)\left({U}_{t}-s\right)}}ds} \text{,} $$ (4) 其中
$ p\left({\xi }_{0},{l}_{0}^{\mathrm{*}}\right) $ 是以$ t $ 为索引的PI线方向的投影,$ {\xi }_{0}=\theta -\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{s}\mathrm{i}\mathrm{n}(t/{R}_{0}) $ ,$ {l}_{0}^{\mathrm{*}}=\displaystyle\frac{{D}_{0}}{{T}_{0}}(x\mathrm{c}\mathrm{o}\mathrm{s}\left({\xi }_{0}\right)+y\mathrm{s}\mathrm{i}\mathrm{n}\left({\xi }_{0}\right)) $ 。1.3 CL几何DBP重建
对于图1所示的CL成像系统,设待重建物体函数为
$ f=f\left(\mathbf{x}\right)=f\left(x,y,z\right) $ ,$ \mathbf{x}=(x,y,z) $ 为重建区域某一点的空间坐标,其投影可用如下公式表示:$$ {p}^{\mathrm{C}\mathrm{L}}\left(\xi ,u,v\right)=\int _{0}^{\mathrm{\infty }}f\left(\boldsymbol{A}\left(\xi \right)+t\boldsymbol{\lambda }\left(\xi ,u,v\right)\right)dt \text{,} $$ (5) 其中
$ \boldsymbol{A}\left(\xi \right)=(-|SO|\mathrm{c}\mathrm{o}\mathrm{s}\left(\alpha \right)\mathrm{c}\mathrm{o}\mathrm{s}\left(\xi \right),-|SO|\mathrm{c}\mathrm{o}\mathrm{s}\left(\alpha \right) \mathrm{s}\mathrm{i}\mathrm{n}\left(\xi \right), -|SO\left|\mathrm{s}\mathrm{i}\mathrm{n}\left(\alpha \right)\right) $ 为射线源空间坐标,$ \mathit{\lambda }\left(\xi ,u,v\right)= \displaystyle\frac{\left(\right|S{O}_{\mathrm{d}}|\mathrm{c}\mathrm{o}\mathrm{s}\left(\alpha \right)\mathrm{c}\mathrm{o}\mathrm{s}\left(\xi \right)+u,|S{O}_{\mathrm{d}}|\mathrm{c}\mathrm{o}\mathrm{s}\left(\alpha \right)\mathrm{s}\mathrm{i}\mathrm{n}\left(\xi \right)+v,|S{O}_{\mathrm{d}}\left|\mathrm{s}\mathrm{i}\mathrm{n}\left(\alpha \right)\right)}{\sqrt{{\left(\left|S{O}_{\mathrm{d}}\right|\mathrm{c}\mathrm{o}\mathrm{s}\left(\alpha \right)\mathrm{c}\mathrm{o}\mathrm{s}\left(\xi \right)+u\right)}^{2}+{\left(\left|S{O}_{\mathrm{d}}\right|\mathrm{c}\mathrm{o}\mathrm{s}\left(\alpha \right)\mathrm{s}\mathrm{i}\mathrm{n}\left(\xi \right)+v\right)}^{2}+{\left(\left|S{O}_{\mathrm{d}}\right|\mathrm{s}\mathrm{i}\mathrm{n}\left(\alpha \right)\right)}^{2}}} $ 为射线的单位方向向量。在图1所示的CL成像系统中,射线源轨道平面与被扫描物体没有交集,且探测器接收的射线与射线源轨道平面之间有很大的锥角,根据Tuy-Smith条件,方形截面视野旋转CL成像系统的投影数据是不完备的,所以该CL成像系统没有对应的精确重建方法。根据Noo等[16]的研究结果,要通过DBP重建方法实现精确重建,需要满足的条件是:当射线源在PI线的两个端点之间移动时,如果完整收集了覆盖PI线与重建区域的交集的投影数据,则可以在 PI线上进行准确的重建。由于图1中CL系统的射线源仅在一个
$ z=-\left|SO\right|\mathrm{cos}(\alpha ) $ 平面沿半径为$ \left|SO\right|\mathrm{s}\mathrm{in}(\alpha ) $ 的圆轨道移动,而重建物体中心位于原点$ O $ ,且与射线源轨道平面没有交集,所以不存在穿过物体的PI线。我们只能采用一种近似的方法进行微分反投影重建。如图3所示,我们在重建区域内取PI线与$ xOy $ 平面平行且与$ x $ 轴的夹角记为$ \theta $ (橙色平面为重建区域内一轴向平面,$ s $ 轴与PI线平行,$ t $ 轴与$ s $ 轴垂直,坐标轴原点位于旋转轴上),并忽略PI线与射线源轨道平面的轴向距离来计算微分反投影。获得微分反投影图像后应用公式(2)的Hilbert逆变换来获得最终重建图像。根据公式(3),微分反投影的计算需要沿垂直中心射线的方向对投影数据进行求导。如图3所示,我们在探测器平面内定义穿过探测器中心
$ {O}_{\mathrm{d}} $ 且与中心射线垂直的方向为$ {v}' $ 轴, 定义$ {u}' $ 轴与$ {v}' $ 轴垂直,可得$ {u}'=u\mathrm{c}\mathrm{o}\mathrm{s}\xi +v\mathrm{s}\mathrm{i}\mathrm{n}\xi $ ,$ {v}'=-u\mathrm{s}\mathrm{i}\mathrm{n}\xi +v\mathrm{c}\mathrm{o}\mathrm{s}\xi $ 。单独考虑CL探测器上的一个点,这个点有一个对应的虚拟等距线探测器$ {v}' $ ,这个虚拟探测器与射线源组成了一个如图2所示的扇束几何,这样,CL探测器的每个点都有各自对应的扇束几何,参考扇束投影的求导计算方式,沿$ {v}^{\text{'}} $ 方向对投影数据进行求导:$$\begin{split} p^{\mathrm{Diff}}(\xi,u^{'*},v^{'*})\triangleq & \frac d{dv'}\left\{\frac D{\sqrt{D^2+v^{'2}+H^2}}p^{\mathrm{CL}}(\xi,u^{'*},v')\right\}\Bigg|_{v'=v^{'*}}= \left(-\sin\xi\frac d{du}\left\{\frac D{\sqrt{D^2+(-u\mathrm{sin}\xi+v\mathrm{cos}\xi)^2+H^2}}\right.\right.\\& \left.\left.p^{\mathrm{CL}}(\xi,u,\nu)\right\}+\cos\xi\frac d{dv}\left\{\frac D{\sqrt{D^2+(-u\mathrm{sin}\xi+v\mathrm{cos}\xi)^2+H^2}}p^{\mathrm{CL}}(\xi,u,\nu)\right\}\right)\Bigg|_{u=u^*, v=v^*},\end{split}$$ (6) 其中,
$ ({u}^{'*},{v}^{'*}) $ 为点$ \mathbf{x} $ 对应的虚拟探测器坐标,$ ({u}^{\mathrm{*}},{v}^{\mathrm{*}}) $ 为实际探测器坐标,$ D=\left|S{O}_{\mathrm{d}}\right| \mathrm{cos}(\alpha)+ u\mathrm{cos}\left(\xi \right)+ v\mathrm{sin}\left(\xi \right) $ 为探测器与虚拟线探测器$ {v}' $ 之间的距离,$ H=\left|S{O}_{\mathrm{d}}\right|\mathrm{sin}(\alpha ) $ 为探测器轨道平面与射线源轨道平面之间的距离。参考锥束DBP重建方法[16]的推导过程,可得CL几何的微分反投影计算公式为,
$$\begin{split} b_\theta(\mathbf{x})= & \frac{1}{2}\int_0^{2\pi}\frac{RD}{T^2}\mathrm{sgn}\left(\sin\left(\xi+\arctan\frac{v^{'*}}{D}-\theta\right)\right) \frac{d}{dv'}\left\{\frac{D}{\sqrt{D^2+v^{'^2}+H^2}}p^{\mathrm{CL}}\big(\xi,u^{'*},v'\big)\right\}\Bigg|_{v'=v^{'*}}\mathrm{d}\xi+\\& \frac{p^{\mathrm{CL}}(\xi_1,u_1^*,v_1^*)}{\|\mathbf{x}-\boldsymbol{A}_1\|}- \frac{p^{\mathrm{CL}}(\xi_2,u_2^*,v_2^*)}{\|\mathbf{x}-A_2\|}=\frac12\int_0^{2\pi}\frac{RD}{T^2}\mathrm{sgn}\left(\sin\left(\xi+\arctan\frac{-u^*\sin\xi+v^*\cos\xi}D- \theta\right)\right)\biggl(-\sin\xi\frac d{du}\\& \left\{\frac D{\sqrt{D^2+(-u\mathrm{sin}\xi+v\mathrm{cos}\xi)^2+H^2}}p^{\mathrm{CL}}(\xi,u,v)\right\}+ \cos\xi\frac d{dv}\Bigg\{\frac D{\sqrt{D^2+(-u\mathrm{sin}\xi+v\mathrm{cos}\xi)^2+H^2}}\\& p^{\mathrm{CL}} (\xi,u,v)\Bigg\}\biggl)\Big|_{v=v^*, u=u^*}\mathrm{d}\xi+ \frac{p^{\mathrm{CL}}(\xi_1,u_1^*,v_1^*)}{\|\mathbf{x}-A_1\|}-\frac{p^{\mathrm{CL}}(\xi_2,u_2^*,v_2^*)}{\|\mathbf{x}-A_2\|}, \end{split}$$ (7) 其中
$ R=\left|SO\right|\mathrm{cos}(\alpha) $ ,$ T=R+x\mathrm{cos}(\xi )+y\mathrm{s}\mathrm{in}(\xi ) $ ,$ {\xi }_{1}= \theta - \mathrm{a}\mathrm{r}\mathrm{c}\mathrm{s}\mathrm{i}\mathrm{n}\left(\right(-x\mathrm{sin}\theta +y\mathrm{cos}\theta )/R) $ ,$ {\xi }_{2}=\pi +2\theta -{\xi }_{1} $ ,$ {\mathit{A}}_{\mathit{i}}= (-R \mathrm{cos}({\xi }_{i}),-R\mathrm{sin}({\xi }_{i}),-|SO\left|\mathrm{sin}(\alpha )\right) $ ,$ {u}^{\mathrm{*}},{v}^{\mathrm{*}} $ 的计算方式如下,$$ {u}^{\mathrm{*}}=\frac{(x+R\mathrm{cos}(\xi ))H}{z+\left|SO\right|\mathrm{s}\mathrm{in}(\alpha )}-\left|S{O}_{\mathrm{d}}\right|\mathrm{cos}(\alpha )\mathrm{cos}\left(\xi \right) , $$ (8) $$ {v}^{\mathrm{*}}=\frac{\left(y+R\mathrm{sin}(\xi )\right)H}{z+\left|SO\right|\mathrm{si}\mathrm{n}(\alpha )}-\left|S{O}_{\mathrm{d}}\right|\mathrm{cos}(\alpha )\mathrm{sin}\left(\xi \right) , $$ (9) 计算出图3中CL系统的微分反投影之后,对
$ z $ 方向不同层进行Hilbert逆变换获得重建结果:$$ {f}^{z}\left(x,y\right)={f}^{z}\left(t{\boldsymbol{\beta }}^{\perp }+s\boldsymbol{\beta }\right)=\frac{-1}{\sqrt{\left(s-{L}_{t}\right)\left({U}_{t}-s\right)}}\times \left(\int _{{L}_{t}}^{{U}_{t}}\sqrt{\left({s}'-{L}_{t}\right)\left({U}_{t}-{s}'\right)}\frac{{{b}_{\theta }}^{z}\left(t{\boldsymbol{\beta }}^{\perp }+s\boldsymbol{\beta }\right)}{-2{\pi }^{2}(s-{s}')}d{s}'+{C}_{t}^{z}\right) \text{,} $$ (10) $$ {C}_{t}^{z}=\displaystyle\frac{-{p}^{z}\left(\theta ,t\right)-\displaystyle\int _{{L}_{t}-{\varepsilon }_{t}}^{{U}_{t}-{\varepsilon }_{t}}\displaystyle\frac{1}{\sqrt{\left(s-{L}_{t}\right)\left({U}_{t}-s\right)}}\int _{{L}_{t}}^{{U}_{t}}\sqrt{\left({s}'-{L}_{t}\right)\left({U}_{t}-{s}'\right)}\displaystyle\frac{{{b}_{\theta }}^{z}\left(t{\boldsymbol{\beta }}^{\perp }+s\boldsymbol{\beta }\right)}{-2{\pi }^{2}\left(s-{s}'\right)}d{s}'ds}{\displaystyle\int _{{L}_{t}-{\varepsilon }_{t}}^{{U}_{t}-{\varepsilon }_{t}}\displaystyle\frac{1}{\sqrt{\left(s-{L}_{t}\right)\left({U}_{t}-s\right)}}ds} 。 $$ (11) 在计算
$ {C}_{t}^{z} $ 时,必须考虑$ {p}^{z}\left(\theta ,t\right) $ 的值,即点$ (x,y,z) $ 沿PI线方向的投影,该点对应$ (t,s,z) $ ,对于固定的$ t $ 和$ z $ ,$ {p}^{z}\left(\theta ,t\right) $ 保持不变。然而,在CL成像系统中,射线源轨道平面与扫描物体不存在交集,因此无法直接从投影数据获得$ {p}^{z}\left(\theta ,t\right) $ 的值。为了估算$ {p}^{z}\left(\theta ,t\right) $ ,我们采用了一种方法,即计算通过PI中点且平行于$ zOs $ 平面的射线的投影来估计$ {p}^{z}\left(\theta ,t\right) $ 的值。通常,存在两条满足这个条件的射线,如图3所示,分别为射线1和射线2。射线1的投影值表示为$ {p}^{\mathrm{C}\mathrm{L}}\left({\xi }_{3},{u}_{3}^{\mathrm{*}},{v}_{3}^{\mathrm{*}}\right) $ ,射线2的投影值表示为$ {p}^{\mathrm{C}\mathrm{L}}\left({\xi }_{4},{u}_{4}^{\mathrm{*}},{v}_{4}^{\mathrm{*}}\right) $ ,我们取这两条射线的投影值的算术平均数作为$ {p}^{z}\left(\theta ,t\right) $ 的近似:$$ {p}^{z}\left(\theta ,t\right)\approx \frac{1}{2}({p}^{\mathrm{C}\mathrm{L}}\left({\xi }_{3},{u}_{3}^{\mathrm{*}},{v}_{3}^{\mathrm{*}}\right)+{p}^{\mathrm{C}\mathrm{L}}\left({\xi }_{4},{u}_{4}^{\mathrm{*}},{v}_{4}^{\mathrm{*}}\right)) \text{,} $$ (12) 其中,
$ {\xi }_{3}=\theta -\mathrm{arcsin}\left(\displaystyle\frac{t}{\left|SO\right|\mathrm{cos}(\alpha )}\right) $ ,$ {\xi }_{4}=\mathrm{\pi }+2\theta -{\xi }_{3} $ ,$ {t}=\mathbf{x}{{\boldsymbol{\beta }}^{\perp }}^{\mathrm{T}} $ ,$ {u}_{i}^{*},{v}_{i}^{*},i=3, 4 $ 以公式(8)和公式(9)计算,至此我们得到了CL成像系统的DBP解析重建方法。根据Yoseob等[22]的研究,在频域观察DBP重建结果,沿矢状方向的DBP重建会导致冠状方向频域数据的缺失,而沿冠状方向的DBP重建会导致矢状方向频域数据的缺失,他们通过频谱混合的方法解决这一问题,冠状方向中缺失的频谱成分可以通过矢状方向的结果来补偿,反之亦然。具体操作是首先将冠状和矢状三维重建结果的轴向图像进行二维傅里叶变换转换为频域,然后应用蝶型频谱加权来抑制缺失的频域信号,并将它们组合在一起。这个过程用数学方法表示为:
$$ {f}^{\mathrm{c}\mathrm{o}\mathrm{m}}={\mathcal{F}}^{-1}\left\{\omega \odot \mathcal{F}\left\{{f}^{\mathrm{c}\mathrm{o}\mathrm{r}}\right\}+(1-\omega )\odot \mathcal{F}\left\{{f}^{\mathrm{s}\mathrm{a}\mathrm{g}}\right\}\right\} , $$ (13) $ \mathcal{F} $ 和$ {\mathcal{F}}^{-1} $ 分别表示正向和反向二维傅里叶变换,$ {f}^{\mathrm{c}\mathrm{o}\mathrm{r}} $ 和$ {f}^{\mathrm{s}\mathrm{a}\mathrm{g}} $ 分别表示沿冠状和矢状方向的DBP重建图像,$ \odot $ 表示逐像素相乘,本文用来抑制缺失的频域信号的蝶型频谱滤波器如图4所示,定义为:$$ \omega (\rho ,\phi )=\left\{\begin{array}{cc}1& \phi \in \left[0,\displaystyle\frac{\pi }{4}\right),\\ 0.5& \phi =\pi /4,\\ 0& \phi \in \left(\displaystyle\frac{\pi }{4},\displaystyle\frac{\pi }{2}\right],\\ \omega (\rho ,\pi -\phi )& \phi \in \left(\displaystyle\frac{\pi }{2},\pi \right],\\ \omega \left(\rho ,2\pi -\phi \right)& \phi \in \left(\pi ,2\pi \right),\end{array}\right. $$ (14) $ (\rho ,\phi ) $ 为频域的极坐标。2. CL重建实验结果
我们用计算机仿真的投影数据和实际CL系统扫描所得的电路板投影数据对我们提出的DBP重建方法进行了验证,并将重建结果与FDK重建方法进行对比。
2.1 仿真数据重建
选取三维离散化SheppLogan模体作为仿真模体,将SheppLogan模体中心置于坐标原点
$ O $ ,按照如表1的CL几何参数进行投影仿真。以仿真投影数据进行FDK和DBP重建,重建体素尺寸为0.0133mm3,重建中心为坐标原点,重建体素数量为30×300×300($ z $ ×$ y $ ×$ x $ )。$ z $ =−0.0065 mm平面的重建结果如图5所示,$ y $ =0.0065 mm平面的重建结果如图6所示,$ x $ =−0.0065 mm平面的重建结果如图7所示,图像的灰度窗设置为[−0.1,1.0]。可以看到FDK和DBP重建图像几何结构和对比度基本一致。表2给出了SheppLogan模体DBP重建和FDK重建的均方根误差(root mean square error,RMSE)和多尺度结构相似性(multi-scale structural similarity,MSSIM)的值。表 1 仿真参数Table 1. Simulation parameters模体像素/mm3 模体体素数量
($ z $×$ y $×$ x $)探测器像素/mm2 探测器数量
($ u $×$ v $)旋转角范围/° 旋转角采
样数量源探距
$ S{O}_{\mathrm{d}} $/mm源轴距
$ SO $/mm倾斜角
$ \alpha $/°0.0065 360×600×600 0.1376 2350×350 0−360 512 269.378 25.058 45 表 2 重建算法RMSE与MSSIM比较Table 2. Comparison of root mean square error (RMSE) and mean structural similarity (MSSIM) for reconstruction algorithms重建算法 DBP FDK RMSE 0.1517 0.1486 MSSIM 0.4108 0.4184 图8是
$ z $ =−0.0065 mm平面重建图像的150行和150列的剖面线,可以看出DBP重建结果与FDK非常接近。按照表1的参数仿真得到的投影数据可以覆盖整个shepplogan模体,我们将表1中的探测器数量从350×350缩小至256×256,其他参数不变,如图9所示,可以得到有明显截断的投影数据。以此截断投影数据进行FDK和DBP重建,其中DBP重建中取PI线与
$ y $ 轴平行。图10所示为$ z $ =−0.0065 mm平面的重建结果,图11所示为图10第150行剖面线,表3为图10中黄色框内区域的RMSE和MSSIM的值。可以看到FDK重建图像在左右两侧截断位置出现了明显的截断伪影,而DBP重建图像则无此明显伪影,DBP重建的黄色框内区域的MSSIM也明显高于FDK。这体现了DBP重建方法对于CL成像在处理截断问题时同样可以发挥一定的优势。表 3 黄色框内区域RMSE与MSSIM比较Table 3. Comparison of root mean square error (RMSE) and mean structural similarity (MSSIM) in the yellow box area重建算法 DBP FDK RMSE 0.1664 0.1662 MSSIM 0.5265 0.4581 2.2 实际电路板扫描数据重建
我们在图12所示的实际CL设备扫描了电路板,扫描参数如表4所示,用投影数据重建了100×
1000 ×1000 ($ z $ ×$ y $ ×$ x $ )的三维体,重建体素尺寸为0.0133mm3,重建中心为坐标原点$ O $ 。FDK方法和DBP方法的重建时间如表5所示。两种算法均使用NVIDIA GeForce RTX3060 Laptop GPU进行加速计算。与FDK方法相比,DBP方法的重建时间更长,主要原因在于DBP方法需要分别在冠状方向和矢状方向进行重建,然后再进行频谱混合。$ z $ =−0.0455 mm平面的重建结果如图13所示,$ y $ =−0.5915 mm平面的重建结果如图14所示,$ x $ =−1.9045 mm平面的重建结果如图15所示,图像的灰度窗设置为[−0.5,1.5]。可以看到FDK和DBP重建图像的结构基本相同,都可以清楚地看到气泡缺陷。由于CL设备获得的投影数据在每个视角下都有截断,不满足有限Hilbert逆变换条件,所以DBP方法的重建图像也存在截断伪影,但DBP重建方法的对投影的微分是局部操作,投影截断造成的边界光晕伪影比FDK弱。由于CL几何扫描数据的不完备性,两种重建方法的层间混叠伪影都比较严重。表 4 扫描参数Table 4. Scan parameters探测器像素/mm2 探测器数量($ u $×$ v $) 旋转角范围/° 旋转角采样数量 源探距$ S{O}_{\mathrm{d}} $/mm 源轴距$ SO $/mm 倾斜角$ \alpha $/° 0.1376 2750×750 0−360 256 256.904 23.898 45 表 5 实际投影重建时间Table 5. Actual projection reconstruction time旋转角采样数量 探测器数量($ u $×$ v $) 重建体素数量($ z $×$ y $×$ x $) DBP重建时间/s FDK重建时间/s 256 750×750 100× 1000 ×1000 30.84 9.87 3. 结论
本文推导了基于DBP的方形截面视野旋转CL的解析重建算法,为这种特殊CL成像提供了一种新的重建方案,我们用计算机仿真的投影数据和实际电路板的投影数据进行了重建方法的验证,并与FDK算法进行了比较。结果表明,与FDK算法相比,DBP重建方法获得的重建图像与之在结构上有很好的一致性,且在处理CL成像投影截断的问题上具有一定优势。在针对PCB的成像场景中,由于分辨率要求高且探测器面积有限,通常只能对感兴趣的局部区域进行成像,因此投影数据不可避免存在截断,这是成像系统物理特性的固有局限。我们提出的DBP方法正是为了有效应对这种截断投影的重建问题。目前的DBP方法中我们仅使用一种简单的方式估计PI线方向的投影,未来我们将继续研究提高PI线方向投影估计的准确性,以期获得比FDK方法更为准确的重建图像。此外,因为CL几何扫描视角有限,重建图像中会出现严重的层间混叠伪影。我们将研究利用深度学习方法解决这一难题。与FDK方法相比,以DBP方法作为深度学习方法的域转换算子,可以有效地减少截断伪影的影响,使深度学习网络能够更专注于学习层间混叠伪影的抑制方法,提高重建图像质量,为PCB等板状物体的检测提供可靠工具。
-
表 1 仿真参数
Table 1 Simulation parameters
模体像素/mm3 模体体素数量
($ z $×$ y $×$ x $)探测器像素/mm2 探测器数量
($ u $×$ v $)旋转角范围/° 旋转角采
样数量源探距
$ S{O}_{\mathrm{d}} $/mm源轴距
$ SO $/mm倾斜角
$ \alpha $/°0.0065 360×600×600 0.1376 2350×350 0−360 512 269.378 25.058 45 表 2 重建算法RMSE与MSSIM比较
Table 2 Comparison of root mean square error (RMSE) and mean structural similarity (MSSIM) for reconstruction algorithms
重建算法 DBP FDK RMSE 0.1517 0.1486 MSSIM 0.4108 0.4184 表 3 黄色框内区域RMSE与MSSIM比较
Table 3 Comparison of root mean square error (RMSE) and mean structural similarity (MSSIM) in the yellow box area
重建算法 DBP FDK RMSE 0.1664 0.1662 MSSIM 0.5265 0.4581 表 4 扫描参数
Table 4 Scan parameters
探测器像素/mm2 探测器数量($ u $×$ v $) 旋转角范围/° 旋转角采样数量 源探距$ S{O}_{\mathrm{d}} $/mm 源轴距$ SO $/mm 倾斜角$ \alpha $/° 0.1376 2750×750 0−360 256 256.904 23.898 45 表 5 实际投影重建时间
Table 5 Actual projection reconstruction time
旋转角采样数量 探测器数量($ u $×$ v $) 重建体素数量($ z $×$ y $×$ x $) DBP重建时间/s FDK重建时间/s 256 750×750 100× 1000 ×1000 30.84 9.87 -
[1] 曹大泉, 王雅霄, 阙介民, 等. 基于SART算法的CL硬化伪影校正方法研究[J/OL]. 原子能科学技术, 2014, 48(7): 1314-1320. DOI: 10.7538/yzk.2014.48.07.1314. CAO D Q, WANG Y X, QUE J M, et al. Research of beam hardening correction method for CL system based on SART algorithm[J/OL]. Atomic Energy Science and Technology, 2014, 48(7): 1314-1320. DOI:10.7538/yzk.2014.48.07.1314. (in Chinese).
[2] YANG M, WANG G, LIU Y. New reconstruction method for x-ray testing of multilayer printed circuit board[J/OL]. Optical Engineering, 2010, 49(5): 056501-056501. DOI: 10. 1117/1. 3430629.
[3] MYAGOTIN A, VOROPAEV A, HELFEN L, et al. Efficient volume reconstruction for parallel-beam computed laminography by filtered backprojection on Multi-Core clusters[J/OL]. IEEE transactions on image processing, 2013, 22(12): 5348-5361. DOI: 10.1109/TIP.2013.2285600.
[4] VOROPAEV A, MYAGOTIN A, HELFEN L, et al. Direct Fourier inversion reconstruction algorithm for computed laminography[J]. IEEE Transactions on Image Processing, 2016, 25(5): 2368-2378. DOI: 10.1109/TIP.2016.2546547.
[5] SUN L, ZHOU G, QIN Z, et al. A reconstruction method for cone-beam computed laminography based on projection transformation[J]. Measurement Science and Technology, 2021, 32(4): 045403.
[6] GORDON R, BENDER R, HERMAN G T. Algebraic reconstruction techniques (ART) for three-dimensional electron microscopy and x-ray photography[J]. Journal of Theoretical Biology, 1970, 29(3): 471, IN1, 477-476, IN2, 481. DOI: 10.1016/0022-5193(70)90109-8.
[7] GILBERT P. Iterative methods for the three-dimensional reconstruction of an object from projections[J]. Journal of Theoretical Biology, 1972, 36(1): 105-117. DOI: 10.1016/0022-5193(72)90180-4.
[8] ANDERSEN A H, KAK A C. Simultaneous algebraic reconstruction technique (SART): A superior implementation of the ART algorithm[J]. Ultrasonic Imaging, 1984, 6(1): 81-94. DOI: 10.1016/0161-7346(84)90008-7.
[9] HEBERT T, LEAHY R. A generalized EM algorithm for 3-D bayesian reconstruction from poisson data using gibbs priors[J/OL]. IEEE Transactions on Medical Imaging, 1989, 8(2): 194-202. DOI: 10.1109/42.24868.
[10] DOBBINSobbins J T, GODFREY D J. Digital X-ray tomosynthesis: Current state of the art and clinical potential[J/OL]. Physics in Medicine & Biology, 2003, 48(19): R65-R106. DOI: 10.1088/0031-9155/48/19/R01.
[11] LAURITSCH G, HAERER W H. Theoretical framework for filtered back projection in tomosynthesis[C]. Proceedings of SPIE: 1998: 1127-1137. DOI: 10.1117/12.310839.
[12] LEVAKHINA Y M, MVLLER J, DUSCHKA R L, et al. Weighted simultaneous algebraic reconstruction technique for tomosynthesis imaging of objects with high-attenuation features[J/OL]. Medical Physics (Lancaster), 2013, 40(3): 031106-n/a. DOI: 10.1118/1.4789592.
[13] ZhAO Y, XU J, LI H, et al. Edge information diffusion-based reconstruction for cone beam computed laminography[J]. IEEE transactions on image processing, 2018, 27(9): 4663-4675. DOI: 10.1109/TIP.2018.2845098.
[14] NOO F, CLACKDOYLE R, PACK J D. A two-step Hilbert transform method for 2D image reconstruction[J/OL]. Physics in Medicine & Biology, 2004, 49(17): 3903-3923. DOI: 10.1088/0031-9155/49/17/006.
[15] ZOU Y, PAN X. Exact image reconstruction on PI-lines from minimum data in helical cone-beam CT[J/OL]. Physics in medicine & biology, 2004, 49(6): 941-959. DOI: 10.1088/0031-9155/49/6/006.
[16] PACK J D, NOO F, CLACKDOYLE R. Cone-beam reconstruction using the backprojection of locally filtered projections[J/OL]. IEEE Transactions on Medical Imaging, 2005, 24(1): 70-85. DOI: 10.1109/TMI.2004.837794.
[17] KUDO H, SUZUKI T, Rashed E A. Image reconstruction for sparse-view CT and interior CT—introduction to compressed sensing and differentiated backprojection[J]. Quantitative Imaging in Medicine and Surgery, 2013, 3(3): 147.
[18] TANG S, TANG X. Radial differential interior tomography and its image reconstruction with differentiated backprojection and projection onto convex sets[J]. Medical Physics, 2013, 40(9): 091914.
[19] HAN Y, YE J C. One network to solve all ROIs: Deep learning CT for any ROI using differentiated backprojection[J]. Medical Physics, 2019, 46(12): e855-e872.
[20] ERNST P, ROSE G, NVRNBERGERA. Sparse view deep differentiated backprojection for circular trajectories in cbct[C]//Proceedings of the 16th Virtual International Meeting on Fully 3D Image Reconstruction in Radiology and Nuclear Medicine. 2021: 463-466.
[21] Zou X, Shi W L, Du M G, et al. A square cross-section FOV rotational CL (SC-CL) and its analytical reconstruction method. 2024.
[22] HAN Y, KIM J, YE J C. Differentiated backprojection domain deep learning for conebeam crtifact removal[J/OL]. IEEE transactions on medical imaging, 2020, 39(11): 3571-3582. DOI:10. 1109/TMI.2020.3000341.