CN109829902A - 一种基于广义S变换和Teager属性的肺部CT图像结节筛选方法 - Google Patents

Abstract

本发明公开了一种基于广义S变换和Teager属性的肺部CT图像结节筛选方法，属于肺部图像处理领域；其包括步骤1：对输入的肺部CT图像进行水平方向和竖直方向上的广义S变换，得到水平方向时频谱和竖直方向时频谱；步骤2：对所述水平方向时频谱和竖直方向时频谱进行Teager主能量的提取，得到水平方向Teager主能量图和竖直方向Teager主能量图；步骤3：将所述水平方向Teager主能量图和竖直方向Teager主能量图进行阈值分割，得到疑似结节；本发明通过广义S变换和计算Teager主能量属性，从时频分析的角度来分析结节和非结节区域的差异，通过寻找非零像素克服肺部CT图像的肺部边界对广义S变换的影响，利用时频谱中的Teager主能量属性的区别从而筛选出疑似结节区域，提高了筛选准确度。

Description

一种基于广义S变换和Teager属性的肺部CT图像结节筛选 方法

技术领域

本发明属于肺部图像处理领域，特别涉及一种基于广义S变换和Teager属性的肺部CT图像结节筛选方法。

背景技术

广义S变换是一种时频分析方法，对非平稳信号而言，时频分析方法可以分析出信号的频率成分和定位频率成分的位置；广义S变换不仅具有S变换的优点，比短时傅里叶变换有更好的时频分辨率，没有Wigner-Ville交叉项的干扰，而且相较于S变换有更广的可调节的频率分辨率。所以，广义S变换有更好的灵活性和更高的时频分辨率。

Teager主能量属性是基于Teager-Kaiser(TK)能量的改进，可以对信号的局部能量变换进行追踪和提取，相较于TK能量只能进行单频计算的缺点，Teager主能量可以多频段的计算；Terger主能量比其他的能量算子有更好的能量聚焦性。

在肺部CT图像处理中，结节筛选的方法分为四类：一、基于阈值的方法，利用结节和非结节区域的灰度值的差异进行筛选；二、基于滤波的方法，利用结节和非结节区域的形态差异进行筛选；三、基于匹配的方法，利用结节的形态信息进行筛选；四、基于聚类的方法，利用灰度值信息进行筛选；虽然四类方法的原理不同，但都是基于结节在空间的灰度信息和形态信息差异来进行筛选，无法完整地体现结节信息，导致假阳性结节个数多，筛选精度不高。因此需要一种筛选方法，从结节的频率信息进行筛选，提供一种高精度的结节筛选方法。

发明内容

本发明的目的在于：针对现有的肺部CT图像中结节筛选的方法中只利用了结节在空间灰度信息和形态信息，并没有利用结节的频率信息的问题，提供一种基于广义S变换和Teager属性的肺部CT图像结节筛选方法，利用广义S变换获得的时频谱，从时频谱中提取Teager主能量信息，为肺部CT图像结节筛选方法提供了新的筛选方法，即时频分析方法，从而通过Teager能量进行结节筛选。

本发明采用的技术方案如下：

一种基于广义S变换和Teager属性的肺部CT图像结节筛选方法，包括以下步骤：

步骤1：对输入的肺部CT图像进行水平方向和竖直方向上的广义S变换，得到水平方向时频谱和竖直方向时频谱；

步骤2：对所述水平方向时频谱和竖直方向时频谱进行Teager主能量的提取，得到水平方向Teager主能量图和竖直方向Teager主能量图；

步骤3：将所述水平方向Teager主能量图和竖直方向Teager主能量图进行阈值分割，得到疑似结节；

所述步骤1中对输入的肺部CT图像进行预处理，所述预处理包括确定肺部CT图像行变量和列变量的取值范围和去除肺部以外的像素，具体如下包括如下：

步骤a：确定输入的肺部CT图像中行变量和列变量的取值范围即非零像素，并计算非零像素所在的最小行mir、最大行mar、最小列mic和最大列mac，计算公式如下：

{(x，y)|I(x，y)＞0}

mir＝min(x)，mar＝max(x)，mic×min(y)，mac＝max(y)

其中，I(x，y)表示肺部CT图像在点(x，y)处的灰度值，x是水平方向位置变量，y是竖直方向位置变量，min表示最小值，max表示最大值，mir表示最小行数值，mar表示最大行数值，mic表示最小列数值，mac表示最大列数值；

步骤b：初始化水平方向行变量i、竖直方向列变量j和广义S变换的窗口调节参数p；

步骤c：提取CT图像水平方向第i行一维信号和第j列一维信号，去除零像素得到水平方向第i行一维信号sigh和竖直方向第j列一维信号sigv。

优选地，所述步骤b具体为：

初始化广义S变换的窗口调节参数p＝0.1；

计算水平方向时频谱时，初始化水平方向行变量i＝mir，i∈mir～mar；

计算竖直方向时频谱时，初始化竖直方向列变量j＝mic，j∈mic～mac。

优选地，所述步骤c中，计算第i行一维信号sigh和竖直方向第j列一维信号sigv的公式如下：

