'''
Description: This script was developed to test the repeatability and reproducibility of the manual thickness
measurements for the MULTIS project. This code grabs thickness xmls from two different directories and compares them
using Bland Altman plots. At this stage, both directories that are selected must have the same locations analyzed
because there is not a check to determine if the thicknesses being compared correlate to each other since they are
sorted by number. Figures showing the differences are included for thicknesses that are more than two standard
deviations away from the mean difference across all samples.
Original Author:
Erica Morrill
Department of Biomedical Engineering
Lerner Research Institute
Cleveland Clinic
Cleveland, OH
morrile2@ccf.org
'''
import matplotlib
# matplotlib.use("TkAgg")
import xml.etree.ElementTree as ET
from pylab import *
import numpy as np
import os
import tkMessageBox
import Tkinter as tk
import tkFileDialog
import matplotlib.cm as cm
def parseXMLinclusion(inclusion):
doc = ET.parse(inclusion)
root = doc.getroot()
thickness_xmls = []
for child in root:
for subChild in child:
XML = subChild.attrib["Anatomical"]
if XML != "None" and XML[-28:-27] == "A":
thickness_xmls.append(XML)
return thickness_xmls
def PlotSum(inclusion_xmls):
# Setup plot formatting for the Bland-Altman plot
fig = plt.figure(figsize=(15,11))
ax = fig.add_subplot(111)
ax.set_xticklabels([])
ax.set_yticklabels([])
ax1 = fig.add_subplot(411)
plt.title("Total")
ax2 = fig.add_subplot(412)
plt.title("Muscle")
ax3 = fig.add_subplot(413)
plt.title("Fat")
ax4 = fig.add_subplot(414)
plt.title("Skin")
# Turn off axis lines and ticks of the big subplot
ax.spines['top'].set_color('none')
ax.spines['bottom'].set_color('none')
ax.spines['left'].set_color('none')
ax.spines['right'].set_color('none')
ax.tick_params(labelcolor='none', top='off', bottom='off', left='off', right='off')
ax.set_axis_bgcolor('none')
# colors = cm.rainbow(np.linspace(0, 1, len(XMLS)))
DATA = []
locations = []
# Find all of the anatomical xml files in one directory at a time
for fname in inclusion_xmls:
dir = os.path.split(fname)[0]
XMLS = parseXMLinclusion(fname)
# Sort list by location number (first three numbers of the xml filename)
XMLS = sorted(XMLS, key=lambda x: x.split('/')[-1])
# Initialize thickness data lists
skin = []
fat = []
muscle = []
tag = []
# Read each xml file and append thickness data to lists
for i, xmlname in enumerate(XMLS):
doc = ET.parse(os.path.join(dir,xmlname))
root = doc.getroot()
subj = root.find('Subject')
src = subj.find('Source')
tag.append(xmlname)
loc = src.find("Frame")
thick = loc.find("Thickness")
skin.append(float(thick.find("Skin").text))
fat.append(float(thick.find("Fat").text))
muscle.append(float(thick.find("Muscle").text))
# Calculate the total thickness of tissue
total = [sum(x) for x in zip(skin, fat, muscle)]
# Store thicknesses and locations in list format
DATA.append([skin, fat, muscle, total])
locations.append(tag)
# Convert data to numpy array and calculate average and difference for each location
DATA = np.array(DATA)
DATA = np.concatenate((np.concatenate((DATA[0], DATA[2], DATA[4], DATA[6], DATA[8]), axis=1),
np.concatenate((DATA[1], DATA[3], DATA[5], DATA[7], DATA[9]), axis=1)))
DATA = DATA.reshape(2, 4, -1)
locations = np.array(locations)
locations = np.concatenate((np.concatenate((locations[0], locations[2], locations[4], locations[6], locations[8])),
np.concatenate((locations[1], locations[3], locations[5], locations[7],
locations[9]))))
locations = locations.reshape(2, -1)
avg = np.mean(np.abs(DATA), axis=0)
# DIFF = np.concatenate((DATA[0]-DATA[1], DATA[2]-DATA[3], DATA[4]-DATA[5], DATA[6]-DATA[7], DATA[8]-DATA[9]), axis=1)
DIFF = DATA[0]-DATA[1]
# Plot data in Bland-Altman form
ax1.scatter(avg[3], DIFF[3])
ax2.scatter(avg[2], DIFF[2])
ax3.scatter(avg[1], DIFF[1])
ax4.scatter(avg[0], DIFF[0])
# Determine average and standard deviation of the paired differences
error_avg = np.mean(DIFF, axis=1)
error_std = np.std(DIFF, axis=1)
error_max = np.max(np.abs(DIFF), axis=1)
error_min = np.min(np.abs(DIFF), axis=1)
print "Absolute Average", np.mean(np.abs(DIFF), axis=1)
print "Absolute Stdev", np.std(np.abs(DIFF), axis=1)
print "MIN", error_min
print "MAX", error_max
# Plot each of the locations that had a difference greater than two standard deviations from the mean difference
for i, err in enumerate(DIFF):
# Create labels for which thickness is responsible for the outlier
if i == 0:
desc = "Skin"
elif i == 1:
desc = "Fat"
elif i == 2:
desc = "Muscle"
elif i == 3:
desc = "Total"
#
# # Plot and save figure and print the locations for verification that they match
# for j, e in enumerate(err):
# if e < error_avg[i]-2*error_std[i] or e > error_avg[i]+2*error_std[i] and desc is not "Total":
# f = plt.figure(figsize=(20,10))
# f.suptitle(desc, fontsize=20)
# f.add_subplot(121)
# # plt.title('Rici')
# plt.title('Directory1')
# plt.imshow(imread(os.path.split(locations[0,j])[0]+'/ThicknessPNG/'+os.path.split(locations[0,j])[1][:-4]+'.png'))
# plt.axis("Off")
# f.add_subplot(122)
# # plt.title('Ahmet')
# plt.title('Directory2')
# plt.imshow(imread(os.path.split(locations[1, j])[0] + '/ThicknessPNG/' + os.path.split(locations[1, j])[1][:-4] + '.png'))
# plt.axis("Off")
# plt.tight_layout()
# # Option to save the comparison of images that have a difference that is greater than 2 standard deviations from the mean
# # f.savefig(os.path.split(os.path.split(locations[0,j])[0])[0]+'/R_A_Comparison/'+os.path.split(locations[0,j])[1][:-25])
# print(i, locations[0, j], locations[1,j])
# Print the average difference and standard deviation of the difference
print("Skin Average Diff = " + str(error_avg[0]) + " mm +/- " + str(error_std[0]) + " mm")
print("Fat Average Diff = " + str(error_avg[1]) + " mm +/- " + str(error_std[1]) + " mm")
print("Muscle Average Diff = " + str(error_avg[2]) + " mm +/- " + str(error_std[2]) + " mm")
print("Total Average Diff = " + str(error_avg[3]) + " mm +/- " + str(error_std[3]) + " mm")
print
# Graph labeling for the Bland-Altman plot, blue line added for mean difference and red dotted lines are added to show 2 standard deviations fromt the mean in either direction
ax.set_xlabel("Mean thickness")
ax.set_ylabel("Difference in thickness between datasets")
ax.yaxis.set_label_coords(-0.2, 0.5)
# start, end = ax1.get_ylim()
# ax1.yaxis.set_ticks(np.linspace(start, end, 3))
start, end = ax1.get_xlim()
start = 0
ax1.plot([start, end], [error_avg[3], error_avg[3]])
ax1.plot([start, end], [error_avg[3]+2*error_std[3], error_avg[3]+2*error_std[3]], 'r--')
ax1.plot([start, end], [error_avg[3] - 2 * error_std[3], error_avg[3] - 2 * error_std[3]], 'r--')
ax1.ticklabel_format(useOffset=False)
ax1.set_xlim(0, end)
print "Total"
for i, d in enumerate(DIFF[3]):
if d > error_avg[3]+2*error_std[3] or d < error_avg[3]-2*error_std[3]:
print locations[0][i]
print "Muscle"
for i, d in enumerate(DIFF[2]):
if d > error_avg[2] + 2 * error_std[2] or d < error_avg[2] - 2 * error_std[2]:
print locations[0][i]
print "Fat"
for i, d in enumerate(DIFF[1]):
if d > error_avg[1] + 2 * error_std[1] or d < error_avg[1] - 2 * error_std[1]:
print locations[0][i]
print "Skin"
for i, d in enumerate(DIFF[0]):
if d > error_avg[0] + 2 * error_std[0] or d < error_avg[0] - 2 * error_std[0]:
print locations[0][i]
# start, end = ax2.get_ylim()
# ax2.yaxis.set_ticks(np.linspace(start, end, 3))
start, end = ax2.get_xlim()
start = 0
ax2.plot([start, end], [error_avg[2], error_avg[2]])
ax2.plot([start, end], [error_avg[2] + 2 * error_std[2], error_avg[2] + 2 * error_std[2]], 'r--')
ax2.plot([start, end], [error_avg[2] - 2 * error_std[2], error_avg[2] - 2 * error_std[2]], 'r--')
ax2.ticklabel_format(useOffset=False)
ax2.set_xlim(0, end)
# start, end = ax3.get_ylim()
# ax3.yaxis.set_ticks(np.linspace(start, end, 3))
start, end = ax3.get_xlim()
start = 0
ax3.plot([start, end], [error_avg[1], error_avg[1]])
ax3.plot([start, end], [error_avg[1] + 2 * error_std[1], error_avg[1] + 2 * error_std[1]], 'r--')
ax3.plot([start, end], [error_avg[1] - 2 * error_std[1], error_avg[1] - 2 * error_std[1]], 'r--')
ax3.ticklabel_format(useOffset=False)
ax3.set_xlim(0, end)
# start, end = ax4.get_ylim()
# ax4.yaxis.set_ticks(np.linspace(start, end, 3))
start, end = ax4.get_xlim()
start = 0
ax4.plot([start, end], [error_avg[0], error_avg[0]])
ax4.plot([start, end], [error_avg[0] + 2 * error_std[0], error_avg[0] + 2 * error_std[0]], 'r--')
ax4.plot([start, end], [error_avg[0] - 2 * error_std[0], error_avg[0] - 2 * error_std[0]], 'r--')
ax4.ticklabel_format(useOffset=False)
ax4.set_xlim(0, end)
ax1.grid("on")
ax2.grid("on")
ax3.grid("on")
ax4.grid("on")
fig.tight_layout()
fig.subplots_adjust(right=0.8)
# plt.savefig('/home/morrile2/Documents/MULTIS Data/MULTIS_test/MULTIS_repeatability_intra02.svg')
plt.show()
return
class FileSelectionApp(tk.Tk):
def __init__(self):
tk.Tk.__init__(self)
# # Search for MULTIS_test directory from the home directory
# home = os.path.expanduser('~')
# for dirname, subdirList, fileList in os.walk(home):
# for dir in subdirList:
# if "MULTIS_test" in dir:
# print(dirname + dir)
# multis_dir = dirname + '/' + dir
# try:
# multis_dir
# except NameError:
# multis_dir = home
# Specify the directories to be compared
# dir1 = tkFileDialog.askdirectory(initialdir = multis_dir, title = "Select first directory")
# dir1 = '/home/morrile2/Documents/MULTIS_test/MULTIS001-1/Analysis_initial'
# dir2 = tkFileDialog.askdirectory(initialdir = multis_dir, title = "Select Second directory")
# dir2 = '/home/morrile2/Documents/MULTIS_test/MULTIS001-1/Analysis'
p1 = ''
p2 = '_2'
subj_ID = 'MULTIS008-1'
xml1 = '/home/morrile2/Documents/MULTIS Data/MULTIS_test/'+ subj_ID +'/TissueThickness/UltrasoundManual/'+ subj_ID +'_TA_inclusion' + p1 +'.xml'
xml2 = '/home/morrile2/Documents/MULTIS Data/MULTIS_test/'+ subj_ID +'/TissueThickness/UltrasoundManual/'+ subj_ID +'_TA_inclusion' + p2 + '.xml'
subj_ID = 'MULTIS023-1'
xml3 = '/home/morrile2/Documents/MULTIS Data/MULTIS_test/'+ subj_ID +'/TissueThickness/UltrasoundManual/'+ subj_ID +'_TA_inclusion' + p1 +'.xml'
xml4 = '/home/morrile2/Documents/MULTIS Data/MULTIS_test/'+ subj_ID +'/TissueThickness/UltrasoundManual/'+ subj_ID +'_TA_inclusion' + p2 + '.xml'
subj_ID = 'MULTIS042-1'
xml5 = '/home/morrile2/Documents/MULTIS Data/MULTIS_test/'+ subj_ID +'/TissueThickness/UltrasoundManual/'+ subj_ID +'_TA_inclusion' + p1 +'.xml'
xml6 = '/home/morrile2/Documents/MULTIS Data/MULTIS_test/'+ subj_ID +'/TissueThickness/UltrasoundManual/'+ subj_ID +'_TA_inclusion' + p2 + '.xml'
subj_ID = 'MULTIS061-1'
xml7 = '/home/morrile2/Documents/MULTIS Data/MULTIS_test/'+ subj_ID +'/TissueThickness/UltrasoundManual/'+ subj_ID +'_TA_inclusion' + p1 +'.xml'
xml8 = '/home/morrile2/Documents/MULTIS Data/MULTIS_test/'+ subj_ID +'/TissueThickness/UltrasoundManual/'+ subj_ID +'_TA_inclusion' + p2 + '.xml'
subj_ID = 'MULTIS087-1'
xml9 = '/home/morrile2/Documents/MULTIS Data/MULTIS_test/'+ subj_ID +'/TissueThickness/UltrasoundManual/'+ subj_ID +'_TA_inclusion' + p1 +'.xml'
xml10 = '/home/morrile2/Documents/MULTIS Data/MULTIS_test/'+ subj_ID +'/TissueThickness/UltrasoundManual/'+ subj_ID +'_TA_inclusion' + p2 + '.xml'
inclusion_xmls = [xml1, xml2, xml3, xml4, xml5, xml6, xml7, xml8, xml9, xml10]
PlotSum(inclusion_xmls)
# plt.show()
if __name__ == "__main__":
app = FileSelectionApp()
app.mainloop()