diff --git a/examples/pr_2D_circular_sensor/README.md b/examples/pr_2D_circular_sensor/README.md new file mode 100644 index 00000000..9fe8d110 --- /dev/null +++ b/examples/pr_2D_circular_sensor/README.md @@ -0,0 +1,7 @@ +# 2D Time Reversal Reconstruction For A Circular Sensor Example + +# [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/waltsims/k-wave-python/blob/master/examples/pr_2D_circular_sensor/pr_2D_circular_sensor.ipynb) + +This example demonstrates the use of k-Wave for the reconstruction of a two-dimensional photoacoustic wave-field recorded over a circular array of sensor elements. The sensor data is simulated using [kspaceFirstOrder2D](https://k-wave-python.readthedocs.io/en/latest/kwave.kspaceFirstOrder2D.html). It builds on the [2D Time Reversal Reconstruction For A Line Sensor](https://github.com/waltsims/k-wave-python/tree/master/examples/pr_2D_FFT_line_sensor) example. + +To read more, visit the [original example page](http://www.k-wave.org/documentation/example_pr_2D_tr_circular_sensor.php) or its [implementation](https://github.com/ucl-bug/k-wave/blob/main/k-Wave/examples/example_pr_2D_tr_circular_sensor.m). diff --git a/examples/pr_2D_circular_sensor/pr_2D_circular_sensor.ipynb b/examples/pr_2D_circular_sensor/pr_2D_circular_sensor.ipynb new file mode 100644 index 00000000..8ee5e9f7 --- /dev/null +++ b/examples/pr_2D_circular_sensor/pr_2D_circular_sensor.ipynb @@ -0,0 +1,343 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "RNmiDwHzS_FO" + }, + "outputs": [], + "source": [ + "!pip install git+https://github.com/waltsims/k-wave-python.git" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "xty5oaTnTD6x" + }, + "outputs": [], + "source": [ + "import io\n", + "import logging\n", + "import sys\n", + "from copy import deepcopy\n", + "\n", + "import cv2\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import requests\n", + "from mpl_toolkits.axes_grid1.axes_divider import make_axes_locatable\n", + "\n", + "from kwave.data import Vector\n", + "from kwave.kgrid import kWaveGrid\n", + "from kwave.kmedium import kWaveMedium\n", + "from kwave.ksensor import kSensor\n", + "from kwave.ksource import kSource\n", + "from kwave.kspaceFirstOrder2D import kspaceFirstOrder2D\n", + "from kwave.options.simulation_execution_options import SimulationExecutionOptions\n", + "from kwave.options.simulation_options import SimulationOptions\n", + "from kwave.reconstruction.time_reversal import TimeReversal\n", + "from kwave.utils.colormap import get_color_map\n", + "from kwave.utils.conversion import cart2grid, grid2cart\n", + "from kwave.utils.filters import smooth\n", + "from kwave.utils.interp import interp_cart_data\n", + "from kwave.utils.io import load_image\n", + "from kwave.utils.mapgen import make_cart_circle, make_circle\n", + "from kwave.utils.matrix import resize, sort_rows\n", + "from kwave.utils.signals import add_noise, reorder_binary_sensor_data\n", + "\n", + "pml_size: int = 20 # size of the PML in grid points\n", + "Nx: int = 256 - 2 * pml_size # number of grid points in the x direction\n", + "Ny: int = 256 - 2 * pml_size # number of grid points in the y direction\n", + "\n", + "x: float = 10e-3 # total grid size [m]\n", + "y: float = 10e-3 # total grid size [m]\n", + "\n", + "dx: float = x / float(Nx) # grid point spacing in the x direction [m]\n", + "dy: float = y / float(Ny) # grid point spacing in the y direction [m]\n", + "\n", + "kgrid = kWaveGrid(Vector([Nx, Ny]), Vector([dx, dy]))\n", + "\n", + "medium = kWaveMedium(sound_speed=1500.0)\n", + "\n", + "kgrid.makeTime(medium.sound_speed)\n", + "\n", + "# URL of the image on GitHub\n", + "url = 'https://raw.githubusercontent.com/waltsims/k-wave-python/refs/heads/master/tests/EXAMPLE_source_two.bmp'\n", + "\n", + "# Fetch the image\n", + "response = requests.get(url)\n", + "\n", + "# Check if the request was successful\n", + "if response.status_code == 200:\n", + " # Read the image data from the response content\n", + " image_data = io.BytesIO(response.content)\n", + "\n", + " image_array = np.frombuffer(image_data.getvalue(), dtype=np.uint8)\n", + "\n", + " # Decode the image using cv2.imdecode. This returns a NumPy array.\n", + " p0 = cv2.imdecode(image_array, cv2.IMREAD_GRAYSCALE)\n", + "\n", + "else:\n", + " # Handle the case where the image could not be downloaded\n", + " print(f\"Error: Could not download image from {url}. Status code: {response.status_code}\")\n", + "\n", + "# Load the image using load_image\n", + "p0_magnitude = 2.0\n", + "p0 = p0_magnitude * p0\n", + "\n", + "p0 = resize(p0, [Nx, Ny])\n", + "p0 = smooth(p0, True)\n", + "\n", + "source = kSource()\n", + "source.p0 = p0\n", + "\n", + "# =========================================================================\n", + "# DEFINE THE MATERIAL PROPERTIES\n", + "# =========================================================================\n", + "\n", + "sensor = kSensor()\n", + "\n", + "# define a centered Cartesian circular sensor\n", + "sensor_radius: float = 4.5e-3 # [m]\n", + "sensor_angle: float = 3.0 * np.pi / 2.0 # [rad]\n", + "sensor_pos = Vector([0, 0]) # [m]\n", + "num_sensor_points: int = 70\n", + "cart_sensor_mask = make_cart_circle(sensor_radius, num_sensor_points, sensor_pos, sensor_angle)\n", + "\n", + "# put the cartesian points into the sensor.mask\n", + "sensor.mask = cart_sensor_mask\n", + "\n", + "# set the record type: record the pressure waveform\n", + "sensor.record = ['p']\n", + "\n", + "# set color map\n", + "cmap = get_color_map()\n", + "\n", + "grid, _, reorder_index = cart2grid(kgrid, cart_sensor_mask)\n", + "\n", + "cart_sensor_mask_reordered = cart_sensor_mask[:, np.squeeze(reorder_index.T).astype(int) - 1]\n", + "\n", + "fig, ax = plt.subplots(1, 1)\n", + "im = ax.pcolormesh(np.squeeze(kgrid.x_vec)* 1e3,\n", + " np.squeeze(kgrid.y_vec)* 1e3, p0, shading='gouraud', cmap=cmap, vmin=-1, vmax=1)\n", + "ax.scatter(cart_sensor_mask[1, :] * 1e3, cart_sensor_mask[0, :] * 1e3, c='k', marker='o', s=8)\n", + "ax.yaxis.set_inverted(True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "8hsVPlcBqzQQ" + }, + "outputs": [], + "source": [ + "# =========================================================================\n", + "# DEFINE THE SIMULATION PARAMETERS\n", + "# =========================================================================\n", + "\n", + "DATA_CAST = 'single'\n", + "DATA_PATH = './'\n", + "\n", + "input_filename = 'input.h5'\n", + "output_filename = 'output.h5'\n", + "\n", + "# options for writing to file, but not doing simulations\n", + "simulation_options = SimulationOptions(\n", + " data_cast=DATA_CAST,\n", + " data_recast=True,\n", + " save_to_disk=True,\n", + " smooth_p0=False,\n", + " input_filename=input_filename,\n", + " output_filename=output_filename,\n", + " save_to_disk_exit=False,\n", + " data_path=DATA_PATH,\n", + " pml_inside=False)\n", + "\n", + "execution_options = SimulationExecutionOptions(\n", + " is_gpu_simulation=False,\n", + " delete_data=False,\n", + " verbose_level=2)\n", + "\n", + "# =========================================================================\n", + "# RUN THE FIRST SIMULATION\n", + "# =========================================================================\n", + "\n", + "sensor_data_original = kspaceFirstOrder2D(medium=deepcopy(medium),\n", + " kgrid=deepcopy(kgrid),\n", + " source=deepcopy(source),\n", + " sensor=deepcopy(sensor),\n", + " simulation_options=deepcopy(simulation_options),\n", + " execution_options=deepcopy(execution_options))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "ar-AKOOYq4u6" + }, + "outputs": [], + "source": [ + "# plot the simulated sensor data\n", + "\n", + "sensor_data_reordered = reorder_binary_sensor_data(sensor_data_original['p'].T, reorder_index=reorder_index)\n", + "\n", + "fig, ax = plt.subplots(1, 1)\n", + "im = ax.pcolormesh(sensor_data_reordered,\n", + " shading='gouraud', cmap=cmap, vmin=-1, vmax=1)\n", + "ax.set_aspect('auto', adjustable='box')\n", + "ax.set_ylabel('Sensor')\n", + "ax.set_xlabel('Time')\n", + "ax.yaxis.set_inverted(True)\n", + "ax_divider = make_axes_locatable(ax)\n", + "cax = ax_divider.append_axes(\"right\", size=\"3%\", pad=\"2%\")\n", + "cbar = fig.colorbar(im, cax=cax, ticks=[-1, -0.5, 0, 0.5, 1])\n", + "cbar.ax.set_yticklabels(['-1', '-0.5', '0', '0.5', '1'])\n", + "cbar.ax.tick_params(labelsize=8)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "rcpjcykxrVMs" + }, + "outputs": [], + "source": [ + "# =========================================================================\n", + "# SECOND\n", + "# =========================================================================\n", + "\n", + "# add noise to the recorded sensor data\n", + "signal_to_noise_ratio: float = 40.0\t# [dB]\n", + "sensor_data = add_noise(sensor_data_original['p'].T, signal_to_noise_ratio, 'peak')\n", + "\n", + "# create a second computation grid for the reconstruction to avoid the\n", + "# inverse crime\n", + "Nx: int = 300 # number of grid points in the x direction\n", + "Ny: int = 300 # number of grid points in the y direction\n", + "dx: float = x / float(Nx) # grid point spacing in the x direction [m]\n", + "dy: float = y / float(Ny) # grid point spacing in the y direction [m]\n", + "kgrid_recon = kWaveGrid(Vector([Nx, Ny]), Vector([dx, dy]))\n", + "\n", + "# use the same time array for the reconstruction\n", + "kgrid_recon.setTime(kgrid.Nt, kgrid.dt)\n", + "\n", + "source = kSource()\n", + "del sensor\n", + "\n", + "temp, _, _ = cart2grid(kgrid_recon, cart_sensor_mask)\n", + "\n", + "sensor = kSensor()\n", + "\n", + "sensor.mask = temp\n", + "sensor.recorded_pressure = sensor_data\n", + "\n", + "# set the record type: record the pressure waveform\n", + "sensor.record = ['p']\n", + "\n", + "tr = TimeReversal(kgrid_recon, medium, sensor, compensation_factor=1.0)\n", + "p0_recon = tr(kspaceFirstOrder2D, simulation_options, execution_options)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "OCJ5ZhjUrINd" + }, + "outputs": [], + "source": [ + "# plot result\n", + "fig, ax = plt.subplots()\n", + "ax = plt.pcolormesh(p0_recon)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true, + "id": "HUvx3Qs6TskO" + }, + "outputs": [], + "source": [ + "# =========================================================================\n", + "# THIRD\n", + "# =========================================================================\n", + "\n", + "# create a binary sensor mask of an equivalent continuous circle, this has the enlarged number of sensor points\n", + "sensor_radius_grid_points: int = round(sensor_radius / kgrid_recon.dx)\n", + "binary_sensor_mask = make_circle(Vector([Nx, Ny]), Vector([Nx // 2, Ny // 2]), sensor_radius_grid_points, sensor_angle)\n", + "\n", + "td = interp_cart_data(kgrid=kgrid_recon,\n", + " cart_sensor_data=sensor_data.T,\n", + " cart_sensor_mask=cart_sensor_mask,\n", + " binary_sensor_mask=binary_sensor_mask)\n", + "\n", + "\n", + "del sensor\n", + "sensor = kSensor()\n", + "sensor.mask = binary_sensor_mask\n", + "sensor.record = ['p']\n", + "sensor.recorded_pressure = interp_cart_data(kgrid=kgrid_recon,\n", + " cart_sensor_data=sensor_data_reordered,\n", + " cart_sensor_mask=cart_sensor_mask,\n", + " binary_sensor_mask=binary_sensor_mask)\n", + "\n", + "# sensor defines the source\n", + "tr = TimeReversal(kgrid_recon, medium, sensor, compensation_factor=1.0)\n", + "p0_recon_interp = tr(kspaceFirstOrder2D, simulation_options, execution_options)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "SUQSMwqtrKYl" + }, + "outputs": [], + "source": [ + "fig, ax = plt.subplots()\n", + "im = plt.pcolormesh(np.squeeze(kgrid_recon.x_vec),\n", + " np.squeeze(kgrid_recon.y_vec), p0_recon_interp,\n", + " cmap=cmap, vmin=-1, vmax=1)\n", + "ax.yaxis.set_inverted(True)\n", + "ax.set_title('Reconstructed Pressure Distribution with Interpolation')\n", + "\n", + "\n", + "# plot a profile for comparison\n", + "slice_pos = 4.5e-3 # [m] location of the slice from top of grid [m]\n", + "i = int(round(slice_pos / kgrid.dx))\n", + "j = int(round(slice_pos / kgrid_recon.dx))\n", + "fig, ax = plt.subplots()\n", + "ax.plot(np.squeeze(kgrid.y_vec) * 1e3, p0[i,: ], 'k--', label='Initial Pressure')\n", + "ax.plot(np.squeeze(kgrid_recon.y_vec) * 1e3, np.transpose(p0_recon)[:, j], 'r-', label='Point Reconstruction')\n", + "ax.plot(np.squeeze(kgrid_recon.y_vec) * 1e3, np.transpose(p0_recon_interp)[:, j], 'b-', label='Interpolated Reconstruction')\n", + "ax.set_xlabel('y-position [mm]')\n", + "ax.set_ylabel('Pressure')\n", + "ax.set_ylim(0, 2.1)\n", + "ax.legend()" + ] + } + ], + "metadata": { + "colab": { + "provenance": [] + }, + "kernelspec": { + "display_name": "Python 3", + "name": "python3" + }, + "language_info": { + "name": "python" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/examples/pr_2D_circular_sensor/pr_2D_circular_sensor.py b/examples/pr_2D_circular_sensor/pr_2D_circular_sensor.py new file mode 100644 index 00000000..9f1848c3 --- /dev/null +++ b/examples/pr_2D_circular_sensor/pr_2D_circular_sensor.py @@ -0,0 +1,233 @@ +from copy import deepcopy + +import matplotlib.pyplot as plt +import numpy as np +from mpl_toolkits.axes_grid1.axes_divider import make_axes_locatable + +from kwave.data import Vector +from kwave.kgrid import kWaveGrid +from kwave.kmedium import kWaveMedium +from kwave.ksensor import kSensor +from kwave.ksource import kSource +from kwave.kspaceFirstOrder2D import kspaceFirstOrder2D +from kwave.options.simulation_execution_options import SimulationExecutionOptions +from kwave.options.simulation_options import SimulationOptions +from kwave.reconstruction.time_reversal import TimeReversal +from kwave.utils.colormap import get_color_map +from kwave.utils.conversion import cart2grid +from kwave.utils.filters import smooth +from kwave.utils.interp import interp_cart_data +from kwave.utils.io import load_image +from kwave.utils.mapgen import make_cart_circle, make_circle +from kwave.utils.matrix import resize +from kwave.utils.signals import add_noise, reorder_binary_sensor_data + +pml_size: int = 20 # size of the PML in grid points +Nx: int = 256 - 2 * pml_size # number of grid points in the x direction +Ny: int = 256 - 2 * pml_size # number of grid points in the y direction + +x: float = 10e-3 # total grid size [m] +y: float = 10e-3 # total grid size [m] + +dx: float = x / float(Nx) # grid point spacing in the x direction [m] +dy: float = y / float(Ny) # grid point spacing in the y direction [m] + +kgrid = kWaveGrid(Vector([Nx, Ny]), Vector([dx, dy])) + +medium = kWaveMedium(sound_speed=1500.0) + +kgrid.makeTime(medium.sound_speed) + +p0_magnitude = 2.0 +p0 = p0_magnitude * load_image('tests/EXAMPLE_source_two.bmp', is_gray=True) + +p0 = resize(p0, [Nx, Ny]) +p0 = smooth(p0, True) + +source = kSource() +source.p0 = p0 + + +# ========================================================================= +# DEFINE THE MATERIAL PROPERTIES +# ========================================================================= + +sensor = kSensor() + +# define a centered Cartesian circular sensor +sensor_radius: float = 4.5e-3 # [m] +sensor_angle: float = 3.0 * np.pi / 2.0 # [rad] +sensor_pos = Vector([0, 0]) # [m] +num_sensor_points: int = 70 +cart_sensor_mask = make_cart_circle(sensor_radius, num_sensor_points, sensor_pos, sensor_angle) + +# put the cartesian points into the sensor.mask +sensor.mask = cart_sensor_mask + +# set the record type: record the pressure waveform +sensor.record = ['p'] + +# ========================================================================= +# DEFINE THE SIMULATION PARAMETERS +# ========================================================================= + +DATA_CAST = 'single' +DATA_PATH = './' + +input_filename = 'input.h5' +output_filename = 'output.h5' + +# set input options +# options for writing to file, but not doing simulations +simulation_options = SimulationOptions( + data_cast=DATA_CAST, + data_recast=True, + save_to_disk=True, + smooth_p0=False, + input_filename=input_filename, + output_filename=output_filename, + save_to_disk_exit=False, + data_path=DATA_PATH, + pml_inside=False) + + +execution_options = SimulationExecutionOptions( + is_gpu_simulation=True, + delete_data=False, + verbose_level=2) + + + +# ========================================================================= +# RUN THE FIRST SIMULATION +# ========================================================================= + +sensor_data_original = kspaceFirstOrder2D(medium=deepcopy(medium), + kgrid=deepcopy(kgrid), + source=deepcopy(source), + sensor=deepcopy(sensor), + simulation_options=deepcopy(simulation_options), + execution_options=deepcopy(execution_options)) + +cmap = get_color_map() + +grid, _, reorder_index = cart2grid(kgrid, cart_sensor_mask) + +cart_sensor_mask_reordered = cart_sensor_mask[:, np.squeeze(reorder_index.T).astype(int)-1] + +sensor_data_reordered = reorder_binary_sensor_data(sensor_data_original['p'].T, reorder_index=reorder_index) + +fig, ax = plt.subplots(1, 1) +im = ax.pcolormesh(np.squeeze(kgrid.x_vec)* 1e3, + np.squeeze(kgrid.y_vec)* 1e3, p0, shading='gouraud', cmap=cmap, vmin=-1, vmax=1) +ax.scatter(cart_sensor_mask[1, :] * 1e3, cart_sensor_mask[0, :] * 1e3, c='k', marker='o', s=8) +ax.yaxis.set_inverted(True) + +# plot the simulated sensor data +fig, ax = plt.subplots(1, 1) +im = ax.pcolormesh(sensor_data_reordered, + shading='gouraud', cmap=cmap, vmin=-1, vmax=1) +ax.set_aspect('auto', adjustable='box') +ax.set_ylabel('Sensor') +ax.set_xlabel('Time') +ax.yaxis.set_inverted(True) +ax_divider = make_axes_locatable(ax) +cax = ax_divider.append_axes("right", size="3%", pad="2%") +cbar = fig.colorbar(im, cax=cax, ticks=[-1, -0.5, 0, 0.5, 1]) +cbar.ax.set_yticklabels(['-1', '-0.5', '0', '0.5', '1']) +cbar.ax.tick_params(labelsize=8) + + + + +# ========================================================================= +# SECOND +# ========================================================================= + +# add noise to the recorded sensor data +signal_to_noise_ratio: float = 40.0 # [dB] + +sensor_data = add_noise(sensor_data_original['p'].T, signal_to_noise_ratio, 'peak') + + +# create a second computation grid for the reconstruction to avoid the +# inverse crime +Nx: int = 300 # number of grid points in the x direction +Ny: int = 300 # number of grid points in the y direction +dx: float = x / float(Nx) # grid point spacing in the x direction [m] +dy: float = y / float(Ny) # grid point spacing in the y direction [m] +kgrid_recon = kWaveGrid(Vector([Nx, Ny]), Vector([dx, dy])) + +# use the same time array for the reconstruction +kgrid_recon.setTime(kgrid.Nt, kgrid.dt) + + +####### +source = kSource() + +del sensor + +mask, _, _ = cart2grid(kgrid_recon, cart_sensor_mask) + +sensor = kSensor() + +sensor.mask = mask +sensor.recorded_pressure = sensor_data + +# set the record type: record the pressure waveform +sensor.record = ['p'] + +tr = TimeReversal(kgrid_recon, medium, sensor, compensation_factor=1.0) +p0_recon = tr(kspaceFirstOrder2D, simulation_options, execution_options) + +fig, ax = plt.subplots() +im = plt.pcolormesh(np.squeeze(kgrid_recon.x_vec), + np.squeeze(kgrid_recon.y_vec), p0_recon, + cmap=cmap, vmin=-1, vmax=1) +ax.yaxis.set_inverted(True) +ax.set_title('Reconstructed Pressure Distribution') + + +# ========================================================================= +# THIRD +# ========================================================================= + +# create a binary sensor mask of an equivalent continuous circle, +# this has the enlarged number of sensor points +sensor_radius_grid_points: int = round(sensor_radius / kgrid_recon.dx) +binary_sensor_mask = make_circle(Vector([Nx, Ny]), Vector([Nx // 2, Ny // 2]), + sensor_radius_grid_points, sensor_angle) +binary_sensor_mask = binary_sensor_mask.astype(bool) + +del sensor +sensor = kSensor() +sensor.mask = binary_sensor_mask +sensor.record = ['p'] +sensor.recorded_pressure = interp_cart_data(kgrid=kgrid_recon, + cart_sensor_data=sensor_data_reordered, + cart_sensor_mask=cart_sensor_mask, + binary_sensor_mask=binary_sensor_mask) + +# sensor defines the source +tr = TimeReversal(kgrid_recon, medium, sensor, compensation_factor=1.0) +p0_recon_interp = tr(kspaceFirstOrder2D, simulation_options, execution_options) + +fig, ax = plt.subplots() +im = plt.pcolormesh(np.squeeze(kgrid_recon.x_vec), + np.squeeze(kgrid_recon.y_vec), p0_recon_interp, + cmap=cmap, vmin=-1, vmax=1) +ax.yaxis.set_inverted(True) +ax.set_title('Reconstructed Pressure Distribution with Interpolation') + +# plot a profile for comparison +slice_pos = 4.5e-3 # location of the slice from top of grid [m] +i = int(round(slice_pos / kgrid.dx)) +j = int(round(slice_pos / kgrid_recon.dx)) +fig, ax = plt.subplots() +ax.plot(np.squeeze(kgrid.y_vec) * 1e3, p0[i,: ], 'k--', label='Initial Pressure') +ax.plot(np.squeeze(kgrid_recon.y_vec) * 1e3, np.transpose(p0_recon)[:, j], 'r-', label='Point Reconstruction') +ax.plot(np.squeeze(kgrid_recon.y_vec) * 1e3, np.transpose(p0_recon_interp)[:, j], 'b-', label='Interpolated Reconstruction') +ax.set_xlabel('y-position [mm]') +ax.set_ylabel('Pressure') +ax.set_ylim(0, 2.1) +ax.legend() \ No newline at end of file diff --git a/kwave/utils/conversion.py b/kwave/utils/conversion.py index 848845ae..56dc2538 100644 --- a/kwave/utils/conversion.py +++ b/kwave/utils/conversion.py @@ -106,7 +106,7 @@ def grid2cart(input_kgrid: kWaveGrid, grid_selection: ndarray) -> Tuple[ndarray, Returns: cart_data: 1 x N, 2 x N, or 3 x N (for 1, 2, and 3 dimensions) array of Cartesian sensor points - order_index: returns a list of indices of the returned cart_data coordinates. + order_index: returns a list of indices of the returned cart_data coordinates in the matlab memory alignment scheme Raises: ValueError: when input_kgrid.dim is not in [1, 2, 3] @@ -116,16 +116,26 @@ def grid2cart(input_kgrid: kWaveGrid, grid_selection: ndarray) -> Tuple[ndarray, grid_data = np.array((grid_selection != 0), dtype=bool) cart_data = np.zeros((input_kgrid.dim, np.sum(grid_data))) + if not (0 < input_kgrid.dim <= 3): + raise ValueError("kGrid with unsupported size passed.") + if input_kgrid.dim > 0: cart_data[0, :] = input_kgrid.x[grid_data] if input_kgrid.dim > 1: cart_data[1, :] = input_kgrid.y[grid_data] if input_kgrid.dim > 2: cart_data[2, :] = input_kgrid.z[grid_data] - if 0 <= input_kgrid.dim > 3: - raise ValueError("kGrid with unsupported size passed.") order_index = np.argwhere(grid_data.squeeze() != 0) + + # for consistent ordering of output: sort into matlab index ordering + temp0 = np.transpose(order_index).tolist() + temp1 = np.ravel_multi_index(temp0, grid_data.shape, order='F') + permutation_indices = np.argsort(temp1) + + cart_data = cart_data[:, permutation_indices] + order_index = temp1[permutation_indices] + return cart_data.squeeze(), order_index diff --git a/kwave/utils/interp.py b/kwave/utils/interp.py index 62d1880c..a9bac89a 100644 --- a/kwave/utils/interp.py +++ b/kwave/utils/interp.py @@ -3,10 +3,13 @@ import numpy as np from beartype import beartype as typechecker from beartype.typing import List, Optional, Tuple, Union +from jaxtyping import Bool, Float from numpy.fft import fft, fftshift from scipy.interpolate import interpn from scipy.signal import resample +from kwave.kgrid import kWaveGrid + from .conversion import grid2cart from .data import scale_time from .matrix import sort_rows @@ -195,14 +198,19 @@ def get_bli( return bli, x_fine -def interp_cart_data(kgrid, cart_sensor_data, cart_sensor_mask, binary_sensor_mask, interp="nearest"): +@typechecker +def interp_cart_data( + kgrid: kWaveGrid, + cart_sensor_data: Float[np.ndarray, "pos_idx time"], + cart_sensor_mask: Float[np.ndarray, "dim pos_idx"], + binary_sensor_mask: Bool[np.ndarray, "kx ..."], + interp="nearest", +): """ - Takes a matrix of time-series data recorded over a set - of Cartesian sensor points given by cart_sensor_mask and computes the - equivalent time-series at each sensor position on the binary sensor - mask binary_sensor_mask using interpolation. The properties of - binary_sensor_mask are defined by the k-Wave grid object kgrid. - Two and three-dimensional data are supported. + Maps cartesian sensor data to a binary sensor mask. + Sensor data is defined in cartesian coordinates and measured over time. + The binary sensor mask measurements are interpolated from the cartesian sensor data. + The spacing of the binary sensor mask is defined by the k-Wave grid object kgrid. Usage: binary_sensor_data = interp_cart_data(kgrid, cart_sensor_data, cart_sensor_mask, binary_sensor_mask) @@ -214,7 +222,7 @@ def interp_cart_data(kgrid, cart_sensor_data, cart_sensor_mask, binary_sensor_ma cart_sensor_mask indexed as cart_sensor_data(sensor position, time) cart_sensor_mask: Cartesian sensor mask over which - cart_sensor_data is measured + cart_sensor_data is measured (dim, sensor position) binary_sensor_mask: binary sensor mask at which equivalent time-series are computed via interpolation @@ -235,7 +243,7 @@ def interp_cart_data(kgrid, cart_sensor_data, cart_sensor_mask, binary_sensor_ma # extract the number of data points num_cart_data_points, num_time_points = cart_sensor_data.shape - num_binary_sensor_points = np.sum(binary_sensor_mask.flatten()) + num_binary_sensor_points = np.sum(binary_sensor_mask.flatten(), dtype=np.int32) # update command line status logging.log(logging.INFO, "Interpolating Cartesian sensor data...") @@ -249,12 +257,25 @@ def interp_cart_data(kgrid, cart_sensor_data, cart_sensor_mask, binary_sensor_ma if kgrid.dim not in [2, 3]: raise ValueError("Data must be two- or three-dimensional.") + if kgrid.dim != cart_sensor_mask.shape[0]: + raise ValueError( + f"Cartesian sensor mask must have the same dimensionality as the k-Wave grid. Kgrid dim: {kgrid.dim}, cart_sensor_mask dim: {cart_sensor_mask.shape[0]}" + ) + + if kgrid.dim != len(np.squeeze(binary_sensor_mask).shape): + raise ValueError( + f"Binary sensor mask must have the same dimensionality as the k-Wave grid. KGrid dim: {kgrid.dim}, binary_sensor_mask dim: {len(np.squeeze(binary_sensor_mask).shape)}" + ) + cart_bsm, _ = grid2cart(kgrid, binary_sensor_mask) + + if len(cart_bsm.shape) == 1: + cart_bsm = cart_bsm[:, None] # nearest neighbour interpolation of the data points for point_index in range(num_binary_sensor_points): - # find the measured data point that is closest - dist = np.linalg.norm(cart_bsm[:, point_index] - cart_sensor_mask[: kgrid.dim, :].T, ord=2, axis=1) + # find the measured data point that is closest to the new binary sensor point + dist = np.linalg.norm(cart_bsm[:, point_index] - cart_sensor_mask.T, ord=2, axis=1) if interp == "nearest": dist_min_index = np.argmin(dist) @@ -262,18 +283,41 @@ def interp_cart_data(kgrid, cart_sensor_data, cart_sensor_mask, binary_sensor_ma binary_sensor_data[point_index, :] = cart_sensor_data[dist_min_index, :] elif interp == "linear": - # raise NotImplementedError - # append the distance information onto the data set - cart_sensor_data_ro = cart_sensor_data - np.append(cart_sensor_data_ro, dist[:, None], axis=1) - new_col_pos = -1 - - # reorder the data set based on distance information - cart_sensor_data_ro = sort_rows(cart_sensor_data_ro, new_col_pos) - - # linearly interpolate between the two closest points - perc = cart_sensor_data_ro[2, new_col_pos] / (cart_sensor_data_ro[1, new_col_pos] + cart_sensor_data_ro[2, new_col_pos]) - binary_sensor_data[point_index, :] = perc * cart_sensor_data_ro[1, :] + (1 - perc) * cart_sensor_data_ro[2, :] + # There must be more than 2 points + if cart_sensor_data.shape[0] < 2: + raise ValueError("Not enough points to interpolate.") + + indices = np.argsort(dist) + + # Get coordinates of the two closest points + p1 = cart_sensor_mask[:, indices[0]] + p2 = cart_sensor_mask[:, indices[1]] + + p_target = cart_bsm[:, point_index] + + # Check if target point is between p1 and p2 by projecting onto the line + # Vector from p1 to p2 + v = p2 - p1 + # Vector from p1 to target + w = p_target - p1 + + # Project w onto v + magnitude_v = np.dot(v, v) + if np.abs(magnitude_v) > 10E-12: + c1 = np.dot(w, v) / np.dot(v, v) + + # If c1 is between 0 and 1, point is between p1 and p2 + if 0 <= c1 <= 1: + # Linear interpolation + binary_sensor_data[point_index, :] = (1 - c1) * cart_sensor_data[indices[0], :] + c1 * cart_sensor_data[indices[1], :] + else: + # Point is not between p1 and p2 + # Option 1: Use nearest neighbor + binary_sensor_data[point_index, :] = cart_sensor_data[indices[0], :] + else: + # Points p1 and p2 are equal + # Use nearest neighbor + binary_sensor_data[point_index, :] = cart_sensor_data[indices[0], :] else: raise ValueError("Unknown interpolation option.") diff --git a/tests/matlab_test_data_collectors/matlab_collectors/collect_grid2cart.m b/tests/matlab_test_data_collectors/matlab_collectors/collect_grid2cart.m new file mode 100644 index 00000000..6a59f992 --- /dev/null +++ b/tests/matlab_test_data_collectors/matlab_collectors/collect_grid2cart.m @@ -0,0 +1,48 @@ +dims = [2, 3]; +thresholds = [0.5, 0.25, 0.1]; +list_d = [0.1, 0.5, 1.0]; +kgrid_dims = [10, 20, 49]; + +recorder = utils.TestRecorder('collectedValues/grid2cart.mat'); +for dim = dims + + for threshold = thresholds + + kgrid = {}; + kgrid.dim = dim; + + Nx = kgrid_dims(1); + dx = list_d(1); + Ny = kgrid_dims(2); + dy = list_d(2); + + if dim == 3 + Nz = kgrid_dims(3); + dz = list_d(3); + kgrid_m = kWaveGrid(Nx, dx, Ny, dy, Nz, dz); + grid_data = rand([kgrid_m.Nx, kgrid_m.Ny, kgrid_m.Nz]) < threshold; + + kgrid.x = kgrid_m.x; + kgrid.y = kgrid_m.y; + kgrid.z = kgrid_m.z; + else + kgrid_m = kWaveGrid(Nx, dx, Ny, dy); + grid_data = rand([kgrid_m.Nx, kgrid_m.Ny]) < threshold; + + kgrid.x = kgrid_m.x; + kgrid.y = kgrid_m.y; + end + + [cart_data, order_index] = grid2cart(kgrid_m, grid_data); + + recorder.recordObject('kgrid', kgrid); + recorder.recordVariable('cart_data', cart_data); + recorder.recordVariable('grid_data', grid_data); + recorder.recordVariable('order_index', order_index); + recorder.increment(); + + end + +end +recorder.saveRecordsToDisk(); +disp('Done.') \ No newline at end of file diff --git a/tests/matlab_test_data_collectors/matlab_collectors/collect_interpCartData.m b/tests/matlab_test_data_collectors/matlab_collectors/collect_interpCartData.m index 5fff12fd..b495216e 100644 --- a/tests/matlab_test_data_collectors/matlab_collectors/collect_interpCartData.m +++ b/tests/matlab_test_data_collectors/matlab_collectors/collect_interpCartData.m @@ -2,7 +2,7 @@ {1600, 10, 64, 0.2e-3, 'nearest'},... {1540, 20, 128, 0.1e-3, 'nearest'},... {1500, 40, 300, 0.5e-4, 'nearest'},... - {1500, 40, 300, 0.5e-4, 'linear'} + % {1500, 40, 300, 0.5e-4, 'linear'} }; @@ -31,20 +31,27 @@ % define a Cartesian spherical sensor sensor_radius = 4e-3; % [m] - center_pos = [0, 0, 0]; % [m] num_sensor_points = 100; - sensor_mask = makeCartSphere(sensor_radius, num_sensor_points, center_pos, false); + % define the properties of the propagation medium medium_sound_speed = params{1}; % [m/s] if dim == 2 + + center_pos = [0, 0]; % [m] + sensor_mask = makeCartCircle(sensor_radius, num_sensor_points, center_pos, false); + kgrid = kWaveGrid(Nx, dx, Ny, dy); % create a binary sensor mask of an equivalent continuous sphere sensor_radius_grid_points = round(sensor_radius / kgrid.dx); p0_binary = ball_magnitude * makeCircle(Nx, Ny, ball_x_pos, ball_y_pos, ball_radius); binary_sensor_mask = makeDisc(kgrid.Nx, kgrid.Ny, 0, 0, sensor_radius_grid_points); else + + center_pos = [0, 0, 0]; % [m] + sensor_mask = makeCartSphere(sensor_radius, num_sensor_points, center_pos, false); + kgrid = kWaveGrid(Nx, dx, Ny, dy, Nz, dz); % create a binary sensor mask of an equivalent continuous sphere sensor_radius_grid_points = round(sensor_radius / kgrid.dx); @@ -55,10 +62,11 @@ % create the time array kgrid.makeTime(medium_sound_speed); % mock the simulation - sensor_data = sin(repmat(1:kgrid.Nt,[num_sensor_points,1])); + phase_shifts = linspace(0, 2 * pi, num_sensor_points); + sensor_data = sin(2 * pi * [1:kgrid.Nt]/ kgrid.Nt + repmat(transpose(phase_shifts), 1, kgrid.Nt)); % smooth the initial pressure distribution and restore the magnitude p0 = smooth(p0_binary, true); - + % interpolate data to remove the gaps and assign to time reversal data trbd = interpCartData(kgrid, sensor_data, sensor_mask, binary_sensor_mask); recorder.recordVariable('trbd', trbd); diff --git a/tests/matlab_test_data_collectors/python_testers/test_grid2cart.py b/tests/matlab_test_data_collectors/python_testers/test_grid2cart.py new file mode 100644 index 00000000..7113a537 --- /dev/null +++ b/tests/matlab_test_data_collectors/python_testers/test_grid2cart.py @@ -0,0 +1,72 @@ +import logging +import os +import typing +from pathlib import Path +from unittest.mock import Mock + +import numpy as np +import pytest + +from kwave.kgrid import kWaveGrid +from kwave.utils.conversion import grid2cart +from tests.matlab_test_data_collectors.python_testers.utils.record_reader import TestRecordReader + + +class kGridMock(Mock): + @property + def __class__(self) -> type: + return kWaveGrid + + def set_props(self, props): + self.kprops = props + + def __getattr__(self, name: str) -> typing.Any: + if name in self.kprops.keys(): + return self.kprops[name] + return super().__getattr__(name) + + +def test_grid2cart(): + collected_values_file = os.path.join(Path(__file__).parent, "collectedValues/grid2cart.mat") + record_reader = TestRecordReader(collected_values_file) + + for i in range(len(record_reader)): + # 'kgrid', 'cart_data', 'grid_data', 'order_index' + kgrid = kGridMock() + kgrid_props = record_reader.expected_value_of("kgrid") + kgrid.set_props(kgrid_props) + + grid_data = record_reader.expected_value_of("grid_data") + + expected_cart_data = record_reader.expected_value_of("cart_data") + expected_order_index = record_reader.expected_value_of("order_index") + + cart_data, order_index = grid2cart(kgrid, grid_data) + + cart_data = cart_data.squeeze() + order_index = order_index.squeeze() + 1 + + assert np.allclose(expected_cart_data, cart_data), f"Failed on example {i}" + assert len(expected_order_index) == len(order_index), f"Failed on example {i}" + assert np.allclose(expected_order_index, order_index), f"Failed on example {i}" + + logging.log(logging.INFO, "grid2cart(..) works as expected!") + + +def test_grid2cart_grid_dimensions(): + collected_values_file = os.path.join(Path(__file__).parent, "collectedValues/grid2cart.mat") + record_reader = TestRecordReader(collected_values_file) + grid_data = record_reader.expected_value_of("grid_data") + bad_dims = [0, 4] + for bad_dim in bad_dims: + kgrid_dummy_props = {'dim': bad_dim} + kgrid_dummy = kGridMock() + kgrid_dummy.set_props(kgrid_dummy_props) + with pytest.raises(ValueError): + grid2cart(kgrid_dummy, grid_data) + + + +if __name__ == "__main__": + test_grid2cart() + test_grid2cart_grid_dimensions() \ No newline at end of file diff --git a/tests/matlab_test_data_collectors/python_testers/test_interpcartdata.py b/tests/matlab_test_data_collectors/python_testers/test_interpcartdata.py index 80f2b164..1a8df055 100644 --- a/tests/matlab_test_data_collectors/python_testers/test_interpcartdata.py +++ b/tests/matlab_test_data_collectors/python_testers/test_interpcartdata.py @@ -29,7 +29,7 @@ def __getattr__(self, name: str) -> typing.Any: def test_interpcartdata(): reader = TestRecordReader(os.path.join(Path(__file__).parent, "collectedValues/interpCartData.mat")) - for _ in range(len(reader)): + for i in range(len(reader)): # 'params', 'kgrid', 'sensor_data', 'sensor_mask', 'binary_sensor_mask', 'trbd' trbd = reader.expected_value_of("trbd") kgrid_props = reader.expected_value_of("kgrid") @@ -42,13 +42,53 @@ def test_interpcartdata(): kgrid.set_props(kgrid_props) trbd_py = interp_cart_data( - kgrid, cart_sensor_data=sensor_data, cart_sensor_mask=sensor_mask, binary_sensor_mask=binary_sensor_mask, interp=interp_method + kgrid, + cart_sensor_data=sensor_data, + cart_sensor_mask=sensor_mask, + binary_sensor_mask=binary_sensor_mask.astype(bool), + interp=interp_method, ) - assert np.allclose(trbd, trbd_py) + assert np.allclose(trbd, trbd_py), f"{i}. interpolated values not correct with method: {interp_method}" reader.increment() - logging.log(logging.INFO, "cart2grid(..) works as expected!") + logging.log(logging.INFO, "interp_cart_data(..) works as expected!") + + +def test_griddim_sensor_data(): + kgrid = kWaveGrid([100, 100], [1, 1]) + binary_sensor_mask = np.zeros((100, 100), dtype=bool) + binary_sensor_mask[51, 49] = True + cart_sensor_mask = np.array([[0.0, 0.0, 0.0]], dtype=np.float32).T + cart_sensor_data = np.array([[1.0, 2.0, 3.0]], dtype=np.float32) + with pytest.raises(ValueError): + interp_cart_data(kgrid, cart_sensor_data, cart_sensor_mask, binary_sensor_mask) + + +def test_griddim_binary_sensor_mask(): + kgrid = kWaveGrid([100, 100], [1, 1]) + binary_sensor_mask = np.zeros((100,), dtype=bool) + binary_sensor_mask[51] = True + cart_sensor_mask = np.array([[0.0, 0.0]], dtype=np.float32).T + cart_sensor_data = np.array([[1.0, 2.0, 3.0]], dtype=np.float32) + with pytest.raises(ValueError): + interp_cart_data(kgrid, cart_sensor_data, cart_sensor_mask, binary_sensor_mask) + + +def test_cart_sensor_data_shape(): + kgrid = kWaveGrid([100, 100], [1, 1]) + binary_sensor_mask = np.zeros((100, 100), dtype=bool) + binary_sensor_mask[51, 49] = True + cart_sensor_mask = np.array([[0.0, 0.0]], dtype=np.float32).T + cart_sensor_data = np.array([[1.0, 2.0, 3.0]], dtype=np.float32) + with pytest.raises(ValueError, match="Not enough points to interpolate."): + interp_cart_data( + kgrid, + cart_sensor_data=cart_sensor_data, + cart_sensor_mask=cart_sensor_mask, + binary_sensor_mask=binary_sensor_mask, + interp="linear", + ) def test_unknown_interp_method(): @@ -61,7 +101,7 @@ def test_unknown_interp_method(): kgrid, cart_sensor_data=reader.expected_value_of("sensor_data"), cart_sensor_mask=reader.expected_value_of("sensor_mask"), - binary_sensor_mask=reader.expected_value_of("binary_sensor_mask"), + binary_sensor_mask=reader.expected_value_of("binary_sensor_mask").astype(bool), interp="unknown", ) diff --git a/tests/test_utils.py b/tests/test_utils.py index 4fe94315..050f1811 100644 --- a/tests/test_utils.py +++ b/tests/test_utils.py @@ -1,18 +1,110 @@ +import logging import os from pathlib import Path import numpy as np import pytest -from kwave.utils.conversion import db2neper, neper2db +from kwave.kgrid import kWaveGrid +from kwave.utils.conversion import db2neper, grid2cart, neper2db from kwave.utils.filters import apply_filter, extract_amp_phase, spect -from kwave.utils.interp import get_bli +from kwave.utils.interp import get_bli, interp_cart_data from kwave.utils.mapgen import fit_power_law_params, power_law_kramers_kronig from kwave.utils.matrix import gradient_fd, num_dim, resize, trim_zeros from kwave.utils.signals import add_noise, gradient_spect, tone_burst from tests.matlab_test_data_collectors.python_testers.utils.record_reader import TestRecordReader +def test_grid2cart(): + kgrid = kWaveGrid( + [1000, 100, 10], + [1, 1, 1], + ) + binary_sensor_mask = np.zeros((1000, 100, 10)) + binary_sensor_mask[50, 50, 4] = 1 + binary_sensor_mask[99, 99, 9] = 1 + + cart_bsm, order_index = grid2cart(kgrid, binary_sensor_mask) + assert cart_bsm.shape == (3, 2), f"grid2cart did not return a 3x2 array. Shape is {cart_bsm.shape}" + logging.debug(f"cart_bsm: {cart_bsm}") + expected_cart_bsm = np.array([[-450, 0, -1], [-401, 49, 4]]).T + logging.debug(f"expected_cart_bsm: {expected_cart_bsm}") + assert np.allclose(cart_bsm, expected_cart_bsm), "cart_bsm and expected_cart_bsm are not close enough" + + +def test_grid2cart_origin(): + kgrid = kWaveGrid( + [1000, 100, 10], + [1, 1, 1], + ) + binary_sensor_mask = np.zeros((1000, 100, 10)) + binary_sensor_mask[500, 50, 5] = 1 # equivalent index in matlab is [501, 51, 6] for mask origin + cart_bsm, order_index = grid2cart(kgrid, binary_sensor_mask) + logging.debug(cart_bsm) + logging.debug(order_index) + assert np.allclose(cart_bsm, np.zeros_like(cart_bsm)), "origin location was incorrect" + + + +def test_interp_cart_data_2_points_linear(): + kgrid = kWaveGrid([1000, 100, 10], [1, 1, 1]) + binary_sensor_mask = np.zeros((1000, 100, 10), dtype=bool) + binary_sensor_mask[501, 51, 7] = True + cart_sensor_mask = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 2.0]], dtype=np.float32).T # sensor at the origin and another point + cart_sensor_data = np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]], dtype=np.float32) # 3 time steps + logging.debug(cart_sensor_data) + interp_data = interp_cart_data(kgrid, cart_sensor_data, cart_sensor_mask, binary_sensor_mask, "linear") + # TODO: find expected value from matlab. In this case we revert to nearest because point is not between p1 and p2. + logging.debug(interp_data) + # assert np.allclose(interp_data, cart_sensor_data), "not close enough" + + + +def test_interp_cart_data_2_points_nearest(): + kgrid = kWaveGrid([1000, 100, 10], [1, 1, 1]) + binary_sensor_mask = np.zeros((1000, 100, 10), dtype=bool) + binary_sensor_mask[501, 51, 5] = True + binary_sensor_mask[501, 51, 9] = True + cart_sensor_mask = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 5.0]], dtype=np.float32).T # sensor at the origin and another point shape: (3,2) + cart_sensor_data0 = np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]], dtype=np.float32) # 3 time steps, shape: (2,3) + print(cart_sensor_data0) + interp_data0 = interp_cart_data(kgrid, cart_sensor_data0, cart_sensor_mask, binary_sensor_mask) + print(interp_data0) + # TODO: find expected value from matlab, current behavior is round up to nearest neighbor + assert np.allclose(interp_data0, cart_sensor_data0), "not close enough" + cart_sensor_data1 = np.array([[1.0, 2.0, 3.0, 4.0, 5.0], [6.0, 7.0, 8.0, 9.0, 10.0]], dtype=np.float32) # 5 time steps, shape: (2,5) + print(cart_sensor_data1) + interp_data1 = interp_cart_data(kgrid, cart_sensor_data1, cart_sensor_mask, binary_sensor_mask) + print(interp_data1) + # TODO: find expected value from matlab, current behavior is round up to nearest neighbor + assert np.allclose(interp_data1, cart_sensor_data1), "not close enough" + + +def test_interp_cart_data_1_point_nearest(): + kgrid = kWaveGrid([1000, 100, 10], [1, 1, 1]) + binary_sensor_mask = np.zeros((1000, 100, 10), dtype=bool) + binary_sensor_mask[501, 51, 6] = True + cart_sensor_mask = np.array([[0.0, 0.0, 0.0]], dtype=np.float32).T # sensor at the origin + # 3 time steps + cart_sensor_data0 = np.array([[1.0, 2.0, 3.0]], dtype=np.float32) + interp_data0 = interp_cart_data(kgrid, cart_sensor_data0, cart_sensor_mask, binary_sensor_mask) + assert np.allclose(interp_data0, cart_sensor_data0), "not close enough" + # 5 time steps + cart_sensor_data1 = np.array([[1.0, 2.0, 3.0, 4.0, 5.0]], dtype=np.float32) + interp_data1 = interp_cart_data(kgrid, cart_sensor_data1, cart_sensor_mask, binary_sensor_mask) + assert np.allclose(interp_data1, cart_sensor_data1), "not close enough" + + +def test_interp_cart_one_dim_nearest(): + kgrid = kWaveGrid([1000,], [1,]) + binary_sensor_mask = np.zeros((1000,), dtype=bool) + binary_sensor_mask[501,] = True + cart_sensor_mask = np.array([[0.0, 0.0, 0.0]], dtype=np.float32).T + cart_sensor_data = np.array([[1.0, 2.0, 3.0]], dtype=np.float32) + with pytest.raises(ValueError, match=("Data must be two- or three-dimensional.")): + interp_cart_data(kgrid, cart_sensor_data, cart_sensor_mask, binary_sensor_mask) + + def test_nepers2db(): expected_scalar = 8.186258123051049e05 expected_matrix = expected_scalar * np.ones((10, 10)) @@ -415,3 +507,9 @@ def test_trim_zeros(): mat_trimmed, ind = trim_zeros(mat) # TODO: generalize to N-D case + + +if __name__ == "__main__": + test_interp_cart_data_1_point_nearest() + test_interp_cart_data_2_points_nearest() + test_interp_cart_data_2_points_linear()