sigh＝{I(x，：)|I(x，：)＞0，I(x，：)∈I(i，：)}

sigv＝{I(：，y)|I(：，y)＞0，I(：，y)∈I(：，j)}

其中，I(i，：)表示肺部CT图像水平方向的第i行一维信号，I(x，：)表示第i行一维信号在x处的灰度值；I(：，j)表示肺部CT图像竖直方向的第j行一维信号，I(：，y)表示第j列一维信号在y处的灰度值；sigh表示水平方向第i行一维信号，sigv表示竖直方向第j列一维信号。

优选地，所述步骤1包括如下步骤：

步骤1.1：初始化水平方向频率变量和竖直方向频率变量：具体如下：

计算水平方向第i行一维信号广义S变换时，初始化水平方向频率变量f_x＝0，f_x取值范围为0～m，m表示水平方向第i行一维信号sigh的长度；

计算竖直方向第j列一维信号广义S变换时，初始化水平方向频率变量f_y＝0，f_y取值范围为0～n，n表示竖直方向第j列一维信号sigv的长度；

步骤1.2：根据步骤1.1的水平方向频率变量和竖直方向频率变量，计算水平方向第i行一维信号sigh和竖直方向第j列一维信号sigv对应的平移后的频谱图；

步骤1.3：根据步骤1.1的水平方向频率变量和竖直方向频率变量，计算频率域高斯窗函数，计算公式如下：

其中，G(f)表示频率域高斯窗函数，f表示水平方向频率变量或者竖直方向频率变量，α表示频率变量，p为窗口调节参数；

步骤1.4：计算步骤1.2所述平移后的频谱图和步骤1.3所述频率域高斯窗函数的乘积；

步骤1.5：对步骤1.4的结果进行反傅里叶变换，分别得到水平方向时频谱和竖直方向时频谱，反傅里叶变换的公式如下：

其中，M(f)表示步骤1.4的乘积结果，f表示水平方向频率变量或者竖直方向频率变量，α表示已知的频率变量。

优选地，所述步骤2包括如下步骤：

步骤2.1：初始化位置参数，具体地：

水平方向初始化时，变量x＝2，x∈2～m-1，m表示水平方向第i行一维信号sigh的长度；

竖直方向初始化时，变量y＝2，y∈2～n-1，n表示竖直方向第j列一维信号sigv的长度；

步骤2.2：根据步骤1.5求得在频率为f时的时频谱S(t，f，p)，计算出TK能量，其计算公式如下：

retk(t，f)＝re(S(t，f，p))²-re(S(t-1，f，p))*re(S(t+1，f，p))

imtk(t，f)＝im(S(t，f，p))²-im(S(t-1，f，p))*im(S(t+1，f，p))

tk(t，f)＝retk(t，f)+imtk(t，f)

其中，re()表示取实部，im()表示取虚部，retk(t，f)表示实部TK能量，imtk(t，f)表示虚部TK能量，tk(t，f)表示TK能量；

步骤2.3：执行步骤2.2后，判断水平方向位置变量和竖直方向位置变量的值是否循环遍历了所有值，若否，则将水平方向位置变量和竖直方向位置变量自增后，继续执行步骤2.2，否则跳出循环后跳至步骤2.4；

步骤2.4：判断水平方向频率变量和竖直方向频率变量的值是否循环遍历了所有值，若否，则将水平方向频率变量和竖直方向频率变量自增后，继续执行步骤1.5到步骤2.3，否则跳出循环后跳至步骤2.5；

步骤2.5：在循环遍历所有频率数值后，计算Teager主能量，其计算公式为：

其中，tm(t，f)表示Teager主能量，tk(t，f)表示步骤2.2求得的TK能量，f表示频率变量；

步骤2.6：执行步骤2.5后，判断水平方向行变量和竖直方向列变量的值是否循环遍历了所有值，若否，则将水平方向行变量和竖直方向列变量自增后，继续执行步骤1.3到步骤2.5，否则跳出循环后跳至步骤2.7；

步骤2.7：循环遍历所有水平方向行变量和竖直方向列变量的值后，得到水平方向Teager主能量图和竖直方向Teager主能量图，计算公式如下：

tm_h(x，y)＝{tm(i，f)，i∈[mir，mar]}

tm_v(x，y)＝{tm(j，f)，j∈[mic，mac]}

其中，tm_h(x，y)表示水平方向Teager主能量图，tm_v(x，y)表示竖直方向Teager主能量图，tm(i，f)表示步骤2.5所求的第i行一维信号sigh的Teager主能量，tm(j，f)表示步骤2.5所求的第j列一维信号sigv的Teager主能量，mir表示最小行数值，mar表示最大行数值，mic表示最小列数值，mac表示最大列数值。

优选地，步骤3具体为：

步骤3.1：对步骤2求得的水平方向Teager主能量图和竖直方向Teager主能量图进行归一化、阈值分割，分别得到水平方向二值化图像NB_x(x，y)和竖直方向二值化图像NB_y(x，y)，阈值T＝0.7，归一化的公式如下：

