CN104867157B - 一种ct探测器像素响应不一致性的校正方法 - Google Patents

一种ct探测器像素响应不一致性的校正方法 Download PDF

Info

Publication number
CN104867157B
CN104867157B CN201510293904.2A CN201510293904A CN104867157B CN 104867157 B CN104867157 B CN 104867157B CN 201510293904 A CN201510293904 A CN 201510293904A CN 104867157 B CN104867157 B CN 104867157B
Authority
CN
China
Prior art keywords
mrow
pixel
response
mover
msup
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.)
Active
Application number
CN201510293904.2A
Other languages
English (en)
Other versions
CN104867157A (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.)
Beijing Wandong Medical Polytron Technologies Inc
Original Assignee
Beijing Wandong Medical Polytron Technologies Inc
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 Beijing Wandong Medical Polytron Technologies Inc filed Critical Beijing Wandong Medical Polytron Technologies Inc
Priority to CN201510293904.2A priority Critical patent/CN104867157B/zh
Publication of CN104867157A publication Critical patent/CN104867157A/zh
Application granted granted Critical
Publication of CN104867157B publication Critical patent/CN104867157B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T3/00—Geometric image transformations in the plane of the image
    • G06T3/40—Scaling of whole images or parts thereof, e.g. expanding or contracting
    • G06T3/4053—Scaling of whole images or parts thereof, e.g. expanding or contracting based on super-resolution, i.e. the output image resolution being higher than the sensor resolution
    • G06T3/4076—Scaling of whole images or parts thereof, e.g. expanding or contracting based on super-resolution, i.e. the output image resolution being higher than the sensor resolution using the original low-resolution images to iteratively correct the high-resolution images
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/80—Geometric correction
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/30—Determination of transform parameters for the alignment of images, i.e. image registration
    • 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)
  • Computer Vision & Pattern Recognition (AREA)
  • Apparatus For Radiation Diagnosis (AREA)

Abstract

本发明公开了一种CT探测器像素响应不一致性的校正方法，该方法包括：将原始的投影正弦图在多个投影角内进行平均，获得初始平均响应曲线；去除初始平均响应曲线中的轮廓线；在去除轮廓线的平均响应曲线中探测响应不一致的像素点，并将响应不一致的像素点分为两类：单双点，以及连续三点及以上；对于单双点的响应不一致像素点使用邻域中值滤波的方法估计偏差；对于连续三点及以上的响应不一致像素点使用基于平滑性约束的填充方法估计偏差；根据估计出的偏差补偿原始的投影正弦图，完成CT探测器像素响应不一致性的校正。本发明公开的方法，可以有效地消除由于探测器像素响应不一致性产生的伪像。

Description

