计算机断层成像(computed tomography, 简称CT)技术是通过对物体进行不同角度下的射线投影测量进而获取物体截面信息的成像技术^[1].CT重建算法以检测目标为中心, 通过旋转采集目标在不同视角下的投射数据, 利用图形重建算法还原出检测目标的断层图像^{[2, 3]}.锥束CT具有空间分辨率高、数据采集时间短、射线利用效率高等特点^{[4, 5]}, 其研发与应用已成为研究的热点.Tuy和Smith^{[6, 7]}给出了CT精确重建的充分条件:即每一个与被检测物体相交的平面上至少包含一个射线源点, 并且指出投影数据的采集和重建算法并不唯一.已经证明, 完整的投影数据可以唯一重建图像^{[8, 9]}.然而受检测环境、物体形状结构等因素影响, 某些角度下射线信号较弱甚至无法穿透, 破坏了投影数据的一致性和完整性, 使得采集到的投影信息部分缺失, 降低了物体的缺陷检测与识别能力^{[10, 11]}.所以, 不同结构的检测对象, 如何经济、高效地采集有效数据, 并对其快速、精确地重建成为解决不完全投影重建的关键.

Zou等人^[12]揭示了待重建图像和投影数据间的某种局部依赖关系, 当投影数据在水平方向有某种截断或投影线只覆盖部分区域时, 只要感兴趣区域数据完整, 图像依然可以重建^[13-15].Guo等人^[16]依据物体几何对称性对截断投影数据进行补齐完成重建.虽然Radon变换作为一种重要的数学工具, 被广泛应用于X射线晶体学、医疗辐射、核磁共振等不同领域, 但这并不意味着针对投影数据不全的有效重建方法已经找到.由于工业物体结构复杂, 可能出现部分位置的投影数据无法获取, 针对投影数据除了传统形式上的顺序采样和固定角度采样以外, 一些学者就已有的投影数据从访问形式上进行了研究.Kong等人^[17]提出了随机访问方式, Herman等人^[18]提出了质数分解法访问方式, Mueller等人^[19]提出了加权距离访问方式, Guan等人^[20]提出了多层访问方式.上述访问方式从重建速度及重建精度上也有自己的局限性.传统的诸如FDK、FBP^[21-24]解析算法在投影数据不完备时, 重建图像存在严重伪影, 而迭代类算法在该方面有着不可替代的优势.随着压缩感知理论的出现, 利用稀疏角度下的投影进行重建变为可能.Demircan^[25]根据少量投影结合离散代数重建算法得到较好的重建结果, Park^[26]、Hashemi^[27]、Bian^[28]将正则项加入到有限投影重建中, 应用范数最小化求解, 取得了欣喜的结果.

针对投影数据不足的重建图像质量问题, 现阶段的解决思路有两种:一是发展新的重建算法, 添加先验约束, 如单材质、低密度等到迭代过程中影响图像质量; 二是从投影数据本身入手研究投影幅数与重建质量之间的关系, 从数据完备性上影响图像.因此, 研究不完全投影基于迭代算法框架下的不同采样策略对重建图像质量的影响, 对不完备数据重建如何通过合理采样获取高质量重建图像有一定指导意义.

1 CT投影过程

工业CT的物理原理与医用CT基本一致, 医用CT通常采用低能射线.射线源电压最高可达80KV, 而工业CT所用的射线能量范围则大得多, 高达450KV.当一定能量的射线穿越物体时, 由于吸收和散射, 射线将产生衰减.研究指出, 射线穿过物质并与之发生作用后, 射线强度将受到射线路径上物质的吸收而衰减, 衰减规律由比尔定律确定^{[29, 30]}.基于圆扫描轨迹的CT投影采集过程, 透射得到的投影数据实际上是被测物体线性衰减系数沿穿越路径的积分^{[31, 32]}, 它是图像重建的原始输入.如图 1所示, 当射线沿某一路径l穿透物体时, 其入射强度I₀发生衰减变为I, 其作用过程可由Beer定律确定:

$I = {I_0}\exp \left( { - \int_l {f(x,y)} {\text{d}}l} \right)$

(1)

