k-Wave Toolbox |

## Product Overview

On this page…
Introduction to the k-Wave Toolbox Overview of functions Time-domain simulation of acoustic wave-fields Simulation example

### Introduction to the k-Wave Toolbox

k-Wave is a freely available third-party toolbox for MATLAB (released under the GNU LGPLv3 license) developed by Bradley Treeby and Ben Cox. The toolbox is designed to make the time-domain simulation of acoustic wave propagation fast and easy. Some of the main features include:

- an advanced time-domain model of acoustic wave propagation that can account for nonlinearity, acoustic heterogeneities, and power law absorption,
- state of the art algorithms for reconstructing photoacoustic images,
- optional input parameters to adjust visualisation and performance, including options to generate movies and to run the simulations on a graphics processing unit (GPU), and
- many simple-to-follow tutorial examples, for example, the simulation of ultrasound beam patterns and B-mode images, and the reconstruction of photoacoustic images with attenuation compensation.

Back to Top

### Overview Of Functions

The MATLAB functions included within the toolbox can be divided into four basic categories:

1. the simulation of acoustic wave-fields,
2. the reconstruction of photoacoustic images,
3. the creation of geometric shapes, and
4. utility and system functions.

The simulation functions compute the time evolution of an acoustic wavefield in a homogeneous or heterogeneous medium in either one, two, or three dimensions. The computations are based on a k-space psuedospectral solution to coupled first-order acoustic equations defined on a regularly spaced Cartesian grid. The photoacoustic image reconstruction functions compute an estimate of initial photoacoustic pressure given a set of time-varying boundary measurements. The two algorithms included in the toolbox are (a) time reversal, which can accept data recorded from an arbitrarily-shaped measurement surface, and (b) FFT-based reconstruction, which can estimate the initial photoacoustic pressure from data recorded over a linear (2D) or planar (3D) measurement surface. The geometry creation functions allow many useful objects to be defined easily, including circles, arcs, discs, spheres, shells, and balls. They can return the shapes either fitted to a known spatial grid (i.e. as a binary matrix of 1's and 0's, with the 1's corresponding to the location of shape) or in exact geometric coordinates. The utility functions are used to perform additional tasks such as grid creation, matrix smoothing, matrix interpolation, file loading, etc.

Back to Top

### Time-domain simulation of acoustic wave-fields

The computational architecture of the k-Wave simulation functions is illustrated below (see `kspaceFirstOrder1D`, `kspaceFirstOrder2D`, and `kspaceFirstOrder3D` for further details).

The functions are given information about the discretisation of the propagation medium (the k-Wave grid), its acoustic properties (sound speed, absorption, etc), the initial (or time varying) pressure or particle velocity distribution, and the location and characteristics of the measurement surface which detects the acoustic wave-field. These properties are assigned as fields within four input structures:

- `kgrid`
- `medium`,
- `source`, and
- `sensor`.

The propagation of the wave-field through the medium is then computed step-by-step, with the pressure (and velocity) values over the sensor stored after each iteration. These values are returned when the time loop has completed.

Back to Top

### Simulation example

A simple example of a 2D initial value problem in a heterogeneous medium is given below (1D and 3D simulations are performed in an analogous fashion).

    % create the computational grid
    Nx = 128;      % number of grid points in the x (row) direction
    Ny = 256;      % number of grid points in the y (column) direction
    dx = 50e-6;    % grid point spacing in the x direction [m]
    dy = 50e-6;    % grid point spacing in the y direction [m]
    kgrid = makeGrid(Nx, dx, Ny, dy);

    % define the medium properties
    medium.sound_speed = 1500*ones(Nx, Ny); % [m/s]
    medium.sound_speed(1:50, :) = 1600;     % [m/s]
    medium.density = 1000*ones(Nx, Ny);     % [kg/m^3]
    medium.density(1:50, :) = 1040;         % [kg/m^3]

    % define the initial pressure using makeDisc
    disc_x_pos = 75;                        % [grid points]
    disc_y_pos = 120;                       % [grid points]
    disc_radius = 8;                        % [grid points]
    disc_mag = 3;                           % [au]
    source.p0 = disc_mag*makeDisc(Nx, Ny, disc_x_pos, disc_y_pos, disc_radius);

    % define a Cartesian sensor mask of a centered circle with 50 sensor elements
    sensor_radius = 2.5e-3;                 % [m]
    num_sensor_points = 50;
    sensor.mask = makeCartCircle(sensor_radius, num_sensor_points);

    % run the simulation
    sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor);

The medium discretisation is created using the utility function `makeGrid`. The size (dx, dy) and number (Nx, Ny) of grid points in each Cartesian direction are used to calculate the Cartesian and k-space discretisations and an object of the kWaveGrid class (called `kgrid` in this example) is returned. The discretisations are calculated to satisfy the requirements of the FFT-based spatial derivatives. This object is used extensively by both the simulation and utility functions within k-Wave. The time steps used in the simulation are defined by `kgrid.t_array` which is set to `'auto'` by default. In this case, the time array is automatically calculated within the simulation functions using the utility function `makeTime` based on the size and properties of the k-Wave grid and sensible stability criterion (a Courant-Friedrichs-Lewy stability value of 0.3).

For a homogeneous medium, `medium.sound_speed` is given simply as a scalar value (the density is not required unless velocity inputs or outputs are used). For a heterogeneous medium, `medium.sound_speed` and `medium.density` are given as Nx by Ny matrices with arbitrary numeric values. The initial acoustic pressure `source.p0` is similarly defined as an Nx by Ny matrix. The measurement surface `sensor.mask` is given either as a binary grid (i.e., an Nx by Ny matrix of 1's and 0's) representing the grid points within the computational grid at which the data will be recorded, or as a 2 by N matrix of Cartesian coordinates where the pressure values are calculated at each time step using interpolation. Note, in 3D the input matrices are instead Nx by Ny by Nz in size, and the Cartesian sensor points are given by a 3 by N matrix. For the example given here, the geometric function `makeDisc` is used to define an initial pressure distribution of a small filled disc, while `makeCartCircle` is used to define a Cartesian sensor mask with a set of evenly spaced points on a circle.

The simulation is invoked by calling `kspaceFirstOrder2D` with the inputs described above. By default, a visualisation of the propagating wave-field and a status bar are displayed, with frame updates every ten time steps. The default k-Wave color map displays positive pressures as yellows through reds to black, zero pressures as white, and negative pressures as light to dark blue-greys. As the function runs, status updates and computational parameters are printed to the command line. When the time loop has completed, the function returns the time series recorded at the detector locations defined by the sensor mask. If the sensor mask is given as a set of Cartesian coordinates, the computed `sensor_data` is returned in the same order. If the sensor mask is given as a binary grid, `sensor_data` is returned using MATLAB's standard column-wise linear index ordering. In both cases, the recorded data is indexed as `sensor_data(sensor_position, time)`. For a binary sensor mask, the pressure values at a particular time can be restored to the sensor positions within the computational grid using the utility function `unmaskSensorData`. The final pressure within the computational domain can also be obtained by calling the simulation functions with two output arguments `[sensor_data, field_data] = kspaceFirstOrder2D(...)`.

Additional properties of the medium, source, and sensor can be assigned using the remaining structure fields. For example, arbitrary power law absorption can be assigned using `medium.alpha_power` and `medium.alpha_coeff`, nonlinear governing equations can be used by setting `medium.BonA`, a time varying pressure source can be used by defining `source.p` and `source.p_mask`, and a time varying velocity source can be used by defining `source.ux`, `source.uy`, `source.uz` and `source.u_mask`.

The default behaviour of the simulation functions within the toolbox can also be controlled through the use of optional input parameters. These are given as `'param'`, `value` pairs following the primary function inputs. For example, the visualisation can be automatically saved as a movie by setting `'RecordMovie'` to `true`; a plot of the initial pressure distribution can be automatically generated by setting `'PlotLayout'` to `true`; the properties of the absorbing boundary condition can be adjusted; the interpolation method used to calculate the pressure values on a Cartesian sensor mask can be set; the smoothing of input data can be controlled; and the particle velocity recorded by the sensor mask can be returned.

For additional help using the k-Wave toolbox, see the list of examples and the individual function references. The examples are provided with source code that can both be viewed and executed from within the help browser.

Back to Top

© 2009, 2010, 2011 Bradley Treeby and Ben Cox.