OneStepDataAnalysis/RadiomicsHistograms.py at master · ImageGuidedTherapyLab/OneStepDataAnalysis · 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
# uses pyradiomics to generate feature maps and prints a histogram of the data values

# Requires pyradiomics version 2.0 or greater

# example usage: python RadiomicsHistograms.py --file BRATSlist.csv --param exampleVoxel.yaml --writeImage

from radiomics import featureextractor
import matplotlib.pyplot as plt
import six
import sys, os, argparse, csv
import numpy as np

import SimpleITK as sitk


# Get params from cmd line
parser = argparse.ArgumentParser(description = "Plot histograms of radiomics feature" +
    "for multiple images. Optionally save feature maps of the form imageName_maskName_labelNumber_featureName.nii.gz. Feature maps are cropped and need to be put back into the image space with c3d")
parser.add_argument("--file", "-f", help = "CSV with filepaths and columns Image, Mask, Label",
      required = True)
#parser.add_argument("--feature", default = None,
#      help = "feature name: i.e. original_firstorder_Mean. See pyradiomics -h for help")
parser.add_argument("--writeImage", "-w", dest="writeImage", action="store_true",
      help = "Include flag to write out a feature map (NIFTI) for each feature for each row. Cropped to the specific ROI")
parser.add_argument("--directory", "-d", default = "",
      help = "Directory to place output files in")
parser.add_argument("--param", "-p", default = "",
      help = "Pyradiomics parameter file to use for feature extraction. Old parameter files will need the VoxelSettings category added in order to be compatible.")

args = parser.parse_args()

# Check that csv file exists
if not os.path.isfile(args.file):
  sys.exit("Error: csv file " + args.file + " not found")

#Read file
inputCsv = csv.DictReader(open(args.file))
if 'Image' not in inputCsv.fieldnames:
  sys.exit("CSV file does not contain 'Image' column")
elif 'Mask' not in inputCsv.fieldnames:
  sys.exit("CSV file does not contain 'Mask' column")
elif 'Label' not in inputCsv.fieldnames:
  sys.exit("CSV file does not contain 'Label' column")

#set up extractor
if not os.path.isfile(args.param):
  sys.exit("No parameter file given")

extractor = featureextractor.RadiomicsFeaturesExtractor(args.param)


#set up flags to add to histogram
newHist  = True
oldImage = False
oldMask  = False

for row in inputCsv:
  imageName = row['Image']
  maskName  = row['Mask']
  label     = int(row['Label'])

  #If same image/mask combination, add to histogram
  if (imageName != oldImage) | (maskName != oldMask):
    newHist = True
    print("Saving new histogram")

  oldImage = imageName
  oldMask  = maskName

  # compute feature maps
  print("Extracting features for:\n" + imageName + "\n" + maskName + "\n" +  "label = " + str(label))
  result = extractor.execute(imageName, maskName, label, voxelBased=True)

  #HACK if more than 1 feature disable plotting on sasme axes since key,val loop over features first
  if len(result.keys()) != 1:
    print("Multiple features being output, disabling overlaid histograms")
    newHist = True

  # extract images
  for key, val in six.iteritems(result):
    outname = args.directory + '_'.join([os.path.basename(imageName).rsplit('.')[0],
                          os.path.basename(maskName).rsplit('.')[0],
                          str(label), key]) + '.nii.gz'

    print("Creating histogram for feature " + key)
    voxelArray = sitk.GetArrayFromImage(val)
    voxelArray = voxelArray[~np.isnan(voxelArray)] #only non-nan voxels
    voxelMean = voxelArray.mean()
    print("Mean: " + str(voxelMean))

    #Print histogram, if different image/mask combo then start new figure
    if newHist:
      plt.figure()
      newHist = False

    n, bins, patched = plt.hist(voxelArray,25,density=True, alpha=0.75, label = str(label))
    plt.xlabel(key)
    plt.ylabel('density')
    plt.title(imageName +"\n"+ maskName +"\n"+str(label))
    plt.grid(True)
    plt.legend()
    plt.savefig(outname.replace('.nii.gz','.png'), bbox_inches='tight')

    if args.writeImage:
      print('Writing Image ' + outname)
      sitk.WriteImage(val, outname)


# Check that feature name is valid


# featureExtractor.enableFeaturesByName()
## key is the feature class name, value is a list of enabled feature names