examples/us_defining_transducer/us_defining_transducer.ipynb

{
 "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.kmedium import kWaveMedium\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",
    "\n",
    "from kwave.options.simulation_execution_options import SimulationExecutionOptions\n",
    "from kwave.utils.dotdictionary import dotdict\n",
    "from kwave.utils.filters import spect\n",
    "from kwave.utils.plot import voxel_plot\n",
    "from kwave.utils.signals import tone_burst"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# simulation settings\n",
    "DATA_CAST = 'single'\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 = 128 - 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])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "# define the medium\n",
    "medium = kWaveMedium(sound_speed=1540,\n",
    "                     density=1000,\n",
    "                     alpha_coeff=0.75,\n",
    "                     alpha_power=1.5,\n",
    "                     BonA=6)\n",
    "\n",
    "# create the time array\n",
    "t_end = 40e-6\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 = 72\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)"
   ]
  },
  {
   "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.zeros((transducer.number_elements, 1))\n",
    "not_transducer.active_elements[21:52] = 1\n",
    "not_transducer.input_signal = input_signal\n",
    "\n",
    "not_transducer = NotATransducer(transducer, kgrid, **not_transducer)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "voxel_plot(np.single(not_transducer.all_elements_mask))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define sensor mask\n",
    "sensor_mask = np.zeros((Nx, Ny, Nz))\n",
    "sensor_mask[Nx//4, Ny//2, Nz//2] = 1\n",
    "sensor_mask[Nx//2, Ny//2, Nz//2] = 1\n",
    "sensor_mask[3*Nx//4, Ny//2, Nz//2] = 1\n",
    "sensor = kSensor(sensor_mask)\n",
    "sensor.record=['p']\n",
    "\n",
    "# SIMULATION\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",
    "execution_options = SimulationExecutionOptions(is_gpu_simulation=True)\n",
    "\n",
    "sensor_data = kspaceFirstOrder3D(medium=medium, kgrid=kgrid, source=not_transducer, sensor=sensor,\n",
    "                                  simulation_options=simulation_options, execution_options=execution_options)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "padded_input_signal = np.concatenate((input_signal, np.zeros((1, 2 * np.shape(input_signal)[1]))), axis=1)\n",
    "f_input, as_input, _ = spect(padded_input_signal, 1/kgrid.dt)\n",
    "_, as_1, _ = spect(sensor_data['p'][:, 0], 1/kgrid.dt)\n",
    "_, as_2, _ = spect(sensor_data['p'][:, 1], 1/kgrid.dt)\n",
    "f, as_3, _ = spect(sensor_data['p'][:, 2], 1/kgrid.dt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(2,1)\n",
    "axes[0].plot(np.arange(0,input_signal.shape[-1]) * kgrid.dt * 1e6, input_signal.T, 'k-')\n",
    "axes[0].set_xlabel('Time [\\mus]')\n",
    "axes[0].set_ylabel('Particle Velocity [m/s]')\n",
    "\n",
    "axes[1].plot(f_input/1e6, np.squeeze(as_input/np.max(as_input)), 'k-')\n",
    "axes[1].plot([tone_burst_freq/1e6, tone_burst_freq/1e6], [0, 1], 'k--')\n",
    "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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Creating a dictionary with the step labels as keys\n",
    "sensor_positions = {\n",
    "    'Sensor Position 1': sensor_data['p'][:,0],\n",
    "    'Sensor Position 2': sensor_data['p'][:,1],\n",
    "    'Sensor Position 3': sensor_data['p'][:,2]\n",
    "}\n",
    "\n",
    "fig, axes = plt.subplots(1,2)\n",
    "fig.set_figwidth(16)\n",
    "fig.tight_layout = True\n",
    "\n",
    "offset = -30e5\n",
    "# Plotting each step using the dictionary\n",
    "for i, (label, data) in enumerate(sensor_positions.items()):\n",
    "    axes[0].plot(kgrid.t_array.squeeze()[:len(data.squeeze())] * 1e6, data.squeeze() + offset * i, label=label)\n",
    "\n",
    "# Set y-ticks and y-labels\n",
    "axes[0].set_yticks([offset * i for i in range(len(sensor_positions.keys()))], list(sensor_positions.keys()))\n",
    "axes[0].set_xlabel('Time [\\u03BCs]')\n",
    "\n",
    "\n",
    "axes[1].plot(f * 1e-6, as_1 / np.max(as_1.flatten()), 'k-', label = 'Sensor Position 1')\n",
    "axes[1].plot(f * 1e-6, as_2 / np.max(as_1.flatten()), 'b-', label = 'Sensor Position 2')\n",
    "axes[1].plot(f * 1e-6, as_3 / np.max(as_1.flatten()), 'r-', label = 'Sensor Position 3')\n",
    "axes[1].legend()\n",
    "axes[1].set_xlabel('Frequency [MHz]')\n",
    "axes[1].set_ylabel('Normalised Amplitude Spectrum [au]')\n",
    "f_max = medium.sound_speed / (2 * dx)\n",
    "axes[1].set_xlim([0, f_max * 1e-6])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "venv",
   "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
}