PSOAS YALLAH IT WORKS

gistfile1.txt

#!/usr/bin/python
import sys, glob
import numpy as np
import os
import matplotlib.pyplot as plt
from matplotlib import colors
from scipy.signal import convolve2d
import scipy
from scipy import ndimage
from scipy.misc import imread

from skimage import io
import morphsnakes
from scipy.signal import decimate
from sklearn.metrics import jaccard_similarity_score
  

from scipy.spatial import Delaunay
from skimage.draw import circle, polygon, ellipse
from scipy.stats import itemfreq
from skimage import measure
from skimage.morphology import convex_hull_image
from scipy.ndimage import binary_fill_holes, map_coordinates, measurements
from scipy.misc import imsave
from scipy.ndimage.filters import median_filter
from skimage import measure
import math

import shapely.geometry as geometry

from shapely.geometry import MultiLineString


from descartes import PolygonPatch
from shapely.geometry import MultiLineString
from shapely.ops import cascaded_union, polygonize

import cv2

from dicom_utils import read_dcm
import random


def color_overlay(ax, mask, color):

  transparent = (0.0, 0.0, 0.0, 0.0)

  cmap = colors.ListedColormap([transparent, color])

  bounds=[0.5]
  norm = colors.BoundaryNorm(bounds, cmap.N)

  ax.imshow(mask, cmap=cmap, norm=norm, alpha=0.5)

def image_overlay(im, layers, colors, show=False, file=None):

  fig, axes = plt.subplots(nrows=1)

  plt.imshow(im, cmap='gray')

  k = 0

  for layer in layers:
    color_overlay(plt, layer, colors[k])
    k += 1

  if show:
    plt.show()
  
  if file != None:
    plt.savefig(file)
  
    plt.close()
    
def show(im, cmap='gray'):
  plt.imshow(im, cmap=cmap)
  plt.show()

def image_in_dir(directory):
  filenames = os.listdir(directory)
  if filenames[0] == '.DS_Store':
    return directory  + '/' + filenames[1]
  else:
    return directory  + '/' + filenames[0]
  
def read_image(image_file):

  im = io.imread(image_file, as_grey=True)
  
  return im
  

def largest_object_mask(mask):

  labeled_mask, number_of_objects = ndimage.label(mask)
  
  labels = itemfreq(labeled_mask.flatten())[:,1]
  
  labels[labeled_mask[0,0]] = 0
  
  max_label = np.where(labels == labels.max())[0][0]
  
  return labeled_mask == max_label

def threshold(im_t):
  
  im = np.copy(im_t)
  im[im > 130] = 0
  im[im < 80] = 0
  
  return im

def align(im1, im2):
   
  # Find size of image1
  sz = im1.shape
 
  warp_matrix = np.eye(2, 3, dtype=np.float32)
 
  # Specify the number of iterations.
  number_of_iterations = 500;
 
  # Specify the threshold of the increment
  # in the correlation coefficient between two iterations
  termination_eps = 1e-10;
 
  # Define termination criteria
  criteria = (cv2.TERM_CRITERIA_EPS | cv2.TERM_CRITERIA_COUNT, number_of_iterations,  termination_eps)
 
  # Run the ECC algorithm. The results are stored in warp_matrix.
  (cc, warp_matrix) = cv2.findTransformECC (im1,im2,warp_matrix, cv2.MOTION_AFFINE, criteria)
 
  return warp_matrix

def warp_image(mask, image, warp_matrix):
  
  sz = mask.shape
    
  warped_mask = cv2.warpAffine(mask, warp_matrix, (sz[1],sz[0]), flags=cv2.INTER_LINEAR + cv2.WARP_INVERSE_MAP);
  warped_image = cv2.warpAffine(image, warp_matrix, (sz[1],sz[0]), flags=cv2.INTER_LINEAR + cv2.WARP_INVERSE_MAP);
  
  warped_image[warped_mask == 0] = np.min(image)
  
  return warped_image
  
def np_to_opencv_image(vis):
  # vis = np.zeros((512, 512), np.float32)
  return cv2.cvtColor(vis, cv2.COLOR_GRAY2BGR)
  
