beamform_image.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "beamform_image.ipynb",
      "provenance": [],
      "toc_visible": true,
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/luuleitner/dasIT/blob/main/beamform_image.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "\n",
        "<h1>Plane Wave Beamforming Tutorial</h1>\n",
        "\n",
        "This is a hands-on introduction to plane wave ultrasound beamforming. During this exercise we will look at the raw ultrasound data acquired and how to convert it into actual images. \n",
        "\n",
        "<br>\n",
        "\n",
        "<h2>Main Take Aways</h2>\n",
        "\n",
        "- understand ultrasonic raw data,\n",
        "- know the necessary processing elements of a plane wave beamformer, and\n",
        "- be able to perform a quality estimation of an ultrasound image.\n",
        "\n",
        "<br>\n",
        "\n",
        "<h2>Supplemental Information</h2>\n",
        "\n",
        "This practical example is part of the ***Graz University of Technology*** lecture series ***Development of Electronic Systems*** and the ***Fundamentals of Biomedical Engineering Laboratory***. \n",
        "\n",
        "This program is free software and licensed under the ***Apache License v2.0*** (see the github LICENSE file for details). It is also citable: ***Leitner, C. (2022). New directions in recording and processing electro-mechanical signals from the human body. Dissertation. Graz University of Technology.***\n",
        "\n",
        "<br>\n",
        "\n",
        "<h2>Uncomplete List of Other Free Ultrasound Ressources:</h2>\n",
        "\n",
        "*   <a href=\"http://www.k-wave.org/\">k-wave ultrasound simulator</a> - free MATLAB and C++ implementations\n",
        "*   <a href=\"https://field-ii.dk/\">field II ultrasound simulator</a> - free MATLAB implementation\n",
        "*   <a href=\"https://github.com/Sergio5714/pybf\">pybf - Python beamformer</a> - optimized for short processing times\n",
        "*   <a href=\"https://www.biomecardio.com/MUST/\">MATLAB ultrasound toolbox</a> - free MATLAB beamformer<br><br>\n",
        "\n",
        "<br><br>\n",
        "\n",
        "I wish everyone a great dive into the topic, and please do not hesitate to <a href=\"mailto:christoph.leitner@tugraz.at\">contact</a> me in case of **any** questions!\n",
        "\n",
        "yours,<br>\n",
        "Christoph<br><br>\n",
        "\n",
        "<br>\n",
        "<br>\n",
        "\n",
        "\n",
        "---\n",
        "\n",
        "\n"
      ],
      "metadata": {
        "id": "oHQTS3hehKVL"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Introduction\n",
        "\n",
        "\n",
        "<br>\n",
        "\n",
        "**But what is beamforming anyway?** Essentially, this involves generating ultrasonic waves (beam patterns) for a transmit event (TX) and reconstructing images from the reflected ultrasonic waves or received signals (RX) on the piezo. \n",
        "\n",
        "<br>\n",
        "\n",
        "**Simplifications:** For ultrasound image reconstruction, we make some additional assumptions, e.g..\n",
        "- sound propagates throughout the tissue at a constant speed of sound,\n",
        "- sound propagates in straight lines,\n",
        "- sound is attenuated by the same amount throughout the tissue.\n",
        "\n",
        "However, all of these assumptions are only model simplifications. Therefore, we must be aware that image distortion can occur if one of these basic assumptions are not met.\n",
        "\n",
        "<br>\n",
        "\n",
        "<h2>Plane wave beamforming principles</h2>\n",
        "\n",
        "This method of beamforming is probably the easiest way to get started. In **plane wave beamforming**, all transducer elements are activated simultaneously, creating a single, unfocused ultrasonic plane wave that propagates in a straight line into the tissue (see figure below). \n",
        "\n",
        "<br>\n",
        "\n",
        "<figure>\n",
        "<center>\n",
        "<img src='https://raw.githubusercontent.com/luuleitner/dasIT/main/example_data/ColabBFfigures/PlaneWaveBF.gif'/ width=\"1200\">\n",
        "</figure>\n",
        "\n",
        "<br>\n",
        "\n",
        "**How does this simplify things?** With a beamformer we have to solve the inverse problem, i.e. the transfer of the ultrasonic signal amplitudes from the time-channel space into an image space. In the case of beamforming plane waves, the solution simplifies to straightforward geometric formulations (we will have a look at them later on).\n",
        "\n",
        "<br>\n",
        "\n",
        "---"
      ],
      "metadata": {
        "collapsed": false,
        "pycharm": {
          "name": "#%% md\n"
        },
        "id": "IfDb8Upqvz_X"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "# 1.Getting Started\n",
        "\n",
        "\n",
        "```\n",
        "[#] Areas shown like this are executable code. Use the mousover play button to run these cells.\n",
        "```"
      ],
      "metadata": {
        "collapsed": false,
        "pycharm": {
          "name": "#%% md\n"
        },
        "id": "xhS9t3Ylvz_Z"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Installation\n",
        "First we need to install the beamformer package from the the <a href=\"https://github.com/luuleitner/dasIT\">GitHub repository</a>:"
      ],
      "metadata": {
        "id": "LPoxjTFywf8T"
      }
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "outputs": [],
      "source": [
        "# Fetch the newest dasIT package from the github repository\n",
        "\n",
        "!pip install git+https://github.com/luuleitner/dasIT"
      ],
      "metadata": {
        "pycharm": {
          "name": "#%%\n"
        },
        "id": "rqSNgrwqvz_a"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Next we load all other necessary libraries into our Colab notebook:"
      ],
      "metadata": {
        "collapsed": false,
        "pycharm": {
          "name": "#%% md\n"
        },
        "id": "VzuCr8Bkvz_b"
      }
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "outputs": [],
      "source": [
        "# import necessary packages\n",
        "import os\n",
        "import numpy as np\n",
        "from datetime import datetime\n",
        "from matplotlib import pyplot as plt\n",
        "\n",
        "# to run dasIT, we need to import the necessary commands from the library\n",
        "from dasIT.data.loader import RFDataloader, TDloader, TGCloader\n",
        "from dasIT.features.transducer import transducer\n",
        "from dasIT.features.medium import medium\n",
        "from dasIT.features.tgc import tg_compensation\n",
        "from dasIT.src.delays import planewave_delays\n",
        "from dasIT.src.apodization import apodization\n",
        "from dasIT.src.das_bf import RXbeamformer\n",
        "from dasIT.features.signal import RFfilter, fftsignal, analytic_signal\n",
        "from dasIT.features.image import interp_lateral\n",
        "from dasIT.visualization.signal_callback import amp_freq_1channel, amp_1channel\n",
        "from dasIT.visualization.image_callback import plot_signal_grid, plot_signal_image\n",
        "from dasIT.features.signal import logcompression"
      ],
      "metadata": {
        "pycharm": {
          "name": "#%%\n"
        },
        "id": "aHXTvjxHvz_b"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Download the example dataset"
      ],
      "metadata": {
        "collapsed": false,
        "pycharm": {
          "name": "#%% md\n"
        },
        "id": "NmxVlWQhvz_b"
      }
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "outputs": [],
      "source": [
        "# Download the dataset from github\n",
        "\n",
        "# The data was obtained in the course of the work for:\n",
        "# Leitner et al. 2020, \"Detection of Motor Endplates in Deep and Pennate Skeletal Muscles in-vivo using Ultrafast Ultrasound\",\n",
        "# 2020 IEEE International Ultrasonics Symposium (IUS).\n",
        "#\n",
        "\n",
        "rfdata_path = '/content/rfdata'\n",
        "\n",
        "if os.path.exists(rfdata_path) == False:\n",
        "  os.mkdir(rfdata_path)\n",
        "  os.chdir(rfdata_path)\n",
        "  !wget -i https://raw.githubusercontent.com/luuleitner/dasIT/main/example_data/CIRSphantom_GE9LD_VVantage/COLABdownload_url.txt\n",
        "\n",
        "os.chdir(rfdata_path)"
      ],
      "metadata": {
        "pycharm": {
          "name": "#%%\n"
        },
        "id": "Pr-F6shMvz_c"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Compilation of the data set\n",
        "Five frames were captured on a **Verasonics Vantage 256** ultrasound research system, using a **GE-9LD** transducer, and a **CIRS General Purpose** ultrasound phantom.\n",
        "\n",
        "The figure below shows the reconstructed ultrasonic image (left) captured on an ultrasound phantom (center) with a GE-9LD transducer (right).\n",
        "\n",
        "<figure>\n",
        "<center>\n",
        "<img src='https://raw.githubusercontent.com/luuleitner/dasIT/main/example_data/ColabBFfigures/ExperimentalSetup.jpg'/ width=\"1200\">\n",
        "</figure>"
      ],
      "metadata": {
        "id": "Rqcy-TxDwBUT"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "# 2.Define Hardware and Imaging Medium\n",
        "\n",
        "To get started, we need to provide our software several static parameters. In particular, we have to define the used ultrasound hardware (e.g. transducer design,..) the emitted ultrasound signal (e.g. used ultrasound frequency,..) and the definition of the medium (e.g. speed of sound, size, pixels,..) below the ultrasound lens.\n",
        "\n",
        "We assume a constant speed of sound (**v = 1540 m/s**) in the tissue. This model assumption allows us to infer the tissue depth of the reflected echo by the relationship v=s/t and the measurement of the delay time of the reflected sound (Puls-Echo method)."
      ],
      "metadata": {
        "id": "HFufHVD70v8C"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Transducer\n",
        "\n",
        "Based on the transducer specs given above fill in the necessary parameters below:"
      ],
      "metadata": {
        "collapsed": false,
        "pycharm": {
          "name": "#%% md\n"
        },
        "id": "OOt8MmWQvz_c"
      }
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "outputs": [],
      "source": [
        "# dasIT transducer\n",
        "physical_transducer = TDloader('transducer.csv')\n",
        "dasIT_transducer = transducer(center_frequency_hz = ,  # <--- FILL IN CENTER FREQUENCY OF THE TRANSDUCER IN [Hz]\n",
        "                              bandwidth_hz=physical_transducer.transducer['bandwidth'].dropna().to_numpy(dtype='float', copy=False),    # [Hz]\n",
        "                              adc_ratio=4,  # [-]\n",
        "                              transducer_elements_nr = , # <--- FILL IN THE NUMBER OF TRANSDUCER ELEMENTS [#]\n",
        "                              element_pitch_m = , # <--- FILL IN THE ELEMENT PITCH IN [m]\n",
        "                              pinmap=physical_transducer.transducer['pinmap'].dropna().to_numpy(dtype='int', copy=False),   # [-]\n",
        "                              pinmapbase=1, # [-]\n",
        "                              elevation_focus=0.028, # [m]\n",
        "                              focus_number=None,\n",
        "                              totalnr_planewaves=1,     # [-]\n",
        "                              planewave_angle_interval=[0,0],   # [rad]\n",
        "                              axial_cutoff_wavelength=5,  # [#]\n",
        "                              speed_of_sound_ms = )  # <--- FILL IN THE SPEED OF SOUND IN [m/s]"
      ],
      "metadata": {
        "pycharm": {
          "name": "#%%\n"
        },
        "id": "h3RhVX6avz_c"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Print the transducer specifications\n",
        "\n",
        "print(f'Transducer properties:')\n",
        "print()\n",
        "vars(dasIT_transducer)"
      ],
      "metadata": {
        "id": "lCmUwz764_-9"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Medium\n",
        "\n",
        "The size of the medium we are imaging is usually specified in \"wavelength numbers\" via the relationship lambda = c/f. We don't worry too much about these definitions and simply load the specifications of the imaging medium."
      ],
      "metadata": {
        "id": "B3QCH7qw5HQ0"
      }
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "outputs": [],
      "source": [
        "# dasIT medium\n",
        "dasIT_medium = medium(speed_of_sound_ms = 1540, # [m/s]\n",
        "                      center_frequency = dasIT_transducer.center_frequency, # [Hz]\n",
        "                      sampling_frequency = dasIT_transducer.sampling_frequency, # [Hz]\n",
        "                      max_depth_wavelength = 177,   # [#]\n",
        "                      lateral_transducer_element_spacing = dasIT_transducer.lateral_transducer_spacing, # [m]\n",
        "                      axial_extrapolation_coef = 1.05,  # [-]\n",
        "                      attenuation_coefficient= 0.75,   # [dB/(MHz^y cm)]\n",
        "                      attenuation_power=1.5   # [-]\n",
        "                      )"
      ],
      "metadata": {
        "pycharm": {
          "name": "#%%\n"
        },
        "id": "zueTe5JNvz_d"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "# 3.Data Wrangling\n",
        "\n",
        "From an ultrasound research system, we usually get more information than we actually need (due to the memory allocation of the system), and the channel data is usually unsorted. Therefore, we need to \"clean up\" before we begin. So we cut off unassigned samples (zero values) and sort our channel data in ascending order from 1-192. "
      ],
      "metadata": {
        "collapsed": false,
        "pycharm": {
          "name": "#%% md\n"
        },
        "id": "qPLrl1IWvz_d"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "<figure>\n",
        "<center>\n",
        "<img src='https://raw.githubusercontent.com/luuleitner/dasIT/main/example_data/ColabBFfigures/DataPreperation.jpg'/ width=\"1200\">\n",
        "</figure>\n"
      ],
      "metadata": {
        "id": "Ny4oxxQ9AQP2"
      }
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "outputs": [],
      "source": [
        "### Load RF Data\n",
        "RFdata = RFDataloader('CIRS_phantom.h5')\n",
        "\n",
        "### Preprocess (Clip and Sort) RF Data\n",
        "# Samples start: at first recorded echo (number of wavelength distance is provided from vendor)\n",
        "# -> null out the rest to not overshadow the real results\n",
        "# Samples end: at penetration depth -> clip rest of samples without data\n",
        "# If necessary sort the transducer channels according to the pin map to get the channels first-last channel\n",
        "RFdata.signal[:dasIT_transducer.start_depth_rec_samples, :, :] = 0\n",
        "RFdata.signal = RFdata.signal[:dasIT_medium.rx_echo_totalnr_samples, dasIT_transducer.transducer_pinmap, :]\n",
        "\n",
        "print(f'Channels of transducer: {RFdata.signal.shape[1]}')\n",
        "print(f'Samples per channel: {RFdata.signal.shape[0]}')\n",
        "print(f'Number of frames: {RFdata.signal.shape[2]}')"
      ],
      "metadata": {
        "pycharm": {
          "name": "#%%\n"
        },
        "id": "gE0pvhsovz_d"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Exploring raw data - Beamforming Principle"
      ],
      "metadata": {
        "id": "kWHSu5Ho4ZFn"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "If we now plot the sorted and trimmed channel data over time (left image), point-like echoes (phantom reflectors in the right image) show up, hyperbolic in the x,t plots.\n",
        "\n",
        "**With our beamformer, we now need to solve the inverse problem and transfer the raw ultrasound data (time and channel space) into an image (width and height space)**.\n",
        "\n",
        "<figure>\n",
        "<center>\n",
        "<img src='https://raw.githubusercontent.com/luuleitner/dasIT/main/example_data/ColabBFfigures/Beamformer_xt_xz.jpg'/ width=\"1200\">\n",
        "</figure>"
      ],
      "metadata": {
        "id": "LA7W8wOwzGvl"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Exploring data of a single channel\n",
        "Let's have a look at the raw data of **Channel 156**. We plot the signal amplitude over time."
      ],
      "metadata": {
        "id": "dR6NrS64Bj0G"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "cutoff_lens = 90\n",
        "cutoff_depth = len(RFdata.signal) - 700\n",
        "\n",
        "channel = 156\n",
        "dbrange = 55\n",
        "signal = logcompression(RFdata.signal[:,:,0], dbrange)\n",
        "\n",
        "fig = plt.figure(figsize=(10, 4), dpi=100)\n",
        "ax_1 = fig.add_subplot(121)\n",
        "ax_1.imshow(signal[90:-700, :],\n",
        "            aspect=1,\n",
        "            cmap='gray')\n",
        "\n",
        "ax_1.set_xlabel('Transducer Element [#]', fontsize=15, fontweight='bold', labelpad=10)\n",
        "ax_1.set_ylabel('Passing Time (Sample [#])', fontsize=15, fontweight='bold', labelpad=10)\n",
        "ax_1.xaxis.tick_top()\n",
        "ax_1.xaxis.set_label_position('top')\n",
        "ax_1.axvline(x=channel, color='red')\n",
        "ax_1.axhline(y=0, color='blue', lw=5)\n",
        "ax_1.axhline(y=cutoff_depth-cutoff_lens-2, color='blue')\n",
        "\n",
        "\n",
        "ax_2 = fig.add_subplot(122)\n",
        "ax_2.plot(RFdata.signal[:,channel,0], 'r')\n",
        "ax_2.set_xlabel('Samples [#]')\n",
        "ax_2.set_ylabel('Signal [V]')\n",
        "ax_2.set_title(f'RF-data channel {channel}')\n",
        "ax_2.axvline(x=cutoff_lens, color='blue',)\n",
        "ax_2.axvline(x=cutoff_depth, color='blue')\n",
        "\n",
        "plt.xticks(fontsize=12)\n",
        "plt.yticks(fontsize=12)\n",
        "plt.tight_layout()\n",
        "plt.show()\n",
        "\n",
        "print('')\n",
        "print('Channel 156 is indicated in red color.\\nThe visible area in the time-signal map on the left is marked with blue lines in the channel plot on the right side.')"
      ],
      "metadata": {
        "id": "v8bS5yk0WJI0"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# 4.Signal Processing"
      ],
      "metadata": {
        "collapsed": false,
        "pycharm": {
          "name": "#%% md\n"
        },
        "id": "b_ayxHWgvz_e"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Time Gain Compensation\n",
        "\n",
        "We assume that acoustic waves are absorbed to the same extent everywhere in the tissue (which is not formally correct, since different tissue types absorb to different amounts). The pressure amplitudes of acoustic waves lose energy as they propagate through tissue due to attenuation. Therefore, reflections from deeper structures appear weaker. For this reason, time gain correction is usually applied to the recorded RF data by adding a time-dependent gain to compensate for the attenuation losses.\n",
        "\n",
        "<figure>\n",
        "<center>\n",
        "<img src='https://raw.githubusercontent.com/luuleitner/dasIT/main/example_data/ColabBFfigures/TGC.jpg'/ width=\"1200\">\n",
        "</figure>"
      ],
      "metadata": {
        "collapsed": false,
        "pycharm": {
          "name": "#%% md\n"
        },
        "id": "9Fuq41FWvz_e"
      }
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "outputs": [],
      "source": [
        "# Load and apply tgc-waveform\n",
        "tgc_cntrl_points = TGCloader(controlpt_path='tgc_cntrl_pt.csv')\n",
        "TGCsignals = tg_compensation(signals=RFdata.signal,\n",
        "                             medium=dasIT_medium,\n",
        "                             center_frequency=dasIT_transducer.center_frequency,\n",
        "                             cntrl_points=tgc_cntrl_points,\n",
        "                             mode='points')"
      ],
      "metadata": {
        "pycharm": {
          "name": "#%%\n"
        },
        "id": "cPGbYGS-vz_e"
      }
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "outputs": [],
      "source": [
        "# Plot channel 156\n",
        "fig = plt.figure(figsize=(5, 3), dpi=100)\n",
        "ax_1 = fig.add_subplot(211)\n",
        "ax_1.plot(RFdata.signal[:,channel,0])\n",
        "ax_1.set_xlabel('Samples [#]')\n",
        "ax_1.set_ylabel('Signal [V]')\n",
        "ax_1.set_title(f'RF-data channel {channel}')\n",
        "\n",
        "ax_2 = fig.add_subplot(212)\n",
        "ax_2.plot(TGCsignals.signals[:,channel,0],'r')\n",
        "ax_2.set_xlabel('Samples [#]')\n",
        "ax_2.set_ylabel('Signal [V]')\n",
        "ax_2.set_title(f'TGC compensated RF-data channel {channel}')\n",
        "\n",
        "plt.tight_layout()\n",
        "plt.show()"
      ],
      "metadata": {
        "pycharm": {
          "name": "#%%\n"
        },
        "id": "YCmSTfbkvz_e"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "###Filtering\n",
        "\n",
        "To eliminate unwanted noise, we filter our raw data with a FIR bandpass filter and Gaussian window. We use the upper and lower limits of the transducer bandwidth as cutoff frequencies."
      ],
      "metadata": {
        "collapsed": false,
        "pycharm": {
          "name": "#%% md\n"
        },
        "id": "GZA8hn7Bvz_e"
      }
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "outputs": [],
      "source": [
        "### Filter RF Data\n",
        "RFdata_filtered = RFfilter(signals=TGCsignals.signals,\n",
        "                           fcutoff_band=dasIT_transducer.bandwidth,\n",
        "                           fsampling=dasIT_transducer.sampling_frequency,\n",
        "                           type='gaussian',\n",
        "                           order=10)"
      ],
      "metadata": {
        "pycharm": {
          "name": "#%%\n"
        },
        "id": "y_-70Gibvz_e"
      }
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "outputs": [],
      "source": [
        "# Plot channel 156\n",
        "fftFil = fftsignal(RFdata_filtered.signal[:,channel,0], dasIT_transducer.sampling_frequency)\n",
        "\n",
        "fig = plt.figure(figsize=(6, 5), dpi=100)\n",
        "ax_1 = fig.add_subplot(311)\n",
        "ax_1.plot(RFdata.signal[:,channel,0])\n",
        "ax_1.set_xlabel('Samples [#]')\n",
        "ax_1.set_ylabel('Signal [V]')\n",
        "ax_1.set_title(f'RF-data channel {channel}')\n",
        "\n",
        "ax_2 = fig.add_subplot(312)\n",
        "ax_2.plot(RFdata_filtered.signal[:,channel,0],'r')\n",
        "ax_2.set_xlabel('Samples [#]')\n",
        "ax_2.set_ylabel('Signal [V]')\n",
        "ax_2.set_title(f'Filtered RF-data channel {channel}')\n",
        " \n",
        "ax_3 = fig.add_subplot(313)\n",
        "ax_3.plot(fftFil[0], fftFil[1], 'r')\n",
        "ax_3.set_xlabel('Frequency [MHz]')\n",
        "ax_3.set_ylabel('Power [W/Hz]')\n",
        "ax_3.set_title(f'FFT channel {channel}')\n",
        "\n",
        "plt.tight_layout()\n",
        "plt.show()"
      ],
      "metadata": {
        "pycharm": {
          "name": "#%%\n"
        },
        "id": "6v_udIY4vz_e"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Convert to analytical signal\n",
        "\n",
        "In the ultrasound domain, we prefer to use signal envelopes to represent data We display the instantaneous energy distributions from fluctuating raw signals. Therefore, we must first calculate the analytical signal of the raw data using the Hilbert transformation."
      ],
      "metadata": {
        "collapsed": false,
        "pycharm": {
          "name": "#%% md\n"
        },
        "id": "Nl-8BN5Yvz_e"
      }
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "outputs": [],
      "source": [
        "####################################################################\n",
        "#------------------------ Analytical Signal -----------------------#\n",
        "\n",
        "### Hilbert Transform\n",
        "RFdata_analytic = analytic_signal(np.squeeze(RFdata_filtered.signal), interp=False)"
      ],
      "metadata": {
        "pycharm": {
          "name": "#%%\n"
        },
        "id": "hbyY3KObvz_f"
      }
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "outputs": [],
      "source": [
        "# Plot channel 156\n",
        "fftFil = fftsignal(RFdata_filtered.signal[:,channel,0], dasIT_transducer.sampling_frequency)\n",
        "\n",
        "fig = plt.figure(figsize=(5, 3), dpi=100)\n",
        "ax_1 = fig.add_subplot(211)\n",
        "ax_1.plot(RFdata.signal[:,channel,0])\n",
        "ax_1.set_xlabel('Samples [#]')\n",
        "ax_1.set_ylabel('Signal [V]')\n",
        "ax_1.set_title(f'RF-data channel {channel}')\n",
        "\n",
        "ax_2 = fig.add_subplot(212)\n",
        "ax_2.plot(RFdata_analytic[:,channel,0].real,'r')\n",
        "ax_2.plot(abs(RFdata_analytic[:,channel,0]),'g', label='envelope')\n",
        "ax_2.set_xlabel('Samples [#]')\n",
        "ax_2.set_ylabel('Signal [V]')\n",
        "ax_2.set_title(f'Analytic signal channel {channel}')\n",
        "ax_2.legend(loc='lower right')\n",
        "\n",
        "plt.tight_layout()\n",
        "plt.show()"
      ],
      "metadata": {
        "pycharm": {
          "name": "#%%\n"
        },
        "id": "Hm15e40kvz_f"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "# 5.Beamforming\n",
        "\n",
        "\n",
        "To solve the beamforming operation, signal amplitudes must be transferred from the (x,time) domain into the (z,x) plane. The figure below depicts the necessary mathematical operations to calculate the time delays needed to determine\n",
        "the signal amplitudes at each pixel point:\n",
        "\n",
        "<figure>\n",
        "<center>\n",
        "<img src='https://raw.githubusercontent.com/luuleitner/dasIT/main/example_data/ColabBFfigures/PWbeamforming.jpg'/ width=\"800\">\n",
        "</figure>\n"
      ],
      "metadata": {
        "collapsed": false,
        "pycharm": {
          "name": "#%% md\n"
        },
        "id": "bEnHRQ0-vz_f"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Acount for element directivity and focus with an apodization table\n",
        "\n",
        "Due to the directivity of the elements, the transducer elements cannot transmit and receive equally in all directions. Therefore, only a partial aperture of the elements contributes to the signal amplitude (returning waves have a hyperbolic shape). This is taken into account in the apodization table."
      ],
      "metadata": {
        "collapsed": false,
        "pycharm": {
          "name": "#%% md\n"
        },
        "id": "aMeoyhRcvz_f"
      }
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "outputs": [],
      "source": [
        "####################################################################\n",
        "#-------------------------- Apodization Table --------------------------#\n",
        "\n",
        "apodization = apodization(delays=None,\n",
        "                          medium=dasIT_medium.medium,\n",
        "                          transducer=dasIT_transducer,\n",
        "                          apo='rec',\n",
        "                          angles=dasIT_transducer.planewave_angles())"
      ],
      "metadata": {
        "pycharm": {
          "name": "#%%\n"
        },
        "id": "ns-0OewAvz_f"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Setup delay table for signal transformation\n",
        "\n",
        "Here we calculate the signal delay between the emission of the ultrasonic pulse, the impact on a reflector target and the back reflection to the transducer element."
      ],
      "metadata": {
        "id": "yv_mMjyut3K2"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "####################################################################\n",
        "#-------------------------- Delay Tables --------------------------#\n",
        "\n",
        "### DAS delay tabels for tilted planewaves\n",
        "delay_table = planewave_delays(medium=dasIT_medium.medium,\n",
        "                               sos=dasIT_medium.speed_of_sound,\n",
        "                               fsampling=dasIT_transducer.sampling_frequency,\n",
        "                               angles=dasIT_transducer.planewave_angles())"
      ],
      "metadata": {
        "id": "cBjzzA5XuFwu"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "###Beamforming\n",
        "\n",
        "Now we apply the apodization and delay tables to the processed raw ultrasound data. In this step, we solve the inverse problem and map the time-channel data to an image space."
      ],
      "metadata": {
        "collapsed": false,
        "pycharm": {
          "name": "#%% md\n"
        },
        "id": "J--ICU8Fvz_f"
      }
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "outputs": [],
      "source": [
        "####################################################################\n",
        "#-------------------------- Beamforming ---------------------------#\n",
        "start_das_timing = datetime.now()\n",
        "\n",
        "# Mask images areas in axial direction which have been included for reconstruction\n",
        "# but are not part of the actual image.\n",
        "RFsignals = RFdata_analytic[:,:,0]\n",
        "\n",
        "RFsignals = np.expand_dims(RFsignals, 2)\n",
        "RFsignals = np.repeat(RFsignals, RFsignals.shape[1], axis=2)\n",
        "RFsignals = np.expand_dims(RFsignals, 3)\n",
        "\n",
        "BFsignals = RXbeamformer(signals=RFsignals,\n",
        "                         delays=delay_table.sample_delays,\n",
        "                         apodization=apodization.table)"
      ],
      "metadata": {
        "pycharm": {
          "name": "#%%\n"
        },
        "id": "pzNe1_ulvz_f"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "# 6.Image Formation\n",
        "\n",
        "The so-called dynamic range defines the logarithmic ratio between the highest (strongest echo near the lens) and the lowest signal amplitude (weak echo from deep targets). Assuming the beamformed data is mapped linearly onto a grayscale image image, significant tissue features are lost because a few echoes with high amplitude echoes overshadow the rest of the image space. Therefore, we balance the dynamic range by logarithmically adjusting the echo amplitudes.\n",
        "\n",
        "That's it. Now we can plot the results graphically:"
      ],
      "metadata": {
        "collapsed": false,
        "pycharm": {
          "name": "#%% md\n"
        },
        "id": "Pj7jjKJ8vz_f"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "####################################################################\n",
        "#------------------------ Image Formation --------------------------\n",
        "\n",
        "# Envelope\n",
        "BFsignals.envelope = abs(BFsignals.frame)\n",
        "\n",
        "# Interpolate over Lateral space\n",
        "BFsignals.interpolated = interp_lateral(signals=BFsignals.envelope,\n",
        "                                        transducer=dasIT_transducer,\n",
        "                                        medium=dasIT_medium,\n",
        "                                        scale=3)\n",
        "\n",
        "\n",
        "# Plot image\n",
        "plot_signal_grid(signals=BFsignals.interpolated.signals_lateral_interp,\n",
        "                 axis_vectors_xz=BFsignals.interpolated.imagegrid_mm,\n",
        "                 axial_clip=[dasIT_transducer.start_depth_rec_m, None],\n",
        "                 compression=True,\n",
        "                 dbrange=58,\n",
        "                 path=rfdata_path)"
      ],
      "metadata": {
        "id": "kINmM40C5nrW"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# 7.Calculate Image Resolution\n",
        "\n",
        "For the indicated point target in the near field (figure below) calculate the axial and lateral resolutions of your beamformed image. \n",
        "\n",
        "<figure>\n",
        "<center>\n",
        "<img src='https://raw.githubusercontent.com/luuleitner/dasIT/main/example_data/ColabBFfigures/Resolution.jpg'/ width=\"300\">\n",
        "</figure>\n",
        "\n",
        "\n",
        "Instructions: \n",
        "- Use a ruler and measure the full width (start and end points) of your reflection target in the lateral and axial directions. \n",
        "- One full width is equals approximately 6 standard deviations (± 3$\\sigma$ = 99.73%).\n",
        "- Divide the full width by 6 to get the approximate value of one sigma.\n",
        "- Calculate the full width at half maximum (FWHM) using a Gaussian distribution: $FWHM = 2 \\sqrt{2 ln2} * \\sigma$\n",
        "- You can do the calculation in <a href=\"https://www.wolframalpha.com/\">https://www.wolframalpha.com/</a>"
      ],
      "metadata": {
        "id": "hvN5S63KU8lh"
      }
    }
  ]
}