且路径l上的线积分等于入射强度与出射强度比的自然对数.所以, 当射线以不同方向和位置穿透检测物体时, 对应的所有路径上的积分构成的集合, 即为扫描过程所获取的投影.

对于给定的二维平面区域R², 设函数f(x, y)为重建图像, 给定射线l∈R², 射线方程由表达式(1)确定:

$x\cos \theta + y\sin \theta = \rho $

(2)

其中, θ表示射线与坐标系x轴的夹角, ρ表示射线到坐标原点的距离, 由射线路径与物体关系可知, 物体衰减系数沿穿越路径l的积分可表示为

$p = \int\limits_{ - \infty }^\infty {f(x,y)\operatorname{d} } l = \int\limits_{ - \infty }^\infty {f(r,\varphi )\operatorname{d} } l = \int\limits_{ - \infty }^\infty {f\left( {\sqrt {{\rho ^2} + {l^2}} ,\theta + \arctan \frac{l}{\rho }} \right)\operatorname{d} } l$

(3)

对沿着投影方向上的每条投影线的密度函数进行积分, 就得到了该射线上物体衰减的投影值^{[33, 34]}, 即该投影方向上的投影函数g_θ(ρ), 其反映了投影方向.如图 2所示.

Radon为CT技术建立了数学理论的基础, 认为:任何物体线性衰减系数分布可由其所有积分集合确定, 模型为

$f(x,y) = - \frac{1}{{2{{\mathtt{π}}^2}}}\mathop {\lim }\limits_{\delta \to 0} \int_\delta ^\infty {\frac{1}{\rho }\int_0^{2{\mathtt{π }}} {(x\cos \theta + y\sin \theta } } + \rho ,\theta ){\text{d}}\theta {\text{d}}\rho $

(4)

由此可知, 只要获取所有积分信息, 即可去求取反映物体断层内部结构和组成的数值信息.对单色窄束X射线而言, 物体中每个体素值将由唯一的衰减系数确定, 获取断层衰减系数分布是断层图像重建的最终目的.若己知μ(x, y, z), 则P可由式(5)得到, 称此问题为正问题, 称P为μ(x, y, z)的投影数据.而CT问题称上述问题的反问题, 即已知P求μ(x, y, z).在实际问题中, 通常是由P的一组采样值来计算μ(x, y, z)的近似值, 即由投影重建图像.图 3给出了投影采样重建流程示意图.

$p = - \ln \left( {\frac{I}{{{I_0}}}} \right) = \int_L {\mu (x,y,z)} {\text{d}}l$

(5)

传统CT重建算法需要完备的投影数据^[35], 然而受医学CT低剂量要求以及工业CT检测物体受射线源能量限制和物体几何形状(横纵比差异大, 结构扭曲异形)的限制, 致使投影正弦图数据缺失, 此时, 这些投影角度范围下的投影数据不能有效反映X射线经过物体的衰减程度.由不同的采样方式得到的投影数据信息有所不同, 因此重建结果各异, 图 4给出了相同物体在圆扫描轨迹下不同扫描策略得到的中间切片投影的正弦图像.该情况下采样得到的不完全投影数据分为稀疏采样投影、分段连续投影^[36-39].稀疏采样投影^[40-42], 即取不连续的若干个扫描角度, 相邻的旋转角之间存有间隔, 如图 4(b)、图 4(c)所示.分段连续角度投影^{[43, 44]}, 如图 4(d)、图 4(e)所示, 其扫描的旋转角限定在连续部分角度范围内.因此, 不完全投影数据采样方式有稀疏均匀采样、稀疏非均匀采样、均匀片段缺失采样、非均匀片段缺失采样等类型.理论上对不完全投影进行重建, 添加约束后问题存在唯一解, 但重建的结果不够稳定^[45].

2 采样策略与重建图像 2.1 范数最小化的迭代类算法

目前, 由投影重建图像问题主要有:解析法和迭代法, 其中解析法在完全投影的条件下可获得高质量的重建图像, 不完全投影条件下重建切片图像的质量较差, 然而迭代法可用于不完全投影条件下图像重建.

