2. Comparing image reconstruction methods

In this section, we will compare several different image reconstruction methods. Several different methods are available in photoacoustic imaging. Each differs in its complexity, speed and physical basis.

PATATO includes two algorithms for image reconstruction, backprojection and a model-based method. Here we will compare the two.

import matplotlib.pyplot as plt

import patato as pat

As an example dataset, we will use a dataset taken from a small-animal photoacoustic scanner (iThera inVision 256-TF, iThera Medical GmbH). We scanned a single slice through a tumour-bearing mouse. Example data can be found in our example data repository here.

pa_data = pat.PAData.from_hdf5("data/so2-timeseries-data.hdf5")[0:1]

First we will define a reconstruction grid to directly compare the reconstruction results.

# Define the reconstruction methods to compare:
nx = 101 # number of pixels
lx = 2.5e-2 # m

We will then set up our backprojection reconstruction pipeline.

# Define a pre-processor for backprojection
pre_bp = pat.MSOTPreProcessor(lp_filter=7e6, hp_filter=5e3)

ts_bp, settings_bp, _ = pre_bp.run(pa_data.get_time_series(), pa_data)


# Delay and sum.
das = pat.Backprojection(field_of_view=(lx, lx, 0),
                             n_pixels=(nx, nx, 1)) # z axis must be specified but is empty in this case.

Next, we define out model-based reconstruction pipeline. Here, the program pre-calculates the forward acoustic model matrix, which is then applied in the following step.

# Apply energy correction to time series data for Model Based reconstruction:
ts_model_based = pa_data.get_time_series().copy()
ts_model_based.raw_data = ts_model_based.raw_data / pa_data.get_overall_correction_factor()[:, :, None, None]

model_based = pat.ModelBasedReconstruction(field_of_view=(lx, lx, 0),
                                           n_pixels=(nx, nx, 1), pa_example=pa_data)

Loading model

/home/cri.camres.org/else01/patato/patato/recon/model_based/model_based.py:219: UserWarning: This class is experimental. If you would like to contribute to further development of this approach, please get in touch.
  warnings.warn("This class is experimental. If you would like to contribute to further development of this "

Now, do the image reconstruction.

rec_backprojection, _, _ = das.run(ts_bp, pa_data, **settings_bp)
rec_modelbased, _, _ = model_based.run(ts_model_based, pa_data)

As before, we can plot the reconstructed images from the two approaches and compare them visually.

fig, axes = plt.subplots(1, 2)
rec_modelbased.imshow(ax = axes[0])
axes[0].set_title("Model based")
rec_backprojection.imshow(ax = axes[1])
axes[1].set_title("Backprojection")
plt.show()