Sim annealing

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
COUNTER = 0

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 > 80
  bone_mask = binary_fill_holes(bone_mask)
  bone_mask = largest_object_mask(bone_mask)
  
  return bone_mask

from skimage.draw import polygon

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 filter_outliers(vfls):
#  
#  print(vfls3)
#
#  vfls3 = np.zeros(len(vfls))
#  
#  vfls1 = ndimage.filters.percentile_filter(vfls, 20, size=100)
#  vfls2 = ndimage.filters.percentile_filter(vfls, 80, size=100)
#  
#  for k in xrange(0, len(vfls)):
#
#    if vfls[k] < vfls1[k] or vfls[k] > vfls2[k]:
#      vfls3[k] = -1
#    else:
#      vfls3[k] = vfls[k]
#    
#    if k == 0 or k == len(vfls) - 1:
#      vfls3[k] = -1
#  
#  
#  #plt.plot(vfls3, 'ro')
#  #plt.plot(vfls, 'b+')
#  #plt.show()
#  
#  return vfls3
def fill_nan(A):
    '''
    interpolate to fill nan values
    '''
    inds = np.arange(A.shape[0])
    good = np.where(np.isfinite(A))
    f = interpolate.interp1d(inds[good], A[good],bounds_error=False)
    B = np.where(np.isfinite(A),A,f(inds))
    return B

from scipy import interpolate

def filter_outliers2(y, threshold):
  
  import statsmodels.api as sm

  x = np.arange(0, len(y)).astype(float)
  
  yc = np.zeros(len(y) + 200)
  yc[0:100] = y[len(y)-100:]
  yc[len(y)+100:] = y[0:100]
  yc[100:len(y)+100] = np.copy(y)
  
  xc = np.arange(0, yc.shape[0])
  
  lowess = sm.nonparametric.lowess(yc, xc, frac=0.2)

  smoothed = lowess[:, 1]
  
  smoothed = smoothed[100:len(y)+100]
  
  yp = np.zeros(len(y))
  
  for i in xrange(0, len(y)):
  
    if np.absolute(y[i] - smoothed[i]) > 25:
      yp[i] = float('NaN')
    else:
      yp[i] = y[i]
      
  yy = fill_nan(yp)
  
  
  #plt.plot(yy, 'r')
  #plt.plot(x, y, 'g+')
  #plt.show()
  
  return yy

def triangulate(points, mask):
  
  poly_mask = np.zeros(mask.shape)
  rr,cc = polygon(points[:,0],points[:,1],poly_mask.shape)
  
  poly_mask[rr,cc] = 1
  
  return poly_mask
  
def filter_outliers(y,perc,size):
  
  
  ys = np.zeros(len(y)+200)
  ys[0:100] = y[-100:]
  ys[-100:] = y[0:100]
  ys[100:-100] = y
  
  yp = np.copy(ys)
  
  #return median_filter(yp, size=100)[100:-100]
  
  return ndimage.filters.percentile_filter(yp,perc,size=size)[100:-100]
  
  print(y)
  
  #x = np.arange(0, len(ys)).astype(float)
  #lowess = sm.nonparametric.lowess(ys, x, frac=0.2)
  #smoothed = lowess[:, 1]
  #
  #smoothed2 = median_filter(ys, size=100)
  
  ys_smooth = median_filter(ys, size=100)
  
  dists = np.absolute(ys - ys_smooth)
  
  dists_smooth = median_filter(dists, size=100)
  yp[dists > dists_smooth] = float('NaN')

  yy = fill_nan(yp)
   
  return yy[100:-100]
  
  
  lowess = sm.nonparametric.lowess(ys, x, frac=0.2)
  smoothed = lowess[:, 1]# median_filter(ys, size=100)
  
  plt.plot(smoothed)
  plt.plot(ys)
  plt.show()
  
  distances = np.absolute(ys - smoothed)
  
  mean = np.mean(distances)
  std = np.std(distances)
  
  yy = np.zeros(len(y)+200)

  yy[0:100] = yy[-100:]
  yy[-100:] = yy[0:100]
  yy[100:-100] = np.copy(y)
  
  yy = yy.astype(float)
  
  yy[distances > 25] = float('NaN')
  
  yy = fill_nan(yy)
  
  plt.plot(y, 'r',  linewidth=1.0)
  plt.plot(yy[100:-100], 'b', linewidth=2.0)
  plt.show()
  
  
  return yy[100:-100]
   
