CN111110259A - 一种识别滤线栅安装状态并进行栅影抑制的方法 - Google Patents

Abstract

本发明公开了一种识别滤线栅安装状态并进行栅影抑制的方法，所述方法包括：获取DR原始图像；从原始图像中抽取关键数据；根据关键数据识别滤线栅安装状态；对已有滤线栅图像计算滤线栅栅影频率；利用滤线栅栅影频率创建滤波器；根据滤波器进行栅影抑制；输出处理后图像。本发明采用的方法可以有效识别出滤线栅安装状态，便于反馈滤线栅是否正确安装；该发明同时对已进行安装的滤线栅装置所产生的栅影进行抑制，从而提高图像质量。

Description

一种识别滤线栅安装状态并进行栅影抑制的方法

技术领域

本发明涉及X射线数字图像处理技术领域，特别是涉及一种识别滤线栅安装状态并进行栅影抑制的方法。

背景技术

在X射线成像过程中，入射X线光子与人体发生相互作用，并产生散射X线光子。散射光子会给图像增加一层随机的灰度信息，降低影像中组织之间的对比度，增加阅片的难度。

滤线栅装置可以滤除散射线，从而提高X光图像的对比度。然而滤线栅本身也会在X光图像上产生相应的影像，这种影像被称作滤线栅影或滤线栅伪影。滤线栅栅影的存在会影响组织的正常显示，增加医生阅片的难度，影响医生对病症的诊断。

目前普遍使用的滤线栅栅影处理方法大多是在已知滤线栅已正确安装的情况下，且滤线栅的频率和方向已知且固定不可改变的情况下进行的。有些方法在计算栅影频率时，仅考虑了高频部分的栅频，没有考虑栅影频率会出现在2个或2个以上栅频段中。在使用频域处理时，大多数方法采用直接对栅影频率截断、或是矩形窗卷积，该方法会导致变换中带进了高频干扰和泄漏，甚至出现负谱现象。有些本发明是完全不知道滤线栅安装状态的情况下，首先进行滤线栅安装状态的识别，如未安装滤线栅、滤线栅安装成反面，或是滤线栅安装为纵向或是横向。然后对确认滤线栅已正确安装的情况，自动计算滤线栅频率，采用滤波器进行滤线栅抑制。

发明内容

本发明的目的是提供一种未知滤线栅安装状态，且未知滤线栅型号和平板探测器像素大小的滤线栅状态识别和栅影抑制方法，该方法能够有效的识别出当前设备安装滤线栅的状态，且对已安装滤线栅的图像进行栅影频段确定，并进行栅影抑制。

本发明提供了一种识别滤线栅安装状态的方法，包括：识别滤线栅是否已安装、滤线栅正反面是否颠倒、滤线栅的安装方向为纵向还是横向的方法。本发明同时提供了一种对已正确安装滤线栅的情况下，进行栅影抑制的方法。

本发明采用如下技术方案：一种识别滤线栅安装状态并进行栅影抑制的方法，包括以下步骤：

A.获取原始图像；

B.从原始数据中抽取关键行和关键列；

C.计算关键行和关键列的频域曲线；

D.根据频域曲线识别滤线栅安装状态；

E.对安装滤线栅的图像，计算滤线栅栅影频率；

F.利用滤线栅栅影频率创建滤波器；

H.利用滤波器进行栅影抑制。

所述步骤B包括：

设定初始需要计算滤线栅频率所占原始图像行数或列数的比率a、b、c，a<b<c；设定计算图像中行方向滤线栅频率和列方向滤线栅频率的个数；

对于行方向和列方向，分别根据设定的个数，在比率[a,b]和[b,c]中插值，使总的行方向比率个数、列方向比率个数等于设定的个数；

用图像总行数乘以根据插值得到的行方向的各个比率，即得到参与滤线栅频率计算的关键行；

用图像总列数乘以根据插值得到的列方向的各个比率，即得到参与滤线栅频率计算的关键列。

所述步骤C包括：

计算行方向栅频率的长度L_{Row_FFT}，计算公式如下：

