Expand interp cart grid testing

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).

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
+}