CN1300938A - X－ct有限角投影数据图象重建方法 - Google Patents

Abstract

X－CT截断投影数据图象重建方法是把信号表示成奇异函数的加权和。从采集到的X－CT有限角截断投影数据中提取用于重建图象的特征信息,并用这些特征信息构造截断奇异函数,再由这些截断奇异函数的线性泛函来重构X－CT图象,达到消除X－CT有限角投影数据图象重建中的截断伪影之目的。不论含噪声的还是无噪声的、对实际的还是对计算机模拟的X－CT有限角投影数据的实验结果,都表明这种方法效果高于现有方法,效果显著。

Description

X-CT有限角投影数据图象重建方法

本发明属于X-CT投影数据图象重建技术领域，涉及截断奇异函数分析的图象重建方法。

计算机析层成象(CT)已广泛地应用于医学成象、地震探测、材料非损伤检测等各种领域。CT成象效果常常依赖于采集CT投影数据条件与理想条件的差距。例如，投影数据和旋转角的采样。在工业应用中，常常遇见被CT检查对象的尺寸比CT扫描系统检测器的宽度大得多，只能用有限角投影数据重建图象。这种有限角投影数据用现有方法重建的图象总有伪影。旋转角只能分布在φ₁～φ₂之间，不能均匀分布于[0,π)，如附图1中的(b)图所示。虽然，在φ₁～φ₂之间可以分割成足够多的角度，投影数据也可以采样成足够多的值，用于投影图象重建。但这样构成的图象重建方程组往往是严重病态的和欠定的，不满足X-CT投影数据重建图象的完备性条件，得到的图象有严重伪影。完全投影数据集的频域分布如附图1中(a)图所示，有限角在φ₁～φ₂的截断投影数据集如附图1中的(b)图所示。目前国内外解决有限角截断投影成象中伪影问题主要且较有效的方法是最大熵和最小交叉熵多目标优化迭代重建方法，但效果都很不理想。

本发明的目的是针对现有的技术问题。应用本方法于有限角投影数据图象重建领域，在消除有限角截断投影成象中伪影的精度和速度方面比现有的方法有重大突破。

为实现上述目的，本发明技术构思关键是按信息论中的信息结构理论，函数可表示为某个函数系的函数的泛函的原理，建立图象的结构模型，从X-CT有限角投影数据中提取的模型参数，再根据结构模型，推出缺损投影数据，最后重建CT图象，达到消除截断伪影，准确重建图象。

本发明的实质性内容

1、截断奇异函数分析数学模型

众所周知，可用泰劳级数、麦克劳林级数、小波级数、傅里叶级数等去代表一个函数。但它们目的都是用来逼近这个函数，求取近似函数值。在此，我们希望找出一种离散函数系中的函数的加权和来表示一个任意的离散信号序列，且有最少个数的非零加权系数。这样，我们就能用尽量少代码表达这个离散信号序列，表示这个离散信号序列的信息。

