k-Wave Toolbox Previous   Next

2D FFT Reconstruction For A Line Sensor Example

Overview

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.

 Back to Top

Defining the initial pressure distribution

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

 Back to Top

Defining the time array

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

 Back to Top

Simulating the sensor data

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.

 Back to Top

Performing the reconstruction

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.

 Back to Top

The limited view problem

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.

 Back to Top

Setting the interpolation mode

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.

 Back to Top


© 2009 Bradley Treeby and Ben Cox.