ntm(x，y)＝(tm(x，y)-min(tm(x，y)))/(max(tm(x，y))-min(tm(x，y))))

其中，ntm(x，y)表示归一化的Teager主能量，tm(x，y)为步骤2所求的Teager主能量，max()表示最大值，min()表示最小值。

步骤3.2：对所述水平方向二值化图像NB_x(x，y)和竖直方向二值化图像NB_y(x，y)进行与操作，筛选出疑似结节。

综上所述，由于采用了上述技术方案，本发明的有益效果是：

1.本发明基于广义S变换，对得到的时频谱进行Teager主能量的提取，然后进行阈值分割，完成对肺部CT图像结节的筛选，通过寻找非零像素克服肺部CT图像的肺部边界对广义S变换的影响，为肺部CT图像结节筛选方法提供了新的筛选方法即时频分析方法，从而通过Teager能量进行结节筛选，降低假阳性结节个数，利于提高结节筛选准确度；

2.本发明为克服肺部CT图像的肺部边界对广义S变换的影响，通过寻找非零像素，去除肺部边界，从而避免了边界异常的变换所带来的异常区域，为后续的操作减少了误判的可能性，进一步利于提高结节筛选准确度；

3.本发明通过广义S变换，将原始的空间信息变换到时频域得到时频谱，为结节提供了一种全新的时频信息；广义S变换相较于其他时频分析方法，时频聚焦性高，无交差项的干扰，获得的频率域信息有更好的区分度；

4.本发明通过在时频域计算Teager主能量，为结节提供了频率域能量信息，从而实现利用时频分析筛选结节，Teager主能量相较于其他能量算子，对信号有更好的能量聚焦性和追踪性；

5.本发明阈值处理阶段的归一化和固定阈值处理对临床复杂的肺部图像有较好的鲁棒性，避免了不同CT图像之间Teager能量变换范围之间的差异而引起的阈值的变换，利于提高结节筛选准确度。

附图说明

为了更清楚地说明本发明实施例的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，应当理解，以下附图仅示出了本发明的某些实施例，因此不应被看作是对范围的限定，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他相关的附图。

图1为本发明的流程图；

图2为本发明输入的肺部CT图像；

图3为本发明的水平方向Teager主能量图；

图4为本发明的竖直方向Teager主能量图；

图5为本发明的水平方向二值化图像；

图6为本发明的竖直方向二值化图像；

图7为本发明最终的处理图像。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明，即所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本发明实施例的组件可以以各种不同的配置来布置和设计。

因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明的实施例，本领域技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。

需要说明的是，术语“第一”和“第二”等之类的关系术语仅仅用来将一个实体或者操作与另一个实体或操作区分开来，而不一定要求或者暗示这些实体或操作之间存在任何这种实际的关系或者顺序。而且，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者设备所固有的要素。在没有更多限制的情况下，由语句“包括一个……”限定的要素，并不排除在包括所述要素的过程、方法、物品或者设备中还存在另外的相同要素。

技术问题：解决现有技术基于结节在空间的灰度信息和形态信息差异来进行筛选，无法完整地体现结节信息，导致假阳性结节个数多，筛选精度不高的问题；

技术手段：一种基于广义S变换和Teager属性的肺部CT图像结节筛选方法，包括以下步骤：

步骤1：对输入的肺部CT图像进行水平方向和竖直方向上的广义S变换，得到水平方向时频谱和竖直方向时频谱；

步骤2：对所述水平方向时频谱和竖直方向时频谱进行Teager主能量的提取，得到水平方向Teager主能量图和竖直方向Teager主能量图；

步骤3：将所述水平方向Teager主能量图和竖直方向Teager主能量图进行阈值分割，得到疑似结节；

步骤1中对输入的肺部CT图像进行预处理，所述预处理包括确定肺部CT图像行变量和列变量的取值范围和去除肺部以外的像素，具体如下包括如下：

步骤a：确定输入的肺部CT图像中行变量和列变量的取值范围即非零像素，并计算非零像素所在的最小行mir、最大行mar、最小列mic和最大列mac，计算公式如下：

{(x，y)|I(x，y)＞0}

mir＝min(x)，mar＝max(x)，mic＝min(y)，mac＝max(y)

其中，I(x，y)表示肺部CT图像在点(x，y)处的灰度值，x是水平方向位置变量，y是竖直方向位置变量，min表示最小值，max表示最大值，mir表示最小行数值，mar表示最大行数值，mic表示最小列数值，mac表示最大列数值；

步骤b：初始化水平方向行变量i、竖直方向列变量j和广义S变换的窗口调节参数p；

步骤c：提取CT图像水平方向第i行一维信号和第j列一维信号，去除零像素得到水平方向第i行一维信号sigh和竖直方向第j列一维信号sigv。

步骤b具体为：

初始化广义S变换的窗口调节参数p＝0.1；

计算水平方向时频谱时，初始化水平方向行变量i＝mir，i∈mir～mar；

计算竖直方向时频谱时，初始化竖直方向列变量j＝mic，j∈mic～mac。

