I am migrating all my code from Matlab to Python. One of the first things I need to do is to have a code that plots a sum of multiple Lorentzian functions (the reason are described here). The first thing to do is to create is to ad the following code to our module (MRI.py
`):
# Collection of MRI function for the Cardenas Laboratory
# Spin Echo for VTR signal
import numpy as np
def Spin_Echo(pars,rep_time):
return pars[0]*np.exp(-rep_time/pars[1])
def Lorentz1(p,sat_offset):
a=p[0]; w=pow(p[1],2)/4; off=pow(sat_offset-p[2],2)
return (a*w)/(w+off)
def Lorentzians(pars,sat_offset):
num_of_pools=len(pars)//3
new_pars=np.reshape(pars,(3,num_of_pools))
L=np.zeros(sat_offset.shape)
for i in range(num_of_pools):
L += Lorentz1(new_pars[:,i],sat_offset)
return L
Next, I added the following code to a Python script to call the MRI.py module and plot the results using Matplotlib:
# Cardenas Lab Library
import MRI
#import what I need
import numpy as np
import matplotlib.pyplot as plt
#import seaborn as sn
# generate data for T1 VTR
rep_time=np.linspace(0,5,20)
pars=[100,2.7]
T1signal=MRI.Spin_Echo(pars,rep_time)
# CESTmri
ppm=np.linspace(-10,10,101)
p=[1,0.5,.7 ,2,1,1 ,0,-3,7]
cest=MRI.Lorentzians(p,ppm)
# Plot
fig = plt.figure()
ax1 = fig.add_subplot(221)
ax1.set_xlabel("Time (sec)")
ax1.plot(rep_time,T1signal, 'r--')
ax1.set_ylabel("Signal (%)")
ax1.set_title("Variable TR")
ax2 = fig.add_subplot(222)
ax2.plot(ppm,cest, 'ob-')
ax2.set_xlabel("Saturation Offset (ppm)")
ax2.set_ylabel("Signal (%)")
plt.tight_layout()
ax1.set_title("CEST MRI")
I am data scientist scientist, passionate about helping people using mathematics, programming and chemistry