其中，W表示图像的宽度，函数log₂(*)表示对2取对数，函数int(*)表示对浮点数向下取整数；

对于参与滤线栅频率计算的关键行，截取所有关键行中间长度为L_{Row_FFT}的像素灰度曲线，并计算其傅里叶变换的频域曲线；所述中间长度为L_{Row_FFT}像素灰度曲线是指从关键行的中点位置分别向该行的左和右取L_{Row_FFT}/2长度，这个范围内的像素点的灰度曲线；

将所有计算的关键行傅里叶频域曲线相加取平均，获得行的频域曲线；

计算列方向栅频率的长度L_{Col_FFT}，计算公式如下：

其中，H表示图像的高度；

对于参与滤线栅频率计算的关键列，截取所有关键列中间长度为L_{Col_FFT}的像素灰度曲线，并计算其傅里叶变换的频域曲线；所述中间长度为L_{Col_FFT}像素灰度曲线是指从关键列的中点位置分别向该列的上和下取L_{Col_FFT}/2长度，这个范围内的像素点的灰度曲线；

将所有计算的关键列傅里叶频域曲线相加取平均，获得列的频域曲线。

所述步骤D包括：

设定滤线栅频率值最大阈值Thd_MaxMag,设定幅值阈值Rat_MaxMag；

对于关键行的频域曲线，在[0.5,1]范围内找到频域幅值最大值V_{W_MaxMag}和该最大值邻域幅值的中间值V_{W_NbMed}，若V_{W_MaxMag}＞Thd_MaxMag并且

则初步认定该图像行方向有栅影存在；

对于关键列的频域曲线，在[0.5,1]范围内找到频域幅值最大值V_{H_MaxMag}和该最大值邻域幅值的中间值V_{H_NbMed}，若V_{H_MaxMag}＞Thd_MaxMag并且

则初步认定该图像列方向有栅影存在；

若初步判定的行方向和列方向均没有栅影，则判定有安装滤线栅，或者滤线栅方向安装失败；

若初步判定行方向和列方向中仅一个方向有栅影，则判定装有滤线栅，且栅影方向同初步判定的栅影方向相同；

若初步判定行方向和列方向都存在栅影，比较行方向频域幅值最大值V_{W_MaxMag}的位置和列方向频域幅值最大值V_{H_MaxMag}的位置，其中较接近频率为1的方向为栅影方向位置；若行方向频域幅值最大位置和列方向频域幅值最大位置相对于频率为1的距离相同，则认为没有安装滤线栅。

步骤E包括以下内容：

根据栅影方向，选择该方向对应的频域曲线，为栅影的频域曲线；

在[0.5,1]频率范围内找到满足

要求的所有频率点；对满足要求的所有频率点，计算某一点与其它依次各点之间的距离，遍历满足要求的所有频率点；若某距离小于阈值，或者有频率点在同一个设定窗口内，则对这些频率进行合并，取频率的平均值；

合并之后的满足

要求的频率点，则为栅影所在的频率，频率个数为N_GridF。

步骤F包括以下内容：

针对每一个栅影所在的频率F_i，设定栅影频域窗口大小Lwin_i，设计Gaussian滤波器如下：

μ_i＝F_i

i＝1,2,......,N_GridF

其中，CutThd_i表示对栅影所在频率F_i前后设定范围频域抑制的阈值；μ_i表示Gaussian函数的位置参数，σ_i表示Gaussian函数的尺度参数；L_FFT表示需要滤波的长度，若行方向有栅影，则L_FFT等于图像行数；若列方向有栅影，则L_FFT等于图像列数；

将所有Gaussian滤波器相乘，获取滤线栅滤波器如下：

步骤G包括以下内容：

若滤线栅栅影方向为行方向，则对图像每一行数据进行频域滤波处理：先对行数据进行一维傅里叶变换，在频域中与滤波器进行相乘，然后通过反傅里叶变换，抑制该行数据中的栅影；

若滤线栅栅影方向为列方向，则对图像每一列数据进行频域滤波处理：先对列数据进行一维傅里叶变换，在频域中与滤波器进行相乘，然后通过反傅里叶变换，抑制该列数据中的栅影。

