k-Wave Toolbox |
![]() ![]() |
On this page… |
---|
Creating the k-space grid and defining the medium properties |
This example provides a simple demonstration of using k-Wave for the simulation and detection of the pressure field generated by an initial pressure distribution within a three-dimensional heterogeneous propagation medium. It builds on the Homogeneous Propagation Medium and Heterogeneous Propagation Medium examples.
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 length [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]
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); source.p0 = ball_1 + ball_2;
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 using cart2grid
and voxelPlot
is given below.
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.
% input arguments input_args = {'PlotLayout', true, 'PlotPML', false, 'DataCast', 'single'}; % run the simulation sensor_data = 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. To remove the PML from the display, the optional input 'PlotPML'
is also set to false
.
As the function runs, status updates and computational parameters are printed to the command line. By setting the optional input 'DataCast'
to 'single'
, the computational time is decreased (see the Optimising k-Wave Performance Example for more details).
Running k-space simulation... start time: 09-Feb-2011 14:00:32 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) maximum supported frequency: 7.5MHz smoothing p0 distribution... casting variables to single type... precomputation completed in 1.5767s starting time loop... estimated simulation time 22.291s... simulation completed in 26.1408s reordering Cartesian measurement data... total computation time 27.75s
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.
As with the two-dimensional simulations, the interpolation mode when using a Cartesian sensor mask may be set to 'linear'
or 'nearest'
(the default in three-dimensions). If 'CartInterp'
is set to 'linear'
, TriScatteredInterp
is used to compute the Delaunay triangulation during the precomputation (note, this was only introduced in MATLAB R2009a). For large 3D matrices, this can signficantly add to the precomputation time (see the output below for the same example run with 'CartInterp'
set to 'linear'
). Note, the inbuilt function TriScatteredInterp
only accepts inputs of type 'double'
, thus 'DataCast'
must be set to 'off'
(the default value) if using linear interpolation.
Running k-space simulation... start time: 09-Feb-2011 14:07:48 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) maximum supported frequency: 7.5MHz smoothing p0 distribution... calculating Delaunay triangulation (TriScatteredInterp)... precomputation completed in 27.6795s starting time loop... estimated simulation time 41.6013s... simulation completed in 47.4548s total computation time 1min 15.163s
In both two- and three- dimensional simulations, if the optional input 'CartInterp'
is set to 'nearest'
, 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.
![]() |
Simulations In One Dimension | Photoacoustic Waveforms in 1D, 2D and 3D | ![]() |
© 2009, 2010, 2011 Bradley Treeby and Ben Cox.