CN112200882B - 基于傅里叶变换微分性质的板状物ct图像快速重建方法 - Google Patents

Abstract

本发明公开了一种基于傅里叶变换微分性质的板状物CT图像快速重建方法，该方法包括：步骤1：利用CT系统扫描被测板状物，获得被测板状物的投影数据；步骤2：为板状物CT图像赋初值，作为板状物CT图像的估计值；步骤3：利用凸集投影算法使板状物CT图像的估计值满足优化模型的约束条件；步骤4：利用傅里叶变换微分性质等算法求解优化模型中的目标函数；步骤5：重复步骤3和步骤4，直到满足终止条件。本发明适用于X射线板状物CT图像重建，不仅可以有效去除板状物CT图像重建的伪影、抑制噪声，并且大大加快板状物CT图像的收敛速度。

Description

基于傅里叶变换微分性质的板状物CT图像快速重建方法

技术领域

本发明涉及X射线板状物CT图像重建技术领域，特别是关于一种基于傅里叶变换微分性质的板状物CT图像快速重建方法。

背景技术

X射线计算机断层成像技术(X ray Computed Tomography,简称X射线CT)可以在不破坏或损伤物体的情况下呈现物体的内部细节，已广泛应用于医学、生物、工业、材料、古化石、航天等众多领域。在工业应用中，常常需要对厚度较薄，而长、宽尺寸比厚度大的多的板状物体进行计算机断层成像(CT)检测，这些板状物如各类芯片、印刷电路板(PCB)、板状化石等。应用CT技术检测板状物时，为了重建出高分辨率的图像，一方面被测物体应尽可能靠近射线源，从而使被测物体在探测器上获得足够的放大比，此时被测物体在旋转过程中可能碰撞射线源，探测器仅能采集到有限角度的投影数据；另一方面，通常采用微焦点或纳米焦点射线源扫描板状物，这类射线源发射功率有限，X射线可能穿不透板状物的长边，或者仅有少部分光子穿透板状物的长边，这时由于探测器光子计数率低导致投影数据中的噪声水平相对提高。传统CT重建算法(ART,SART等)用有限角度的投影数据重建出的图像存在严重的有限角伪影。

在投影数据不理想的情况下，CT重建问题在数学上对应于病态反问题的求解问题，此时通常将其转化为包含数据项和正则化项的最优化模型进行求解。近年来，基于压缩感知理论，有多种正则化方法可用于图像重建，如离散梯度变换、非局部滤波、小波紧框架、低秩约束、先验图像约束的压缩感知(PICCS)、字典学习、基于深度学习的图像正则化等。其中非局部滤波会引起图像中与扫描方向相关的伪影扩散，而PICCS、字典学习、深度学习等需要先验图像或训练样本。而基于离散梯度域稀疏的先验约束简单有效，且不需要先验图像，被广泛应用于稀疏角、有限角、低剂量CT、圆轨迹锥束CT等重建问题。研究者提出了各项同性的TV算法(TV)、各向异性的TV算法(ATV)、再加权的各项异性TV算法、保边界扩散和光滑的梯度l₀算法等一系列有限角CT重建算法，均在一定程度上改善了有限角CT重建图像的质量。但在板状物CT重建问题中，这些算法不能快速得到高质量的板状物CT图像，限制了该类算法在实际中的应用。

发明内容

本发明的目的在于提供一种基于傅里叶变换微分性质的板状物CT图像快速重建方法，可以由采集到的板状物CT的投影数据快速重建出板状物的图像。

为实现上述目的，本发明提供一种基于傅里叶变换微分性质的板状物CT图像快速重建方法，该方法包括：

步骤1：利用CT系统扫描被测板状物，获得被测板状物的投影数据；

Af＝p (1)