设一有限长离散信号f(n),n=0,1,...,N-1的完整频谱为 $\hat{f}\left(u\right),u=\mathrm{0,1},...,N-1;{\hat{R}}_{u_1-u_2}\left(u\right)$ 是单位离散矩形函数，即 ${\hat{R}}_{u_1-u_2}\left(u\right)={\{}_{0,u\∉\{u_1,u_1+1,...,u_2-1\}}^{1,u\∈\{u_1,u_1+1,...,u_2-1\}}$ 其中，u₁是截断频谱的起始频率，u₂是截断频谱的终止频率。则截断频谱信号 $f_{u_1-u_2}\left(n\right)$ 可表示为： $f_{u_1-u_2}\left(n\right)=\mathrm{IDFT}\[\hat{f}\left(u\right){\hat{R}}_{u_1-u_2}\left(u\right)\]$ 即 $f_{u_1-u_2}\left(n\right)=f\left(n\right)\⊗r_{u_1-u_2}\left(n\right)\·\·\·\·\left(1\right)$ 这里，IDFT[·]表示逆离散傅里叶变换，是圆周卷积符号， $r_{u_1-u_2}\left(n\right)=\mathrm{IDFT}\[{\hat{R}}_{u_1-u_2}\left(u\right)\],$ 即 $r_{u_i-u_2}\left(n\right)=\frac{1}{N}\frac{e^{j2\mathrm{\π}{u}_{1}n/N}-e^{j2\mathrm{\π}{u}_2/N}}{1-e^{j2\mathrm{\πn}/N}}$ 离散信号的截断频谱信号重建问题可以表述为从 $f_{u_1-u_2}\left(n\right)$ 中恢复出原来的信号f(n),n=0,1….,N-1。在中国发明专利申请号为98107580.0的专利申请文件中，我们定义了奇异空间{W_k(n),k,n为整数}，其中W_k(n)为 $W_k\left(n\right)=\underset{l=0}{\overset{n}{\mathrm{\Σ}}}\mathrm{\δ}(l-k),n,k=0,1,...,N-1,$ 其中，δ(n)为离散单位脉冲函数。任何有限长离散信号f(n),n=0,1,…,N-1可表示为奇异空间{W_k(n),k,n∈0,1,…,N-1}中的线性泛函。 $f\left(n\right)=\underset{k=0}{\overset{N-1}{\mathrm{\Σ}}}{a}_k{W}_k\left(n\right)\·\·\·\·\left(2\right)$ 对上式两边同时进行 $r_{u_1-u_2}\left(n\right)$ 的圆周卷积得： $f_{u_1-u_2}\left(n\right)=\underset{k=0}{\overset{N-1}{\mathrm{\Σ}}}{a}_k{W}_k\left(n\right)\⊗r_{u_1-u_2}\left(n\right)\·\·\·\left(3\right)$ 我们定义 $W_{k,u_1-u_2}\left(n\right)=W_k\left(n\right)\⊗r_{u_1-u_2}\left(n\right)$ 为截断奇异函数。诸截断奇异函数 $\left\{W_{k,u_1-u_2}\left(n\right)\right|k=0,1,...,N-1\}$ 构成N维截断奇异函数空间。式(3)表明：截断频谱信号 $f_{u_1-u_2}\left(n\right)$ 可以表达为N维截断奇异函数空间中线性泛函数。即 若信号中仅有Q个奇异点(在a_k,k=0,1,…,N-1中仅有Q个不为0)，则截断频谱信号 $f_{u_1-u_2}\left(n\right)$ 可由 $\left\{W_{k,u_1-u_2}\left(n\right)\right|k=0,1,...,N-1\}$ 中对应的Q个截断奇异函数的线性泛函表示。即 $f_{u_1-u_2}\left(n\right)=\underset{i=1}{\overset{Q}{\mathrm{\Σ}}}{a}_{\mathrm{bi}}{W}_{\mathrm{bi},u_1-u_2}\left(n\right),a_{\mathrm{bi}}\≠0,i=1,...,Q\·\·\·\·\left(5\right)$ 其中 $W_{k,u_1-u_2}\left(n\right)$ 是以bi位置为奇异点的截断奇异函数，a_bi是 $W_{k,u_1-u_2}\left(n\right)$ 的加权系数。当u₁=0,u₂=N-1时，显而易见应有f_0-(N-1)(n)=f(n)，这时式(5)就变为：

所以，获得信号的奇异点及其奇异谱分析加权系数就等价于获取信号的频谱数据。从 $f_{u_1-u_2}\left(n\right)$ 重建f(n)的关键问题是由 $f_{u_1-u_2}\left(n\right)$ 求出奇异点bi和相应的加权系数a_bi,I=1,2…,Q。2、模型参数的求取

根据投影定理，X-CT投影数据p_y′(x′)与g(x,Y)的频谱G(ω,φ)之间存在以下关系： $G(\mathrm{\ω},\mathrm{\φ})={\∫}_{-\∞}^{\∞}\mathrm{Py}\′(x\′)e^{-\mathrm{jaxx}\′}\mathrm{dx}\′---\left(7\right)$ 对于有限角范围投影数据集，我们可以由投影定理获得有限角范围内的CT截断频谱(简称为有限角截断频谱)，如附图1中的(b)图所示。即，X-CT有限角截断投影图象重建问题可以转变为有限角截断频谱的图象重建问题。

对于有限角截断频谱的图象重建问题，由于信号能量损失严重，小波系数误差太大，故不能用类似磁共振截断频谱图象重建方法的思路，用小波变换方法预选信号的奇异点。为此，我们根据奇异函数分析理论，运用层析法筛选信号的奇异点，确定奇异谱分析的加权系数。

(1)奇异点检测：我们把 $f_{u_1-u_2}\left(n\right)$ 看成周期为N的周期函数，则其差分函数 $d_{u_1-u_2}\left(n\right)$ 为

由于 $\mathrm{\δ}(n-\mathrm{bi})\⊗u_{u_1-u_2}\left(n\right)=r_{u_1-u_2}(n-\mathrm{bi})$ ，即 $d_{u_1-u_2}\left(n\right)$ 立刻成为： $d_{u_1-u_2}\left(n\right)=\underset{i=1}{\overset{Q}{\mathrm{\Σ}}}{a}_{\mathrm{bi}}{r}_{u_1-u_2}(n-\mathrm{bi})\·\·\·\·\left(8\right)$ 式(8)表明， $d_{u_1-u_2}\left(n\right)$ 是Q个 $r_{u_1-u_2}\left(n\right)$ 分别位于bi,i=1,2,…,Q,上的函数的加权和。当u₂=N-u₁，0＜u₁＜N/2时, $r_{u_1-u_2}(n-\mathrm{bi})$ 可写为 $r_{u_1-u_2}(n-\mathrm{bi})=\frac{1}{N}{e}^{\mathrm{jx}(1-n/N)}\frac{\sin (2{\mathrm{\πu}}_1(n-\mathrm{bi})/N)}{\sin (\mathrm{\π}(n-\mathrm{bi})/N)}$ $r_{u_1-u_2}(n-\mathrm{bi})$ 的幅值是一个振荡衰减序列。 $r_{u_1-u_2}(n-\mathrm{bi})$ 的幅值必然在bi,i=1,2…,Q各点分别取得局部最大值。据此，可检测出信号的各奇异点bi,i=1,2…,Q。为此，我们提出找到最大奇异强度的奇异点，便从含截断伪影的差分信号 $d_{u_1-u_2}\left(j\right)$ 中，从大到小的顺序检测奇异点。

层析法检测检测奇异点算法如下：

第一步，初始化。置奇异点集B为空。奇异点集元素个数计数器Q归零。

第二步，找出 $a=\max \left\{\right|d_{u_1-u_2}\left(0\right)|,|d_{u_1-u_2}\left(1\right)|,...,|d_{u_1-u_2}(N-1)\left|\right\}$ 的最大值点b_i；如果 $b_i\∉B$ ，则b_i加入到

奇异点集B中，并且计数器Q加1。

第三步，估算b_i的奇异强度 $a_{b_i}=\frac{N}{u_2-u_1}a$

第四步，从 $d_{u_1-u_2}\left(j\right)$ 减法正弦振荡衰减信号 $a_{b_i}{r}_{u_1-u_2}(j-b_i)$ ，即 $d_{u_1-u_2}\left(j\right)=d_{u_1-u_2}\left(j\right)-a_{b_i}{r}_{u_1-u_2}(j-b_i),j=\mathrm{0,1},...,N-1$

第五步，如果α＞T(预先给定的阀值，一般取初始 $d_{u_1-u_2}\left(j\right)$ 最大值的1.0％)，则转回到第二步。

第六步，输出奇异点集B及奇异点个数计数器值Q。

第七步，结束。

(2)权系数确定：

设已测得的信号f(n),i=0,1,…,N-1的奇异点为bi，i=1,2…,Q。根据式(5)，可由 $f_{u_1-u_2}\left(n\right)$ 建立方程组：

记y=[y_b1,y_b2,…,y_bQ]^T为Q维加权系数构成的矢量， $f_{u_1-u_2}=\[f_{u_1-u_2}\left(0\right),f_{u_1-u_2}\left(1\right),...,f_{u_1-u_2}(N-1){\]}^T$ 式(8)的矩阵形式为： $W_{u_1-u_2}\·y=f_{u_1-u_2}\·\·\·\left(9\right)$ 其中

是由奇异点bi、频率u₁和u₂确定的、已知的矩阵。方程组(8)的解为： $y=W_{u_1-u_2}^{+}{f}_{u_1-u_2}\·\·\·\·\left(10\right)$ 其中，

的伪逆矩阵。

(3)X-CT有限角截断投影数据奇异谱分析图象重建算法

第一步，根据投影定理，把原始有限角截断投影数据转变为有限角截断频谱数据，如附图1(b)所示。

第二步，按列计算截断频率(参见附图2：有限角截断频谱数据截断频率计算图)：u₁=N/2^*tgφ,u₂=-N/2^*tgφ

第三步，对第u₁～N/2和u₂～-N/2各行进行付里叶反变换。

第四步，从0到N-1列分别执行：

1)运用层析法检测奇异点。