from skimage.filter import canny
from skimage.measure import regionprops
def draw_contour_lines(directory, im, contour, center, params):

  global COUNTER
  # 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 = []
  
  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 % 1 == 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)
  #mat = np.gradient(mat)[1]
  
  
  #print(mat)
  #plt.imshow(mat)
  #plt.show()
  #
  x = []
  y = []
  
  last = 0
  
  for i in xrange(0, mat.shape[0]):
    found = False
    for j in xrange(0+int(params[0]), mat.shape[1]):
      if mat[i,j] > params[1]:
        y.append(j)
        found = True
        last = j
        break
    if not found:
      y.append(last)
    
    x.append(i)
  
  #lowess = sm.nonparametric.lowess(y, x, frac=0.01)
  #vfls = lowess[:, 1] 
  #plt.plot(y)
  #plt.show()
  
  vfls = filter_outliers(y, params[2], params[3])
  
  old_points_out = np.copy(y)
  #old_points_out[old_points_out == None] = vfls[old_points_out == None]
  points_out = np.copy(vfls)
  
  plt.imshow(mat.T, cmap='gray')
  
  for k in xrange(0, len(vfls)-1):
    plt.plot([k, k+1], [points_out[k],points_out[k+1]], 'y-', linewidth=1)
  

  plt.savefig('./cuts/' + directory + '_cut'+str(COUNTER)+'.png')
  plt.close()
  #plt.show()
  
  return [orig, units, points_out]
  
  #######
  
  x = []
  y = []
  last = None
  
  for i in xrange(0, mat.shape[0]):
    
    start = int(old_points_out[i])
    found = False
    for j in xrange(int(start+params[4]), mat.shape[1]):
      if mat[i,j] < params[5] or mat[i,j] > 120:
          y.append(j)
          last = j
          found = True
          break
    if not found:
      y.append(last)
    x.append(i)
  
  old_points_in = np.copy(y)
  vfls = filter_outliers(y, params[6], params[7])
  
  points_in = np.copy(vfls)
  
  #######
  
  x = []
  y = []
  last = 0.0 #####
  
  for i in xrange(0, mat.shape[0]):
    start = int(old_points_in[i])
    found = False
    for j in xrange(int(start+params[8]), mat.shape[1]):
      if mat[i,j] > params[9]:
          y.append(j)
          last = j
          found = True
          break
    if not found:
      y.append(last)
    x.append(i)
  
  vfls = y
  vfls = filter_outliers(y, params[10], params[11])
  points_in_in_in = np.copy(vfls)

  old_points_in_in_in = np.copy(y)
  #######
  
  x = []
  y = []
  last = 0.0 ###3
  
  for i in xrange(0, mat.shape[0]):
    start = int(old_points_in_in_in[i])
    found = False
    for j in xrange(int(start+params[12]), mat.shape[1]):
      if mat[i,j] < params[13] or mat[i,j] > 100:
          y.append(j)
          last = j
          found = True
          break
    if not found:
      y.append(last)
    x.append(i)
  
  vfls = y
  vfls = filter_outliers(y, params[14], params[15])
  points_in_in_in_in = np.copy(vfls)
  
  
  #######
  
  # [10, -30, 30, 100, 
  # 25, -30, 50, 100, 
  # 40, -30, 50, 100, 
  # 50, -50, 75, 100]
  
  # Plot outer points
  go_plot = True
  if go_plot:
    plt.imshow(mat.T, cmap='gray')
    
    for k in xrange(0, len(vfls)-1):
      plt.plot([k, k+1], [points_out[k],points_out[k+1]], 'y-', linewidth=1)
    
    for k in xrange(0, len(vfls)-1):
      plt.plot([k, k+1], [points_in[k],points_in[k+1]], 'g-', linewidth=1)
    
    for k in xrange(0, len(vfls)-1):
      plt.plot([k, k+1], [points_in_in_in[k],points_in_in_in[k+1]], 'c-', linewidth=1)
      
    for k in xrange(0, len(vfls)-1):
      plt.plot([k, k+1], [points_in_in_in_in[k],points_in_in_in_in[k+1]], 'm-', linewidth=1)
    

    plt.savefig('./cuts/' + directory + '_cut'+str(COUNTER)+'.png')
    plt.close()
    #plt.show()

  return [orig, units, points_out, points_in, points_in_in_in, points_in_in_in_in]