式(1)中，f是离散化的二维数字图像，其像素总个数为N＝K×K，K表示二维数字图像的行数或列数，其元素f_n表示被测物体在第n个像素处的线性衰减系数；p为探测器在不同投影角下测量的投影数据，M表示CT扫描的总的射线数；A为正投影矩阵，其元素a_m,n表示第m条射线与第n个像素的交线长；

步骤2：为板状物CT图像赋初值，作为板状物CT图像的估计值；

步骤3：利用凸集投影算法使板状物CT图像的估计值满足优化模型的约束条件；其中，所述优化模型包括：

(一)目标函数，其表示为式(7)：

(二)约束条件，其包括：

数据保真项约束，其表示为式(8)：

非负约束，其表示为式(9)：

f≥0 (9)边界约束，其表示为式(10)：

上述各式中，f^*表示优化模型最终求出的最优解，D_x为沿x方向的梯度算子，

为满足约束条件的板状物CT图像的估计值，

为待求图像沿x方向的近似梯度

其中κ为阈值参数；R_s(f)为正则化项，用于去除求解过程中引入的伪结构；W为所述步骤1中投影数据的权重矩阵；

表示二范数的平方；ε为与投影数据噪声水平相关的数，实际应用时是优化方程的误差门限；Ω表示整幅二维图像；

步骤4：利用傅里叶变换微分性质等算法求解优化模型中的目标函数；

步骤5：重复步骤3和步骤4，直到满足终止条件。

进一步地，所述步骤3包括：

步骤3.1，利用凸集投影算法使板状物CT图像的估计值满足式(8)表示的数据保真项约束，其包括：

采用凸集投影方法，使重建结果满足约束条件的要求，假定f^k表示第k次迭代的结果，当前图像f^k在

上的投影∏(f^k)表示为：

∏(f^k)＝f^k+τA^T(AA^T)^-1W(p-Af^k) (13)

其中，参数τ＞0；

步骤3.2：根据式(9)表示的非负约束，利用式(15)和(16)求解

令

其中，f^k _i,j表示当前图像f^k中的第i行，第j列的像素。

进一步地，采用如式(14)所示的逐射线投影方法求解式(13)：

符号:＝为赋值运算，即将该符号右边的结果赋给左边的变量。

进一步地，

其中，

表示第p个探测器在第q个扫描角度下对应的投影数据的标准差，N_0p表示在不放物体的情况下第p个探测器探测到的光子数，

是探测器上电子噪声的方差，

表示第p个探测器在第q个扫描角度下下对应的投影数据的值。

进一步地，所述步骤4中，求解优化模型的目标函数时，将式(6)表示的优化问题转化为如下两个子问题，以更新f^k+1：

子问题1：

子问题2：

式中，λ₁表示用于平衡

之间的关系的参数，λ₂表示用于平衡R_s(f)和

之间的关系的参数。

进一步地，子问题1中，求解式(17)表示的优化问题，新的迭代点为：

根据傅里叶变换微分性质，式(19)转化为：

的计算公式如式(21)：

其中，参数κ根据x方向的梯度图像可见边界处的灰度值确定；

令

由(22)式算出的g，经过(23)式修正后的结果作为新的g：

g:＝g-repmat(g(:,1),1,K) (23)

其中，g(:,1)表示图像g的第一列，repmat(g(:,1),1,K)表示将图像g的第一列复制K份，组成一个K×K大小的图像。

进一步地，子问题1求解过程中，当ω_x＝0时，令ω_x＝10^-10，并通过式(24)，使g的结果与

的结果灰度值量级相同：

令参数λ₁随迭代次数由小到大逐渐变化，

表示f^k在(i,j)处的像素值，g_i,j表示g在(i,j)处的像素值；

子问题1的解表示为式(25)：

进一步地，在板状物CT重建时，子问题2的解采用将式(11)或式(12)所示的R_s(f)的表达式代入(18)式后分别对应的求解算法，或者采用沿y方向的中值滤波；