迭代类算法可以解决解析算法不能处理的离散数据问题, 其对不完全数据重建有着不可替代的优势^[46].诸如代数重建算法(algebraic reconstruction technique, 简称ART)^[47]、联合代数重建算法(simultaneous algebraic reconstruction technique, 简称SART)^[48]、联合迭代重建算法(simultaneous iterative reconstruction technique, 简称SIRT)^[49]等, 压缩感知(compressed sensing, 简称CS)理论^{[50, 51]}的出现, 将迭代类算法带入新的领域, 如基于CS框架的全变分正则化的重建算法(total variation augmented lagrangian alternating direction algorithm, 简称TVAL3).根据CS理论, 在迭代重建中增加关于稀疏变换后图像的稀疏性约束, 凭借其少于奈奎斯特采样数的投影观测值重建出原图像, 并且CS理论不受采样定理影响, 为解决少视角等不完备投影数据的重建问题提供了理论支撑.压缩感知理论的核心问题是信息重构, 即如何将数据从低维观测向量中恢复出来.其处理的模型如下:

$\left\{ \begin{array}{l} {\rm{minimize}}{\left\| {\Psi x} \right\|_p}\\ {\rm{s}}{\rm{.t}}{\rm{. }}\Theta x = b \end{array} \right.$

(6)

其中, x∈Rⁿ为待重建图像, b∈R^m为采样数据, Θ∈R^m×n为描述X射线投影测量的观测矩阵, ${\left\| \cdot \right\|_p}$表示向量的l_p范数.

由于l_p范数正则化的非光滑行对于求解优化问题具有一定的挑战^[52], 一般采用交替方向乘子法(alternating direction method of multipliers, 简称ADMM)求解优化问题.ADMM利用的是扩展拉格朗日函数来对数据进行迭代更新, 解决的模型如下:

$L(x,y,u) = \left\| x \right\|_2^2 + \frac{\rho }{2}\left\| {\Theta x - y} \right\|_2^2 + u(\Theta x - y)$

(7)

式中, u为对偶约束参数.ADMM算法在每次迭代时包括以下更新:

$\left. \begin{gathered} {x^{k + 1}}: = \arg \min {L_\rho }(x,{y^k},{u^k}) \hfill \\ {y^{k + 1}}: = \arg \min {L_\rho }({x^{k + 1}},y,{u^k}) \hfill \\ {u^{k + 1}}: = {u^k} + \rho (\Theta {x^{k + 1}} - {y^{k + 1}}) \hfill \\ \end{gathered} \right\}$

(8)

利用L(x, y, u)最小化问题和对偶可行性条件:$\frac{{\partial L}}{{\partial x}} = 0,\frac{{\partial L}}{{\partial y}} = 0,\frac{{\partial L}}{{\partial u}} = 0,$则主要迭代过程可表示为

$\left. \begin{matrix} {{x}^{k+1}}:={{(I+\rho {{\Theta }^{T}}\Theta )}^{-1}}(y+{{\Theta }^{T}}\mu +\rho {{\Theta }^{T}}(y-u)) \\ {{y}^{k+1}}:={{S}_{{\gamma }/{\rho }\;}}(\Theta {{x}^{k+1}}+{{u}^{k}}/\rho ) \\ {{u}^{k+1}}:={{u}^{k}}+\rho (\Theta {{x}^{k+1}}-{{y}^{k+1}}) \\ \end{matrix} \right\}$

(9)

其中, S_γ/ρ函数定义为

${{S}_{{\gamma }/{\rho }\;}}(x)={\rm{sgn}} (x)\cdot \max (\left| x \right|-{\gamma }/{\rho }\;,0)$

(10)

对于模型中参数λ、ρ的选取没有固定准则, 根据经验, 通常选取λ=0.2, ρ=0.1^[53].

由于不同的迭代算法模型各异, 性能不同, 重建结果也不同, 本文主要探讨投影采样策略对应的不同数据类型对重建质量的影响, 在算法本身对重建图像质量产生的影响方面未作比较.下文中不同采样模式的投影重建均采用基于CS构架的TVAL3重建算法^[54]来实现.

2.2 稀疏欠采样投影 2.2.1 均匀和非均匀采样策略

