k-Wave Toolbox Previous   Next

Homogeneous Propagation Medium Example

Overview

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.

 Back to Top

Creating the k-space grid

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).

 Back to Top

Defining the medium properties

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.

 Back to Top

Defining the initial pressure distribution

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;

 Back to Top

Defining the sensor mask

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.

 Back to Top

Running the simulation

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.

 Back to Top


© 2009 Bradley Treeby and Ben Cox.