其中，f_x,y表示在点(x,y)处f的值；#{}是计数算子，在式(11)中计算满足|f_x,y-f_x,y-1|≠0的像素点的数目，在式(11)中计算满足|f_x,y-f_x-1,y|+|f_x,y-f_x,y-1|≠0的像素点的数目。

进一步地，为板状物CT图像赋初值f⁰＝0。

由于采取以上技术方案，本发明适用于X射线板状物CT图像重建，不仅可以有效去除板状物CT图像重建的伪影、抑制噪声，并且大大加快板状物CT图像的收敛速度。

附图说明

图1为本发明提出的一种基于傅里叶变换微分性质的板状物CT图像快速重建方法的流程图。

图2a为本发明一个实施例用作测试模体的三维模拟PCB模体。

图2b为本发明一个实施例用作测试模体的二维模拟PCB模体。

图3为本发明一个实施例采集到的二维模拟PCB模体在120°扫描角度下得到的投影数据图像。

图4为本发明一个实施例采用SART方法重建的结果图像。

图5为本发明一个实施例采用本发明方法重建的结果图像。

具体实施方式

下面结合附图和实施例对本发明进行详细的描述。

参阅图1，为本发明提出的一种基于傅里叶变换微分性质的板状物CT图像快速重建方法的流程图。由图1可以看出，本发明方法包括以下步骤：

步骤1：利用CT系统扫描被测板状物，获得被测板状物的投影数据；

步骤2：为板状物CT图像赋初值，作为板状物CT图像的估计值；

步骤3：利用凸集投影算法使板状物CT图像的估计值满足优化模型的约束条件；

步骤4：利用傅里叶变换微分性质等算法求解优化模型的目标函数；

步骤5：重复步骤3和步骤4，直到满足终止条件。

在一个实施例中，步骤1所述的获得被测板状物的投影数据的过程在数学上可以用如下线性方程组表示：

Af＝p (1)

其中：

f＝(f₁,f₂,…,f_n,…,f_N)^T (2)

f表示一个二维函数f(x,y)，此处的f为f(x,y)的简写，其是离散化的二维数字图像，像素总个数为N＝K×K，K表示二维数字图像的行数或列数，其元素f_n表示被测物体在第n个像素处的线性衰减系数，f(x,y)＝f_(i-1)*K+j。通常，一幅图像的左上角为坐标原点即(x,y)＝(0,0)，沿图像水平方向从左到右为x轴正方向，沿图像竖直方向从上到下为y轴正方向，水平或竖直相邻像素间距为1。但也可以根据实际需求，调整原点位置以及x轴和y轴的位置及其正方向。下文中的采用(i,j)表示图像中像素点的位置，i表示横坐标，j表示纵坐标。

p＝(p₁,p₂,…,p_m,…,p_M)^T (3)

p为探测器在不同投影角下测量的投影数据，M表示CT扫描的总的射线数(比如X射线数)。

A为正投影矩阵，适用于多种CT扫描模式，包括平行束CT、扇束CT和锥束CT等，其元素a_m,n表示第m条射线与第n个像素的交线长。此时，CT重建问题可以描述为由已知的正投影矩阵A和投影数据p求解图像f。

在一个实施例中，步骤2中的为板状物CT图像赋初值，是以0作为初始图像，即图像中的每个像素均为0。

在一个实施例中，以一维函数f(x)为例，所述步骤4中，傅里叶变换的微分性质给出了原函数与导函数之间的等式关系：

其中，F(·)为傅里叶变换算子，

为f沿x方向的导数，

为f(x)的傅里叶变换，ω为频率。式(4)变形可得：

则有式(6)：

其中，F^-1(·)为傅里叶逆变换算子。假设本发明实施例待求CT图像为f(x)，那么沿x方向的梯度为

在扫描PCB等板状物时，被测物体材质差异较大，可以得到相对准确的