步骤c中，计算第i行一维信号sigh和竖直方向第j列一维信号sigv的公式如下：

sigh＝{I(x，:)|I(x，:)＞0，I(x，:)∈I(i，:)}

sigv＝{I(:，y)|I(:，y)＞0，I(:，y)∈I(:，j)}

其中，I(i，：)表示肺部CT图像水平方向的第i行一维信号，I(x，：)表示第i行一维信号在处的灰度值；I(：，j)表示肺部CT图像竖直方向的第j行一维信号，I(：，y)表示第j列一维信号在y处的灰度值；sigh表示水平方向第i行一维信号，sigv表示竖直方向第j列一维信号。

步骤1包括如下步骤：

步骤1.1：初始化水平方向频率变量和竖直方向频率变量：具体如下：

计算水平方向第i行一维信号广义S变换时，初始化水平方向频率变量f_x＝0，f_x取值范围为0～m，m表示水平方向第i行一维信号sigh的长度；

计算竖直方向第j列一维信号广义S变换时，初始化水平方向频率变量f_y＝0，f_y取值范围为0～n，n表示竖直方向第j列一维信号sigv的长度；

步骤1.2：根据步骤1.1的水平方向频率变量和竖直方向频率变量，计算水平方向第i行一维信号sigh和竖直方向第j列一维信号sigv对应的平移后的频谱图；

步骤1.3：根据步骤1.1的水平方向频率变量和竖直方向频率变量，计算频率域高斯窗函数，计算公式如下：

其中，G(f)表示频率域高斯窗函数，f表示水平方向频率变量或者竖直方向频率变量，α表示频率变量，p为窗口调节参数；

步骤1.4：计算步骤1.2所述平移后的频谱图和步骤1.3所述频率域高斯窗函数的乘积；

步骤1.5：对步骤1.4的结果进行反傅里叶变换，分别得到水平方向时频谱和竖直方向时频谱，反傅里叶变换的公式如下：

其中，M(f)表示步骤1.4的乘积结果，f表示水平方向频率变量或者竖直方向频率变量，α表示已知的频率变量。

步骤2包括如下步骤：

步骤2.1：初始化位置参数，具体地：

水平方向初始化时，变量x＝2，x∈2～m-1，m表示水平方向第i行一维信号sigh的长度；

竖直方向初始化时，变量y＝2，y∈2～n-1，n表示竖直方向第j列一维信号sigv的长度；

步骤2.2：根据步骤1.5求得在频率为f时的时频谱S(t，f，p)，计算出TK能量，其计算公式如下：

retk(t，f)＝re(S(t，f，p))²-re(S(t-1，f，p))*re(S(t+1，f，p))

imtk(t，f)＝im(S(t，f，p))²-im(S(t-1，f，p))*im(S(t+1，f，p))

tk(t，f)＝retk(t，f)+imtk(t，f)

其中，re()表示取实部，im()表示取虚部，retk(t，f)表示实部TK能量，imtk(t，f)表示虚部TK能量，tk(t，f)表示TK能量；

步骤2.3：执行步骤2.2后，判断水平方向位置变量和竖直方向位置变量的值是否循环遍历了所有值，若否，则将水平方向位置变量和竖直方向位置变量自增后，继续执行步骤2.2，否则跳出循环后跳至步骤2.4；

步骤2.4：判断水平方向频率变量和竖直方向频率变量的值是否循环遍历了所有值，若否，则将水平方向频率变量和竖直方向频率变量自增后，继续执行步骤1.5到步骤2.3，否则跳出循环后跳至步骤2.5；

步骤2.5：在循环遍历所有频率数值后，计算Teager主能量，其计算公式为：

其中，tm(t，f)表示Teager主能量，tk(t，f)表示步骤2.2求得的TK能量，f表示频率变量；

步骤2.6：执行步骤2.5后，判断水平方向行变量和竖直方向列变量的值是否循环遍历了所有值，若否，则将水平方向行变量和竖直方向列变量自增后，继续执行步骤1.3到步骤2.5，否则跳出循环后跳至步骤2.7；

步骤2.7：循环遍历所有水平方向行变量和竖直方向列变量的值后，得到水平方向Teager主能量图和竖直方向Teager主能量图，计算公式如下：

tm_h(x，y)＝{tm(i，f)，i∈[mir，mar]}

tm_v(x，y)＝{tm(j，f)，j∈[mic，mac]}

其中，tm_h(x，y)表示水平方向Teager主能量图，tm_v(x，y)表示竖直方向Teager主能量图，tm(i，f)表示步骤2.5所求的第i行一维信号sigh的Teager主能量，tm(j，f)表示步骤2.5所求的第j列一维信号sigv的Teager主能量，mir表示最小行数值，mar表示最大行数值，mic表示最小列数值，mac表示最大列数值。

步骤3具体为：

步骤3.1：对步骤2求得的水平方向Teager主能量图和竖直方向Teager主能量图进行归一化、阈值分割，分别得到水平方向二值化图像NB_x(x，y)和竖直方向二值化图像NB_y(x，y)，阈值T＝0.7，归一化的公式如下：

