CN105869137A - 一种基于多尺度技术的肿瘤纹理特征的提取方法 - Google Patents

Abstract

本发明公开了一种基于多尺度技术的肿瘤纹理特征的提取方法，包括：从临床获取肿瘤PET图像数据，勾画理想肿瘤治疗区域；按照勾画的靶区边界分割并提取大体肿瘤体积GTV，获得只包含GTV存在的PET图像；对每幅PET图像进行灰度级量化，对量化后图像进行二维双树复小波变换，每幅图像得到若干个子带图像；将所有PET图像分解得到的各子带图像重新组合成新立体GTV；计算二维双树复小波变换后的细节系数矩阵的均值和方差，并计算肿瘤纹理特征值；分别建立纹理特征值与标准摄取值的最大值SUVmax、肿瘤TNM分期和浸润深度之间的关系，对量化肿瘤的异质性的纹理特征值进行筛选和优化。

Description

一种基于多尺度技术的肿瘤纹理特征的提取方法

技术领域

本发明涉及图像处理领域，尤其涉及一种基于多尺度技术的肿瘤纹理特征的提取方法。

背景技术

纹理特征是图像最基本的视觉特征，纹理分析技术一直是医学图像处理领域的研究热点。用F-18标记的氟脱氧葡萄糖正电子断层显像(PET)已被广泛应用于预测肿瘤患者放化疗疗效、指导个体治疗及预后中。近些年的研究表明，PET图像中肿瘤的灰度纹理变化特征可以有效评价肿瘤内FDG摄取的高低和非均一程度，从而可以有效量化肿瘤的异质性程度，获得的异质性参数可以用来预测放化疗疗效或指导放疗。由于PET图像灰度变化较复杂，其中隐藏有大量的纹理细节信息，因此寻找合适的纹理特征提取方法与多纹理特征相结合来分析肿瘤的异质性程度，从而帮助临床医生制定更为科学合理的放疗计划是非常必要的。

灰度共生矩阵是分析图像纹理特征的经典统计方法，但它是对整幅图像的总体进行分析处理，缺少对局部细节纹理信息的利用，与人类视觉模型脱节，难于研究纹理尺度像素间的依赖关系。鉴于纹理本身具有多尺度的特性，小波变换以其良好的多分辨率多尺度特性为纹理特征的提取赢得了较好的空间。传统的金字塔小波变换在各分解级仅对低频部分进行分解，在纹理特征提取时仅利用了纹理图像低频子带的信息，对于高频子带中的重要纹理特征信息无法利用。因此，针对PET图像的肿瘤内纹理特征，有必要研究在不同尺度和不同子带上的丰富纹理细节，这对于实时精确量化肿瘤的异质性至关重要。

发明内容

本发明提出了一种基于多尺度技术的肿瘤纹理特征的提取方法，该方法能有效地提取PET图像中多尺度、多方向性的纹理特征，对PET图像进行多层次的紧致描述，同时兼顾了纹理特征的局部随机性和整体规律性的特点。

为实现上述目的，本发明采用如下技术方案：

一种基于多尺度技术的肿瘤纹理特征的提取方法，包括以下步骤：

步骤(1)：从临床获取肿瘤PET图像数据，勾画理想肿瘤治疗区域；

步骤(2)：将步骤(1)中勾画好理想肿瘤治疗区域的PET图像读入MATLAB中，按照勾画的靶区边界分割并提取大体肿瘤体积GTV，并精确提取肿瘤区域GTV，获得只包含GTV存在的PET图像；GTV指临床可见，通过诊断检查证实的肿瘤范围；

步骤(3)：对每幅PET图像进行灰度级量化，对量化后图像进行若干层的二维双树复小波变换，每幅图像得到若干个子带图像；将所有PET图像分解得到的各子带图像按照不同层次不同方向进行分类，同层次同方向的子带图像按原图像编号顺序重新组合成新立体GTV；

步骤(4)：计算二维双树复小波变换后的细节系数矩阵的均值和方差，并选取步骤(3)形成的新立体GTV内三维空间若干个不同方向来统计灰度共生矩阵GLCM，并计算肿瘤纹理特征值；

步骤(5)：分别建立纹理特征值与标准摄取值的最大值SUVmax、肿瘤TNM分期和浸润深度之间的关系，对量化肿瘤的异质性的纹理特征值进行筛选和优化。

所述步骤(1)中勾画肿瘤治疗区域的具体过程，包括：

步骤(1.1)：PET图像读入MATLAB中；

步骤(1.2)：获取步骤(1.1)中PET图像对应的包含SUV值的图像，先获取整个图像中SUV值最大的点值为初始点，然后基于初始点进行区域增长，直至包含所有医生建议阈值的区域，获得肿瘤区域；

步骤(1.3)：勾画步骤(1.2)所得图像中标准摄取值SUV值大于医生建议阈值的区域，获得肿瘤区域；

步骤(1.4)：将步骤(1.3)所获取的肿瘤区域图像再读入MATLAB中，再经过步骤(1.2)和步骤(1.3)获得理想肿瘤治疗区域。

所述步骤(3)的具体过程为：

采用二维双树复小波变换对PET图像进行多尺度小波分解，将每幅PET图像首先进行三层二维双树复小波变换，每层变换得到六个方向的高频子带图像，所述高频子带图像所在的方向分别是±15°、±45°和±75°；三层分解后每幅图像可以得到3×6＝18个高频子带图像。

所述二维双树复小波变换的小波基函数为：

式中，表示母小波；s表示尺度因子，s为整数；k₀表示位移因子。

所述的二维双树复小波域小波函数为：

ψ(x,y)＝[ψ_h(x)+jψ_g(x)][ψ_h(y)+jψ_g(y)]

＝ψ_h(x)ψ_h(y)－ψ_g(x)ψ_g(y)+j[ψ_g(x)ψ_h(y)+ψ_h(x)ψ_g(y)]

其中，ψ_h(x),ψ_g(x)表示树a中正交的实小波，且能形成Hilbert变换对；ψ_h(y)，ψ_g(y)表示树b中正交的小波，且能形成Hilbert变换对。

所述步骤(4)中求取细节系数矩阵的均值和方差的表达式为：

${\mathrm{\μ}}_i=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}i\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{P}_d(i,j)$

${\mathrm{\σ}}_i=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{(i-{\mathrm{\μ}}_i)}^2\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{P}_d(i,j)$

其中，μ_i和σ_i分别表示灰度为i的子带图像的均值和方差，M代表图像的长度，N代表图像的宽度，P_d(i,j)是细节系数矩阵。

所述灰度共生矩阵GLCM用P(i,j|d,θ)表示，灰度共生矩阵GLCM为从灰度为i的像素统计与其距离为d，方向为θ上，灰度为j的像素对(i,j)出现的联合概率密度。

所述纹理特征值，包括：对比度、相关度、角二阶距和熵，其中，

对比度CON的表达式为：

$\mathrm{CON}=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{|i-j|}^{2}P(i,j)$

式中，i,j表示像素的灰度，P(i,j)表示灰度共生矩阵，其中，P(i,j)＝P(i,j|d,θ)。

相关度COR的表达式为：

$\mathrm{COR}=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\frac{(i-{\mathrm{\μ}}_i)(j-{\mathrm{\μ}}_j)P(i,j)}{{\mathrm{\σ}}_i{\mathrm{\σ}}_j}$

式中，μ_i和μ_j分别表示像素灰度为i和j的细节系数矩阵的均值；σ_i和σ_j分别表示像素灰度为i和j的细节系数矩阵的方差。

角二阶距ASM的表达式为：

$\mathrm{ASM}=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}P{(i,j)}^2$

熵ENT的表达式为：

$\mathrm{ENT}=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}P(i,j)\mathrm{lgP}(i,j).$

本发明的有益效果是：

(1)本发明选用双树复小波变换作为PET图像的分解小波，可以全面满足小波变换多尺度多分辨率的要求，更好的体现纹理特征的全方位性；

(2)通过重组分解细化得到的同尺度同方向高频子带图像，从而提取常用的纹理特征。

附图说明

图1为本发明提取纹理特征方法流程图；

图2为本发明的双树复小波分解方向图；

图3为本发明的新GTV重组形成图；

图4为本发明的三维空间中13个不同的邻域方向图。

具体实施方式

下面结合附图与实施对本发明内容做进一步说明：

本发明中提出一种基于多尺度技术的肿瘤纹理特征的提取方法，从而获取更为有效的纹理特征值以帮助临床医生量化肿瘤的异质性、辅助优化放疗剂量、为病人制定更为合理有效的放疗计划。

一种基于多尺度技术的肿瘤纹理特征的提取方法，具体包括以下步骤：

步骤(1)：从临床获取PET图像并去除其床位伪影及噪声成分，并由临床医生结合CT图像在SUV值大于2.5区域再次编辑勾画确定出肿瘤治疗区域GTV；

步骤(2)：勾画好GTV的PET图像导入MATLAB软件，对勾画区域进行分割提取并保存只含GTV区域的图像；

步骤(3)：将每幅PET图像进行16级灰度等级量化，量化后的PET图像采用二维双树复小波变换进行三层小波分解，每层都得到六个方向的高频子带图像；

