Run_IBSI_Benchmark.py

#!/usr/bin/env python
import logging
import os

import numpy as np
import pandas as pd
import SimpleITK as sitk
import six

import radiomics
from radiomics import setVerbosity, featureextractor

if not os.path.isdir("results"):
  os.mkdir("results")

CASES = six.moves.range(1, 6)
TYPES = ['', '_Combined']
IBSI_BINNING = False
IBSI_RESAMPLING = False
GRAY_VALUE_ROUNDING = True

#"""
rLogger = radiomics.logger
logHandler = logging.FileHandler(filename='results/IBSIlog.log', mode='w')
logHandler.setLevel(logging.DEBUG)
logHandler.setFormatter(logging.Formatter('%(levelname)-.1s: %(name)s: %(message)s'))

rLogger.setLevel(logging.DEBUG)
rLogger.addHandler(logHandler)

setVerbosity(logging.INFO)
# """

ibsiLogger = logging.getLogger('radiomics.ibsi')


def run_phantom():
  global TYPES
  ibsiLogger.info('################################### Extracting Phantom #############################')
  extractor = featureextractor.RadiomicsFeaturesExtractor()

  image = sitk.ReadImage(os.path.join('data', 'digital_phantom', 'Phantom.nrrd'))
  mask = sitk.ReadImage(os.path.join('data', 'digital_phantom', 'Phantom-label.nrrd'))

  result_series = pd.Series()
  extraction_mode = ('2D', '3D')
  for e in extraction_mode:
    ibsiLogger.info('######################### MODE %s ####################' % e)
    for t in TYPES:
      ibsiLogger.info('######################### TYPE %s ####################' % t)
      params = os.path.join('configuration', 'Phantom%s.yml' % t)
      if not os.path.isfile(params):
        continue
      extractor.loadParams(params)
      extractor.addProvenance(provenance_on=(t == ''))
      extractor.settings['force2D'] = e == '2D'

      fv = pd.Series(extractor.execute(image, mask))

      if t != '':
        fv = fv.add_prefix(t[1:] + '_')

      if extractor.settings.get('force2D', False):
        fv = fv.add_prefix('2D_')
      else:
        fv = fv.add_prefix('3D_')
      result_series = result_series.append(fv)

    result_series.name = 'phantom'
    #result_series.to_csv('resCase%d.csv' % case_idx)  # Uncomment to enable saving intermediate results
  return result_series


def run_case(case_idx, image, mask):
  global TYPES, GRAY_VALUE_ROUNDING, ibsiLogger

  ibsiLogger.info('################################### Extracting Case %d #############################' % case_idx)
  extractor = featureextractor.RadiomicsFeaturesExtractor()

  result_series = pd.Series()
  for t in TYPES:
    ibsiLogger.info('######################### TYPE %s ####################' % t)

    params = os.path.join('configuration', 'case%d%s.yml' % (case_idx, t))
    if not os.path.isfile(params):
      continue

    extractor.loadParams(params)
    if GRAY_VALUE_ROUNDING:
      extractor.settings['grayValuePrecision'] = 0  # round to nearest integer when using IBSI resampling
    extractor.addProvenance(provenance_on=(t == ''))

    fv = pd.Series(extractor.execute(image, mask))

    if t != '':
      fv = fv.add_prefix(t[1:] + '_')

    if extractor.settings.get('force2D', False):
      fv = fv.add_prefix('2D_')
    else:
      fv = fv.add_prefix('3D_')

    result_series = result_series.append(fv)

  result_series.name = case_idx
  return result_series


def IBSI_binning(parameterValues, **kwargs):
  global ibsiLogger
  binWidth = kwargs.get('binWidth', 25)
  binCount = kwargs.get('binCount')
  resegmentRange = kwargs.get('resegmentRange')
  resegmentMode = kwargs.get('resegmentMode', 'absolute')

  if binCount is not None:
    binEdges = np.histogram(parameterValues, binCount)[1]
    binEdges[-1] += 1  # Ensures that the maximum value is included in the topmost bin when using numpy.digitize
  else:
    minimum = min(parameterValues)
    maximum = max(parameterValues)

    # Start binning form the first value lesser than or equal to the minimum value and evenly dividable by binwidth
    lowBound = minimum - (minimum % binWidth)
    # Add + binwidth to ensure the maximum value is included in the range generated by numpy.arange
    highBound = maximum + binWidth

    # #####################################
    # IBSI difference
    # #####################################
    if resegmentRange is not None:
      if resegmentMode == 'absolute':
        lowBound = min(resegmentRange)
      elif resegmentMode == 'sigma':
        lowBound = minimum
    # #####################################

    binEdges = np.arange(lowBound, highBound, binWidth)

    # if min(parameterValues) % binWidth = 0 and min(parameterValues) = max(parameterValues), binEdges will only contain
    # 1 value. If this is the case (flat region) ensure that numpy.histogram creates 1 bin (requires 2 edges). For
    # numpy.histogram, a binCount (1) would also suffice, however, this is not accepted by numpy.digitize, which also uses
    # binEdges calculated by this function.
    if len(binEdges) == 1:  # Flat region, ensure that there is 1 bin
      binEdges = [binEdges[0] - .5, binEdges[0] + .5]  # Simulates binEdges returned by numpy.histogram if bins = 1

    ibsiLogger.debug('Calculated %d bins for bin width %g with edges: %s)', len(binEdges) - 1, binWidth, binEdges)

  return binEdges  # numpy.histogram(parameterValues, bins=binedges)


