diff --git a/.gitignore b/.gitignore index 18a37dc4..68cc2f63 100644 --- a/.gitignore +++ b/.gitignore @@ -4,6 +4,7 @@ # Virtual environments .venv/ +venv/ # Compiled files *.pyc diff --git a/examples/README.md b/examples/README.md index 77238d05..6927bfac 100644 --- a/examples/README.md +++ b/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)) @@ -25,9 +26,10 @@ 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)) +- [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)) + ## Contributing new examples @@ -38,10 +40,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="^%" ` 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] \". Thanks for contributing to k-wave-python! diff --git a/examples/sd_focussed_detector_2D/README.md b/examples/sd_focussed_detector_2D/README.md new file mode 100644 index 00000000..9a05f008 --- /dev/null +++ b/examples/sd_focussed_detector_2D/README.md @@ -0,0 +1,7 @@ +# Focussed Detector In 2D 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_2D/sd_focussed_detector_2D.ipynb) + +This example shows how k-Wave-python can be used to model the output of a focussed semicircular 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_2D.php). diff --git a/examples/sd_focussed_detector_2D/sd_focussed_detector_2D.ipynb b/examples/sd_focussed_detector_2D/sd_focussed_detector_2D.ipynb new file mode 100644 index 00000000..99115624 --- /dev/null +++ b/examples/sd_focussed_detector_2D/sd_focussed_detector_2D.ipynb @@ -0,0 +1,207 @@ +{ + "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 2D Example\n", + "This example shows how k-Wave-python can be used to model the output of a focused semicircular detector, where the directionality arises from spatially averaging across the detector surface. Unlike the original example in k-Wave, this example does not visualize the simulation, as this functionality is not intrinsically supported by the accelerated binaries." + ] + }, + { + "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.kspaceFirstOrder2D import kspaceFirstOrder2DC\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_disc, make_circle\n", + "from kwave.utils.data import scale_SI\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# create the computational grid\n", + "grid_size = Vector([180, 180]) # [grid points]\n", + "grid_spacing = Vector([0.1e-3, 0.1e-3]) # [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 a sensor as part of a circle centred on the grid\n", + "sensor_radius = 65 # [grid points]\n", + "arc_angle = np.pi # [rad]\n", + "sensor_mask = make_circle(grid_size, grid_size // 2 + 1, sensor_radius, arc_angle)\n", + "sensor = kSensor(sensor_mask)\n", + "\n", + "# define the array of temporal points\n", + "t_end = 11e-6 # [s]\n", + "_ = kgrid.makeTime(medium.sound_speed, t_end=t_end)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Define simulation parameters" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "input_filename = \"example_sd_focused_2d_input.h5\"\n", + "pathname = gettempdir()\n", + "input_file_full_path = os.path.join(pathname, input_filename)\n", + "simulation_options = SimulationOptions(save_to_disk=True, input_filename=input_filename, data_path=pathname)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Run simulation with first source" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# place a disc-shaped source near the focus of the detector\n", + "source = kSource()\n", + "source.p0 = 2 * make_disc(grid_size, grid_size / 2, 4)\n", + "\n", + "# run the simulation\n", + "sensor_data1 = kspaceFirstOrder2DC(\n", + " medium=medium,\n", + " kgrid=kgrid,\n", + " source=deepcopy(source),\n", + " sensor=sensor,\n", + " simulation_options=simulation_options,\n", + " execution_options=SimulationExecutionOptions(),\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "sensor_data1['p'].shape\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Run simulation with second source" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# place a disc-shaped source horizontally shifted from the focus of the detector\n", + "source.p0 = 2 * make_disc(grid_size, grid_size / 2 + [0, 20], 4)\n", + "\n", + "sensor_data2 = kspaceFirstOrder2DC(\n", + " medium=medium,\n", + " kgrid=kgrid,\n", + " source=deepcopy(source),\n", + " sensor=sensor,\n", + " simulation_options=simulation_options,\n", + " execution_options=SimulationExecutionOptions(),\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Visualize recorded data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "sensor_output1 = np.sum(sensor_data1['p'], axis=1) / np.sum(sensor.mask)\n", + "sensor_output2 = np.sum(sensor_data2['p'], axis=1) / np.sum(sensor.mask)\n", + "\n", + "t_sc, t_scale, t_prefix, _ = scale_SI(t_end)\n", + "t_array = kgrid.t_array.squeeze() * t_scale\n", + "\n", + "plt.plot(t_array, sensor_output1, 'k')\n", + "plt.plot(t_array, sensor_output2, 'r')\n", + "\n", + "plt.xlabel('Time [' + t_prefix + 's]')\n", + "plt.ylabel('Average Pressure Measured Over Detector [au]')\n", + "plt.legend([\n", + " f\"Source on focus, sum(output^2) = {round(np.sum(sensor_output1**2) * 100) / 100}\",\n", + " f\"Source off focus, sum(output^2) = {round(np.sum(sensor_output2**2) * 100) / 100}\"\n", + "])\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_2D/sd_focussed_detector_2D.py b/examples/sd_focussed_detector_2D/sd_focussed_detector_2D.py new file mode 100644 index 00000000..eba51063 --- /dev/null +++ b/examples/sd_focussed_detector_2D/sd_focussed_detector_2D.py @@ -0,0 +1,108 @@ +# # Focussed Detector In 2D Example +# This example shows how k-Wave-python can be used to model the output of a focused semicircular detector, where the directionality arises from spatially averaging across the detector surface. Unlike the original example in k-Wave, this example does not visualize the simulation, as this functionality is not intrinsically supported by the accelerated binaries. + +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.kspaceFirstOrder2D import kspaceFirstOrder2DC +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_disc, make_circle +from kwave.utils.data import scale_SI + + +# In[3]: + + +# create the computational grid +grid_size = Vector([180, 180]) # [grid points] +grid_spacing = Vector([0.1e-3, 0.1e-3]) # [m] +kgrid = kWaveGrid(grid_size, grid_spacing) + +# define the properties of the propagation medium +medium = kWaveMedium(sound_speed=1500) + +# define a sensor as part of a circle centred on the grid +sensor_radius = 65 # [grid points] +arc_angle = np.pi # [rad] +sensor_mask = make_circle(grid_size, grid_size // 2 + 1, sensor_radius, arc_angle) +sensor = kSensor(sensor_mask) + +# define the array of temporal points +t_end = 11e-6 # [s] +_ = kgrid.makeTime(medium.sound_speed, t_end=t_end) + + +# ## Define simulation parameters +input_filename = "example_sd_focused_2d_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) + + +# ## Run simulation with first source +# place a disc-shaped source near the focus of the detector +source = kSource() +source.p0 = 2 * make_disc(grid_size, grid_size / 2, 4) + +# run the simulation +sensor_data1 = kspaceFirstOrder2DC( + medium=medium, + kgrid=kgrid, + source=deepcopy(source), + sensor=sensor, + simulation_options=simulation_options, + execution_options=SimulationExecutionOptions(), +) + + +sensor_data1["p"].shape + + +# ## Run simulation with second source + + +# place a disc-shaped source horizontally shifted from the focus of the detector +source.p0 = 2 * make_disc(grid_size, grid_size / 2 + [0, 20], 4) + +sensor_data2 = kspaceFirstOrder2DC( + medium=medium, + kgrid=kgrid, + source=deepcopy(source), + sensor=sensor, + simulation_options=simulation_options, + execution_options=SimulationExecutionOptions(), +) + + +# ## Visualize recorded data + + +sensor_output1 = np.sum(sensor_data1["p"], axis=1) / np.sum(sensor.mask) +sensor_output2 = np.sum(sensor_data2["p"], axis=1) / np.sum(sensor.mask) + +t_sc, t_scale, t_prefix, _ = scale_SI(t_end) +t_array = kgrid.t_array.squeeze() * t_scale + +plt.plot(t_array, sensor_output1, "k") +plt.plot(t_array, sensor_output2, "r") + +plt.xlabel("Time [" + t_prefix + "s]") +plt.ylabel("Average Pressure Measured Over Detector [au]") +plt.legend( + [ + f"Source on focus, sum(output^2) = {round(np.sum(sensor_output1**2) * 100) / 100}", + f"Source off focus, sum(output^2) = {round(np.sum(sensor_output2**2) * 100) / 100}", + ] +) + +plt.show() diff --git a/examples/us_beam_patterns/README.md b/examples/us_beam_patterns/README.md new file mode 100644 index 00000000..e88c42a0 --- /dev/null +++ b/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). \ No newline at end of file diff --git a/examples/us_beam_patterns/us_beam_patterns.ipynb b/examples/us_beam_patterns/us_beam_patterns.ipynb new file mode 100644 index 00000000..0a1abcfd --- /dev/null +++ b/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 +} diff --git a/examples/us_beam_patterns/us_beam_patterns.py b/examples/us_beam_patterns/us_beam_patterns.py new file mode 100644 index 00000000..bcb62f1f --- /dev/null +++ b/examples/us_beam_patterns/us_beam_patterns.py @@ -0,0 +1,274 @@ +#!/usr/bin/env python +# coding: utf-8 + +# In[ ]: + + +import matplotlib.pyplot as plt +import numpy as np + +from kwave.data import Vector +from kwave.kgrid import kWaveGrid +from kwave.kspaceFirstOrder3D import kspaceFirstOrder3D +from kwave.ksensor import kSensor +from kwave.ktransducer import kWaveTransducerSimple, NotATransducer +from kwave.kWaveSimulation import SimulationOptions +from kwave.kmedium import kWaveMedium +from kwave.utils.filters import spect + +from kwave.options.simulation_execution_options import SimulationExecutionOptions +from kwave.utils.dotdictionary import dotdict +from kwave.utils.math import find_closest +from kwave.utils.signals import tone_burst + + +# In[ ]: + + +# simulation settings +DATA_CAST = "single" +# set to 'xy' or 'xz' to generate the beam pattern in different planes +MASK_PLANE = "xy" +# set to true to compute the rms or peak beam patterns, set to false to compute the harmonic beam patterns +USE_STATISTICS = True + +# define the grid +PML_X_SIZE = 20 +PML_Y_SIZE = 10 +PML_Z_SIZE = 10 +Nx = 128 - 2 * PML_X_SIZE +Ny = 64 - 2 * PML_Y_SIZE +Nz = 64 - 2 * PML_Z_SIZE +x = 40e-3 +dx = x / Nx +dy = dx +dz = dx + +kgrid = kWaveGrid([Nx, Ny, Nz], [dx, dy, dz]) + + +# In[ ]: + + +# define the medium +medium = kWaveMedium(sound_speed=1540, density=1000, alpha_coeff=0.75, alpha_power=1.5, BonA=6) + +# create the time array - using a time == time to travel the hypot of the grid +t_end = 45e-6 # alternatively, use np.sqrt(kgrid.x_size ** 2 + kgrid.y_size ** 2) / medium.sound_speed + +kgrid.makeTime(medium.sound_speed, t_end=t_end) + + +# In[ ]: + + +# define the input signal +source_strength = 1e6 +tone_burst_freq = 0.5e6 +tone_burst_cycles = 5 +input_signal = tone_burst(1 / kgrid.dt, tone_burst_freq, tone_burst_cycles) +input_signal = (source_strength / (medium.sound_speed * medium.density)) * input_signal + + +# In[ ]: + + +# define the transducer +transducer = dotdict() +transducer.number_elements = 32 +transducer.element_width = 1 +transducer.element_length = 12 +transducer.element_spacing = 0 +transducer.radius = np.inf + +# calculate the width of the transducer in grid points +transducer_width = transducer.number_elements * transducer.element_width + (transducer.number_elements - 1) * transducer.element_spacing + +# use this to position the transducer in the middle of the computational grid +transducer.position = np.round([1, Ny / 2 - transducer_width / 2, Nz / 2 - transducer.element_length / 2]) +transducer = kWaveTransducerSimple(kgrid, **transducer) + + +# In[ ]: + + +not_transducer = dotdict() +not_transducer.sound_speed = medium.sound_speed # sound speed [m/s] +not_transducer.focus_distance = 20e-3 # focus distance [m] +not_transducer.elevation_focus_distance = 19e-3 # focus distance in the elevation plane [m] +not_transducer.steering_angle = 0 # steering angle [degrees] +not_transducer.transmit_apodization = "Rectangular" +not_transducer.receive_apodization = "Rectangular" +not_transducer.active_elements = np.ones((transducer.number_elements, 1)) +not_transducer.input_signal = input_signal + +not_transducer = NotATransducer(transducer, kgrid, **not_transducer) + + +# In[ ]: + + +sensor_mask = np.zeros((Nx, Ny, Nz)) + +if MASK_PLANE == "xy": + sensor_mask[:, :, Nz // 2] = 1 + # store y axis properties + Nj = Ny + j_vec = kgrid.y_vec + j_label = "y" +elif MASK_PLANE == "xz": + sensor_mask[:, Ny // 2, :] = 1 + # store z axis properties + Nj = Nz + j_vec = kgrid.z_vec + j_label = "z" + + +# In[ ]: + + +sensor = kSensor(sensor_mask) +if USE_STATISTICS: + sensor.record = ["p_rms", "p_max"] + +simulation_options = SimulationOptions( + pml_inside=False, + pml_size=Vector([PML_X_SIZE, PML_Y_SIZE, PML_Z_SIZE]), + data_cast=DATA_CAST, + save_to_disk=True, +) + +if not USE_STATISTICS: + simulation_options.stream_to_disk = "harmonic_data.h5" + +execution_options = SimulationExecutionOptions(is_gpu_simulation=True) + +sensor_data = kspaceFirstOrder3D( + medium=medium, + kgrid=kgrid, + source=not_transducer, + sensor=sensor, + simulation_options=simulation_options, + execution_options=execution_options, +) + + +# In[ ]: + + +if USE_STATISTICS: + fig, axes = plt.subplots(1, 2) + fig.set_figwidth(8) + fig.set_figheight(6) + + for ind, measures in enumerate(sensor.record): + im1 = axes[ind].imshow( + sensor_data[measures].reshape([Nj, Nx]).T * 1e-6, + extent=[ + min(j_vec * 1e3), + max(j_vec * 1e3), + min((kgrid.x_vec - min(kgrid.x_vec)) * 1e3), + max((kgrid.x_vec - min(kgrid.x_vec)) * 1e3), + ], + aspect="auto", + cmap="jet", + ) + axes[ind].set_xlabel("y-position (mm)") + axes[ind].set_ylabel("x-position (mm)") + title_text = f"Total Beam Pattern Using {str.upper(measures.split('_')[1])} \n of Recorded Pressure" + axes[ind].set_title(title_text) + axes[ind].set_yticks(axes[ind].get_yticks().tolist()) + axes[ind].set_yticklabels(axes[ind].get_yticklabels()[::-1]) + axes[ind].set_xlim([-10, 10]) + fig.colorbar(im1, label="Pressure [MPa]", orientation="vertical") + + plt.tight_layout() + plt.show() + + +# In[ ]: + + +if not USE_STATISTICS: + sensor_data_array = np.reshape(sensor_data["p"], [kgrid.Nt, kgrid.Ny, kgrid.Nx]).transpose(2, 1, 0) + # compute the amplitude spectrum + [freq, amp_spect, _] = spect(sensor_data_array, 1 / kgrid.dt, dim=2) + + # compute the index at which the source frequency and its harmonics occur + [f1_value, f1_index] = find_closest(freq, tone_burst_freq) + [f2_value, f2_index] = find_closest(freq, 2 * tone_burst_freq) + + # extract the amplitude at the source frequency and store + beam_pattern_f1 = np.squeeze(amp_spect[:, :, f1_index]) + + # extract the amplitude at the second harmonic and store + beam_pattern_f2 = np.squeeze(amp_spect[:, :, f2_index]) + + # extract the integral of the total amplitude spectrum + beam_pattern_total = np.squeeze(np.sum(amp_spect, axis=2)) + + +# In[ ]: + + +if not USE_STATISTICS: + fig, axes = plt.subplots(1, 3) + fig.set_figwidth(8) + fig.set_figheight(6) + + for ind, measure in enumerate([beam_pattern_f1, beam_pattern_f2, beam_pattern_total]): + im1 = axes[ind].imshow( + np.squeeze(measure) * 1e-6, + extent=[ + min(j_vec * 1e3), + max(j_vec * 1e3), + min((kgrid.x_vec - min(kgrid.x_vec)) * 1e3), + max((kgrid.x_vec - min(kgrid.x_vec)) * 1e3), + ], + aspect="auto", + cmap="jet", + ) + axes[ind].set_xlabel("y-position (mm)") + axes[ind].set_ylabel("x-position (mm)") + + axes[ind].set_yticks(axes[ind].get_yticks().tolist()) + axes[ind].set_yticklabels(axes[ind].get_yticklabels()[::-1]) + + axes[0].set_title("Beam Pattern \nAt Source Fundamental") + axes[1].set_title("Beam Pattern \nAt Second Harmonic") + axes[2].set_title("Total Beam Pattern\nUsing Integral Of Recorded Pressure") + + plt.tight_layout() + plt.show() + + +# In[ ]: + + +if not USE_STATISTICS: + # Compute the directivity at each of the harmonics + directivity_f1 = beam_pattern_f1[round(not_transducer.focus_distance / dx), :] + directivity_f2 = beam_pattern_f2[round(not_transducer.focus_distance / dx), :] + + # Normalize the directivity + directivity_f1 /= np.max(directivity_f1) + directivity_f2 /= np.max(directivity_f2) + + # Compute relative angles from the transducer + if MASK_PLANE == "xy": + horz_axis = ((np.arange(Ny) + 1) - Ny / 2) * dy + else: + horz_axis = ((np.arange(Nz) + 1) - Nz / 2) * dz + + angles = 180 * np.arctan2(horz_axis, not_transducer.focus_distance) / np.pi + + # Plot the directivity + plt.figure() + plt.plot(angles, directivity_f1, "k-", label="Fundamental") + plt.plot(angles, directivity_f2, "k--", label="Second Harmonic") + plt.axis("tight") + plt.xlabel("Angle [deg]") + plt.ylabel("Normalized Amplitude") + plt.legend(loc=2) + plt.xlim([-25, 25]) + plt.show() diff --git a/examples/us_defining_transducer/README.md b/examples/us_defining_transducer/README.md index 3efd92f2..7472dfd6 100644 --- a/examples/us_defining_transducer/README.md +++ b/examples/us_defining_transducer/README.md @@ -1,7 +1,8 @@ -# Simulating B-mode Ultrasound Images Example +# Defining An Ultrasound Transducer 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_defining_transducer/us_defining_transducer.ipynb) -This example illustrates how this object is created and can be used to simulate the field produced by an ultrasound transducer. +This example illustrates how kWaveTransducerSimple can be used to simulate the field produced by an ultrasound transducer. + To read more, visit the [original example page](http://www.k-wave.org/documentation/example_us_defining_transducer.php). diff --git a/examples/us_defining_transducer/us_defining_transducer.ipynb b/examples/us_defining_transducer/us_defining_transducer.ipynb index 5b9f080e..f62ced92 100644 --- a/examples/us_defining_transducer/us_defining_transducer.ipynb +++ b/examples/us_defining_transducer/us_defining_transducer.ipynb @@ -21,6 +21,7 @@ "\n", "from kwave.data import Vector\n", "from kwave.kgrid import kWaveGrid\n", + "from kwave.kmedium import kWaveMedium\n", "from kwave.kspaceFirstOrder3D import kspaceFirstOrder3D\n", "from kwave.ksensor import kSensor\n", "from kwave.ktransducer import kWaveTransducerSimple, NotATransducer\n", @@ -200,8 +201,7 @@ "axes[1].set_xlabel('Frequency [MHz]')\n", "axes[1].set_ylabel('Amplitude Spectrum [au]')\n", "f_max = medium.sound_speed / (2*dx)\n", - "axes[1].set_xlim([0, f_max/1e6])\n", - "\n" + "axes[1].set_xlim([0, f_max/1e6])\n" ] }, { diff --git a/examples/us_defining_transducer/us_defining_transducer.py b/examples/us_defining_transducer/us_defining_transducer.py index c77e7696..ca5b6ffc 100644 --- a/examples/us_defining_transducer/us_defining_transducer.py +++ b/examples/us_defining_transducer/us_defining_transducer.py @@ -1,4 +1,3 @@ -# %% import matplotlib.pyplot as plt import numpy as np @@ -16,7 +15,6 @@ from kwave.utils.plot import voxel_plot from kwave.utils.signals import tone_burst -# %% # simulation settings DATA_CAST = "single" @@ -34,7 +32,6 @@ kgrid = kWaveGrid([Nx, Ny, Nz], [dx, dy, dz]) -# %% # define the medium medium = kWaveMedium(sound_speed=1540, density=1000, alpha_coeff=0.75, alpha_power=1.5, BonA=6) @@ -42,7 +39,6 @@ t_end = 40e-6 kgrid.makeTime(medium.sound_speed, t_end=t_end) -# %% # define the input signal source_strength = 1e6 tone_burst_freq = 0.5e6 @@ -50,8 +46,6 @@ input_signal = tone_burst(1 / kgrid.dt, tone_burst_freq, tone_burst_cycles) input_signal = (source_strength / (medium.sound_speed * medium.density)) * input_signal - -# %% # define the transducer transducer = dotdict() transducer.number_elements = 72 @@ -60,14 +54,10 @@ transducer.element_spacing = 0 transducer.radius = np.inf -# calculate the width of the transducer in grid points -transducer_width = transducer.number_elements * transducer.element_width + (transducer.number_elements - 1) * transducer.element_spacing - # use this to position the transducer in the middle of the computational grid -transducer.position = np.round([1, Ny / 2 - transducer_width / 2, Nz / 2 - transducer.element_length / 2]) +transducer.position = np.round([1, Ny / 2 - Nx // 2, Nz / 2 - transducer.element_length / 2]) transducer = kWaveTransducerSimple(kgrid, **transducer) -# %% not_transducer = dotdict() not_transducer.sound_speed = medium.sound_speed # sound speed [m/s] not_transducer.focus_distance = 20e-3 # focus distance [m] @@ -81,10 +71,8 @@ not_transducer = NotATransducer(transducer, kgrid, **not_transducer) -# %% voxel_plot(np.single(not_transducer.all_elements_mask)) -# %% # define sensor mask sensor_mask = np.zeros((Nx, Ny, Nz)) sensor_mask[Nx // 4, Ny // 2, Nz // 2] = 1 @@ -112,16 +100,12 @@ execution_options=execution_options, ) - -# %% - padded_input_signal = np.concatenate((input_signal, np.zeros((1, 2 * np.shape(input_signal)[1]))), axis=1) f_input, as_input, _ = spect(padded_input_signal, 1 / kgrid.dt) _, as_1, _ = spect(sensor_data["p"][:, 0], 1 / kgrid.dt) _, as_2, _ = spect(sensor_data["p"][:, 1], 1 / kgrid.dt) f, as_3, _ = spect(sensor_data["p"][:, 2], 1 / kgrid.dt) -# %% fig, axes = plt.subplots(2, 1) axes[0].plot(np.arange(0, input_signal.shape[-1]) * kgrid.dt * 1e6, input_signal.T, "k-") axes[0].set_xlabel("Time [\mus]") @@ -134,8 +118,6 @@ f_max = medium.sound_speed / (2 * dx) axes[1].set_xlim([0, f_max / 1e6]) - -# %% # Creating a dictionary with the step labels as keys sensor_positions = { "Sensor Position 1": sensor_data["p"][:, 0], diff --git a/kwave/executor.py b/kwave/executor.py index 1865b177..8456fe5b 100644 --- a/kwave/executor.py +++ b/kwave/executor.py @@ -32,16 +32,12 @@ def _make_binary_executable(self): raise e def run_simulation(self, input_filename: str, output_filename: str, options: str): - command = ( - f"{self.execution_options.system_string} " - f"{self.execution_options.binary_path} " - f"-i {input_filename} " - f"-o {output_filename} " - f"{options}" - ) + command = [str(self.execution_options.binary_path), "-i", input_filename, "-o", output_filename, options] try: - with subprocess.Popen(command, stdout=subprocess.PIPE, stderr=subprocess.PIPE, shell=True, text=True) as proc: + with subprocess.Popen( + command, env=self.execution_options.env_vars, stdout=subprocess.PIPE, stderr=subprocess.PIPE, text=True + ) as proc: stdout, stderr = "", "" if self.execution_options.show_sim_log: # Stream stdout in real-time diff --git a/kwave/options/simulation_execution_options.py b/kwave/options/simulation_execution_options.py index d689ffd3..3dcf42c6 100644 --- a/kwave/options/simulation_execution_options.py +++ b/kwave/options/simulation_execution_options.py @@ -1,9 +1,9 @@ from dataclasses import dataclass +import os from typing import Optional, Union from kwave import PLATFORM, BINARY_PATH from kwave.ksensor import kSensor -from kwave.utils.checks import is_unix @dataclass @@ -117,29 +117,22 @@ def get_options_string(self, sensor: kSensor) -> str: return options_string @property - def system_string(self): - # set OS string for setting environment variables - if is_unix(): - env_set_str = "" - sys_sep_str = " " - else: - env_set_str = "set " - sys_sep_str = " & " + def env_vars(self) -> dict: + env = os.environ - # set system string to define domain for thread migration - system_string = env_set_str + "OMP_PLACES=cores" + sys_sep_str + if PLATFORM != "darwin": + env.update({"OMP_PLACES": "cores"}) if self.thread_binding is not None: + if PLATFORM == "darwin": + raise ValueError("Thread binding is not supported in MacOS.") # read the parameters and update the system options if self.thread_binding: - system_string = system_string + " " + env_set_str + "OMP_PROC_BIND=SPREAD" + sys_sep_str + env.update({"OMP_PROC_BIND": "SPREAD"}) else: - system_string = system_string + " " + env_set_str + "OMP_PROC_BIND=CLOSE" + sys_sep_str + env.update({"OMP_PROC_BIND": "CLOSE"}) else: - # set to round-robin over places - system_string = system_string + " " + env_set_str + "OMP_PROC_BIND=SPREAD" + sys_sep_str - - if self.system_call: - system_string = system_string + " " + self.system_call + sys_sep_str + if PLATFORM != "darwin": + env.update({"OMP_PROC_BIND": "SPREAD"}) - return system_string + return env diff --git a/kwave/utils/filters.py b/kwave/utils/filters.py index 50aab667..bbd7c31f 100644 --- a/kwave/utils/filters.py +++ b/kwave/utils/filters.py @@ -98,18 +98,15 @@ def spect( # automatically set dimension to first non - singleton dimension if dim == "auto": - dim_index = 0 - while dim_index <= len(sz): - if sz[dim_index] > 1: - dim = dim_index - break - dim_index = dim_index + 1 + dim = np.argmax(np.array(sz) > 1) + if sz[dim] <= 1: + raise ValueError("All dimensions are singleton; unable to determine valid dimension.") # assign the number of points being analysed func_length = sz[dim] # set the length of the FFT - if not fft_len > func_length: + if fft_len <= 0 or fft_len < func_length: if power_two: # find an appropriate FFT length of the form 2 ^ N that is equal to or # larger than the length of the input signal @@ -119,34 +116,32 @@ def spect( fft_len = func_length # window the signal, reshaping the window to be in the correct direction - win, coherent_gain = get_win(func_length, window, symmetric=False) - win = np.reshape(win, tuple(([1] * dim + [func_length] + [1] * (len(sz) - 2)))) + win, coherent_gain = get_win(func_length, type_=window, symmetric=False) + win_shape = [1] * len(sz) + win_shape[dim] = func_length + win = np.reshape(win, tuple(win_shape)) func = win * func # compute the fft using the defined FFT length, if fft_len > # func_length, the input signal is padded with zeros - func_fft = fft(func, fft_len, dim) + func_fft = np.fft.fft(func, n=fft_len, axis=dim) # correct for the magnitude scaling of the FFT and the coherent gain of the # window(note that the correction is equal to func_length NOT fft_len) - func_fft = func_fft / (func_length * coherent_gain) + epsilon = 1e-10 # Small value to prevent division by zero + func_fft = func_fft / (func_length * coherent_gain + epsilon) # reduce to a single sided spectrum where the number of unique points for # even numbered FFT lengths is given by N / 2 + 1, and for odd(N + 1) / 2 num_unique_pts = int(np.ceil((fft_len + 1) / 2)) - if dim == 0: - func_fft = func_fft[0:num_unique_pts] - elif dim == 1: - func_fft = func_fft[:, 0:num_unique_pts] - elif dim == 2: - func_fft = func_fft[:, :, 0:num_unique_pts] - elif dim == 3: - func_fft = func_fft[:, :, :, 0:num_unique_pts] + slicing = [slice(None)] * len(sz) + slicing[dim] = slice(0, num_unique_pts) + func_fft = func_fft[tuple(slicing)] func_fft = single_sided_correction(func_fft, fft_len, dim) # create the frequency axis variable - f = np.arange(0, func_fft.shape[dim]) * Fs / fft_len + f = np.arange(0, num_unique_pts) * Fs / fft_len # calculate the amplitude spectrum func_as = np.abs(func_fft) @@ -156,7 +151,7 @@ def spect( # unwrap the phase spectrum if required if unwrap_phase: - func_ps = scipy.unwrap(func_ps, [], dim) + func_ps = np.unwrap(func_ps, axis=dim) return f, func_as, func_ps diff --git a/pyproject.toml b/pyproject.toml index 9ff02103..7c5ddd1b 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -30,7 +30,7 @@ dependencies = [ "deepdiff==8.0.1", "matplotlib==3.9.2", "numpy>=1.22.2,<1.27.0", - "beartype==0.18.5", + "beartype==0.19.0", "jaxtyping==0.2.34" ] @@ -42,17 +42,17 @@ Bug-tracker = "https://github.com/waltsims/k-wave-python/issues" [project.optional-dependencies] test = ["pytest", - "coverage==7.6.1", + "coverage==7.6.4", "phantominator", "testfixtures==8.3.0", "requests==2.32.3"] example = ["gdown==5.2.0"] docs = [ "sphinx-mdinclude==0.6.2", "sphinx-copybutton==0.5.2", - "sphinx-tabs==3.4.5", + "sphinx-tabs==3.4.7", "sphinx-toolbox==3.8.0", "furo==2024.8.6"] -dev = ["pre-commit==3.8.0"] +dev = ["pre-commit==4.0.1"] [tool.hatch.version] path = "kwave/__init__.py" diff --git a/tests/test_executor.py b/tests/test_executor.py index e0a6b209..f65f37d8 100644 --- a/tests/test_executor.py +++ b/tests/test_executor.py @@ -99,9 +99,11 @@ def test_run_simulation_success(self): sensor_data = executor.run_simulation("input.h5", "output.h5", "options") normalized_path = os.path.normpath(self.execution_options.binary_path) - expected_command = f"{self.execution_options.system_string} " f"{normalized_path} " f"-i input.h5 " f"-o output.h5 " f"options" + expected_command = [normalized_path, "-i", "input.h5", "-o", "output.h5", "options"] - self.mock_popen.assert_called_once_with(expected_command, stdout=subprocess.PIPE, stderr=subprocess.PIPE, shell=True, text=True) + self.mock_popen.assert_called_once_with( + expected_command, env=self.execution_options.env_vars, stdout=subprocess.PIPE, stderr=subprocess.PIPE, text=True + ) self.mock_proc.communicate.assert_called_once() self.assertEqual(sensor_data, dotdict()) diff --git a/tests/test_simulation_execution_options.py b/tests/test_simulation_execution_options.py new file mode 100644 index 00000000..2da09e72 --- /dev/null +++ b/tests/test_simulation_execution_options.py @@ -0,0 +1,71 @@ +import unittest +from kwave.options.simulation_execution_options import SimulationExecutionOptions +from kwave.ksensor import kSensor +from unittest.mock import patch + + +class TestSimulationExecutionOptions(unittest.TestCase): + def setUp(self): + # Set up a default kSensor object for testing + self.sensor = kSensor() + self.sensor.record = ["p", "u"] + self.sensor.record_start_index = 10 + + def test_default_initialization(self): + with patch("kwave.options.simulation_execution_options.PLATFORM", "linux"): + options = SimulationExecutionOptions() + self.assertFalse(options.is_gpu_simulation) + self.assertEqual(options.binary_name, "kspaceFirstOrder-OMP") + self.assertTrue(options.delete_data) + self.assertEqual(options.num_threads, None) # TODO: confusing logic here + self.assertEqual(options.verbose_level, 0) + + def test_gpu_simulation_initialization(self): + with patch("kwave.options.simulation_execution_options.PLATFORM", "linux"): + options = SimulationExecutionOptions(is_gpu_simulation=True) + self.assertTrue(options.is_gpu_simulation) + self.assertEqual(options.binary_name, "kspaceFirstOrder-CUDA") + + def test_binary_name_extension_on_windows(self): + with patch("kwave.options.simulation_execution_options.PLATFORM", "windows"): + options = SimulationExecutionOptions() + self.assertTrue(options.binary_name.endswith(".exe")) + + def test_get_options_string(self): + options = SimulationExecutionOptions() + options_string = options.get_options_string(self.sensor) + self.assertIn("--p_raw", options_string) + self.assertIn("--u_raw", options_string) + self.assertIn("-s 10", options_string) + + def test_env_vars_linux(self): + with patch("kwave.options.simulation_execution_options.PLATFORM", "linux"): + options = SimulationExecutionOptions() + env_vars = options.env_vars + self.assertIn("OMP_PLACES", env_vars) + self.assertEqual(env_vars["OMP_PLACES"], "cores") + self.assertIn("OMP_PROC_BIND", env_vars) + self.assertEqual(env_vars["OMP_PROC_BIND"], "SPREAD") + + def test_thread_binding_linux(self): + with patch("kwave.options.simulation_execution_options.PLATFORM", "linux"): + options = SimulationExecutionOptions(thread_binding=True) + env_vars = options.env_vars + self.assertEqual(env_vars["OMP_PROC_BIND"], "SPREAD") + + def test_thread_binding_darwin(self): + with patch("kwave.options.simulation_execution_options.PLATFORM", "darwin"): + options = SimulationExecutionOptions(thread_binding=True) + with self.assertRaises(ValueError, msg="Thread binding is not supported in MacOS."): + _ = options.env_vars + + def test_env_vars_darwin(self): + with patch("kwave.options.simulation_execution_options.PLATFORM", "darwin"): + options = SimulationExecutionOptions() + env_vars = options.env_vars + self.assertNotIn("OMP_PLACES", env_vars) + self.assertNotIn("OMP_PROC_BIND", env_vars) + + +if __name__ == "__main__": + unittest.main()