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 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 and density are set as scalars
in SI units. Power law absorption can also be optionally set by
assigning values to medium.alpha_power
and
medium.alpha_coeff
.
% define the properties of the propagation medium medium.sound_speed = 1500; % [m/s] medium.density = 1000; % [kg/m^3] medium.alpha_power = 1.5; % [dB/(MHz^y cm)] medium.alpha_coeff = 0.75; % [dB/(MHz^y cm)]
It is useful to note that for a homogeneous medium, the wave propagation
is not dependent on the ambient density. However, an arbitrary value for
medium.density
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. This is assigned
to the p0
field of the source
structure.
% 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); source.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
).
% run the simulation sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor);
As the function runs, status updates and computational parameters are printed to the command line.
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) maximum supported frequency: 7.5MHz smoothing p0 distribution... calculating Delaunay triangulation... precomputation completed in 1.3966s starting time loop... computation completed in 6.9013s
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, 2010 Bradley Treeby and Ben Cox.