PATLAB: A graphical computational software package for photoacoustic computed tomography research.

Photoacoustic tomography (PAT) provides high resolution optical images of tissue at depths of up to several centimetres. This modality has been of interest to researchers for at least 30 years and is still the subject of intensive research. However, PAT researchers lack access to a comprehensive open-source graphical simulation and reconstruction software package. In this article, we introduce PATLAB, an open-source MATLAB-based graphical software package that can perform both PAT simulation and image reconstruction. PATLAB is simple to use yet is capable of complex PAT data processing tasks and offers advanced users a framework to build and test new methods.

Introduction

Photoacoustic computed tomography (PAT) is an imaging modality with the ability to provide in-depth visualization of soft tissue [1], [2]. A typical PAT system uses light pulses to illuminate and penetrate the tissue. In response, absorbers present in the tissue (e.g. blood-filled vessels) produce acoustic signals. The signals are captured with detectors, recorded in digital format, processed, and reconstructed into PAT images [3]. The hybrid nature of PAT imparts the advantages of high resolution due to acoustic readout and high contrast due to the differences in optical absorption between tissue types [4]. As with any tomographic imaging modality, the image reconstruction process plays a key role. There have been many developments in PAT reconstruction algorithms over the past three decades mostly adapted from other imaging modalities such as ultrasound computed tomography (USCT) and x-ray computed tomography (CT). The algorithms have evolved over time with advancements such as signal denoising, accounting for spatial and temporal responses of detectors, handling of incomplete data, faster computation, and accounting for acoustically heterogeneous media [5]. Researchers seeking reconstruction software to evaluate hardware systems often face challenges, which include finding an appropriate algorithm, adapting the algorithm to their imaging system, and visualizing the data and images. It can be extremely time-consuming and difficult to meet these challenges, especially for researchers with no prior experience with image reconstruction. Fortunately, in recent years there has been significant development of software packages that have been made openly available to PAT researchers. For example, k-Wave is a MATLAB toolbox capable of performing photoacoustic simulation and image reconstruction [6]. Recently, another toolbox has been published for reconstruction and analysis of multispectral PAT [7]. These two toolboxes are both open source and offer researchers a comprehensive set of PAT data and image processing tools. However, each toolbox lacks a graphical interface and requires significant coding experience to use. Another published software package, 3D PHOVIS, incorporates a graphical interface, however it is designed for photoacoustic microscopy and lacks many of the tools needed to perform tomographic image reconstruction. Therefore, our goal was to develop an easy-to-use software package for PAT data processing and image reconstruction that does not require prior coding experience and allows users to work within a familiar graphical interface. The resultant PATLAB software was developed within the MATLAB environment, which we believe will motivate others to actively extend the functionality over time.

PATLAB is a simulation and reconstruction tool for PAT experiments that explore the development of PAT systems and PAT reconstruction algorithms. The goal of PATLAB was to provide a GUI-based software tool capable of reconstructing both two-dimensional (2D) and three-dimensional (3D) images from sets of photoacoustic time series. PATLAB includes an intuitive graphical interface and requires no prior coding experience. PATLAB contains a multitude of algorithms to simulate, pre-process, and reconstruct tomographic photoacoustic data. PATLAB has methods for loading raw data in various formats, for examining and trimming the data, for image reconstruction, and for image visualization. Furthermore, PATLAB provides a graphical user interface for interacting with the popular k-Wave toolbox, so that users can simulate PAT signals without coding. The PATLAB application and code is available on GitHub as an open-source project and includes a step-by-step user manual.

Methods

PATLAB was written in the popular MATLAB numerical environment. The 'App Designer' toolkit provided with MATLAB 2021a was used to develop the graphical user interface (GUI) included in PATLAB. Within the PATLAB documentation (including this article) several common terms are used to avoid ambiguity. First, a time series can be thought of as a series of indexed data points in time that, in the context of PAT, represent the set of samples of the acoustic wave that is recorded from an acoustic detector. A detector location is the cartesian location of an acoustic detector in a global coordinate system. The coordinate of any location can be represented by an X and Y set of values for 2D scenarios, and an X, Y, and Z set of values for 3D scenarios. A detector array or a detector configuration represent a set of detector locations and are used interchangeably. The detector array in a PAT system can be translated or rotated as necessary to scan an object. In PATLAB, a scan is referred to as a frame.