再根据公式(6)，即可快速求出图像f(x)的近似结果。由于CT图像的可见边界处的梯度值是近似准确，而不是完全准确，所以本发明方法采用迭代算法，根据步骤3求得的满足约束条件的板状物CT图像的估计值

逐步修正板状物CT重建图像f^k+1。

在一个实施例中，所述步骤3中，优化模型包括：

(一)目标函数，其表示为式(7)：

(二)约束条件，其包括：

数据保真项约束，其表示为式(8)：

非负约束，其表示为式(9)：

f≥0 (9)

边界约束，其表示为式(10)：

其中，f^*表示优化模型最终求出的最优解，D_x为沿x方向的梯度算子，D_xf相当于

为满足约束条件的板状物CT图像的估计值，R_s(f)为正则化项，用于去除求解过程中引入的伪结构。W为投影数据的权重矩阵。

表示二范数的平方；ε为与投影数据噪声水平相关的数，实际应用时是优化方程的误差门限比如0.0001；

为待求图像沿x方向的近似梯度

其中κ为阈值参数。

在投影数据获取时，由于模型离散化、噪声和散射等因素的影响，本发明实施例很难满足式(1)所示的等式条件，因此本发明方法将数据保真项约束定义为一个不等式约束。此外，如前所述，板状物CT扫描时，不同角度下的光子计数差别很大，相应地不同角度下的统计噪声幅度也差别很大。因此不同角度下的投影数据可靠性不同，需要对可靠性不同的投影数据加权使用。为此，数据保真项要求由重建图像生成的投影数据与实际投影数据的加权距离的平方和在给定的距离ε²内。优化模型的非负约束f≥0表示待重建图像的灰度值大于0。优化模型的边界约束

表示待求图像f(x)的边界为0，在不超视野扫描的情况下，符合CT图像的特点。

投影矩阵的权重矩阵W与投影数据的噪声模型相关，不同的噪声模型推导出的权重矩阵对板状物CT图像噪声的抑制效果不同。本发明方法数据保真项权重W选取时，可有多种选取方式，如泊松噪声模型、非平稳高斯噪声模型、考虑系统电子噪声的噪声模型等。

R_s(f)有多种选择方式，如沿y方向的梯度l₀范数(11)、各向同性的梯度l₀范数(12)、或者沿y方向的中值滤波等。

其中，f_x,y表示在点(x,y)处f的值；#{}是计数算子，在式(11)中计算满足|f_x,y-f_x,y-1|≠0的像素点的数目，在式(11)中计算满足|f_x,y-f_x-1,y|+|f_x,y-f_x,y-1|≠0的像素点的数目。

在一个实施例中，所述步骤4中，傅里叶变换微分性质等算法，包括傅里叶变换微分性质和将式(11)或式(12)所示的R_s(f)的表达式代入(18)式后分别对应的求解算法，或者采用沿y方向的中值滤波。

在一个实施例中，所述步骤3包括：

步骤3.1，利用凸集投影算法使板状物CT图像的估计值满足式(8)表示的数据保真项约束，其包括：

假定f^k表示第k次迭代的结果，采用凸集投影方法，使当前图像f^k满足数据保真项约束的要求，当前图像f^k在

上的投影Π(f^k)可表示为：

∏(f^k)＝f^k+τA^T(AA^T)^-1W(p-Af^k) (13)

其中，参数τ＞0,CT重建问题的系统矩阵A的维度很大，直接求解式(13)较为困难，本发明方法求解时，采用如式(14)所示的逐射线投影方法求解式(13)

其中，符号:＝为赋值运算，即将该符号右边的结果赋给左边的变量；A_i表示第i条射线与待重建图像各个像素的交线长，A_i是矩阵A的第i行；W_i,i表示矩阵W第i行、第i列的元素值，W_i,i是矩阵W中的一个元素；p_i表示第i条射线穿过被测物体后得到的投影值，p_i是p中的一个元素。

步骤3.2：根据式(9)表示的非负约束，利用式(15)和(16)求解