本发明具有以下优点及有益效果：

1.本发明所述方法通过软件来检测和识别滤线栅安装状态，不需要依靠电子或其他辅助装置进行检测。

2.本发明所述方法可以在不知道滤线栅参数和平板探测器像素参数的情况，很好的计算出栅频率，从而进行栅影抑制。

3.本发明所述方法可以计算出多组栅频，从而在栅影抑制时，更好的抑制栅影。

4.本发明所述方法采用的滤波器，不会引起高频干扰和泄漏，也不会引起负谱情况。

附图说明

图1为滤线栅装置与平板探测器放置的正确位置关系示意图；

图2为滤线栅正反面颠倒的示意图；

图3a为滤线栅横向放置的俯视示意图；

图3b为滤线栅纵向放置的俯视示意图；

图4为本发明的流程图；

图5为从原始数据中抽取关键行和关键列的流程图；

图6为根据初步判定栅影方向计算滤线栅安装状态的分支分析图；

图7为某方向有栅影存在时的频域曲线图；

图8为某方向不存在栅影时的频域曲线图；

图9为存在有滤线栅的临床图像；

图10为采用本发明进行滤线栅状态识别并进行栅影抑制的处理后图像。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面将结合本发明具体实施例及相应的附图对本发明技术方案进行清楚、完整地描述。显然，所描述的实施例仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

下面结合附图及实施例对本发明做进一步的详细说明。

目前DR设备中使用的滤线栅装置普遍采用的聚焦型的滤线栅装置，使用的时候放置在平板探测器上方，尽量贴近平板探测器，放置的相对位置可参考图1所示，滤线栅正面朝上。滤线栅通过高衰减材料(通常是铅条)滤掉经过X线通过人体产生的散射线，从而提高图像质量。由于高衰减材料会遮挡X线穿过物体，因此会在图像上留下明暗交替的周期变化的栅影。若将滤线栅正反面放置颠倒，参考图2所示，滤线栅不仅会滤掉更多非散射的射线，还会在图像上形成大面积的伪影，严重影响图像质量。结合上述特点，本发明首先实施从数字图像上计算频域的方法进行识别滤线栅的安装状态。

对于正确安放滤线栅的状态，也有两种放置方式，可以使水平放置，也可以使垂直放置，参考图3a、图3b所示。获取到的数字图像会有滤线栅栅影存在，分析其频域特征，发现栅影频率会在其放置方向上的高频区域，且在中间频域带也会存在一部分栅影频率，因此本发明针对该情况采用合适的滤波器，对滤线栅栅影进行有效抑制。

本发明所述的一种识别滤线栅安装状态并进行栅影抑制的方法，包括以下步骤：

A.获取原始图像；

B.从原始数据中抽取关键行和关键列；

C.计算关键行和关键列的频域曲线；

D.根据频域曲线识别滤线栅安装状态；

E.对正确安装滤线栅的图像，计算滤线栅栅影频率；

F.利用滤线栅栅影频率创建滤波器；

G.利用滤波器进行栅影抑制；

H.获取处理后图像。

本发明具体的流程图如图4所示：

步骤402是从原始数据中抽取关键行和关键列。对关键行和关键列的提取分析需要分开进行，但两者之间不需要有先后顺序。本实施例中以行为例，介绍关键行的提取方法，方法流程图见图5。

子步骤501：设定初始需要计算滤线栅频率所占行数比率，比率的取值范围为[0,1]。本步骤的实施方案之一为设定三个比率值，分别为0.5,0.25,0.75，但具体实施时不局限于使用该三个比率值；

子步骤502：设定计算图像中行方向滤线栅频率的个数，个数要不小于3；

其中，子步骤501和子步骤502没有先后顺序。

子步骤503：根据设定的计算个数，对子步骤设定的行数比率值，进行插值，保证总的比率个数等于设定的计算频率个数；该步骤在子步骤501的实施方案一的插值方法，即在比率[0.25,0.5]和[0.5,0.75]中插值，采用线性插值来保证。当然，插值方式无需特殊说明，仅要求保证总的比率个数等于设定值即可。