根据短扫描方式下重建条件, 只要保证扫描角度为180+γ(γ为锥形束CT的锥角)就可获得重建所需的完整数据^[55].因此, 在[0, 240°]范围内进行均匀和非均匀采样获取投影(如图 5所示), 对其数据进行重建.针对圆盘含噪声仿真模型, 在圆轨迹扫描下, 分别获取均匀采样的90幅、60幅、30幅投影信息及90幅、60幅投影下随机去掉30幅构成的非均匀采样投影等1~5这5种模式, 进行重建切片获取.

图 6所示为各采样策略下的仿真圆盘的sinogram信息分布.

按照上述采集方式得到的投影数据进行重建, 仿真参数如下:射线源到探测器距离600mm, 射线源到旋转中心距离500mm, 采样范围[0, 240°], 重建图像分辨率为256×256.根据真实的扫描投影背景区域的幅值变化, 仿真投影过程中增加N(μ=0, σ²=10000), 噪声幅值约400dBw的高斯噪声, 以下仿真投影设置一致.

图 7给出了圆盘的重建结果, 其中, 图 7(a)~图 7(e)为多孔圆盘对应采样策略1~策略5, 即均匀采样90幅、均匀采样60幅、稀疏非均匀采样60幅、均匀采样30幅、稀疏非均匀采样30幅的含噪声投影不同层切片信息, 图像灰度显示窗口为[0, 0.035].图 7(a1~a3)~图 7(e1~e3)为第120层切片图像对比, _1对应切片的横断面切片(X-Y), _2对应切片的矢状面切片(Y-Z), _3对应切片的冠状面切片(X-Z).显然, 在均匀采样90幅、60幅、30幅的重建结果中投影数越多, 重建结果越好.将图 7(b1~b3)与图 7(c1~c3)、图 7(d1~d3)与图 7(e1~e3)图像进行对比, 可以看出, 相同采样投影数情况下, 均匀采样重建结果要比非均匀采样的重建结果好很多, 在轮廓边缘上效果明显.图 8给出了含噪声投影重建图像沿直线的灰度分布, 其中均匀采样与非均匀采样重建结果有明显的不同.

图 9给出了不同幅数策略下多材质圆盘均匀与非均匀投影的重建结果, 图像灰度显示窗口为[0, 0.023].图 9(a)为均匀采样90幅的重建结果, 图 9(b)、图 9(c)分别为非均匀采样120幅、150幅的重建结果.同样地, _1对应切片的横断面切片(X-Y), _2对应切片的矢状面切片(Y-Z), _3对应切片的冠状面切片(X-Z).结果表明, 均匀采样策略重建图像在少于非均匀采样投影角度的情况下, 图像质量依然较高.因此, 获取数据时应优先选择均匀采样策略, 对于某些情况下无法进行均匀采样的应尽可能使采样角度趋于均匀分布.由于投影空间的对称性, [0, 360°]范围的情况完全类似, 这里不再赘述.

为了验证投影方式对重建切片影响规律的一致性, 对真实的铝质工件#A1获取投影.实验数据的CBCT系统为德国Yxlon公司的Y.TU 450-D02, 探测器为Varian公司的PaxScan 2520, 投影采集过程中扫描电压为200KV, 曝光量0.4mA·s, DSO为925.041 7mm, DOD为286.925 1mm, 锥形束CT的锥角为5°, 圆周扫描采集360幅投影, 在[0, 240°]范围内均匀选取60幅、30幅投影和非均匀的90幅、60幅投影进行重建, 重建图像分辨率为512×512.

图 10为真实铝件投影重建结果, 图像灰度显示窗口为[0, 0.03].图 10(a)~图 10(d)分别给出了试件及其不同位置切片.其中, 图 10(a1)~图 10(a4), 图 10(b1)~图 10(b4)和图 10(c1)~图 10(c4)分别为第50层、第100层和第140层切片对比, 依次对应均匀60幅、均匀30幅投影、非均匀90幅和非均匀60幅投影.对比后不难发现, 相同采样投影数情况下, 均匀采样重建结果比非均匀采样的重建结果要好, 比如, 图 10(_1)的60幅均匀投影与图 10(_4)的60幅非均匀投影.而且均匀采样策略重建图像在少于非均匀采样投影角度的情况下, 图像质量依然较高, 比如图 10(_1)的60幅均匀投影与图 10(_3)的90幅非均匀投影.因此, 真实投影数据重建结果与仿真得到的结论相一致.

