SemanticFeature-Analysis/get_radiomic.py at main · hsu-lab/SemanticFeature-Analysis · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
#!/usr/bin/env python
# coding: utf-8
# %%
import numpy as np
import pandas as pd
import radiomics
import six
import SimpleITK as sitk
from radiomics import featureextractor
import sys
import os
import argparse
import nrrd
import pickle
import matplotlib.pyplot as plt
import nibabel as nib


# %%

def read_file(file_path):
#     print(file_path)
    img = nib.load(file_path)
    data = img.get_fdata()
    header = img.header
    return data, header

def get_radiomics(extractor, image, label):
    feats = {}
    featureVector = extractor.execute(image, label)
    for (key, val) in six.iteritems(featureVector):
        if key.startswith("original") or key.startswith("log") or key.startswith("wavelet") :
            feats[key] = val
    return feats

# (-1000, 500) for NLST dataset, for UCLA 0-1500
def scale_HU(image_array, lower_bound=-1000, upper_bound=500):
    image_array[image_array > upper_bound] = upper_bound
    image_array[image_array < lower_bound] = lower_bound
    return image_array

def get_radiomics_feat(image_dir, masks, ids, scale_HU_value=True):
    all_nodules_feats_dict = {'pid':[]}
    for idx, mask_pth in enumerate(masks):
        pid = ids[idx]
#         print(pid)
#         print(image_dir[idx])
#         print(masks[idx])
        # read image
        data, header = read_file(image_dir[idx])
        if scale_HU_Value:
            data = scale_HU(data) # scale to 0-1500 range
        sitk_image = sitk.GetImageFromArray(data)
        # read mask
        mask = nib.load(mask_pth)
        mask_data = mask.get_fdata()
        mask_header = mask.header
        assert data.shape == mask_data.shape
        sitk_mask = sitk.GetImageFromArray(mask_data)

        # extract features
        feats = get_radiomics(extractor, sitk_image, sitk_mask)
        feats_values = [float(x) for x in feats.values()]

        # write header
        if idx == 0:
            for header in feats.keys():
                all_nodules_feats_dict[header] = []
        # write values
        for feat_key in feats:
            all_nodules_feats_dict[feat_key].append(float(feats[feat_key]))
        all_nodules_feats_dict['pid'].append(pid)
        print(f"PID:{pid} complete!")
    return all_nodules_feats_dict

if __name__ == '__main__':
    # setting
    parser = argparse.ArgumentParser()
    parser.add_argument('--input_csv', type=str, help = "Input csv file path")
    parser.add_argument('--output_csv', type=str, help = "Output csv file path")

    args = parser.parse_args()
    input_csv = args.input_csv
    output_dset = args.output_csv
    if not os.path.exists('/'.join(output_dset.split('/')[:-1])):
        os.makedirs('/'.join(output_dset.split('/')[:-1]))

    # additional options
    scale_HU_Value = True
    paramPath = './CT.yaml'

    # read radiomics configuration
#     print('Parameter file, absolute path:', os.path.abspath(paramPath))
    extractor = featureextractor.RadiomicsFeatureExtractor(paramPath)

    input_df = pd.read_csv(input_csv)
    pids = input_df["PID"].values
    print('num of cases:', len(pids))

#     for idx, cnd in enumerate(cnds):
#     print(f"counter:{idx}")
    imgs = input_df["img_path"].values
    masks = input_df["seg_path"].values
#         print('Reading cnd:', path_to_read)
    df_dict = get_radiomics_feat(imgs, masks, pids, scale_HU_Value)
    df_dict = pd.DataFrame(df_dict)
    df_dict = df_dict.rename(columns = {'pid':'PID'})
    df_dict.to_csv(output_dset, index=False)

    print(f"Feature extraction complete. Saved feature dataset to {output_dset}")


# %%