k-Wave Toolbox |

# kspaceFirstOrder2D

2D time-domain simulation of wave propagation

## Syntax

    sensor_data = kspaceFirstOrder2D(p0, kgrid, c, rho, t_array, sensor_mask)
    sensor_data = kspaceFirstOrder2D(p0, kgrid, c, rho, t_array, sensor_mask, ...)

## Description

`kspaceFirstOrder2D` simulates the time-domain propagation of linear compressional waves through a two-dimensional homogeneous or heterogeneous acoustic medium defined by `c` and `rho` given the initial pressure distribution `p0`. The size and discretisation of the acoustic domain are defined by the k-space grid structure `kgrid`. At each time-step the pressure at the positions defined by `sensor_mask` are recorded and stored.

The computation is based on a first-order k-space model which allows a heterogeneous sound speed and density. An absorbing boundary condition (in this case a perfectly matched layer) 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 (PML). 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.

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 the same size as `p0`) representing the pixels 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 2 by N matrix corresponding to the x and z positions, respectively.

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

The code may also be used for time-reversal image reconstruction by setting the optional input `'TimeRev'` to `true`. This enforces the pressure given by `p0` as a time varying Dirichlet boundary condition over the sensor mask. In this mode, the input pressure `p0` must be indexed as `p0(sensor position, time)`. If `sensor_mask` is given as a set of Cartesian coordinates then `p0` 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 then `p0` must be given as an array ordered using MATLAB's standard column-wise linear matrix indexing.

Note To run the inverse crime, the `sensor_data` returned when `'TimeRev'` is set to `false` can be used without modification as `p0` when `'TimeRev'` is set to `true`.

## Inputs

`p0` | map of the initial pressure within the medium over the discretisation given by `kgrid` (or the time varying pressure across the sensor mask if `'TimeRev'` is set to `true`)
`kgrid` | k-space grid structure returned by `makeGrid`
`c` | sound speed within the acoustic medium [m/s]
`rho` | density within the acoustic medium [kg/m^3]
`t_array` | evenly spaced array of time values [s] (`t_array` can alternatively be set to `'auto'` to automatically generate the array using `makeTime`)
`sensor_mask` | binary grid or a set of Cartesian points where the pressure is recorded at each time-step

Note For heterogeneous medium parameters, `c` and `rho` must be given in matrix form with the same dimensions as `p0`. For homogeneous medium parameters, `c` and `rho` can be given as single numeric values.

## Optional Inputs

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

Input | Valid Settings | Default | Description
--- | --- | --- | ---
`'AdaptThresh'` | (numeric scalar) | `0.005` | Adaptive boundary condition threshold used when `'TimeRev'` is set to `3`.
`'CartInterp'` | `'linear'` `'nearest'` | `'linear'` | 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 pixel, duplicated data points are discarded and `sensor_data` will be returned with less points than that specified by `sensor_mask`.
`'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 `fft2` and `ifft2` 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'`.
`'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 a four panel plot of the initial simulation layout is produced (initial pressure, sensor mask, sound speed, density).
`'PlotScale'` | (numeric two element vector) | `[-1, 1]` | [min, max] values used to control the scaling for `imagesc` (visualisation) and `im2frame` (movie capture).
`'PlotSim'` | (boolean scalar) | `true` | Boolean controlling whether the simulation iterations are progressively plotted.
`'PMLAlpha'` | (numeric scalar or two element vector) | `4` | Attenuation in Nepers per m of the absorption within the perfectly matched layer.
`'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 two element vector) | `20` | Size of the perfectly matched layer in pixels. By default, the PML is added evenly to all sides of the grid, however, both `PMLSize` and `PMLAlpha` can be given as two element arrays to specify the x and z properties, respectively. To remove the PML, set the appropriate `PMLAlpha` to zero rather than forcing the PML to be of zero size.
`'MovieName'` | (string) | `'date-time'` | Name of the movie produced when `'RecordMovie'` is set to `true`.
`'MovieArgs'` | (string cell array) | `{}` | Settings for `movie2avi`. Parameters must be given as {param, value, ...} pairs within a cell array.
`'RecordMovie'` | (boolean scalar) | `false` | Boolean controlling whether the displayed image frames are captured using `im2frame` and stored as a movie using `movie2avi`.
`'Smooth'` | (boolean scalar or three element vector) | `true` | Boolean controlling whether `p0`, `c`, and `rho` are smoothed using `smooth` before computation. `'Smooth'` can also be given as a 3 element array to control the smoothing of `p0`, `c`, and `rho`, respectively.
`'TimeRev'` | (boolean scalar or integer numeric scalar) | `false` | Boolean controlling whether the code is used in time-reversal mode. If set to `true` (or `1`), the time-reversal is computed by enforcing the pressure values given by `p0` over the sensor surface at each time-step as a time varying Dirichlet boundary condition (conventional time-reversal). If set to `2`, the time-reversal is computed by introducing the pressure values given by `p0` over the sensor surface as a source term. If set to `3`, the time-reversal is computed by using an adaptive boundary condition with the threshold set by `'AdaptThresh'`.

## Outputs

`sensor_data` | array of pressure time-series recorded at the sensor positions given by `sensor_mask`

## Examples

- Homogeneous Propagation Medium Example
- Using A Binary Sensor Mask Example
- Loading External Image Maps Example
- Heterogeneous Propagation Medium Example
- Saving Movie Files Example
- Controlling The Absorbing Boundary Layer Example
- 2D FFT Reconstruction For A Line Sensor Example
- 2D time-reversal Reconstruction For A Line Sensor
- 2D time-reversal Reconstruction For A Circular Sensor

## See Also

`fft2`, `ifft2`, `im2frame`, `imagesc`, `kspaceFirstOrder1D`, `kspaceFirstOrder3D`, `makeGrid`, `makeTime`, `movie2avi`, `smooth`, `unmaskSensorData`

© 2009 Bradley Treeby and Ben Cox.