2.2.2 奇偶投影数采样策略

虽然迭代类算法可以有效对投影数据进行重建, 但基于算法本身的属性, 需要更多的投影数据补充到迭代过程中来修正重建带来的伪影对图像质量的影响.圆轨迹采样奇偶数对重建结果也存在影响, 合理的采样策略对重建结果至关重要.对于[0, 180°]范围内的采样幅数来说, 无论奇数投影数还是偶数投影数, 对应在[0, 360°]范围内均为偶数情况, 所以仅讨论0~360°范围内的情况.这里仅对含噪声投影图像进行讨论, 噪声幅度与前面的设置一致.图 11(a)所示为在[0, 360°]范围内均匀获取2N(N=1, 2, 3, …)幅含噪声投影, 图 11(b)所示为在[0, 360°]范围内均匀获取2N-1(N=1, 2, 3, …)幅含噪声投影.

对本节奇偶投影采样投影进行重建, 图 12给出了多孔圆盘和Phantom头模型奇/偶数投影的切片结果, 图 12(a)、图 12(b)对应不同模型多层重建结果.重建图像分辨率为512×512, 图像灰度显示窗口为[0, 0.035].其中, 图 12(A1)~图 12(F1)表示偶数幅投影重建, 分别对应均匀采样60幅、30幅、10幅; 图 12(a1)~图 12(f1)表示奇数幅投影重建, 分别对应均匀采样59幅、29幅、9幅, 第一、二列表示多孔圆盘, 第三、四列表示Phantom头模型.图 13给出了第128层切片沿直线CD/EF的灰度变化曲线对比情况, 其中, 图 13(a)~图 13(c)对应多孔圆盘, 图 13(d)~图 13(f)对应头模型.可以看出:当采样数较多的情况下, 奇、偶数投影重建结果差别较小, 但是奇数采样还是优于偶数采样, 随着采样幅数的减少, 奇数投影数重建的结果在边界轮廓上较偶数投影数要好, 内部点重建质量平滑.原因在于偶数投影数情况下, 镜像采样位置的射线与物体的权值关系关于两位置中心连线对称, 即约束条件受限, 而奇数投影数可以避免共轭位置存在的问题.

对重建结果进行整体评价, 评价方法采用均方误差(mean square error, 简称MSE)^[56]、峰值信噪比(peak signal to noise ratio, 简称PSNR)^[57]、归一化平均绝对距离(normal average absolute distance, 简称NAAD)、归一化均方距离(normalized mean square distance, 简称NMSD)^[58].NAAD反映某些点产生较大误差的情况, NMSD反映许多点均有小误差的情况, 两者很相似, 其距离越大表明重建图像的偏差越大.表 1给出了数值结果.

由图 12、图 13和表 1可以看出, 圆周扫描下的奇、偶数投影, 随着投影数目的增加, 图像间的差别逐渐降低, 由于增加信息量相当于重建时的约束增加, 致使重建结果更接近真实图像.但就内部信息而言, 奇数幅采样的重建结果优于偶数幅采样.由其采样空间对称性可知, [0, 360°]范围内的2N(N=1, 2, 3, …)幅采样在α和180+α (α∈[0, 180°])的位置射线与体素相交权值关系相同, 约束出现了重复, 相比之下, 相同水平下的2N-1(N=1, 2, 3, …)奇数幅采样, 体素对射线的贡献关系约束量多, 致使内部重建结果更优.而对于轮廓而言, 边界信息的保持需要有切向射线经过, 而没有切向射线经过的边缘相对模糊^{[59, 60]}.

同样, 图 14所示为真实投影重建的结果, 其中, 图 14(a)对应钛材工件#A2, 图 14(b)对应铝质工件#A1.图 14(a1)~图 14(a4)依次对应30幅、29幅、10幅和9幅投影重建的第128层切片结果, 图 14(b1)~图 14(b4)、图 14(c1)~图 14(c4)依次对应30幅、29幅、10幅和9幅投影重建的第60层和第128层切片结果, 重建图像分辨率为512×512, 图像灰度显示窗口为[0, 0.03].

