k-Wave Toolbox |

# kspaceFirstOrder3D

3D time-domain simulation of wave propagation

## Syntax

    sensor_data = kspaceFirstOrder3D(kgrid, medium, source, sensor)
    sensor_data = kspaceFirstOrder3D(kgrid, medium, source, sensor, ...)

    [sensor_data, field_data] = kspaceFirstOrder3D(kgrid, medium, source, sensor)
    [sensor_data, field_data] = kspaceFirstOrder3D(kgrid, medium, source, sensor, ...)

    kspaceFirstOrder2D(kgrid, medium, source, [])
    kspaceFirstOrder2D(kgrid, medium, source, [], ...)

## Description

`kspaceFirstOrder3D` simulates the time-domain propagation of linear compressional waves through a three-dimensional homogeneous or heterogeneous acoustic medium given four input structures: `kgrid`, `medium`, `source`, and `sensor`. The computation is based on a first-order k-space model which accounts for power law absorption and a heterogeneous sound speed and density. If `medium.BonA` is specified, cumulative nonlinear effects are also modelled. At each time-step (defined by `kgrid.t_array`), the pressure at the positions defined by `sensor.mask` are recorded and stored. If `kgrid.t_array` is set to `'auto'`, this array is automatically generated using `makeTime`. An anisotropic absorbing boundary condition called a perfectly matched layer (PML) is implemented to prevent waves that leave one side of the domain being reintroduced from the opposite side (a consequence of using the FFT to compute the spatial derivatives in the wave equation). This allows infinite domain simulations to be computed using small computational grids.

For a homogeneous medium the formulation is exact and the time-steps are only limited by the effectiveness of the perfectly matched layer. For a heterogeneous medium, the solution represents a leap-frog pseudospectral method with a Laplacian correction that improves the accuracy of computing the temporal derivatives. This allows larger time-steps to be taken without instability compared to conventional pseudospectral time-domain methods. The computational grids are staggered both spatially and temporally.

An initial pressure distribution can be specified by assigning a matrix (the same size as the computational grid) of arbitrary numeric values to `source.p0`. A time varying pressure source can similarly be specified 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 grid points that form part of the source. The time varying input signals are 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. A time varying velocity source can be specified in an analogous fashion, where the source location is specified by `source.u_mask`, and the time varying input velocity is assigned to `source.ux`, `source.uy`, and `source.uz`.

The pressure is returned as an array of time series at the sensor locations defined by `sensor.mask`. This can be given either as a binary grid (i.e., a matrix of 1's and 0's with the same dimensions as the computational grid) representing the grid points within the computational grid that will collect the data, or as a series of arbitrary Cartesian coordinates within the grid at which the pressure values are calculated at each time step via interpolation. The Cartesian points must be given as a 3 by N matrix corresponding to the x, y, and z positions, respectively. The final pressure field over the complete computational grid can also be obtained using the output `field_data`. If no output is required, the `sensor` input can be replaced with an empty array `[]`. Both the `source` and `sensor` inputs can also be replaced by an object of the `kWaveTransducer` class created using `makeTransducer`

If `sensor.mask` is given as a set of Cartesian coordinates, the computed `sensor_data` is returned in the same order. If `sensor.mask` is given as a binary grid, `sensor_data` is returned using MATLAB's standard column-wise linear matrix index ordering. In both cases, the recorded data is indexed as `sensor_data(sensor_position, time)`. For a binary sensor mask, the pressure values at a particular time can be restored to the sensor positions within the computation grid using `unmaskSensorData`.

By default, the recorded pressure field is passed directly to the output arguments `sensor_data` and `field_data`. However, the particle velocity can also be recorded by setting the optional input parameter `'ReturnVelocity'` to `true`. In this case, the output arguments `sensor_data` and `field_data` are returned as structures with the pressure and particle velocity appended as separate fields. In two dimensions, these fields are given by `sensor_data.p`, `sensor_data.ux`, `sensor_data.uy`, and `sensor_data.uz`, and `field_final.p`, `field_final.ux`, `field_final.uy`, and `field_final.uz`.

`kspaceFirstOrder3D` may also be used for time reversal image reconstruction by assigning the time varying pressure recorded over an arbitrary sensor surface to the input field `sensor.time_reversal_boundary_data`. This data is then enforced in time reversed order as a time varying Dirichlet boundary condition over the sensor surface given by `sensor.mask`. The boundary data must be indexed as `sensor.time_reversal_boundary_data(sensor_position, time)`. If `sensor.mask` is given as a set of Cartesian coordinates, the boundary data must be given in the same order. An equivalent binary sensor mask (computed using nearest neighbour interpolation) is then used to place the pressure values into the computational grid at each time step. If `sensor.mask` is given as a binary grid of sensor points, the boundary data must be ordered using MATLAB's standard column-wise linear matrix indexing. If no additional inputs are required, the `source` input can be replaced with an empty array `[]`.

Acoustic attenuation compensation can also be included during time reversal image reconstruction by assigning the absorption parameters `medium.alpha_coeff` and `medium.alpha_power` and reversing the sign of the absorption term by setting `medium.alpha_sign = [-1, 1]`. This forces the propagating waves to grow according to the absorption parameters instead of decay. The reconstruction should then be regularised by assigning a filter to `medium.alpha_filter` (this can be created using `getAlphaFilter`).

Note To run a simple reconstruction example using time reversal (that commits the 'inverse crime' of using the same numerical parameters and model for data simulation and image reconstruction), the `sensor_data` returned from a k-Wave simulation can be passed directly to `sensor.time_reversal_boundary_data` with the input fields `source.p0` and `source.p` removed or set to zero.

## Inputs

The minimum fields that must be assigned to run an initial value problem (for example, a photoacoustic forward simulation) are marked with a *.

`kgrid*` | k-Wave grid structure returned by `makeGrid` containing Cartesian and k-space grid fields
`kgrid.t_array*` | evenly spaced array of time values [s] (set to `'auto'` by `makeGrid`)
`medium.sound_speed*` | sound speed distribution within the acoustic medium [m/s]
`medium.sound_speed_ref` | reference sound speed used within the k-space operator (phase correction term) [m/s]
`medium.density*` | density distribution within the acoustic medium [kg/m^3]
`medium.BonA` | parameter of nonlinearity
`medium.alpha_power` | power law absorption exponent
`medium.alpha_coeff` | power law absorption coefficient [dB/(MHz^y cm)]
`medium.alpha_mode` | optional input to force either the absorption or dispersion terms in the equation of state to be excluded; valid inputs are `'no_absorption'` or `'no_dispersion'`
`medium.alpha_filter` | frequency domain filter applied to the absorption and dispersion terms in the equation of state
`medium.alpha_sign` | two element array used to control the sign of absorption and dispersion terms in the equation of state
`source.p0*` | initial pressure within the acoustic medium
`source.p` | time varying pressure at each of the source positions given by `source.p_mask`
`source.p_mask` | binary grid specifying the positions of the time varying pressure source distribution
`source.p_mode` | optional input to control whether the input pressure is injected as a mass source or enforced as a dirichlet boundary condition; valid inputs are `'additive'` (the default) or `'dirichlet'`
`source.ux` | time varying particle velocity in the x-direction at each of the source positions given by `source.u_mask`
`source.uy` | time varying particle velocity in the y-direction at each of the source positions given by `source.u_mask`
`source.uz` | time varying particle velocity in the z-direction at each of the source positions given by `source.u_mask`
`source.u_mask` | binary grid specifying the positions of the time varying particle velocity distribution
`source.u_mode` | optional input to control whether the input velocity is applied as a force source or enforced as a dirichlet boundary condition; valid inputs are `'additive'` (the default) or `'dirichlet'`
`sensor.mask*` | binary grid or a set of Cartesian points where the pressure is recorded at each time-step
`sensor.time_reversal_boundary_data` | time varying pressure enforced as a Dirichlet boundary condition over `sensor.mask`
`sensor.frequency_response` | two element array specifying the center frequency and percentage bandwidth of a frequency domain Gaussian filter applied to the `sensor_data`
`sensor.record_mode` | optional input to save the statistics of the wave field at the sensor elements rather than the time series; valid inputs are `'statistics'` or `'time_history'` (the default)

