k-Wave Toolbox |
![]() ![]() |
On this page… |
---|
Defining the medium properties |
This example provides a simple demonstration of using k-Wave for the simulation and detection of photoacoustic waves within a two-dimensional homogeneous propagation medium.
The computational architecture of k-Wave is illustrated below (see kspaceFirstOrder2D
for further details). The k-space propagation code is given information about the discretisation of the propagation medium (the k-space grid), the acoustic properties of the medium, the initial photoacoustic pressure, and the location and characteristics of the sensor which detects the photoacoustic wave-field. Additional information such as the time-step and length of the simulation can be either computed automatically or given explicitly. The propagation of the wave-field through the medium is then computed step-by-step using a k-space pseudo-spectral time-domain computation. At each time-step, the pressure values over the sensor mask are stored. These are returned when the time-loop has completed.
The medium discretisation is performed by makeGrid
. Both the size (dx, dz) and number (Nx, Nz) of pixels in each Cartesian direction are used to compute the discretisation, and a k-Wave grid structure is returned.
% create the computational grid Nx = 128; % number of pixels in the x (column) direction Nz = 128; % number of pixels in the z (row) direction dx = 0.1e-3; % pixel width [m] dz = 0.1e-3; % pixel height [m] kgrid = makeGrid(Nx, dx, Nz, dz);
As the k-space computation is based heavily on the FFT, the computation will be fastest for grids where the number of pixels is given by a power of two (see fft2
).
For a homogeneous medium, the sound speed (c
) and density (rho
) are set as scalars in SI units.
% define the properties of the propagation medium c = 1500; % [m/s] rho = 1000; % [kg/m^3]
It is useful to note that for a homogeneous medium, the wave propagation is not dependent on the ambient density (rho
). However, an arbitrary value for rho
must still be given due to code's inherent ability to deal with heterogeneous media.
The initial pressure distribution is set as a matrix of arbitrary numeric values over the medium discretisation given by the k-space grid. Several functions are included in the toolbox for the creation of simple geometric shapes. Here makeDisc
is used to create an initial pressure distribution of two small discs with different diameters.
% create initial pressure distribution using makeDisc disc_magnitude = 5; disc_x_pos = 50; % pixels disc_z_pos = 50; % pixels disc_radius = 8; % pixels disc_1 = disc_magnitude*makeDisc(Nx, Nz, disc_x_pos, disc_z_pos, disc_radius); disc_magnitude = 3; disc_x_pos = 60; % pixels disc_z_pos = 80; % pixels disc_radius = 5; % pixels disc_2 = disc_magnitude*makeDisc(Nx, Nz, disc_x_pos, disc_z_pos, disc_radius); p0 = disc_1 + disc_2;
The sensor mask defines the locations where the pressure field is recorded at each time-step. The sensor mask can be given either as a binary grid of sensor positions (see the Using A Binary Sensor Mask Example), or as a set of arbitrary Cartesian coordinates (given in 2D as a 2 x N matrix) within the dimensions of the centered origin k-space grid. Here makeCartCircle
is used to set the sensor mask to the Cartesian coordinates of a set of evenly spaced points on a circle.
% define a centered circular sensor sensor_radius = 4e-3; % [m] num_sensor_points = 50; sensor_mask = makeCartCircle(sensor_radius, num_sensor_points);
A plot of the initial pressure distribution and the sensor mask (using nearest neighbour interpolation) is given below. By default, the pressure at the Cartesian points in 2D is computed using linear interpolation.
The computation is invoked by calling kspaceFirstOrder2D
with the inputs defined 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 to reds to black, and negative pressures as light to dark blue-greys (see getColorMap
). In this example the time-array t_array
is set to 'auto'
, which sets kspaceFirstOrder2D
to select a suitable array based on the size and properties of the k-space grid, and sensible stability criterion (see makeTime
).
% run the simulation sensor_data = kspaceFirstOrder2D(p0, kgrid, c, rho, 'auto', sensor_mask);
As the function runs, status updates and computational parameters are printed to the command line. This output can be saved by using the inbuilt Matlab function diary
.
Running k-space simulation... dt: 20ns, t_end: 12.06us, time steps: 604 input grid size: 128 by 128 pixels (12.8 by 12.8mm) smoothing p0 distribution... calculating Delaunay triangulation... precomputation completed in 0.99242s starting time loop... computation completed in 5.5464s
When the time loop has completed, the function returns the recorded time series at each of sensor locations defined by sensor_mask
. A visualisation is given below.
% plot the simulated sensor data figure; imagesc(sensor_data, [-1, 1]); colormap(getColorMap); ylabel('Sensor Position'); xlabel('Time Step'); colorbar;
Running the same example with num_sensor_points = 500;
records the wave-field with a much higher spatial resolution.
![]() |
Simulating Photoacoustic Wave Propagation | Using A Binary Sensor Mask | ![]() |
© 2009 Bradley Treeby and Ben Cox.