Software architecture

PATLAB consists of four different modules designed to work together as a multi-stage pipeline (Fig. 1). The Loading/Simulating Inputs module provides toolboxes for loading time series and detector data from disk and toolboxes to simulate both time series data and detector configurations. The Pre-processing module provides toolboxes for manipulating time series data prior to image reconstruction. The Image Reconstruction toolbox provides access to a wide variety of PAT image reconstruction methods. The Display & Save module provides methods to plot and save time series data, detector configurations, and reconstructed images. The following sub-sections explain the functionality of each module and the toolboxes available within each module.

Fig. 1

PATLAB architectural layout consisting of four modules: loading/simulating inputs, pre-processing, image reconstruction, and display & save.

Loading/simulating Inputs

This module takes detector locations and photoacoustic time series as inputs. Each of these inputs can be loaded from a file or simulated within PATLAB. The load option uses pre-saved data in the MATLAB array format. This module also enables simulation of photoacoustic time series, which can be saved in MATLAB array format. The first toolbox in the module enables users to specify sets of detector locations, i.e., a detector array or a detector configuration. Settings are available to create a computational grid map that encompasses the detector array, as well as settings to configure 2D and 3D scenarios. Fig. 2 provides an overview of the settings available within the Design Detector Locations toolbox.

Fig. 2

An overview of the available settings within the GUI for designing sets of detector locations.

The Simulate Time Series Data toolbox provides an interface to the popular k-Wave MATLAB toolbox [6]. The k-Wave toolbox allows users to simulate and reconstruct photoacoustic wave fields and extract the raw time-series signals. The toolbox allows users to define a computational grid, specify properties of the medium, incorporate detector properties, and define initial pressure distributions for a pre-defined 2D or 3D detector configuration. A flow chart describing the settings provided in the Simulate Time Series Data toolbox can be found in Fig. 3.

Fig. 3

A flow chart showing the functions available in the Simulate Time Series Data toolbox.

Pre-processing

The purpose of the Pre-Processing module is to modify raw time-series signals so that they are prepared for image reconstruction. There are five toolboxes in this module that are useful for processing the raw time-series data (Fig. 4). (i) The data cleaning toolbox provides methods to remove time series related to faulty detectors; trim samples from the beginning and end of time series; zero pad time series; and threshold time-series data. (ii) The normalization toolbox provides methods for normalizing the time-series data, which is useful to correct for offset differences between detectors and variations in laser power over time. Normalization can result in images with higher signal-to-noise ratio (SNR). (iii) The filtering toolbox can be used to apply low pass, high pass, band pass, and band stop filters to time series to enhance or suppress contributions from specific spectral ranges. For example, systematic electronic noise that degrades the final image can be removed prior to reconstruction with a band-pass filter designed to preserve the PAT signal. Another example is where the time series contain frequencies above the Nyquist frequency resulting in aliasing artifacts in the final image. This can be overcome to some degree by low pass filtering the time series prior to image reconstruction. (iv) The down-sampling toolbox offers users the ability to apply two methods of decimation to time series. Down-sampling is advantageous in situations where the acquired time-series signal is sampled at a higher rate compared to the rate needed to achieve a specific image resolution. In standard decimation, every Mth sample is picked, and the intervening samples are discarded (M is an integer value and called the decimation factor). Processing times can be reduced significantly by decimating the time-series data, since the amount of data is reduced. Also, time series that are oversampled well beyond the required target Nyquist rate can be down-sampled by this method prior to follow-up processing. The second decimation method called moving average decimation applies averaging (a low-pass filtering) and decimation at the same time. (v) The final toolbox in the pre-processing module is matched filtering. Matched filtering correlates a template signal with a time series, then convolves the template signal with the local maxima found by correlation. Matched filtering can improve the contrast and SNR of reconstructed images by a factor of 2 at the expense of reduction in spatial resolution [8]. Examples of template signals are the frequency response of a detector, a standard N-shaped photoacoustic signal, or even a segment of time series data. It is worth mentioning that the pre-processing tools can be used sequentially or in any order.

Fig. 4

A flowchart showing the toolboxes and settings available in the Pre-Processing module.