从真实投影重建结果来看, 根据图 13的真实投影重建结果, 图 14(a2)、图 14(b2)、图 14(c2)对应的29幅投影重建的内部灰度信息比图 14(a1)、图 14(b1)、图 14(c1)对应的30幅投影重建结果要好.同样, 图 14(a4)、图 14(b4)、图 14(c4)对应的9幅投影重建结果显然比图 14(a3)、图 14(b3)、图 14(c3)对应的10幅投影数据轮廓结构好很多, 且内部灰度分布较均匀.上述结果表明, 投影数据的奇偶性对真实数据同样有效, 与仿真结论相符.

2.3 分段连续投影

基于工业上零件的高密度性局部角度射线无法穿透的问题以及医学上人体辐射的低剂量问题, 通常导致信息获取时局部范围采样信息缺失, 继而出现投影信息分段.分段连续信息采样往往是物体投影无法获取情况下的数据结构^[61-66], 如图 15所示, 可以根据投影数据缺失的位置分布将采样策略分为两种类型.

均匀分段连续投影的过程是采样点和非采样点均匀连续分布的叠加, 如图 15(a)所示.非均匀分段连续投影的过程是采样点和非采样点非均匀连续分布的叠加, 如图 15(b)所示.为了说明不同采样策略的优劣性, 在[0, 360°]范围内每间隔15°进行连续采样(每幅/3°), 作为构成均匀分段连续的投影类型; 在[0, 360°]范围内, 不同间隔获取连续投影, 构成非均匀分段连续(每幅/3°).表 2给出了不同数据采集策略, 图 16为按照表 2模式得到的心脏模型含噪声sinogram图像, 其中, 图 16(a)~图 16(d)分别对应表 2中的模式1~模式4.

对获得的叶片和心脏模型分段连续投影进行重建, 图 17给出了不同高度处的切片结果, 图像显示窗口为[0, 0.045].

图 17(a)~图 17(d)分别给出了叶片和心脏模型不同位置处的切片图像.图 17(a1)~图 17(a4)、图 17(b1)~图 17(b4)分别对应叶片第80层和第140层不同分段形式投影的重建结果; 图 17(c1)~图 17(c4)、图 17(d1)~图 17(d4)分别对应心脏第95层和第130层不同分段形式投影的重建结果.可以看到, 均匀分段连续投影重建结果中(如图 17(a1)~图 17(d1)、图 17(a2)~图 17(d2)所示), 缺失角度较少情况下重建图像伪影较少, 缺失角度较多情况下伪影较多, 且投影缺失位置伪影严重, 轮廓结构性差.非均匀分段连续投影重建结果中(如图 17(a3)~图 17(d3)、图 17(a4)~图 17(d4)所示), 同样在投影缺失位置伪影严重, 边缘轮廓差, 且投影缺失越多情况下质量越差.

上述分段连续投影采样策略重建结果均在缺失投影较少情况下伪影不是特别明显, 但是局部有偏差; 当缺失投影较多时, 重建结果均较差.由于缺失角度对应位置无投影, 该角度位置的约束力弱, 其仅由其他采样位置的射线来补偿, 即采样位置的数据反馈给缺失位置体素的信息过少, 因此局部轮廓位置重建结果较差.

2.4 结果分析

迭代重建算法相比解析算法最主要的区别就是将重建图像模型化.迭代算法中, 重建图像被离散化, 因此迭代重建算法的实质就是解线性方程组, 其将真实的成像几何结构与成像物理效应模型化.通过对上述不同类型的采样策略下获得的投影信息应用迭代算法进行重建, 得到了不同质量的重建图像.

从图 7、图 9所示稀疏均匀和非均匀采样的重建结果中不难发现, 均匀采样重建结果要比非均匀采样的重建结果好很多, 且轮廓边缘效果明显.由投影数据的访问方式(即扫描得到的投影数据按投影角度或其他规则排列次序)可知, 迭代重建算法以采样次序所得投影进行反投影.然而, 迭代类算法的重建性能依赖于系统矩阵、投影序列等因素较大, 且不完全投影数据所带来的最大问题是使采样得到的线性方程组高度欠定, 因此带来的投影数据噪声严重, 伪影较大.