将所有子带图像按不同尺度不同方向子带分类，同层次同方向子带图像分到一组共形成18组，每一组子带图像按原图形编号顺序排列构成一个新的三维立体GTV，选取新GTV内任一体素周围的13个互不重复的方向统计灰度共生矩阵，计算纹理参数，并求得各个方向上的平均值；

步骤(4)：计算二维双树复小波变换后的细节系数矩阵的均值和方差，将得到的数据归一化；在纹理特征提取时，在各个新GTV内d取1,选择任一体素周围的13个互不重复的方向统计三维灰度共生矩阵，计算各方向的纹理参数，并取各方向平均值作为异质性参数；

本发明使用的均值和方差的计算公式为：

${\mathrm{\μ}}_i=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}i\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{P}_d(i,j)$

${\mathrm{\σ}}_i=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{(i-{\mathrm{\μ}}_i)}^2\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{P}_d(i,j)$

其中，μ_i和σ_i分别表示灰度为i的子带图像的均值和方差，M代表图像的长度，N代表图像的宽度，P_d(i,j)是细节系数矩阵。

步骤(5)：使用SPSS统计分析方法分别建立纹理特征值与标准摄取值的最大值SUVmax、肿瘤TNM分期和浸润深度之间的关系，将各方向提取的纹理特征值取平均值作为异质性参数，并与临床参数相结合来进行统计分析，对特征值进行筛选和优化。

所述步骤(3)中所用到的二维双树复小波变换是一个多尺度多分辨率的工具，选用如下小波基函数：

式中，表示母小波；s表示尺度因子，s为整数；k₀表示位移因子。每经过一级分解后，小波基函数就把输入图像分解为一组不同方向的函数。

双树复小波变换采用两颗平行的小波树(Tree a和Tree b)作为复数的实部和虚部，滤波器H_0a、H_1a与H_0b、H_1b之间的延迟作为一个采样值的间隔，来确保采样过程中能够全面取到所需采样点。

二维双树复小波变换复数域小波函数表达式如下：

ψ(x,y)＝[ψ_h(x)+jψ_g(x)][ψ_h(y)+jψ_g(y)]

＝ψ_h(x)ψ_h(y)－ψ_g(x)ψ_g(y)+j[ψ_g(x)ψ_h(y)+ψ_h(x)ψ_g(y)]

实部小波函数为：

ψ_i(x，y)＝ψ_1，i(x，y)－ψ_2，i(x，y)

ψ_i+3(x，y)＝ψ_1，i(x，y)+ψ_2，i(x，y)

式中，i＝1，2，3。

其中，ψ_1，i(x，y)和ψ_2，i(x，y)分别为六个方向的小波函数，定义为:

ψ_1，1(x，y)＝φ_h(x)ψ_h(y)

ψ_1，2(x，y)＝ψ_h(x)φ_h(y)

ψ_1，3(x，y)＝ψ_h(x)ψ_h(y)

ψ_2，1(x，y)＝φ_g(x)ψ_g(y)

ψ_2，2(x，y)＝ψ_g(x)φ_g(y)

ψ_2，3(x，y)＝ψ_g(x)ψ_g(y)

虚部小波函数为：

ψ_i(x，y)＝ψ_3，i(x，y)+ψ_4，i(x，y)

ψ_i+3(x，y)＝ψ_3，i(x，y)－ψ_4，i(x，y)

其中，ψ_3，i(x，y)和ψ_4，i(x，y)分别为六个方向的小波函数，定义为

ψ_3，1(x，y)＝φ_h(x)ψ_h(y)

ψ_3，2(x，y)＝ψ_h(x)φ_h(y)

ψ_3，3(x，y)＝ψ_h(x)ψ_h(y)

ψ_4，1(x，y)＝φ_g(x)ψ_g(y)

ψ_4，2(x，y)＝ψ_g(x)φ_g(y)

ψ_4，3(x，y)＝ψ_g(x)ψ_g(y)

按照上式对图像进行双树复小波分解，可以得到多个分解层，每个分解层有六个方向的细节子带图像，分别用I_k,α表示，其中，k表示分解层数，α表示分解方向。

本实施例选用3层DT-CWT对PET图像进行分解，并得到18个方向上的高频子带图像作为研究对象，将获得的同层次同方向的子带图像按原图形编号顺序重新组合成一个新的立体GTV，共得到18个新GTV。