Image reconstruction

Image reconstruction is the primary task of PATLAB. The reconstruction module utilizes user provided settings to configure and run one of the several reconstruction algorithms. Settable parameters relate to the 2D or 3D grid map defining the imaging plane; detector array properties, frames and samples (time-series data) to include in the reconstruction, pre-processing steps; and the reconstruction algorithm. Image reconstruction algorithms fall into four categories: back projection (BP), time reversal (TR), Fourier transform (FT), and model-based. All reconstruction algorithms in PATLAB assume homogeneous media. Each image reconstruction algorithm implemented in PATLAB is described in more detail below.

Back Projection: With BP, an image reconstruction domain is defined as a cartesian space representing the location of each light absorber in relation to the detectors. Then, the algorithm determines the time-of-flight using the location of the acoustic detectors, the speed-of-sound, and the location of the image grid points relative to the detectors. Using the principle of time-of-flight and knowledge of the sampling of data acquisition, each grid point is assigned to a set of samples in each recorded acoustic signal. Finally, the image is reconstructed by adding up the sets of samples from all available time series signals at each grid point.

Derivative algorithms of BP mainly differ in how they handle the acoustic signals. The simplest algorithm is delay-and-sum (DAS), which computes a summation of the time points associated with each correspond grid point. The algorithm assumes that the detectors have an omnidirectional sensitivity and the signal acquisition is an integration over circular curves for 2D and spherical surfaces for 3D images. However, due to the limited view of physical detectors, it is preferable to assign a respective weight to each grid point for each detector based on the directional axis of the detector. The weighting method is effective at reducing reconstruction artifacts such as streaking. The DAS algorithm also incorporates an approximation where the detectors are modelled as point-like acoustic detectors. In reality, detectors have a finite sensing surface and the detector signal is the integration of the received energy over the surface. With this approximation, the spatial resolution will also be impaired. By increasing the number of views during data collection, the DAS algorithm can provide improved image resolution, but may not preserve the high-frequency components of the light absorbers. An approach called filtered BP (FBP) addresses this issue by applying a high-pass filter (e.g. a ramp filter) in the Fourier domain before back projection to the image space [9]. Universal BP (UBP) is the most advanced BP algorithm as it accounts for the phase sensitivity of the PA signals and preserves both low and high frequency components simultaneously [10].

Time Reversal: The time reversal algorithm solves the inversion problem by simulating the acoustic wave propagation from each transducer back into the Cartesian image space [4], [11]. The method is independent of the detector configuration. Also, since the method iteratively computes the propagating wave field, it can consider the heterogeneity of the medium to perform highly accurate image reconstructions; however, it is computationally expensive [6]. For specific geometries such as planar and linear arrays, it is possible to perform the reconstruction in just one step. Knowing the geometry of the detector array allows the algorithm to map the measured signal into a pressure map of absorbers. This reconstruction algorithm has been named as Fourier transform (FT) reconstruction in k-Wave and PATLAB.

Model-based: To increase image reconstruction performance, model-based image reconstruction methods have been introduced. Using these techniques, an initial image of the absorber structure is converted into photoacoustic signals by means of a mathematical model [2]. The converted signals are then compared with the original signals to optimize the original image iteratively. Currently, PATLAB does not have any model-based image reconstruction methods, but could be extend to include them in the future.

Currently, PATLAB has built-in implementations of the DAS, UBP, FBP, TR, and FT algorithms. Future versions may include model-based algorithms. PATLAB enables the use of directional detector sensitivity and view angle to be specified during image reconstruction. PATLAB provides several ways to define detector directivity based on the focal point concept. For example, a global unit vector can be specified to define a common focal point towards which all detectors point. Alternatively, users can load an external file containing vector information about the directional axis of sensitivity corresponding to each detector. In addition to detector directivity, PATLAB contains a noise reduction technique named coherence factor (CF), which can decrease artifacts in the reconstructed image through suppression of the detector sidelobe sensitivity [12]. This technique of noise reduction can be used along with image reconstruction algorithms such as DAS, UBP, and FBP.

Display and save