分析认为, 数据采样过程实际上是射线与重建体素相交的过程, 射线的重建体素之间的观测关系可以由一系列线性方程组表示, 如图 18所示, 模型如下:

$\left\{ \begin{gathered} {\omega _{11}}{x_1} + {\omega _{12}}{x_2} + {\omega _{13}}{x_3} + ... + {\omega _{1n}}{x_n} = {p_1} \hfill \\ {\omega _{21}}{x_1} + {\omega _{22}}{x_2} + {\omega _{23}}{x_3} + ... + {\omega _{2n}}{x_n} = {p_2} \hfill \\ {\omega _{31}}{x_1} + {\omega _{32}}{x_2} + {\omega _{33}}{x_3} + ... + {\omega _{3n}}{x_n} = {p_3} \hfill \\ \begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {}&{}&{} \end{array}}&{}&{}& \vdots \end{array} \hfill \\ {\omega _{m1}}{x_1} + {\omega _{m2}}{x_2} + {\omega _{m3}}{x_3} + ... + {\omega _{mn}}{x_n} = {p_m} \hfill \\ \end{gathered} \right.$

(11)

由于迭代算法中投影方法与重建模型中的权系数w_ij以射线与体素相交长度来体现, 所以相邻角度下的射线束所对应的方程彼此相关, 相关性的大小取决于相邻射束相交体素的相关程度, 如图 18所示.对于均匀采样来说, 采样位置反映出的射线与体素的权值关系相近, 而对于非均匀采样, 投影角度在某些范围内出现聚集现象, 其相关性较大, 使得均匀采样策略下的均匀相关性失衡.即均匀采样策略使每个待求体素与其他体素的关联相当, 而非均匀采样策略下, 局部射线方程聚集, 在某种程度上只是部分体素发生关系, 体素对射线的贡献出现不均, 缺失角度下体素之间的模型关系丢失, 导致体素之间关系失衡.因此, 非均匀策略下的投影重建在图像局部边缘轮廓处存在啃食现象, 且内部灰度分布不均匀.

整个影响过程可以由迭代算法的思想来反映:对于给定的投影数据, 将所得方程组看成是N个超平面, 在给定初始图像估计值x⁽⁰⁾的条件下, 将其投影到第1个方程所代表的超平面, 得到x⁽¹⁾, 再将x⁽¹⁾投影到第2个超平面, 得到x⁽²⁾, 以此类推, 则将x^(i–1)投影到第i个方程所表示的超平面时得到x⁽ⁱ⁾, 当投影进行到最后, 得到x^(M), 如图 19所示, 其表达如下:

$x_j^{\left( k \right)} = x_j^{\left( {k - 1} \right)} + \frac{{\left( {{p_i} - \sum\limits_{n = 1}^N {x_{}^{(k - 1)}{w_{in}}} } \right)}}{{\sum\limits_{n = 1}^N {w_{in}^2} }}{w_{ij}},j = 1,2,...,N$

(12)

从表达式(12)可以看出, 第i条射线遍历的体素迭代的过程以残差最小化的方式在逼近.对于没有进行投影的位置, 没有残差更新.而非均匀采样情况下, 不同位置体素更新的次数不同, 缺失投影位置体素迭代更新较少, 所以相对采样点来说, 缺失投影位置伪影比较严重.

图 20给出了含有两个超平面的二维情形, 显然, 超平面夹角越大, 收敛得越快(如图 20(a)和图 20(b)所示).当超平面互相垂直时, 收敛得最快(如图 20(c)所示).虽然收敛变快, 但是图像的细节信息忽略得较多.N维空间中, 第i个和第j个超平面构成的角度可以表示为

$\cos \theta = \frac{{{\omega _i} \cdot {\omega _j}}}{{\left| {{\omega _i}} \right| \cdot \left| {{\omega _j}} \right|}}$

(13)