ntm(x，y)＝(tm(x，y)-min(tm(x，y)))/(max(tm(x，y))-min(tm(x，y))))

其中，ntm(x，y)表示归一化的Teager主能量，tm(x，y)为步骤2所求的Teager主能量，max()表示最大值，min()表示最小值。

步骤3.2：对所述水平方向二值化图像NB_x(x，y)和竖直方向二值化图像NB_y(x，y)进行与操作，筛选出疑似结节。

技术效果：本发明基于广义S变换，对得到的时频谱进行Teager主能量的提取，然后进行阈值分割，完成对肺部CT图像结节的筛选，通过寻找非零像素克服肺部CT图像的肺部边界对广义S变换的影响，为肺部CT图像结节筛选方法提供了新的筛选方法即时频分析方法，从而通过Teager能量进行结节筛选，降低假阳性结节个数，利于提高结节筛选准确度；为克服肺部CT图像的肺部边界对广义S变换的影响，通过寻找非零像素，去除肺部边界，从而避免了边界异常的变换所带来的异常区域，为后续的操作减少了误判的可能性，进一步利于提高结节筛选准确度；通过广义S变换，将原始的空间信息变换到时频域得到时频谱，为结节提供了一种全新的时频信息；广义S变换相较于其他时频分析方法，时频聚焦性高，无交差项的干扰，获得的频率域信息有更好的区分度；通过在时频域计算Teager主能量，为结节提供了频率域能量信息，从而实现利用时频分析筛选结节，Teager主能量相较于其他能量算子，对信号有更好的能量聚焦性和追踪性；阈值处理阶段的归一化和固定阈值处理对临床复杂的肺部图像有较好的鲁棒性，避免了不同CT图像之间Teager能量变换范围之间的差异而引起的阈值的变换，利于提高结节筛选准确度。

以下结合实施例对本发明的特征和性能作进一步的详细描述。

实施例1

如图1所示，一种基于广义S变换和Teager属性的肺部CT图像结节筛选方法，包括以下步骤：

步骤1：对输入的肺部CT图像进行水平方向和竖直方向上的广义S变换，得到水平方向时频谱和竖直方向时频谱；

步骤2：对所述水平方向时频谱和竖直方向时频谱进行Teager主能量的提取，得到水平方向Teager主能量图和竖直方向Teager主能量图；

步骤3：将所述水平方向Teager主能量图和竖直方向Teager主能量图进行阈值分割，得到疑似结节；

步骤1中对输入的肺部CT图像进行预处理，所述预处理包括确定肺部CT图像行变量和列变量的取值范围和去除肺部以外的像素，具体如下包括如下：

步骤a：确定输入的肺部CT图像中行变量和列变量的取值范围即非零像素，并计算非零像素所在的最小行mir、最大行mar、最小列mic和最大列mac，计算公式如下：

{(x，y)|I(x，y)＞0}

mir＝min(x)，mar＝max(x)，mic＝min(y)，mac＝max(y)

其中，I(x，y)表示肺部CT图像在点(x，y)处的灰度值，x是水平方向位置变量，y是竖直方向位置变量，min表示最小值，max表示最大值，mir表示最小行数值，mar表示最大行数值，mic表示最小列数值，mac表示最大列数值；

步骤b：初始化水平方向行变量i、竖直方向列变量j和广义S变换的窗口调节参数p；

步骤c：提取CT图像水平方向第i行一维信号和第j列一维信号，去除零像素得到水平方向第i行一维信号sigh和竖直方向第j列一维信号sigv。

步骤b具体为：

初始化广义S变换的窗口调节参数p＝0.1；

计算水平方向时频谱时，初始化水平方向行变量i＝mir，i∈mir～mar；

计算竖直方向时频谱时，初始化竖直方向列变量j＝mic，j∈mic～mac。

步骤c中，计算第i行一维信号sigh和竖直方向第j列一维信号sigv的公式如下：

sigh＝{I(x，:)|I(x，:)＞0，I(x，:)∈I(i，:)}

sigv＝{I(:，y)|I(:，y)＞0，I(:，y)∈I(:，j)}

其中，I(i，：)表示肺部CT图像水平方向的第i行一维信号，I(x，：)表示第i行一维信号在处的灰度值；I(：，j)表示肺部CT图像竖直方向的第j行一维信号，I(：，y)表示第j列一维信号在y处的灰度值；sigh表示水平方向第i行一维信号，sigv表示竖直方向第j列一维信号。

步骤3具体为：

步骤3.1：对步骤2求得的水平方向Teager主能量图和竖直方向Teager主能量图进行归一化、阈值分割，分别得到水平方向二值化图像NB_x(x，y)和竖直方向二值化图像NB_y(x，y)，阈值T＝0.7，归一化的公式如下：

ntm(x，y)＝(tm(x，y)-min(tm(x，y)))/(max(tm(x，y))-min(tm(x，y))))

其中，ntm(x，y)表示归一化的Teager主能量，tm(x，y)为步骤2所求的Teager主能量，max()表示最大值，min()表示最小值；

