k-Wave Toolbox |
![]() ![]() |
On this page… |
---|
Introduction to the k-Wave Toolbox |
k-Wave is a freely available third-party MATLAB toolbox developed by Bradley Treeby and Ben Cox (released under the GNU LGPLv3 license). The toolbox is designed to make the time-domain simulation of photoacoustic and ultrasonic wave fields easy and fast. The simulation codes are based on a k-space pseudo-spectral time domain solution to coupled first-order wave equations for heterogeneous media in one, two, and three dimensions.
The toolbox includes:
The functions included within the toolbox can be divided into four categories:
The simulation functions compute the time evolution of an acoustic wave-field within homogeneous or heterogeneous media in either one, two, and three dimensions. The computations are based on a k-space solution to coupled first-order wave equations. These functions can also be used for time reversal image reconstruction. The additional image reconstruction functions allow the initial photoacoustic pressure to be estimated from data recorded over a linear (2D) or planar (3D) measurement surface. The geometry creation functions allow both Cartesian and grid based geometries to be defined, including circles, arcs, discs, spheres, shells, and balls. The Cartesian based functions return the geometric coordinates of the particular shape, while the grid based functions return a binary matrix (i.e., matrix of 1's and 0's) where the 1's correspond to the location of shape. The utility functions are used to perform additional tasks such as grid creation, matrix smoothing, matrix interpolation, file loading etc.
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-space grid), its acoustic properties, the initial (or time varying) pressure 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
,sensor
.The propagation of the wave-field through the medium is then computed step-by-step, with the pressure values over the sensor stored after each iteration. These values are returned when the time loop has completed.
A simple example of a 2D simulation in a heterogeneous layered medium is given below (1D and 3D simulations are performed in an analogous fashion).
% create the computational grid Nx = 256; Nz = 128; dx = 50e-6; dz = 50e-6; kgrid = makeGrid(Nx, dx, Nz, dz); % define the medium properties medium.sound_speed = 1500*ones(Nz, Nx); medium.sound_speed(1:50, :) = 1600; medium.density = 1000*ones(Nz, Nx); medium.density(1:50, :) = 1040; % define the initial pressure disc_x_pos = 120; disc_z_pos = 75; disc_radius = 8; disc_mag = 3; source.p0 = disc_mag*makeDisc(Nx, Nz, disc_x_pos, disc_z_pos, disc_radius); % define a centered circular sensor sensor_radius = 2.5e-3; 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 performed using the utility function
makeGrid
. The size (dx, dz) and
number (Nx, Nz) of pixels in each Cartesian direction are used to calculate
the Cartesian and k-space discretisations and a k-Wave grid structure kgrid
encapsulating this information is returned. The discretisations are calculated
to satisfy the requirements of the FFT-based spatial derivatives. This
structure 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 makeGrid
.
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-space grid and sensible stability
criterion (a Courant-Friedrichs-Lewy stability value of 0.3).
For a homogeneous medium, medium.sound_speed
and medium.density
are given simply as scalar values. For a heterogeneous medium, these are given as
Nz by Nx matrices with arbitrary numeric values. The initial photoacoustic
pressure distribution source.p0
is similarly defined as an Nz by Nx matrix.
The measurement surface sensor.mask
is given either as a binary grid
(i.e., an Nz by Nx matrix of 1's and 0's) representing the pixels within the
computational grid that will collect the data, 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 Nz by Nx by Ny 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
.
Additional properties of the medium, source, and sensor can be similarly
assigned using the remaining structure fields. For example, arbitrary power
law absorption can be assigned using medium.alpha_power
and
medium.alpha_coeff
, sensor directivity can be specified using
sensor.directivity_angle
, sensor.directivity_pattern
and sensor.directivity_size
, and a time varying source can be used
by defining source.p
and source.p_mask
.
The default behaviour of the simulation functions within the toolbox can 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
. Similarly, a plot of the initial
pressure distribution can be automatically generated by setting 'PlotLayout'
to true
; the properties of the perfectly matched layer can be adjusted;
the interpolation method used to calculate the pressure values on a Cartesian sensor mask
can be set; and the smoothing of input matrices can be controlled.
For additional help using the k-Wave toolbox, see the list of examples and the individual function references. These examples are provided with source code that can both be viewed and executed from the help files.
![]() |
Getting Started | Modelling Wave Propagation | ![]() |
© 2009, 2010 Bradley Treeby and Ben Cox.