Skip to content
Snippets Groups Projects
Commit 1c0fd89d authored by Paul McCarthy's avatar Paul McCarthy :mountain_bicyclist:
Browse files

Merge branch 'master' into 'master' 

Fixed weird bug w.r.t. ipynbs created in VS Code, switched ode to odeint

See merge request pytreat-practicals-2020!26
parents d74a8c2c 6e5f1b9a
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
talks/matlab_vs_python/bloch/bloch.ipynb
+ 20
− 43
View file @ 1c0fd89d
%% Cell type:markdown id: tags:
Imports
%% Cell type:code id: tags:
``` python
import numpy as np
from scipy.integrate import ode
from scipy.integrate import ode, solve_ivp
import matplotlib.pyplot as plt
```
%% Cell type:markdown id: tags:
# Define the Bloch equation
$$\frac{d\vec{M}}{dt} = \vec{M}\times \vec{B} - \frac{M_x + M_y}{T2} - \frac{M_z - M_0}{T1}$$
In this section:
- define a function
- numpy functions like np.cross
%% Cell type:code id: tags:
``` python
# define bloch equation
def bloch_ode(t, M, T1, T2):
def bloch(t, M, T1, T2):
# get effective B-field at time t
B = B_eff(t)
# cross product of M and B, add T1 and T2 relaxation terms
return np.array([M[1]*B[2] - M[2]*B[1] - M[0]/T2,
M[2]*B[0] - M[0]*B[2] - M[1]/T2,
M[0]*B[1] - M[1]*B[0] - (M[2]-1)/T1])
# alternatively
#return np.cross(M,B) - np.array([M[0]/T2, M[1]/T2, (M[2]-1)/T1])
```
%% Cell type:markdown id: tags:
# Define the pulse sequence
We work in the rotating frame, so we only need the amplitude envelope of the RF pulse
Typically, B1 excitation fields point in the x- and/or y-directions
Gradient fields point in the z-direction
In this simple example, a simple sinc-pulse excitation pulse is applied for 1 ms along the x-axis
then a gradient is turned on for 1.5 ms after that.
In this section:
- constants such as np.pi
- functions like np.sinc
%% Cell type:code id: tags:
``` python
# define effective B-field
def B_eff(t):
# Do nothing for 0.25 ms
if t < 0.25:
return np.array([0, 0, 0])
# Sinc RF along x-axis and slice-select gradient on for 1.00 ms
elif t < 1.25:
return np.array([np.pi*np.sinc((t-0.75)*4), 0, np.pi])
# Do nothing for 0.25 ms
elif t < 1.50:
return np.array([0, 0, 0])
# Slice refocusing gradient on for 1.50 ms
# Half the area of the slice-select gradient lobe
elif t < 3.00:
return np.array([0, 0, -(1/3)*np.pi])
# Pulse sequence finished
else:
return np.array([0, 0, 0])
```
%% Cell type:markdown id: tags:
# Plot the pulse sequence
In this section:
- unpacking return values
- unwanted return values
- list comprehension
- basic plotting
%% Cell type:code id: tags:
``` python
# Create 2 vertical subplots
_, ax = plt.subplots(2, 1, figsize=(12,12))
# Get pulse sequence B-fields from 0 - 5 ms
pulseq = [B_eff(t) for t in np.linspace(0,5,1000)]
pulseq = np.array(pulseq)
# Plot RF
ax[0].plot(pulseq[:,0])
ax[0].set_ylabel('B1')
# Plot gradient
ax[1].plot(pulseq[:,2])
ax[1].set_ylabel('Gradient')
```
%% Cell type:markdown id: tags:
# Integrate ODE
This uses a Runge-Kutta variant called the "Dormand-Prince method"
This uses a Runge-Kutta variant called the "Dormand-Prince method" by default. Other ode integration methods are available.
In this section:
- list of arrays
- ode solvers
- list appending
- lambdas (anonymous functions)
%% Cell type:code id: tags:
``` python
# Set the initial conditions
# time (t) = 0
# equilibrium magnetization (M) = (0, 0, 1)
t = [0]
M = [np.array([0, 0, 1])]
M = [0, 0, 1]
# Set integrator time-step
dt= 0.005
# Set time interval for integration
t = [0, 5]
# Set up ODE integrator object
r = ode(bloch_ode)
# Choose the integrator method
r.set_integrator('dopri5')
# Pass in initial values
r.set_initial_value(M[0], t[0])
# Set max step size
dt = 0.005
# Set T1 and T2
T1, T2 = 1500, 50
r.set_f_params(T1, T2)
# Integrate by looping over time, moving dt by step size each iteration
# Append new time point and Magnetisation vector at every step to t and M
while r.successful() and r.t < 5:
t.append(r.t + dt)
M.append(r.integrate(t[-1]))
# Integrate ODE
# In Scipy 1.2.0, the first argument to solve_ivp must be a function that takes exactly 2 arguments
sol = solve_ivp(lambda t, M : bloch(t, M, T1, T2), t, M, max_step=dt)
# Convert M to 2-D numpy array from list of arrays
M = np.array(M)
# Grab output
t = sol.t
M = sol.y
```
%% Cell type:markdown id: tags:
# Plot Results
In this section:
- more plotting
%% Cell type:code id: tags:
``` python
# Create single axis
_, ax = plt.subplots(figsize=(12,12))
# Plot x, y and z components of Magnetisation
ax.plot(t, M[:,0], label='Mx')
ax.plot(t, M[:,1], label='My')
ax.plot(t, M[:,2], label='Mz')
ax.plot(t, M[0,:], label='Mx')
ax.plot(t, M[1,:], label='My')
ax.plot(t, M[2,:], label='Mz')
# Add legend and grid
ax.legend()
ax.grid()
```
......
talks/matlab_vs_python/partial_fourier/partial_fourier.ipynb
+ 2
− 14
View file @ 1c0fd89d
......@@ -2,9 +2,7 @@
"cells": [
{
"cell_type": "markdown",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"Imports"
]
......@@ -22,9 +20,7 @@
},
{
"cell_type": "markdown",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Load data\n",
"\n",
......@@ -57,9 +53,7 @@
},
{
"cell_type": "markdown",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# 6/8 Partial Fourier sampling\n",
"\n",
......@@ -91,9 +85,7 @@
},
{
"cell_type": "markdown",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Estimate phase\n",
"\n",
......@@ -130,9 +122,7 @@
},
{
"cell_type": "markdown",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# POCS reconstruction\n",
"\n",
......@@ -192,9 +182,7 @@
},
{
"cell_type": "markdown",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Display error and plot reconstruction\n",
"\n",
......@@ -252,9 +240,9 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.3-final"
"version": "3.7.3"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
\ No newline at end of file
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment