k-Wave Toolbox |

## Simulations In Three Dimensions Example

On this page…
Overview Creating the k-space grid and defining the medium properties Defining the time varying pressure source Defining the sensor mask Running the simulation

### Overview

This example provides a simple demonstration of using k-Wave for the simulation and detection of a time varying pressure source within a three-dimensional heterogeneous propagation medium. It builds on the Point Source In A Homogeneous Propagation Medium and Simulations In Three Dimensions 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
    medium.sound_speed = 1500*ones(Nz, Nx, Ny);	% [m/s]
    medium.sound_speed(1:Nz/2, :, :) = 1800;    % [m/s]
    medium.density = 1000*ones(Nz, Nx, Ny);     % [kg/m^3]
    medium.density(:, Nx/4:end, :) = 1200;      % [kg/m^3]

Back to Top

### Defining the time varying pressure source

As in one- and two-dimensions, a time varying pressure source is defined by assigning a binary matrix (i.e., a matrix of 1's and 0's with the same dimensions as the computational grid) to `source.p_mask` where the 1's represent the pixels that form part of the source. The time varying input signal is then assigned to `source.p`. This must be the same length as `kgrid.t_array` and can be a single time series (in which case it is applied to all source elements), or a matrix of time series following the source elements using MATLAB's standard column-wise linear matrix index ordering. Here a sinusoidal input is assigned to a square source element. To avoid numerical stabilities, the input should first be filtered using `filterTimeSeries` (see the Filtering A Delta Function Input Signal Example for more information).

    % create the time array
    [kgrid.t_array dt] = makeTime(kgrid, medium.sound_speed);

    % define a square source element
    source_radius = 5;
    source.p_mask = zeros(Nz, Nx, Ny);
    source.p_mask(Nz/4, Nx/2 - source_radius:Nx/2 + source_radius, Ny/2 - source_radius:Ny/2 + source_radius) = 1;

    % define a time varying sinusoidal source
    source_freq = 2e6;
    source_mag = 1;
    source.p = source_mag*sin(2*pi*source_freq*kgrid.t_array);

    % smooth the source
    source.p = filterTimeSeries(kgrid, medium, source.p);

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). To allow visualisation of the source elements within the grid, the source mask is assigned to the optional input `'DisplayMask'`. This mask is overlayed onto the plot during the simulation.

    % input arguments
    input_args = {'PlotLayout', true, 'PMLSize', 10, 'DisplayMask', source.p_mask};

    % run the simulation
    [sensor_data p_final] = kspaceFirstOrder3D(kgrid, medium, source, sensor, 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. A plot of display during simulation is shown below.

The effective visualisation of three-dimensional matrix data remains an important part of data exploration and presentation. An animated slice-by-slice visualisation of the final pressure field can be viewed using `flyThrough`.

    % view final pressure field slice by slice
    flyThrough(p_final);

Back to Top

© 2009, 2010 Bradley Treeby and Ben Cox.