k-Wave Toolbox Previous   Next

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 p_final] = kspaceFirstOrder3D(kgrid, medium, source, sensor)
[sensor_data p_final] = kspaceFirstOrder3D(kgrid, medium, source, sensor, ...) 

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 allows a heterogeneous sound speed and density and power law absorption. 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 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.

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 the grid defined by kgrid) 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 3 by N matrix corresponding to the x, y, 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 assigning the boundary data to the input field sensor.time_reversal_boundary_data. This data is then enforced as a time varying Dirichlet boundary condition over the 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 then 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 then the boundary data must be given as an array ordered using MATLAB's standard column-wise linear matrix indexing.

Inputs

The fields that are required for a photoacoustic forward simulation are marked with a *.

kgrid*

k-space 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.density*

density distribution within the acoustic medium [kg/m^3]

medium.alpha_power

power law absorption exponent [dB/(MHz^y cm)]

medium.alpha_coeff

power law absorption coefficient [dB/(MHz^y cm)]

 

 

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

 

 

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

adaptive boundary condition threshold

Optional Inputs

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

Input Valid Settings Default Description

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

'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)

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.

'MovieName'

(string)

'date-time-kspaceFirstOrder3D'

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.

'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 three plots are produced of the initial simulation layout (initial pressure, sound speed, density).

'PlotScale'

(numeric two element vector)

[-1, 1]

[min, max] values used to control the scaling for imagesc (visualisation).

'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. Currently only true is supported in 3D.

'PMLSize'

(integer numeric scalar or three 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 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 using getframe and stored as a movie using movie2avi.

'Smooth'

(boolean scalar or three element vector)

true

Boolean controlling whether source.p0, medium.sound_speed, and medium.density are smoothed using smooth before computation. 'Smooth' can also be given as a 3 element array to control the smoothing of source.p0, medium.sound_speed, and medium.density, respectively.

Outputs

sensor_data

array of pressure time-series recorded at the sensor positions given by sensor.mask

p_final

final pressure field

Examples

See Also

fftn, ifftn, imagesc, kspaceFirstOrder1D, kspaceFirstOrder2D, makeGrid, makeTime, smooth, unmaskSensorData


© 2009, 2010 Bradley Treeby and Ben Cox.