diff --git a/codes/CS.py b/codes/CS.py index 1d0d2bd..85a1a1c 100644 --- a/codes/CS.py +++ b/codes/CS.py @@ -14,487 +14,403 @@ rank = comm.Get_rank() # 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=norm: + P = pdf(ks, k+1, klo, q) + if np.sum(P) >= norm: break - - P = np.fft.fftshift(P) - - - #if np.mod(n,2)!=0: - # P = np.concatenate(([1],P),axis=None) - - + 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 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): -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) + 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: + 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 != 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 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 - + 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)) - + 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 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])) - + + 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 + 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 - + + 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 - + 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) +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 0) + s = s + 1*(s < tt) + ss = ss/s return ss*X -def CSMETHOD(ITOT,R): - - ''' Compressed Function. +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 + 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) + MASK = GeneratePattern(ITOT.shape, R) + CS1 = np.zeros(ITOT.shape, dtype=complex) + + nit = 0 + nit_tot = (numt2-1)/20 + + if np.ndim(ITOT) == 3: - 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) + 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) + 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) + 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 - + 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: - + + 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)) - #if np.mod(t,nit_tot)<1: - # sys.stdout.write('\r') - # # Progress bar - # if numt2==3: - # sys.stdout.write("{4d-CS} [%-6s] %d%%" % ('=='*nit, 100*t/numt2)) - # else: - # sys.stdout.write("{4d-CS} [%-40s] %d%%" % ('=='*nit, 100*t/numt2)) - # sys.stdout.flush() - # nit = nit +1 - + 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 + 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 + 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]) + 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) + 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) + 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) + 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 - + 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 - - - - + + CS1[:, :, :, t] = img + return CS1 -def CSMETHOD_SENSE(ITOT,R,R_SENSE): - ''' Compressed Function. +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 + 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]) + [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) - 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)) + 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 + 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 + 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]) + 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) + 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) + 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) + 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 + 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 + ITOTCS[rs][:, :, :, t] = img - - return [ITOTCS[0],ITOTCS[1]] - -def CSMETHOD_peaksystole(ITOT,R,tstar): - - ''' Compressed Function. + return [ITOTCS[0], ITOTCS[1]] - Args: - ITOT: a numpy matrix with the full sampled (3D or 4D) dynamical data - R: the acceleration factor - tstar: the time when the flux in the inlet it's maximum - ''' - # Method parameters - ninner = 5 - nbreg = 10 - lmbda = 4 - mu = 20 - gam = 1 - - - - [row,col,dep,numt2] = ITOT.shape - - MASK = GeneratePattern(ITOT.shape,R) - CS1 = np.zeros([row,col,dep],dtype=complex) - - for t in range(tstar,tstar+1): - Kdata = np.fft.fftn(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,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[:,:,:] = img - - - return CS1 - -def phase_contrast(M1,M0,VENC,scantype='0G'): - param = 1 - if scantype=='-G+G': - param = 0.5 +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'): + +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 @@ -504,202 +420,178 @@ def GenerateMagnetization(Sq,VENC,noise,scantype='0G'): 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)) + 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 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]) - - if np.ndim(Sq)==4: + 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]) + 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) + 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 == '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 + 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 + 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): - return [RHO0,RHO1] + R = options['cs']['R'] -def undersampling(Sqx,Sqy,Sqz,options,savepath): - - R = options['cs']['R'] - for r in R: - - if rank==0: + + 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']) - + + [M0, M1] = GenerateMagnetization( + Sqx, options['cs']['VENC'], options['cs']['noise']) print('\n Component x of M0') - M0_cs = CSMETHOD(M0,r) + M0_cs = CSMETHOD(M0, r) print('\n Component x of M1') - M1_cs = CSMETHOD(M1,r) + M1_cs = CSMETHOD(M1, r) - Sqx_cs = phase_contrast(M1_cs,M0_cs,options['cs']['VENC']) - del M0,M1 + 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']) - + + [M0, M1] = GenerateMagnetization( + Sqy, options['cs']['VENC'], options['cs']['noise']) print('\n Component y of M0') - M0_cs = CSMETHOD(M0,r) + M0_cs = CSMETHOD(M0, r) print('\n Component y of M1') - M1_cs = CSMETHOD(M1,r) - + M1_cs = CSMETHOD(M1, r) - Sqy_cs = phase_contrast(M1_cs,M0_cs,options['cs']['VENC']) - - del M0,M1 + 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: + + [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: + M0_cs = CSMETHOD(M0, r) + if rank == 0: print('\n Component z of M1') - M1_cs = CSMETHOD(M1,r) - if rank==0: + M1_cs = CSMETHOD(M1, r) + if rank == 0: print(' ') - - Sqz_cs = phase_contrast(M1_cs,M0_cs,options['cs']['VENC']) - - - if rank==0: + + 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_peakpv(Sqx,Sqy,Sqz,options,R): - - Sqx_cs = {} - Sqy_cs = {} - Sqz_cs = {} - [Mx0,Mx1] = GenerateMagnetization(Sqx,options['cs']['VENC'],options['cs']['noise'],scantype='0G') - [My0,My1] = GenerateMagnetization(Sqy,options['cs']['VENC'],options['cs']['noise'],scantype='0G') - [Mz0,Mz1] = GenerateMagnetization(Sqz,options['cs']['VENC'],options['cs']['noise'],scantype='0G') - + 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) - Mx0_cs = CSMETHOD(Mx0,R) - Mx1_cs = CSMETHOD(Mx1,R) - My0_cs = CSMETHOD(My0,R) - My1_cs = CSMETHOD(My1,R) - Mz0_cs = CSMETHOD(Mz0,R) - Mz1_cs = CSMETHOD(Mz1,R) + del Sqx_cs, Sqy_cs, Sqz_cs - Sqx_cs = phase_contrast(Mx1_cs,Mx0_cs,options['cs']['VENC'],scantype='0G') - Sqy_cs = phase_contrast(My1_cs,My0_cs,options['cs']['VENC'],scantype='0G') - Sqz_cs = phase_contrast(Mz1_cs,Mz0_cs,options['cs']['VENC'],scantype='0G') - - - - return [Sqx_cs,Sqy_cs,Sqz_cs] - -def undersampling_short(Mx,My,Mz,options): - - R = options['cs']['R'] - savepath = options['cs']['savepath'] +def undersampling_short(Mx, My, Mz, options): + R = options['cs']['R'] + savepath = options['cs']['savepath'] - R_SENSE = 1 + R_SENSE = 1 if 'R_SENSE' in options['cs']: - R_SENSE = options['cs']['R_SENSE'][0] + R_SENSE = options['cs']['R_SENSE'][0] - for r in R: - if rank==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: + 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) + 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: + 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 + 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 index 91d2685..89bf6ed 100644 --- a/codes/Graphics.py +++ b/codes/Graphics.py @@ -586,25 +586,15 @@ def CenterComparison(A,B,C,R,center): def VelocityChannel(M,repeat,recorder): [row,col,numt2] = M.shape - - [X,Y] = np.meshgrid(np.linspace(0,col,col),np.linspace(0,row,row)) - plt.ion() - rown = row*6/row coln = col*6/row - fig = plt.figure()#figsize=(4, 6) , dpi=200) ax = fig.add_subplot(111, projection='3d') - - for rr in range(-1,repeat): - for t in range(numt2): - V = M[:,:,t] - #ax.plot_surface(X, Y, V, cmap=plt.cm.magma, vmin=-30, vmax=50, linewidth=0, antialiased=False) #ax.plot_wireframe(X, Y, V,rcount=20,ccount=20,linewidth=0.5) ax.plot_surface(X, Y, V,cmap='magma',vmin=-30, vmax=100, linewidth=0, antialiased=False) @@ -615,7 +605,6 @@ def VelocityChannel(M,repeat,recorder): ax.set_xlabel('$x$') ax.set_ylabel('$y$') ax.set_zlabel('$v(r)$') - ax.set_title('phase ' + str(t)) plt.pause(0.001) plt.draw() @@ -634,116 +623,7 @@ def PlotTri(tri,pos,p): #ax.plot_trisurf(pos[:,0], pos[:,1], pos[:,2], triangles=tri.simplices, cmap=plt.cm.Spectral,linewidth=0.1,edgecolors='k') #ax.tripcolor(pos[:,0], pos[:,1], pos[:,2], triangles=tri.simplices, facecolors=p2.T, edgecolor='black') plt.show() - -def PlotPressureDrop(mode): - - barye2mmHg = 1/1333.22387415 - - CT = np.loadtxt('Pressure/DROPS/pd_NS_coarse.txt') - CT2 = np.loadtxt('Pressure/DROPS/pd_NS_coarse2.txt') - - ref_ao = np.loadtxt('Pressure/DROPS2/refPPE1_coarse.txt') - ref = np.loadtxt('Pressure/DROPS2/refPPE1_coarse_leo1.txt') - CS2 = np.loadtxt('Pressure/DROPS2/CSPPE2_coarse_leo1.txt') - - - PPE_CT = np.loadtxt('Pressure/DROPS/pd_PPE_NS_coarse.txt') - PPE_CT_leo = np.loadtxt('Pressure/DROPS/pd_PPE_NS_coarse_leo1.txt') - - - STE_CT = np.loadtxt('Pressure/DROPS/test_STE_R1_coarse.txt') - test_STE_CT = np.loadtxt('Pressure/DROPS/test2_STE_R1_coarse_leo1.txt') - - STEint_CT = np.loadtxt('Pressure/DROPS/test_STEint_R1_coarse.txt') - test_STEint_CT = np.loadtxt('Pressure/DROPS/test2_STEint_R1_coarse_leo1.txt') - - - # UNDERSAMPLING - PPE_R1_leo1 = np.loadtxt('Pressure/DROPS/test_PPE_R1_coarse_leo1.txt') - - PPE_KT_R2_leo1 = np.loadtxt('Pressure/DROPS/KT_PPE_R2_coarse_leo1.txt') - PPE_KT_R4_leo1 = np.loadtxt('Pressure/DROPS/KT_PPE_R4_coarse_leo1.txt') - PPE_KT_R8_leo1 = np.loadtxt('Pressure/DROPS/KT_PPE_R8_coarse_leo1.txt') - PPE_KT_R12_leo1 = np.loadtxt('Pressure/DROPS/KT_PPE_R12_coarse_leo1.txt') - PPE_KT_R20_leo1 = np.loadtxt('Pressure/DROPS/KT_PPE_R20_coarse_leo1.txt') - - PPE_CS_R2_leo1 = np.loadtxt('Pressure/DROPS/CS_PPE_R2_coarse_leo1.txt') - PPE_CS_R4_leo1 = np.loadtxt('Pressure/DROPS/CS_PPE_R4_coarse_leo1.txt') - PPE_CS_R8_leo1 = np.loadtxt('Pressure/DROPS/CS_PPE_R8_coarse_leo1.txt') - PPE_CS_R12_leo1 = np.loadtxt('Pressure/DROPS/CS_PPE_R12_coarse_leo1.txt') - PPE_CS_R20_leo1 = np.loadtxt('Pressure/DROPS/CS_PPE_R20_coarse_leo1.txt') - #CT_fine = np.loadtxt('Pressure/DROPS/pd_NS_fine.txt') - #CT_fine_T = np.loadtxt('Pressure/DROPS/pd_NS_fine_T.txt') - #PPE_CT_fine = np.loadtxt('Pressure/DROPS/pd_PPE_NS_fine.txt') - #PPE_CT_fine_leo3 = np.loadtxt('Pressure/DROPS/pd_PPE_NS_fine_leo3.txt') - - - tvec2 = np.linspace(0, 2.5 ,PPE_CT_leo.size) - tvec = np.linspace(tvec2[1]*0.5 ,2.5+tvec2[1]*0.5 , CT.size) - #tvec_T = np.linspace(0,2.5,CT_fine_T.size) - - fig = plt.figure() - - if mode=='KT': - plt.plot(tvec,CT,'-k',linewidth=2,label='$ref$') - plt.plot(tvec2,PPE_R1_leo1,'xkcd:red',linewidth=2,linestyle='-' , label='$R = 1$') - plt.plot(tvec2,PPE_KT_R2_leo1,'xkcd:blue',linewidth=2,linestyle='-' ,label='$R = 2$') - plt.plot(tvec2,PPE_KT_R4_leo1,'xkcd:green',linewidth=2,linestyle='-',label='$R = 4$') - plt.plot(tvec2,PPE_KT_R8_leo1,'xkcd:orange',linewidth=2,linestyle='-',label='$R = 8$') - plt.plot(tvec2,PPE_KT_R20_leo1,'xkcd:magenta',linewidth=2,linestyle='-', label='$R = 20$') - plt.title('$kt-BLAST$',fontsize=20) - - if mode=='CS': - plt.plot(tvec,CT,'-k',linewidth=2,label='$ref$') - plt.plot(tvec2,ref_ao,'xkcd:red',linewidth=2,linestyle='--' , label='$aorta$') - plt.plot(tvec2,ref,'xkcd:red',linewidth=2,linestyle='-' , label='$leo$') - plt.plot(tvec2,CS2,'xkcd:blue',linewidth=2,linestyle='-' ,label='$R = 2$') - - #plt.plot(tvec2,PPE_R1_leo1,'xkcd:red',linewidth=2,linestyle='-' , label='$R = 1$') - #plt.plot(tvec2,PPE_CS_R2_leo1,'xkcd:blue',linewidth=2,linestyle='-' ,label='$R = 2$') - #plt.plot(tvec2,PPE_CS_R4_leo1,'xkcd:green',linewidth=2,linestyle='-',label='$R = 4$') - #plt.plot(tvec2,PPE_CS_R8_leo1,'xkcd:orange',linewidth=2,linestyle='-',label='$R = 8$') - #plt.plot(tvec2,PPE_CS_R20_leo1,'xkcd:magenta',linewidth=2,linestyle='-', label='$R = 20$') - plt.title('$Compressed \ \ Sensing$',fontsize=20) - - - - if mode=='STE': - plt.plot(tvec,CT,'-k',linewidth=2,label='$ref$') - plt.plot(tvec2,STE_CT , 'xkcd:purple' , label='STE-CT aorta') - plt.plot(tvec2,test_STE_CT , 'xkcd:purple' , linestyle='--', marker='o', label='STE-CT leo') - - - if mode=='STEint': - plt.plot(tvec,CT,'-k',linewidth=2,label='$ref$') - plt.plot(tvec2,STEint_CT, 'xkcd:aquamarine', label='STEint-CT aorta') - plt.plot(tvec2,test_STEint_CT, 'xkcd:aquamarine', linestyle='--', marker='o',label='STEint-CT aorta') - - - plt.ylim([-2,7]) - plt.xlabel(r'$time \ \ \ (s)$',fontsize=20) - plt.ylabel(r'$\delta p \ \ \ (mmHg) $',fontsize=20) - plt.legend(fontsize=16) - ############################################################################################################################## - - - #fig = plt.figure() - #plt.plot(tvec2,CT_fine,'-ok',label='CT') - ##plt.plot(tvec_T,CT_fine_T,'--k',label='CT T') - #plt.plot(tvec2 ,PPE_CT_fine ,'-om',label='PPE-CT aorta') - #plt.plot(tvec2 ,PPE_CT_fine_leo3,'-oc',label='PPE-CT leo3') - #plt.title('aorta fine') - #plt.xlabel(r'$time \ \ (s) $',fontsize=20) - #plt.ylabel(r'$\delta p \ \ \ (mmHg) $',fontsize=20) - #plt.legend(fontsize=14) - - - - - - plt.show() - def Plot_flux(masterpath,meshpath,options,mode,R): mesh = Mesh() @@ -1090,8 +970,7 @@ def Plot_dP(masterpath,options,mode,R): tcat[l-2] = float(row[0]) tcat = tcat+shift_t ax.plot(tcat[cstar[0]:tcat.size-cstar[1]],catheter[cstar[0]:tcat.size-cstar[1]],'white',linewidth=linesize,linestyle='--',label='$catheter$') - - + tm = np.linspace(0,Dt*tend,tend) ax.set_xlim([-0.05,0.81]) @@ -1131,226 +1010,7 @@ def Plot_dP(masterpath,options,mode,R): plt.show() - -def Plot_peaksystole(datapath,options,meshes,dt,R): - import pickle - barye2mmHg = 1/1333.22387415 - - - for mesh_size in meshes: - t_star = 6 - - - PPE_MEAN = np.zeros([len(R)]) - PPE_STD = np.zeros([len(R)]) - STE_MEAN = np.zeros([len(R)]) - STE_STD = np.zeros([len(R)]) - V_MEAN = np.zeros([len(R)]) - V_STD = np.zeros([len(R)]) - - ref_P = 0 - ref_V = 0 - - - - - if mesh_size=='Ucoarse': - ref_V = 326.95828118309191 - if mesh_size=='Ufine': - ref_V = 232.95021682714497 - if mesh_size=='Uffine': - ref_V = 234.66445211879045 - - - - - for l in range(len(R)): - if R[l]==0: - ref = np.loadtxt('/home/yeye/N_MRI/codes/pressure_drop/'+mesh_size+'/dt' + str(dt) + '/ref_'+mesh_size+'.txt') - PPE0_raw = open('/home/yeye/N_MRI/codes/pressure_drop/'+mesh_size+'/dt' + str(dt) + '/R0/pdrop_PPE_impl_stan.dat','rb') - STE0_raw = open('/home/yeye/N_MRI/codes/pressure_drop/'+mesh_size+'/dt' + str(dt) + '/R0/pdrop_STE_impl_stan.dat','rb') - PPE0 = pickle.load(PPE0_raw)['pdrop']*(-barye2mmHg) - STE0 = pickle.load(STE0_raw)['pdrop']*(-barye2mmHg) - - curpath = datapath + 'sequences/aorta_'+mesh_size+'.npz' - p = np.load(curpath) - px = p['x'] - py = p['y'] - pz = p['z'] - v = np.sqrt(px[:,:,:,t_star]**2 + py[:,:,:,t_star]**2 + pz[:,:,:,t_star]**2) - max0 = np.where(v==np.max(v)) - - ref_P = ref[6] - V_MEAN[l] = ref_V - V_STD[l] = 0 - PPE_MEAN[l] = PPE0[6] - PPE_STD[l] = 0 - STE_MEAN[l] = STE0[6] - STE_STD[l] = 0 - else: - PPE_MEAN[l] = np.loadtxt(datapath + 'maxpv/'+mesh_size + '/ppemean_R'+ str(R[l]) +'.txt') - PPE_STD[l] = np.loadtxt(datapath + 'maxpv/'+mesh_size + '/ppestd_R'+ str(R[l]) +'.txt') - STE_MEAN[l] = np.loadtxt(datapath + 'maxpv/'+mesh_size + '/stemean_R'+ str(R[l]) +'.txt') - STE_STD[l] = np.loadtxt(datapath + 'maxpv/'+mesh_size + '/stestd_R'+ str(R[l]) +'.txt') - V_MEAN[l] = np.loadtxt(datapath + 'maxpv/'+mesh_size + '/vmean_R'+ str(R[l]) +'.txt') - V_STD[l] = np.loadtxt(datapath + 'maxpv/'+mesh_size + '/vstd_R'+ str(R[l]) +'.txt') - - - - plt.figure(figsize=(10, 6), dpi=100) - Rvec = np.linspace(-10,R[-1]+5,100) - hline = Rvec*0+ref_P - plt.subplot(1,2,1) - plt.plot([0],[ref_P],color='k',marker='o',label= '$ref$') - plt.plot(Rvec,hline,color='k',linestyle='--') - plt.errorbar(R,PPE_MEAN, yerr=PPE_STD,color=options['ppecol'],marker='o',label= '$PPE$') - plt.errorbar(R,STE_MEAN, yerr=STE_STD,color=options['stecol'],marker='o',label= '$STE$') - plt.xlim([-2,R[-1]+5]) - plt.ylim([-7,18]) - plt.xlabel(r'$R$',fontsize=20) - plt.title('$Peak \ Systole \ pressure$',fontsize=18) - plt.ylabel(r'$ max \ \big ( \delta P \big ) \ \ mmHg$',fontsize=16) - - - plt.subplot(1,2,2) - #plt.plot([0],[ref_V],color='k',marker='o',label= '$ref$') - Rvec = np.linspace(-10,R[-1]+5,100) - hline = Rvec*0+ref_V - plt.plot(Rvec,hline,color='k',linestyle='--') - plt.errorbar(R,V_MEAN, yerr=V_STD,color='royalblue',marker='o',label= '$' + mesh_size + '$') - plt.xlim([-2,R[-1]+5]) - plt.ylim([0,350]) - plt.xlabel(r'$R$',fontsize=20) - plt.ylabel(r'$ max \ \big ( v \big ) \ \ cm/s$',fontsize=16) - plt.title('$Peak \ Systole \ velocity$',fontsize=18) - #plt.legend(fontsize=15) - - plt.annotate('$'+mesh_size+'$', xy=(-0.2, 1.1), xycoords='axes fraction',fontsize=15) - - - plt.show() - -def Plot_peaksystole_flux(datapath,options,meshes,dt,R): - import pickle - barye2mmHg = 1/1333.22387415 - - - for mesh_size in meshes: - t_star = 6 - PPE_MEAN = np.zeros([len(R)]) - PPE_STD = np.zeros([len(R)]) - STE_MEAN = np.zeros([len(R)]) - STE_STD = np.zeros([len(R)]) - V_MEAN = np.zeros([len(R)]) - V_STD = np.zeros([len(R)]) - Q_MEAN = np.zeros([len(R)]) - Q_STD = np.zeros([len(R)]) - - ref_V = 0 - - if mesh_size=='Ucoarse': - ref_V = 326.95828118309191 - if mesh_size=='Ufine': - ref_V = 232.95021682714497 - if 'Uffine' in mesh_size: - ref_V = 226.93523675462458 - maxv1 = 184.05091675316763 - ref_Q = 454.77517472437495 - maxq1 = 429.01393994253556 - - - for l in range(len(R)): - if R[l]==0: - ref = np.loadtxt('/home/yeye/N_MRI/codes/pressure_drop/'+mesh_size+'/dt' + str(dt) + '/ref_'+mesh_size+'.txt') - PPE0_raw = open('/home/yeye/N_MRI/codes/pressure_drop/'+mesh_size+'/dt' + str(dt) + '/R0/pdrop_PPE_impl_stan.dat','rb') - STE0_raw = open('/home/yeye/N_MRI/codes/pressure_drop/'+mesh_size+'/dt' + str(dt) + '/R0/pdrop_STE_impl_stan.dat','rb') - PPE0 = pickle.load(PPE0_raw)['pdrop']*(-barye2mmHg) - STE0 = pickle.load(STE0_raw)['pdrop']*(-barye2mmHg) - - curpath = datapath + 'sequences/aorta_'+mesh_size+'.npz' - p = np.load(curpath) - px = p['x'] - py = p['y'] - pz = p['z'] - v = np.sqrt(px[:,:,:,t_star]**2 + py[:,:,:,t_star]**2 + pz[:,:,:,t_star]**2) - max0 = np.where(v==np.max(v)) - V_MEAN[l] = ref_V - V_STD[l] = 0 - PPE_MEAN[l] = PPE0[6] - PPE_STD[l] = 0 - STE_MEAN[l] = STE0[6] - STE_STD[l] = 0 - elif R[l]==1: - PPE_MEAN[l] = 0 - PPE_STD[l] = 0 - STE_MEAN[l] = 0 - STE_STD[l] = 0 - V_MEAN[l] = maxv1 - V_STD[l] = 0 - Q_MEAN[l] = maxq1 - Q_STD[l] = 0 - - else: - PPE_MEAN[l] = np.loadtxt(datapath + 'maxpv/'+mesh_size + '/ppemean_R'+ str(R[l]) +'.txt') - PPE_STD[l] = np.loadtxt(datapath + 'maxpv/'+mesh_size + '/ppestd_R'+ str(R[l]) +'.txt') - STE_MEAN[l] = np.loadtxt(datapath + 'maxpv/'+mesh_size + '/stemean_R'+ str(R[l]) +'.txt') - STE_STD[l] = np.loadtxt(datapath + 'maxpv/'+mesh_size + '/stestd_R'+ str(R[l]) +'.txt') - V_MEAN[l] = np.loadtxt(datapath + 'maxpv/'+mesh_size + '/vmean_R'+ str(R[l]) +'.txt') - V_STD[l] = np.loadtxt(datapath + 'maxpv/'+mesh_size + '/vstd_R'+ str(R[l]) +'.txt') - Q_MEAN[l] = -np.loadtxt(datapath + 'maxpv/'+mesh_size + '/Q2_R'+ str(R[l]) +'.txt') - Q_STD[l] = -np.loadtxt(datapath + 'maxpv/'+mesh_size + '/Qstd2_R'+ str(R[l]) +'.txt') - - - - plt.figure(figsize=(15, 5), dpi=100) - plt.annotate('$'+mesh_size+'$', xy=(-0.2, 1.1), xycoords='axes fraction',fontsize=15) - - - Rvec = np.linspace(-10,R[-1]+5,100) - plt.subplot(1,3,1) - #plt.plot(Rvec,hline,color='k',linestyle='--') - plt.errorbar(R,PPE_MEAN, yerr=PPE_STD,color=options['ppecol'],marker='o',label= '$PPE$') - plt.errorbar(R,STE_MEAN, yerr=STE_STD,color=options['stecol'],marker='o',label= '$STE$') - plt.xlim([-2,R[-1]+5]) - plt.ylim([-0.1,2.0]) - plt.xlabel(r'$R$',fontsize=20) - plt.legend(fontsize=12,loc=2) - plt.title('$ l_2 \ error \ in \ Peak \ pressure$',fontsize=16) - plt.xticks(R) - plt.ylabel(r'$ || p - p^{ref} ||/|| p^{ref} ||_{l2} $',fontsize=16) - - - plt.subplot(1,3,2) - #plt.plot([0],[ref_V],color='k',marker='o',label= '$ref$') - Rvec = np.linspace(-10,R[-1]+7,100) - hline = Rvec*0+ref_V - plt.plot(Rvec,hline,color='royalblue',linestyle='--') - plt.errorbar(R,V_MEAN, yerr=V_STD,color='royalblue',marker='o',label= '$' + mesh_size + '$') - plt.xlim([-2,R[-1]+5]) - plt.ylim([0,350]) - plt.xlabel(r'$R$',fontsize=20) - plt.ylabel(r'$ v \ \ cm/s$',fontsize=16) - plt.xticks(R) - plt.title('$Peak \ velocity$',fontsize=16) - #plt.legend(fontsize=15) - - plt.subplot(1,3,3) - Rvec = np.linspace(-10,R[-1]+7,100) - hline = Rvec*0+ref_Q - plt.plot(Rvec,hline,color='mediumvioletred',linestyle='--') - plt.errorbar(R,Q_MEAN, yerr=Q_STD,color='mediumvioletred',marker='o',label= '$' + mesh_size + '$') - plt.xlim([-2,R[-1]+5]) - plt.xlabel(r'$R$',fontsize=20) - plt.ylabel(r'$ Q \ \ ml/s$',fontsize=16) - plt.xticks(R) - plt.ylim([150,550]) - plt.title('$Peak \ Flux $',fontsize=16) - - - - - plt.show() - + def CLOCK(t1, t2): tot_time = np.round(t2 - t1, 2) if tot_time < 60: @@ -1411,6 +1071,9 @@ def ROUTINE(options): if options['Error-curves']['apply']: print('--- Error-curves analysis ---') ratio = False + fac=1 + if '11mm' in options['Error-curves']['subfolders']: + fac=100 for types in options['Error-curves']['type']: if types=='mean_ratio': types = 'mean' @@ -1418,6 +1081,9 @@ def ROUTINE(options): if types=='max_ratio': types = 'max' ratio = True + if types=='norm2_m': + types='norm2' + meas=True nc = 0 if len(options['Error-curves']['subfolders'])==0: ucomp = [] @@ -1445,46 +1111,262 @@ def ROUTINE(options): for subf in options['Error-curves']['subfolders']: ucomp = [] wcomp = [] - colorset = options['Error-curves']['colors'] + if 'colors' in options['Error-curves']: + print('colors setted...') + colorset = options['Error-curves']['colors'] + colorsetted = True + else: + colorsetted = False + styles = options['Error-curves']['styles'] labelset = options['Error-curves']['labels'] path = options['Error-curves']['folder'] + subf + '/' - try: - ucomp = np.loadtxt(path + '/u'+types+'.txt') - wcomp = np.loadtxt(path + '/w'+types+'.txt') - times = np.loadtxt(path + '/times.txt') - except IOError: - print('no cError-curves for ' + subf) + wcomp = np.loadtxt(path + 'w'+types+'.txt') + times = np.loadtxt(path + 'times.txt') + if types != 'grad': + ucomp = np.loadtxt(path + 'u'+types+'.txt') - if not ratio: - plt.plot( - times, ucomp, color=colorset[nc], linestyle='-', label= '$u'+ subf +'$' ) - plt.plot( - times, wcomp, color=colorset[nc], linestyle='--', label='$w'+subf+'$') + if colorsetted: + if not ratio: + if types not in ['grad']: + if meas: + plt.plot(times, fac*ucomp, color=colorset[nc], linestyle='--', label= '$ u \ '+ labelset[nc] +'$',linewidth=2.5 ) + plt.plot(times, fac*wcomp, color=colorset[nc], linestyle=styles[nc], label= '$ w \ '+ labelset[nc] +'$',linewidth=2.5) + else: + wu = wcomp/ucomp + plt.plot( + times, wu, color=colorset[nc], linestyle=styles[nc], label= '$'+ labelset[nc] +'$' ) + nc +=1 else: - wu = wcomp/ucomp - plt.plot( - times, wu, color=colorset[nc], linestyle='-', label= '$'+ labelset[nc] +'$' ) - nc +=1 + if not ratio: + if types not in ['grad']: + if meas: + plt.plot(times, fac*ucomp, linestyle='--', label= '$'+ labelset[nc] +'$' ) + else: + plt.plot(times, fac*ucomp, linestyle='--', label= '$'+ labelset[nc] +'$' ) + plt.plot( + times, fac*wcomp, linestyle=styles[nc], label= '$'+ labelset[nc] +'$') + else: + wu = wcomp/ucomp + plt.plot( + times, wu, linestyle=styles[nc], label= '$'+ labelset[nc] +'$' ) + nc +=1 + - #plt.ylim([0, 170]) plt.xlabel('$time \ \ (s)$', fontsize=18) - plt.legend(fontsize=16) + plt.legend(fontsize=14) if options['Error-curves']['title']: plt.title(options['Error-curves']['title'], fontsize=18) if not ratio: - plt.ylabel('$velocity \ \ (cm/s)$', fontsize=18) + if types == 'grad': + plt.ylim([0, 230]) + plt.ylabel('$ |grad(w)|_2 \ \ (1/s)$', fontsize=18) + elif types == 'norm2': + plt.ylim([0, 52]) + plt.ylabel('$ |w|_2 \ \ (cm/s)$', fontsize=18) + else: + plt.ylabel('$velocity \ \ (cm/s)$', fontsize=18) plt.savefig(options['Error-curves']['outpath'] + types + '.png', dpi=500, bbox_inches='tight') else: - plt.ylabel('$w/u$', fontsize=18) - if 'max' in types: - plt.ylim([0, 1.8]) + plt.ylabel('$|w/u|$', fontsize=18) + #if 'max' in types: + #plt.ylim([0, 1.8]) plt.savefig(options['Error-curves']['outpath'] + types + '_ratio.png', dpi=500, bbox_inches='tight') plt.show() + if 'Histograms_checkpoint' in options: + if options['Histograms_checkpoint']['apply']: + print('--- Histograms ---') + path = options['Histograms_checkpoint']['path'] + freq = np.loadtxt(path + 'hist_freq.txt') + edges = np.loadtxt(path + 'hist_edges.txt') + fig, ax = plt.subplots() + #ax.bar(edges[:-1], freq, width=np.diff(edges), edgecolor="black", align="edge") + ax.bar(edges[:-1], freq, width=np.diff(edges), align="edge") + plt.title(options['Histograms_checkpoint']['title'], fontsize=18) + plt.xlim([0 , 50]) + plt.ylim([0 , 1.8]) + plt.savefig(path + 'hist.png', dpi=500, bbox_inches='tight') + plt.show() + if 'Pressure_drops' in options: + if options['Pressure_drops']['apply']: + print('--- Pressure_drops ---') + import pickle + nc = 0 + tommhg = options['Pressure_drops']['convertor'] + + for subf in options['Pressure_drops']['subfolders']: + ucomp = [] + wcomp = [] + if 'colors' in options['Pressure_drops']: + colorset = options['Pressure_drops']['colors'] + colorsetted = True + else: + colorsetted = False + styles = options['Pressure_drops']['styles'] + labelset = options['Pressure_drops']['labels'] + path = options['Pressure_drops']['folder'] + subf + '/' + dataname = 'pdrop_COR_impl_stan.dat' + if 'STE' in path: + dataname = 'pdrop_STE_impl_stan.dat' + if labelset[nc]=='ref': + dataname = 'pdrop.dat' + data = open(path+dataname, 'rb') + p_drop = pickle.load(data)['pdrop']/tommhg + data = open(path+dataname, 'rb') + times = pickle.load(data)['time'] + + if labelset[nc] == 'ref': + plt.plot(times, p_drop,color='black', linestyle='-', label= '$ref$' ) + else: + + if colorsetted: + plt.plot( + times, p_drop, color=colorset[nc], linestyle=styles[nc], label= '$'+ labelset[nc] +'$' ) + else: + plt.plot(times, p_drop, linestyle=styles[nc], label= '$'+ labelset[nc] +'$' ) + + if options['Pressure_drops']['catheter']: + + c_path = '/home/yeye/Desktop/PhD/MEDICAL_DATA/Study_David/catheter_data/catheter_'+ labelset[nc]+'_rest_stats.csv' + + + with open(c_path, 'r') as csvfile: + mylist = [row[0] for row in csv.reader(csvfile, delimiter=';')] + + Values = np.array(mylist) + catheter = np.zeros([len(Values)-2]) + tcat = np.zeros([len(Values)-2]) + for l in range(2,len(Values)): + row = Values[l].split(',') + catheter[l-2] = float(row[5]) + tcat[l-2] = float(row[0]) + + + if '11mm' in subf: + tdelay = 0.015 + elif '9mm' in subf: + tdelay = -0.12 + elif '13mm' in subf: + tdelay = 0.1 + elif 'Normal' in subf: + tdelay = -0.01 + + plt.plot(tcat+tdelay,catheter,color=colorset[nc],linestyle='--')# ,label='$cat' + subf + '$') + + nc +=1 + + + #plt.ylim([0, 170]) + plt.xlabel('$time \ \ (s)$', fontsize=18) + plt.legend(fontsize=14) + if options['Pressure_drops']['title']: + plt.title(options['Pressure_drops']['title'], fontsize=18) + plt.ylabel('$ \delta P \ \ (mmHg)$', fontsize=18) + plt.savefig(options['Pressure_drops']['outpath'] + 'pressure_drops.png', dpi=500, bbox_inches='tight') + plt.show() + + if 'l2_comp' in options: + if options['l2_comp']['apply']: + print('--- L2 component analysis ---') + fig, ax = plt.subplots() + for subf in options['l2_comp']['subfolder']: + path = options['l2_comp']['folder'] + subf + '/' + colors = options['l2_comp']['colors'] + mode = options['l2_comp']['mode']['type'] + + if mode in ['gain','gain_compressed']: + gain = True + path_to_comp = options['l2_comp']['mode']['comp'] + wx = np.loadtxt(path_to_comp + '/wx.txt') + wy = np.loadtxt(path_to_comp + '/wy.txt') + wz = np.loadtxt(path_to_comp + '/wz.txt') + else: + gain = False + wx = np.loadtxt(path + '/wx.txt') + wy = np.loadtxt(path + '/wy.txt') + wz = np.loadtxt(path + '/wz.txt') + + times = np.loadtxt(path + '/times.txt') + if subf != 'SNRinfV120' and gain: + varux = np.loadtxt(path + '/varux.txt') + varuy = np.loadtxt(path + '/varuy.txt') + varuz = np.loadtxt(path + '/varuz.txt') + if 'SNRinfV120' in subf: + lsty = '--' + labels = ['','','','',''] + else: + lsty = '-' + lsty2 = '--' + labels = ['$wx$','$wy$','$wz$','$div$'] + labels2 = ['$\delta u_x$','$\delta u_y$','$\delta u_z$'] + labels3 = ['$x$','$y$','$z$'] + + if mode == 'gain': + plt.plot(times, varux, color = colors[0], linestyle=lsty2 , label= labels2[0] ) + plt.plot(times, varuy, color = colors[1], linestyle=lsty2 , label= labels2[1] ) + plt.plot(times, varuz, color = colors[2], linestyle=lsty2 , label= labels2[2] ) + plt.plot(times, wx, color = colors[0], linestyle=lsty, label= labels[0] ) + plt.plot(times, wy, color = colors[1], linestyle=lsty, label= labels[1] ) + plt.plot(times, wz, color = colors[2], linestyle=lsty, label= labels[2] ) + elif mode == 'gain_compressed': + plt.plot(times, varux-wx, color = colors[0], linestyle=lsty, label= labels3[0] ) + plt.plot(times, varuy-wy, color = colors[1], linestyle=lsty, label= labels3[1] ) + plt.plot(times, varuz-wz, color = colors[2], linestyle=lsty, label= labels3[2] ) + else: + plt.plot(times, wx, color = colors[0], linestyle=lsty, label= labels[0] ) + plt.plot(times, wy, color = colors[1], linestyle=lsty, label= labels[1] ) + plt.plot(times, wz, color = colors[2], linestyle=lsty, label= labels[2] ) + + plt.xlabel('$time \ \ (s)$', fontsize=18) + + if options['l2_comp']['div']: + div = np.loadtxt(path + 'div.txt') + div_rescaled = div*0.01 + div_rescaled = div_rescaled + 0.1 + plt.plot(times, div_rescaled, color = 'indigo', linestyle=lsty, label= labels[3] ) + + if 'title' in options['l2_comp']: + if options['l2_comp']['title']: + plt.title(options['l2_comp']['title'], fontsize=18) + + + + if options['l2_comp']['aliasing']: + if 'Poiseuille' in options['l2_comp']['folder']: + print('adding alaising color in Poiseuille') + import matplotlib.transforms as mtransforms + trans = mtransforms.blended_transform_factory(ax.transData, ax.transAxes) + if 'SNR10' in subf: + time_al = [0.15,0.78] + elif 'SNRinf' in subf: + time_al = [0.21,0.78] + pm_al = 0.5*(time_al[1] + time_al[0]) + r_al = 0.5*(time_al[1] - time_al[0]) + ax.fill_between(times, 0, 1, where=np.abs(times -pm_al)<= r_al ,facecolor='gold', alpha=0.4, transform=trans) + if mode == 'gain_compressed': + ax.text(pm_al/times[-1], 0.55, '$aliasing$', horizontalalignment='center',verticalalignment='center', transform=ax.transAxes,fontsize=17) + else: + ax.text(pm_al/times[-1], 0.82, '$aliasing$', horizontalalignment='center',verticalalignment='center', transform=ax.transAxes,fontsize=17) + + leg = plt.legend(fancybox=True,fontsize=16) + leg.get_frame().set_linewidth(0.0) + ax.tick_params(axis='both', which='major', labelsize=14) + #plt.ylim([-0.005 , 0.235]) + + if mode == 'gain_compressed': + plt.ylim([-0.25 , 1.75]) + plt.ylabel('$|| \delta u ||_{L2} - || w ||_{L2}$', fontsize=18) + if not gain: + plt.ylim([-0.005 , 0.5]) + plt.ylabel('$ \sqrt{\int w^2 dx /| \Omega|}/ venc$', fontsize=18) + + plt.savefig(options['l2_comp']['folder'] + options['l2_comp']['name'] + '.png', dpi=500, bbox_inches='tight') + plt.show() + diff --git a/codes/MRI.py b/codes/MRI.py index a1eab4b..9707027 100644 --- a/codes/MRI.py +++ b/codes/MRI.py @@ -3,13 +3,10 @@ # Workspace for MRI analysis of the 4Dflow data # # written by Jeremias Garay L: j.e.garay.labra@rug.nl -# Fernanda te amo # for autoreload in ipython3 # %load_ext autoreload # %autoreload 2 ################################################################ -import h5py -from dPdirectestim.dPdirectestim import * from dolfin import * import dolfin import numpy as np @@ -999,14 +996,8 @@ def CLOCK(rank,t1,t2): print('Total time: ' + str(time_hour) + ' hrs, ' + str(time_min) + ' min, ' + str(time_sec) + ' s') - def SCANNER(options): - ######################################## - # - # Basic Tools - # - ######################################## if 'kspace_cib' in options: if options['kspace_cib']['apply']: print('--- kspace from CIB data ---') @@ -1416,24 +1407,18 @@ def SCANNER(options): else: CIBtoH5(path_to_cib,times,dt,outpath,flip=flip) - ######################################## - # - # Undersampling - # - ######################################## 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 options['cs']['short']: - [Mx,My,Mz] = LOADsequences(options['cs']['Mpath']) - CS.undersampling_short(Mx,My,Mz,options) + 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']) @@ -1457,12 +1442,6 @@ def SCANNER(options): print('saving the sequences' + options['SENSE']['savepath']) np.savez_compressed(options['SENSE']['savepath'], x=MxS,y=MyS,z=MzS) - ######################################## - # - # Writing Checkpoint from Sequence - # - ######################################## - if 'create_checkpoint' in options: if options['create_checkpoint']['apply']: print('--- Create Checkpoint ---') @@ -1529,7 +1508,7 @@ def SCANNER(options): comm = MESH['mesh'].mpi_comm() dt = options['create_checkpoint']['dt'] if options['create_checkpoint']['xdmf']: - xdmf_u = XDMFFile(options['create_checkpoint']['savepath']+'u.xdmf') + xdmf_u = XDMFFile(options['create_checkpoint']['savepath']+ 'R' + str(r) + '/u.xdmf') for l in range(len(vel_seq)): if rank == 0: @@ -1545,11 +1524,6 @@ def SCANNER(options): inout.write_HDF5_data( comm, path + '/u.h5', vel_seq[l], '/u', t=l*dt) - ######################################## - # - # Relative Pressure Estimators - # - ######################################## if 'reference' in options: if options['reference']['apply']: if rank == 0: @@ -1688,157 +1662,6 @@ def SCANNER(options): if rank==0: print(' ') - if 'peak_pv' in options: - - if options['peak_pv']['apply']: - import CS - import pickle - import sys - import logging - # DPestim - logging.getLogger().setLevel(logging.INFO) - parameters['form_compiler']['optimize'] = True - parameters['form_compiler']['cpp_optimize'] = True - parameters['form_compiler']['cpp_optimize_flags'] = '-O3 -ffast-math -march=native' - infile_dp = options['peak_pv']['infile_dp'] - estimator = DPDirectEstim(infile_dp) - - barye2mmHg = 1/1333.22387415 - t_star = 0.185 - if rank==0: - print('Computing Velocity and Pressure at Peak Systole') - print('The inlet max occurs at ' + str(t_star) + ' sec') - - - [BOX,AORTA,LEO] = CREATEmeshes(options) - - if options['peak_pv']['mesh_into']=='leo': - MESH = LEO - del AORTA - if options['peak_pv']['mesh_into']=='aorta': - MESH = AORTA - del LEO - - if options['peak_pv']['p_comp']=='error': - if rank==0: - print('Reading Pressure Reference') - P0_PPE = READcheckpoint(MESH,'p',options,options['checkpoint_path'],'p_PPE_impl_stan') - P0_STE = READcheckpoint(MESH,'p',options,options['checkpoint_path'],'p_STE_impl_stan') - - [Sqx,Sqy,Sqz] = LOADsequences(options['peak_pv']['orig_seq']) - Nite = options['peak_pv']['N_CS'] - [row,col,dep,numt] = Sqz.shape - frms_num = 3 - peakslice = options['peak_pv']['peak_slice'] - R = options['peak_pv']['R'] - - v0 = np.sqrt(Sqx[:,:,:,peakslice]**2 + Sqy[:,:,:,peakslice]**2 + Sqz[:,:,:,peakslice]**2) - max0 = np.where(v0==np.max(v0)) - - qqvec = {} - for ss in options['peak_pv']['flux_bnd']: - qqvec[ss] = np.zeros([Nite]) - - ppemax = np.zeros([Nite]) - stemax = np.zeros([Nite]) - vmax = np.zeros([Nite]) - slrdmax = peakslice - # Selecting around the max only - slrdmax = 1 - Sqx = Sqx[:,:,:,peakslice-1:peakslice+2] - Sqy = Sqy[:,:,:,peakslice-1:peakslice+2] - Sqz = Sqz[:,:,:,peakslice-1:peakslice+2] - - for l in range(len(R)): - if rank==0: - print('Peak Systole velocity and pressure at R = ' + str(R[l])) - - sx_cs = np.zeros([row,col,dep,frms_num,Nite]) - sy_cs = np.zeros([row,col,dep,frms_num,Nite]) - sz_cs = np.zeros([row,col,dep,frms_num,Nite]) - - for k in range(Nite): - if rank==0: - print('CS iteration number ' + str(k+1)) - [xcs,ycs,zcs] = CS.undersampling_peakpv(Sqx,Sqy,Sqz,options,R[l]) - sx_cs[:,:,:,:,k] = xcs - sy_cs[:,:,:,:,k] = ycs - sz_cs[:,:,:,:,k] = zcs - vk = np.sqrt(sx_cs[:,:,:,1,k]**2 + sy_cs[:,:,:,1,k]**2 + sz_cs[:,:,:,1,k]**2) - vmax[k] = vk[max0] - - if rank==0: - print('\n CS done') - - # To write the checkpoints - vel_seq = SqtoH5(BOX,MESH,sx_cs[:,:,:,:,k],sy_cs[:,:,:,:,k],sz_cs[:,:,:,:,k]) - comm = MESH['mesh'].mpi_comm() - - # Computing the Fluxes - if rank==0: - print('\n Computing the Flux') - QQ = Fluxes(MESH,vel_seq,options,options['peak_pv']['flux_bnd']) - - for ss in options['peak_pv']['flux_bnd']: - qqvec[ss][k] = QQ[ss][slrdmax] - - if rank==0: - print('\n Writing checkpoints') - - for ns in range(len(vel_seq)): - pathss = options['peak_pv']['savepath'] + 'H5/checkpoint/{i}/'.format(i=ns) - if l<10 and l>0: - pathss = options['peak_pv']['savepath'] + 'H5/checkpoint/0{i}/'.format(i=ns) - inout.write_HDF5_data(comm, pathss + '/u.h5', vel_seq[ns], '/u', t=0) - if rank==0: - print('\n The checkpoints were wrote') - - # Computing the Pressure Drop - estimator.estimate() - # Reading the results - if options['peak_pv']['p_comp']=='peak': - ppe_raw = open(options['peak_pv']['savepath'] + '/H5/pdrop_PPE_impl_stan.dat','rb') - ste_raw = open(options['peak_pv']['savepath'] + '/H5/pdrop_STE_impl_stan.dat','rb') - ppe = pickle.load(ppe_raw)['pdrop']*(-barye2mmHg) - ste = pickle.load(ste_raw)['pdrop']*(-barye2mmHg) - p1max[k] = ppe[slrdmax] - p2max[k] = ste[slrdmax] - elif options['peak_pv']['p_comp']=='error': - PPE = READcheckpoint(MESH,'p',options,options['peak_pv']['savepath']+'H5/checkpoint/','p_PPE_impl_stan') - STE = READcheckpoint(MESH,'p',options,options['peak_pv']['savepath']+'H5/checkpoint/','p_STE_impl_stan') - ppe_vec_0 = P0_PPE[peakslice].vector().get_local() - P0_PPE[peakslice].vector().get_local()[0] - ppe_vec = PPE[slrdmax].vector().get_local() - PPE[slrdmax].vector().get_local()[0] - ste_vec_0 = P0_STE[peakslice].vector().get_local() - P0_STE[peakslice].vector().get_local()[0] - ste_vec = STE[slrdmax].vector().get_local() - STE[slrdmax].vector().get_local()[0] - ppemax[k] = np.linalg.norm(ppe_vec_0 - ppe_vec)/np.linalg.norm(ppe_vec_0) - stemax[k] = np.linalg.norm(ste_vec_0 - ste_vec)/np.linalg.norm(ste_vec_0) - else: - raise Exception('Pressure computation not recognize!') - - - # VELOCITIES - vmean = np.mean(vmax) - vstd = np.std(vmax) - # PRESSURES - ppemean = np.mean(ppemax) - stemean = np.mean(stemax) - ppestd = np.std(ppemax) - stestd = np.std(stemax) - - - if options['peak_pv']['save']: - if rank==0: - print('\n saving the files in ' + options['peak_pv']['savepath']) - for ss in options['peak_pv']['flux_bnd']: - np.savetxt( options['peak_pv']['savepath'] + 'Q'+str(ss) +'_R'+str(R[l])+'.txt', [np.mean(qqvec[ss])]) - np.savetxt( options['peak_pv']['savepath'] + 'Qstd'+str(ss) +'_R'+str(R[l])+'.txt', [np.std(qqvec[ss])]) - np.savetxt( options['peak_pv']['savepath'] + 'ppemean_R'+str(R[l])+'.txt', [ppemean] ) - np.savetxt( options['peak_pv']['savepath'] + 'stemean_R'+str(R[l])+'.txt', [stemean] ) - np.savetxt( options['peak_pv']['savepath'] + 'vmean_R'+str(R[l])+'.txt', [vmean] ) - np.savetxt( options['peak_pv']['savepath'] + 'ppestd_R'+str(R[l])+'.txt', [ppestd] ) - np.savetxt( options['peak_pv']['savepath'] + 'stestd_R'+str(R[l])+'.txt', [stestd] ) - np.savetxt( options['peak_pv']['savepath'] + 'vstd_R'+str(R[l])+'.txt', [vstd] ) - if 'change_mesh' in options: if options['change_mesh']['apply']: if rank == 0: @@ -1861,9 +1684,9 @@ def SCANNER(options): changed = {} #W = LEO['FEM'].sub(0).collapse() comm = MESH_out['mesh'].mpi_comm() - v2 = Function(MESH_out['FEM']) for k in range(len(list(origin))): + v2 = Function(MESH_out['FEM']) if rank == 0: print('CHANGING: index', k) LagrangeInterpolator.interpolate(v2, origin[k]) @@ -1878,7 +1701,7 @@ def SCANNER(options): if rank == 0: print('saving checkpoint', l) path = options['change_mesh']['savepath'] + \ - 'R1/checkpoint/{i}/'.format(i=l) + '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) diff --git a/codes/PostCheck.py b/codes/PostCheck.py index f3d75e4..0187b37 100644 --- a/codes/PostCheck.py +++ b/codes/PostCheck.py @@ -109,14 +109,10 @@ def WORKcheck(MESH, mode, output_path, checkpoint_path, filename, outname, optio if mode == 'p' or mode == 'p_cib': xdmf_p = XDMFFile(output_path+outname+'.xdmf') - for k in range(0, len(indexes), 1): + 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' - - if filename == 'p': - if k < 10 and k > 0: - path = checkpoint_path + '0' + \ - str(indexes[k]) + '/'+filename+'.h5' - p = Function(W) if mode == 'p': barye2mmHg = 1/1333.22387415 @@ -126,36 +122,10 @@ def WORKcheck(MESH, mode, output_path, checkpoint_path, filename, outname, optio 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) + #p1vec = p1vec - np.mean(p1vec) p.vector()[:] = p1vec*barye2mmHg xdmf_p.write(p, te) - te = te + dt - numt = numt + 1 - - 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' - - if indexes[k] < 10 and indexes[k] > 0: - path = checkpoint_path + '0' + \ - str(indexes[k]) + '/'+filename+'.h5' - - v = Function(V) - dv = Function(W) - - dv.rename('div', outname) - comm = MPI.COMM_WORLD - 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, te) - te = te + dt - numt = numt + 1 def READcheckpoint(MESH, mode, output_path, checkpoint_path, filename, outname, options, flow=False, bnds=None): @@ -203,7 +173,7 @@ def READcheckpoint(MESH, mode, output_path, checkpoint_path, filename, outname, freq,edges = np.histogram(ValuesPeak, bins=80, density=True) #Saving the histogram print('Saving at ' + output_path) - np.savetxt(output_path + 'hist_frew.txt', freq) + np.savetxt(output_path + 'hist_freq.txt', freq) np.savetxt(output_path + 'hist_edges.txt', edges) if mode == 'perturbation': @@ -235,6 +205,20 @@ def READcheckpoint(MESH, mode, output_path, checkpoint_path, filename, outname, 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') @@ -244,9 +228,7 @@ def READcheckpoint(MESH, mode, output_path, checkpoint_path, filename, outname, uvec = u.vector().get_local() if Phase_Contrast: - ufactor = options['Perturbation']['type']['phase_contrast']/100 - VENC = np.max(np.abs(uvec))*ufactor - gamma = 267.513e6 # rad/Tesla/sec Gyromagnetic ratio for H nuclei + #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 @@ -474,7 +456,11 @@ def ERRORmap(MESH, mode, outpath, reference_path, checkpoint_path, refname,comna 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() @@ -486,6 +472,80 @@ def ERRORmap(MESH, mode, outpath, reference_path, checkpoint_path, refname,comna 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: @@ -495,27 +555,54 @@ def ERRORmap(MESH, mode, outpath, reference_path, checkpoint_path, refname,comna 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 = [] - dt = options['Error-curves']['dt'] - for k in range(1,len(indexes)): - path_w = checkpoint_path + str(indexes[k]) + '/'+comname+'.h5' - path_u = reference_path + str(indexes0[k]) + '/'+refname+'.h5' + + 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 comname: + if 'w' in comname2: hdf_w.read(w, 'w/vector_0') else: hdf_w.read(w, 'u/vector_0') - hdf_u = HDF5File(MESH['mesh'].mpi_comm(), path_u, 'r') - hdf_u.read(u, 'u/vector_0') - hdf_u.close() + 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() - u_vec = u.vector().get_local() w_vec = w.vector().get_local() print('computing error in timestep numer',k) @@ -525,59 +612,53 @@ def ERRORmap(MESH, mode, outpath, reference_path, checkpoint_path, refname,comna 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) - - np.savetxt(outpath+'u' +typs+'.txt',ucomp) - np.savetxt(outpath+'w' +typs+'.txt',wcomp) + + 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 == 'histogram': - from pathlib import Path - import pickle - u = Function(V) - w = Function(V) - errors = {} - - for k in range(len(indexes)): - path_w = checkpoint_path + str(indexes[k]) + '/'+comname+'.h5' - path_u = reference_path + str(indexes0[k]) + '/'+refname+'.h5' - - if indexes0[k] < 10: - path_u = reference_path + '0' + \ - str(indexes0[k]) + '/'+refname+'.h5' - - u.rename('meas', 'meas') - w.rename('w', 'w') - 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() - errors[k] = np.zeros(u_vec.size) - - for l in range(len(errors[k])): - #errors[k][l] = np.nan - if u_vec[l]<1e-9: - errors[k][l] = -1 - else: - eta = np.abs(w_vec[l]/u_vec[l]) - if np.abs(eta)>50: - errors[k][l] = -1 - else: - errors[k][l] = eta - - - write_path = Path(outpath) - fpath = write_path.joinpath('errors.dat') - pickle.dump(errors, fpath.open('wb')) - if mode == 'h5': xdmf_u = XDMFFile(outpath+'meas.xdmf') #xdmf_ucor = XDMFFile(output_path+'ucor.xdmf') @@ -592,10 +673,6 @@ def ERRORmap(MESH, mode, outpath, reference_path, checkpoint_path, refname,comna path_w = checkpoint_path + str(indexes[k]) + '/'+comname+'.h5' path_u = reference_path + str(indexes0[k]) + '/'+refname+'.h5' - if indexes0[k] < 10: - path_u = reference_path + '0' + \ - str(indexes0[k]) + '/'+refname+'.h5' - u.rename('meas', 'meas') w.rename('w', 'w') #ucor.rename('ucor', 'ucor') @@ -612,34 +689,173 @@ def ERRORmap(MESH, mode, outpath, reference_path, checkpoint_path, refname,comna #ucor.vector()[:] = u_vec + w_vec xdmf_u.write(u, k) - #xdmf_ucor.write(ucor, k) xdmf_w.write(w, k) - if mode == 'colormap': - colormap = XDMFFile(outpath+'colormap.xdmf') + 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(W) - dt = options['Corrector']['dt'] + cm = Function(DGs) + vs = TestFunction(DGs) for k in range(len(indexes)): - print('making the colormap in the time',np.round(k*dt,2)) 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('color','color') + 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() - uvec = u.vector().get_local() - wvec = w.vector().get_local() - cm.vector()[:] = np.sqrt((uvec - wvec)**2) - colormap.write(cm, k*dt) + 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') @@ -874,6 +1090,11 @@ 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] @@ -897,22 +1118,25 @@ def ESTIMpressure(MESH, outpath, checkpoint_path, filename, options): 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) - 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) - p = Function(W) - - one_mesh = interpolate(Constant(1), W) - LagrangeInterpolator.interpolate(s1, one_mesh) - LagrangeInterpolator.interpolate(s2, one_mesh) - vol1 = assemble(s1*dx) - vol2 = assemble(s1*dx) - - dt = options['Estim_Pressure']['dt'] - p_drop_lst = [] - time_ = [] + 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' @@ -920,11 +1144,14 @@ def ESTIMpressure(MESH, outpath, checkpoint_path, filename, options): hdf.read(p, 'p/vector_0') hdf.close() - LagrangeInterpolator.interpolate(s1, p) - LagrangeInterpolator.interpolate(s2, p) - - P1 = assemble(s1*dx)/vol1 - P2 = assemble(s2*dx)/vol2 + 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) @@ -1012,41 +1239,60 @@ def OUTLETwind(MESH, output_path, checkpoint_path, filename, bnds): 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 ---') - mesh_path = options['Outlet_Wind']['mesh_path'] + 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(mesh_path) + MESH = LOADmesh(meshpath) OUTLETwind(MESH, out_path, checkpoint, filename, bnds) - if 'Corrector' in options: - if options['Corrector']['apply']: - if rank == 0: - print('Applying Corrector') + if 'Colormap' in options: + if options['Colormap']['apply']: - MESH = LOADmesh(options['Corrector']['mesh_path']) - u_path = options['Corrector']['u_path'] + 'checkpoint/' - w_path = options['Corrector']['w_path'] + 'checkpoint/' - wname = options['Corrector']['wname'] - uname = options['Corrector']['uname'] - mode = options['Corrector']['mode'] - outpath = options['Corrector']['outpath'] + for mode in options['Colormap']['mode']: + if rank == 0: + print('Applying Colormap') + print('Choosen mode: ' + mode) - ERRORmap(MESH, mode, outpath, u_path, - w_path, uname, wname, options) + 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']['mesh_path']) + MESH = LOADmesh(options['Velocity']['meshpath']) filename = options['Velocity']['filename'] checkpoint_path = options['Velocity']['checkpoint'] + 'checkpoint/' outpath = options['Velocity']['checkpoint'] @@ -1067,7 +1313,7 @@ def ROUTINE(options): if options['W_map']['apply']: if rank == 0: print('Applying W--MAP') - MESH = LOADmesh(options['mesh_path']) + MESH = LOADmesh(options['meshpath']) filename = options['W_map']['filename'] outname = options['W_map']['out_name'] mode = 'w' @@ -1078,7 +1324,7 @@ def ROUTINE(options): if options['GradW_map']['apply']: if rank == 0: print('Applying Grad W--MAP') - MESH = LOADmesh(options['mesh_path']) + MESH = LOADmesh(options['meshpath']) filename = options['GradW_map']['filename'] outname = options['GradW_map']['out_name'] mode = 'gradw' @@ -1088,21 +1334,21 @@ def ROUTINE(options): if 'Pressure' in options: if options['Pressure']['apply']: if rank == 0: - print('Applying Pressure-MAP') + print('---Applying Pressure---') - MESH = LOADmesh(options['mesh_path']) + MESH = LOADmesh(options['meshpath']) + outpath = options['Pressure']['checkpoint'] + checkpoint_path = options['Pressure']['checkpoint'] + 'checkpoint/' filename = options['Pressure']['filename'] - outname = options['Pressure']['out_name'] - mode = 'p' - WORKcheck(MESH, mode, output_path, checkpoint_path, - filename, outname, options) + 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['mesh_path']) + MESH = LOADmesh(options['meshpath']) refname = options['Error_P']['refname'] reference_path = options['Error_P']['refpath'] filename = options['Error_P']['filename'] @@ -1115,10 +1361,13 @@ def ROUTINE(options): if options['Estim_Pressure']['apply']: if rank == 0: print('Applying Pressure Estimator') + print('Method choosed: '+options['Estim_Pressure']['method']) - MESH = LOADmesh(options['mesh_path']) + MESH = LOADmesh(options['Estim_Pressure']['meshpath']) filename = options['Estim_Pressure']['filename'] - outpath = options['Estim_Pressure']['outpath'] + outpath = options['Estim_Pressure']['checkpath'] + checkpoint_path = options['Estim_Pressure']['checkpath'] + 'checkpoint/' + ESTIMpressure(MESH, outpath, checkpoint_path, filename, options) if 'Error-curves' in options: @@ -1162,7 +1411,11 @@ def ROUTINE(options): MESH = LOADmesh(options['Perturbation']['meshpath']) checkpath = options['Perturbation']['checkpath'] + 'checkpoint/' - name_folder = 'SNR' +options['Perturbation']['type']['SNR'] + 'V' + str(options['Perturbation']['type']['phase_contrast']) + 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) diff --git a/codes/__pycache__/CS.cpython-36.pyc b/codes/__pycache__/CS.cpython-36.pyc index c0eede0..a3acd0f 100644 Binary files a/codes/__pycache__/CS.cpython-36.pyc and b/codes/__pycache__/CS.cpython-36.pyc differ diff --git a/codes/__pycache__/Graphics.cpython-36.pyc b/codes/__pycache__/Graphics.cpython-36.pyc index d0cf857..1cad2a7 100644 Binary files a/codes/__pycache__/Graphics.cpython-36.pyc and b/codes/__pycache__/Graphics.cpython-36.pyc differ diff --git a/codes/monolithic.py b/codes/monolithic.py index 193d96e..8576d8a 100644 --- a/codes/monolithic.py +++ b/codes/monolithic.py @@ -1,5 +1,4 @@ from dolfin import * -import matplotlib.pyplot as plt import numpy as np from common import inout from mpi4py import MPI @@ -23,7 +22,6 @@ def solv_NavierStokes(options): theta = Constant(options['theta']) Tf = options['Tf'] dt = options['dt'] - dt_w = options['dt_write'] # CREATING THE FILES xdmf_u = XDMFFile(options['savepath']+'u.xdmf') @@ -53,7 +51,6 @@ def solv_NavierStokes(options): BCS = [] bc = options['boundary_conditions'] # For Windkessel implementation - flows = {} pii0 = {} pii = {} press = {} @@ -63,6 +60,23 @@ def solv_NavierStokes(options): 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': @@ -82,43 +96,27 @@ def solv_NavierStokes(options): windkessel = True if rank==0: print('Adding Windkessel BC at boundary ',nid) - [R_p,C,R_d] = bc[nbc]['value'] + [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: - flows[nid] = Constant(0) - press[nid] = Constant(R_p*0) + print('no capacitance C for boundary',nid) + press[nid] = Constant(0) else: - flows[nid] = Constant(0) pii0[nid] = Constant(bc[nbc]['p0'][0]*bc[nbc]['p0'][1]) - pii[nid] = Constant(alpha[nid]*pii0[nid] + beta[nid]*flows[nid]) - press[nid] = Constant(gamma[nid]*flows[nid] + alpha[nid]*pii0[nid]) + pii[nid] = Constant(bc[nbc]['p0'][0]*bc[nbc]['p0'][1]) + press[nid] = Constant(bc[nbc]['p0'][0]*bc[nbc]['p0'][1]) - Windkvar[nid] = dt*press[nid]*dot(v,n)*ds(nid) + Windkvar[nid] = press[nid]*dot(v,n)*ds(nid) - 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 - ) - # Stabilization Terms if options['stabilization']['temam']: if rank==0: - print('Addint Temam stabilization term') + print('Adding Temam stabilization term') F += 0.5*rho*div(u0)*dot(u_, v)*dx if pspg: @@ -132,7 +130,7 @@ def solv_NavierStokes(options): for nid in options['stabilization']['forced_normal']['boundaries']: if rank==0: print('Forcing normal velocity in border ', nid) - ut = u - n*dot(u,n) + ut = u - n*dot(u0,n) vt = v - n*dot(v,n) F += gparam*dot(ut,vt)*ds(nid) @@ -145,52 +143,67 @@ def solv_NavierStokes(options): print('adding backflow stabilization in border number:' + str(nk)) F -= 0.5*rho*abs_n(dot(u0, n))*dot(u_, v)*ds(nk) - a = lhs(F) - L = rhs(F) if windkessel: for nid in Windkvar.keys(): - L = L - Windkvar[nid] + 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['checkpoint_dt']/options['dt_write']) + writecicle = int(options['xdmf_dt']/options['dt']) while t<=Tf+dt: - if windkessel: - for k in flows.keys(): - flows[k].assign(assemble(dot(u0,n)*ds(k))) - pii0[k].assign(pii[k]) - pii[k].assign(alpha[k]*pii0[k] + beta[k]*flows[k]) - press[k].assign(gamma[k]*flows[k] + alpha[k]*pii0[k]) - # To solve assemble(a, tensor=A) b = assemble(L) [bcs.apply(A, b) for bcs in BCS] - solve(A, w.vector(), b , 'lu') + 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 = u.function_space().mesh().mpi_comm() - inout.write_HDF5_data(comm, checkpath + '/u.h5', u, '/u', t=t) - inout.write_HDF5_data(comm, checkpath + '/p.h5', p, '/p', t=t) + 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) - inflow.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__': diff --git a/input/Graphics_input.yaml b/input/Graphics_input.yaml index a57b9da..1d99157 100755 --- a/input/Graphics_input.yaml +++ b/input/Graphics_input.yaml @@ -3,12 +3,6 @@ # Input file for Graphics # ################################################# -ppecol: 'coral' -stecol: 'yellowgreen' -daecol: 'cyan' -corcol: 'indigo' -refcol: 'black' - dP: apply: false data: '/home/yeye/Desktop/dP/results/11AoCoPhantomRest2/R1/testBF/' @@ -23,7 +17,7 @@ dP: save: True name: '/home/yeye/Desktop/11AoCoPhantomRest2_BF.png' -Histograms: +Histograms_meshes: apply: false outpath: '/home/yeye/Desktop/' colors: @@ -44,20 +38,45 @@ HistCorrector: 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 - norm: 2 - folder: '/home/yeye/Desktop/First/' - resol: ['H2'] - mode: ['meas','corrs'] - times: ['30ms'] - colors: ['orangered','lime','dodgerblue'] - outpath: '/home/yeye/Desktop/' + 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/' -Components: +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 - folder: '/home/yeye/Desktop/Corrector_2019/figures/Errorcurves/' - subfolders: ['11mm','13mm','Normal'] - type: ['mean','max'] - colors: ['orangered','lime','dodgerblue'] - outpath: '/home/yeye/Desktop/' \ No newline at end of file + 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 index 5efb529..c44d2c3 100755 --- a/input/PostCheck_input.yaml +++ b/input/PostCheck_input.yaml @@ -3,63 +3,95 @@ # Input file for the checkpoint postprocessing # ################################################# -Corrector: +Algebra: apply: false - mode: 'colormap' # h5 , histogram, colormap - outpath: '/home/yeye/Desktop/Corrector_2019/figure3/11mm/' - mesh_path: '/home/yeye/Desktop/Corrector_2019/meshes/11AoCoPhantomRest2.h5' - u_path: '/home/yeye/Desktop/Corrector_2019/First/meas/11mmAoCo2/' - w_path: '/home/yeye/Desktop/Corrector_2019/First/corrector/11mmAoCo2/' - wname: 'w_COR_impl_stan' + 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' - dt: 0.03072 + wname: 'w_COR_impl_stan' Velocity: apply: false - mesh_path: '/home/yeye/Desktop/Corrector_2019/meshes/NormalAoCoPhantomRest2.h5' - checkpoint: '/home/yeye/Desktop/Corrector_2019/First/meas/NormalAoCo2/' + meshpath: '/home/yeye/Desktop/Corrector_2019/Poiseuille/poiseuille.h5' + checkpoint: '/home/yeye/Desktop/Corrector_2019/Poiseuille/COR/SNR15V120/' undersampling: 1 - dt: 0.03072 - filename: 'u' + 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' - outpath: '/home/yeye/Desktop/dP/results/11AoCoPhantomRest2/R1/testBF/' + checkpath: '/home/yeye/Desktop/Corrector_2019/AoCo/9mm_pspg/' method: spheres # slices, boundaries, spheres dt: 0.03072 boundaries: [2,3] spheres: - - center: [0.09, 0.093, 0.08] # 11mm - #- center : [0.08, 0.09, 0.08] - #- center : [0.093, 0.09, 0.085] # Normal + - 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.134, 0.135, 0.026] # 11mm - #- center : [0.125, 0.135, 0.025] # 13 mm - #- center : [0.13, 0.14, 0.03] # Normal + - 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/NavierStokes/examples/meshes/aorta.h5' - checkpoint: '/home/yeye/NavierStokes/examples/results/aorta/' + 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/NavierStokes/examples/results/aorta/' + out_path: '/home/yeye/Desktop/aorta/results/mono/' Error-curves: apply: true - dt: 0.03072 - type: ['mean','max'] - meshpath: '/home/yeye/Desktop/Corrector_2019/meshes/NormalAoCoPhantomRest1.4.h5' - ref_check: '/home/yeye/Desktop/Corrector_2019/I/meas/NormalRest1.4/R1/' + 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/Corrector_2019/I/corrector/NormalRest1.4/' + com_check: '/home/yeye/Desktop/Poiseuille/Corrector/leo2mm/SNR10V120/' com_name: 'w_COR_impl_stan' - outpath: '/home/yeye/Desktop/Corrector_2019/figures/Errorcurves/Normal/' + 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 @@ -77,4 +109,15 @@ Temporal-Interpolation: xdmf: true dt: 0.03072 dt_new: 0.015 - out_check: '/home/yeye/Desktop/Corrector_2019/I/meas/11mmRest1.4/dt15ms/' \ No newline at end of file + 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/cfd_input.yaml b/input/cfd_input.yaml deleted file mode 100755 index c7ef27c..0000000 --- a/input/cfd_input.yaml +++ /dev/null @@ -1,50 +0,0 @@ -################################################# -# -# Input file for CFD Monolithic Solver -# -################################################# -mesh_path : '/home/p283370/Desktop/PhD/AORTA/MESH/aorta_fine/aorta_fine_marked.h5' -write_xdmf : False -savepath : '/home/yeye/Desktop/PhD/AORTA/DATA/monolithic/' -write_checkpoint : False -checkpoint_path : '/home/yeye/Desktop/PhD/AORTA/DATA/monolithic/checkpoint/' -write_outlets : False -outlets_path : '/home/yeye/Desktop/PhD/AORTA/DATA/monolithic/outlets/' - -Tf : 0.9 # seg -density : 1.2 # gr/cm3 -dynamic_viscosity : 0.035 -dt : 0.005 -dt_write : 0.01 -checkpoint_dt : 0.02 - - -stab: - pspg : 1 - epsilon : 0.01 - backflow : 1 - back_id : [3,4,5,6] - temam : 1 - supg : 0 - pssu : 0 - -param: - U : 60 - period : 0.7 - Nvel : 1 - Npress : 1 - theta : 1 - -windkessel: - id : [3,4,5,6] - R_p : {'3':100 , '4':250, '5':250 , '6':250 } - R_d : {'3':1000, '4':8000, '5':8000 , '6':8000} - C : {'3':0.01, '4':0.0001, '5':0.00001, '6':0.0001} - alpha : {'3':0 , '4':0 , '5':0 , '6':0} - beta : {'3':0 , '4':0 , '5':0 , '6':0} - gamma : {'3':0 , '4':0 , '5':0 , '6':0} - - - - - diff --git a/input/mri_input.yaml b/input/mri_input.yaml index ffbb89b..a2bf246 100755 --- a/input/mri_input.yaml +++ b/input/mri_input.yaml @@ -5,42 +5,43 @@ ################################################# phase_contrast: - apply : False - dealiased : True - VENC : 400 - meas0 : '/home/yeye/Desktop/PhD/AORTA/DATA/sequences/9AoCoPhantom0.9/Mag9AoCo0.90.npz' - measG : '/home/yeye/Desktop/PhD/AORTA/DATA/sequences/9AoCoPhantom0.9/Mag9AoCo0.9G.npz' - output : '/home/yeye/Desktop/PhD/AORTA/DATA/sequences/9AoCoPhantom0.9/u_R1.npz' + apply: False + dealiased: True + VENC: 400 + meas0: '/home/yeye/Desktop/PhD/AORTA/DATA/sequences/9AoCoPhantom0.9/Mag9AoCo0.90.npz' + measG: '/home/yeye/Desktop/PhD/AORTA/DATA/sequences/9AoCoPhantom0.9/Mag9AoCo0.9G.npz' + output: '/home/yeye/Desktop/PhD/AORTA/DATA/sequences/9AoCoPhantom0.9/u_R1.npz' resize: - apply : False - dim : [120,120,84] - loadseq : '/home/yeye/Desktop/PhD/AORTA/DATA/sequences/Uffine/aoreal.npz' - save : True - savepath : '/home/yeye/Desktop/PhD/AORTA/DATA/sequences/Uffine/aoreal_rs.npz' + apply: False + dim: [120,120,84] + loadseq: '/home/yeye/Desktop/PhD/AORTA/DATA/sequences/Uffine/aoreal.npz' + save: True + savepath: '/home/yeye/Desktop/PhD/AORTA/DATA/sequences/Uffine/aoreal_rs.npz' magnetization_CIB: - apply : False - loadpath : '/home/yeye/Desktop/PhD/MEDICAL_DATA/Phantom/With_AoCo_9mm/MATLAB_FILES/' - savepath : '/home/yeye/Desktop/PhD/AORTA/DATA/sequences/9AoCoPhantom0.9/Mag9AoCo0.9.npz' + apply: False + loadpath: '/home/yeye/Desktop/PhD/MEDICAL_DATA/Phantom/With_AoCo_9mm/MATLAB_FILES/' + savepath: '/home/yeye/Desktop/PhD/AORTA/DATA/sequences/9AoCoPhantom0.9/Mag9AoCo0.9.npz' magnetization: - apply : False - loadseq : '/home/yeye/Desktop/PhD/AORTA/DATA/sequences/Uffine/aoreal.npz' - VENC : 250 - save : True - savepath : '/home/yeye/Desktop/PhD/AORTA/DATA/sequences/MAG/Mag9AoCoE1.npz' + apply: False + loadseq: '/home/yeye/Desktop/PhD/AORTA/DATA/sequences/Uffine/aoreal.npz' + VENC: 250 + save: True + savepath: '/home/yeye/Desktop/PhD/AORTA/DATA/sequences/MAG/Mag9AoCoE1.npz' sequence: - apply: False - checkpoint_path: '/home/yeye/Desktop/First/H1/dt5ms/coaorta/' - meshpath: '/home/yeye/Desktop/PhD/AORTA/MESH/coaorta/H1/coaortaH1.h5' + apply: false + checkpoint_path: '/home/yeye/Desktop/Poiseuille/Meas/' + meshpath: '/home/yeye/Desktop/Poiseuille/Meas/poiseuille.h5' sampling_rate: 1 boxtype: 'fix_resolution' - ranges: {'x': [10.8,21] , 'y': [13,19] , 'z': [-0.1,11.5] } - resol: [0.14,0.14,0.14] - boxsize: [80,80,80] - name: 'SH1_5ms' + #ranges: {'x': [10.8,21] , 'y': [13,19] , 'z': [-0.1,11.5] } + ranges: {'x': [-1,1] , 'y': [-1,1] , 'z': [0,6] } + resol: [0.2,0.2,0.2] + boxsize: [100,100,100] + name: 'seq' cs: apply: false @@ -54,14 +55,14 @@ cs: name: 'ucs' kt-BLAST: - apply : False - VENC : 3.5 # cm/s - mode : 'kxky' - R : [4] # Acceleration factors - noise : False - save : True - savepath : '/home/yeye/Desktop/PhD/AORTA/DATA/sequences/9AoCoPhantomRest2E1/' - name : 'KT_kxky' + apply: False + VENC: 3.5 # cm/s + mode: 'kxky' + R: [4] # Acceleration factors + noise: False + save: True + savepath: '/home/yeye/Desktop/PhD/AORTA/DATA/sequences/9AoCoPhantomRest2E1/' + name: 'KT_kxky' reference: apply : False @@ -76,11 +77,6 @@ reference: save : True name : 'ref' -norms: - apply : False - field_path : '/home/yeye/Desktop/PhD/AORTA/DATA/ct/aorta_coarse/' - plot : True - create_checkpoint: apply: false loadseq: '/home/yeye/Desktop/Corrector_2019/Second/Vol2/dicom/mod/' @@ -100,28 +96,15 @@ create_checkpoint: savepath: '/home/yeye/Desktop/Corrector_2019/Second/Vol2/dicom/' xdmf: true -peak_pv: - apply : False - orig_seq : '/home/yeye/Desktop/PhD/AORTA/DATA/sequences/aorta_Ucoarse.npz' - peak_slice : 6 - R : [1,2] - p_comp : 'error' # 'peak' or 'error' - flux_bnd : [2] - N_CS : 1 - mesh_into : 'leo' - save : True - savepath : '/home/yeye/Desktop/PhD/AORTA/DATA/maxpv/Ucoarse_test/' - infile_dp : '/home/yeye/Desktop/dP/input/all_points.yaml' - change_mesh: apply: false mode: 'u' - dt: 0.005 - checkpoint_path: '/home/yeye/Desktop/R1/checkpoint/' + dt: 0.03 + checkpoint_path: '/home/yeye/Desktop/Poiseuille/Meas/checkpoint/' under_rate: 1 - mesh_in: '/home/yeye/Desktop/PhD/AORTA/MESH/coaorta/H1/coaortaH1.h5' - mesh_out: '/home/yeye/Desktop/meshes/leo_1.42.h5' - savepath: '/home/yeye/Desktop/test/' + mesh_in: '/home/yeye/Desktop/Poiseuille/Meas/poiseuille.h5' + mesh_out: '/home/yeye/Desktop/Poiseuille/Meas_leo/poiseuille.h5' + savepath: '/home/yeye/Desktop/Poiseuille/Meas_leo/' xdmf: True SENSE: @@ -133,8 +116,9 @@ SENSE: create_leo: apply: false - velseq: '/home/yeye/Desktop/Corrector_2019/Second/Vol2/dicom/mod_nocoro/u_R1.mat' - resol: [0.09,0.09,0.09] + velseq: '/home/yeye/Desktop/Corrector_2019/Meshes/SH3_2.0_R1.npz' + resol: [0.2,0.2,0.2] + #ranges: {'x': [-1.2,1.2] , 'y': [-1.2,1.2] , 'z': [-0.2,6.2] } ranges: {'x': [10.8,21] , 'y': [13,19] , 'z': [-0.1,11.5] } masked: false segmentation: '/home/yeye/Desktop/PhD/MEDICAL_DATA/Phantom/With_AoCo_9mm/MATLAB_FILES/Segmented_Aorta.mat' @@ -148,13 +132,13 @@ kspace_cib: CIBtoH5: apply: false - data_path: '/home/yeye/Desktop/PhD/MEDICAL_DATA/AoCo/9AoCoPhantomRest0.9/' - interpolate: true - flip: false - dt: 0.03072 - mesh_path: '/home/yeye/Desktop/PhD/MEDICAL_DATA/AoCo/9AoCoPhantomRest1.4/' - 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] - outpath: '/home/yeye/Desktop/9mmRest1.4/R1/' + data_path: '/home/yeye/Desktop/PhD/MEDICAL_DATA/Study_David/Patients/Study_14_GM/velmesh_from_velmat/' + interpolate: false + flip: false # Set it True for AoCo = 13mm + dt: 0.03831417624521074 + mesh_path: '/home/yeye/Desktop/PhD/MEDICAL_DATA/Study_David/Patients/Study_14_GM/velmesh_from_velmat/' + 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,26,27] + outpath: '/home/yeye/Desktop/Patient_GM/'