aloma/SimVascular_report_PostProcess
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Compare commits

base: aloma:8f27505ea4274676cdffae9586accf7e9d52cb18
Branches Tags
aloma:master
..
compare: aloma:c3c570764bd15cee2219b3173deda8351583526a
Branches Tags
aloma:master

No commits in common. "8f27505ea4274676cdffae9586accf7e9d52cb18" and "c3c570764bd15cee2219b3173deda8351583526a" have entirely different histories.
8f27505ea4 ... c3c570764b

2 changed files with 34 additions and 139 deletions
Show Stats Expand all files Collapse all files

72
functions.py
View File

@ -15,10 +15,11 @@ from math import floor
from statistics import mean from statistics import mean
from scipy.signal import find_peaks from scipy.signal import find_peaks
def cycle(folder,dt,N_ts,save_path): def cycle(folder,dt,save_path):
pressure = np.loadtxt(folder + '/PHistRCR.dat',skiprows=1) pressure = np.loadtxt(folder + '/PHistRCR.dat',skiprows=2)
N_ts = pressure.shape[0]
T = round(dt*N_ts,3) T = round(dt*N_ts,3)
time = np.linspace(0,T,N_ts+1) time = np.linspace(0,T,N_ts)
# Find peaks, keep only the maximums # Find peaks, keep only the maximums
peaks, _ = find_peaks(pressure[:,-1]) peaks, _ = find_peaks(pressure[:,-1])
idx = np.where(pressure[peaks,-1] >= np.mean(pressure)) idx = np.where(pressure[peaks,-1] >= np.mean(pressure))
@ -91,8 +92,8 @@ def error_plot(folder,t_step,r_criteria,save_path):
def periodicity(project,folder,dt,T_cyc,n_cyc,save_path): def periodicity(project,folder,dt,T_cyc,n_cyc,save_path):
pressure = np.loadtxt(folder+'/PHistRCR.dat',skiprows=1,) pressure = np.loadtxt(folder+'/PHistRCR.dat',skiprows=2,)
time = np.linspace(0,T_cyc*n_cyc,round(T_cyc/dt*n_cyc)+1) time = np.linspace(0,T_cyc*n_cyc,round(T_cyc/dt*n_cyc))
peak_P = [] peak_P = []
peak_P_pos = [] peak_P_pos = []
Nc = round(T_cyc/dt) Nc = round(T_cyc/dt)
@ -118,13 +119,10 @@ def periodicity(project,folder,dt,T_cyc,n_cyc,save_path):
if (peak_Pdiff[-1]<=1): if (peak_Pdiff[-1]<=1):
print('The numerical simulation \'{0}\' has achieve periodicity!\nSystolic Blood Pressure (SBP):\nsecond-last cycle = {1:.2f} mmHg,\nlast cycle = {2:.2f} mmHg,\n\u0394mmHg = {3:.2f} mmHg'.format(project,peak_P[-2],peak_P[-1],peak_Pdiff[-1])) print('The numerical simulation \'{0}\' has achieve periodicity!\nSystolic Blood Pressure (SBP):\nsecond-last cycle = {1:.2f} mmHg,\nlast cycle = {2:.2f} mmHg,\n\u0394mmHg = {3:.2f} mmHg'.format(project,peak_P[-2],peak_P[-1],peak_Pdiff[-1]))
txt = ['The numerical simulation \'{0}\' has achieved periodicity.'.format(project), 'Systolic Blood Pressure (SBP):', 'Second-last cycle = {0:.2f} mmHg'.format(peak_P[-2]), 'Last cycle = {0:.2f} mmHg'.format(peak_P[-1]), 'Delta_mmHg = {0:.2f} mmHg'.format(peak_Pdiff[-1])] txt = ['The numerical simulation \'{0}\' has achieved periodicity.'.format(project), 'Systolic Blood Pressure (SBP):', 'Second-last cycle = {0:.2f} mmHg'.format(peak_P[-2]), 'Last cycle = {0:.2f} mmHg'.format(peak_P[-1]), 'Delta_mmHg = {0:.2f} mmHg'.format(peak_Pdiff[-1])]
else:
print('The numerical simulation \'{0}\' has not achieve periodicity...\nSystolic Blood Pressure (SBP):\nsecond-last cycle = {1:.2f} mmHg,\nlast cycle = {2:.2f} mmHg,\n\u0394mmHg = {3:.2f} mmHg'.format(project,peak_P[-2],peak_P[-1],peak_Pdiff[-1]))
txt = ['The numerical simulation \'{0}\' has achieved periodicity.'.format(project), 'Systolic Blood Pressure (SBP):', 'Second-last cycle = {0:.2f} mmHg'.format(peak_P[-2]), 'Last cycle = {0:.2f} mmHg'.format(peak_P[-1]), 'Delta_mmHg = {0:.2f} mmHg'.format(peak_Pdiff[-1])]
return txt return txt
def pressure(folder,N_ts,T_cyc,dt,n_cyc,save_path): def pressure(folder,N_ts,T_cyc,dt,n_cyc,save_path):
pressure = np.loadtxt(folder+'/PHistRCR.dat',skiprows=1,) pressure = np.loadtxt(folder+'/PHistRCR.dat',skiprows=2,)
Nc = round(T_cyc/dt) Nc = round(T_cyc/dt)
time = np.linspace(0,T_cyc,Nc) time = np.linspace(0,T_cyc,Nc)
fig, ax = plt.subplots() fig, ax = plt.subplots()
@ -149,14 +147,14 @@ def pressure(folder,N_ts,T_cyc,dt,n_cyc,save_path):
def flow(folder,N_ts,T_cyc,dt,n_cyc,save_path): def flow(folder,N_ts,T_cyc,dt,n_cyc,save_path):
flow = np.loadtxt(folder+'/QHistRCR.dat',skiprows=1,) flow = np.loadtxt(folder+'/QHistRCR.dat',skiprows=2,)
Nc = round(T_cyc/dt) Nc = round(T_cyc/dt)
time = np.linspace(0,T_cyc,Nc) time = np.linspace(0,T_cyc,Nc)
fig, ax = plt.subplots() fig, ax = plt.subplots()
Q = np.empty(flow.shape[1]) Q = np.empty(flow.shape[1])
for i in range(0,flow.shape[1]): for i in range(0,flow.shape[1]):
ax.plot(time,flow[N_ts-Nc+1:N_ts+1,i],label='ROI-'+str(i+2)) ax.plot(time,flow[N_ts-Nc:N_ts,i],label='ROI-'+str(i+2))
Q[i] = (mean(flow[N_ts-Nc+1:N_ts+1,i])) Q[i] = (mean(flow[N_ts-Nc:N_ts,i]))
ax.set(xlabel='time [s]', ylabel='Flow [mL/s]', ax.set(xlabel='time [s]', ylabel='Flow [mL/s]',
title='Flow at each outlet') title='Flow at each outlet')
@ -171,10 +169,7 @@ def flow(folder,N_ts,T_cyc,dt,n_cyc,save_path):
def inlet_flow_waveform(project_folder,t_btw_rst,N_ts,dt,T_cyc,n_cyc,save_path): def inlet_flow_waveform(project_folder,t_btw_rst,N_ts,dt,T_cyc,n_cyc,save_path):
x = np.loadtxt(project_folder+'/ROI-1.flow') x = np.loadtxt(project_folder+'/ROI-1.flow')
t = x[:,0] t = x[:,0]
if (x[4,1] < 0): Q = -x[:,1]
Q = -x[:,1]
else:
Q = x[:,1]
Nt_pts = np.linspace(t_btw_rst,N_ts,int(N_ts/t_btw_rst)) Nt_pts = np.linspace(t_btw_rst,N_ts,int(N_ts/t_btw_rst))
t_pts = Nt_pts*dt t_pts = Nt_pts*dt
# Put all the time values on a single cardiac cylce # Put all the time values on a single cardiac cylce
@ -229,11 +224,10 @@ def barPlot(project_folder,DBP,MBP,SBP,PP,Q_avg,save_path):
aimed_all = [1] #No li agrada fer append en una llista buida, no recordo perquè i segur que es pot fer millor aimed_all = [1] #No li agrada fer append en una llista buida, no recordo perquè i segur que es pot fer millor
f = open(project_folder+'/aimed.txt', 'r') #obrir el fitxer f = open(project_folder+'/aimed.txt', 'r') #obrir el fitxer
for line in f: for line in f:
aimed_all.append(line.split()) #fas una llista amb cada valor de la línia carregan-te els espais buits aimed_all = line.split() #fas una llista amb cada valor de la línia carregan-te els espais buits
del aimed_all[0] #esborrar el primer valor que has posat perquè no estigues buida la llista del aimed_all[0] #esborrar el primer valor que has posat perquè no estigues buida la llista
def plot_bar(CFD,aimed,name,save_path,unit,decimals): def plot_bar(CFD,aimed,name,save_path,unit):
labels = [] labels = []
for i in range(0,len(CFD)): for i in range(0,len(CFD)):
labels.append('ROI-'+str(i+2)) labels.append('ROI-'+str(i+2))
@ -249,19 +243,19 @@ def barPlot(project_folder,DBP,MBP,SBP,PP,Q_avg,save_path):
ax.legend() ax.legend()
high = max(max(aimed,CFD)) high = max(max(aimed,CFD))
ax.set_ylim(top = 1.7*high) ax.set_ylim(top = 1.7*high)
def autolabel(rects,decimals): def autolabel(rects):
"""Attach a text label above each bar in *rects*, displaying its height.""" """Attach a text label above each bar in *rects*, displaying its height."""
for rect in rects: for rect in rects:
height = rect.get_height() height = rect.get_height()
ax.annotate(decimals.format(height), ax.annotate('{}'.format(height),
xy=(rect.get_x() + rect.get_width() / 2, height), xy=(rect.get_x() + rect.get_width() / 2, height),
xytext=(0,3), # 3 points vertical offset xytext=(0,3), # 3 points vertical offset
textcoords="offset points", textcoords="offset points",
ha='center', va='bottom', ha='center', va='bottom',
fontsize = 7.5) fontsize = 7.5)
autolabel(rects1,decimals) autolabel(rects1)
autolabel(rects2,decimals) autolabel(rects2)
# Create labels # Create labels
err = [] err = []
zip_object = zip(CFD, aimed) zip_object = zip(CFD, aimed)
@ -278,26 +272,26 @@ def barPlot(project_folder,DBP,MBP,SBP,PP,Q_avg,save_path):
fig.savefig(save_path + '/'+name+'_bar.jpg',dpi=150) fig.savefig(save_path + '/'+name+'_bar.jpg',dpi=150)
# DBP # DBP
CFD = [float(i) for i in DBP] CFD = [round(float(i),1) for i in DBP]
aimed = [float(i) for i in aimed_all[0]] aimed = [round(float(i),1) for i in aimed_all[0]]
plot_bar(CFD,aimed,'DBP',save_path,' [mmHg]','{0:.1f}') plot_bar(CFD,aimed,'DBP',save_path,' [mmHg]')
# MBP # MBP
CFD = [float(i) for i in MBP] CFD = [round(float(i),1) for i in MBP]
aimed = [float(i) for i in aimed_all[1]] aimed = [round(float(i),1) for i in aimed_all[1]]
plot_bar(CFD,aimed,'MBP',save_path,' [mmHg]','{0:.1f}') plot_bar(CFD,aimed,'MBP',save_path,' [mmHg]')
# SBP # SBP
CFD = [float(i) for i in SBP] CFD = [round(float(i),1) for i in SBP]
aimed = [float(i) for i in aimed_all[2]] aimed = [round(float(i),1) for i in aimed_all[2]]
plot_bar(CFD,aimed,'SBP',save_path,' [mmHg]','{0:.1f}') plot_bar(CFD,aimed,'SBP',save_path,' [mmHg]')
# PP # PP
CFD = [float(i) for i in PP] CFD = [round(float(i),1) for i in PP]
aimed = [float(i) for i in aimed_all[3]] aimed = [round(float(i),1) for i in aimed_all[3]]
plot_bar(CFD,aimed,'PP',save_path,' [mmHg]','{0:.1f}') plot_bar(CFD,aimed,'PP',save_path,' [mmHg]')
# Q_avg # DBP
CFD = [float(i) for i in Q_avg] CFD = [round(float(i),1) for i in Q_avg]
aimed = [float(i) for i in aimed_all[4]] aimed = [round(float(i),1) for i in aimed_all[4]]
plot_bar(CFD,aimed,'Q_avg',save_path,' [mL/s]','{0:.2f}') plot_bar(CFD,aimed,'Q_avg',save_path,' [mL/s]')