The Display and Save module provides users with the ability to visualize and save raw data, detector locations, pre-processed data, and reconstructed images. These functions offer users the ability to review data processing at each step and save data for later use. Tools are provided for selecting specific samples, detectors, and frames for visualization and saving. Visualization of reconstructed images can be performed on individual images or image stacks. Controls are provided to adjust the brightness and contrast of images. The module also provides a post-processing tool to reduce the effect of negativity artifacts caused by limited view and limited bandwidth of the detectors [13].

Results

PATLAB has been in use by our research group for several years to reconstruct PA images from different PAT systems. We designed PATLAB for ease of use. The software handles data in a flexible manner allowing users to define the organization of their data files. For example, loaded time-series data is typically presented in a matrix format that has several dimensions, which may include time-series samples, detectors, and frames (in cases where there are scans in the data acquisition). Within PATLAB, the meaning of each dimension is assigned after the matrix has been loaded, eliminating the need for users to pre-format data in advance. Similar design principles have been used to enable users to load information about detector locations. PATLAB can also be used for educational purposes. Having a 100% illustrative nature, it visually guides users through complex data processing steps, which may be useful as a method for teaching data handling and image reconstruction principles related to PAT. Algorithm developers can use PATLAB to improve their current pre-processing and reconstruction methods and to make fast comparisons. For example, developers can insert new functions into the code and link them to the PATLAB easily. Additionally, PATLAB can be used simply as a front-end to k-Wave to reduce the learning curve of this sophisticated MATLAB toolbox.

Below, we provide image reconstruction results for three separate datasets processed using PATLAB. The first two datasets were obtained from simulations using the k-Wave toolbox. These simulations were carried out to demonstrate the reconstruction algorithms included in PATLAB. Data for the third set was acquired from star phantom using a PAT system available in our lab. The PAT system was capable of capturing a near full view dataset. The third dataset was used to demonstrate the different pre-processing steps available in PATLAB.

Image reconstruction with simulated data

With the Design Detector Locations toolbox, we generated two detector configurations in PATLAB. In the first configuration, 512 detectors were arranged in a circle with a radius of 20 mm, while in the second design, 512 detectors were arranged in a line with a length of 40 mm. We then used the Simulate Time Series Data toolbox to perform PAT simulations using each detector configuration. A modified Shepp-Logan phantom 20 mm in width was used as the initial pressure distribution. Also, a computational grid with the size of 41 × 41 mm2 was defined to include both the phantom and detector array. The simulations were conducted using a sampling rate of 50 MHz. Given the sampling rate and a speed of sound of 1500 m/s, the grid pitch was 0.03 mm in both the X and Y directions resulting in 1367 * 1367 grid points. Each detector was configured with an omnidirectional sensitivity profile.

Different image reconstruction algorithms provided by PATLAB were used to reconstruct the simulated time series data. A note must be made that the FT algorithm could only be used for linear and planar configurations and was therefore not used for the ring detector configuration. However, with the FT algorithm there was no need to load the detector locations, only some prior information about the array size and the detector pitch was required. In total, five algorithms were used to reconstruct the linear array data and four algorithms were used to reconstruct the ring array data as shown in Fig. 5. Based on the comparison of the output from the different algorithms, it was evident that each provides different results for a given dataset. Some algorithms provide higher contrast compared to others, but perhaps at the expense of a global property such as image resolution. With PATLAB, users can easily study the performance differences between the available reconstruction algorithms for their specific data and system configuration. Different reconstruction algorithms provide different dynamic ranges, which are influenced by pre- and post-processing. For example, matched-filtering and normalization in pre-processing and coherence factor and directivity weighting in post-processing can affect the dynamic range. For this reason, the images (ground truth and reconstructed) in Fig. 5 have been normalized between 0 and 1 to compare the quality of the results obtained by different reconstruction algorithms.

Fig. 5

Comparison of five PAT reconstruction algorithms for two datasets obtained with Simulate Time Series Data toolbox provided in PATLAB. (a) Ring array consisting of 512 detectors (dashed red line), Shepp-Logan phantom used as the initial pressure distribution, and two trajectories (yellow arrows) used to plot line profiles in horizontal and vertical directions. (b-e) DAS, UBP, FBP, and TR reconstruction algorithms used to reconstruct time series data simulated from (a), respectively. (f) Linear array consisting of 512 detectors (dashed red line) and Shepp-Logan phantom. (g-k) images reconstructed using the DAS, UBP, FBP, TR, and FT algorithms using time series data simulated from (f), respectively. (l & m) Line profiles extracted from (a-e) in the vertical and horizontal directions, respectively. (n & o) Line profiles extracted from (f-k) in the vertical and horizontal directions, respectively. In (l-o), the vertical axis represents the amplitude (a.u.) and the horizontal axis represents pixel position along each trajectory.

3D reconstruction using experimental data

Next, we used PATLAB to process and reconstruct experimental data collected with a PAT system available in our laboratory [14]. The PAT system provided near full view 3D PAT imaging of objects. The system incorporated an array consisting of 64 detectors arranged in two parallel rings (Fig. 6a). Each ring contained two types of detectors evenly distributed but alternating in order. The first type of detector had a central frequency of 0.4 MHz. The second detector type had a central frequency of 0.9 MHz. The detector array was attached to a robot arm with six degrees of freedom, which enabled precise rotation and translation of the array. For the experiment, the following scanning procedure was used. (i) With the axial center of the array at an XY grid position, the array was incremented in rotation 37 steps about its Z-axis (i.e. axis of symmetry in Z-direction for the array), where each step was 4.5°. (ii) The array was incremented in a single rotational step of 10° about its axial axis (Fig. 6a). (iii) like (i), the array was incremented 37 steps in rotation about its Z-axis, where each step was 4.5°. (iv) The array was homed to its starting position in both rotation about its Z-axis and its axial axis. As shown in Fig. 6b, the rotation steps (i) through (iii) resulted in nearly a complete sampling of detector locations on a sphere. (v) The detector array was incremented in XY to a new grid position. (vi) Steps (i)-(v) were repeated to sample a 7 by 7 grid in the XY plane with a step size of 10 mm in both the X- and Y-directions (Fig. 6c). The resulting dense scan minimized signal loss in space while taking less than 20 min to complete (Fig. 6d). A phantom made of 2% w/v agarose and 2% v/v Intralipid™ containing absorber sections made of 0.015% v/v India ink was used as the imaging target (Fig. 7). The phantom was approximately 125 mm in diameter and 11 mm in height and included alternating wedges of absorbing and non-absorbing material. Spheres 10 mm in diameter were embedded in each wedge such that the sphere was made from the same material as the wedge, but with a complimentary absorption property.

Fig. 6

The detector configuration used to collect the experimental PAT data. (a) The detector array consisted of two rings of 32 elements each. (b) A 330-degree sphere of detectors formed by 74 rotational positions of the array shown in (a). (c) Grid points in XY, where each represents the center of rotation for the scan described in (b). (d) The resultant dense scan showing the combined detector locations for 49 scans described in (c).

Fig. 7

Phantom used as the imaging target for the experimental dataset. The phantom consisted of agarose gel doped with Intralipid™ and/or India ink to simulate scattering and absorption of light in tissue, respectively.

To highlight the versatility of the pre-processing functions built into PATLAB, the time series data collected from the star phantom (Fig. 7) was reconstructed with UBP after applying different pre-processing algorithms (Fig. 8). Each row of Fig. 8 shows the result for each type of pre-processing algorithm. Within each row are plotted examples of pre-processed data (a), reconstructed images (b-f), and line profiles extracted from the image data (g). Within each row, the five images (b-f) represent 2D image slices obtained from the 3D image, where the image slices are 3 mm apart and ordered from top (b) to bottom (f). The image size was 130 mm by 130 mm for each slice and the grid spacing was 0.5 mm in the X and Y directions (260 ×260 pixels per slice). As shown in Fig. 8a (row 1), the raw time-series contained electronic noise. Notably, there was a large signal recorded at the beginning of the time series. The large signal was caused by photoacoustic loading of the detector by the pulsed laser. Detector loading is often considered an artifact and does not convey valuable information about the target. We used the PATLAB data cleaning tool to remove the loading signal and threshold the time series to reduce electronic noise. The results show that the artifact introduced by detector loading and the electronic noise before and after the PA signals could be artificially removed (Fig. 8a (row 2)). The resultant images had higher contrast which was apparent by the relative feature sizes shown in the line profiles.

Fig. 8

