k-Wave Toolbox |

## Homogeneous Propagation Medium Example

On this page…
Overview Creating the k-space grid Defining the medium properties Defining the initial pressure distribution Defining the sensor mask Running the simulation

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

- open the file in the MATLAB Editor
- run the file in MATLAB

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.