#real_dcm, dcm = read_dcm(input_file)

input_dcm_file = './exports/ADAMS_A/image.dcm'
input_dcm_image, input_dcm_converted = read_dcm(input_dcm_file)

# export DYLD_LIBRARY_PATH=/usr/lib/opencv
root = './exports/'

i = 0

# export DYLD_FALLBACK_LIBRARY_PATH=/Users/louis-mullie/Code/opencv/build/lib:$DYLD_FALLBACK_LIBRARY_PATH 
# export PYTHONPATH=/Users/louis-mullie/Code/opencv/build/lib/python2.7/site-packages:$PYTHONPATH 
from skimage.feature import blob_dog

FAT_LOWER_THRESHOLD = -190

def extract_body_mask(dcm_image):
  body_mask = dcm_image > FAT_LOWER_THRESHOLD
  body_mask = binary_fill_holes(body_mask)
  body_mask = largest_object_mask(body_mask)

  return body_mask.astype(np.uint8) * 255

def extract_mask(image, mask):
  
  tmp = np.copy(image)
  tmp[mask == 0] = 0
  
  return tmp.astype(np.uint8)

def extract_bone_mask(body_image):

  bone_mask = body_image > 120
  bone_mask = binary_fill_holes(bone_mask)
  bone_mask = largest_object_mask(bone_mask)
  
  return bone_mask
  
def extract_perispinal_mask(aligned_dcm):
  
  bone_mask = extract_bone_mask(np.copy(aligned_dcm))
  
  center = measurements.center_of_mass(bone_mask)
  
  mask = np.copy(aligned_dcm)
  mask[bone_mask == 1] = 0
  
  mask = mask[center[0] - 75:center[0] + 125,
              center[1] - 125:center[1] + 125]
  
  mask[mask < -30] = 0
  mask[mask > 150] = 0
  
  mask = mask > 0

  paraspinal_muscle_mask = np.zeros((512, 512))
  paraspinal_muscle_mask[center[0] - 75:center[0] + 125,
              center[1] - 125:center[1] + 125] = mask
  
  paraspinal_mask = median_filter(paraspinal_muscle_mask, size=1)
  
  
  paraspinal_image = aligned_dcm[center[0] - 75:center[0] + 125,
              center[1] - 125:center[1] + 125]
  
  return [center, paraspinal_image, paraspinal_mask.astype(np.uint8)]
  

def triangulate(mask, alpha=1000):

  contours = measure.find_contours(mask, 0.8)

  points = []

  for n, contour in enumerate(contours):

    for m in xrange(0, len(contour[:, 0])):
      y = contour[:, 0][m]
      x = contour[:, 1][m]
      points.append([y, x])
      
  points = np.asarray(points)

  tri = Delaunay(points)
  edges = set()
  edge_points = []

  def add_edge(i, j):
    if (i, j) in edges or (j, i) in edges: return
    edges.add( (i, j) )
    edge_points.append(points[ [i, j] ])

  for ia, ib, ic in tri.vertices:
    try:
      pa = points[ia]
      pb = points[ib]
      pc = points[ic]
      # Lengths of sides of triangle
      a = math.sqrt((pa[0]-pb[0])**2 + (pa[1]-pb[1])**2)
      b = math.sqrt((pb[0]-pc[0])**2 + (pb[1]-pc[1])**2)
      c = math.sqrt((pc[0]-pa[0])**2 + (pc[1]-pa[1])**2)
      # Semiperimeter of triangle
      s = (a + b + c)/2.0
      # Area of triangle by Heron's formula
      area = math.sqrt(s*(s-a)*(s-b)*(s-c))
      circum_r = a*b*c/(4.0*area)
      # Here's the radius filter.
      if circum_r < alpha:
          add_edge(ia, ib)
          add_edge(ib, ic)
          add_edge(ic, ia)
    except:
      pass
      #print('Tri error')

  m = MultiLineString(edge_points)
  triangles = list(polygonize(m))

  poly = PolygonPatch(cascaded_union(triangles), alpha=0.5)

  vertices = poly.get_path().vertices

  rr, cc = polygon(vertices[:,0], vertices[:,1])

  img = np.zeros(mask.shape)
  img[rr, cc] = 1

  return img