令

其中，f^k _i,j表示当前图像f^k中的第i行，第j列的像素，

表示步骤3.2的输出结果。

在一个实施例中，所述步骤4中，求解优化模型的目标函数时。更新f^k+1可以将式(7)表示的优化问题转化为如下两个子问题：

子问题1：

子问题2：

式中，f^k+1/2表示满足式(17)代表的优化方程的最优解，

代表求解该符号左边表达式值最小的时候的f的值，

代表D_xf与

的距离的二范数的平方，

代表f与

的距离的二范数的平方。f^k+1表示满足(18)式代表的优化方程的最优解，R_s(f)表示与f相关的正则化项，

代表f与f^k+1/2的距离的二范数的平方。λ₁表示用于平衡

和

之间的关系的参数，取值方法下文有叙述，参数λ₁越大，算法收敛越慢，为了加速优化问题(20)的收敛速度，本发明方法采用Fixed pointcontinuation(FPC)技术，令参数λ₁随迭代次数由小到大逐渐变化。λ₂表示用于平衡R_s(f)和

之间的关系的参数，其具体数值通常是根据实验现象来调节，选取一个合适的数值以达到想要的结果。

子问题1用来利用可见边界处的梯度值快速恢复板状物的结构，但在子问题1求解时，是对图像沿x方向逐行进行求解，此时有可能会打乱图像相邻行灰度值之间的连续性，从而引入沿x方向的条状伪影。子问题2用来消除子问题1求解过程中引入的条状伪影。上述子问题拆分方式及求解算法，是基于目标函数的实际意义，而不是根据严格的数学推导得来的。

在一个实施例中，子问题1的求解方法如下：

求解(17)式表示的优化问题，新的迭代点为：

根据傅里叶变换微分性质，式(19)可转化为：

其中，ω_x表示ω沿x方向的频率，在一维的情况下，ω_x与ω等价。

本文算法的关键是选取合适的参数κ，以求得相对准确的

计算公式如下，

其中，参数κ根据x方向的梯度图像可见边界处的灰度值确定。

令

g不能满足式(10)表示的边界处灰度值为0的约束，为了满足该约束条件，数值实现时，(22)式的计算结果g，经过(23)式修正后的结果作为新的g，即每行分别减去g每行左边界处的灰度值，作为(22)式等号右端项的结果：

g:＝g-repmat(g(:,1),1,K) (23)

其中，g(:,1)表示图像g的第一列。repmat()为向量广播运算符，repmat(g(:,1),1,K)表示将图像g的第一列复制K份，组成一个K×K大小的图像。

在计算式(20)时，当ω_x取0频率时，式子没有意义，为避免这一情况，本文算法实现时，当ω_x取0频率时，令ω_x＝10^-10。这会使g的结果与

的结果的灰度值近似为倍数关系，本发明实施例通过(24)式对g进行简单处理，使g的结果与

的结果灰度值量级相同，再进行后续计算g。

式中，

表示f^k在(i,j)处的像素值，g_i,j表示g在(i,j)处的像素值。

此时，子问题1的解可表示如下：

参数λ₁越大，算法收敛越慢，为了加速优化问题(20)的收敛速度，本发明方法采用Fixed point continuation(FPC)技术，令参数λ₁随迭代次数由小到大逐渐变化。

子问题2的求解方法如下：根据如上所述，正则化项R_s(f)的选取方式不同，对应不同的子问题2的求解算法。在板状物CT重建时，用沿y方向的中值滤波的结果作为子问题2的解，相较于其它几种正则化方式，更能有效地去除子问题1引入的条状伪影。

f^k+1＝medfilter1d(f^k+1/2) (26)其中，medfilter1d(f^k+1/2)表示对f^k+1/2做沿y方向的中值滤波。

以下通过一个具体实施例来说明本发明提出的一种基于傅里叶变换微分性质的板状物CT图像快速重建方法的具体实现过程。

图2a为本具体实施例中用作测试模体的三维模拟PCB模体，截取该数据的一层作为本实施例中用作测试模体的二维模拟PCB模体，如图2b所示。

实验扫描参数设置如下：射线源到转台中心的距离为600mm，射线源到探测器的距离为1000mm。线探测器由512个探测器单元组成，每个探测器单元的尺寸为0.25mm。本实施例中，令α为射线源和转台中心的连线与转台水平方向的夹角，本发明实施例将二维模拟PCB模体沿α∈[π/6,5π/6]进行扫描，投影数据用光线投射算法作用于模体数据上生成，如图3所示。

具体实施步骤如下：

步骤1：为板状物CT图像赋初值f⁰＝0；

步骤2：假定经过k次迭代后，有板状物CT图像的估计值f^k。利用公式(14)、(15)和(16)使板状物CT图像的估计值满足优化模型的约束条件。W有多种选取方式，如泊松噪声模型、非平稳高斯噪声模型、考虑系统电子噪声的噪声模型等。本实施例中，数据保真项权重W选取时，选择考虑系统电子噪声的噪声模型，即

其中，

表示第p个探测器在第q个扫描角度下对应的投影数据的标准差，N_0p表示在不放物体的情况下第p个探测器探测到的光子数，

是探测器上电子噪声的方差,

表示第p个探测器在第q个扫描角度下下对应的投影数据的值。其中计算所需的CT系统空扫时探测器的计数N_0p根据实际扫描条件获得，电子噪声方差较难获得精确值，本实施例中，本发明实施例取

步骤3：利用式(25)和沿y方向的中值滤波求解优化模型的目标函数；

步骤4：判断终止条件是否满足，若满足则终止迭代，否则转向步骤2。

图4为采用SART方法重建的结果图像。本发明方法进行5次迭代得到的板状物CT图像如图5所示，灰度窗设置为[0,0.15]。由图5结果可看出，本发明方法可以由采集到的板状物CT的投影数据快速重建出被测板状物的CT图像。

本领域普通技术人员可以理解：附图只是一个实施例的示意图，附图中的模块或流程并不一定是实施本发明所必须的。

最后需要指出的是：以上实施例仅用以说明本发明的技术方案，而非对其限制。本领域的普通技术人员应当理解：可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。

Claims (9)

1.一种基于傅里叶变换微分性质的板状物CT图像快速重建方法，其特征在于，包括：

步骤1：利用CT系统扫描被测板状物，获得被测板状物的投影数据；

Af＝p (1)

式(1)中，f是离散化的二维数字图像，其像素总个数为N＝K×K，K表示二维数字图像的行数或列数，其元素f_n表示被测物体在第n个像素处的线性衰减系数；p为探测器在不同投影角下测量的投影数据，M表示CT扫描的总的射线数；A为正投影矩阵，其元素a_m,n表示第m条射线与第n个像素的交线长；

步骤2：为板状物CT图像赋初值，作为板状物CT图像的估计值；

步骤3：利用凸集投影算法使板状物CT图像的估计值满足优化模型的约束条件；其中，所述优化模型包括：

(一)目标函数，其表示为式(7)：

(二)约束条件，其包括：

数据保真项约束，其表示为式(8)：

非负约束，其表示为式(9)：

f≥0 (9)

边界约束，其表示为式(10)：

上述各式中，f^*表示优化模型最终求出的最优解，D_x为沿x方向的梯度算子，

为满足约束条件的板状物CT图像的估计值，

为待求图像沿x方向的近似梯度

其中κ为阈值参数；R_s(f)为正则化项，用于去除求解过程中引入的伪结构；W为所述步骤1中投影数据的权重矩阵；

表示二范数的平方；ε为与投影数据噪声水平相关的数，实际应用时是优化方程的误差门限；Ω表示整幅二维图像；