子步骤504：结合图像总行数(即图像的高度，单位为像素)，根据插值计算的比率，用总行数乘以各个比率，即得到参与滤线栅频率计算的关键行；

关键列的计算方式可以等同或者类似于子步骤501、502、503，仅将步骤504中的计算改为，用总列数乘以各个比率，得到参与滤线栅频率计算的关键列。

步骤403：计算关键行列的频域曲线。

对于明暗交替的周期栅影，采用单一行或者单一列计算的频域，栅影频率不会太明显，因此将步骤402计算的多组关键行或者多组关键列进行分别计算，然后求和取平均计算栅影频率，更能凸出栅影频率的特点。

下列以计算行方向的栅影频率为例，进行方法说明。

首先，计算行方向栅频率的长度L_{Row_FFT}，计算公式如下：

其中，W表示图像的宽度，函数log₂(*)表示对2取对数，函数int(*)表示对浮点数向下取整数；计算之后的栅频率的长度小于图像宽度W，且是2的整数次幂。

其次，根据步骤402计算得到的滤线栅频率关键行，截取该行中间长度为L_{Row_FFT}的像素灰度曲线，并计算其傅里叶变换的频域曲线；这里不从改行的左右两侧开始取L_{Row_FFT}长度的原因，是由于图像中四周可能会有限束器成像区域，该部分无栅影信息。

最后，将所有计算的关键行傅里叶频域曲线相加取平均，获得行方向的频域曲线；

同样的方法，将所有计算的关键列傅里叶频域曲线相加取平均，获得列的频域曲线。

步骤404：根据频域曲线识别滤线栅安装状态。

滤线栅安装状态分为四种：未安装滤线栅、正反面颠倒安装滤线栅(如图2所示)，水平方向正确安装滤线栅(如图1和图3a所示)、垂直方向安装滤线栅(如图1和图3b所示)。

该步骤根据步骤403计算得到的行方向频域曲线和列方向频域曲线，识别滤线栅的四种安装状态，具体步骤如下：

首先，设定栅频幅值最大阈值Thd_MaxMag,设定最大幅值与邻域中值的倍数Rat_MaxMag；

其次，对于关键行的频域曲线，在[0.5,1]范围内找到频域幅值最大值V_{W_MaxMag}和该最大值邻域的中值V_{W_NbMed}，若V_{W_MaxMag}＞Thd_MaxMag并且

则初步认定该图像行方向有栅影存在；对于关键列的频域曲线，在[0.5,1]范围内找到频域幅值最大值V_{H_MaxMag}和该最大值邻域的中值V_{H_NbMed}，若V_{H_MaxMag}＞Thd_MaxMag并且

则初步认定该图像列方向有栅影存在；

最后，根据初步认定结果进行栅影方向的最终确认，具体分析见图6。若初步判定的行方向和列方向均没有栅影，则判定该系统没有安装滤线栅，或者滤线栅方向安装正反颠倒；若初步判定行方向和列方向仅有一个方向有栅影，则判定该系统安装有滤线栅，且栅影方向同初步判定的栅影方向相同；若初步判定行方向和列方向都存在栅影，比较行方向幅值最大值V_{W_MaxMag}的位置和列方向幅值最大值V_{H_MaxMag}的位置，较接近频率为1的方向为栅影方向位置。若两位位置相对于频率为1的距离一致，则认为没有安装滤线栅。

步骤405：判断滤线栅是否正确安装。若正确安装，即行方向有栅影或者列方向有栅影两种情况。

步骤406：对已有滤波器图像计算滤线栅伪影频率。对于已经正确安装滤线栅的设备，采集到的图像会看到明暗交替的栅影。本发明通过频域分析，提取栅影所在的最大频率和频段。

首先，根据步骤404确认栅影计算方向，选择该方向对应的频域曲线，为栅影的频域曲线；

其次，由于栅影在频域曲线中峰值不仅为1个，需要在[0.5,1]频率范围内找到满足