阈值分割的公式如下：

步骤3.2：对所述水平方向二值化图像NB_x(x，y)和竖直方向二值化图像NB_y(x，y)进行与操作，筛选出疑似结节。阈值处理阶段的归一化和固定阈值处理对临床复杂的肺部图像有较好的鲁棒性，避免了不同CT图像之间Teager能量变换范围之间的差异而引起的阈值的变换，利于提高结节筛选准确度；

克服的技术问题：肺部CT图像的肺部边界对广义S变换有影响，通过寻找非零像素，非零像素确定行和列变量的取值范围，去除肺部边界，从而避免了边界异常的变换所带来的异常区域，为后续的操作减少了误判的可能性，进一步利于提高结节筛选准确度。

如图2所示，为输入的肺部CT图像，从图中可看出有部分大面积高亮区域为结节区域，通过本申请进行筛选获取如图7所示的筛选结果图像，本申请具体处理如图3-6所示，图7所示的筛选结果图像不仅突出了结节区域，还将结节区域的非影响因素进行了抑制，提供了新的筛选方法即时频分析方法，提高了结节筛选准确度；综上，本发明基于广义S变换，对得到的时频谱进行Teager主能量的提取，然后进行阈值分割，完成对肺部CT图像结节的筛选，通过寻找非零像素克服肺部CT图像的肺部边界对广义S变换的影响，为肺部CT图像结节筛选方法提供了新的筛选方法即时频分析方法，从而通过Teager能量进行结节筛选，降低假阳性结节个数，利于提高结节筛选准确度；为克服肺部CT图像的肺部边界对广义S变换的影响，通过寻找非零像素，去除肺部边界，从而避免了边界异常的变换所带来的异常区域，如图3-4所示，肺部边界无异常值(数值较大值)，为后续的操作减少了误判的可能性，进一步利于提高结节筛选准确度。

实施例2

基于实施例1，步骤1：对输入的肺部CT图像进行水平方向和竖直方向上的广义S变换，得到水平方向时频谱和竖直方向时频谱，包括如下步骤：

步骤1.1：初始化水平方向频率变量和竖直方向频率变量：具体如下：

计算水平方向第i行一维信号广义S变换时，初始化水平方向频率变量f_x＝0，f_x取值范围为0～m，m表示水平方向第i行一维信号sigh的长度；

计算竖直方向第j列一维信号广义S变换时，初始化水平方向频率变量f_y＝0，f_y取值范围为0～n，n表示竖直方向第j列一维信号sigv的长度；

步骤1.2：根据步骤1.1的水平方向频率变量和竖直方向频率变量，计算水平方向第i行一维信号sigh和竖直方向第j列一维信号sigv对应的平移后的频谱图，傅里叶频谱计算的公式如下：

其中，H(α+f)表示得到的平移后的频谱图，h(x)表示水平方向第i行一维信号和竖直方向第j列一维信号，f表示水平方向或者竖直方向的频率数值，α表示频率变量；

步骤1.3：根据步骤1.1的水平方向频率变量和竖直方向频率变量，计算频率域高斯窗函数，计算公式如下：

其中，G(f)表示频率域高斯窗函数，f表示水平方向频率变量或者竖直方向频率变量，α表示频率变量，p为窗口调节参数；

步骤1.4：计算步骤1.2所述平移后的频谱图和步骤1.3所述频率域高斯窗函数的乘积；

步骤1.5：对步骤1.4的结果进行反傅里叶变换，分别得到水平方向时频谱和竖直方向时频谱，反傅里叶变换的公式如下：

其中，M(f)表示步骤1.4的乘积结果，f表示水平方向频率变量或者竖直方向频率变量，α表示已知的频率变量。

通过广义S变换，将CT图像从空域变换到时频域，从而获得结节的时频信息，为结节的筛选提供了一种新颖的频率信息。

实施例3

基于实施例1，步骤2：对所述水平方向时频谱和竖直方向时频谱进行Teager主能量的提取，得到水平方向Teager主能量图和竖直方向Teager主能量图，包括如下步骤：

步骤2.1：初始化位置参数，具体地：

水平方向初始化时，变量x＝2，x∈2～m-1，m表示水平方向第i行一维信号sigh的长度；

竖直方向初始化时，变量y＝2，y∈2～n-1，n表示竖直方向第j列一维信号sigv的长度；

步骤2.2：根据步骤1.5求得在频率为f时的时频谱S(t，f，p)，计算出TK能量，其计算公式如下：

retk(t，f)＝re(S(t，f，p))²-re(S(t-1，f，p))*re(S(t+1，f，p))

imtk(t，f)＝im(S(t，f，p))²-im(S(t-1，f，p))*im(S(t+1，f，p))

tk(t，f)＝retk(t，f)+imtk(t，f)

其中，re()表示取实部，im()表示取虚部，retk(t，f)表示实部TK能量，imtk(t，f)表示虚部TK能量，tk(t，f)表示TK能量；

步骤2.3：执行步骤2.2后，判断水平方向位置变量和竖直方向位置变量的值是否循环遍历了所有值，若否，则将水平方向位置变量和竖直方向位置变量自增后，继续执行步骤2.2，否则跳出循环后跳至步骤2.4；

