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 ultrasound and photoacoustic wave fields easy and fast. The simulation codes are based on a time domain solution to linear acoustic equations for heterogeneous media with power law absorption 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 wavefield in a homogeneous or heterogeneous linear medium in either one, two, or three dimensions. The computations are based on a k-space solution to coupled first-order acoustic equations using 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.
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 (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
, andsensor
.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.
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 = 256; % number of pixels in the x (column) direction Nz = 128; % number of pixels in the z (row) direction dx = 50e-6; % pixel width [m] dz = 50e-6; % pixel height [m] kgrid = makeGrid(Nx, dx, Nz, dz); % define the medium properties medium.sound_speed = 1500*ones(Nz, Nx); % [m/s] medium.sound_speed(1:50, :) = 1600; % [m/s] medium.density = 1000*ones(Nz, Nx); % [kg/m^3] medium.density(1:50, :) = 1040; % [kg/m^3] % define the initial pressure using makeDisc disc_x_pos = 120; % [pixels] disc_z_pos = 75; % [pixels] disc_radius = 8; % [pixels] disc_mag = 3; % [au] source.p0 = disc_mag*makeDisc(Nx, Nz, disc_x_pos, disc_z_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 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 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-space 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 Nz by Nx matrices with arbitrary numeric
values. The initial acoustic pressure 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 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
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
.
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
, 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.
![]() |
Troubleshooting | Modelling Wave Propagation | ![]() |
© 2009, 2010, 2011 Bradley Treeby and Ben Cox.