Dicom_parser.py

# -*- coding: utf-8 -*-
"""Copy of Dicom_Parser.ipynb

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/1zi5G9sqmvT0mTafxmcgr2MTJXbeDx-yv
"""

!nvidia-smi

!lscpu

!pip install plotly==4.11.0

import os

from google.colab import drive
drive.mount('/content/drive')

os.chdir('/content/drive/My Drive/Dicom_Leg')
!ls

#!unzip A-20201011T051721Z-001.zip

!pip install pydicom

!pip install chart_studio

# common packages 
import numpy as np 
import os
import copy
from math import *
import matplotlib.pyplot as plt
from functools import reduce
# reading in dicom files
import pydicom
# skimage image processing packages
from skimage import measure, morphology
from skimage.morphology import ball, binary_closing
from skimage.measure import label, regionprops
# scipy linear algebra functions 
from scipy.linalg import norm
import scipy.ndimage
# ipywidgets for some interactive plots
from ipywidgets.widgets import * 
import ipywidgets as widgets
# plotly 3D interactive graphs 
import plotly
from plotly.graph_objs import *
import chart_studio.plotly as py
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
from plotly.tools import FigureFactory as FF
from plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot

def load_scan(path):
  slices = []
  for root, dirs, files in os.walk(path):
    path = root.split(os.sep)
    for file in files:
        slices.append(pydicom.dcmread('/'.join(path) + '/' + file))
  
  slices = [s for s in slices if 'SliceLocation' in s]
  slices.sort(key = lambda x: int(x.InstanceNumber))
  try:
    slice_thickness = np.abs(slices[0].ImagePositionPatient[2] - slices[1].ImagePositionPatient[2])
  except:
    slice_thickness = np.abs(slices[0].SliceLocation - slices[1].SliceLocation)
  for s in slices:
      s.SliceThickness = slice_thickness
  return slices

def get_pixels_hu(scans):
  print(len(scans))
  image = np.stack([s.pixel_array for s in scans if len(s.pixel_array) == 512])
  image = image.astype(np.int16)
  # Set outside-of-scan pixels to 0
  # The intercept is usually -1024, so air is approximately 0
  image[image == -2000] = 0
  
  # Convert to Hounsfield units (HU)
  intercept = scans[0].RescaleIntercept
  slope = scans[0].RescaleSlope
  
  if slope != 1:
      image = slope * image.astype(np.float64)
      image = image.astype(np.int16)
      
  image += np.int16(intercept)
  
  return np.array(image, dtype=np.int16)

# set path and load files 
path = './A/A/B'
path1 = './A/A/A'
path2= './A/A/C'

print('DEBUG : Loading Files.')
patient_dicom = load_scan(path)
print('DEBUG : Files Loaded.')

from collections import Counter
Counter([len(s.pixel_array) for s in patient_dicom])

patient_pixels = get_pixels_hu(patient_dicom)

file_used= patient_pixels
imgs_to_process = file_used.astype(np.float64) 

plt.hist(imgs_to_process.flatten(), bins=50, color='c')
plt.xlabel("Hounsfield Units (HU)")
plt.ylabel("Frequency")
plt.show()

id = 0
imgs_to_process = patient_pixels

def sample_stack(stack, rows=6, cols=6, start_with=10, show_every=3):
    fig,ax = plt.subplots(rows,cols,figsize=[12,12])
    for i in range(rows*cols):
        ind = start_with + i*show_every
        ax[int(i/rows),int(i % rows)].set_title('slice %d' % ind)
        ax[int(i/rows),int(i % rows)].imshow(stack[ind],cmap='gray')
        ax[int(i/rows),int(i % rows)].axis('off')
    plt.show()

sample_stack(imgs_to_process)

print(f"Slice Thickness:{patient_dicom[0].SliceThickness}")
print(f"Pixel Spacing (row, col): ({patient_dicom[0].PixelSpacing[0]}, {patient_dicom[0].PixelSpacing[1]}). ")