步骤2.4：判断水平方向频率变量和竖直方向频率变量的值是否循环遍历了所有值，若否，则将水平方向频率变量和竖直方向频率变量自增后，继续执行步骤1.5到步骤2.3，否则跳出循环后跳至步骤2.5；

步骤2.5：在循环遍历所有频率数值后，计算Teager主能量，其计算公式为：

其中，tm(t，f)表示Teager主能量，tk(t，f)表示步骤2.2求得的TK能量，f表示频率变量；

步骤2.6：执行步骤2.5后，判断水平方向行变量和竖直方向列变量的值是否循环遍历了所有值，若否，则将水平方向行变量和竖直方向列变量自增后，继续执行步骤1.3到步骤2.5，否则跳出循环后跳至步骤2.7；

步骤2.7：循环遍历所有水平方向行变量和竖直方向列变量的值后，得到水平方向Teager主能量图和竖直方向Teager主能量图，计算公式如下：

tm_h(x，y)＝{tm(i，f)，i∈[mir，mar]}

tm_v(x，y)＝{tm(j，f)，j∈[mic，mac]}

其中，tm_h(x，y)表示水平方向Teager主能量图，tm_v(x，y)表示竖直方向Teager主能量图，tm(i，f)表示步骤2.5所求的第i行一维信号sigh的Teager主能量，tm(j，f)表示步骤2.5所求的第j列一维信号sigv的Teager主能量，mir表示最小行数值，mar表示最大行数值，mic表示最小列数值，mac表示最大列数值。

为量化和凸显结节的时频信息所引起的能量异常，通过Teager能量算子，计算时频域的Teager主能量，从而凸显结节的异常和压制部分细小血管的响应；水平方向Teager主能量如图3所示，竖直方向Teager主能量如图4所示，从图中可以看出，Teager主能量不仅检测到了结节的能量异常，而且抑制了部分血管的响应，利于提高筛选准确度。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (6)

1.一种基于广义S变换和Teager属性的肺部CT图像结节筛选方法，其特征在于，包括以下步骤：

步骤1：对输入的肺部CT图像进行水平方向和竖直方向上的广义S变换，得到水平方向时频谱和竖直方向时频谱；

步骤2：对所述水平方向时频谱和竖直方向时频谱进行Teager主能量的提取，得到水平方向Teager主能量图和竖直方向Teager主能量图；

步骤3：将所述水平方向Teager主能量图和竖直方向Teager主能量图进行阈值分割，得到疑似结节；

所述步骤1中对输入的肺部CT图像进行预处理，所述预处理包括确定肺部CT图像行变量和列变量的取值范围和去除肺部以外的像素，具体如下包括如下：

步骤a：确定输入的肺部CT图像中行变量和列变量的取值范围即非零像素，并计算非零像素所在的最小行mir、最大行mar、最小列mic和最大列mac，计算公式如下：

{(x，y)|I(x，y)＞0}

mir＝min(x)，mar＝max(x)，mic＝min(y)，mac＝max(y)

其中，I(x，y)表示肺部CT图像在点(x，y)处的灰度值，x是水平方向位置变量，y是竖直方向位置变量，min表示最小值，max表示最大值，mir表示最小行数值，mar表示最大行数值，mic表示最小列数值，mac表示最大列数值；

步骤b：初始化水平方向行变量i、竖直方向列变量j和广义S变换的窗口调节参数p；

步骤c：提取CT图像水平方向第i行一维信号和第j列一维信号，去除零像素得到水平方向第i行一维信号sigh和竖直方向第j列一维信号sigv。

2.如权利要求1所述的一种基于广义S变换和Teager属性的肺部CT图像结节筛选方法，其特征在于，所述步骤b具体为：

初始化广义S变换的窗口调节参数p＝0.1；

计算水平方向时频谱时，初始化水平方向行变量i＝mir，i∈mir～mar；

计算竖直方向时频谱时，初始化竖直方向列变量j＝mic，j∈mic～mac。

3.如权利要求1所述的一种基于广义S变换和Teager属性的肺部CT图像结节筛选方法，其特征在于，所述步骤c中，计算第i行一维信号sigh和竖直方向第j列一维信号sigv的公式如下：

sigh＝{I(x，：)|I(x，：)＞0，I(x，：)∈I(i，：)}

sigv＝{I(：,y)|I(：,y)>O,I(：,y)∈I(：,j)}

其中，I(i，：)表示肺部CT图像水平方向的第i行一维信号，I(x，：)表示第i行一维信号在处的灰度值；I(：,j)表示肺部CT图像竖直方向的第j行一维信号，I(：，y)表示第j列一维信号在y处的灰度值；sigh表示水平方向第i行一维信号，sigv表示竖直方向第j列一维信号。

4.如权利要求1或者2或者3所述的一种基于广义S变换和Teager属性的肺部CT图像结节筛选方法，其特征在于，所述步骤1包括如下步骤：

