使用pyradiomics提取影像组学特征【详细】

最近由于项目需求要使用pyradiomics提取影像组学特征，网上阅读了很多别人的博客，学到一些，然后去查看了pyradiomics的官方文档，最后自己实现了特征的提取，写下此文记录，方便日后查看。

首先放上官方文档：https://pyradiomics.readthedocs.io/
和pyradiomics源代码地址：https://github.com/AIM-Harvard/pyradiomics

第一步、安装

在已经安装好的Python环境下，执行：

pip install pyradiomics

第二步、提取特征

通过自定义提取特征可以根据自己的需求来设置并提取。

• 首先，可以指定用于提取特征的图像类型。在pyradiomics包中为我们提供了许多可以使用的滤波器，所以我们可以使用原始图像及经过各种滤波器之后的图像，如下方表格，具体可以参考官方：图像类型

图像类型	解释
Original	原始图像
Wavelet	小波变换。产生在三个维度中每个维度分别使用高通、低通滤波器的所有组合（LLH、LHL、LHH、HLL、HLH、HHL、HHH、LLL）
LoG	Laplacian of Gaussian filter高斯滤波器的拉普拉斯算子，是一种边缘增强滤波器。使用它需要指定参数sigma，低 sigma 强调精细纹理，高 sigma 值强调粗糙纹理
Square	平方。取原始像素的平方并将它们线性缩放回原始范围
SquareRoot	平方根。取绝对图像强度的平方根并将它们缩放回原始范围
Logarithm	对数。取绝对强度 + 1 的对数，值缩放到原始范围
Exponential	指数。采用e^(绝对强度)获取强度的指数值，值被缩放到原始范围
Gradient	梯度。返回局部梯度的大小
LocalBinaryPattern2D	在每一片中进行的本地二进制模式
LocalBinaryPattern3D	在3d中进行的本地二进制模式

指定方式：设置特征提取器后可以在下边指定（如何设置特征提取器：extractor）
可以直接使用：enableAllImageTypes()启用所有类型
也可以使用：enableImageTypeByName( imageType , enabled=True , customArgs=None )启用你想用的类型
例如：

# 所有类型
extractor.enableAllFeatures()
# 指定使用LoG和Wavelet滤波器
extractor.enableImageTypeByName('LoG')
extractor.enableImageTypeByName('Wavelet')

• 然后可以指定所要提取的特征。pyradiomics包也为我们提供了很多种可选的特征，如下方表格，具体特征类型详细资料可参考：radiomics-features

特征类型	特征	解释
First Order Features（共19个）	Energy	能量
	Total Energy	总能量
	Entropy	熵
	Minimum	最小特征值
	10Percentile	特征值的百分之10的值
	90Percentile	特征值得百分之90的值
	Maximum	最大特征值
	Mean	均值
	Median	中位数
	InterquartileRange	四分位距离
	Range	灰度值范围
	MeanAbsoluteDeviation	(MAD)平均绝对误差
	RobustMeanAbsoluteDeviation	(rMAD) 鲁棒平均绝对偏差
	RootMeanSquared	(RMS)均方根误差
	StandardDeviation	标准差。测量平均值的变化或离散量，默认不启用，因为与方差相关
	Skewness	偏度。测量值的分布关于平均值的不对称性
	Kurtosis	峰度。是图像 ROI 中值分布的“峰值”的量度
	Variance	方差。是每个强度值与平均值的平方距离的平均值
	Uniformity	均匀度。是每个强度值的平方和

特征类型	特征	解释
Shape Features (3D) （共17个）	MeshVolume	网格体积
	VoxelVolume	体素体积
	SurfaceArea	表面积
	SurfaceVolumeRatio	表面积体积比
	Sphericity	球形度
	Compactness1	紧凑度1。衡量肿瘤形状相对于球形的紧密程度，由于与球形度、紧凑度2、球形不对称相关，默认不启用
	Compactness2	紧凑度2。衡量肿瘤形状相对于球形的紧密程度，与1的计算公式不同，默认不启用
	SphericalDisproportion	球形不对称。是肿瘤区域的表面积与与肿瘤区域具有相同体积的球体的表面积之比，是球形度的倒数，默认不启用
	Maximum3DDiameter	最大3D直径
	Maximum2DDiameterSlice	最大二维直径（切片），定义为行列（通常是轴向）平面中肿瘤表面网格顶点之间的最大成对欧几里得距离。
	Maximum2DDiameterColumn	最大2D直径（列）。定义为行切片（通常是冠状）平面中肿瘤表面网格顶点之间的最大成对欧几里得距离。
	Maximum2DDiameterRow	最大二维直径（行）。定义为列切片（通常为矢状）平面中肿瘤表面网格顶点之间的最大成对欧几里得距离。
	MajorAxisLength	长轴长度
	MinorAxisLength	短轴长度
	LeastAxisLength	最小轴长度
	Elongation	伸长率。表示ROI形状中两个最大的主成分之间的关系
	Flatness	平整度。表示ROI形状中最大和最小主成分之间的关系

特征类型	特征	解释
Shape Features (2D)	MeshSurface	网格面。由三角形网格定义，计算ROI的表面
	PixelSurface	像素面。ROI 中的像素数乘以单个像素的表面积
	Perimeter	周长
	PerimeterSurfaceRatio	周长比。值越低越紧凑
	Sphericity	球形度
	SphericalDisproportion	球形不均匀度。是球形度的倒数 ，与球形度相关，默认不启用
	MaximumDiameter	最大二维直径。定义为肿瘤表面网格顶点之间的最大成对欧几里得距离
	MajorAxisLength	长轴长度
	MinorAxisLength	短轴长度
	Elongation	伸长率。显示了 ROI 形状中两个最大主成分之间的关​​系

特征类型	特征	解释
Gray Level Co-occurrence Matrix (GLCM) Features灰度共生矩阵特征（共24个）	Autocorrelation	自相关。是纹理细度和粗糙度大小的量度
	JointAverage	联合平均
	ClusterProminence	集群显著性。衡量 GLCM 偏度和不对称性的指标
	ClusterShade	集群阴影。衡量 GLCM 偏度和均匀度的指标
	ClusterTendency	集群趋势
	Contrast	对比度
	Correlation	相关性
	DifferenceAverage	平均差
	DifferenceEntropy	差分熵
	DifferenceVariance	差异方差。是异质性的一种度量
	Dissimilarity	差异性。已弃用
	JointEnergy	联合能量
	JointEntropy	联合熵。是邻域强度值的随机性/可变性的度量
	Imc1	相关性信息度量 (IMC) 1
	Imc2	相关性信息度量 (IMC) 2
	Idm	逆差矩
	MCC	最大相关系数
	Idmn	逆差分矩归一化
	Id	逆差
	Idn	逆差分归一化
	InverseVariance	逆方差
	MaximumProbability	最大概率
	SumAverage	求和平均
	SumEntropy	总熵
	SumSquares	平方和

特征类型	特征	解释
Gray Level Size Zone Matrix (GLSZM) Features灰度级大小区域矩阵特征 （共16个）	SmallAreaEmphasis	是小尺寸区域分布的量度，数值越大表示区域越小，纹理越精细
	LargeAreaEmphasis	是大面积区域分布的度量，值越大表示区域越大，纹理越粗糙。
	GrayLevelNonUniformity	灰度非均匀性。测量图像中灰度级强度值的可变性，值越低表示强度值越均匀。
	GrayLevelNonUniformityNormalized	灰度非均匀归一化。测量图像中灰度强度值的可变性，值越低表示强度值的相似性越大
	SizeZoneNonUniformity	尺寸区域非均匀性。测量图像中大小区域体积的可变性，值越低表示大小区域体积越均匀
	SizeZoneNonUniformityNormalized	尺寸区域非均匀归一化。测量整个图像中大小区域体积的可变性，值越低表示图像中区域大小体积之间的同质性越高
	ZonePercentage	区域百分比。获取 ROI 中的区域数和体素数的比率来测量纹理的粗糙度
	GrayLevelVariance	灰度方差
	ZoneVariance	区域方差
	ZoneEntropy	区域熵
	LowGrayLevelZoneEmphasis	低灰度区域强调
	HighGrayLevelZoneEmphasis	高灰度区域强调
	SmallAreaLowGrayLevelEmphasis	小面积低灰度区域强调
	SmallAreaHighGrayLevelEmphasis	小区域高灰度强调
	LargeAreaLowGrayLevelEmphasis	大面积低灰度强调
	LargeAreaHighGrayLevelEmphasis	大面积高灰度强调

