From 251d4a1e97e6320c633acbec0fa062181383ca2f Mon Sep 17 00:00:00 2001 From: David Sinden Date: Mon, 24 Nov 2025 09:58:33 +0100 Subject: [PATCH 01/11] Add fullwave_simple_plane_wave simulation example This script replicates the fullwave25 simple_plane_wave example using k-wave, including methods for mapping coordinates and generating a Gaussian-modulated sinusoidal signal. It sets up the simulation parameters, defines the acoustic medium properties, initializes the pressure source, runs the simulation, and visualizes the results. --- .../fullwave_simple_plane_wave.py | 307 ++++++++++++++++++ 1 file changed, 307 insertions(+) create mode 100644 examples/fullwave_simple_plane_wave/fullwave_simple_plane_wave.py diff --git a/examples/fullwave_simple_plane_wave/fullwave_simple_plane_wave.py b/examples/fullwave_simple_plane_wave/fullwave_simple_plane_wave.py new file mode 100644 index 00000000..97baeea0 --- /dev/null +++ b/examples/fullwave_simple_plane_wave/fullwave_simple_plane_wave.py @@ -0,0 +1,307 @@ +""" +k-wave python example which replicates the fullwave25 simple_plane_wave example, directly +copying methods (maps-to_coords and gaussian_modulated_sinusoidal_signal) from the package. +""" + +from copy import deepcopy + +from matplotlib import animation +import matplotlib.pyplot as plt +import numpy as np + +from mpl_toolkits.axes_grid1 import make_axes_locatable +from numpy.typing import NDArray +from tqdm import tqdm + +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.utils.colormap import get_color_map +from kwave.kspaceFirstOrder2D import kspaceFirstOrder2D + +from kwave.options.simulation_execution_options import SimulationExecutionOptions +from kwave.options.simulation_options import SimulationOptions + + + +def map_to_coords( + map_data: NDArray[np.float64 | np.int64 | np.bool], + *, + export_as_xyz: bool = False, +) -> NDArray[np.int64]: + """Map the mask map to coordinates. + + Returns: + NDArray[np.int64]: An array of coordinates corresponding to non-zero elements in the mask. + + """ + is_3d = map_data.ndim == 3 + # indices = np.where(map_data.T != 0) + indices = np.where(map_data != 0) + if is_3d: + # out = np.array([indices[2], indices[1], indices[0]]).T + # out = np.array([indices[2], indices[1], indices[0]]).T + out = np.array([*indices]).T + + if export_as_xyz: + out = np.stack([out[:, 2], out[:, 1], out[:, 0]], axis=1) + else: + out = np.array([*indices]).T + if export_as_xyz: + out = np.stack([out[:, 1], out[:, 0]], axis=1) + return out + + +def gaussian_modulated_sinusoidal_signal( + nt: int, + duration: float, + ncycles: int, + drop_off: int, + f0: float, + p0: float, + delay_sec: float = 0.0, + i_layer: int | None = None, + dt_for_layer_delay: float | None = None, + cfl_for_layer_delay: float | None = None, +) -> NDArray[np.float64]: + """Generate a pulse signal based on input parameters. + + Parameters + ---------- + nt: int + Number of time samples of the simulation. + duration: float + Total duration of the simulation. + ncycles: int + Number of cycles in the pulse. + drop_off: int + Controls the pulse decay. + f0: float + Frequency of the pulse. + p0: float + Amplitude scaling factor. + delay_sec: float + Delay in seconds. Default is 0.0. + i_layer: int + Index of the layer where the source is located. Default is None. + This variable is used to shift the pulse signal in time + so that the signal is emmitted within the transducer layer correctly. + dt_for_layer_delay: float + Time step of the simulation. Default is None. + This variable is used to shift the pulse signal in time + so that the signal is emmitted within the transducer layer correctly. + cfl_for_layer_delay: float + Courant-Friedrichs-Lewy number. Default is None. + This variable is used to shift the pulse signal in time + so that the signal is emmitted within the transducer layer correctly. + + Returns + ------- + NDArray[np.float64]: The generated pulse signal. + + """ + + t = (np.arange(0, nt)) / nt * duration - ncycles / f0 + t = t - delay_sec + + if i_layer: + assert dt_for_layer_delay, "dt must be provided if i_layer is provided" + assert cfl_for_layer_delay, "cfl must be provided if i_layer is provided" + t = t - (dt_for_layer_delay / cfl_for_layer_delay) * i_layer + + omega0 = 2 * np.pi * f0 + return ( + np.multiply( + np.exp( + -((1.05 * t * omega0 / (ncycles * np.pi)) ** (2 * drop_off)), + ), + np.sin(t * omega0), + ) + * p0 + ) + + + +domain_size = (3e-2, 2e-2) # meters + +f0 = 3.0e6 +c0 = 1540.0 + +duration = domain_size[0] / c0 * 2 + +cfl = 0.2 +ppw = 12 + +wavelength = c0 / f0 + +dx = wavelength / ppw +dy = dx +dt = cfl * dx / c0 + +Nx = int(np.round(domain_size[0] / dx)) +Ny = int(np.round(domain_size[1] / dy)) + +grid = kWaveGrid(Vector([Nx, Ny]), Vector([dx, dy])) + +duration = domain_size[0] / c0 * 2 +Nt = int(np.round(duration / dt)) +grid.setTime(Nt, dt) + +# +# --- define the acoustic medium properties --- +# +# Define the base 2D medium arrays +sound_speed_map = 1540 * np.ones((grid.Nx, grid.Ny)) # m/s +density_map = 1000 * np.ones((grid.Nx, grid.Ny)) # kg/m^3 +alpha_coeff_map = 0.5 * np.ones((grid.Nx, grid.Ny)) # dB/(MHz^y cm) + +# embed an object with different properties in the center of the medium +obj_x_start = grid.Nx // 3 +obj_x_end = 2 * grid.Nx // 3 +obj_y_start = grid.Ny // 3 +obj_y_end = 2 * grid.Ny // 3 + +sound_speed_map[obj_x_start:obj_x_end, obj_y_start:obj_y_end] = 1600 # m/s +density_map[obj_x_start:obj_x_end, obj_y_start:obj_y_end] = 1100 # kg/m^3 +alpha_coeff_map[obj_x_start:obj_x_end, obj_y_start:obj_y_end] = 0.75 # dB/(MHz^y cm) + +# setup the Medium instance +medium = kWaveMedium(sound_speed=sound_speed_map, + density=density_map, + alpha_coeff=alpha_coeff_map, + alpha_power=1.1) + +# initialize the pressure source mask +p_mask = np.zeros((grid.Nx, grid.Ny), dtype=bool) + +# set the source location at the top rows of the grid with specified thickness +element_thickness_px = 3 +p_mask[0:element_thickness_px, :] = True + +# define the pressure source [n_sources, nt] +p0 = np.zeros((p_mask.sum(), grid.Nt)) + +p_coordinates = map_to_coords(p_mask) + +for i_thickness in range(element_thickness_px): + # create a gaussian-modulated sinusoidal pulse as the source signal with layer delay + p0_vec = gaussian_modulated_sinusoidal_signal(nt=grid.Nt, # number of time steps + f0=f0, # center frequency [Hz] + duration=duration, # duration [s] + ncycles=2, # number of cycles + drop_off=2, # drop off factor + p0=1e5, # maximum amplitude [Pa] + i_layer=i_thickness, + dt_for_layer_delay=grid.dt, + cfl_for_layer_delay=cfl) + + # assign the source signal to the corresponding layer + n_y = p_coordinates.shape[0] // element_thickness_px + p0[n_y * i_thickness : n_y * (i_thickness + 1), :] = p0_vec.copy() + +# create the kSource instance +source = kSource() +source.p_mask = p_mask +source.p = p0 + +# setup the Sensor instance +sensor_mask = np.ones((grid.Nx, grid.Ny), dtype=bool) +sensor = kSensor() +sensor.mask = sensor_mask +sensor.record = ["p"] + + +# -------------------- +# SIMULATION +# -------------------- + +simulation_options = SimulationOptions(pml_auto=True, data_recast=True, save_to_disk=True, save_to_disk_exit=False, pml_inside=False) + +execution_options = SimulationExecutionOptions(is_gpu_simulation=True, delete_data=False, verbose_level=2) + +sensor_data = kspaceFirstOrder2D(kgrid=deepcopy(grid), + source=deepcopy(source), + sensor=deepcopy(sensor), + medium=deepcopy(medium), + simulation_options=simulation_options, + execution_options=execution_options) + +# -------------------- +# VISUALIZATION +# -------------------- + +propagation_map = np.reshape(sensor_data["p"], (grid.Nt, grid.Nx, grid.Ny), order='F') # Store the recorded pressure data + +p_max_plot = np.abs(propagation_map).max().item() / 4 + +time_step = propagation_map.shape[0] // 3 + + + +cmap = get_color_map() + +fig, ax = plt.subplots(1, 1) +im = ax.imshow(propagation_map[time_step, :, :], + extent=[grid.y_vec.min() * 1e3, grid.y_vec.max() * 1e3, grid.x_vec.min() * 1e3, grid.x_vec.max() * 1e3], + vmin=-p_max_plot, + vmax=p_max_plot, + cmap=cmap) +title = "Snapshot" +ax.set_title(title) +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, figsize=(4,6)) + +num_plot_image: int = 50 +skip_every_n_frame =int(grid.Nt / num_plot_image) + +c_map = medium.sound_speed +rho_map = medium.density + +start = 0 +end = None +z_map = c_map * rho_map +z_map = (z_map - np.min(z_map)) / (np.max(z_map) - np.min(z_map) + 1e-9) + +z_map_offset = p_max_plot * 0.8 + +animation_list = [] + +# propagation_map = propagation_map.transpose(2, 0, 1) +for i, p_map_i in tqdm( + enumerate(propagation_map[::skip_every_n_frame, start:end, start:end]), + total=len(propagation_map[::skip_every_n_frame, start:end, start:end]), + desc="plotting animation"): + + processed_p_map = p_map_i + z_map_offset * (z_map) + + image2 = ax.imshow( + processed_p_map, + vmin=-p_max_plot, + vmax=p_max_plot, + interpolation="nearest" + ) + # set text to show the current time step + text = ax.text(0.5, 1.05, f"t = {i * skip_every_n_frame} / {propagation_map.shape[0]}", fontsize=4, + ha="center", + animated=True, + transform=ax.transAxes) + animation_list.append([image2, text]) + +animation_data = animation.ArtistAnimation( + fig, + animation_list, + interval=150, + blit=True, + repeat_delay=500, +) +animation_data.save("fullwave_plane_wave.mp4", writer="ffmpeg", dpi=300) +plt.close("all") From 945bf9bc741a221de1cb2ca16974a3466e814dcb Mon Sep 17 00:00:00 2001 From: David Sinden Date: Mon, 24 Nov 2025 10:05:04 +0100 Subject: [PATCH 02/11] Add README for Simple Plane Wave example Added README for Simple Plane Wave example showcasing compatibility with fullwave25. --- examples/fullwave_simple_plane_wave/README.md | 6 ++++++ 1 file changed, 6 insertions(+) create mode 100644 examples/fullwave_simple_plane_wave/README.md diff --git a/examples/fullwave_simple_plane_wave/README.md b/examples/fullwave_simple_plane_wave/README.md new file mode 100644 index 00000000..923baa5c --- /dev/null +++ b/examples/fullwave_simple_plane_wave/README.md @@ -0,0 +1,6 @@ +# Simple Plane Wave + + [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/waltsims/k-wave-python/blob/HEAD/examples/fullwave_simple_plane_wave.ipynb) + +This example showcases the compatibility between [fullwave25](https://github.com/pinton-lab/fullwave25) and k-wave-python. It directly copies two methods over. + From ef8f7f855ad270790c903f0a8e0ec2107040d773 Mon Sep 17 00:00:00 2001 From: David Sinden Date: Mon, 24 Nov 2025 10:22:20 +0100 Subject: [PATCH 03/11] Add ipynb upload notebook --- .../fullwave_simple_plane_wave.ipynb | 471 ++++++++++++++++++ 1 file changed, 471 insertions(+) create mode 100644 examples/fullwave_simple_plane_wave/fullwave_simple_plane_wave.ipynb diff --git a/examples/fullwave_simple_plane_wave/fullwave_simple_plane_wave.ipynb b/examples/fullwave_simple_plane_wave/fullwave_simple_plane_wave.ipynb new file mode 100644 index 00000000..4aa32997 --- /dev/null +++ b/examples/fullwave_simple_plane_wave/fullwave_simple_plane_wave.ipynb @@ -0,0 +1,471 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "5884fede-6f1f-4cc6-9ae9-d19e03951883", + "metadata": {}, + "source": [ + "# Fullwave25 Plane wave example\n", + "\n", + "k-wave python example which replicates the fullwave25 simple_plane_wave example.\n", + "\n", + "Begin with installation and imports" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f8c5e68e-7e51-404a-838d-f219b38b8fb6", + "metadata": {}, + "outputs": [], + "source": [ + "%%capture\n", + "!pip install k-wave-python" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "36d103c4-d931-4a53-aa5e-4a70a794b95e", + "metadata": {}, + "outputs": [], + "source": [ + "from copy import deepcopy\n", + "\n", + "from matplotlib import animation\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "\n", + "from mpl_toolkits.axes_grid1 import make_axes_locatable\n", + "from numpy.typing import NDArray\n", + "from tqdm import tqdm\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.utils.colormap import get_color_map\n", + "from kwave.kspaceFirstOrder2D import kspaceFirstOrder2D\n", + "\n", + "from kwave.options.simulation_execution_options import SimulationExecutionOptions\n", + "from kwave.options.simulation_options import SimulationOptions" + ] + }, + { + "cell_type": "markdown", + "id": "fca4f639-728e-40ea-a98a-42e48e8c7cc0", + "metadata": {}, + "source": [ + "Now need two methods from the fullwave25 library:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4cbc5395-5a34-4f24-a5c0-8fc74158388f", + "metadata": {}, + "outputs": [], + "source": [ + "def map_to_coords(\n", + " map_data: NDArray[np.float64 | np.int64 | np.bool],\n", + " *,\n", + " export_as_xyz: bool = False,\n", + ") -> NDArray[np.int64]:\n", + " \"\"\"Map the mask map to coordinates.\n", + "\n", + " Returns:\n", + " NDArray[np.int64]: An array of coordinates corresponding to non-zero elements in the mask.\n", + "\n", + " \"\"\"\n", + " is_3d = map_data.ndim == 3\n", + " # indices = np.where(map_data.T != 0)\n", + " indices = np.where(map_data != 0)\n", + " if is_3d:\n", + " # out = np.array([indices[2], indices[1], indices[0]]).T\n", + " # out = np.array([indices[2], indices[1], indices[0]]).T\n", + " out = np.array([*indices]).T\n", + "\n", + " if export_as_xyz:\n", + " out = np.stack([out[:, 2], out[:, 1], out[:, 0]], axis=1)\n", + " else:\n", + " out = np.array([*indices]).T\n", + " if export_as_xyz:\n", + " out = np.stack([out[:, 1], out[:, 0]], axis=1)\n", + " return out" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "aee4f3b9-2d62-4fce-8ac6-94156bda9764", + "metadata": {}, + "outputs": [], + "source": [ + "def gaussian_modulated_sinusoidal_signal(\n", + " nt: int,\n", + " duration: float,\n", + " ncycles: int,\n", + " drop_off: int,\n", + " f0: float,\n", + " p0: float,\n", + " delay_sec: float = 0.0,\n", + " i_layer: int | None = None,\n", + " dt_for_layer_delay: float | None = None,\n", + " cfl_for_layer_delay: float | None = None,\n", + ") -> NDArray[np.float64]:\n", + " \"\"\"Generate a pulse signal based on input parameters.\n", + "\n", + " Parameters\n", + " ----------\n", + " nt: int\n", + " Number of time samples of the simulation.\n", + " duration: float\n", + " Total duration of the simulation.\n", + " ncycles: int\n", + " Number of cycles in the pulse.\n", + " drop_off: int\n", + " Controls the pulse decay.\n", + " f0: float\n", + " Frequency of the pulse.\n", + " p0: float\n", + " Amplitude scaling factor.\n", + " delay_sec: float\n", + " Delay in seconds. Default is 0.0.\n", + " i_layer: int\n", + " Index of the layer where the source is located. Default is None.\n", + " This variable is used to shift the pulse signal in time\n", + " so that the signal is emmitted within the transducer layer correctly.\n", + " dt_for_layer_delay: float\n", + " Time step of the simulation. Default is None.\n", + " This variable is used to shift the pulse signal in time\n", + " so that the signal is emmitted within the transducer layer correctly.\n", + " cfl_for_layer_delay: float\n", + " Courant-Friedrichs-Lewy number. Default is None.\n", + " This variable is used to shift the pulse signal in time\n", + " so that the signal is emmitted within the transducer layer correctly.\n", + "\n", + " Returns\n", + " -------\n", + " NDArray[np.float64]: The generated pulse signal.\n", + "\n", + " \"\"\"\n", + " t = (np.arange(0, nt)) / nt * duration - ncycles / f0\n", + " t = t - delay_sec\n", + "\n", + " if i_layer:\n", + " assert dt_for_layer_delay, \"dt must be provided if i_layer is provided\"\n", + " assert cfl_for_layer_delay, \"cfl must be provided if i_layer is provided\"\n", + " t = t - (dt_for_layer_delay / cfl_for_layer_delay) * i_layer\n", + "\n", + " omega0 = 2 * np.pi * f0\n", + " return (\n", + " np.multiply(\n", + " np.exp(\n", + " -((1.05 * t * omega0 / (ncycles * np.pi)) ** (2 * drop_off)),\n", + " ),\n", + " np.sin(t * omega0),\n", + " )\n", + " * p0\n", + " )" + ] + }, + { + "cell_type": "markdown", + "id": "dfb55809-d0a5-4f0d-849b-d2d3b9f9a915", + "metadata": {}, + "source": [ + "Now set up the computation grid. This requires the domain size in metres and the CFL number and points per wavelength (and hence the wavelength). " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3e117208-b260-493e-b988-13d1afdaa5af", + "metadata": {}, + "outputs": [], + "source": [ + "domain_size = (3e-2, 2e-2) # meters\n", + "\n", + "f0 = 3.0e6\n", + "c0 = 1540.0\n", + "wavelength = c0 / f0\n", + "\n", + "cfl = 0.2\n", + "ppw = 12\n", + "\n", + "dx = wavelength / ppw\n", + "dy = dx\n", + "dt = cfl * dx / c0 \n", + "\n", + "Nx = int(np.round(domain_size[0] / dx))\n", + "Ny = int(np.round(domain_size[1] / dy))\n", + "\n", + "grid = kWaveGrid(Vector([Nx, Ny]), Vector([dx, dy]))\n", + "\n", + "# set the time steps\n", + "duration = domain_size[0] / c0 * 2\n", + "Nt = int(np.round(duration / dt))\n", + "grid.setTime(Nt, dt)" + ] + }, + { + "cell_type": "markdown", + "id": "bf3e081e-9954-4021-969e-df16ed58c75c", + "metadata": {}, + "source": [ + "Now the medium" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cbbdbd64-cd15-4a17-b725-e1b51e60bb20", + "metadata": {}, + "outputs": [], + "source": [ + "# Define the base 2D medium arrays\n", + "sound_speed_map = 1540 * np.ones((grid.Nx, grid.Ny)) # m/s\n", + "density_map = 1000 * np.ones((grid.Nx, grid.Ny)) # kg/m^3\n", + "alpha_coeff_map = 0.5 * np.ones((grid.Nx, grid.Ny)) # dB/(MHz^y cm)\n", + "\n", + "# embed an object with different properties in the center of the medium\n", + "obj_x_start = grid.Nx // 3\n", + "obj_x_end = 2 * grid.Nx // 3\n", + "obj_y_start = grid.Ny // 3\n", + "obj_y_end = 2 * grid.Ny // 3\n", + "\n", + "sound_speed_map[obj_x_start:obj_x_end, obj_y_start:obj_y_end] = 1600 # m/s\n", + "density_map[obj_x_start:obj_x_end, obj_y_start:obj_y_end] = 1100 # kg/m^3\n", + "alpha_coeff_map[obj_x_start:obj_x_end, obj_y_start:obj_y_end] = 0.75 # dB/(MHz^y cm)\n", + "\n", + "# setup the Medium instance\n", + "medium = kWaveMedium(sound_speed=sound_speed_map,\n", + " density=density_map,\n", + " alpha_coeff=alpha_coeff_map,\n", + " alpha_power=1.1)" + ] + }, + { + "cell_type": "markdown", + "id": "656a524a-eeaa-41e4-a796-de17fee2a482", + "metadata": {}, + "source": [ + "Now set up the `kSource` instance" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "923c0eee-b8cb-48c2-8e0e-bb1f9f89a9a1", + "metadata": {}, + "outputs": [], + "source": [ + "# initialize the pressure source mask\n", + "p_mask = np.zeros((grid.Nx, grid.Ny), dtype=bool)\n", + "\n", + "# set the source location at the top rows of the grid with specified thickness\n", + "element_thickness_px = 3\n", + "p_mask[0:element_thickness_px, :] = True\n", + "\n", + "# define the pressure source [n_sources, nt]d\n", + "p0 = np.zeros((p_mask.sum(), grid.Nt)) # [n_sources, nt]\n", + "\n", + "p_coordinates = map_to_coords(p_mask)\n", + "\n", + "for i_thickness in range(element_thickness_px):\n", + " # create a gaussian-modulated sinusoidal pulse as the source signal with layer delay\n", + " p0_vec = gaussian_modulated_sinusoidal_signal(nt=grid.Nt, # number of time steps\n", + " f0=f0, # center frequency [Hz]\n", + " duration=duration, # duration [s]\n", + " ncycles=2, # number of cycles\n", + " drop_off=2, # drop off factor\n", + " p0=1e5, # maximum amplitude [Pa]\n", + " i_layer=i_thickness,\n", + " dt_for_layer_delay=grid.dt,\n", + " cfl_for_layer_delay=cfl)\n", + "\n", + " # assign the source signal to the corresponding layer\n", + " n_y = p_coordinates.shape[0] // element_thickness_px\n", + " p0[n_y * i_thickness : n_y * (i_thickness + 1), :] = p0_vec.copy()\n", + "\n", + "source = kSource()\n", + "source.p_mask = p_mask\n", + "source.p = p0" + ] + }, + { + "cell_type": "markdown", + "id": "e5b5bfed-4dd5-4cfa-aa96-cddb9188beae", + "metadata": {}, + "source": [ + "Set up the `kSensor` instance" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "810885a9-4ed5-470d-b150-2c1d5055c3ad", + "metadata": {}, + "outputs": [], + "source": [ + "sensor_mask = np.ones((grid.Nx, grid.Ny), dtype=bool)\n", + "sensor = kSensor()\n", + "sensor.mask = sensor_mask\n", + "sensor.record = [\"p\"]" + ] + }, + { + "cell_type": "markdown", + "id": "cea0ca1e-4561-4bb9-b8ef-6108b0f1c190", + "metadata": {}, + "source": [ + "Perform the simulation" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8ef5a369-abca-4544-a29b-b81c14e21608", + "metadata": {}, + "outputs": [], + "source": [ + "simulation_options = SimulationOptions(pml_auto=True, data_recast=True, save_to_disk=True, save_to_disk_exit=False, pml_inside=False)\n", + "\n", + "execution_options = SimulationExecutionOptions(is_gpu_simulation=True, delete_data=False, verbose_level=2)\n", + "\n", + "sensor_data = kspaceFirstOrder2D(kgrid=deepcopy(grid),\n", + " source=deepcopy(source),\n", + " sensor=deepcopy(sensor),\n", + " medium=deepcopy(medium),\n", + " simulation_options=simulation_options,\n", + " execution_options=execution_options)\n", + "\n", + "propagation_map = np.reshape(sensor_data[\"p\"], (grid.Nt, grid.Nx, grid.Ny), order='F') \n" + ] + }, + { + "cell_type": "markdown", + "id": "256f9054-fd27-416c-bafc-42deb8a7eeb9", + "metadata": {}, + "source": [ + "Visualise the data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ec2a06a1-2f5f-43f9-8aae-df4c70eb1022", + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "\n", + "time_step = propagation_map.shape[0] // 3\n", + "\n", + "cmap = get_color_map()\n", + "\n", + "fig, ax = plt.subplots(1, 1)\n", + "im = ax.imshow(propagation_map[time_step, :, :],\n", + " extent=[grid.y_vec.min() * 1e3, grid.y_vec.max() * 1e3, grid.x_vec.min() * 1e3, grid.x_vec.max() * 1e3],\n", + " vmin=-p_max_plot,\n", + " vmax=p_max_plot,\n", + " cmap=cmap)\n", + "title = \"Snapshot\"\n", + "ax.set_title(title)\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()" + ] + }, + { + "cell_type": "markdown", + "id": "28bbe254-5275-44f6-a61d-0d0e43912433", + "metadata": {}, + "source": [ + "An animation" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ccb611ae-5059-4033-aa95-1691ddf386b6", + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(1, 1, figsize=(4,6))\n", + "\n", + "num_plot_image: int = 50\n", + "skip_every_n_frame = int(grid.Nt / num_plot_image)\n", + "\n", + "c_map = medium.sound_speed\n", + "rho_map = medium.density \n", + "\n", + "start = 0\n", + "end = None\n", + "z_map = c_map * rho_map\n", + "z_map = (z_map - np.min(z_map)) / (np.max(z_map) - np.min(z_map) + 1e-9)\n", + "\n", + "z_map_offset = p_max_plot * 0.8\n", + "\n", + "animation_list = []\n", + "\n", + "# propagation_map = propagation_map.transpose(2, 0, 1)\n", + "for i, p_map_i in tqdm(\n", + " enumerate(propagation_map[::skip_every_n_frame, start:end, start:end]),\n", + " total=len(propagation_map[::skip_every_n_frame, start:end, start:end]),\n", + " desc=\"plotting animation\"):\n", + "\n", + " processed_p_map = p_map_i + z_map_offset * (z_map)\n", + "\n", + " image2 = ax.imshow(\n", + " processed_p_map,\n", + " vmin=-p_max_plot,\n", + " vmax=p_max_plot,\n", + " interpolation=\"nearest\"\n", + " )\n", + " # set text to show the current time step\n", + " text = ax.text(0.5, 1.05, f\"t = {i * skip_every_n_frame} / {propagation_map.shape[0]}\", \n", + " fontsize=4,\n", + " ha=\"center\",\n", + " animated=True,\n", + " transform=ax.transAxes)\n", + " animation_list.append([image2, text])\n", + "\n", + "animation_data = animation.ArtistAnimation(\n", + " fig,\n", + " animation_list,\n", + " interval=150,\n", + " blit=True,\n", + " repeat_delay=500,\n", + ")\n", + "animation_data.save(\"fullwave_plane_wave.mp4\", writer=\"ffmpeg\", dpi=300)\n", + "plt.close(\"all\")" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "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", + "version": "3.13.9" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} From 877a4561aa7ca77a591382bd5c9b9f1d05f0b7ca Mon Sep 17 00:00:00 2001 From: David Sinden Date: Mon, 24 Nov 2025 10:28:01 +0100 Subject: [PATCH 04/11] Update README to include k-wave-python link --- examples/fullwave_simple_plane_wave/README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/fullwave_simple_plane_wave/README.md b/examples/fullwave_simple_plane_wave/README.md index 923baa5c..05d3c025 100644 --- a/examples/fullwave_simple_plane_wave/README.md +++ b/examples/fullwave_simple_plane_wave/README.md @@ -2,5 +2,5 @@ [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/waltsims/k-wave-python/blob/HEAD/examples/fullwave_simple_plane_wave.ipynb) -This example showcases the compatibility between [fullwave25](https://github.com/pinton-lab/fullwave25) and k-wave-python. It directly copies two methods over. +This example showcases the compatibility between [fullwave25](https://github.com/pinton-lab/fullwave25) and [k-wave-python](https://github.com/waltsims/k-wave-python). It directly copies two methods over. From ef81f9f172a92a33d4a1cbf2e2dd212889ab5e9b Mon Sep 17 00:00:00 2001 From: David Sinden Date: Mon, 24 Nov 2025 10:33:03 +0100 Subject: [PATCH 05/11] Add HTML display for animation in notebook --- .../fullwave_simple_plane_wave.ipynb | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/examples/fullwave_simple_plane_wave/fullwave_simple_plane_wave.ipynb b/examples/fullwave_simple_plane_wave/fullwave_simple_plane_wave.ipynb index 4aa32997..46e92241 100644 --- a/examples/fullwave_simple_plane_wave/fullwave_simple_plane_wave.ipynb +++ b/examples/fullwave_simple_plane_wave/fullwave_simple_plane_wave.ipynb @@ -32,6 +32,7 @@ "source": [ "from copy import deepcopy\n", "\n", + "from IPython.display import HTML", "from matplotlib import animation\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", @@ -361,6 +362,7 @@ "source": [ "%matplotlib inline\n", "\n", + "p_max_plot = np.abs(propagation_map).max().item() / 4", "time_step = propagation_map.shape[0] // 3\n", "\n", "cmap = get_color_map()\n", @@ -443,7 +445,8 @@ " repeat_delay=500,\n", ")\n", "animation_data.save(\"fullwave_plane_wave.mp4\", writer=\"ffmpeg\", dpi=300)\n", - "plt.close(\"all\")" + "plt.close(\"all\")", + "HTML(animation_data.to_html5_video())" ] } ], From c0e2fad60f494af82c8937aab095a88d3f27c582 Mon Sep 17 00:00:00 2001 From: David Sinden Date: Mon, 24 Nov 2025 10:35:36 +0100 Subject: [PATCH 06/11] Fix import statement formatting in Jupyter notebook --- .../fullwave_simple_plane_wave/fullwave_simple_plane_wave.ipynb | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/fullwave_simple_plane_wave/fullwave_simple_plane_wave.ipynb b/examples/fullwave_simple_plane_wave/fullwave_simple_plane_wave.ipynb index 46e92241..5a89fe26 100644 --- a/examples/fullwave_simple_plane_wave/fullwave_simple_plane_wave.ipynb +++ b/examples/fullwave_simple_plane_wave/fullwave_simple_plane_wave.ipynb @@ -32,7 +32,7 @@ "source": [ "from copy import deepcopy\n", "\n", - "from IPython.display import HTML", + "from IPython.display import HTML\n", "from matplotlib import animation\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", From 4064d079d34f319d72823d7a891a97fee9faed31 Mon Sep 17 00:00:00 2001 From: David Sinden Date: Mon, 24 Nov 2025 10:37:15 +0100 Subject: [PATCH 07/11] Remove commented out code for clarity --- .../fullwave_simple_plane_wave/fullwave_simple_plane_wave.ipynb | 2 -- 1 file changed, 2 deletions(-) diff --git a/examples/fullwave_simple_plane_wave/fullwave_simple_plane_wave.ipynb b/examples/fullwave_simple_plane_wave/fullwave_simple_plane_wave.ipynb index 5a89fe26..15651e09 100644 --- a/examples/fullwave_simple_plane_wave/fullwave_simple_plane_wave.ipynb +++ b/examples/fullwave_simple_plane_wave/fullwave_simple_plane_wave.ipynb @@ -83,8 +83,6 @@ " # indices = np.where(map_data.T != 0)\n", " indices = np.where(map_data != 0)\n", " if is_3d:\n", - " # out = np.array([indices[2], indices[1], indices[0]]).T\n", - " # out = np.array([indices[2], indices[1], indices[0]]).T\n", " out = np.array([*indices]).T\n", "\n", " if export_as_xyz:\n", From 645ed942317f687a0b40c157f8d97f26c9e3fe4d Mon Sep 17 00:00:00 2001 From: David Sinden Date: Mon, 24 Nov 2025 10:52:35 +0100 Subject: [PATCH 08/11] Fix type hint and formatting in Jupyter notebook --- .../fullwave_simple_plane_wave.ipynb | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/fullwave_simple_plane_wave/fullwave_simple_plane_wave.ipynb b/examples/fullwave_simple_plane_wave/fullwave_simple_plane_wave.ipynb index 15651e09..5779733a 100644 --- a/examples/fullwave_simple_plane_wave/fullwave_simple_plane_wave.ipynb +++ b/examples/fullwave_simple_plane_wave/fullwave_simple_plane_wave.ipynb @@ -69,7 +69,7 @@ "outputs": [], "source": [ "def map_to_coords(\n", - " map_data: NDArray[np.float64 | np.int64 | np.bool],\n", + " map_data: NDArray[np.float64 | np.int64 | bool],\n", " *,\n", " export_as_xyz: bool = False,\n", ") -> NDArray[np.int64]:\n", @@ -360,7 +360,7 @@ "source": [ "%matplotlib inline\n", "\n", - "p_max_plot = np.abs(propagation_map).max().item() / 4", + "p_max_plot = np.abs(propagation_map).max().item() / 4\n", "time_step = propagation_map.shape[0] // 3\n", "\n", "cmap = get_color_map()\n", @@ -443,7 +443,7 @@ " repeat_delay=500,\n", ")\n", "animation_data.save(\"fullwave_plane_wave.mp4\", writer=\"ffmpeg\", dpi=300)\n", - "plt.close(\"all\")", + "plt.close(\"all\")\n", "HTML(animation_data.to_html5_video())" ] } From 5aa36b020ccafffeda89d2b47229c50be2b24631 Mon Sep 17 00:00:00 2001 From: David Sinden Date: Mon, 24 Nov 2025 10:53:54 +0100 Subject: [PATCH 09/11] Remove animation save and plot close commands Removed saving of animation to MP4 and closing of plots. --- .../fullwave_simple_plane_wave/fullwave_simple_plane_wave.ipynb | 2 -- 1 file changed, 2 deletions(-) diff --git a/examples/fullwave_simple_plane_wave/fullwave_simple_plane_wave.ipynb b/examples/fullwave_simple_plane_wave/fullwave_simple_plane_wave.ipynb index 5779733a..a2b7a0d6 100644 --- a/examples/fullwave_simple_plane_wave/fullwave_simple_plane_wave.ipynb +++ b/examples/fullwave_simple_plane_wave/fullwave_simple_plane_wave.ipynb @@ -442,8 +442,6 @@ " blit=True,\n", " repeat_delay=500,\n", ")\n", - "animation_data.save(\"fullwave_plane_wave.mp4\", writer=\"ffmpeg\", dpi=300)\n", - "plt.close(\"all\")\n", "HTML(animation_data.to_html5_video())" ] } From a418f1390bab2c7f3d66d0e7d2f49ddeaa4d7b97 Mon Sep 17 00:00:00 2001 From: David Sinden Date: Mon, 24 Nov 2025 11:12:43 +0100 Subject: [PATCH 10/11] Reorganize imports in fullwave_simple_plane_wave.py --- .../fullwave_simple_plane_wave.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/examples/fullwave_simple_plane_wave/fullwave_simple_plane_wave.py b/examples/fullwave_simple_plane_wave/fullwave_simple_plane_wave.py index 97baeea0..cfcb0007 100644 --- a/examples/fullwave_simple_plane_wave/fullwave_simple_plane_wave.py +++ b/examples/fullwave_simple_plane_wave/fullwave_simple_plane_wave.py @@ -7,9 +7,8 @@ from matplotlib import animation import matplotlib.pyplot as plt -import numpy as np - from mpl_toolkits.axes_grid1 import make_axes_locatable +import numpy as np from numpy.typing import NDArray from tqdm import tqdm @@ -18,11 +17,10 @@ from kwave.kmedium import kWaveMedium from kwave.ksensor import kSensor from kwave.ksource import kSource -from kwave.utils.colormap import get_color_map from kwave.kspaceFirstOrder2D import kspaceFirstOrder2D - from kwave.options.simulation_execution_options import SimulationExecutionOptions from kwave.options.simulation_options import SimulationOptions +from kwave.utils.colormap import get_color_map From 7ac735cc8d2f8280a09c17c279042b073e1a9350 Mon Sep 17 00:00:00 2001 From: David Sinden Date: Mon, 24 Nov 2025 11:15:41 +0100 Subject: [PATCH 11/11] Fix import order and update type hint for map_data --- .../fullwave_simple_plane_wave.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/fullwave_simple_plane_wave/fullwave_simple_plane_wave.py b/examples/fullwave_simple_plane_wave/fullwave_simple_plane_wave.py index cfcb0007..2f13e1aa 100644 --- a/examples/fullwave_simple_plane_wave/fullwave_simple_plane_wave.py +++ b/examples/fullwave_simple_plane_wave/fullwave_simple_plane_wave.py @@ -5,10 +5,10 @@ from copy import deepcopy -from matplotlib import animation import matplotlib.pyplot as plt -from mpl_toolkits.axes_grid1 import make_axes_locatable import numpy as np +from matplotlib import animation +from mpl_toolkits.axes_grid1 import make_axes_locatable from numpy.typing import NDArray from tqdm import tqdm @@ -25,7 +25,7 @@ def map_to_coords( - map_data: NDArray[np.float64 | np.int64 | np.bool], + map_data: NDArray[np.float64 | np.int64 | np.bool_ | bool], *, export_as_xyz: bool = False, ) -> NDArray[np.int64]: