diff --git a/examples/README.md b/examples/README.md index 6927bfac..16fada4d 100644 --- a/examples/README.md +++ b/examples/README.md @@ -28,7 +28,9 @@ Every example has a short readme.md file which briefly describes the purpose of - [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)) -- [Focussed Detector In 2D Example](at_focused_annular_array_3D/) ([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)) +- [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)) + +- [Focussed Detector In 3D Example](sd_focussed_detector_3D/) ([original example](http://www.k-wave.org/documentation/example_sd_focussed_detector_3D.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_3D/sd_focussed_detector_3D.ipynb)) ## Contributing new examples diff --git a/examples/sd_focussed_detector_3D/README.md b/examples/sd_focussed_detector_3D/README.md new file mode 100644 index 00000000..0243ef37 --- /dev/null +++ b/examples/sd_focussed_detector_3D/README.md @@ -0,0 +1,7 @@ +# Focussed Detector In 3D 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/sd_focussed_detector_3D/sd_focussed_detector_3D.ipynb) + +This example shows how k-Wave can be used to model the output of a focussed bowl detector where the directionality arises from spatially averaging across the detector surface. + +To read more, visit the [original example page](http://www.k-wave.org/documentation/example_sd_focussed_detector_3D.php). diff --git a/examples/sd_focussed_detector_3D/sd_focussed_detector_3D.ipynb b/examples/sd_focussed_detector_3D/sd_focussed_detector_3D.ipynb new file mode 100644 index 00000000..b45df3a7 --- /dev/null +++ b/examples/sd_focussed_detector_3D/sd_focussed_detector_3D.ipynb @@ -0,0 +1,276 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%pip install git+https://github.com/waltsims/k-wave-python " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Focussed Detector In 3D Example\n", + "This example shows how k-Wave can be used to model the output of a focussed bowl detector where the directionality arises from spatially averaging across the detector surface." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "\n", + "import os\n", + "from copy import deepcopy\n", + "from tempfile import gettempdir\n", + "\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "\n", + "from kwave.data import Vector\n", + "from kwave.kgrid import kWaveGrid\n", + "from kwave.kmedium import kWaveMedium\n", + "from kwave.ksource import kSource\n", + "from kwave.ktransducer import kSensor\n", + "from kwave.options.simulation_execution_options import SimulationExecutionOptions\n", + "from kwave.options.simulation_options import SimulationOptions\n", + "from kwave.utils.mapgen import make_bowl\n", + "from kwave.utils.data import scale_SI\n", + "\n", + "from kwave.kspaceFirstOrder3D import kspaceFirstOrder3D\n", + "from kwave.utils.filters import filter_time_series" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# create the computational grid\n", + "grid_size = Vector([64, 64, 64]) # [grid points]\n", + "grid_spacing_single = 100e-3 / grid_size.x\n", + "grid_spacing = Vector([grid_spacing_single, grid_spacing_single, grid_spacing_single]) # [m]\n", + "kgrid = kWaveGrid(grid_size, grid_spacing)\n", + "\n", + "# define the properties of the propagation medium\n", + "medium = kWaveMedium(sound_speed=1500)\n", + "\n", + "# define the array of temporal points\n", + "_ = kgrid.makeTime(medium.sound_speed)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Define simulation parameters" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "input_filename = \"example_sd_focused_3d_input.h5\"\n", + "pathname = gettempdir()\n", + "input_file_full_path = os.path.join(pathname, input_filename)\n", + "simulation_options = SimulationOptions(\n", + " save_to_disk=True, \n", + " input_filename=input_filename, \n", + " data_path=pathname, \n", + " pml_size=10, \n", + " data_cast='single'\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Define a concave sensor" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# create a concave sensor\n", + "sphere_offset = 10\n", + "diameter = grid_size.x / 2 + 1\n", + "radius = grid_size.x / 2\n", + "bowl_pos = Vector([1 + sphere_offset, grid_size.y / 2, grid_size.z / 2])\n", + "focus_pos = grid_size / 2\n", + "\n", + "sensor_mask = make_bowl(grid_size, bowl_pos, radius, diameter, focus_pos)\n", + "sensor = kSensor(sensor_mask)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Run simulation with first source" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# define a time varying sinusoidal source\n", + "source = kSource()\n", + "\n", + "source_freq = 0.25e6 # [Hz]\n", + "source_mag = 1 # [Pa]\n", + "source.p = source_mag * np.sin(2 * np.pi * source_freq * kgrid.t_array)\n", + "source.p = filter_time_series(kgrid, medium, source.p)\n", + "\n", + "# place the first point source near the focus of the detector\n", + "source1 = np.zeros(grid_size)\n", + "source1[int(sphere_offset + radius), grid_size.y // 2 - 1, grid_size.z // 2 - 1] = 1\n", + "\n", + "# run the first simulation\n", + "source.p_mask = source1\n", + "\n", + "sensor_data1 = kspaceFirstOrder3D(\n", + " medium=medium,\n", + " kgrid=kgrid,\n", + " source=deepcopy(source),\n", + " sensor=deepcopy(sensor),\n", + " simulation_options=simulation_options,\n", + " execution_options=SimulationExecutionOptions(is_gpu_simulation=False),\n", + ")\n", + "\n", + "# average the data recorded at each grid point to simulate the measured signal from a single element focused detector\n", + "sensor_data1 = np.sum(sensor_data1['p'], axis=1)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Run simulation with second source" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# place the second point source off axis\n", + "source2 = np.zeros(grid_size)\n", + "source2[int(1 + sphere_offset + radius), grid_size.y // 2 + 5, grid_size.z // 2 + 5] = 1\n", + "\n", + "\n", + "# run the second simulation\n", + "source.p_mask = source2\n", + "sensor_data2 = kspaceFirstOrder3D(\n", + " medium=medium,\n", + " kgrid=kgrid,\n", + " source=source,\n", + " sensor=sensor,\n", + " simulation_options=simulation_options,\n", + " execution_options=SimulationExecutionOptions(is_gpu_simulation=False),\n", + ")\n", + "\n", + "# average the data recorded at each grid point to simulate the measured signal from a single element focused detector\n", + "sensor_data2 = np.sum(sensor_data2['p'], axis=1)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Visualize detector and sources" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Combine arrays as in MATLAB: sensor.mask + source1 + source2\n", + "combined_array = sensor_mask + source1 + source2\n", + "\n", + "# Find the indices of non-zero elements\n", + "x, y, z = np.nonzero(combined_array)\n", + "\n", + "# Enable interactive mode\n", + "plt.ion()\n", + "\n", + "# Create an interactive 3D plot with matplotlib\n", + "fig = plt.figure(figsize=(8, 6))\n", + "ax = fig.add_subplot(111, projection='3d')\n", + "sc = ax.scatter(x, y, z, c='blue', marker='s', s=100, depthshade=True) # 's' for square-like markers\n", + "\n", + "# Set the view angle to mimic MATLAB's `view([130, 40])`\n", + "ax.view_init(elev=40, azim=130)\n", + "\n", + "# Customize the axes\n", + "ax.set_xlabel('X')\n", + "ax.set_ylabel('Y')\n", + "ax.set_zlabel('Z')\n", + "ax.set_title('Interactive 3D Voxel Plot')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Visualize recorded data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "t_sc, t_scale, t_prefix, _ = scale_SI(kgrid.t_array[-1])\n", + "t_array = kgrid.t_array.squeeze() * t_scale\n", + "\n", + "plt.figure()\n", + "plt.plot(t_array, sensor_data1)\n", + "plt.plot(t_array, sensor_data2, 'r')\n", + "\n", + "plt.xlabel('Time [' + t_prefix + 's]')\n", + "plt.ylabel('Average Pressure Measured By Focussed Detector [Pa]')\n", + "plt.legend(['Source on axis', 'Source off axis'])\n", + "\n", + "plt.show()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "env_kwave", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/examples/sd_focussed_detector_3D/sd_focussed_detector_3D.py b/examples/sd_focussed_detector_3D/sd_focussed_detector_3D.py new file mode 100644 index 00000000..b6ee4df2 --- /dev/null +++ b/examples/sd_focussed_detector_3D/sd_focussed_detector_3D.py @@ -0,0 +1,134 @@ +# Focussed Detector In 3D Example +# This example shows how k-Wave can be used to model the output of a focussed bowl detector where the directionality arises from spatially averaging across the detector surface. + +import os +from copy import deepcopy +from tempfile import gettempdir + +import numpy as np +import matplotlib.pyplot as plt + +from kwave.data import Vector +from kwave.kgrid import kWaveGrid +from kwave.kmedium import kWaveMedium +from kwave.ksource import kSource +from kwave.ktransducer import kSensor +from kwave.options.simulation_execution_options import SimulationExecutionOptions +from kwave.options.simulation_options import SimulationOptions +from kwave.utils.mapgen import make_bowl +from kwave.utils.data import scale_SI + +from kwave.kspaceFirstOrder3D import kspaceFirstOrder3D +from kwave.utils.filters import filter_time_series + + +# create the computational grid +grid_size = Vector([64, 64, 64]) # [grid points] +grid_spacing_single = 100e-3 / grid_size.x +grid_spacing = Vector([grid_spacing_single, grid_spacing_single, grid_spacing_single]) # [m] +kgrid = kWaveGrid(grid_size, grid_spacing) + +# define the properties of the propagation medium +medium = kWaveMedium(sound_speed=1500) + +# define the array of temporal points +_ = kgrid.makeTime(medium.sound_speed) + +input_filename = "example_sd_focused_3d_input.h5" +pathname = gettempdir() +input_file_full_path = os.path.join(pathname, input_filename) +simulation_options = SimulationOptions( + save_to_disk=True, input_filename=input_filename, data_path=pathname, pml_size=10, data_cast="single" +) + +# create a concave sensor +sphere_offset = 10 +diameter = grid_size.x / 2 + 1 +radius = grid_size.x / 2 +bowl_pos = Vector([1 + sphere_offset, grid_size.y / 2, grid_size.z / 2]) +focus_pos = grid_size / 2 + +sensor_mask = make_bowl(grid_size, bowl_pos, radius, diameter, focus_pos) +sensor = kSensor(sensor_mask) + +# define a time varying sinusoidal source +source = kSource() + +source_freq = 0.25e6 # [Hz] +source_mag = 1 # [Pa] +source.p = source_mag * np.sin(2 * np.pi * source_freq * kgrid.t_array) +source.p = filter_time_series(kgrid, medium, source.p) + +# place the first point source near the focus of the detector +source1 = np.zeros(grid_size) +source1[int(sphere_offset + radius), grid_size.y // 2 - 1, grid_size.z // 2 - 1] = 1 + +# run the first simulation +source.p_mask = source1 + +sensor_data1 = kspaceFirstOrder3D( + medium=medium, + kgrid=kgrid, + source=deepcopy(source), + sensor=deepcopy(sensor), + simulation_options=simulation_options, + execution_options=SimulationExecutionOptions(is_gpu_simulation=False), +) + +# average the data recorded at each grid point to simulate the measured signal from a single element focused detector +sensor_data1 = np.sum(sensor_data1["p"], axis=1) + +# place the second point source off axis +source2 = np.zeros(grid_size) +source2[int(1 + sphere_offset + radius), grid_size.y // 2 + 5, grid_size.z // 2 + 5] = 1 + + +# run the second simulation +source.p_mask = source2 +sensor_data2 = kspaceFirstOrder3D( + medium=medium, + kgrid=kgrid, + source=source, + sensor=sensor, + simulation_options=simulation_options, + execution_options=SimulationExecutionOptions(is_gpu_simulation=False), +) + +# average the data recorded at each grid point to simulate the measured signal from a single element focused detector +sensor_data2 = np.sum(sensor_data2["p"], axis=1) + +# Combine arrays as in MATLAB: sensor.mask + source1 + source2 +combined_array = sensor_mask + source1 + source2 + +# Find the indices of non-zero elements +x, y, z = np.nonzero(combined_array) + +# Enable interactive mode +plt.ion() + +# Create an interactive 3D plot with matplotlib +fig = plt.figure(figsize=(8, 6)) +ax = fig.add_subplot(111, projection="3d") +sc = ax.scatter(x, y, z, c="blue", marker="s", s=100, depthshade=True) # 's' for square-like markers + +# Set the view angle to mimic MATLAB's `view([130, 40])` +ax.view_init(elev=40, azim=130) + +# Customize the axes +ax.set_xlabel("X") +ax.set_ylabel("Y") +ax.set_zlabel("Z") +ax.set_title("Interactive 3D Voxel Plot") + +t_sc, t_scale, t_prefix, _ = scale_SI(kgrid.t_array[-1]) +t_array = kgrid.t_array.squeeze() * t_scale + +plt.figure() +plt.plot(t_array, sensor_data1) +plt.plot(t_array, sensor_data2, "r") + +plt.xlabel("Time [" + t_prefix + "s]") +plt.ylabel("Average Pressure Measured By Focussed Detector [Pa]") +plt.legend(["Source on axis", "Source off axis"]) + +plt.show()