特征类型	特征	解释
Gray Level Run Length Matrix (GLRLM) Features灰度游程矩阵 (GLRLM) 特征 （共16个）	ShortRunEmphasis	是衡量短运行长度分布的指标，值越大表示运行长度越短，纹理纹理越精细。
	LongRunEmphasis	是对长行程分布的度量，值越大表示行程越长，结构纹理越粗糙。
	GrayLevelNonUniformity	灰度非均匀性
	GrayLevelNonUniformityNormalized	灰度非均匀归一化
	RunLengthNonUniformity	运行长度非均匀性
	RunLengthNonUniformityNormalized	运行长度非均匀归一化
	RunPercentage	运行百分比。通过获取 ROI 中的运行次数和体素数的比率来测量纹理的粗糙度
	GrayLevelVariance	灰度级方差
	RunVariance	运行方差
	RunEntropy	运行熵
	LowGrayLevelRunEmphasis	测量低灰度值的分布，值越高表示图像中低灰度值的集中度越高
	HighGrayLevelRunEmphasis	测量较高灰度值的分布，值越高表示图像中高灰度值的集中度越高
	ShortRunLowGrayLevelEmphasis	测量具有较低灰度值的较短运行长度的联合分布
	ShortRunHighGrayLevelEmphasis	测量具有较高灰度值的较短运行长度的联合分布
	LongRunLowGrayLevelEmphasis	测量具有较低灰度值的长运行长度的联合分布
	LongRunHighGrayLevelEmphasis	测量具有较高灰度值的长运行长度的联合分布

特征类型	特征	解释
Neighbouring Gray Tone Difference Matrix (NGTDM) Features相邻灰度差矩阵（共5个） (NGTDM) 特征	Coarseness	粗糙度
	Contrast	对比度
	Busyness	繁忙度。从像素到其邻居的变化的度量。繁忙度的高值表示“繁忙”图像，像素与其邻域之间的强度快速变化
	Complexity	复杂度
	Strength	强度

特征类型	特征	解释
Gray Level Dependence Matrix (GLDM) Features 灰度依赖矩阵 (GLDM) 特征（共14个）	SmallDependenceEmphasis	小依赖分布的度量，值越大表示依赖越小，纹理越不均匀
	LargeDependenceEmphasis	大依赖分布的度量，值越大表示依赖越大，纹理越均匀
	GrayLevelNonUniformity	灰度非均匀性。测量图像中灰度级强度值的相似性
	DependenceNonUniformity	衡量整个图像中依赖关系的相似性，值越低表示图像中依赖关系之间的同质性越高
	DependenceNonUniformityNormalized	依赖非均匀归一化。衡量整个图像中依赖关系的相似性，值越低表示图像中依赖关系之间的同质性越高
	GrayLevelVariance	灰度方差
	DependenceVariance	依赖方差
	DependenceEntropy	依赖熵
	LowGrayLevelEmphasis	测量低灰度值的分布，值越大表示图像中低灰度值的集中度越高
	HighGrayLevelEmphasis	测量较高灰度值的分布，值越高表示图像中高灰度值的集中度越高
	SmallDependenceLowGrayLevelEmphasis	测量具有较低灰度值的小相关性的联合分布
	SmallDependenceHighGrayLevelEmphasis	测量具有较高灰度值的小相关性的联合分布
	LargeDependenceLowGrayLevelEmphasis	测量具有较低灰度值的大相关性的联合分布
	LargeDependenceHighGrayLevelEmphasis	测量具有较高灰度值的大相关性的联合分布

注意：除了形状特征类外，其他特征都可以在原始图像和滤波后的图像上进行计算

具体的指定方法：
可以直接使用：enableAllFeatures( )启用所有类型
也可以使用：enableFeatureClassByName(featureClass, enabled=True)启用你想用的类型
例如：

# 设置一阶特征
extractor.enableFeatureClassByName('firstorder')
# 设置只提取一阶特征的'Mean'和'Skewness'
extractor.enableFeaturesByName(firstorder=['Mean', 'Skewness'])

也可以直接通过特征类别来启用具体的特征，示例见：https://github.com/AIM-Harvard/pyradiomics/blob/master/examples/helloFeatureClass.py

• 第三步可以对特征提取器进行一些设置。是否对原图归一化、是否重采样，原文档见：feature-extractor-level

图像归一化：

• normalize：默认为false。设置为true时进行图像归一化。
• normalizeScale：确定图像归一化后的比例。默认为1。
• removeOutliers：定义要从图像中删除的异常值。默认为0。

图像/mask重采样：

• resampledPixelSpacing：设置重采样时的体素大小。默认无。
• interpolator：设置用于重采样的插值器。仅适用于重采样图像，sitkNearestNeighbor始终用于重采样掩码以保留标签值。可选的插值器：
  
• padDistance：设置重采样期间裁剪肿瘤体时的体素补充数量。

例如：

