# Iterative Image Improvement Using Time Reversal Example

This example demonstrates how photoacoustic image reconstruction may be improved iteratively using time reversal and a positivity condition. This example builds on the 2D Time Reversal Reconstruction For A Line Sensor Example.

- Open the file in the MATLAB Editor
- Run the file in MATLAB

## Contents

- Simulating the sensor data
- Image reconstruction using time reversal
- Iterating to improve the image reconstruction

## Simulating the sensor data

The example begins by defining an initial pressure distribution using a pre-defined image of the word 'k-Wave'. The acoustic pressure time series measured at each of the sensor positions in an L-shaped sensor array are then simulated using `kspaceFirstOrder2D`. The reason for using an L-shaped detector is so there will be a 'visible region' (the convex hull) in which the iterations below will converge to the true solution.

    % load an image for the initial pressure distribution
    p0_image = loadImage('EXAMPLE_k-Wave.png');

    % make it binary
    p0_image = double(p0_image > 0);

    % smooth the initial pressure distribution
    p0 = smooth(p0_image, true);

    % assign to the source structure
    source.p0 = p0;

## Image reconstruction using time reversal

The simulated data is then used in a time reversal image reconstruction. One way to setup time reversal is described in 2D Time Reversal Reconstruction For A Line Sensor Example. Here, a different way to setup time reversal is described, in which the time series are explicitly reversed in time and the `source.p_mode = 'dirichlet'` option is used to enforce the values at the sensor points.

    % use the sensor points as sources in time reversal
    source.p_mask = sensor.mask;

    % time reverse and assign the data
    source.p = fliplr(sensor_data);

    % enforce, rather than add, the time-reversed pressure values
    source.p_mode = 'dirichlet';

    % set the simulation to record the final image (at t = 0)
    sensor.record = {'p_final'};

    % run the time reversal reconstruction
    p0_estimate = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});

    % apply a positivity condition
    p0_estimate.p_final = p0_estimate.p_final .* (p0_estimate.p_final > 0);

Note that a positivity condition is applied to the reconstructed image because, from physical considerations, the initial acoustic pressure distribution must be zero or greater. The initial pressure distribution and image resulting from this first reconstruction are shown below.

## Iterating to improve the image reconstruction

Beginning with the image obtained using time reversal, an iterative procedure is then used to improve it. First, the image from the first step is used as a photoacoustic source in a forward model. The resulting time series are subtracted from the original measured time series to give a residual.

    % set the initial pressure to be the latest estimate of p0
    source.p0 = p0_estimate.p_final;

    % set the simulation to record the time series
    sensor = rmfield(sensor, 'record');

    % calculate the time series using the latest estimate of p0
    sensor_data2 = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});

    % calculate the error in the estimated time series
    data_difference = sensor_data - sensor_data2;

An image is then formed from this residual data using time reversal reconstruction in the usual way. The resulting image is then added to the original image to give an improved estimate.

    % assign the data_difference as a time-reversal source
    source.p_mask = sensor.mask;
    source.p = fliplr(data_difference);
    source = rmfield(source,'p0');
    source.p_mode = 'dirichlet';

    % set the simulation to record the final image (at t = 0)
    sensor.record = {'p_final'};

    % run the time reversal reconstruction
    p0_update = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});

    % add the update to the latest image
    p0_estimate.p_final = p0_estimate.p_final + p0_update.p_final;

    % apply a positivity condition
    p0_estimate.p_final = p0_estimate.p_final .* (p0_estimate.p_final > 0);

This simulation-reconstruction loop can then be repeated until there is no further noticeable improvement in the image. The image after three iterations is shown below.

Notice that within the region in which the visibility condition is satisfied (above the dashed line) the amplitude has recovered the correct value and the artefacts have reduced, but outside this region the amplitude remains low and there are still artefacts. This could be improved by recording data in more positions around the object, thereby increasing the visible region, or making other a priori assumptions about the image. In the latter case, an alternative approach to image reconstruction may be required. See Iterative Image Reconstruction Using the Adjoint Example.