要求的所有位置。对满足要求的所有位置相对距离进行判定，若距离较近，则进行频率合并，取频率的平均值。实施方案一可以设定一个窗口大小，若满足要求的位置在同一个窗口内，则进行合并。实施方案二可以设定相对距离阈值，若小于该阈值，则进行位置合并。

最后，位置合并之后的满足要求的位置，则为栅影所在的频率，频率个数为N_GridF。如图7为计算的一个频域曲线，该频域中有两个栅影所在的频率，该频率所在的栅影频段已在图7中用椭圆图形标记出来。图7、图8的横坐标为归一化的频率值，纵坐标为频域赋值。

步骤407：利用滤线栅伪影频率创建滤波器。

针对每一个栅影所在的频率F_i，设定栅影频域窗口大小Lwin_i，设计Gaussian滤波器如下：

μ_i＝F_i

i＝1,2,......,N_GridF

其中，CutThd表示对栅影所在频率F_i附近频域抑制的阈值，取值范围在(0,1]；μ_i表示Gaussian函数的位置参数，σ_i表示Gaussian函数的尺度参数。

将所有Gaussian滤波器相乘，获取进行滤线栅滤波器如下：

Gaussian滤波器的优势在于不会引起高频干扰和泄漏，也不会引起负谱情况。

步骤408：利用滤波器进行伪影抑制。若滤线栅栅影方向为行方向，则对图像每一行数据进行频域滤波处理：一维傅里叶变换，在频域中使用滤波器进行滤波，然后通过反傅里叶变换，抑制该行数据中的栅影。对所有行进行重复频域处理。若滤线栅栅影方向为列方向，则对图像每一列数据进行频域滤波处理：一维傅里叶变换，在频域中使用滤波器进行滤波，然后通过反傅里叶变换，抑制该列数据中的栅影。对所有列进行重复频域处理。

图9存在有滤线栅的临床图像，图10采用本发明进行滤线栅状态识别并进行栅影抑制的处理后图像，可以看出该方法可以很好的对栅影进行抑制。

本发明实施例中所述的方法可以是或者也可以不是物理分开的，可以根据实际的需要选择其中的部分或全部方法来达到本实施例方案的目的。按照本发明的上述步骤，本领域的普通技术人员在不付出创造性劳动的情况下，即可以理解并实施。

Claims (7)

1.一种识别滤线栅安装状态并进行栅影抑制的方法，其特征在于包括以下步骤：

A.获取原始图像；

B.从原始数据中抽取关键行和关键列；

C.计算关键行和关键列的频域曲线；

D.根据频域曲线识别滤线栅安装状态；

E.对安装滤线栅的图像，计算滤线栅栅影频率；

F.利用滤线栅栅影频率创建滤波器；

H.利用滤波器进行栅影抑制。

2.根据权利要求1所述的一种识别滤线栅安装状态并进行栅影抑制的方法，其特征在于，所述步骤B包括：

设定初始需要计算滤线栅频率所占原始图像行数或列数的比率a、b、c，a<b<c；设定计算图像中行方向滤线栅频率和列方向滤线栅频率的个数；

对于行方向和列方向，分别根据设定的个数，在比率[a,b]和[b,c]中插值，使总的行方向比率个数、列方向比率个数等于设定的个数；

用图像总行数乘以根据插值得到的行方向的各个比率，即得到参与滤线栅频率计算的关键行；

用图像总列数乘以根据插值得到的列方向的各个比率，即得到参与滤线栅频率计算的关键列。

3.根据权利要求1所述的一种识别滤线栅安装状态并进行栅影抑制的方法，其特征在于，所述步骤C包括：

计算行方向栅频率的长度L_{Row_FFT}，计算公式如下：

其中，W表示图像的宽度，函数log₂(*)表示对2取对数，函数int(*)表示对浮点数向下取整数；

对于参与滤线栅频率计算的关键行，截取所有关键行中间长度为L_{Row_FFT}的像素灰度曲线，并计算其傅里叶变换的频域曲线；所述中间长度为L_{Row_FFT}像素灰度曲线是指从关键行的中点位置分别向该行的左和右取L_{Row_FFT}/2长度，这个范围内的像素点的灰度曲线；

