| 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 array of sensor elements. 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 used for the reconstruction is generated using kspaceFirstOrder2D in the same way as in the Initial Value Problems examples. (The reconstruction of experimental data can approached in the same fashion by substituting the numerical sensor data with experimental measurements.) It is important to note that, by default, the function kspaceFirstOrder2D spatially smooths the initial pressure distribution using smooth (see the Source Smoothing Example). 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 modified via the optional input parameter 'Smooth'. If given as a single Boolean value, this setting controls the smoothing of the initial pressure along with the sound speed and density distributions. To control the smoothing of these distributions individually, 'Smooth' can also be given as a three element array controlling the smoothing of the initial pressure, the sound speed, and the density, respectively.
Although kspaceFirstOrder2D
can be used to automatically define the simulation length and time-step (by default kgrid.t_array is set to 'auto' by makeGrid), this information is also required by kspaceLineRecon and thus in this example kgrid.t_array must be explicitly created. This can be easily achieved by directly calling makeTime. This is the same function that is called internally by the first-order simulation functions 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. A sensor with this 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 x = const).
% define a binary line sensor sensor.mask = zeros(Nx, Ny); sensor.mask(1, :) = 1;
The initial pressure distribution (along with the binary sensor mask) and the recorded sensor data returned after running the simulation are shown below.
The reconstruction is invoked by calling kspaceLineRecon with the recorded sensor data, as well as the properties of the acoustic medium and the sampling parameters. By default, the sensor data input must be indexed as p(time, sensor_position). This means the simulated sensor data returned by kspaceFirstOrder2D which is indexed as sensor_data(sensor_position, time) must first be rotated. Alternatively, the optional input parameter 'DataOrder' can be set to 'yt' (the default settings is 'ty'). By setting the optional input parameter 'Plot' to true, a plot of the reconstructed initial pressure distribution is also produced.
% reconstruct the initial pressure p_xy = kspaceLineRecon(sensor_data.', dy, dt, medium.sound_speed, 'Plot', true);
As the reconstruction runs, the size of the recorded data and the time to compute the reconstruction are both printed to the command line.
Running k-Wave line reconstruction... grid size: 216 by 778 grid points interpolation mode: *nearest computation completed in 0.11074s
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., x = const). The resolution of the reconstruction in the y-direction is defined by the physical location and spacing of the sensor elements while the resolution in the x-direction is defined by the sampling rate at which the pressure field is recorded (i.e., dt). Consequently, the reconstructed initial pressure map typically will have a much finer discretisation in the x- (time) direction. 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.
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_xy [Nx_recon, Ny_recon] = size(p_xy); kgrid_recon = makeGrid(Nx_recon, dt*medium.sound_speed, Ny_recon, dy); % resample p_xy to be the same size as source.p0 p_xy_rs = interp2(kgrid_recon.y, kgrid_recon.x - min(kgrid_recon.x(:)), p_xy, kgrid.y, kgrid.x - min(kgrid.x(:)));
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 x = 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 '*linear' or '*cubic', and re-running the simulation. This 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.
In practice, '*linear' interpolation provides a good balance between reconstruction speed and image artefacts.
| |
Photoacoustic Image Reconstruction | 3D FFT Reconstruction For A Planar Sensor | |
© 2009-2014 Bradley Treeby and Ben Cox.
