P283370/NuMRI
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

codes copied

Browse Source
This commit is contained in:
jeremias 2020-01-17 16:24:10 +01:00
parent 14cffd5876
commit 2da253f278
21 changed files with 7186 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

15
codes/.vscode/launch.json vendored Normal file
View File

@ -0,0 +1,15 @@
{
// Use IntelliSense to learn about possible attributes.
// Hover to view descriptions of existing attributes.
// For more information, visit: https://go.microsoft.com/fwlink/?linkid=830387
"version": "0.2.0",
"configurations": [
{
"name": "Python: Current File",
"type": "python",
"request": "launch",
"program": "${file}",
"console": "integratedTerminal"
}
]
}

12
codes/.vscode/tasks.json vendored Normal file
View File

@ -0,0 +1,12 @@
{
// See https://go.microsoft.com/fwlink/?LinkId=733558
// for the documentation about the tasks.json format
"version": "2.0.0",
"tasks": [
{
"label": "echo",
"type": "shell",
"command": "echo Hello"
}
]
}

705
codes/CS.py Normal file
View File

@ -0,0 +1,705 @@
import numpy as np
from numpy import linalg as LA
import sys
from mpi4py import MPI
comm = MPI.COMM_WORLD
size = comm.Get_size()
rank = comm.Get_rank()
# COMPRESSED SENSING: LINEAR BREGMAN METHOD
# Translated and adapted into python from tinycs
#
# *tinycs* is a minimal compressed sensing (CS) toolkit designed
# to allow MR imaging scientists to design undersampled
# acquisitions and reconstruct the resulting data with CS without
# needing to be a CS expert.
#
# The Cartesian reconstruction is based on the split Bregman
# code written by Tom Goldstein, originally available here:
# <http://tag7.web.rice.edu/Split_Bregman.html>
def pdf(k,kw,klo,q):
p = (np.abs(k)/kw)**(-q)
p[np.where(k==0)] = 0
p[np.where(np.abs(k)<=kw)] = 1
p[np.where(k<klo)] = 0
return p
def mask_pdf_1d(n,norm,q,pf):
ks = np.arange(0,n) - np.ceil(n/2) - 1
kmax = np.floor(n/2)
npf = np.round(pf*n)
klo = ks[n-npf]
for k in range(int(kmax)):
P = pdf(ks,k+1,klo,q)
if np.sum(P)>=norm:
break
P = np.fft.fftshift(P)
#if np.mod(n,2)!=0:
# P = np.concatenate(([1],P),axis=None)
return P
def mask_pdf_2d(dims,norm,q,pf):
nz = dims[1]
ny = dims[0]
yc = round(ny/2)
zc = round(nz/2)
rmax = np.sqrt((ny-yc)**2 + (nz-zc)**2)
[Z,Y] = np.meshgrid(np.arange(0,nz),np.arange(0,ny))
RR = np.sqrt( (Y-yc)**2 + (Z-zc)**2)
Z = np.abs(Z - nz/2 - 0.5)
Y = np.abs(Y - ny/2 - 0.5)
for rw in range(1,int(rmax)+1):
P = np.ones([ny,nz])/pf
C = np.logical_and( Z <= rw , Y <= rw)
W = np.logical_or( Z > rw , Y > rw)
P[W] = (RR[W]/rw)**(-q)
if np.sum(P) >= norm:
break
return [P,C]
def GeneratePattern(dim,R):
# 3D CASE
if np.size(dim)==3:
nro = dim[0]
npe = dim[1]
nacq = round(npe/R)
q = 1
pf = 1
P = mask_pdf_1d(npe, nacq, q, pf)
while True:
M = np.random.rand(npe)
M = 1*(M<=P)
if np.sum(M)==nacq:
break
# remove partial Fourier plane and compensate sampling density
M = M!=0
M = np.tile(M,[nro,1]);
#M = M.T
# 4D CASE
if np.size(dim)==4:
nro = dim[0]
npe1 = dim[1]
npe2 = dim[2]
nacq = round(npe1*npe2/R)
q = 1
pf = 1
[P,C] = mask_pdf_2d([npe1,npe2], nacq, q, pf)
RR = np.random.rand(npe1,npe2)
M = (RR <= P)
nchosen = np.sum(M)
if nchosen > nacq: # Correct for inexact number chosen
#outerOn = np.logical_and( M , P!=1 )
outerOn = np.where((M)*(P!=1))
numToFlip = nchosen-nacq
idxs = np.random.permutation(outerOn[0].size)
idxx = outerOn[0][idxs[0:numToFlip]]
idxy = outerOn[1][idxs[0:numToFlip]]
M[idxx,idxy] = False
elif nchosen < nacq:
outerOff = np.where(~M)
idxs = np.random.permutation(outerOff[0].size)
numToFlip = nacq - nchosen
idxx = outerOff[0][idxs[0:numToFlip]]
idxy = outerOff[1][idxs[0:numToFlip]]
M[idxx,idxy] = True
M = np.rollaxis(np.tile(np.rollaxis(M,1),[nro,1,1]),2)
M = np.fft.ifftshift(M)
M = M.transpose((1,0,2))
return M
def get_norm_factor(MASK,uu):
UM = MASK==1
return UM.shape[0]/LA.norm(uu)
def Dxyzt(X):
if np.ndim(X)==3:
dd0 = X[:,:,0]
dd1 = X[:,:,1]
DA = dd0 - np.vstack((dd0[1::,:],dd0[0,:]))
DB = dd1 - np.hstack((dd1[:,1::],dd1[:,0:1]))
return DA + DB
if np.ndim(X)==4:
dd0 = X[:,:,:,0]
dd1 = X[:,:,:,1]
dd2 = X[:,:,:,2]
DA = dd0 - np.vstack((dd0[1::,:,:],dd0[0,:,:][np.newaxis,:,:]))
DB = dd1 - np.hstack((dd1[:,1::,:],dd1[:,0,:][:,np.newaxis,:]))
DC = dd2 - np.dstack((dd2[:,:,1::],dd2[:,:,0][:,:,np.newaxis]))
return DA + DB + DC
def Dxyz(u):
if np.ndim(u)==2:
dx = u[:,:]- np.vstack((u[-1,:],u[0:-1,:]))
dy = u[:,:]- np.hstack((u[:,-1:],u[:,0:-1]))
D = np.zeros([dx.shape[0],dx.shape[1],2],dtype=complex)
D[:,:,0] = dx
D[:,:,1] = dy
return D
if np.ndim(u)==3:
dx = u[:,:,:]- np.vstack((u[-1,:,:][np.newaxis,:,:],u[0:-1,:,:]))
dy = u[:,:,:]- np.hstack((u[:,-1,:][:,np.newaxis,:],u[:,0:-1,:]))
dz = u[:,:,:]- np.dstack((u[:,:,-1][:,:,np.newaxis],u[:,:,0:-1]))
D = np.zeros([dx.shape[0],dx.shape[1],dx.shape[2],3],dtype=complex)
D[:,:,:,0] = dx
D[:,:,:,1] = dy
D[:,:,:,2] = dz
return D
def shrink(X,pgam):
p = 1
s = np.abs(X)
tt = pgam/(s)**(1-p)
# t = pgam/np.sqrt(s)
ss = s-tt
ss = ss*(ss>0)
s = s + 1*(s<tt)
ss = ss/s
return ss*X
def CSMETHOD(ITOT,R):
''' Compressed Function.
Args:
ITOT: a numpy matrix with the full sampled (3D or 4D) dynamical data
R: the acceleration factor
'''
# Method parameters
ninner = 5
nbreg = 10
lmbda = 4
mu = 20
gam = 1
if np.ndim(ITOT)==3:
[row,col,numt2] = ITOT.shape
elif np.ndim(ITOT)==4:
[row,col,dep,numt2] = ITOT.shape
else:
raise Exception('Dynamical data is requested')
MASK = GeneratePattern(ITOT.shape,R)
CS1 = np.zeros(ITOT.shape,dtype=complex)
nit = 0
nit_tot = (numt2-1)/20
if np.ndim(ITOT)==3:
for t in range(numt2):
if rank==0:
print('{3D COMPRESSED SENSING} t = ',t)
Kdata = np.fft.fft2(ITOT[:,:,t])*MASK
data_ndims = Kdata.ndim
mask = Kdata!=0 # not perfect, but good enough
# normalize the data so that standard parameter values work
norm_factor = get_norm_factor(mask, Kdata)
Kdata = Kdata*norm_factor
# Reserve memory for the auxillary variables
Kdata0 = Kdata
img = np.zeros([row,col],dtype=complex)
X = np.zeros([row,col, data_ndims])
B = np.zeros([row,col, data_ndims])
# Build Kernels
scale = np.sqrt(row*col)
murf = np.fft.ifft2(mu*mask*Kdata)*scale
uker = np.zeros([row,col])
uker[0,0] = 4
uker[0,1] = -1 ; uker[1,0] = -1
uker[-1,0] = -1 ; uker[0,-1] = -1
uker = 1/(mu*mask + lmbda*np.fft.fftn(uker) + gam)
# Do the reconstruction
for outer in range(nbreg):
for inner in range(ninner):
# update u
rhs = murf + lmbda*Dxyzt(X-B) + gam*img
img = np.fft.ifft2(np.fft.fft2(rhs)*uker)
# update x and y
A = Dxyz(img) + B
X = shrink(A, 1/lmbda)
# update bregman parameters
B = A - X
Kdata = Kdata + Kdata0 - mask*np.fft.fftn(img)/scale
murf = np.fft.ifftn(mu*mask*Kdata)*scale
# undo the normalization so that results are scaled properly
img = img / norm_factor / scale
CS1[:,:,t] = img
if np.ndim(ITOT)==4:
for t in range(numt2):
if rank==0:
print('[4D CS] R = {re} t = {te}/{tef}'.format(re=R,te=t,tef=numt2))
#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
Kdata_0 = np.fft.fftn(ITOT[:,:,:,t])
Kdata = Kdata_0*MASK
data_ndims = Kdata.ndim
mask = Kdata!=0 # not perfect, but good enough
# normalize the data so that standard parameter values work
norm_factor = get_norm_factor(mask, Kdata)
Kdata = Kdata*norm_factor
# Reserve memory for the auxillary variables
Kdata0 = Kdata
img = np.zeros([row,col,dep],dtype=complex)
X = np.zeros([row,col,dep, data_ndims])
B = np.zeros([row,col,dep, data_ndims])
# Build Kernels
scale = np.sqrt(row*col*dep)
murf = np.fft.ifftn(mu*mask*Kdata)*scale
uker = np.zeros([row,col,dep])
uker[0,0,0] = 8
uker[1,0,0] = -1 ; uker[0,1,0] = -1 ; uker[0,0,1] = -1
uker[-1,0,0] = -1 ; uker[0,-1,0] = -1 ; uker[0,0,-1] = -1
uker = 1/(mu*mask + lmbda*np.fft.fftn(uker) + gam)
# Do the reconstruction
for outer in range(nbreg):
for inner in range(ninner):
# update u
rhs = murf + lmbda*Dxyzt(X-B) + gam*img
img = np.fft.ifft2(np.fft.fft2(rhs)*uker)
# update x and y
A = Dxyz(img) + B
X = shrink(A, 1/lmbda)
# update bregman parameters
B = A - X
Kdata = Kdata + Kdata0 - mask*np.fft.fftn(img)/scale
murf = np.fft.ifftn(mu*mask*Kdata)*scale
# undo the normalization so that results are scaled properly
img = img / norm_factor / scale
CS1[:,:,:,t] = img
return CS1
def CSMETHOD_SENSE(ITOT,R,R_SENSE):
''' Compressed 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
[row,col,dep,numt2] = ITOT.shape
MASK = {}
ITOTCS = {}
MASK[0] = GeneratePattern([row,int(np.ceil(col/2)),dep,numt2],R)
MASK[1] = GeneratePattern([row,int(np.ceil(col/2)),dep,numt2],R)
SenseMAP = {}
[SenseMAP[0],SenseMAP[1]] = Sensitivity_Map([row,col,dep])
col = int(np.ceil(col/2))
ITOTCS[0] = np.zeros([row,col,dep,numt2],dtype=complex)
ITOTCS[1] = np.zeros([row,col,dep,numt2],dtype=complex)
for rs in range(R_SENSE):
for t in range(numt2):
if rank==0:
print('[4D CS] R = {re} t = {te}/{tef}'.format(re=R,te=t,tef=numt2))
Kdata_0 = np.fft.fftn(ITOT[:,:,:,t])
Kdata_0 = Kdata_0*SenseMAP[rs]
Kdata_0 = Kdata_0[:,0::R_SENSE,:]
Kdata = Kdata_0*MASK[rs]
data_ndims = Kdata.ndim
mask = Kdata!=0 # not perfect, but good enough
# normalize the data so that standard parameter values work
norm_factor = get_norm_factor(mask, Kdata)
Kdata = Kdata*norm_factor
# Reserve memory for the auxillary variables
Kdata0 = Kdata
img = np.zeros([row,col,dep],dtype=complex)
X = np.zeros([row,col,dep, data_ndims])
B = np.zeros([row,col,dep, data_ndims])
# Build Kernels
scale = np.sqrt(row*col*dep)
murf = np.fft.ifftn(mu*mask*Kdata)*scale
uker = np.zeros([row,col,dep])
uker[0,0,0] = 8
uker[1,0,0] = -1 ; uker[0,1,0] = -1 ; uker[0,0,1] = -1
uker[-1,0,0] = -1 ; uker[0,-1,0] = -1 ; uker[0,0,-1] = -1
uker = 1/(mu*mask + lmbda*np.fft.fftn(uker) + gam)
# Do the reconstruction
for outer in range(nbreg):
for inner in range(ninner):
# update u
rhs = murf + lmbda*Dxyzt(X-B) + gam*img
img = np.fft.ifft2(np.fft.fft2(rhs)*uker)
# update x and y
A = Dxyz(img) + B
X = shrink(A, 1/lmbda)
# update bregman parameters
B = A - X
Kdata = Kdata + Kdata0 - mask*np.fft.fftn(img)/scale
murf = np.fft.ifftn(mu*mask*Kdata)*scale
# undo the normalization so that results are scaled properly
img = img / norm_factor / scale
ITOTCS[rs][:,:,:,t] = img
return [ITOTCS[0],ITOTCS[1]]
def CSMETHOD_peaksystole(ITOT,R,tstar):
''' Compressed Function.
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
return VENC*param*(np.angle(M1) - np.angle(M0))/np.pi
def GenerateMagnetization(Sq,VENC,noise,scantype='0G'):
# MRI PARAMETERS
gamma = 267.513e6 # rad/Tesla/sec Gyromagnetic ratio for H nuclei
B0 = 1.5 # Tesla Magnetic Field Strenght
TE = 5e-3 # Echo-time
PHASE0 = np.zeros(Sq.shape)
PHASE1 = np.zeros(Sq.shape)
RHO0 = np.zeros(Sq.shape, dtype=complex)
RHO1 = np.zeros(Sq.shape, dtype=complex)
if np.ndim(Sq)==3:
[row,col,numt2] = Sq.shape
[X,Y] = np.meshgrid(np.linspace(0,col,col),np.linspace(0,row,row))
for k in range(numt2):
if noise:
Drho = np.random.normal(0,0.2,[row,col])
Drho2 = np.random.normal(0,0.2,[row,col])
else:
Drho = np.zeros([row,col])
Drho2 = np.zeros([row,col])
varPHASE0 = np.random.randint(-10,11,size=(row,col))*np.pi/180*(np.abs(Sq[:,:,k])<0.001) #Hugo's observation
modulus = 0.5 + 0.5*(np.abs(Sq[:,:,k])>0.001)
if scantype=='0G':
PHASE0[:,:,k] = (gamma*B0*TE+0.01*X)*(np.abs(Sq[:,:,k])>0.001) + 10*varPHASE0
PHASE1[:,:,k] = (gamma*B0*TE+0.01*X)*(np.abs(Sq[:,:,k])>0.001) + 10*varPHASE0 + np.pi*Sq[:,:,k]/VENC
if scantype=='-G+G':
PHASE0[:,:,k] = gamma*B0*TE*np.ones([row,col]) + 10*varPHASE0 - np.pi*Sq[:,:,k]/VENC
PHASE1[:,:,k] = gamma*B0*TE*np.ones([row,col]) + 10*varPHASE0 + np.pi*Sq[:,:,k]/VENC
RHO0[:,:,k] = modulus*np.cos(PHASE0[:,:,k]) + Drho + 1j*modulus*np.sin(PHASE0[:,:,k]) + 1j*Drho2
RHO1[:,:,k] = modulus*np.cos(PHASE1[:,:,k]) + Drho + 1j*modulus*np.sin(PHASE1[:,:,k]) + 1j*Drho2
if np.ndim(Sq)==4:
[row, col, dep, numt2] = Sq.shape
[X, Y, Z] = np.meshgrid(np.linspace(0, col, col), np.linspace(
0, row, row), np.linspace(0, dep, dep))
for k in range(numt2):
if noise:
Drho = np.random.normal(0, 0.2, [row, col, dep])
Drho2 = np.random.normal(0, 0.2, [row, col, dep])
else:
Drho = np.zeros([row, col, dep])
Drho2 = np.zeros([row, col, dep])
varPHASE0 = np.random.randint(-10, 11, size=(row, col, dep)) * \
np.pi/180*(np.abs(Sq[:, :, :, k]) < 0.001)
modulus = 0.5 + 0.5*(np.abs(Sq[:, :, :, k]) > 0.001)
if scantype == '0G':
PHASE0[:, :, :, k] = (gamma*B0*TE+0.01*X) * \
(np.abs(Sq[:, :, :, k]) > 0.001) + 10*varPHASE0
PHASE1[:, :, :, k] = (gamma*B0*TE+0.01*X)*(np.abs(Sq[:, :, :, k])
> 0.001) + 10*varPHASE0 + np.pi*Sq[:, :, :, k]/VENC
if scantype == '-G+G':
PHASE0[:, :, :, k] = gamma*B0*TE * \
np.ones([row, col, dep]) + varPHASE0 - np.pi*Sq[:, :, :, k]/VENC
PHASE1[:, :, :, k] = gamma*B0*TE * \
np.ones([row, col, dep]) + varPHASE0 + np.pi*Sq[:, :, :, k]/VENC
RHO0[:, :, :, k] = modulus*np.cos(PHASE0[:, :, :, k]) + \
Drho + 1j*modulus*np.sin(PHASE0[:, :, :, k]) + 1j*Drho2
RHO1[:, :, :, k] = modulus*np.cos(PHASE1[:, :, :, k]) + \
Drho + 1j*modulus*np.sin(PHASE1[:, :, :, k]) + 1j*Drho2
return [RHO0,RHO1]
def undersampling(Sqx,Sqy,Sqz,options,savepath):
R = options['cs']['R']
for r in R:
if rank==0:
print('Using Acceleration Factor R = ' + str(r))
print('Component x of M0')
[M0,M1] = GenerateMagnetization(Sqx,options['cs']['VENC'],options['cs']['noise'])
print('\n Component x of M0')
M0_cs = CSMETHOD(M0,r)
print('\n Component x of M1')
M1_cs = CSMETHOD(M1,r)
Sqx_cs = phase_contrast(M1_cs,M0_cs,options['cs']['VENC'])
del M0,M1
del M0_cs, M1_cs
[M0,M1] = GenerateMagnetization(Sqy,options['cs']['VENC'],options['cs']['noise'])
print('\n Component y of M0')
M0_cs = CSMETHOD(M0,r)
print('\n Component y of M1')
M1_cs = CSMETHOD(M1,r)
Sqy_cs = phase_contrast(M1_cs,M0_cs,options['cs']['VENC'])
del M0,M1
del M0_cs, M1_cs
[M0,M1] = GenerateMagnetization(Sqz,options['cs']['VENC'],options['cs']['noise'])
if rank==0:
print('\n Component z of M0')
M0_cs = CSMETHOD(M0,r)
if rank==0:
print('\n Component z of M1')
M1_cs = CSMETHOD(M1,r)
if rank==0:
print(' ')
Sqz_cs = phase_contrast(M1_cs,M0_cs,options['cs']['VENC'])
if rank==0:
print('saving the sequences in ' + savepath)
seqname = options['cs']['name'] +'_R' + str(r) + '.npz'
print('sequence name: ' + seqname)
np.savez_compressed( savepath + seqname, x=Sqx_cs, y=Sqy_cs,z=Sqz_cs)
del Sqx_cs,Sqy_cs,Sqz_cs
def undersampling_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')
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)
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']
R_SENSE = 1
if 'R_SENSE' in options['cs']:
R_SENSE = options['cs']['R_SENSE'][0]
for r in R:
if rank==0:
print('Using Acceleration Factor R = ' + str(r))
if R_SENSE==2:
[MxS0_cs,MxS1_cs] = CSMETHOD_SENSE(Mx,r,2)
[MyS0_cs,MyS1_cs] = CSMETHOD_SENSE(My,r,2)
[MzS0_cs,MzS1_cs] = CSMETHOD_SENSE(Mz,r,2)
if rank==0:
print('saving the sequences in ' + savepath)
seqname_s0 = options['cs']['name'] +'S0_R' + str(r) + '.npz'
seqname_s1 = options['cs']['name'] +'S1_R' + str(r) + '.npz'
print('sequence name: ' + seqname_s0)
np.savez_compressed( savepath + seqname_s0, x=MxS0_cs, y=MyS0_cs,z=MzS0_cs)
print('sequence name: ' + seqname_s1)
np.savez_compressed( savepath + seqname_s1, x=MxS1_cs, y=MyS1_cs,z=MzS1_cs)
del MxS0_cs, MyS0_cs, MzS0_cs
del MxS1_cs, MyS1_cs, MzS1_cs
elif R_SENSE==1:
Mx_cs = CSMETHOD(Mx,r)
My_cs = CSMETHOD(My,r)
Mz_cs = CSMETHOD(Mz,r)
if rank==0:
print('saving the sequences in ' + savepath)
seqname = options['cs']['name'] +'_R' + str(r) + '.npz'
print('sequence name: ' + seqname)
np.savez_compressed( savepath + seqname, x=Mx_cs, y=My_cs,z=Mz_cs)
del Mx_cs,My_cs,Mz_cs
else:
raise Exception('Only implemented for 2-fold SENSE!!')
# THE END