其中, ω_i和ω_j表示采样矩阵向量, 重建模型中, ω_ij被取为射线与网格相交长度.投影过程中, 相同角度下的相邻射线束对应的超平面方程彼此相关, 相关性的大小取决于与这两条相邻射束同时相交的像素的相交程度.一般情况下, 相邻角度射线束所对应的方程彼此相关性较大, 则出现局部轮廓约束信息不足, 它们之间的角度越大, 相关性越小.所以, 为了减小被访问投影角度的相关性, 采样的投影角度应该在视角范围内尽可能均匀地分布, 且投影角度避免出现聚集.

根据上述图 14所示铝件切片结果, 图 21(a1)、图 21(a2)、图 21(b1)、图 21(b2)中的切片局部细节信息, 奇数幅投影重建比偶数幅投影重建所得图像边缘信息较好, 内部灰度也较均匀.由于奇偶采样位置的特殊性因素, 偶数幅投影采样(如图 21(a)所示), 共轭位置角度的投影信息相关性较大, 如图中黑色阴影部分, 其所含约束信息相对减少; 奇数数幅投影采样(如图 21(b)所示), 不存在共轭采样位置, 投影相关性较小, 包含的约束信息会较多, 因此结果较好.同样地, 根据图 17所示切片结果, 图 21(c1)、图 21(c2)、图 21(d1)、图 21(d2)中的切片局部细节信息由于投影信息聚集, 星状伪影引起的灰度分布均匀性差是制约轮廓信息的关键.对于分段连续采样而言(如图 21(c)、图 21(d)所示), 投影的缺失位置导致了约束信息较少, 使得该位置轮廓存在偏差.

基于上述分析, 针对不同采样策略获得的投影信息必须发展新的算法, 进而改进边缘轮廓和完成内部灰度更新.从而使得在相同投影数据结构下进行重构, 实现数据挖掘.

3 总结与讨论

本文主要针对有限角度类型下的不完全投影数据重建结果探讨采样策略对其质量的影响, 研究了均匀采样和非均匀采样情况下投影数据的特点.分别对仿真投影数据和真实投影数据按照不同采样策略进行获取, 并结合迭代类算法的优势对不同采样策略下的离散投影进行重建, 分析了不同投影类型对重建质量的影响.结果表明, 圆轨迹扫描下的均匀采样策略较非均匀采样策略重建质量要好, 且少量均匀数据比多非均匀数据结果更优; 再者就是对于满足重建要求的投影采集来说, 奇数幅投影相对偶数幅投影重建结果要好一些.因此, 对于数据采样, 可以得出以下总结:圆周采样时, 应优先选择均匀采样策略; 对于无法实施均匀采样策略的应尽可能地使采样角度逼近均匀分布; 优先选择奇数投影数进行重建.同时, 为了减小被访问投影角度的相关性, 采样的投影角度应该在视角范围内尽可能均匀分布, 且投影角度避免出现聚集.实验结果包括了不同结构类型的仿真对象和真实铝质工件#A1、钛质工件#A2扫描对象, 结果和分析表明了该结论的可靠性和针对不同对象具有的广泛适用性.所以, 在现有的大量文献中, 学者们研究着重于低剂量稀疏采样重建^[67-70]算法对重建结果的影响, 而针对本文涉及的分段连续采样形式对应的投影数据而言, 如何发展新的重建算法, 成为后续研究的着落点.比如, 可以在缺失位置进行投影差值, 构建均匀投影分布, 或者针对缺失位置投影进行边界约束和投影数据挖掘等, 都将是重建算法研究的方向.

基于工业对象异形弯曲等特点, 某些角度范围内X射线无法有效穿透物体导致该位置投影缺失, 出现投影数据分段, 所以研究有效算法、降低由于投影数据结构引入的伪影影响, 无疑有着更广的应用价值.当然, 投影信息的保真度是后续高质量重建的保障, 例如, 重建之前对投影的去散射^{[71, 72]}、射束硬化处理等.

通过对稀疏采样、奇偶性采样、分段连续采样等不同采样策略下的重建结果进行分析探讨, 为相关研究学者建立了圆轨迹投影扫描策略的整体框架, 方便相关技术人员和学者结合具体的检测对象选择合适的采样方法, 并为其算法研究提供依据.借鉴压缩感知框架下的数据压缩与挖掘思想, 为不完全投影数据重建高质量图像等研究提供突破点, 充分利用迭代类算法的优势, 推动数据恢复算法与采样方式的进程.