Test of different pre-processing methods provided in PATLAB for an experimentally obtained dataset. All images were reconstructed using the UBP algorithm. Each row (i-vi) shows results from the same time series dataset, except the dataset was pre-processed prior to image reconstruction (pre-processing algorithm indicated at the beginning of each row). (a) Example of time series data in which the vertical axes represent amplitude, and the horizontal axes representing sample points. b-f) Five different slices of the imaged volume spaced 3 mm apart. Each image is 130 mm by 130 mm and comprised of 260 × 260 pixels. g) Line profiles for the indicated yellow arrows in column (c). The blue dashed lines are the expected locations of the edges of features in the imaged target. For this column all vertical axes indicate amplitude (A.U.), and horizontal axes indicate pixels.

For PAT systems, detector sensitivity is an important for obtaining high quality images. Many PAT systems contain an array of detectors, where the sensitivity of each detector may vary with respect to its neighbours. The differences in detector sensitivity can affect the final reconstructed image since some regions of an object may be sampled with higher sensitivity compared to other regions. One method to mitigate this problem is to normalize the time series data before reconstruction. The third row in Fig. 8 shows the results for a standard z-score normalization as a pre-processing step. The reconstructed image results demonstrate improved contrast compared to the images reconstructed from the raw data.

Filtering of PA time series to reduce electronic noise is often used prior to image reconstruction. Filters can be designed to select and remove a specific frequency range. In the example shown in row (iv) of Fig. 8., we applied a Butterworth bandpass IIR filter with cut off frequencies of 0.2 MHz and 1 MHz in order of 100. In the filtered time series data there was a noticeable reduction in high frequency noise, but possibly at the expense of high frequency features. Row (v) of Fig. 8, shows an example of matched filtering, where we used the frequency response of the detectors as the template signal. Matched filtering resulted in time series data with higher signal-to-noise compared to the raw time series data. The related reconstructed images and the line profile had sharp contrast changes at the edges for the top layer of the target. However, intermediate image slices were less defined than comparable images processed with other filter methods.

The final example of pre-processing relates to down-sampling of time series data (Fig. 8 row (vi)). We applied a moving average with a decimation factor of 5 and an offset of 0. The algorithm resulted in a smaller time series data set, but produced image reconstruction results almost identical to those obtained from the raw data.

We demonstrate additional features of PATLAB in Fig. 9 using the same raw dataset used for creating Fig. 8. As mentioned earlier, PATLAB allows users to select among loaded detectors, frames, and samples for the reconstruction algorithms. In some cases, these features can be used to remove defective data or select a subset of detectors to prevent oversampling. Within the dataset, we had a total of 64 × 2 × 37 × 49 = 232,064 effective detectors. While this provided a dense, near-full view 3D scan of the star phantom, the density of the scan was not uniform about the phantom. In particular, regions near the top and bottom of the array featured greater scanning density. In principle, these oversampled regions would contribute more greatly to the final image compared to regions were there was less dense scanning (e.g. sides). To partially compensate for this oversampling effect, we omitted eight of the detectors at the top and bottom poles of the ring array and reconstructed the PAT images (first row of Fig. 9). The resultant reconstructed images had less background signal near the centre of the field of view compared to the images reconstructed from the raw dataset. In the second row of Fig. 9, we present the effect of applying the coherence factor gaining method to the dataset prior to image reconstruction. The method resulted in images with better background flatness and higher contrast features; however, images representing interior slices of the phantom had lower signal intensity. In the third row of Fig. 9, we included a 10° directivity pattern for each detector, where the directivity was symmetric about the normal vector of each detector. The reconstructed images appeared to have features that were more uniform in appearance and higher in contrast compared to images reconstructed from datasets pre-processed in other ways (see Fig. 8).

Fig. 9

Demonstration of PATLAB image reconstruction with the ability to select a limited number of detectors (first row), to include adaptive beamforming with coherence factor (second row), and to define and use detector directivity (third row). (a-e) Five different image slices from a 3D reconstruction of the star phantom (Fig. 7). Each slice is separated by 3 mm. f) Line profiles extracted from each image in (b) for the indicated yellow trajectory. The blue dashed lines indicate the locations of the edges in the imaged target. The vertical axes indicate the amplitude (a.u.) and the horizontal axes indicate the pixels along the trajectory.