Note For heterogeneous medium parameters, `medium.sound_speed` and `medium.density` must be given in matrix form with the same dimensions as `kgrid`. For homogeneous medium parameters, these can be given as single numeric values. If the medium is homogeneous and velocity inputs or outputs are not required, it is not necessary to specify `medium.density`.

## Optional Inputs

Optional 'string', value pairs that may be used to modify the default computational settings.

Input | Valid Settings | Default | Description
--- | --- | --- | ---
`'CartInterp'` | `'linear'` `'nearest'` | `'nearest'` | Interpolation mode used to extract the pressure when a Cartesian sensor mask is given. If set to `'nearest'` and more than one Cartesian point maps to the same grid point, duplicated data points are discarded and `sensor_data` will be returned with less points than that specified by `sensor.mask`.
`'CreateLog'` | (Boolean scalar) | `false` | Boolean controlling whether the command line output is saved using the diary function with a date and time stamped filename.
`'DataCast'` | (string of data type) | `'off'` | String input of the data type that variables are cast to before computation. For example, setting to `'single'` will speed up the computation time (due to the improved efficiency of `fftn` and `ifftn` for this data type) at the expense of a loss in precision. This variable is also useful for utilising GPU parallelisation through libraries such as GPUmat or AccelerEyesJacket by setting `'DataCast'` to `'GPUsingle'` or `'gsingle'`.
`'DisplayMask'` | (binary matrix) or `'off'` | `sensor.mask` | Binary matrix overlayed onto the animated simulation display. Elements set to 1 within the display mask are set to black within the display.
`'LogScale'` | (Boolean scalar) or (numeric scalar) | `false` | Boolean controlling whether the pressure field is log compressed before display. The data is compressed by scaling both the positive and negative values between 0 and 1 (truncating the data to the given plot scale), adding a scalar value (compression factor) and then using the corresponding portion of a log10 plot for the compression (the negative parts are remapped to be negative thus the default color scale will appear unchanged). The amount of compression can be controlled by adjusting the compression factor which can be given in place of the Boolean input. The closer the compression factor is to zero, the steeper the corresponding part of the log10 plot used, and the greater the compression (the default compression factor is 0.02).
`'MovieArgs'` | (string cell array) | `{}` | Settings for `movie2avi`. Parameters must be given as {param, value, ...} pairs within a cell array.
`'MovieName'` | (string) | `'date-time-kspaceFirstOrder2D'` | Name of the movie produced when `'RecordMovie'` is set to `true`.
`'PlotFreq'` | (integer numeric scalar) | `10` | The number of iterations which must pass before the simulation plot is updated.
`'PlotLayout'` | (Boolean scalar) | `false` | Boolean controlling whether three plots are produced of the initial simulation layout (initial pressure, sound speed, density).
`'PlotPML'` | (Boolean scalar) | `true` | Boolean controlling whether the perfectly matched layer is shown in the simulation plots. If set to `false`, the PML is not displayed.
`'PlotScale'` | (numeric two element vector) or `'auto'` | `[-1, 1]` | [min, max] values used to control the scaling for `imagesc` (visualisation). If set to `'auto'`, a symmetric plot scale is chosen automatically for each plot frame.
`'PlotSim'` | (Boolean scalar) | `true` | Boolean controlling whether the simulation iterations are progressively plotted.
`'PMLAlpha'` | (numeric scalar or three element vector) | `2` | Absorption within the perfectly matched layer in Nepers per metre.
`'PMLInside'` | (Boolean scalar) | `true` | Boolean controlling whether the perfectly matched layer is inside or outside the grid. If set to `false`, the input grids are enlarged by `PMLSize` before running the simulation.
`'PMLSize'` | (integer numeric scalar or three element vector) | `10` | Size of the perfectly matched layer in grid points. By default, the PML is added evenly to all sides of the grid, however, both `PMLSize` and `PMLAlpha` can be given as three element arrays to specify the x, y, and z properties, respectively. To remove the PML, set the appropriate `PMLAlpha` to zero rather than forcing the PML to be of zero size.
`'RecordMovie'` | (Boolean scalar) | `false` | Boolean controlling whether the displayed image frames are captured and stored as a movie using `movie2avi`.
`'ReturnVelocity'` | (Boolean scalar) | `false` | Boolean controlling whether the acoustic particle velocity at the positions defined by `sensor.mask` are also returned. If set to `true`, the output argument `sensor_data` is returned as a structure with the pressure and particle velocity appended as separate fields. In three dimensions, these fields are given by `sensor_data.p`, `sensor_data.ux`, `sensor_data.uy`, and `sensor_data.uz`. The field data is similarly returned as `field_data.p`, `field_data.ux`, `field_data.uy`, and `field_data.uz`.
`'Smooth'` | (Boolean scalar or three element vector) | `[true, false, false]` | Boolean controlling whether `source.p0`, `medium.sound_speed`, and `medium.density` are smoothed using `smooth` before computation. `'Smooth'` can either be given as a single Boolean value or as a 3 element array to control the smoothing of `source.p0`, `medium.sound_speed`, and `medium.density`, independently.
`'StreamToDisk'` | (Boolean scalar or integer numeric scalar) | `false` | Boolean controlling whether `sensor_data` is periodically saved to disk to avoid storing the complete matrix in memory. `StreamToDisk` may also be given as an integer which specifies the number of times steps that are taken before the data is saved to disk (default = 200).

## Outputs

If `'ReturnVelocity'` is set to `false`:

`sensor_data` | time varying pressure recorded at the sensor positions given by `sensor.mask`
`field_data` | final pressure field

If `'ReturnVelocity'` is set to `true`:

`sensor_data.p` | time varying pressure recorded at the sensor positions given by `sensor.mask`
`sensor_data.ux` | time varying particle velocity in the x-direction recorded at the sensor positions given by `sensor.mask`
`sensor_data.uy` | time varying particle velocity in the y-direction recorded at the sensor positions given by `sensor.mask`
`sensor_data.uz` | time varying particle velocity in the z-direction recorded at the sensor positions given by `sensor.mask`
`field_data.p` | final pressure field
`field_data.ux` | final field of the particle velocity in the x-direction
`field_data.uy` | final field of the particle velocity in the y-direction
`field_data.uz` | final field of the particle velocity in the z-direction

## Examples

- Simulations In Three Dimensions
- Photoacoustic Waveforms in 1D, 2D and 3D
- Simulations In Three-Dimensions
- Focussed Detector in 3D
- Modelling Sensor Directivity in 3D
- 3D FFT Reconstruction For A Planar Sensor
- 3D Time Reversal For A Planar Sensor
- 3D Time Reversal For A Spherical Sensor
- Defining An Ultrasound Transducer
- Simulating Ultrasound Beam Patterns
- Using An Ultrasound Transducer As A Sensor
- Simulating B-mode Ultrasound Images
- Optimising k-Wave Performance

## See Also

`benchmark`, `fftn`, `ifftn`, `imagesc`, `kspaceFirstOrder1D`, `kspaceFirstOrder2D`, `makeGrid`, `makeTime`, `makeTransducer`, `smooth`, `unmaskSensorData`

© 2009-2012 Bradley Treeby and Ben Cox.