def find_fat_boundary(directory, im, body_mask, bone_mask, params):
  
  pp = np.copy(params)
  
  #params = [1,0,50,25,1,0,50,25,1,0,50,25,1,0,50,25]
  
  #params[0:4] = pp
  
  bone_center = measurements.center_of_mass(bone_mask)
  bone_center = [bone_center[1], regionprops(bone_mask)[0].bbox[0]]
  
  #print(bone_center)
  
  contour = measure.find_contours(body_mask, 0.8)[0]
  
  midpoint = int(contour.shape[0] / 2)
  
  contour2 = np.zeros(contour.shape)
  
  for k in xrange(0, contour.shape[0]):
    
    if k < (contour.shape[0]-midpoint):
      contour2[k,:] = contour[k+midpoint,:]
    else:
      contour2[k,:] = contour[(k-(contour.shape[0]-midpoint)),:]
    
  contour = contour2
  
  orig, units, points_out = draw_contour_lines(directory, im, contour, bone_center, params)

  return [orig, units, points_out]

import statsmodels.api as sm


def find_outer_wall(directory, im, body_mask, bone_mask, params):

  orig, units, points_out = find_fat_boundary(directory, im, body_mask, bone_mask, params)
  
  mask = np.zeros(im.shape)
  
  go_plot = False
  
  if go_plot: plt.imshow(im, cmap='gray')

  i = 0
  
  points = []
  #print(len(points_in))
  #print(len(units))
  for unit, length in units:
    
    if i >= len(units)-2: break
    if points_out[i] < 10:
      i += 1
      continue
    p = orig + unit * (length - points_out[i] / 500.0 * length)
    if math.isnan(p[0][0]) or math.isnan(p[1][0]):
      i += 1
      continue
    mask[p[1][0], p[0][0]] = 1
    points.append((p[1][0],p[0][0]))
    if go_plot: plt.plot(p[0][0], p[1][0], 'r+')
    i += 1
    
  #plt.show()
  points = np.asarray(points)
  
  mask_out = triangulate(points, mask) #, 100)
  
  return mask_out
  
  bone_center = measurements.center_of_mass(bone_mask)
  bone_top = regionprops(bone_mask)[0].bbox[0]
  
  mask = np.zeros(im.shape)

  i = 0
  
  points = []
  
  for unit, length in units:
    
    if i >= len(units)-2: break
    
    if points_in[i] < 0:
      i += 1
      continue
    
      
    diffx = (p[0][0] - bone_center[1])
    p = orig + unit * (length - points_in[i] / 500.0 * length)
    
    if np.absolute(diffx) < 100 and p[1][0] > bone_top:
      p = orig + unit * (length - points_in_in[i] / 500.0 * length)
  
      
    if math.isnan(p[0][0]) or math.isnan(p[1][0]):
      i += 1
      continue
      
    mask[p[1][0], p[0][0]] = 1
    points.append((p[1][0],p[0][0]))
    if go_plot: plt.plot(p[0][0], p[1][0], 'b+')
    i += 1
    
  #plt.show()
  
  points = np.asarray(points)
  mask_in = triangulate(points,mask) #ndimage.binary_erosion(triangulate(points,mask),iterations=4)
  
  # Show outer wall mask over image
  
  mask = np.zeros(im.shape)
  points = []
  i = 0
  
  for unit, length in units:
    
    if i >= len(units)-2: break
    vfl1 = points_out[i]
    vfl2 = points_in[i]
    
    p = orig + unit * np.mean([length - vfl1 / 500.0 * length, \
                               length - vfl2 / 500.0 * length])
    
    if math.isnan(p[0][0]) or math.isnan(p[1][0]):
      i += 1
      continue
    
    i += 1
    
    mask[p[1][0], p[0][0]] = 1
    points.append((p[1][0],p[0][0]))
    if go_plot: plt.plot(p[0][0], p[1][0], 'y+')

  if go_plot: plt.show()
  
  points = np.asarray(points)
  
  mask_middle = triangulate(points,mask)
  
  ############ PSOAS MUSCLE 
  
  i = 0
  
  mask1 = np.zeros(im.shape)
  mask2 = np.zeros(im.shape)
  points = []
  
  for unit, length in units:
    
    if i >= len(units)-2: break
    
    if points_in[i] < 0:
      i += 1
      continue
    
    p1 = orig + unit * (length - points_in_in_in_in[i] / 500.0 * length)
    p2 = orig + unit * (length - points_in_in[i] / 500.0 * length)
    
    diffx1 = (p1[0][0] - bone_center[1])
    diffx2 = (p2[0][0] - bone_center[1])
    
    if (np.absolute(diffx1) > 80) or \
       (np.absolute(diffx2) > 80 or np.absolute(diffx2) < 30 or p2[1][0] < bone_top):
      i += 1
      continue
    
    if math.isnan(p1[0][0]) or math.isnan(p1[1][0]) or \
       math.isnan(p2[0][0]) or math.isnan(p2[1][0]):
      i += 1
      continue
      
    if diffx2 < 0:
      mask1[p1[1][0], p1[0][0]] = 1.0
      mask1[p2[1][0], p2[0][0]] = 1.0
    else:
      mask2[p1[1][0], p1[0][0]] = 1.0
      mask2[p2[1][0], p2[0][0]] = 1.0
    
    if go_plot:
      plt.plot(p1[0][0], p1[1][0], 'b+')
    
    i += 1
    
  #plt.imshow(mask2)
  #plt.show()
  
  mask_lpsoas = convex_hull_image(mask1)
  mask_rpsoas = convex_hull_image(mask2)
  
  #image_overlay(im, [mask_lpsoas, mask_rpsoas], ['red', 'blue'], show=True)
  
  #image_overlay(im, [mask_out], ['red'], show=True)
  #image_overlay(im, [mask_in], ['red'], show=True)
  #image_overlay(im, [mask_middle], ['red'], show=True)
  
  mask_psoas = mask_lpsoas + mask_rpsoas
  
  return [mask_out, mask_in, mask_middle, mask_psoas.astype(bool)]

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

