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.

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

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

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

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

 Back to Top


© 2009, 2010 Bradley Treeby and Ben Cox.