2)频建立诸奇异点对应的奇异谱函数的线性泛函，并由原始截断频谱数据确定其中的加权系数。

第五步，通过奇异函数分析公式(6)计算无截断伪影信号：

第六步，把各列信号综合成X-CT图象。

本方法的关键点是建立截断奇异函数分析模型，难点是奇异点的检测。但是奇异点的检测有相当的鲁棒性，即如检测到的奇异点集B中，含有非奇异点bi，则必有y_bi=0。W_bi(n)在f(n)中贡献为零。但是过多的非奇异包含在奇异点集B中，必将加重求 $W_{u_1-u_2}^{+}$ 的计算量和造成方程 $W_{u_1-u_2}\·y=f_{u_1-u_2}$ 欠定(欠定方程的伪逆解，是权向量模｜y｜最小的解)。我们一般选取奇异点B集的元素数量Q不大于被截取的频率数u₂-u₁，即Q≤u₂-u₁。当信号的奇异点数确实大于u₁-u₂时(如含噪声的情况)，我仅取局部极较大的前u₂-u₁个奇异点。一般地说(如果各奇异点的彼此间隔大于5)，局部极值较大的奇异点对应着较大的加权系数，局部极值较小的奇异点对应着较小的加权系数。因而，按 $f\left(n\right)=\underset{i=1}{\overset{Q}{\mathrm{\Σ}}}{y}_{\mathrm{bi}}{W}_{\mathrm{bi}}\left(n\right)$ ，局部极值较大奇异点对信号的贡献大，忽略局部极值较小的奇异点重建信号的误差相对较小。参照附图说明本发明的效果