1581
codes/Graphics.py Normal file
View File

File diff suppressed because it is too large Load Diff

57
codes/MATLAB/createU.m Normal file
View File

@ -0,0 +1,57 @@
clear all; close all
folder_name = uigetdir([],'Load Folder...');
data = load(strcat(folder_name,'/data.mat'));
SEG = load(strcat(folder_name,'/SEG.mat'));
data = data.data;
SEG = SEG.SEG;
VENC = data.VENC;
VoxelSize = data.voxel_MR;
vel_AP = data.MR_PCA_AP;
vel_RL = data.MR_PCA_RL;
vel_FH = data.MR_PCA_FH;
SEG2 = permute(SEG,[2,3,1]);
SEG2 = SEG2(:,:,:);
vel_AP_seg = vel_AP.*SEG2(2:end-1,2:end-1,2:end-1);
vel_RL_seg = vel_RL.*SEG2(2:end-1,2:end-1,2:end-1);
vel_FH_seg = vel_FH.*SEG2(2:end-1,2:end-1,2:end-1);
u_R1 = [] ;
u_R1.x = vel_FH_seg;
u_R1.y = vel_AP_seg;
u_R1.z = vel_RL_seg;
u_R1.VoxelSize = VoxelSize;
save('/home/yeye/Desktop/u_R1.mat','u_R1');
disp('data saved')
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FIGURES
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure
size_vel = size(vel_FH);
for n=1:size_vel(3)
imshow(squeeze(vel_FH_seg(:,:,n,8)),[-100,100],'InitialMagnification',300);
colormap(gca);
pause(0.1)
end
%%
size_seg2 = size(SEG2);
for n=1:size_seg2(3)
imshow(squeeze(SEG2(:,:,n)),'InitialMagnification',300);
colormap(gca);
pause(0.1)
end

