基于对称重复的模板校正锥束CT系统几何失真的方法

摘要

一种基于对称重复的模板校正锥束CT系统几何失真的方法，属于图像处理领域里的计算机层析成像技术领域。其特征是本发明建立了探测器角度偏差与纹理图像的秩之间的联系，并将探测器角度偏差的求解问题等效为一定约束下纹理图像的秩最小化问题，根据求得的角度偏差对探测器进行调节，消除角度偏差后再利用特定的几何关系对探测器位移偏差进行求解。本发明的效果和益处是，利用含有纹理低秩性的模板进行CT成像系统的角度偏差求取，仅需要采集单角度投影图像，而位移偏差的求取也只需移动旋转台一次，获取变化前后的两幅投影图像，该方法模板制作容易，测量方法简单有效。

主权项

一种基于对称重复的模板校正锥束CT系统几何失真的方法，首先建立了探测器角度偏差与纹理图像的秩之间的联系，并将探测器角度偏差的求解问题等效为一定约束下纹理图像的秩最小化问题，通过求解秩最小化问题，得到CT成像系统中探测器的角度偏差，其特征在于，步骤如下：S1在所采集的投影图像中选取包含低秩纹理的一块矩形区域作为输入，记为I；将参数进行初始化，角度参数初始值θ＝0,β＝0,η＝0，归一化位移参数的初始值t′₁＝0,t'₂＝0,t'₃＝0；S2将I和初始化的转换参数τ作为初始条件带入，对公式(8)的优化问题进行求解：Iоτ运算被转换为：将局部最优解τ对应的投影变换作用在图像的每一个点上，这样便转化成了一幅新的图像；将所得的新图像重新赋值给I(u_ij(τ),v_ij(τ))，图像转换过程中利用公式(5)；$\begin{array}{c}{u}_2=\frac{\frac{D}{p_x}\left(r_{11}{u}_1+r_{12}{v}_1+\frac{D}{p_x·R}{t}_1\right)}{r_{31}{u}_1+r_{32}{v}_1+\frac{D}{p_x·R}·t_3}=\frac{\frac{D}{p_x}\left(r_{11}{u}_1+r_{12}{v}_1+\frac{D}{p_x}{t^{\prime }}_1\right)}{r_{31}{u}_1+r_{32}{v}_1+\frac{D}{p_x}·{t^{\prime }}_3}\\ v_2=\frac{\frac{D}{p_x}\left(r_{21}{u}_1+r_{22}{v}_1+\frac{D}{p_x·R}{t}_2\right)}{r_{31}{u}_1+r_{32}{v}_1+\frac{D}{p_x·R}·t_3}=\frac{\frac{D}{p_x}\left(r_{21}{u}_1+r_{22}{v}_1+\frac{D}{p_x}{t^{\prime }}_2\right)}{r_{31}{u}_1+r_{32}{v}_1+\frac{D}{p_x}·{t^{\prime }}_3}\end{array}---\left(5\right)$其中，利用探测器绕X轴、Y轴、和Z_W轴的旋转角θ,β,η表示旋转矩阵r，具体表达式为：$r=R_x{R}_y{R}_z=\left[\begin{array}{ccc}{r}_{11}& r_{12}& r_{13}\\ r_{21}& r_{22}& r_{23}\\ r_{31}& r_{32}& r_{33}\end{array}\right]$$R_x=\left[\begin{array}{ccc}1& 0& 0\\ 0& cos\theta & sin\theta \\ 0& -\sin \theta & \cos \theta \end{array}\right]$$R_y=\left[\begin{array}{ccc}cos\beta & 0& -\sin \beta \\ 0& 1& 0\\ sin\beta & 0& \cos \beta \end{array}\right]$$R_z=\left[\begin{array}{ccc}cos\eta & sin\eta & 0\\ -sin\eta & cos\eta & 0\\ 0& 0& 1\end{array}\right]$$r=\left[\begin{array}{ccc}\cos \beta \cos \eta & \cos \beta \sin \eta & -\sin \beta \\ -\cos \theta \sin \eta +\sin \theta \sin \beta \cos \eta & \cos \theta \cos \eta +\sin \theta \sin \beta \sin \eta & \sin \theta \cos \beta \\ \sin \theta \sin \eta +\cos \theta \sin \beta cos\eta & -\sin \theta \cos \eta +\cos \theta \sin \beta \sin \eta & \cos \theta \cos \beta \end{array}\right]---\left(2\right)$式(8)求解的具体步骤为：S2.1将式(8)中的目标函数进行凸优化，约束条件进行线性化，得到式(10)：其中▽I为雅克比矩阵，是一个三维张量，转换函数τ中各转换参数的增量是其中的元素，在线性化过程中，由于是用线性问题近似表示局部非线性问题，因此，在求解该问题时，需要先求解各局部解的最优值，然后进行迭代以求得全局最优值，使用迭代凸优化算法对该问题进行求解，步骤如下：S2.2计算雅克比矩阵首先要对图像进行归一化处理，$I\left(u_{ij}\left(\tau \right),v_{ij}\left(\tau \right)\right)\leftarrow \frac{I\left(u_{ij}\left(\tau \right),v_{ij}\left(\tau \right)\right)}{{||I\left(u_{ij}\left(\tau \right),v_{ij}\left(\tau \right)\right)||}_F},{▿}_{\tau }I\left(u_{ij}\left(\tau \right),v_{ij}\left(\tau \right)\right)\leftarrow \frac{\partial }{\partial \tau }\left(\frac{I\left(u_{ij}\left(\tau \right),v_{ij}\left(\tau \right)\right)}{{||I\left(u_{ij}\left(\tau \right),v_{ij}\left(\tau \right)\right)||}_F}\right);$S2.3利用交替方向迭代法求解下式：$\left(I^0,E,\Delta \tau \right)={\mathrm{argmin}}_{I^0,E,\Delta \tau }{L}_{\mu }\left(I^0,E,\Delta \tau ,Y\right)$其中L_μ(I⁰,E,Δτ,Y)＝f(I,E)+＜Y,为增广拉格朗日方程，Y表示拉格朗日乘数矩阵，μ＞0表示惩罚量，＜.,.＞表示矩阵的内积，f(I⁰,E)＝||I⁰||_*+λ||E||₁表示目标方程，R(I⁰,E,Δτ)表示约束条件的函数式，R(I⁰,E,Δτ)＝Iоτ+▽IΔτ‑I⁰‑E；具体求解步骤如下：S2.3.1初始化参数设置为：Y₀＝0,E₀＝0,Δτ₀＝0,μ₀＞0,ρ＞1,k＝0S2.3.2迭代循环，循环主要过程为使用交替方向法逐个求取各个变量的优化解：$I_{k+1}^0\leftarrow U_k{S}_{{\mu }_k^{-1}}\left[{\Sigma }_k\right]V_k^T$μ_k+1＝ρμ_k其中，▽I⁺是▽I的广义逆矩阵，S_μ[·]代表一个参数为μ的收缩算子，其一般操作为S_μ[x]＝sign(x)·(|x|‑μ)，μ表示收缩阈值的标量，为非负实数；判断目标函数值||Iоτ+▽IΔτ‑I⁰‑E||_F是否收敛，若不收敛，则用循环求得的参数值重新初始化变量，然后再进入循环求最优值，若收敛，继续执行下一步；S2.3.3更新变换参数τ＝τ+Δτ_k+1，判断目标函数值是否全局收敛，若不收敛，则返回S2.2继续执行，若收敛，继续执行下一步；S3输出低秩图像τ，随机噪声的局部最优解E＝E_K+1作为最终解；根据τ的计算结果得到探测器的角度偏差θ,β,η；根据求得的探测器角度偏差的值θ,β,η调节探测器，由于旋转顺序为先绕Z_W轴旋转，再绕Y轴旋转，最后绕X轴旋转，因此，再求得角度偏差后，将探测器先绕X轴旋转‑θ，再绕Y轴旋转‑β，最后绕Z_W旋转‑η，此时探测器的角度偏差已经调节完毕。