本发明采用二维双树复小波变换来进行图像变换，图像经双树复小波变换后在每一尺度上都产生±15°、±45°、±75°六个方向的高频复系数子带图像和一个低频子带图像。各高频子带是对原图滤波后抽取的不同方向上的对应纹理信息，并组成各自只包含此特定方向信息的子带图像。低频子带图像看作是原图形的近似，作为下一层分解的输入。因此每进行一层分解后就将图像纹理信息细化了一层，这样分解后的高频子带图像在局部结构上较传统小波变换具有了更多的方向和细节信息。

所述步骤(4)中计算灰度共生矩阵GLCM：

灰度共生矩阵GLCM为从灰度为i的像素统计与其距离为d，方向为θ上，灰度为j的像素对(i,j)出现的联合概率密度。为全面体现肿瘤纹理灰度的空间变化，在新GTV内d取1，选任一体素周围的13个互不重复的方向统计灰度共生矩阵。

使用医学图像分析中较为常用的具有代表性的纹理参数：对比度(CON)、相关度(COR)、角二阶距(ASM)、熵(ENT)四个纹理特征值来量化肿瘤的异质性。

四个参数所代表的意义及计算公式如下：

(1)对比度CON：反应图像特定位置关系下的像素对的灰度差，即清晰度。对比度越大，图像中纹理基元对比越强烈，视觉效果越清晰。

因此，对比度值的大小反映了图像纹理的粗细程度。对于粗纹理，对比度值较小，对于细纹理，对比度值较大，肿瘤摄取异质性越高对比度值越大。

$\mathrm{CON}=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{|i-j|}^{2}P(i,j)$

式中，i,j表示像素的灰度，P(i,j)表示灰度共生矩阵，其中，P(i,j)＝P(i,j|d,θ)。

(2)相关度COR：反映图像纹理区域在某个方向上的相似度，是用来度量图像局部灰度线性相关性的，肿瘤摄取异质性越高时，其值越小。

$\mathrm{COR}=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\frac{(i-{\mathrm{\μ}}_i)(j-{\mathrm{\μ}}_j)P(i,j)}{{\mathrm{\σ}}_i{\mathrm{\σ}}_j}$

式中，μ_i和μ_j分别表示像素灰度为i和j的细节系数矩阵的均值；σ_i和σ_j分别表示像素灰度为i和j的细节系数矩阵的方差。

(3)角二阶距ASM：当矩阵中少数的元素值较大即特定的像素对较多时，图像有较好的一致性，灰度分布均匀，能量较大，它反映了图像灰度分布的一致性程度，肿瘤异质性程度越高时其值越小。

$\mathrm{ASM}=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}P{(i,j)}^2$

(4)熵ENT：反映图像中纹理复杂程度，图像中纹理密集时熵值大，纹理稀疏时熵值小，肿瘤摄取异质性越高时其值越大。

$\mathrm{ENT}=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}P(i,j)\mathrm{lgP}(i,j)$

双树复小波变换具有保持有限数据冗余、计算效率高等特点，使用双树复小波变换来进行图像分解，保证了样本在空间的分布均匀性及近似的平移不变性。双树复小波变换在保持传统小波变换具有的良好时频局部化分析能力的同时，还具有良好的方向分析能力、近似的平移不变性、有限的数据冗余性和能够完全重建的能力。

上述内容虽然结合附图对本发明的具体实施方式过程进行了描述，但并不是对本发明保护范围的限制，所属领域技术人员明白，在本发明技术方案基础上，本领域技术人员可以不付出创造性劳动就能够做出的各种修改或变形仍在本发明的保护范围以内。

Claims (10)

1.一种基于多尺度技术的肿瘤纹理特征的提取方法，其特征在于，包括以下步骤：

步骤(1)：从临床获取肿瘤PET图像数据，勾画理想肿瘤治疗区域；

步骤(2)：将步骤(1)中勾画好理想肿瘤治疗区域的PET图像读入MATLAB中，按照勾画的靶区边界分割并提取大体肿瘤体积GTV，获得只包含GTV存在的PET图像；

步骤(3)：对每幅PET图像进行灰度级量化，对量化后图像进行若干层的二维双树复小波变换，每幅图像得到若干个子带图像；将所有PET图像分解得到的各子带图像按照不同层次不同方向进行分类，同层次同方向的子带图像按原图像编号顺序重新组合成新立体GTV；

步骤(4)：计算二维双树复小波变换后的细节系数矩阵的均值和方差，并选取步骤(3)形成的新立体GTV内三维空间若干个不同方向来统计灰度共生矩阵GLCM，并计算肿瘤纹理特征值；

步骤(5)：分别建立纹理特征值与标准摄取值的最大值SUVmax、肿瘤TNM分期和浸润深度之间的关系，对量化肿瘤的异质性的纹理特征值进行筛选和优化。