X-CT有限角截断投影数据奇异谱分析图象重建算法的测试是在一台PIII的微机上进行的。为了充分认识奇异谱分析成象方法的效果，我们用对有噪声的、无噪声的，实际的和模拟的X-CT有限角截断投影数据进行了成象算法测试，并把本发明方法重建的图象和滤波反投影重建的图象进行如下比较。

(1)计算机模拟有限角(27°-153°)投影数据成象的图象比较

图3表示无噪声计算机模拟X-CT有限角投影数据图象重建情况。其中图(a)是128×128个象素的无噪声样本图象，图(b)和(c)是图(a)的有限角(27°-153°)投影数据分别用滤波反投影和奇异谱分析方法重建的图象。把图(b)、(c)分别和图(a)比较得，图(b)比图(c)有多得多的伪影，这说明奇异谱成象法能完全消除截断伪影，由表1知奇异谱成象法的保真指标比滤波反投影高出三个数量级，其中保真指标的计算公式为：

归一化均方误差，即 $\mathrm{\σ}=\frac{\underset{I=1}{\overset{M}{\mathrm{\Σ}}}\underset{J=1}{\overset{N}{\mathrm{\Σ}}}(x_{I,J}^o-x_{I,J}{)}^2}{\underset{I=1}{\overset{M}{\mathrm{\Σ}}}\underset{J=1}{\overset{N}{\mathrm{\Σ}}}{(x_{I,J}^o-\mathrm{\μ})}^2}\·\·\·\·\left(11\right)$ 归一化绝对误差，即 $E=\frac{\underset{I=1}{\overset{M}{\mathrm{\Σ}}}\underset{J=1}{\overset{N}{\mathrm{\Σ}}}|x_{I,J}^o-x_{I,J}|}{\underset{I=1}{\overset{M}{\mathrm{\Σ}}}\underset{J=1}{\overset{N}{\mathrm{\Σ}}}|x_{I,J}^o-\mathrm{\μ}|}\·\·\·\·\left(12\right)$ 其中 $X_{I,J}^o,X_{I,J}$ 分别是样本图象和重建图象的第I行第J列的图象素灰度值。μ是样本图象的灰度均值。M是图象的行数，N是图象的列数。

图4表示含5％噪声计算机模拟X-CT有限角投影数据图象重建情况。其中图(a)是128×128个象素的含5％噪声样本图象。图(b)和(c)是图(a)的有限角(27°-153°)投影数据分别用滤波反投影和奇异谱分析方法重建的图象。把图(b)、(c)分别和图(a)比较得，图(b)比图(c)有多得多的伪影。由表2知奇异谱成象法的保真指标比滤波反投影高出二个数量级。这说明奇异谱成象法在有噪声情况下也消除截断伪影，但由表1、表2比较可见算法对噪声敏感，这是因为由于噪声的引入使得奇异点数量Q急剧上升，频谱数据量n小于奇异点数量Q或方程(9)成为不相容，使方程(9)只能得一伪逆解。

(2)实际X-CT投影数据重建算法的比较

