notebook.community

Edit and run 





In [1]:



    


%load_ext autoreload
%autoreload 2
from __future__ import division
# %pylab
import numpy as np
import epgcpmg as epg



    











In [2]:



    


np.set_printoptions(suppress=True, precision=4, linewidth=300)



    











In [35]:



    


T1 = 1000e-3
T2 = 100e-3
TE = 5e-3



P_z = np.array([[0],[0],[1]])
P_xy = np.array([[1],[1],[0]])

def eRF(P, a):
        return epg.rf(P, np.pi/180 * a)
def eTE(P, a):
        return epg.FSE_TE(P, a * np.pi/180, TE, T1, T2)
def eRelax(P, T):
        return epg.relax(P, T, T1, T2)
def eGrad(P):
        return epg.grad(P)



    











In [36]:



    


Q1 = eGrad(P_xy)
print Q1



    


















    






[[0 1]
 [0 0]
 [0 0]]

















In [37]:



    


Q2 = eRF(Q1, 120)
print Q2



    


















    






[[ 0.     0.25 ]
 [ 0.     0.75 ]
 [ 0.    -0.433]]

















In [38]:



    


Q3 = eGrad(Q2)
print Q3



    


















    






[[ 0.75   0.     0.25 ]
 [ 0.75   0.     0.   ]
 [ 0.    -0.433  0.   ]]

















In [39]:



    


Q4 = eGrad(eRF(eGrad(Q3), 120))
print Q4



    


















    






[[ 0.9375  0.     -0.1875  0.      0.0625]
 [ 0.9375  0.      0.1875  0.      0.    ]
 [ 0.     -0.1083  0.     -0.1083  0.    ]]

















In [40]:



    


Q5 = eGrad(eRF(eGrad(Q4), 100))
print Q5



    


















    






[[ 0.7342  0.      0.3908  0.     -0.1841  0.      0.0258]
 [ 0.7342  0.     -0.0034  0.      0.0367  0.      0.    ]
 [ 0.     -0.3505  0.      0.1111  0.     -0.0308  0.    ]]

















In [41]:



    


Q6 = eGrad(eRF(eGrad(Q5), 120))
print Q6



    


















    






[[ 0.8534  0.     -0.1226  0.      0.2214  0.     -0.0727  0.      0.0065]
 [ 0.8534  0.      0.206   0.     -0.1114  0.      0.0194  0.      0.    ]
 [ 0.     -0.1442  0.     -0.2089  0.      0.0951  0.     -0.0112  0.    ]]

















In [ ]:

sadasd





In [42]:



    


P1 = eRF(P_z, 90)
P1



    


















    
Out[42]:








array([[ 1.],
       [-1.],
       [ 0.]])

















In [43]:



    


P2 = eTE(P_xy, 120)
P2



    


















    
Out[43]:








array([[ 0.7134,  0.0021,  0.2378],
       [ 0.7134,  0.    ,  0.    ],
       [ 0.0013, -0.4213,  0.    ]])

















In [44]:



    


P3 = eTE(P2, 120)
P3



    


















    
Out[44]:








array([[ 0.8639,  0.0032, -0.1853,  0.0005,  0.0566],
       [ 0.8639,  0.0015,  0.1697,  0.    ,  0.    ],
       [ 0.0006, -0.091 , -0.0009, -0.1002,  0.    ]])

















In [45]:



    


P4 = eRelax(P3, TE/2)
P4



    


















    
Out[45]:








array([[ 0.8426,  0.0031, -0.1807,  0.0005,  0.0552],
       [ 0.8426,  0.0015,  0.1655,  0.    ,  0.    ],
       [ 0.0031, -0.0907, -0.0009, -0.0999,  0.    ]])

















In [46]:



    


P5 = eGrad(P4)
P5



    


















    
Out[46]:








array([[ 0.0015,  0.8426,  0.0031, -0.1807,  0.0005,  0.0552],
       [ 0.0015,  0.1655,  0.    ,  0.    ,  0.    ,  0.    ],
       [ 0.0031, -0.0907, -0.0009, -0.0999,  0.    ,  0.    ]])

















In [47]:



    


P6 = eRF(P5, 120)
P6



    


















    
Out[47]:








array([[ 0.0042,  0.2562,  0.    , -0.1317,  0.0001,  0.0138],
       [-0.0012,  0.7519,  0.0031, -0.049 ,  0.0004,  0.0414],
       [-0.0016, -0.2478, -0.0009,  0.1282, -0.0002, -0.0239]])

















In [48]:



    


P7 = eGrad(P6)
P7



    


















    
Out[48]:








array([[ 0.7519,  0.0042,  0.2562,  0.    , -0.1317,  0.0001,  0.0138],
       [ 0.7519,  0.0031, -0.049 ,  0.0004,  0.0414,  0.    ,  0.    ],
       [-0.0016, -0.2478, -0.0009,  0.1282, -0.0002, -0.0239,  0.    ]])

















In [49]:



    


P8 = eRelax(P7, TE/2)
P8



    


















    
Out[49]:








array([[ 0.7333,  0.0041,  0.2498,  0.    , -0.1285,  0.0001,  0.0134],
       [ 0.7333,  0.003 , -0.0478,  0.0004,  0.0403,  0.    ,  0.    ],
       [ 0.0009, -0.2472, -0.0009,  0.1279, -0.0002, -0.0238,  0.    ]])

















In [50]:



    


P9 = eTE(P8, 100)
P9



    


















    
Out[50]:








array([[ 0.6274,  0.0062,  0.0247,  0.0009,  0.2432, -0.0002, -0.0733,  0.    ,  0.0053],
       [ 0.6274,  0.0033,  0.0328,  0.0002, -0.0489,  0.0001,  0.0075,  0.    ,  0.    ],
       [ 0.0019, -0.3315, -0.0016, -0.1225,  0.    ,  0.0657, -0.0001, -0.0064,  0.    ]])

















In [51]:



    


P10 = eTE(P9, 100)
P10



    


















    
Out[51]:








array([[ 0.6807,  0.0073, -0.0527,  0.001 , -0.1349,  0.0004,  0.1627, -0.0001, -0.035 ,  0.    ,  0.0021],
       [ 0.6807,  0.0051,  0.1119,  0.0005,  0.0758, -0.0001, -0.0348,  0.    ,  0.003 ,  0.    ,  0.    ],
       [ 0.0017, -0.2276, -0.0026, -0.0141, -0.0004, -0.1243,  0.0001,  0.0362, -0.    , -0.0025,  0.    ]])

















In [52]:



    


P11 = eTE(P10, 100)
P11



    


















    
Out[52]:








array([[ 0.642 ,  0.0089,  0.1119,  0.0007,  0.0081, -0.0001, -0.1915,  0.0003,  0.1003, -0.0001, -0.0162,  0.    ,  0.0008],
       [ 0.642 ,  0.0068,  0.0139,  0.0009,  0.0301,  0.0002,  0.0573, -0.0001, -0.0171,  0.    ,  0.0012,  0.    ,  0.    ],
       [ 0.0018, -0.2332, -0.0028,  0.064 , -0.0004,  0.0695, -0.0002, -0.0828,  0.0001,  0.0172, -0.    , -0.001 ,  0.    ]])

















In [131]:



    


np.exp(-TE/2 / T1)



    


















    
Out[131]:








0.99750312239746008

















In [150]:



    


T = 250



    











In [151]:



    


M1 = 20
M2 = 100
M3 = 100
M4 = T - M1 - M2 - M3
a0 = np.r_[np.linspace(120, 50, M1), np.linspace(50, 80, M2), np.linspace(80, 90, M3), np.linspace(90, 120, M4)]
plt.plot(a0)
plt.show()



    











In [152]:



    


a1 = a0.copy()
a2 = np.array([120] * T)
a3 = np.array([180] * T)



    











In [153]:



    


T1 = 1
T21 = 100e-3
T22 = 40e-3



    











In [154]:



    


x11 = epg.FSE_signal(np.pi/180*a1, TE, T1, T21).ravel()
x12 = epg.FSE_signal(np.pi/180*a2, TE, T1, T21).ravel()
x13 = epg.FSE_signal(np.pi/180*a3, TE, T1, T21).ravel()

x21 = epg.FSE_signal(np.pi/180*a1, TE, T1, T22).ravel()
x22 = epg.FSE_signal(np.pi/180*a2, TE, T1, T22).ravel()
x23 = epg.FSE_signal(np.pi/180*a3, TE, T1, T22).ravel()



    











In [ ]:



    






    











In [155]:



    


plt.figure()
plt.plot(x11)
plt.plot(x21)
plt.plot(x12)
plt.plot(x22)
plt.plot(x13)
plt.plot(x23)



    


















    
Out[155]:








[<matplotlib.lines.Line2D at 0x1157ce810>]

















In [156]:



    


plt.show()



    











In [157]:



    


print np.linalg.norm(x11)
print np.linalg.norm(x12)
print np.linalg.norm(x13)
print np.linalg.norm(x21)
print np.linalg.norm(x22)
print np.linalg.norm(x23)



    


















    






2.73910716239
2.83252417209
3.08355832516
1.71779149367
1.72557076802
1.87638260069

















In [158]:



    


print np.dot(x11, x21)
print np.dot(x12, x22)
print np.dot(x13, x23)



    


















    






4.3442598635
4.43013589734
5.22886160947

















In [159]:



    


print np.linalg.norm(x11 - x21)**2
print np.linalg.norm(x12 - x22)**2
print np.linalg.norm(x13 - x23)**2



    


















    






1.76499593578
2.14051586621
2.5714203899

















In [160]:



    


def normalized_loss(theta1, theta2, angles_rad, TE, TR):
    T = len(angles_rad)
    x1 = epg.FSE_signal(angles_rad, TE, theta1['T1'], theta1['T2']) * (1 - np.exp(-(TR - T * TE)/theta1['T1']))
    x2 = epg.FSE_signal(angles_rad, TE, theta2['T1'], theta2['T2']) * (1 - np.exp(-(TR - T * TE)/theta2['T1']))

    x1 = x1 / np.linalg.norm(x1, ord=2)
    x2 = x2 / np.linalg.norm(x2, ord=2)

    return - np.dot(x1.ravel(), x2.ravel())



    











In [161]:



    


TR = 1.4
theta1 = {'T1': T1, 'T2': T21}
theta2 = {'T1' : T1, 'T2': T22}
print normalized_loss(theta1, theta2, a1, TE, TR)
print normalized_loss(theta1, theta2, a2, TE, TR)
print normalized_loss(theta1, theta2, a3, TE, TR)



    


















    






-0.938086863237
-0.958632990652
-0.90561248144

















In [ ]:

Content source: jtamir/mri-sim-py

Similar notebooks:

notebook.community | gallery | about