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

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

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

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

第一步、安装

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

pip install pyradiomics

第二步、提取特征

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

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

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

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

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

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

具体的指定方法：
可以直接使用：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)'] = 25settings['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_functionimport siximport os  # needed navigate the system to get the input dataimport numpy as npimport radiomicsimport pandas as pdimport SimpleITK as sitkfrom radiomics import featureextractor  # This module is used for interaction with pyradiomicsimport argparsedef 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:],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 = idname_df.columns = namewriter = pd.ExcelWriter('NSCLC-Radiomics-features.xlsx')name_df.to_excel(writer)writer.save()

最后结果如下：

来源地址：https://blog.csdn.net/weixin_46428351/article/details/123592586