k-Wave Toolbox |

## 2D FFT Reconstruction For A Line Sensor Example

On this page…
Overview Defining the initial pressure distribution Defining the time array Simulating the sensor data Performing the reconstruction The limited view problem Setting the interpolation mode

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

- open the file in the MATLAB Editor
- run the file in MATLAB

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.