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

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

Info

Publication number
CN112200882A
CN112200882A CN202011077364.1A CN202011077364A CN112200882A CN 112200882 A CN112200882 A CN 112200882A CN 202011077364 A CN202011077364 A CN 202011077364A CN 112200882 A CN112200882 A CN 112200882A
Authority
CN
China
Prior art keywords
image
plate
formula
fourier transform
shaped object
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202011077364.1A
Other languages
English (en)
Other versions
CN112200882B (zh
Inventor
赵树森
夏迪梦
赵星
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Capital Normal University
Original Assignee
Capital Normal University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Capital Normal University filed Critical Capital Normal University
Priority to CN202011077364.1A priority Critical patent/CN112200882B/zh
Publication of CN112200882A publication Critical patent/CN112200882A/zh
Application granted granted Critical
Publication of CN112200882B publication Critical patent/CN112200882B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/00—2D [Two Dimensional] image generation
    • G06T11/003—Reconstruction from projections, e.g. tomography
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G06F17/12—Simultaneous equations, e.g. systems of linear equations
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/14—Fourier, Walsh or analogous domain transformations, e.g. Laplace, Hilbert, Karhunen-Loeve, transforms
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/60—Analysis of geometric attributes
    • G06T7/64—Analysis of geometric attributes of convexity or concavity
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/04—Constraint-based CAD
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10072—Tomographic images
    • G06T2207/10081—Computed x-ray tomography [CT]

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Data Mining & Analysis (AREA)
  • General Engineering & Computer Science (AREA)
  • Geometry (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Evolutionary Computation (AREA)
  • Computer Hardware Design (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Operations Research (AREA)
  • Analysing Materials By The Use Of Radiation (AREA)
  • Apparatus For Radiation Diagnosis (AREA)
  • Image Analysis (AREA)

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算法、保边界扩散和光滑的梯度l0算法等一系列有限角CT重建算法，均在一定程度上改善了有限角CT重建图像的质量。但在板状物CT重建问题中，这些算法不能快速得到高质量的板状物CT图像，限制了该类算法在实际中的应用。
发明内容
本发明的目的在于提供一种基于傅里叶变换微分性质的板状物CT图像快速重建方法，可以由采集到的板状物CT的投影数据快速重建出板状物的图像。
为实现上述目的，本发明提供一种基于傅里叶变换微分性质的板状物CT图像快速重建方法，该方法包括：
步骤1：利用CT系统扫描被测板状物，获得被测板状物的投影数据；
Af＝p (1)
式(1)中，f是离散化的二维数字图像，其像素总个数为N＝K×K，K表示二维数字图像的行数或列数，其元素fn表示被测物体在第n个像素处的线性衰减系数；p为探测器在不同投影角下测量的投影数据，M表示CT扫描的总的射线数；A为正投影矩阵，其元素am,n表示第m条射线与第n个像素的交线长；
步骤2：为板状物CT图像赋初值，作为板状物CT图像的估计值；
步骤3：利用凸集投影算法使板状物CT图像的估计值满足优化模型的约束条件；其中，所述优化模型包括：
(一)目标函数，其表示为式(7)：
Figure BDA0002717721590000021
(二)约束条件，其包括：
数据保真项约束，其表示为式(8)：
Figure BDA0002717721590000022
非负约束，其表示为式(9)：
f≥0 (9)边界约束，其表示为式(10)：
Figure BDA0002717721590000025
上述各式中，f*表示优化模型最终求出的最优解，Dx为沿x方向的梯度算子，
Figure BDA0002717721590000023
为满足约束条件的板状物CT图像的估计值，
Figure BDA0002717721590000024
为待求图像沿x方向的近似梯度
Figure BDA0002717721590000031
其中κ为阈值参数；Rs(f)为正则化项，用于去除求解过程中引入的伪结构；W为所述步骤1中投影数据的权重矩阵；
Figure BDA0002717721590000032
表示二范数的平方；ε为与投影数据噪声水平相关的数，实际应用时是优化方程的误差门限；Ω表示整幅二维图像；
步骤4：利用傅里叶变换微分性质等算法求解优化模型中的目标函数；
步骤5：重复步骤3和步骤4，直到满足终止条件。
进一步地，所述步骤3包括：
步骤3.1，利用凸集投影算法使板状物CT图像的估计值满足式(8)表示的数据保真项约束，其包括：
采用凸集投影方法，使重建结果满足约束条件的要求，假定fk表示第k次迭代的结果，当前图像fk在
Figure BDA0002717721590000033
上的投影∏(fk)表示为：
∏(fk)＝fk+τAT(AAT)-1W(p-Afk) (13)
其中，参数τ＞0；
步骤3.2：根据式(9)表示的非负约束，利用式(15)和(16)求解
Figure BDA0002717721590000034
Figure BDA0002717721590000035
令
Figure BDA0002717721590000036
其中，fk i,j表示当前图像fk中的第i行，第j列的像素。
进一步地，采用如式(14)所示的逐射线投影方法求解式(13)：
Figure BDA0002717721590000037
符号:＝为赋值运算，即将该符号右边的结果赋给左边的变量。
进一步地，
Figure BDA0002717721590000038
Figure BDA0002717721590000041
其中，
Figure BDA0002717721590000042
表示第p个探测器在第q个扫描角度下对应的投影数据的标准差，N0p表示在不放物体的情况下第p个探测器探测到的光子数，
Figure BDA0002717721590000043
是探测器上电子噪声的方差，
Figure BDA0002717721590000044
表示第p个探测器在第q个扫描角度下下对应的投影数据的值。
进一步地，所述步骤4中，求解优化模型的目标函数时，将式(6)表示的优化问题转化为如下两个子问题，以更新fk+1：
子问题1：
Figure BDA0002717721590000045
子问题2：
Figure BDA0002717721590000046
式中，λ1表示用于平衡
Figure BDA0002717721590000047
之间的关系的参数，λ2表示用于平衡Rs(f)和
Figure BDA0002717721590000048
之间的关系的参数。
进一步地，子问题1中，求解式(17)表示的优化问题，新的迭代点为：
Figure BDA0002717721590000049
根据傅里叶变换微分性质，式(19)转化为：
Figure BDA00027177215900000410
Figure BDA00027177215900000411
的计算公式如式(21)：
Figure BDA00027177215900000412
其中，参数κ根据x方向的梯度图像可见边界处的灰度值确定；
令
Figure BDA0002717721590000051
由(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的结果与
Figure BDA0002717721590000052
的结果灰度值量级相同：
Figure BDA0002717721590000053
令参数λ1随迭代次数由小到大逐渐变化，
Figure BDA0002717721590000054
表示fk在(i,j)处的像素值，gi,j表示g在(i,j)处的像素值；
子问题1的解表示为式(25)：
Figure BDA0002717721590000055
进一步地，在板状物CT重建时，子问题2的解采用将式(11)或式(12)所示的Rs(f)的表达式代入(18)式后分别对应的求解算法，或者采用沿y方向的中值滤波；
Figure BDA0002717721590000056
Figure BDA0002717721590000061
其中，fx,y表示在点(x,y)处f的值；#{}是计数算子，在式(11)中计算满足|fx,y-fx,y-1|≠0的像素点的数目，在式(11)中计算满足|fx,y-fx-1,y|+|fx,y-fx,y-1|≠0的像素点的数目。
进一步地，为板状物CT图像赋初值f0＝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＝(f1,f2,…,fn,…,fN)T (2)
f表示一个二维函数f(x,y)，此处的f为f(x,y)的简写，其是离散化的二维数字图像，像素总个数为N＝K×K，K表示二维数字图像的行数或列数，其元素fn表示被测物体在第n个像素处的线性衰减系数，f(x,y)＝f(i-1)*K+j。通常，一幅图像的左上角为坐标原点即(x,y)＝(0,0)，沿图像水平方向从左到右为x轴正方向，沿图像竖直方向从上到下为y轴正方向，水平或竖直相邻像素间距为1。但也可以根据实际需求，调整原点位置以及x轴和y轴的位置及其正方向。下文中的采用(i,j)表示图像中像素点的位置，i表示横坐标，j表示纵坐标。
p＝(p1,p2,…,pm,…,pM)T (3)
p为探测器在不同投影角下测量的投影数据，M表示CT扫描的总的射线数(比如X射线数)。
A为正投影矩阵，适用于多种CT扫描模式，包括平行束CT、扇束CT和锥束CT等，其元素am,n表示第m条射线与第n个像素的交线长。此时，CT重建问题可以描述为由已知的正投影矩阵A和投影数据p求解图像f。
在一个实施例中，步骤2中的为板状物CT图像赋初值，是以0作为初始图像，即图像中的每个像素均为0。
在一个实施例中，以一维函数f(x)为例，所述步骤4中，傅里叶变换的微分性质给出了原函数与导函数之间的等式关系：
Figure BDA0002717721590000071
其中，F(·)为傅里叶变换算子，
Figure BDA0002717721590000072
为f沿x方向的导数，
Figure BDA0002717721590000073
为f(x)的傅里叶变换，ω为频率。式(4)变形可得：
Figure BDA0002717721590000074
则有式(6)：
Figure BDA0002717721590000075
其中，F-1(·)为傅里叶逆变换算子。假设本发明实施例待求CT图像为f(x)，那么沿x方向的梯度为
Figure BDA0002717721590000081
在扫描PCB等板状物时，被测物体材质差异较大，可以得到相对准确的
Figure BDA0002717721590000082
再根据公式(6)，即可快速求出图像f(x)的近似结果。由于CT图像的可见边界处的梯度值是近似准确，而不是完全准确，所以本发明方法采用迭代算法，根据步骤3求得的满足约束条件的板状物CT图像的估计值
Figure BDA00027177215900000810
逐步修正板状物CT重建图像fk+1。
在一个实施例中，所述步骤3中，优化模型包括：
(一)目标函数，其表示为式(7)：
Figure BDA0002717721590000083
(二)约束条件，其包括：
数据保真项约束，其表示为式(8)：
Figure BDA0002717721590000084
非负约束，其表示为式(9)：
f≥0 (9)
边界约束，其表示为式(10)：
Figure BDA00027177215900000811
其中，f*表示优化模型最终求出的最优解，Dx为沿x方向的梯度算子，Dxf相当于
Figure BDA0002717721590000085
Figure BDA0002717721590000086
为满足约束条件的板状物CT图像的估计值，Rs(f)为正则化项，用于去除求解过程中引入的伪结构。W为投影数据的权重矩阵。
Figure BDA0002717721590000087
表示二范数的平方；ε为与投影数据噪声水平相关的数，实际应用时是优化方程的误差门限比如0.0001；
Figure BDA0002717721590000088
为待求图像沿x方向的近似梯度
Figure BDA0002717721590000089
其中κ为阈值参数。
在投影数据获取时，由于模型离散化、噪声和散射等因素的影响，本发明实施例很难满足式(1)所示的等式条件，因此本发明方法将数据保真项约束定义为一个不等式约束。此外，如前所述，板状物CT扫描时，不同角度下的光子计数差别很大，相应地不同角度下的统计噪声幅度也差别很大。因此不同角度下的投影数据可靠性不同，需要对可靠性不同的投影数据加权使用。为此，数据保真项要求由重建图像生成的投影数据与实际投影数据的加权距离的平方和在给定的距离ε2内。优化模型的非负约束f≥0表示待重建图像的灰度值大于0。优化模型的边界约束
Figure BDA0002717721590000091
表示待求图像f(x)的边界为0，在不超视野扫描的情况下，符合CT图像的特点。
投影矩阵的权重矩阵W与投影数据的噪声模型相关，不同的噪声模型推导出的权重矩阵对板状物CT图像噪声的抑制效果不同。本发明方法数据保真项权重W选取时，可有多种选取方式，如泊松噪声模型、非平稳高斯噪声模型、考虑系统电子噪声的噪声模型等。
Rs(f)有多种选择方式，如沿y方向的梯度l0范数(11)、各向同性的梯度l0范数(12)、或者沿y方向的中值滤波等。
Figure BDA0002717721590000092
Figure BDA0002717721590000093
其中，fx,y表示在点(x,y)处f的值；#{}是计数算子，在式(11)中计算满足|fx,y-fx,y-1|≠0的像素点的数目，在式(11)中计算满足|fx,y-fx-1,y|+|fx,y-fx,y-1|≠0的像素点的数目。
在一个实施例中，所述步骤4中，傅里叶变换微分性质等算法，包括傅里叶变换微分性质和将式(11)或式(12)所示的Rs(f)的表达式代入(18)式后分别对应的求解算法，或者采用沿y方向的中值滤波。
在一个实施例中，所述步骤3包括：
步骤3.1，利用凸集投影算法使板状物CT图像的估计值满足式(8)表示的数据保真项约束，其包括：
假定fk表示第k次迭代的结果，采用凸集投影方法，使当前图像fk满足数据保真项约束的要求，当前图像fk在
Figure BDA0002717721590000094
上的投影Π(fk)可表示为：
∏(fk)＝fk+τAT(AAT)-1W(p-Afk) (13)
其中，参数τ＞0,CT重建问题的系统矩阵A的维度很大，直接求解式(13)较为困难，本发明方法求解时，采用如式(14)所示的逐射线投影方法求解式(13)
Figure BDA0002717721590000101
其中，符号:＝为赋值运算，即将该符号右边的结果赋给左边的变量；Ai表示第i条射线与待重建图像各个像素的交线长，Ai是矩阵A的第i行；Wi,i表示矩阵W第i行、第i列的元素值，Wi,i是矩阵W中的一个元素；pi表示第i条射线穿过被测物体后得到的投影值，pi是p中的一个元素。
步骤3.2：根据式(9)表示的非负约束，利用式(15)和(16)求解
Figure BDA0002717721590000102
Figure BDA0002717721590000103
令
Figure BDA0002717721590000104
其中，fk i,j表示当前图像fk中的第i行，第j列的像素，
Figure BDA0002717721590000105
表示步骤3.2的输出结果。
在一个实施例中，所述步骤4中，求解优化模型的目标函数时。更新fk+1可以将式(7)表示的优化问题转化为如下两个子问题：
子问题1：
Figure BDA0002717721590000106
子问题2：
Figure BDA0002717721590000107
式中，fk+1/2表示满足式(17)代表的优化方程的最优解，
Figure BDA0002717721590000108
代表求解该符号左边表达式值最小的时候的f的值，
Figure BDA0002717721590000109
代表Dxf与
Figure BDA00027177215900001010
的距离的二范数的平方，
Figure BDA0002717721590000111
代表f与
Figure BDA0002717721590000112
的距离的二范数的平方。fk+1表示满足(18)式代表的优化方程的最优解，Rs(f)表示与f相关的正则化项，
Figure BDA0002717721590000113
代表f与fk+1/2的距离的二范数的平方。λ1表示用于平衡
Figure BDA0002717721590000114
和
Figure BDA0002717721590000115
之间的关系的参数，取值方法下文有叙述，参数λ1越大，算法收敛越慢，为了加速优化问题(20)的收敛速度，本发明方法采用Fixed pointcontinuation(FPC)技术，令参数λ1随迭代次数由小到大逐渐变化。λ2表示用于平衡Rs(f)和
Figure BDA0002717721590000116
之间的关系的参数，其具体数值通常是根据实验现象来调节，选取一个合适的数值以达到想要的结果。
子问题1用来利用可见边界处的梯度值快速恢复板状物的结构，但在子问题1求解时，是对图像沿x方向逐行进行求解，此时有可能会打乱图像相邻行灰度值之间的连续性，从而引入沿x方向的条状伪影。子问题2用来消除子问题1求解过程中引入的条状伪影。上述子问题拆分方式及求解算法，是基于目标函数的实际意义，而不是根据严格的数学推导得来的。
在一个实施例中，子问题1的求解方法如下：
求解(17)式表示的优化问题，新的迭代点为：
Figure BDA0002717721590000117
根据傅里叶变换微分性质，式(19)可转化为：
Figure BDA0002717721590000118
其中，ωx表示ω沿x方向的频率，在一维的情况下，ωx与ω等价。
本文算法的关键是选取合适的参数κ，以求得相对准确的
Figure BDA0002717721590000119
Figure BDA00027177215900001110
计算公式如下，
Figure BDA00027177215900001111
其中，参数κ根据x方向的梯度图像可见边界处的灰度值确定。
令
Figure BDA0002717721590000121
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的结果与
Figure BDA0002717721590000122
的结果的灰度值近似为倍数关系，本发明实施例通过(24)式对g进行简单处理，使g的结果与
Figure BDA0002717721590000123
的结果灰度值量级相同，再进行后续计算g。
Figure BDA0002717721590000124
式中，
Figure BDA0002717721590000125
表示fk在(i,j)处的像素值，gi,j表示g在(i,j)处的像素值。
此时，子问题1的解可表示如下：
Figure BDA0002717721590000126
参数λ1越大，算法收敛越慢，为了加速优化问题(20)的收敛速度，本发明方法采用Fixed point continuation(FPC)技术，令参数λ1随迭代次数由小到大逐渐变化。
子问题2的求解方法如下：根据如上所述，正则化项Rs(f)的选取方式不同，对应不同的子问题2的求解算法。在板状物CT重建时，用沿y方向的中值滤波的结果作为子问题2的解，相较于其它几种正则化方式，更能有效地去除子问题1引入的条状伪影。
fk+1＝medfilter1d(fk+1/2) (26)其中，medfilter1d(fk+1/2)表示对fk+1/2做沿y方向的中值滤波。
以下通过一个具体实施例来说明本发明提出的一种基于傅里叶变换微分性质的板状物CT图像快速重建方法的具体实现过程。
图2a为本具体实施例中用作测试模体的三维模拟PCB模体，截取该数据的一层作为本实施例中用作测试模体的二维模拟PCB模体，如图2b所示。
实验扫描参数设置如下：射线源到转台中心的距离为600mm，射线源到探测器的距离为1000mm。线探测器由512个探测器单元组成，每个探测器单元的尺寸为0.25mm。本实施例中，令α为射线源和转台中心的连线与转台水平方向的夹角，本发明实施例将二维模拟PCB模体沿α∈[π/6,5π/6]进行扫描，投影数据用光线投射算法作用于模体数据上生成，如图3所示。
具体实施步骤如下：
步骤1：为板状物CT图像赋初值f0＝0；
步骤2：假定经过k次迭代后，有板状物CT图像的估计值fk。利用公式(14)、(15)和(16)使板状物CT图像的估计值满足优化模型的约束条件。W有多种选取方式，如泊松噪声模型、非平稳高斯噪声模型、考虑系统电子噪声的噪声模型等。本实施例中，数据保真项权重W选取时，选择考虑系统电子噪声的噪声模型，即
Figure BDA0002717721590000131
Figure BDA0002717721590000132
其中，
Figure BDA0002717721590000133
表示第p个探测器在第q个扫描角度下对应的投影数据的标准差，N0p表示在不放物体的情况下第p个探测器探测到的光子数，
Figure BDA0002717721590000134
是探测器上电子噪声的方差,
Figure BDA0002717721590000135
表示第p个探测器在第q个扫描角度下下对应的投影数据的值。其中计算所需的CT系统空扫时探测器的计数N0p根据实际扫描条件获得，电子噪声方差较难获得精确值，本实施例中，本发明实施例取
Figure BDA0002717721590000141
步骤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表示二维数字图像的行数或列数，其元素fn表示被测物体在第n个像素处的线性衰减系数；p为探测器在不同投影角下测量的投影数据，M表示CT扫描的总的射线数；A为正投影矩阵，其元素am,n表示第m条射线与第n个像素的交线长；
步骤2：为板状物CT图像赋初值，作为板状物CT图像的估计值；
步骤3：利用凸集投影算法使板状物CT图像的估计值满足优化模型的约束条件；其中，所述优化模型包括：
(一)目标函数，其表示为式(7)：
Figure FDA0002717721580000011
(二)约束条件，其包括：
数据保真项约束，其表示为式(8)：
Figure FDA0002717721580000012
非负约束，其表示为式(9)：
f≥0 (9)
边界约束，其表示为式(10)：
Figure FDA0002717721580000017
上述各式中，f*表示优化模型最终求出的最优解，Dx为沿x方向的梯度算子，
Figure FDA0002717721580000013
为满足约束条件的板状物CT图像的估计值，
Figure FDA0002717721580000014
为待求图像沿x方向的近似梯度
Figure FDA0002717721580000015
其中κ为阈值参数；Rs(f)为正则化项，用于去除求解过程中引入的伪结构；W为所述步骤1中投影数据的权重矩阵；
Figure FDA0002717721580000016
表示二范数的平方；ε为与投影数据噪声水平相关的数，实际应用时是优化方程的误差门限；Ω表示整幅二维图像；
步骤4：利用傅里叶变换微分性质等算法求解优化模型中的目标函数；
步骤5：重复步骤3和步骤4，直到满足终止条件。
2.如权利要求1所述的基于傅里叶变换微分性质的板状物CT图像快速重建方法，其特征在于，所述步骤3包括：
步骤3.1，利用凸集投影算法使板状物CT图像的估计值满足式(8)表示的数据保真项约束，其包括：
采用凸集投影方法，使重建结果满足约束条件的要求，假定fk表示第k次迭代的结果，当前图像fk在
Figure FDA0002717721580000021
上的投影∏(fk)表示为：
∏(fk)＝fk+τAT(AAT)-1W(p-Afk) (13)
其中，参数τ＞0；
步骤3.2：根据式(9)表示的非负约束，利用式(15)和(16)求解
Figure FDA0002717721580000022
Figure FDA0002717721580000023
令
Figure FDA0002717721580000024
其中，fk i,j表示当前图像fk中的第i行，第j列的像素。
3.如权利要求2所述的基于傅里叶变换微分性质的板状物CT图像快速重建方法，其特征在于，采用如式(14)所示的逐射线投影方法求解式(13)：
Figure FDA0002717721580000025
符号:＝为赋值运算，即将该符号右边的结果赋给左边的变量。
4.如权利要求2所述的基于傅里叶变换微分性质的板状物CT图像快速重建方法，其特征在于，
Figure FDA0002717721580000026
Figure FDA0002717721580000027
其中，
Figure FDA0002717721580000031
表示第p个探测器在第q个扫描角度下对应的投影数据的标准差，N0p表示在不放物体的情况下第p个探测器探测到的光子数，
Figure FDA0002717721580000032
是探测器上电子噪声的方差，
Figure FDA0002717721580000033
表示第p个探测器在第q个扫描角度下下对应的投影数据的值。
5.如权利要求1至4中任一项所述的基于傅里叶变换微分性质的板状物CT图像快速重建方法，其特征在于，所述步骤4中，求解优化模型的目标函数时，将式(6)表示的优化问题转化为如下两个子问题，以更新fk+1：
子问题1：
Figure FDA0002717721580000034
子问题2：
Figure FDA0002717721580000035
式中，λ1表示用于平衡
Figure FDA0002717721580000036
和
Figure FDA0002717721580000037
之间的关系的参数，λ2表示用于平衡Rs(f)和
Figure FDA0002717721580000038
之间的关系的参数。
6.如权利要求5所述的基于傅里叶变换微分性质的板状物CT图像快速重建方法，其特征在于，子问题1中，求解式(17)表示的优化问题，新的迭代点为：
Figure FDA0002717721580000039
根据傅里叶变换微分性质，式(19)转化为：
Figure FDA00027177215800000310
Figure FDA00027177215800000311
的计算公式如式(21)：
Figure FDA0002717721580000041
其中，参数κ根据x方向的梯度图像可见边界处的灰度值确定；
令
Figure FDA0002717721580000042
由(22)式算出的g，经过(23)式修正后的结果作为新的g：
g:＝g-repmat(g(:,1),1,K) (23)
其中，g(:,1)表示图像g的第一列，repmat(g(:,1),1,K)表示将图像g的第一列复制K份，组成一个K×K大小的图像。
7.如权利要求6所述的基于傅里叶变换微分性质的板状物CT图像快速重建方法，其特征在于，子问题1求解过程中，当ωx＝0时，令ωx＝10-10，并通过式(24)，使g的结果与
Figure FDA0002717721580000043
的结果灰度值量级相同：
Figure FDA0002717721580000044
令参数λ1随迭代次数由小到大逐渐变化，
Figure FDA0002717721580000045
表示fk在(i,j)处的像素值，gi,j表示g在(i,j)处的像素值；
子问题1的解表示为式(25)：
Figure FDA0002717721580000046
8.如权利要求5所述的基于傅里叶变换微分性质的板状物CT图像快速重建方法，其特征在于，在板状物CT重建时，子问题2的解采用将式(11)或式(12)所示的Rs(f)的表达式代入(18)式后分别对应的求解算法，或者采用沿y方向的中值滤波；
Figure FDA0002717721580000051
Figure FDA0002717721580000052
其中，fx,y表示在点(x,y)处f的值；#{}是计数算子，在式(11)中计算满足|fx,y-fx,y-1|≠0的像素点的数目，在式(11)中计算满足|fx,y-fx-1,y|+|fx,y-fx,y-1|≠0的像素点的数目。
9.如权利要求5所述的基于傅里叶变换微分性质的板状物CT图像快速重建方法，其特征在于，为板状物CT图像赋初值f0＝0。
CN202011077364.1A 2020-10-10 2020-10-10 基于傅里叶变换微分性质的板状物ct图像快速重建方法 Active CN112200882B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011077364.1A CN112200882B (zh) 2020-10-10 2020-10-10 基于傅里叶变换微分性质的板状物ct图像快速重建方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011077364.1A CN112200882B (zh) 2020-10-10 2020-10-10 基于傅里叶变换微分性质的板状物ct图像快速重建方法

Publications (2)

Publication Number Publication Date
CN112200882A true CN112200882A (zh) 2021-01-08
CN112200882B CN112200882B (zh) 2022-11-25

Family

ID=74013933

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011077364.1A Active CN112200882B (zh) 2020-10-10 2020-10-10 基于傅里叶变换微分性质的板状物ct图像快速重建方法

Country Status (1)

Country Link
CN (1) CN112200882B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112712572A (zh) * 2021-01-11 2021-04-27 明峰医疗系统股份有限公司 Ct扫描设备的低信号噪声的抑制方法、系统及计算机可读存储介质

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20080273651A1 (en) * 2007-05-05 2008-11-06 Franz Edward Boas Methods and apparatus for reducing artifacts in computed tomography images
US20110007864A1 (en) * 2007-04-18 2011-01-13 Agency For Science, Technology And Research Method and apparatus for reorientated resconstruction of computed tomography images of planar objects
CN105023282A (zh) * 2014-04-30 2015-11-04 华中科技大学 一种基于压缩感知的稀疏投影超声ct图像重建方法
CN106846427A (zh) * 2017-01-25 2017-06-13 浙江大学 一种基于重加权各向异性全变分的有限角度ct重建方法
CN108010099A (zh) * 2017-12-04 2018-05-08 首都师范大学 一种x射线多能谱ct有限角扫描和图像迭代重建方法
CN109085190A (zh) * 2018-08-10 2018-12-25 首都师范大学 一种针对板状物的x射线三维ct数据扫描系统及其扫描方法

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20110007864A1 (en) * 2007-04-18 2011-01-13 Agency For Science, Technology And Research Method and apparatus for reorientated resconstruction of computed tomography images of planar objects
US20080273651A1 (en) * 2007-05-05 2008-11-06 Franz Edward Boas Methods and apparatus for reducing artifacts in computed tomography images
CN105023282A (zh) * 2014-04-30 2015-11-04 华中科技大学 一种基于压缩感知的稀疏投影超声ct图像重建方法
CN106846427A (zh) * 2017-01-25 2017-06-13 浙江大学 一种基于重加权各向异性全变分的有限角度ct重建方法
CN108010099A (zh) * 2017-12-04 2018-05-08 首都师范大学 一种x射线多能谱ct有限角扫描和图像迭代重建方法
CN109085190A (zh) * 2018-08-10 2018-12-25 首都师范大学 一种针对板状物的x射线三维ct数据扫描系统及其扫描方法

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
F. ZHANG ET AL.: "An image reconstruction strategy for truncated projections of planar object in cone-beam CT", 《2015 12TH INTERNATIONAL CONFERENCE ON FUZZY SYSTEMS AND KNOWLEDGE DISCOVERY (FSKD)》 *
LIU, TONG AND XU, JIAN: "Differential Reconstruction for Planar Object in Computed Tomography", 《JOURNAL OF X-RAY SCIENCE AND TECHNOLOGY》 *
XIAO XUE,SHUSEN ZHAO,YUNSONG ZHAO,PENG ZHANG. MEASUREMENT SCIEN: "Image reconstruction for limited-angle computed tomography with curvature constraint", 《MEASUREMENT SCIENCE AND TECHNOLOGY》 *
刘华,潘晋孝,陈平: "板状物成像方法研究", 《核电子学与探测技术》 *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112712572A (zh) * 2021-01-11 2021-04-27 明峰医疗系统股份有限公司 Ct扫描设备的低信号噪声的抑制方法、系统及计算机可读存储介质
CN112712572B (zh) * 2021-01-11 2023-10-24 明峰医疗系统股份有限公司 Ct扫描设备的低信号噪声的抑制方法、系统及计算机可读存储介质

Also Published As

Publication number Publication date
CN112200882B (zh) 2022-11-25

Similar Documents

Publication Publication Date Title
Calvetti et al. Hypermodels in the Bayesian imaging framework
EP3814759B1 (en) Item inspection by radiation imaging using an iterative projection-matching approach
CN107328798B (zh) 一种新型icl系统及实现方法
CN112819723B (zh) 一种高能x射线图像盲复原方法及系统
CN102590248B (zh) 平移式微焦ct检测装置在线检测电子元件的方法
Xu et al. Deep residual learning in CT physics: scatter correction for spectral CT
Shtok et al. Sparsity-based sinogram denoising for low-dose computed tomography
Zhao et al. Edge information diffusion-based reconstruction for cone beam computed laminography
CN112200882B (zh) 基于傅里叶变换微分性质的板状物ct图像快速重建方法
US20210262947A1 (en) Method and device for acquiring tomographic image data by oversampling, and control program
CN110415307A (zh) 一种基于张量补全的多能ct成像方法、装置及其存储设备
CN112102428A (zh) Ct锥形束扫描图像重建方法、扫描系统计及存储介质
Dabravolski et al. Dynamic angle selection in x-ray computed tomography
EP3951434A1 (en) Estimating background radiation from unknown sources
Schut et al. Top-ct: Trajectory with overlapping projections x-ray computed tomography
Bappy et al. Modified simultaneous iterative reconstruction technique for fast, high‐quality CT reconstruction
CN110717959B (zh) 基于曲率约束的x射线有限角ct图像重建方法和装置
CN105608719B (zh) 一种基于两阶段投影调整的快速ct图像重建方法
CN115359049B (zh) 基于非线性扩散模型的有限角ct图像重建方法及装置
Venkatakrishnan et al. Model-based iterative reconstruction for neutron laminography
US20140169651A1 (en) Shape-Extraction Method and Shape-Extraction System
Vengrinovich Bayesian image and pattern reconstruction from incomplete and noisy data
CN116188615A (zh) 一种基于正弦域和图像域的稀疏角度ct重建方法
JP5493072B2 (ja) Ｃｔ装置、ｃｔ装置における画像再構成方法、及び電子回路部品
Ziabari et al. X-ray ct reconstruction of additively manufactured parts using 2.5 d deep learning mbir

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant