From 4874253d6baef7179e995ee53ac94b23b9752177 Mon Sep 17 00:00:00 2001 From: Fabian Schneider <56137202+faberno@users.noreply.github.com> Date: Tue, 24 Dec 2024 03:34:35 +0100 Subject: [PATCH] FFT + Time Reversal Reconstruction (#475) --- examples/README.md | 5 + examples/pr_2D_FFT_line_sensor/README.md | 7 + .../pr_2D_FFT_line_sensor.ipynb | 295 +++++++++++++++ .../pr_2D_FFT_line_sensor.py | 150 ++++++++ examples/pr_2D_TR_line_sensor/README.md | 7 + .../pr_2D_TR_line_sensor.ipynb | 337 ++++++++++++++++++ .../pr_2D_TR_line_sensor.py | 173 +++++++++ examples/pr_3D_FFT_planar_sensor/README.md | 7 + .../pr_3D_FFT_planar_sensor.ipynb | 310 ++++++++++++++++ .../pr_3D_FFT_planar_sensor.py | 163 +++++++++ examples/pr_3D_TR_planar_sensor/README.md | 7 + .../pr_3D_TR_planar_sensor.ipynb | 302 ++++++++++++++++ .../pr_3D_TR_planar_sensor.py | 176 +++++++++ kwave/executor.py | 62 ++-- kwave/kWaveSimulation.py | 24 +- kwave/kspaceLineRecon.py | 136 +++++++ kwave/kspacePlaneRecon.py | 136 +++++++ tests/test_executor.py | 68 ++++ tests/test_simulation.py | 51 ++- 19 files changed, 2376 insertions(+), 40 deletions(-) create mode 100644 examples/pr_2D_FFT_line_sensor/README.md create mode 100644 examples/pr_2D_FFT_line_sensor/pr_2D_FFT_line_sensor.ipynb create mode 100644 examples/pr_2D_FFT_line_sensor/pr_2D_FFT_line_sensor.py create mode 100644 examples/pr_2D_TR_line_sensor/README.md create mode 100644 examples/pr_2D_TR_line_sensor/pr_2D_TR_line_sensor.ipynb create mode 100644 examples/pr_2D_TR_line_sensor/pr_2D_TR_line_sensor.py create mode 100644 examples/pr_3D_FFT_planar_sensor/README.md create mode 100644 examples/pr_3D_FFT_planar_sensor/pr_3D_FFT_planar_sensor.ipynb create mode 100644 examples/pr_3D_FFT_planar_sensor/pr_3D_FFT_planar_sensor.py create mode 100644 examples/pr_3D_TR_planar_sensor/README.md create mode 100644 examples/pr_3D_TR_planar_sensor/pr_3D_TR_planar_sensor.ipynb create mode 100644 examples/pr_3D_TR_planar_sensor/pr_3D_TR_planar_sensor.py create mode 100644 kwave/kspaceLineRecon.py create mode 100644 kwave/kspacePlaneRecon.py diff --git a/examples/README.md b/examples/README.md index bf5da62e..664a1a91 100644 --- a/examples/README.md +++ b/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)) diff --git a/examples/pr_2D_FFT_line_sensor/README.md b/examples/pr_2D_FFT_line_sensor/README.md new file mode 100644 index 00000000..950cef92 --- /dev/null +++ b/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). \ No newline at end of file diff --git a/examples/pr_2D_FFT_line_sensor/pr_2D_FFT_line_sensor.ipynb b/examples/pr_2D_FFT_line_sensor/pr_2D_FFT_line_sensor.ipynb new file mode 100644 index 00000000..9ddcc80a --- /dev/null +++ b/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 +} diff --git a/examples/pr_2D_FFT_line_sensor/pr_2D_FFT_line_sensor.py b/examples/pr_2D_FFT_line_sensor/pr_2D_FFT_line_sensor.py new file mode 100644 index 00000000..48143df5 --- /dev/null +++ b/examples/pr_2D_FFT_line_sensor/pr_2D_FFT_line_sensor.py @@ -0,0 +1,150 @@ +import matplotlib.pyplot as plt +from mpl_toolkits.axes_grid1 import make_axes_locatable +import numpy as np +from scipy.interpolate import RegularGridInterpolator + +from kwave.data import Vector +from kwave.kgrid import kWaveGrid +from kwave.kmedium import kWaveMedium +from kwave.ksensor import kSensor +from kwave.ksource import kSource +from kwave.kspaceFirstOrder2D import kspaceFirstOrder2D +from kwave.kspaceLineRecon import kspaceLineRecon +from kwave.options.simulation_execution_options import SimulationExecutionOptions +from kwave.options.simulation_options import SimulationOptions +from kwave.utils.colormap import get_color_map +from kwave.utils.mapgen import make_disc +from kwave.utils.filters import smooth + + +# 2D FFT Reconstruction For A Line Sensor Example + +# 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 +# and reconstructed using kspaceLineRecon. It builds on the Homogeneous +# Propagation Medium and Heterogeneous Propagation Medium examples. + + +def main(): + + # -------------------- + # SIMULATION + # -------------------- + + # create the computational grid + PML_size = 20 # size of the PML in grid points + N = Vector([128, 256]) - 2 * PML_size # number of grid points + d = Vector([0.1e-3, 0.1e-3]) # grid point spacing [m] + kgrid = kWaveGrid(N, d) + + # define the properties of the propagation medium + medium = kWaveMedium(sound_speed=1500) # [m/s] + + # create initial pressure distribution using makeDisc + disc_magnitude = 5 # [Pa] + disc_pos = Vector([60, 140]) + disc_radius = 5 + disc_2 = disc_magnitude * make_disc(N, disc_pos, disc_radius) + + disc_pos = Vector([30, 110]) + disc_radius = 8 + disc_1 = disc_magnitude * make_disc(N, disc_pos, disc_radius) + + # smooth the initial pressure distribution and restore the magnitude + p0 = disc_1 + disc_2 + p0 = smooth(p0, restore_max=True) + + source = kSource() + source.p0 = p0 + + # define a binary line sensor + sensor = kSensor() + sensor.mask = np.zeros(N) + sensor.mask[0] = 1 + + # create the time array + kgrid.makeTime(medium.sound_speed) + + # set the input arguments: force the PML to be outside the computational grid + simulation_options = SimulationOptions( + save_to_disk=True, + pml_inside=False, + pml_size=PML_size, + smooth_p0=False, + ) + execution_options = SimulationExecutionOptions(is_gpu_simulation=True) + + # run the simulation + sensor_data = kspaceFirstOrder2D(kgrid, source, sensor, medium, simulation_options, execution_options) + sensor_data = sensor_data['p'].T + + # reconstruct the initial pressure + p_xy = kspaceLineRecon(sensor_data.T, dy=d[1], dt=kgrid.dt.item(), c=medium.sound_speed.item(), + pos_cond=True, interp='linear') + + # define a second k-space grid using the dimensions of p_xy + N_recon = Vector(p_xy.shape) + d_recon = Vector([kgrid.dt.item() * medium.sound_speed.item(), kgrid.dy]) + kgrid_recon = kWaveGrid(N_recon, d_recon) + + # resample p_xy to be the same size as source.p0 + interp_func = RegularGridInterpolator( + (kgrid_recon.x_vec[:, 0] - kgrid_recon.x_vec.min(), kgrid_recon.y_vec[:, 0]), + p_xy, method='linear' + ) + query_points = np.stack((kgrid.x - kgrid.x.min(), kgrid.y), axis=-1) + p_xy_rs = interp_func(query_points) + + + # -------------------- + # VISUALIZATION + # -------------------- + + cmap = get_color_map() + + # plot the initial pressure and sensor distribution + fig, ax = plt.subplots(1, 1) + im = ax.imshow(p0 + sensor.mask[PML_size:-PML_size, PML_size:-PML_size] * disc_magnitude, + extent=[kgrid.y_vec.min() * 1e3, kgrid.y_vec.max() * 1e3, kgrid.x_vec.max() * 1e3, kgrid.x_vec.min() * 1e3], + vmin=-disc_magnitude, vmax=disc_magnitude, cmap=cmap) + divider = make_axes_locatable(ax) + cax = divider.append_axes("right", size="3%", pad="2%") + ax.set_ylabel('x-position [mm]') + ax.set_xlabel('y-position [mm]') + fig.colorbar(im, cax=cax) + plt.show() + + fig, ax = plt.subplots(1, 1) + im = ax.imshow(sensor_data, vmin=-1, vmax=1, cmap=cmap, aspect='auto') + divider = make_axes_locatable(ax) + cax = divider.append_axes("right", size="3%", pad="2%") + ax.set_ylabel('Sensor Position') + ax.set_xlabel('Time Step') + fig.colorbar(im, cax=cax) + plt.show() + + # plot the reconstructed initial pressure + fig, ax = plt.subplots(1, 1) + im = ax.imshow(p_xy_rs, + extent=[kgrid.y_vec.min() * 1e3, kgrid.y_vec.max() * 1e3, kgrid.x_vec.max() * 1e3, kgrid.x_vec.min() * 1e3], + vmin=-disc_magnitude, vmax=disc_magnitude, cmap=cmap) + divider = make_axes_locatable(ax) + cax = divider.append_axes("right", size="3%", pad="2%") + ax.set_ylabel('x-position [mm]') + ax.set_xlabel('y-position [mm]') + fig.colorbar(im, cax=cax) + plt.show() + + # plot a profile for comparison + plt.plot(kgrid.y_vec[:, 0] * 1e3, p0[disc_pos[0], :], 'k-', label='Initial Pressure') + plt.plot(kgrid.y_vec[:, 0] * 1e3, p_xy_rs[disc_pos[0], :], 'r--', label='Reconstructed Pressure') + plt.xlabel('y-position [mm]') + plt.ylabel('Pressure') + plt.legend() + plt.axis('tight') + plt.ylim([0, 5.1]) + plt.show() + +if __name__ == '__main__': + main() \ No newline at end of file diff --git a/examples/pr_2D_TR_line_sensor/README.md b/examples/pr_2D_TR_line_sensor/README.md new file mode 100644 index 00000000..92384a05 --- /dev/null +++ b/examples/pr_2D_TR_line_sensor/README.md @@ -0,0 +1,7 @@ +# 2D Time Reversal 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_TR_line_sensor/pr_2D_TR_line_sensor.ipynb) + +This example demonstrates the use of k-Wave for the time-reversal reconstruction of a two-dimensional photoacoustic wave-field recorded over a linear array of sensor elements. The sensor data is simulated and then time-reversed using [kspaceFirstOrder2D](https://k-wave-python.readthedocs.io/en/latest/kwave.kspaceFirstOrder2D.html). It builds on the [2D FFT Reconstruction For A Line Sensor Example](../pr_2D_TR_line_sensor/). + +To read more, visit the [original example page](http://www.k-wave.org/documentation/example_pr_2D_tr_line_sensor.php) or its [implementation](https://github.com/ucl-bug/k-wave/blob/main/k-Wave/examples/example_pr_2D_TR_line_sensor.m). \ No newline at end of file diff --git a/examples/pr_2D_TR_line_sensor/pr_2D_TR_line_sensor.ipynb b/examples/pr_2D_TR_line_sensor/pr_2D_TR_line_sensor.ipynb new file mode 100644 index 00000000..2d0eef2a --- /dev/null +++ b/examples/pr_2D_TR_line_sensor/pr_2D_TR_line_sensor.ipynb @@ -0,0 +1,337 @@ +{ + "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": "4d18dbe95c4a69b0", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": [ + "## 2D Time Reversal Reconstruction For A Line Sensor Example\n", + "\n", + "This example demonstrates the use of k-Wave for the time-reversal\n", + "reconstruction of a two-dimensional photoacoustic wave-field recorded\n", + "over a linear array of sensor elements. The sensor data is simulated and\n", + "then time-reversed using kspaceFirstOrder2D. It builds on the 2D FFT \n", + "Reconstruction For A Line Sensor Example." + ], + "id": "88cb3b08164514b0" + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": "### SIMULATION", + "id": "76db584c54e13248" + }, + { + "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": "6eeb7e260b43f43", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "cell_type": "code", + "source": [ + "# define the properties of the propagation medium\n", + "medium = kWaveMedium(\n", + " sound_speed=1500,\t# [m/s]\n", + ")" + ], + "id": "7f72a8c783169161", + "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]) # [grid points]\n", + "disc_radius = 5 # [grid points]\n", + "disc_2 = disc_magnitude * make_disc(N, disc_pos, disc_radius)\n", + "\n", + "disc_pos = Vector([30, 110]) # [grid points]\n", + "disc_radius = 8 # [grid points]\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 = smooth(disc_1 + disc_2, restore_max=True)\n", + "\n", + "# assign to the source structure\n", + "source = kSource()\n", + "source.p0 = p0" + ], + "id": "619c6b4d5a08c7a5", + "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": "59bbb876871f155e", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "cell_type": "code", + "source": [ + "%%capture\n", + "\n", + "# create the time array\n", + "kgrid.makeTime(medium.sound_speed)" + ], + "id": "a23b0c47baf1bde2", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "cell_type": "code", + "source": [ + "# set the input arguements: force the PML to be outside the computational\n", + "# grid; switch off p0 smoothing within kspaceFirstOrder2D\n", + "simulation_options = SimulationOptions(\n", + " pml_inside=False,\n", + " pml_size=PML_size,\n", + " smooth_p0=False,\n", + " save_to_disk=True,\n", + " data_cast=\"single\",\n", + ")\n", + "execution_options = SimulationExecutionOptions(is_gpu_simulation=True)" + ], + "id": "c85b2aacc4cc6bcc", + "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" + ], + "id": "db85c9943ed5f91a", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "cell_type": "code", + "source": [ + "# reset the initial pressure and sensor\n", + "source = kSource()\n", + "sensor = kSensor()\n", + "sensor.mask = np.zeros(N)\n", + "sensor.mask[0] = 1\n", + "\n", + "# assign the time reversal data\n", + "sensor.time_reversal_boundary_data = sensor_data\n", + "\n", + "# run the time reversal reconstruction\n", + "p0_recon = kspaceFirstOrder2D(kgrid, source, sensor, medium, simulation_options, execution_options)\n", + "p0_recon = p0_recon['p_final'].T\n", + "\n", + "# add first order compensation for only recording over a half plane\n", + "p0_recon = 2 * p0_recon" + ], + "id": "fc22514bed1e2a98", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "cell_type": "code", + "source": [ + "# repeat the FFT reconstruction for comparison\n", + "p_xy = kspaceLineRecon(sensor_data.T, dy=d[1], dt=kgrid.dt.item(), c=medium.sound_speed.item(), 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 * 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": "431b7bb0ee525b5e", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": "### VISUALIZATION", + "id": "7305f8eba3374df2" + }, + { + "metadata": {}, + "cell_type": "code", + "source": "cmap = get_color_map()", + "id": "7c8e399f00971c60", + "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 * 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": "b159649563e8767e", + "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(p0_recon,\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": "9bc5d6d4b6761d4", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "cell_type": "code", + "source": [ + "# Apply a positivity condition\n", + "p0_recon[p0_recon < 0] = 0\n", + "\n", + "# Plot the reconstructed initial pressure with positivity condition\n", + "fig, ax = plt.subplots(1, 1)\n", + "im = ax.imshow(p0_recon,\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": "4b70ab1746f0ed88", + "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='FFT Reconstruction')\n", + "plt.plot(kgrid.y_vec[:, 0] * 1e3, p0_recon[disc_pos[0], :], 'b:', label='Time Reversal')\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": "a95d4b11c4111e59", + "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 +} diff --git a/examples/pr_2D_TR_line_sensor/pr_2D_TR_line_sensor.py b/examples/pr_2D_TR_line_sensor/pr_2D_TR_line_sensor.py new file mode 100644 index 00000000..44630804 --- /dev/null +++ b/examples/pr_2D_TR_line_sensor/pr_2D_TR_line_sensor.py @@ -0,0 +1,173 @@ +import matplotlib.pyplot as plt +from mpl_toolkits.axes_grid1 import make_axes_locatable +import numpy as np +from scipy.interpolate import RegularGridInterpolator + +from kwave.data import Vector +from kwave.kgrid import kWaveGrid +from kwave.kmedium import kWaveMedium +from kwave.ksensor import kSensor +from kwave.ksource import kSource +from kwave.kspaceFirstOrder2D import kspaceFirstOrder2D +from kwave.kspaceLineRecon import kspaceLineRecon +from kwave.options.simulation_execution_options import SimulationExecutionOptions +from kwave.options.simulation_options import SimulationOptions +from kwave.utils.colormap import get_color_map +from kwave.utils.mapgen import make_disc +from kwave.utils.filters import smooth + + +# 2D Time Reversal Reconstruction For A Line Sensor Example + +# This example demonstrates the use of k-Wave for the time-reversal +# reconstruction of a two-dimensional photoacoustic wave-field recorded +# over a linear array of sensor elements. The sensor data is simulated and +# then time-reversed using kspaceFirstOrder2D. It builds on the 2D FFT +# Reconstruction For A Line Sensor Example. + + +def main(): + + # -------------------- + # SIMULATION + # -------------------- + + # create the computational grid + PML_size = 20 # size of the PML in grid points + N = Vector([128, 256]) - 2 * PML_size # number of grid points + d = Vector([0.1e-3, 0.1e-3]) # grid point spacing [m] + kgrid = kWaveGrid(N, d) + + # define the properties of the propagation medium + medium = kWaveMedium( + sound_speed=1500, # [m/s] + ) + + # create initial pressure distribution using makeDisc + disc_magnitude = 5 # [Pa] + disc_pos = Vector([60, 140]) # [grid points] + disc_radius = 5 # [grid points] + disc_2 = disc_magnitude * make_disc(N, disc_pos, disc_radius) + + disc_pos = Vector([30, 110]) # [grid points] + disc_radius = 8 # [grid points] + disc_1 = disc_magnitude * make_disc(N, disc_pos, disc_radius) + + # smooth the initial pressure distribution and restore the magnitude + p0 = smooth(disc_1 + disc_2, restore_max=True) + + # assign to the source structure + source = kSource() + source.p0 = p0 + + # define a binary line sensor + sensor = kSensor() + sensor.mask = np.zeros(N) + sensor.mask[0] = 1 + + # set the input arguements: force the PML to be outside the computational + # grid; switch off p0 smoothing within kspaceFirstOrder2D + simulation_options = SimulationOptions( + pml_inside=False, + pml_size=PML_size, + smooth_p0=False, + save_to_disk=True, + data_cast="single", + ) + execution_options = SimulationExecutionOptions(is_gpu_simulation=True) + + # run the simulation + sensor_data = kspaceFirstOrder2D(kgrid, source, sensor, medium, simulation_options, execution_options) + sensor_data = sensor_data['p'].T + + # reset the initial pressure and sensor + source = kSource() + sensor = kSensor() + sensor.mask = np.zeros(N) + sensor.mask[0] = 1 + + # assign the time reversal data + sensor.time_reversal_boundary_data = sensor_data + + # run the time reversal reconstruction + p0_recon = kspaceFirstOrder2D(kgrid, source, sensor, medium, simulation_options, execution_options) + p0_recon = p0_recon['p_final'].T + + # add first order compensation for only recording over a half plane + p0_recon = 2 * p0_recon + + # repeat the FFT reconstruction for comparison + p_xy = kspaceLineRecon(sensor_data.T, dy=d[1], dt=kgrid.dt.item(), c=medium.sound_speed.item(), pos_cond=True, interp='linear') + + # define a second k-space grid using the dimensions of p_xy + N_recon = Vector(p_xy.shape) + d_recon = Vector([kgrid.dt * medium.sound_speed.item(), kgrid.dy]) + kgrid_recon = kWaveGrid(N_recon, d_recon) + + # resample p_xy to be the same size as source.p0 + interp_func = RegularGridInterpolator( + (kgrid_recon.x_vec[:, 0] - kgrid_recon.x_vec.min(), kgrid_recon.y_vec[:, 0]), + p_xy, method='linear' + ) + query_points = np.stack((kgrid.x - kgrid.x.min(), kgrid.y), axis=-1) + p_xy_rs = interp_func(query_points) + + + # -------------------- + # VISUALIZATION + # -------------------- + + cmap = get_color_map() + + # Plot the initial pressure and sensor distribution + fig, ax = plt.subplots(1, 1) + im = ax.imshow(p0 + sensor.mask * disc_magnitude, + extent=[kgrid.y_vec.min() * 1e3, kgrid.y_vec.max() * 1e3, kgrid.x_vec.max() * 1e3, kgrid.x_vec.min() * 1e3], + vmin=-disc_magnitude, vmax=disc_magnitude, cmap=cmap) + divider = make_axes_locatable(ax) + cax = divider.append_axes("right", size="3%", pad="2%") + ax.set_ylabel('x-position [mm]') + ax.set_xlabel('y-position [mm]') + fig.colorbar(im, cax=cax) + plt.show() + + # Plot the reconstructed initial pressure + fig, ax = plt.subplots(1, 1) + im = ax.imshow(p0_recon, + extent=[kgrid.y_vec.min() * 1e3, kgrid.y_vec.max() * 1e3, kgrid.x_vec.max() * 1e3, kgrid.x_vec.min() * 1e3], + vmin=-disc_magnitude, vmax=disc_magnitude, cmap=cmap) + divider = make_axes_locatable(ax) + cax = divider.append_axes("right", size="3%", pad="2%") + ax.set_ylabel('x-position [mm]') + ax.set_xlabel('y-position [mm]') + fig.colorbar(im, cax=cax) + plt.show() + + # Apply a positivity condition + p0_recon[p0_recon < 0] = 0 + + # Plot the reconstructed initial pressure with positivity condition + fig, ax = plt.subplots(1, 1) + im = ax.imshow(p0_recon, + extent=[kgrid.y_vec.min() * 1e3, kgrid.y_vec.max() * 1e3, kgrid.x_vec.max() * 1e3, kgrid.x_vec.min() * 1e3], + vmin=-disc_magnitude, vmax=disc_magnitude, cmap=cmap) + divider = make_axes_locatable(ax) + cax = divider.append_axes("right", size="3%", pad="2%") + ax.set_ylabel('x-position [mm]') + ax.set_xlabel('y-position [mm]') + fig.colorbar(im, cax=cax) + plt.show() + + # Plot a profile for comparison + plt.plot(kgrid.y_vec[:, 0] * 1e3, p0[disc_pos[0], :], 'k-', label='Initial Pressure') + plt.plot(kgrid.y_vec[:, 0] * 1e3, p_xy_rs[disc_pos[0], :], 'r--', label='FFT Reconstruction') + plt.plot(kgrid.y_vec[:, 0] * 1e3, p0_recon[disc_pos[0], :], 'b:', label='Time Reversal') + plt.xlabel('y-position [mm]') + plt.ylabel('Pressure') + plt.legend() + plt.axis('tight') + plt.ylim([0, 5.1]) + plt.show() + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/examples/pr_3D_FFT_planar_sensor/README.md b/examples/pr_3D_FFT_planar_sensor/README.md new file mode 100644 index 00000000..98ec2623 --- /dev/null +++ b/examples/pr_3D_FFT_planar_sensor/README.md @@ -0,0 +1,7 @@ +# 3D FFT Reconstruction For A Planar 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_3D_FFT_planar_sensor/pr_3D_FFT_planar_sensor.ipynb) + +This example demonstrates the use of k-Wave for the reconstruction of a three-dimensional photoacoustic wave-field recorded over a planar array of sensor elements. The sensor data is simulated using [kspaceFirstOrder3D](https://k-wave-python.readthedocs.io/en/latest/kwave.kspaceFirstOrder3D.html) and reconstructed using kspacePlaneRecon. It builds on the Simulations In Three Dimensions and [2D FFT Reconstruction For A Line Sensor](../pr_3D_FFT_planar_sensor/) examples. + +To read more, visit the [original example page](http://www.k-wave.org/documentation/example_pr_3D_fft_planar_sensor.php) or its [implementation](https://github.com/ucl-bug/k-wave/blob/main/k-Wave/examples/example_pr_3D_FFT_planar_sensor.m). \ No newline at end of file diff --git a/examples/pr_3D_FFT_planar_sensor/pr_3D_FFT_planar_sensor.ipynb b/examples/pr_3D_FFT_planar_sensor/pr_3D_FFT_planar_sensor.ipynb new file mode 100644 index 00000000..dbb4bf97 --- /dev/null +++ b/examples/pr_3D_FFT_planar_sensor/pr_3D_FFT_planar_sensor.ipynb @@ -0,0 +1,310 @@ +{ + "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": [ + "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.kspaceFirstOrder3D import kspaceFirstOrder3D\n", + "from kwave.kspacePlaneRecon import kspacePlaneRecon\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_ball\n", + "from kwave.utils.filters import smooth\n", + "\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "from scipy.interpolate import RegularGridInterpolator" + ], + "id": "611488037b162a21", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": [ + "## 3D FFT Reconstruction For A Planar Sensor Example\n", + "\n", + "This example demonstrates the use of k-Wave for the reconstruction of a\n", + "three-dimensional photoacoustic wave-field recorded over a planar array\n", + "of sensor elements. The sensor data is simulated using kspaceFirstOrder3D\n", + "and reconstructed using kspacePlaneRecon. It builds on the Simulations In\n", + "Three Dimensions and 2D FFT Reconstruction For A Line Sensor examples." + ], + "id": "d7d832f9f7f7caf3" + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": "### SIMULATION", + "id": "6e5153f4bb9a8e8b" + }, + { + "metadata": {}, + "cell_type": "code", + "source": "scale = 1", + "id": "5ef3d745e587ea7d", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "cell_type": "code", + "source": [ + "# create the computational grid\n", + "PML_size = 10 # size of the PML in grid points\n", + "N = Vector([32, 64, 64]) * scale - 2 * PML_size # number of grid points\n", + "d = Vector([0.2e-3, 0.2e-3, 0.2e-3]) / scale # grid point spacing [m]\n", + "kgrid = kWaveGrid(N, d)" + ], + "id": "62a1b78d509423b1", + "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": "9606aea53baa8a1", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "cell_type": "code", + "source": [ + "# create initial pressure distribution using makeBall\n", + "ball_magnitude = 10 # [Pa]\n", + "ball_radius = 3 * scale # [grid points]\n", + "p0 = ball_magnitude * make_ball(N, N / 2, ball_radius)\n", + "p0 = smooth(p0, restore_max=True)\n", + "\n", + "source = kSource()\n", + "source.p0 = p0" + ], + "id": "a5ab847f521ec9b9", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "cell_type": "code", + "source": [ + "# define a binary planar sensor\n", + "sensor = kSensor()\n", + "sensor_mask = np.zeros(N)\n", + "sensor_mask[0] = 1\n", + "sensor.mask = sensor_mask" + ], + "id": "80a91cfee6c4b79f", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "cell_type": "code", + "source": [ + "%%capture\n", + "\n", + "# create the time array\n", + "kgrid.makeTime(medium.sound_speed)" + ], + "id": "39e4252f2d55caae", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "cell_type": "code", + "source": [ + "# set the input arguments\n", + "simulation_options = SimulationOptions(\n", + " save_to_disk=True,\n", + " pml_size=PML_size,\n", + " pml_inside=False,\n", + " smooth_p0=False,\n", + " data_cast='single'\n", + ")\n", + "\n", + "execution_options = SimulationExecutionOptions(is_gpu_simulation=True)" + ], + "id": "d036c62f78598f2d", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "cell_type": "code", + "source": [ + "# run the simulation\n", + "sensor_data = kspaceFirstOrder3D(kgrid, source, sensor, medium, simulation_options, execution_options)\n", + "sensor_data = sensor_data['p'].T\n", + "\n", + "# reshape sensor data to y, z, t\n", + "sensor_data_rs = sensor_data.reshape(N[1], N[2], kgrid.Nt)\n", + "\n", + "# reconstruct the initial pressure\n", + "p_xyz = kspacePlaneRecon(sensor_data_rs, kgrid.dy, kgrid.dz, kgrid.dt.item(),\n", + " medium.sound_speed.item(), data_order='yzt', pos_cond=True)" + ], + "id": "e74cd39b11526d81", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "cell_type": "code", + "source": [ + "# define a k-space grid using the dimensions of p_xyz\n", + "N_recon = Vector(p_xyz.shape)\n", + "d_recon = Vector([kgrid.dt.item() * medium.sound_speed.item(), kgrid.dy, kgrid.dz])\n", + "kgrid_recon = kWaveGrid(N_recon, d_recon)\n", + "\n", + "# define a k-space grid with the same z-spacing as p0\n", + "kgrid_interp = kWaveGrid(N, d)\n", + "\n", + "# resample the p_xyz to be the same size as p0\n", + "interp_func = RegularGridInterpolator(\n", + " (kgrid_recon.x_vec[:, 0] - kgrid_recon.x_vec[:, 0].min(),\n", + " kgrid_recon.y_vec[:, 0] - kgrid_recon.y_vec[:, 0].min(),\n", + " kgrid_recon.z_vec[:, 0] - kgrid_recon.z_vec[:, 0].min()),\n", + " p_xyz, method='linear'\n", + ")\n", + "query_points = np.stack((kgrid_interp.x - kgrid_interp.x.min(),\n", + " kgrid_interp.y - kgrid_interp.y.min(),\n", + " kgrid_interp.z - kgrid_interp.z.min()),\n", + " axis=-1)\n", + "p_xyz_rs = interp_func(query_points)" + ], + "id": "c0d05511c962ede8", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": "### VISUALIZATION", + "id": "c63150dfb5bc2f2f" + }, + { + "metadata": {}, + "cell_type": "code", + "source": [ + "# plot the initial pressure and sensor surface in voxel form\n", + "# from kwave.utils.plot import voxel_plot\n", + "# voxel_plot(np.single((p0 + sensor_mask) > 0)) # todo: needs unsmoothed po + plot not working" + ], + "id": "d612fbaa36e986a9", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "cell_type": "code", + "source": [ + "# plot the initial pressure\n", + "plot_scale = [-10, 10]\n", + "fig, axs = plt.subplots(2, 2)\n", + "axs[0, 0].imshow(p0[:, :, N[2] // 2],\n", + " extent=[kgrid_interp.y_vec.min() * 1e3, kgrid_interp.y_vec.max() * 1e3,\n", + " kgrid_interp.x_vec.max() * 1e3, kgrid_interp.x_vec.min() * 1e3],\n", + " vmin=plot_scale[0], vmax=plot_scale[1], cmap=get_color_map())\n", + "axs[0, 0].set_title('x-y plane')\n", + "\n", + "axs[0, 1].imshow(p0[:, N[1] // 2, :],\n", + " extent=[kgrid_interp.z_vec.min() * 1e3, kgrid_interp.z_vec.max() * 1e3,\n", + " kgrid_interp.x_vec.max() * 1e3, kgrid_interp.x_vec.min() * 1e3],\n", + " vmin=plot_scale[0], vmax=plot_scale[1], cmap=get_color_map())\n", + "axs[0, 1].set_title('x-z plane')\n", + "\n", + "axs[1, 0].imshow(p0[N[0] // 2, :, :],\n", + " extent=[kgrid_interp.z_vec.min() * 1e3, kgrid_interp.z_vec.max() * 1e3,\n", + " kgrid_interp.y_vec.max() * 1e3, kgrid_interp.y_vec.min() * 1e3],\n", + " vmin=plot_scale[0], vmax=plot_scale[1], cmap=get_color_map())\n", + "axs[1, 0].set_title('y-z plane')\n", + "\n", + "axs[1, 1].axis('off')\n", + "axs[1, 1].set_title('(All axes in mm)')\n", + "plt.show()" + ], + "id": "2db2c76fa8f0b234", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "cell_type": "code", + "source": [ + "# plot the reconstructed initial pressure\n", + "fig, axs = plt.subplots(2, 2)\n", + "axs[0, 0].imshow(p_xyz_rs[:, :, N[2] // 2],\n", + " extent=[kgrid_interp.y_vec.min() * 1e3, kgrid_interp.y_vec.max() * 1e3,\n", + " kgrid_interp.x_vec.max() * 1e3, kgrid_interp.x_vec.min() * 1e3],\n", + " vmin=plot_scale[0], vmax=plot_scale[1], cmap=get_color_map())\n", + "axs[0, 0].set_title('x-y plane')\n", + "\n", + "axs[0, 1].imshow(p_xyz_rs[:, N[1] // 2, :],\n", + " extent=[kgrid_interp.z_vec.min() * 1e3, kgrid_interp.z_vec.max() * 1e3,\n", + " kgrid_interp.x_vec.max() * 1e3, kgrid_interp.x_vec.min() * 1e3],\n", + " vmin=plot_scale[0], vmax=plot_scale[1], cmap=get_color_map())\n", + "axs[0, 1].set_title('x-z plane')\n", + "\n", + "axs[1, 0].imshow(p_xyz_rs[N[0] // 2, :, :],\n", + " extent=[kgrid_interp.z_vec.min() * 1e3, kgrid_interp.z_vec.max() * 1e3,\n", + " kgrid_interp.y_vec.max() * 1e3, kgrid_interp.y_vec.min() * 1e3],\n", + " vmin=plot_scale[0], vmax=plot_scale[1], cmap=get_color_map())\n", + "axs[1, 0].set_title('y-z plane')\n", + "\n", + "axs[1, 1].axis('off')\n", + "axs[1, 1].set_title('(All axes in mm)')\n", + "plt.show()" + ], + "id": "d5f54008217d6532", + "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 +} diff --git a/examples/pr_3D_FFT_planar_sensor/pr_3D_FFT_planar_sensor.py b/examples/pr_3D_FFT_planar_sensor/pr_3D_FFT_planar_sensor.py new file mode 100644 index 00000000..aca56fa6 --- /dev/null +++ b/examples/pr_3D_FFT_planar_sensor/pr_3D_FFT_planar_sensor.py @@ -0,0 +1,163 @@ +from kwave.data import Vector +from kwave.kgrid import kWaveGrid +from kwave.kmedium import kWaveMedium +from kwave.ksensor import kSensor +from kwave.ksource import kSource +from kwave.kspaceFirstOrder3D import kspaceFirstOrder3D +from kwave.kspacePlaneRecon import kspacePlaneRecon +from kwave.options.simulation_execution_options import SimulationExecutionOptions +from kwave.options.simulation_options import SimulationOptions +from kwave.utils.colormap import get_color_map +from kwave.utils.mapgen import make_ball +from kwave.utils.filters import smooth + +import matplotlib.pyplot as plt +import numpy as np +from scipy.interpolate import RegularGridInterpolator + + +# 3D FFT Reconstruction For A Planar Sensor Example + +# This example demonstrates the use of k-Wave for the reconstruction of a +# three-dimensional photoacoustic wave-field recorded over a planar array +# of sensor elements. The sensor data is simulated using kspaceFirstOrder3D +# and reconstructed using kspacePlaneRecon. It builds on the Simulations In +# Three Dimensions and 2D FFT Reconstruction For A Line Sensor examples. + + +def main(): + + # -------------------- + # SIMULATION + # -------------------- + + scale = 1 + + # create the computational grid + PML_size = 10 # size of the PML in grid points + N = Vector([32, 64, 64]) * scale - 2 * PML_size # number of grid points + d = Vector([0.2e-3, 0.2e-3, 0.2e-3]) / scale # grid point spacing [m] + kgrid = kWaveGrid(N, d) + + # define the properties of the propagation medium + medium = kWaveMedium(sound_speed=1500) # [m/s] + + # create initial pressure distribution using makeBall + ball_magnitude = 10 # [Pa] + ball_radius = 3 * scale # [grid points] + p0 = ball_magnitude * make_ball(N, N / 2, ball_radius) + p0 = smooth(p0, restore_max=True) + + source = kSource() + source.p0 = p0 + + # define a binary planar sensor + sensor = kSensor() + sensor_mask = np.zeros(N) + sensor_mask[0] = 1 + sensor.mask = sensor_mask + + # set the input arguments + simulation_options = SimulationOptions( + save_to_disk=True, + pml_size=PML_size, + pml_inside=False, + smooth_p0=False, + data_cast='single' + ) + + execution_options = SimulationExecutionOptions(is_gpu_simulation=True) + + # run the simulation + sensor_data = kspaceFirstOrder3D(kgrid, source, sensor, medium, simulation_options, execution_options) + sensor_data = sensor_data['p'].T + + # reshape sensor data to y, z, t + sensor_data_rs = sensor_data.reshape(N[1], N[2], kgrid.Nt) + + # reconstruct the initial pressure + p_xyz = kspacePlaneRecon(sensor_data_rs, kgrid.dy, kgrid.dz, kgrid.dt.item(), + medium.sound_speed.item(), data_order='yzt', pos_cond=True) + + # define a k-space grid using the dimensions of p_xyz + N_recon = Vector(p_xyz.shape) + d_recon = Vector([kgrid.dt.item() * medium.sound_speed.item(), kgrid.dy, kgrid.dz]) + kgrid_recon = kWaveGrid(N_recon, d_recon) + + # define a k-space grid with the same z-spacing as p0 + kgrid_interp = kWaveGrid(N, d) + + # resample the p_xyz to be the same size as p0 + interp_func = RegularGridInterpolator( + (kgrid_recon.x_vec[:, 0] - kgrid_recon.x_vec[:, 0].min(), + kgrid_recon.y_vec[:, 0] - kgrid_recon.y_vec[:, 0].min(), + kgrid_recon.z_vec[:, 0] - kgrid_recon.z_vec[:, 0].min()), + p_xyz, method='linear' + ) + query_points = np.stack((kgrid_interp.x - kgrid_interp.x.min(), + kgrid_interp.y - kgrid_interp.y.min(), + kgrid_interp.z - kgrid_interp.z.min()), + axis=-1) + p_xyz_rs = interp_func(query_points) + + + # -------------------- + # VISUALIZATION + # -------------------- + + # plot the initial pressure and sensor surface in voxel form + # from kwave.utils.plot import voxel_plot + # voxel_plot(np.single((p0 + sensor_mask) > 0)) # todo: needs unsmoothed po + plot not working + + # plot the initial pressure + plot_scale = [-10, 10] + fig, axs = plt.subplots(2, 2) + axs[0, 0].imshow(p0[:, :, N[2] // 2], + extent=[kgrid_interp.y_vec.min() * 1e3, kgrid_interp.y_vec.max() * 1e3, + kgrid_interp.x_vec.max() * 1e3, kgrid_interp.x_vec.min() * 1e3], + vmin=plot_scale[0], vmax=plot_scale[1], cmap=get_color_map()) + axs[0, 0].set_title('x-y plane') + + axs[0, 1].imshow(p0[:, N[1] // 2, :], + extent=[kgrid_interp.z_vec.min() * 1e3, kgrid_interp.z_vec.max() * 1e3, + kgrid_interp.x_vec.max() * 1e3, kgrid_interp.x_vec.min() * 1e3], + vmin=plot_scale[0], vmax=plot_scale[1], cmap=get_color_map()) + axs[0, 1].set_title('x-z plane') + + axs[1, 0].imshow(p0[N[0] // 2, :, :], + extent=[kgrid_interp.z_vec.min() * 1e3, kgrid_interp.z_vec.max() * 1e3, + kgrid_interp.y_vec.max() * 1e3, kgrid_interp.y_vec.min() * 1e3], + vmin=plot_scale[0], vmax=plot_scale[1], cmap=get_color_map()) + axs[1, 0].set_title('y-z plane') + + axs[1, 1].axis('off') + axs[1, 1].set_title('(All axes in mm)') + plt.show() + + # plot the reconstructed initial pressure + fig, axs = plt.subplots(2, 2) + axs[0, 0].imshow(p_xyz_rs[:, :, N[2] // 2], + extent=[kgrid_interp.y_vec.min() * 1e3, kgrid_interp.y_vec.max() * 1e3, + kgrid_interp.x_vec.max() * 1e3, kgrid_interp.x_vec.min() * 1e3], + vmin=plot_scale[0], vmax=plot_scale[1], cmap=get_color_map()) + axs[0, 0].set_title('x-y plane') + + axs[0, 1].imshow(p_xyz_rs[:, N[1] // 2, :], + extent=[kgrid_interp.z_vec.min() * 1e3, kgrid_interp.z_vec.max() * 1e3, + kgrid_interp.x_vec.max() * 1e3, kgrid_interp.x_vec.min() * 1e3], + vmin=plot_scale[0], vmax=plot_scale[1], cmap=get_color_map()) + axs[0, 1].set_title('x-z plane') + + axs[1, 0].imshow(p_xyz_rs[N[0] // 2, :, :], + extent=[kgrid_interp.z_vec.min() * 1e3, kgrid_interp.z_vec.max() * 1e3, + kgrid_interp.y_vec.max() * 1e3, kgrid_interp.y_vec.min() * 1e3], + vmin=plot_scale[0], vmax=plot_scale[1], cmap=get_color_map()) + axs[1, 0].set_title('y-z plane') + + axs[1, 1].axis('off') + axs[1, 1].set_title('(All axes in mm)') + plt.show() + + +if __name__ == '__main__': + main() \ No newline at end of file diff --git a/examples/pr_3D_TR_planar_sensor/README.md b/examples/pr_3D_TR_planar_sensor/README.md new file mode 100644 index 00000000..49c3a43c --- /dev/null +++ b/examples/pr_3D_TR_planar_sensor/README.md @@ -0,0 +1,7 @@ +# 3D Time Reversal Reconstruction For A Planar 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_3D_TR_planar_sensor/pr_3D_TR_planar_sensor.ipynb) + +This example demonstrates the use of k-Wave for the reconstruction of a three-dimensional photoacoustic wave-field recorded over a planar array of sensor elements. The sensor data is simulated and then time-reversed using [kspaceFirstOrder3D](https://k-wave-python.readthedocs.io/en/latest/kwave.kspaceFirstOrder3D.html). It builds on the [3D FFT Reconstruction For A Planar Sensor](../pr_3D_FFT_planar_sensor/) and [2D Time Reversal Reconstruction For A Line Sensor](../pr_2D_TR_line_sensor) examples. + +To read more, visit the [original example page](http://www.k-wave.org/documentation/example_pr_3D_tr_planar_sensor.php) or its [implementation](https://github.com/ucl-bug/k-wave/blob/main/k-Wave/examples/example_pr_3D_TR_planar_sensor.m). \ No newline at end of file diff --git a/examples/pr_3D_TR_planar_sensor/pr_3D_TR_planar_sensor.ipynb b/examples/pr_3D_TR_planar_sensor/pr_3D_TR_planar_sensor.ipynb new file mode 100644 index 00000000..31d8a5d4 --- /dev/null +++ b/examples/pr_3D_TR_planar_sensor/pr_3D_TR_planar_sensor.ipynb @@ -0,0 +1,302 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "initial_id", + "metadata": {}, + "outputs": [], + "source": [ + "%%capture\n", + "!pip install k-wave-python " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6a6cdd4e64365d3e", + "metadata": {}, + "outputs": [], + "source": [ + "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.kspaceFirstOrder3D import kspaceFirstOrder3D\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_ball\n", + "from kwave.utils.filters import smooth\n", + "\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np" + ] + }, + { + "cell_type": "markdown", + "id": "e7d1f7d7d6481a2c", + "metadata": {}, + "source": [ + "## 3D Time Reversal Reconstruction For A Planar Sensor Example\n", + "\n", + "This example demonstrates the use of k-Wave for the reconstruction of a three-dimensional photoacoustic wave-field recorded over a planar array of sensor elements. The sensor data is simulated and then time-reversed using kspaceFirstOrder3D. It builds on the 3D FFT Reconstruction For A Planar Sensor and 2D Time Reversal Reconstruction For A Line Sensor examples." + ] + }, + { + "cell_type": "markdown", + "id": "4a805d6488afc4b7", + "metadata": {}, + "source": [ + "### SIMULATION" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b8f9b1f3d5bef716", + "metadata": {}, + "outputs": [], + "source": [ + "# change scale to 2 to reproduce the higher resolution figures used in the\n", + "# help file\n", + "scale = 1" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "249a1bd033f3abe7", + "metadata": {}, + "outputs": [], + "source": [ + "# create the computational grid\n", + "PML_size = 10 # size of the PML in grid points\n", + "N = Vector([32, 64, 64]) * scale - 2 * PML_size # number of grid points\n", + "d = Vector([0.2e-3, 0.2e-3, 0.2e-3]) / scale # grid point spacing [m]\n", + "kgrid = kWaveGrid(N, d)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "25d647c865e569b", + "metadata": {}, + "outputs": [], + "source": [ + "# define the properties of the propagation medium\n", + "medium = kWaveMedium(sound_speed=1500) # [m/s]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cb75ff29974130cf", + "metadata": {}, + "outputs": [], + "source": [ + "# create initial pressure distribution using makeBall\n", + "ball_magnitude = 10 # [Pa]\n", + "ball_radius = 3 * scale # [grid points]\n", + "p0 = ball_magnitude * make_ball(N, N / 2, ball_radius)\n", + "\n", + "# smooth the initial pressure distribution and restore the magnitude\n", + "p0 = smooth(p0, True)\n", + "\n", + "# assign to the source structure\n", + "source = kSource()\n", + "source.p0 = p0" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9eedddfe4204cea5", + "metadata": {}, + "outputs": [], + "source": [ + "# define a binary planar sensor\n", + "sensor = kSensor()\n", + "sensor.mask = np.zeros(N)\n", + "sensor.mask[0] = 1" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7ea0924151a96a6f", + "metadata": {}, + "outputs": [], + "source": [ + "%%capture\n", + "\n", + "# create the time array\n", + "kgrid.makeTime(medium.sound_speed)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7ad7415e71716b05", + "metadata": {}, + "outputs": [], + "source": [ + "# set the input arguements\n", + "simulation_options = SimulationOptions(\n", + "save_to_disk=True,\n", + "pml_size=PML_size,\n", + "pml_inside=False,\n", + "smooth_p0=False,\n", + "data_cast='single',\n", + ")\n", + "execution_options = SimulationExecutionOptions(is_gpu_simulation=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "192d9eeefd209077", + "metadata": {}, + "outputs": [], + "source": [ + "# run the simulation\n", + "sensor_data = kspaceFirstOrder3D(kgrid, source, sensor, medium, simulation_options, execution_options)\n", + "sensor_data = sensor_data['p'].T" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "790676732e5f078c", + "metadata": {}, + "outputs": [], + "source": [ + "# reset the initial pressure\n", + "source = kSource()\n", + "sensor = kSensor()\n", + "sensor.mask = np.zeros(N)\n", + "sensor.mask[0] = 1\n", + "\n", + "# assign the time reversal data\n", + "sensor.time_reversal_boundary_data = sensor_data\n", + "\n", + "# run the time-reversal reconstruction\n", + "p0_recon = kspaceFirstOrder3D(kgrid, source, sensor, medium, simulation_options, execution_options)\n", + "p0_recon = p0_recon['p_final'].T\n", + "\n", + "# add first order compensation for only recording over a half plane\n", + "p0_recon = 2 * p0_recon\n", + "\n", + "# apply a positivity condition\n", + "p0_recon[p0_recon < 0] = 0" + ] + }, + { + "cell_type": "markdown", + "id": "7bc3a5e3b2f4571", + "metadata": {}, + "source": [ + "### VISUALIZATION" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7f16ca4a29e4b671", + "metadata": {}, + "outputs": [], + "source": [ + "cmap = get_color_map()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ac432a6b8e284bda", + "metadata": {}, + "outputs": [], + "source": [ + "plot_scale = [-10, 10]\n", + "\n", + "# plot the initial pressure\n", + "fig, axs = plt.subplots(2, 2)\n", + "axs[0, 0].imshow(p0[:, :, N[2] // 2],\n", + " extent=[kgrid.y_vec.min() * 1e3, kgrid.y_vec.max() * 1e3, kgrid.x_vec.max() * 1e3, kgrid.x_vec.min() * 1e3],\n", + " vmin=plot_scale[0], vmax=plot_scale[1], cmap=cmap)\n", + "axs[0, 0].set_title('x-y plane')\n", + "axs[0, 0].axis('image')\n", + "\n", + "axs[0, 1].imshow(p0[:, N[1] // 2, :],\n", + " extent=[kgrid.z_vec.min() * 1e3, kgrid.z_vec.max() * 1e3, kgrid.x_vec.max() * 1e3, kgrid.x_vec.min() * 1e3],\n", + " vmin=plot_scale[0], vmax=plot_scale[1], cmap=cmap)\n", + "axs[0, 1].set_title('x-z plane')\n", + "axs[0, 1].axis('image')\n", + "\n", + "axs[1, 0].imshow(p0[N[0] // 2, :, :],\n", + " extent=[kgrid.z_vec.min() * 1e3, kgrid.z_vec.max() * 1e3, kgrid.y_vec.max() * 1e3, kgrid.y_vec.min() * 1e3],\n", + " vmin=plot_scale[0], vmax=plot_scale[1], cmap=cmap)\n", + "axs[1, 0].set_title('y-z plane')\n", + "axs[1, 0].axis('image')\n", + "\n", + "axs[1, 1].axis('off')\n", + "axs[1, 1].set_title('(All axes in mm)')\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e946817f00b55e9e", + "metadata": {}, + "outputs": [], + "source": [ + "# plot the reconstructed initial pressure\n", + "fig, axs = plt.subplots(2, 2)\n", + "axs[0, 0].imshow(p0_recon[:, :, N[2] // 2],\n", + " extent=[kgrid.y_vec.min() * 1e3, kgrid.y_vec.max() * 1e3, kgrid.x_vec.max() * 1e3, kgrid.x_vec.min() * 1e3],\n", + " vmin=plot_scale[0], vmax=plot_scale[1], cmap=cmap)\n", + "axs[0, 0].set_title('x-y plane')\n", + "axs[0, 0].axis('image')\n", + "\n", + "axs[0, 1].imshow(p0_recon[:, N[1] // 2, :],\n", + " extent=[kgrid.z_vec.min() * 1e3, kgrid.z_vec.max() * 1e3, kgrid.x_vec.max() * 1e3, kgrid.x_vec.min() * 1e3],\n", + " vmin=plot_scale[0], vmax=plot_scale[1], cmap=cmap)\n", + "axs[0, 1].set_title('x-z plane')\n", + "axs[0, 1].axis('image')\n", + "\n", + "axs[1, 0].imshow(p0_recon[N[0] // 2, :, :],\n", + " extent=[kgrid.z_vec.min() * 1e3, kgrid.z_vec.max() * 1e3, kgrid.y_vec.max() * 1e3, kgrid.y_vec.min() * 1e3],\n", + " vmin=plot_scale[0], vmax=plot_scale[1], cmap=cmap)\n", + "axs[1, 0].set_title('y-z plane')\n", + "axs[1, 0].axis('image')\n", + "\n", + "axs[1, 1].axis('off')\n", + "axs[1, 1].set_title('(All axes in mm)')\n", + "\n", + "plt.show()" + ] + } + ], + "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" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/examples/pr_3D_TR_planar_sensor/pr_3D_TR_planar_sensor.py b/examples/pr_3D_TR_planar_sensor/pr_3D_TR_planar_sensor.py new file mode 100644 index 00000000..82b4e110 --- /dev/null +++ b/examples/pr_3D_TR_planar_sensor/pr_3D_TR_planar_sensor.py @@ -0,0 +1,176 @@ +from kwave.data import Vector +from kwave.kgrid import kWaveGrid +from kwave.kmedium import kWaveMedium +from kwave.ksensor import kSensor +from kwave.ksource import kSource +from kwave.kspaceFirstOrder3D import kspaceFirstOrder3D +from kwave.options.simulation_execution_options import SimulationExecutionOptions +from kwave.options.simulation_options import SimulationOptions +from kwave.utils.colormap import get_color_map +from kwave.utils.mapgen import make_ball +from kwave.utils.filters import smooth + +import matplotlib.pyplot as plt +import numpy as np + + +# 3D Time Reversal Reconstruction For A Planar Sensor Example + +# This example demonstrates the use of k-Wave for the reconstruction +# of a three-dimensional photoacoustic wave-field recorded over a planar +# array of sensor elements. The sensor data is simulated and then +# time-reversed using kspaceFirstOrder3D. It builds on the 3D FFT +# Reconstruction For A Planar Sensor and 2D Time Reversal Reconstruction +# For A Line Sensor examples. + + +def main(): + # SIMULATION + + # change scale to 2 to reproduce the higher resolution figures used in the + # help file + scale = 1 + + # create the computational grid + PML_size = 10 # size of the PML in grid points + N = Vector([32, 64, 64]) * scale - 2 * PML_size # number of grid points + d = Vector([0.2e-3, 0.2e-3, 0.2e-3]) / scale # grid point spacing [m] + kgrid = kWaveGrid(N, d) + + # define the properties of the propagation medium + medium = kWaveMedium(sound_speed=1500) # [m/s] + + # create initial pressure distribution using makeBall + ball_magnitude = 10 # [Pa] + ball_radius = 3 * scale # [grid points] + p0 = ball_magnitude * make_ball(N, N / 2, ball_radius) + + # smooth the initial pressure distribution and restore the magnitude + p0 = smooth(p0, True) + + # assign to the source structure + source = kSource() + source.p0 = p0 + + # define a binary planar sensor + sensor = kSensor() + sensor.mask = np.zeros(N) + sensor.mask[0] = 1 + + # create the time array + kgrid.makeTime(medium.sound_speed) + + # set the input arguements + simulation_options = SimulationOptions( + save_to_disk=True, + pml_size=PML_size, + pml_inside=False, + smooth_p0=False, + data_cast="single", + ) + execution_options = SimulationExecutionOptions(is_gpu_simulation=True) + + # run the simulation + sensor_data = kspaceFirstOrder3D(kgrid, source, sensor, medium, simulation_options, execution_options) + sensor_data = sensor_data["p"].T + + # reset the initial pressure + source = kSource() + sensor = kSensor() + sensor.mask = np.zeros(N) + sensor.mask[0] = 1 + + # assign the time reversal data + sensor.time_reversal_boundary_data = sensor_data + + # run the time-reversal reconstruction + p0_recon = kspaceFirstOrder3D(kgrid, source, sensor, medium, simulation_options, execution_options) + p0_recon = p0_recon["p_final"].T + + # add first order compensation for only recording over a half plane + p0_recon = 2 * p0_recon + + # apply a positivity condition + p0_recon[p0_recon < 0] = 0 + + # VISUALIZATION + cmap = get_color_map() + plot_scale = [-10, 10] + + # plot the initial pressure + fig, axs = plt.subplots(2, 2) + axs[0, 0].imshow( + p0[:, :, N[2] // 2], + extent=[kgrid.y_vec.min() * 1e3, kgrid.y_vec.max() * 1e3, kgrid.x_vec.max() * 1e3, kgrid.x_vec.min() * 1e3], + vmin=plot_scale[0], + vmax=plot_scale[1], + cmap=cmap, + ) + axs[0, 0].set_title("x-y plane") + axs[0, 0].axis("image") + + axs[0, 1].imshow( + p0[:, N[1] // 2, :], + extent=[kgrid.z_vec.min() * 1e3, kgrid.z_vec.max() * 1e3, kgrid.x_vec.max() * 1e3, kgrid.x_vec.min() * 1e3], + vmin=plot_scale[0], + vmax=plot_scale[1], + cmap=cmap, + ) + axs[0, 1].set_title("x-z plane") + axs[0, 1].axis("image") + + axs[1, 0].imshow( + p0[N[0] // 2, :, :], + extent=[kgrid.z_vec.min() * 1e3, kgrid.z_vec.max() * 1e3, kgrid.y_vec.max() * 1e3, kgrid.y_vec.min() * 1e3], + vmin=plot_scale[0], + vmax=plot_scale[1], + cmap=cmap, + ) + axs[1, 0].set_title("y-z plane") + axs[1, 0].axis("image") + + axs[1, 1].axis("off") + axs[1, 1].set_title("(All axes in mm)") + + plt.show() + + # plot the reconstructed initial pressure + fig, axs = plt.subplots(2, 2) + axs[0, 0].imshow( + p0_recon[:, :, N[2] // 2], + extent=[kgrid.y_vec.min() * 1e3, kgrid.y_vec.max() * 1e3, kgrid.x_vec.max() * 1e3, kgrid.x_vec.min() * 1e3], + vmin=plot_scale[0], + vmax=plot_scale[1], + cmap=cmap, + ) + axs[0, 0].set_title("x-y plane") + axs[0, 0].axis("image") + + axs[0, 1].imshow( + p0_recon[:, N[1] // 2, :], + extent=[kgrid.z_vec.min() * 1e3, kgrid.z_vec.max() * 1e3, kgrid.x_vec.max() * 1e3, kgrid.x_vec.min() * 1e3], + vmin=plot_scale[0], + vmax=plot_scale[1], + cmap=cmap, + ) + axs[0, 1].set_title("x-z plane") + axs[0, 1].axis("image") + + axs[1, 0].imshow( + p0_recon[N[0] // 2, :, :], + extent=[kgrid.z_vec.min() * 1e3, kgrid.z_vec.max() * 1e3, kgrid.y_vec.max() * 1e3, kgrid.y_vec.min() * 1e3], + vmin=plot_scale[0], + vmax=plot_scale[1], + cmap=cmap, + ) + axs[1, 0].set_title("y-z plane") + axs[1, 0].axis("image") + + axs[1, 1].axis("off") + axs[1, 1].set_title("(All axes in mm)") + + plt.show() + + +if __name__ == "__main__": + main() diff --git a/kwave/executor.py b/kwave/executor.py index de8cd4cb..77e10803 100644 --- a/kwave/executor.py +++ b/kwave/executor.py @@ -58,36 +58,58 @@ def run_simulation(self, input_filename: str, output_filename: str, options: lis raise sensor_data = self.parse_executable_output(output_filename) + if not self.simulation_options.pml_inside: + self._crop_pml(sensor_data) return sensor_data + def _crop_pml(self, sensor_data: dotdict): + Nx = sensor_data["Nx"].item() + Ny = sensor_data["Ny"].item() + Nz = sensor_data["Nz"].item() + pml_x_size = sensor_data["pml_x_size"].item() + pml_y_size = sensor_data["pml_y_size"].item() + pml_z_size = 0 if Nz <= 1 else sensor_data["pml_z_size"].item() + axisymmetric = sensor_data["axisymmetric_flag"].item() + + # if the PML is outside, set the index variables to remove the pml + # from the _all and _final variables + x1 = pml_x_size + x2 = Nx - pml_x_size + y1 = 0 if axisymmetric else pml_y_size + y2 = Ny - pml_y_size + z1 = pml_z_size + z2 = Nz - pml_z_size + + possible_fields = [ + "p_final", + "p_max_all", + "p_min_all", + "ux_max_all", + "uy_max_all", + "uz_max_all", + "ux_min_all", + "uy_min_all", + "uz_min_all", + "ux_final", + "uy_final", + "uz_final", + ] + for field in possible_fields: + if field in sensor_data: + if sensor_data[field].ndim == 2: + sensor_data[field] = sensor_data[field][y1:y2, x1:x2] + else: + sensor_data[field] = sensor_data[field][z1:z2, y1:y2, x1:x2] + @staticmethod def parse_executable_output(output_filename: str) -> dotdict: - # Load the simulation and pml sizes from the output file - # with h5py.File(output_filename, 'r') as output_file: - # Nx, Ny, Nz = output_file['/Nx'][0].item(), output_file['/Ny'][0].item(), output_file['/Nz'][0].item() - # pml_x_size, pml_y_size = output_file['/pml_x_size'][0].item(), output_file['/pml_y_size'][0].item() - # pml_z_size = output_file['/pml_z_size'][0].item() if Nz > 1 else 0 - - # # Set the default index variables for the _all and _final variables - # x1, x2 = 1, Nx - # y1, y2 = ( - # 1, 1 + pml_y_size) if self.simulation_options.simulation_type is not SimulationType.AXISYMMETRIC else ( - # 1, Ny) - # z1, z2 = (1 + pml_z_size, Nz - pml_z_size) if Nz > 1 else (1, Nz) - # - # # Check if the PML is set to be outside the computational grid - # if self.simulation_options.pml_inside: - # x1, x2 = 1 + pml_x_size, Nx - pml_x_size - # y1, y2 = (1, Ny) if self.simulation_options.simulation_type is SimulationType.AXISYMMETRIC else ( - # 1 + pml_y_size, Ny - pml_y_size) - # z1, z2 = 1 + pml_z_size, Nz - pml_z_size if Nz > 1 else (1, Nz) - # Load the C++ data back from disk using h5py with h5py.File(output_filename, "r") as output_file: sensor_data = {} for key in output_file.keys(): sensor_data[key] = output_file[f"/{key}"][:].squeeze() + # if self.simulation_options.cuboid_corners: # sensor_data = [output_file[f'/p/{index}'][()] for index in range(1, len(key['mask']) + 1)] # diff --git a/kwave/kWaveSimulation.py b/kwave/kWaveSimulation.py index f4f7d9c7..676a436a 100644 --- a/kwave/kWaveSimulation.py +++ b/kwave/kWaveSimulation.py @@ -44,19 +44,33 @@ def __init__( # FLAGS WHICH DEPEND ON USER INPUTS (THESE SHOULD NOT BE MODIFIED) # ========================================================================= # flags which control the characteristics of the sensor + #: Whether the sensor sensor mask is defined by cuboid corners #: Whether time reversal simulation is enabled + # check if the sensor mask is defined as a list of cuboid corners + if self.sensor.mask is not None and self.sensor.mask.shape[0] == (2 * self.kgrid.dim): + self.userarg_cuboid_corners = True + else: + self.userarg_cuboid_corners = False + # check if performing time reversal, and replace inputs to explicitly use a # source with a dirichlet boundary condition if self.sensor.time_reversal_boundary_data is not None: # define a new source structure - source = {"p_mask": self.sensor.p_mask, "p": np.flip(self.sensor.time_reversal_boundary_data, 2), "p_mode": "dirichlet"} + self.source = kSource() + self.source.p_mask = self.sensor.mask + self.source.p = np.flip(self.sensor.time_reversal_boundary_data, 1) + self.source.p_mode = "dirichlet" # define a new sensor structure Nx, Ny, Nz = self.kgrid.Nx, self.kgrid.Ny, self.kgrid.Nz - sensor = kSensor(mask=np.ones((Nx, Ny, max(1, Nz))), record=["p_final"]) + self.sensor = kSensor(mask=np.ones((Nx, Ny, max(1, Nz))), record=["p_final"]) + # set time reversal flag self.userarg_time_rev = True + + self.record = Recorder() + self.record.p = False else: # set time reversal flag self.userarg_time_rev = False @@ -69,12 +83,6 @@ def __init__( #: Whether the sensor.mask is binary self.binary_sensor_mask = True - # check if the sensor mask is defined as a list of cuboid corners - if self.sensor.mask is not None and self.sensor.mask.shape[0] == (2 * self.kgrid.dim): - self.userarg_cuboid_corners = True - else: - self.userarg_cuboid_corners = False - #: If tse sensor is an object of the kWaveTransducer class self.transducer_sensor = False diff --git a/kwave/kspaceLineRecon.py b/kwave/kspaceLineRecon.py new file mode 100644 index 00000000..e464a7f4 --- /dev/null +++ b/kwave/kspaceLineRecon.py @@ -0,0 +1,136 @@ +import logging + +from kwave.data import Vector +from kwave.kgrid import kWaveGrid + +import numpy as np +from numpy.fft import fftn, ifftn, fftshift, ifftshift +from scipy.interpolate import RegularGridInterpolator + + +def kspaceLineRecon( + p: np.ndarray, + dy: float, + dt: float, + c: float, + data_order: str = 'ty', + interp: str = 'nearest', + pos_cond: bool = False +) -> np.ndarray: + """ + kspaceLineRecon takes an acoustic pressure time-series p_ty recorded + over an evenly spaced array of sensor points on a line, and + constructs an estimate of the initial acoustic pressure distribution + that gave rise to those measurements using an algorithm based on the + FFT. The input p_ty must be indexed p_ty(time step, sensor y + position), where the sensor spacing is given by dy, the temporal + spacing given by dt, and the sound speed in the propagation medium + (which is assumed to be acoustically homogeneous) is given by c. The + output p_xy is indexed as p_xy(x position, y position). + + The code uses a k-space algorithm which performs (1) a Fourier + transform on the data p_ty along both t and y dimensions (into + wavenumber-frequency space, where the wavenumber component is along + the detector line), (2) a mapping, based on the dispersion relation + for a plane wave in an acoustically homogeneous medium, from + wavenumber-frequency space to wavenumber-wavenumber space (where the + second component is orthogonal to the detector line), and finally (3) + an inverse Fourier transform back from the wavenumber domain to the + spatial domain. The result is an estimate of the initial acoustic + pressure distribution from which the acoustic waves originated. + + Steps (1) and (3) can be performed efficiently using the fast Fourier + transform (FFT); they are therefore fastest when the number of + samples and number of detector points are both powers of 2. The + mapping in step (2) requires an interpolation of the data from an + evenly spaced grid of points in the wavenumber-frequency domain to + an evenly-spaced grid of points in the wavenumber-wavenumber domain. + The option 'Interp' may be used to choose the interpolation method. + + The physics of photoacoustics requires that the acoustic pressure is + initially non-negative everywhere. The estimate of the initial + pressure distribution generated by this code may have negative + regions due to artefacts arising from differences between the + assumed model and the real situation, e.g., homogeneous medium vs. + real, somewhat heterogeneous, medium; infinite measurement surface + vs. finite-sized region-of-detection, etc. A positivity (or + non-negativity) condition can be enforced by setting the optional + 'PosCond' to true which simply sets any negative parts of the final + image to zero. + """ + p = p.copy() + + # reorder the data if needed (p_ty) + if data_order == 'yt': + p = p.T + + # mirror the time domain data about t = 0 to allow the cosine transform to + # be computed using an FFT (p_ty) + p = np.vstack((np.flipud(p), p[1:])) + + # extract the size of mirrored input data + Nt, Ny = p.shape + + # update command line status + logging.log(logging.INFO, "Running k-Wave line reconstruction...\n" + f"grid size: {Ny} by {(Nt + 1) / 2} grid points\n" + f"interpolation mode: {interp}") + + # create a computational grid that is evenly spaced in w and ky, where + # Nx = Nt and dx = dt*c + N = Vector([Nt, Ny]) + d = Vector([dt * c, dy]) + kgrid = kWaveGrid(N, d) + + # from the grid for kx, create a computational grid for w using the + # relation dx = dt*c; this represents the initial sampling of p(w, ky) + w = c * kgrid.kx + + # remap the computational grid for kx onto w using the dispersion + # relation w/c = (kx^2 + ky^2)^1/2. This gives an w grid that is + # evenly spaced in kx. This is used for the interpolation from p(w, ky) + # to p(kx, ky). Only real w is taken to force kx (and thus x) to be + # symmetrical about 0 after the interpolation. + w_new = c * kgrid.k + + # calculate the scaling factor using the value of kx, where + # kx = sqrt( (w/c).^2 - kgrid.ky.^2 ) and then manually + # replacing the DC value with its limit (otherwise NaN results) + with np.errstate(divide='ignore', invalid='ignore'): + sf = c ** 2 * np.emath.sqrt((w / c) ** 2 - kgrid.ky ** 2) / (2 * w) + sf[(w == 0) & (kgrid.ky == 0)] = c / 2 + + # compute the FFT of the input data p(t, y) to yield p(w, ky) and scale + p = sf * fftshift(fftn(ifftshift(p))) + + # exclude the inhomogeneous part of the wave + p[np.abs(w) < np.abs(c * kgrid.ky)] = 0 + + # compute the interpolation from p(w, ky) to p(kx, ky)and then force to be + # symmetrical + interp_func = RegularGridInterpolator( + (w[:, 0], kgrid.ky[0]), + p, bounds_error=False, fill_value=0, method=interp + ) + query_points = np.stack((w_new, kgrid.ky), axis=-1) + p = interp_func(query_points) + + # compute the inverse FFT of p(kx, ky) to yield p(x, y) + p = np.real(fftshift(ifftn(ifftshift(p)))) + + # remove the left part of the mirrored data which corresponds to the + # negative part of the mirrored time data + p = p[(Nt // 2):, :] + + # correct the scaling - the forward FFT is computed with a spacing of dt + # and the reverse requires a spacing of dy = dt*c, the reconstruction + # assumes that p0 is symmetrical about y, and only half the plane collects + # data (first approximation to correcting the limited view problem) + p = 2 * 2 * p / c + + # enforce positivity condition + if pos_cond: + logging.log(logging.INFO, 'applying positivity condition...') + p[p < 0] = 0 + + return p diff --git a/kwave/kspacePlaneRecon.py b/kwave/kspacePlaneRecon.py new file mode 100644 index 00000000..b84ea435 --- /dev/null +++ b/kwave/kspacePlaneRecon.py @@ -0,0 +1,136 @@ +import logging + +from kwave.data import Vector +from kwave.kgrid import kWaveGrid + +import numpy as np +from numpy.fft import fftn, ifftn, fftshift, ifftshift +from scipy.interpolate import RegularGridInterpolator + + +def kspacePlaneRecon( + p: np.ndarray, + dy: float, + dz: float, + dt: float, + c: float, + data_order: str = 'tyz', + interp: str = 'nearest', + pos_cond: bool = False +): + """ + kspacePlaneRecon takes an acoustic pressure time-series p_tyz + recorded over an uniform array of sensor points on a plane, and + constructs an estimate of the initial acoustic pressure distribution + that gave rise to those measurements using an algorithm based on the + FFT. The input p_tyz must be indexed p_tyz(time step, sensor y + position, sensor z position), where the sensor spacing is given by dy + and dz, the temporal spacing given by dt, and the sound speed in the + propagation medium (which is assumed to be acoustically homogeneous) + is given by c. The output p_xyz is indexed as p_xyz(x position, y + position, z position). + + The code uses a k-space algorithm which performs (1) a Fourier + transform on the data p_tyz along both t, y, and z dimensions (into + wavenumber-frequency space), (2) a mapping, based on the dispersion + relation for a plane wave in an acoustically homogeneous medium, from + wavenumber-frequency space to wavenumber-wavenumber space, and + finally (3) an inverse Fourier transform back from the wavenumber + domain to the spatial domain. The result is an estimate of the + initial acoustic pressure distribution from which the acoustic waves + originated. + + Steps (1) and (3) can be performed efficiently using the fast Fourier + transform (FFT); they are therefore fastest when the number of + samples and number of detector points are both powers of 2. The + mapping in step (2) requires an interpolation of the data from an + evenly spaced grid of points in the wavenumber-frequency domain to + an evenly-spaced grid of points in the wavenumber-wavenumber + domain. The option 'Interp' may be used to choose the interpolation + method. + + The physics of photoacoustics requires that the acoustic pressure is + initially non-negative everywhere. The estimate of the initial + pressure distribution generated by this code may have negative + regions due to artefacts arising from differences between the + assumed model and the real situation, e.g., homogeneous medium vs. + real, somewhat heterogeneous, medium; infinite measurement surface + vs. finite-sized region-of-detection, etc. A positivity (or + non-negativity) condition can be enforced by setting the optional + 'PosCond' to true which simply sets any negative parts of the final + image to zero. + """ + p = p.copy() + + # reorder the data to p(t, y, z) if needed + if data_order == 'yzt': + p = np.transpose(p, (2, 0, 1)) + + # mirror the time domain data about t = 0 to allow the cosine transform in + # the t direction to be computed using an FFT + p = np.concatenate((np.flip(p, axis=0), p[1:, :, :]), axis=0) + + # extract the size of mirrored input data + Nt, Ny, Nz = p.shape + + # update command line status + logging.log(logging.INFO, "Running k-Wave planar reconstruction...\n" + f"grid size: {(Nt + 1) // 2} by {Ny} by {Nz} grid points\n" + f"interpolation mode: {interp}") + + # create a computational grid that is evenly spaced in w, ky, and kz, where + # Nx = Nt and dx = dt*c + N = Vector([Nt, Ny, Nz]) + d = Vector([dt * c, dy, dz]) + kgrid = kWaveGrid(N, d) + + # from the grid for kx, create a computational grid for w using the + # relation dx = dt*c; this represents the initial sampling of p(w, ky, kz) + w = c * kgrid.kx + + # remap the computational grid for kx onto w using the dispersion + # relation w/c = (kx^2 + ky^2 + kz^2)^1/2. This gives an w grid that is + # evenly spaced in kx. This is used for the interpolation from p(w, ky, kz) + # to p(kx, ky, kz). Only real w is taken to force kx (and thus x) to be + # symmetrical about 0 after the interpolation. + w_new = c * kgrid.k + + # calculate the scaling factor using the value of kx, where + # kx = sqrt( (w/c)^2 - kgrid.ky^2 - kgrid.kz^2 ) and then manually + # replacing the DC value with its limit (otherwise NaN results) + with np.errstate(divide='ignore', invalid='ignore'): + sf = c ** 2 * np.emath.sqrt((w / c) ** 2 - kgrid.ky ** 2 - kgrid.kz ** 2) / (2 * w) + sf[(w == 0) & (kgrid.ky == 0) & (kgrid.kz == 0)] = c / 2 + + # compute the FFT of the input data p(t, y, z) to yield p(w, ky, kz) and scale + p = sf * fftshift(fftn(ifftshift(p))) + + # exclude the inhomogeneous part of the wave + p[np.abs(w) < (c * np.sqrt(kgrid.ky ** 2 + kgrid.kz ** 2))] = 0 + + # compute the interpolation from p(w, ky, kz) to p(kx, ky, kz) + interp_func = RegularGridInterpolator( + (w[:, 0, 0], kgrid.ky[0, :, 0], kgrid.kz[0, 0, :]), + p, bounds_error=False, fill_value=0, method=interp) + query_points = np.stack((w_new, kgrid.ky, kgrid.kz), axis=-1) + p = interp_func(query_points) + + # compute the inverse FFT of p(kx, ky, kz) to yield p(x, y, z) + p = np.real(fftshift(ifftn(ifftshift(p)))) + + # remove the left part of the mirrored data which corresponds to the + # negative part of the mirrored time data + p = p[(Nt // 2):, ] + + # correct the scaling - the forward FFT is computed with a spacing of dt + # and the reverse requires a spacing of dz = dt*c, the reconstruction + # assumes that p0 is symmetrical about z, and only half the plane collects + # data (first approximation to correcting the limited view problem) + p = 2 * 2 * p / c + + # enforce positivity condition + if pos_cond: + logging.log(logging.INFO, 'applying positivity condition...') + p[p < 0] = 0 + + return p diff --git a/tests/test_executor.py b/tests/test_executor.py index c1826c4f..b8c1116e 100644 --- a/tests/test_executor.py +++ b/tests/test_executor.py @@ -21,6 +21,7 @@ def setUp(self): self.execution_options.binary_path = Path("/fake/path/to/binary") self.execution_options.system_string = "fake_system" self.execution_options.show_sim_log = False + self.simulation_options.pml_inside = True # Mock for stat result self.mock_stat_result = Mock() @@ -43,6 +44,28 @@ def setUp(self): self.mock_proc.communicate.return_value = ("stdout content", "stderr content") self.mock_popen.return_value.__enter__.return_value = self.mock_proc + # Mock return dictionary + N = np.array([20, 20, 20]) + pml = np.array([2, 2, 2]) + two_d_output = MagicMock() + two_d_output.ndim = 2 + three_d_output = MagicMock() + three_d_output.ndim = 3 + self.mock_dict_values = { + "Nx": N[0], + "Ny": N[1], + "Nz": N[2], + "pml_x_size": pml[0], + "pml_y_size": pml[1], + "pml_z_size": pml[2], + "axisymmetric_flag": np.array(False), + "p_final": two_d_output, + "p_max_all": three_d_output, + } + self.mock_dict = MagicMock() + self.mock_dict.__getitem__.side_effect = self.mock_dict_values.__getitem__ + self.mock_dict.__contains__.side_effect = self.mock_dict_values.__contains__ + def tearDown(self): # Stop patchers self.patcher_stat.stop() @@ -151,6 +174,51 @@ def test_parse_executable_output(self): self.assertIn("data", result) self.assertTrue(np.all(result["data"] == np.ones((10, 10)))) + def test_sensor_data_cropping_with_pml_outside(self): + """If pml is outside, fields like p_final and p_max_all should be cropped.""" + self.mock_proc.returncode = 0 + self.simulation_options.pml_inside = False + + # Instantiate the Executor + executor = Executor(self.execution_options, self.simulation_options) + + # Mock the parse_executable_output method + with patch.object(executor, "parse_executable_output", return_value=self.mock_dict): + sensor_data = executor.run_simulation("input.h5", "output.h5", ["options"]) + + # because pml is outside, the output should be cropped + two_d_output = sensor_data["p_final"] + two_d_output.__getitem__.assert_called_once_with((slice(2, 18), slice(2, 18))) + three_d_output = sensor_data["p_max_all"] + three_d_output.__getitem__.assert_called_once_with((slice(2, 18), slice(2, 18), slice(2, 18))) + + # check that the other fields are changed + for field in self.mock_dict_values.keys(): + if field not in ["p_final", "p_max_all"]: + self.assertEqual(sensor_data[field], self.mock_dict_values[field]) + + def test_sensor_data_cropping_with_pml_inside(self): + """If pml is inside, no field should be cropped.""" + self.mock_proc.returncode = 0 + + # Instantiate the Executor + executor = Executor(self.execution_options, self.simulation_options) + + # Mock the parse_executable_output method + with patch.object(executor, "parse_executable_output", return_value=self.mock_dict): + sensor_data = executor.run_simulation("input.h5", "output.h5", ["options"]) + + # because pml is inside, the output should not be cropped + two_d_output = sensor_data["p_final"] + two_d_output.__getitem__.assert_not_called() + three_d_output = sensor_data["p_max_all"] + three_d_output.__getitem__.assert_not_called() + + # check that the other fields are changed + for field in self.mock_dict_values.keys(): + if field not in ["p_final", "p_max_all"]: + self.assertEqual(sensor_data[field], self.mock_dict_values[field]) + def test_executor_file_not_found_on_non_darwin(self): # Configure the mock path object mock_binary_path = MagicMock(spec=Path) diff --git a/tests/test_simulation.py b/tests/test_simulation.py index bff0cde1..64d220f0 100644 --- a/tests/test_simulation.py +++ b/tests/test_simulation.py @@ -50,37 +50,64 @@ def setUp(self): # define simulation options self.simulation_options = SimulationOptions( - pml_inside=False, pml_size=self.PML_size, data_cast=self.DATA_CAST, - save_to_disk=True, data_recast=True, smooth_p0=False, + pml_inside=False, + pml_size=self.PML_size, + data_cast=self.DATA_CAST, + save_to_disk=True, + data_recast=True, + smooth_p0=False, ) # create the time array self.kgrid.makeTime(self.medium.sound_speed) def test_record_final_pressure(self): - self.sensor.record = ['p_final'] - k_sim = kWaveSimulation(kgrid=self.kgrid, source=self.source, sensor=self.sensor, - medium=self.medium, simulation_options=self.simulation_options) + self.sensor.record = ["p_final"] + k_sim = kWaveSimulation( + kgrid=self.kgrid, source=self.source, sensor=self.sensor, medium=self.medium, simulation_options=self.simulation_options + ) k_sim.input_checking("kspaceFirstOrder2D") recorder = k_sim.record.__dict__ for key, val in recorder.items(): - if key == 'p_final': + if key == "p_final": assert val - elif key.startswith(('p', 'u', 'I')): + elif key.startswith(("p", "u", "I")): + assert not val + + def test_time_reversal(self): + self.sensor.time_reversal_boundary_data = np.zeros((len(self.sensor.mask[0]), self.kgrid.Nt)) + + k_sim = kWaveSimulation( + kgrid=self.kgrid, source=self.source, sensor=self.sensor, medium=self.medium, simulation_options=self.simulation_options + ) + k_sim.input_checking("kspaceFirstOrder2D") + + # source and sensor are replaced when time-reversal is enabled + assert not (k_sim.source is self.source) + assert not (k_sim.sensor is self.sensor) + + assert k_sim.userarg_time_rev + + recorder = k_sim.record.__dict__ + for key, val in recorder.items(): + if key == "p_final": + assert val + elif key.startswith(("p", "u", "I")): assert not val def test_record_pressure(self): - self.sensor.record = ['p'] - k_sim = kWaveSimulation(kgrid=self.kgrid, source=self.source, sensor=self.sensor, - medium=self.medium, simulation_options=self.simulation_options) + self.sensor.record = ["p"] + k_sim = kWaveSimulation( + kgrid=self.kgrid, source=self.source, sensor=self.sensor, medium=self.medium, simulation_options=self.simulation_options + ) k_sim.input_checking("kspaceFirstOrder2D") recorder = k_sim.record.__dict__ for key, val in recorder.items(): - if key == 'p': + if key == "p": assert val - elif key.startswith(('p', 'u', 'I')): + elif key.startswith(("p", "u", "I")): assert not val