2.如权利要求1所述的一种基于多尺度技术的肿瘤纹理特征的提取方法，其特征在于，所述步骤(1)中勾画肿瘤治疗区域的具体过程，包括：

步骤(1.1)：PET图像读入MATLAB中；

步骤(1.2)：获取步骤(1.1)中PET图像对应的包含SUV值的图像，先获取整个图像中SUV值最大的点值为初始点，然后基于初始点进行区域增长，直至包含所有医生建议阈值的区域，获得肿瘤区域；

步骤(1.3)：勾画步骤(1.2)所得图像中标准摄取值SUV值大于医生建议阈值的区域，获得肿瘤区域；

步骤(1.4)：将步骤(1.3)所获取的肿瘤区域图像再读入MATLAB中，再经过步骤(1.2)和步骤(1.3)获得理想肿瘤治疗区域。

3.如权利要求1所述的一种基于多尺度技术的肿瘤纹理特征的提取方法，其特征在于，所述步骤(3)的具体过程为：

采用二维双树复小波变换对PET图像进行多尺度小波分解，将每幅PET图像首先进行三层二维双树复小波变换，每层变换得到六个方向的高频子带图像，三层分解后每幅图像可以得到3×6＝18个高频子带图像。

4.如权利要求3所述的一种基于多尺度技术的肿瘤纹理特征的提取方法，其特征在于，所述高频子带图像所在的方向分别是±15⁰、±45⁰和±75⁰。

5.如权利要求3所述的一种基于多尺度技术的肿瘤纹理特征的提取方法，其特征在于，所述二维双树复小波变换的小波基函数为：

式中，表示母小波；s表示尺度因子，s为整数；k₀表示位移因子。

6.如权利要求3所述的一种基于多尺度技术的肿瘤纹理特征的提取方法，其特征在于，所述的二维双树复小波域小波函数为：

ψ(x,y)＝[ψ_h(x)+jψ_g(x)][ψ_h(y)+jψ_g(y)]

＝ψ_h(x)ψ_h(y)－ψ_g(x)ψ_g(y)+j[ψ_g(x)ψ_h(y)+ψ_h(x)ψ_g(y)]

其中，ψ_h(x),ψ_g(x)表示树a中正交的实小波，且能形成Hilbert变换对；ψ_h(y)，ψ_g(y)表示树b中正交的小波，且能形成Hilbert变换对。

7.如权利要求1所述的一种基于多尺度技术的肿瘤纹理特征的提取方法，其特征在于，所述步骤(4)中求取细节系数矩阵的均值和方差的表达式为：

${\mathrm{\μ}}_i=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}i\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{P}_d(i,j)$

${\mathrm{\σ}}_i=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{(i-{\mathrm{\μ}}_i)}^2\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{P}_d(i,j)$

其中，μ_i和σ_i分别表示灰度为i的子带图像的均值和方差，M代表图像的长度，N代表图像的宽度，P_d(i,j)是细节系数矩阵。

8.如权利要求1所述的一种基于多尺度技术的肿瘤纹理特征的提取方法，其特征在于，所述灰度共生矩阵GLCM用P(i,j|d,θ)表示，灰度共生矩阵GLCM为从灰度为i的像素统计与其距离为d，方向为θ上，灰度为j的像素对(i,j)出现的联合概率密度。

9.如权利要求1所述的一种基于多尺度技术的肿瘤纹理特征的提取方法，其特征在于，所述纹理特征值，包括：对比度和相关度，其中，

对比度CON的表达式为：

$\mathrm{CON}=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{|i-j|}^{2}P(i,j)$

式中，i,j表示像素的灰度，P(i,j)表示灰度共生矩阵，其中，P(i,j)＝P(i,j|d,θ)；

相关度COR的表达式为：

$\mathrm{COR}=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\frac{(i-{\mathrm{\μ}}_i)(j-{\mathrm{\μ}}_j)P(i,j)}{{\mathrm{\σ}}_i{\mathrm{\σ}}_j}$

式中，μ_i和μ_j分别表示像素灰度为i和j的细节系数矩阵的均值；σ_i和σ_j分别表示像素灰度为i和j的细节系数矩阵的方差。

10.如权利要求1所述的一种基于多尺度技术的肿瘤纹理特征的提取方法，其特征在于，所述纹理特征值，还包括：角二阶距ASM和熵ENT，其中，

角二阶距ASM的表达式为：

$\mathrm{ASM}=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}P{(i,j)}^2$

熵ENT的表达式为：

$\mathrm{ENT}=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}P(i,j)\mathrm{lgP}(i,j).$