CT(Computed Tomography)技术在生活中发挥着日益重要的作用。CT可以在不破坏样品的情况下, 利用样品对射线能量的吸收特性对生物组织和工程材料进行断层成像, 由此获取样品内部的结构信息。

目前在利用Matlab将扫描数据进行图像重建时需要用到iradon^[1]函数, 利用X射线接收器的位置与得到的衰减后的能量进行Radon逆变换得到原物质图像。康克军^[2]等人基于Radon变换导数对三维CT图像重建进行了研究, 对目前最新研究成果作了简要的介绍和概括。Cirrone^[3]等利用滤波反投影重建算法实现了高能量的质子束计算机断层成像。Kim^[4]等提出了结合有序子集和动量算法加速X射线CT图像重建的算法, 利用可分二次规划的有序子集方法具有可平行化, 有效地提高了算法的收敛速度。Pelt和Batenburg^[5]设计了一种与数据相关的滤波器以减少投影误差, 提高滤波反投影重建算法的重建质量。叶海霞^[6]研究了工业CT窄角扇形束滤波反投影并行图像重建, 将窄角扇形束扫描看作多探测器的平行束扫描, 将窄角扇形束重排成平行束, 再由平行束滤波反投影重建算法实现图像重建。阎春生^[7]等介绍了层析成像技术的图像重建算法, 并从正向问题数学模型的简化和反向问题数学模型的映射结构的角度比较了各种算法的特点和优劣。高红霞^[8]等为了解决不完备投影数据的重建问题而提出了基于粒子群优化的随机稀疏重建算法, 保证了算法效率, 也提高了重建质量和算法稳健性。张铁^[9]等基于Radon变换对改进的Fourier进行了误差分析。虽然CT图像重建算法有很大的改进, 但是利用Matlab中的iradon函数进行图像重建得到的图像位置以及形状等信息与原物质还存在一定的差距, 需要利用已知物质进行参数标定，从而确定重建图像位置相对原物质偏差状况, 吻合效果有待提高。

本文基于Radon变换, 利用包络线对Matlab中CT图像重建的过程进行优化, 准确计算出代还原物体的形状、角度及其位置。

1 Radon算法图像重建

一种典型的二维CT系统如图 1^[10]所示：平行入射的X射线垂直于探测器平面, 每个探测器单元看成一个接收点, 且等距排列。X射线发射器和探测器相对位置固定不变, 整个发射-接收系统绕某固定的旋转中心逆时针旋转180次。对每一个X射线方向, 在具有512个等距单元的探测器上测量经位置固定的二维待检测介质吸收衰减后的射线能量, 并经过处理后得到180组接收信息, 部分数据如表 1所示。经Matlab可视化后所构建的吸收信息图像如图 2所示。

射线经过物质的衰减规律^[11]为：

$ I = {I_0}{e^{ - \int\limits_L {u\left( {x,y} \right)dl} }} $

(1)

其中, I₀为初始射线强度, $\int\limits_L {} u\left( {x, y} \right)dl $为对射线l的积分, 在此情景中射线的角度随着时间变化而变化, 故考虑利用Radon变换^[12]对问题的求解域进行转换, 如图 3所示, 变换形式为：

$ {R_u} = \left( {p,\theta } \right)\int\limits_L {u\left( {x,y} \right)dl} $

(2)

其中(x, y)为积分直线上点坐标, $ {\vec r}$为旋转中心到积分直线上一点(x′, y′)的向量, $ {\vec p}$为旋转中心到积分直线的垂直向量, p为旋转中心到积分直线的距离, θ为(x, y)到旋转中心的连线与x轴正方形之间的夹角。直线l:y=kx+b可表示为：

$ p - x\cos \theta - y\sin \theta = 0 $

(3)

如果借助Delta函数, 上述积分还可写为：

$ \begin{array}{l}{R_{u}(p, \theta)=\int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} f d x d y} \\ {f=u(x, y) \delta(p-x \cos \theta-y \sin \theta)}\end{array} $

(4)

其中, R_u(p, θ)定义为沿由p和θ定义的直线的线积分。在R_u(p, θ)空间中, 利用Radon逆变换给出从投影重建的解u(x, y), 其中u(x, y)可表示为^[11]：

$ \begin{array}{l}{u(x, y)=\frac{1}{2 \pi^{2}} \int_{0}^{\pi} \int_{-\infty}^{+\infty} g d p d \theta} \\ {g=\frac{1}{x \cos \theta+y \sin \theta-p} \frac{\partial R_{u}}{\partial p}}\end{array} $

(5)

若要解出某定点的吸收率u(x, y), 首先要知道R_u(p, θ)具体函数表达式, 但由于只知道部分p与θ对应的R_u(p, θ)的离散值, 故可以考虑通过以下两种方式求解u(x, y)。

1) 对大量的R_u(p, θ)进行二元函数拟合^[13], 得到二元拟合函数R_u(p, θ)。此方式有以下缺点：一是导致最终的结果偏差太大, 二是拟合耗时过长。

2) 第二种方法退而求其次, 不再求解R_u(p, θ)的具体表达式, 而是利用已知的R_u(p, θ)值进行求解近似积分值^[14]：