步骤4：利用傅里叶变换微分性质等算法求解优化模型中的目标函数；

步骤5：重复步骤3和步骤4，直到满足终止条件。

2.如权利要求1所述的基于傅里叶变换微分性质的板状物CT图像快速重建方法，其特征在于，所述步骤3包括：

步骤3.1，利用凸集投影算法使板状物CT图像的估计值满足式(8)表示的数据保真项约束，其包括：

采用凸集投影方法，使重建结果满足约束条件的要求，假定f^k表示第k次迭代的结果，当前图像f^k在

上的投影Π(f^k)表示为：

Π(f^k)＝f^k+τA^T(AA^T)^-1W(p-Af^k) (13)

其中，参数τ>0；

步骤3.2：根据式(9)表示的非负约束，利用式(15)和(16)求解

令

其中，f^k _i,j表示当前图像f^k中的第i行，第j列的像素。

3.如权利要求2所述的基于傅里叶变换微分性质的板状物CT图像快速重建方法，其特征在于，采用如式(14)所示的逐射线投影方法求解式(13)：

符号:＝为赋值运算，即将该符号右边的结果赋给左边的变量。

4.如权利要求2所述的基于傅里叶变换微分性质的板状物CT图像快速重建方法，其特征在于，

其中，

表示第p个探测器在第q个扫描角度下对应的投影数据的标准差，N_0p表示在不放物体的情况下第p个探测器探测到的光子数，

是探测器上电子噪声的方差，

表示第p个探测器在第q个扫描角度下对应的投影数据的值。

5.如权利要求1至4中任一项所述的基于傅里叶变换微分性质的板状物CT图像快速重建方法，其特征在于，所述步骤4中，求解优化模型的目标函数时，将式(7)表示的优化问题转化为如下两个子问题，以更新f^k+1：

子问题1：

子问题2：

式中，λ₁表示用于平衡

和

之间的关系的参数，λ₂表示用于平衡R_s(f)和

之间的关系的参数。

6.如权利要求5所述的基于傅里叶变换微分性质的板状物CT图像快速重建方法，其特征在于，子问题1中，求解式(17)表示的优化问题，新的迭代点为：

根据傅里叶变换微分性质，式(19)转化为：

的计算公式如式(21)：

其中，参数κ根据x方向的梯度图像可见边界处的灰度值确定；

令

由(22)式算出的g，经过(23)式修正后的结果作为新的g：

g:＝g-repmat(g(:,1),1,K) (23)

其中，g(:,1)表示图像g的第一列，repmat(g(:,1),1,K)表示将图像g的第一列复制K份，组成一个K×K大小的图像，ω_x表示ω沿x方向的频率。

7.如权利要求6所述的基于傅里叶变换微分性质的板状物CT图像快速重建方法，其特征在于，子问题1求解过程中，当ω_x＝0时，令ω_x＝10^-10，并通过式(24)，使g的结果与

的结果灰度值量级相同：

令参数λ₁随迭代次数由小到大逐渐变化，

表示f^k在(i,j)处的像素值，g_i,j表示g在(i,j)处的像素值；

子问题1的解表示为式(25)：

8.如权利要求5所述的基于傅里叶变换微分性质的板状物CT图像快速重建方法，其特征在于，在板状物CT重建时，子问题2的解采用将式(11)或式(12)所示的R_s(f)的表达式代入(18)式后分别对应的求解算法，或者采用沿y方向的中值滤波；

其中，f_x,y表示在点(x,y)处f的值；#{}是计数算子，在式(11)中计算满足|f_x,y-f_x,y-1|≠0的像素点的数目，在式(12)中计算满足|f_x,y-f_x-1,y|+|f_x,y-f_x,y-1|≠0的像素点的数目。

9.如权利要求5所述的基于傅里叶变换微分性质的板状物CT图像快速重建方法，其特征在于，为板状物CT图像赋初值f⁰＝0。

Images