k-Wave Toolbox |
![]() ![]() |
On this page… |
---|
This example demonstrates the use of k-Wave for the reconstruction of a two-dimensional photoacoustic wave-field recorded over a linear sensor array. The sensor data is simulated using kspaceFirstOrder2D
and reconstructed using kspaceLineRecon
. It builds on the Homogeneous Propagation Medium and Heterogeneous Propagation Medium examples.
The sensor data is simulated using kspaceFirstOrder2D
in the same way as in the preceeding simulation examples. By default, this
function applies a k-space Hanning filter (using smooth
)
to both the initial pressure, and the sound speed and density matrices (if they are heterogeneous).
This filtering removes high spatial frequencies (generated by sharp edges)
that cause numerical aliasing in the simulation (see the Filtering A Delta Function Input Signal Example
for more information).
To provide a valid comparison, the output of the reconstruction code kspaceLineRecon
must be compared to the smoothed version of the initial pressure used in the simulation. This can be achieved by explicitly using smooth
before passing the initial pressure to kspaceFirstOrder2D
.
% smooth the initial pressure distribution and restore the magnitude source.p0 = smooth(kgrid, source.p0, true);
The default smoothing within kspaceFirstOrder2D
can be switched off by setting the optional input 'Smooth'
to false
. To control the smoothing of the initial pressure and heterogeneity maps individually, 'Smooth'
can also be given as a three element array (controlling smoothing of the initial pressure, the heterogeneous sound speed map, and the heterogeneous density map, respectively).
Although kspaceFirstOrder2D
can be used to automatically define the simulation length and time-step
(kgrid.t_array
is set to 'auto'
by makeGrid
), this information
is also required by kspaceLineRecon
and thus in this example, t_array
must be explicitly created.
This can be easily achieved by directly calling makeTime
,
the same function that is called internally from within the first-order simulation codes when
kgrid.t_array
is set to 'auto'
.
% create the time array [kgrid.t_array dt] = makeTime(kgrid, medium.sound_speed);
The FFT reconstruction function kspaceLineRecon
requires data recorded along an equally spaced line-shaped array of sensor points. This sensor shape can be created by defining a binary sensor mask matrix with a line of 1's along the first matrix row (i.e., a line along z = const).
% define a binary line sensor sensor.mask = zeros(kgrid.Nz, kgrid.Nx); sensor.mask(1, :) = 1;
The initial pressure distribution (along with the binary sensor mask) and the recorded time series are shown below.
The reconstruction is invoked by calling kspaceLineRecon
with the recorded time-series data, as well as the properties of the acoustic medium and the sampling parameters. The time-series data input must be indexed as p_tx(time, sensor position)
so the simulated sensor_data
returned by kspaceFirstOrder2D
must first be rotated. By setting the optional input parameter 'PlotRecon'
to true
, a plot of the reconstructed initial pressure distribution is also produced. Regardless of the physical alignment of the sensor within the acoustic medium, the reconstruction is always returned as if the sensor was located along the first matrix row (i.e., z = const).
% reconstruct the initial pressure p_zx = kspaceLineRecon(sensor_data.', kgrid.dx, dt, c, 'PlotRecon', true);
The size of the recorded data and the time to compute the reconstruction are both printed to the command line.
Running k-space line reconstruction... grid size: 216 by 778 pixels interpolation mode: *nearest computation completed in 0.18423s
By default, the reconstructed initial pressure distribution is not re-scaled or thresholded in any way. However, a positivity condition can be automatically enforced by setting the optional input parameter 'PosCond'
to true
.
The resolution of the reconstruction in the x-direction is defined by the physical location and spacing of the sensor elements. The resolution in the z-direction is defined by sample rate (i.e., dt) at which the pressure field is recorded. The reconstructed initial pressure map will thus typically have a much finer discretisation in the z- (time) direction.
To directly compare the initial pressure distribution with that produced from the reconstruction, it is convenient to interpolate the reconstructed pressure onto a k-space grid with the same dimensions as source.p0
.
% define a second k-space grid using the dimensions of p_zx [Nz_recon Nx_recon] = size(p_zx); kgrid_recon = makeGrid(Nx_recon, kgrid.dx, Nz_recon, dt*medium.sound_speed); % resample p_zx to be the same size as source.p0 p_zx_rs = interp2(kgrid_recon.x_off, kgrid_recon.z_off, p_zx, kgrid.x_off, kgrid.z_off);
The interpolated reconstructed initial pressure distribution (using a positivity condition) with the same plot scaling as the plot of source.p0
above is shown below. A z = const slice through the center of the larger disc is also shown for comparison. The reconstructed pressure magnitude is decreased due to the limited view problem; the reconstruction is only exact if the data is collected over an infinite line, here a finite length sensor is used.
The FFT reconstruction algorithm relies on the interpolation between a temporal and a spatial domain coordinate with different inherent spacings. Both the speed and accuracy of the reconstruction is thus dependent on the method used for this interpolation. This can be controlled via the optional input parameter 'Interp'
which is passed directly to interp2
. By default, this is set to '*nearest'
which optimises the interpolation for speed. The accuracy of the interpolation can be improved by setting 'Interp'
to '*cubic'
(for example), and re-running the simulation.
Running k-space line reconstruction... grid size: 216 by 778 pixels interpolation mode: *cubic applying positivity condition... computation completed in 0.595s
This particular interpolation setting considerably increases the time to compute the reconstruction, however, the noise in the reconstruction is noticeably improved, and the error in the magnitude of the reconstruction is also slightly reduced.
![]() |
Photoacoustic Image Reconstruction | 3D FFT Reconstruction For A Planar Sensor | ![]() |
© 2009, 2010 Bradley Treeby and Ben Cox.