Merge branch 'master' into 'master'

Added matlab vs python comparisons, bloch simulator and partial fourier

See merge request pytreat-practicals-2020!2

talks/matlab_vs_python/bloch/bloch.ipynb

+{
+ "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
+}

talks/matlab_vs_python/bloch/bloch.m

+%% 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

talks/matlab_vs_python/bloch/bloch.py

+#%% 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

+{
+ "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
+}

talks/matlab_vs_python/partial_fourier/partial_fourier.m

+%% 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

talks/matlab_vs_python/partial_fourier/partial_fourier.py

+#%% 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