def do_image(params, directory):
  
  global root

  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
  psoas_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
  
  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)
  
  aligned_dcm_copy = np.copy(aligned_dcm)
  aligned_dcm_copy[skin_mask == 0] = np.min(aligned_dcm)
  
  fat_mask = (aligned_dcm_copy <= -30) - (aligned_dcm_copy < -190)
  muscle_mask = (aligned_dcm_copy <= 150) - (aligned_dcm_copy < -60)
  
  muscle_dcm = np.copy(aligned_dcm)
  muscle_dcm[muscle_mask == 0] = -100
  
  bone_mask = extract_bone_mask(np.copy(aligned_dcm))
  bone_center = measurements.center_of_mass(bone_mask)
  
  #try:
  outer_wall_mask = find_outer_wall(directory, aligned_dcm_copy, aligned_body_mask, bone_mask, params)
  
  #except:
  #  print('Error')
  #  continue
    
  #image_overlay(aligned_dcm, [inner_wall_mask], ['red'], show=True)
    
  fat_mask = (aligned_dcm <= -30) - (aligned_dcm < -190)
  sc_fat_mask = np.copy(fat_mask)
  #vc_fat_mask = np.copy(fat_mask)
  
  sc_fat_mask[outer_wall_mask == 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
  #
  #psoas_mask = (aligned_dcm < 150) - (aligned_dcm < -29)
  #
  #psoas_mask[mask_psoas == False] = 0
  #psoas_mask[muscle_mask == 1] = 0
  
  
  #psoas_mask2 = (aligned_dcm < 150) - (aligned_dcm < -29)
  #psoas_mask2[ndimage.binary_dilation(psoas_mask, iterations=4) == 0] = 0
  #
  #psoas_mask = psoas_mask2
  #
  #psoas_mask[bone_mask == 1] = 0
  
  ### Show results
  
  #image_overlay(aligned_dcm, [vc_fat_mask, sc_fat_mask, muscle_mask, psoas_mask], ['orange', 'yellow', 'cyan', 'red'], show=False, file='./output/' + directory + '.png')
  
  e1 = 1-jaccard_similarity_score(aligned_sfat_mask_ref, sc_fat_mask)
  
  #e1 = jaccard_similarity_score(aligned_vfat_mask_ref, vc_fat_mask)
  
  #e2 = jaccard_similarity_score(aligned_sfat_mask_ref, sc_fat_mask)
  #e3 = jaccard_similarity_score(aligned_muscle_mask_ref, muscle_mask)
  
  #err = (1 - e1) + (1 - e2) + (1 - e3)
  
  print('Step - err ' + str(e1) + ', params: ' + str(params) )
  
  return e1
  
   
  #image_overlay(aligned_dcm, [np.logical_xor(sc_fat_mask, aligned_sfat_mask_ref)], ['red'], show=True)
  
  #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))  )
  
  #show(aligned_psoas_mask_ref)
  #show(psoas_mask)
  
  print('Psoas, Jaccard = %f, Error = %f' % ( jaccard_similarity_score(aligned_psoas_mask_ref, psoas_mask),
                                                  100 * ( np.absolute(np.sum(aligned_psoas_mask_ref) - 
                                                   np.sum(psoas_mask)) )  / np.sum(aligned_psoas_mask_ref) ) )
                                                   
  return err                                               

def do_images(x):

  global COUNTER
  global root
  
  COUNTER += 1
  
  errors = []
  i = 0
  
  for directory in os.listdir(root):
  
    i += 1
  
    if i > 8: continue
  
    if not os.path.isdir(root + directory): continue
  
    print(directory + ' --------------------------------')
    errors.append(do_image(x,directory))
  
  print('===================== counter ' + str(COUNTER) + ', mean ' + str(np.mean(errors)))
  print("\n\n")
  return np.mean(errors)
  
def printf(text):
   print(text)

  
f = lambda x: do_images(x)
c = lambda x: printf("    " + str(x) + "\n")

x0 = [10, -30, 30, 100, 25, -30, 50, 100, 40, -30, 50, 100, 50, -50, 75, 100]

#!
x0 = [10, -30, 30, 100] #,1,0,50,25,1,0,50,25,1,0,50,25]
# [1,-5,10,25,5,-5,10,25,5,-5,10,25,5,-5,10,25]

opt_info = scipy.optimize.anneal(f, x0, lower=[1,-100,5,50], upper=[15,0,95,200]) # , callback=c

xopt = opt_info[0]