Browse Source

delete all

master
J.E. Garay Labra 2 years ago
parent
commit
5d6198ee9c
  1. 15
      codes/.vscode/launch.json
  2. 12
      codes/.vscode/tasks.json
  3. 597
      codes/CS.py
  4. 1412
      codes/Graphics.py
  5. 57
      codes/MATLAB/createU.m
  6. 14
      codes/MATLAB/leo/CREATE_MESH.m
  7. 19
      codes/MATLAB/leo/maskFEM.m
  8. 169
      codes/MATLAB/leo/meshStructTess.m
  9. 97
      codes/MATLAB/leo/writemesh.m
  10. 126
      codes/MATLAB/load_dicom.m
  11. 1755
      codes/MRI.py
  12. 1452
      codes/PostCheck.py
  13. 115
      codes/SENSE.py
  14. BIN
      codes/__pycache__/CS.cpython-36.pyc
  15. BIN
      codes/__pycache__/Graphics.cpython-36.pyc
  16. BIN
      codes/__pycache__/MRI.cpython-36.pyc
  17. BIN
      codes/__pycache__/SENSE.cpython-36.pyc
  18. BIN
      codes/__pycache__/ktBLAST.cpython-36.pyc
  19. 870
      codes/ktBLAST.py
  20. 71
      codes/mesh_generator.py
  21. 229
      codes/monolithic.py
  22. 82
      input/Graphics_input.yaml
  23. 123
      input/PostCheck_input.yaml
  24. 138
      input/aorta.yaml
  25. 131
      input/aorta_roukf.yaml
  26. 147
      input/mri_input.yaml

15
codes/.vscode/launch.json

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

12
codes/.vscode/tasks.json

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

597
codes/CS.py

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

1412
codes/Graphics.py
File diff suppressed because it is too large
View File

57
codes/MATLAB/createU.m

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

14
codes/MATLAB/leo/CREATE_MESH.m

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

19
codes/MATLAB/leo/maskFEM.m

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

169
codes/MATLAB/leo/meshStructTess.m

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

97
codes/MATLAB/leo/writemesh.m