id = 0
imgs_to_process = patient_pixels
def resample(image, scan, new_spacing=[1,1,1]):
  # Determine current pixel spacing
  spacing = map(float, ([scan[0].SliceThickness] + list(scan[0].PixelSpacing)))
  spacing = np.array(list(spacing))

  resize_factor = spacing / new_spacing
  new_real_shape = image.shape * resize_factor
  new_shape = np.round(new_real_shape)
  real_resize_factor = new_shape / image.shape
  new_spacing = spacing / real_resize_factor
  
  image = scipy.ndimage.interpolation.zoom(image, real_resize_factor)
  
  return image, new_spacing

print("Shape before resampling\t", imgs_to_process.shape)
imgs_after_resamp, spacing = resample(imgs_to_process, patient_dicom, [1,1,1])
#imgs_after_resamp = imgs_to_process
print("Shape after resampling\t", imgs_after_resamp.shape)

def make_mesh(image, threshold=226, step_size=1):

  print("Transposing surface")
  p = image.transpose(2,1,0)
  
  print("Calculating surface")
  verts, faces, norm, val = measure.marching_cubes_lewiner(p, threshold, step_size=step_size, allow_degenerate=True) 
  return verts, faces
import plotly.graph_objects as go
def plotly_3d(verts, faces):
  x,y,z = zip(*verts) 
  
  print("Drawing")
  
  # Make the colormap single color since the axes are positional not intensity. 
#    colormap=['rgb(255,105,180)','rgb(255,255,51)','rgb(0,191,255)']
  colormap=['rgb(236, 236, 212)','rgb(236, 236, 212)']
  
  fig = FF.create_trisurf(x=x,
                      y=y, 
                      z=z, 
                      plot_edges=False,
                      colormap=colormap,
                      simplices=faces,
                      backgroundcolor='rgb(64, 64, 64)',
                      title="Interactive Visualization")
  iplot(fig)

  #fig = go.Figure(data=go.Isosurface(
  #    x=x,
  #    y=y,
  #    z=z,
  #    value=faces,
  #    isomin=10,
  #    isomax=50,
  #    surface_count=5, # number of isosurfaces, 2 by default: only min and max
  #    colorbar_nticks=5, # colorbar ticks correspond to isosurface values
  #    caps=dict(x_show=False, y_show=False)
  #    ))
  #fig.show()

def plt_3d(verts, faces):
  print("Drawing")
  x,y,z = zip(*verts) 
  fig = plt.figure(figsize=(10, 10))
  ax = fig.add_subplot(111, projection='3d')

  # Fancy indexing: `verts[faces]` to generate a collection of triangles
  mesh = Poly3DCollection(verts[faces], linewidths=0.05, alpha=1)
  face_color = [1, 1, 0.9]
  mesh.set_facecolor(face_color)
  ax.add_collection3d(mesh)

  ax.set_xlim(0, max(x))
  ax.set_ylim(0, max(y))
  ax.set_zlim(0, max(z))
  ax.set_fc((0.7, 0.7, 0.7))
  plt.show()

# v, f = make_mesh(imgs_after_resamp, 300)
# plt_3d(v, f)
# v, f = make_mesh(imgs_after_resamp,250,2)
# plotly_3d(v, f)

v, f = make_mesh(imgs_after_resamp,226,2)

plotly_3d(v, f)

pip install numpy-stl

import numpy as np
from stl import mesh

body = mesh.Mesh(np.zeros(f.shape[0], dtype=mesh.Mesh.dtype))
for i, fc in enumerate(f):
    for j in range(3):
        body.vectors[i][j] = v[fc[j],:]

body.save('body.stl')

import csv
file = open('f.csv','w')
writer = csv.writer(file, delimiter=',')
writer.writerows(f)
file.close()

#v.to_csv('data_v.csv')

"""## Refs
- https://www.raddq.com/dicom-processing-segmentation-visualization-in-python/
- https://medium.com/@hengloose/a-comprehensive-starter-guide-to-visualizing-and-analyzing-dicom-images-in-python-7a8430fcb7ed
- https://plotly.com/python/3d-isosurface-plots/
- https://pypi.org/project/numpy-stl/
- https://vtk.org/vtk-in-action/#image-gallery
"""