14
codes/MATLAB/leo/CREATE_MESH.m Executable file
View File

@ -0,0 +1,14 @@
% Program to create a structured mesh using the codes of Leo Sok
clear all; close all
nodes = load('LEO_files/nodes.txt');
ux = load('LEO_files/ux.txt') ;
uy = load('LEO_files/uy.txt') ;
uz = load('LEO_files/uz.txt') ;
u = sqrt(ux.^2 + uy.^2 + uz.^2);
resol = load('LEO_files/resol.txt') ;
dx = resol(1); dy = resol(2) ; dz = resol(3);
nodes_masked = maskFEM(nodes,u);
[N,tets,faces] = meshStructTess(nodes_masked,dx,dy,dz,0,0);
writemesh('/home/yeye/Desktop/leomesh',N,tets,faces)

19
codes/MATLAB/leo/maskFEM.m Executable file
View File

@ -0,0 +1,19 @@
function nodes2 = maskFEM(nodes,vel)
a = [];
b = [];
c = [];
ind = 1;
for i=1:length(nodes)
if vel(i)>0
a(ind) = nodes(i,1);
b(ind) = nodes(i,2);
c(ind) = nodes(i,3);
ind = ind +1;
end
end
nodes2 = [a', b', c'];

169
codes/MATLAB/leo/meshStructTess.m Executable file
View File

@ -0,0 +1,169 @@
function [nodes, tets, faces, P] = meshStructTess(nodes, dx, dy, dz, check_mesh, plot_mesh)
%% [nodes, tets, faces] = meshStructTess(nodes, dx, dy, dz, check_mesh, plot_mesh)
% Generate a tessalation from a list of structured nodes.
% input: nodes: n times 3 matrix with on the rows the coordinates of
% the n points in the mesh
% dx, dy, dz: the mesh-size in the directions x, y and z
% check_mesh: if true, then it solves a Poisson problem
% plot_mesh: if true, then it plots the mesh
% output: nodes: m times 3 matrix with on the rows the coordinates of
% the m <= n points in the triangulationedi
% tets: l times 4 matrix with on the rows the tetrahedra
% faces: k times 3 matrix with on the rows the triangles of the
% boundary of the mesh
% P: Transformation matrix from input nodes to output nodes.
% Useful also for transforming node-valued functions on
% the input nodes to node-valued functions on the output
% nodes
%
% The triangulation can be plotted using tetramesh(tets,nodes)
% compute the minimum and number of points in each direction
if size(nodes,1) < 4
error('Triangulation needs at least 4 points')
end
mn = min(nodes);
xmin = mn(1);
ymin = mn(2);
zmin = mn(3);
mn = max(nodes);
xmax = mn(1);
ymax = mn(2);
zmax = mn(3);
nx = round((xmax-xmin)/dx +1);
ny = round((ymax-ymin)/dy +1);
nz = round((zmax-zmin)/dz +1);
Nnodes = size(nodes,1);
% Define tensor which consist of nodes indices, used for the creation of
% the tetrahedra
nodes3d = zeros(nx,ny,nz); % preallocate
for i=1:Nnodes
nodes3d(round((nodes(i,1)-xmin)/dx)+1,round((nodes(i,2)-ymin)/dy)+1,round((nodes(i,3)-zmin)/dz)+1)=i;
end
disp('Creating Tetrahedra')
% create tetrahedral mesh in cube, which we will reuse.
ii = 1;
X = zeros(8,3);
for i=0:1
for j=0:1
for k=0:1
X(ii,:) = [i,j,k];
ii = ii+1;
end
end
end
cubetet = delaunay(X);
% Run through the mesh
el = 1;
Tetrahedra = zeros(6*(nnz(nodes3d)),4); % preallocate
for i=1:nx-1
for j=1:ny-1
for k=1:nz-1
% take [i:i+1,j:j+1,k:k+1] as cube
nod = zeros(1,8); % perallocate
for l = 1:8
% nod is vector with node indices of cube
nod(l) = nodes3d(i + X(l,1), j + X(l,2), k + X(l,3));
end
if nnz(nod) == 8 % then the cube is inside the mesh
tet = nod(cubetet);
else % then there is at least one point of the cube outside the mesh
Xs = X(logical(nod),:); % take only nodes inside the mesh
nodx = nod(logical(nod));
if nnz(nod) == 4 % 4 nodes, check if points are coplanar
C = cross(Xs(2,:)-Xs(1,:), Xs(3,:)-Xs(1,:));
cop = logical(dot(C,Xs(4,:)-Xs(1,:)));
% if cop = 0, then points are coplanar end thus no
% tetrahedra exists.
end
if (nnz(nod)>4) || (nnz(nod) == 4 && cop)
% create tetrahedra
tet1 = delaunay(Xs);
tet = nodx(tet1);
else % no tetrahedra exists
tet = [];
end
end
% add new tetrahedra to list
Tetrahedra(el:el+size(tet,1)-1,:) = tet;
el = el+size(tet,1);
end
end
end
tets = Tetrahedra(1:el-1,:); % Delete extra preallocated rows.
clear Tetrahedra
disp([num2str(size(tets,1)), ' tetrahedra created'])
% Delete nodes which are not in any tetrahedra.
disp('Update mesh')
contr = zeros(size(nodes,1),1);
for i=1:size(tets,1)
for j=1:4
contr(tets(i,j))=1;
end
end
nodes = nodes(logical(contr),:);
% compute P
P = speye(Nnodes);
P = P(logical(contr),:);
disp([num2str(nnz(~contr)), ' unused nodes in triangulation deleted.'])
disp('Update tetrahedra')
% make tetrahedra compatible with new node indices
cumcon = cumsum(~contr)';
tets = tets - cumcon(tets);
% create triangles
if size(tets,1) == 0
warning('No tetrahedra created')
faces = zeros(0,3);
else
disp('Create Triangles')
faces = freeBoundary(triangulation(tets,nodes));
disp([num2str(size(faces,1)), ' triangles created'])
end
% checking the mesh by solving a Poisson problem
if check_mesh
% Builds the P1 stiffness matrix from tets and nodes
[A,volumes]=stifness_matrixP1_3D(tets,nodes);
% Check if element volumes may be negative
if any(volumes<=0)
warning('Some elements have zero or negative volume')
end
% solve the Poisson problem with Dirichlet BC
A(2:end,2:end)\ones(size(A(2:end,2:end),1),1);
disp('If there are no warnings, it probably means that the mesh is fine')
end
% Plots mesh
if plot_mesh
tetramesh(tets,nodes)
xlabel('x')
ylabel('y')
zlabel('z')
end
end

97
codes/MATLAB/leo/writemesh.m Executable file
View File

@ -0,0 +1,97 @@
function writemesh(varargin)
%% writemesh(path, mesh)
% Save triangulation as path.xml and path.msh
% mesh is a struct with fields Pts, Tet, Tri
% alernatively one can use writemesh(path, Pts, Tet, Tri)
% Pts should by a n times 3 matrix consisting points of the mesh
% Tet is the m times 4 matrix consisting the tetrahedra
% Tri is the l times 3 matrix consisting the triangles at the boundary
if nargin > 3
mesh.Pts=varargin{2};
mesh.Tet=varargin{3};
mesh.Tri=varargin{4};
writemesh(varargin{1},mesh,varargin(nargin));
elseif isstruct(varargin{2})
rootMeshFile = varargin{1};
% NEW FILE
obj = [rootMeshFile,'.msh'];
meshfile = fopen(obj,'w');
obj2 = [rootMeshFile,'.xml'];
xmlfile = fopen(obj2,'w');
% MESH
fprintf(meshfile,['$MeshFormat','\n']);
fprintf(meshfile,['2.2 0 8','\n']);
fprintf(meshfile,['$EndMeshFormat','\n']);
fprintf(xmlfile,['<?xml version="1.0" encoding="UTF-8"?>','\n']);
fprintf(xmlfile,'\n');
fprintf(xmlfile,['<dolfin xmlns:dolfin="http://www.fenicsproject.org">','\n']);
mesh = varargin{2};
Nodes = mesh.('Pts');
mesh = rmfield(mesh,'Pts');
Nodes = [(1:size(Nodes,1))' Nodes(:,1:3)];
% POINTS
if ~strcmp(varargin{nargin},'mute')
disp('Write Points')
end
fprintf(meshfile,['$Nodes','\n']);
fprintf(meshfile,['%i','\n'],size(Nodes,1));
fprintf(xmlfile,[' <mesh celltype="tetrahedron" dim="3">','\n']);
fprintf(xmlfile,[' <vertices size="%i">','\n'],size(Nodes,1));
fprintf(meshfile,'%i %13.6f %13.6f %13.6f\n',Nodes');
Nodes(:,1) = Nodes(:,1) - 1;
fprintf(xmlfile,' <vertex index="%i" x="%0.16e" y="%0.16e" z="%0.16e"/>\n',Nodes');
fprintf(meshfile,['$EndNodes','\n']);
fprintf(meshfile,['$Elements','\n']);
fprintf(meshfile,['%i','\n'],size(mesh.Tet,1)+size(mesh.Tri,1));
fprintf(xmlfile,[' </vertices>','\n']);
fprintf(xmlfile,[' <cells size="%i">','\n'],size(mesh.Tet,1));
% Triangles
if ~strcmp(varargin{nargin},'mute')
disp('Write Triangles')
end
tri = mesh.('Tri');
tri = [(1:size(tri,1))' 2*ones(size(tri,1),1) 2*ones(size(tri,1),1) zeros(size(tri,1),1) 2*ones(size(tri,1),1) tri(:,1:3)];
fprintf(meshfile,'%i %i %i %i %i %i %i %i\n',tri');
% Tetrahedra
if ~strcmp(varargin{nargin},'mute')
disp('Write Tetrahedra')
end
tet = mesh.('Tet');
tet = [(size(tri,1)+1:size(tri,1)+size(tet,1))' 4*ones(size(tet,1),1) 2*ones(size(tet,1),1) zeros(size(tet,1),1) ones(size(tet,1),1) tet(:,1:4)];
fprintf(meshfile,'%i %i %i %i %i %i %i %i %i\n',tet');
tet = mesh.('Tet');
tet = [(0:size(tet,1)-1)' (tet(:,1:4)-1)];
fprintf(xmlfile,' <tetrahedron index="%i" v0="%i" v1="%i" v2="%i" v3="%i"/>\n',tet');
fprintf(meshfile,['$EndElements','\n']);
fprintf(xmlfile,' </cells>\n </mesh>\n</dolfin>\n');
fclose('all');
end

126
codes/MATLAB/load_dicom.m Normal file
View File

@ -0,0 +1,126 @@
clear all ; close all
% Load dicom
name = 'Ronald' ;
if strcmp(name, 'Ronald')
path_all = [
'/home/yeye/Desktop/PhD/MEDICAL_DATA/DatosSEPT2019/20190909_Ronald/FH/DICOM/IM_0001',
'/home/yeye/Desktop/PhD/MEDICAL_DATA/DatosSEPT2019/20190909_Ronald/AP/DICOM/IM_0001',
'/home/yeye/Desktop/PhD/MEDICAL_DATA/DatosSEPT2019/20190909_Ronald/RL/DICOM/IM_0001'
] ;
end
if strcmp(name, 'Jeremias')
path_all = [
'/home/yeye/Desktop/PhD/MEDICAL_DATA/DatosSEPT2019/20190909_Jeremias/FH/DICOM/IM_0001',
'/home/yeye/Desktop/PhD/MEDICAL_DATA/DatosSEPT2019/20190909_Jeremias/AP/DICOM/IM_0001',
'/home/yeye/Desktop/PhD/MEDICAL_DATA/DatosSEPT2019/20190909_Jeremias/RL/DICOM/IM_0001'
] ;
end
if strcmp(name, 'Hugo')
path_all = [
'/home/yeye/Desktop/PhD/MEDICAL_DATA/DatosSEPT2019/20190924_Hugo/Dicom/DICOM/IM_0013',
'/home/yeye/Desktop/PhD/MEDICAL_DATA/DatosSEPT2019/20190924_Hugo/Dicom/DICOM/IM_0009',
'/home/yeye/Desktop/PhD/MEDICAL_DATA/DatosSEPT2019/20190924_Hugo/Dicom/DICOM/IM_0005'
] ;
end
for i=1:3
if i==1
%path = '/home/yeye/Desktop/PhD/MEDICAL_DATA/DatosSEPT2019/20190924_Paloma/Dicom/DICOM/IM_0013'
disp('Reading the FH component from ...')
path = path_all(1,:)
end
if i==2
%path = '/home/yeye/Desktop/PhD/MEDICAL_DATA/DatosSEPT2019/20190924_Paloma/Dicom/DICOM/IM_0009' ;
disp('Reading the AP component from ...')
path = path_all(2,:)
end
if i==3
%path = '/home/yeye/Desktop/PhD/MEDICAL_DATA/DatosSEPT2019/20190924_Paloma/Dicom/DICOM/IM_0005' ;
disp('Reading the RL component from ...')
path = path_all(3,:)
end
I_info = dicominfo(path);
I = double(dicomread(path));
VENC = eval(['I_info.PerFrameFunctionalGroupsSequence.Item_1.MRVelocityEncodingSequence.Item_1.VelocityEncodingMaximumValue']) ;
heart_rate = eval(['I_info.PerFrameFunctionalGroupsSequence.Item_1.Private_2005_140f.Item_1.HeartRate']);
MAG = zeros(size(I,1),size(I,2),I_info.Private_2001_1018,I_info.Private_2001_1017);
PHASE = zeros(size(I,1),size(I,2),I_info.Private_2001_1018,I_info.Private_2001_1017);
for n=1:size(I,4)
RI = eval(['I_info.PerFrameFunctionalGroupsSequence.Item_',num2str(n),'.Private_2005_140f.Item_1.RescaleIntercept']); % intercept
RS = eval(['I_info.PerFrameFunctionalGroupsSequence.Item_',num2str(n),'.Private_2005_140f.Item_1.RescaleSlope']); % slope
cp = eval(['I_info.PerFrameFunctionalGroupsSequence.Item_',num2str(n),'.Private_2005_140f.Item_1.Private_2001_1008']); %cp
slc = eval(['I_info.PerFrameFunctionalGroupsSequence.Item_',num2str(n),'.Private_2005_140f.Item_1.Private_2001_100a']); %scl
id = eval(['I_info.PerFrameFunctionalGroupsSequence.Item_',num2str(n),'.Private_2005_140f.Item_1.Private_2005_106e']); % PCA o FFE
if strcmp(id,'FFE')==1
MAG(:,:,slc,cp) = I(:,:,1,n)*RS + RI;
else
PHASE(:,:,slc,cp) = I(:,:,1,n)*RS + RI;
end
end
MASK = double(abs((PHASE==PHASE(1,1,1,1))-1));
PHASE = PHASE.*MASK;
if i==1
MR_FFE_FH = MAG;
MR_PCA_FH = VENC*PHASE/pi/100;
end
if i==2
MR_FFE_AP = MAG;
MR_PCA_AP = VENC*PHASE/pi/100;
end
if i==3
MR_FFE_RL = MAG;
MR_PCA_RL = VENC*PHASE/pi/100;
end
end
disp('Saving the data ...')
spaceslices = eval(['I_info.PerFrameFunctionalGroupsSequence.Item_1.PixelMeasuresSequence.Item_1.SpacingBetweenSlices']);
pixelspacing = eval(['I_info.PerFrameFunctionalGroupsSequence.Item_1.PixelMeasuresSequence.Item_1.PixelSpacing']);
disp('voxel-size recognized:')
voxel_MR = [pixelspacing(1),pixelspacing(1),spaceslices]
data = [];
data.MR_FFE_AP = MR_FFE_AP;
data.MR_FFE_RL = MR_FFE_RL;
data.MR_FFE_FH = MR_FFE_FH;
data.MR_PCA_AP = MR_PCA_AP;
data.MR_PCA_RL = MR_PCA_RL;
data.MR_PCA_FH = MR_PCA_FH;
data.type = 'DAT';
data.VENC = VENC ;
data.voxel_MR = voxel_MR;
data.heart_rate = heart_rate;
save('/home/yeye/Desktop/data.mat','data','-v7.3');
disp('data saved')

1936
codes/MRI.py Normal file
View File

File diff suppressed because it is too large Load Diff

1072
codes/PostCheck.py Normal file
View File

File diff suppressed because it is too large Load Diff

115
codes/SENSE.py Normal file
View File

@ -0,0 +1,115 @@
import numpy as np
from numpy import linalg as LA
import sys
from mpi4py import MPI
comm = MPI.COMM_WORLD
size = comm.Get_size()
rank = comm.Get_rank()
# SENSE: Simulation of SENSitive Encoding algorithm proposed by K. Pruessmann, et. al. in:
# "SENSE: Sensitivity Enconding for Fast MRI" Mag. Res. in Medicine 42. (1999)
# written by Jeremias Garay (j.e.garay.labra@rug.nl)
def Sensitivity_Map(shape):
[Nx,Ny,Nz] = shape
[X,Y,Z] = np.meshgrid(np.linspace(0,Ny,Ny),np.linspace(0,Nx,Nx),np.linspace(0,Nz,Nz))
Xsense1 = (X/(Nx*2)-1)**2
Xsense2 = ((Nx-X)/(Nx*2)-1)**2
S_MAPS = [np.fft.fftshift(Xsense1),np.fft.fftshift(Xsense2)]
return S_MAPS
def SENSE_recon(S1,M1,S2,M2):
[Nx,Ny,Nz,Nt] = M1.shape
M = np.zeros([Nx,int(2*Ny),Nz,Nt],dtype=complex)
sm1 = np.fft.fftshift(S1)[:,:,0]
sm2 = np.fft.fftshift(S2)[:,:,0]
for j in range(Ny):
for k in range(Nx):
l1 = M1[k,j,:,:]; a1 = sm1[k,j]; a2 = sm1[k,j+Ny]
l2 = M2[k,j,:,:]; b1 = sm2[k,j]; b2 = sm2[k,j+Ny]
B = (l1*b1 - l2*a1)/(a2*b1 - b2*a1)
A = (l1*b2 - l2*a2)/(a1*b2 - a2*b1)
M[k,j,:,:] = A
M[k,j+Ny,:,:] = B
return M
def SENSE_recon2(S1,M1,S2,M2):
# With matrices as in the original paper!
[Nx,Ny,Nz,Nt] = M1.shape
M = np.zeros([Nx,int(2*Ny),Nz,Nt],dtype=complex)
sm1 = np.fft.fftshift(S1)[:,:,0]
sm2 = np.fft.fftshift(S2)[:,:,0]
sigma2 = 0.049**2
sigma2 = 1
Psi = np.diagflat(np.array([sigma2,sigma2])) # Error matrix Psi
Psi_inv = np.linalg.inv(Psi)
for j in range(Ny):
for k in range(Nx):
l1 = M1[k,j,:,:]; a1 = sm1[k,j]; a2 = sm1[k,j+Ny]
l2 = M2[k,j,:,:]; b1 = sm2[k,j]; b2 = sm2[k,j+Ny]
S = np.array([[a1,a2],[b1,b2]])
U = np.linalg.inv((np.transpose(S)*Psi_inv*S))*np.transpose(S)*Psi_inv
a = np.array([l1,l2])
a_resized = np.resize(a,(2,Nz*Nt))
v_resized = np.dot(U,a_resized)
v = np.resize(v_resized,(2,Nz,Nt))
M[k,j,:,:] = v[0,:,:]
M[k,j+Ny,:,:] = v[1,:,:]
return M
def SENSE_METHOD(Seq,R):
'''
Args:
ITOT: a numpy matrix with the full sampled (3D or 4D) dynamical data
R: the acceleration factor
'''
[row,col,dep,numt2] = Seq.shape
Seq_red = {}
SenseMAP = {}
[SenseMAP[0],SenseMAP[1]] = Sensitivity_Map([row,col,dep])
col2 = int(np.ceil(col/2))
for rs in range(R):
Seq_red[rs] = np.zeros([row,col2,dep,numt2],dtype=complex)
for t in range(numt2):
Kdata_0 = np.fft.fftn(Seq[:,:,:,t])
Kdata_0 = Kdata_0*SenseMAP[rs]
Kdata_0 = Kdata_0[:,0::R,:]
Seq_red[rs][:,:,:,t] = np.fft.ifftn(Kdata_0)
Seq_recon = SENSE_recon2(SenseMAP[0],Seq_red[0],SenseMAP[1],Seq_red[1])
return Seq_recon
def undersampling(Mx,My,Mz,options):
R = options['SENSE']['R']
for r in R:
if rank==0:
print('Using Acceleration Factor R = ' + str(r))
print('applying into x component')
Mx_s = SENSE_METHOD(Mx,r)
if rank==0:
print('applying into y component')
My_s = SENSE_METHOD(My,r)
if rank==0:
print('applying into z component')
Mz_s = SENSE_METHOD(Mz,r)
return [Mx_s,My_s,Mz_s]

BIN
codes/__pycache__/CS.cpython-36.pyc Normal file
View File

Binary file not shown.

BIN
codes/__pycache__/Graphics.cpython-36.pyc Normal file
View File

Binary file not shown.

BIN
codes/__pycache__/MRI.cpython-36.pyc Normal file
View File

Binary file not shown.

BIN
codes/__pycache__/SENSE.cpython-36.pyc Normal file
View File

Binary file not shown.

BIN
codes/__pycache__/ktBLAST.cpython-36.pyc Normal file
View File

Binary file not shown.

327
codes/cfd_mono.py Normal file
View File

@ -0,0 +1,327 @@
from dolfin import *
import matplotlib.pyplot as plt
import numpy as np
import dolfin
from common import inout
from mpi4py import MPI
import sys
import os
#
# NAVIER STOKES PROBLEM IN THE AORTA with a MONOLITHIC SOLVER
# THIS SCRIPT INCLUDE THE 0-WINDKESSEL BOUNDARY CONDITION
#
# Written by Jeremias Garay L: j.e.garay.labra@rug.nl
#
if '2017' in dolfin.__version__:
class MPI(MPI):
comm_world = mpi_comm_world()
set_log_active(False)
else:
parameters["std_out_all_processes"] = False
def BACH():
import os
tempof = 0.7
semicorchea = 0.1*tempof
corchea = 0.2*tempof
negra = 0.4*tempof
blanca = 0.8*tempof
LA3 = 220
SIb = 233.08
SI = 246.94
DO = 261
REb = 277.18
RE = 293.66
FA = 349.23
MIb = 311.13
MI = 329.63
SOL = 391
LA = 440
os.system('play -V0 --no-show-progress --null --channels 1 synth %s sine %f' % (negra, DO))
os.system('play -V0 --no-show-progress --null --channels 1 synth %s sine %f' % (negra, LA))
os.system('play -V0 --no-show-progress --null --channels 1 synth %s sine %f' % (negra, SOL))
os.system('play -V0 --no-show-progress --null --channels 1 synth %s sine %f' % (negra, FA))
os.system('play -V0 --no-show-progress --null --channels 1 synth %s sine %f' % (semicorchea, MI))
os.system('play -V0 --no-show-progress --null --channels 1 synth %s sine %f' % (semicorchea, FA))
os.system('play -V0 --no-show-progress --null --channels 1 synth %s sine %f' % (corchea, SOL))
os.system('play -V0 --no-show-progress --null --channels 1 synth %s sine %f' % (blanca, FA))
def flux(u,n,ds,idd):
Q = 0
for k in idd:
Q += assemble(dot(u,n)*ds(k))
return Q
def save_outlet_data(options,Qin,Qout3,Qout4,Qout5,Qout6,Pin,Pout3,Pout4,Pout5,Pout6,dt):
# saving the fluxes
np.savetxt(options['outlets_path'] + 'fluxes/' + 'q_in_dt' + str(dt) + '.txt', Qin)
np.savetxt(options['outlets_path'] + 'fluxes/' + 'q3_dt' + str(dt) + '.txt', Qout3)
np.savetxt(options['outlets_path'] + 'fluxes/' + 'q4_dt' + str(dt) + '.txt', Qout4)
np.savetxt(options['outlets_path'] + 'fluxes/' + 'q5_dt' + str(dt) + '.txt', Qout5)
np.savetxt(options['outlets_path'] + 'fluxes/' + 'q6_dt' + str(dt) + '.txt', Qout6)
# saving the pressures
np.savetxt(options['outlets_path'] + 'pressures/' + 'p_in_dt' + str(dt) + '.txt', Pin)
np.savetxt(options['outlets_path'] + 'pressures/' + 'p3_dt' + str(dt) + '.txt', Pout3)
np.savetxt(options['outlets_path'] + 'pressures/' + 'p4_dt' + str(dt) + '.txt', Pout4)
np.savetxt(options['outlets_path'] + 'pressures/' + 'p5_dt' + str(dt) + '.txt', Pout5)
np.savetxt(options['outlets_path'] + 'pressures/' + 'p6_dt' + str(dt) + '.txt', Pout6)
def windkessel_update(u0,n,ds,fluxes,press,pii,pii0,windkessel):
# Updating the time dependent windkessel parameters
if windkessel['C']:
for nk in windkessel['id']:
k = str(nk)
fluxes[k].assign(assemble(dot(u0,n)*ds(nk)))
pii0[k].assign(pii[k])
pii[k].assign(windkessel['alpha'][k]*pii0[k] + windkessel['beta'][k]*fluxes[k])
press[k].assign(windkessel['gamma'][k]*fluxes[k] + windkessel['alpha'][k]*pii0[k])
else:
for nk in windkessel['id']:
k = str(nk)
fluxes[k] = assemble(dot(u0,n)*ds(nk))
press[k].assign( windkessel['R_p'][k]*assemble(dot(u0,n)*ds(nk)))
def solv_NavierStokes(options):
# Assign physical parameters
rho = Constant(options['density'])
mu = Constant(options['dynamic_viscosity'])
otheta = Constant(1-options['param']['theta'])
theta = Constant(options['param']['theta'])
Tf = options['Tf']
dt = options['dt']
dt_w = options['dt_write']
Qin = np.zeros([int(Tf/dt)])
Qout3 = np.zeros([int(Tf/dt)])
Qout4 = np.zeros([int(Tf/dt)])
Qout5 = np.zeros([int(Tf/dt)])
Qout6 = np.zeros([int(Tf/dt)])
Pin = np.zeros([int(Tf/dt)])
Pout3 = np.zeros([int(Tf/dt)])
Pout4 = np.zeros([int(Tf/dt)])
Pout5 = np.zeros([int(Tf/dt)])
Pout6 = np.zeros([int(Tf/dt)])
barye2mmHg = 1/1333.22387415
# CREATING THE FILES
xdmf_u = XDMFFile(options['savepath']+'u.xdmf')
xdmf_p = XDMFFile(options['savepath']+'p.xdmf')
xdmf_u.parameters['rewrite_function_mesh'] = False
xdmf_p.parameters['rewrite_function_mesh'] = False
# LOADING THE MESH
#mesh = Mesh(options['mesh_path'])
#boundaries = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
#boundaries = MarkBoundaries(boundaries)
mesh = Mesh()
hdf = HDF5File(mesh.mpi_comm(), options['mesh_path'] , 'r')
hdf.read(mesh, '/mesh', False)
boundaries = MeshFunction('size_t', mesh , mesh.topology().dim() - 1)
hdf.read(boundaries, '/boundaries')
# To save the boundaries information
#path2 = '/home/p283370/Desktop/marked_fine/boundaries.xdmf'
#XDMFFile(path2).write(boundaries)
# DEFINE FUNCTION SPACES
V = VectorElement('Lagrange', mesh.ufl_cell(), options['param']['Nvel'])
Q = FiniteElement('Lagrange', mesh.ufl_cell(), options['param']['Npress'])
TH = V * Q
W = FunctionSpace(mesh, TH)
# No-slip boundary condition for velocity on walls
noslip = Constant((0, 0, 0))
inflow = Expression(('0','0','(-0.5*U*fabs(sin(w*t)) - 0.5*U*sin(w*t))*(t<0.7)'), degree=3,t=0,U=options['param']['U'], w=2*DOLFIN_PI/options['param']['period'])
bc_inflow = DirichletBC(W.sub(0), inflow, boundaries, 2)
bc_walls = DirichletBC(W.sub(0), noslip, boundaries, 1)
# COLLECTING THE BOUNDARY CONDITIONS
BCS = [bc_inflow,bc_walls]
u, p = TrialFunctions(W)
v, q = TestFunctions(W)
w = Function(W)
n = FacetNormal(mesh)
ds = Measure("ds", subdomain_data=boundaries)
h = CellDiameter(mesh)
f = Constant((0,0,0))
u0, p0 = w.split()
u_ = theta*u + (1 - theta)*u0
theta_p = theta
p_ = theta_p*p + (1 - theta_p)*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
)
if options['stab']['temam']:
F += 0.5*rho*div(u0)*dot(u_, v)*dx
if options['stab']['pspg']:
eps = Constant(options['stab']['epsilon'])
F += eps/mu*h**2*inner(grad(p_), grad(q))*dx
if options['stab']['backflow']:
def abs_n(x):
return 0.5*(x - abs(x))
for nk in options['stab']['back_id']:
if rank==0:
print('adding backflow stabilization in border number:' + str(nk))
F -= 0.5*rho*abs_n(dot(u0, n))*dot(u_, v)*ds(nk)
a = lhs(F)
L = rhs(F)
# Initialization of Windkessel Boundaries
if options['windkessel']['id']:
windkessel = options['windkessel']
# Coeficients
fluxes = {str(k):[] for k in windkessel['id']}
pii0 = {str(k):[] for k in windkessel['id']}
pii = {str(k):[] for k in windkessel['id']}
press = {str(k):[] for k in windkessel['id']}
if windkessel['C']:
for nk in windkessel['id']:
k = str(nk)
if rank==0:
print('Capacitance of windkessel model is: ', windkessel['C'][k])
# computing coeficients
windkessel['alpha'][k] = windkessel['R_d'][k]*windkessel['C'][k]/(windkessel['R_d'][k]*windkessel['C'][k] + options['dt'])
windkessel['beta'][k] = windkessel['R_d'][k]*(1-windkessel['alpha'][k])
windkessel['gamma'][k] = windkessel['R_p'][k] + windkessel['beta'][k]
if rank==0:
print('Using 0-Windkessel complete at outlet number: ' + str(k))
# setting initial values for flux, distal pressure and proximal pressure
fluxes[k] = Constant(0)
pii0[k] = Constant(47/barye2mmHg)
pii[k] = Constant(windkessel['alpha'][k]*pii0[k] + windkessel['beta'][k]*fluxes[k])
press[k] = Constant(windkessel['gamma'][k]*fluxes[k] + windkessel['alpha'][k]*pii0[k])
# Adding to RHS
L = L - dt*press[k]*dot(v,n)*ds(nk)
else:
for nk in windkessel['id']:
k = str(nk)
if rank==0:
print('Using 0-Windkessel reduced at outlet number: ' + str(nk))
fluxes[k] = Constant(0)
press[k] = Constant(windkessel['R_p'][k]*0)
# Adding to RHS
L = L - dt*press[k]*dot(v,n)*ds(nk)
# The static part of the matrix
A = assemble(a)
#b = assemble(L)
[bc.apply(A) for bc in BCS]
#[bc.apply(b) for bc in BCS]
#solv = LUSolver()
#solv.set_operator(A)
#solv.parameters['linear_solver'] = 'mumps'
#solv.parameters['reuse_factorization'] = True
u, p = w.split()
u.rename('velocity', 'u')
p.rename('pressure', 'p')
ind = 0
t = dt
ones = interpolate(Constant(1),W.sub(1).collapse())
A2 = assemble(ones*ds(2))
A3 = assemble(ones*ds(3))
A4 = assemble(ones*ds(4))
A5 = assemble(ones*ds(5))
A6 = assemble(ones*ds(6))
checkcicle = int(options['checkpoint_dt']/options['dt'])
writecicle = int(options['checkpoint_dt']/options['dt_write'])
while t<=Tf+dt:
if options['windkessel']['id']:
windkessel_update(u,n,ds,fluxes,press,pii,pii0,windkessel)
# To solve
assemble(a, tensor=A)
b = assemble(L)
[bc.apply(A, b) for bc in BCS]
solve(A, w.vector(), b)
Qin[ind] = flux(u,n,ds,[2])
Qout3[ind] = flux(u,n,ds,[3])
Qout4[ind] = flux(u,n,ds,[4])
Qout5[ind] = flux(u,n,ds,[5])
Qout6[ind] = flux(u,n,ds,[6])
Pin[ind] = barye2mmHg*assemble(p*ds(2))/A2
Pout3[ind] = barye2mmHg*assemble(p*ds(3))/A3
Pout4[ind] = barye2mmHg*assemble(p*ds(4))/A4
Pout5[ind] = barye2mmHg*assemble(p*ds(5))/A5
Pout6[ind] = barye2mmHg*assemble(p*ds(6))/A6
if rank==0:
print('t = ',t)
# print('|u|:', norm(u0))
# print('|p|:', norm(p0))
# print('div(u):', assemble(div(u0)*dx))
print('Dp = ',np.round(Pin[ind]-Pout3[ind],3),'mmHg')
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)
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)
inflow.t = t
t += dt
# saving the data at outlets: fluxes and pressures
if options['write_outlets']:
if rank==0:
save_outlet_data(options,Qin,Qout3,Qout4,Qout5,Qout6,Pin,Pout3,Pout4,Pout5,Pout6,options['dt'])
if __name__ == '__main__':
comm = MPI.COMM_WORLD
size = comm.Get_size()
rank = comm.Get_rank()
if len(sys.argv) > 1:
if os.path.exists(sys.argv[1]):
inputfile = sys.argv[1]
if rank==0:
print('Found input file ' + inputfile)
else:
raise Exception('Command line arg given but input file does not exist:'
' {}'.format(sys.argv[1]))
else:
raise Exception('An input file is required as argument!')
if 'Zion' in os.popen('hostname').read():
user = 'yeye'
np.set_printoptions(threshold=5)
if 'fwn-bborg-5-166' in os.popen('hostname').read():
user = 'p283370'
if rank==0:
print('Welcome user {uss}'.format(uss=user))
options = inout.read_parameters(inputfile)
solv_NavierStokes(options)

870
codes/ktBLAST.py Normal file
View File

@ -0,0 +1,870 @@
import numpy as np
import scipy as sc
from scipy import signal
from mpi4py import MPI
comm = MPI.COMM_WORLD
size = comm.Get_size()
rank = comm.Get_rank()
# kt-BLAST (NO DC TERM) method for reconstruction of undersampled MRI image based on
# l2 minimization.
def EveryAliased3D2(i,j,k,PP,Nx,Ny,Nz,BB,R):
ivec = [i,j,k]
Nvec = [Nx,Ny,Nz]
[ktot,ltot] = PP.shape
Ptot = np.zeros([ktot**ltot,ltot])
PP2 = np.zeros([ktot**ltot,ltot])
tt = -1
for kk in range(Ptot.shape[0]):
nn = int(np.mod(kk,3))
mm = int(np.mod(np.floor(kk/3),3))
if np.mod(kk,9)==0:
tt+=1
Ptot[kk,0] = PP[tt,0] + ivec[0]
Ptot[kk,1] = PP[mm,1] + ivec[1]
Ptot[kk,2] = PP[nn,2] + ivec[2]
for kk in range(Ptot.shape[0]):
for ll in range(Ptot.shape[1]):
if Ptot[kk,ll]<0:
Ptot[kk,ll] = Ptot[kk,ll] + Nvec[ll]
if Ptot[kk,ll]>=Nvec[ll]:
Ptot[kk,ll] = Ptot[kk,ll] - Nvec[ll]
CC = np.zeros([3,Ptot.shape[0]+1])
YY = np.array([ [i] , [j], [k] ])
CC[0,0] = i
CC[1,0] = j
CC[2,0] = k
psel = 0
for l in range(1,Ptot.shape[0]+1):
CC[0,l] = int(Ptot[l-1,0])
CC[1,l] = int(Ptot[l-1,1])
CC[2,l] = int(Ptot[l-1,2])
if CC[0,l]==YY[0,psel] and CC[1,l]==YY[1,psel] and CC[2,l]==YY[2,psel] and BB[int(CC[1,l]),int(CC[2,l]),int(CC[0,l])]!=0:
pass
else:
War = False
for ww in range(psel):
if CC[0,l]==YY[0,ww] and CC[1,l]==YY[1,ww] and CC[2,l]==YY[2,ww] and BB[int(CC[1,l]),int(CC[2,l]),int(CC[0,l])]!=0:
War = True
if not War:
psel += 1
CCC = np.array([ [CC[0,l] ] , [CC[1,l]] , [CC[2,l]]])
YY = np.concatenate( ( YY, CCC ) ,axis=1 )
return YY.astype(int)
def EveryAliased3D(i,j,k,DP,Nx,Ny,Nz,BB,R,SPREAD=None):
ivec = [i,j,k]
Nvec = [Nx,Ny,Nz]
[ktot,ltot] = DP.shape
DPN = np.zeros([ktot,ltot])
if SPREAD is not None: # WITH SPREAD FUNCTIONS FORMALISM
Maux = np.zeros([Ny,Nz,Nx])
Maux[j,k,i] = 1
SP2 = SPREAD[::-1,::-1,::-1]
MS = R*sc.signal.convolve(Maux,SP2, mode='same')
ms = np.abs(MS)
Ims = 1*(ms>np.max(ms)*0.405)
Pas = np.where(Ims==1)
PP = np.array(Pas[:])
PEA = PP[::-1,:]
for ll in range(PEA.shape[1]):
if PEA[0,ll]>=Nx:
PEA[0,ll] = PEA[0,ll] - Nx
if PEA[1,ll]>=Ny:
PEA[1,ll] = PEA[1,ll] - Ny
if PEA[2,ll]>=Nz:
PEA[2,ll] = PEA[2,ll] - Nz
Ntot = PEA.shape[1]
ind = 0
PEAnew = PEA
for ll in range(Ntot):
if BB[PEA[1,ll],PEA[2,ll],PEA[0,ll]]!=0:
PEAnew = np.delete(PEAnew,(ll-ind),axis=1)
ind +=1
return PEA
else:
for kk in range(DPN.shape[0]):
for l in range(DPN.shape[1]):
DPN[kk,l] = DP[kk,l] + ivec[l]
if DPN[kk,l]<0:
DPN[kk,l] = DPN[kk,l] + Nvec[l]
if DPN[kk,l]>=Nvec[l]:
DPN[kk,l] = DPN[kk,l] - Nvec[l]
CC = np.zeros([3,ktot+1])
YY = np.array([ [i] , [j], [k] ])
CC[0,0] = i
CC[1,0] = j
CC[2,0] = k
for l in range(1,ktot+1):
CC[0,l] = DPN[l-1,0]
CC[1,l] = DPN[l-1,1]
CC[2,l] = DPN[l-1,2]
if CC[0,l]!=CC[0,l-1] and CC[1,l]!=CC[1,l-1] and CC[2,l]!=CC[2,l-1] and BB[int(CC[1,l]),int(CC[2,l]),int(CC[0,l])]==0:
CCC = np.array([ [CC[0,l] ] , [CC[1,l]] , [CC[2,l]]])
YY = np.concatenate( ( YY, CCC ) ,axis=1 )
return YY.astype(int)
def EveryAliased(i,j,DP,Nx,Ny,BB,R,mode):
if mode==1: # USING GEOMETRICAL ASSUMPTIONS
ivec = [i,j]
Nvec = [Nx,Ny]
DPN = 0*DP
[ktot,ltot] = DP.shape
for k in range(ktot):
for l in range(ltot):
DPN[k,l] = DP[k,l] + ivec[l]
if DPN[k,l]<0:
#DPN[k,l] = ivec[l]
DPN[k,l] = DPN[k,l] + Nvec[l]
if DPN[k,l]>=Nvec[l]:
#DPN[k,l] = ivec[l]
DPN[k,l] = DPN[k,l] - Nvec[l]
CC = np.zeros([2,ktot+1])
YY = np.array([ [i] , [j] ])
CC[0,0] = i
CC[1,0] = j
for l in range(1,ktot+1):
CC[0,l] = DPN[l-1,0]
CC[1,l] = DPN[l-1,1]
if CC[0,l]!=CC[0,l-1] and CC[1,l]!=CC[1,l-1] and BB[int(CC[1,l]),int(CC[0,l])]==0:
CCC = np.array([ [CC[0,l] ] , [CC[1,l]] ])
YY = np.concatenate( ( YY, CCC ) ,axis=1 )
return YY.astype(int)
if mode=='spread': # WITH SPREAD FUNCTIONS FORMALISM
Maux = np.zeros([row,numt2])
Maux[l,k] = 1
MS = R*ConvolSP(Maux,SPREAD)
ms = np.abs(MS)
Ims = 1*(ms>np.max(ms)*0.405)
Pas = np.where(Ims==1)
PP = np.array(Pas[:])
return PP[::-1,:]
def GetSymmetric(M):
[row,numt2] = M.shape
S = np.zeros(M.shape,dtype=complex)
aux = np.zeros([1,row])
nmid = 0.5*(numt2+1)
for k in range(int(nmid)):
aux = 0.5*( M[:,k] + M[:,numt2-k-1] )
S[:,k] = aux
S[:,numt2-k-1] = aux
return S
def UNDER(A,mode,R,k):
start = np.mod(k,R)
I1B = np.zeros(A.shape,dtype=complex)
# Not quite efficient ! better to work with vectors
if mode=='ky':
for k in range(start,A.shape[0],R):
I1B[k,:,:] = A[k,:,:]
if mode=='kxky':
for k in range(start,A.shape[0],R):
for l in range(start,A.shape[1],R):
I1B[k,l,:] = A[k,l,:]
if mode=='kxkykz':
for k in range(start,A.shape[0],R):
for l in range(start,A.shape[2],R):
for r in range(start,A.shape[1],R):
I1B[k,r,l] = A[k,r,l]
return I1B
def FilteringHigh(M,fac):
if M.ndim==2:
[row,col] = M.shape
inx = np.linspace(0,col-1,col)
MF = np.zeros(M.shape,dtype=complex)
for k in range(row):
vecfou = np.fft.fft(M[k,:])
window = signal.tukey(2*col,fac)
vecfou2 = vecfou*window[col:2*col]
MF[k,:] = np.fft.ifft(vecfou2)
return MF
if M.ndim==3:
[row,col,dep] = M.shape
MF = np.zeros(M.shape,dtype=complex)
inx = np.linspace(0,col-1,col)
for l in range(dep):
for k in range(row):
vecfou = np.fft.fft(M[k,:,l])
window = signal.tukey(2*col,fac)
vecfou2 = vecfou*window[col:2*col]
MF[k,:,l] = np.fft.ifft(vecfou2)
return MF
def InterpolateM(M,numt2):
if M.ndim==2:
[row,numt] = M.shape
MNew = np.zeros([row,numt2],dtype=complex)
xdat = np.linspace(0,numt,numt)
xdat2 = np.linspace(0,numt,numt2)
nstar = int(0.5*(numt2-numt))
for t in range(nstar,nstar+numt):
MNew[:,t] = M[:,t-nstar]
for l in range(row):
ydat = M[l,:]
fdat = sc.interpolate.interp1d(xdat,ydat,kind='cubic')
MNew[l,1:nstar] = fdat(xdat2)[1:nstar]
MNew[l,nstar+numt:numt2] = fdat(xdat2)[nstar+numt:numt2]
if M.ndim==3:
[row,col,numt] = M.shape
MNew = np.zeros([row,col,numt2],dtype=complex)
xdat = np.linspace(0,numt,numt)
xdat2 = np.linspace(0,numt,numt2)
nstar = int(0.5*(numt2-numt))
for t in range(nstar,nstar+numt):
MNew[:,:,t] = M[:,:,t-nstar]
for c in range(col):
for l in range(row):
ydat = M[l,c,:]
fdat = sc.interpolate.interp1d(xdat,ydat,kind='cubic')
MNew[l,c,1:nstar] = fdat(xdat2)[1:nstar]
MNew[l,c,nstar+numt:numt2] = fdat(xdat2)[nstar+numt:numt2]
return MNew
def KTT(M,scol):
#Maux = M[:,scol,:]
#Maux = np.fft.ifftshift(Maux,axes=0)
#Maux = np.fft.ifft(Maux,axis=0)
#Maux = np.fft.ifft(Maux,axis=1)
#Maux = np.fft.fftshift(Maux,axes=1)
# TAO STYLE
Maux = np.zeros(M.shape,dtype=complex)
for k in range(M.shape[2]):
Maux[:,:,k] = np.fft.ifftshift(M[:,:,k])
Maux[:,:,k] = np.fft.ifft2(Maux[:,:,k])
Maux = Maux[:,scol,:]
Maux = np.fft.ifft(Maux,axis=1)
Maux = np.fft.fftshift(Maux,axes=1)
return Maux
def IKTT(M):
#Maux = np.fft.ifftshift(M,axes=1)
#Maux = np.fft.ifft(Maux,axis=1)
#Maux = np.fft.fft(Maux,axis=0)
#Maux = np.fft.fftshift(Maux,axes=0)
# TAO STYLE
Maux = np.fft.ifftshift(M,axes=1)
Maux = np.fft.fft(Maux,axis=1)
return Maux
def KTT3D(M,sdep):
Maux = np.zeros(M.shape,dtype=complex)
for k in range(M.shape[3]):
Maux[:,:,:,k] = np.fft.ifftshift(M[:,:,:,k])
Maux[:,:,:,k] = np.fft.ifftn(Maux[:,:,:,k])
Maux = Maux[:,:,sdep,:]
Maux = np.fft.ifft(Maux,axis=2)
Maux = np.fft.fftshift(Maux,axes=2)
return Maux
def IKTT3D(M):
Maux = np.fft.ifftshift(M,axes=2)
Maux = np.fft.fft(Maux,axis=2)
return Maux
def get_points4D(row,col,dep,numt2,R,mode):
bb = np.ceil(row/R)
cc = np.ceil(col/R)
aa = np.ceil(numt2/R)
points = R+1
kmid = int(R/2)
if mode=='kxky':
PC = [np.ceil(numt2/2),np.ceil(row/2),np.ceil(col/2)]
PP = np.zeros([points,3])
DP = np.zeros([points-1,3])
for k in range(points):
PP[k,0] = numt2-aa*(k) + 1
PP[k,1] = bb*(k)
PP[k,2] = cc*(k)
if PP[k,0]>=numt2:
PP[k,0] -= 1
if PP[k,0]<0:
PP[k,0] += 1
if PP[k,1]>=row:
PP[k,1] -= 1
if PP[k,1]<0:
PP[k,1] += 1
if PP[k,2]>=col:
PP[k,2] -= 1
if PP[k,2]<0:
PP[k,2] += 1
if k<kmid:
DP[k,0] = PP[k,0] - PC[0]
DP[k,1] = PP[k,1] - PC[1]
DP[k,2] = PP[k,2] - PC[2]
if k>kmid:
DP[k-1,0] = PP[k,0] - PC[0]
DP[k-1,1] = PP[k,1] - PC[1]
DP[k-1,2] = PP[k,2] - PC[2]
kmax = int((PP[kmid,0] + PP[kmid-1,0])/2 )
kmin = int((PP[kmid,0] + PP[kmid+1,0])/2 )
cmax = int((PP[kmid,1] + PP[kmid-1,1])/2 )
cmin = int((PP[kmid,1] + PP[kmid+1,1])/2 )
#DP2 = np.zeros([DP.shape[0]**DP.shape[1],DP.shape[1]])
#DP2[0,0] = DP[0,0]; DP2[0,1] = DP[0,1] ; DP2[0,2] = DP[0,2]
#DP2[1,0] = DP[0,0]; DP2[1,1] = DP[0,1] ; DP2[1,2] = DP[1,2]
#DP2[2,0] = DP[0,0]; DP2[2,1] = DP[1,1] ; DP2[2,2] = DP[0,2]
#DP2[3,0] = DP[0,0]; DP2[3,1] = DP[1,1] ; DP2[3,2] = DP[1,2]
#DP2[4,0] = DP[1,0]; DP2[4,1] = DP[0,1] ; DP2[4,2] = DP[0,2]
#DP2[5,0] = DP[1,0]; DP2[5,1] = DP[0,1] ; DP2[5,2] = DP[1,2]
#DP2[6,0] = DP[1,0]; DP2[6,1] = DP[1,1] ; DP2[6,2] = DP[0,2]
#P2[7,0] = DP[1,0]; DP2[7,1] = DP[1,1] ; DP2[7,2] = DP[1,2]
return [kmin,kmax,PP,DP]
if mode=='ky':
PC = [np.ceil(numt2/2),np.ceil(row/2)]
PP = np.zeros([points,2])
DP = np.zeros([points-1,2])
for k in range(points):
PP[k,0] = numt2-(aa-1)*(k)
PP[k,1] = bb*(k)
if k<kmid:
DP[k,0] = PP[k,0] - PC[0]
DP[k,1] = PP[k,1] - PC[1]
if k>kmid:
DP[k-1,0] = PP[k,0] - PC[0]
DP[k-1,1] = PP[k,1] - PC[1]
kmax = int((PP[kmid,0] + PP[kmid-1,0])/2 )
kmin = int((PP[kmid,0] + PP[kmid+1,0])/2 )
return [kmin,kmax,PP,DP]
def SpreadPoint3D(M,R,sdep):
[row,col,dep,numt2] = M.shape
PS = np.zeros([row,col,dep,numt2],dtype=complex)
inx = 0
iny = 0
for k in range(np.mod(inx,R),row,R):
for ss in range(np.mod(iny,R),col,R):
PS[k,ss,:,:] = 1
iny = iny + 1
inx = inx + 1
for k in range(numt2):
PS[:,:,:,k] = np.fft.ifftn(PS[:,:,:,k])
PS[:,:,:,k] = np.fft.fftshift(PS[:,:,:,k])
SPREAD = PS[:,:,sdep,:]
SPREAD = np.fft.ifft(SPREAD,axis=2)
SPREAD = np.fft.fftshift(SPREAD,axes=2)
return SPREAD
def SpreadPoint(M,R,scol):
[row,col,numt2] = M.shape
PS = np.zeros([row,col,numt2],dtype=complex)
inx = 0
for l in range(0,numt2):
for k in range(np.mod(inx,R),row,R):
PS[k,:,l] = 1
inx = inx + 1
#PS = 1*(M!=0)
#PS = 0*M + 1
#SPREAD = KTT(PS,0)
for k in range(numt2):
PS[:,:,k] = np.fft.ifft2(PS[:,:,k])
PS[:,:,k] = np.fft.fftshift(PS[:,:,k])
SPREAD = PS[:,scol,:]
SPREAD = np.fft.ifft(SPREAD,axis=1)
SPREAD = np.fft.fftshift(SPREAD,axes=1)
return SPREAD
def ConvolSP(M1,M2):
M2 = M2[::-1,::-1]
M3 = sc.signal.convolve2d(M1,M2, boundary='wrap', mode='same')
return M3
def KTBLASTMETHOD_4D_kxky(ITOT,R,mode):
###################################################################
# Training Stage #
###################################################################
[row,col,dep,numt2] = ITOT.shape
# INPUT PARAMETERS
iteshort = 1
tetest = int(dep/2)
numt = int(numt2)
Dyy = int(row*0.1)
rmid = int(row/2)
cmid = int(col/2)
TKdata = np.zeros(ITOT.shape,dtype=complex)
UKdata = np.zeros(ITOT.shape,dtype=complex)
Kdata = np.zeros(ITOT.shape,dtype=complex)
Kdata_NEW0 = np.zeros(ITOT.shape,dtype=complex)
Kdata_NEW = np.zeros(ITOT.shape,dtype=complex)
KTBLAST0 = np.zeros(ITOT.shape,dtype=complex)
KTBLAST = np.zeros(ITOT.shape,dtype=complex)
for k in range(numt2):
# THE FULL KSPACE
Kdata[:,:,:,k] = np.fft.fftn(ITOT[:,:,:,k])
Kdata[:,:,:,k] = np.fft.fftshift(Kdata[:,:,:,k])
# UNDERSAMPLING STEP AND FILLED WITH ZEROS THE REST
UKdata[:,:,:,k] = UNDER(Kdata[:,:,:,k],mode,R,k)
# GENERATING THE TRAINING DATA WITH SUBSAMPLING the Center IN KX , KY
for k in range(numt):
TKdata[rmid-Dyy:rmid+Dyy+1,cmid-Dyy:cmid+Dyy+1,:,k] = Kdata[rmid-Dyy:rmid+Dyy+1,cmid-Dyy:cmid+Dyy+1,:,k]
[kmin,kmax,PP,DP] = get_points4D(row,col,dep,numt2,R,mode)
if iteshort==1:
print(PP)
print(DP)
print('range of k = ',kmin,kmax)
SPREAD = SpreadPoint3D(UKdata,R,tetest)
###################################################################
# RECONSTRUCTION #
###################################################################
ZE1 = iteshort + (tetest-1)*(iteshort)
ZE2 = (tetest+1)*(iteshort) + (dep)*(1-iteshort)
for zi in range(ZE1,ZE2):
if rank==0:
print('4D KTBLAST: R = ' + str(R) + ' and z = ' + str(zi)+'/'+str(dep))
### CONSTRUCT THE REFERENCE M_TRAINING
B2 = KTT3D(TKdata,zi)
B2 = FilteringHigh(B2,0.3)
M2 = 4*np.abs(B2)**2
#M2 = GetSymmetric(M2)
### INTERPOLATE IF NUMT<NUMT2 IS REQUIRED
if numt<numt2:
M2 = InterpolateM(M2,numt2)
#M2 = GetSymmetric(M2)
###################################################################
# Adquisition Stage #
###################################################################
VLund = KTT3D(UKdata,zi)
VLref = KTT3D(Kdata,zi)
VLnew = np.zeros(VLund.shape,dtype=complex)
#for k in range(kmin*0+12,kmax*0+13):
for k in range(kmin,kmax):
#for k in range(0,numt2):
for l in range(row):
for cc in range(col):
PEA = EveryAliased3D2(k,l,cc,PP,numt2,row,col,VLnew,R)
NAA = PEA.shape[1]
Mvec = np.zeros([1,NAA],dtype=complex)
un = np.ones([1,NAA],dtype=complex)
Mvec[0,:] = M2[PEA[1,:],PEA[2,:],PEA[0,:]]
ralia = NAA*VLund[l,cc,k]
MDIAG = np.diagflat(Mvec)
rnew = np.transpose(np.inner(un,MDIAG))/np.sum(np.inner(un,MDIAG))*ralia
VLnew[PEA[1,:],PEA[2,:],PEA[0,:]] = rnew[:,0]
#VLnew[PEA[1,:],PEA[2,:],PEA[0,:]] += 1
KTBLAST[:,:,zi,:] = IKTT3D(VLnew)
KTBLAST = np.fft.ifftshift(KTBLAST,axes=3)
if iteshort==1:
D0 = np.abs(M2)
D1 = np.abs(VLund)
D2 = np.abs(VLnew)
D3 = np.abs(VLref)
return [D0,D1,D2,D3]
if iteshort==0:
return KTBLAST
def KTBLASTMETHOD_4D_ky(ITOT,R,mode):
###################################################################
# Training Stage #
###################################################################
[row,col,dep,numt2] = ITOT.shape
# INPUT PARAMETERS
iteshort = 0
tetest = int(col/2)
zetest = int(3*dep/5)
numt = int(numt2/1)
Dyy = int(row*0.1)
rmid = int(row/2)
cmid = int(col/2)
TKdata = np.zeros(ITOT.shape,dtype=complex)
UKdata = np.zeros(ITOT.shape,dtype=complex)
Kdata = np.zeros(ITOT.shape,dtype=complex)
KTBLAST0 = np.zeros(ITOT.shape,dtype=complex)
KTBLAST = np.zeros(ITOT.shape,dtype=complex)
for k in range(numt2):
# THE FULL KSPACE
Kdata[:,:,:,k] = np.fft.fft2(ITOT[:,:,:,k], axes=(0,1))
Kdata[:,:,:,k] = np.fft.fftshift(Kdata[:,:,:,k])
# UNDERSAMPLING STEP AND FILLED WITH ZEROS THE REST
UKdata[:,:,:,k] = UNDER(Kdata[:,:,:,k],mode,R,k)
# GENERATING THE TRAINING DATA WITH ONLY SUBSAMPLING IN KY
for k in range(numt):
TKdata[rmid-Dyy:rmid+Dyy+1,:,:,k] = Kdata[rmid-Dyy:rmid+Dyy+1,:,:,k]
[kmin,kmax,PP,DP] = get_points4D(row,col,dep,numt2,R,mode)
if iteshort==1:
print(PP)
print(DP)
print('range of y = ',kmin,kmax)
#SPREAD = SpreadPoint(UKdata,R,tetest)
###################################################################
# RECONSTRUCTION #
###################################################################
TE1 = iteshort + (tetest-1)*(iteshort)
TE2 = (tetest+1)*(iteshort) + (col)*(1-iteshort)
ZE1 = iteshort + (zetest-1)*(iteshort)
ZE2 = (zetest+1)*(iteshort) + (dep)*(1-iteshort)
for zi in range(ZE1,ZE2):
if rank==0:
print('4D KTBLAST: z = ',zi)
for te in range(TE1,TE2):
### CONSTRUCT THE REFERENCE M_TRAINING
B2 = KTT(TKdata[:,:,zi,:],te)
B2 = FilteringHigh(B2,0.3)
M2 = 4*np.abs(B2)**2
M2 = GetSymmetric(M2)
### INTERPOLATE IF NUMT<NUMT2 IS REQUIRED
if numt<numt2:
M2 = InterpolateM(M2,numt2)
M2 = GetSymmetric(M2)
###################################################################
# Adquisition Stage #
###################################################################
VLund = KTT(UKdata[:,:,zi,:],te)
VLref = KTT(Kdata[:,:,zi,:],te)
VLnew = np.zeros(VLund.shape,dtype=complex)
for k in range(kmin,kmax):
for l in range(row):
PEA = EveryAliased(k,l,DP,numt2,row,VLnew,R,1)
NAA = PEA.shape[1]
Mvec = np.zeros([1,NAA],dtype=complex)
un = np.ones([1,NAA],dtype=complex)
Mvec[0,:] = M2[PEA[1,:],PEA[0,:]]
ralia = NAA*VLund[l,k]
MDIAG = np.diagflat(Mvec)
rnew = np.transpose(np.inner(un,MDIAG))/np.sum(np.inner(un,MDIAG))*ralia
VLnew[PEA[1,:],PEA[0,:]] = rnew[:,0]
KTBLAST[:,te,zi,:] = IKTT(VLnew)
KTBLAST = np.fft.ifftshift(KTBLAST,axes=2)
if iteshort==1:
D0 = np.abs(M2)
D1 = np.abs(VLund)
D2 = np.abs(VLnew)
D3 = np.abs(VLref)
return [D0,D1,D2,D3]
if iteshort==0:
return KTBLAST
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'):
# MRI PARAMETERS
gamma = 267.513e6 # rad/Tesla/sec Gyromagnetic ratio for H nuclei
B0 = 1.5 # Tesla Magnetic Field Strenght
TE = 5e-3 # Echo-time
M1 = np.pi/(gamma*VENC)
PHASE0 = np.zeros(Sq.shape)
PHASE1 = np.zeros(Sq.shape)
RHO0 = np.zeros(Sq.shape,dtype=complex)
RHO1 = np.zeros(Sq.shape,dtype=complex)
# THE K DATA
KRHO0 = np.zeros(Sq.shape,dtype=complex)
KRHO1 = np.zeros(Sq.shape,dtype=complex)
if np.ndim(Sq)==3:
[row,col,numt2] = Sq.shape
[X,Y] = np.meshgrid(np.linspace(0,col,col),np.linspace(0,row,row))
for k in range(numt2):
if noise:
Drho = np.random.normal(0,0.2,[row,col])
Drho2 = np.random.normal(0,0.2,[row,col])
else:
Drho = np.zeros([row,col])
Drho2 = np.zeros([row,col])
varPHASE0 = np.random.randint(-10,11,size=(row,col))*np.pi/180*(np.abs(Sq[:,:,k])<0.001) #Hugo's observation
modulus = 0.5 + 0.5*(np.abs(Sq[:,:,k])>0.001)
if scantype=='0G':
PHASE0[:,:,k] = (gamma*B0*TE+0.01*X)*(np.abs(Sq[:,:,k])>0.001) + 10*varPHASE0
PHASE1[:,:,k] = (gamma*B0*TE+0.01*X)*(np.abs(Sq[:,:,k])>0.001) + 10*varPHASE0 + np.pi*Sq[:,:,k]/VENC
if scantype=='-G+G':
PHASE0[:,:,k] = gamma*B0*TE*np.ones([row,col]) + 10*varPHASE0 - np.pi*Sq[:,:,k]/VENC
PHASE1[:,:,k] = gamma*B0*TE*np.ones([row,col]) + 10*varPHASE0 + np.pi*Sq[:,:,k]/VENC
RHO0[:,:,k] = modulus*np.cos(PHASE0[:,:,k]) + Drho + 1j*modulus*np.sin(PHASE0[:,:,k]) + 1j*Drho2
RHO1[:,:,k] = modulus*np.cos(PHASE1[:,:,k]) + Drho + 1j*modulus*np.sin(PHASE1[:,:,k]) + 1j*Drho2
if np.ndim(Sq)==4:
[row,col,dep,numt2] = Sq.shape
[X,Y,Z] = np.meshgrid(np.linspace(0,col,col),np.linspace(0,row,row),np.linspace(0,dep,dep))
for k in range(numt2):
if noise:
Drho = np.random.normal(0,0.2,[row,col,dep])
Drho2 = np.random.normal(0,0.2,[row,col,dep])
else:
Drho = np.zeros([row,col,dep])
Drho2 = np.zeros([row,col,dep])
varPHASE0 = np.random.randint(-10,11,size=(row,col,dep))*np.pi/180*(np.abs(Sq[:,:,:,k])<0.001)
modulus = 0.5 + 0.5*(np.abs(Sq[:,:,:,k])>0.001)
if scantype=='0G':
PHASE0[:,:,:,k] = (gamma*B0*TE+0.01*X)*(np.abs(Sq[:,:,:,k])>0.001) + 10*varPHASE0
PHASE1[:,:,:,k] = (gamma*B0*TE+0.01*X)*(np.abs(Sq[:,:,:,k])>0.001) + 10*varPHASE0 + np.pi*Sq[:,:,:,k]/VENC
if scantype=='-G+G':
PHASE0[:,:,:,k] = gamma*B0*TE*np.ones([row,col,dep]) + varPHASE0 - np.pi*Sq[:,:,:,k]/VENC
PHASE1[:,:,:,k] = gamma*B0*TE*np.ones([row,col,dep]) + varPHASE0 + np.pi*Sq[:,:,:,k]/VENC
RHO0[:,:,:,k] = modulus*np.cos(PHASE0[:,:,:,k]) + Drho + 1j*modulus*np.sin(PHASE0[:,:,:,k]) + 1j*Drho2
RHO1[:,:,:,k] = modulus*np.cos(PHASE1[:,:,:,k]) + Drho + 1j*modulus*np.sin(PHASE1[:,:,:,k]) + 1j*Drho2
return [RHO0,RHO1]
def undersampling(Sqx,Sqy,Sqz,options,savepath):
R = options['kt-BLAST']['R']
mode = options['kt-BLAST']['mode']
transpose = True
for r in R:
if rank==0:
print('Using Acceleration Factor R = ' + str(r))
print('Component x of M0')
[M0,M1] = GenerateMagnetization(Sqx,options['kt-BLAST']['VENC'],options['kt-BLAST']['noise'],scantype='0G')
if transpose:
M0 = M0.transpose((0,2,1,3))
M1 = M1.transpose((0,2,1,3))
if mode=='ky':
M0_kt = KTBLASTMETHOD_4D_ky(M0,r,mode)
if mode=='kxky':
M0_kt = KTBLASTMETHOD_4D_kxky(M0,r,mode)
if rank==0:
print('\n Component x of M1')
if mode=='ky':
M1_kt = KTBLASTMETHOD_4D_ky(M1,r,mode)
if mode=='kxky':
M1_kt = KTBLASTMETHOD_4D_kxky(M1,r,mode)
Sqx_kt = phase_contrast(M1_kt,M0_kt,options['kt-BLAST']['VENC'],scantype='0G')
del M0,M1
del M0_kt, M1_kt
[M0,M1] = GenerateMagnetization(Sqy,options['kt-BLAST']['VENC'],options['kt-BLAST']['noise'],scantype='0G')
if transpose:
M0 = M0.transpose((0,2,1,3))
M1 = M1.transpose((0,2,1,3))
if rank==0:
print('\n Component y of M0')
if mode=='ky':
M0_kt = KTBLASTMETHOD_4D_ky(M0,r,mode)
if mode=='kxky':
M0_kt = KTBLASTMETHOD_4D_kxky(M0,r,mode)
if rank==0:
print('\n Component y of M1')
if mode=='ky':
M1_kt = KTBLASTMETHOD_4D_ky(M1,r,mode)
if mode=='kxky':
M1_kt = KTBLASTMETHOD_4D_kxky(M1,r,mode)
Sqy_kt = phase_contrast(M1_kt,M0_kt,options['kt-BLAST']['VENC'],scantype='0G')
del M0,M1
del M0_kt, M1_kt
[M0,M1] = GenerateMagnetization(Sqz,options['kt-BLAST']['VENC'],options['kt-BLAST']['noise'],scantype='0G')
if transpose:
M0 = M0.transpose((0,2,1,3))
M1 = M1.transpose((0,2,1,3))
if rank==0:
print('\n Component z of M0')
if mode=='ky':
M0_kt = KTBLASTMETHOD_4D_ky(M0,r,mode)
if mode=='kxky':
M0_kt = KTBLASTMETHOD_4D_kxky(M0,r,mode)
if rank==0:
print('\n Component z of M1')
if mode=='ky':
M1_kt = KTBLASTMETHOD_4D_ky(M1,r,mode)
if mode=='kxky':
M1_kt = KTBLASTMETHOD_4D_kxky(M1,r,mode)
if rank==0:
print(' ')
Sqz_kt = phase_contrast(M1_kt,M0_kt,options['kt-BLAST']['VENC'],scantype='0G')
if transpose:
Sqx_kt = Sqx_kt.transpose((0,2,1,3))
Sqy_kt = Sqy_kt.transpose((0,2,1,3))
Sqz_kt = Sqz_kt.transpose((0,2,1,3))
if options['kt-BLAST']['save']:
if rank==0:
print('saving the sequences in ' + savepath)
seqname = options['kt-BLAST']['name'] +'_R' + str(r) + '.npz'
print('sequence name: ' + seqname)
np.savez_compressed( savepath + seqname, x=Sqx_kt, y=Sqy_kt,z=Sqz_kt)
del Sqx_kt,Sqy_kt,Sqz_kt

71
codes/mesh_generator.py Executable file
View File

@ -0,0 +1,71 @@
from dolfin import *
mesh_in = '/home/yeye/Desktop/leomesh.xml'
mesh_out = '/home/yeye/Desktop/aorta.h5'
mesh = Mesh(mesh_in)
hdf = HDF5File(mesh.mpi_comm(), mesh_out, 'w')
boundaries = MeshFunction('size_t', mesh,2)
marked = 1
testmesh = 0
hdf.write(mesh, '/mesh')
if marked==1:
class Inlet(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and between(x[0],(0.1975,0.1989)) and between(x[2],(0.07,0.1))
class Outlet(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and between(x[0],(0.1975,0.1989)) and between(x[2],(0,0.04))
class Walls(SubDomain):
def inside(self, x, on_boundary):
return on_boundary
outlet = Outlet()
inlet = Inlet()
walls = Walls()
boundaries.set_all(0)
walls.mark(boundaries,1)
outlet.mark(boundaries,3)
inlet.mark(boundaries,2)
hdf.write(boundaries, '/boundaries')
hdf.close()
if testmesh:
print('Testing Mesh...')
meshname = mesh_out
pathtoB = '/home/yeye/Desktop/boundaries.xdmf'
mesh = Mesh()
hdf = HDF5File(mesh.mpi_comm(), meshname , 'r')
hdf.read(mesh, '/mesh', False)
boundaries = MeshFunction('size_t', mesh , mesh.topology().dim() - 1)
hdf.read(boundaries, '/boundaries')
# To save the boundaries information
XDMFFile(pathtoB).write(boundaries)
print('Boundary info printed in ' + pathtoB)