def extract_intensity_line(im, x0, x1, y0, y1, num=500):
  
  x, y = np.linspace(x0, x1, num), np.linspace(y0, y1, num)
  zi = scipy.ndimage.map_coordinates(np.transpose(im), np.vstack((x,y)))
  
  return x, y, zi

from skimage import exposure


def gaussian(x, mu, sig):
  return np.exp(-np.power(x - mu, 2.) / (2 * np.power(sig, 2.)))
 

def mexican_hat(t, sig=1):

  return np.absolute(2 / (math.sqrt(3*sig)*math.pi**0.25) * (1-t**2/sig**2) * math.exp(-t**2/(2*sig**2)))
 
from skimage.filter import canny
from skimage.measure import regionprops
def draw_contour_lines(im, contour, center, offset=0, skip_first_max=False):

  # p2, p98 = np.percentile(im, (2, 98))
  # img_rescale = exposure.rescale_intensity(im, in_range=(p2, p98))
  # im = img_rescale
  
  tick = 0

  points_out = []
  points_in = []

  mat = []
  units = []
  lengths = []
  bone_mask = extract_bone_mask(np.copy(im))
  bone_center = measurements.center_of_mass(bone_mask)
  bone_top = regionprops(bone_mask)[0].bbox[0]
  
  started = False
  
  contour = contour.tolist()
  cent = None
  
  #k = 0
  
  for x, y in contour:
    
    tick += 1
    
    do_plot = False
    
    if tick % 100 == 0:
      do_plot = True
    
    orig = np.asarray([[center[0]], [center[1]]])
    
    if tick % 2 == 0:
      
      start = (y, x)
      
      stop = center

      x0, y0 = start[0], start[1]
      x1, y1 = stop[0], stop[1]
      
      x, yy, z = extract_intensity_line(im, x0, x1, y0, y1)

      wt, ht = im.shape[1], im.shape[0]

      dx, dy = [x1 - x0], [((ht - y1) - (ht - y0)) * -1]
      unit = np.asarray([dx, dy])
      unit /= (-1 * np.linalg.norm(unit))
      length = np.sqrt((x1-x0)**2+(y1-y0)**2)
      
      units.append([unit, length])
      mat.append(z)

  orig_mat = np.copy(mat)
  mat = median_filter(mat, size=10)
  mat = np.asarray(mat)
  
  # Plot transformed image
  #plt.imshow(mat,cmap='gray')
  #plt.show()
  
  mat2 = np.zeros((mat.shape[0] + 100, mat.shape[1]))
  
  mat2[0:50,:] = mat[mat.shape[0]-50:,:]
  mat2[50:mat.shape[0]+50,:] = mat 
  mat2[mat.shape[0]+50:,:] = mat[0:50,:]
  
  mat = mat2
  
  units2 = [0 for i in xrange(0, len(units) + 100)]
  units2[50:len(units)+50] = units
  units2[0:50] = units[len(units)-50:]
  units2[len(units2)-50:] = units[0:50]

  
  x = []
  y = []
  
  for i in xrange(0, mat.shape[0]):
    for j in xrange(0, mat.shape[1]):
      if mat[i,j] > -53:
        y.append(j)
        break
    x.append(i)
  
  lowess = sm.nonparametric.lowess(y, x, frac=0.1)
  y = lowess[:, 1]

  
  vfls = y
  
  # Plot smoothed point
  # plt.plot(x, y, 'r')
  # plt.plot(x, vfls, 'b')
  # plt.show()
  
  # Plot outer points
  
  # plt.imshow(mat, cmap='gray')
  # 
  # for k in xrange(0, len(vfls)):
  #   plt.plot(vfls[k], k, 'go')
  # 
  # plt.show()
  
  points_out = np.copy(vfls)
  
  #######
  
  x = []
  y = []
  
  for i in xrange(0, mat.shape[0]):
    start = int(points_out[i])
    for j in xrange(int(start+25), mat.shape[1]):
      if mat[i,j] < -53:
        y.append(j)
        break
          
    x.append(i)
  
  #lowess = sm.nonparametric.lowess(y, x, frac=0.1)
  #vfls = lowess[:, 1]
  #vfls = ndimage.filters.percentile_filter(y, 90, size=20)#, size=size)
  
  vfls = y
  
  # Plot inner wall points
  
  #plt.imshow(mat, cmap='gray')
  #
  #for k in xrange(0, len(vfls)):
  #  plt.plot(vfls[k], k, 'go')
  #
  #plt.show()
  
  points_in = np.copy(vfls)

  mat = mat[50:len(mat)-50,:]
  
  points_in = points_in[50:len(points_in)-50]
  points_out = points_out[50:len(points_out)-50]
  
  return [orig, units, points_out, points_in]

