Compare commits
No commits in common. "8f27505ea4274676cdffae9586accf7e9d52cb18" and "c3c570764bd15cee2219b3173deda8351583526a" have entirely different histories.
8f27505ea4
...
c3c570764b
74
functions.py
74
functions.py
|
|
@ -15,10 +15,11 @@ from math import floor
|
|||
from statistics import mean
|
||||
from scipy.signal import find_peaks
|
||||
|
||||
def cycle(folder,dt,N_ts,save_path):
|
||||
pressure = np.loadtxt(folder + '/PHistRCR.dat',skiprows=1)
|
||||
def cycle(folder,dt,save_path):
|
||||
pressure = np.loadtxt(folder + '/PHistRCR.dat',skiprows=2)
|
||||
N_ts = pressure.shape[0]
|
||||
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
|
||||
peaks, _ = find_peaks(pressure[:,-1])
|
||||
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):
|
||||
pressure = np.loadtxt(folder+'/PHistRCR.dat',skiprows=1,)
|
||||
time = np.linspace(0,T_cyc*n_cyc,round(T_cyc/dt*n_cyc)+1)
|
||||
pressure = np.loadtxt(folder+'/PHistRCR.dat',skiprows=2,)
|
||||
time = np.linspace(0,T_cyc*n_cyc,round(T_cyc/dt*n_cyc))
|
||||
peak_P = []
|
||||
peak_P_pos = []
|
||||
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):
|
||||
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])]
|
||||
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
|
||||
|
||||
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)
|
||||
time = np.linspace(0,T_cyc,Nc)
|
||||
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):
|
||||
flow = np.loadtxt(folder+'/QHistRCR.dat',skiprows=1,)
|
||||
flow = np.loadtxt(folder+'/QHistRCR.dat',skiprows=2,)
|
||||
Nc = round(T_cyc/dt)
|
||||
time = np.linspace(0,T_cyc,Nc)
|
||||
fig, ax = plt.subplots()
|
||||
Q = np.empty(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))
|
||||
Q[i] = (mean(flow[N_ts-Nc+1:N_ts+1,i]))
|
||||
ax.plot(time,flow[N_ts-Nc:N_ts,i],label='ROI-'+str(i+2))
|
||||
Q[i] = (mean(flow[N_ts-Nc:N_ts,i]))
|
||||
|
||||
ax.set(xlabel='time [s]', ylabel='Flow [mL/s]',
|
||||
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):
|
||||
x = np.loadtxt(project_folder+'/ROI-1.flow')
|
||||
t = x[:,0]
|
||||
if (x[4,1] < 0):
|
||||
Q = -x[:,1]
|
||||
else:
|
||||
Q = x[:,1]
|
||||
Q = -x[:,1]
|
||||
Nt_pts = np.linspace(t_btw_rst,N_ts,int(N_ts/t_btw_rst))
|
||||
t_pts = Nt_pts*dt
|
||||
# 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
|
||||
f = open(project_folder+'/aimed.txt', 'r') #obrir el fitxer
|
||||
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 = []
|
||||
for i in range(0,len(CFD)):
|
||||
labels.append('ROI-'+str(i+2))
|
||||
|
|
@ -249,19 +243,19 @@ def barPlot(project_folder,DBP,MBP,SBP,PP,Q_avg,save_path):
|
|||
ax.legend()
|
||||
high = max(max(aimed,CFD))
|
||||
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."""
|
||||
for rect in rects:
|
||||
height = rect.get_height()
|
||||
ax.annotate(decimals.format(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,decimals)
|
||||
autolabel(rects2,decimals)
|
||||
autolabel(rects1)
|
||||
autolabel(rects2)
|
||||
# Create labels
|
||||
err = []
|
||||
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)
|
||||
|
||||
# DBP
|
||||
CFD = [float(i) for i in DBP]
|
||||
aimed = [float(i) for i in aimed_all[0]]
|
||||
plot_bar(CFD,aimed,'DBP',save_path,' [mmHg]','{0:.1f}')
|
||||
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 = [float(i) for i in MBP]
|
||||
aimed = [float(i) for i in aimed_all[1]]
|
||||
plot_bar(CFD,aimed,'MBP',save_path,' [mmHg]','{0:.1f}')
|
||||
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 = [float(i) for i in SBP]
|
||||
aimed = [float(i) for i in aimed_all[2]]
|
||||
plot_bar(CFD,aimed,'SBP',save_path,' [mmHg]','{0:.1f}')
|
||||
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 = [float(i) for i in PP]
|
||||
aimed = [float(i) for i in aimed_all[3]]
|
||||
plot_bar(CFD,aimed,'PP',save_path,' [mmHg]','{0:.1f}')
|
||||
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]')
|
||||
|
||||
# Q_avg
|
||||
CFD = [float(i) for i in Q_avg]
|
||||
aimed = [float(i) for i in aimed_all[4]]
|
||||
plot_bar(CFD,aimed,'Q_avg',save_path,' [mL/s]','{0:.2f}')
|
||||
# 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]')
|
||||
|
|
@ -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)
|
||||
Loading…
Reference in New Issue
Block a user