FFT + Time Reversal Reconstruction (#475)

examples/README.md

@@ -27,6 +27,11 @@ Every example has a short readme.md file which briefly describes the purpose of
 - [Focused bowl transducer (axisymmetric)](at_focused_bowl_AS/) ([original example](http://www.k-wave.org/documentation/example_at_piston_and_bowl_transducers.php#heading6), [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/waltsims/k-wave-python/blob/master/examples/at_focused_bowl_AS/at_focused_bowl_AS.ipynb))
 
 - [Focused annular array](at_focused_annular_array_3D/) ([original example](http://www.k-wave.org/documentation/example_at_piston_and_bowl_transducers.php#heading7), [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/waltsims/k-wave-python/blob/master/examples/at_focused_annular_array_3D/at_focused_annular_array_3D.ipynb))
+- [2D FFT Reconstruction For A Line Sensor](pr_2D_FFT_line_sensor/) ([original example](http://www.k-wave.org/documentation/example_pr_2D_fft_line_sensor.php), [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/waltsims/k-wave-python/blob/master/examples/pr_2D_FFT_line_sensor/pr_2D_FFT_line_sensor.ipynb))
+- [3D FFT Reconstruction For A Planar Sensor](pr_3D_FFT_planar_sensor/) ([original example](http://www.k-wave.org/documentation/example_pr_3D_fft_planar_sensor.php), [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/waltsims/k-wave-python/blob/master/examples/pr_3D_FFT_planar_sensor/pr_3D_FFT_planar_sensor.ipynb))
+- [2D Time Reversal Reconstruction For A Line Sensor](pr_2D_TR_line_sensor/) ([original example](http://www.k-wave.org/documentation/example_pr_2D_tr_line_sensor.php), [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/waltsims/k-wave-python/blob/master/examples/pr_2D_TR_line_sensor/pr_2D_TR_line_sensor.ipynb))
+- [3D Time Reversal Reconstruction For A Planar Sensor](pr_3D_TR_planar_sensor/) ([original example](http://www.k-wave.org/documentation/example_pr_3D_tr_planar_sensor.php), [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/waltsims/k-wave-python/blob/master/examples/pr_3D_TR_planar_sensor/pr_3D_TR_planar_sensor.ipynb))
+
 
 - [Focussed Detector In 2D Example](sd_focussed_detector_2D/) ([original example](http://www.k-wave.org/documentation/example_sd_focussed_detector_2D.php), [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/waltsims/k-wave-python/blob/master/examples/sd_focussed_detector_2D/sd_focussed_detector_2D.ipynb))
 

examples/pr_2D_FFT_line_sensor/README.md

@@ -0,0 +1,7 @@
+# 2D FFT Reconstruction For A Line Sensor Example 
+
+[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/waltsims/k-wave-python/blob/master/examples/pr_2D_FFT_line_sensor/pr_2D_FFT_line_sensor.ipynb)
+
+This example demonstrates the use of k-Wave for the reconstruction of a two-dimensional photoacoustic wave-field recorded over a linear array of sensor elements. The sensor data is simulated using [kspaceFirstOrder2D](https://k-wave-python.readthedocs.io/en/latest/kwave.kspaceFirstOrder2D.html) and reconstructed using kspaceLineRecon. It builds on the Homogeneous Propagation Medium and Heterogeneous Propagation Medium examples.
+
+To read more, visit the [original example page](http://www.k-wave.org/documentation/example_pr_2D_fft_line_sensor.php) or its [implementation](https://github.com/ucl-bug/k-wave/blob/main/k-Wave/examples/example_pr_2D_FFT_line_sensor.m).

examples/pr_2D_FFT_line_sensor/pr_2D_FFT_line_sensor.ipynb

@@ -0,0 +1,295 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "id": "initial_id",
+   "metadata": {
+    "collapsed": true
+   },
+   "source": [
+    "%%capture\n",
+    "!pip install k-wave-python "
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "metadata": {},
+   "cell_type": "code",
+   "source": [
+    "import matplotlib.pyplot as plt\n",
+    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
+    "import numpy as np\n",
+    "from scipy.interpolate import RegularGridInterpolator\n",
+    "\n",
+    "from kwave.data import Vector\n",
+    "from kwave.kgrid import kWaveGrid\n",
+    "from kwave.kmedium import kWaveMedium\n",
+    "from kwave.ksensor import kSensor\n",
+    "from kwave.ksource import kSource\n",
+    "from kwave.kspaceFirstOrder2D import kspaceFirstOrder2D\n",
+    "from kwave.kspaceLineRecon import kspaceLineRecon\n",
+    "from kwave.options.simulation_execution_options import SimulationExecutionOptions\n",
+    "from kwave.options.simulation_options import SimulationOptions\n",
+    "from kwave.utils.colormap import get_color_map\n",
+    "from kwave.utils.mapgen import make_disc\n",
+    "from kwave.utils.filters import smooth"
+   ],
+   "id": "48677f5b96e1511e",
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "metadata": {},
+   "cell_type": "markdown",
+   "source": [
+    "## 2D FFT Reconstruction For A Line Sensor Example\n",
+    "\n",
+    "This example demonstrates the use of k-Wave for the reconstruction of a\n",
+    "two-dimensional photoacoustic wave-field recorded  over a linear array of\n",
+    "sensor elements  The sensor data is simulated using kspaceFirstOrder2D\n",
+    "and reconstructed using kspaceLineRecon. It builds on the Homogeneous\n",
+    "Propagation Medium and Heterogeneous Propagation Medium examples. "
+   ],
+   "id": "262efa49a13b9574"
+  },
+  {
+   "metadata": {},
+   "cell_type": "markdown",
+   "source": "### SIMULATION",
+   "id": "310bcfe8541609c2"
+  },
+  {
+   "metadata": {},
+   "cell_type": "code",
+   "source": [
+    "# create the computational grid\n",
+    "PML_size = 20  # size of the PML in grid points\n",
+    "N = Vector([128, 256]) - 2 * PML_size  # number of grid points\n",
+    "d = Vector([0.1e-3, 0.1e-3])  # grid point spacing [m]\n",
+    "kgrid = kWaveGrid(N, d)"
+   ],
+   "id": "6d7c5328ae5dd3ee",
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "metadata": {},
+   "cell_type": "code",
+   "source": [
+    "# define the properties of the propagation medium\n",
+    "medium = kWaveMedium(sound_speed=1500)  # [m/s]"
+   ],
+   "id": "bd5eece140453f0f",
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "metadata": {},
+   "cell_type": "code",
+   "source": [
+    "# create initial pressure distribution using makeDisc\n",
+    "disc_magnitude = 5  # [Pa]\n",
+    "disc_pos = Vector([60, 140])\n",
+    "disc_radius = 5\n",
+    "disc_2 = disc_magnitude * make_disc(N, disc_pos, disc_radius)\n",
+    "\n",
+    "disc_pos = Vector([30, 110])\n",
+    "disc_radius = 8\n",
+    "disc_1 = disc_magnitude * make_disc(N, disc_pos, disc_radius)\n",
+    "\n",
+    "# smooth the initial pressure distribution and restore the magnitude\n",
+    "p0 = disc_1 + disc_2\n",
+    "p0 = smooth(p0, restore_max=True)\n",
+    "\n",
+    "source = kSource()\n",
+    "source.p0 = p0"
+   ],
+   "id": "2ce3b8059a9ec319",
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "metadata": {},
+   "cell_type": "code",
+   "source": [
+    "# define a binary line sensor\n",
+    "sensor = kSensor()\n",
+    "sensor.mask = np.zeros(N)\n",
+    "sensor.mask[0] = 1"
+   ],
+   "id": "cf6961e373c6f8e8",
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "metadata": {},
+   "cell_type": "code",
+   "source": [
+    "%%capture\n",
+    "\n",
+    "# create the time array\n",
+    "kgrid.makeTime(medium.sound_speed)"
+   ],
+   "id": "5a278b1e9b253295",
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "metadata": {},
+   "cell_type": "code",
+   "source": [
+    "# set the input arguments: force the PML to be outside the computational grid\n",
+    "simulation_options = SimulationOptions(\n",
+    "    save_to_disk=True,\n",
+    "    pml_inside=False,\n",
+    "    pml_size=PML_size,\n",
+    "    smooth_p0=False,\n",
+    ")\n",
+    "execution_options = SimulationExecutionOptions(is_gpu_simulation=True)"
+   ],
+   "id": "a3802e3e482dc77",
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "metadata": {},
+   "cell_type": "code",
+   "source": [
+    "# run the simulation\n",
+    "sensor_data = kspaceFirstOrder2D(kgrid, source, sensor, medium, simulation_options, execution_options)\n",
+    "sensor_data = sensor_data['p'].T\n",
+    "\n",
+    "# reconstruct the initial pressure\n",
+    "p_xy = kspaceLineRecon(sensor_data.T, dy=d[1], dt=kgrid.dt.item(), c=medium.sound_speed.item(),\n",
+    "                       pos_cond=True, interp='linear')\n",
+    "\n",
+    "# define a second k-space grid using the dimensions of p_xy\n",
+    "N_recon = Vector(p_xy.shape)\n",
+    "d_recon = Vector([kgrid.dt.item() * medium.sound_speed.item(), kgrid.dy])\n",
+    "kgrid_recon = kWaveGrid(N_recon, d_recon)\n",
+    "\n",
+    "# resample p_xy to be the same size as source.p0\n",
+    "interp_func = RegularGridInterpolator(\n",
+    "    (kgrid_recon.x_vec[:, 0] - kgrid_recon.x_vec.min(), kgrid_recon.y_vec[:, 0]),\n",
+    "    p_xy, method='linear'\n",
+    ")\n",
+    "query_points = np.stack((kgrid.x - kgrid.x.min(), kgrid.y), axis=-1)\n",
+    "p_xy_rs = interp_func(query_points)"
+   ],
+   "id": "154c372469cd6063",
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "metadata": {},
+   "cell_type": "markdown",
+   "source": "### VISUALIZATION\n",
+   "id": "b07f65d05feaf902"
+  },
+  {
+   "metadata": {},
+   "cell_type": "code",
+   "source": "cmap = get_color_map()",
+   "id": "30edd081e39aa21b",
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "metadata": {},
+   "cell_type": "code",
+   "source": [
+    "# plot the initial pressure and sensor distribution\n",
+    "fig, ax = plt.subplots(1, 1)\n",
+    "im = ax.imshow(p0 + sensor.mask[PML_size:-PML_size, PML_size:-PML_size] * disc_magnitude,\n",
+    "               extent=[kgrid.y_vec.min() * 1e3, kgrid.y_vec.max() * 1e3, kgrid.x_vec.max() * 1e3, kgrid.x_vec.min() * 1e3],\n",
+    "               vmin=-disc_magnitude, vmax=disc_magnitude, cmap=cmap)\n",
+    "divider = make_axes_locatable(ax)\n",
+    "cax = divider.append_axes(\"right\", size=\"3%\", pad=\"2%\")\n",
+    "ax.set_ylabel('x-position [mm]')\n",
+    "ax.set_xlabel('y-position [mm]')\n",
+    "fig.colorbar(im, cax=cax)\n",
+    "plt.show()"
+   ],
+   "id": "d02e73401dee274c",
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "metadata": {},
+   "cell_type": "code",
+   "source": [
+    "fig, ax = plt.subplots(1, 1)\n",
+    "im = ax.imshow(sensor_data, vmin=-1, vmax=1, cmap=cmap, aspect='auto')\n",
+    "divider = make_axes_locatable(ax)\n",
+    "cax = divider.append_axes(\"right\", size=\"3%\", pad=\"2%\")\n",
+    "ax.set_ylabel('Sensor Position')\n",
+    "ax.set_xlabel('Time Step')\n",
+    "fig.colorbar(im, cax=cax)\n",
+    "plt.show()"
+   ],
+   "id": "fc203d6e2a22dec6",
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "metadata": {},
+   "cell_type": "code",
+   "source": [
+    "# plot the reconstructed initial pressure\n",
+    "fig, ax = plt.subplots(1, 1)\n",
+    "im = ax.imshow(p_xy_rs,\n",
+    "               extent=[kgrid.y_vec.min() * 1e3, kgrid.y_vec.max() * 1e3, kgrid.x_vec.max() * 1e3, kgrid.x_vec.min() * 1e3],\n",
+    "               vmin=-disc_magnitude, vmax=disc_magnitude, cmap=cmap)\n",
+    "divider = make_axes_locatable(ax)\n",
+    "cax = divider.append_axes(\"right\", size=\"3%\", pad=\"2%\")\n",
+    "ax.set_ylabel('x-position [mm]')\n",
+    "ax.set_xlabel('y-position [mm]')\n",
+    "fig.colorbar(im, cax=cax)\n",
+    "plt.show()"
+   ],
+   "id": "2da3a2a311fd1683",
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "metadata": {},
+   "cell_type": "code",
+   "source": [
+    "# plot a profile for comparison\n",
+    "plt.plot(kgrid.y_vec[:, 0] * 1e3, p0[disc_pos[0], :], 'k-', label='Initial Pressure')\n",
+    "plt.plot(kgrid.y_vec[:, 0] * 1e3, p_xy_rs[disc_pos[0], :], 'r--', label='Reconstructed Pressure')\n",
+    "plt.xlabel('y-position [mm]')\n",
+    "plt.ylabel('Pressure')\n",
+    "plt.legend()\n",
+    "plt.axis('tight')\n",
+    "plt.ylim([0, 5.1])\n",
+    "plt.show()"
+   ],
+   "id": "f7ae5b0a2f5e147e",
+   "outputs": [],
+   "execution_count": null
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 2
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython2",
+   "version": "2.7.6"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}