步骤1.1：初始化水平方向频率变量和竖直方向频率变量：具体如下：

计算水平方向第i行一维信号广义S变换时，初始化水平方向频率变量f_x＝0，f_x取值范围为0～m，m表示水平方向第i行一维信号sigh的长度；

计算竖直方向第j列一维信号广义S变换时，初始化水平方向频率变量f_y＝0，f_y取值范围为0～n，n表示竖直方向第j列一维信号sigv的长度；

步骤1.2：根据步骤1.1的水平方向频率变量和竖直方向频率变量，计算水平方向第i行一维信号sigh和竖直方向第j列一维信号sigv对应的平移后的频谱图；

步骤1.3：根据步骤1.1的水平方向频率变量和竖直方向频率变量，计算频率域高斯窗函数，计算公式如下：

其中，G(f)表示频率域高斯窗函数，f表示水平方向频率变量或者竖直方向频率变量，α表示频率变量，p为窗口调节参数；

步骤1.4：计算步骤1.2所述平移后的频谱图和步骤1.3所述频率域高斯窗函数的乘积；

步骤1.5：对步骤1.4的结果进行反傅里叶变换，分别得到水平方向时频谱和竖直方向时频谱，反傅里叶变换的公式如下：

其中，M(f)表示步骤1.4的乘积结果，f表示水平方向频率变量或者竖直方向频率变量，α表示已知的频率变量。

5.如权利要求4所述的一种基于广义S变换和Teager属性的肺部CT图像结节筛选方法，其特征在于，所述步骤2包括如下步骤：

步骤2.1：初始化位置参数，具体地：

水平方向初始化时，变量x＝2，x∈2～m-1，m表示水平方向第i行一维信号sigh的长度；

竖直方向初始化时，变量y＝2，y∈2～n-1，n表示竖直方向第j列一维信号sigv的长度；

步骤2.2：根据步骤1.5求得在频率为f时的时频谱S(t，f，p)，计算出TK能量，其计算公式如下：

retk(t，f)＝re(S(t，f，p))²-re(S(t-1，f，p))*re(S(t+1，f，p))

imtk(t，f)＝im(S(t，f，p))²-im(S(t-1，f，p))*im(S(t+1，f，p))

tk(t，f)＝retk(t，f)+imtk(t，f)

其中，re()表示取实部，im()表示取虚部，retk(t，f)表示实部TK能量，imtk(t，f)表示虚部TK能量，tk(t，f)表示TK能量；

步骤2.3：执行步骤2.2后，判断水平方向位置变量和竖直方向位置变量的值是否循环遍历了所有值，若否，则将水平方向位置变量和竖直方向位置变量自增后，继续执行步骤2.2，否则跳出循环后跳至步骤2.4；

步骤2.4：判断水平方向频率变量和竖直方向频率变量的值是否循环遍历了所有值，若否，则将水平方向频率变量和竖直方向频率变量自增后，继续执行步骤1.5到步骤2.3，否则跳出循环后跳至步骤2.5；

步骤2.5：在循环遍历所有频率数值后，计算Teager主能量，其计算公式为：

其中，tm(t，f)表示Teager主能量，tk(t，f)表示步骤2.2求得的TK能量，f表示频率变量；

步骤2.6：执行步骤2.5后，判断水平方向行变量和竖直方向列变量的值是否循环遍历了所有值，若否，则将水平方向行变量和竖直方向列变量自增后，继续执行步骤1.3到步骤2.5，否则跳出循环后跳至步骤2.7；

步骤2.7：循环遍历所有水平方向行变量和竖直方向列变量的值后，得到水平方向Teager主能量图和竖直方向Teager主能量图，计算公式如下：

tm_h(x，y)＝{tm(i，f)，i∈[mir，mar]}

tm_v(x，y)＝{tm(j，f)，j∈[mic，mac]}

其中，tm_h(x，y)表示水平方向Teager主能量图，tm_v(x，y)表示竖直方向Teager主能量图，tm(i，f)表示步骤2.5所求的第i行一维信号sigh的Teager主能量，tm(j，f)表示步骤2.5所求的第j列一维信号sigv的Teager主能量，mir表示最小行数值，mar表示最大行数值，mi表示最小列数值，mac表示最大列数值。

6.如权利要求1或者5所述的一种基于广义S变换和Teager属性的肺部CT图像结节筛选方法，其特征在于，步骤3具体为：

步骤3.1：对步骤2求得的水平方向Teager主能量图和竖直方向Teager主能量图进行归一化、阈值分割，分别得到水平方向二值化图像NB_x(x，y)和竖直方向二值化图像NB_y(x，y)，阈值T＝0.7，归一化的公式如下：

ntm(x，y)＝(tm(x，y)-min(tm(x，y)))/(max(tm(x，y))-min(tm(x，y))))

其中，ntm(x，y)表示归一化的Teager主能量，tm(x，y)为步骤2所求的Teager主能量，max()表示最大值，min()表示最小值。

步骤3.2：对所述水平方向二值化图像NB_x(x，y)和竖直方向二值化图像NB_y(x，y)进行与操作，筛选出疑似结节。