Delay And Sum

Read data

In this example, we are going to go through a Delay And Sum (DAS) routine using ultraspy on raw RFs simulated data from the PICMUS challenge (you can download them directly from their website, the one used here is the resolution_distorsion_simu_dataset_rf.hdf5 data file. Those Radio-Frequencies were simulated using FieldII on a few scatterers in an empty homogeneous medium. We provide a Reader object to read saved data.

1from ultraspy.io.reader import Reader
2
3reader = Reader('/path/to/picmus_simu_rfs.hdf5', 'picmus')

If successfully loaded, an overview of the data should have been printed:

====== Loaded: ======
Data looking like:
    data_shape: (1, 75, 128, 1527)
    data_type: float32
    is_iq: False
Acquired with the settings:
    sampling_freq: 20832000.0
    t0: 0.0
    prf: 100.0
    signal_duration: 5e-07
    delays: (75, 128) (numpy.ndarray - float64) -> [...]
    sound_speed: 1540.0
    sequence_elements
      emitted: (75, 128) (numpy.ndarray - int32) -> [...]
      received: (75, 128) (numpy.ndarray - int32) -> [...]
Using the probe:
====== Probe: ======
name: L11-4v, with a central freq of 5208000.0Hz
Geometry type: linear
Probe elements:
    lateral (x): [-0.019049999999999997 ; 0.019049999999999997]
    elevational (y): [0.0 ; 0.0]
    axial (z): [0.0 ; 0.0]

Here, we’ve used the ‘picmus’ option in the reader, meaning that it would work only with data from the PICMUS challenge, but the extracted information are in a specific format. Four components have been loaded:

  • data: The data that has been loaded, a numpy array of shape (nb_frames, nb_transmissions, nb_probe_elements, nb_time_samples). Meaning that if we’ve recorded 25 acquisitions of a 11 plane-waves transmission using a 128 linear probe, the format will be (25, 11, 128, nb_time_samples).

  • data_info: A dictionary with three information: the shape of the data, its data type and a boolean telling if the loaded data are I/Qs or RFs.

  • acquisition_info: A dictionary with the main information regarding the acquisition, including the sampling frequency, delays, PRF, etc. Those are detailed more in depth in the other tutorial about how to use your own data.

  • probe: This is a probe object with the information of the used probe (geometry of the elements, central frequency, etc). Check out the probe class description in the dedicated section to pick yours.

If you want to use your own data, you will need to respect this format, or to create your own reader option in the dedicated class. A detailed tutorial on how to do this is available.

Initialize the beamformer

Now we’re going to introduce the DAS beamformer and set the parameters / options for our specific application:

1from ultraspy.beamformers.das import DelayAndSum
2
3# DAS Beamformer
4beamformer = DelayAndSum()
5
6# Set the information about the acquisition (probe, angles, ...) and the
7# options to use (here the transmit method, which needs to be centered for
8# PICMUS data).
9beamformer.automatic_setup(reader.acquisition_info, reader.probe)

Note that by default, ultraspy will try to initialize everything to work on GPU unless you define it otherwise. You can manually initialize the CPU mode by using beamformer = DelayAndSum(on_gpu=False) or, if the beamformer is already initialized, by updating its mode with beamformer.set_on_gpu(False). However, if there is no GPU available on your system (or if cupy hasn’t been installed), it’ll automatically switch back to the CPU mode. If you want to now the current mode of the beamformer, and also its setup, you can print it (print(beamformer)):

====== Beamformer: ======
name: das
(on RFs, with GPU)
    emitted_probe: (3, 75, 128) (cupy.ndarray - float32) -> [...]
    received_probe: (3, 75, 128) (cupy.ndarray - float32) -> [...]
    emitted_thetas: (75, 128) (cupy.ndarray - float32) -> [...]
    received_thetas: (75, 128) (cupy.ndarray - float32) -> [...]
    delays: (75, 128) (cupy.ndarray - float32) -> [...]
    sound_speed: 1540.0
    f_number: 1.0
    t0: 0.0
    signal_duration: 5e-07
    fixed_t0: -2.5e-07
    sampling_freq: 20832000.0
    central_freq: 5208000.0
    bandwidth: 67.0
    prf: 100.0

The first line specifies if the code is running on CPU or GPU, and if it is expecting RFs data or I/Qs (more on this in the doppler tutorial). The next lines are listing all the setups for the beamforming, that you can update any time by using:

1# Additional parameters
2beamformer.update_setup('f_number', 1.75)

A list of the detailed setup and options can be found in the beamformer section.

Define a Scan (region of interest)

We also have to define the area of interest in our medium. Here, the probe is ~4cm wide, centered at 0, and scatterers are spread between 1 and 4.5cm depth, so we are choosing an area covering it all:

 1import numpy as np
 2from ultraspy.scan import GridScan
 3
 4# Zone of interest (in m)
 5x = np.linspace(-20, 20, 500) * 1e-3
 6z = np.linspace(5, 50, 1000) * 1e-3
 7
 8# Or, if you want to define the grid using the lateral / axial steps you
 9# want based on f0 to respect the Shannon-Nyquist theorem
10#wavelength = beamformer.sound_speed / reader.probe.central_freq
11#x = np.arange(-20e-3, 20e-3, wavelength / 4)
12#z = np.arange(5e-3, 50e-3, wavelength / 4)
13
14# Then define a regular grid for the beamforming
15scan = GridScan(x, z)

Beamforming

Time to perform the actual beamforming! Let’s do it on the first frame of PICMUS (there is only one anyway):

 1import ultraspy as us
 2
 3# Get the first frame
 4first_frame = reader.data[0]
 5
 6# On GPU mode, the data needs to be sent to the memory of the GPU.
 7# If is also possible to send directly the numpy array, then the
 8# beamformer will handle it. However, if you want to apply many
 9# beamformers / setups or whatever on the same data, this will
10# prevent it to be stored multiple times on GPU
11import cupy as cp
12d_data = cp.asarray(first_frame, np.float32)
13
14# Actual beamforming operation, then we compute its envelope
15d_output = beamformer.beamform(d_data, scan)
16d_envelope = beamformer.compute_envelope(d_output, scan)

The d_ prefix stands for device, it is a good practice to define your variables this way in order to keep track of what is on GPU or in CPU. Also note that the beamformed signals are still on GPU for now, so we now can use them for something else if needed (Doppler imaging for example).

Now that we have the envelope, we just want to convert it to B-Mode, then recover the data from the GPU memory to the CPU, in order to display it using matplotlib.

1# Get the B-Mode to display
2us.to_b_mode(d_envelope)
3
4# Get the B-Mode back on CPU memory
5b_mode = d_envelope.get()

All good, time to display the B-Mode:

1import matplotlib.pyplot as plt
2extent = [x * 1e3 for x in [x[0], x[-1], z[-1], z[0]]]  # In mm
3plt.imshow(b_mode.T, extent=extent, cmap='gray', clim=[-60, 0])
4plt.title('DAS on PICMUS - 75 plane waves')
5plt.xlabel('Axial (mm)')
6plt.ylabel('Depth (mm)')
7plt.show()
../_images/das_bmode.png

And voilà! Feel free to have a look to the next tutorial (about Doppler) if you want to see another application of ultraspy.