For the results presented in Fig. 8 and Fig. 9 UBP was used, but other algorithms such as DAS, FBP, and TR can also be used. In general, UBP, DAS, and FBP use similar computational methods with differences in how they handle raw signals, but TR uses a completely different procedure. The TR algorithm propagates acoustic waves from each detector into the field of view until the signals converge to the shape of the initial wave generator. This procedure needs high computational resources because it requires a grid that includes all detectors and the desired field of view. In the case of our experimental PAT system presented in Fig. 6d, a reconstruction volume of 400 * 400 * 300 mm3 was needed to perform a complete TR reconstruction considering all available detectors. This was not possible with the computer resources we had on hand. In PATLAB, it is possible to decrease the computational grid size by considering only a portion of the detectors, for example those in a limited z range. In Fig. 10, an example of the TR algorithm along with UBP for detectors with z ranging from − 5 to + 5 mm is shown. The reconstructed images span the coordinates of − 65 mm < x,y < 65 mm, and z = 0. There is no doubt that the TR algorithm can outperform the UBP, but at the cost of more resources and time. As an indicator, TR for this example took almost 8 h, which is considerably longer than UBP, which took approximately 1 min.

Fig. 10

A demonstration of TR reconstruction (a) and UBP reconstruction (b) using a limited number of detectors (c) for our experimental system (Fig. 6) and the star phantom (Fig. 7).

Discussion

PATLAB is a platform for facilitating photoacoustic data simulation and image reconstruction, accelerating workflows, and serving both novice and expert researchers. This software has the potential to expedite research in PAT because it organizes numerous algorithms and tools for each stage of data processing and image reconstruction. The graphical interface makes it useful for all PAT researchers, regardless of expertise in coding and algorithms. PATLAB is currently the only open-access program that makes PAT simulation and reconstruction accessible with a graphical user interface.

Even with all the advantages that come with PATLAB, it is not without flaws. One of the main downsides of the software is its dependence on CPU-based algorithms. Optimization and transfer of these algorithms to GPU processing would result in substantial reductions in computation time. In addition, PATLAB was designed for offline computations and currently does not support real-time imaging systems.

The limitations described above could be addressed by taking advantage of the built-in extensibility of PATLAB. Algorithms designed for GPU processing could be added as modules in future updates. The dramatic reductions in computation time provided by GPU processing could be the basis for developing real-time image data loading, data pre-processing, and image reconstruction modules. Furthermore, PATLAB could become more feature complete if modules for model-based and AI-enabled image reconstructions were included in future updates.

Conclusion

PATLAB is a new computational software package for PAT image reconstruction with a graphical user interface. The software has been designed to be useful for researchers in the field of PAT imaging with the mission of accelerating research in this field. A comprehensive user manual is provided with the software to guide the users through each processing step. Several pre-processing methods are available to improve the quality of time series data before reconstruction and several reconstruction algorithms are available for obtaining comparative results. Working directly with PATLAB, we conducted PAT simulations using the included Simulate Time Series Data toolbox, which has been developed to extend the k-Wave toolbox. We then demonstrated several image reconstruction techniques using the simulated datasets. An experimentally obtained dataset was also used to illustrate the potential of PATLAB for pre-processing noisy time series data and reconstructing 3D PAT images of a large and extended target.

Declaration of Competing Interest

The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:Jeffrey J.L. Carson reports financial support was provided by Natural Sciences and Engineering Research Council of Canada. Jeffrey J.L. Carson reports financial support was provided by Canadian Institutes of Health Research.

Metadata

Nr	Code metadata description
C1	Current code version	v1.0
C2	Permanent link to code/repository used for this code version	https://github.com/Lawson-Optics-Lab/PATLAB
C3	Permanent link to reproducible capsule	NA
C4	Legal code license	LGPL-3.0
C5	Code versioning system used	git
C6	Software code languages, tools and services used	MATLAB
C7	Compilation requirements, operating environments and dependencies	MATLAB R2021b or later
C8	If available, link to developer documentation/manual	https://github.com/Lawson-Optics-Lab/PATLAB/tree/main/docs
C9	Support email for questions	jcarson@lawsonimaging.ca