@ -1,97 +0,0 @@
function writemesh(varargin)
%% writemesh(path, mesh)
% Save triangulation as path.xml and path.msh
% mesh is a struct with fields Pts, Tet, Tri
% alernatively one can use writemesh(path, Pts, Tet, Tri)
% Pts should by a n times 3 matrix consisting points of the mesh
% Tet is the m times 4 matrix consisting the tetrahedra
% Tri is the l times 3 matrix consisting the triangles at the boundary
if nargin > 3
mesh.Pts=varargin{2};
mesh.Tet=varargin{3};
mesh.Tri=varargin{4};
writemesh(varargin{1},mesh,varargin(nargin));
elseif isstruct(varargin{2})
rootMeshFile = varargin{1};
% NEW FILE
obj = [rootMeshFile,'.msh'];
meshfile = fopen(obj,'w');
obj2 = [rootMeshFile,'.xml'];
xmlfile = fopen(obj2,'w');
% MESH
fprintf(meshfile,['$MeshFormat','\n']);
fprintf(meshfile,['2.2 0 8','\n']);
fprintf(meshfile,['$EndMeshFormat','\n']);
fprintf(xmlfile,['<?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

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

1755
codes/MRI.py
File diff suppressed because it is too large
View File

1452
codes/PostCheck.py
File diff suppressed because it is too large
View File

115
codes/SENSE.py

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

BIN
codes/__pycache__/CS.cpython-36.pyc

BIN
codes/__pycache__/Graphics.cpython-36.pyc

BIN
codes/__pycache__/MRI.cpython-36.pyc

BIN
codes/__pycache__/SENSE.cpython-36.pyc

BIN
codes/__pycache__/ktBLAST.cpython-36.pyc

870
codes/ktBLAST.py

@ -1,870 +0,0 @@
import numpy as np
import scipy as sc
from scipy import signal
from mpi4py import MPI
comm = MPI.COMM_WORLD
size = comm.Get_size()
rank = comm.Get_rank()
# kt-BLAST (NO DC TERM) method for reconstruction of undersampled MRI image based on
# l2 minimization.
def EveryAliased3D2(i,j,k,PP,Nx,Ny,Nz,BB,R):
ivec = [i,j,k]
Nvec = [Nx,Ny,Nz]
[ktot,ltot] = PP.shape
Ptot = np.zeros([ktot**ltot,ltot])
PP2 = np.zeros([ktot**ltot,ltot])
tt = -1
for kk in range(Ptot.shape[0]):
nn = int(np.mod(kk,3))
mm = int(np.mod(np.floor(kk/3),3))
if np.mod(kk,9)==0:
tt+=1
Ptot[kk,0] = PP[tt,0] + ivec[0]
Ptot[kk,1] = PP[mm,1] + ivec[1]
Ptot[kk,2] = PP[nn,2] + ivec[2]
for kk in range(Ptot.shape[0]):
for ll in range(Ptot.shape[1]):
if Ptot[kk,ll]<0:
Ptot[kk,ll] = Ptot[kk,ll] + Nvec[ll]
if Ptot[kk,ll]>=Nvec[ll]:
Ptot[kk,ll] = Ptot[kk,ll] - Nvec[ll]
CC = np.zeros([3,Ptot.shape[0]+1])
YY = np.array([ [i] , [j], [k] ])
CC[0,0] = i
CC[1,0] = j
CC[2,0] = k
psel = 0
for l in range(1,Ptot.shape[0]+1):
CC[0,l] = int(Ptot[l-1,0])
CC[1,l] = int(Ptot[l-1,1])
CC[2,l] = int(Ptot[l-1,2])
if CC[0,l]==YY[0,psel] and CC[1,l]==YY[1,psel] and CC[2,l]==YY[2,psel] and BB[int(CC[1,l]),int(CC[2,l]),int(CC[0,l])]!=0:
pass
else:
War = False
for ww in range(psel):
if CC[0,l]==YY[0,ww] and CC[1,l]==YY[1,ww] and CC[2,l]==YY[2,ww] and BB[int(CC[1,l]),int(CC[2,l]),int(CC[0,l])]!=0:
War = True
if not War:
psel += 1
CCC = np.array([ [CC[0,l] ] , [CC[1,l]] , [CC[2,l]]])
YY = np.concatenate( ( YY, CCC ) ,axis=1 )
return YY.astype(int)
def EveryAliased3D(i,j,k,DP,Nx,Ny,Nz,BB,R,SPREAD=None):
ivec = [i,j,k]
Nvec = [Nx,Ny,Nz]
[ktot,ltot] = DP.shape
DPN = np.zeros([ktot,ltot])
if SPREAD is not None: # WITH SPREAD FUNCTIONS FORMALISM
Maux = np.zeros([Ny,Nz,Nx])
Maux[j,k,i] = 1
SP2 = SPREAD[::-1,::-1,::-1]
MS = R*sc.signal.convolve(Maux,SP2, mode='same')
ms = np.abs(MS)
Ims = 1*(ms>np.max(ms)*0.405)
Pas = np.where(Ims==1)
PP = np.array(Pas[:])
PEA = PP[::-1,:]
for ll in range(PEA.shape[1]):
if PEA[0,ll]>=Nx:
PEA[0,ll] = PEA[0,ll] - Nx
if PEA[1,ll]>=Ny:
PEA[1,ll] = PEA[1,ll] - Ny
if PEA[2,ll]>=Nz:
PEA[2,ll] = PEA[2,ll] - Nz
Ntot = PEA.shape[1]
ind = 0
PEAnew = PEA
for ll in range(Ntot):
if BB[PEA[1,ll],PEA[2,ll],PEA[0,ll]]!=0:
PEAnew = np.delete(PEAnew,(ll-ind),axis=1)
ind +=1
return PEA
else:
for kk in range(DPN.shape[0]):
for l in range(DPN.shape[1]):
DPN[kk,l] = DP[kk,l] + ivec[l]
if DPN[kk,l]<0:
DPN[kk,l] = DPN[kk,l] + Nvec[l]
if DPN[kk,l]>=Nvec[l]:
DPN[kk,l] = DPN[kk,l] - Nvec[l]
CC = np.zeros([3,ktot+1])
YY = np.array([ [i] , [j], [k] ])
CC[0,0] = i
CC[1,0] = j
CC[2,0] = k
for l in range(1,ktot+1):
CC[0,l] = DPN[l-1,0]
CC[1,l] = DPN[l-1,1]
CC[2,l] = DPN[l-1,2]
if CC[0,l]!=CC[0,l-1] and CC[1,l]!=CC[1,l-1] and CC[2,l]!=CC[2,l-1] and BB[int(CC[1,l]),int(CC[2,l]),int(CC[0,l])]==0:
CCC = np.array([ [CC[0,l] ] , [CC[1,l]] , [CC[2,l]]])
YY = np.concatenate( ( YY, CCC ) ,axis=1 )
return YY.astype(int)
def EveryAliased(i,j,DP,Nx,Ny,BB,R,mode):
if mode==1: # USING GEOMETRICAL ASSUMPTIONS
ivec = [i,j]
Nvec = [Nx,Ny]
DPN = 0*DP
[ktot,ltot] = DP.shape
for k in range(ktot):
for l in range(ltot):
DPN[k,l] = DP[k,l] + ivec[l]
if DPN[k,l]<0:
#DPN[k,l] = ivec[l]
DPN[k,l] = DPN[k,l] + Nvec[l]
if DPN[k,l]>=Nvec[l]:
#DPN[k,l] = ivec[l]
DPN[k,l] = DPN[k,l] - Nvec[l]
CC = np.zeros([2,ktot+1])
YY = np.array([ [i] , [j] ])
CC[0,0] = i
CC[1,0] = j
for l in range(1,ktot+1):
CC[0,l] = DPN[l-1,0]