使用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)
'] = 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()

最后结果如下：