图5表示为实际X-CT有限角投影数据图象重建情况。其中图(a)实际X-CT图象，图(b)(c)为图(a)的有限角(27°-153°)投影数据分别用滤波反投影和奇异谱分析方法重建的图象。实际X-CT有限角(27°-153°)投影数据成象情况，奇异谱成象法也能很好地消除截断伪影。

5、结论

在奇异点个数较少的情况下，奇异谱分析法重建图象的误差仅由计算工具精度决定，即理论上讲缺损频率分量可以准确恢复。但对于含噪声情况下，由于在预测奇异点时采用层析法，得到重建图象的优化解。使得奇异谱分析成象，不论含噪声的还是无噪声的、对实际的还是对计算机模拟的X-CT有限角投影数据的实验结果，都表明了奇异谱分析成象是一种适合于有限角投影数据图象重建的高精度重建方法，保证能消除重建的图象的截断伪影，图象的质量大大优于用传统方法重建的图象。

表1．无噪声时的算法重建误差、重建时间比较表

算法	归一化均方误差	归一化绝对误差
			滤波反投影	0.2308772	0.5937076
奇异谱分析成象法	0.0002016	0.0006317

表2．含5％高斯噪声时的算法重建误差、重建时间比较表

算法	归一化均方误差	归一化绝对误差
			滤波反投影	0.3033662	0.5854027
奇异谱分析成象法	0.0040634	0.0058211

附图说明：

图1为投影数据示意图，其中(a)完整投影数据频谱图，(b)有限角投影数据频谱图

图2为有限角截断频谱数据截断频率计算图。

图3无噪声计算机模拟X-CT有限角投影数据图象例比较，其中(a)是128×128个象素的无噪声样本图象，(b)为有限角(27°-153°)投影数据分别用滤波反投影法重建图象，(c)是本方法重建的图象。

图4有噪声计算机模拟X-CT有限角投影数据图象例比较，其中(a)是128×128个象素的无噪声样本图象，(b)为有限角(27°-153°)投影数据分别用滤波反投影法重建图象，(c)是本方法重建的图象。

图5实际X-CT有限角投影数据图象例比较，其中(a)是256×256个象素的样本图象，(b)为有限角(27°-153°)投影数据分别用滤波反投影法重建图象，(c)是本方法重建的图象。

本发明受国家自然科学基金资助，批准号码：39870211和39970219

Claims (1)

1、这种X-CT有限角投影数据的图象重建方法，其特征在于采用下列操作步骤：

(1)根据投影定理，把原始有限角截断投影数据转变为有限角截断频谱数据；

(2)按列计算截断频率：u₁=N/2^*tgφ,u₂=-N/2^*tgφ；

(3)对第u₁～N/2和u₂～-N/2各行进行付里叶反变换；

(4)从0到N-1列分别执行：

1)运用层析法检测奇异点。即执行下列①初始化，置奇异点集B为空。奇异点集元素个数计数器Q归零。②找出 $a=\max \left\{\right|d_{u_1-u_2}\left(0\right)|,|d_{u_1-u_2}\left(1\right)|,...,|d_{u_1-u_2}(N-1)\left|\right\}$ 的最大值点b_i，如果 $b_i\∉B$ ，则b_i加入到奇异点集B中，并且计数器Q加1。③估算b_i的奇异强度 $a_{b_i}=\frac{N}{u_2-u_1}a;$ ④从 $d_{u_1-u_2}\left(j\right)$ 减法正弦振荡衰减信号 $a_{b_i}{r}_{u_1-u_2}(j-b_i)$ ，即 $d_{u_1-u_2}\left(j\right)=d_{u_1-u_2}\left(j\right)-a_{b_i}{r}_{u_1-u_2}(j-b_i),j=0,1,...,N-1;$ ⑤如果a＞T(预先给定的阀值，一般取初始 $d_{u_1-u_2}\left(j\right)$ 最大值的1.0％)，则转回到第二步。⑥输出奇异点集B及奇异点个数计数器值Q。

2)频建立诸奇异点对应的奇异谱函数的线性泛函，并由原始截断频谱数据确定其中的加权系数；

(5)按下式计算无截断伪影信号： $f\left(n\right)=\underset{i=1}{\overset{Q}{\mathrm{\Σ}}}{a}_{\mathrm{bi}}{W}_{\mathrm{bi}}\left(n\right),$ $a_{\mathrm{bi}}\≠0,i=1,...,Q,n=0,1,...,N-1;$

(6)把各列信号综合成X-CT图象；

(7)结束。

Images