$ \begin{array}{l} u(x,y) = \frac{1}{{2{\pi ^2}}}\int_0^\pi {\int_{ - \infty }^{ + \infty } h } dpd\theta \\ \approx \frac{1}{{2{\pi ^2}}}\sum\limits_{i = 1}^{511} {\sum\limits_{j = 1}^{179} k } \left( {{\theta _{j + 1}} - {\theta _j}} \right)\\ h = \frac{1}{{x\cos \theta + y\sin \theta - p}}\frac{{\partial {R_u}}}{{\partial p}}\\ k = \frac{{{R_u}\left( {{p_{i + 1}},{\theta _j}} \right) - {R_u}\left( {{p_i},{\theta _j}} \right)}}{{x\cos {\theta _j} + y\sin {\theta _j} - {p_i}}} \end{array} $

(6)

其中, p_i代表第i条X射线到旋转中心的距离, θ_i代表第i个旋转角, R_u(p_i, θ_j)为(p_i, θ_j)处的接收值, 由此可计算出指定点外的吸收率u(x, y)。

2 包络图像函数优化

为了在特殊条件下准确还原未知物质的图像, 在此引入物质的包络线——某一物质的外围轮廓线。下面将利用包络线对未知物质经扫描得出的数据进行图像精确重建。

2.1 切点处探测单元的确定

CT系统每一时刻扫描时的X射线对物体边缘均有两个切点：一个为X射线与物体上切面切线的切点, 另一个为X射线与物体下切面切线的切点, 如图 4中m、n两点。若数据中第i列为全零元素, 第i+1列不为全零元素, 可得此时刻第i个探测器发出的X射线与物体相切。同理可得另外一条与物体相切的X射线。其余时刻的情况以此类推, 可得到每一时刻发出相切射线的探测单元序号, 如表 2所示。

2.2 上下切线截距的确定

本例中由CT扫描系统旋转中心的坐标(-9.104 7, 5.793 9)能够得出经过旋转中心的X射线的函数公式为：

$ y-y_{0}=-\cot \theta_{i}\left(x-x_{0}\right) $

(7)

其中, (x₀, y₀)是CT旋转中心坐标, θ_i为X射线在竖直方向上的夹角, 如图 5所示。

由(7)式可得上下两条切线的方程为：

$ y=-\cot \theta_{i} x+\cot \theta_{i} x_{0}+y_{0}+\Delta c $

(8)

$ \Delta c=\frac{\Delta s}{\sin \theta_{i}} $

(9)

$ \Delta s_{1}=(256-1) \delta $

(10)

$ \Delta s_{2}=(256-512) \delta $

(11)

其中, Δc为两条与位置物质相切的X射线在竖直方向上的距离, Δs=Δs₁+Δs₂, Δs₁, Δs₂为第i条X射线到旋转中心的距离。

易求得180次照射各个切线相对于第256条射线的增加的截距Δc以及上下切线的截距, 如表 3所示。

由斜率与Δc易得上下切线方程, 依次分别联立旋转前后两组切线相邻方程, 得到两组切线交点, 连接交点进而得出图形的未知物质的包络线和位置。如图 6所示。

由于只求包络线的思想存在误差, 为了更加精确地对位置和形状进行确定, 此处采用二次取中点的方法进行校正：依次连接相邻两个交点, 求出各个交点组成线段的中间坐标；连接相邻中间坐标, 得到经过二次取中点进行校正后的包络线和位置。如图 7所示。

此包络图像确定了原未知物质的形状轮廓与位置, 为后续图像的准确还原提供了客观依据。

3 计算机仿真分析 3.1 传统方法图像还原

由X射线能量接收数据u(x, y)运用Matlab中的iradon函数得到初步重建图像, 如图 8所示。

由于这只是经过iradon函数得出的初步重建图像^[14], 并不知道其方向与位置的准确性, 下面进行验证。

将初步重建图像与原物质对比, 如图 8、图 9所示, 可见重建图像与原物质之间存在着角度上的误差。

若想将所得到的初步图像消除位置上的误差, 必须利用已知模板进行参数标定, 得到图像与模板之间的对应关系后再经过人工观察进行逐次实验, 从而将误差逐步减小^[11]。如图 10所示, 过程相对来说是很麻烦的。

由此可见, 仅经过iradon函数处理得到的重建图像与原物质的位置之间产生了较为明显的误差, 且误差需要人为观察后逐步进行修正。

3.2 运用包络图像

根据以上验证, 直接运用iradon函数进行图像重建, 会产生位置以及角度上的误差。根据图 7所示原未知物质的角度与位置的信息, 对初步重建图像的位置、形状等参数进行人工调整^[15], 最终将得到优化后的重建图像, 与原未知物质对比效果如图 11、图 12所示。

通过对比可以看出, 利用包络思想进行图像重建所得到的效果较为理想, 基本吻合原未知物质的形状、位置以及角度等特征信息。

综上, 在实际操作时, 可以先通过包络思想, 利用包络线近似计算图像的形状与位置, 可以简化图像重建流程, 进一步得到未知物质的准确还原图像。

4 结束语

本文阐明了利用Matlab中iradon函数处理CT系统接收数据所得重建图像在角度、位置等反面的不足, 并且提供了包络线这一概念, 来提前确定未知介质的形状和角度, 简化了分析过程, 在算法与模型上进行了优化。