def IBSI_resampling(image, mask, **kwargs):
  global ibsiLogger
  # resample image to new spacing, align centers of both resampling grids.
  spacing = kwargs.get('resampledPixelSpacing')
  grayValuePrecision = kwargs.get('grayValuePrecision')
  interpolator = kwargs.get('interpolator', sitk.sitkLinear)

  try:
    if isinstance(interpolator, six.string_types):
      interpolator = getattr(sitk, interpolator)
  except Exception:
    ibsiLogger.warning('interpolator "%s" not recognized, using sitkLinear', interpolator)
    interpolator = sitk.sitkLinear

  im_spacing = np.array(image.GetSpacing(), dtype='float')
  im_size = np.array(image.GetSize(), dtype='float')

  spacing = np.where(np.array(spacing) == 0, im_spacing, spacing)

  spacingRatio = im_spacing / spacing
  newSize = np.ceil(im_size * spacingRatio)

  # Calculate center in real-world coordinates
  im_center = image.TransformContinuousIndexToPhysicalPoint((im_size - 1) / 2)

  new_origin = tuple(np.array(image.GetOrigin()) + 0.5 * ((im_size - 1) * im_spacing - (newSize - 1) * spacing))

  ibsiLogger.info('Resampling from %s to %s (size %s to %s), aligning Centers', im_spacing, spacing, im_size, newSize)

  rif = sitk.ResampleImageFilter()
  rif.SetOutputOrigin(new_origin)
  rif.SetSize(np.array(newSize, dtype='int').tolist())
  rif.SetOutputDirection(image.GetDirection())
  rif.SetOutputSpacing(spacing)

  rif.SetOutputPixelType(sitk.sitkFloat32)
  rif.SetInterpolator(interpolator)
  res_im = rif.Execute(sitk.Cast(image, sitk.sitkFloat32))

  # Round to n decimals (0 = to nearest integer)
  if grayValuePrecision is not None:
    ibsiLogger.debug('Rounding Image Gray values to %d decimals', grayValuePrecision)
    im_arr = sitk.GetArrayFromImage(res_im)
    im_arr = np.round(im_arr, grayValuePrecision)
    round_im = sitk.GetImageFromArray(im_arr)
    round_im.CopyInformation(res_im)
    res_im = round_im

  # Sanity check: Compare Centers!
  new_center = res_im.TransformContinuousIndexToPhysicalPoint((newSize - 1) / 2)
  ibsiLogger.debug("diff centers: %s" % np.abs(np.array(im_center) - np.array(new_center)))

  rif.SetOutputPixelType(sitk.sitkFloat32)
  rif.SetInterpolator(sitk.sitkLinear)
  res_ma = rif.Execute(sitk.Cast(mask, sitk.sitkFloat32))
  res_ma = sitk.BinaryThreshold(res_ma, lowerThreshold=0.5)

  return res_im, res_ma


def index_func(series, *args, **kwargs):
  if 'idx' not in series or np.isnan(series['idx']):
    return series

  idx_loc = series.index.get_loc('idx')

  vals = []
  for val in series.iloc[idx_loc+1:]:
    vals.append(val)

  for val_idx, value in enumerate(vals):
    try:
      # new_val = eval(value)
      series.values[int(idx_loc + 1 + val_idx)] = value[int(series['idx'])]
    except:
      pass
  return series


def correct_kurtosis(series, *args, **kwargs):
  if 'pyradiomics_feature' not in series:
    return series

  if 'firstorder_Kurtosis' not in str(series['pyradiomics_feature']):
    return series

  idx_loc = series.index.get_loc('pyradiomics_feature')

  if 'idx' in series:
    idx_loc += 1

  vals = []
  for val in series.iloc[idx_loc+1:]:
    vals.append(val)

  for val_idx, value in enumerate(vals):
    try:
      # new_val = eval(value)
      series.values[int(idx_loc + 1 + val_idx)] = value - 3
    except:
      pass
  return series


if __name__ == '__main__':
  if IBSI_BINNING:
    radiomics.imageoperations.getBinEdges = IBSI_binning
  if IBSI_RESAMPLING:
    radiomics.imageoperations.resampleImage = IBSI_resampling

  mapping_file = os.path.join('mapping', 'mapping_phantom.csv')
  mapping = pd.read_csv(mapping_file)

  results = mapping.join(run_phantom(), on='pyradiomics_feature', how='left')
  results = results.apply(correct_kurtosis, axis=1)

  results.sort_index(inplace=True)
  results.to_csv('results/results_phantom.csv')

  im = sitk.ReadImage(os.path.join('data', 'patient_cases', 'PAT1', 'PAT1.nrrd'))
  ma = sitk.ReadImage(os.path.join('data', 'patient_cases', 'PAT1', 'GTV.nrrd'))

  for case in CASES:
    mapping_file = os.path.join('mapping', 'mapping_case%s.csv' % case)
    mapping = pd.read_csv(mapping_file)

    results = mapping.join(run_case(case, im, ma), on='pyradiomics_feature', how='left')
    results = results.apply(index_func, axis=1)
    results = results.apply(correct_kurtosis, axis=1)

    results.sort_index(inplace=True)
    results.to_csv('results/results_case%s.csv' % case)
  exit(0)