def find_fat_boundary(im, body_mask, center=False, offset=0, skip_first_max=False):
  
  bone_mask = extract_bone_mask(np.copy(im))
  
  bone_center = measurements.center_of_mass(bone_mask)
  bone_center = [bone_center[0], bone_center[1] - 100]

  contour = measure.find_contours(body_mask, 0.8)[0]
  
  orig, units, points_out, points_in = draw_contour_lines(im, contour, bone_center, \
    offset, skip_first_max=skip_first_max)

  return [orig, units, points_out, points_in]

import statsmodels.api as sm

def find_outer_wall(im, body_mask, center=False):

  orig, units, points_out, points_in = find_fat_boundary(im, body_mask, center=center)

  #mask = np.zeros(im.shape)
  
  mask = np.zeros(im.shape)

  #plt.imshow(im, cmap='gray')

  i = 0
  
  for unit, length in units:
    p = orig + unit * (length - points_out[i] / 500.0 * length)
    mask[p[1][0], p[0][0]] = 1
    #plt.plot(p[0][0], p[1][0], 'go')
    i += 1
    
  #plt.show()

  mask_out = triangulate(mask, 100)
  
  # Show outer wall mask over image
  # image_overlay(im, [mask_out], ['red'], show=True)

  bone_mask = extract_bone_mask(np.copy(im))
  bone_center = measurements.center_of_mass(bone_mask)
  bone_top = regionprops(bone_mask)[0].bbox[0]
  
  mask = np.zeros(im.shape)

  plt.imshow(im, cmap='gray')

  i = 0
  for unit, length in units:
    
    p = orig + unit * (length - points_in[i] / 500.0 * length)
    diffx = (p[0][0] - bone_center[1])
    diffy = (p[1][0] - (bone_top-50))
    
    if np.absolute(diffx) < 10:
      p[1][0] = bone_top
    elif np.absolute(diffx) < 100 and diffy > 0:
      p[1][0] = (bone_top-20) * 1.0/gaussian(diffx, 0, 200) + np.absolute(diffx)**0.75
    
    mask[p[1][0], p[0][0]] = 1
    plt.plot(p[0][0], p[1][0], 'go')
    i += 1

  plt.show()
  
  mask_in = triangulate(mask, 50)
  
  mask = np.zeros(im.shape)

  #plt.imshow(im, cmap='gray')
  

  i = 0
  for unit, length in units:
    vfl1 = points_out[i]
    vfl2 = points_in[i]

    p = orig + unit * (length - np.mean([vfl1, vfl2]) / 500.0 * length)
    #p = orig + unit * (length - (np.mean(vfl1)) / 500.0 * length)
    
    mask[p[1][0], p[0][0]] = 1
    i += 1
    #plt.plot(p[0][0], p[1][0], 'ro')

  #plt.show()
  
  mask_middle = triangulate(mask, 100)
  
  return [mask_out, mask_in, mask_middle]

def extract_perispinal_mask(aligned_dcm):
  
  bone_mask = extract_bone_mask(np.copy(aligned_dcm))
  
  center = measurements.center_of_mass(bone_mask)
  
  mask = np.copy(aligned_dcm)
  mask[bone_mask == 1] = 0
  
  mask = mask[center[0] - 75:center[0] + 125,
              center[1] - 125:center[1] + 125]
  
  mask[mask < -30] = 0
  mask[mask > 150] = 0
  
  mask = mask > 0

  paraspinal_muscle_mask = np.zeros((512, 512))
  paraspinal_muscle_mask[center[0] - 75:center[0] + 125,
              center[1] - 125:center[1] + 125] = mask
  
  paraspinal_mask = median_filter(paraspinal_muscle_mask, size=1)
  
  
  paraspinal_image = aligned_dcm[center[0] - 75:center[0] + 125,
              center[1] - 125:center[1] + 125]
  
  return [center, paraspinal_image, paraspinal_mask.astype(np.uint8)]

reference_body_mask = extract_body_mask(input_dcm_image) 
reference_body_image = extract_mask(input_dcm_image, reference_body_mask)

_, _, reference_paraspinal = extract_perispinal_mask(input_dcm_image)

reference_muscle2_image = read_image('./aligned/mean_vfat.png')
reference_muscle_image = read_image('./aligned/mean_sfat.png')

reference_paraspinal2 =  read_image('./aligned/mean_muscle.png')
  
mean_vfat222_image = np.zeros((512, 512))
mean_sfat222_image = np.zeros((512, 512))

