diff --git a/talks/matlab_vs_python/bloch/bloch.m b/talks/matlab_vs_python/bloch/bloch.m index 64fdc0afc2f25e8e63f8bc506fafb6a00fba14ba..9af943967c9df78e862d99f52cf8e67a642e05e1 100644 --- a/talks/matlab_vs_python/bloch/bloch.m +++ b/talks/matlab_vs_python/bloch/bloch.m @@ -18,7 +18,7 @@ subplot(2,1,2); plot(pulseq(:,3)); ylabel('Gradient'); -% Integrate ODE +%% Integrate ODE T1 = 1500; T2 = 50; t0 = 0; @@ -27,20 +27,20 @@ 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 +%% 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)); +plot(t, M(:,1), 'linewidth', 2); +plot(t, M(:,2), 'linewidth', 2); +plot(t, M(:,3), 'linewidth', 2); % add legend and grid legend({'Mx','My','Mz'}); grid on; -% define the bloch equation +%% define the bloch equation function dM = bloch_ode(t, M, T1, T2) % get effective B-field at time t B = B_eff(t); @@ -51,7 +51,7 @@ function dM = bloch_ode(t, M, T1, T2) M(1)*B(2) - M(2)*B(1) - (M(3)-1)/T1]; end -% define effective B-field +%% define effective B-field function b = B_eff(t) % Do nothing for 0.25 ms if t < 0.25 diff --git a/talks/matlab_vs_python/partial_fourier/partial_fourier.ipynb b/talks/matlab_vs_python/partial_fourier/partial_fourier.ipynb index 87f826f3d7d200f19fb6ea5241a6a349230f1361..b0424ee7e416b38424aa0974ef2e06e88c4be628 100644 --- a/talks/matlab_vs_python/partial_fourier/partial_fourier.ipynb +++ b/talks/matlab_vs_python/partial_fourier/partial_fourier.ipynb @@ -64,7 +64,8 @@ "In this section:\n", "- np.random.randn\n", "- np.fft\n", - "- 0-based indexing" + "- 0-based indexing\n", + "- image plotting" ] }, { @@ -80,7 +81,11 @@ "y = np.fft.fftshift(np.fft.fft2(img), axes=0) + n\n", "\n", "# set bottom 24/96 lines to 0\n", - "y[72:,:] = 0" + "y[72:,:] = 0\n", + "\n", + "# show sampling\n", + "_, ax = plt.subplots()\n", + "ax.imshow(np.log(np.abs(np.fft.fftshift(y, axes=1))))" ] }, { @@ -96,7 +101,8 @@ "In this section:\n", "- np.pad\n", "- np.hanning\n", - "- reshaping 1D array to 2D array using np.newaxis (or None)" + "- reshaping 1D array to 2D array using np.newaxis (or None)\n", + "- subplots with titles" ] }, { @@ -117,7 +123,14 @@ "low = np.fft.ifft2(np.fft.ifftshift(y*filt, axes=0))\n", "\n", "# get phase image\n", - "phs = np.exp(1j*np.angle(low))" + "phs = np.exp(1j*np.angle(low))\n", + "\n", + "# show phase estimate alongside true phase\n", + "_, ax = plt.subplots(1,2)\n", + "ax[0].imshow(np.angle(img))\n", + "ax[0].set_title('True image phase')\n", + "ax[1].imshow(np.angle(phs))\n", + "ax[1].set_title('Estimated phase')" ] }, { @@ -190,7 +203,7 @@ "\n", "In this section:\n", "- print formatted strings to standard output\n", - "- plotting, with min/max scales\n", + "- 2D subplots with min/max scales, figure size\n", "- np.sum sums over all elements by default" ] }, @@ -212,15 +225,17 @@ "print(f'RMSE for POCS recon: {err_pocs}')\n", "\n", "# plot both recons side-by-side\n", - "_, ax = plt.subplots(1,2,figsize=(16,16))\n", + "_, ax = plt.subplots(2,2,figsize=(16,16))\n", "\n", "# plot zero-filled\n", - "ax[0].imshow(np.abs(zf), vmin=0, vmax=1)\n", - "ax[0].set_title('Zero-filled')\n", + "ax[0,0].imshow(np.abs(zf), vmin=0, vmax=1)\n", + "ax[0,0].set_title('Zero-filled')\n", + "ax[1,0].plot(np.abs(zf[:,47]))\n", "\n", "# plot POCS\n", - "ax[1].imshow(est, vmin=0, vmax=1)\n", - "ax[1].set_title('POCS recon')" + "ax[0,1].imshow(est, vmin=0, vmax=1)\n", + "ax[0,1].set_title('POCS recon')\n", + "ax[1,1].plot(np.abs(est[:,47]))" ] } ], diff --git a/talks/matlab_vs_python/partial_fourier/partial_fourier.m b/talks/matlab_vs_python/partial_fourier/partial_fourier.m index 81e1cc01ceb625f36d565bae12d5e72c3d6a3907..d8b2d084951b43c9a08cf54cbfec9a6a628811b3 100644 --- a/talks/matlab_vs_python/partial_fourier/partial_fourier.m +++ b/talks/matlab_vs_python/partial_fourier/partial_fourier.m @@ -15,6 +15,10 @@ y = fftshift(fft2(img),1) + n; % set bottom 24/96 lines to 0 y(73:end,:) = 0; +% show sampling +figure(); +imshow(log(abs(fftshift(y,2))), [], 'colormap', jet) + %% Estimate phase % create zero-padded hanning filter for ky-filtering filt = padarray(hann(48),24); @@ -25,6 +29,16 @@ low = ifft2(ifftshift(y.*filt,1)); % get phase image phs = exp(1j*angle(low)); +% show phase estimate alongside true phase +figure(); +subplot(1,2,1); +imshow(angle(img), [-pi,pi], 'colormap', hsv); +title('True image phase'); + +subplot(1,2,2); +imshow(angle(phs), [-pi,pi], 'colormap', hsv) +title('Estimated phase'); + %% POCS reconstruction % initialise image estimate to be zeros est = zeros(96); @@ -74,11 +88,15 @@ fprintf(1, 'RMSE for POCS recon: %f\n', err_pocs); figure(); % plot zero-filled -subplot(1,2,1); +subplot(2,2,1); imshow(abs(zf), [0 1]); title('Zero-Filled'); +subplot(2,2,3); +plot(abs(zf(:,48)), 'linewidth', 2); % plot POCS -subplot(1,2,2); +subplot(2,2,2); imshow(est, [0 1]); title('POCS recon'); +subplot(2,2,4); +plot(abs(est(:,48)), 'linewidth', 2);