| 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.
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 p0 = smooth(p0, kgrid, 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 (by setting t_array to 'auto'), 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.
% create the time array [t_array dt] = makeTime(kgrid, c);
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 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*c); % resample p_zx to be the same size as 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 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 will thus be 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 Bradley Treeby and Ben Cox.