diff --git a/talks/matlab_vs_python/bloch/bloch.ipynb b/talks/matlab_vs_python/bloch/bloch.ipynb new file mode 100644 index 0000000000000000000000000000000000000000..a96717a06ac7f261330974fdc39c8e68b2348208 --- /dev/null +++ b/talks/matlab_vs_python/bloch/bloch.ipynb @@ -0,0 +1,127 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Imports" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "from scipy.integrate import ode\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Define bloch and B_eff functions" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def bloch_ode(t,M,T1,T2):\n", + " B = B_eff(t)\n", + " return np.array([M[1]*B[2] - M[2]*B[1] - M[0]/T2,\n", + " M[2]*B[0] - M[0]*B[2] - M[1]/T2,\n", + " M[0]*B[1] - M[1]*B[0] - (M[2]-1)/T1])\n", + "\n", + "def B_eff(t):\n", + " if t < 0.25:\n", + " return np.array([0, 0, 0])\n", + " elif t < 1.25:\n", + " return np.array([1.8*np.sinc(t-0.75), 0, 0])\n", + " elif t < 1.50:\n", + " return np.array([0, 0, 0])\n", + " elif t < 3.00:\n", + " return np.array([0, 0, 2*np.pi])\n", + " else:\n", + " return np.array([0, 0, 0])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Integrate ODE" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "t = np.array([0])\n", + "M = np.array([[0, 0, 1]])\n", + "dt= 0.005\n", + "r = ode(bloch_ode)\n", + "r.set_integrator('dopri5')\n", + "r.set_initial_value(M[0],t[0])\n", + "r.set_f_params(1500, 50)\n", + "while r.successful() and r.t < 5:\n", + " t = np.append(t,r.t+dt)\n", + " M = np.append(M, np.array([r.integrate(t[-1])]),axis=0)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Plot Results" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "_, ax = plt.subplots(figsize=(12,12))\n", + "ax.plot(t,M[:,0], label='Mx')\n", + "ax.plot(t,M[:,1], label='My')\n", + "ax.plot(t,M[:,2], label='Mz')\n", + "ax.legend()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.7.3" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/talks/matlab_vs_python/bloch/bloch.m b/talks/matlab_vs_python/bloch/bloch.m new file mode 100644 index 0000000000000000000000000000000000000000..324fcc85ebf9ff8bd6ea0631f50fca34b6c761e2 --- /dev/null +++ b/talks/matlab_vs_python/bloch/bloch.m @@ -0,0 +1,32 @@ +%% 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 diff --git a/talks/matlab_vs_python/bloch/bloch.py b/talks/matlab_vs_python/bloch/bloch.py new file mode 100644 index 0000000000000000000000000000000000000000..873798d44709e21eb70e3294ede8aee10b118705 --- /dev/null +++ b/talks/matlab_vs_python/bloch/bloch.py @@ -0,0 +1,40 @@ +#%% 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 diff --git a/talks/matlab_vs_python/partial_fourier/data.mat b/talks/matlab_vs_python/partial_fourier/data.mat new file mode 100644 index 0000000000000000000000000000000000000000..c63c984d088a3280bafb2a681a3a8a50853531cf Binary files /dev/null and b/talks/matlab_vs_python/partial_fourier/data.mat differ diff --git a/talks/matlab_vs_python/partial_fourier/partial_fourier.ipynb b/talks/matlab_vs_python/partial_fourier/partial_fourier.ipynb new file mode 100644 index 0000000000000000000000000000000000000000..ea2bb3e782730d4a5ad488077d4ba7a330f35e65 --- /dev/null +++ b/talks/matlab_vs_python/partial_fourier/partial_fourier.ipynb @@ -0,0 +1,136 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Imports" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import h5py\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Load Data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "h = h5py.File('data.mat','r')\n", + "img = np.transpose(h.get('img')['real']+1j*h.get('img')['imag'])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "6/8 Partial Fourier sampling" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "n = np.random.randn(96,96) + 1j*np.random.randn(96,96)\n", + "y = np.fft.fftshift(np.fft.fft2(img),axes=0) + n\n", + "y[72:,:] = 0" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Estimate phase" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "pad = np.pad(np.hanning(48),24,'constant')[:,None]\n", + "phs = np.exp(1j*np.angle(np.fft.ifft2(np.fft.ifftshift(y*pad,axes=0))))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "POCS reconstruction" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "est = np.zeros((96,96))\n", + "iters = 10\n", + "for i in range(iters):\n", + " est = np.fft.fftshift(np.fft.fft2(est*phs),axes=0)\n", + " est[:72,:] = y[:72,:]\n", + " est = np.maximum(np.real(np.fft.ifft2(np.fft.ifftshift(est,axes=0))*np.conj(phs)),0)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Plot reconstruction" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "_, ax = plt.subplots(1,2,figsize=(16,16))\n", + "ax[0].imshow(np.abs(np.fft.ifft2(np.fft.ifftshift(y,axes=0))), vmin=0, vmax=1)\n", + "ax[0].set_title('Zero-filled')\n", + "ax[1].imshow(est, vmin=0, vmax=1)\n", + "ax[1].set_title('POCS recon')" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.7.3" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/talks/matlab_vs_python/partial_fourier/partial_fourier.m b/talks/matlab_vs_python/partial_fourier/partial_fourier.m new file mode 100644 index 0000000000000000000000000000000000000000..b364ab7a192beafbeac83d9dba59ac0fbd655425 --- /dev/null +++ b/talks/matlab_vs_python/partial_fourier/partial_fourier.m @@ -0,0 +1,32 @@ +%% 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 diff --git a/talks/matlab_vs_python/partial_fourier/partial_fourier.py b/talks/matlab_vs_python/partial_fourier/partial_fourier.py new file mode 100644 index 0000000000000000000000000000000000000000..2ceb7ece4f5a288fba11ae2345cb225cc2e1642b --- /dev/null +++ b/talks/matlab_vs_python/partial_fourier/partial_fourier.py @@ -0,0 +1,32 @@ +#%% 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