settings = {}
settings['![binWidth](https://img-blog.csdnimg.cn/c9b0896a5eea4eaf8217d0ed7f23e92b.png)
'] = 25
settings['resampledPixelSpacing'] = [3,3,3]  # [3,3,3] is an example for defining resampling (voxels with size 3x3x3mm)
settings['interpolator'] = sitk.sitkBSpline

还有其他的一些设置，根据自己的需要修改即可。
最后举一个完整的例子，我要处理的是CT数据，所以我参考了pyradiomics的github源码中ct参数的设置，最后我的代码如下：

from __future__ import print_function
import six
import os  # needed navigate the system to get the input data
import numpy as np
import radiomics
import pandas as pd
import SimpleITK as sitk
from radiomics import featureextractor  # This module is used for interaction with pyradiomics
import argparse
def catch_features(imagePath,maskPath):
    if imagePath is None or maskPath is None:  # Something went wrong, in this case PyRadiomics will also log an error
        raise Exception('Error getting testcase!')  # Raise exception to prevent cells below from running in case of "run all"
    settings = {}
    settings['binWidth'] = 25  # 5
    settings['sigma'] = [3, 5]
    settings['Interpolator'] = sitk.sitkBSpline
    settings['resampledPixelSpacing'] = [1, 1, 1]  # 3,3,3
    settings['voxelArrayShift'] = 1000  # 300
    settings['normalize'] = True
    settings['normalizeScale'] = 100
    extractor = featureextractor.RadiomicsFeatureExtractor(**settings)
    #extractor = featureextractor.RadiomicsFeatureExtractor()
    print('Extraction parameters:\n\t', extractor.settings)

    extractor.enableImageTypeByName('LoG')
    extractor.enableImageTypeByName('Wavelet')
    extractor.enableAllFeatures()
    extractor.enableFeaturesByName(firstorder=['Energy', 'TotalEnergy', 'Entropy', 'Minimum', '10Percentile', '90Percentile', 'Maximum', 'Mean', 'Median', 'InterquartileRange', 'Range', 'MeanAbsoluteDeviation', 'RobustMeanAbsoluteDeviation', 'RootMeanSquared', 'StandardDeviation', 'Skewness', 'Kurtosis', 'Variance', 'Uniformity'])
    extractor.enableFeaturesByName(shape=['VoxelVolume', 'MeshVolume', 'SurfaceArea', 'SurfaceVolumeRatio', 'Compactness1', 'Compactness2', 'Sphericity', 'SphericalDisproportion','Maximum3DDiameter','Maximum2DDiameterSlice','Maximum2DDiameterColumn','Maximum2DDiameterRow', 'MajorAxisLength', 'MinorAxisLength', 'LeastAxisLength', 'Elongation', 'Flatness'])
# 上边两句我将一阶特征和形状特征中的默认禁用的特征都手动启用，为了之后特征筛选
    print('Enabled filters:\n\t', extractor.enabledImagetypes)
    feature_cur = []
    feature_name = []
    result = extractor.execute(imagePath, maskPath, label=255)
        for key, value in six.iteritems(result):
        print('\t', key, ':', value)
        feature_name.append(key)
        feature_cur.append(value)
    print(len(feature_cur[37:]))
    name = feature_name[37:]
    name = np.array(name)
    '''
    flag=1
    if flag:
        name = np.array(feature_name)
        name_df = pd.DataFrame(name)
        writer = pd.ExcelWriter('key.xlsx')
        name_df.to_excel(writer)
        writer.save()
        flag = 0
    '''
    for i in range(len(feature_cur[37:])):
        #if type(feature_cur[i+22]) != type(feature_cur[30]):
        feature_cur[i+37] = float(feature_cur[i+37])
    return feature_cur[37:],name

image_dir = r'D:\study\extract_radiomics_feature\linshi'
mask_dir = r'D:\study\extract_radiomics_feature\linshi'
patient_list = os.listdir(image_dir)
save_file = np.empty(shape=[1,1051])
id = []
for patient in patient_list:
    print(patient)
    for file in os.listdir(os.path.join(image_dir,patient)):
        if file =='image.nii':
            imagePath = os.path.join(image_dir,patient,file)
        if file =='mask_gtv-1.nii':
            maskPath = os.path.join(mask_dir,patient,file)
    print(imagePath)
    print(maskPath)
    save_curdata,name = catch_features(imagePath,maskPath)
    save_curdata = np.array(save_curdata)
    save_curdata = save_curdata.reshape([1, 1051])
    id.append(patient.split('.')[0])
    np.concatenate((patient,save_curdata),axis=1)
    save_file = np.append(save_file,save_curdata,axis=0)
    print(save_file.shape)
save_file = np.delete(save_file,0,0)
#save_file = save_file.transpose()
#print(save_file.shape)
id_num = len(id)
id = np.array(id)
name_df = pd.DataFrame(save_file)
name_df.index = id
name_df.columns = name
writer = pd.ExcelWriter('NSCLC-Radiomics-features.xlsx')
name_df.to_excel(writer)
writer.save()

最后结果如下：