From a0040a70909bf3834ac21ca5b41afda347221ac2 Mon Sep 17 00:00:00 2001 From: Aloma Blanch Date: Thu, 23 Jul 2020 11:06:35 -0500 Subject: [PATCH] Added plot bars for DBP, MBP, SBP, Q --- functions.py | 71 +++++++++++++++++++++++++++++++++++++++++++++++--- generatePDF.py | 15 ++++++++--- main.py | 7 +++-- 3 files changed, 85 insertions(+), 8 deletions(-) diff --git a/functions.py b/functions.py index 5577723..e71da69 100644 --- a/functions.py +++ b/functions.py @@ -5,6 +5,7 @@ Created on Thu Jul 16 15:08:27 2020 @author: Aloma Blanch """ +import matplotlib import numpy as np import pandas as pd import matplotlib.pyplot as plt @@ -47,7 +48,7 @@ def cycle(folder,dt,save_path): ax.legend(loc=0) # plt.show() plt.savefig(save_path + '/Pressure_max_n_cycles.pdf') - return (floor(mean(peak_tdiff)*1000)/1000,len(t_pk_loc)) + return (floor(mean(peak_tdiff[1:])*1000)/1000,len(t_pk_loc)) def error_plot(folder,t_step,r_criteria,save_path): # Load iterations and residual error @@ -105,7 +106,7 @@ def periodicity(project,folder,dt,T_cyc,n_cyc,save_path): fig, ax = plt.subplots() ax.plot(time,pressure[:,-1]/1333.22,'b') - ax.plot(time[peak_P_pos], peak_P,'ro',label='Cylce pike') + ax.plot(time[peak_P_pos], peak_P,'ro',label='Cycle pike') ax.set(xlabel='time [s]', ylabel='Pressure [mmHg]', title='Pressure at last outlet') ax.spines['right'].set_visible(False) @@ -219,4 +220,68 @@ def rcr(project_folder): Rc_C_Rd[i][2] = float(rcr_lines[4+i*6][:-1]) return Rc_C_Rd - \ No newline at end of file +def barPlot(project_folder,DBP,MBP,SBP,PP,Q_avg,save_path): + import csv + with open(project_folder+'/aimed.csv', 'r') as f: + aimed_all = list(csv.reader(f, delimiter=';')) + + def plot_bar(CFD,aimed,name,save_path,unit): + labels = [] + for i in range(0,len(CFD)): + labels.append('ROI-'+str(i+2)) + x = np.arange(len(labels)) # the label locations + width = 0.38 # the width of the bars + fig, ax = plt.subplots() + rects1 = ax.bar(x - width/2, CFD, width, label='CFD',color='k') + rects2 = ax.bar(x + width/2, aimed, width, label='Aimed',color='white', edgecolor='black', hatch='/') + ax.set_ylabel(name + unit) + ax.set_title(name) + ax.set_xticks(x) + ax.set_xticklabels(labels) + ax.legend() + high = max(max(aimed,CFD)) + ax.set_ylim(top = 1.7*high) + def autolabel(rects): + """Attach a text label above each bar in *rects*, displaying its height.""" + for rect in rects: + height = rect.get_height() + ax.annotate('{}'.format(height), + xy=(rect.get_x() + rect.get_width() / 2, height), + xytext=(0,3), # 3 points vertical offset + textcoords="offset points", + ha='center', va='bottom', + fontsize = 7.5) + + autolabel(rects1) + autolabel(rects2) + fig.tight_layout() + # plt.show() + plt.savefig(save_path + '/'+name+'_bar.pdf') + fig.savefig(save_path + '/'+name+'_bar.jpg',dpi=150) + + # DBP + CFD = [round(float(i),1) for i in DBP] + aimed = [round(float(i),1) for i in aimed_all[0]] + plot_bar(CFD,aimed,'DBP',save_path,' [mmHg]') + + # MBP + CFD = [round(float(i),1) for i in MBP] + aimed = [round(float(i),1) for i in aimed_all[1]] + plot_bar(CFD,aimed,'MBP',save_path,' [mmHg]') + + # SBP + CFD = [round(float(i),1) for i in SBP] + aimed = [round(float(i),1) for i in aimed_all[2]] + plot_bar(CFD,aimed,'SBP',save_path,' [mmHg]') + + # PP + CFD = [round(float(i),1) for i in PP] + aimed = [round(float(i),1) for i in aimed_all[3]] + plot_bar(CFD,aimed,'PP',save_path,' [mmHg]') + + # DBP + CFD = [round(float(i),1) for i in Q_avg] + aimed = [round(float(i),1) for i in aimed_all[4]] + plot_bar(CFD,aimed,'Q_avg',save_path,' [mL/s]') + + \ No newline at end of file diff --git a/generatePDF.py b/generatePDF.py index 7bd0ba5..e47046a 100644 --- a/generatePDF.py +++ b/generatePDF.py @@ -135,14 +135,23 @@ def generatePDF(save_path,project,DBP,MBP,SBP,PP,Q_avg,txt1,txt2,Rc_C_Rd,T_cyc,n pdf.cell(width_cell[0],10,'ROI-'+ str(i+2),1,0,'C') pdf.cell(width_cell[1],10,str(round(Q_avg[i],3)),1,1,'C') + pdf.print_section(4, 'Checking CFD values vs aimed') + pdf.image(save_path +'/DBP_bar.jpg', x = 40, y = 50, w = 140, h = 100, type = '', link = '') + pdf.image(save_path +'/MBP_bar.jpg', x = 40, y = 170, w = 140, h = 100, type = '', link = '') + pdf.print_section(4, 'Checking CFD values vs aimed') + pdf.image(save_path +'/SBP_bar.jpg', x = 40, y = 50, w = 140, h = 100, type = '', link = '') + pdf.image(save_path +'/PP_bar.jpg', x = 40, y = 170, w = 140, h = 100, type = '', link = '') + pdf.print_section(4, 'Checking CFD values vs aimed') + pdf.image(save_path +'/Q_avg_bar.jpg', x = 40, y = 50, w = 140, h = 100, type = '', link = '') + pdf.cell(0, 10,'', 0, 1) - pdf.print_section(4, 'Checking Inlet Flow Waveform and time steps saved') + pdf.print_section(5, 'Checking Inlet Flow Waveform and time steps saved') pdf.image(save_path +'/inlet_waveform.jpg', x = 40, y = 50, w = 140, h = 100, type = '', link = '') pdf.set_font('Arial', '', 10) pdf.cell(0, 5,txt2, 0, 1) width_cell=[30,20,110]; - pdf.print_section(5, 'Outlet Boundary Condition Applied') + pdf.print_section(6, 'Outlet Boundary Condition Applied') pdf.set_font('Arial', '', 10) pdf.cell(0, 10,'3-Element Windkessel model outlet boundary condition applied in SimVascular', 0, 1) pdf.cell(0, 5,' ', 0, 1) @@ -162,7 +171,7 @@ def generatePDF(save_path,project,DBP,MBP,SBP,PP,Q_avg,txt1,txt2,Rc_C_Rd,T_cyc,n pdf.cell(width_cell[0],10,str(Rc_C_Rd[i][1]),1,0,'C') # Third column of row 1 pdf.cell(width_cell[0],10,str(Rc_C_Rd[i][2]),1,1,'C') # Third column of row 1 - pdf.print_section(6, 'Solver Parameters') + pdf.print_section(7, 'Solver Parameters') pdf.set_font('Arial', '', 10) pdf.cell(0, 10,'Cardiac cycle period - ' + str(T_cyc) + 's', 0, 1) pdf.cell(0, 10,'Number of cycles - ' + str(n_cyc), 0, 1) diff --git a/main.py b/main.py index 92bfef2..a535abe 100644 --- a/main.py +++ b/main.py @@ -20,7 +20,7 @@ import os.path # from fpdf import FPDF -from functions import error_plot, periodicity, pressure, flow, inlet_flow_waveform, rcr, cycle +from functions import error_plot, periodicity, pressure, flow, inlet_flow_waveform, rcr, cycle, barPlot from generatePDF import generatePDF @@ -53,7 +53,7 @@ Rc_C_Rd = rcr(project_folder) # Extracting number of cycles and period of cardiac cycle (T_cyc,n_cyc) = cycle(folder,dt,save_path) - +T_cyc=0.477 # Cehcking convergency and periodicity error_plot(folder,dt,rc,save_path) txt1 = periodicity(project,folder,dt,T_cyc,n_cyc,save_path) @@ -67,5 +67,8 @@ txt1 = periodicity(project,folder,dt,T_cyc,n_cyc,save_path) # Inlet Flow Rate + and t saved txt2 = inlet_flow_waveform(project_folder,t_btw_rst,N_ts,dt,T_cyc,n_cyc,save_path) +# Extract bar plots CFD vs. aimed +barPlot(project_folder,DBP,MBP,SBP,PP,Q_avg,save_path) + # Create PDF report generatePDF(save_path,project,DBP,MBP,SBP,PP,Q_avg,txt1,txt2,Rc_C_Rd,T_cyc,n_cyc,N_ts,dt)