将所有计算的关键行傅里叶频域曲线相加取平均，获得行的频域曲线；

计算列方向栅频率的长度L_{Col_FFT}，计算公式如下：

其中，H表示图像的高度；

对于参与滤线栅频率计算的关键列，截取所有关键列中间长度为L_{Col_FFT}的像素灰度曲线，并计算其傅里叶变换的频域曲线；所述中间长度为L_{Col_FFT}像素灰度曲线是指从关键列的中点位置分别向该列的上和下取L_{Col_FFT}/2长度，这个范围内的像素点的灰度曲线；

将所有计算的关键列傅里叶频域曲线相加取平均，获得列的频域曲线。

4.根据权利要求1所述的一种识别滤线栅安装状态并进行栅影抑制的方法，其特征在于，所述步骤D包括：

设定滤线栅频率值最大阈值Thd_MaxMag,设定幅值阈值Rat_MaxMag；

对于关键行的频域曲线，在[0.5,1]范围内找到频域幅值最大值V_{W_MaxMag}和该最大值邻域幅值的中间值V_{W_NbMed}，若V_{W_MaxMag}＞Thd_MaxMag并且

则初步认定该图像行方向有栅影存在；

对于关键列的频域曲线，在[0.5,1]范围内找到频域幅值最大值V_{H_MaxMag}和该最大值邻域幅值的中间值V_{H_NbMed}，若V_{H_MaxMag}＞Thd_MaxMag并且

则初步认定该图像列方向有栅影存在；

若初步判定的行方向和列方向均没有栅影，则判定有安装滤线栅，或者滤线栅方向安装失败；

若初步判定行方向和列方向中仅一个方向有栅影，则判定装有滤线栅，且栅影方向同初步判定的栅影方向相同；

若初步判定行方向和列方向都存在栅影，比较行方向频域幅值最大值V_{W_MaxMag}的位置和列方向频域幅值最大值V_{H_MaxMag}的位置，其中较接近频率为1的方向为栅影方向位置；若行方向频域幅值最大位置和列方向频域幅值最大位置相对于频率为1的距离相同，则认为没有安装滤线栅。

5.根据权利要求1所述的一种识别滤线栅安装状态并进行栅影抑制的方法，其特征在于，步骤E包括以下内容：

根据栅影方向，选择该方向对应的频域曲线，为栅影的频域曲线；

在[0.5,1]频率范围内找到满足

要求的所有频率点；对满足要求的所有频率点，计算某一点与其它依次各点之间的距离，遍历满足要求的所有频率点；若某距离小于阈值，或者有频率点在同一个设定窗口内，则对这些频率进行合并，取频率的平均值；

合并之后的满足

要求的频率点，则为栅影所在的频率，频率个数为N_GridF。

6.根据权利要求1所述的一种识别滤线栅安装状态并进行栅影抑制的方法，其特征在于，步骤F包括以下内容：

针对每一个栅影所在的频率F_i，设定栅影频域窗口大小Lwin_i，设计Gaussian滤波器如下：

μ_i＝F_i

其中，CutThd_i表示对栅影所在频率F_i前后设定范围频域抑制的阈值；μ_i表示Gaussian函数的位置参数，σ_i表示Gaussian函数的尺度参数；L_FFT表示需要滤波的长度，若行方向有栅影，则L_FFT等于图像行数；若列方向有栅影，则L_FFT等于图像列数；

将所有Gaussian滤波器相乘，获取滤线栅滤波器如下：

7.根据权利要求1所述的一种识别滤线栅安装状态并进行栅影抑制的方法，其特征在于，步骤G包括以下内容：

若滤线栅栅影方向为行方向，则对图像每一行数据进行频域滤波处理：先对行数据进行一维傅里叶变换，在频域中与滤波器进行相乘，然后通过反傅里叶变换，抑制该行数据中的栅影；

若滤线栅栅影方向为列方向，则对图像每一列数据进行频域滤波处理：先对列数据进行一维傅里叶变换，在频域中与滤波器进行相乘，然后通过反傅里叶变换，抑制该列数据中的栅影。

Images