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

Updated bloch/partial_fourier, more comments

parent 363f308a
No related branches found
No related tags found
1 merge request!13Updated bloch/partial_fourier, more comments
Hide whitespace changes
Inline Side-by-side
Showing
with 412 additions and 147 deletions
talks/matlab_vs_python/bloch/bloch.ipynb
+ 144
− 24
View file @ 54c1373b
%% Cell type:markdown id: tags:
Imports
%% Cell type:code id: tags:
``` python
import numpy as np
from scipy.integrate import ode
import matplotlib.pyplot as plt
```
%% Cell type:markdown id: tags:
Define bloch and B_eff functions
# 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
def bloch_ode(t,M,T1,T2):
# define bloch equation
def bloch_ode(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([1.8*np.sinc(t-0.75), 0, 0])
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, 2*np.pi])
return np.array([0, 0, -(1/3)*np.pi])
# Pulse sequence finished
else:
return np.array([0, 0, 0])
```
%% Cell type:markdown id: tags:
Integrate ODE
# 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"
In this section:
- list of arrays
- ode solvers
- list appending
%% Cell type:code id: tags:
``` python
t = np.array([0])
M = np.array([[0, 0, 1]])
# Set the initial conditions
# time (t) = 0
# equilibrium magnetization (M) = (0, 0, 1)
t = [0]
M = [np.array([0, 0, 1])]
# Set integrator time-step
dt= 0.005
# Set up ODE integrator object
r = ode(bloch_ode)
# Choose the integrator method
r.set_integrator('dopri5')
r.set_initial_value(M[0],t[0])
r.set_f_params(1500, 50)
# Pass in initial values
r.set_initial_value(M[0], t[0])
# 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 = np.append(t,r.t+dt)
M = np.append(M, np.array([r.integrate(t[-1])]),axis=0)
t.append(r.t + dt)
M.append(r.integrate(t[-1]))
# Convert M to 2-D numpy array from list of arrays
M = np.array(M)
```
%% Cell type:markdown id: tags:
Plot Results
# Plot Results
In this section:
- more plotting
%% Cell type:code id: tags:
``` python
# Create single axis
_, ax = plt.subplots(figsize=(12,12))
ax.plot(t,M[:,0], label='Mx')
ax.plot(t,M[:,1], label='My')
ax.plot(t,M[:,2], label='Mz')
ax.legend()
```
%% Cell type:code id: tags:
# 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')
``` python
# Add legend and grid
ax.legend()
ax.grid()
```
......
talks/matlab_vs_python/bloch/bloch.m
+ 66
− 25
View file @ 54c1373b
%% Imports
%% Integrate ODE
[t, M] = ode45(@(t,M)bloch_ode(t,M,1500,50),linspace(0,5,1000),[0;0;1]);
%% Plot Results
clf();hold on;
plot(t,M(:,1));
plot(t,M(:,2));
plot(t,M(:,3));
%% Define bloch and b_eff functions
function dM = bloch_ode(t,M,T1,T2)
B = B_eff(t); % B-effective
dM = [M(2)*B(3) - M(3)*B(2) - M(1)/T2; % dMx/dt
M(3)*B(1) - M(1)*B(3) - M(2)/T2; % dMy/dt
M(1)*B(2) - M(2)*B(1) - (M(3)-1)/T1]; % dMz/dt
% Plot the pulse sequence
% create figure
figure();
% get pulse sequence B-fields from 0-5 ms
pulseq = zeros(1000,3);
for i = 1:1000
pulseq(i,:) = B_eff(i*0.005);
end
% plot RF
subplot(2,1,1);
plot(pulseq(:,1));
ylabel('B1');
% plot gradient
subplot(2,1,2);
plot(pulseq(:,3));
ylabel('Gradient');
% Integrate ODE
T1 = 1500;
T2 = 50;
t0 = 0;
t1 = 5;
dt = 0.005;
M0 = [0; 0; 1];
[t, M] = ode45(@(t,M)bloch_ode(t, M, T1, T2), linspace(t0, t1, (t1-t0)/dt), M0);
% Plot Results
% create figure
figure();hold on;
% plot x, y and z components of Magnetisation
plot(t, M(:,1));
plot(t, M(:,2));
plot(t, M(:,3));
% add legend and grid
legend({'Mx','My','Mz'});
grid on;
% define the bloch equation
function dM = bloch_ode(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
dM = [M(2)*B(3) - M(3)*B(2) - M(1)/T2;
M(3)*B(1) - M(1)*B(3) - M(2)/T2;
M(1)*B(2) - M(2)*B(1) - (M(3)-1)/T1];
end
% define effective B-field
function b = B_eff(t)
if t < 0.25 % No B-field
% Do nothing for 0.25 ms
if t < 0.25
b = [0, 0, 0];
elseif t < 1.25 % 1-ms excitation around x-axis
b = [1.8*sinc(t-0.75), 0, 0];
elseif t < 1.50 % No B-field
% Sinc RF along x-axis and slice-select gradient on for 1.00 ms
elseif t < 1.25
b = [pi*sinc((t-0.75)*4), 0, pi];
% Do nothing for 0.25 ms
elseif t < 1.50
b = [0, 0, 0];
elseif t < 3.00 % Gradient in y-direction
b = [0, 0, 2*pi];
else % No B-field
% Slice refocusing gradient on for 1.50 ms
% Half the area of the slice-select gradient lobe
elseif t < 3.00
b = [0, 0, -(1/3)*pi];
% pulse sequence finished
else
b = [0, 0, 0];
end
end
\ No newline at end of file
end
talks/matlab_vs_python/bloch/bloch.py deleted 100644 → 0
+ 0
− 40
View file @ 363f308a
#%% Imports
import numpy as np
from scipy.integrate import ode
import matplotlib.pyplot as plt
#%% Define bloch and B_eff functions
def bloch_ode(t,M,T1,T2):
B = B_eff(t)
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])
def B_eff(t):
if t < 0.25:
return np.array([0, 0, 0])
elif t < 1.25:
return np.array([1.8*np.sinc(t-0.75), 0, 0])
elif t < 1.50:
return np.array([0, 0, 0])
elif t < 3.00:
return np.array([0, 0, 2*np.pi])
else:
return np.array([0, 0, 0])
#%% Integrate ODE
t = np.array([0])
M = np.array([[0, 0, 1]])
dt= 0.005
r = ode(bloch_ode)\
.set_integrator('dopri5')\
.set_initial_value(M[0],t[0])\
.set_f_params(1500, 50)
while r.successful() and r.t < 5:
t = np.append(t,r.t+dt)
M = np.append(M, np.array([r.integrate(t[-1])]),axis=0)
#%% Plot Results
plt.plot(t,M[:,0])
plt.plot(t,M[:,1])
plt.plot(t,M[:,2])
\ No newline at end of file
talks/matlab_vs_python/partial_fourier/partial_fourier.ipynb
+ 138
− 14
View file @ 54c1373b
%% Cell type:markdown id: tags:
Imports
%% Cell type:code id: tags:
``` python
import h5py
import matplotlib.pyplot as plt
import numpy as np
```
%% Cell type:markdown id: tags:
Load Data
# Load data
Load complex image data from MATLAB mat-file (v7.3 or later), which is actually an HDF5 format
Complex data is loaded as a (real, imag) tuple, so it neds to be explicitly converted to complex double
In this section:
- using h5py module
- np.transpose
- 1j as imaginary constant
%% Cell type:code id: tags:
``` python
# get hdf5 object for the mat-file
h = h5py.File('data.mat','r')
img = np.transpose(h.get('img')['real']+1j*h.get('img')['imag'])
# get img variable from the mat-file
dat = h.get('img')
# turn array of (real, imag) tuples into an array of complex doubles
# transpose to keep data in same orientation as MATLAB
img = np.transpose(dat['real'] + 1j*dat['imag'])
```
%% Cell type:markdown id: tags:
6/8 Partial Fourier sampling
# 6/8 Partial Fourier sampling
Fourier transform the image to get k-space data, and add complex Gaussian noise
To simulate 6/8 Partial Fourier sampling, zero out the bottom 1/4 of k-space
In this section:
- np.random.randn
- np.fft
- 0-based indexing
%% Cell type:code id: tags:
``` python
# generate normally-distributed complex noise
n = np.random.randn(96,96) + 1j*np.random.randn(96,96)
y = np.fft.fftshift(np.fft.fft2(img),axes=0) + n
# Fourier transform the image and add noise
y = np.fft.fftshift(np.fft.fft2(img), axes=0) + n
# set bottom 24/96 lines to 0
y[72:,:] = 0
```
%% Cell type:markdown id: tags:
Estimate phase
# Estimate phase
Filter the k-space data and extract a low-resolution phase estimate
Filtering can help reduce ringing in the phase image
In this section:
- np.pad
- np.hanning
- reshaping 1D array to 2D array using np.newaxis (or None)
%% Cell type:code id: tags:
``` python
pad = np.pad(np.hanning(48),24,'constant')[:,None]
phs = np.exp(1j*np.angle(np.fft.ifft2(np.fft.ifftshift(y*pad,axes=0))))
# create zero-padded hanning filter for ky-filtering
filt = np.pad(np.hanning(48),24,'constant')
# reshape 1D array into 2D array
filt = filt[:,np.newaxis]
# or
# filt = filt[:,None]
# generate low-res image with inverse Fourier transform
low = np.fft.ifft2(np.fft.ifftshift(y*filt, axes=0))
# get phase image
phs = np.exp(1j*np.angle(low))
```
%% Cell type:markdown id: tags:
POCS reconstruction
# POCS reconstruction
Perform the projection-onto-convex-sets (POCS) partial Fourier reconstruction method.
POCS is an iterative scheme estimates the reconstructed image as any element in the intersection of the following two (convex) sets:
1. Set of images consistent with the measured data
2. Set of images that are non-negative real
This requires prior knowledge of the image phase (hence the estimate above), and it works because although we have less than a full k-space of measurements, we're now only estimating half the number of free parameters (real values only, instead of real + imag), and we're no longer under-determined. Equivalently, consider the fact that real-valued images have conjugate symmetric k-spaces, so we only require half of k-space to reconstruct our image.
In this section:
- np.zeros
- range() builtin
- point-wise multiplication (*)
- np.fft operations default to last axis, not first
- np.maximum vs np.max
%% Cell type:code id: tags:
``` python
# initialise image estimate to be zeros
est = np.zeros((96,96))
# set the number of iterations
iters = 10
# each iteration cycles between projections
for i in range(iters):
est = np.fft.fftshift(np.fft.fft2(est*phs),axes=0)
# projection onto data-consistent set:
# use a-priori phase to get complex image
est = est*phs
# Fourier transform to get k-space
est = np.fft.fftshift(np.fft.fft2(est), axes=0)
# replace data with measured lines
est[:72,:] = y[:72,:]
est = np.maximum(np.real(np.fft.ifft2(np.fft.ifftshift(est,axes=0))*np.conj(phs)),0)
# inverse Fourier transform to get back to image space
est = np.fft.ifft2(np.fft.ifftshift(est, axes=0))
# projection onto non-negative reals:
# remove a-priori phase
est = est*np.conj(phs)
# get real part
est = np.real(est)
# ensure output is non-negative
est = np.maximum(est, 0)
```
%% Cell type:markdown id: tags:
Plot reconstruction
# Display error and plot reconstruction
The POCS reconstruction is compared to a zero-filled reconstruction (i.e., where the missing data is zeroed prior to inverse Fourier Transform)
In this section:
- print formatted strings to standard output
- plotting, with min/max scales
- np.sum sums over all elements by default
%% Cell type:code id: tags:
``` python
# compute zero-filled recon
zf = np.fft.ifft2(np.fft.ifftshift(y, axes=0))
# compute rmse for zero-filled and POCS recon
err_zf = np.sqrt(np.sum(np.abs(zf - img)**2))
err_pocs = np.sqrt(np.sum(np.abs(est*phs - img)**2))
# print errors
print(f'RMSE for zero-filled recon: {err_zf}')
print(f'RMSE for POCS recon: {err_pocs}')
# plot both recons side-by-side
_, ax = plt.subplots(1,2,figsize=(16,16))
ax[0].imshow(np.abs(np.fft.ifft2(np.fft.ifftshift(y,axes=0))), vmin=0, vmax=1)
# plot zero-filled
ax[0].imshow(np.abs(zf), vmin=0, vmax=1)
ax[0].set_title('Zero-filled')
# plot POCS
ax[1].imshow(est, vmin=0, vmax=1)
ax[1].set_title('POCS recon')
```
......
talks/matlab_vs_python/partial_fourier/partial_fourier.m
+ 64
− 12
View file @ 54c1373b
%% Imports
%% Load Data
% get matfile object
h = matfile('data.mat');
% get img variable from the mat-file
img = h.img;
%% 6/8 Partial Fourier sampling
% generate normally-distributed complex noise
n = randn(96) + 1j*randn(96);
y = fftshift(fft2(img),1) + 0*n;
%y(73:end,:) = 0;
% Fourier transform the image and add noise
y = fftshift(fft2(img),1) + n;
% set bottom 24/96 lines to 0
y(73:end,:) = 0;
%% Estimate phase
pad = padarray(hann(48),24);
phs = exp(1j*angle(ifft2(ifftshift(y.*pad,1))));
% create zero-padded hanning filter for ky-filtering
filt = padarray(hann(48),24);
% generate low-res image with inverse Fourier transform
low = ifft2(ifftshift(y.*filt,1));
% get phase image
phs = exp(1j*angle(low));
%% POCS reconstruction
% initialise image estimate to be zeros
est = zeros(96);
% set number of iterations
iters = 10;
% each iteration cycles between projections
for i = 1:iters
est = fftshift(fft2(est.*phs),1);
% projection onto data-consistent set:
% use a-priori phase to get complex image
est = est.*phs;
% Fourier transform to get k-space
est = fftshift(fft2(est), 1);
% replace data with measured lines
est(1:72,:) = y(1:72,:);
est = max(real(ifft2(ifftshift(est,1)).*conj(phs)),0);
% inverse Fourier transform to get back to image space
est = ifft2(ifftshift(est, 1));
% projection onto non-negative reals
% remove a-priori phase
est = est.*conj(phs);
% get real part
est = real(est);
% ensure output is non-negative
est = max(est, 0);
end
%% Plot reconstruction
%% Display error and plot reconstruction
% compute zero-filled recon
zf = ifft2(ifftshift(y, 1));
% compute rmse for zero-filled and POCS recon
err_zf = norm(zf(:) - img(:));
err_pocs = norm(est(:).*phs(:) - img(:));
% print errors
fprintf(1, 'RMSE for zero-filled recon: %f\n', err_zf);
fprintf(1, 'RMSE for POCS recon: %f\n', err_pocs);
% plot both recons side-by-side
figure();
% plot zero-filled
subplot(1,2,1);
imshow(abs(ifft2(ifftshift(y,1))),[0 1],'colormap',jet);
imshow(abs(zf), [0 1]);
title('Zero-Filled');
% plot POCS
subplot(1,2,2);
imshow(est,[0 1],'colormap',jet);
title('POCS recon');
\ No newline at end of file
imshow(est, [0 1]);
title('POCS recon');
talks/matlab_vs_python/partial_fourier/partial_fourier.py deleted 100644 → 0
+ 0
− 32
View file @ 363f308a
#%% Imports
import h5py
import matplotlib.pyplot as plt
import numpy as np
#%% Load Data
h = h5py.File('data.mat','r')
img = np.transpose(h.get('img')['real']+1j*h.get('img')['imag'])
#%% 6/8 Partial Fourier sampling
n = np.random.randn(96,96) + 1j*np.random.randn(96,96)
y = np.fft.fftshift(np.fft.fft2(img),axes=0) + n
y[72:,:] = 0
#%% Estimate phase
pad = np.pad(np.hanning(48),24,'constant')[:,None]
phs = np.exp(1j*np.angle(np.fft.ifft2(np.fft.ifftshift(y*pad,axes=0))))
#%% POCS reconstruction
est = np.zeros((96,96))
iters = 10
for i in range(iters):
est = np.fft.fftshift(np.fft.fft2(est*phs),axes=0)
est[:72,:] = y[:72,:]
est = np.maximum(np.real(np.fft.ifft2(np.fft.ifftshift(est,axes=0))*np.conj(phs)),0)
#%% Plot reconstruction
_, ax = plt.subplots(1,2)
ax[0].imshow(np.abs(np.fft.ifft2(np.fft.ifftshift(y,axes=0))),vmin=0,vmax=1)
ax[0].set_title('Zero-filled')
ax[1].imshow(est, vmin=0, vmax=1)
ax[1].set_title('POCS recon')
\ 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