Commit 6cf3e2af authored by Mark Chiew's avatar Mark Chiew
Browse files

Added some matlab vs python comparisons, bloch simulator and partial

fourier
parent a550709d
%% 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
%% Cell type:code id: tags:
``` python
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])
```
%% Cell type:markdown id: tags:
Integrate ODE
%% Cell type:code id: tags:
``` python
t = np.array([0])
M = np.array([[0, 0, 1]])
dt= 0.005
r = ode(bloch_ode)
r.set_integrator('dopri5')
r.set_initial_value(M[0],t[0])
r.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)
```
%% Cell type:markdown id: tags:
Plot Results
%% Cell type:code id: tags:
``` python
_, 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:
``` python
```
%% 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
end
function b = B_eff(t)
if t < 0.25 % No B-field
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
b = [0, 0, 0];
elseif t < 3.00 % Gradient in y-direction
b = [0, 0, 2*pi];
else % No B-field
b = [0, 0, 0];
end
end
\ No newline at end of file
#%% 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
%% 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
%% Cell type:code id: tags:
``` python
h = h5py.File('data.mat','r')
img = np.transpose(h.get('img')['real']+1j*h.get('img')['imag'])
```
%% Cell type:markdown id: tags:
6/8 Partial Fourier sampling
%% Cell type:code id: tags:
``` python
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
```
%% Cell type:markdown id: tags:
Estimate phase
%% 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))))
```
%% Cell type:markdown id: tags:
POCS reconstruction
%% Cell type:code id: tags:
``` python
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)
```
%% Cell type:markdown id: tags:
Plot reconstruction
%% Cell type:code id: tags:
``` python
_, 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)
ax[0].set_title('Zero-filled')
ax[1].imshow(est, vmin=0, vmax=1)
ax[1].set_title('POCS recon')
```
%% Imports
%% Load Data
h = matfile('data.mat');
img = h.img;
%% 6/8 Partial Fourier sampling
n = randn(96) + 1j*randn(96);
y = fftshift(fft2(img),1) + 0*n;
%y(73:end,:) = 0;
%% Estimate phase
pad = padarray(hann(48),24);
phs = exp(1j*angle(ifft2(ifftshift(y.*pad,1))));
%% POCS reconstruction
est = zeros(96);
iters = 10;
for i = 1:iters
est = fftshift(fft2(est.*phs),1);
est(1:72,:) = y(1:72,:);
est = max(real(ifft2(ifftshift(est,1)).*conj(phs)),0);
end
%% Plot reconstruction
figure();
subplot(1,2,1);
imshow(abs(ifft2(ifftshift(y,1))),[0 1],'colormap',jet);
title('Zero-Filled');
subplot(1,2,2);
imshow(est,[0 1],'colormap',jet);
title('POCS recon');
\ No newline at end of file
#%% 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
Supports Markdown
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