#441 Add us beams example (#445)

* add us_defining_transducer

* add to list of examples

* remove unrequired imports

* add beam pattern example upto broken spect call

* burst freq to 500Hz

* working beam pattern example

* add README.md

* fixes to us_beam_patterns

* fix linting issues raised by ruff

* update us_beam_pattern to use spect

* remove useless expression

* update notebook

* update color map and default values

---------

Co-authored-by: Walter Simson <[email protected]>
Co-authored-by: Walter Simson <[email protected]>

.gitignore

@@ -4,6 +4,7 @@
 
 # Virtual environments
 .venv/
+venv/
 
 # Compiled files
 *.pyc

examples/README.md

@@ -14,6 +14,7 @@ Every example has a short readme.md file which briefly describes the purpose of
 - [Controlling the PML](na_controlling_the_pml/)
 ([original example](http://www.k-wave.org/documentation/example_na_controlling_the_pml.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/na_controlling_the_pml/na_controlling_the_pml.ipynb))
 - [Defining An Ultrasound Transducer Example](us_defining_transducer) ([original example](http://www.k-wave.org/documentation/example_us_defining_transducer.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/us_defining_transducer/us_defining_transducer.ipynb))
+- [Simulating Ultrasound Beam Patterns](us_beam_patterns/)([original example](http://www.k-wave.org/documentation/example_us_beam_patterns), [![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/us_beam_patterns/us_beam_patterns.ipynb))
 - [Linear transducer B-mode](us_bmode_linear_transducer/) ([original example](http://www.k-wave.org/documentation/example_us_bmode_linear_transducer.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/us_bmode_linear_transducer/us_bmode_linear_transducer.ipynb))
 - [Phased array B-mode](us_bmode_phased_array/)
 ([original example](http://www.k-wave.org/documentation/example_us_bmode_phased_array.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/us_bmode_phased_array/us_bmode_phased_array.ipynb))
@@ -38,10 +39,10 @@ When adding a new example notebook, follow these steps:
 3. Fork and clone the repository and create a branch for your new example.
 2. Create an example sub-directory using the name from the hyperlink of the original k-wave example if it exists (e.g. for http://www.k-wave.org/documentation/example_ivp_loading_external_image.php name the directory "ivp_loading_external_image).
 3. Add your example notebook to your example directory.
-4. Create a Python script that mirrors your example notebook in the same directory.
+4. Create a Python script that mirrors your example notebook in the same directory. Using `nbconvert` is an easy way to convert to a Python script using ` jupyter nbconvert --to python --RegexRemovePreprocessor.patterns="^%" <notebook_path>`
 5. Add a README.md file to your example directory briefly describing the concept or principle the example is meant to display and linking to the origonal k-wave example page if it exists.
-6. Include a link in the readme.md in the examples directory to a colab notebook for your example.
-7. Add a your example to the list on this readme.md and add a colab badge [using html](https://openincolab.com/) OR copy the pure markdown version above.
+6. Include a link in the `README.md` in the examples directory to a colab notebook for your example.
+7. Add a your example to the list on this `README.md` and add a colab badge [using html](https://openincolab.com/) OR copy the pure markdown version above.
 8. Open a pull request that [closes the open issue](https://docs.github.com/en/issues/tracking-your-work-with-issues/linking-a-pull-request-to-an-issue) from your forked example branch and name pull request "[Example] \<name of your example\>".
 
 Thanks for contributing to k-wave-python!

examples/us_beam_patterns/README.md

@@ -0,0 +1,7 @@
+# Simulating Ultrasound Beam Patterns 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/us_beam_patterns/us_beam_patterns.ipynb)
+
+This example shows how the nonlinear beam pattern from an ultrasound transducer can be modelled. It builds on the Defining An Ultrasound Transducer and Simulating Transducer Field Patterns examples.
+
+To read more, visit the [original example page](http://www.k-wave.org/documentation/example_us_beam_patterns.php).

examples/us_beam_patterns/us_beam_patterns.ipynb

@@ -0,0 +1,346 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%%capture\n",
+    "# !pip install k-wave-python"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
+    "\n",
+    "from kwave.data import Vector\n",
+    "from kwave.kgrid import kWaveGrid\n",
+    "from kwave.kspaceFirstOrder3D import kspaceFirstOrder3D\n",
+    "from kwave.ksensor import kSensor\n",
+    "from kwave.ktransducer import kWaveTransducerSimple, NotATransducer\n",
+    "from kwave.kWaveSimulation import SimulationOptions\n",
+    "from kwave.kmedium import kWaveMedium\n",
+    "from kwave.utils.filters import spect\n",
+    "\n",
+    "from kwave.options.simulation_execution_options import SimulationExecutionOptions\n",
+    "from kwave.utils.dotdictionary import dotdict\n",
+    "from kwave.utils.math import find_closest\n",
+    "from kwave.utils.signals import tone_burst"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# simulation settings\n",
+    "DATA_CAST = \"single\"\n",
+    "# set to 'xy' or 'xz' to generate the beam pattern in different planes\n",
+    "MASK_PLANE = \"xy\"\n",
+    "# set to true to compute the rms or peak beam patterns, set to false to compute the harmonic beam patterns\n",
+    "USE_STATISTICS = True\n",
+    "\n",
+    "# define the grid\n",
+    "PML_X_SIZE = 20\n",
+    "PML_Y_SIZE = 10\n",
+    "PML_Z_SIZE = 10\n",
+    "Nx = 128 - 2 * PML_X_SIZE\n",
+    "Ny = 64 - 2 * PML_Y_SIZE\n",
+    "Nz = 64 - 2 * PML_Z_SIZE\n",
+    "x = 40e-3\n",
+    "dx = x / Nx\n",
+    "dy = dx\n",
+    "dz = dx\n",
+    "\n",
+    "kgrid = kWaveGrid([Nx, Ny, Nz], [dx, dy, dz])\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# define the medium\n",
+    "medium = kWaveMedium(sound_speed=1540, density=1000, alpha_coeff=0.75, alpha_power=1.5, BonA=6)\n",
+    "\n",
+    "# create the time array - using a time == time to travel the hypot of the grid\n",
+    "t_end = 45e-6  # alternatively, use np.sqrt(kgrid.x_size ** 2 + kgrid.y_size ** 2) / medium.sound_speed\n",
+    "\n",
+    "kgrid.makeTime(medium.sound_speed, t_end=t_end)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# define the input signal\n",
+    "source_strength = 1e6\n",
+    "tone_burst_freq = 0.5e6\n",
+    "tone_burst_cycles = 5\n",
+    "input_signal = tone_burst(1 / kgrid.dt, tone_burst_freq, tone_burst_cycles)\n",
+    "input_signal = (source_strength / (medium.sound_speed * medium.density)) * input_signal\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# define the transducer\n",
+    "transducer = dotdict()\n",
+    "transducer.number_elements = 32\n",
+    "transducer.element_width = 1\n",
+    "transducer.element_length = 12\n",
+    "transducer.element_spacing = 0\n",
+    "transducer.radius = np.inf\n",
+    "\n",
+    "# calculate the width of the transducer in grid points\n",
+    "transducer_width = transducer.number_elements * transducer.element_width + (transducer.number_elements - 1) * transducer.element_spacing\n",
+    "\n",
+    "# use this to position the transducer in the middle of the computational grid\n",
+    "transducer.position = np.round([1, Ny / 2 - transducer_width / 2, Nz / 2 - transducer.element_length / 2])\n",
+    "transducer = kWaveTransducerSimple(kgrid, **transducer)\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "not_transducer = dotdict()\n",
+    "not_transducer.sound_speed = medium.sound_speed  # sound speed [m/s]\n",
+    "not_transducer.focus_distance = 20e-3  # focus distance [m]\n",
+    "not_transducer.elevation_focus_distance = 19e-3  # focus distance in the elevation plane [m]\n",
+    "not_transducer.steering_angle = 0  # steering angle [degrees]\n",
+    "not_transducer.transmit_apodization = \"Rectangular\"\n",
+    "not_transducer.receive_apodization = \"Rectangular\"\n",
+    "not_transducer.active_elements = np.ones((transducer.number_elements, 1))\n",
+    "not_transducer.input_signal = input_signal\n",
+    "\n",
+    "not_transducer = NotATransducer(transducer, kgrid, **not_transducer)\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sensor_mask = np.zeros((Nx, Ny, Nz))\n",
+    "\n",
+    "if MASK_PLANE == \"xy\":\n",
+    "    sensor_mask[:, :, Nz // 2] = 1\n",
+    "    # store y axis properties\n",
+    "    Nj = Ny\n",
+    "    j_vec = kgrid.y_vec\n",
+    "    j_label = \"y\"\n",
+    "elif MASK_PLANE == \"xz\":\n",
+    "    sensor_mask[:, Ny // 2, :] = 1\n",
+    "    # store z axis properties\n",
+    "    Nj = Nz\n",
+    "    j_vec = kgrid.z_vec\n",
+    "    j_label = \"z\"\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sensor = kSensor(sensor_mask)\n",
+    "if USE_STATISTICS:\n",
+    "    sensor.record = [\"p_rms\", \"p_max\"]\n",
+    "\n",
+    "simulation_options = SimulationOptions(\n",
+    "    pml_inside=False,\n",
+    "    pml_size=Vector([PML_X_SIZE, PML_Y_SIZE, PML_Z_SIZE]),\n",
+    "    data_cast=DATA_CAST,\n",
+    "    save_to_disk=True,\n",
+    ")\n",
+    "\n",
+    "if not USE_STATISTICS:\n",
+    "    simulation_options.stream_to_disk = \"harmonic_data.h5\"\n",
+    "\n",
+    "execution_options = SimulationExecutionOptions(is_gpu_simulation=True)\n",
+    "\n",
+    "sensor_data = kspaceFirstOrder3D(\n",
+    "    medium=medium,\n",
+    "    kgrid=kgrid,\n",
+    "    source=not_transducer,\n",
+    "    sensor=sensor,\n",
+    "    simulation_options=simulation_options,\n",
+    "    execution_options=execution_options,\n",
+    ")\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "if USE_STATISTICS:\n",
+    "    fig, axes = plt.subplots(1, 2)\n",
+    "    fig.set_figwidth(8)\n",
+    "    fig.set_figheight(6)\n",
+    "\n",
+    "    for ind, measures in enumerate(sensor.record):\n",
+    "        im1 = axes[ind].imshow(\n",
+    "            sensor_data[measures].reshape([Nj, Nx]).T * 1e-6,\n",
+    "            extent=[\n",
+    "                min(j_vec * 1e3),\n",
+    "                max(j_vec * 1e3),\n",
+    "                min((kgrid.x_vec - min(kgrid.x_vec)) * 1e3),\n",
+    "                max((kgrid.x_vec - min(kgrid.x_vec)) * 1e3),\n",
+    "            ],\n",
+    "            aspect=\"auto\",\n",
+    "            cmap=\"jet\",\n",
+    "        )\n",
+    "        axes[ind].set_xlabel(\"y-position (mm)\")\n",
+    "        axes[ind].set_ylabel(\"x-position (mm)\")\n",
+    "        title_text = f\"Total Beam Pattern Using {str.upper(measures.split('_')[1])} \\n of Recorded Pressure\"\n",
+    "        axes[ind].set_title(title_text)\n",
+    "        axes[ind].set_yticks(axes[ind].get_yticks().tolist())\n",
+    "        axes[ind].set_yticklabels(axes[ind].get_yticklabels()[::-1])\n",
+    "        axes[ind].set_xlim([-10, 10])\n",
+    "        fig.colorbar(im1, label=\"Pressure [MPa]\", orientation=\"vertical\")\n",
+    "\n",
+    "    plt.tight_layout()\n",
+    "    plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "if not USE_STATISTICS:\n",
+    "    sensor_data_array = np.reshape(sensor_data[\"p\"], [kgrid.Nt, kgrid.Ny, kgrid.Nx]).transpose(2, 1, 0)\n",
+    "    # compute the amplitude spectrum\n",
+    "    [freq, amp_spect, _] = spect(sensor_data_array, 1 / kgrid.dt, dim=2)\n",
+    "\n",
+    "    # compute the index at which the source frequency and its harmonics occur\n",
+    "    [f1_value, f1_index] = find_closest(freq, tone_burst_freq)\n",
+    "    [f2_value, f2_index] = find_closest(freq, 2 * tone_burst_freq)\n",
+    "\n",
+    "    # extract the amplitude at the source frequency and store\n",
+    "    beam_pattern_f1 = np.squeeze(amp_spect[:, :, f1_index])\n",
+    "\n",
+    "    # extract the amplitude at the second harmonic and store\n",
+    "    beam_pattern_f2 = np.squeeze(amp_spect[:, :, f2_index])\n",
+    "\n",
+    "    # extract the integral of the total amplitude spectrum\n",
+    "    beam_pattern_total = np.squeeze(np.sum(amp_spect, axis=2))\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "if not USE_STATISTICS:\n",
+    "    fig, axes = plt.subplots(1, 3)\n",
+    "    fig.set_figwidth(8)\n",
+    "    fig.set_figheight(6)\n",
+    "\n",
+    "    for ind, measure in enumerate([beam_pattern_f1, beam_pattern_f2, beam_pattern_total]):\n",
+    "        im1 = axes[ind].imshow(\n",
+    "            np.squeeze(measure) * 1e-6,\n",
+    "            extent=[\n",
+    "                min(j_vec * 1e3),\n",
+    "                max(j_vec * 1e3),\n",
+    "                min((kgrid.x_vec - min(kgrid.x_vec)) * 1e3),\n",
+    "                max((kgrid.x_vec - min(kgrid.x_vec)) * 1e3),\n",
+    "            ],\n",
+    "            aspect=\"auto\",\n",
+    "            cmap=\"jet\",\n",
+    "        )\n",
+    "        axes[ind].set_xlabel(\"y-position (mm)\")\n",
+    "        axes[ind].set_ylabel(\"x-position (mm)\")\n",
+    "\n",
+    "        axes[ind].set_yticks(axes[ind].get_yticks().tolist())\n",
+    "        axes[ind].set_yticklabels(axes[ind].get_yticklabels()[::-1])\n",
+    "\n",
+    "    axes[0].set_title(\"Beam Pattern \\nAt Source Fundamental\")\n",
+    "    axes[1].set_title(\"Beam Pattern \\nAt Second Harmonic\")\n",
+    "    axes[2].set_title(\"Total Beam Pattern\\nUsing Integral Of Recorded Pressure\")\n",
+    "\n",
+    "    plt.tight_layout()\n",
+    "    plt.show()\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "if not USE_STATISTICS:\n",
+    "    # Compute the directivity at each of the harmonics\n",
+    "    directivity_f1 = beam_pattern_f1[round(not_transducer.focus_distance / dx), :]\n",
+    "    directivity_f2 = beam_pattern_f2[round(not_transducer.focus_distance / dx), :]\n",
+    "\n",
+    "    # Normalize the directivity\n",
+    "    directivity_f1 /= np.max(directivity_f1)\n",
+    "    directivity_f2 /= np.max(directivity_f2)\n",
+    "\n",
+    "    # Compute relative angles from the transducer\n",
+    "    if MASK_PLANE == \"xy\":\n",
+    "        horz_axis = ((np.arange(Ny) + 1) - Ny / 2) * dy\n",
+    "    else:\n",
+    "        horz_axis = ((np.arange(Nz) + 1) - Nz / 2) * dz\n",
+    "\n",
+    "    angles = 180 * np.arctan2(horz_axis, not_transducer.focus_distance) / np.pi\n",
+    "\n",
+    "    # Plot the directivity\n",
+    "    plt.figure()\n",
+    "    plt.plot(angles, directivity_f1, \"k-\", label=\"Fundamental\")\n",
+    "    plt.plot(angles, directivity_f2, \"k--\", label=\"Second Harmonic\")\n",
+    "    plt.axis(\"tight\")\n",
+    "    plt.xlabel(\"Angle [deg]\")\n",
+    "    plt.ylabel(\"Normalized Amplitude\")\n",
+    "    plt.legend(loc=2)\n",
+    "    plt.xlim([-25, 25])\n",
+    "    plt.show()"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}