一种CT探测器像素响应不一致性的校正方法
技术领域
本发明涉及医学图像处理技术领域，尤其涉及一种CT探测器像素响应不一致性的校正方法。
背景技术
如图1所示，一台螺旋CT机通常由如下部件构成：X线管球、CT探测器、机架、数据采集系统。在扫描过程中，机架带动X线管球和CT探测器进行旋转的同时，X线管球进行曝光，CT探测器采集扫描物体在各个角度下的投影数据，并利用工作站重建获得表示物体X线衰减系数的断层图像。
CT探测器由光电阵列模块在支撑结构件上拼接而成。支撑结构件成弧形，弧形的圆心是管球的焦点，每个光电模块通过插槽和螺钉固定在支撑结构上，从而连接成一个弧形的探测器。光电阵列模块的结构如图2所示，主要由滤线器，闪烁体和硅光电二极管阵列构成。滤线器由一系列薄的钨片排列成梳子状的结构，用于滤除由其他方向传播过来的散射光子。闪烁体的作用是将X线光子转化为可见光子。入射X光子与闪烁体作用，激发出可见光子，硅光电二极管将可见光转换为电流信号，数据采集系统再对电流信号进行采集。
在理想情况下，探测器的各个像素对X光子的响应是一致的，即在相同的X光子入射的情况下，探测器各个像素的输出不会随着X光子入射探测器的位置而改变。但是在实际情况中，探测器响应的一致性是无法达到的。一方面在探测器制造过程中，不可能制造出完全一致的探测器像素。例如，由于机械加工和安装的误差，滤线器钨片之间的距离不能完全一致。光电阵列模块每个像素的光电二级管的响应也可能有所差别。另一方面，在长时间的X线照射情况下，探测器的某些像素会严重退化，会在探测器像素的输出信号上产生一个不稳定的偏置，或是探测器像素的响应与输入的X光子数呈非线性关系。这些都会导致一些探测器像素的响应与其他像素的响应明显不同。
探测器像素响应的不一致性通常会导致重建图像中的环带状伪像。现有的校正方法主要分为如下两类：
一类校正是在原始的投影正弦图数据中进行。首先标记出探测器中这些响应非正常的像素，使用模体测量出探测器在不同X光子入射量下的投影误差。误差根据一定的模型估计出当前投影数据的误差，对这些探测器像素的响应进行补偿。这类方法的缺点是，需要经常性的对像素响应的偏差进行标定，并且对于像素响应偏差不稳定的情况不能有效的补偿。
另一类校正方法是在重建图像中进行。在重建图像中检测同心的环或者弧，测定这些环的强度，并从重建图像中减去。这类方法的缺点在于，有时环状伪像和真实的解剖结构的边界是难于区分的，在去除环状伪像的过程中，也会损失真实边界信息。
发明内容
本发明的目的是提供一种CT探测器像素响应不一致性的校正方法，可以有效地消除由于探测器像素响应不一致性产生的伪像。
本发明的目的是通过以下技术方案实现的：
一种CT探测器像素响应不一致性的校正方法，该方法包括：
将原始的投影正弦图在多个投影角内进行平均，获得初始平均响应曲线；
去除初始平均响应曲线中的轮廓线；
在去除轮廓线的平均响应曲线中探测响应不一致的像素点，并将响应不一致的像素点分为两类：单双点，以及连续三点及以上；
对于单双点的响应不一致像素点使用邻域中值滤波的方法估计偏差；
对于连续三点及以上的响应不一致像素点使用基于平滑性约束的填充方法估计偏差；
根据估计出的偏差补偿原始的投影正弦图，完成CT探测器像素响应不一致性的校正。
所述将原始的投影正弦图在多个投影角内进行平均，获得初始平均响应曲线包括：
沿投影数n方向将原始的投影正弦图P(m,n)进行平均滤波，得到再将沿投影数n方向进行下采样得到初始平均响应曲线
其中，m表示探测器像素点坐标，n表示投影数。
所述去除初始平均响应曲线中的轮廓线包括：
沿探测器像素点坐标m方向使用savitzky-golay滤波器进行滤波获取轮廓线Z(m,n)，表示为：
将获取轮廓线Z(m,n)从初始平均响应曲线中减去，再沿投影数n方向进行一次低频滤波，获得去除轮廓线的平均响应曲线R，其中的像素点记为R(m,n)，表示为：
所述在去除轮廓线的平均响应曲线中探测响应不一致的像素点包括：
利用每一像素点与其领域像素点平均值的差异来探测响应不一致的像素点，表示为：
其中，R(l,n)表示R(m,n)的相邻的2W个像素点，2W+1表示领域像素点的窗宽；
当T(m,n)大于阈值Th，则对应像素点为响应不一致的像素点。
所述对于单双点的响应不一致像素点使用邻域中值滤波的方法估计偏差包括：
对于单像素点响应不一致以及连续两个像素点响应不一致的情况，使用7点中值滤波的方法对进行滤波，计算出响应不一致像素点的响应期望值，表示为：
其中，p表示响应不一致像素点的坐标，n表示投影数，median表示取数组的中值；
再计算对应的偏差值，表示为：
所述对于连续三点及以上的响应不一致像素点使用基于平滑性约束的填充方法估计偏差包括：
在初始平均响应曲线中标记连续三点及以上的响应不一致像素点，并使用线性差值对所标记的响应不一致像素点进行填充；
将标记的响应不一致像素点作为坏像素区Φ，与之相邻的5个像素点作为过渡区Ψ，并使用之前的线性差值填充时的结果作为初始值，求解下面的最优化问题，获得响应不一致像素点的响应期望值：
其中，I(p,n)表示最优化问题的变量，▽I(p,n)表示对I(p,n)求梯度，p表示响应不一致像素点的坐标，n表示投影数，β为常数；
再计算对应的偏差值，表示为：
求解最优化问题时，使用了最速梯度下降法，f(I(p,n))的梯度的表达式为：
所述根据估计出的偏差补偿原始的投影正弦图包括：
将获得的偏差值进行1：8的线性插值，并叠加回原始的投影正弦图。
由上述本发明提供的技术方案可以看出，将响应不一致的像素点分为两类，根据类别不同使用不同的方法来准确计算其偏差值，再对原始投影数据进行补偿校正，不仅可以消除环形伪影，也可以消除像素响应不一致带来的其它伪影；并且由于在处理过程中不需要对探测器的异常像素进行标定，所以该校正方法对于响应偏差稳定和不稳定的探测器像素都是有效的。
附图说明
为了更清楚地说明本发明实施例的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域的普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他附图。
图1为本发明背景技术提供的螺旋CT机的结构示意图；
图2为本发明背景技术提供的光电阵列模块的结构示意图；
图3为本发明实施例提供的一种CT探测器像素响应不一致性的校正方法的流程图；
图4为本发明实施例提供的去除轮廓线前后的平均响响应曲线的示意图；
图5为本发明实施例提供的将平均响应曲线图分成三类区域的示意图；
图6为本发明实施例提供的响应不一致点表现形式的示意图。
具体实施方式
下面结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明的保护范围。
本发明实施例使用的校正方法是针对投影正弦图进行处理的方法。投影正弦图中，每个像素的值表示了对应X线传播路线上物质的衰减系数的积分。在相同的X光入射的情况下，响应不一致像素的输出的信号往往会显著低于或者高于相邻的像素，所以可以用像素和相邻像素的差异来检测出这些响应不一致的像素。但是由于噪声和扫描物体的结构的影响，直接在正弦图中检测这些异常像素是非常困难的。
本发明实施例中，将探测器的响应在多个投影角范围内进行平均，可以获得一个平均响应曲线(mean curve)。在平均响应曲线中，噪声和扫描物体结构的影响可以被减弱，从而突出探测器像素响应的不一致性。这样做的好处是可以将结构中的强边缘造成的响应差异和像素自身引入的响应偏差区分开来。对于单个或者连续两个响应不一致的异常像素，可以使用邻域像素中值滤波来估计像素在部分投影角范围内的响应偏差。在获得了这些像素响应偏差后，可以补偿原始投影正弦图，实现无伪影的图像重建。
对于连续3个及3个以上的不一致像素使用中值滤波的方法估计像素期望值是不准确的。一方面，由于不一致像素较多，本身就会使中值滤波的结果引入一定的偏差。另一方面，在靠近中心的位置像素，是由远离该点的像素估计古来的，估计的准确性也大为下降。所以在本发明实施例中使用基于平滑性约束的填充的方法来估计不一致像素的期望值，进而获取像素偏差。
由于本方法针对原始投影数据进行校正，所以不仅可以消除环形伪影，也可以消除像素响应不一致带来的其它伪影。并且由于在处理过程中不需要对探测器的异常像素进行标定，所以本校正方法对于响应偏差稳定和不稳定的探测器像素都是有效的。
实施例
本发明实施例提供一种CT探测器像素响应不一致性的校正方法。如图3所示，该方法主要包括如下步骤：
步骤11、将原始的投影正弦图在多个投影角内进行平均，获得初始平均响应曲线；
步骤12、去除初始平均响应曲线中的轮廓线；
步骤13、在去除轮廓线的平均响应曲线中探测响应不一致的像素点，并将响应不一致的像素点分为两类：单双点，以及连续三点及以上；
步骤14、对于单双点的响应不一致像素点使用邻域中值滤波的方法估计偏差；
步骤15、对于连续三点及以上的响应不一致像素点使用基于平滑性约束的填充方法估计偏差；
步骤16、根据估计出的偏差补偿原始的投影正弦图，完成CT探测器像素响应不一致性的校正。
在上述步骤12中，为了准确的检测出响应不一致点，使用了去除轮廓线的方法。由于认为每一点和相邻点的像素偏差是一个随时间缓慢变化的量，所以在估计偏差前，先对每个像素在投影数方向上作平均滤波获得初始的平均响应曲线。由于扫描物体通常是近似一个椭圆形结构，在两侧衰减厚度小，在中间衰减大。在初始的平均响应曲线中，像素值有一个快速上升和下降的趋势，称之为轮廓线。而像素响应的偏差非常小，会被这种像素值的迅速变化趋势覆盖掉。因此在检测不一致的像素点时，首先将平均响应曲线中的轮廓线减去。
在求轮廓线的过程中，使用了savitzky-golay滤波器对初始的平均响应曲线进行滤波获取轮廓线。图4a-图4b显示了去除轮廓线之前(图4a)和之后(图4b)的平均响响应曲线。从图中可以看出来，像素点的非一致性造成的偏差被明显的凸显出来，方便进行随后的检测。在随后的非一致性像素的检测中，通过的简单的阈值方法就能将这些点检测出来。
在上述步骤14中，对于单点及双点的响应不一致点，其像素的响应总是明显高于或者低于相邻的像素，相当于产生了一个“突变”。本实施例中用相邻点的中值滤波值作为该响应不一致点的响应期望值可以很好的消除此类“突变”。
但是多于连续三点及以上的不一致点，相邻点的中值滤波就不能准确的估计出响应期望值。此时的中值滤波的结果仍然是有偏差的。连续多点的响应不一致点在重建图像中产生宽带状的环状伪影。在上述步骤15中，使用了基于平滑性约束像素填充的技术。如图5所示，将平均响应曲线图分成三类区域，一类是需要估计偏差值的坏像素区，记为Φ，坏像素与好像素过度区，记为Ψ，好像素区，记为Ω。像素点的坐标记为(x,y)，需要估计像素点的灰度值记为I(x,y)，像素的原始灰度值记为I0(x,y).估计期望值的算法可以表示为一个最小化估计问题：
上述最小化估计问题包含了两个目标函数。一个目标函数是希望在坏像素区Φ和过渡区Ψ的梯度和最小，另一个目标函数是希望在过渡区Ψ和原始图I0的差最小。通过最优化过程来使填充的坏像素区Φ自然过渡到好像素区Ω，并且Φ区也是平滑的。上述的参数β是一个常数(系数)，用来调节上面的两个目标函数的加权比例，求解过程使用迭代方法来达到最优化。填充的结果避免了坏像素区Φ与相邻像素连接时出现不连续的情况，同时也相当于进行了一次平滑的滤波，消除了Φ区域内响应不一致造成的像素值突变。将Φ填充的结果I(x,y)作为像素响应的期望值，然后将I(x,y)-I0(x,y)作为坏像素区Φ的偏差。最终偏差被叠加到原始的投影图，完成对投影图的补偿。
为了便于理解本发明，下面结合一具体的示例来做进一步的说明。需要说明的是，在下述示例中所涉及的某些具体的参数仅为举例说明并非构成限制，用户还可以根据实际需求来更换为其他合适的数值。
本示例中，处理步骤主要包括：
步骤1、沿投影数n方向将原始的投影正弦图P(m,n)进行平均滤波，得到表示为：
其中，m表示探测器像素点坐标，n表示投影数。
本示例中，平均滤波器的长度(M+1)选为65；也可以根据实际情况进行调整。
步骤2、将沿投影数n方向进行下采样得到初始平均响应曲线
因为在投影数n方向是响应是低频的，下采样基本不会影响后续运算的精度，同时可以极大的降低后续的运算量，提高运算速度。本示例使用的下采样率为8：1也可根据具体情况调整。
步骤3、获取初始平均响应曲线中的轮廓线。
本实例中，沿探测器像素点坐标m方向使用savitzky-golay滤波器进行滤波获取轮廓线Z(m,n)，表示为：
其中，w表示滤波器数组的的第w个系数，如9个系数的滤波的数组，w的取值从-4到4。
savitzky-golay(SG)滤波器是一种基于局域多项式最小二乘拟合的滤波方法，这种滤波的最大优势是滤除噪声的同时保持信号的形状、结构和宽度等分布特性。经过SG滤波后，保留了平均响应曲线的轮廓和结构，消除轮廓线Z(m,n)后主要保留的就是噪声和像素响应的“突变”。
本示例中使用的SG滤波器长度为9点，滤波器的系数表示为：
SG＝[0.0350 -0.1282 0.0699 0.3147 0.4172 0.3147 0.0699 -0.12820.0350]
步骤4、去除轮廓线。
将轮廓线Z(m,n)从初始平均响应曲线中减去；去除轮廓线后，平均曲线中保留的主要就是噪声和像素偏差；为了进一步降低噪声对偏差的影响，可以将再沿投影数n方向再进行一次低频滤波，获得去除轮廓线的平均响应曲线R，其中的像素点记为R(m,n)，其噪声被进一步的平滑，响应不一致造成的响应偏差就被突显出来了，便于后续对这些点的检测。
上述处理表示为：
步骤5、探测响应不一致像素点。
在去除了轮廓线的响应R(m,n)中，响应不一致点表现为如下图6a-图6c所示的三种形式：如图6a所示，为单点的“突变”，通常是该像素的暗电流异常造成的。如图6b所示，为连续两点的响应不一致，这种情况通常是由滤线器中像素分隔的机械误差造成的。如图6c所示，为连续三点以上的响应不一致，这通常是由于探测器模块的退化造成像素响应的非线性造成的。上述三类点响应的在R(m,n)都表示为显著的高于或低于领域的平均值。因此我们可以用每一像素点与其领域像素点平均值的差异来检测这些相应不一致的点；定义如下判定参数：
当T(m,n)大于阈值Th，则对应像素点为响应不一致的像素点；其中，R(l,n)表示R(m,n)的相邻的2W个像素点，2W+1表示领域像素点的窗宽(窗宽表示了邻域包含的像素点的个数)，本示例中选为11点；阈值Th选为领域像素点平均值的3倍。
步骤6、对于单像素点响应不一致以及连续两个像素点响应不一致的情况，使用7点中值滤波的方法对进行滤波，计算出响应不一致像素点的响应期望值及偏差值。对于不一致像素点R(p,n)，其响应期望值计算公式为：
其中，p表示响应不一致像素点的坐标，median表示取数组的中值；
再计算对应的偏差值，其公式为：
步骤7、对于连续三点及以上的响应不一致像素点，在初始平均响应曲线中标记连续三点及以上的响应不一致像素点，并使用线性差值对所标记的响应不一致像素点进行填充。
步骤8、将标记的响应不一致像素点作为坏像素区Φ，与之相邻的5个像素点作为过渡区Ψ，并使用步骤7的线性差值填充时的结果作为初始值(可极大的加速迭代的收敛过程)，求解下面的最优化问题，获得响应不一致像素点的响应期望值：
其中，I(p,n)表示最优化问题的变量，I(p,n)的目标解使函数f(I(p,n))达到最小值，▽I(p,n)表示对I(p,n)求梯度。本示例中，β为常数，其值可选为0.5；
优选的，求解上述最优化问题时，使用了最速梯度下降法，f(I(p,n))的梯度的表达式为：
再计算对应的偏差值，表示为：
步骤9、根据步骤6与步骤8计算出的偏差值补偿原始的投影正弦图。
例如，将获得的偏差值进行1：8的线性插值，并叠加回原始的投影正弦图。
本发明实施例将响应不一致的像素点分为两类，根据类别不同使用不同的方法来准确计算其偏差值，再对原始投影数据进行补偿校正，不仅可以消除环形伪影，也可以消除像素响应不一致带来的其它伪影；并且由于在处理过程中不需要对探测器的异常像素进行标定，所以该校正方法对于响应偏差稳定和不稳定的探测器像素都是有效的。
通过以上的实施方式的描述，本领域的技术人员可以清楚地了解到上述实施例可以通过软件实现，也可以借助软件加必要的通用硬件平台的方式来实现。基于这样的理解，上述实施例的技术方案可以以软件产品的形式体现出来，该软件产品可以存储在一个非易失性存储介质(可以是CD-ROM，U盘，移动硬盘等)中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本发明各个实施例所述的方法。
以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明披露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求书的保护范围为准。