99
report_main.py
View File

@ -1,99 +0,0 @@
#!/usr/bin/env python
"""
Created on Thu Jul 16 15:08:27 2020
@author: Aloma Blanch
"""
# import numpy as np
# import matplotlib.pyplot as plt
# import matplotlib.ticker as mtick
from tkinter import Tk
# from tkinter.filedialog import askopenfilename, asksaveasfile, askdirectory
from tkinter.filedialog import askdirectory
# import pandas as pd
# import tkinter as tk
import os.path
# from scipy.signal import find_peaks_cwt
# from scipy import signal
# import statistics
# from fpdf import FPDF
from functions import error_plot, periodicity, pressure, flow, inlet_flow_waveform, rcr, cycle, barPlot
from generatePDF import generatePDF
# Selct dir
Tk().withdraw()
folder = askdirectory()
project_folder = os.path.dirname(folder)
project = os.path.basename(project_folder)
save_path = folder+'/'+project+'-report'
os.mkdir(save_path)
save_pdf = save_path + '/' + project + '-report.pdf'
# Read important parameters from - solver.inp file
mylines = [] # Declare an empty list named mylines.
with open (project_folder + '/solver.inp', 'rt') as myfile: # Open lorem.txt for reading text data.
for myline in myfile: # For each line, stored as myline,
mylines.append(myline) # add its contents to mylines.
# Number of Timesteps - idx 3
# Idea: remove the text to extract the number, the text part will be the same no matter the simulation
N_ts = int(mylines[3][20:-1])
# Time Step Size - idx 4
dt = float(mylines[4][16:-1])
# Residual criteria - idx 4
rc = float(mylines[26][18:-1])
# Imesteps between Restarts - idx 6
t_btw_rst = int(mylines[6][37:-1])
# Extracting Outlet Boundary Conditions from - rcrt.dat fle
Rc_C_Rd = rcr(project_folder)
# Extracting number of cycles and period of cardiac cycle
(T_cyc,n_cyc) = cycle(folder,dt,N_ts,save_path)
T_cyc=0.4769
n_cyc=4
# Cehcking convergency and periodicity
error_plot(folder,dt,rc,save_path)
txt1 = periodicity(project,folder,dt,T_cyc,n_cyc,save_path)
# Pressure - Outlets
(DBP,MBP,SBP,PP) = pressure(folder,N_ts,T_cyc,dt,n_cyc,save_path)
# Flow Rate - Outlets
(Q_avg) = flow(folder,N_ts,T_cyc,dt,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)
# # Flow comparison specific for patient 120.
# # Only ROI-5,6,8 because that is the PC-MRA data that we have.
# def flow_comparison(folder,N_ts,T_cyc,dt,n_cyc,save_path):
# flow_sim = np.loadtxt(folder+'/QHistRCR.dat',skiprows=2,)
# flow_doc = np.loadtxt(os.path.dirname(os.path.dirname(os.path.dirname(os.path.dirname(os.path.dirname(folder)))))+'/Data/flow_waveforms.txt',skiprows=2,)
# Nc = round(T_cyc/dt)
# time_sim = np.linspace(0,T_cyc,Nc)
# time_doc = np.linspace(0,T_cyc,flow_doc.shape[0])
# Q = np.empty(flow_sim.shape[1])
# for i in range(0,flow_sim.shape[1]):
# fig, ax = plt.subplots()
# ax.plot(time_sim,flow_sim[N_ts-Nc:N_ts,i],label='sim ROI-'+str(i+2))
# ax.plot(time_doc,flow_doc[:,1+i],label='doc ROI-'+str(i+2))
# ax.set(xlabel='time [s]', ylabel='Flow [mL/s]',
# title='Flow at each outlet')
# ax.spines['right'].set_visible(False)
# ax.spines['top'].set_visible(False)
# ax.legend(loc=0)
# # plt.show()
# plt.savefig(save_path + '/flow_comparison.pdf')
# fig.savefig(save_path + '/flow_comparison.jpg',dpi=150)
# 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)