k-Wave Toolbox |

## 2D FFT Reconstruction For A Line Sensor Example

On this page…
Overview Defining a smooth 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 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.

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

Back to Top

### Defining a smooth initial pressure distribution

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.

Back to Top

### Defining the time array

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

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

Back to Top

### Performing the reconstruction

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

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

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

Back to Top

© 2009-2014 Bradley Treeby and Ben Cox.