Claims (8)

1.一种CT探测器像素响应不一致性的校正方法，其特征在于，该方法包括：
将原始的投影正弦图在多个投影角内进行平均，获得初始平均响应曲线；
去除初始平均响应曲线中的轮廓线；
在去除轮廓线的平均响应曲线中探测响应不一致的像素点，并将响应不一致的像素点分为两类：单双点，以及连续三点及以上；
对于单双点的响应不一致像素点使用邻域中值滤波的方法估计偏差；
对于连续三点及以上的响应不一致像素点使用基于平滑性约束的填充方法估计偏差；
根据估计出的偏差补偿原始的投影正弦图，完成CT探测器像素响应不一致性的校正。
2.根据权利要求1所述的方法，其特征在于，所述将原始的投影正弦图在多个投影角内进行平均，获得初始平均响应曲线包括：
沿投影数n方向将原始的投影正弦图P(m,n)进行平均滤波，得到再将沿投影数n方向进行下采样得到初始平均响应曲线
其中，m表示探测器像素点坐标，n表示投影数。
3.根据权利要求2所述的方法，其特征在于，所述去除初始平均响应曲线中的轮廓线包括：
沿探测器像素点坐标m方向使用savitzky-golay滤波器进行滤波获取轮廓线Z(m,n)，表示为：
<mrow> <mi>Z</mi> <mrow> <mo>(</mo> <mi>m</mi> <mo>,</mo> <mi>n</mi> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>w</mi> <mo>=</mo> <mo>-</mo> <mi>L</mi> </mrow> <mi>L</mi> </munderover> <mover> <mi>P</mi> <mo>^</mo> </mover> <mrow> <mo>(</mo> <mi>m</mi> <mo>-</mo> <mi>w</mi> <mo>,</mo> <mi>n</mi> <mo>)</mo> </mrow> <mi>S</mi> <mi>G</mi> <mrow> <mo>(</mo> <mi>w</mi> <mo>)</mo> </mrow> <mo>;</mo> </mrow>
其中，SG(w)表示使用savitzky-golay滤波器第w个系数进行滤波，平均滤波器的长度为M+1；
将获取轮廓线Z(m,n)从初始平均响应曲线中减去，再沿投影数n方向进行一次低频滤波，获得去除轮廓线的平均响应曲线R，其中的像素点的值记为R(m,n)，表示为：
<mrow> <mover> <mi>R</mi> <mo>^</mo> </mover> <mrow> <mo>(</mo> <mi>m</mi> <mo>,</mo> <mi>n</mi> <mo>)</mo> </mrow> <mo>=</mo> <mover> <mi>P</mi> <mo>^</mo> </mover> <mrow> <mo>(</mo> <mi>m</mi> <mo>,</mo> <mi>n</mi> <mo>)</mo> </mrow> <mo>-</mo> <mi>Z</mi> <mrow> <mo>(</mo> <mi>m</mi> <mo>,</mo> <mi>n</mi> <mo>)</mo> </mrow> <mo>;</mo> </mrow>
<mrow> <mi>R</mi> <mrow> <mo>(</mo> <mi>m</mi> <mo>,</mo> <mi>n</mi> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mi>n</mi> <mo>-</mo> <mi>M</mi> <mo>/</mo> <mn>2</mn> </mrow> <mrow> <mi>n</mi> <mo>+</mo> <mi>M</mi> <mo>/</mo> <mn>2</mn> </mrow> </munderover> <mover> <mi>R</mi> <mo>^</mo> </mover> <mrow> <mo>(</mo> <mi>m</mi> <mo>,</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>.</mo> </mrow>
4.根据权利要求3所述的方法，其特征在于，所述在去除轮廓线的平均响应曲线中探测响应不一致的像素点包括：
利用每一像素点与其邻域像素点平均值的差异来探测响应不一致的像素点，表示为：
<mrow> <mi>T</mi> <mrow> <mo>(</mo> <mi>m</mi> <mo>,</mo> <mi>n</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <mrow> <mn>2</mn> <mi>W</mi> <mo>+</mo> <mn>1</mn> </mrow> </mfrac> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>l</mi> <mo>=</mo> <mi>m</mi> <mo>-</mo> <mi>W</mi> </mrow> <mrow> <mi>l</mi> <mo>=</mo> <mi>m</mi> <mo>+</mo> <mi>W</mi> </mrow> </munderover> <mo>|</mo> <mi>R</mi> <mrow> <mo>(</mo> <mi>l</mi> <mo>,</mo> <mi>n</mi> <mo>)</mo> </mrow> <mo>|</mo> <mo>-</mo> <mo>|</mo> <mi>R</mi> <mrow> <mo>(</mo> <mi>m</mi> <mo>,</mo> <mi>n</mi> <mo>)</mo> </mrow> <mo>|</mo> <mo>;</mo> </mrow>
其中，R(l,n)表示R(m,n)的相邻的2W个像素点的值，2W+1表示邻域像素点的窗宽；
当T(m,n)大于阈值Th，则对应像素点为响应不一致的像素点。
5.根据权利要求1所述的方法，其特征在于，所述对于单双点的响应不一致像素点使用邻域中值滤波的方法估计偏差包括：
对于单像素点响应不一致以及连续两个像素点响应不一致的情况，使用7点中值滤波的方法进行滤波，计算出响应不一致像素点的响应期望值，表示为：
<mrow> <mover> <mi>R</mi> <mo>~</mo> </mover> <mrow> <mo>(</mo> <mi>p</mi> <mo>,</mo> <mi>n</mi> <mo>)</mo> </mrow> <mo>=</mo> <mi>m</mi> <mi>e</mi> <mi>d</mi> <mi>i</mi> <mi>a</mi> <mi>n</mi> <mrow> <mo>(</mo> <mo>&amp;lsqb;</mo> <mi>R</mi> <mo>(</mo> <mrow> <mi>p</mi> <mo>-</mo> <mn>3</mn> <mo>,</mo> <mi>n</mi> </mrow> <mo>)</mo> <mo>,</mo> <mi>R</mi> <mo>(</mo> <mrow> <mi>p</mi> <mo>-</mo> <mn>2</mn> <mo>,</mo> <mi>n</mi> </mrow> <mo>)</mo> <mo>...</mo> <mo>,</mo> <mi>R</mi> <mo>(</mo> <mrow> <mi>p</mi> <mo>+</mo> <mn>3</mn> <mo>,</mo> <mi>n</mi> </mrow> <mo>)</mo> <mo>&amp;rsqb;</mo> <mo>)</mo> </mrow> <mo>;</mo> </mrow>
其中，p表示响应不一致像素点的坐标，n表示投影数，median表示取数组的中值；
再计算对应的偏差值，表示为：
<mrow> <mi>E</mi> <mrow> <mo>(</mo> <mi>p</mi> <mo>,</mo> <mi>n</mi> <mo>)</mo> </mrow> <mo>=</mo> <mi>R</mi> <mrow> <mo>(</mo> <mi>p</mi> <mo>,</mo> <mi>n</mi> <mo>)</mo> </mrow> <mo>-</mo> <mover> <mi>R</mi> <mo>~</mo> </mover> <mrow> <mo>(</mo> <mi>p</mi> <mo>,</mo> <mi>n</mi> <mo>)</mo> </mrow> <mo>;</mo> </mrow>
其中，R(p,n)为不一致像素点的值。
6.根据权利要求3所述的方法，其特征在于，所述对于连续三点及以上的响应不一致像素点使用基于平滑性约束的填充方法估计偏差包括：
在初始平均响应曲线中标记连续三点及以上的响应不一致像素点，并使用线性差值对所标记的响应不一致像素点进行填充；
将标记的响应不一致像素点作为坏像素区Φ，与之相邻的5个像素点作为过渡区Ψ，并使用之前的线性差值填充时的结果作为初始值，求解下面的最优化问题，获得响应不一致像素点的响应期望值
<mrow> <mtable> <mtr> <mtd> <mrow> <mi>f</mi> <mrow> <mo>(</mo> <mi>I</mi> <mo>(</mo> <mrow> <mi>p</mi> <mo>,</mo> <mi>n</mi> </mrow> <mo>)</mo> <mo>)</mo> </mrow> <mo>=</mo> <munder> <mo>&amp;Sigma;</mo> <mrow> <mo>(</mo> <mi>p</mi> <mo>,</mo> <mi>n</mi> <mo>)</mo> <mo>&amp;Element;</mo> <mi>&amp;Phi;</mi> <mo>&amp;cup;</mo> <mi>&amp;Psi;</mi> </mrow> </munder> <mo>|</mo> <mo>&amp;dtri;</mo> <mi>I</mi> <mo>(</mo> <mrow> <mi>p</mi> <mo>,</mo> <mi>n</mi> </mrow> <mo>)</mo> <mo>|</mo> <mo>)</mo> <mo>+</mo> <mi>&amp;beta;</mi> <munder> <mo>&amp;Sigma;</mo> <mrow> <mo>(</mo> <mi>p</mi> <mo>,</mo> <mi>n</mi> <mo>)</mo> <mo>&amp;Element;</mo> <mi>&amp;Psi;</mi> </mrow> </munder> <msup> <mrow> <mo>(</mo> <mi>I</mi> <mo>(</mo> <mrow> <mi>p</mi> <mo>,</mo> <mi>n</mi> </mrow> <mo>)</mo> <mo>-</mo> <mover> <mi>P</mi> <mo>^</mo> </mover> <mo>(</mo> <mrow> <mi>p</mi> <mo>,</mo> <mi>n</mi> </mrow> <mo>)</mo> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mover> <mi>P</mi> <mo>&amp;OverBar;</mo> </mover> <msub> <mrow> <mo>(</mo> <mi>p</mi> <mo>,</mo> <mi>n</mi> <mo>)</mo> </mrow> <mrow> <mo>(</mo> <mi>p</mi> <mo>,</mo> <mi>n</mi> <mo>)</mo> <mo>&amp;Element;</mo> <mi>&amp;Phi;</mi> </mrow> </msub> <mo>=</mo> <mi>min</mi> <mrow> <mo>(</mo> <mi>f</mi> <mo>(</mo> <mrow> <mi>I</mi> <mrow> <mo>(</mo> <mrow> <mi>p</mi> <mo>,</mo> <mi>n</mi> </mrow> <mo>)</mo> </mrow> </mrow> <mo>)</mo> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> </mtable> <mo>;</mo> </mrow>
其中，I(p,n)表示最优化问题的变量，表示对I(p,n)求梯度，p表示响应不一致像素点的坐标，n表示投影数，β为常数；
再计算对应的偏差值，表示为：
<mrow> <mi>E</mi> <msup> <mrow> <mo>(</mo> <mi>p</mi> <mo>,</mo> <mi>n</mi> <mo>)</mo> </mrow> <mo>&amp;prime;</mo> </msup> <mo>=</mo> <mover> <mi>P</mi> <mo>^</mo> </mover> <mrow> <mo>(</mo> <mi>p</mi> <mo>,</mo> <mi>n</mi> <mo>)</mo> </mrow> <mo>-</mo> <mover> <mi>P</mi> <mo>&amp;OverBar;</mo> </mover> <mrow> <mo>(</mo> <mi>p</mi> <mo>,</mo> <mi>n</mi> <mo>)</mo> </mrow> <mo>.</mo> </mrow>
7.根据权利要求6所述的方法，其特征在于，求解最优化问题时，使用了最速梯度下降法，f(I(p,n))的梯度的表达式为：
<mrow> <mtable> <mtr> <mtd> <mrow> <mfrac> <mrow> <mo>&amp;part;</mo> <mi>f</mi> <mrow> <mo>(</mo> <mrow> <mi>I</mi> <mrow> <mo>(</mo> <mrow> <mi>p</mi> <mo>,</mo> <mi>n</mi> </mrow> <mo>)</mo> </mrow> </mrow> <mo>)</mo> </mrow> </mrow> <mrow> <mo>&amp;part;</mo> <mi>I</mi> <mrow> <mo>(</mo> <mrow> <mi>p</mi> <mo>,</mo> <mi>n</mi> </mrow> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>=</mo> <mfrac> <mrow> <mrow> <mo>(</mo> <mrow> <mi>I</mi> <mrow> <mo>(</mo> <mrow> <mi>p</mi> <mo>,</mo> <mi>n</mi> </mrow> <mo>)</mo> </mrow> <mo>-</mo> <mi>I</mi> <mrow> <mo>(</mo> <mrow> <mi>p</mi> <mo>-</mo> <mn>1</mn> <mo>,</mo> <mi>n</mi> </mrow> <mo>)</mo> </mrow> </mrow> <mo>)</mo> </mrow> <mo>+</mo> <mrow> <mo>(</mo> <mrow> <mi>I</mi> <mrow> <mo>(</mo> <mrow> <mi>p</mi> <mo>,</mo> <mi>n</mi> </mrow> <mo>)</mo> </mrow> <mo>-</mo> <mi>I</mi> <mrow> <mo>(</mo> <mrow> <mi>p</mi> <mo>,</mo> <mi>n</mi> <mo>-</mo> <mn>1</mn> </mrow> <mo>)</mo> </mrow> </mrow> <mo>)</mo> </mrow> </mrow> <msqrt> <mrow> <msup> <mrow> <mo>(</mo> <mrow> <mi>I</mi> <mrow> <mo>(</mo> <mrow> <mi>p</mi> <mo>,</mo> <mi>n</mi> </mrow> <mo>)</mo> </mrow> <mo>-</mo> <mi>I</mi> <mrow> <mo>(</mo> <mrow> <mi>p</mi> <mo>-</mo> <mn>1</mn> <mo>,</mo> <mi>n</mi> </mrow> <mo>)</mo> </mrow> </mrow> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>+</mo> <msup> <mrow> <mo>(</mo> <mrow> <mi>I</mi> <mrow> <mo>(</mo> <mrow> <mi>p</mi> <mo>,</mo> <mi>n</mi> </mrow> <mo>)</mo> </mrow> <mo>-</mo> <mi>I</mi> <mrow> <mo>(</mo> <mrow> <mi>p</mi> <mo>,</mo> <mi>n</mi> <mo>-</mo> <mn>1</mn> </mrow> <mo>)</mo> </mrow> </mrow> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow> </msqrt> </mfrac> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mo>+</mo> <mfrac> <mrow> <mo>(</mo> <mrow> <mi>I</mi> <mrow> <mo>(</mo> <mrow> <mi>p</mi> <mo>,</mo> <mi>n</mi> </mrow> <mo>)</mo> </mrow> <mo>-</mo> <mi>I</mi> <mrow> <mo>(</mo> <mrow> <mi>p</mi> <mo>+</mo> <mn>1</mn> <mo>,</mo> <mi>n</mi> </mrow> <mo>)</mo> </mrow> </mrow> <mo>)</mo> </mrow> <msqrt> <mrow> <msup> <mrow> <mo>(</mo> <mrow> <mi>I</mi> <mrow> <mo>(</mo> <mrow> <mi>p</mi> <mo>+</mo> <mn>1</mn> <mo>,</mo> <mi>n</mi> </mrow> <mo>)</mo> </mrow> <mo>-</mo> <mi>I</mi> <mrow> <mo>(</mo> <mrow> <mi>p</mi> <mo>,</mo> <mi>n</mi> </mrow> <mo>)</mo> </mrow> </mrow> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>+</mo> <msup> <mrow> <mo>(</mo> <mrow> <mi>I</mi> <mrow> <mo>(</mo> <mrow> <mi>p</mi> <mo>+</mo> <mn>1</mn> <mo>,</mo> <mi>n</mi> </mrow> <mo>)</mo> </mrow> <mo>-</mo> <mi>I</mi> <mrow> <mo>(</mo> <mrow> <mi>p</mi> <mo>+</mo> <mn>1</mn> <mo>,</mo> <mi>n</mi> <mo>-</mo> <mn>1</mn> </mrow> <mo>)</mo> </mrow> </mrow> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow> </msqrt> </mfrac> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mo>+</mo> <mfrac> <mrow> <mo>(</mo> <mrow> <mi>I</mi> <mrow> <mo>(</mo> <mrow> <mi>p</mi> <mo>,</mo> <mi>n</mi> </mrow> <mo>)</mo> </mrow> <mo>-</mo> <mi>I</mi> <mrow> <mo>(</mo> <mrow> <mi>p</mi> <mo>,</mo> <mi>n</mi> <mo>+</mo> <mn>1</mn> </mrow> <mo>)</mo> </mrow> </mrow> <mo>)</mo> </mrow> <msqrt> <mrow> <msup> <mrow> <mo>(</mo> <mrow> <mi>I</mi> <mrow> <mo>(</mo> <mrow> <mi>p</mi> <mo>,</mo> <mi>n</mi> <mo>+</mo> <mn>1</mn> </mrow> <mo>)</mo> </mrow> <mo>-</mo> <mi>I</mi> <mrow> <mo>(</mo> <mrow> <mi>p</mi> <mo>,</mo> <mi>n</mi> </mrow> <mo>)</mo> </mrow> </mrow> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>+</mo> <msup> <mrow> <mo>(</mo> <mrow> <mi>I</mi> <mrow> <mo>(</mo> <mrow> <mi>p</mi> <mo>,</mo> <mi>n</mi> <mo>+</mo> <mn>1</mn> </mrow> <mo>)</mo> </mrow> <mo>-</mo> <mi>I</mi> <mrow> <mo>(</mo> <mrow> <mi>p</mi> <mo>-</mo> <mn>1</mn> <mo>,</mo> <mi>n</mi> <mo>+</mo> <mn>1</mn> </mrow> <mo>)</mo> </mrow> </mrow> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow> </msqrt> </mfrac> <mo>+</mo> <mn>2</mn> <mi>&amp;beta;</mi> <mrow> <mo>(</mo> <mrow> <mi>I</mi> <mrow> <mo>(</mo> <mrow> <mi>p</mi> <mo>,</mo> <mi>n</mi> </mrow> <mo>)</mo> </mrow> <mo>-</mo> <mover> <mi>P</mi> <mo>^</mo> </mover> <mrow> <mo>(</mo> <mrow> <mi>p</mi> <mo>,</mo> <mi>n</mi> </mrow> <mo>)</mo> </mrow> </mrow> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> </mtable> <mo>.</mo> </mrow>
8.根据权利要求6所述的方法，其特征在于，所述根据估计出的偏差补偿原始的投影正弦图包括：
将获得的偏差值进行1：8的线性插值，并叠加回原始的投影正弦图。
CN201510293904.2A 2015-06-01 2015-06-01 一种ct探测器像素响应不一致性的校正方法 Active CN104867157B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201510293904.2A CN104867157B (zh) 2015-06-01 2015-06-01 一种ct探测器像素响应不一致性的校正方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201510293904.2A CN104867157B (zh) 2015-06-01 2015-06-01 一种ct探测器像素响应不一致性的校正方法

Publications (2)

Publication Number Publication Date
CN104867157A CN104867157A (zh) 2015-08-26
CN104867157B true CN104867157B (zh) 2018-02-02

Family

ID=53912972

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201510293904.2A Active CN104867157B (zh) 2015-06-01 2015-06-01 一种ct探测器像素响应不一致性的校正方法

Country Status (1)

Country Link
CN (1) CN104867157B (zh)

Families Citing this family (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2017066248A1 (en) * 2015-10-16 2017-04-20 Varian Medical Systems, Inc. Iterative image reconstruction in image-guided radiation therapy
CN106910164A (zh) * 2015-12-23 2017-06-30 通用电气公司 一种对ct投影数据进行滤波的方法及装置
CN105738073B (zh) * 2016-02-03 2019-02-26 中国科学院国家空间科学中心 一种在空间频率域进行像素响应函数测量的方法
CN105841925B (zh) * 2016-03-22 2019-02-26 中国科学院国家空间科学中心 一种基于探测器像素响应傅里叶谱获取的图像重建方法
CN109754446B (zh) * 2018-12-11 2023-01-24 北京纳米维景科技有限公司 一种探测器模块之间的拼接缝宽度估计方法及系统
CN109697476B (zh) * 2019-02-01 2023-06-23 重庆大学 一种基于深度学习的x射线光子计数探测器一致性校准方法
CN117310789B (zh) * 2023-11-30 2024-03-15 赛诺威盛科技(北京)股份有限公司 探测器通道响应线性校正方法、装置、设备和存储介质

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5721427A (en) * 1996-12-19 1998-02-24 Hughes Electronics Scene-based nonuniformity correction processor incorporating motion triggering
CN102300057B (zh) * 2011-06-14 2013-05-01 北京空间机电研究所 线阵ccd像元响应不一致性校正方法
CN104156917A (zh) * 2014-07-30 2014-11-19 天津大学 基于双能谱的x射线ct图像增强方法

Also Published As

Publication number Publication date
CN104867157A (zh) 2015-08-26

Similar Documents

Publication Publication Date Title
CN104867157B (zh) 一种ct探测器像素响应不一致性的校正方法
US8965144B2 (en) Method and system utilizing parameter-less filter for substantially reducing streak and or noise in computer tomography (CT) images
La Riviere et al. Reduction of noise-induced streak artifacts in X-ray computed tomography through spline-based penalized-likelihood sinogram smoothing
US10304217B2 (en) Method and system for generating image using filtered backprojection with noise weighting and or prior in
Tang et al. Cone beam volume CT image artifacts caused by defective cells in x‐ray flat panel imagers and the artifact removal using a wavelet‐analysis‐based algorithm
Chen et al. Three-dimensional point spread function measurement of cone-beam computed tomography system by iterative edge-blurring algorithm
JP4799053B2 (ja) Ｘ線画像における画像乱れの補償方法およびｘ線装置
US20050249431A1 (en) Method for post- reconstructive correction of images of a computer tomograph
CN101690418A (zh) 促进修正图像重建中的增益波动的方法和系统
US20070140416A1 (en) X-ray attenuation correction method, image generating apparatus, x-ray ct apparatus, and image generating method
JP2013085960A (ja) 画像再構成方法及び画像再構成システム
CN109801343A (zh) 基于重建前后图像的环形伪影校正方法、ct控制系统
Anzengruber et al. Convergence rates for Morozov's discrepancy principle using variational inequalities
US6944258B2 (en) Beam hardening post-processing method and X-ray CT apparatus
JPWO2012173205A1 (ja) 画像のノイズレベルを推定する方法
CN101416485A (zh) 温度伪影校正
WO2019085433A1 (zh) 一种基于全变分的图像复原方法及系统
US20110075907A1 (en) X-ray computed tomography apparatus and image processing apparatus
CN114998399B (zh) 一种异源光学遥感卫星影像立体像对预处理方法
US7676074B2 (en) Systems and methods for filtering data in a medical imaging system
US7668360B2 (en) Image processing method and X-ray CT system
KR20110125696A (ko) 엑스선 씨티에서의 링 아티팩트 제거방법 및 장치
Dekker et al. Stray light in cone beam optical computed tomography: II. Reduction using a convergent light source
Yang et al. Post-processing method for the removal of mixed ring artifacts in CT images
Kharfi et al. Spatial resolution limit study of a CCD camera and scintillator based neutron imaging system according to MTF determination and analysis

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
EXSB Decision made by sipo to initiate substantive examination
SE01 Entry into force of request for substantive examination
CB02 Change of applicant information

Address after: 100016, No. 3, building 9, Jiuxianqiao East Road, Beijing, Chaoyang District

Applicant after: Beijing Wandong Medical Polytron Technologies Inc

Address before: 100016 Jiuxianqiao East Road, Beijing, No. 9 A3

Applicant before: Huarun Wandong Medical Equipment Co., Ltd.

CB02 Change of applicant information
GR01 Patent grant
GR01 Patent grant