https://github.com/jennyhelyanwe/post_MI_postprocessing
Tip revision: 48c36cc0c5a8e853a4d97b18309493957e608091 authored by Jenny on 03 December 2024, 11:43:08 UTC
Cell model copyright and clean version
Cell model copyright and clean version
Tip revision: 48c36cc
ECGPV_visualisation.py
import os, sys
import numpy
import matplotlib
#matplotlib.use('tkagg') # For use on CSCS daint
from matplotlib import pyplot as plt
from matplotlib.gridspec import GridSpec
import matplotlib.animation as animation
import time
import argparse
import pickle
from scipy.signal import lfilter, resample
from scipy import io
class ECGPV_visualisation:
def __init__(self, CL):
self.CL = CL
self.beat_fig_size = [5, 5]
self.pv_fig_size = [6, 7]
def read_ecg_pv(self, name, dir):
meshname = dir+name
print(meshname)
if os.path.exists(dir+'ecgs.pl'):
ecgs = pickle.load(open(dir+'ecgs.pl', 'rb'))
else:
ecgs = self._read_ECG(meshname)
pickle.dump(ecgs, open(dir+'ecgs.pl', 'wb'))
if os.path.exists(dir+'pvs.pl'):
pvs = pickle.load(open(dir+'pvs.pl', 'rb'))
else:
pvs = self._read_PV(meshname)
pickle.dump(pvs, open(dir+'pvs.pl', 'wb'))
# Save as .mat file for delineation using Karlsruhe code:
raw_leads = numpy.vstack([ecgs['I']/ecgs['max_all_leads'],ecgs['II']/ecgs['max_all_leads'],ecgs['III']/ecgs['max_all_leads'],ecgs['aVL']/ecgs['max_all_leads'],ecgs['aVR']/ecgs['max_all_leads'],ecgs['aVF']/ecgs['max_all_leads'],ecgs['V1']/ecgs['max_all_leads'],ecgs['V2']/ecgs['max_all_leads'],ecgs['V3']/ecgs['max_all_leads'],ecgs['V4']/ecgs['max_all_leads'],ecgs['V5']/ecgs['max_all_leads'],ecgs['V6']/ecgs['max_all_leads']]).T
raw_leads_resampled = numpy.zeros((int(ecgs['t'][-1]/0.002), 12))
# Down sample to 500 Hz
for i in range(0, 12):
nsample = int(ecgs['t'][-1]/0.002)
raw_leads_resampled[:,i] = resample(raw_leads[:,i], nsample)
io.savemat('raw_ecg_leads.mat', {'signal':raw_leads_resampled})
io.savemat('ecgs.mat', {'ecgs': ecgs})
io.savemat('pvs.mat', {'pvs': pvs})
return ecgs, pvs
def _read_ECG(self, meshname):
filename = meshname + '.exm.vin'
with open(filename,'r') as f:
data = f.readlines()
LA = RA= LL = RL = V1 = V2 = V3 = V4 = V5 = V6 = t= []
# First 7 lines are header lines.
if len(data) > 7:
LA = numpy.zeros(len(data)-7)
RA = numpy.zeros(len(data)-7)
LL = numpy.zeros(len(data)-7)
RL = numpy.zeros(len(data)-7)
V1 = numpy.zeros(len(data)-7)
V2 = numpy.zeros(len(data)-7)
V3 = numpy.zeros(len(data)-7)
V4 = numpy.zeros(len(data)-7)
V5 = numpy.zeros(len(data)-7)
V6 = numpy.zeros(len(data)-7)
t = numpy.zeros(len(data)-7)
for i in range(7, len(data)):
t[i-7] = float(data[i].split()[-12])
LA[i-7] = float(data[i].split()[-10])
RA[i-7] = float(data[i].split()[-9])
LL[i-7] = float(data[i].split()[-8])
RL[i-7] = float(data[i].split()[-7])
V1[i-7] = float(data[i].split()[-6])
V2[i-7] = float(data[i].split()[-5])
V3[i-7] = float(data[i].split()[-4])
V4[i-7] = float(data[i].split()[-3])
V5[i-7] = float(data[i].split()[-2])
V6[i-7] = float(data[i].split()[-1])
sort_i = numpy.argsort(t)
t = t[sort_i]
LA = LA[sort_i]
RA = RA[sort_i]
LL = LL[sort_i]
RL = RL[sort_i]
V1 = V1[sort_i]
V2 = V2[sort_i]
V3 = V3[sort_i]
V4 = V4[sort_i]
V5 = V5[sort_i]
V6 = V6[sort_i]
# Ealuate Wilson's central terminal
VW = 1.0/3.0*(RA + LA + LL)
# Evaluate simulated ECG lead traces
V1 = V1 - VW
V2 = V2 - VW
V3 = V3 - VW
V4 = V4 - VW
V5 = V5 - VW
V6 = V6 - VW
I = LA - RA
II = LL - RA
III = LL - LA
aVL = LA - (RA + LL)/2.0
aVF = LL - (LA + RA)/2.0
aVR = RA - (LA + LL)/2.0
all_leads = numpy.concatenate((V1,V2,V3,V4,V5,V6,I,II,III,aVR,aVL,aVF))
precord_leads = numpy.concatenate((V1,V2,V3,V4,V5,V6))
limb_leads = numpy.concatenate((I,II,III,aVR,aVL,aVF))
max_all_leads = max(abs(all_leads))
max_precord_leads = max(abs(precord_leads))
max_limb_leads = max(abs(limb_leads))
# Divide into beats
ts, V1s = self._divide_signal(t, V1, self.CL)
ts, V2s = self._divide_signal(t, V2, self.CL)
ts, V3s = self._divide_signal(t, V3, self.CL)
ts, V4s = self._divide_signal(t, V4, self.CL)
ts, V5s = self._divide_signal(t, V5, self.CL)
ts, V6s = self._divide_signal(t, V6, self.CL)
ts, aVRs = self._divide_signal(t, aVR, self.CL)
ts, aVLs = self._divide_signal(t, aVL, self.CL)
ts, aVFs = self._divide_signal(t, aVF, self.CL)
ts, Is = self._divide_signal(t, I, self.CL)
ts, IIs = self._divide_signal(t, II, self.CL)
ts, IIIs = self._divide_signal(t, III, self.CL)
output_dict = {'t':t,'ts':ts,'V1':V1,'V2':V2,'V3':V3,'V4':V4,'V5':V5,'V6':V6,
'aVR':aVR,'aVL':aVL,'aVF':aVF,'I':I,'II':II,'III':III,
'V1s':V1s,'V2s':V2s,'V3s':V3s,'V4s':V4s,'V5s':V5s,'V6s':V6s,
'aVRs':aVRs,'aVLs':aVLs,'aVFs':aVFs,'Is':Is,'IIs':IIs,'IIIs':IIIs,
'max_all_leads':max_all_leads,'max_precord_leads':max_precord_leads,
'max_limb_leads':max_limb_leads}
return output_dict
def _read_PV(self,meshname):
filename = meshname + '-cardiac-cycle.sld.res'
with open(filename, 'r') as f:
data = f.readlines()
if (len(data)>18):
pl = numpy.zeros(len(data)-17)
pr = numpy.zeros(len(data)-17)
vl= numpy.zeros(len(data)-17)
vr = numpy.zeros(len(data)-17)
curtime = numpy.zeros(len(data)-17)
phasel = numpy.zeros(len(data)-17)
phaser = numpy.zeros(len(data)-17)
vl[0] = float(data[18].split()[5])
vr[0] = float(data[18].split()[-3])
for i in range(18, len(data)):
pl[i-17] = float(data[i].split()[6])/10000
pr[i-17] = float(data[i].split()[-2])/10000
vl[i-17] = float(data[i].split()[5])
vr[i-17] = float(data[i].split()[-3])
curtime[i-17] = float(data[i].split()[1])
phasel[i-17] = float(data[i].split()[4])
phaser[i-17] = float(data[i].split()[10])
sort_i = numpy.argsort(curtime)
curtime = curtime[sort_i]
pl = pl[sort_i]
pr = pr[sort_i]
vl = vl[sort_i]
vr = vr[sort_i]
phasel = phasel[sort_i]
phaser = phaser[sort_i]
# Divide into beats
ts, pls = self._divide_signal(curtime, pl, self.CL)
ts, vls = self._divide_signal(curtime, vl, self.CL)
ts, prs = self._divide_signal(curtime, pr, self.CL)
ts, vrs = self._divide_signal(curtime, vr, self.CL)
ts, phasels = self._divide_signal(curtime, phasel, self.CL)
ts, phasers = self._divide_signal(curtime, phaser, self.CL)
output_dict = {'t':curtime, 'ts':ts, 'pl':pl, 'vl':vl, 'pr':pr, 'vr':vr, 'phasel':phasel, 'phaser':phaser,
'pls':pls, 'vls':vls, 'prs':prs, 'vrs':vrs, 'phasels':phasels, 'phasers':phasers,
'plabel':'Pressure (kPa)', 'vlabel':'Volume (mL)'}
return output_dict
def _divide_signal(self, curtime, signal, CL):
if (curtime.max() > CL):
t_offsets = []
signals = []
idx_start = []
idx_end = []
for i in range(0, int(curtime.max()/CL)+1):
idx_start.append(numpy.where(curtime>=i*CL)[0][0])
if curtime.max() > (i+1)*CL:
idx_end.append(numpy.where(curtime>=(i+1)*CL)[0][0])
else:
idx_end.append(len(curtime)-1)
t_offsets.append(curtime[idx_start[i]:idx_end[i]]-curtime[idx_start[i]])
signals.append(signal[idx_start[i]:idx_end[i]])
else:
t_offsets = [curtime]
signals = [signal]
return t_offsets, signals
def plot_ecgpv_live(self, ecgs, pvs, title, show, ecgs2=[], pvs2=[]):
print('Plotting ECG PV live')
matplotlib.rcParams.update({'font.size':'11'})
matplotlib.rcParams.update({'text.color':'black'})
matplotlib.rcParams.update({'lines.linewidth':'1'})
fig = plt.figure(tight_layout=True, figsize=[15,7])
fig.suptitle(title)
gs = GridSpec(3,6)
axs = []
axs.append(fig.add_subplot(gs[:,0]))
axs.append(fig.add_subplot(gs[0,1]))
axs.append(fig.add_subplot(gs[1,1]))
axs.append(fig.add_subplot(gs[2,1]))
axs.append(fig.add_subplot(gs[0,2]))
axs.append(fig.add_subplot(gs[1,2]))
axs.append(fig.add_subplot(gs[2,2]))
axs.append(fig.add_subplot(gs[0,3]))
axs.append(fig.add_subplot(gs[1,3]))
axs.append(fig.add_subplot(gs[2,3]))
axs.append(fig.add_subplot(gs[0,4]))
axs.append(fig.add_subplot(gs[1,4]))
axs.append(fig.add_subplot(gs[2,4]))
axs.append(fig.add_subplot(gs[0,5]))
axs.append(fig.add_subplot(gs[1,5]))
axs.append(fig.add_subplot(gs[2,5]))
def animate_single(i):
# Plot PV
axs[0].clear()
axs[0].plot(pvs['vl'], pvs['pl'], pvs['vr'], pvs['pr'])
axs[1].clear()
axs[1].plot(pvs['t'], pvs['pl'], pvs['t'], pvs['pr'])
axs[2].clear()
axs[2].plot(pvs['t'], pvs['vl'], pvs['t'], pvs['vr'])
axs[4].clear()
axs[3].plot(pvs['t'], pvs['phasel'], pvs['t'], pvs['phaser'])
# Plot ECGs:
axs[4].clear()
axs[4].plot(ecgs['t'], ecgs['I']/ecgs['max_limb_leads'])
axs[5].clear()
axs[5].plot(ecgs['t'], ecgs['II']/ecgs['max_limb_leads'])
axs[6].clear()
axs[6].plot(ecgs['t'], ecgs['III']/ecgs['max_limb_leads'])
axs[7].clear()
axs[7].plot(ecgs['t'], ecgs['aVR']/ecgs['max_limb_leads'])
axs[8].clear()
axs[8].plot(ecgs['t'], ecgs['aVL']/ecgs['max_limb_leads'])
axs[9].clear()
axs[9].plot(ecgs['t'], ecgs['aVF']/ecgs['max_limb_leads'])
axs[10].clear()
axs[10].plot(ecgs['t'], ecgs['V1']/ecgs['max_precord_leads'])
axs[11].clear()
axs[11].plot(ecgs['t'], ecgs['V2']/ecgs['max_precord_leads'])
axs[12].clear()
axs[12].plot(ecgs['t'], ecgs['V3']/ecgs['max_precord_leads'])
axs[13].clear()
axs[13].plot(ecgs['t'], ecgs['V4']/ecgs['max_precord_leads'])
axs[14].clear()
axs[14].plot(ecgs['t'], ecgs['V5']/ecgs['max_precord_leads'])
axs[15].clear()
axs[15].plot(ecgs['t'], ecgs['V6']/ecgs['max_precord_leads'])
axs[0].set_xlabel('Volume (mL)')
axs[1].set_xlabel('Time (ms)')
axs[2].set_xlabel('Time (ms)')
axs[3].set_xlabel('Time (ms)')
axs[0].set_ylabel('Pressure (kPa)')
axs[1].set_ylabel('Pressure (kPa)')
axs[2].set_ylabel('Volume (mL)')
axs[3].set_ylabel('Phase')
axs[3].set_ylim([0,5])
axs[1].set_title('Pressure transient')
axs[2].set_title('Volume transient')
axs[3].set_title('Cardiac cycle')
axs[4].set_xlabel('Time (s)')
axs[5].set_xlabel('Time (s)')
axs[6].set_xlabel('Time (s)')
axs[7].set_xlabel('Time (s)')
axs[8].set_xlabel('Time (s)')
axs[9].set_xlabel('Time (s)')
axs[10].set_xlabel('Time (s)')
axs[11].set_xlabel('Time (s)')
axs[12].set_xlabel('Time (s)')
axs[13].set_xlabel('Time (s)')
axs[14].set_xlabel('Time (s)')
axs[15].set_xlabel('Time (s)')
axs[4].set_ylabel('Normalised ECG')
axs[5].set_ylabel('Normalised ECG')
axs[6].set_ylabel('Normalised ECG')
axs[7].set_ylabel('Normalised ECG')
axs[8].set_ylabel('Normalised ECG')
axs[9].set_ylabel('Normalised ECG')
axs[10].set_ylabel('Normalised ECG')
axs[11].set_ylabel('Normalised ECG')
axs[12].set_ylabel('Normalised ECG')
axs[13].set_ylabel('Normalised ECG')
axs[14].set_ylabel('Normalised ECG')
axs[15].set_ylabel('Normalised ECG')
axs[4].set_ylim(-1,1)
axs[5].set_ylim(-1,1)
axs[6].set_ylim(-1,1)
axs[7].set_ylim(-1,1)
axs[8].set_ylim(-1,1)
axs[9].set_ylim(-1,1)
axs[10].set_ylim(-1,1)
axs[11].set_ylim(-1,1)
axs[12].set_ylim(-1,1)
axs[13].set_ylim(-1,1)
axs[14].set_ylim(-1,1)
axs[15].set_ylim(-1,1)
axs[4].set_title('I')
axs[5].set_title('II')
axs[6].set_title('III')
axs[7].set_title('aVR')
axs[8].set_title('aVL')
axs[9].set_title('aVF')
axs[10].set_title('V1')
axs[11].set_title('V2')
axs[12].set_title('V3')
axs[13].set_title('V4')
axs[14].set_title('V5')
axs[15].set_title('V6')
def animate_double(i):
# Plot PV
axs[0].clear()
axs[0].plot(pvs['vl'], pvs['pl'], pvs['vr'], pvs['pr'], pvs2['vl'],pvs2['pl'],'--', pvs2['vr'], pvs2['pr'], '--')
axs[1].clear()
axs[1].plot(pvs['t'], pvs['pl'], pvs['t'], pvs['pr'],pvs2['t'], pvs2['pl'],'--', pvs2['t'], pvs2['pr'], '--')
axs[2].clear()
axs[2].plot(pvs['t'], pvs['vl'], pvs['t'], pvs['vr'],pvs2['t'], pvs2['vl'], '--', pvs2['t'], pvs2['vr'], '--')
axs[4].clear()
axs[3].plot(pvs['t'], pvs['phasel'], pvs['t'], pvs['phaser'],pvs2['t'], pvs2['phasel'], '--', pvs2['t'], pvs2['phaser'], '--')
axs[0].set_xlabel('Volume (mL)')
axs[1].set_xlabel('Time (ms)')
axs[2].set_xlabel('Time (ms)')
axs[3].set_xlabel('Time (ms)')
axs[0].set_ylabel('Pressure (kPa)')
axs[1].set_ylabel('Pressure (kPa)')
axs[2].set_ylabel('Volume (mL)')
axs[3].set_ylabel('Phase')
axs[3].set_ylim([0,5])
axs[1].set_title('Pressure transient')
axs[2].set_title('Volume transient')
axs[3].set_title('Cardiac cycle')
# Plot ECGs:
max_all_leads = max([ecgs['max_all_leads'], ecgs2['max_all_leads']])
max_precord_leads = max([ecgs['max_precord_leads'], ecgs2['max_precord_leads']])
max_limb_leads = max([ecgs['max_limb_leads'], ecgs2['max_limb_leads']])
axs[4].clear()
axs[4].plot(ecgs['t'], ecgs['I']/max_limb_leads, ecgs2['t'], ecgs2['I']/max_limb_leads, '--')
axs[5].clear()
axs[5].plot(ecgs['t'], ecgs['II']/max_limb_leads, ecgs2['t'],ecgs2['II']/max_limb_leads, '--')
axs[6].clear()
axs[6].plot(ecgs['t'], ecgs['III']/max_limb_leads, ecgs2['t'], ecgs2['III']/max_limb_leads, '--')
axs[7].clear()
axs[7].plot(ecgs['t'], ecgs['aVR']/max_limb_leads, ecgs2['t'], ecgs2['aVR']/max_limb_leads, '--')
axs[8].clear()
axs[8].plot(ecgs['t'], ecgs['aVL']/max_limb_leads, ecgs2['t'], ecgs2['aVL']/max_limb_leads, '--')
axs[9].clear()
axs[9].plot(ecgs['t'], ecgs['aVF']/max_limb_leads, ecgs2['t'], ecgs2['aVF']/max_limb_leads, '--')
axs[10].clear()
axs[10].plot(ecgs['t'], ecgs['V1']/max_precord_leads, ecgs2['t'], ecgs2['V1']/max_precord_leads, '--')
axs[11].clear()
axs[11].plot(ecgs['t'], ecgs['V2']/max_precord_leads, ecgs2['t'], ecgs2['V2']/max_precord_leads, '--')
axs[12].clear()
axs[12].plot(ecgs['t'], ecgs['V3']/max_precord_leads, ecgs2['t'], ecgs2['V3']/max_precord_leads, '--')
axs[13].clear()
axs[13].plot(ecgs['t'], ecgs['V4']/max_precord_leads, ecgs2['t'], ecgs2['V4']/max_precord_leads, '--')
axs[14].clear()
axs[14].plot(ecgs['t'], ecgs['V5']/max_precord_leads, ecgs2['t'], ecgs2['V5']/max_precord_leads, '--')
axs[15].clear()
axs[15].plot(ecgs['t'], ecgs['V6']/max_precord_leads, ecgs2['t'], ecgs2['V6']/max_precord_leads, '--')
axs[4].set_xlabel('Time (s)')
axs[5].set_xlabel('Time (s)')
axs[6].set_xlabel('Time (s)')
axs[7].set_xlabel('Time (s)')
axs[8].set_xlabel('Time (s)')
axs[9].set_xlabel('Time (s)')
axs[10].set_xlabel('Time (s)')
axs[11].set_xlabel('Time (s)')
axs[12].set_xlabel('Time (s)')
axs[13].set_xlabel('Time (s)')
axs[14].set_xlabel('Time (s)')
axs[15].set_xlabel('Time (s)')
axs[4].set_ylabel('Normalised ECG')
axs[5].set_ylabel('Normalised ECG')
axs[6].set_ylabel('Normalised ECG')
axs[7].set_ylabel('Normalised ECG')
axs[8].set_ylabel('Normalised ECG')
axs[9].set_ylabel('Normalised ECG')
axs[10].set_ylabel('Normalised ECG')
axs[11].set_ylabel('Normalised ECG')
axs[12].set_ylabel('Normalised ECG')
axs[13].set_ylabel('Normalised ECG')
axs[14].set_ylabel('Normalised ECG')
axs[15].set_ylabel('Normalised ECG')
axs[4].set_ylim(-1,1)
axs[5].set_ylim(-1,1)
axs[6].set_ylim(-1,1)
axs[7].set_ylim(-1,1)
axs[8].set_ylim(-1,1)
axs[9].set_ylim(-1,1)
axs[10].set_ylim(-1,1)
axs[11].set_ylim(-1,1)
axs[12].set_ylim(-1,1)
axs[13].set_ylim(-1,1)
axs[14].set_ylim(-1,1)
axs[15].set_ylim(-1,1)
axs[4].set_title('I')
axs[5].set_title('II')
axs[6].set_title('III')
axs[7].set_title('aVR')
axs[8].set_title('aVL')
axs[9].set_title('aVF')
axs[10].set_title('V1')
axs[11].set_title('V2')
axs[12].set_title('V3')
axs[13].set_title('V4')
axs[14].set_title('V5')
axs[15].set_title('V6')
# if ecgs2:
# ani = animation.FuncAnimation(fig, animate_double, fargs=( axs, ecgs, pvs, ecgs2, pvs2), interval=1000)
# else:
ani = animation.FuncAnimation(fig, animate_single, interval=1000)
if show:
plt.show(block=True)
def _set_ecg_ticks(self,ax, t_end, CL):
t_start = 0
t_end = numpy.ceil(t_end/CL) * CL
t_end = int(numpy.ceil(t_end*10.0/2.0)*2)/10.0
#t_end = numpy.round(t_end*10.0)/10.0
#t_end = numpy.round(t_end)
#t_end = self.CL
minor_ticks = numpy.arange(0, t_end+0.04, 0.04)
major_ticks = numpy.arange(0, t_end+0.2, 0.2)
#minor_ticks = numpy.linspace(0, t_end, int((t_end-t_start)/0.04) + 1)
#major_ticks = numpy.linspace(0, t_end, int((t_end-t_start)/0.2) + 1)
ax.set_xticks(minor_ticks, minor=True)
ax.set_xticks(major_ticks)
minor_ticks = numpy.linspace(-1, 1, 21)
major_ticks = numpy.linspace(-1, 1, 5)
ax.set_yticks(minor_ticks, minor=True)
ax.set_yticks(major_ticks)
ax.grid(which="minor", color='r', linestyle='-', linewidth=1, alpha=0.5)
ax.grid(which="major", color='r', linestyle='-', linewidth=2, alpha=0.5)
def _set_full_ecg_ticks(self,ax, t_end, v_max, v_min, CL):
t_start = 0
t_end = numpy.ceil(t_end/CL) * CL
t_end = int(numpy.ceil(t_end*10.0/2.0)*2)/10.0
minor_ticks = numpy.arange(0, t_end+0.04, 0.04)
major_ticks = numpy.arange(0, t_end+0.2, 0.2)
ax.set_xticks(minor_ticks, minor=True)
ax.set_xticks(major_ticks)
v_max = numpy.ceil(v_max)
minor_ticks = numpy.arange(v_min, v_max+0.1, 0.1)
major_ticks = numpy.arange(v_min, v_max+0.5, 0.5)
ax.set_yticks(minor_ticks, minor=True)
ax.set_yticks(major_ticks)
ax.grid(which="minor", color='r', linestyle='-', linewidth=1, alpha=0.3)
ax.grid(which="major", color='r', linestyle='-', linewidth=2, alpha=0.3)
ax.spines['right'].set_visible(False)
ax.spines['left'].set_visible(False)
ax.spines['top'].set_visible(False)
ax.spines['bottom'].set_visible(False)
for tic in ax.xaxis.get_major_ticks():
tic.tick1line.set_visible(False)
tic.tick2line.set_visible(False)
tic.label1.set_visible(False)
tic.label2.set_visible(False)
for tic in ax.xaxis.get_minor_ticks():
tic.tick1line.set_visible(False)
tic.tick2line.set_visible(False)
tic.label1.set_visible(False)
tic.label2.set_visible(False)
for tic in ax.yaxis.get_major_ticks():
tic.tick1line.set_visible(False)
tic.tick2line.set_visible(False)
tic.label1.set_visible(False)
tic.label2.set_visible(False)
for tic in ax.yaxis.get_minor_ticks():
tic.tick1line.set_visible(False)
tic.tick2line.set_visible(False)
tic.label1.set_visible(False)
tic.label2.set_visible(False)
def _set_pv_ticks(self,ax):
x_ticks = numpy.arange(50, 250, 50)
y_ticks = numpy.arange(0, 140, 20)
ax.set_xticks(x_ticks)
ax.set_yticks(y_ticks)
#ax.grid(which="major", color='k', linestyle='--', linewidth=1, alpha=0.5)
def plot_ecg_lead(self, ecgs, lead_name, filename,show, beat=0, ecg2=[]):
print('Plotting ECG lead '+str(lead_name))
matplotlib.rcParams.update({'font.size':'24'})
matplotlib.rcParams.update({'text.color':'black'})
matplotlib.rcParams.update({'lines.linewidth':'3'})
if ecg2:
if beat > 0:
max_all_leads = max([ecgs['max_all_leads'], ecgs2['max_all_leads']])
self._plot_double(self.beat_fig_size,
ecgs['ts'][beat-1],ecgs[lead_name+'s'][beat-1]/max_all_leads,'-',
ecgs2['ts'][beat-1], ecgs[lead_name+'s'][beat-1]/max_all_leads,'--',
'Time (s)', lead_name, filename+'_beat'+str(beat)+'_comparison.png',
show,ecg_grid=True)
else:
figsize=[self.beat_fig_size[0]*(int(ecgs['t'][-1]/1.0)+1), self.beat_fig_size[1]]
self._plot_double(figsize, ecgs['t'], ecgs[lead_name]/max_all_leads,'-',
ecgs2['t'], ecgs2[lead_name]/max_all_leads,'--'
'Time (s)', lead_name, filename+'_full_comparison.png',show, ecg_grid=True)
else:
if beat > 0:
self._plot_single(self.beat_fig_size, ecgs['ts'][beat-1],
ecgs[lead_name+'s'][beat-1]/ecgs['max_all_leads'],
'Time (s)', lead_name, filename+'_beat'+str(beat)+'.png',
show, ecg_grid=True)
else:
figsize=[self.beat_fig_size[0]*(int(ecgs['t'][-1]/1.0)+1), self.beat_fig_size[1]]
self._plot_single(figsize, ecgs['t'], ecgs[lead_name]/ecgs['max_all_leads'],
'Time (s)', lead_name, filename+'_full.png', show,ecg_grid=True)
def plot_ecg_all_leads(self, ecgs, show, beat=0, ecgs2=[]):
print('Plotting all ECG leads')
matplotlib.rcParams.update({'font.size':'24'})
matplotlib.rcParams.update({'text.color':'black'})
matplotlib.rcParams.update({'lines.linewidth':'1'})
if ecgs2:
print('Full ECG comparison not yet implemented...')
else:
if beat > 0:
# t_padding
t_padding = 0.2
v_padding = 0.5
# Generate calibrating step function
t_res = ecgs['t'][1] - ecgs['t'][0]
t_calib = numpy.arange(0, 0.2+t_res, t_res)
v_calib = numpy.zeros(numpy.shape(t_calib))
for i in range(0, len(t_calib)):
if t_calib[i] < 0.04:
v_calib[i] = 0.0
elif t_calib[i] > t_calib[-1]-0.04:
v_calib[i] = 0.0
else:
v_calib[i] = 1.0
# Concatenate the leads together to plot in a single figure
t_end = ecgs['ts'][beat-1][-1]
full_t = numpy.concatenate([t_calib, ecgs['ts'][beat-1]+t_calib[-1], ecgs['ts'][beat-1]+t_end+t_calib[-1], ecgs['ts'][beat-1]+2*t_end+t_calib[-1], ecgs['ts'][beat-1]+3*t_end+t_calib[-1]]) + t_padding
scale = 1
fig_size = [numpy.ceil(full_t[-1]+t_padding)*0.5/0.2*scale,6*scale]
fig = plt.figure(tight_layout=True, figsize=fig_size)
gs = GridSpec(1,1)
ax = fig.add_subplot(gs[0,0])
# Bottom row: III, aVF, V3, V6
bottom_V = numpy.concatenate([v_calib, ecgs['IIIs'][beat-1]/ecgs['max_limb_leads'], ecgs['aVFs'][beat-1]/ecgs['max_limb_leads'], ecgs['V3s'][beat-1]/ecgs['max_precord_leads'], ecgs['V6s'][beat-1]/ecgs['max_precord_leads']])
ax.plot(full_t, bottom_V, 'k')
# Middle row: II, aVL, V2, V5
offset = 2
midrow_V = numpy.concatenate([v_calib+offset, ecgs['IIs'][beat-1]/ecgs['max_limb_leads']+offset, ecgs['aVLs'][beat-1]/ecgs['max_limb_leads']+offset, ecgs['V2s'][beat-1]/ecgs['max_precord_leads']+offset, ecgs['V5s'][beat-1]/ecgs['max_precord_leads']+offset])
ax.plot(full_t, midrow_V, 'k')
# Top row: I, aVR, V1, V4
toprow_V = numpy.concatenate([v_calib+offset*2, ecgs['Is'][beat-1]/ecgs['max_limb_leads']+offset*2, ecgs['aVRs'][beat-1]/ecgs['max_limb_leads']+offset*2, ecgs['V1s'][beat-1]/ecgs['max_precord_leads']+offset*2, ecgs['V4s'][beat-1]/ecgs['max_precord_leads']+offset*2])
ax.plot(full_t, toprow_V, 'k')
# ax.set_xlim([0, t_end*4])
# ax.set_ylim([-1, offset*2+1])
#
# self._set_full_ecg_ticks(ax, full_t[-1]+t_padding, offset*2+1, self.CL)
# Add ECG red grid
self._set_full_ecg_ticks(ax, full_t[-1]+t_padding*2, offset*2 + 1 + v_padding, -1-v_padding, self.CL)
ax.set_xlim([0, full_t[-1]+t_padding*2])
ax.set_ylim([-1-v_padding, offset*2 + 1 +v_padding])
# Label the leads
label_fontsize = 'large'
label_t_offset = t_calib[-1] + t_padding
label_v_offset = 0.1
plt.text(0+label_t_offset, -1+label_v_offset, 'III', fontsize=label_fontsize)
plt.text(t_end+label_t_offset, -1+label_v_offset, 'aVF', fontsize=label_fontsize)
plt.text(2*t_end+label_t_offset, -1+label_v_offset, 'V3', fontsize=label_fontsize)
plt.text(3*t_end+label_t_offset, -1+label_v_offset, 'V6', fontsize=label_fontsize)
plt.text(0+label_t_offset, -1+offset+label_v_offset, 'II', fontsize=label_fontsize)
plt.text(t_end+label_t_offset, -1+offset+label_v_offset, 'aVL', fontsize=label_fontsize)
plt.text(2*t_end+label_t_offset, -1+offset+label_v_offset, 'V2', fontsize=label_fontsize)
plt.text(3*t_end+label_t_offset, -1+offset+label_v_offset, 'V5', fontsize=label_fontsize)
plt.text(0+label_t_offset, -1+offset*2+label_v_offset, 'I', fontsize=label_fontsize)
plt.text(t_end+label_t_offset, -1+offset*2+label_v_offset, 'aVR', fontsize=label_fontsize)
plt.text(2*t_end+label_t_offset, -1+offset*2+label_v_offset, 'V1', fontsize=label_fontsize)
plt.text(3*t_end+label_t_offset, -1+offset*2+label_v_offset, 'V4', fontsize=label_fontsize)
plt.savefig('full_ecg_beat_'+str(beat)+'.png')
if show:
plt.show()
else:
# Padding
t_padding = 0.2
v_padding = 0.5
# Generate calibrating step function
t_res = ecgs['t'][1] - ecgs['t'][0]
t_calib = numpy.arange(0, 0.2+t_res, t_res)
v_calib = numpy.zeros(numpy.shape(t_calib))
for i in range(0, len(t_calib)):
if t_calib[i] < 0.04:
v_calib[i] = 0.0
elif t_calib[i] > t_calib[-1]-0.04:
v_calib[i] = 0.0
else:
v_calib[i] = 1.0
# Concatenate the leads together to plot in a single figure
t_end = ecgs['t'][-1]
full_t = numpy.concatenate([t_calib, ecgs['t']+t_calib[-1], ecgs['t']+t_end+t_calib[-1], ecgs['t']+2*t_end+t_calib[-1], ecgs['t']+3*t_end+t_calib[-1]]) + t_padding
scale = 0.8
fig_size = [numpy.ceil(full_t[-1]+t_padding)*0.5/0.2*scale,6*scale]
fig = plt.figure(tight_layout=True, figsize=fig_size)
gs = GridSpec(1,1)
ax = fig.add_subplot(gs[0,0])
# Bottom row: III, aVF, V3, V6
bottom_V = numpy.concatenate([v_calib, ecgs['III']/ecgs['max_limb_leads'], ecgs['aVF']/ecgs['max_limb_leads'], ecgs['V3']/ecgs['max_precord_leads'], ecgs['V6']/ecgs['max_precord_leads']])
ax.plot(full_t, bottom_V, 'k')
# Middle row: II, aVL, V2, V5
offset = 2
midrow_V = numpy.concatenate([v_calib+offset, ecgs['II']/ecgs['max_limb_leads']+offset, ecgs['aVL']/ecgs['max_limb_leads']+offset, ecgs['V2']/ecgs['max_precord_leads']+offset, ecgs['V5']/ecgs['max_precord_leads']+offset])
ax.plot(full_t, midrow_V, 'k')
# Top row: I, aVR, V1, V4
toprow_V = numpy.concatenate([v_calib+offset*2, ecgs['I']/ecgs['max_limb_leads']+offset*2, ecgs['aVR']/ecgs['max_limb_leads']+offset*2, ecgs['V1']/ecgs['max_precord_leads']+offset*2, ecgs['V4']/ecgs['max_precord_leads']+offset*2])
ax.plot(full_t, toprow_V, 'k')
# Add ECG red grid
self._set_full_ecg_ticks(ax, full_t[-1]+t_padding*2, offset*2 + 1 + v_padding, -1-v_padding, self.CL)
ax.set_xlim([0, full_t[-1]+t_padding*2])
ax.set_ylim([-1-v_padding, offset*2 + 1 +v_padding])
# Label the leads
label_fontsize = 'xx-small'
label_t_offset = t_calib[-1] + t_padding
label_v_offset = 0.1
plt.text(0+label_t_offset, -1+label_v_offset, 'III', fontsize=label_fontsize)
plt.text(t_end+label_t_offset, -1+label_v_offset, 'aVF', fontsize=label_fontsize)
plt.text(2*t_end+label_t_offset, -1+label_v_offset, 'V3', fontsize=label_fontsize)
plt.text(3*t_end+label_t_offset, -1+label_v_offset, 'V6', fontsize=label_fontsize)
plt.text(0+label_t_offset, -1+offset+label_v_offset, 'II', fontsize=label_fontsize)
plt.text(t_end+label_t_offset, -1+offset+label_v_offset, 'aVL', fontsize=label_fontsize)
plt.text(2*t_end+label_t_offset, -1+offset+label_v_offset, 'V2', fontsize=label_fontsize)
plt.text(3*t_end+label_t_offset, -1+offset+label_v_offset, 'V5', fontsize=label_fontsize)
plt.text(0+label_t_offset, -1+offset*2+label_v_offset, 'I', fontsize=label_fontsize)
plt.text(t_end+label_t_offset, -1+offset*2+label_v_offset, 'aVR', fontsize=label_fontsize)
plt.text(2*t_end+label_t_offset, -1+offset*2+label_v_offset, 'V1', fontsize=label_fontsize)
plt.text(3*t_end+label_t_offset, -1+offset*2+label_v_offset, 'V4', fontsize=label_fontsize)
plt.savefig('full_ecg_all_beats.png')
if show:
plt.show()
def plot_pv_signal(self, pvs, signal_name, filename, show, beat=0, pvs2=[]):
matplotlib.rcParams.update({'font.size':'28'})
matplotlib.rcParams.update({'text.color':'black'})
matplotlib.rcParams.update({'lines.linewidth':'3'})
if pvs2:
if beat > 0:
if (signal_name == 'p') | (signal_name == 'v'):
self._plot_quadruple(self.pv_fig_size,
pvs['ts'][beat-1], pvs[signal_name+'ls'][beat-1],'b--',
pvs['ts'][beat-1], pvs[signal_name+'rs'][beat-1],'g--',
pvs2['ts'][beat-1], pvs2[signal_name+'ls'][beat-1],'b-',
pvs2['ts'][beat-1], pvs2[signal_name+'rs'][beat-1],'g-',
'Time (s)',pvs[signal_name+'label'], filename+'_beat'+str(beat)+'_compare.png', show)
elif (signal_name == 'pv'):
self._plot_quadruple(self.pv_fig_size,
pvs['vls'][beat-1],pvs['pls'][beat-1]*7.5,'b--',
pvs['vrs'][beat-1],pvs['prs'][beat-1]*7.5,'g--',
pvs2['vls'][beat-1],pvs2['pls'][beat-1]*7.5,'b-',
pvs2['vrs'][beat-1],pvs2['prs'][beat-1]*7.5,'g-',
'Volume (mL)', 'Pressure (mmHg)',
filename+'_beat'+str(beat)+'_compare.png', show)
else:
if (signal_name == 'p') | (signal_name == 'v'):
figsize=[self.pv_fig_size[0]*(int(pvs['t'][-1]/1.0)+1), self.pv_fig_size[1]]
self._plot_quadruple(figsize,
pvs['t'], pvs[signal_name+'l'],'b--',
pvs['t'], pvs[signal_name+'r'],'g--',
pvs2['t'], pvs2[signal_name+'l'],'b-',
pvs2['t'], pvs2[signal_name+'r'],'g-','Time (s)',
pvs[signal_name+'label'], filename+'_full_compare.png', show)
elif (signal_name == 'pv'):
self._plot_quadruple(self.pv_fig_size,
pvs['vl'], pvs['pl']*7.5,'b--',
pvs['vr'], pvs['pr']*7.5, 'g--',
pvs2['vl'], pvs2['pl']*7.5,'b-',
pvs2['vr'], pvs2['pr']*7.5, 'g-',
'Volume (mL)','Pressure (mmHg)', filename+'_full_compare.png', show)
else:
if beat > 0:
if (signal_name == 'p') | (signal_name == 'v'):
self._plot_double(self.pv_fig_size,
pvs['ts'][beat-1],pvs[signal_name+'ls'][beat-1],'b-',
pvs['ts'][beat-1],pvs[signal_name+'rs'][beat-1],'g-',
'Time (s)',pvs[signal_name+'label'], filename+'_beat'+str(beat)+'.png',show)
elif (signal_name == 'pv'):
self._plot_double(self.pv_fig_size,
pvs['vls'][beat-1],pvs['pls'][beat-1], 'b-',
pvs['vrs'][beat-1],pvs['prs'][beat-1], 'g-',
'Volume (mL)', 'Pressure (kPa)',
filename+'_beat'+str(beat)+'.png',show)
else:
if (signal_name == 'p') | (signal_name == 'v'):
figsize=[self.pv_fig_size[0]*(int(pvs['t'][-1]/1.0)+1), self.pv_fig_size[1]]
self._plot_double(figsize, pvs['t'], pvs[signal_name+'l'],'b-',
pvs['t'], pvs[signal_name+'r'],'g-',
'Time (s)',pvs[signal_name+'label'], filename+'_full.png',show)
elif (signal_name == 'pv'):
self._plot_double(self.pv_fig_size, pvs['vl'], pvs['pl'],'b-',
pvs['vr'], pvs['pr'], 'r-',
'Volume (mL)','Pressure (kPa)', filename+'_full.png',show)
def _plot_single(self, fig_size, x, y, xlabel, ylabel, filename, show, ecg_grid=False):
fig = plt.figure(tight_layout=True, figsize=fig_size)
gs = GridSpec(1,1)
ax = fig.add_subplot(gs[0,0])
ax.plot(x, y)
ax.set_xlabel(xlabel)
ax.set_ylabel(ylabel)
if ecg_grid:
print('setting ecg grids')
self._set_ecg_ticks(ax, x[-1], self.CL)
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.spines['bottom'].set_visible(False)
ax.spines['left'].set_visible(False)
ax.tick_params(left=False, right=False, top=False, bottom=False, labelleft=False, labelbottom=False)
ax.set_xlim([0, numpy.ceil(x[-1]/self.CL)*CL])
if show:
plt.show(block=True)
plt.savefig(filename)
def _plot_double(self, fig_size, x, y, linestyle, x1, y1, linestyle1, xlabel, ylabel, filename, show, ecg_grid=False):
fig = plt.figure(tight_layout=True, figsize=fig_size)
gs = GridSpec(1,1)
ax = fig.add_subplot(gs[0,0])
ax.plot(x, y, linestyle, x1, y1, linestyle1)
ax.set_xlabel(xlabel)
ax.set_ylabel(ylabel)
if ecg_grid:
self._set_ecg_ticks(ax, x[-1],self.CL)
if show:
plt.show(block=True)
plt.savefig(filename)
def _plot_quadruple(self, fig_size, x, y, linestyle, x1, y1,linestyle1, x2, y2, linestyle2, x3, y3,linestyle3, xlabel, ylabel, filename, show, ecg_grid=False):
fig = plt.figure(tight_layout=True, figsize=fig_size)
gs = GridSpec(1,1)
ax = fig.add_subplot(gs[0,0])
ax.plot(x, y, linestyle,x1, y1,linestyle1, x2, y2, linestyle2, x3, y3,linestyle3)
ax.set_xlabel(xlabel)
ax.set_ylabel(ylabel)
ax.set_ylim([0,120])
self._set_pv_ticks(ax)
if ecg_grid:
self._set_ecg_ticks(ax, x[-1], self.CL)
if show:
plt.show(block=True)
plt.savefig(filename)
def read_single_cell(self, dir, material):
# Read membrane potentials
filename = dir+'torord_ode.m'+str(material)+'c3.csv'
with open(filename, 'r') as f:
data = f.readlines()
data = data[1:]
epi = numpy.zeros((len(data),5))
for i in range(0, len(epi)): epi[i,1] = float(data[i].split()[2])*1000000.0 # Calcium transient
for i in range(0, len(epi)): epi[i,2] = float(data[i].split()[3]) # Membrane potential
for i in range(0, len(epi)): epi[i,0] = float(data[i].split()[0])/1000.0 # real time
for i in range(0, len(epi)): epi[i,3] = float(data[i].split()[1])/1000.0 # each beat time
for i in range(0, len(epi)): epi[i,4] = float(data[i].split()[4]) # Ta
filename = dir+'torord_ode.m'+str(material)+'c2.csv'
with open(filename, 'r') as f:
data = f.readlines()
data = data[1:]
mid = numpy.zeros((len(data),5))
for i in range(0, len(mid)): mid[i,1] = float(data[i].split()[2])*1000000.0 # Calcium transient
for i in range(0, len(mid)): mid[i,2] = float(data[i].split()[3]) # Membrane potential
for i in range(0, len(mid)): mid[i,0] = float(data[i].split()[0])/1000.0 # real time
for i in range(0, len(mid)): mid[i,3] = float(data[i].split()[1])/1000.0 # each beat time
for i in range(0, len(mid)): mid[i,4] = float(data[i].split()[4]) # Ta
filename = dir+'torord_ode.m'+str(material)+'c1.csv'
with open(filename, 'r') as f:
data = f.readlines()
data = data[1:]
endo = numpy.zeros((len(data),5))
for i in range(0, len(endo)): endo[i,1] = float(data[i].split()[2])*1000000.0 # Calcium transient
for i in range(0, len(endo)): endo[i,2] = float(data[i].split()[3]) # Membrane potential
for i in range(0, len(endo)): endo[i,0] = float(data[i].split()[0])/1000.0 # real time
for i in range(0, len(endo)): endo[i,3] = float(data[i].split()[1])/1000.0 # each beat time
for i in range(0, len(endo)): endo[i,4] = float(data[i].split()[4]) # Ta
return epi, mid, endo
def plot_single_cell(self, epi, mid, endo, material,analysis, show, epi2=[], mid2=[], endo2=[]):
# Set up subplots
fig = plt.figure(tight_layout=True, figsize=[15,5])
gs = GridSpec(3,7)
ax1 = fig.add_subplot(gs[0,0])
ax2 = fig.add_subplot(gs[1,0])
ax3 = fig.add_subplot(gs[2,0])
ax4 = fig.add_subplot(gs[:,1:3])
ax5 = fig.add_subplot(gs[:,3:5])
ax6 = fig.add_subplot(gs[:,5:7])
if len(epi2)>0:
ax1.plot(epi[:,0], epi[:,2], mid[:,0], mid[:,2], endo[:,0], endo[:,2],
epi2[:,0], epi2[:,2],'--', mid2[:,0], mid2[:,2],'--', endo2[:,0], endo2[:,2], '--')
ax2.plot(epi[:,0], epi[:,1],mid[:,0],mid[:,1],endo[:,0],endo[:,1],
epi2[:,0], epi2[:,1],'--',mid2[:,0],mid2[:,1],'--',endo2[:,0],endo2[:,1],'--')
ax3.plot(epi[:,0], epi[:,4],mid[:,0],mid[:,4],endo[:,0],endo[:,4],
epi2[:,0], epi2[:,4],'--',mid2[:,0],mid2[:,4],'--',endo2[:,0],endo2[:,4],'--')
dt = epi[1,0] - epi[0,0]
idx = int(self.CL / dt)
t_s_epi = epi[-idx,0]
t_s_mid = mid[-idx,0]
t_s_endo = endo[-idx,0]
t_s_epi2 = epi2[-idx,0]
t_s_mid2 = mid2[-idx,0]
t_s_endo2 = endo2[-idx,0]
ax4.plot(epi[-idx:,0]-t_s_epi, epi[-idx:,2], mid[-idx:,0]-t_s_mid, mid[-idx:,2], endo[-idx:,0]-t_s_endo, endo[-idx:,2],
epi2[-idx:,0]-t_s_epi2, epi2[-idx:,2],'--', mid2[-idx:,0]-t_s_mid2, mid2[-idx:,2],'--', endo2[-idx:,0]-t_s_endo2, endo2[-idx:,2], '--')
ax5.plot(epi[-idx:,0]-t_s_epi, epi[-idx:,1], mid[-idx:,0]-t_s_mid, mid[-idx:,1], endo[-idx:,0]-t_s_endo, endo[-idx:,1],
epi2[-idx:,0]-t_s_epi2, epi2[-idx:,1],'--',mid2[-idx:,0]-t_s_mid2, mid2[-idx:,1],'--', endo2[-idx:,0]-t_s_endo2, endo2[-idx:,1],'--')
ax6.plot(epi[-idx:,0]-t_s_epi, epi[-idx:,4], mid[-idx:,0]-t_s_mid, mid[-idx:,4], endo[-idx:,0]-t_s_endo, endo[-idx:,4],
epi2[-idx:,0]-t_s_epi2, epi2[-idx:,4],'--',mid2[-idx:,0]-t_s_mid2, mid2[-idx:,4],'--', endo2[-idx:,0]-t_s_endo2, endo2[-idx:,4],'--')
if analysis:
print ('1: Diastolic calcium (nM) epi: '+str(epi[-1,1])+', mid: '+str(mid[-1,1])+', endo: '+str(endo[-1,1]))
print ('1: Calcium amplitude (nM) epi: '+str(max(epi[-idx:,1]))+', mid: '+str(max(mid[-idx:,1]))+', endo: '+str(max(endo[-idx:,1])))
print ('1: Diastolic active tension (kPa) epi: '+str(epi[-1,4])+', mid: '+str(mid[-1,4])+', endo: '+str(endo[-1,4]))
print ('1: Active tension amplitude (kPa) epi: '+str(max(epi[-idx:,4]))+', mid: '+str(max(mid[-idx:,4]))+', endo: '+str(max(endo[-idx:,4])))
print ('2: Diastolic calcium (nM) epi: '+str(epi2[-1,1])+', mid: '+str(mid2[-1,1])+', endo: '+str(endo2[-1,1]))
print ('2: Calcium amplitude (nM) epi: '+str(max(epi2[-idx:,1]))+', mid: '+str(max(mid2[-idx:,1]))+', endo: '+str(max(endo2[-idx:,1])))
print ('2: Diastolic active tension (kPa) epi: '+str(epi2[-1,4])+', mid: '+str(mid2[-1,4])+', endo: '+str(endo2[-1,4]))
print ('2: Active tension amplitude (kPa) epi: '+str(max(epi2[-idx:,4]))+', mid: '+str(max(mid2[-idx:,4]))+', endo: '+str(max(endo2[-idx:,4])))
else:
ax1.plot(epi[:,0], epi[:,2], mid[:,0], mid[:,2], endo[:,0], endo[:,2])
ax2.plot(epi[:,0],epi[:,1],mid[:,0],mid[:,1],endo[:,0],endo[:,1])
ax3.plot(epi[:,0],epi[:,4],mid[:,0],mid[:,4],endo[:,0],endo[:,4])
dt = epi[1,0] - epi[0,0]
idx = int(self.CL / dt)
t_s_epi = epi[-idx,0]
t_s_mid = mid[-idx,0]
t_s_endo = endo[-idx,0]
ax4.plot(epi[-idx:,0]-t_s_epi, epi[-idx:,2], mid[-idx:,0]-t_s_mid, mid[-idx:,2], endo[-idx:,0]-t_s_endo, endo[-idx:,2])
ax5.plot(epi[-idx:,0]-t_s_epi, epi[-idx:,1], mid[-idx:,0]-t_s_mid, mid[-idx:,1], endo[-idx:,0]-t_s_endo, endo[-idx:,1])
ax6.plot(epi[-idx:,0]-t_s_epi, epi[-idx:,4], mid[-idx:,0]-t_s_mid, mid[-idx:,4], endo[-idx:,0]-t_s_endo, endo[-idx:,4])
if analysis:
print ('Diastolic calcium (nM) epi: {:.2f}, mid: {:.2f}, endo: {:.2f}'.format(epi[-1,1], mid[-1,1], endo[-1,1]))
print ('Calcium amplitude (nM) epi: {:.2f}, mid: {:.2f}, endo: {:.2f}'.format(max(epi[-idx:,1]), max(mid[-idx:,1]), max(endo[-idx:,1])))
print ('Diastolic active tension (kPa) epi: {:.2f}, mid: {:.2f}, endo: {:.2f}'.format(epi[-1,4], mid[-1,4], endo[-1,4]))
print ('Active tension amplitude (kPa) epi: {:.2f}, mid: {:.2f}, endo: {:.2f}'.format(max(epi[-idx:,4]), max(mid[-idx:,4]), max(endo[-idx:,4])))
ax1.set_ylabel('(mV)')
ax1.set_xlabel('Time (s)')
#ax1.set_title('Action potentials')
ax1.legend(['Epi', 'Mid', 'Endo'])
#ax2.set_ylabel('Intracellular Calcium (nM)')
ax2.set_xlabel('Time (s)')
ax2.set_ylabel('CaT (nM)')
ax3.set_xlabel('Time (s)')
ax3.set_ylabel('Active tension (kPa)')
dt = epi[1,0] - epi[0,0]
idx = int(0.8 / dt)
ax4.set_ylabel('Vm (mV)')
ax4.set_xlabel('Time (s)')
ax4.set_ylim([-100, 40])
ax4.grid()
ax4.set_xticks(numpy.arange(0, self.CL+0.2, 0.2))
ax5.set_ylabel('CaT (nM)')
ax5.set_xlabel('Time (s)')
ax5.set_ylim([0, 1000])
ax5.grid()
ax5.set_xticks(numpy.arange(0, self.CL+0.2, 0.2))
ax6.set_ylabel('Active tension (kPa)')
ax6.set_xlabel('Time (s)')
ax6.set_ylim([0, 70])
ax6.grid()
ax6.set_xticks(numpy.arange(0, self.CL+0.2, 0.2))
if show:
plt.show(block=True)
plt.savefig('single_cell_'+str(material)+'.png')
def analysis_PV(self, pvs, beat=0):
if beat > 0:
EDVL = max(pvs['vls'][beat-1])
EDVR = max(pvs['vrs'][beat-1])
ESVL = min(pvs['vls'][beat-1])
ESVR = min(pvs['vrs'][beat-1])
PmaxL = int(max(pvs['pls'][beat-1]))
PmaxR = int(max(pvs['prs'][beat-1]))
LVEF = int((EDVL-ESVL)/EDVL*100)
RVEF = int((EDVR-ESVR)/EDVR*100)
SVL = int(EDVL - ESVL)
SVR = int(EDVR - ESVR)
else:
LVEF = []
RVEF = []
EDVL = []
EDVR = []
ESVL = []
ESVR = []
PmaxL = []
PmaxR = []
SVL = []
SVR = []
for i in range(0, len(pvs['vls'])):
EDVL.append(max(pvs['vls'][i]))
ESVL.append(min(pvs['vls'][i]))
LVEF.append(int((EDVL[i]-ESVL[i])/EDVL[i] * 100))
PmaxL.append(int(max(pvs['pls'][i])))
SVL.append(int(EDVL[i]-ESVL[i]))
EDVR.append(max(pvs['vrs'][i]))
ESVR.append(min(pvs['vrs'][i]))
RVEF.append(int((EDVR[i]-ESVR[i])/EDVR[i] * 100))
PmaxR.append(int(max(pvs['prs'][i])))
SVR.append(int(EDVR[i]-ESVR[i]))
print('EDVL:'+str(int(EDVL))+' mL,\tEDVR: '+str(int(EDVR))+' mL')
print('ESVL: '+str(int(ESVL))+' mL,\tESVR: '+str(int(ESVR))+' mL')
print('LVEF: '+str(LVEF)+' %,\tRVEF: '+str(RVEF)+' %')
print('PmaxL: '+str(PmaxL)+' kPa,\tPmaxR: '+str(PmaxR)+' kPa')
print('SVL: '+str(SVL)+' mL,\tSVR: '+str(SVR)+' mL')
def analysis_ECG(self,ecgs,beat,show):
fig = plt.figure(tight_layout=True, figsize=[11,9])
gs = GridSpec(3,4)
ax5 = fig.add_subplot(gs[0,0])
ax6 = fig.add_subplot(gs[1,0])
ax7 = fig.add_subplot(gs[2,0])
ax8 = fig.add_subplot(gs[0,1])
ax9 = fig.add_subplot(gs[1,1])
ax10 = fig.add_subplot(gs[2,1])
ax11 = fig.add_subplot(gs[0,2])
ax12 = fig.add_subplot(gs[1,2])
ax13 = fig.add_subplot(gs[2,2])
ax14 = fig.add_subplot(gs[0,3])
ax15 = fig.add_subplot(gs[1,3])
ax16 = fig.add_subplot(gs[2,3])
ecg_biomarkers = numpy.zeros((12,6))
ecg_biomarkers[0, :] = self._plot_with_landmarks(ax5, ecgs['ts'][beat-1], ecgs['max_all_leads'], ecgs['Is'][beat-1], 'I', CL)
ecg_biomarkers[1, :] = self._plot_with_landmarks(ax6, ecgs['ts'][beat-1], ecgs['max_all_leads'], ecgs['IIs'][beat-1], 'II', CL)
ecg_biomarkers[2, :] = self._plot_with_landmarks(ax7, ecgs['ts'][beat-1], ecgs['max_all_leads'], ecgs['IIIs'][beat-1], 'III', CL)
ecg_biomarkers[3, :] = self._plot_with_landmarks(ax8, ecgs['ts'][beat-1], ecgs['max_all_leads'], ecgs['aVRs'][beat-1], 'aVR', CL)
ecg_biomarkers[4, :] = self._plot_with_landmarks(ax9, ecgs['ts'][beat-1], ecgs['max_all_leads'], ecgs['aVLs'][beat-1], 'aVL', CL)
ecg_biomarkers[5, :] = self._plot_with_landmarks(ax10, ecgs['ts'][beat-1], ecgs['max_all_leads'], ecgs['aVFs'][beat-1], 'aVF', CL)
ecg_biomarkers[6, :] = self._plot_with_landmarks(ax11, ecgs['ts'][beat-1], ecgs['max_all_leads'], ecgs['V1s'][beat-1], 'V1', CL)
ecg_biomarkers[7, :] = self._plot_with_landmarks(ax12, ecgs['ts'][beat-1], ecgs['max_all_leads'], ecgs['V2s'][beat-1], 'V2', CL)
ecg_biomarkers[8, :] = self._plot_with_landmarks(ax13, ecgs['ts'][beat-1], ecgs['max_all_leads'], ecgs['V3s'][beat-1], 'V3', CL)
ecg_biomarkers[9, :] = self._plot_with_landmarks(ax14, ecgs['ts'][beat-1], ecgs['max_all_leads'], ecgs['V4s'][beat-1], 'V4', CL)
ecg_biomarkers[10, :] = self._plot_with_landmarks(ax15, ecgs['ts'][beat-1], ecgs['max_all_leads'], ecgs['V5s'][beat-1], 'V5', CL)
ecg_biomarkers[11, :] = self._plot_with_landmarks(ax16, ecgs['ts'][beat-1], ecgs['max_all_leads'], ecgs['V6s'][beat-1], 'V6', CL)
if show:
plt.show()
plt.savefig('ecg_analysis.png')
print('QRS dur, T dur, T pe, T op, T amp, QT dur')
print (ecg_biomarkers)
QRS_mean = numpy.average(ecg_biomarkers[:,0])
QRS_std = numpy.std(ecg_biomarkers[:,0])
T_dur_mean = numpy.average(ecg_biomarkers[:,1])
T_dur_std = numpy.std(ecg_biomarkers[:,1])
T_pe_mean = numpy.average(ecg_biomarkers[:,2])
T_pe_std = numpy.std(ecg_biomarkers[:,2])
T_op_mean = numpy.average(ecg_biomarkers[:,3])
T_op_std = numpy.std(ecg_biomarkers[:,3])
QT_dur_mean = numpy.average(ecg_biomarkers[:,5])
QT_dur_std = numpy.std(ecg_biomarkers[:,5])
QT_dispersion_6 = ecg_biomarkers[6:12,5].max() - ecg_biomarkers[6:12,5].min()
QT_dispersion_12 = ecg_biomarkers[:,5].max() - ecg_biomarkers[:,5].min()
print('QRS: {:.3f} +- {:.3f}'.format(QRS_mean, QRS_std))
print('T duration: {:.3f} +- {:.3f}'.format(T_dur_mean, T_dur_std))
print('T pe: {:.3f} +- {:.3f}'.format(T_pe_mean, T_pe_std))
print('T op: {:.3f} +- {:.3f}'.format(T_op_mean, T_op_std))
print('QT duration: {:.3f} +- {:.3f}'.format(QT_dur_mean, QT_dur_std))
print('QT dispersion (precordial): {:.3f}'.format(QT_dispersion_6))
print('QT dispersion (12 leads): {:.3f}'.format(QT_dispersion_12))
def _plot_with_landmarks(self, ax, t, max_all_leads, V, name, CL):
ax.clear()
qrs_dur, t_dur, t_amp, t_ep, t_op, qt_dur, landmarks =self._measurements(V, 0.0, CL, t, 2e-5, 0.02)
ax.plot(t, V/max_all_leads, landmarks[:,0], landmarks[:,1]/max_all_leads, '*')
ax.set_title(name + ' '+str(int(qrs_dur*1000))+' '+str(int(qt_dur*1000))+'\n'+str(int(t_dur*1000))+' '+str(int(t_amp/max_all_leads)))
ax.set_xlim(0.0,CL)
ax.set_ylim(-1,1)
ax.set_xlabel('Time (s)')
ax.set_ylabel('Normalised ECG')
return qrs_dur, t_dur, t_amp, t_ep, t_op, qt_dur
def _measurements(self, V, start_t, end_t, t, t_tol, v_tol):
idx_start = numpy.where(abs(numpy.array(t-start_t)) < t_tol)[0][0]
try:
idx_end = numpy.where(abs(numpy.array(t-end_t)) < t_tol)[0][0]
except:
idx_end = len(V)-1
# Offset voltage using t_end
n = 100
b = [1.0 / n] * n
V = V - V[idx_end]
V = lfilter(b,1, V)
dV = numpy.gradient(V)
dV[0:2] = 0.0 # Remove gradient artefacts
n = 500
b = [1.0 / n] * n
dV = lfilter(b, 1, dV)
dV = abs(dV)/abs(dV).max()
ddV = numpy.gradient(dV)
ddV[0:2] = 0.0 # Remove gradient artefacts
ddV = lfilter(b, 1, ddV)
ddV = abs(ddV)/abs(ddV).max()
TOL = v_tol
# Find Q start
for i in range(idx_start, len(V)):
if (dV[i] > TOL) & (i > 10):
break
q_start_idx = i
q_start_t = t[i]
# Find T end
for i in range(idx_end-1, idx_start, -1):
if (dV[i] > TOL):
break
t_end_idx = i
t_end_t = t[i]
# Find QRS end
window = 5000
dV = abs(dV)
max_window = numpy.zeros(t_end_idx-q_start_idx+1)
for i in range(q_start_idx, t_end_idx):
max_window[i-q_start_idx] = dV[i:(i+window)].max()
print('evaluated max window')
min = max_window[0]
for i in range(0, len(max_window)):
if (min >= max_window[i]):
min = max_window[i]
elif (min < (max_window.min()+0.05)):
break
qrs_end_idx = i + q_start_idx
qrs_end_t = t[qrs_end_idx]
qrs_dur = qrs_end_t - q_start_t
# Find T peak timing and amplitude
segment = V[qrs_end_idx:(t_end_idx+1)]
t_amplitude = abs(segment).max()
t_peak_idx = numpy.where(abs(segment) == t_amplitude)[0][0] + qrs_end_idx
t_sign = numpy.sign(segment[t_peak_idx-qrs_end_idx])
t_amplitude = t_sign * t_amplitude
t_peak_t = t[t_peak_idx]
# Find T start
segment = ddV[qrs_end_idx:t_peak_idx]
min_dd_idx = numpy.where(segment == segment.min())[0][0] + qrs_end_idx
for i in range(min_dd_idx, t_peak_idx):
if (abs(ddV[i]) > 0.0005):
break
t_start_idx = i
t_start_t = t[t_start_idx]
t_dur = t_end_t - t_start_t
qt_dur = t_end_t - q_start_t
t_pe = t_end_t - t_peak_t
t_op = t_peak_t - t_start_t
landmarks = numpy.array([[q_start_t, V[q_start_idx]], [qrs_end_t, V[qrs_end_idx]], [t_start_t, V[t_start_idx]], [t_peak_t, V[t_peak_idx]], [t_end_t, V[t_end_idx]]])
return qrs_dur, t_dur, t_pe, t_op, t_amplitude, qt_dur, landmarks
if __name__=='__main__':
false = ['False', 'false', 'F', 'f']
true = ['True', 'true', 'T', 't']
parser = argparse.ArgumentParser()
parser.add_argument("--type", help='Either ecg or pv, or single cell')
parser.add_argument("--name", help='Name of simulation', default='heart_remeshed_3D')
parser.add_argument("--refresh", help='Delete previous reading of ECG and PV', default=False)
parser.add_argument("--lead", help='Name of the signal lead to plot')
parser.add_argument("--CL", help='Cycle length of simulation, (s)',default=0.8)
parser.add_argument("--beat", help='Beat number for single beat plots', default=0)
parser.add_argument("--figure_title", help='Title of the figure', default = '')
parser.add_argument("--analysis", help='Title of the figure', default =False)
parser.add_argument("--compare", help='Directory to compare with', default='')
parser.add_argument("--material", help='Directory to compare with', default=1)
parser.add_argument("--show", help='Toggle whether to show plot', default=True)
args = parser.parse_args()
plot_type = args.type
name = args.name
refresh = bool(args.refresh)
if refresh in true:
os.system('rm ecgs.pl pvs.pl')
lead_name = args.lead
CL = float(args.CL)
beat = int(args.beat)
figure_title = args.figure_title
analysis = args.analysis
compare = args.compare
material = int(args.material)
show = args.show
if show in false:
show = False
if analysis in true:
analysis = True
# Set up visualisation
a = ECGPV_visualisation(CL)
if not (plot_type=='cell'):
ecgs, pvs = a.read_ecg_pv(name, './')
epi, mid, endo = a.read_single_cell('./', material)
if compare:
ecgs2, pvs2 = a.read_ecg_pv(name, compare)
epi2, mid2, endo2 = a.read_single_cell(compare, material)
if analysis:
if (plot_type=='live')|(plot_type=='p')|(plot_type=='v')|(plot_type=='pv'):
a.analysis_PV(pvs, beat)
print('second one...')
a.analysis_PV(pvs2, beat)
elif (plot_type=='ecg'):
a.analysis_ECG(ecgs,beat,show)
a.analysis_ECG(ecgs2,beat,show)
if plot_type == 'live':
a.plot_ecgpv_live(ecgs, pvs, figure_title, show, ecgs2, pvs2)
elif plot_type == 'ecg':
if lead_name:
a.plot_ecg_lead(ecgs, lead_name, lead_name, show, beat, ecgs2)
else:
a.plot_ecg_all_leads(ecgs, show, beat, ecgs2)
elif plot_type == 'p':
a.plot_pv_signal(pvs, 'p', 'Pt', show, beat, pvs2)
elif plot_type == 'v':
a.plot_pv_signal(pvs, 'v', 'Vt', show, beat, pvs2)
elif plot_type == 'pv':
a.plot_pv_signal(pvs, 'pv', 'PV_loop', show, beat, pvs2)
elif plot_type == 'cell':
a.plot_single_cell(epi, mid, endo, material,analysis, show, epi2, mid2, endo2)
else:
if analysis:
if (plot_type=='live')|(plot_type=='p')|(plot_type=='v')|(plot_type=='pv'):
a.analysis_PV(pvs, beat)
elif (plot_type=='ecg'):
a.analysis_ECG(ecgs,beat,show)
if plot_type == 'live':
a.plot_ecgpv_live(ecgs, pvs, figure_title, show)
elif plot_type == 'ecg':
if lead_name:
assert lead_name, 'Please enter lead name using option: lead='
a.plot_ecg_lead(ecgs, lead_name, lead_name, show, beat)
else:
a.plot_ecg_all_leads(ecgs, show, beat)
elif plot_type == 'p':
a.plot_pv_signal(pvs, 'p', 'Pt', show, beat)
elif plot_type == 'v':
a.plot_pv_signal(pvs, 'v', 'Vt', show, beat)
elif plot_type == 'pv':
a.plot_pv_signal(pvs,'pv', 'PV_loop', show, beat)
elif plot_type == 'cell':
a.plot_single_cell(epi, mid, endo, material, analysis, show)
