Skip to content
Snippets Groups Projects
Commit 6e5f1b9a authored by Mark Chiew's avatar Mark Chiew
Browse files

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

parent d74a8c2c
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 @ 6e5f1b9a
%% 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 @ 6e5f1b9a
......@@ -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