dicom_utils.py

import os
import pydicom
import numpy as np
import scipy

# Read in a DICOM file.
def read_dcm(file_name):
  
  dcm = pydicom.dcmread(file_name)
  return dcm

def find_dcm_files(path):

  files = []
  
  for r, d, f in os.walk(path):
    if 'mask' not in r: # exclude masks
        for file in f:
          if '.dcm' in file:
            files.append(os.path.join(r, file))
  
  files = sorted(files)
  print(files)
  return files

def read_dcm_archive(path):
  
  filenames = find_dcm_files(path)
  print('Reading %d files in %s' % (len(filenames), path))
  slices = [read_dcm(filename) for filename in filenames]
  print('Finished reading DICOM archives.')
  
  # Get rid of "cover" images generated by PACS, which 
  # don't belong in the current series. 
  # This is typically a single coronal slice at the 
  # 1st position in the axial series.
  
  if len(slices) > 2:
    slice_1, slice_2 = slices[0], slices[1]
    slice_1_iop = slice_1[1]['image_orientation_patient']
    slice_2_iop = slice_2[1]['image_orientation_patient']
    if slice_1_iop != slice_2_iop:
      files.remove(slice_1)
  
  # Get the first slice and length
  first_slice_data, first_slice_attrs = slices[0]
  print('Slice thickness: %d mm' % first_slice_attrs['slice_thickness'])
  
  image_orientation_patient = first_slice_attrs['image_orientation_patient']
  patient_position = first_slice_attrs['patient_position']
  
  # Calculate Z vector by cross product of X and Y vectors
  u = [
    float(image_orientation_patient[0]),
    float(image_orientation_patient[1]),
    float(image_orientation_patient[2])
  ]
  
  v = [
    float(image_orientation_patient[3]),
    float(image_orientation_patient[4]),
    float(image_orientation_patient[5])
  ]
  
  ab = [u[1]*v[2] - u[2]*v[1], 
        u[2]*v[0] - u[0]*v[2], 
        u[0]*v[1] - u[1]*v[0]]
  
  # Parse out imagePositionPatient for each slice
  ipp = []
  
  for slice_n in slices:
    slice_attrs = slice_n[1]
    image_position_patient = slice_attrs['image_position_patient']
    ipp.append([float(image_position_patient[0]), 
                float(image_position_patient[1]), 
                float(image_position_patient[2])])

  # Sort slices by scalar product between vector and position
  # which gives the distance along the Z vector
  distances = []
  slice_nums = []
  
  for i in range(0, len(slices)):
    distances.append(ipp[i][0] * ab[0] + 
    					       ipp[i][1] * ab[1] +
    						     ipp[i][2] * ab[2])
    slice_nums.append(i)
  
  slice_nums.sort(key=lambda x: distances[x])
  sorted_slices = [slices[i] for i in slice_nums]
  
  if patient_position == "FFDR" or \
     patient_position == "FFDL" or \
     patient_position == "FFP":
     sorted_slices = reversed(sorted_slices)
	
  series = [np.asarray(s[0]) for s in sorted_slices]
  
  return series

def window_level(arr, window_center, window_width, lut_min=0, lut_max=255):

    # Basic sanity checking.
    if np.isreal(arr).sum() != arr.size: raise ValueError
    if lut_max != 255: raise ValueError
    if arr.dtype != np.float64: arr = arr.astype(np.float64)

    # Get window information.
    window_width = max(1, window_width)
    wc, ww = np.float64(window_center), np.float64(window_width)
    lut_range = np.float64(lut_max) - lut_min

    # Transform the image.
    minval = wc - 0.5 - (ww - 1.0) / 2.0
    maxval = wc - 0.5 + (ww - 1.0) / 2.0
    min_mask = (minval >= arr)
    to_scale = (arr > minval) & (arr < maxval)
    max_mask = (arr >= maxval)
    if min_mask.any(): arr[min_mask] = lut_min

    # Scale the image to the right proportions.
    if to_scale.any(): arr[to_scale] = \
      ((arr[to_scale] - (wc - 0.5)) /
      (ww - 1.0) + 0.5) * lut_range + lut_min
    if max_mask.any(): arr[max_mask] = lut_max

    arr = np.rint(arr).astype(np.uint8)

    return arr

def resample_image(img, pixel_spacing):
    
    new_spacing = [1, 1]
    spacing = np.array(pixel_spacing, dtype=np.float32)

    resize_factor = spacing / new_spacing
    new_real_shape = img.shape * resize_factor
    new_shape = np.round(new_real_shape)
    real_resize_factor = new_shape / img.shape
    new_spacing = spacing / real_resize_factor

    img_resampled = scipy.ndimage.interpolation.zoom(img, real_resize_factor)
    
    #print('Before resampling: %dx%d' % (img.shape[1], img.shape[0]))
    #print('After resampling: %dx%d' % (img_resampled.shape[1], img_resampled.shape[0]))
    
    return img_resampled

def resample_image_3d(itk_image, out_spacing=[2.0, 2.0, 2.0], is_label=False):
    
    # Resample images to 2mm spacing with SimpleITK
    original_spacing = itk_image.GetSpacing()
    original_size = itk_image.GetSize()

    out_size = [
        int(np.round(original_size[0] * (original_spacing[0] / out_spacing[0]))),
        int(np.round(original_size[1] * (original_spacing[1] / out_spacing[1]))),
        int(np.round(original_size[2] * (original_spacing[2] / out_spacing[2])))]

    resample = sitk.ResampleImageFilter()
    resample.SetOutputSpacing(out_spacing)
    resample.SetSize(out_size)
    resample.SetOutputDirection(itk_image.GetDirection())
    resample.SetOutputOrigin(itk_image.GetOrigin())
    resample.SetTransform(sitk.Transform())
    resample.SetDefaultPixelValue(itk_image.GetPixelIDValue())

    if is_label:
        resample.SetInterpolator(sitk.sitkNearestNeighbor)
    else:
        resample.SetInterpolator(sitk.sitkBSpline)

    return resample.Execute(itk_image)

def read_dcm_as_grayscale(filename):

  pixel_arr, attrs = read_dcm(filename)

  return pixel_arr

def read_ct_dcm_as_grayscale(filename, window_center=None, window_width=None, resample=True):
  
  pixel_arr, attrs = read_dcm(filename)
  
  if (window_center is None) or (window_width is None):
    
    window_center = (pixel_arr.max() + pixel_arr.min()) / 2.0
    window_width = pixel_arr.max() - pixel_arr.min() + 1.0
  
    if ('window_center' in attrs) and ('window_width' in attrs):
      window_center = attrs['window_center']
      window_width = attrs['window_width']
      try: window_center = window_center[0]
      except: pass
      try: window_width = window_width[0]
      except: pass
      window_center = int(window_center)
      window_width = int(window_width)
  
  img = window_level(pixel_arr, window_center, window_width)
  
  if resample:
    img = resample_image(img, attrs['pixel_spacing'])
  
  return img

def read_dcm_as_mask(filename, resample=True):
  
  pixel_arr, attrs = read_dcm(filename)
  
  is_empty = np.sum(pixel_arr) == 0
  
  if resample:
    pixel_arr = resample_image(pixel_arr, attrs['pixel_spacing'])
  
  threshold = threshold_mean(pixel_arr)
  
  if is_empty: return np.zeros(pixel_arr.shape).astype(np.uint8)
  
  return (pixel_arr > threshold).astype(np.uint8)