diff --git a/codes/.vscode/launch.json b/codes/.vscode/launch.json deleted file mode 100644 index 17e15f2..0000000 --- a/codes/.vscode/launch.json +++ /dev/null @@ -1,15 +0,0 @@ -{ - // Use IntelliSense to learn about possible attributes. - // Hover to view descriptions of existing attributes. - // For more information, visit: https://go.microsoft.com/fwlink/?linkid=830387 - "version": "0.2.0", - "configurations": [ - { - "name": "Python: Current File", - "type": "python", - "request": "launch", - "program": "${file}", - "console": "integratedTerminal" - } - ] -} \ No newline at end of file diff --git a/codes/.vscode/tasks.json b/codes/.vscode/tasks.json deleted file mode 100644 index 6c28385..0000000 --- a/codes/.vscode/tasks.json +++ /dev/null @@ -1,12 +0,0 @@ -{ - // See https://go.microsoft.com/fwlink/?LinkId=733558 - // for the documentation about the tasks.json format - "version": "2.0.0", - "tasks": [ - { - "label": "echo", - "type": "shell", - "command": "echo Hello" - } - ] -} \ No newline at end of file diff --git a/codes/CS.py b/codes/CS.py deleted file mode 100644 index 85a1a1c..0000000 --- a/codes/CS.py +++ /dev/null @@ -1,597 +0,0 @@ -import numpy as np -from numpy import linalg as LA -import sys -from mpi4py import MPI -comm = MPI.COMM_WORLD -size = comm.Get_size() -rank = comm.Get_rank() - - -# COMPRESSED SENSING: LINEAR BREGMAN METHOD -# Translated and adapted into python from tinycs -# -# *tinycs* is a minimal compressed sensing (CS) toolkit designed -# to allow MR imaging scientists to design undersampled -# acquisitions and reconstruct the resulting data with CS without -# needing to be a CS expert. -# -# The Cartesian reconstruction is based on the split Bregman -# code written by Tom Goldstein, originally available here: -# - - -def pdf(k, kw, klo, q): - - p = (np.abs(k)/kw)**(-q) - p[np.where(k == 0)] = 0 - p[np.where(np.abs(k) <= kw)] = 1 - p[np.where(k < klo)] = 0 - - return p - -def mask_pdf_1d(n, norm, q, pf): - - ks = np.arange(0, n) - np.ceil(n/2) - 1 - kmax = np.floor(n/2) - npf = np.round(pf*n) - klo = ks[n-npf] - - for k in range(int(kmax)): - P = pdf(ks, k+1, klo, q) - if np.sum(P) >= norm: - break - - P = np.fft.fftshift(P) - - return P - -def mask_pdf_2d(dims, norm, q, pf): - - nz = dims[1] - ny = dims[0] - yc = round(ny/2) - zc = round(nz/2) - rmax = np.sqrt((ny-yc)**2 + (nz-zc)**2) - [Z, Y] = np.meshgrid(np.arange(0, nz), np.arange(0, ny)) - RR = np.sqrt((Y-yc)**2 + (Z-zc)**2) - Z = np.abs(Z - nz/2 - 0.5) - Y = np.abs(Y - ny/2 - 0.5) - - for rw in range(1, int(rmax)+1): - P = np.ones([ny, nz])/pf - C = np.logical_and(Z <= rw, Y <= rw) - W = np.logical_or(Z > rw, Y > rw) - P[W] = (RR[W]/rw)**(-q) - if np.sum(P) >= norm: - break - - return [P, C] - -def GeneratePattern(dim, R): - - # 3D CASE - if np.size(dim) == 3: - - nro = dim[0] - npe = dim[1] - nacq = round(npe/R) - q = 1 - pf = 1 - P = mask_pdf_1d(npe, nacq, q, pf) - while True: - M = np.random.rand(npe) - M = 1*(M <= P) - if np.sum(M) == nacq: - break - # remove partial Fourier plane and compensate sampling density - M = M != 0 - M = np.tile(M, [nro, 1]) - #M = M.T - - # 4D CASE - if np.size(dim) == 4: - - nro = dim[0] - npe1 = dim[1] - npe2 = dim[2] - nacq = round(npe1*npe2/R) - q = 1 - pf = 1 - [P, C] = mask_pdf_2d([npe1, npe2], nacq, q, pf) - RR = np.random.rand(npe1, npe2) - M = (RR <= P) - nchosen = np.sum(M) - if nchosen > nacq: # Correct for inexact number chosen - #outerOn = np.logical_and( M , P!=1 ) - outerOn = np.where((M)*(P != 1)) - numToFlip = nchosen-nacq - idxs = np.random.permutation(outerOn[0].size) - idxx = outerOn[0][idxs[0:numToFlip]] - idxy = outerOn[1][idxs[0:numToFlip]] - M[idxx, idxy] = False - - elif nchosen < nacq: - outerOff = np.where(~M) - idxs = np.random.permutation(outerOff[0].size) - numToFlip = nacq - nchosen - idxx = outerOff[0][idxs[0:numToFlip]] - idxy = outerOff[1][idxs[0:numToFlip]] - M[idxx, idxy] = True - - M = np.rollaxis(np.tile(np.rollaxis(M, 1), [nro, 1, 1]), 2) - M = np.fft.ifftshift(M) - M = M.transpose((1, 0, 2)) - - return M - -def get_norm_factor(MASK, uu): - UM = MASK == 1 - return UM.shape[0]/LA.norm(uu) - -def Dxyzt(X): - - if np.ndim(X) == 3: - dd0 = X[:, :, 0] - dd1 = X[:, :, 1] - DA = dd0 - np.vstack((dd0[1::, :], dd0[0, :])) - DB = dd1 - np.hstack((dd1[:, 1::], dd1[:, 0:1])) - - return DA + DB - - if np.ndim(X) == 4: - dd0 = X[:, :, :, 0] - dd1 = X[:, :, :, 1] - dd2 = X[:, :, :, 2] - - DA = dd0 - np.vstack((dd0[1::, :, :], dd0[0, :, :][np.newaxis, :, :])) - DB = dd1 - np.hstack((dd1[:, 1::, :], dd1[:, 0, :][:, np.newaxis, :])) - DC = dd2 - np.dstack((dd2[:, :, 1::], dd2[:, :, 0][:, :, np.newaxis])) - - return DA + DB + DC - -def Dxyz(u): - - if np.ndim(u) == 2: - - dx = u[:, :] - np.vstack((u[-1, :], u[0:-1, :])) - dy = u[:, :] - np.hstack((u[:, -1:], u[:, 0:-1])) - D = np.zeros([dx.shape[0], dx.shape[1], 2], dtype=complex) - D[:, :, 0] = dx - D[:, :, 1] = dy - - return D - - if np.ndim(u) == 3: - - dx = u[:, :, :] - \ - np.vstack((u[-1, :, :][np.newaxis, :, :], u[0:-1, :, :])) - dy = u[:, :, :] - \ - np.hstack((u[:, -1, :][:, np.newaxis, :], u[:, 0:-1, :])) - dz = u[:, :, :] - \ - np.dstack((u[:, :, -1][:, :, np.newaxis], u[:, :, 0:-1])) - - D = np.zeros([dx.shape[0], dx.shape[1], dx.shape[2], 3], dtype=complex) - - D[:, :, :, 0] = dx - D[:, :, :, 1] = dy - D[:, :, :, 2] = dz - - return D - -def shrink(X, pgam): - p = 1 - s = np.abs(X) - tt = pgam/(s)**(1-p) - # t = pgam/np.sqrt(s) - ss = s-tt - ss = ss*(ss > 0) - s = s + 1*(s < tt) - ss = ss/s - return ss*X - -def CSMETHOD(ITOT, R): - ''' Compressed Sensing Function. - - Args: - ITOT: a numpy matrix with the full sampled (3D or 4D) dynamical data - R: the acceleration factor - ''' - # Method parameters - ninner = 5 - nbreg = 10 - lmbda = 4 - mu = 20 - gam = 1 - - if np.ndim(ITOT) == 3: - [row, col, numt2] = ITOT.shape - elif np.ndim(ITOT) == 4: - [row, col, dep, numt2] = ITOT.shape - else: - raise Exception('Dynamical data is requested') - - MASK = GeneratePattern(ITOT.shape, R) - CS1 = np.zeros(ITOT.shape, dtype=complex) - - nit = 0 - nit_tot = (numt2-1)/20 - - if np.ndim(ITOT) == 3: - - for t in range(numt2): - - if rank == 0: - print('{3D COMPRESSED SENSING} t = ', t) - - Kdata = np.fft.fft2(ITOT[:, :, t])*MASK - - data_ndims = Kdata.ndim - mask = Kdata != 0 # not perfect, but good enough - # normalize the data so that standard parameter values work - norm_factor = get_norm_factor(mask, Kdata) - Kdata = Kdata*norm_factor - # Reserve memory for the auxillary variables - Kdata0 = Kdata - img = np.zeros([row, col], dtype=complex) - X = np.zeros([row, col, data_ndims]) - B = np.zeros([row, col, data_ndims]) - # Build Kernels - scale = np.sqrt(row*col) - murf = np.fft.ifft2(mu*mask*Kdata)*scale - uker = np.zeros([row, col]) - uker[0, 0] = 4 - uker[0, 1] = -1 - uker[1, 0] = -1 - uker[-1, 0] = -1 - uker[0, -1] = -1 - - uker = 1/(mu*mask + lmbda*np.fft.fftn(uker) + gam) - - # Do the reconstruction - for outer in range(nbreg): - for inner in range(ninner): - # update u - rhs = murf + lmbda*Dxyzt(X-B) + gam*img - img = np.fft.ifft2(np.fft.fft2(rhs)*uker) - # update x and y - A = Dxyz(img) + B - X = shrink(A, 1/lmbda) - # update bregman parameters - B = A - X - Kdata = Kdata + Kdata0 - mask*np.fft.fftn(img)/scale - murf = np.fft.ifftn(mu*mask*Kdata)*scale - - # undo the normalization so that results are scaled properly - img = img / norm_factor / scale - - CS1[:, :, t] = img - - if np.ndim(ITOT) == 4: - - for t in range(numt2): - - if rank == 0: - print( - '[4D CS] R = {re} t = {te}/{tef}'.format(re=R, te=t, tef=numt2)) - - Kdata_0 = np.fft.fftn(ITOT[:, :, :, t]) - Kdata = Kdata_0*MASK - data_ndims = Kdata.ndim - mask = Kdata != 0 # not perfect, but good enough - # normalize the data so that standard parameter values work - norm_factor = get_norm_factor(mask, Kdata) - Kdata = Kdata*norm_factor - # Reserve memory for the auxillary variables - Kdata0 = Kdata - img = np.zeros([row, col, dep], dtype=complex) - X = np.zeros([row, col, dep, data_ndims]) - B = np.zeros([row, col, dep, data_ndims]) - # Build Kernels - scale = np.sqrt(row*col*dep) - murf = np.fft.ifftn(mu*mask*Kdata)*scale - uker = np.zeros([row, col, dep]) - uker[0, 0, 0] = 8 - uker[1, 0, 0] = -1 - uker[0, 1, 0] = -1 - uker[0, 0, 1] = -1 - uker[-1, 0, 0] = -1 - uker[0, -1, 0] = -1 - uker[0, 0, -1] = -1 - uker = 1/(mu*mask + lmbda*np.fft.fftn(uker) + gam) - - # Do the reconstruction - for outer in range(nbreg): - for inner in range(ninner): - # update u - rhs = murf + lmbda*Dxyzt(X-B) + gam*img - img = np.fft.ifft2(np.fft.fft2(rhs)*uker) - # update x and y - A = Dxyz(img) + B - X = shrink(A, 1/lmbda) - # update bregman parameters - B = A - X - Kdata = Kdata + Kdata0 - mask*np.fft.fftn(img)/scale - murf = np.fft.ifftn(mu*mask*Kdata)*scale - - # undo the normalization so that results are scaled properly - img = img / norm_factor / scale - - CS1[:, :, :, t] = img - - return CS1 - -def CSMETHOD_SENSE(ITOT, R, R_SENSE): - ''' Compressed sense algorith with SENSE... in contruction!. - - Args: - ITOT: a numpy matrix with the full sampled (3D or 4D) dynamical data - R: the acceleration factor - ''' - # Method parameters - ninner = 5 - nbreg = 10 - lmbda = 4 - mu = 20 - gam = 1 - - [row, col, dep, numt2] = ITOT.shape - MASK = {} - ITOTCS = {} - MASK[0] = GeneratePattern([row, int(np.ceil(col/2)), dep, numt2], R) - MASK[1] = GeneratePattern([row, int(np.ceil(col/2)), dep, numt2], R) - SenseMAP = {} - [SenseMAP[0], SenseMAP[1]] = Sensitivity_Map([row, col, dep]) - - col = int(np.ceil(col/2)) - ITOTCS[0] = np.zeros([row, col, dep, numt2], dtype=complex) - ITOTCS[1] = np.zeros([row, col, dep, numt2], dtype=complex) - - for rs in range(R_SENSE): - for t in range(numt2): - if rank == 0: - print( - '[4D CS] R = {re} t = {te}/{tef}'.format(re=R, te=t, tef=numt2)) - - Kdata_0 = np.fft.fftn(ITOT[:, :, :, t]) - Kdata_0 = Kdata_0*SenseMAP[rs] - Kdata_0 = Kdata_0[:, 0::R_SENSE, :] - - Kdata = Kdata_0*MASK[rs] - data_ndims = Kdata.ndim - mask = Kdata != 0 # not perfect, but good enough - # normalize the data so that standard parameter values work - norm_factor = get_norm_factor(mask, Kdata) - Kdata = Kdata*norm_factor - # Reserve memory for the auxillary variables - Kdata0 = Kdata - img = np.zeros([row, col, dep], dtype=complex) - X = np.zeros([row, col, dep, data_ndims]) - B = np.zeros([row, col, dep, data_ndims]) - # Build Kernels - scale = np.sqrt(row*col*dep) - murf = np.fft.ifftn(mu*mask*Kdata)*scale - uker = np.zeros([row, col, dep]) - uker[0, 0, 0] = 8 - uker[1, 0, 0] = -1 - uker[0, 1, 0] = -1 - uker[0, 0, 1] = -1 - uker[-1, 0, 0] = -1 - uker[0, -1, 0] = -1 - uker[0, 0, -1] = -1 - uker = 1/(mu*mask + lmbda*np.fft.fftn(uker) + gam) - - # Do the reconstruction - for outer in range(nbreg): - for inner in range(ninner): - # update u - rhs = murf + lmbda*Dxyzt(X-B) + gam*img - img = np.fft.ifft2(np.fft.fft2(rhs)*uker) - # update x and y - A = Dxyz(img) + B - X = shrink(A, 1/lmbda) - # update bregman parameters - B = A - X - Kdata = Kdata + Kdata0 - mask*np.fft.fftn(img)/scale - murf = np.fft.ifftn(mu*mask*Kdata)*scale - - # undo the normalization so that results are scaled properly - img = img / norm_factor / scale - - ITOTCS[rs][:, :, :, t] = img - - return [ITOTCS[0], ITOTCS[1]] - -def phase_contrast(M1, M0, VENC, scantype='0G'): - param = 1 - if scantype == '-G+G': - param = 0.5 - return VENC*param*(np.angle(M1) - np.angle(M0))/np.pi - -def GenerateMagnetization(Sq, VENC, noise, scantype='0G'): - ''' Simulation of a typical magnetization. A x-dependent plane is added into the - reference phase. - ''' - # MRI PARAMETERS - gamma = 267.513e6 # rad/Tesla/sec Gyromagnetic ratio for H nuclei - B0 = 1.5 # Tesla Magnetic Field Strenght - TE = 5e-3 # Echo-time - PHASE0 = np.zeros(Sq.shape) - PHASE1 = np.zeros(Sq.shape) - RHO0 = np.zeros(Sq.shape, dtype=complex) - RHO1 = np.zeros(Sq.shape, dtype=complex) - - if np.ndim(Sq) == 3: - [row, col, numt2] = Sq.shape - [X, Y] = np.meshgrid(np.linspace(0, col, col), - np.linspace(0, row, row)) - for k in range(numt2): - if noise: - Drho = np.random.normal(0, 0.2, [row, col]) - Drho2 = np.random.normal(0, 0.2, [row, col]) - else: - Drho = np.zeros([row, col]) - Drho2 = np.zeros([row, col]) - - varPHASE0 = np.random.randint(-10, 11, size=(row, col))*np.pi/180*( - np.abs(Sq[:, :, k]) < 0.001) # Hugo's observation - modulus = 0.5 + 0.5*(np.abs(Sq[:, :, k]) > 0.001) - - if scantype == '0G': - PHASE0[:, :, k] = (gamma*B0*TE+0.01*X) * \ - (np.abs(Sq[:, :, k]) > 0.001) + 10*varPHASE0 - PHASE1[:, :, k] = (gamma*B0*TE+0.01*X)*(np.abs(Sq[:, :, k]) - > 0.001) + 10*varPHASE0 + np.pi*Sq[:, :, k]/VENC - - if scantype == '-G+G': - PHASE0[:, :, k] = gamma*B0*TE * \ - np.ones([row, col]) + 10*varPHASE0 - np.pi*Sq[:, :, k]/VENC - PHASE1[:, :, k] = gamma*B0*TE * \ - np.ones([row, col]) + 10*varPHASE0 + np.pi*Sq[:, :, k]/VENC - - RHO0[:, :, k] = modulus*np.cos(PHASE0[:, :, k]) + \ - Drho + 1j*modulus*np.sin(PHASE0[:, :, k]) + 1j*Drho2 - RHO1[:, :, k] = modulus*np.cos(PHASE1[:, :, k]) + \ - Drho + 1j*modulus*np.sin(PHASE1[:, :, k]) + 1j*Drho2 - - if np.ndim(Sq) == 4: - [row, col, dep, numt2] = Sq.shape - [X, Y, Z] = np.meshgrid(np.linspace(0, col, col), np.linspace( - 0, row, row), np.linspace(0, dep, dep)) - - for k in range(numt2): - - if noise: - Drho = np.random.normal(0, 0.2, [row, col, dep]) - Drho2 = np.random.normal(0, 0.2, [row, col, dep]) - else: - Drho = np.zeros([row, col, dep]) - Drho2 = np.zeros([row, col, dep]) - - varPHASE0 = np.random.randint(-10, 11, size=(row, col, dep)) * \ - np.pi/180*(np.abs(Sq[:, :, :, k]) < 0.001) - modulus = 0.5 + 0.5*(np.abs(Sq[:, :, :, k]) > 0.001) - - if scantype == '0G': - PHASE0[:, :, :, k] = (gamma*B0*TE+0.01*X) * \ - (np.abs(Sq[:, :, :, k]) > 0.001) + 10*varPHASE0 - PHASE1[:, :, :, k] = (gamma*B0*TE+0.01*X)*(np.abs(Sq[:, :, :, k]) - > 0.001) + 10*varPHASE0 + np.pi*Sq[:, :, :, k]/VENC - - if scantype == '-G+G': - PHASE0[:, :, :, k] = gamma*B0*TE * \ - np.ones([row, col, dep]) + varPHASE0 - \ - np.pi*Sq[:, :, :, k]/VENC - PHASE1[:, :, :, k] = gamma*B0*TE * \ - np.ones([row, col, dep]) + varPHASE0 + \ - np.pi*Sq[:, :, :, k]/VENC - - RHO0[:, :, :, k] = modulus*np.cos(PHASE0[:, :, :, k]) + \ - Drho + 1j*modulus*np.sin(PHASE0[:, :, :, k]) + 1j*Drho2 - RHO1[:, :, :, k] = modulus*np.cos(PHASE1[:, :, :, k]) + \ - Drho + 1j*modulus*np.sin(PHASE1[:, :, :, k]) + 1j*Drho2 - - return [RHO0, RHO1] - -def undersampling(Sqx, Sqy, Sqz, options, savepath): - - R = options['cs']['R'] - - for r in R: - - if rank == 0: - print('Using Acceleration Factor R = ' + str(r)) - print('Component x of M0') - - [M0, M1] = GenerateMagnetization( - Sqx, options['cs']['VENC'], options['cs']['noise']) - - print('\n Component x of M0') - M0_cs = CSMETHOD(M0, r) - print('\n Component x of M1') - M1_cs = CSMETHOD(M1, r) - - Sqx_cs = phase_contrast(M1_cs, M0_cs, options['cs']['VENC']) - del M0, M1 - del M0_cs, M1_cs - - [M0, M1] = GenerateMagnetization( - Sqy, options['cs']['VENC'], options['cs']['noise']) - - print('\n Component y of M0') - M0_cs = CSMETHOD(M0, r) - print('\n Component y of M1') - M1_cs = CSMETHOD(M1, r) - - Sqy_cs = phase_contrast(M1_cs, M0_cs, options['cs']['VENC']) - - del M0, M1 - del M0_cs, M1_cs - - [M0, M1] = GenerateMagnetization( - Sqz, options['cs']['VENC'], options['cs']['noise']) - - if rank == 0: - print('\n Component z of M0') - M0_cs = CSMETHOD(M0, r) - if rank == 0: - print('\n Component z of M1') - M1_cs = CSMETHOD(M1, r) - if rank == 0: - print(' ') - - Sqz_cs = phase_contrast(M1_cs, M0_cs, options['cs']['VENC']) - - if rank == 0: - print('saving the sequences in ' + savepath) - seqname = options['cs']['name'] + '_R' + str(r) + '.npz' - print('sequence name: ' + seqname) - np.savez_compressed(savepath + seqname, - x=Sqx_cs, y=Sqy_cs, z=Sqz_cs) - - del Sqx_cs, Sqy_cs, Sqz_cs - -def undersampling_short(Mx, My, Mz, options): - - R = options['cs']['R'] - savepath = options['cs']['savepath'] - - R_SENSE = 1 - if 'R_SENSE' in options['cs']: - R_SENSE = options['cs']['R_SENSE'][0] - - for r in R: - if rank == 0: - print('Using Acceleration Factor R = ' + str(r)) - - if R_SENSE == 2: - [MxS0_cs, MxS1_cs] = CSMETHOD_SENSE(Mx, r, 2) - [MyS0_cs, MyS1_cs] = CSMETHOD_SENSE(My, r, 2) - [MzS0_cs, MzS1_cs] = CSMETHOD_SENSE(Mz, r, 2) - if rank == 0: - print('saving the sequences in ' + savepath) - seqname_s0 = options['cs']['name'] + 'S0_R' + str(r) + '.npz' - seqname_s1 = options['cs']['name'] + 'S1_R' + str(r) + '.npz' - print('sequence name: ' + seqname_s0) - np.savez_compressed(savepath + seqname_s0, - x=MxS0_cs, y=MyS0_cs, z=MzS0_cs) - print('sequence name: ' + seqname_s1) - np.savez_compressed(savepath + seqname_s1, - x=MxS1_cs, y=MyS1_cs, z=MzS1_cs) - del MxS0_cs, MyS0_cs, MzS0_cs - del MxS1_cs, MyS1_cs, MzS1_cs - elif R_SENSE == 1: - Mx_cs = CSMETHOD(Mx, r) - My_cs = CSMETHOD(My, r) - Mz_cs = CSMETHOD(Mz, r) - if rank == 0: - print('saving the sequences in ' + savepath) - seqname = options['cs']['name'] + '_R' + str(r) + '.npz' - print('sequence name: ' + seqname) - np.savez_compressed(savepath + seqname, - x=Mx_cs, y=My_cs, z=Mz_cs) - del Mx_cs, My_cs, Mz_cs - else: - raise Exception('Only implemented for 2-fold SENSE!!') - - -# THE END diff --git a/codes/Graphics.py b/codes/Graphics.py deleted file mode 100644 index 89bf6ed..0000000 --- a/codes/Graphics.py +++ /dev/null @@ -1,1412 +0,0 @@ -import h5py -from dolfin import * -import numpy as np -from matplotlib import pyplot as plt -import scipy as sc -import scipy.io as sio -from mpl_toolkits.mplot3d import Axes3D -from mpl_toolkits.mplot3d import proj3d -import os, sys -import csv -from common import inout -import time - -from matplotlib import rc -#rc('font',**{'family':'sans-serif','sans-serif':['Helvetica']}) -rc('text', usetex=True) - - -def MATmedical(path): - - - datatype = None - - if 'ODV' in path: - datatype = 'odv' - print('The DATA file has the next structure: [x,y,z,coil,CardiacPhase,gradient,NSA] where x \ - and y are the spatial dimensions for the image. z is one because is only one slice. Coil is \ - the number of coils (5), CardiacPhase is the time dimension. Gradient is the number of gradients used in the \ - adquisition (2). NSA is the number of signal average.') - - if 'Phantom' in path: - datatype = 'cib_phantom' - #print('CIB phantom data which contains magnitude and velocity measurements.') - - - if not datatype=='cib_phantom': - STRUC = sio.loadmat(path) - DAT0 = STRUC['data'] - if datatype=='odv': - KDAT0 = STRUC['kdata'] - [row,col,coil,numt2,p4,p5] = DAT0.shape - else: - KDAT0 = STRUC['k_data'] - [row,col,slicesel,coil,numt2,p4,p5] = DAT0.shape - - - z0 = 0 - NSA = 0 - coilnum = 3 - - DAT1 = np.zeros([row,col,numt2],dtype=complex) - DAT2 = np.zeros([row,col,numt2],dtype=complex) - KDAT1 = np.zeros([row,col,numt2],dtype=complex) - KDAT2 = np.zeros([row,col,numt2],dtype=complex) - - if datatype=='odv': - for k in range(numt2): - DAT1[:,:,k] = DAT0[:,:,coilnum,k,0,NSA] - DAT2[:,:,k] = DAT0[:,:,coilnum,k,1,NSA] - KDAT1[:,:,k] = KDAT0[:,:,coilnum,k,0,NSA] - KDAT2[:,:,k] = KDAT0[:,:,coilnum,k,1,NSA] - - return [DAT1,DAT2,KDAT1,KDAT2] - else: - - UK1 = np.zeros([row,col,numt2],dtype=complex) - UK2 = np.zeros([row,col,numt2],dtype=complex) - - for k in range(numt2): - DAT1[:,:,k] = DAT0[:,:,z0,coilnum,k,0,NSA] - DAT2[:,:,k] = DAT0[:,:,z0,coilnum,k,1,NSA] - KDAT1[:,:,k] = KDAT0[:,:,z0,coilnum,k,0,NSA] - KDAT2[:,:,k] = KDAT0[:,:,z0,coilnum,k,1,NSA] - UK1[:,:,k] = np.fft.fftshift(np.fft.fft2(DAT1[:,:,k])) - UK2[:,:,k] = np.fft.fftshift(np.fft.fft2(DAT2[:,:,k])) - - return [DAT1,DAT2,KDAT1,KDAT2] - - elif datatype=='cib_phantom': - - lp = len(path) - for k in range(lp): - if path[lp-1-k]=='/': - name_var = path[lp-k:lp] - break - - name_var = name_var.replace('.mat','') - - if name_var in ['VENC','Segmented_Aorta','voxel_MR']: - STRUC = sio.loadmat(path) - return STRUC[name_var] - else: - #import h5py - arrays = {} - f = h5py.File(path) - for k, v in f.items(): - arrays[k] = np.array(v) - - return arrays[name_var].transpose((1,2,3,0)) - -def Plot_Inlet(): - ''' To Show the flow given at the inlet in the FEM simulation ''' - - Tf = 0.9 - dt = 0.0005 - t = np.linspace(0,Tf,Tf/dt) - U = 60 - Th = 0.37 - DOLFIN_PI = 3.1416 - y = U*np.sin(DOLFIN_PI*t/Th)*(t<=Th) + (Th0: - plt.ylim(ywin) - plt.xlim(xwin) - ax.set_aspect('auto') - - plt.xticks([]) - plt.yticks([]) - #plt.xlabel('$x$',Fontsize=20) - #plt.ylabel('$y$',Fontsize=20) - plt.title('phase ' + str(k)) - plt.show() - - if recorder: - plt.savefig('results/pictures/frame'+str(k),bbox_inches='tight') - - plt.pause(0.101) - plt.clf() - -def PrintSequence2(A,X,Y,Z,x,y,z): - - plt.ion() - fig = plt.figure(figsize=(7, 7), dpi=100) - ax = fig.add_subplot(1,1,1) - [row,col,dep] = A.shape - - - xi = np.linspace(np.min(x),np.max(x),row) - yi = np.linspace(np.min(y),np.max(y),col) - yi2 = np.min(yi) + np.max(yi) - yi - - from matplotlib.colors import LinearSegmentedColormap - basic_cols = ['#ff4e2b', '#000000', '#1893db'] - my_cmap = LinearSegmentedColormap.from_list('mycmap', basic_cols) - ccmap = mpl.cm.Spectral - - xt = x - yt = y - zt = z - yt2 = (np.min(y) + np.max(y)) - y - Dz = (np.max(z) - np.min(z))/dep - - def shiftedColorMap(cmap, start=0, midpoint=0.5, stop=1.0, name='shifted'): - ''' - Function to offset the "center" of a colormap. Useful for - data with a negative min and positive max and you want the - middle of the colormap's dynamic range to be at zero. - - Input - ----- - cmap : The matplotlib colormap to be altered - start : Offset from lowest point in the colormap's range. - Defaults to 0.0 (no lower offset). Should be between - 0.0 and `midpoint`. - midpoint : The new center of the colormap. Defaults to - 0.5 (no shift). Should be between 0.0 and 1.0. In - general, this should be 1 - vmax / (vmax + abs(vmin)) - For example if your data range from -15.0 to +5.0 and - you want the center of the colormap at 0.0, `midpoint` - should be set to 1 - 5/(5 + 15)) or 0.75 - stop : Offset from highest point in the colormap's range. - Defaults to 1.0 (no upper offset). Should be between - `midpoint` and 1.0. - ''' - - cdict = { - 'red': [], - 'green': [], - 'blue': [], - 'alpha': [] - } - - # regular index to compute the colors - reg_index = np.linspace(start, stop, 257) - - # shifted index to match the data - shift_index = np.hstack([ - np.linspace(0.0, midpoint, 128, endpoint=False), - np.linspace(midpoint, 1.0, 129, endpoint=True) - ]) - - for ri, si in zip(reg_index, shift_index): - r, g, b, a = cmap(ri) - - cdict['red'].append((si, r, r)) - cdict['green'].append((si, g, g)) - cdict['blue'].append((si, b, b)) - cdict['alpha'].append((si, a, a)) - - newcmap = mpl.colors.LinearSegmentedColormap(name, cdict) - plt.register_cmap(cmap=newcmap) - - return newcmap - - ccmap = mpl.cm.Spectral - shifted_ccmap = shiftedColorMap(ccmap, midpoint=0.75, name='shifted') - - for k in range(dep): - - zmin = np.min(z) + k*Dz - zmax = zmin + Dz - xs = np.array([]) - ys = np.array([]) - zs = np.array([]) - - for j in range(z.size): - if z[j]>= zmin and z[j] 1: - if os.path.exists(sys.argv[1]): - inputfile = sys.argv[1] - print('Found input file ' + inputfile) - else: - raise Exception('Command line arg given but input file does not exist:' - ' {}'.format(sys.argv[1])) - else: - raise Exception('An input file is required as argument!') - - - start_time = time.time() - - - options = inout.read_parameters(inputfile) - ROUTINE(options) - end_time = time.time() - CLOCK(start_time,end_time) - - - -# END - - - - - - - diff --git a/codes/MATLAB/createU.m b/codes/MATLAB/createU.m deleted file mode 100644 index 67eb34e..0000000 --- a/codes/MATLAB/createU.m +++ /dev/null @@ -1,57 +0,0 @@ -clear all; close all - -folder_name = uigetdir([],'Load Folder...'); - -data = load(strcat(folder_name,'/data.mat')); -SEG = load(strcat(folder_name,'/SEG.mat')); - -data = data.data; -SEG = SEG.SEG; - - -VENC = data.VENC; -VoxelSize = data.voxel_MR; - -vel_AP = data.MR_PCA_AP; -vel_RL = data.MR_PCA_RL; -vel_FH = data.MR_PCA_FH; - -SEG2 = permute(SEG,[2,3,1]); -SEG2 = SEG2(:,:,:); - - -vel_AP_seg = vel_AP.*SEG2(2:end-1,2:end-1,2:end-1); -vel_RL_seg = vel_RL.*SEG2(2:end-1,2:end-1,2:end-1); -vel_FH_seg = vel_FH.*SEG2(2:end-1,2:end-1,2:end-1); - - - - -u_R1 = [] ; -u_R1.x = vel_FH_seg; -u_R1.y = vel_AP_seg; -u_R1.z = vel_RL_seg; -u_R1.VoxelSize = VoxelSize; -save('/home/yeye/Desktop/u_R1.mat','u_R1'); -disp('data saved') -%% -%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% FIGURES -%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -figure -size_vel = size(vel_FH); -for n=1:size_vel(3) - imshow(squeeze(vel_FH_seg(:,:,n,8)),[-100,100],'InitialMagnification',300); - colormap(gca); - pause(0.1) -end -%% -size_seg2 = size(SEG2); -for n=1:size_seg2(3) - imshow(squeeze(SEG2(:,:,n)),'InitialMagnification',300); - colormap(gca); - pause(0.1) -end - - diff --git a/codes/MATLAB/leo/CREATE_MESH.m b/codes/MATLAB/leo/CREATE_MESH.m deleted file mode 100755 index 47bab09..0000000 --- a/codes/MATLAB/leo/CREATE_MESH.m +++ /dev/null @@ -1,14 +0,0 @@ -% Program to create a structured mesh using the codes of Leo Sok -clear all; close all - -nodes = load('LEO_files/nodes.txt'); -ux = load('LEO_files/ux.txt') ; -uy = load('LEO_files/uy.txt') ; -uz = load('LEO_files/uz.txt') ; -u = sqrt(ux.^2 + uy.^2 + uz.^2); -resol = load('LEO_files/resol.txt') ; -dx = resol(1); dy = resol(2) ; dz = resol(3); - -nodes_masked = maskFEM(nodes,u); -[N,tets,faces] = meshStructTess(nodes_masked,dx,dy,dz,0,0); -writemesh('/home/yeye/Desktop/leomesh',N,tets,faces) \ No newline at end of file diff --git a/codes/MATLAB/leo/maskFEM.m b/codes/MATLAB/leo/maskFEM.m deleted file mode 100755 index b7dca59..0000000 --- a/codes/MATLAB/leo/maskFEM.m +++ /dev/null @@ -1,19 +0,0 @@ -function nodes2 = maskFEM(nodes,vel) - -a = []; -b = []; -c = []; -ind = 1; - -for i=1:length(nodes) -if vel(i)>0 - a(ind) = nodes(i,1); - b(ind) = nodes(i,2); - c(ind) = nodes(i,3); - ind = ind +1; -end -end - -nodes2 = [a', b', c']; - - diff --git a/codes/MATLAB/leo/meshStructTess.m b/codes/MATLAB/leo/meshStructTess.m deleted file mode 100755 index 4e23f07..0000000 --- a/codes/MATLAB/leo/meshStructTess.m +++ /dev/null @@ -1,169 +0,0 @@ -function [nodes, tets, faces, P] = meshStructTess(nodes, dx, dy, dz, check_mesh, plot_mesh) -%% [nodes, tets, faces] = meshStructTess(nodes, dx, dy, dz, check_mesh, plot_mesh) -% Generate a tessalation from a list of structured nodes. -% input: nodes: n times 3 matrix with on the rows the coordinates of -% the n points in the mesh -% dx, dy, dz: the mesh-size in the directions x, y and z -% check_mesh: if true, then it solves a Poisson problem -% plot_mesh: if true, then it plots the mesh -% output: nodes: m times 3 matrix with on the rows the coordinates of -% the m <= n points in the triangulationedi -% tets: l times 4 matrix with on the rows the tetrahedra -% faces: k times 3 matrix with on the rows the triangles of the -% boundary of the mesh -% P: Transformation matrix from input nodes to output nodes. -% Useful also for transforming node-valued functions on -% the input nodes to node-valued functions on the output -% nodes -% -% The triangulation can be plotted using tetramesh(tets,nodes) - - -% compute the minimum and number of points in each direction -if size(nodes,1) < 4 - error('Triangulation needs at least 4 points') -end -mn = min(nodes); -xmin = mn(1); -ymin = mn(2); -zmin = mn(3); - -mn = max(nodes); -xmax = mn(1); -ymax = mn(2); -zmax = mn(3); - -nx = round((xmax-xmin)/dx +1); -ny = round((ymax-ymin)/dy +1); -nz = round((zmax-zmin)/dz +1); - -Nnodes = size(nodes,1); - - -% Define tensor which consist of nodes indices, used for the creation of -% the tetrahedra - -nodes3d = zeros(nx,ny,nz); % preallocate -for i=1:Nnodes - nodes3d(round((nodes(i,1)-xmin)/dx)+1,round((nodes(i,2)-ymin)/dy)+1,round((nodes(i,3)-zmin)/dz)+1)=i; -end - - -disp('Creating Tetrahedra') - -% create tetrahedral mesh in cube, which we will reuse. -ii = 1; -X = zeros(8,3); -for i=0:1 - for j=0:1 - for k=0:1 - X(ii,:) = [i,j,k]; - ii = ii+1; - end - end -end -cubetet = delaunay(X); - -% Run through the mesh -el = 1; -Tetrahedra = zeros(6*(nnz(nodes3d)),4); % preallocate - -for i=1:nx-1 - for j=1:ny-1 - for k=1:nz-1 - % take [i:i+1,j:j+1,k:k+1] as cube - nod = zeros(1,8); % perallocate - - for l = 1:8 - % nod is vector with node indices of cube - nod(l) = nodes3d(i + X(l,1), j + X(l,2), k + X(l,3)); - end - - if nnz(nod) == 8 % then the cube is inside the mesh - tet = nod(cubetet); - else % then there is at least one point of the cube outside the mesh - Xs = X(logical(nod),:); % take only nodes inside the mesh - nodx = nod(logical(nod)); - if nnz(nod) == 4 % 4 nodes, check if points are coplanar - C = cross(Xs(2,:)-Xs(1,:), Xs(3,:)-Xs(1,:)); - cop = logical(dot(C,Xs(4,:)-Xs(1,:))); - % if cop = 0, then points are coplanar end thus no - % tetrahedra exists. - end - if (nnz(nod)>4) || (nnz(nod) == 4 && cop) - % create tetrahedra - tet1 = delaunay(Xs); - tet = nodx(tet1); - else % no tetrahedra exists - tet = []; - end - end - - % add new tetrahedra to list - Tetrahedra(el:el+size(tet,1)-1,:) = tet; - el = el+size(tet,1); - end - end -end - -tets = Tetrahedra(1:el-1,:); % Delete extra preallocated rows. -clear Tetrahedra - -disp([num2str(size(tets,1)), ' tetrahedra created']) - -% Delete nodes which are not in any tetrahedra. -disp('Update mesh') -contr = zeros(size(nodes,1),1); -for i=1:size(tets,1) - for j=1:4 - contr(tets(i,j))=1; - end -end - -nodes = nodes(logical(contr),:); - -% compute P -P = speye(Nnodes); -P = P(logical(contr),:); - -disp([num2str(nnz(~contr)), ' unused nodes in triangulation deleted.']) - -disp('Update tetrahedra') - -% make tetrahedra compatible with new node indices -cumcon = cumsum(~contr)'; -tets = tets - cumcon(tets); - -% create triangles -if size(tets,1) == 0 - warning('No tetrahedra created') - faces = zeros(0,3); -else - disp('Create Triangles') - faces = freeBoundary(triangulation(tets,nodes)); - disp([num2str(size(faces,1)), ' triangles created']) -end - -% checking the mesh by solving a Poisson problem -if check_mesh - % Builds the P1 stiffness matrix from tets and nodes - [A,volumes]=stifness_matrixP1_3D(tets,nodes); - % Check if element volumes may be negative - if any(volumes<=0) - warning('Some elements have zero or negative volume') - end - % solve the Poisson problem with Dirichlet BC - A(2:end,2:end)\ones(size(A(2:end,2:end),1),1); - disp('If there are no warnings, it probably means that the mesh is fine') -end - -% Plots mesh -if plot_mesh - tetramesh(tets,nodes) - xlabel('x') - ylabel('y') - zlabel('z') -end - -end - diff --git a/codes/MATLAB/leo/writemesh.m b/codes/MATLAB/leo/writemesh.m deleted file mode 100755 index a79ad0a..0000000 --- a/codes/MATLAB/leo/writemesh.m +++ /dev/null @@ -1,97 +0,0 @@ -function writemesh(varargin) -%% writemesh(path, mesh) -% Save triangulation as path.xml and path.msh -% mesh is a struct with fields Pts, Tet, Tri -% alernatively one can use writemesh(path, Pts, Tet, Tri) -% Pts should by a n times 3 matrix consisting points of the mesh -% Tet is the m times 4 matrix consisting the tetrahedra -% Tri is the l times 3 matrix consisting the triangles at the boundary - -if nargin > 3 - mesh.Pts=varargin{2}; - mesh.Tet=varargin{3}; - mesh.Tri=varargin{4}; - writemesh(varargin{1},mesh,varargin(nargin)); - -elseif isstruct(varargin{2}) - rootMeshFile = varargin{1}; - - % NEW FILE - obj = [rootMeshFile,'.msh']; - meshfile = fopen(obj,'w'); - - obj2 = [rootMeshFile,'.xml']; - xmlfile = fopen(obj2,'w'); - - % MESH - fprintf(meshfile,['$MeshFormat','\n']); - fprintf(meshfile,['2.2 0 8','\n']); - fprintf(meshfile,['$EndMeshFormat','\n']); - - fprintf(xmlfile,['','\n']); - fprintf(xmlfile,'\n'); - fprintf(xmlfile,['','\n']); - - mesh = varargin{2}; - - Nodes = mesh.('Pts'); - mesh = rmfield(mesh,'Pts'); - - Nodes = [(1:size(Nodes,1))' Nodes(:,1:3)]; - - % POINTS - if ~strcmp(varargin{nargin},'mute') - disp('Write Points') - end - fprintf(meshfile,['$Nodes','\n']); - fprintf(meshfile,['%i','\n'],size(Nodes,1)); - fprintf(xmlfile,[' ','\n']); - fprintf(xmlfile,[' ','\n'],size(Nodes,1)); - - - fprintf(meshfile,'%i %13.6f %13.6f %13.6f\n',Nodes'); - - Nodes(:,1) = Nodes(:,1) - 1; - - fprintf(xmlfile,' \n',Nodes'); - - fprintf(meshfile,['$EndNodes','\n']); - fprintf(meshfile,['$Elements','\n']); - fprintf(meshfile,['%i','\n'],size(mesh.Tet,1)+size(mesh.Tri,1)); - fprintf(xmlfile,[' ','\n']); - fprintf(xmlfile,[' ','\n'],size(mesh.Tet,1)); - - % Triangles - - if ~strcmp(varargin{nargin},'mute') - disp('Write Triangles') - end - - tri = mesh.('Tri'); - tri = [(1:size(tri,1))' 2*ones(size(tri,1),1) 2*ones(size(tri,1),1) zeros(size(tri,1),1) 2*ones(size(tri,1),1) tri(:,1:3)]; - fprintf(meshfile,'%i %i %i %i %i %i %i %i\n',tri'); - - - - % Tetrahedra - if ~strcmp(varargin{nargin},'mute') - disp('Write Tetrahedra') - end - - tet = mesh.('Tet'); - tet = [(size(tri,1)+1:size(tri,1)+size(tet,1))' 4*ones(size(tet,1),1) 2*ones(size(tet,1),1) zeros(size(tet,1),1) ones(size(tet,1),1) tet(:,1:4)]; - fprintf(meshfile,'%i %i %i %i %i %i %i %i %i\n',tet'); - - tet = mesh.('Tet'); - tet = [(0:size(tet,1)-1)' (tet(:,1:4)-1)]; - fprintf(xmlfile,' \n',tet'); - - - - fprintf(meshfile,['$EndElements','\n']); - fprintf(xmlfile,' \n \n\n'); - - fclose('all'); -end - - diff --git a/codes/MATLAB/load_dicom.m b/codes/MATLAB/load_dicom.m deleted file mode 100644 index 35eb488..0000000 --- a/codes/MATLAB/load_dicom.m +++ /dev/null @@ -1,126 +0,0 @@ -clear all ; close all -% Load dicom - - -name = 'Ronald' ; - -if strcmp(name, 'Ronald') - path_all = [ - '/home/yeye/Desktop/PhD/MEDICAL_DATA/DatosSEPT2019/20190909_Ronald/FH/DICOM/IM_0001', - '/home/yeye/Desktop/PhD/MEDICAL_DATA/DatosSEPT2019/20190909_Ronald/AP/DICOM/IM_0001', - '/home/yeye/Desktop/PhD/MEDICAL_DATA/DatosSEPT2019/20190909_Ronald/RL/DICOM/IM_0001' - ] ; -end - -if strcmp(name, 'Jeremias') - path_all = [ - '/home/yeye/Desktop/PhD/MEDICAL_DATA/DatosSEPT2019/20190909_Jeremias/FH/DICOM/IM_0001', - '/home/yeye/Desktop/PhD/MEDICAL_DATA/DatosSEPT2019/20190909_Jeremias/AP/DICOM/IM_0001', - '/home/yeye/Desktop/PhD/MEDICAL_DATA/DatosSEPT2019/20190909_Jeremias/RL/DICOM/IM_0001' - ] ; -end - -if strcmp(name, 'Hugo') - path_all = [ - '/home/yeye/Desktop/PhD/MEDICAL_DATA/DatosSEPT2019/20190924_Hugo/Dicom/DICOM/IM_0013', - '/home/yeye/Desktop/PhD/MEDICAL_DATA/DatosSEPT2019/20190924_Hugo/Dicom/DICOM/IM_0009', - '/home/yeye/Desktop/PhD/MEDICAL_DATA/DatosSEPT2019/20190924_Hugo/Dicom/DICOM/IM_0005' - ] ; -end - -for i=1:3 - - if i==1 - %path = '/home/yeye/Desktop/PhD/MEDICAL_DATA/DatosSEPT2019/20190924_Paloma/Dicom/DICOM/IM_0013' - disp('Reading the FH component from ...') - path = path_all(1,:) - end - - if i==2 - %path = '/home/yeye/Desktop/PhD/MEDICAL_DATA/DatosSEPT2019/20190924_Paloma/Dicom/DICOM/IM_0009' ; - disp('Reading the AP component from ...') - path = path_all(2,:) - end - - if i==3 - %path = '/home/yeye/Desktop/PhD/MEDICAL_DATA/DatosSEPT2019/20190924_Paloma/Dicom/DICOM/IM_0005' ; - disp('Reading the RL component from ...') - path = path_all(3,:) - end - - -I_info = dicominfo(path); -I = double(dicomread(path)); -VENC = eval(['I_info.PerFrameFunctionalGroupsSequence.Item_1.MRVelocityEncodingSequence.Item_1.VelocityEncodingMaximumValue']) ; -heart_rate = eval(['I_info.PerFrameFunctionalGroupsSequence.Item_1.Private_2005_140f.Item_1.HeartRate']); - - -MAG = zeros(size(I,1),size(I,2),I_info.Private_2001_1018,I_info.Private_2001_1017); -PHASE = zeros(size(I,1),size(I,2),I_info.Private_2001_1018,I_info.Private_2001_1017); - -for n=1:size(I,4) - - RI = eval(['I_info.PerFrameFunctionalGroupsSequence.Item_',num2str(n),'.Private_2005_140f.Item_1.RescaleIntercept']); % intercept - RS = eval(['I_info.PerFrameFunctionalGroupsSequence.Item_',num2str(n),'.Private_2005_140f.Item_1.RescaleSlope']); % slope - cp = eval(['I_info.PerFrameFunctionalGroupsSequence.Item_',num2str(n),'.Private_2005_140f.Item_1.Private_2001_1008']); %cp - slc = eval(['I_info.PerFrameFunctionalGroupsSequence.Item_',num2str(n),'.Private_2005_140f.Item_1.Private_2001_100a']); %scl - id = eval(['I_info.PerFrameFunctionalGroupsSequence.Item_',num2str(n),'.Private_2005_140f.Item_1.Private_2005_106e']); % PCA o FFE - - if strcmp(id,'FFE')==1 - MAG(:,:,slc,cp) = I(:,:,1,n)*RS + RI; - else - PHASE(:,:,slc,cp) = I(:,:,1,n)*RS + RI; - end - -end - - -MASK = double(abs((PHASE==PHASE(1,1,1,1))-1)); -PHASE = PHASE.*MASK; - - -if i==1 - MR_FFE_FH = MAG; - MR_PCA_FH = VENC*PHASE/pi/100; -end - -if i==2 - MR_FFE_AP = MAG; - MR_PCA_AP = VENC*PHASE/pi/100; -end -if i==3 - MR_FFE_RL = MAG; - MR_PCA_RL = VENC*PHASE/pi/100; -end - - -end - - -disp('Saving the data ...') - -spaceslices = eval(['I_info.PerFrameFunctionalGroupsSequence.Item_1.PixelMeasuresSequence.Item_1.SpacingBetweenSlices']); -pixelspacing = eval(['I_info.PerFrameFunctionalGroupsSequence.Item_1.PixelMeasuresSequence.Item_1.PixelSpacing']); - -disp('voxel-size recognized:') -voxel_MR = [pixelspacing(1),pixelspacing(1),spaceslices] - - -data = []; -data.MR_FFE_AP = MR_FFE_AP; -data.MR_FFE_RL = MR_FFE_RL; -data.MR_FFE_FH = MR_FFE_FH; -data.MR_PCA_AP = MR_PCA_AP; -data.MR_PCA_RL = MR_PCA_RL; -data.MR_PCA_FH = MR_PCA_FH; -data.type = 'DAT'; -data.VENC = VENC ; -data.voxel_MR = voxel_MR; -data.heart_rate = heart_rate; - -save('/home/yeye/Desktop/data.mat','data','-v7.3'); -disp('data saved') - - - - \ No newline at end of file diff --git a/codes/MRI.py b/codes/MRI.py deleted file mode 100644 index 9707027..0000000 --- a/codes/MRI.py +++ /dev/null @@ -1,1755 +0,0 @@ -################################################################ -# -# Workspace for MRI analysis of the 4Dflow data -# -# written by Jeremias Garay L: j.e.garay.labra@rug.nl -# for autoreload in ipython3 -# %load_ext autoreload -# %autoreload 2 -################################################################ -from dolfin import * -import dolfin -import numpy as np -import sys, os -from mpi4py import MPI -from common import inout -import time - - -PRINT_BARS = False -parameters["std_out_all_processes"] = False - -def Etcetera(mode): - if mode=='happy': - print(r""" - /\ /\ - //\\_//\\ ____ - \_ _/ / / - / ^ ^ \ /^^^] - \_\O/_/ [ ] - / \_ [ / - \ \_ / / - [ [ / \/ _/ - _[ [ \ /_/ - """) - - if mode=='serious': - print(r""" - |\__/| - / \ - /_.- -,_\ - \@/ - """) - -def dealiased(VENC,vx,vy,vz): - - vx2 = np.zeros(vx.shape) - vy2 = np.zeros(vy.shape) - vz2 = np.zeros(vz.shape) - - badx_pos = np.where(vx > VENC) - badx_neg = np.where(vx < -VENC) - bady_pos = np.where(vy > VENC) - bady_neg = np.where(vy < -VENC) - badz_pos = np.where(vz > VENC) - badz_neg = np.where(vz < -VENC) - # More than VENC - vx2[badx_pos] = -(2*VENC - vx[badx_pos]) - vy2[bady_pos] = -(2*VENC - vy[bady_pos]) - vz2[badz_pos] = -(2*VENC - vz[badz_pos]) - # Less than VENC - vx2[badx_neg] = 2*VENC + vx[badx_neg] - vy2[bady_neg] = 2*VENC + vy[bady_neg] - vz2[badz_neg] = 2*VENC + vz[badz_neg] - # The rest - otherx = np.where(vx2 == 0) - othery = np.where(vy2 == 0) - otherz = np.where(vz2 == 0) - vx2[otherx] = vx[otherx] - vy2[othery] = vy[othery] - vz2[otherz] = vz[otherz] - - return [vx2,vy2,vz2] - -def checkpoint_norm(MESH,fieldpath,mode): - V = MESH['FEM'] - v = Function(V) - path = fieldpath + 'checkpoint/' - # Checkpoints folders - unsort_indexes = os.listdir(path) - indexes = [int(x) for x in unsort_indexes] - indexes.sort() - norm = np.zeros([len(indexes)]) - j = 0 - comm = MPI.COMM_WORLD - - if mode=='ct': - # CT DATA - for l in indexes: - path1 = path + str(l) + '/u.h5' - hdf = HDF5File(MESH['mesh'].mpi_comm(),path1,'r') - hdf.read(v, 'u/vector_0') - hdf.close() - norm[j] = np.linalg.norm(v.vector().get_local()) - j = j+1 - - if mode=='measurement': - # MEASUREMENT - path2 = fieldpath + 'measurements/' - j = 0 - for l in indexes: - meapath = path2 + 'u' + str(l) + '.h5' - comm = MPI.COMM_WORLD - hdf = HDF5File(MESH['mesh'].mpi_comm(),meapath,'r') - hdf.read(v, 'u/vector_0') - hdf.close() - norm[j] = np.linalg.norm(v.vector().get_local()) - j = j+1 - - - if mode=='roukf': - indexes2 = indexes[1:-1] - j = 1 - for l in indexes2: - path3 = fieldpath + 'roukf/checkpoint/' + str(l) + '/X0.h5' - hdf = HDF5File(MESH['mesh'].mpi_comm(),path3,'r') - hdf.read(v, 'X') - hdf.close() - norm[j] = np.linalg.norm(v.vector().get_local()) - j = j+1 - - - return norm - -def LOADmesh(meshtype,resol=None,boxsize=None,meshpath=None, ranges=None): - - if meshtype=='box_zeros': - # The size of the mesh is given - # The spacements are given - # The x range goes from 0 - [Lx, Ly, Lz] = boxsize - [ddx, ddy, ddz] = resol - x0 = 0 - y0 = 0 - z0 = 0 - x1 = (Lx-1)*ddx - y1 = (Ly-1)*ddy - z1 = (Lz-1)*ddz - mesh0 = BoxMesh(Point(x0, y0, z0), Point(x1, y1, z1), Lx-1, Ly-1, Lz-1) - P1 = FunctionSpace(mesh0, "Lagrange", 1) - # This need to be saved from Leo's construction - Xt = P1.tabulate_dof_coordinates() - #boundaries = MeshFunction('size_t', mesh0 , mesh0.topology().dim() - 1) - #XDMFFile(path2).write(boundaries) - if '2017' in dolfin.__version__: - Xt = np.reshape(Xt, (round(Xt.size/3), 3)) - SqXt = np.zeros(Xt.shape) - for k in range(Xt.shape[0]): - xk = round((Xt[k][0]-x0)/ddx) - yk = round((Xt[k][1]-y0)/ddy) - zk = round((Xt[k][2]-z0)/ddz) - SqXt[k, 0] = int(xk) - SqXt[k, 1] = int(yk) - SqXt[k, 2] = int(zk) - if rank == 0: - print('BOX real created: H = ' + str(mesh0.hmin()) + - ' H_max = ' + str(mesh0.hmax())) - - BOX = {} - BOX['mesh'] = mesh0 - BOX['FEM'] = P1 - BOX['nodes'] = Xt - BOX['seq'] = SqXt.astype(int) - BOX['dimension'] = [Lx,Ly,Lz] - return BOX - - if meshtype=='box_ranges': - # The size of the mesh is given - # The X ranges are given - # The spacements are calculated - [Lx,Ly,Lz] = boxsize - x0 = ranges['x'][0] - x1 = ranges['x'][1] - y0 = ranges['y'][0] - y1 = ranges['y'][1] - z0 = ranges['z'][0] - z1 = ranges['z'][1] - ddx = (x1-x0)/(Lx-1) - ddy = (y1-y0)/(Ly-1) - ddz = (z1-z0)/(Lz-1) - mesh0 = BoxMesh(Point(x0, y0, z0), Point(x1, y1, z1), Lx-1, Ly-1, Lz-1) - P1 = FunctionSpace(mesh0, "Lagrange", 1) - Xt = P1.tabulate_dof_coordinates() - if '2017' in dolfin.__version__: - Xt = np.reshape(Xt, (round(Xt.size/3), 3)) - SqXt = np.zeros(Xt.shape) - for k in range(Xt.shape[0]): - xk = round((Xt[k][0]-x0)/ddx) - yk = round((Xt[k][1]-y0)/ddy) - zk = round((Xt[k][2]-z0)/ddz) - SqXt[k, 0] = int(xk) - SqXt[k, 1] = int(yk) - SqXt[k, 2] = int(zk) - if rank == 0: - print('BOX real created: H = ' + str(mesh0.hmin()) + - ' H_max = ' + str(mesh0.hmax())) - - BOX = {} - BOX['mesh'] = mesh0 - BOX['FEM'] = P1 - BOX['nodes'] = Xt - BOX['seq'] = SqXt.astype(int) - BOX['dimension'] = [Lx, Ly, Lz] - return BOX - - if meshtype=='fix_resolution': - # The size of the mesh is given - # The spacements are given - # The X ranges are given - # The X ranges are redefined from the oldest - [Lx,Ly,Lz] = boxsize - [ddx,ddy,ddz] = resol - x0 = ranges['x'][0] - x1 = ranges['x'][1] - y0 = ranges['y'][0] - y1 = ranges['y'][1] - z0 = ranges['z'][0] - z1 = ranges['z'][1] - deltax = (ddx*(Lx-1) - (x1 - x0))/2 - deltay = (ddy*(Ly-1) - (y1 - y0))/2 - deltaz = (ddz*(Lz-1) - (z1 - z0))/2 - x0_new = x0 - deltax - x1_new = x1 + deltax - y0_new = y0 - deltay - y1_new = y1 + deltay - z0_new = z0 - deltaz - z1_new = z1 + deltaz - - mesh0 = BoxMesh(Point(x0_new,y0_new,z0_new),Point(x1_new,y1_new,z1_new),Lx-1,Ly-1,Lz-1) - P1 = FunctionSpace(mesh0, "Lagrange", 1) - Xt = P1.tabulate_dof_coordinates() # This need to be saved from Leo's construction - #boundaries = MeshFunction('size_t', mesh0 , mesh0.topology().dim() - 1) - #XDMFFile('/home/yeye/Desktop/bound.xdmf').write(boundaries) - if '2017' in dolfin.__version__: - Xt = np.reshape(Xt,( round(Xt.size/3) ,3)) - SqXt = np.zeros(Xt.shape) - for k in range(Xt.shape[0]): - xk = round((Xt[k][0]-x0_new)/ddx) - yk = round((Xt[k][1]-y0_new)/ddy) - zk = round((Xt[k][2]-z0_new)/ddz) - SqXt[k, 0] = int(xk) - SqXt[k, 1] = int(yk) - SqXt[k, 2] = int(zk) - - if rank==0: - print('BOX real created: H = ' + str(mesh0.hmin()) + ' H_max = ' + str(mesh0.hmax()) ) - - BOX = {} - BOX['mesh'] = mesh0 - BOX['FEM'] = P1 - BOX['nodes'] = Xt - BOX['seq'] = SqXt.astype(int) - BOX['dimension'] = [Lx, Ly, Lz] - - return BOX - - if meshtype=='cib': - ddx = 2.00032/1000; ddy = 2.00032/1000; ddz = 2.00526316/1000 - x0 = ddx*17 ; y0 = ddy*17; z0 = ddz - Lx = 75; Ly = 85; Lz = 57 - x1 = x0+ddx*(Lx-1); y1 = y0+ddy*(Ly-1); z1 = z0+ddz*(Lz-1) - - - mesh0 = BoxMesh(Point(x0,y0,z0),Point(x1,y1,z1),Lx-1,Ly-1,Lz-1) - P1 = FunctionSpace(mesh0, "Lagrange", 1) - Xt = P1.tabulate_dof_coordinates() - #XDMFFile('mesh_hex.xdmf').write(mesh0) - # MESH HEXA - #mesh0 = BoxMesh.create([Point(x0,y0,z0), Point(x1,y1,z1)], [Lx, Ly, Lz], CellType.Type.hexahedron) - #P1 = FunctionSpace(mesh0, "DG", 0) - #t = P1.tabulate_dof_coordinates() - - if '2017' in dolfin.__version__: - Xt = np.reshape(Xt,( round(Xt.size/3) ,3)) - SqXt = np.zeros(Xt.shape) - - for k in range(Xt.shape[0]): - xk = round((Xt[k][0]-x0)/ddx) - yk = round((Xt[k][1]-y0)/ddy) - zk = round((Xt[k][2]-z0)/ddz) - SqXt[k,0] = int(xk) - SqXt[k,1] = int(yk) - SqXt[k,2] = int(zk) - if rank==0: - print('BOX cib created: H = ' + str(mesh0.hmin()) + ' H_max = ' + str(mesh0.hmax()) ) - - BOX = {} - BOX['mesh'] = mesh0 - BOX['FEM'] = P1 - BOX['nodes'] = Xt - BOX['seq'] = SqXt.astype(int) - BOX['name'] = 'cib' - BOX['dimension'] = [Lx,Ly,Lz] - return BOX - - if meshtype=='cib_RT': - ddx = 1.8/1000; ddy = 1.8/1000; ddz = 1.8/1000 - x0 = ddx*33 ; y0 = ddy*43; z0 = ddz*13 - Lx = 103; Ly = 63; Lz = 47 - x1 = x0+ddx*(Lx-1); y1 = y0+ddy*(Ly-1); z1 = z0+ddz*(Lz-1) - - mesh0 = BoxMesh(Point(x0,y0,z0),Point(x1,y1,z1),Lx-1,Ly-1,Lz-1) - P1 = FunctionSpace(mesh0, "Lagrange", 1) - Xt = P1.tabulate_dof_coordinates() - #XDMFFile('mesh_hex.xdmf').write(mesh0) - # MESH HEXA - #mesh0 = BoxMesh.create([Point(x0,y0,z0), Point(x1,y1,z1)], [Lx, Ly, Lz], CellType.Type.hexahedron) - #P1 = FunctionSpace(mesh0, "DG", 0) - #t = P1.tabulate_dof_coordinates() - - if '2017' in dolfin.__version__: - Xt = np.reshape(Xt,( round(Xt.size/3) ,3)) - SqXt = np.zeros(Xt.shape) - - for k in range(Xt.shape[0]): - xk = round((Xt[k][0]-x0)/ddx) - yk = round((Xt[k][1]-y0)/ddy) - zk = round((Xt[k][2]-z0)/ddz) - SqXt[k,0] = int(xk) - SqXt[k,1] = int(yk) - SqXt[k,2] = int(zk) - if rank==0: - print('BOX cib created: H = ' + str(mesh0.hmin()) + ' H_max = ' + str(mesh0.hmax()) ) - - BOX = {} - BOX['mesh'] = mesh0 - BOX['FEM'] = P1 - BOX['nodes'] = Xt - BOX['seq'] = SqXt.astype(int) - BOX['name'] = 'cib' - BOX['dimension'] = [Lx,Ly,Lz] - return BOX - - if meshtype=='cib_GM': - ddx = 1.8/1000; ddy = 1.8/1000; ddz = 1.8/1000 - x0 = ddx*20 ; y0 = ddy*43; z0 = -ddz*5 - Lx = 110; Ly = 70; Lz = 43 - x1 = x0+ddx*(Lx-1); y1 = y0+ddy*(Ly-1); z1 = z0+ddz*(Lz-1) - - mesh0 = BoxMesh(Point(x0,y0,z0),Point(x1,y1,z1),Lx-1,Ly-1,Lz-1) - P1 = FunctionSpace(mesh0, "Lagrange", 1) - Xt = P1.tabulate_dof_coordinates() - #XDMFFile('mesh_hex.xdmf').write(mesh0) - # MESH HEXA - #mesh0 = BoxMesh.create([Point(x0,y0,z0), Point(x1,y1,z1)], [Lx, Ly, Lz], CellType.Type.hexahedron) - #P1 = FunctionSpace(mesh0, "DG", 0) - #t = P1.tabulate_dof_coordinates() - - if '2017' in dolfin.__version__: - Xt = np.reshape(Xt,( round(Xt.size/3) ,3)) - SqXt = np.zeros(Xt.shape) - - for k in range(Xt.shape[0]): - xk = round((Xt[k][0]-x0)/ddx) - yk = round((Xt[k][1]-y0)/ddy) - zk = round((Xt[k][2]-z0)/ddz) - SqXt[k,0] = int(xk) - SqXt[k,1] = int(yk) - SqXt[k,2] = int(zk) - if rank==0: - print('BOX cib created: H = ' + str(mesh0.hmin()) + ' H_max = ' + str(mesh0.hmax()) ) - - BOX = {} - BOX['mesh'] = mesh0 - BOX['FEM'] = P1 - BOX['nodes'] = Xt - BOX['seq'] = SqXt.astype(int) - BOX['name'] = 'cib' - BOX['dimension'] = [Lx,Ly,Lz] - return BOX - - if meshtype=='aorta': - - if '.xml' in meshpath: - mesh1 = Mesh(meshpath) - Evel1 = VectorElement('Lagrange', mesh1.ufl_cell(), 1) - V1 = FunctionSpace(mesh1, Evel1) - boundaries1 = MeshFunction( - "size_t", mesh1, mesh1.topology().dim()-1) - if '.h5' in meshpath: - mesh1 = Mesh() - hdf1 = HDF5File(mesh1.mpi_comm(), meshpath, 'r') - hdf1.read(mesh1, '/mesh', False) - boundaries1 = MeshFunction( - 'size_t', mesh1, mesh1.topology().dim() - 1) - hdf1.read(boundaries1, '/boundaries') - Evel1 = VectorElement('Lagrange', mesh1.ufl_cell(), 1) - V1 = FunctionSpace(mesh1, Evel1) - - if rank == 0: - print('AORTA created: h_min = ' + str(mesh1.hmin()) + - ' h_max = ' + str(mesh1.hmax())) - - AORTA = {} - AORTA['mesh'] = mesh1 - AORTA['boundaries'] = boundaries1 - AORTA['FEM'] = V1 - AORTA['name'] = 'aorta' - - return AORTA - - if meshtype == 'leo': - - # LEO - mesh = Mesh(meshpath) - Evel = VectorElement('Lagrange', mesh.ufl_cell(), 1) - #Evel = VectorElement('DG',mesh2.ufl_cell(),0) - V = FunctionSpace(mesh, Evel) - boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1) - marked = True - - if marked: - class Outlet(SubDomain): - def inside(self, x, on_boundary): - return on_boundary and between(x[2], (11.18, 11.4)) - - class Inlet(SubDomain): - def inside(self, x, on_boundary): - return on_boundary and between(x[2], (7.56, 7.8)) and between(x[0], (10, 16)) - - outlet = Outlet() - inlet = Inlet() - boundaries.set_all(0) - outlet.mark(boundaries, 3) - inlet.mark(boundaries, 2) - - if rank == 0: - print('LEO created: H = ' + str(mesh.hmin()) + - ' H_max = ' + str(mesh.hmax())) - - LEO = {} - LEO['mesh'] = mesh - LEO['FEM'] = V - LEO['boundaries'] = boundaries - - return LEO - - if meshtype=='other': - - if '.xml' in meshpath: - mesh1 = Mesh(meshpath) - Evel1 = VectorElement('Lagrange', mesh1.ufl_cell(), 1) - V1 = FunctionSpace(mesh1, Evel1) - boundaries1 = MeshFunction( - "size_t", mesh1, mesh1.topology().dim()-1) - if '.h5' in meshpath: - mesh1 = Mesh() - hdf1 = HDF5File(mesh1.mpi_comm(), meshpath, 'r') - hdf1.read(mesh1, '/mesh', False) - boundaries1 = MeshFunction( - 'size_t', mesh1, mesh1.topology().dim() - 1) - hdf1.read(boundaries1, '/boundaries') - Evel1 = VectorElement('Lagrange', mesh1.ufl_cell(), 1) - V1 = FunctionSpace(mesh1, Evel1) - - if rank == 0: - print('MESH created: h_min = ' + str(mesh1.hmin()) + - ' h_max = ' + str(mesh1.hmax())) - - MESH = {} - MESH['mesh'] = mesh1 - MESH['boundaries'] = boundaries1 - MESH['FEM'] = V1 - - return MESH - -def CREATEsequences(BOX,AORTA,options): - # Checkpoints folders - checkpoint_path = options['sequence']['checkpoint_path']+'checkpoint/' - unsort_indexes = os.listdir(checkpoint_path) - indexes = [int(x) for x in unsort_indexes] - indexes.sort() - [Lx, Ly, Lz] = BOX['dimension'] - - Lt = len(indexes) - Sqx = np.zeros([Lx, Ly, Lz, Lt]) - Sqy = np.zeros([Lx, Ly, Lz, Lt]) - Sqz = np.zeros([Lx, Ly, Lz, Lt]) - - v1x = Function(BOX['FEM']) - v1y = Function(BOX['FEM']) - v1z = Function(BOX['FEM']) - v2 = Function(AORTA['FEM']) - N_sp = int(options['sequence']['sampling_rate']) - - if rank==0: - print('Total number of checkpoints folders: ',len(indexes)) - - if options['sequence']['sampling_rate']>1: - if rank==0: - print('Applying an undersampling of ' + str(N_sp)) - print('list of indexes',indexes) - - for k in range(0, len(indexes), N_sp): - if rank == 0: - print('READING CHECKPOINT FOLDER: ' + str(indexes[k])) - - hdf = HDF5File(AORTA['mesh'].mpi_comm( - ), checkpoint_path + str(indexes[k])+'/u.h5', 'r') - hdf.read(v2, 'u/vector_0') - hdf.close() - - LagrangeInterpolator.interpolate(v1x, v2.sub(0)) - LagrangeInterpolator.interpolate(v1y, v2.sub(1)) - LagrangeInterpolator.interpolate(v1z, v2.sub(2)) - - v1_gather = v1x.vector().gather( - np.array(range(BOX['FEM'].dim()), 'intc')) - v2_gather = v1y.vector().gather( - np.array(range(BOX['FEM'].dim()), 'intc')) - v3_gather = v1z.vector().gather( - np.array(range(BOX['FEM'].dim()), 'intc')) - - SqXt = BOX['seq'] - #SqXt = SqXt.astype(int) - # To Gather the Sorting-Sequence in parallel with Gatherv - loc_X = SqXt[:, 0] - loc_Y = SqXt[:, 1] - loc_Z = SqXt[:, 2] - sendbuf_X = np.array(loc_X) - sendbuf_Y = np.array(loc_Y) - sendbuf_Z = np.array(loc_Z) - sendcounts_X = np.array(comm.gather(len(sendbuf_X), 0)) - sendcounts_Y = np.array(comm.gather(len(sendbuf_Y), 0)) - sendcounts_Z = np.array(comm.gather(len(sendbuf_Z), 0)) - - if rank == 0: - SStot_X = np.zeros([Lx*Ly*Lz], dtype=int) - SStot_Y = np.zeros([Lx*Ly*Lz], dtype=int) - SStot_Z = np.zeros([Lx*Ly*Lz], dtype=int) - else: - SStot_X = None - SStot_Y = None - SStot_Z = None - - comm.Gatherv(sendbuf=sendbuf_X, recvbuf=( - SStot_X, sendcounts_X), root=0) - comm.Gatherv(sendbuf=sendbuf_Y, recvbuf=( - SStot_Y, sendcounts_Y), root=0) - comm.Gatherv(sendbuf=sendbuf_Z, recvbuf=( - SStot_Z, sendcounts_Z), root=0) - - if rank == 0: - for l in range(Lx*Ly*Lz): - Sqx[SStot_X[l], SStot_Y[l], SStot_Z[l], k] = v1_gather[l] - Sqy[SStot_X[l], SStot_Y[l], SStot_Z[l], k] = v2_gather[l] - Sqz[SStot_X[l], SStot_Y[l], SStot_Z[l], k] = v3_gather[l] - - if options['sequence']['sampling_rate']>1: - Sqx = Sqx[:, :, :, 0::N_sp] - Sqy = Sqy[:, :, :, 0::N_sp] - Sqz = Sqz[:, :, :, 0::N_sp] - - if rank == 0: - print('saving the sequences' + options['sequence']['checkpoint_path']) - savepath = options['sequence']['checkpoint_path'] + options['sequence']['name'] + '.npz' - np.savez_compressed(savepath, x=Sqx, y=Sqy, z=Sqz) - - - return [Sqx,Sqy,Sqz] - -def LOADsequences(loadpath): - # for loading existing sequences alocated in loadpath - if rank==0: - print('{reading} ' + loadpath) - - if '.mat' in loadpath: - import scipy.io as sio - S = sio.loadmat(loadpath) - S = S['u_R1'] - Sqx = S['x'][0][0] - Sqy = S['y'][0][0] - Sqz = S['z'][0][0] - else: - S = np.load(loadpath) - Sqx = S['x'] - Sqy = S['y'] - Sqz = S['z'] - - return [Sqx,Sqy,Sqz] - -def SqtoH5(BOX,MESH,Sqx,Sqy,Sqz): - - Xt = BOX['nodes'] - SqXt = BOX['seq'] - P1 = BOX['FEM'] - V = MESH['FEM'] - velocity = {} - Lt = Sqx.shape[3] - VX = V.sub(0).collapse() - VY = V.sub(1).collapse() - VZ = V.sub(2).collapse() - vv_x = Function(VX) - vv_y = Function(VY) - vv_z = Function(VZ) - - if rank == 0: - print('{SQtoH5} total number of timesteps: ' + str(Lt)) - - nit = 0 - nit_tot = (Lt-1)/20 - - for t in range(Lt): - if rank == 0: - if PRINT_BARS: - if np.mod(t, nit_tot) < 1: - sys.stdout.write('\r') - # Progress bar - sys.stdout.write( - "{SQtoH5} [%-40s] %d%%" % ('=='*nit, 100*t/Lt)) - sys.stdout.flush() - nit = nit + 1 - else: - print('iteration number:', t) - - v1_x = Function(P1) - v1_y = Function(P1) - v1_z = Function(P1) - - values_x = np.zeros(v1_x.vector().get_local().shape) - values_y = np.zeros(v1_y.vector().get_local().shape) - values_z = np.zeros(v1_z.vector().get_local().shape) - S0x = Sqx[:, :, :, t] - S0y = Sqy[:, :, :, t] - S0z = Sqz[:, :, :, t] - - for k in range(Xt.shape[0]): - values_x[k] = S0x[SqXt[k, 0], SqXt[k, 1], SqXt[k, 2]] - values_y[k] = S0y[SqXt[k, 0], SqXt[k, 1], SqXt[k, 2]] - values_z[k] = S0z[SqXt[k, 0], SqXt[k, 1], SqXt[k, 2]] - - v1_x.vector()[:] = values_x - v1_y.vector()[:] = values_y - v1_z.vector()[:] = values_z - - LagrangeInterpolator.interpolate(vv_x, v1_x) - LagrangeInterpolator.interpolate(vv_y, v1_y) - LagrangeInterpolator.interpolate(vv_z, v1_z) - - v = Function(V) - assign(v.sub(0), vv_x) - assign(v.sub(1), vv_y) - assign(v.sub(2), vv_z) - - velocity[t] = v - - - return velocity - -def MASKED(S_REF,S_RAW): - - [Sq0x,Sq0y,Sq0z] = S_REF - [Sq1x,Sq1y,Sq1z] = S_RAW - - if Sq0z.shape != Sq1z.shape: - print('Original',Sq0z.shape) - print('To mask',Sq1z.shape) - - raise Exception('Dimension of masking sequences not match!') - - # masked Sq1 according to Sq0. Both sequences must have the same sizes. - Sq2x = np.zeros(Sq0x.shape) - Sq2y = np.zeros(Sq0y.shape) - Sq2z = np.zeros(Sq0z.shape) - - nonzerox = np.where(Sq0x!=0) - nonzeroy = np.where(Sq0y!=0) - nonzeroz = np.where(Sq0z!=0) - - Sq2x[nonzerox] = Sq1x[nonzerox] - Sq2y[nonzeroy] = Sq1y[nonzeroy] - Sq2z[nonzeroz] = Sq1z[nonzeroz] - - return [Sq2x,Sq2y,Sq2z] - -def READcheckpoint(MESH,mode,ckpath,under_rate,name,zeros=True): - - V = MESH['FEM'] - W = MESH['FEM'].sub(0).collapse() - # Checkpoints folders - unsort_indexes = os.listdir(ckpath) - indexes = [int(x) for x in unsort_indexes] - indexes.sort() - - numt = 0 - nit = 0 - nit_tot = (len(indexes)-1)/20 - - from dolfin import HDF5File - if mode == 'u': - velocity = {} - for k in range(0, len(indexes), under_rate): - path = ckpath + str(indexes[k]) + '/'+name+'.h5' - - if zeros: - if indexes[k] < 10 and indexes[k] > 0: - path = ckpath + '0' + str(indexes[k]) + '/'+name+'.h5' - - if rank == 0: - if PRINT_BARS: - if np.mod(numt, nit_tot) < 1: - sys.stdout.write('\r') - sys.stdout.write( - "{vel checkpoint} [%-40s] %d%%" % ('=='*nit, 100*numt/(len(indexes)-1))) - sys.stdout.flush() - nit = nit + 1 - - #if options['P_measurement']=='p2': - #v = Function(V) - #vmix = Function(V) - #v,p = vmix.split() - #from dolfin import HDF5File - #hdf = HDF5File(AORTA['mesh'].mpi_comm(), path, 'r') - #hdf.read(v, 'u') - #hdf.close() - #udat = h5py.File(path, 'r') - #v.vector()[:] = udat['u']['vector_0'] - #else: - - v = Function(V) - comm = MPI.COMM_WORLD - hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r') - hdf.read(v, 'u/vector_0') - hdf.close() - velocity[numt] = v - numt = numt + 1 - - return velocity - - if mode == 'p': - pressure = {} - for k in range(0, len(indexes), under_rate): - path = ckpath + str(indexes[k]) + '/'+name+'.h5' - if zeros: - if indexes[k] < 10 and indexes[k] > 0: - path = ckpath + '0' + str(indexes[k]) + '/'+name+'.h5' - - if rank == 0: - if PRINT_BARS: - if np.mod(numt, nit_tot) < 1: - sys.stdout.write('\r') - sys.stdout.write( - "{press checkpoint} [%-40s] %d%%" % ('=='*nit, 100*numt/(len(indexes)-1))) - sys.stdout.flush() - nit = nit + 1 - p = Function(W) - hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r') - hdf.read(p, 'p/vector_0') - hdf.close() - pressure[numt] = p - numt = numt + 1 - - return pressure - -def CHANGEmesh(MESH,field,mode): - # To change between meshes the same result - field2 = {} - if mode=='velocity': - V = MESH['FEM'] - elif mode=='pressure': - V = MESH['FEM'].sub(0).collapse() - else: - raise Exception('You need to assign a mode for the changing!') - for k in range(len(list(field))): - if rank==0: - print('CHANGING MESH: index',k) - v2 = Function(V) - LagrangeInterpolator.interpolate(v2,field[k]) - field2[k] = v2 - - - return field2 - -def GradP(MESH,pressure,options): - - W = MESH['FEM'].sub(0).collapse() - boundaries = MESH['boundaries'] - barye2mmHg = 1/1333.22387415 - - Tf = Constant(options['Tf']) - dt = Constant(options['Tf']/len(list(pressure))) - NT = int(float(Tf)/float(dt)) - - ind = 0 - t = float(dt) - ds = Measure("ds", subdomain_data=boundaries) - ones = interpolate(Constant(1),W) - DELTA_p = np.zeros([NT]) - - p0 = Function(W) - nit = 0 - nit_tot = (NT-1)/20 - - zero_point = None - if options['reference']['fix_const']: - zero_point = options['reference']['zero_point'] - zero_point = np.array(zero_point) - ndim = W.mesh().topology().dim() - zero_point = W.tabulate_dof_coordinates().reshape((-1, ndim))[0] - - #################### THE TIME LOOP ######################### - while ind < NT: - p0.assign(pressure[ind]) - - if rank == 0: - if PRINT_BARS: - if np.mod(ind, nit_tot) < 1: - sys.stdout.write('\r') - # Progress bar - sys.stdout.write("{ref} [%-40s] %f seg" % - ('=='*nit, round(t, 2))) - sys.stdout.flush() - nit = nit + 1 - - if options['reference']['fix_const']: - p0.vector()[:] -= p0(zero_point) - - i_pin = assemble(p0*ds(2)) - i_pout = assemble(p0*ds(3)) - i_s2 = assemble(ones*ds(2)) - i_s3 = assemble(ones*ds(3)) - - DELTA_p[ind] = barye2mmHg*(i_pin/i_s2-i_pout/i_s3) - ind = ind + 1 - t = t + float(dt) - - if options['reference']['save']: - - if rank==0: - np.savetxt(options['reference']['savepath'],DELTA_p) - print('\n Pressure drop save in ' + options['reference']['savepath']) - -def Fluxes(MESH,velocity,options,bnd): - - mesh = MESH['mesh'] - W = MESH['FEM'] - boundaries = MESH['boundaries'] - n = FacetNormal(mesh) - Tf = Constant(options['Tf']) - rho = Constant(options['density']) - dt = Constant(options['Tf']/len(list(velocity))) - NT = int(float(Tf)/float(dt)) - ind = 0 - t = float(dt) - ds = Measure("ds", subdomain_data=boundaries) - QQ = {} - for k in bnd: - QQ[k] = np.zeros([NT]) - - u = Function(W) - nit = 0 - nit_tot = (NT-1)/20 - #################### THE TIME LOOP ######################### - while indVENC) - INDy = np.any(np.abs(vy)>VENC) - INDz = np.any(np.abs(vz)>VENC) - - if INDx or INDy or INDz: - if not options['phase_contrast']['dealiased']: - if rank==0: - raise Exception('There is some velocity greater or lower than VENC. You should turn dealised true!') - - - if options['phase_contrast']['dealiased']: - [vx_d,vy_d,vz_d] = dealiased(VENC,vx,vy,vz) - else: - vx_d = vx; vy_d = vy; vz_d = vz - - if rank==0: - print('saving velocity as ' + options['phase_contrast']['output']) - np.savez_compressed(options['phase_contrast']['output'], x=vx_d, y=vy_d, z=vz_d) - - if 'resize' in options: - if options['resize']['apply']: - if rank == 0: - print('Changing the size of the sequence') - - [Nx, Ny, Nz] = options['resize']['dim'] - [Sqx, Sqy, Sqz] = LOADsequences(options['resize']['loadseq']) - [Nx0, Ny0, Nz0, Nt] = Sqx.shape - - Sqx_new = np.zeros([Nx, Ny, Nz, Nt]) - Sqy_new = np.zeros([Nx, Ny, Nz, Nt]) - Sqz_new = np.zeros([Nx, Ny, Nz, Nt]) - - ax = int((Nx - Nx0)/2) - bx = ax + Nx0 - ay = int((Ny - Ny0)/2) - by = ay + Ny0 - az = int((Nz - Nz0)/2) - bz = az + Nz0 - - for t in range(Nt): - Sqx_new[ax:bx, ay:by, az:bz, t] = Sqx[:, :, :, t] - Sqy_new[ax:bx, ay:by, az:bz, t] = Sqy[:, :, :, t] - Sqz_new[ax:bx, ay:by, az:bz, t] = Sqz[:, :, :, t] - - if options['resize']['save']: - if rank == 0: - print('saving the sequences' + - options['resize']['savepath']) - np.savez_compressed( - options['resize']['savepath'], x=Sqx_new, y=Sqy_new, z=Sqz_new) - - if 'magnetization_CIB' in options: - if options['magnetization_CIB']['apply']: - if rank == 0: - print('-- Generating k-space for CIB phantom data-- ') - print('Using 0-G gradients') - - import Graphics - masterpath = options['magnetization_CIB']['loadpath'] - mask = Graphics.MATmedical(masterpath+'Segmented_Aorta.mat') - mask = mask.transpose((2, 0, 1)) - VENC = Graphics.MATmedical(masterpath+'VENC.mat') - VENC = VENC[0][0] - if rank == 0: - print('Recognized VENC = ' + str(VENC) + ' cm/s') - - # Magnitudes - magnitude_AP = Graphics.MATmedical(masterpath+'MR_FFE_AP.mat') - magnitude_FH = Graphics.MATmedical(masterpath+'MR_FFE_FH.mat') - magnitude_RL = Graphics.MATmedical(masterpath+'MR_FFE_RL.mat') - # Velocities - vel_AP = Graphics.MATmedical(masterpath+'MR_PCA_AP.mat') - vel_FH = Graphics.MATmedical(masterpath+'MR_PCA_FH.mat') - vel_RL = Graphics.MATmedical(masterpath+'MR_PCA_RL.mat') - - # Cropping images - a_crop = 48 - b_crop = 178 - mask = mask[:, :, a_crop:b_crop] - magnitude_AP = magnitude_AP[:, :, a_crop:b_crop, :] - vel_AP = vel_AP[:, :, a_crop:b_crop, :] - magnitude_RL = magnitude_RL[:, :, a_crop:b_crop, :] - vel_RL = vel_RL[:, :, a_crop:b_crop, :] - magnitude_FH = magnitude_FH[:, :, a_crop:b_crop, :] - vel_FH = vel_FH[:, :, a_crop:b_crop, :] - - # Reference phase - [Nx, Ny, Nz, Nt] = vel_AP.shape - gamma = 267.513e6 # rad/Tesla/sec Gyromagnetic ratio for H nuclei - B0 = 1.5 # Tesla Magnetic Field Strenght - TE = 5e-3 # Echo-time - [X, Y, Z] = np.meshgrid(np.linspace( - 0, Ny, Ny), np.linspace(0, Nx, Nx), np.linspace(0, Nz, Nz)) - frecX = 0.21 - p_ref = np.zeros([Nx, Ny, Nz]) - p_flow_AP = np.zeros([Nx, Ny, Nz]) - p_flow_RL = np.zeros([Nx, Ny, Nz]) - p_flow_FH = np.zeros([Nx, Ny, Nz]) - #Sx2G = np.zeros([Nx,Ny,Nz,Nt],dtype=complex) - #Sy2G = np.zeros([Nx,Ny,Nz,Nt],dtype=complex) - #Sz2G = np.zeros([Nx,Ny,Nz,Nt],dtype=complex) - Sx20 = np.zeros([Nx, Ny, Nz, Nt], dtype=complex) - Sy20 = np.zeros([Nx, Ny, Nz, Nt], dtype=complex) - Sz20 = np.zeros([Nx, Ny, Nz, Nt], dtype=complex) - - treshold = np.mean(magnitude_AP) - #p_ref = (gamma*B0*TE+frecX*X)*(np.abs(magnitude_AP[:,:,:,7])>treshold) - p_ref = (gamma*B0*TE+frecX*X)*(mask > 0.5) - - for t in range(Nt): - p_flow_AP = p_ref + np.pi*vel_AP[:, :, :, t]/VENC - p_flow_RL = p_ref + np.pi*vel_RL[:, :, :, t]/VENC - p_flow_FH = p_ref + np.pi*vel_FH[:, :, :, t]/VENC - - Sx20[:, :, :, t] = magnitude_AP[:, :, :, t] * \ - np.cos(p_ref) + 1j*magnitude_AP[:, :, :, t]*np.sin(p_ref) - Sy20[:, :, :, t] = magnitude_RL[:, :, :, t] * \ - np.cos(p_ref) + 1j*magnitude_RL[:, :, :, t]*np.sin(p_ref) - Sz20[:, :, :, t] = magnitude_FH[:, :, :, t] * \ - np.cos(p_ref) + 1j*magnitude_FH[:, :, :, t]*np.sin(p_ref) - - #Sx2G[:,:,:,t] = magnitude_AP[:,:,:,t]*np.cos(p_flow_AP) + 1j*magnitude_AP[:,:,:,t]*np.sin(p_flow_AP) - #Sy2G[:,:,:,t] = magnitude_RL[:,:,:,t]*np.cos(p_flow_RL) + 1j*magnitude_RL[:,:,:,t]*np.sin(p_flow_RL) - #Sz2G[:,:,:,t] = magnitude_FH[:,:,:,t]*np.cos(p_flow_FH) + 1j*magnitude_FH[:,:,:,t]*np.sin(p_flow_FH) - - if rank == 0: - print('saving the sequences' + - options['magnetization_CIB']['savepath']) - savepath = options['magnetization_CIB']['savepath'] - basepath = savepath.replace('.npz', '') - savepath0 = basepath + '0.npz' - savepathG = basepath + 'G.npz' - np.savez_compressed(savepath0, x=Sx20, y=Sy20, z=Sz20) - #np.savez_compressed(savepathG, x=Sx2G,y=Sy2G,z=Sz2G) - - if 'magnetization' in options: - if options['magnetization']['apply']: - if rank==0: - print('-- Generating k-space -- ') - print('Using 0-G gradients') - - [Sqx,Sqy,Sqz] = LOADsequences(options['magnetization']['loadseq']) - [Nx,Ny,Nz,Nt] = Sqx.shape - VENC = options['magnetization']['VENC'] - gamma = 267.513e6 # rad/Tesla/sec Gyromagnetic ratio for H nuclei - B0 = 1.5 # Tesla Magnetic Field Strenght - TE = 5e-3 # Echo-time - [X,Y,Z] = np.meshgrid(np.linspace(0,Ny,Ny),np.linspace(0,Nx,Nx),np.linspace(0,Nz,Nz)) - frecX = 0.21 # 0.21 without leaps: 0.01 - - Mx = np.zeros(Sqx.shape) - My = np.zeros(Sqx.shape) - Mz = np.zeros(Sqx.shape) - Sx2G = np.zeros(Sqx.shape,dtype=complex) - Sy2G = np.zeros(Sqx.shape,dtype=complex) - Sz2G = np.zeros(Sqx.shape,dtype=complex) - Sx20 = np.zeros(Sqx.shape,dtype=complex) - Sy20 = np.zeros(Sqx.shape,dtype=complex) - Sz20 = np.zeros(Sqx.shape,dtype=complex) - Fx0 = np.zeros(Sqx.shape) - Fy0 = np.zeros(Sqx.shape) - Fz0 = np.zeros(Sqx.shape) - FxG = np.zeros(Sqx.shape) - FyG = np.zeros(Sqx.shape) - FzG = np.zeros(Sqx.shape) - - # Creating M - for t in range(Nt): - Mx[:,:,:,t] = 0.05 + 0.9*(np.abs(Sqx[:,:,:,t])>0.001) + 0.05*np.sqrt(np.abs(Sqx[:,:,:,t])) - My[:,:,:,t] = 0.05 + 0.9*(np.abs(Sqy[:,:,:,t])>0.001) + 0.05*np.sqrt(np.abs(Sqy[:,:,:,t])) - Mz[:,:,:,t] = 0.05 + 0.9*(np.abs(Sqz[:,:,:,t])>0.001) + 0.05*np.sqrt(np.abs(Sqz[:,:,:,t])) - # Organs - Mx[:,:,:,t] = Mx[:,:,:,t] + 1*( (X-Ny/2)**2 + 0.75*(Y-Nx/2)**2 < 35**2 ) * ( (X-Ny/2)**2 + 0.75*(Y-Nx/2)**2 > 33**2 ) - Mx[:,:,:,t] = Mx[:,:,:,t] + 0.7*( 0.3*(X-Ny/2+9)**2 + 0.7*(Y-Nx/2-4)**2 < 3**2 )* ( np.abs(Z-Nz/2)<25 ) - My[:,:,:,t] = My[:,:,:,t] + 1*( (X-Ny/2)**2 + 0.75*(Y-Nx/2)**2 < 35**2 ) * ( (X-Ny/2)**2 + 0.75*(Y-Nx/2)**2 > 33**2 ) - My[:,:,:,t] = My[:,:,:,t] + 0.7*( 0.3*(X-Ny/2+9)**2 + 0.7*(Y-Nx/2-4)**2 < 3**2 )* ( np.abs(Z-Nz/2)<25 ) - Mz[:,:,:,t] = Mz[:,:,:,t] + 1*( (X-Ny/2)**2 + 0.75*(Y-Nx/2)**2 < 35**2 ) * ( (X-Ny/2)**2 + 0.75*(Y-Nx/2)**2 > 33**2 ) - Mz[:,:,:,t] = Mz[:,:,:,t] + 0.7*( 0.3*(X-Ny/2+9)**2 + 0.7*(Y-Nx/2-4)**2 < 3**2 )* ( np.abs(Z-Nz/2)<25 ) - FxG[:,:,:,t] = (gamma*B0*TE+frecX*X)*(np.abs(Mz[:,:,:,t])>0.4) + np.pi*Sqx[:,:,:,t]/VENC - FyG[:,:,:,t] = (gamma*B0*TE+frecX*X)*(np.abs(Mz[:,:,:,t])>0.4) + np.pi*Sqy[:,:,:,t]/VENC - FzG[:,:,:,t] = (gamma*B0*TE+frecX*X)*(np.abs(Mz[:,:,:,t])>0.4) + np.pi*Sqz[:,:,:,t]/VENC - Fx0[:,:,:,t] = (gamma*B0*TE+frecX*X)*(np.abs(Mz[:,:,:,t])>0.4) - Fy0[:,:,:,t] = (gamma*B0*TE+frecX*X)*(np.abs(Mz[:,:,:,t])>0.4) - Fz0[:,:,:,t] = (gamma*B0*TE+frecX*X)*(np.abs(Mz[:,:,:,t])>0.4) - - Sx20[:,:,:,t] = Mx[:,:,:,t]*np.cos(Fx0[:,:,:,t]) + 1j*Mx[:,:,:,t]*np.sin(Fx0[:,:,:,t]) - Sy20[:,:,:,t] = My[:,:,:,t]*np.cos(Fy0[:,:,:,t]) + 1j*My[:,:,:,t]*np.sin(Fy0[:,:,:,t]) - Sz20[:,:,:,t] = Mz[:,:,:,t]*np.cos(Fz0[:,:,:,t]) + 1j*Mz[:,:,:,t]*np.sin(Fz0[:,:,:,t]) - Sx2G[:,:,:,t] = Mx[:,:,:,t]*np.cos(FxG[:,:,:,t]) + 1j*Mx[:,:,:,t]*np.sin(FxG[:,:,:,t]) - Sy2G[:,:,:,t] = My[:,:,:,t]*np.cos(FyG[:,:,:,t]) + 1j*My[:,:,:,t]*np.sin(FyG[:,:,:,t]) - Sz2G[:,:,:,t] = Mz[:,:,:,t]*np.cos(FzG[:,:,:,t]) + 1j*Mz[:,:,:,t]*np.sin(FzG[:,:,:,t]) - - - # Adding noise by cartesian components - for t in range(Nt): - noise1 = np.random.normal(-0.0007,0.049,[Nx,Ny,Nz]) - noise2 = np.random.normal(-0.0007,0.049,[Nx,Ny,Nz]) - Sx20[:,:,:,t] += noise1 + 1j*noise2 - Sy20[:,:,:,t] += noise1 + 1j*noise2 - Sz20[:,:,:,t] += noise1 + 1j*noise2 - - noise1g = np.random.normal(-0.0007,0.049,[Nx,Ny,Nz]) - noise2g = np.random.normal(-0.0007,0.049,[Nx,Ny,Nz]) - Sx2G[:,:,:,t] += noise1g + 1j*noise2g - Sy2G[:,:,:,t] += noise1g + 1j*noise2g - Sz2G[:,:,:,t] += noise1g + 1j*noise2g - - - if rank==0: - print('saving the sequences' + options['magnetization']['savepath']) - savepath = options['magnetization']['savepath'] - basepath = savepath.replace('.npz','') - savepath0 = basepath + '0.npz' - savepathG = basepath + 'G.npz' - np.savez_compressed(savepath0, x=Sx20,y=Sy20,z=Sz20) - np.savez_compressed(savepathG, x=Sx2G,y=Sy2G,z=Sz2G) - - if 'sequence' in options: - if options['sequence']['apply']: - if rank == 0: - print('--- Creating Sequences ---') - - resol = options['sequence']['resol'] - boxsize = options['sequence']['boxsize'] - ranges = options['sequence']['ranges'] - BOX = LOADmesh(options['sequence']['boxtype'], resol=resol, boxsize=boxsize, ranges=ranges) - meshpath = options['sequence']['meshpath'] - AORTA = LOADmesh('aorta', meshpath=meshpath) - [Sqx, Sqy, Sqz] = CREATEsequences(BOX, AORTA, options) - - if 'norms' in options: - if options['norms']['apply']: - [BOX,AORTA,LEO] = CREATEmeshes(options) - if rank==0: - print('Computing Norms between fields') - ct_norm = checkpoint_norm(AORTA,options['norms']['field_path'],'ct') - kal_norm = checkpoint_norm(AORTA,options['norms']['field_path'],'roukf') - - if options['norms']['plot']: - import matplotlib.pyplot as plt - plt.plot(ct_norm, color='blue',marker='.',label='ct') - #plt.plot(meas_norm,label='meas') - plt.plot(kal_norm,color='green',label='kalman') - #plt.ylim([0,8500]) - plt.legend(fontsize=14) - plt.show() - - if 'CIBtoH5' in options: - print('---Converting CIB files into H5---') - if options['CIBtoH5']['apply']: - path_to_cib = options['CIBtoH5']['data_path'] - outpath = options['CIBtoH5']['outpath'] - times = options['CIBtoH5']['times'] - dt = options['CIBtoH5']['dt'] - flip = options['CIBtoH5']['flip'] - - if options['CIBtoH5']['interpolate']: - mesh_path = options['CIBtoH5']['mesh_path'] - CIBtoH5(path_to_cib,times,dt,outpath,interpolate=True,mesh_path=mesh_path,flip=flip) - else: - CIBtoH5(path_to_cib,times,dt,outpath,flip=flip) - - if 'cs' in options: - if options['cs']['apply']: - if rank==0: - print('Applying Compressed Sensing') - - [Sqx, Sqy, Sqz] = LOADsequences(options['cs']['seqpath']) - - import CS - if 'short' in options['cs']: - if options['cs']['short']: - [Mx,My,Mz] = LOADsequences(options['cs']['Mpath']) - CS.undersampling_short(Mx,My,Mz,options) - else: - CS.undersampling(Sqx,Sqy,Sqz,options,options['cs']['savepath']) - - if 'kt-BLAST' in options: - if options['kt-BLAST']['apply']: - if rank==0: - print('Applying kt-BLAST') - import ktBLAST - ktBLAST.undersampling(Sqx,Sqy,Sqz,options,options['kt-BLAST']['savepath']) - - if 'SENSE' in options: - if options['SENSE']['apply']: - if rank==0: - print('-- SENSE reconstruction --') - import SENSE - - [Mx,My,Mz] = LOADsequences(options['SENSE']['Mpath']) - [MxS,MyS,MzS] = SENSE.undersampling(Mx,My,Mz,options) - - if rank==0: - print('saving the sequences' + options['SENSE']['savepath']) - np.savez_compressed(options['SENSE']['savepath'], x=MxS,y=MyS,z=MzS) - - if 'create_checkpoint' in options: - if options['create_checkpoint']['apply']: - print('--- Create Checkpoint ---') - - - mesh_into = options['create_checkpoint']['mesh_into'] - - - if options['create_checkpoint']['extension'] == '.mat': - import scipy.io as sio - seqpath = options['create_checkpoint']['loadseq'] - seqpath = seqpath + options['create_checkpoint']['seqname'] - seqpath = seqpath + '_R1.mat' - S = sio.loadmat(seqpath) - [ddx,ddy,ddz] = S['u_R1']['VoxelSize'][0][0][0] - [ddx,ddy,ddz] = [ddx/10,ddy/10,ddz/10] #mm to cm - S = S['u_R1'] - Sqx = S['x'][0][0] - [Lx,Ly,Lz,Lt] = Sqx.shape - BOX = LOADmesh('box_zeros', resol=[ddx,ddy,ddz], boxsize=[Lx,Ly,Lz]) - del S, Sqx, seqpath - else: - # Box mesh definition - boxsize = options['create_checkpoint']['boxsize'] - resol = options['create_checkpoint']['resol'] - ranges = options['create_checkpoint']['ranges'] - boxtype = options['create_checkpoint']['boxtype'] - if boxtype == 'fix_resolution': - BOX = LOADmesh(boxtype, resol=resol, boxsize=boxsize, ranges=ranges) - - - MESH = LOADmesh('other', meshpath=mesh_into) - - - for r in options['create_checkpoint']['Rseq']: - if rank == 0: - print('Writing checkpoint from sequence ' + - options['create_checkpoint']['seqname'] + '_R' + str(r) + '.npz') - - seqpath = options['create_checkpoint']['loadseq'] + options['create_checkpoint']['seqname'] + \ - '_R' + str(r) + options['create_checkpoint']['extension'] - - [Sx, Sy, Sz] = LOADsequences(seqpath) - - if options['create_checkpoint']['under_rate']>1: - print('Undersampling with factor: ', - options['create_checkpoint']['under_rate']) - Sx = Sx[:, :, :, 0::options['create_checkpoint']['under_rate']] - Sy = Sy[:, :, :, 0::options['create_checkpoint']['under_rate']] - Sz = Sz[:, :, :, 0::options['create_checkpoint']['under_rate']] - - if options['create_checkpoint']['masked']: - [S0x, S0y, S0z] = LOADsequences( - options['create_checkpoint']['masked_seq']) - [Sx, Sy, Sz] = MASKED([S0x, S0y, S0z], [Sx, Sy, Sz]) - del S0x, S0y, S0z - - vel_seq = SqtoH5(BOX, MESH, -Sx, Sy, -Sz) - #vel_seq = SqtoH5(BOX, MESH, Sx, Sy, Sz) - - if rank == 0: - print(' ') - - comm = MESH['mesh'].mpi_comm() - dt = options['create_checkpoint']['dt'] - if options['create_checkpoint']['xdmf']: - xdmf_u = XDMFFile(options['create_checkpoint']['savepath']+ 'R' + str(r) + '/u.xdmf') - - for l in range(len(vel_seq)): - if rank == 0: - print('saving checkpoint', l) - - path = options['create_checkpoint']['savepath'] + \ - 'R' + str(r) + '/checkpoint/{i}/'.format(i=l) - - if options['create_checkpoint']['xdmf']: - vel_seq[l].rename('velocity', 'u') - xdmf_u.write(vel_seq[l], l*dt) - - inout.write_HDF5_data( - comm, path + '/u.h5', vel_seq[l], '/u', t=l*dt) - - if 'reference' in options: - if options['reference']['apply']: - if rank == 0: - print('---Applying Reference Pressure Gradient---') - - ckpath = options['reference']['checkpoint_path'] - under_rate = options['reference']['under_rate'] - MESH = LOADmesh('other', meshpath=options['reference']['meshpath']) - zeros = not 'ct' in ckpath - pressure = READcheckpoint( - MESH, 'p', ckpath, under_rate, 'p', zeros=zeros) - GradP(MESH, pressure, options) - - if 'ppe' in options: - - if options['ppe']['apply']: - [BOX,AORTA,LEO] = CREATEmeshes(options) - if rank==0: - print('Applying PPE') - if options['ppe']['mesh_type']=='aorta': - velocity_model = READcheckpoint(AORTA,'v',options,options['checkpoint_path'],'u') - if rank==0: - print(' ') - PPE(AORTA,velocity_model,options,options['ppe']['name']+ '.txt') - if rank==0: - print(' ') - - if options['ppe']['mesh_type']=='leo': - if not options['cs']['apply'] and not options['ppe']['read_check']: - raise Exception('You need first apply a CS sequence') - if options['ppe']['read_check']: - if options['meshname']=='leo': - velocity = READcheckpoint(LEO,'v',options,options['checkpoint_path'],'u') - MESH = LEO - if rank==0: - print(' ') - PPE(MESH,velocity,options,options['ppe']['name']+ '.txt') - if rank==0: - print(' ') - - - if not options['ppe']['read_check']: - for r in options['cs']['R']: - velocity[r] = SqtoH5(BOX,LEO,Sqx_cs[r],Sqy_cs[r],Sqz_cs[r]) - if rank==0: - print(' ') - PPE(LEO,velocity[r],options,options['ppe']['name']+options['ppe']['under_type']+'_R'+str(r)+'.txt') - if rank==0: - print(' ') - - if options['ppe']['mesh_type']=='leoct': - velocity_model = READcheckpoint(AORTA,'v',options,options['checkpoint_path'],'u') - if rank==0: - print(' ') - velocity_model_leo = CHANGEvel(LEO,velocity_model) - if rank==0: - print(' ') - PPE(LEO,velocity_model_leo,options,options['ppe']['name']+ '.txt') - if rank==0: - print(' ') - - if 'ste' in options: - - if options['ste']['apply']: - [BOX,AORTA,LEO] = CREATEmeshes(options) - if rank==0: - print('Applying STE') - - - if options['ste']['mesh_type']=='aorta': - if velocity_model: - STE(AORTA,velocity_model,options,options['ste']['name']+ '.txt') - else: - velocity_model = READcheckpoint(AORTA,'v',options,options['checkpoint_path'],'u') - - if rank==0: - print('') - - - if options['ste']['int']: - STEint(AORTA,velocity_model,options,options['ste']['name']+ '.txt') - else: - STE(AORTA,velocity_model,options,options['ste']['name']+ '.txt') - - if options['ste']['mesh_type']=='leo': - - if not options['cs']['apply'] and not options['ste']['read_check']: - raise Exception('You need first apply a CS sequence') - if options['ste']['read_check']: - if options['meshname']=='leo': - velocity = READcheckpoint(LEO,'v',options,options['checkpoint_path'],'u') - MESH = LEO - if rank==0: - print(' ') - - if options['ste']['int']: - STEint(MESH,velocity,options,options['ste']['name']+ '.txt') - else: - STE(MESH,velocity,options,options['ste']['name']+ '.txt') - if rank==0: - print(' ') - - if not options['ste']['read_check']: - for r in options['cs']['R']: - if velocity[r]: - STE(LEO,velocity[r],options,options['ste']['name']+options['ste']['under_type']+'_R'+str(r)+'.txt') - if rank==0: - print(' ') - else: - velocity[r] = SqtoH5(BOX,LEO,Sqx_cs[r],Sqy_cs[r],Sqz_cs[r]) - if rank==0: - print(' ') - if options['ste']['int']: - STEint(LEO,velocity[r],options,options['ste']['name']+options['ste']['under_type']+'_R'+str(r)+'.txt') - else: - STE(LEO,velocity[r],options,options['ste']['name']+options['ste']['under_type']+'_R'+str(r)+'.txt') - if rank==0: - print(' ') - - if options['ste']['mesh_type']=='leoct': - if velocity_model_leo: - STE(LEO,velocity_model_leo,options,options['ste']['name']+ '.txt') - if rank==0: - print(' ') - else: - velocity_model = READcheckpoint(AORTA,'v',options,options['checkpoint_path'],'u') - if rank==0: - print(' ') - velocity_model_leo = CHANGEvel(LEO,velocity_model) - if rank==0: - print(' ') - if options['ste']['int']: - STEint(LEO,velocity_model_leo,options,options['ste']['name']+ '.txt') - else: - STE(LEO,velocity_model_leo,options,options['ste']['name']+ '.txt') - if rank==0: - print(' ') - - if 'change_mesh' in options: - if options['change_mesh']['apply']: - if rank == 0: - print('Changing '+options['change_mesh']['mode'] + ' from: ') - print('Mesh input: ' + options['change_mesh']['mesh_in']) - print('Mesh output: ' + options['change_mesh']['mesh_out']) - - MESH_in = LOADmesh( - 'other', meshpath=options['change_mesh']['mesh_in']) - MESH_out = LOADmesh( - 'other', meshpath=options['change_mesh']['mesh_out']) - - ckpath = options['change_mesh']['checkpoint_path'] - under_rate = options['change_mesh']['under_rate'] - mode = options['change_mesh']['mode'] - - origin = READcheckpoint(MESH_in, mode, ckpath, under_rate, mode, zeros=False) - - # To change between meshes the same result - changed = {} - #W = LEO['FEM'].sub(0).collapse() - comm = MESH_out['mesh'].mpi_comm() - - for k in range(len(list(origin))): - v2 = Function(MESH_out['FEM']) - if rank == 0: - print('CHANGING: index', k) - LagrangeInterpolator.interpolate(v2, origin[k]) - changed[k] = v2 - print(max(origin[k].vector().get_local()),max(v2.vector().get_local())) - - - dt = options['change_mesh']['dt'] - if options['create_checkpoint']['xdmf']: - xdmf_u = XDMFFile(options['change_mesh']['savepath']+'u.xdmf') - for l in range(len(changed)): - if rank == 0: - print('saving checkpoint', l) - path = options['change_mesh']['savepath'] + \ - 'checkpoint/{i}/'.format(i=l) - writepath = path + '/'+options['change_mesh']['mode']+'.h5' - inout.write_HDF5_data( - comm, writepath, changed[l] ,'/'+options['change_mesh']['mode'], t=l*dt) - - if options['create_checkpoint']['xdmf']: - changed[l].rename('velocity', 'u') - xdmf_u.write(changed[l], l*dt) - - - -if __name__ == '__main__': - - comm = MPI.COMM_WORLD - size = comm.Get_size() - rank = comm.Get_rank() - - if len(sys.argv) > 1: - if os.path.exists(sys.argv[1]): - inputfile = sys.argv[1] - if rank==0: - print('Found input file ' + inputfile) - else: - Etcetera('serious') - raise Exception('Command line arg given but input file does not exist:' - ' {}'.format(sys.argv[1])) - else: - Etcetera('serious') - raise Exception('An input file is required as argument!') - - - user = 'p283370' # Default's user - if 'Zion' in os.popen('hostname').read(): - user = 'yeye' - np.set_printoptions(threshold=5) - - if 'fwn-bborg-5-166' in os.popen('hostname').read(): - user = 'p283370' - - if rank==0: - print('Welcome user {uss}'.format(uss=user)) - Etcetera('happy') - - - start_time = time.time() - - options = inout.read_parameters(inputfile) - SCANNER(options) - - end_time = time.time() - CLOCK(rank,start_time,end_time) - diff --git a/codes/PostCheck.py b/codes/PostCheck.py deleted file mode 100644 index 0187b37..0000000 --- a/codes/PostCheck.py +++ /dev/null @@ -1,1452 +0,0 @@ -from dolfin import * -import numpy as np -import sys -import os -from mpi4py import MPI -from common import inout - - -def LOADmesh(pathtomesh): - - mesh = Mesh() - hdf = HDF5File(mesh.mpi_comm(), pathtomesh, 'r') - hdf.read(mesh, '/mesh', False) - boundaries = MeshFunction('size_t', mesh, mesh.topology().dim() - 1) - hdf.read(boundaries, '/boundaries') - - #Evel = VectorElement('P',mesh.ufl_cell(),1) - #V = FunctionSpace(mesh,Evel) - - #V1 = VectorElement('P', mesh.ufl_cell(), 1) - #V = FunctionSpace(mesh, V1) - P1 = FiniteElement('P', mesh.ufl_cell(), 1) - V1 = VectorElement('P', mesh.ufl_cell(), 1) - V = FunctionSpace(mesh, V1*P1) - - if rank == 0: - print('MESH ' + ' created: H = ' + - str(mesh.hmin()) + ' H_max = ' + str(mesh.hmax())) - - MESH = {} - MESH['mesh'] = mesh - MESH['FEM'] = V - MESH['boundaries'] = boundaries - - return MESH - -def LOADsequences(loadpath): - # for loading existing sequences alocated in loadpath - if rank == 0: - print('{reading} ' + loadpath) - S = np.load(loadpath) - Sqx = S['x'] - Sqy = S['y'] - Sqz = S['z'] - return [Sqx, Sqy, Sqz] - -def WORKcheck(MESH, mode, output_path, checkpoint_path, filename, outname, options): - - #from dolfin import HDF5File - V = MESH['FEM'] - W = MESH['FEM'].sub(0).collapse() - # Checkpoints folders - unsort_indexes = os.listdir(checkpoint_path) - indexes = [int(x) for x in unsort_indexes] - indexes.sort() - - - if mode == 'u': - v = Function(W) - xdmf_u = XDMFFile(output_path+outname+'.xdmf') - dt = options['Velocity']['dt'] - for k in range(0,len(indexes),options['Velocity']['undersampling']): - te = k*dt - path = checkpoint_path + str(indexes[k]) + '/'+filename+'.h5' - v.rename('velocity', outname) - hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r') - if 'w' in filename: - hdf.read(v, 'w/vector_0') - else: - hdf.read(v, 'u/vector_0') - hdf.close() - xdmf_u.write(v, te) - - if mode == 'w': - xdmf_u = XDMFFile(output_path+outname+'.xdmf') - for k in range(0, len(indexes), 1): - path = checkpoint_path + str(indexes[k]) + '/'+filename+'.h5' - v = Function(V) - v.rename('waux', outname) - comm = MPI.COMM_WORLD - hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r') - hdf.read(v, 'w/vector_0') - hdf.close() - wvec = v.vector().get_local() - wnorm = wvec - np.mean(wvec) - v.vector()[:] = wnorm - xdmf_u.write(v, te) - te = te + dt - numt = numt + 1 - - if mode == 'gradw': - xdmf_u = XDMFFile(output_path+outname+'.xdmf') - for k in range(0, len(indexes), 1): - path = checkpoint_path + str(indexes[k]) + '/'+filename+'.h5' - v = Function(V) - gw = Function(W) - gw.rename(outname, outname) - comm = MPI.COMM_WORLD - hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r') - hdf.read(v, 'w/vector_0') - hdf.close() - wvec = v.vector().get_local() - wnorm = wvec - np.mean(wvec) - v.vector()[:] = wnorm - gw.assign(project(sqrt(inner(grad(v), grad(v))), W)) - xdmf_u.write(gw, te) - te = te + dt - numt = numt + 1 - - if mode == 'p' or mode == 'p_cib': - xdmf_p = XDMFFile(output_path+outname+'.xdmf') - dt = options['Pressure']['dt'] - for k in range(0, len(indexes), options['Pressure']['undersampling']): - te = k*dt - path = checkpoint_path + str(indexes[k]) + '/'+filename+'.h5' - p = Function(W) - if mode == 'p': - barye2mmHg = 1/1333.22387415 - if mode == 'p_cib': - barye2mmHg = -0.00750062 - p.rename('pressure', outname) - hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r') - hdf.read(p, 'p/vector_0') - hdf.close() - p1vec = p.vector().get_local() - #p1vec = p1vec - np.mean(p1vec) - p.vector()[:] = p1vec*barye2mmHg - xdmf_p.write(p, te) - -def READcheckpoint(MESH, mode, output_path, checkpoint_path, filename, outname, options, flow=False, bnds=None): - - #from dolfin import HDF5File - V = MESH['FEM'] - W = MESH['FEM'].sub(0).collapse() - # Checkpoints folders - unsort_indexes = os.listdir(checkpoint_path) - indexes = [int(x) for x in unsort_indexes] - indexes.sort() - - if flow: - QQ = {} - ds = Measure('ds', domain=MESH['mesh'], - subdomain_data=MESH['boundaries']) - n = FacetNormal(MESH['mesh']) - for bb in bnds: - QQ[bb] = [] - - if mode == 'Histograms': - u = Function(W) - if 'w_COR' in filename: - h5name = 'w/vector_0' - else: - h5name = 'u/vector_0' - - - ValuesTime = np.zeros([len(indexes),len(u.vector().get_local())]) - ValuesPeak = np.zeros([len(u.vector().get_local())]) - ix = int(0) - - for k in indexes: - path = checkpoint_path + str(k) + '/'+filename+'.h5' - hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r') - print('Reading checkpoint ',indexes[k]) - hdf.read(u, h5name) - time = hdf.attributes(h5name).to_dict()['timestamp'] - hdf.close() - uvec = u.vector().get_local() - ValuesTime[ix,:] = np.abs(uvec) - ix +=1 - - CoordPeak = np.where(ValuesTime==np.max(ValuesTime)) - ValuesPeak = np.abs(ValuesTime[CoordPeak[0],:]) - freq,edges = np.histogram(ValuesPeak, bins=80, density=True) - #Saving the histogram - print('Saving at ' + output_path) - np.savetxt(output_path + 'hist_freq.txt', freq) - np.savetxt(output_path + 'hist_edges.txt', edges) - - if mode == 'perturbation': - u = Function(W) - unew = Function(W) - u.rename('velocity', outname) - unew.rename('velocity', outname) - if options['Perturbation']['xdmf']: - xdmf_u = XDMFFile(output_path+'u.xdmf') - - if not options['Perturbation']['type']['SNR']=='inf': - Noise = True - def Add_Noise(signal,SNR): - Psignal = signal**2 - Psignal_av = np.mean(Psignal) - Psignal_av_db = 10*np.log10(Psignal_av) - Pnoise_av_db = Psignal_av_db - SNR - Pnoise_av = 10**(Pnoise_av_db/10) - noise_std = np.sqrt(Pnoise_av) - noise = np.random.normal(0,noise_std,len(signal)) - return signal + noise - else: - Noise = False - - if not options['Perturbation']['type']['phase_contrast']==0: - Phase_Contrast = True - else: - Phase_Contrast = False - - noise_in_coil = options['Perturbation']['type']['coil'] - - umaxima = [] - for k in indexes: - path = checkpoint_path + str(k) + '/'+filename+'.h5' - hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r') - hdf.read(u, 'u/vector_0') - time = hdf.attributes('u/vector_0').to_dict()['timestamp'] - hdf.close() - uvec = u.vector().get_local() - umaxima.append(np.max(uvec)) - - ufactor = options['Perturbation']['type']['phase_contrast']/100 - VENC = np.floor(np.ceil(np.max(umaxima))*ufactor) - print('VENC chosen = ',VENC) - - for k in indexes: - path = checkpoint_path + str(k) + '/'+filename+'.h5' - hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r') - hdf.read(u, 'u/vector_0') - time = hdf.attributes('u/vector_0').to_dict()['timestamp'] - hdf.close() - uvec = u.vector().get_local() - - if Phase_Contrast: - #gamma = 267.513e6 # rad/Tesla/sec Gyromagnetic ratio for H nuclei - B0 = 1.5 # Tesla Magnetic Field Strenght - TE = 5e-3 # Echo-time - Phi1 = 1*B0*TE + 0*uvec - Phi2 = 1*B0*TE + np.pi*uvec/VENC - M1 = np.cos(Phi1) + 1j*np.sin(Phi1) - M2 = np.cos(Phi2) + 1j*np.sin(Phi2) - - if noise_in_coil: - SNR = options['Perturbation']['type']['SNR'] - m1x_new = Add_Noise(np.real(M1),SNR) - m1y_new = Add_Noise(np.imag(M1),SNR) - m2x_new = Add_Noise(np.real(M2),SNR) - m2y_new = Add_Noise(np.imag(M2),SNR) - M1_new = m1x_new + 1j*m1y_new - M2_new = m2x_new + 1j*m2y_new - uvec = (np.angle(M2_new) - np.angle(M1_new))*VENC/np.pi - unew.vector()[:] = uvec - else: - uvec = (np.angle(M2) - np.angle(M1))*VENC/np.pi - else: - if noise_in_coil: - raise Exception('In order to perturb in coils some PC should be selected') - - - if not noise_in_coil: - if Noise: - SNR = options['Perturbation']['type']['SNR'] - unew.vector()[:] = Add_Noise(uvec,SNR) - else: - unew.vector()[:] = uvec - - print('Writing checkpoint number ',k) - write_path = output_path + 'checkpoint/{i}/'.format(i=k) - hdf2 = HDF5File(MESH['mesh'].mpi_comm(), write_path + 'u.h5', 'w') - hdf2.write(unew, '/u', time) - hdf2.close() - - if options['Perturbation']['xdmf']: - xdmf_u.write(unew, time) - - if mode == 'interpolation': - dt_new = options['Temporal-Interpolation']['dt_new'] - dt = options['Temporal-Interpolation']['dt'] - Nfolder = len(indexes) - Tf = dt*(Nfolder-1) - Nfolder_new = np.int(Tf/dt_new + 1) - indexes_new = [k for k in range(1,Nfolder_new+1)] - - u = Function(W) - unew = Function(W) - u.rename('velocity', outname) - unew.rename('velocity', outname) - - if options['Temporal-Interpolation']['xdmf']: - xdmf_u = XDMFFile(output_path+'u.xdmf') - - uvec = {} - # Reading first all the points for every original time-steps - for k in indexes: - path = checkpoint_path + str(k) + '/'+filename+'.h5' - print('Reading timestep number: ' + str(k)) - hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r') - if 'w' in filename: - hdf.read(u, 'w/vector_0') - else: - hdf.read(u, 'u/vector_0') - hdf.close() - uvec[k] = u.vector().get_local() - - times_old = np.linspace(0,Tf,Nfolder) - times_new = np.linspace(0,Tf,Nfolder_new) - velnodes = np.zeros(times_old.size) - velnodes_new = np.zeros(times_new.size) - uvec_new = {} - - for k in indexes_new: - uvec_new[k] = np.zeros(uvec[1].size) - - from scipy.interpolate import interp1d - print('Interpolating in every node across time') - # FOR every single node!!! - for l in range(len(uvec[1])): - for k in range(len(velnodes)): - velnodes[k] = uvec[indexes[k]][l] - inter_f = interp1d(times_old, velnodes , kind='cubic') - velnodes_new = inter_f(times_new) - for k in range(len(velnodes_new)): - uvec_new[indexes_new[k]][l] = velnodes_new[k] - print('Interpolation done') - - for k in indexes_new: - print('Writing timestep number: ' + str(k) ) - - unew.vector()[:] = uvec_new[k][:] - - write_path = output_path + 'checkpoint/{i}/'.format(i=k) - hdf2 = HDF5File(MESH['mesh'].mpi_comm(), write_path + 'u.h5', 'w') - hdf2.write(unew, '/u', float((k-1)*dt_new)) - hdf2.close() - - if options['Temporal-Interpolation']['xdmf']: - xdmf_u.write(unew, (k-1)*dt_new) - - if mode == 'average': - N_av = options['Temporal-Average']['subsampling_rate'] - dt = options['Temporal-Average']['dt'] - ks = 1 - mean_fill = False - u = Function(W) - usub = Function(W) - u.rename('velocity', outname) - usub.rename('velocity', outname) - if options['Temporal-Average']['xdmf']: - xdmf_u = XDMFFile(output_path+'test.xdmf') - - for k in range(len(indexes)): - path = checkpoint_path + str(indexes[k]) + '/'+filename+'.h5' - print('Reading timestep number: ' + str(indexes[k])) - hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r') - if 'w' in filename: - hdf.read(u, 'w/vector_0') - else: - hdf.read(u, 'u/vector_0') - hdf.close() - - if not mean_fill: - mean_u = u.vector().get_local()/N_av - mean_fill = True - else: - mean_u += u.vector().get_local()/N_av - - if np.mod(k+1,N_av) == 0: - mean_fill = False - usub.vector()[:] = mean_u - print('saving timestep number: ' , ks) - write_path = output_path + 'checkpoint/{i}/'.format(i=ks) - hdf2 = HDF5File(MESH['mesh'].mpi_comm(), write_path + 'u.h5', 'w') - hdf2.write(usub, '/u', float((k-N_av+1)*dt)) - hdf2.close() - if options['Temporal-Average']['xdmf']: - xdmf_u.write(usub, (k-N_av+1)*dt) - ks +=1 - - if mode == 'u': - xdmf_u = XDMFFile(output_path+outname+'.xdmf') - for k in range(0, len(indexes), 1): - path = checkpoint_path + str(indexes[k]) + '/'+filename+'.h5' - u = Function(V) - u.rename('velocity', outname) - hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r') - hdf.read(u, 'u/vector_0') - hdf.close() - if flow: - for bb in bnds: - QQ[bb].append(assemble(dot(u, n)*ds(bb))) - - xdmf_u.write(u, k) - te = te + dt - - if flow: - for bb in bnds: - np.savetxt(output_path+'flowrate'+str(bb)+'.txt', QQ[bb]) - - if mode == 'w': - xdmf_u = XDMFFile(output_path+outname+'.xdmf') - for k in range(0, len(indexes), 1): - path = checkpoint_path + str(indexes[k]) + '/'+filename+'.h5' - v = Function(V) - v.rename('waux', outname) - hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r') - hdf.read(v, 'w/vector_0') - hdf.close() - wvec = v.vector().get_local() - wnorm = wvec - np.mean(wvec) - v.vector()[:] = wnorm - xdmf_u.write(v, k) - - if mode == 'gradw': - xdmf_u = XDMFFile(output_path+outname+'.xdmf') - for k in range(0, len(indexes), 1): - path = checkpoint_path + str(indexes[k]) + '/'+filename+'.h5' - v = Function(V) - gw = Function(W) - gw.rename(outname, outname) - hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r') - hdf.read(v, 'w/vector_0') - hdf.close() - wvec = v.vector().get_local() - wnorm = wvec - np.mean(wvec) - v.vector()[:] = wnorm - gw.assign(project(sqrt(inner(grad(v), grad(v))), W)) - xdmf_u.write(gw, k) - - if mode == 'p' or mode == 'p_cib': - xdmf_p = XDMFFile(output_path+outname+'.xdmf') - for k in range(0, len(indexes), 1): - path = checkpoint_path + str(indexes[k]) + '/'+filename+'.h5' - p = Function(W) - if mode == 'p': - barye2mmHg = -1/1333.22387415 - if mode == 'p_cib': - barye2mmHg = -0.00750062 - p.rename('pressure', outname) - hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r') - hdf.read(p, 'p/vector_0') - hdf.close() - p1vec = p.vector().get_local()*barye2mmHg - p.vector()[:] = p1vec - xdmf_p.write(p, k) - - if mode == 'divu': - xdmf_u = XDMFFile(output_path+outname+'.xdmf') - for k in range(0, len(indexes), 1): - path = checkpoint_path + str(indexes[k]) + '/'+filename+'.h5' - v = Function(V) - dv = Function(W) - dv.rename('div', outname) - hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r') - hdf.read(v, 'u/vector_0') - hdf.close() - dv.assign(project(div(v), W)) - xdmf_u.write(dv, k) - -def ERRORmap(MESH, mode, outpath, reference_path, checkpoint_path, refname,comname,options): - - from dolfin import HDF5File - V = MESH['FEM'] - V1 = MESH['FEM'].sub(1).collapse() - W = MESH['FEM'].sub(0).collapse() - DGs = FunctionSpace( MESH['mesh'], 'DG',0) - cellv = CellVolume(MESH['mesh']) - h = CellDiameter(MESH['mesh']) - unsort_indexes = os.listdir(checkpoint_path) - indexes = [int(x) for x in unsort_indexes] - indexes.sort() - unsort_indexes0 = os.listdir(reference_path) - indexes0 = [int(x) for x in unsort_indexes0] - indexes0.sort() - - - if len(indexes)!=len(indexes0): - raise Exception('The lengh of the checkpoints are not the same!') - - if mode == 'algebra': - outname = options['Algebra']['outname'] - output = XDMFFile(outpath+outname+'.xdmf') - checkpoint = options['Algebra']['checkpoint'] - v1 = Function(W) - v2 = Function(W) - vout = Function(W) - vout.rename('velocity','velocity') - vout_vec = np.zeros(vout.vector().get_local().size) - times_vec = np.zeros(len(indexes0)) - vdict1 = {} - vdict2 = {} - if options['Algebra']['mode'] == 'aliasing': - print('Assuming the corrector in the second path entered') - - # Reading all the timestamps first - for k in range(len(indexes)): - path_v1 = reference_path + str(indexes0[k]) + '/'+refname+'.h5' - path_v2 = checkpoint_path + str(indexes[k]) + '/'+comname+'.h5' - hdf_v1 = HDF5File(MESH['mesh'].mpi_comm(), path_v1, 'r') - - if 'w' in refname: - hdf_v1.read(v1, 'w/vector_0') - time = hdf_v1.attributes('w/vector_0').to_dict()['timestamp'] - else: - hdf_v1.read(v1, 'u/vector_0') - time = hdf_v1.attributes('u/vector_0').to_dict()['timestamp'] - times_vec[k] = time - hdf_v2 = HDF5File(MESH['mesh'].mpi_comm(), path_v2, 'r') - if 'w' in comname: - hdf_v2.read(v2, 'w/vector_0') - else: - hdf_v2.read(v2, 'u/vector_0') - - print('computing algebra for the time',time) - hdf_v1.close() - hdf_v2.close() - vdict1[k] = v1.vector().get_local() - vdict2[k] = v2.vector().get_local() - - - for k in range(len(indexes)): - if options['Algebra']['mode'] == '+': - vout.vector()[:] = vdict1[k] + vdict2[k] - elif options['Algebra']['mode'] == '-': - vout.vector()[:] = vdict1[k] - vdict2[k] - elif options['Algebra']['mode'] == 'aliasing': - VENC = options['Algebra']['VENC'] - aliasing = False - for l in range(len(vout_vec)): - if k>0: - #mean1 = abs(np.mean(vdict2[0][:])) - #mean2 = abs(np.mean(vdict2[1][:])) - #mean3 = abs(np.mean(vdict2[2][:])) - #treshold = 10*max(mean1,mean2,mean3) - #if abs(np.mean(vdict2[k][:]))>treshold: - # vout_vec[l] = 0 - if vdict1[k][l]-vdict1[k-1][l] < -VENC: - vdict1[k][l] = vdict1[k][l] + 2*VENC - vout_vec[l] = vdict1[k][l] - else: - vout_vec[l] = vdict1[k][l] - else: - vout_vec[l] = vdict1[k][l] # first case is suppouse to be aliasing-free - - vout.vector()[:] = vout_vec - else: - raise Exception('Not supported operation between vectors!') - - output.write(vout, times_vec[k]) - if checkpoint: - path = outpath + 'checkpoint/' + str(indexes[k]) + '/' + 'u.h5' - inout.write_HDF5_data(MESH['mesh'].mpi_comm(), path , vout, '/u', t=times_vec[k]) - - if mode =='curves': - - if options['Error-curves']['undersampling']>1: - print('undersampling in the reference measurements!') - n_under = options['Error-curves']['undersampling'] - indexes0 = indexes0[0::n_under] - - u = Function(W) - w = Function(W) - ones = interpolate(Constant(1), V.sub(1).collapse()) - L_sys = assemble(ones*dx) - VENC = options['Error-curves']['VENC'] - - for typs in options['Error-curves']['type']: - ucomp = [] - wcomp = [] - div_array = [] - varu = [] - times = [] - - if typs == 'l2_comp': - wy = Function(V1) - wz = Function(V1) - comname2 = comname - wx_array = [] - wy_array = [] - wz_array = [] - - if typs == 'utrue-uobs': - u_path = options['Error-curves']['true_check'] + 'checkpoint/' - w_path = options['Error-curves']['ref_check'] + 'checkpoint/' - comname2 = 'u' - varux_array = [] - varuy_array = [] - varuz_array = [] - else: - u_path = reference_path - w_path = checkpoint_path - - dt = options['Error-curves']['dt'] - - for k in range(1,len(indexes)): - path_w = w_path + str(indexes[k]) + '/'+comname2+'.h5' - path_u = u_path + str(indexes0[k]) + '/'+refname+'.h5' - u.rename('meas', 'meas') - w.rename('w', 'w') - hdf_w = HDF5File(MESH['mesh'].mpi_comm(),path_w,'r') - if 'w' in comname2: - hdf_w.read(w, 'w/vector_0') - else: - hdf_w.read(w, 'u/vector_0') - if typs != 'l2_comp': - hdf_u = HDF5File(MESH['mesh'].mpi_comm(), path_u, 'r') - hdf_u.read(u, 'u/vector_0') - hdf_u.close() - u_vec = u.vector().get_local() - hdf_w.close() - w_vec = w.vector().get_local() - print('computing error in timestep numer',k) - - if typs == 'mean': - ucomp.append(np.mean(abs(u_vec))) - wcomp.append(np.mean(abs(w_vec))) - elif typs == 'max': - ucomp.append(np.max(abs(u_vec))) - wcomp.append(np.max(abs(w_vec))) - elif typs == 'norm2': - ucomp.append(np.sqrt(assemble(dot(u,u)*dx)/L_sys)) - wcomp.append(np.sqrt(assemble(dot(w,w)*dx)/L_sys)) - elif typs == 'grad': - wcomp.append(np.sqrt(assemble(inner(grad(w),grad(w))*dx)/L_sys)) - elif typs == 'l2_comp': - wx = w.sub(0,deepcopy=True) - wy = w.sub(1,deepcopy=True) - wz = w.sub(2,deepcopy=True) - wx_array.append(np.sqrt(assemble(wx*wx*dx)/L_sys)/VENC) - wy_array.append(np.sqrt(assemble(wy*wy*dx)/L_sys)/VENC) - wz_array.append(np.sqrt(assemble(wz*wz*dx)/L_sys)/VENC) - elif typs == 'div': - div_array.append(np.sqrt(assemble(div(u)**2/h*dx)/L_sys)/VENC) - elif typs == 'utrue-uobs': - wx = w.sub(0,deepcopy=True) - wy = w.sub(1,deepcopy=True) - wz = w.sub(2,deepcopy=True) - ux = u.sub(0,deepcopy=True) - uy = u.sub(1,deepcopy=True) - uz = u.sub(2,deepcopy=True) - varux_array.append(np.sqrt(assemble((ux-wx)**2*dx)/L_sys)/VENC) - varuy_array.append(np.sqrt(assemble((uy-wy)**2*dx)/L_sys)/VENC) - varuz_array.append(np.sqrt(assemble((uz-wz)**2*dx)/L_sys)/VENC) - else: - raise Exception('No defined type for curve printing!') - - times.append(k*dt) - - if typs == 'grad': - np.savetxt(outpath+'w' +typs+'.txt',wcomp) - elif typs == 'l2_comp': - np.savetxt(outpath+'wx.txt',wx_array) - np.savetxt(outpath+'wy.txt',wy_array) - np.savetxt(outpath+'wz.txt',wz_array) - elif typs == 'div': - np.savetxt(outpath+'div.txt',div_array) - elif typs == 'utrue-uobs': - np.savetxt(outpath+'varux.txt',varux_array) - np.savetxt(outpath+'varuy.txt',varuy_array) - np.savetxt(outpath+'varuz.txt',varuz_array) - else: - np.savetxt(outpath+'u' +typs+'.txt',ucomp) - np.savetxt(outpath+'w' +typs+'.txt',wcomp) - - np.savetxt(outpath+'times.txt',times) - - if mode == 'h5': - xdmf_u = XDMFFile(outpath+'meas.xdmf') - #xdmf_ucor = XDMFFile(output_path+'ucor.xdmf') - xdmf_w = XDMFFile(outpath+'w.xdmf') - ds = Measure("ds", subdomain_data=MESH['boundaries']) - - u = Function(W) - w = Function(W) - #ucor = Function(W) - - for k in range(0, len(indexes), 1): - path_w = checkpoint_path + str(indexes[k]) + '/'+comname+'.h5' - path_u = reference_path + str(indexes0[k]) + '/'+refname+'.h5' - - u.rename('meas', 'meas') - w.rename('w', 'w') - #ucor.rename('ucor', 'ucor') - hdf_w = HDF5File(MESH['mesh'].mpi_comm(),path_w,'r') - hdf_w.read(w, 'w/vector_0') - hdf_u = HDF5File(MESH['mesh'].mpi_comm(), path_u, 'r') - hdf_u.read(u, 'u/vector_0') - - hdf_u.close() - hdf_w.close() - - #u_vec = u.vector().get_local() - #w_vec = w.vector().get_local() - #ucor.vector()[:] = u_vec + w_vec - - xdmf_u.write(u, k) - xdmf_w.write(w, k) - - if mode == 'u': - colormap = XDMFFile(outpath+'u.xdmf') - #ds = Measure("ds", subdomain_data=MESH['boundaries']) - u = Function(W) - #P2 = FunctionSpace( MESH['mesh'], 'P',1) - cm = Function(DGs) - vs = TestFunction(DGs) - - for k in range(len(indexes)): - path_u = reference_path + str(indexes0[k]) + '/'+refname+'.h5' - u.rename('meas', 'meas') - cm.rename('vel','vel') - hdf_u = HDF5File(MESH['mesh'].mpi_comm(), path_u, 'r') - hdf_u.read(u, 'u/vector_0') - time = hdf_u.attributes('u/vector_0').to_dict()['timestamp'] - print('making the colormap in the time',time) - hdf_u.close() - r = assemble(dot(u,u)*vs/cellv*dx).get_local() - cm.vector()[:] = np.sqrt(np.abs(r)) - #cm.vector()[:] = np.sqrt(r) - colormap.write(cm, time) - - if mode == 'w': - colormap = XDMFFile(outpath+'w.xdmf') - #ds = Measure("ds", subdomain_data=MESH['boundaries']) - w = Function(W) - #P2 = FunctionSpace( MESH['mesh'], 'P',2) - cm = Function(DGs) - vs = TestFunction(DGs) - - for k in range(len(indexes)): - path_w = checkpoint_path + str(indexes0[k]) + '/'+comname+'.h5' - w.rename('cor', 'cor') - cm.rename('w','w') - hdf_w = HDF5File(MESH['mesh'].mpi_comm(), path_w, 'r') - hdf_w.read(w, 'w/vector_0') - time = hdf_w.attributes('w/vector_0').to_dict()['timestamp'] - print('making the colormap in the time',time) - hdf_w.close() - r = assemble(dot(w,w)*vs/cellv*dx).get_local() - cm.vector()[:] = np.sqrt(np.abs(r)) - colormap.write(cm, time) - - if mode == 'w/u': - colormap = XDMFFile(outpath+'w_u.xdmf') - #ds = Measure("ds", subdomain_data=MESH['boundaries']) - u = Function(W) - w = Function(W) - cm = Function(DGs) - vs = TestFunction(DGs) - - for k in range(len(indexes)): - path_w = checkpoint_path + str(indexes[k]) + '/'+comname+'.h5' - path_u = reference_path + str(indexes0[k]) + '/'+refname+'.h5' - u.rename('meas', 'meas') - w.rename('w', 'w') - cm.rename('w_u','w_u') - hdf_w = HDF5File(MESH['mesh'].mpi_comm(),path_w,'r') - hdf_w.read(w, 'w/vector_0') - hdf_u = HDF5File(MESH['mesh'].mpi_comm(), path_u, 'r') - hdf_u.read(u, 'u/vector_0') - time = hdf_u.attributes('u/vector_0').to_dict()['timestamp'] - print('making the colormap in the time',time) - hdf_u.close() - hdf_w.close() - rw = assemble(dot(w,w)*vs/cellv*dx).get_local() - ru = assemble(dot(u,u)*vs/cellv*dx).get_local() - cm.vector()[:] = np.sqrt(rw)/np.sqrt(ru) - colormap.write(cm, time) - - if mode == 'dot': - colormap = XDMFFile(outpath+'dot.xdmf') - #ds = Measure("ds", subdomain_data=MESH['boundaries']) - u = Function(W) - w = Function(W) - cm = Function(V1) - u.rename('meas', 'meas') - w.rename('w', 'w') - - for k in range(len(indexes)): - path_w = checkpoint_path + str(indexes[k]) + '/'+comname+'.h5' - path_u = reference_path + str(indexes0[k]) + '/'+refname+'.h5' - hdf_w = HDF5File(MESH['mesh'].mpi_comm(),path_w,'r') - hdf_w.read(w, 'w/vector_0') - hdf_u = HDF5File(MESH['mesh'].mpi_comm(), path_u, 'r') - hdf_u.read(u, 'u/vector_0') - time = hdf_u.attributes('u/vector_0').to_dict()['timestamp'] - print('making the colormap in the time',time) - hdf_u.close() - hdf_w.close() - cm = project(inner(u,w),V1) - cm.rename('dot','dot') - colormap.write(cm, time) - - if mode == 'div(u)': - - colormap = XDMFFile(outpath+'divu.xdmf') - #ds = Measure("ds", subdomain_data=MESH['boundaries']) - u = Function(W) - cm = Function(DGs) - vs = TestFunction(DGs) - - for k in range(len(indexes)): - path_u = reference_path + str(indexes0[k]) + '/'+refname+'.h5' - u.rename('meas', 'meas') - cm.rename('divu','divu') - hdf_u = HDF5File(MESH['mesh'].mpi_comm(), path_u, 'r') - hdf_u.read(u, 'u/vector_0') - time = hdf_u.attributes('u/vector_0').to_dict()['timestamp'] - print('making the colormap in the time',time) - hdf_u.close() - r = assemble(div(u)**2*vs/cellv*dx).get_local() - cm.vector()[:] = np.sqrt(r) - colormap.write(cm, time) - - if mode == 'div(u)/u': - - colormap = XDMFFile(outpath+'divu_u.xdmf') - #ds = Measure("ds", subdomain_data=MESH['boundaries']) - u = Function(W) - cm = Function(DGs) - vs = TestFunction(DGs) - - surface = np.zeros(cm.vector().get_local().size) - volume = np.zeros(cm.vector().get_local().size) - Cells = MESH['mesh'].cells() - Points = MESH['mesh'].coordinates() - for k in range(MESH['mesh'].num_cells()): - A = Points[Cells[k]][0] - B = Points[Cells[k]][1] - C = Points[Cells[k]][2] - D = Points[Cells[k]][3] - # volume - Vol = 1/6*np.abs(np.inner(A-D, np.cross(B-D, C-D))) - volume[k] = Vol - # surface area - SA = 0 - for l in range(4): - if l == 0: - [P1, P2, P3, P4] = [A, B, C, D] - if l == 1: - [P1, P2, P3, P4] = [D, A, B, C] - if l == 2: - [P1, P2, P3, P4] = [C, D, A, B] - if l == 3: - [P1, P2, P3, P4] = [B, C, D, A] - - SA += 0.5*np.linalg.norm(np.cross(P2-P1, P3-P1)) - surface[k] = SA - - - - for k in range(len(indexes)): - path_u = reference_path + str(indexes0[k]) + '/'+refname+'.h5' - u.rename('meas', 'meas') - cm.rename('divu','divu') - hdf_u = HDF5File(MESH['mesh'].mpi_comm(), path_u, 'r') - hdf_u.read(u, 'u/vector_0') - time = hdf_u.attributes('u/vector_0').to_dict()['timestamp'] - print('making the colormap in the time',time) - hdf_u.close() - rdiv = assemble(div(u)**2*vs*dx).get_local() - ru = assemble(dot(u,u)*vs*dx).get_local() - cm.vector()[:] = np.sqrt(rdiv)/np.sqrt(ru)*volume/surface/np.sqrt(120) - colormap.write(cm, time) - - if mode == 'error_u': - xdmf_u = XDMFFile(output_path+'error_u.xdmf') - for k in range(0, len(indexes), 1): - path = checkpoint_path + str(indexes[k]) + '/'+comname+'.h5' - path0 = reference_path + str(indexes0[k]) + '/'+refname+'.h5' - - if k < 10: - path = checkpoint_path + '0'+str(indexes[k]) + '/u.h5' - path0 = reference_path + '0'+str(indexes[k]) + '/u.h5' - - u = Function(V) - u0 = Function(V) - - eu = Function(W) - eu.rename('error_u', 'error_u') - - comm = MPI.COMM_WORLD - hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r') - hdf.read(u, 'u/vector_0') - hdf0 = HDF5File(MESH['mesh'].mpi_comm(), path0, 'r') - hdf0.read(u0, 'u/vector_0') - hdf0.close() - hdf.close() - - eu.assign(project(sqrt(inner(u-u0, u-u0)), W)) - - xdmf_u.write(eu, te) - te = te + dt - numt = numt + 1 - - if mode == 'error_p' or mode == 'error_p_cib': - xdmf_p = XDMFFile(output_path+outname+'.xdmf') - for k in range(0, len(indexes0), 1): - path = checkpoint_path + str(indexes[k]) + '/'+comname+'.h5' - path0 = reference_path + str(indexes0[k]) + '/'+refname+'.h5' - - if refname == 'p': - if k < 10 and k > 0: - path0 = reference_path + '0' + \ - str(indexes0[k]) + '/'+refname+'.h5' - - #zero_point = [12.79353, 14.32866, 6.51101] - #zero_point = np.array(zero_point) - #ndim = W.mesh().topology().dim() - #zero_point = W.tabulate_dof_coordinates().reshape((-1, ndim))[0] - - if 'cib' in mode: - barye2mmHg = 0.00750062 - else: - barye2mmHg = 1/1333.22387415 - p = Function(W) - p0 = Function(W) - - p.rename('error_p', outname) - - hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r') - hdf.read(p, 'p/vector_0') - hdf0 = HDF5File(MESH['mesh'].mpi_comm(), path0, 'r') - hdf0.read(p0, 'p/vector_0') - - hdf.close() - hdf0.close() - - #p0.vector()[:] -= p0(zero_point) - pvec = p.vector().get_local() - p0vec = p0.vector().get_local() - - pvec = pvec - np.mean(pvec) - p0vec = p0vec - np.mean(p0vec) - - errvec = np.sqrt((pvec - p0vec)**2)*barye2mmHg - #errvec = np.abs((pvec - p0vec)/(p0vec)) - - p.vector()[:] = errvec - #p.vector()[:] = p0vec - - xdmf_p.write(p, te) - te = te + dt - numt = numt + 1 - -def SEQCIBH5(pathtocib, pathtocibseq, ndata): - - import csv - - times = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, - 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25] - #times = [8] - - if ndata == 1: - pathtocib += '9_mm_AoCo_phantom_rest_vox0.9mm/' - pathtocibseq += '9AoCo_rest_vox0.9/' - elif ndata == 2: - pathtocib += '9_mm_AoCo_phantom_stress_vox0.9mm/' - pathtocibseq += '9AoCo_stress_vox0.9/' - elif ndata == 3: - pathtocib += '11_mm_AoCo_phantom_rest_vox0.9mm/' - pathtocibseq += '11AoCo_rest_vox0.9/' - elif ndata == 4: - pathtocib += '11_mm_AoCo_phantom_stress_vox0.9mm/' - pathtocibseq += '11AoCo_stress_vox0.9/' - else: - raise Exception('Data file not recognize!') - - xdmf_v = XDMFFile(pathtocibseq + 'vel.xdmf') - - mesh = Mesh(pathtocib+'aorta.xml') - Evel = VectorElement('Lagrange', mesh.ufl_cell(), 1) - V = FunctionSpace(mesh, Evel) - boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1) - - VX = V.sub(0).collapse() - VY = V.sub(1).collapse() - VZ = V.sub(2).collapse() - vv_x = Function(VX) - vv_y = Function(VY) - vv_z = Function(VZ) - v = Function(V) - - v.rename('velocity', 'vel') - mx = dof_to_vertex_map(VX) - my = dof_to_vertex_map(VY) - mz = dof_to_vertex_map(VZ) - - for ti in times: - print('reading aorta '+str(ti) + '.csv') - with open(pathtocib + 'aorta'+str(ti)+'.csv', 'r') as csvfile: - mylist = [row[0] for row in csv.reader(csvfile, delimiter=';')] - - Values = np.array(mylist) - file_x = np.zeros([len(Values)]) - file_y = np.zeros([len(Values)]) - file_z = np.zeros([len(Values)]) - - for l in range(len(Values)): - row = Values[l].split(',') - file_x[l] = float(row[0]) - file_y[l] = float(row[1]) - file_z[l] = float(row[2]) - - # print(np.max(file_z)) - #print(np.where(file_z==np.max(file_z)) ) - # print(np.max(file_x[np.where(file_z==np.max(file_z))])) - # print(np.max(file_y[np.where(file_z==np.max(file_z))])) - - #S = np.load(pathtocibseq) - #ux = S['x'] - #uy = S['y'] - #uz = S['z'] - #ux = ux.transpose((0,2,1,3)) - #uy = uy.transpose((0,2,1,3)) - #uz = uz.transpose((0,2,1,3)) - - vv_x.vector()[:] = file_x[mx] - vv_y.vector()[:] = file_y[my] - vv_z.vector()[:] = file_z[mz] - - assign(v.sub(0), vv_x) - assign(v.sub(1), vv_y) - assign(v.sub(2), vv_z) - xdmf_v.write(v, ti) - # LagrangeInterpolator.interpolate(vv_x,v1_x) - # LagrangeInterpolator.interpolate(vv_y,v1_y) - # LagrangeInterpolator.interpolate(vv_z,v1_z) - - xdmf_v.close() - -def SqtoH5(BOX, MESH, Sqx, Sqy, Sqz, output_path, uname): - - xdmf_u = XDMFFile(output_path+uname+'.xdmf') - Xt = BOX['nodes'] - SqXt = BOX['seq'] - P1 = BOX['FEM'] - V = MESH['FEM'] - [Lx, Ly, Lz, Lt] = Sqx.shape - - VX = V.sub(0).collapse() - VY = V.sub(1).collapse() - VZ = V.sub(2).collapse() - - vv_x = Function(VX) - vv_y = Function(VY) - vv_z = Function(VZ) - - m = dof_to_vertex_map(P1) - - if rank == 0: - print('{SQtoH5} total number of timesteps: ' + str(Lt)) - - v = Function(V) - v.rename('velocity', uname) - dt = (Lt-1) - te = 0 - - for t in range(Lt): - - if rank == 0: - print('timestep number', t) - v1_x = Function(P1) - v1_y = Function(P1) - v1_z = Function(P1) - - values_x = np.zeros(v1_x.vector().get_local().shape) - values_y = np.zeros(v1_y.vector().get_local().shape) - values_z = np.zeros(v1_z.vector().get_local().shape) - - S0x = Sqx[:, :, :, t] - S0y = Sqy[:, :, :, t] - S0z = Sqz[:, :, :, t] - for k in range(Xt.shape[0]): - values_x[k] = S0x[SqXt[k, 0], SqXt[k, 1], SqXt[k, 2]] - values_y[k] = S0y[SqXt[k, 0], SqXt[k, 1], SqXt[k, 2]] - values_z[k] = S0z[SqXt[k, 0], SqXt[k, 1], SqXt[k, 2]] - - v1_x.vector()[:] = values_x - v1_y.vector()[:] = values_y - v1_z.vector()[:] = values_z - - LagrangeInterpolator.interpolate(vv_x, v1_x) - LagrangeInterpolator.interpolate(vv_y, v1_y) - LagrangeInterpolator.interpolate(vv_z, v1_z) - - assign(v.sub(0), vv_x) - assign(v.sub(1), vv_y) - assign(v.sub(2), vv_z) - - xdmf_u.write(v, te) - te = te+dt - -def ESTIMpressure(MESH, outpath, checkpoint_path, filename, options): - - #from dolfin import HDF5File - V = MESH['FEM'] - W = MESH['FEM'].sub(1).collapse() - p = Function(W) - one_mesh = interpolate(Constant(1), W) - dt = options['Estim_Pressure']['dt'] - p_drop_lst = [] - time_ = [] - # Checkpoints folders - unsort_indexes = os.listdir(checkpoint_path) - indexes = [int(x) for x in unsort_indexes] - indexes.sort() - - from pathlib import Path - import pickle - - # Create the Spheres - if options['Estim_Pressure']['method'] == 'spheres': - import mshr - Npts0 = options['Estim_Pressure']['spheres'][0]['Npts'] - Npts1 = options['Estim_Pressure']['spheres'][1]['Npts'] - center0 = options['Estim_Pressure']['spheres'][0]['center'] - center1 = options['Estim_Pressure']['spheres'][1]['center'] - radius0 = options['Estim_Pressure']['spheres'][0]['radius'] - radius1 = options['Estim_Pressure']['spheres'][1]['radius'] - x0 = Point(np.array(center0, dtype=float)) - x1 = Point(np.array(center1, dtype=float)) - sphere0 = mshr.Sphere(x0, radius0) - sphere1 = mshr.Sphere(x1, radius1) - mesh_s1 = mshr.generate_mesh(sphere0, Npts0) - mesh_s2 = mshr.generate_mesh(sphere1, Npts1) - # New elements - VS1 = FunctionSpace(mesh_s1, FiniteElement('P', mesh_s1.ufl_cell(), 1)) - VS2 = FunctionSpace(mesh_s2, FiniteElement('P', mesh_s2.ufl_cell(), 1)) - s1 = Function(VS1) - s2 = Function(VS2) - LagrangeInterpolator.interpolate(s1, one_mesh) - LagrangeInterpolator.interpolate(s2, one_mesh) - vol1 = assemble(s1*dx) - vol2 = assemble(s1*dx) - - if options['Estim_Pressure']['method'] == 'boundaries': - boundaries = MESH['boundaries'] - ds = Measure("ds", subdomain_data=boundaries) - bnd1 = options['Estim_Pressure']['boundaries'][0] - bnd2 = options['Estim_Pressure']['boundaries'][1] - area1 = assemble(one_mesh*ds(bnd1)) - area2 = assemble(one_mesh*ds(bnd2)) - - - - for k in range(0, len(indexes)): - path = checkpoint_path + str(indexes[k]) + '/'+filename+'.h5' - hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r') - hdf.read(p, 'p/vector_0') - hdf.close() - - if options['Estim_Pressure']['method'] == 'spheres': - LagrangeInterpolator.interpolate(s1, p) - LagrangeInterpolator.interpolate(s2, p) - P1 = assemble(s1*dx)/vol1 - P2 = assemble(s2*dx)/vol2 - if options['Estim_Pressure']['method'] == 'boundaries': - P1 = assemble(p*ds(bnd1))/area1 - P2 = assemble(p*ds(bnd2))/area2 - - p_drop_lst.append(P2-P1) - time_.append(k*dt) - if rank == 0: - print('Pressure drop :', P2-P1) - - # Saving the Result - write_path = Path(outpath) - write_path.mkdir(exist_ok=True) - methods = filename[2:] - data = { - 'pdrop': np.array(p_drop_lst), - 'time': np.array(time_) - } - fpath = write_path.joinpath('pdrop_' + methods + '.dat') - pickle.dump(data, fpath.open('wb')) - -def OUTLETwind(MESH, output_path, checkpoint_path, filename, bnds): - - #from dolfin import HDF5File - V = MESH['FEM'].sub(0).collapse() - W = MESH['FEM'].sub(1).collapse() - # Checkpoints folders - unsort_indexes = os.listdir(checkpoint_path) - indexes = [int(x) for x in unsort_indexes] - indexes.sort() - QQ = {} - PP = {} - ds = Measure('ds', domain=MESH['mesh'], subdomain_data=MESH['boundaries']) - n = FacetNormal(MESH['mesh']) - ones = interpolate(Constant(1), W) - for bb in bnds: - PP[bb] = [] - QQ[bb] = [] - - for k in range(0, len(indexes)): - path_u = checkpoint_path + str(indexes[k]) + '/'+filename[0]+'.h5' - path_p = checkpoint_path + str(indexes[k]) + '/'+filename[1]+'.h5' - - u = Function(V) - p = Function(W) - #comm = MPI.COMM_WORLD - hdf_u = HDF5File(MESH['mesh'].mpi_comm(), path_u, 'r') - hdf_u.read(u, 'u/vector_0') - hdf_u.close() - hdf_p = HDF5File(MESH['mesh'].mpi_comm(), path_p, 'r') - hdf_p.read(p, 'p/vector_0') - hdf_p.close() - - print('Computing flows and pressures at timestep number ' + str(indexes[k])) - for bb in bnds: - if bb == 2: - QQ[bb].append(-assemble(dot(u, n)*ds(bb))) - else: - QQ[bb].append(assemble(dot(u, n)*ds(bb))) - - PP[bb].append(assemble(p*ds(bb))/1333.22387415 / - assemble(ones*ds(bb))) - - print('saving the figure at' + output_path) - from matplotlib import pyplot as plt - from matplotlib import rc - rc('text', usetex=True) - fig = plt.figure(figsize=(12, 5), dpi=100) - t = np.linspace(0, len(indexes), len(indexes)) - plt.subplot(1, 2, 2) - for bb in bnds: - plt.plot(t, QQ[bb], linewidth=1.2, linestyle='-', - label='$ Outlet: '+str(bb)+'$') - plt.xlabel(r'$frames $', fontsize=20) - plt.ylabel(r'$ Q $', fontsize=20) - plt.xlim([0, 1.05*max(t)]) - plt.legend(fontsize=14) - - plt.subplot(1, 2, 1) - for bb in bnds: - plt.plot(t, PP[bb], linewidth=1.2, linestyle='-', - label='$ Outlet: '+str(bb)+'$') - plt.xlabel(r'$frames $', fontsize=20) - plt.ylabel(r'$ P \ \ mmHg$', fontsize=20) - plt.xlim([0, 1.05*max(t)]) - plt.legend(fontsize=14) - - fig.savefig(output_path + 'Outlet_Windkessel', dpi=500) - -def ROUTINE(options): - - if 'Algebra' in options: - if options['Algebra']['apply']: - - if rank == 0: - print('Applying Algebra') - print('Choosen mode: ' + options['Algebra']['mode']) - - MESH = LOADmesh(options['Algebra']['meshpath']) - v1path = options['Algebra']['v1']['path'] + 'checkpoint/' - v1name = options['Algebra']['v1']['name'] - v2path = options['Algebra']['v2']['path'] + 'checkpoint/' - v2name = options['Algebra']['v2']['name'] - outpath = options['Algebra']['outpath'] - - ERRORmap(MESH, 'algebra', outpath, v1path, - v2path, v1name, v2name, options) - - if 'Outlet_Wind' in options: - if options['Outlet_Wind']['apply']: - if rank == 0: - print('--- Outlet Windkessel ---') - - meshpath = options['Outlet_Wind']['meshpath'] - out_path = options['Outlet_Wind']['out_path'] - filename = options['Outlet_Wind']['filename'] - checkpoint = options['Outlet_Wind']['checkpoint'] + 'checkpoint/' - bnds = options['Outlet_Wind']['bnds'] - MESH = LOADmesh(meshpath) - OUTLETwind(MESH, out_path, checkpoint, filename, bnds) - - if 'Colormap' in options: - if options['Colormap']['apply']: - - for mode in options['Colormap']['mode']: - if rank == 0: - print('Applying Colormap') - print('Choosen mode: ' + mode) - - MESH = LOADmesh(options['Colormap']['meshpath']) - upath = options['Colormap']['upath'] + 'checkpoint/' - wpath = options['Colormap']['wpath'] + 'checkpoint/' - wname = options['Colormap']['wname'] - uname = options['Colormap']['uname'] - outpath = options['Colormap']['outpath'] - - ERRORmap(MESH, mode, outpath, upath, - wpath, uname, wname, options) - - if 'Velocity' in options: - if options['Velocity']['apply']: - if rank == 0: - print('--- Reading Velocity ---') - - MESH = LOADmesh(options['Velocity']['meshpath']) - filename = options['Velocity']['filename'] - checkpoint_path = options['Velocity']['checkpoint'] + 'checkpoint/' - outpath = options['Velocity']['checkpoint'] - WORKcheck(MESH, 'u', outpath, checkpoint_path, - filename, filename, options) - - if 'Vel_from_check' in options: - if options['Vel_from_check']['apply']: - if rank == 0: - print('Applying Velocity--MAP from sequence') - - [Sqx, Sqy, Sqz] = LOADsequences(options['Velocity']['pathtoseq']) - BOX = LOADmesh('cib') - AORTA = LOADmesh(pathtomesh) - SqtoH5(BOX, AORTA, Sqx, Sqy, Sqz, output_path, uname) - - if 'W_map' in options: - if options['W_map']['apply']: - if rank == 0: - print('Applying W--MAP') - MESH = LOADmesh(options['meshpath']) - filename = options['W_map']['filename'] - outname = options['W_map']['out_name'] - mode = 'w' - WORKcheck(MESH, mode, output_path, checkpoint_path, - filename, outname, options) - - if 'GradW_map' in options: - if options['GradW_map']['apply']: - if rank == 0: - print('Applying Grad W--MAP') - MESH = LOADmesh(options['meshpath']) - filename = options['GradW_map']['filename'] - outname = options['GradW_map']['out_name'] - mode = 'gradw' - WORKcheck(MESH, mode, output_path, checkpoint_path, - filename, outname, options) - - if 'Pressure' in options: - if options['Pressure']['apply']: - if rank == 0: - print('---Applying Pressure---') - - MESH = LOADmesh(options['meshpath']) - outpath = options['Pressure']['checkpoint'] - checkpoint_path = options['Pressure']['checkpoint'] + 'checkpoint/' - filename = options['Pressure']['filename'] - WORKcheck(MESH, 'p', outpath, checkpoint_path, - filename, filename, options) - - if 'Error_P' in options: - if options['Error_P']['apply']: - if rank == 0: - print('Applying L2 error to Pressure') - - MESH = LOADmesh(options['meshpath']) - refname = options['Error_P']['refname'] - reference_path = options['Error_P']['refpath'] - filename = options['Error_P']['filename'] - outname = options['Error_P']['out_name'] - mode = 'error_p' - ERRORmap(MESH, mode, output_path, reference_path, - checkpoint_path, refname, filename, outname, options) - - if 'Estim_Pressure' in options: - if options['Estim_Pressure']['apply']: - if rank == 0: - print('Applying Pressure Estimator') - print('Method choosed: '+options['Estim_Pressure']['method']) - - MESH = LOADmesh(options['Estim_Pressure']['meshpath']) - filename = options['Estim_Pressure']['filename'] - outpath = options['Estim_Pressure']['checkpath'] - checkpoint_path = options['Estim_Pressure']['checkpath'] + 'checkpoint/' - - ESTIMpressure(MESH, outpath, checkpoint_path, filename, options) - - if 'Error-curves' in options: - if options['Error-curves']['apply']: - if rank == 0: - print('--- Error curves ---' ) - - MESH = LOADmesh(options['Error-curves']['meshpath']) - ref_check = options['Error-curves']['ref_check'] + 'checkpoint/' - com_check = options['Error-curves']['com_check'] + 'checkpoint/' - outpath = options['Error-curves']['outpath'] - refname = options['Error-curves']['ref_name'] - comname = options['Error-curves']['com_name'] - ERRORmap(MESH, 'curves', outpath, ref_check, com_check, - refname, comname, options) - - if 'Temporal-Average' in options: - if options['Temporal-Average']['apply']: - if rank == 0: - print('--- Temporal Average ---') - - MESH = LOADmesh(options['Temporal-Average']['meshpath']) - ref_check = options['Temporal-Average']['original_check'] + 'checkpoint/' - out_check = options['Temporal-Average']['out_check'] - READcheckpoint(MESH,'average', out_check,ref_check,'u','u',options) - - if 'Temporal-Interpolation' in options: - if options['Temporal-Interpolation']['apply']: - if rank == 0: - print('--- Temporal Interpolation ---') - - MESH = LOADmesh(options['Temporal-Interpolation']['meshpath']) - ref_check = options['Temporal-Interpolation']['original_check'] + 'checkpoint/' - out_check = options['Temporal-Interpolation']['out_check'] - READcheckpoint(MESH,'interpolation', out_check,ref_check,'u','u',options) - - if 'Perturbation' in options: - if options['Perturbation']['apply']: - if rank==0: - print('--- Perturbation in measurements ---') - - MESH = LOADmesh(options['Perturbation']['meshpath']) - checkpath = options['Perturbation']['checkpath'] + 'checkpoint/' - if type(options['Perturbation']['type']['SNR'])=='str': - name_folder = 'SNR' +options['Perturbation']['type']['SNR'] + 'V' + str(options['Perturbation']['type']['phase_contrast']) - else: - name_folder = 'SNR' +str(options['Perturbation']['type']['SNR']) + 'V' + str(options['Perturbation']['type']['phase_contrast']) - - out_check = options['Perturbation']['checkpath'] + name_folder + '/' - READcheckpoint(MESH,'perturbation', out_check,checkpath,'u','u',options) - - if 'Histograms' in options: - if options['Histograms']['apply']: - if rank==0: - print('--- Histograms in measurements ---') - - MESH = LOADmesh(options['Histograms']['meshpath']) - checkpath = options['Histograms']['checkpath'] + 'checkpoint/' - outpath = options['Histograms']['checkpath'] - field_name = options['Histograms']['field_name'] - READcheckpoint(MESH,'Histograms', outpath, checkpath,field_name,field_name,options) - - -if __name__ == '__main__': - - comm = MPI.COMM_WORLD - size = comm.Get_size() - rank = comm.Get_rank() - - if len(sys.argv) > 1: - if os.path.exists(sys.argv[1]): - inputfile = sys.argv[1] - if rank == 0: - print('Found input file ' + inputfile) - else: - raise Exception('Command line arg given but input file does not exist:' - ' {}'.format(sys.argv[1])) - else: - raise Exception('An input file is required as argument!') - - options = inout.read_parameters(inputfile) - ROUTINE(options) diff --git a/codes/SENSE.py b/codes/SENSE.py deleted file mode 100644 index 83b82f2..0000000 --- a/codes/SENSE.py +++ /dev/null @@ -1,115 +0,0 @@ -import numpy as np -from numpy import linalg as LA -import sys -from mpi4py import MPI -comm = MPI.COMM_WORLD -size = comm.Get_size() -rank = comm.Get_rank() - -# SENSE: Simulation of SENSitive Encoding algorithm proposed by K. Pruessmann, et. al. in: -# "SENSE: Sensitivity Enconding for Fast MRI" Mag. Res. in Medicine 42. (1999) -# written by Jeremias Garay (j.e.garay.labra@rug.nl) - -def Sensitivity_Map(shape): - - [Nx,Ny,Nz] = shape - [X,Y,Z] = np.meshgrid(np.linspace(0,Ny,Ny),np.linspace(0,Nx,Nx),np.linspace(0,Nz,Nz)) - Xsense1 = (X/(Nx*2)-1)**2 - Xsense2 = ((Nx-X)/(Nx*2)-1)**2 - S_MAPS = [np.fft.fftshift(Xsense1),np.fft.fftshift(Xsense2)] - - return S_MAPS - -def SENSE_recon(S1,M1,S2,M2): - - [Nx,Ny,Nz,Nt] = M1.shape - M = np.zeros([Nx,int(2*Ny),Nz,Nt],dtype=complex) - sm1 = np.fft.fftshift(S1)[:,:,0] - sm2 = np.fft.fftshift(S2)[:,:,0] - - for j in range(Ny): - for k in range(Nx): - l1 = M1[k,j,:,:]; a1 = sm1[k,j]; a2 = sm1[k,j+Ny] - l2 = M2[k,j,:,:]; b1 = sm2[k,j]; b2 = sm2[k,j+Ny] - B = (l1*b1 - l2*a1)/(a2*b1 - b2*a1) - A = (l1*b2 - l2*a2)/(a1*b2 - a2*b1) - M[k,j,:,:] = A - M[k,j+Ny,:,:] = B - - - return M - -def SENSE_recon2(S1,M1,S2,M2): - # With matrices as in the original paper! - - [Nx,Ny,Nz,Nt] = M1.shape - M = np.zeros([Nx,int(2*Ny),Nz,Nt],dtype=complex) - sm1 = np.fft.fftshift(S1)[:,:,0] - sm2 = np.fft.fftshift(S2)[:,:,0] - sigma2 = 0.049**2 - sigma2 = 1 - Psi = np.diagflat(np.array([sigma2,sigma2])) # Error matrix Psi - Psi_inv = np.linalg.inv(Psi) - - for j in range(Ny): - for k in range(Nx): - l1 = M1[k,j,:,:]; a1 = sm1[k,j]; a2 = sm1[k,j+Ny] - l2 = M2[k,j,:,:]; b1 = sm2[k,j]; b2 = sm2[k,j+Ny] - S = np.array([[a1,a2],[b1,b2]]) - U = np.linalg.inv((np.transpose(S)*Psi_inv*S))*np.transpose(S)*Psi_inv - a = np.array([l1,l2]) - a_resized = np.resize(a,(2,Nz*Nt)) - v_resized = np.dot(U,a_resized) - v = np.resize(v_resized,(2,Nz,Nt)) - M[k,j,:,:] = v[0,:,:] - M[k,j+Ny,:,:] = v[1,:,:] - - - return M - -def SENSE_METHOD(Seq,R): - ''' - Args: - ITOT: a numpy matrix with the full sampled (3D or 4D) dynamical data - R: the acceleration factor - ''' - - [row,col,dep,numt2] = Seq.shape - Seq_red = {} - SenseMAP = {} - [SenseMAP[0],SenseMAP[1]] = Sensitivity_Map([row,col,dep]) - - col2 = int(np.ceil(col/2)) - - for rs in range(R): - Seq_red[rs] = np.zeros([row,col2,dep,numt2],dtype=complex) - for t in range(numt2): - Kdata_0 = np.fft.fftn(Seq[:,:,:,t]) - Kdata_0 = Kdata_0*SenseMAP[rs] - Kdata_0 = Kdata_0[:,0::R,:] - Seq_red[rs][:,:,:,t] = np.fft.ifftn(Kdata_0) - - Seq_recon = SENSE_recon2(SenseMAP[0],Seq_red[0],SenseMAP[1],Seq_red[1]) - - return Seq_recon - -def undersampling(Mx,My,Mz,options): - - R = options['SENSE']['R'] - - for r in R: - if rank==0: - print('Using Acceleration Factor R = ' + str(r)) - print('applying into x component') - Mx_s = SENSE_METHOD(Mx,r) - if rank==0: - print('applying into y component') - My_s = SENSE_METHOD(My,r) - if rank==0: - print('applying into z component') - Mz_s = SENSE_METHOD(Mz,r) - - return [Mx_s,My_s,Mz_s] - - - diff --git a/codes/__pycache__/CS.cpython-36.pyc b/codes/__pycache__/CS.cpython-36.pyc deleted file mode 100644 index a3acd0f..0000000 Binary files a/codes/__pycache__/CS.cpython-36.pyc and /dev/null differ diff --git a/codes/__pycache__/Graphics.cpython-36.pyc b/codes/__pycache__/Graphics.cpython-36.pyc deleted file mode 100644 index 1cad2a7..0000000 Binary files a/codes/__pycache__/Graphics.cpython-36.pyc and /dev/null differ diff --git a/codes/__pycache__/MRI.cpython-36.pyc b/codes/__pycache__/MRI.cpython-36.pyc deleted file mode 100644 index eaa8bcc..0000000 Binary files a/codes/__pycache__/MRI.cpython-36.pyc and /dev/null differ diff --git a/codes/__pycache__/SENSE.cpython-36.pyc b/codes/__pycache__/SENSE.cpython-36.pyc deleted file mode 100644 index ff25320..0000000 Binary files a/codes/__pycache__/SENSE.cpython-36.pyc and /dev/null differ diff --git a/codes/__pycache__/ktBLAST.cpython-36.pyc b/codes/__pycache__/ktBLAST.cpython-36.pyc deleted file mode 100644 index 7cba033..0000000 Binary files a/codes/__pycache__/ktBLAST.cpython-36.pyc and /dev/null differ diff --git a/codes/ktBLAST.py b/codes/ktBLAST.py deleted file mode 100644 index 7dadb9c..0000000 --- a/codes/ktBLAST.py +++ /dev/null @@ -1,870 +0,0 @@ -import numpy as np -import scipy as sc -from scipy import signal -from mpi4py import MPI -comm = MPI.COMM_WORLD -size = comm.Get_size() -rank = comm.Get_rank() - - - -# kt-BLAST (NO DC TERM) method for reconstruction of undersampled MRI image based on -# l2 minimization. - -def EveryAliased3D2(i,j,k,PP,Nx,Ny,Nz,BB,R): - - ivec = [i,j,k] - Nvec = [Nx,Ny,Nz] - [ktot,ltot] = PP.shape - Ptot = np.zeros([ktot**ltot,ltot]) - PP2 = np.zeros([ktot**ltot,ltot]) - tt = -1 - - for kk in range(Ptot.shape[0]): - nn = int(np.mod(kk,3)) - mm = int(np.mod(np.floor(kk/3),3)) - if np.mod(kk,9)==0: - tt+=1 - - Ptot[kk,0] = PP[tt,0] + ivec[0] - Ptot[kk,1] = PP[mm,1] + ivec[1] - Ptot[kk,2] = PP[nn,2] + ivec[2] - - for kk in range(Ptot.shape[0]): - for ll in range(Ptot.shape[1]): - if Ptot[kk,ll]<0: - Ptot[kk,ll] = Ptot[kk,ll] + Nvec[ll] - if Ptot[kk,ll]>=Nvec[ll]: - Ptot[kk,ll] = Ptot[kk,ll] - Nvec[ll] - - - CC = np.zeros([3,Ptot.shape[0]+1]) - YY = np.array([ [i] , [j], [k] ]) - CC[0,0] = i - CC[1,0] = j - CC[2,0] = k - psel = 0 - - - for l in range(1,Ptot.shape[0]+1): - CC[0,l] = int(Ptot[l-1,0]) - CC[1,l] = int(Ptot[l-1,1]) - CC[2,l] = int(Ptot[l-1,2]) - - - - if CC[0,l]==YY[0,psel] and CC[1,l]==YY[1,psel] and CC[2,l]==YY[2,psel] and BB[int(CC[1,l]),int(CC[2,l]),int(CC[0,l])]!=0: - pass - else: - War = False - for ww in range(psel): - if CC[0,l]==YY[0,ww] and CC[1,l]==YY[1,ww] and CC[2,l]==YY[2,ww] and BB[int(CC[1,l]),int(CC[2,l]),int(CC[0,l])]!=0: - War = True - - if not War: - psel += 1 - CCC = np.array([ [CC[0,l] ] , [CC[1,l]] , [CC[2,l]]]) - YY = np.concatenate( ( YY, CCC ) ,axis=1 ) - - return YY.astype(int) - -def EveryAliased3D(i,j,k,DP,Nx,Ny,Nz,BB,R,SPREAD=None): - - ivec = [i,j,k] - Nvec = [Nx,Ny,Nz] - [ktot,ltot] = DP.shape - DPN = np.zeros([ktot,ltot]) - - if SPREAD is not None: # WITH SPREAD FUNCTIONS FORMALISM - Maux = np.zeros([Ny,Nz,Nx]) - Maux[j,k,i] = 1 - SP2 = SPREAD[::-1,::-1,::-1] - MS = R*sc.signal.convolve(Maux,SP2, mode='same') - ms = np.abs(MS) - Ims = 1*(ms>np.max(ms)*0.405) - Pas = np.where(Ims==1) - PP = np.array(Pas[:]) - PEA = PP[::-1,:] - - for ll in range(PEA.shape[1]): - if PEA[0,ll]>=Nx: - PEA[0,ll] = PEA[0,ll] - Nx - if PEA[1,ll]>=Ny: - PEA[1,ll] = PEA[1,ll] - Ny - if PEA[2,ll]>=Nz: - PEA[2,ll] = PEA[2,ll] - Nz - - - Ntot = PEA.shape[1] - ind = 0 - PEAnew = PEA - - - for ll in range(Ntot): - if BB[PEA[1,ll],PEA[2,ll],PEA[0,ll]]!=0: - PEAnew = np.delete(PEAnew,(ll-ind),axis=1) - ind +=1 - - - return PEA - else: - - for kk in range(DPN.shape[0]): - for l in range(DPN.shape[1]): - DPN[kk,l] = DP[kk,l] + ivec[l] - if DPN[kk,l]<0: - DPN[kk,l] = DPN[kk,l] + Nvec[l] - if DPN[kk,l]>=Nvec[l]: - DPN[kk,l] = DPN[kk,l] - Nvec[l] - - - - CC = np.zeros([3,ktot+1]) - YY = np.array([ [i] , [j], [k] ]) - CC[0,0] = i - CC[1,0] = j - CC[2,0] = k - - - for l in range(1,ktot+1): - CC[0,l] = DPN[l-1,0] - CC[1,l] = DPN[l-1,1] - CC[2,l] = DPN[l-1,2] - - if CC[0,l]!=CC[0,l-1] and CC[1,l]!=CC[1,l-1] and CC[2,l]!=CC[2,l-1] and BB[int(CC[1,l]),int(CC[2,l]),int(CC[0,l])]==0: - CCC = np.array([ [CC[0,l] ] , [CC[1,l]] , [CC[2,l]]]) - YY = np.concatenate( ( YY, CCC ) ,axis=1 ) - - return YY.astype(int) - -def EveryAliased(i,j,DP,Nx,Ny,BB,R,mode): - - if mode==1: # USING GEOMETRICAL ASSUMPTIONS - ivec = [i,j] - Nvec = [Nx,Ny] - DPN = 0*DP - [ktot,ltot] = DP.shape - - for k in range(ktot): - for l in range(ltot): - - DPN[k,l] = DP[k,l] + ivec[l] - - if DPN[k,l]<0: - #DPN[k,l] = ivec[l] - DPN[k,l] = DPN[k,l] + Nvec[l] - - if DPN[k,l]>=Nvec[l]: - #DPN[k,l] = ivec[l] - DPN[k,l] = DPN[k,l] - Nvec[l] - - CC = np.zeros([2,ktot+1]) - YY = np.array([ [i] , [j] ]) - CC[0,0] = i - CC[1,0] = j - - for l in range(1,ktot+1): - CC[0,l] = DPN[l-1,0] - CC[1,l] = DPN[l-1,1] - if CC[0,l]!=CC[0,l-1] and CC[1,l]!=CC[1,l-1] and BB[int(CC[1,l]),int(CC[0,l])]==0: - CCC = np.array([ [CC[0,l] ] , [CC[1,l]] ]) - YY = np.concatenate( ( YY, CCC ) ,axis=1 ) - - return YY.astype(int) - - - - if mode=='spread': # WITH SPREAD FUNCTIONS FORMALISM - Maux = np.zeros([row,numt2]) - Maux[l,k] = 1 - MS = R*ConvolSP(Maux,SPREAD) - ms = np.abs(MS) - Ims = 1*(ms>np.max(ms)*0.405) - Pas = np.where(Ims==1) - PP = np.array(Pas[:]) - return PP[::-1,:] - -def GetSymmetric(M): - [row,numt2] = M.shape - S = np.zeros(M.shape,dtype=complex) - aux = np.zeros([1,row]) - nmid = 0.5*(numt2+1) - for k in range(int(nmid)): - aux = 0.5*( M[:,k] + M[:,numt2-k-1] ) - S[:,k] = aux - S[:,numt2-k-1] = aux - - return S - -def UNDER(A,mode,R,k): - - start = np.mod(k,R) - I1B = np.zeros(A.shape,dtype=complex) - - # Not quite efficient ! better to work with vectors - - if mode=='ky': - for k in range(start,A.shape[0],R): - I1B[k,:,:] = A[k,:,:] - - if mode=='kxky': - for k in range(start,A.shape[0],R): - for l in range(start,A.shape[1],R): - I1B[k,l,:] = A[k,l,:] - - if mode=='kxkykz': - for k in range(start,A.shape[0],R): - for l in range(start,A.shape[2],R): - for r in range(start,A.shape[1],R): - I1B[k,r,l] = A[k,r,l] - - - return I1B - -def FilteringHigh(M,fac): - - if M.ndim==2: - - [row,col] = M.shape - inx = np.linspace(0,col-1,col) - MF = np.zeros(M.shape,dtype=complex) - - for k in range(row): - vecfou = np.fft.fft(M[k,:]) - window = signal.tukey(2*col,fac) - vecfou2 = vecfou*window[col:2*col] - MF[k,:] = np.fft.ifft(vecfou2) - - return MF - - - if M.ndim==3: - [row,col,dep] = M.shape - MF = np.zeros(M.shape,dtype=complex) - inx = np.linspace(0,col-1,col) - for l in range(dep): - for k in range(row): - vecfou = np.fft.fft(M[k,:,l]) - window = signal.tukey(2*col,fac) - vecfou2 = vecfou*window[col:2*col] - MF[k,:,l] = np.fft.ifft(vecfou2) - - - return MF - -def InterpolateM(M,numt2): - - if M.ndim==2: - - [row,numt] = M.shape - MNew = np.zeros([row,numt2],dtype=complex) - xdat = np.linspace(0,numt,numt) - xdat2 = np.linspace(0,numt,numt2) - nstar = int(0.5*(numt2-numt)) - - for t in range(nstar,nstar+numt): - MNew[:,t] = M[:,t-nstar] - - for l in range(row): - ydat = M[l,:] - fdat = sc.interpolate.interp1d(xdat,ydat,kind='cubic') - MNew[l,1:nstar] = fdat(xdat2)[1:nstar] - MNew[l,nstar+numt:numt2] = fdat(xdat2)[nstar+numt:numt2] - - - if M.ndim==3: - [row,col,numt] = M.shape - MNew = np.zeros([row,col,numt2],dtype=complex) - xdat = np.linspace(0,numt,numt) - xdat2 = np.linspace(0,numt,numt2) - nstar = int(0.5*(numt2-numt)) - for t in range(nstar,nstar+numt): - MNew[:,:,t] = M[:,:,t-nstar] - - for c in range(col): - for l in range(row): - ydat = M[l,c,:] - fdat = sc.interpolate.interp1d(xdat,ydat,kind='cubic') - MNew[l,c,1:nstar] = fdat(xdat2)[1:nstar] - MNew[l,c,nstar+numt:numt2] = fdat(xdat2)[nstar+numt:numt2] - - - return MNew - -def KTT(M,scol): - - #Maux = M[:,scol,:] - #Maux = np.fft.ifftshift(Maux,axes=0) - #Maux = np.fft.ifft(Maux,axis=0) - #Maux = np.fft.ifft(Maux,axis=1) - #Maux = np.fft.fftshift(Maux,axes=1) - - # TAO STYLE - Maux = np.zeros(M.shape,dtype=complex) - - for k in range(M.shape[2]): - Maux[:,:,k] = np.fft.ifftshift(M[:,:,k]) - Maux[:,:,k] = np.fft.ifft2(Maux[:,:,k]) - Maux = Maux[:,scol,:] - Maux = np.fft.ifft(Maux,axis=1) - Maux = np.fft.fftshift(Maux,axes=1) - return Maux - -def IKTT(M): - - #Maux = np.fft.ifftshift(M,axes=1) - #Maux = np.fft.ifft(Maux,axis=1) - #Maux = np.fft.fft(Maux,axis=0) - #Maux = np.fft.fftshift(Maux,axes=0) - - # TAO STYLE - Maux = np.fft.ifftshift(M,axes=1) - Maux = np.fft.fft(Maux,axis=1) - - return Maux - -def KTT3D(M,sdep): - - Maux = np.zeros(M.shape,dtype=complex) - for k in range(M.shape[3]): - Maux[:,:,:,k] = np.fft.ifftshift(M[:,:,:,k]) - Maux[:,:,:,k] = np.fft.ifftn(Maux[:,:,:,k]) - Maux = Maux[:,:,sdep,:] - Maux = np.fft.ifft(Maux,axis=2) - Maux = np.fft.fftshift(Maux,axes=2) - return Maux - -def IKTT3D(M): - - Maux = np.fft.ifftshift(M,axes=2) - Maux = np.fft.fft(Maux,axis=2) - - return Maux - -def get_points4D(row,col,dep,numt2,R,mode): - - bb = np.ceil(row/R) - cc = np.ceil(col/R) - aa = np.ceil(numt2/R) - points = R+1 - kmid = int(R/2) - - if mode=='kxky': - PC = [np.ceil(numt2/2),np.ceil(row/2),np.ceil(col/2)] - PP = np.zeros([points,3]) - DP = np.zeros([points-1,3]) - for k in range(points): - - PP[k,0] = numt2-aa*(k) + 1 - PP[k,1] = bb*(k) - PP[k,2] = cc*(k) - - if PP[k,0]>=numt2: - PP[k,0] -= 1 - if PP[k,0]<0: - PP[k,0] += 1 - - if PP[k,1]>=row: - PP[k,1] -= 1 - if PP[k,1]<0: - PP[k,1] += 1 - - if PP[k,2]>=col: - PP[k,2] -= 1 - if PP[k,2]<0: - PP[k,2] += 1 - - if kkmid: - DP[k-1,0] = PP[k,0] - PC[0] - DP[k-1,1] = PP[k,1] - PC[1] - DP[k-1,2] = PP[k,2] - PC[2] - - kmax = int((PP[kmid,0] + PP[kmid-1,0])/2 ) - kmin = int((PP[kmid,0] + PP[kmid+1,0])/2 ) - cmax = int((PP[kmid,1] + PP[kmid-1,1])/2 ) - cmin = int((PP[kmid,1] + PP[kmid+1,1])/2 ) - - - #DP2 = np.zeros([DP.shape[0]**DP.shape[1],DP.shape[1]]) - #DP2[0,0] = DP[0,0]; DP2[0,1] = DP[0,1] ; DP2[0,2] = DP[0,2] - #DP2[1,0] = DP[0,0]; DP2[1,1] = DP[0,1] ; DP2[1,2] = DP[1,2] - #DP2[2,0] = DP[0,0]; DP2[2,1] = DP[1,1] ; DP2[2,2] = DP[0,2] - #DP2[3,0] = DP[0,0]; DP2[3,1] = DP[1,1] ; DP2[3,2] = DP[1,2] - #DP2[4,0] = DP[1,0]; DP2[4,1] = DP[0,1] ; DP2[4,2] = DP[0,2] - #DP2[5,0] = DP[1,0]; DP2[5,1] = DP[0,1] ; DP2[5,2] = DP[1,2] - #DP2[6,0] = DP[1,0]; DP2[6,1] = DP[1,1] ; DP2[6,2] = DP[0,2] - #P2[7,0] = DP[1,0]; DP2[7,1] = DP[1,1] ; DP2[7,2] = DP[1,2] - - - return [kmin,kmax,PP,DP] - - if mode=='ky': - PC = [np.ceil(numt2/2),np.ceil(row/2)] - PP = np.zeros([points,2]) - DP = np.zeros([points-1,2]) - for k in range(points): - PP[k,0] = numt2-(aa-1)*(k) - PP[k,1] = bb*(k) - - if kkmid: - DP[k-1,0] = PP[k,0] - PC[0] - DP[k-1,1] = PP[k,1] - PC[1] - - kmax = int((PP[kmid,0] + PP[kmid-1,0])/2 ) - kmin = int((PP[kmid,0] + PP[kmid+1,0])/2 ) - return [kmin,kmax,PP,DP] - -def SpreadPoint3D(M,R,sdep): - - [row,col,dep,numt2] = M.shape - PS = np.zeros([row,col,dep,numt2],dtype=complex) - inx = 0 - iny = 0 - - for k in range(np.mod(inx,R),row,R): - for ss in range(np.mod(iny,R),col,R): - PS[k,ss,:,:] = 1 - iny = iny + 1 - inx = inx + 1 - - - for k in range(numt2): - PS[:,:,:,k] = np.fft.ifftn(PS[:,:,:,k]) - PS[:,:,:,k] = np.fft.fftshift(PS[:,:,:,k]) - - - SPREAD = PS[:,:,sdep,:] - SPREAD = np.fft.ifft(SPREAD,axis=2) - SPREAD = np.fft.fftshift(SPREAD,axes=2) - - return SPREAD - -def SpreadPoint(M,R,scol): - - [row,col,numt2] = M.shape - PS = np.zeros([row,col,numt2],dtype=complex) - inx = 0 - - for l in range(0,numt2): - for k in range(np.mod(inx,R),row,R): - PS[k,:,l] = 1 - inx = inx + 1 - - #PS = 1*(M!=0) - #PS = 0*M + 1 - #SPREAD = KTT(PS,0) - - for k in range(numt2): - PS[:,:,k] = np.fft.ifft2(PS[:,:,k]) - PS[:,:,k] = np.fft.fftshift(PS[:,:,k]) - - - SPREAD = PS[:,scol,:] - SPREAD = np.fft.ifft(SPREAD,axis=1) - SPREAD = np.fft.fftshift(SPREAD,axes=1) - - return SPREAD - -def ConvolSP(M1,M2): - M2 = M2[::-1,::-1] - M3 = sc.signal.convolve2d(M1,M2, boundary='wrap', mode='same') - return M3 - -def KTBLASTMETHOD_4D_kxky(ITOT,R,mode): - - - ################################################################### - # Training Stage # - ################################################################### - [row,col,dep,numt2] = ITOT.shape - # INPUT PARAMETERS - iteshort = 1 - tetest = int(dep/2) - numt = int(numt2) - Dyy = int(row*0.1) - rmid = int(row/2) - cmid = int(col/2) - - TKdata = np.zeros(ITOT.shape,dtype=complex) - UKdata = np.zeros(ITOT.shape,dtype=complex) - Kdata = np.zeros(ITOT.shape,dtype=complex) - Kdata_NEW0 = np.zeros(ITOT.shape,dtype=complex) - Kdata_NEW = np.zeros(ITOT.shape,dtype=complex) - KTBLAST0 = np.zeros(ITOT.shape,dtype=complex) - KTBLAST = np.zeros(ITOT.shape,dtype=complex) - - - for k in range(numt2): - # THE FULL KSPACE - Kdata[:,:,:,k] = np.fft.fftn(ITOT[:,:,:,k]) - Kdata[:,:,:,k] = np.fft.fftshift(Kdata[:,:,:,k]) - # UNDERSAMPLING STEP AND FILLED WITH ZEROS THE REST - UKdata[:,:,:,k] = UNDER(Kdata[:,:,:,k],mode,R,k) - - # GENERATING THE TRAINING DATA WITH SUBSAMPLING the Center IN KX , KY - for k in range(numt): - TKdata[rmid-Dyy:rmid+Dyy+1,cmid-Dyy:cmid+Dyy+1,:,k] = Kdata[rmid-Dyy:rmid+Dyy+1,cmid-Dyy:cmid+Dyy+1,:,k] - - [kmin,kmax,PP,DP] = get_points4D(row,col,dep,numt2,R,mode) - - if iteshort==1: - print(PP) - print(DP) - print('range of k = ',kmin,kmax) - - SPREAD = SpreadPoint3D(UKdata,R,tetest) - ################################################################### - # RECONSTRUCTION # - ################################################################### - ZE1 = iteshort + (tetest-1)*(iteshort) - ZE2 = (tetest+1)*(iteshort) + (dep)*(1-iteshort) - - for zi in range(ZE1,ZE2): - - if rank==0: - print('4D KTBLAST: R = ' + str(R) + ' and z = ' + str(zi)+'/'+str(dep)) - - ### CONSTRUCT THE REFERENCE M_TRAINING - B2 = KTT3D(TKdata,zi) - B2 = FilteringHigh(B2,0.3) - M2 = 4*np.abs(B2)**2 - #M2 = GetSymmetric(M2) - ### INTERPOLATE IF NUMT0.001) - - if scantype=='0G': - PHASE0[:,:,k] = (gamma*B0*TE+0.01*X)*(np.abs(Sq[:,:,k])>0.001) + 10*varPHASE0 - PHASE1[:,:,k] = (gamma*B0*TE+0.01*X)*(np.abs(Sq[:,:,k])>0.001) + 10*varPHASE0 + np.pi*Sq[:,:,k]/VENC - - if scantype=='-G+G': - PHASE0[:,:,k] = gamma*B0*TE*np.ones([row,col]) + 10*varPHASE0 - np.pi*Sq[:,:,k]/VENC - PHASE1[:,:,k] = gamma*B0*TE*np.ones([row,col]) + 10*varPHASE0 + np.pi*Sq[:,:,k]/VENC - - RHO0[:,:,k] = modulus*np.cos(PHASE0[:,:,k]) + Drho + 1j*modulus*np.sin(PHASE0[:,:,k]) + 1j*Drho2 - RHO1[:,:,k] = modulus*np.cos(PHASE1[:,:,k]) + Drho + 1j*modulus*np.sin(PHASE1[:,:,k]) + 1j*Drho2 - - - if np.ndim(Sq)==4: - [row,col,dep,numt2] = Sq.shape - [X,Y,Z] = np.meshgrid(np.linspace(0,col,col),np.linspace(0,row,row),np.linspace(0,dep,dep)) - - for k in range(numt2): - - if noise: - Drho = np.random.normal(0,0.2,[row,col,dep]) - Drho2 = np.random.normal(0,0.2,[row,col,dep]) - else: - Drho = np.zeros([row,col,dep]) - Drho2 = np.zeros([row,col,dep]) - - varPHASE0 = np.random.randint(-10,11,size=(row,col,dep))*np.pi/180*(np.abs(Sq[:,:,:,k])<0.001) - modulus = 0.5 + 0.5*(np.abs(Sq[:,:,:,k])>0.001) - - if scantype=='0G': - PHASE0[:,:,:,k] = (gamma*B0*TE+0.01*X)*(np.abs(Sq[:,:,:,k])>0.001) + 10*varPHASE0 - PHASE1[:,:,:,k] = (gamma*B0*TE+0.01*X)*(np.abs(Sq[:,:,:,k])>0.001) + 10*varPHASE0 + np.pi*Sq[:,:,:,k]/VENC - - if scantype=='-G+G': - PHASE0[:,:,:,k] = gamma*B0*TE*np.ones([row,col,dep]) + varPHASE0 - np.pi*Sq[:,:,:,k]/VENC - PHASE1[:,:,:,k] = gamma*B0*TE*np.ones([row,col,dep]) + varPHASE0 + np.pi*Sq[:,:,:,k]/VENC - - RHO0[:,:,:,k] = modulus*np.cos(PHASE0[:,:,:,k]) + Drho + 1j*modulus*np.sin(PHASE0[:,:,:,k]) + 1j*Drho2 - RHO1[:,:,:,k] = modulus*np.cos(PHASE1[:,:,:,k]) + Drho + 1j*modulus*np.sin(PHASE1[:,:,:,k]) + 1j*Drho2 - - - - - return [RHO0,RHO1] - -def undersampling(Sqx,Sqy,Sqz,options,savepath): - - R = options['kt-BLAST']['R'] - mode = options['kt-BLAST']['mode'] - transpose = True - - for r in R: - - if rank==0: - print('Using Acceleration Factor R = ' + str(r)) - print('Component x of M0') - - [M0,M1] = GenerateMagnetization(Sqx,options['kt-BLAST']['VENC'],options['kt-BLAST']['noise'],scantype='0G') - if transpose: - M0 = M0.transpose((0,2,1,3)) - M1 = M1.transpose((0,2,1,3)) - - if mode=='ky': - M0_kt = KTBLASTMETHOD_4D_ky(M0,r,mode) - if mode=='kxky': - M0_kt = KTBLASTMETHOD_4D_kxky(M0,r,mode) - - if rank==0: - print('\n Component x of M1') - - if mode=='ky': - M1_kt = KTBLASTMETHOD_4D_ky(M1,r,mode) - if mode=='kxky': - M1_kt = KTBLASTMETHOD_4D_kxky(M1,r,mode) - - - Sqx_kt = phase_contrast(M1_kt,M0_kt,options['kt-BLAST']['VENC'],scantype='0G') - - del M0,M1 - del M0_kt, M1_kt - - [M0,M1] = GenerateMagnetization(Sqy,options['kt-BLAST']['VENC'],options['kt-BLAST']['noise'],scantype='0G') - if transpose: - M0 = M0.transpose((0,2,1,3)) - M1 = M1.transpose((0,2,1,3)) - if rank==0: - print('\n Component y of M0') - - if mode=='ky': - M0_kt = KTBLASTMETHOD_4D_ky(M0,r,mode) - if mode=='kxky': - M0_kt = KTBLASTMETHOD_4D_kxky(M0,r,mode) - - - if rank==0: - print('\n Component y of M1') - - if mode=='ky': - M1_kt = KTBLASTMETHOD_4D_ky(M1,r,mode) - if mode=='kxky': - M1_kt = KTBLASTMETHOD_4D_kxky(M1,r,mode) - - Sqy_kt = phase_contrast(M1_kt,M0_kt,options['kt-BLAST']['VENC'],scantype='0G') - - del M0,M1 - del M0_kt, M1_kt - - [M0,M1] = GenerateMagnetization(Sqz,options['kt-BLAST']['VENC'],options['kt-BLAST']['noise'],scantype='0G') - if transpose: - M0 = M0.transpose((0,2,1,3)) - M1 = M1.transpose((0,2,1,3)) - if rank==0: - print('\n Component z of M0') - - if mode=='ky': - M0_kt = KTBLASTMETHOD_4D_ky(M0,r,mode) - - if mode=='kxky': - M0_kt = KTBLASTMETHOD_4D_kxky(M0,r,mode) - - if rank==0: - print('\n Component z of M1') - - if mode=='ky': - M1_kt = KTBLASTMETHOD_4D_ky(M1,r,mode) - if mode=='kxky': - M1_kt = KTBLASTMETHOD_4D_kxky(M1,r,mode) - - if rank==0: - print(' ') - - Sqz_kt = phase_contrast(M1_kt,M0_kt,options['kt-BLAST']['VENC'],scantype='0G') - - - if transpose: - Sqx_kt = Sqx_kt.transpose((0,2,1,3)) - Sqy_kt = Sqy_kt.transpose((0,2,1,3)) - Sqz_kt = Sqz_kt.transpose((0,2,1,3)) - - - if options['kt-BLAST']['save']: - if rank==0: - print('saving the sequences in ' + savepath) - seqname = options['kt-BLAST']['name'] +'_R' + str(r) + '.npz' - print('sequence name: ' + seqname) - np.savez_compressed( savepath + seqname, x=Sqx_kt, y=Sqy_kt,z=Sqz_kt) - - del Sqx_kt,Sqy_kt,Sqz_kt - - - - diff --git a/codes/mesh_generator.py b/codes/mesh_generator.py deleted file mode 100755 index 5207cc8..0000000 --- a/codes/mesh_generator.py +++ /dev/null @@ -1,71 +0,0 @@ -from dolfin import * - - -mesh_in = '/home/yeye/Desktop/leomesh.xml' -mesh_out = '/home/yeye/Desktop/aorta.h5' -mesh = Mesh(mesh_in) - -hdf = HDF5File(mesh.mpi_comm(), mesh_out, 'w') -boundaries = MeshFunction('size_t', mesh,2) - -marked = 1 -testmesh = 0 - -hdf.write(mesh, '/mesh') - -if marked==1: - - class Inlet(SubDomain): - def inside(self, x, on_boundary): - return on_boundary and between(x[0],(0.1975,0.1989)) and between(x[2],(0.07,0.1)) - - class Outlet(SubDomain): - def inside(self, x, on_boundary): - return on_boundary and between(x[0],(0.1975,0.1989)) and between(x[2],(0,0.04)) - - class Walls(SubDomain): - def inside(self, x, on_boundary): - return on_boundary - - - outlet = Outlet() - inlet = Inlet() - walls = Walls() - - boundaries.set_all(0) - walls.mark(boundaries,1) - outlet.mark(boundaries,3) - inlet.mark(boundaries,2) - - - - -hdf.write(boundaries, '/boundaries') -hdf.close() - - - -if testmesh: - print('Testing Mesh...') - meshname = mesh_out - pathtoB = '/home/yeye/Desktop/boundaries.xdmf' - mesh = Mesh() - hdf = HDF5File(mesh.mpi_comm(), meshname , 'r') - hdf.read(mesh, '/mesh', False) - boundaries = MeshFunction('size_t', mesh , mesh.topology().dim() - 1) - hdf.read(boundaries, '/boundaries') - # To save the boundaries information - XDMFFile(pathtoB).write(boundaries) - print('Boundary info printed in ' + pathtoB) - - - - - - - - - - - - diff --git a/codes/monolithic.py b/codes/monolithic.py deleted file mode 100644 index 8576d8a..0000000 --- a/codes/monolithic.py +++ /dev/null @@ -1,229 +0,0 @@ -from dolfin import * -import numpy as np -from common import inout -from mpi4py import MPI -import sys -import os -# -# NAVIER STOKES PROBLEM IN THE AORTA with a MONOLITHIC SOLVER -# THIS SCRIPT INCLUDE THE 0-WINDKESSEL BOUNDARY CONDITION -# -# Written by Jeremias Garay L: j.e.garay.labra@rug.nl -# - -parameters["std_out_all_processes"] = False - -def solv_NavierStokes(options): - - # Assign physical parameters - rho = Constant(options['density']) - mu = Constant(options['dynamic_viscosity']) - otheta = Constant(1-options['theta']) - theta = Constant(options['theta']) - Tf = options['Tf'] - dt = options['dt'] - - # CREATING THE FILES - xdmf_u = XDMFFile(options['savepath']+'u.xdmf') - xdmf_p = XDMFFile(options['savepath']+'p.xdmf') - # LOADING THE MESH - mesh = Mesh() - hdf = HDF5File(mesh.mpi_comm(), options['mesh_path'] , 'r') - hdf.read(mesh, '/mesh', False) - bnds = MeshFunction('size_t', mesh , mesh.topology().dim() - 1) - hdf.read(bnds, '/boundaries') - - # DEFINE FUNCTION SPACES - if options['fem_space'] == 'p2p1': - V = VectorElement('Lagrange', mesh.ufl_cell(), 2) - Q = FiniteElement('Lagrange', mesh.ufl_cell(), 1) - pspg = False - elif options['fem_space'] == 'p1p1': - V = VectorElement('Lagrange', mesh.ufl_cell(), 1) - Q = FiniteElement('Lagrange', mesh.ufl_cell(), 1) - pspg = True - TH = V * Q - W = FunctionSpace(mesh, TH) - n = FacetNormal(mesh) - ds = Measure("ds", subdomain_data=bnds) - v, q = TestFunctions(W) - # Boundary Conditions - BCS = [] - bc = options['boundary_conditions'] - # For Windkessel implementation - pii0 = {} - pii = {} - press = {} - alpha = {} - beta = {} - gamma = {} - Windkvar = {} - windkessel = False - - u, p = TrialFunctions(W) - w = Function(W) - h = CellDiameter(mesh) - - u0, p0 = w.split() - u_ = theta*u + otheta*u0 - p_ = theta*p + otheta*p0 - # The variational formulation - # Mass matrix - F = ( - (rho/dt)*dot(u - u0, v)*dx - + mu*inner(grad(u_), grad(v))*dx - - p_*div(v)*dx + q*div(u)*dx - + rho*dot(grad(u_)*u0, v)*dx - ) - - - for nbc in range(len(bc)): - nid = bc[nbc]['id'] - if bc[nbc]['type'] == 'dirichlet': - if rank==0: - print('Adding Dirichlet BC at boundary ', nid) - - val = bc[nbc]['value'] - if isinstance(val, (int, float)): - val = Constant(val) - BCS.append(DirichletBC(W.sub(0), val, bnds, nid)) - else: - params = bc[nbc]['parameters'] if 'parameters' in bc[nbc] else dict() - inflow = Expression(val, degree=3, **params) - BCS.append(DirichletBC(W.sub(0), inflow, bnds, nid)) - - if bc[nbc]['type'] == 'windkessel': - windkessel = True - if rank==0: - print('Adding Windkessel BC at boundary ',nid) - [R_p,R_d,C] = bc[nbc]['value'] - # coeficients - alpha[nid] = R_d*C/(R_d*C + dt) - beta[nid] = R_d*(1-alpha[nid]) - gamma[nid] = R_p + beta[nid] - # dynamical terms - if C ==0: - print('no capacitance C for boundary',nid) - press[nid] = Constant(0) - else: - pii0[nid] = Constant(bc[nbc]['p0'][0]*bc[nbc]['p0'][1]) - pii[nid] = Constant(bc[nbc]['p0'][0]*bc[nbc]['p0'][1]) - press[nid] = Constant(bc[nbc]['p0'][0]*bc[nbc]['p0'][1]) - - Windkvar[nid] = press[nid]*dot(v,n)*ds(nid) - - - # Stabilization Terms - if options['stabilization']['temam']: - if rank==0: - print('Adding Temam stabilization term') - F += 0.5*rho*div(u0)*dot(u_, v)*dx - - if pspg: - if rank==0: - print('Adding PSPG stabilization term') - eps = Constant(options['stabilization']['eps']) - F += eps/mu*h**2*inner(grad(p_), grad(q))*dx - - if options['stabilization']['forced_normal']['enabled']: - gparam = options['stabilization']['forced_normal']['gamma'] - for nid in options['stabilization']['forced_normal']['boundaries']: - if rank==0: - print('Forcing normal velocity in border ', nid) - ut = u - n*dot(u0,n) - vt = v - n*dot(v,n) - F += gparam*dot(ut,vt)*ds(nid) - - - if len(options['stabilization']['backflow_boundaries'])>0: - def abs_n(x): - return 0.5*(x - abs(x)) - for nk in options['stabilization']['backflow_boundaries']: - if rank==0: - print('adding backflow stabilization in border number:' + str(nk)) - F -= 0.5*rho*abs_n(dot(u0, n))*dot(u_, v)*ds(nk) - - - if windkessel: - for nid in Windkvar.keys(): - F += Windkvar[nid] - - a = lhs(F) - L = rhs(F) - - # The static part of the matrix - A = assemble(a) - u, p = w.split() - uwrite = Function(W.sub(0).collapse()) - pwrite = Function(W.sub(1).collapse()) - uwrite.rename('velocity', 'u') - pwrite.rename('pressure', 'p') - u.rename('velocity', 'u') - p.rename('pressure', 'p') - ind = 0 - t = dt - checkcicle = int(options['checkpoint_dt']/options['dt']) - writecicle = int(options['xdmf_dt']/options['dt']) - - while t<=Tf+dt: - # To solve - assemble(a, tensor=A) - b = assemble(L) - [bcs.apply(A, b) for bcs in BCS] - print('solving for t = ' + str(np.round(t,4))) - solve(A, w.vector(), b ) - ind += 1 - - if options['write_xdmf']: - if np.mod(ind,writecicle)<0.1 or ind==1: - xdmf_u.write(u, t) - xdmf_p.write(p, t) - - assign(uwrite, w.sub(0)) - assign(pwrite, w.sub(1)) - - if np.mod(ind,checkcicle)<0.1 or ind==1: - if options['write_checkpoint']: - checkpath = options['savepath'] +'checkpoint/{i}/'.format(i=ind) - comm = uwrite.function_space().mesh().mpi_comm() - inout.write_HDF5_data(comm, checkpath + '/u.h5', uwrite, '/u', t=t) - inout.write_HDF5_data(comm, checkpath + '/p.h5', pwrite, '/p', t=t) - - t += dt - inflow.t = t - assign(u0, w.sub(0)) - - if windkessel: - for nid in Windkvar.keys(): - qq = assemble(dot(u0,n)*ds(nid)) - pii0[nid].assign(pii[nid]) - pii[nid].assign(Constant(alpha[nid]*float(pii0[nid]) + beta[nid]*qq)) - pp = Constant(gamma[nid]*qq + alpha[nid]*float(pii0[nid])) - press[nid].assign(pp) - - - - - - -if __name__ == '__main__': - - comm = MPI.COMM_WORLD - size = comm.Get_size() - rank = comm.Get_rank() - - if len(sys.argv) > 1: - if os.path.exists(sys.argv[1]): - inputfile = sys.argv[1] - if rank==0: - print('Found input file ' + inputfile) - else: - raise Exception('Command line arg given but input file does not exist:' - ' {}'.format(sys.argv[1])) - else: - raise Exception('An input file is required as argument!') - - - options = inout.read_parameters(inputfile) - solv_NavierStokes(options) - diff --git a/input/Graphics_input.yaml b/input/Graphics_input.yaml deleted file mode 100755 index 1d99157..0000000 --- a/input/Graphics_input.yaml +++ /dev/null @@ -1,82 +0,0 @@ -################################################# -# -# Input file for Graphics -# -################################################# -dP: - apply: false - data: '/home/yeye/Desktop/dP/results/11AoCoPhantomRest2/R1/testBF/' - R: [1] - #catheter: '/home/yeye/Desktop/PhD/MEDICAL_DATA/DatosSEPT2019/Phantom_catheter/cat9mm_Stress/pressure_drop.npz' - catheter: 'None' - #factors: [130,90,1.58,-0.03,0] # 9mm Rest [x0,x1,a,b,theta] - #factors: [52,35,0.9,0.05,0] # 9mm Stress [x0,x1,a,b,theta] - #factors: [162,90,1.67,0.03,0] # 11mm Rest [x0,x1,a,b,theta] - mode: 'cib' - estim: ['PPE','STE','COR'] - save: True - name: '/home/yeye/Desktop/11AoCoPhantomRest2_BF.png' - -Histograms_meshes: - apply: false - outpath: '/home/yeye/Desktop/' - colors: - 0: 'orangered' - 1: 'lime' - 2: 'dodgerblue' - meshnames: - 0: 'coaortaH1' - 1: 'coaortaRH2' - 2: 'coaortaRH3' - paths: - 0: '/home/yeye/Desktop/PhD/AORTA/MESH/coaorta/H1/' - 1: '/home/yeye/Desktop/PhD/AORTA/MESH/coaorta/H2/' - 2: '/home/yeye/Desktop/PhD/AORTA/MESH/coaorta/H3/' - -HistCorrector: - apply: false - errors: '/home/yeye/Desktop/testH/H2/errors.dat' - outpath: '/home/yeye/Desktop/h2_160.png' - -Histograms_checkpoint: - apply: false - path: '/home/yeye/Desktop/Corrector_2019/mono2/dt10ms_SUPGcon/' - title: '$dt = 10 \ ms$' - -Error-curves: - apply: false - folder: '/home/yeye/Desktop/Corrector_2019/Perturbation/' - type: ['norm2_m'] - subfolders: ['SNRinfV120','SNR10V120','SNRinfV80','SNR10V80'] - labels: ['SNRinfV120','SNR10V120','SNRinfV80','SNR10V80'] - #subfolders: ['BA_p2p1','BAA','BBA'] - #labels: ['BA_{P2P1}','BAA','BBA'] - colors: ['indigo','limegreen','dodgerblue','orangered','yellow'] - styles: ['-','-','-','-','-','-','-','-'] - title: '$leo2.0$' - outpath: '/home/yeye/Desktop/Corrector_2019/Perturbation/' - -Pressure_drops: - apply: false - folder: '/home/yeye/Desktop/Corrector_2019/Perturbation/STE/' - convertor: -1333.22387415 - #convertor: -133.322 - catheter: false - subfolders: ['SNRinfV120','SNR15V120','SNRinfV80','SNR15V80','dt10ms'] - labels: ['SNRinfV120','SNR15V120','SNRinfV80','SNR15V80','ref'] - colors: ['indigo','limegreen','dodgerblue','orangered','yellow'] - styles: ['-','-','-','-','-','-','-','-','-','-'] - title: '$STE \ leo2.0$' - outpath: '/home/yeye/Desktop/Corrector_2019/Perturbation/STE/' - -l2_comp: - apply: True - colors: ['dodgerblue','orangered','limegreen'] - div: false - aliasing: false - folder: '/home/yeye/Desktop/Poiseuille/curves/' - subfolder: ['SNR10V80'] - name: 'SNR10V120_gain' - mode: - type: 'gain_compressed' - comp: '/home/yeye/Desktop/Poiseuille/curves/SNR10V120/' \ No newline at end of file diff --git a/input/PostCheck_input.yaml b/input/PostCheck_input.yaml deleted file mode 100755 index c44d2c3..0000000 --- a/input/PostCheck_input.yaml +++ /dev/null @@ -1,123 +0,0 @@ -################################################# -# -# Input file for the checkpoint postprocessing -# -################################################# -Algebra: - apply: false - mode: '+' - VENC: 97 - outname: 'ucor' - outpath: '/home/yeye/Desktop/Corrector_2019/Perturbation/gain/SNR10V60/' - checkpoint: true - #meshpath: '/home/yeye/Desktop/Poiseuille/Meas_leo/poiseuille.h5' - meshpath: '/home/yeye/Desktop/Corrector_2019/Meshes/leoH3_2.0.h5' - v1: - name: 'u' - path: '/home/yeye/Desktop/Corrector_2019/Perturbation/Meas/SNR10V60/' - v2: - name: 'w_COR_impl_stan' - path: '/home/yeye/Desktop/Corrector_2019/Perturbation/COR/SNR10V60/' - -Colormap: - apply: false - #mode: ['u','w','w/u','div(u)','div(u)/u'] - mode: ['dot'] - outpath: '/home/yeye/Desktop/Poiseuille/H5/leo2mm/SNRinfV80/' - meshpath: '/home/yeye/Desktop/Poiseuille/Meas_leo/poiseuille.h5' - upath: '/home/yeye/Desktop/Poiseuille/Meas_leo/SNRinfV80/' - wpath: '/home/yeye/Desktop/Poiseuille/Corrector/leo2mm/SNRinfV80/' - uname: 'u' - wname: 'w_COR_impl_stan' - -Velocity: - apply: false - meshpath: '/home/yeye/Desktop/Corrector_2019/Poiseuille/poiseuille.h5' - checkpoint: '/home/yeye/Desktop/Corrector_2019/Poiseuille/COR/SNR15V120/' - undersampling: 1 - dt: 0.03 - filename: 'w_COR_impl_stan' - -Estim_Pressure: - apply: false - meshpath: '/home/yeye/Desktop/Corrector_2019/Meshes/9mmRest2.0.h5' - filename: 'p_COR_impl_stan' - checkpath: '/home/yeye/Desktop/Corrector_2019/AoCo/9mm_pspg/' - method: spheres # slices, boundaries, spheres - dt: 0.03072 - boundaries: [2,3] - spheres: - - center: [0.0980092, 0.0909768, 0.0802258] # 9mm - #- center: [0.110266, 0.086805, 0.0794744] # 11mm - #- center : [0.0940797, 0.0766444, 0.080433] #13mm - #- center: [0.0870168, 0.0901715, 0.0883529] # Normal - radius: 0.005 - Npts: 32 - - center: [0.0980092, 0.135047, 0.0252659] # 9mm - #- center: [0.110266, 0.133002, 0.0263899] # 11mm - #- center : [0.0940573, 0.123321, 0.0260755] # 13 mm - #- center: [0.0869949, 0.12838, 0.0269889] # Normal - radius: 0.005 - Npts: 32 - -Outlet_Wind: - apply: false - mesh_path: '/home/yeye/Desktop/aorta/mesh/aorta_coarse_marked.h5' - checkpoint: '/home/yeye/Desktop/aorta/results/mono/' - filename: ['u','p'] - bnds: [3,4,5,6] - out_path: '/home/yeye/Desktop/aorta/results/mono/' - -Error-curves: - apply: true - dt: 0.03 - VENC: 204 #194/129/113/97 for aorta(120,80,70,60)/ 204/136 - type: ['utrue-uobs'] - #type: ['l2_comp','div'] - meshpath: '/home/yeye/Desktop/Poiseuille/Meas_leo/poiseuille.h5' - #meshpath: '/home/yeye/Desktop/Corrector_2019/Meshes/leoH3_2.0.h5' - true_check: '/home/yeye/Desktop/Poiseuille/Meas_leo/SNRinfV120/' - true_name: 'u' - ref_check: '/home/yeye/Desktop/Poiseuille/Meas_leo/SNR10V120/' - ref_name: 'u' - undersampling: 1 - com_check: '/home/yeye/Desktop/Poiseuille/Corrector/leo2mm/SNR10V120/' - com_name: 'w_COR_impl_stan' - outpath: '/home/yeye/Desktop/Poiseuille/curves/SNR10V120/' - -Histograms: - apply: false - type: ['normal','grad'] - meshpath: '/home/yeye/Desktop/Corrector_2019/meshes/11AoCoPhantomRest1.4.h5' - checkpath: '/home/yeye/Desktop/Corrector_2019/I/corrector/11mmRest1.4/' - field_name: 'w_COR_impl_stan' - outpath: '/home/yeye/Desktop/Corrector_2019/I/corrector/11mmRest1.4/' - -Temporal-Average: - apply: false - meshpath: '/home/yeye/Desktop/PhD/AORTA/MESH/coaorta/H1/coaortaH1.h5' - original_check: '/home/yeye/Desktop/First/H1/dt5ms/coaorta/' - xdmf: false - dt: 0.005 - subsampling_rate: 6 - out_check: '/home/yeye/Desktop/First/H1/test/' - -Temporal-Interpolation: - apply: false - meshpath: '/home/yeye/Desktop/Corrector_2019/meshes/11AoCoPhantomRest1.4.h5' - original_check: '/home/yeye/Desktop/Corrector_2019/I/meas/11mmRest1.4/R1/' - xdmf: true - dt: 0.03072 - dt_new: 0.015 - out_check: '/home/yeye/Desktop/Corrector_2019/I/meas/11mmRest1.4/dt15ms/' - -Perturbation: - apply: false - undersampling: 1 - type: - SNR: 10 # dB signal to noise ratio - coil: true - phase_contrast: 80 # venc % respect u max - meshpath: '/home/yeye/Desktop/Poiseuille/Meas_leo/poiseuille.h5' - checkpath: '/home/yeye/Desktop/Poiseuille/Meas_leo/original/' - xdmf: true \ No newline at end of file diff --git a/input/aorta.yaml b/input/aorta.yaml deleted file mode 100755 index 71314ed..0000000 --- a/input/aorta.yaml +++ /dev/null @@ -1,138 +0,0 @@ -# Set of default parameters for steady Navier-Stokes -mesh: '/home/yeye/Desktop/PhD/AORTA/MESH/coaorta/H1/coaortaH1.h5' -density: 1.119 -dynamic_viscosity: 0.0483 -stokes: False - -io: - write_hdf5: True - write_hdf5_timeseries: False - write_xdmf: True - write_path: '/home/yeye/Desktop/coaorta' - restart: - path: '' # './projects/nse_coa3d/results/test_restart2/' - time: 0 - write_checkpoints: true - write_velocity: 'update' # tentative - log: False - - -boundary_conditions: - - id: 2 - type: 'dirichlet' - value: ['0','0','U*sin(DOLFIN_PI*t/Th)*(t<=Th) + (Th