专利名称：X－ct有限角投影数据图象重建方法
技术领域：
本发明属于X-CT投影数据图象重建技术领域，涉及截断奇异函数分析的图象重建方法。
计算机析层成象(CT)已广泛地应用于医学成象、地震探测、材料非损伤检测等各种领域。CT成象效果常常依赖于采集CT投影数据条件与理想条件的差距。例如，投影数据和旋转角的采样。在工业应用中，常常遇见被CT检查对象的尺寸比CT扫描系统检测器的宽度大得多，只能用有限角投影数据重建图象。这种有限角投影数据用现有方法重建的图象总有伪影。旋转角只能分布在φ1～φ2之间，不能均匀分布于
]>是单位离散矩形函数，即R^u1-u2(u)={0,u∉{u1,u1+1,...,u2-1}1,u∈{u1,u1+1,...,u2-1}]]>其中，u1是截断频谱的起始频率，u2是截断频谱的终止频率。则截断频谱信号fu1-u2(n)]]>可表示为fu1-u2(n)=IDFT[f^(u)R^u1-u2(u)]]]>即fu1-u2(n)=f(n)⊗ru1-u2(n)····(1)]]>这里，IDFT[·]表示逆离散傅里叶变换，是圆周卷积符号，ru1-u2(n)=IDFT[R^u1-u2(u)],]]>即rui-u2(n)=1Nej2πu1n/N-ej2πu2/N1-ej2πn/N]]>离散信号的截断频谱信号重建问题可以表述为从fu1-u2(n)]]>中恢复出原来的信号f(n),n=0,1….,N-1。在中国发明专利申请号为98107580.0的专利申请文件中，我们定义了奇异空间{Wk(n),k,n为整数}，其中Wk(n)为Wk(n)=Σl=0nδ(l-k),n,k=0,1,...,N-1,]]>其中，δ(n)为离散单位脉冲函数。任何有限长离散信号f(n),n=0,1,…,N-1可表示为奇异空间{Wk(n),k,n∈0,1,…,N-1}中的线性泛函。f(n)=Σk=0N-1akWk(n)····(2)]]>对上式两边同时进行ru1-u2(n)]]>的圆周卷积得fu1-u2(n)=Σk=0N-1akWk(n)⊗ru1-u2(n)···(3)]]>我们定义Wk,u1-u2(n)=Wk(n)⊗ru1-u2(n)]]>为截断奇异函数。诸截断奇异函数{Wk,u1-u2(n)|k=0,1,...,N-1}]]>构成N维截断奇异函数空间。式(3)表明截断频谱信号fu1-u2(n)]]>可以表达为N维截断奇异函数空间中线性泛函数。即fu1-u2(n)=Σk=0N-1akWk,u1-u2(n)····(4)]]>若信号中仅有Q个奇异点(在ak,k=0,1,…,N-1中仅有Q个不为0)，则截断频谱信号fu1-u2(n)]]>可由{Wk,u1-u2(n)|k=0,1,...,N-1}]]>中对应的Q个截断奇异函数的线性泛函表示。即fu1-u2(n)=Σi=1QabiWbi,u1-u2(n),abi≠0,i=1,...,Q····(5)]]>其中Wk,u1-u2(n)]]>是以bi位置为奇异点的截断奇异函数，abi是Wk,u1-u2(n)]]>的加权系数。当u1=0,u2=N-1时，显而易见应有f0-(N-1)(n)=f(n)，这时式(5)就变为 所以，获得信号的奇异点及其奇异谱分析加权系数就等价于获取信号的频谱数据。从fu1-u2(n)]]>重建f(n)的关键问题是由fu1-u2(n)]]>求出奇异点bi和相应的加权系数abi,I=1,2…,Q。2、模型参数的求取根据投影定理，X-CT投影数据py′(x′)与g(x,Y)的频谱G(ω,φ)之间存在以下关系G(ω,φ)=∫-∞∞Py′(x′)e-jaxx′dx′---(7)]]>对于有限角范围投影数据集，我们可以由投影定理获得有限角范围内的CT截断频谱(简称为有限角截断频谱)，如附图1中的(b)图所示。即，X-CT有限角截断投影图象重建问题可以转变为有限角截断频谱的图象重建问题。
对于有限角截断频谱的图象重建问题，由于信号能量损失严重，小波系数误差太大，故不能用类似磁共振截断频谱图象重建方法的思路，用小波变换方法预选信号的奇异点。为此，我们根据奇异函数分析理论，运用层析法筛选信号的奇异点，确定奇异谱分析的加权系数。
(1)奇异点检测我们把fu1-u2(n)]]>看成周期为N的周期函数，则其差分函数du1-u2(n)]]>为 由于δ(n-bi)⊗uu1-u2(n)=ru1-u2(n-bi)]]>，即du1-u2(n)]]>立刻成为du1-u2(n)=Σi=1Qabiru1-u2(n-bi)····(8)]]>式(8)表明，du1-u2(n)]]>是Q个ru1-u2(n)]]>分别位于bi,i=1,2,…,Q,上的函数的加权和。当u2=N-u1，0＜u1＜N/2时,ru1-u2(n-bi)]]>可写为ru1-u2(n-bi)=1Nejx(1-n/N)sin(2πu1(n-bi)/N)sin(π(n-bi)/N)]]>ru1-u2(n-bi)]]>的幅值是一个振荡衰减序列。ru1-u2(n-bi)]]>的幅值必然在bi,i=1,2…,Q各点分别取得局部最大值。据此，可检测出信号的各奇异点bi,i=1,2…,Q。为此，我们提出找到最大奇异强度的奇异点，便从含截断伪影的差分信号du1-u2(j)]]>中，从大到小的顺序检测奇异点。
层析法检测检测奇异点算法如下第一步，初始化。置奇异点集B为空。奇异点集元素个数计数器Q归零。
第二步，找出a=max{|du1-u2(0)|,|du1-u2(1)|,...,|du1-u2(N-1)|}]]>的最大值点bi；如果bi∉B]]>，则bi加入到奇异点集B中，并且计数器Q加1。
第三步，估算bi的奇异强度abi=Nu2-u1a]]>第四步，从du1-u2(j)]]>减法正弦振荡衰减信号abiru1-u2(j-bi)]]>，即du1-u2(j)=du1-u2(j)-abiru1-u2(j-bi),j=0,1,...,N-1]]>第五步，如果α＞T(预先给定的阀值，一般取初始du1-u2(j)]]>最大值的1.0％)，则转回到第二步。
第六步，输出奇异点集B及奇异点个数计数器值Q。
第七步，结束。
(2)权系数确定设已测得的信号f(n),i=0,1,…,N-1的奇异点为bi，i=1,2…,Q。根据式(5)，可由fu1-u2(n)]]>建立方程组 记y=[yb1,yb2,…,ybQ]T为Q维加权系数构成的矢量，fu1-u2=[fu1-u2(0),fu1-u2(1),...,fu1-u2(N-1)]T]]>式(8)的矩阵形式为Wu1-u2·y=fu1-u2···(9)]]>其中 是由奇异点bi、频率u1和u2确定的、已知的矩阵。方程组(8)的解为y=Wu1-u2+fu1-u2····(10)]]>其中， 的伪逆矩阵。
(3)X-CT有限角截断投影数据奇异谱分析图象重建算法第一步，根据投影定理，把原始有限角截断投影数据转变为有限角截断频谱数据，如附图1(b)所示。
第二步，按列计算截断频率(参见附图2有限角截断频谱数据截断频率计算图)u1=N/2*tgφ,u2=-N/2*tgφ第三步，对第u1～N/2和u2～-N/2各行进行付里叶反变换。
第四步，从0到N-1列分别执行1)运用层析法检测奇异点。
2)频建立诸奇异点对应的奇异谱函数的线性泛函，并由原始截断频谱数据确定其中的加权系数。
第五步，通过奇异函数分析公式(6)计算无截断伪影信号 第六步，把各列信号综合成X-CT图象。
本方法的关键点是建立截断奇异函数分析模型，难点是奇异点的检测。但是奇异点的检测有相当的鲁棒性，即如检测到的奇异点集B中，含有非奇异点bi，则必有ybi=0。Wbi(n)在f(n)中贡献为零。但是过多的非奇异包含在奇异点集B中，必将加重求Wu1-u2+]]>的计算量和造成方程Wu1-u2·y=fu1-u2]]>欠定(欠定方程的伪逆解，是权向量模｜y｜最小的解)。我们一般选取奇异点B集的元素数量Q不大于被截取的频率数u2-u1，即Q≤u2-u1。当信号的奇异点数确实大于u1-u2时(如含噪声的情况)，我仅取局部极较大的前u2-u1个奇异点。一般地说(如果各奇异点的彼此间隔大于5)，局部极值较大的奇异点对应着较大的加权系数，局部极值较小的奇异点对应着较小的加权系数。因而，按f(n)=Σi=1QybiWbi(n)]]>，局部极值较大奇异点对信号的贡献大，忽略局部极值较小的奇异点重建信号的误差相对较小。参照
本发明的效果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知奇异谱成象法的保真指标比滤波反投影高出三个数量级，其中保真指标的计算公式为归一化均方误差，即σ=ΣI=1MΣJ=1N(xI,Jo-xI,J)2ΣI=1MΣJ=1N(xI,Jo-μ)2····(11)]]>归一化绝对误差，即E=ΣI=1MΣJ=1N|xI,Jo-xI,J|ΣI=1MΣJ=1N|xI,Jo-μ|····(12)]]>其中XI,Jo,XI,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．无噪声时的算法重建误差、重建时间比较表
表2．含5％高斯噪声时的算法重建误差、重建时间比较表
图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
权利要求
1.这种X-CT有限角投影数据的图象重建方法，其特征在于采用下列操作步骤(1)根据投影定理，把原始有限角截断投影数据转变为有限角截断频谱数据；(2)按列计算截断频率u1=N/2*tgφ,u2=-N/2*tgφ；(3)对第u1～N/2和u2～-N/2各行进行付里叶反变换；(4)从0到N-1列分别执行1)运用层析法检测奇异点。即执行下列①初始化，置奇异点集B为空。奇异点集元素个数计数器Q归零。②找出a=max{|du1-u2(0)|,|du1-u2(1)|,...,|du1-u2(N-1)|}]]>的最大值点bi，如果bi∉B]]>，则bi加入到奇异点集B中，并且计数器Q加1。③估算bi的奇异强度abi=Nu2-u1a;]]>④从du1-u2(j)]]>减法正弦振荡衰减信号abiru1-u2(j-bi)]]>，即du1-u2(j)=du1-u2(j)-abiru1-u2(j-bi),j=0,1,...,N-1;]]>⑤如果a＞T(预先给定的阀值，一般取初始du1-u2(j)]]>最大值的1.0％)，则转回到第二步。⑥输出奇异点集B及奇异点个数计数器值Q。2)频建立诸奇异点对应的奇异谱函数的线性泛函，并由原始截断频谱数据确定其中的加权系数；(5)按下式计算无截断伪影信号f(n)=Σi=1QabiWbi(n),]]>abi≠0,i=1,...,Q,n=0,1,...,N-1;]]>(6)把各列信号综合成X-CT图象；(7)结束。
全文摘要
X－CT截断投影数据图象重建方法是把信号表示成奇异函数的加权和。从采集到的X－CT有限角截断投影数据中提取用于重建图象的特征信息,并用这些特征信息构造截断奇异函数,再由这些截断奇异函数的线性泛函来重构X－CT图象,达到消除X－CT有限角投影数据图象重建中的截断伪影之目的。不论含噪声的还是无噪声的、对实际的还是对计算机模拟的X－CT有限角投影数据的实验结果,都表明这种方法效果高于现有方法,效果显著。
文档编号G01N23/04GK1300938SQ9912687
公开日2001年6月27日 申请日期1999年12月22日 优先权日1999年12月22日
发明者骆建华 申请人:上海交通大学