使用pyradiomics提取影像组学特征【详细】
使用提取影像组学特征【详细】
最近由于项目需求要使用提取影像组学特征，网上阅读了很多别人的博客，学到一些，然后去查看了的官方文档，最后自己实现了特征的提取，写下此文记录，方便日后查看。
首先放上官方文档：
和源代码地址：
第一步、安装
在已经安装好的环境下，执行：
pip install pyradiomics
第二步、提取特征
通过自定义提取特征可以根据自己的需求来设置并提取。
图像类型解释
原始图像
小波变换。产生在三个维度中每个维度分别使用高通、低通滤波器的所有组合（LLH、LHL、LHH、HLL、HLH、HHL、HHH、LLL）
LoG
of 高斯滤波器的拉普拉斯算子，是一种边缘增强滤波器。使用它需要指定参数sigma，低 sigma 强调精细纹理，高 sigma 值强调粗糙纹理
平方。取原始像素的平方并将它们线性缩放回原始范围
平方根。取绝对图像强度的平方根并将它们缩放回原始范围
对数。取绝对强度 + 1 的对数，值缩放到原始范围
指数。采用e^(绝对强度)获取强度的指数值，值被缩放到原始范围
梯度。返回局部梯度的大小
在每一片中进行的本地二进制模式
在3d中进行的本地二进制模式
指定方式：设置特征提取器后可以在下边指定（如何设置特征提取器：）
可以直接使用：()启用所有类型
也可以使用：e( , =True , =None )启用你想用的类型
例如：
# 所有类型
extractor.enableAllFeatures()
# 指定使用LoG和Wavelet滤波器
extractor.enableImageTypeByName('LoG')
extractor.enableImageTypeByName('Wavelet')
特征类型特征解释
First Order （共19个）
能量
Total
总能量
熵
最小特征值
特征值的百分之10的值
特征值得百分之90的值
最大特征值
Mean
均值
中位数
四分位距离
Range
灰度值范围
n
(MAD)平均绝对误差
(rMAD) 鲁棒平均绝对偏差
(RMS)均方根误差
标准差。测量平均值的变化或离散量，默认不启用，因为与方差相关
偏度。测量值的分布关于平均值的不对称性
峰度。是图像 ROI 中值分布的“峰值”的量度
方差。是每个强度值与平均值的平方距离的平均值
均匀度。是每个强度值的平方和
特征类型特征解释
Shape (3D) （共17个）
网格体积
体素体积
表面积
表面积体积比
球形度
紧凑度1。衡量肿瘤形状相对于球形的紧密程度，由于与球形度、紧凑度2、球形不对称相关，默认不启用
紧凑度2。衡量肿瘤形状相对于球形的紧密程度，与1的计算公式不同，默认不启用
on
球形不对称。是肿瘤区域的表面积与与肿瘤区域具有相同体积的球体的表面积之比，是球形度的倒数，默认不启用
最大3D直径
ce
最大二维直径（切片），定义为行列（通常是轴向）平面中肿瘤表面网格顶点之间的最大成对欧几里得距离。
umn
最大2D直径（列）。定义为行切片（通常是冠状）平面中肿瘤表面网格顶点之间的最大成对欧几里得距离。
最大二维直径（行）。定义为列切片（通常为矢状）平面中肿瘤表面网格顶点之间的最大成对欧几里得距离。
长轴长度
短轴长度
最小轴长度
伸长率。表示ROI形状中两个最大的主成分之间的关系
平整度。表示ROI形状中最大和最小主成分之间的关系
特征类型特征解释
Shape (2D)
网格面。由三角形网格定义，计算ROI的表面
像素面。ROI 中的像素数乘以单个像素的表面积
周长
o
周长比。值越低越紧凑
球形度
on
球形不均匀度。是球形度的倒数 ，与球形度相关，默认不启用
最大二维直径。定义为肿瘤表面网格顶点之间的最大成对欧几里得距离
长轴长度
短轴长度
伸长率。显示了 ROI 形状中两个最大主成分之间的关系
特征类型特征解释
Gray Level Co- (GLCM) 灰度共生矩阵特征（共24个）
自相关。是纹理细度和粗糙度大小的量度
联合平均
集群显著性。衡量 GLCM 偏度和不对称性的指标
集群阴影。衡量 GLCM 偏度和均匀度的指标
集群趋势
对比度
相关性
平均差
差分熵
差异方差。是异质性的一种度量
差异性。已弃用
联合能量
联合熵。是邻域强度值的随机性/可变性的度量
Imc1
相关性信息度量 (IMC) 1
Imc2
相关性信息度量 (IMC) 2
Idm
逆差矩
MCC
最大相关系数
Idmn
逆差分矩归一化
Id
逆差
Idn
逆差分归一化
逆方差
最大概率
求和平均
总熵
平方和
特征类型特征解释
Gray Level Size Zone (GLSZM) 灰度级大小区域矩阵特征 （共16个）
是小尺寸区域分布的量度，数值越大表示区域越小，纹理越精细
是大面积区域分布的度量，值越大表示区域越大，纹理越粗糙。
ty
灰度非均匀性。测量图像中灰度级强度值的可变性，值越低表示强度值越均匀。
灰度非均匀归一化。测量图像中灰度强度值的可变性，值越低表示强度值的相似性越大
y
尺寸区域非均匀性。测量图像中大小区域体积的可变性，值越低表示大小区域体积越均匀
尺寸区域非均匀归一化。测量整个图像中大小区域体积的可变性，值越低表示图像中区域大小体积之间的同质性越高
区域百分比。获取 ROI 中的区域数和体素数的比率来测量纹理的粗糙度
灰度方差
区域方差
区域熵
asis
低灰度区域强调
hasis
高灰度区域强调
小面积低灰度区域强调
小区域高灰度强调
大面积低灰度强调
大面积高灰度强调
特征类型特征解释
Gray Level Run (GLRLM) 灰度游程矩阵 (GLRLM) 特征 （共16个）
是衡量短运行长度分布的指标，值越大表示运行长度越短，纹理纹理越精细。
是对长行程分布的度量，值越大表示行程越长，结构纹理越粗糙。
ty
灰度非均匀性
灰度非均匀归一化
ty
运行长度非均匀性
运行长度非均匀归一化
运行百分比。通过获取 ROI 中的运行次数和体素数的比率来测量纹理的粗糙度
灰度级方差
运行方差
运行熵
sis
测量低灰度值的分布，值越高表示图像中低灰度值的集中度越高
asis
测量较高灰度值的分布，值越高表示图像中高灰度值的集中度越高
测量具有较低灰度值的较短运行长度的联合分布
测量具有较高灰度值的较短运行长度的联合分布
测量具有较低灰度值的长运行长度的联合分布
测量具有较高灰度值的长运行长度的联合分布
特征类型特征解释
Gray Tone (NGTDM) 相邻灰度差矩阵（共5个） (NGTDM) 特征
粗糙度
对比度
繁忙度。从像素到其邻居的变化的度量。繁忙度的高值表示“繁忙”图像，像素与其邻域之间的强度快速变化
复杂度
强度
特征类型特征解释
Gray Level (GLDM) 灰度依赖矩阵 (GLDM) 特征（共14个）
sis
小依赖分布的度量，值越大表示依赖越小，纹理越不均匀
sis
大依赖分布的度量，值越大表示依赖越大，纹理越均匀
ty
灰度非均匀性。测量图像中灰度级强度值的相似性
ity
衡量整个图像中依赖关系的相似性，值越低表示图像中依赖关系之间的同质性越高
依赖非均匀归一化。衡量整个图像中依赖关系的相似性，值越低表示图像中依赖关系之间的同质性越高
灰度方差
依赖方差
依赖熵
测量低灰度值的分布，值越大表示图像中低灰度值的集中度越高
s
测量较高灰度值的分布，值越高表示图像中高灰度值的集中度越高
测量具有较低灰度值的小相关性的联合分布
测量具有较高灰度值的小相关性的联合分布
测量具有较低灰度值的大相关性的联合分布
测量具有较高灰度值的大相关性的联合分布
注意：除了形状特征类外，其他特征都可以在原始图像和滤波后的图像上进行计算
具体的指定方法：
可以直接使用：( )启用所有类型
也可以使用：Name(, =True)启用你想用的类型
例如：
# 设置一阶特征
extractor.enableFeatureClassByName('firstorder')
# 设置只提取一阶特征的'Mean'和'Skewness'
extractor.enableFeaturesByName(firstorder=['Mean', 'Skewness'])
也可以直接通过特征类别来启用具体的特征，示例见：
图像归一化：
图像/mask重采样：
例如：
settings = {}
settings['
'] = 25
settings['resampledPixelSpacing'] = [3,3,3] # [3,3,3] is an example for defining resampling (voxels with size 3x3x3mm)
settings['interpolator'] = sitk.sitkBSpline
还有其他的一些设置，根据自己的需要修改即可。
最后举一个完整的例子，我要处理的是CT数据，所以我参考了的源码中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 errorraise Exception('Error getting testcase!') # Raise exception to prevent cells below from running in case of "run all"settings = {}settings['binWidth'] = 25 # 5settings['sigma'] = [3, 5]settings['Interpolator'] = sitk.sitkBSplinesettings['resampledPixelSpacing'] = [1, 1, 1] # 3,3,3settings['voxelArrayShift'] = 1000 # 300settings['normalize'] = Truesettings['normalizeScale'] = 100extractor = 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=1if 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:],nameimage_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()
最后结果如下：