for directory in os.listdir(root):
  
  i += 1
  
  #print(directory)
  
  if not os.path.isdir(root + directory): continue
  
  dcm_image_file = root + directory + '/image.dcm'
  if dcm_image_file == input_dcm_file: continue
  
  png_imagefile = root + directory + '/image.png'
  vfat_image_file = image_in_dir(root + directory + '/vfat_slices')
  sfat_image_file = image_in_dir(root + directory + '/sfat_slices')
  wallm_image_file = image_in_dir(root + directory + '/wallmuscle_slices')
  lpsoas_image_file = image_in_dir(root + directory + '/leftpsoas_slices')
  rpsoas_image_file = image_in_dir(root + directory + '/rightpsoas_slices')
  
  dcm_image, png_image = read_dcm(dcm_image_file)
  
  vfat_image_mask = read_image(vfat_image_file)
  sfat_image_mask = read_image(sfat_image_file)
  wallm_image_mask = read_image(wallm_image_file)
  lpsoas_image_mask = read_image(lpsoas_image_file)
  rpsoas_image_mask = read_image(rpsoas_image_file)
  muscle_mask = wallm_image_mask + lpsoas_image_mask + rpsoas_image_mask
  
  body_mask = extract_body_mask(dcm_image)
  body_image = extract_mask(dcm_image, body_mask)
  
  warp_matrix = align(reference_body_mask, body_mask)
  aligned_body_mask = warp_image(body_mask, body_mask, warp_matrix)
  aligned_dcm = warp_image(body_mask, dcm_image, warp_matrix)
  aligned_png = warp_image(body_mask, png_image, warp_matrix)
  
  bone_mask = extract_bone_mask(np.copy(aligned_dcm))
  bone_center = measurements.center_of_mass(bone_mask)
  
  aligned_vfat_mask_ref = warp_image(vfat_image_mask, vfat_image_mask, warp_matrix) > 0
  aligned_sfat_mask_ref = warp_image(sfat_image_mask, sfat_image_mask, warp_matrix) > 0
  psoas_mask = lpsoas_image_mask + rpsoas_image_mask
  aligned_psoas_mask_ref = warp_image(psoas_mask, psoas_mask, warp_matrix) > 0
  aligned_muscle_mask_ref = warp_image(muscle_mask, muscle_mask, warp_matrix) > 0
  
  #
  skin_mask = ndimage.binary_erosion(aligned_body_mask, iterations=10)
  #
  #stripped_skin = np.copy(aligned_dcm)
  #stripped_skin[skin_mask == 0] = np.min(stripped_skin)
  
  aligned_dcm_copy = np.copy(aligned_dcm)
  aligned_dcm_copy[skin_mask == 0] = np.min(aligned_dcm)
  
  ####
  
  bone_mask = extract_bone_mask(np.copy(aligned_dcm))
  bone_center = measurements.center_of_mass(bone_mask)
  
  #try:
  outer_wall_mask, inner_wall_mask, mask_middle = find_outer_wall(aligned_dcm_copy, aligned_body_mask)
  
  image_overlay(aligned_dcm, [inner_wall_mask], ['red'], show=True)
    
  fat_mask = (aligned_dcm <= -30) - (aligned_dcm < -150)
  sc_fat_mask = np.copy(fat_mask)
  vc_fat_mask = np.copy(fat_mask)
  
  sc_fat_mask[mask_middle == 1] = 0
  vc_fat_mask[mask_middle == 0] = 0
  
  muscle_mask = (aligned_dcm < 150) - (aligned_dcm < -29)
  muscle_mask[inner_wall_mask == 1] = 0
  muscle_mask[outer_wall_mask == 0] = 0
  muscle_mask[bone_mask == 1] = 0
  
  image_overlay(aligned_dcm, [vc_fat_mask, sc_fat_mask, muscle_mask], ['yellow', 'cyan', 'red', 'purple'], show=True)
  
  image_overlay(aligned_dcm, [np.logical_xor(vc_fat_mask, aligned_vfat_mask_ref)], ['red'], show=True)
  
  print(directory + ' --------------------------------')
  
  print('VC fat, Jaccard = %f, Error = %f' % ( jaccard_similarity_score(aligned_vfat_mask_ref, vc_fat_mask),
                                                  100 * ( np.absolute(np.sum(aligned_vfat_mask_ref) - 
                                                   np.sum(vc_fat_mask)) )  / np.sum(aligned_vfat_mask_ref) )  )
  
  
  print('SC fat, Jaccard = %f, Error = %f' % ( jaccard_similarity_score(aligned_sfat_mask_ref, sc_fat_mask),
                                                  100 * ( np.absolute(np.sum(aligned_sfat_mask_ref) - 
                                                   np.sum(sc_fat_mask)) )  / np.sum(aligned_sfat_mask_ref) ) )
  
  
  print('Muscle, Jaccard = %f, Error = %f' % ( jaccard_similarity_score(aligned_muscle_mask_ref, muscle_mask),
                                                  100 * ( np.absolute(np.sum(aligned_muscle_mask_ref) - 
                                                   np.sum(muscle_mask)) )  / np.sum(aligned_muscle_mask_ref) )  )
  
  
  #print(directory + ',' + str(jaccard_similarity_score(aligned_vfat_mask_ref, vc_fat_mask)) + ',' +
  #                        str(jaccard_similarity_score(aligned_sfat_mask_ref, sc_fat_mask)) + ',' +
  #                        str(jaccard_similarity_score(aligned_muscle_mask_ref, muscle_mask)) + ',' +
  #                        str(100.0 * ( np.absolute(np.sum(aligned_vfat_mask_ref) - 
  #                            np.sum(vc_fat_mask)) )  / np.sum(aligned_vfat_mask_ref) ) + ',' +
  #                        str(100.0 * ( np.absolute(np.sum(aligned_sfat_mask_ref) - 
  #                            np.sum(sc_fat_mask)) )  / np.sum(aligned_sfat_mask_ref) ) + ',' +
  #                        str(100.0 * ( np.absolute(np.sum(aligned_muscle_mask_ref) - 
  #                            np.sum(muscle_mask)) )  / np.sum(aligned_muscle_mask_ref)))
  
  
  #print(jaccard_similarity_score(aligned_vfat_mask_ref, vc_fat_mask))
  #print(jaccard_similarity_score(aligned_muscle_mask_ref, wall_muscle_mask))
  
  
  #print( 100 * ( np.absolute(np.sum(aligned_vfat_mask_ref) - np.sum(vc_fat_mask)) )  / np.sum(aligned_vfat_mask_ref))
  #print( 100 * ( np.absolute(np.sum(aligned_muscle_mask_ref) - np.sum(wall_muscle_mask)) )  / np.sum(aligned_muscle_mask_ref))
  
  #image_overlay(aligned_dcm, [vc_fat_mask, sc_fat_mask, muscle_mask], ['yellow', 'cyan', 'red', 'purple'], show=True)
  
  
  #if i > 3: exit()