k-Wave Toolbox |

## Simulations In Three Dimensions Example

On this page…
Overview Creating the k-space grid and defining the medium properties Defining the initial pressure distribution Defining the sensor mask Running the simulation A note on Cartesian sensor masks

### Overview

This example provides a demonstration of using k-Wave for the simulation and detection of photoacoustic waves within a three-dimensional heterogeneous propagation medium. It builds on the Homogeneous Propagation Medium, Heterogeneous Propagation Medium, and Controlling The Absorbing Boundary Layer examples.

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

Back to Top

### Creating the k-space grid and defining the medium properties

Simulations in three-dimensions are performed in an analogous fashion to those in one- and two-dimensions. The medium discretisation is again performed by `makeGrid` with two additional inputs for the y-dimension (note the ordering of the inputs). Similarly, matrices for the properties of the acoustic propagation medium also have an extra dimension. The matrix indexing convention (z, x, y) is an extension to that used in two-dimenions (z, x) where matrix row and column data correspond to the z- and x- axes, respectively (this convention yields predictable results when using `imagesc` to plot image matrices).

    % create the computational grid
    Nx = 64;            % number of pixels in the x direction
    Ny = 64;            % number of pixels in the y direction
    Nz = 64;            % number of pixels in the z direction
    dx = 0.1e-3;        % pixel width [m]
    dy = 0.1e-3;        % pixel width [m]
    dz = 0.1e-3;        % pixel height [m]
    kgrid = makeGrid(Nx, dx, Ny, dy, Nz, dz);

    % define the properties of the propagation medium
    c0 = 1500;                      % [m/s]
    c = c0*ones(Nz, Nx, Ny);        % [m/s]
    c(1:Nz/2, :, :) = c0*1.2;

    rho0 = 1000;                    % [kg/m^3]
    rho = rho0*ones(Nz, Nx, Ny);    % [kg/m^3]
    rho(:, Nx/4:end, :) = rho0*1.2;

Back to Top

### Defining the initial pressure distribution

As in one- and two-dimensions, the initial pressure distribution is simply a matrix of arbitrary numeric values over the medium discretisation given by the k-space grid. Here `makeBall` is used to create an initial pressure distribution of two small filled balls with different diameters.

    % create initial pressure distribution using makeBall
    ball_magnitude = 10;
    ball_x_pos = 36;    % pixels
    ball_y_pos = 36;    % pixels
    ball_z_pos = 36;    % pixels
    ball_radius = 5;    % pixels
    ball_1 = ball_magnitude*makeBall(Nx, Ny, Nz, ball_x_pos, ball_y_pos, ball_z_pos, ball_radius);

    ball_magnitude = 10;
    ball_x_pos = 20;    % pixels
    ball_y_pos = 20;    % pixels
    ball_z_pos = 20;    % pixels
    ball_radius = 3;    % pixels
    ball_2 = ball_magnitude*makeBall(Nx, Ny, Nz, ball_x_pos, ball_y_pos, ball_z_pos, ball_radius);

    p0 = ball_1 + ball_2;

Back to Top

### Defining the sensor mask

Again both Cartesian and binary sensor masks can be used to define the locations where the pressure-field is recorded at each time-step. Several functions for three-dimensional geometry creation are included in the k-Wave toolbox, including `makeSphere` which returns a binary map of single pixel spherical shell created using the midpoint circle algorithm, and `makeCartSphere` which returns the Cartesian location of an evenly distributed array of an arbitrary number of points on the surface of a sphere using the Golden Section Spiral method. Here a linear sensor array is explicitly created.

    % define a series of Cartesian points to collect the data
    z = (-22:2:22)*dz;          % [m]
    x = (-22:2:22)*dx;          % [m]
    y = 22*dy*ones(size(z));    % [m]
    sensor_mask = [x; y; z];

A visualisation of the initial pressure distribution and the sensor mask is given below.

Back to Top

### Running the simulation

The computation is invoked by running `kspaceFirstOrder3D` with the inputs defined above. By default, a visualisation of the propagating wave-field (through the horizonal, median, and frontal planes) and a status bar are displayed. As the size of the computation grid used in this example is particularly small (64 pixels in each direction), the default PML size is reduced to 10 pixels to prevent the sensor mask from being within the absorbing boundary layer (see the Controlling The Absorbing Boundary Layer Example; the PML will not work as effectively with this reduced size).

    % input arguments
    input_args = {'PlotLayout', true, 'PMLSize', 10};

    % run the simulation
    sensor_data = kspaceFirstOrder3D(p0, kgrid, c, rho, 'auto', sensor_mask, input_args{:});

By setting the optional input parameter `'PlotLayout'` to `true`, a plot of the initial pressure, and sound speed and density (if heterogenous) are produced.

As the function runs, status updates and computational parameters are printed to the command line.

    Running k-space simulation...
      dt: 16.6667ns, t_end: 6.15us, time steps: 370
      input grid size: 64 by 64 by 64 pixels (6.4 by 6.4 by 6.4mm)
      smoothing p0 distribution...
      smoothing c distribution...
      smoothing rho distribution...
      precomputation completed in 1.9101s
      starting time loop...
      reordering Cartesian measurement data...
      computation completed in 1min 4.1399s

When the time loop has completed, the function returns the recorded time series at each of sensor locations defined by `sensor_mask`. The ordering is again dependent on whether a Cartesian or binary sensor mask is used. A visualisation is given below.

Back to Top

### A note on Cartesian sensor masks

As with the two-dimensional simulations, a Cartesian sensor mask may be used to define the location of the sensor points in Cartesian space (in place of the binary sensor mask used in this example). Currently, in three-dimensions, only nearest neighbour interpolation is supported. In both two- and three- dimensional simulations, if the optional input `'CartInterp'` is set to `'nearest'` (the default in three-dimensions), the Cartesian points are internally mapped onto a binary grid. If more than one Cartesian sensor point maps to the same location, the duplicated points are discarded. Consequently, if the Cartesian sensor mask has a dense array of points, the number of returned time-series may not correspond to the original number of Cartesian sensor points. Details of discared duplicate points are printed to the command line. This situation can be avoided by increasing the resolution of the computational grid, or by reducing the packing of the Cartesian sensor points.

Back to Top

© 2009 Bradley Treeby and Ben Cox.