Photoacoustic tomography (PAT), also referred to as opto-acoustic tomography, is a cross-sectional or three-dimensional (3D) biomedical imaging technique. The physical foundation for PAT is the photoacoustic effect discovered in 1880 by Alexander Graham Bell^[1], who found that modulated light can excite sound waves in materials. Modern PAT typically uses a short-pulsed laser to illuminate biological tissues. The laser energy is absorbed by tissue chromophores and converted into heat. Under the conditions of stress confinement and thermal confinement, acoustic signals are generated due to the thermoelastic effect and are recorded by ultrasound detectors to recover the optical absorption property of tissues^[2–5]. By combining optical excitation and acoustic detection, PAT has the advantages of rich optical contrast and high ultrasonic resolution in deep tissues and has found unique applications in a range of biomedical fields^[6–9].

PAT has three major embodiments: photoacoustic computed tomography (PACT), photoacoustic microscopy (PAM), and photoacoustic endoscopy (PAE). Among them, PACT uses diffused light to illuminate biological tissues and can achieve hundred-micron-level imaging resolution at centimeter-level tissue depth^[10–13]. Analogous to X-ray computed tomography, high-quality image formation in PACT relies on sophisticated image reconstruction procedures, without which the details in an image may not be able to be properly resolved. In 1981, Bowen studied thermoacoustic imaging of soft tissues using ionizing radiation (high-energy electrons, X-ray photons, neutrons, and other charged particles) and non-ionizing radiation (radio waves, microwaves, and ultrasonic waves)^[14–16] and presented one-dimensional (1D) depth-resolved thermoacoustic signals from a human upper arm using radio-wave excitation. However, the results are only 1D depth-resolved and have no lateral resolution. In the mid-1990s, Oraevsky^[17–19] and Kruger^[20,21] independently studied laser-induced photoacoustic imaging of biological tissues and presented experimental 1D signals and two-dimensional (2D) scan images (see Fig. 1). The results are depth-resolved but poorly laterally resolved due to the lack of image reconstruction procedures. To obtain cross-sectional or 3D images resolved both axially and laterally, Kruger et al. in 1995 developed an approximate filtered back projection (FBP) image reconstruction algorithm for laser-induced PACT analogous to the FBP algorithm used in X-ray computed tomography^[22]. The results suggest that the approximate FBP algorithm helps resolve imaging targets both in the axial and lateral directions (see Fig. 1).

High-performance image reconstruction is critical to image formation in modern PACT imaging systems. Ever since the publication of the approximate FBP algorithm by Kruger and colleagues, several novel image reconstruction approaches have been proposed to reconstruct high-quality PACT images. According to whether the emerging deep learning technique is involved, image reconstruction approaches in PACT can be divided into two main categories: conventional image reconstruction methods and deep learning (DL)-based image reconstruction methods. Conventional image reconstruction methods do not involve deep learning and mainly include five classes of algorithms: FBP, delay and sum (DAS), series expansion (SE), time reversal (TR), and iterative reconstruction (IR).

FBP is a class of analytical image reconstruction algorithms originating from X-ray computed tomography that project filtered photoacoustic signals back to target regions according to the time of flight (TOF) of photoacoustic signals. Kruger adopted this idea and proposed an approximate version of FBP for image reconstruction in PACT in 1995, as mentioned previously^[22]. However, this algorithm can only achieve approximate image reconstruction. To reduce reconstruction error, Finch et al. developed an exact FBP algorithm for image reconstruction in a spherical detection geometry in 2004^[23], and Xu and Wang developed a universal FBP algorithm for exact image reconstruction in infinite planar, cylindrical, and spherical detection geometries in 2005^[24]. FBP algorithms are simple yet efficient and have been widely used for image reconstruction in PACT.

DAS-based reconstruction is essentially a simplified version of the FBP algorithm, which directly projects measured photoacoustic signals back to target regions without filtering. It is one of the most widely used beamforming methods in ultrasound imaging and can potentially be adapted for image reconstruction in PACT. In 1998, Hoelen et al. first used a simple weighted DAS algorithm for image reconstruction in PACT and achieved good reconstruction results^[25,26]. Several modified algorithms were subsequently developed to improve the reconstruction performance of DAS^[27–30]. DAS-based reconstruction algorithms are typically fast and easy to implement. However, they are intuitive and empirical image reconstruction techniques and can only achieve approximate image reconstruction.

SE is a class of image reconstruction algorithms that use mathematical series to represent the image to be reconstructed. In 2001, Köstli et al. derived an exact SE-based inverse formula in the Fourier domain for image reconstruction in a planar detection geometry^[31]. Kunyansky took a step further and derived SE reconstruction formulas for image reconstruction in circular, spherical, cylindrical, and cubic geometries^[32,33]. Compared with other image reconstruction algorithms in PACT, SE algorithms may be computationally efficient for certain detection geometries (e.g., planar geometry) due to the application of the fast Fourier transform (FFT) algorithm during computation^[31].

TR algorithms recover photoacoustic images by running a forward numerical acoustic propagation model backward and re-transmitting the photoacoustic signals measured by each detector in a temporally reversed order. In 2004, Xu and Wang presented an analytical TR model for image reconstruction in arbitrary closed detection geometries^[34]. Burgholzer et al. developed a numerical TR algorithm for image reconstruction in arbitrary detection geometries based on the finite difference method^[35]. To improve computational efficiency and image quality, Treeby et al. further developed a numerical k-space pseudospectral TR algorithm for image reconstruction in heterogeneous media^[34]. The TR algorithms can couple the acoustic properties of media (e.g., SOS, density, dispersion, and absorption) and can be used for image reconstruction in arbitrary closed detection geometries. Therefore, they are regarded as the “least restrictive” image reconstruction algorithms in PACT.

The last class of conventional image reconstruction methods is IR. It solves the image reconstruction problem by iteratively seeking the optimal solution that minimizes the error between measured projection data and the estimate from constructed mathematical models. In 2002, Paltauf et al. proposed the first IR method to improve image reconstruction quality under non-ideal imaging conditions and achieved 3D image reconstruction with reduced artifacts^[36]. After that, much research was conducted by different groups to improve the performance of IR, such as improving the computational accuracy of the system matrix in the model^[37,38], compensating for the response of detectors^[39,40], coupling the acoustic property of the media^[41], and accelerating the reconstruction process^[42,43]. Compared with other algorithms, IR algorithms are typically slow but can yield high-quality imaging results when only limited projection data are available. Therefore, they are more suitable for non-ideal imaging scenarios, such as limited-view imaging, sparse-view imaging, and imaging in acoustically heterogeneous media.

The second major category of image reconstruction approaches in PACT is fast-developing DL-based methods, which are inspired by the successful use of DL in a range of fields. DL-based methods typically train neural networks and use them to automatically transform input data into output photoacoustic images. Compared with conventional image reconstruction approaches, DL methods are more efficient and can handle more complicated scenarios. In 2018, Antholzer et al. first proposed a deep convolutional neural network (CNN)-based method for PACT image reconstruction under sparse-view sampling and opened new opportunities for intelligent image reconstruction in PACT^[44]. After that, a series of DL-based methods was developed for image reconstruction using both simulated and experimental datasets^[45–47]. State-of-the-art DL methods are powerful enough to achieve preprocessing in the data domain, postprocessing in the image domain, hybrid processing in both the data and image domains, learned IR, and direct image reconstruction from the data domain to the image domain^[48–50]. They have been used to address a range of image reconstruction problems in PACT, such as detector bandwidth expansion^[51,52], resolution enhancement^[53,54], artifact removal^[55], ultralow-laser-energy imaging^[56], reconstruction acceleration^[57], and reconstruction enhancement in sparse-view and limited-view imaging^{[44,46,47,58,59]}.

Over the past three decades, great achievements have been made in PACT image reconstruction. To date, there have been a few excellent reviews devoted to summarizing the research progress in this field^[60–63]. However, most review papers in the literature are from more than 10 years ago and cannot reflect the latest research advances^[60–62]. Moreover, some review papers cover only particular subjects of the field, such as IR reconstruction^[63,64] or DL-based reconstruction^{[48–50,65–68]}. There is a pressing need to prepare new work to systematically review the recent achievements in this field. For the above reasons, we prepared this review, which differs from other works in the following aspects. First, the review was prepared from a historical perspective. We surveyed the development of each type of reconstruction algorithm and put them into the context of the entire history of PACT. Second, the review is comprehensive. It not only contains the five classes of conventional reconstruction algorithms (FBP, DAS, SE, TR, and IR) but also covers the state-of-the-art DL-based reconstruction algorithms. Third, the review contains comparative studies. A dedicated part (Sec. 5) was prepared to compare the performance of each type of image reconstruction algorithm in terms of image reconstruction quality, reconstruction speed, and memory footprint. Finally, the review is expected to be beginner-friendly. The entire review contains 51 figures, 14 tables, 126 mathematical equations, and comparative studies of each algorithm and is easy for novices to understand.

The remainder of this review is organized as follows. Section 2 describes the forward problem of PACT, including photoacoustic signal generation, propagation, detection, and the mathematical foundation (i.e., the Radon transform) for image reconstruction. Section 3 reviews the basic principles of the five classes of conventional image reconstruction algorithms, namely, FBP, DAS, SE, TR, and IR. Section 4 surveys the DL-based image reconstruction methods from the perspective of preprocessing, postprocessing, hybrid processing, learned IR, and direct image reconstruction. Section 5 provides comparative studies of major image reconstruction algorithms in PACT. Section 6 highlights the major challenges and future directions for PACT image reconstruction. Finally, Sec. 7 offers concluding remarks. Figure 2, Table 1, and Table 2 list major topics, abbreviations, and symbols used in this review, respectively.

It is worth pointing out that image reconstruction in PACT involves two aspects: acoustic inversion for determining the distribution of initial acoustic pressure and optical inversion for determining the map of optical absorption^[61,62,69]. This review focuses on the acoustic inversion problem, which is important for image formation in PACT. Moreover, the mathematical foundation for image reconstruction in PACT is equivalent to that in thermoacoustic tomography^[70–73]. Therefore, the image reconstruction algorithms discussed in this review can in principle be used for thermoacoustic tomography.

To achieve efficient photoacoustic signal excitation and generation, two conditions, i.e., thermal confinement and stress confinement, should be satisfied. The two conditions are related to two timescales, i.e., the thermal relaxation time τth and the stress relaxation time τs. The thermal relaxation time τth characterizes the diffusion of heat from the region heated by laser pulses, while the stress relaxation time τs describes the propagation of acoustic waves from the heated region. τth is given by^[74]τth=dc2αth,(1)where dc is the characteristic dimension of the heated region and αth is the thermal diffusivity (m2/s). τs is given by^[74]τs=dcv0,(2)where v0 is the speed of sound (SOS). If the laser pulse duration is less than the thermal relaxation time, i.e., τ<τth, the thermal diffusion is negligible during laser heating, in which case the thermal confinement condition can be considered to be satisfied. Similarly, if the laser pulse duration is less than the stress relaxation time, i.e., τ<τs, the pressure propagation can be ignored during pulse heating, in which case the stress confinement condition can be considered to be satisfied. In PACT, the thermal and stress confinement conditions should be satisfied simultaneously. For instance, assuming that in soft tissues, the SOS is 1500 m/s, the thermal diffusivity is 0.14mm2/s, and the dimension of the heated region is 50 µm, the thermal relaxation time τth is calculated to be 18 ms, and the stress relaxation time τs is 33 ns. This indicates that the laser used for imaging should have a pulse duration of less than 33 ns to achieve efficient photoacoustic signal excitation and generation.

In PACT, the local fractional volume expansion ΔV/V induced by the photoacoustic effect can be described as^[74]ΔVV=βΔT−κΔp,(3)where ΔT and Δp represent the changes in temperature and pressure, respectively, β is the thermal coefficient of volume expansion, and κ is the isothermal compressibility, which can be written as κ=Cpρv02Cv,(4)where ρ is the mass density and Cp and Cv are the specific heat capacities at constant pressure and volume, respectively.

Under the conditions of thermal and stress confinement, the fractional volume expansion in Eq. (3) is negligible (i.e., ΔV/V=0). The initial photoacoustic pressure p0=Δp and can be written as p0=βΔTκ,(5)where the temperature rise can be calculated by^[75]ΔT=ηthHρCv.(6)

Here ηth is the percentage of the laser energy converted into heat, and H is the heat energy deposited per unit volume, which is defined as the product of the absorption coefficient μa and the optical fluence F (i.e., H=μaF). The initial acoustic pressure p0 can thus be rewritten as p0=(βv02Cp)ηthH=ΓηthH,(7)where Γ is the Grüneisen parameter, a dimensionless constant representing the efficiency of the conversion of heat to pressure. For water at body temperature (37°C), Γ≈0.2, which indicates that 20% of the thermal energy deposited by a laser in water couples into acoustic energy.

The generation of photoacoustic signals can be illustrated using a numerical example. Assume that the laser employed for imaging has a fluence of 10mJ/cm2 (F=10mJ/cm2), which is within the American National Standards Institute (ANSI) safety limit^[76], and the biological tissue being imaged has the following physical parameters: Γ=0.2, μa=0.1cm−1, ρ=1.0g/cm3, and Cv=4.0Jg−1K−1. The factor ηth can be set to 1 because the fluorescence and nonthermal absorption of biological tissues are typically weak. In this way, the temperature rise ΔT is calculated to be 0.25 mK and the initial acoustic pressure p0 is 200 Pa, which indicates that a 1 mK temperature rise can produce an acoustic pressure rise of 800 Pa.

The generation and propagation of photoacoustic waves can be mathematically modeled by the photoacoustic wave equation^[77]. Generally, three first-order equations, including the linearized equation of motion, the linearized equation of continuity, and the thermal elastic equation, can be used to model the properties of acoustic wave generation and propagation^[78]. If the medium is lossless and the thermal conductivity can be ignored, the three equations can be written as^[79]∂∂tu(r,t)=−1ρ(r)∇p(r,t),(8)∇·u(r,t)=−1ρ(r)v02(r)∂∂tp(r,t)+β∂∂tΔT(r,t),(9)ρ(r)Cp∂∂tT(r,t)=H(r,t),(10)where u(r,t) is the particle velocity, v0(r) is the SOS, ρ(r) is the mass density, p(r,t) is the acoustic pressure at position r and time t, T(r,t) and ΔT(r,t) are the temperature and temperature rise, respectively, and H(r,t) is the heating function that represents the thermal energy deposited per unit volume and per unit time. The second-order photoacoustic wave equation can be obtained by eliminating the variable u(r,t) from the three equations as^[78,79]{ρ(r)∇·[1ρ(r)∇]−1v02(r)∂2∂t2}p(r,t)=−βCp∂H(r,t)∂t.(11)

Equation (11) describes the relationship between the acoustic pressure p(r,t) and the photoacoustic source associated with the heating function H(r,t). The source term [right side of Eq. (11)] is proportional to the first derivative of the heating function H(r,t), which indicates that the heating function H(r,t) should be time-varying to achieve efficient photoacoustic signal generation.

When the medium is acoustically homogeneous, the photoacoustic wave equation in Eq. (11) reduces to^[78,79](∇2−1v02∂2∂t2)p(r,t)=−βCp∂H(r,t)∂t.(12)

The forward solution for the wave equation in a homogeneous medium in Eq. (12) can be obtained using Green’s function^[75,80] as p(rd,t)=βCp∫−∞+∞∫VG(rd,t;rs,ts)∂H(rs,ts)∂tsdrsdts,(13)where p(rd,t) is the acoustic pressure detected at position rd and time t and rs and ts represent the location and time of the photoacoustic source, respectively. The Green’s function G(rd,t;rs,ts) in an infinite 3D space has the following form: G(rd,t;rs,ts)=δ(t−ts−|rd−rs|/v0)4π|rd−rs|,(14)where δ is the Dirac delta function.

Under the condition of stress confinement, the heating function can be decomposed as H(rs,ts)=H(rs)δ(ts), where H(rs) is the heat deposited per unit volume. In this way, Eq. (13) becomes p(rd,t)=βCp∫−∞+∞∫VG(rd,t;rs,ts)H(rs)δ′(ts)drsdts,(15)where δ′ is the derivative of the Dirac delta function. Using the property ∫δ′(t−t0)f(t)dt=−f′(t0), Eq. (15) becomes p(rd,t)=βCp∫VH(rs)∂G(rd,t;rs,ts)∂tdrs|ts=0.(16)

Applying the relation in Eq. (7) (ηth is set to 1), p0(rs)=ΓH(rs)=(βv02Cp)H(rs),(17)and substituting the Green’s function [Eq. (14)] into Eq. (16), we obtain the general forward solution of the photoacoustic wave equation, i.e., p(rd,t)=14πv02∂∂t∫Vp0(rs)|rd−rs|δ(t−|rd−rs|v0)drs.(18)

The measured pressure p(rd,t) is associated with the velocity potential ϕ(rd,t) via^[81]p(rd,t)=−ρ∂ϕ(rd,t)∂t.(19)

The velocity potential can be written as ϕ(rd,t)=−14πv02ρ∫Vp0(rs)|rd−rs|δ(t−|rd−rs|v0)drs.(20)

Equations (18) and (20) indicate that the measured pressure p(rd,t) and velocity potential ϕ(rd,t) are linearly related to the initial acoustic pressure p0(rs) and are inversely proportional to the distance between the detector and the source, i.e., |rdrs|.

The concept of acoustic pressure and velocity potential can be illustrated through a simple example^[81]. Suppose that the photoacoustic source is a uniform sphere and that the detector is a point detector. We examine the velocity potential and acoustic pressure received by the point detector in two situations. In the first situation, the point detector is located at the center of the photoacoustic source, as shown in Fig. 3(a). The velocity potential can be obtained as [Eq. (20)] ϕ(rd,t)={−p0(rs)v0ρt,ifv0t≤ra,0,ifv0t>ra,(21)where ra is the radius of the spherical source. In the second situation, the point detector is located outside the photoacoustic source, as shown in Fig. 3(d). The velocity potential can be calculated as ϕ(rd,t)={−p0(rs)4v0ρd[(d−v0t)2−ra2],ifd−ra≤v0t≤d+ra,0,otherwise,(22)where d is the distance between the point detector and the center of the photoacoustic source. The acoustic pressure p(rd,t) in these two cases can be accordingly calculated using the relationship in Eq. (19).

To quantify the velocity potential and acoustic pressure measured by the point detector, we make the following assumptions: the radius of the spherical source is 4 mm (ra=4mm), the mass density of the medium is 1000kg/m3 (ρ=1000kg/m3), the SOS is 1480 m/s (v0=1480m/s), the Grüneisen parameter at room temperature is 0.12 (Γ=0.12), the optical absorption coefficient μa is 100cm−1 (μa=100cm−1), and the optical fluence is 10mJ/cm2 (F=10mJ/cm2). Therefore, the heat energy H deposited at time zero is 106J/m3, and the initial acoustic pressure p0(rs) [Eq. (7)] is 1.2×105Pa.

Figures 3(b) and 3(c) show the velocity potential and acoustic pressure measured by the point detector in the first case (the point detector is located at the center of the source). The results show that the negative velocity potential −ϕ(rd,t) measured by the point detector first linearly increases as the shell of the integration [Eq. (20)] increases with time and then suddenly drops to zero due to no energy deposition outside the source. Consequently, the acoustic pressure p(rd,t) (i.e., the first derivative of the velocity potential) is a non-zero constant at the beginning and becomes zero with a negative impulse caused by the sudden drop of velocity potential to zero. Figures 3(e) and 3(f) show the velocity potential and acoustic pressure measured by the point detector in the second case (the point detector is located 10 mm away from the center of the source). The results show that the negative velocity potential −ϕ(rd,t) measured by the point detector first increases as the incremental hemispherical shell of the integration [Eq. (20)] advances to the center of the source, and then decreases with time as the shell advances to the rear of the source. Consequently, the acoustic pressure p(rd,t) is initially positive, passes through zero, then becomes negative, and has an N-shaped waveform.

Based on the example in Fig. 3, we can also deduce that the size of a photoacoustic source impacts the characteristics of the photoacoustic signals received by a detector in the time domain. To illustrate this, Fig. 4(a) shows three time-domain photoacoustic signals measured by a point detector 3 mm away from the centers of three uniform spherical sources with diameters of 1 mm, 200 µm, and 50 µm. The time-domain photoacoustic signals all have an N-shaped waveform as expected. However, the amplitude and duration of the photoacoustic signals are distinct. The photoacoustic signals produced by sources with a smaller size have smaller amplitudes and shorter durations, which correspond to higher frequencies and broader bandwidths in the frequency domain, as shown in Fig. 4(b). The center frequency of the photoacoustic signals fc generated by a spherical source can be estimated by fc=v0/3ra, where ra is the radius of the spherical source^[82]. This indicates that detectors with a high center frequency should be employed for high-frequency imaging.

Numerical simulation of the forward propagation of photoacoustic fields helps visualize the sound fields in complex media and solve the acoustic inversion problem in PACT. Numerical simulation can be implemented using the powerful open-source k-Wave toolbox^[83], which was developed by the Photoacoustic Imaging Group at University College London (UCL). In the k-Wave toolbox, the propagation of photoacoustic fields is modeled by three coupled first-order partial differential equations^[84], which are fundamentally the same as the three equations in Eqs. (8)–(10). The propagation model considers the acoustic properties of media, such as acoustic speed, dispersion, and absorption, and can characterize the acoustic propagation problem in heterogeneous media. The k-Wave toolbox solves the propagation model via a k-space pseudospectral method^[83,85,86], which can perform fast and accurate computations with reduced memory. Figure 5 is an example showing the propagation of the photoacoustic fields of a 2D disk and a 3D sphere using the k-Wave toolbox. The results visualize the instant characteristics of the photoacoustic fields during propagation.

The photoacoustic signals propagating outward from a photoacoustic source need to be captured by ultrasound detectors for image formation. Ideally, for photoacoustic signal detection, detector arrays that can perfectly record original signals in time and space should be used. However, perfect detector arrays never exist in reality. The characteristics of a practical detector array, such as detector aperture, detector bandwidth, detector number, and view angle, impact detected photoacoustic signals and final images.

Aperture and bandwidth are two fundamental characteristics of an ultrasound detector and have important impacts on measured photoacoustic signals. An ideal ultrasound detector should have a point-like aperture, in which way it has an omnidirectional response. Moreover, an ideal ultrasound detector should also have an infinite bandwidth so that it can respond to all frequency contents of a signal. Neither of the two conditions, however, is attainable in practice. A practical ultrasound detector always has a finite aperture size and a finite bandwidth. The finite aperture averages photoacoustic signals in space, resulting in a smoothed spatial impulse response (SIR), while the finite bandwidth affects the conversion of photoacoustic signals to electrical signals, leading to a degraded electrical impulse response (EIR). Taking the SIR and EIR of an ultrasound detector into account, the photoacoustic signal measured by a practical ultrasound detector can be mathematically formulated as^[4]p(rd,t)=pideal(rd,t)∗hSIR(rd,rs,t)∗hEIR(t),(23)where pideal(rd,t) is the ideal photoacoustic signal, * represents temporal convolution, and hSIR(rd,rs,t) and hEIR(t) represent the SIR and EIR of a detector, respectively. The non-ideal SIR and EIR of an ultrasound detector distort photoacoustic signals measured in the time domain, as illustrated in Fig. 6, which eventually degrades the image quality of photoacoustic images.

In addition to the aperture and bandwidth of a detector, the detector number and view angle of a detector array also play important roles in photoacoustic signal detection. Ideally, the number of detectors in a detector array used for photoacoustic signal measurement should meet the spatial Nyquist sampling theorem for complete signal acquisition^[4]. The view angle of a detector array should be 4π steradian (full 3D view) to record complete photoacoustic signals in 3D space. However, these two conditions are actually unattainable in practice due to the high fabrication cost of a large number of detectors and the configuration constraints of an imaging system (e.g., a separate space is required for laser illumination). Violating these conditions leads to the problems of image reconstruction from sparse-view and limited-view projections, which are mathematically challenging and will be discussed later in this review. Figure 7 showcases commonly used 2D and 3D detection geometries in PACT^[87,88], which have limited view angles except for the closed spherical array.

The forward signal propagation and detection processes in PACT can be described by the well-known Radon transform^[89], which is the mathematical foundation of computed tomography (CT).

Before discussing the Radon transform in PACT, we first introduce the linear Radon transform in X-ray CT, which is defined as the integral of a function along a straight line. Specifically, the linear Radon transform in X-ray CT is written as^[90]g(s,θ)=∫−∞∞∫−∞∞f(x,y)δ(xcosθ+ysinθ−s)dxdy,(24)where f(x,y) is the original function, g(s,θ) is the sinogram or projection data of the function f(x,y) along the straight line defined in the Dirac delta function δ, and (s,θ) are two parameters of the normal equation of the straight line in the delta function. The inverse Radon transform reverses the forward process and recovers the original function f(x,y) from measured sinograms g(s,θ). A schematic representation of the linear Radon transform is shown in Fig. 8(a).

In PACT and thermoacoustic tomography, the Radon transform integrates a function along a circle (2D) or a sphere (3D) instead of a straight line. The 2D circular Radon transform can be written as^[91]g(s,θ)=∫−∞∞∫−∞∞f(x,y)δ[s−(x−rcosθ)2+(y−rsinθ)2]dxdy,(25)where (r,θ) represents the position of the detector in a polar coordinate system and s is the radius of the integral circle. Similarly, the 3D spherical Radon transform g(s,θ) can be mathematically written as^[92]g(s,θ)=∫−∞∞∫−∞∞∫−∞∞f(x,y,z)δ[s−(x−rsinφcosθ)2+(y−rsinφsinθ)2+(z−rcosφ)2]dxdydz,(26)where (r, φ, θ) represents the position of the detector in a spherical coordinate system and s is the radius of the spherical shell to be integrated. Schematic representations of the circular and spherical Radon transforms are shown in Figs. 8(b) and 8(c).

Using the definition in Eq. (26), the spherical Radon transform in PACT can be obtained from Eq. (18) and has the following form: g(rd,t)=∫p0(rs)δ(t−|rd−rs|v0)drs.(27)

Equation (27) is in the form of the spherical Radon transform [Eq. (26)], representing the integral over a spherical shell with a radius of v0t centered at the detector position rd. The projection data g(rd,t) are related to the measured photoacoustic signal p(rd,t) via g(rd,t)=4πv02|rd−rs|∫0tp(rd,τ)dτ.(28)

The projection data g(rd,t) have a similar definition as the velocity potential in Eq. (19) and can be calculated from the measured photoacoustic signal p(rd,t). The task of image reconstruction in PACT is to find the inverse spherical Radon transform and recover the original function f(x,y) or p0(rs) from the measured photoacoustic signal p(rd,t), which is the focus of the following sections.

Delay and sum (DAS)-based beamforming is a commonly used image reconstruction technique in ultrasound imaging^[93]. Due to similar image formation processes, DAS-based beamforming is also widely used in PACT, where it reconstructs an image by summing the delayed raw photoacoustic signals of each detector. The time delay between each channel is calculated according to the acoustic TOF of the photoacoustic signals from the point of interest to each ultrasound detector. To yield high-quality images, preprocessing of raw photoacoustic signals and/or postprocessing of reconstructed images are usually needed. Here we review the five most commonly used DAS-type image reconstruction algorithms: DAS, delay multiply and sum (DMAS), short-lag spatial coherence (SLSC), minimum variance (MV), and coherence factor (CF), which adopt different preprocessing and/or postprocessing strategies. The basic workflows of these methods are first illustrated in Fig. 9, and detailed descriptions are presented below.

DAS is the most basic beamforming algorithm in ultrasound imaging due to its simplicity, speed, and robustness^[94–96]. Hoelen et al. introduced a DAS method for 3D PACT imaging of blood vessels based on a planar detection geometry in 1998^[26,97]. Feng et al. applied a DAS method in linear-scanning thermoacoustic tomography in 2001^[98]. The DAS method can be mathematically formulated as SDAS(x,z)=∑i=1Msi(t),(29)where SDAS is the reconstructed image, si(t) is the photoacoustic signal measured by the ith detector at time t, M is the total number of detectors, and (x,z) represents the position in a coordinate system. The variable t in si(t) denotes the TOF of the photoacoustic signals from position (x,z) to the ith detector. The workflow and principle of the DAS algorithm are illustrated in Figs. 9(a) and 10, respectively.

The performance of the DAS algorithm was evaluated using a designed numerical phantom with multiple point sources equally distributed along the longitudinal centerline [Fig. 11(a)]. In the evaluation, a linear detector array with a width of 200 mm and 128 elements was placed at the top of the phantom to receive the photoacoustic signals generated from the sources. Figure 11(b) shows the image reconstructed by the DAS method. The results show that DAS can successfully reconstruct the structural information of photoacoustic sources but fails to reproduce the correct amplitude. The reconstructed image has both positive and negative components and is bipolar in magnitude. A unipolar image can be obtained by finding the envelope of the bipolar image, as shown in Fig. 11(c). Log transformation can be further used to improve the contrast of the reconstructed image, as shown in Fig. 11(d). Since DAS treats the delayed photoacoustic signals si(t) from different detectors equally and is non-adaptive, significant side lobes are present in the reconstructed image in Fig. 11(d), which degrades the spatial resolution of the image.

Delay multiply and sum (DMAS) is a variant of the DAS algorithm that can provide improved contrast, signal-to-noise ratio (SNR), and resolution. Similar to DAS, DMAS also sums photoacoustic signals measured by different detectors according to calculated time delays. However, the measured photoacoustic signals in DMAS need to be combinatorically coupled and multiplied before summation. Therefore, DMAS is essentially a nonlinear algorithm.

The DMAS algorithm was first proposed by Lim et al. for confocal microwave detection of breast cancer in 2008 and showed improved identification of embedded malignant tumors in a variety of numerical breast phantoms compared with DAS^[99]. In 2015, Matrone et al. modified this method for B-mode ultrasound imaging and demonstrated that DMAS can provide higher contrast resolution than DAS^[100]. Inspired by the success of DMAS in confocal microwave imaging and ultrasound imaging, several research groups worldwide have conducted in-depth studies of DMAS-based photoacoustic image reconstruction. For example, Alshaya et al. introduced DMAS to the field of PACT in 2016 and proposed a filter DMAS to improve the SNR of reconstructed images^[101]. To further improve the performance of DMAS, in 2018, Mozaffarzadeh et al. developed a double-stage DMAS algorithm that can produce images with improved quality compared with DAS and DMAS but at the expense of greater computational cost^[27]. In the same year, Kirchner et al. proposed a signed DMAS algorithm that can provide linear image reconstruction with increased image quality^[102]. In 2022, Mulani et al. presented a high-order DMAS method, in which multi-term (e.g., three, four, or five terms) multiplication is used to replace two-term multiplication in the original DMAS algorithm^[103].

Generally, the original DMAS algorithm can be mathematically written as^[100]SDMAS(x,z)=∑i=1M−1∑j=i+1Mxij(t),(30)where SDMAS is the reconstructed image, M is the total number of detectors, and xij(t) is given by xij(t)=sign[si(ti)sj(tj)]|si(ti)sj(tj)|,(31)where sign(·) denotes the signum function. The workflow of the DMAS algorithm is illustrated in Fig. 9(b). If the center frequency of the photoacoustic signals si(t) and sj(t) is f0, the multiplication operation in Eq. (31) shifts the center frequency of the resultant signal to 0 and 2f0. As a result, the output signal SDMAS needs to be filtered by a bandpass filter centered at 2f0 to extract the second harmonic component 2f0 while removing the direct current (DC) component. Therefore, DMAS is also called filtered-DMAS. Compared with DAS, DMAS uses both the amplitude and spatial correlation of delayed photoacoustic signals from different detectors for image reconstruction. For this reason, it can partially overcome the drawbacks of DAS and reconstruct images with improved spatial resolution and reduced side lobes^[100].

A downside of DMAS is that the combinational multiplication in the algorithm increases the computational complexity. To solve this problem, in 2019, Jeon et al. proposed a DMAS algorithm with a modified CF, which avoids combinatorial multiplication in DMAS and significantly reduces the total number of multiplication operations^[104]. In 2022, Paul et al. proposed a simplified-delay-multiply-and-standard-deviation (SDMASD) method^[105], which is based on the measurement of the standard deviation of delayed and multiplied signals instead of normal delayed signals. Compared with native DAS and DMAS algorithms, the SDMASD algorithm can achieve real-time imaging using graphics processing units (GPUs) and produce improved image quality.

Short-lag spatial coherence (SLSC) is a beamforming technique that was initially developed by Lediju et al. for ultrasound imaging in 2011^[106]. SLSC reconstructs an ultrasound image by calculating the spatial coherence of measured signals, and the reconstructed image is thus independent of the amplitude of the signals. As a result, SLSC can eliminate adverse effects caused by the different strengths of scatterers in ultrasound imaging. It has been demonstrated that compared with the conventional DAS beamforming algorithm, SLSC can achieve image reconstruction with considerable improvements in terms of resolution, contrast-to-noise ratio (CNR), and SNR.

In PACT, Muyinatu Bell et al. applied the SLSC algorithm to achieve imaging of prostate brachytherapy seeds^[28,107] and demonstrated that the SLSC algorithm can enhance photoacoustic image quality compared with DAS, especially when the intensity of laser illumination is insufficient. Graham and Muyinatu Bell later developed a spatial coherence theory based on the van Cittert–Zernike theorem, a classical theorem in statistical optics, to explore the strengths and limitations of the SLSC algorithm^[108,109]. In 2021, Mora et al. combined SLSC with DMAS and proposed a generalized spatial coherence algorithm for PACT, which can preserve relative signal magnitudes and improve the CNR and SNR of a reconstructed image^[110].

Similar to DMAS, the photoacoustic signals in SLSC are also first delayed according to the TOF from the point of interest to each ultrasound detector and then combinatorically coupled and multiplied before summation. The SLSC algorithm can be formulated as^[106]R^p(l)=1M−l∑i=1M−l∑t=t1t2si(t)si+l(t)∑t=t1t2si2(t)∑t=t1t2si+l2(t),(32)where R^p(l) is the normalized spatial coherence of the signals measured by a detector, M is the total number of detectors, and l represents the number of intervals between two detectors used to calculate the spatial coherence and is called the lag. The final SLSC image is obtained by summing all R^p(l) terms, i.e., SSLSC(x,z)=∑l=1LR^p(l),(33)where L is the total number of lags. A large L helps improve the lateral resolution but decreases the CNR and SNR. Therefore, the value of L needs to be elaborately selected to achieve the best tradeoff between key image quality metrics, such as lateral spatial resolution, CNR, and SNR. The workflow of the SLSC algorithm is illustrated in Fig. 9(c).

Minimum variance (MV) is a weighted DAS method that was originally devised by Capon in 1969 for narrowband signal processing applications such as sonar, radar, and wireless communication^[111]. In 2002, Mann and Walker used a constrained adaptive beamformer (the MV method) for medical ultrasound imaging^[112] and demonstrated its effectiveness in spatial resolution and contrast enhancement. Due to the remarkable improvement in spatial resolution, a large number of MV-based methods have been developed for ultrasound imaging^[113–116].

MV-based image reconstruction has also been studied in PACT in recent years^{[30,117–120]}. The MV reconstruction formula can be written as^[114,116]SMV(x,z)=∑i=1Mwi(t)si(t)=w(t)Hs(t),(34)where M is the total number of detectors, wi(t) is the optimal weight for the photoacoustic signal si(t) measured by the ith detector, s(t) is a vector containing delayed photoacoustic signals from all detectors, w(t) is a vector containing optimal weights for the delayed photoacoustic signals in s(t), and the superscript H denotes the conjugate transpose. The workflow of the MV algorithm is illustrated in Fig. 9(d). The weight w(t) can be adaptively calculated by minimizing the variance of the output SMV while maintaining the unit signal gain at the focal imaging point. The optimization problem can be mathematically written as^[114]minwwHRw,subjecttowHa=1,(35)where R=E[ssH] is the spatial covariance matrix of the delayed photoacoustic signals s(t), E denotes the expectation, and a is the equivalent of the steering vector in narrowband applications^[114,116]. When the photoacoustic signals of each detector are delayed, a is a simple all-one vector. The optimization problem in Eq. (35) can be solved by the method of Lagrange multipliers^[121], which gives w=R−1aaHR−1a.(36)

By substituting Eq. (36) into Eq. (34), MV-based image reconstruction can be achieved.

To improve the robustness of the MV method, the covariance matrix R can be calculated based on a spatial smoothing strategy, where a detector array is divided into a group of overlapping subarrays, and the covariance matrices of all subarrays are calculated and averaged to form the final covariance matrix^[114]. In this way, the covariance matrix R is given as^[114]R=1Nd−L+1∑l=1Nd−L+1sl(t)sl(t)H,(37)where L is the number of detectors in a subarray and l is the index of the detector in the current subarray. To estimate the covariance matrix more accurately and enhance the contrast of MV, temporal averaging of multiple samples can be used together with spatial smoothing^[116]. Once the covariance matrix R is estimated, the optimal weights w can be obtained by Eq. (36). The final image reconstructed by the MV algorithm can be written as^[114]SMV(x,z)=1Nd−L+1∑l=1Nd−L+1w(t)Hsl(t).(38)

Although the MV method can produce narrower main lobes (i.e., the lobe located at the target) and higher spatial resolution than algorithms such as DAS and DMAS, its performance in terms of side lobe (i.e., the lobes adjacent to the main lobe) suppression and image contrast enhancement is limited. Therefore, many studies have been devoted to the development of MV-based hybrid beamforming algorithms. In 2008, Park et al. imposed additional CF weights on MV and achieved enhanced spatial resolution, contrast, and side lobe suppression^[117]. In 2018, Mozaffarzadeh et al. developed an MV-based DMAS method^[118] and an eigenspace-based MV method combined with DMAS for resolution improvement and side lobe reduction^[119]. In 2021, Paul et al. proposed an adaptive-weighting-based MV to address the side lobe issue in MV beamformed images. It was demonstrated that the weighted MV approach can improve SNR while reducing main lobe width and side lobe intensity and has the potential for use in PACT imaging systems with a limited number of ultrasound detectors^[30].

The coherence factor (CF) is a pixel-level adaptive weighting factor that can improve the performance of DAS-based beamforming methods in side lobe suppression and spatial resolution improvement. CF was originally proposed by Mallart and Fink in 1994 for the evaluation of phase aberration correction techniques in scattering media^[122] and was later used as an image quality metric for ultrasound imaging by Hollman et al.^[123]. In 2003, Li et al. presented a generalized CF for ultrasound beamforming in heterogeneous media and showed that the combination of generalized CF and DAS-based beamformers could yield improved image quality^[124].

In PACT, Liao et al. incorporated CF into DAS in 2004 and demonstrated the superiority of the CF-weighted DAS method in improving spatial resolution and SNR compared with DAS^[29]. Wang and Li considered the local SNR in the formulation of CF and developed an SNR-dependent CF for ultrasound and photoacoustic imaging in 2014^[125]. Wang et al. integrated CF with a focal-line-based image formation method to improve the contrast and elevational resolution of 3D photoacoustic imaging in 2016^[126]. Mozaffarzadeh et al. proposed a high-resolution CF weighting technique and achieved improved resolution, SNR, and side lobe suppression in 2019^[127]. Paul et al. considered the noise level variations of raw ultrasound data in the formulation of CF and achieved improvements in image resolution, side lobe reduction, SNR, and contrast in 2021^[128]. In the same year, Mukaddim and Varghese extended CF from the spatial domain to the spatiotemporal domain^[129]. This extension helps cancel out signals with low spatial and temporal coherence and results in higher background noise cancellation while preserving the main features of interest in reconstructed images.

Mathematically, the CF is defined as the ratio of the coherent sum of photoacoustic signals across detectors to the incoherent sum and can be formulated as^[123]CF(x,z)=|∑i=1Msi(t)|2M∑i=1M|si(t)|2,(39)where si(t) is the photoacoustic signal measured by the ith detector at time t and M is the total number of detectors. The value of the CF ranges from zero to one. A large value means that the signals at that point are highly focused and that the point can be reconstructed with high quality. In contrast, a small CF value indicates that the signals are weakly focused and will result in lower image quality. The CF-weighted DAS method can be expressed as SCF-DAS(x,z)=CF(x,z)SDAS(x,z),(40)where SCF-DAS(x,z) is the reconstructed image. The term SDAS(x,z) can be replaced by the outputs of other beamforming methods, such as DMAS, SLSC, and MV. The workflow of the CF algorithm is illustrated in Fig. 9(e).

Each beamforming method mentioned above has advantages and disadvantages. To improve the performance of image reconstruction in PACT, they can be combined to yield hybrid beamforming methods such as DMASD plus DAS/DMAS^[105], SLSC plus filter DMAS^[110], MV plus DMAS^[118,119], CF plus DMAS^[104,130], and CF plus MV^{[117,127,131]}.

Representative work on DAS-based image reconstruction in PACT is summarized in Table 3. For completeness, Table 3 also includes relevant work in ultrasound or microwave imaging.

DAS-type algorithms can achieve approximate photoacoustic image reconstruction and are inexact reconstruction techniques. To achieve exact image reconstruction, advanced algorithms are needed. Filtered back projection (FBP) is a class of algorithms that are based on the Radon transform (see Sec. 2.4). It achieves image reconstruction by first filtering measured photoacoustic signals and then back-projecting the filtered signals into the image domain. The back-projection operation in FBP is similar to the reconstruction procedure in DAS-type algorithms. Therefore, the native DAS algorithm can be regarded as a simplified version of FBP, which achieves image reconstruction by directly back-projecting original photoacoustic signals into the image domain.

Early FBP algorithms were developed based on the condition of far-field approximation^[35], which states that if the distance between a detector and a photoacoustic source is much greater than the size of the photoacoustic source itself (i.e., |rs−rd|≫d, d: source size), the approximation |rd−rs|≈|rd−rs·(rd/|rd|)| holds (see Fig. 12). Under the condition of far-field approximation, the photoacoustic integral over a spherical shell can be approximated by the integral over its tangential plane. Consequently, image reconstruction in PACT can be achieved by inverting the linear Radon transform. Kruger et al. proposed an approximate FBP algorithm in 1995^[22], which is the first formal image reconstruction method in PACT. Xu and Wang developed approximate FBP algorithms from a more rigorous perspective in 2002^[134,135] and 2003^[136]. They deduced exact solutions to the image reconstruction problem and proposed approximate time-domain FBP algorithms for circular^[134], spherical^[135], planar, and cylindrical detection geometries^[136].

The preceding approximate FBP algorithms were generally developed for specific detection geometries. In 2007, Burgholzer and Matt extended the approximate FBP algorithms for image reconstruction in arbitrarily closed detection geometries under the assumption of the far-field approximation. The extended FBP reconstruction formula can be written as^[35]p0(rs)≅14π∫Sb(rd,t)δ(t−|rs−rd|v0)dΩ,(41)where p0(rs) is the reconstructed photoacoustic image, δ is the Dirac delta function, b(rd,t) is the back-projection term given by b(rd,t)=−2t∂p(rd,t)∂t,(42)and dΩ is the solid angle subtended by the element dσ of a detection surface and can be calculated by dΩ=dσ|rs−rd|2(nd·rs−rd|rs−rd|).(43)Here nd is the unit normal vector of the detector surface pointing to the region of interest (ROI).

The time-domain first derivative in Eq. (42) can be interpreted in the frequency domain as ∂p(rd,t)∂t=F−1{iωF{p(rd,t)}},(44)where F and F−1 denote the forward and inverse Fourier transforms, respectively; ω is a ramp filter, which suppresses the low-frequency contents of the measured photoacoustic signal p(rd,t) and amplifies high-frequency information. Since the value of ω extends from −∞ to ∞, the ramp filter is not integrable, and the inverse Fourier transform in Eq. (44) is undefined. To solve this problem, the ramp filter can be band-limited by a window function and Eq. (44) becomes ∂p0(rd,t)∂t=F−1{iωW(ω)F{p0(rd,t)}},(45)where W(ω) represents the window function. In practice, a smooth window such as the Hanning function is preferred over a box function because the latter may introduce undesirable ringing artifacts in the image domain.

One benefit of the far-field approximation in FBP is that it allows for simplified calculation of the solid angle dΩ. In other words, under the far-field approximation, dΩ in Eq. (43) reduces to dΩ≅−dσ|rd|2(nd·rd|rd|).(46)

Compared with Eq. (43), Eq. (46) involves only the detector position rd when calculating dΩ and is independent of the source position rs, thereby reducing the computational complexity. For example, assuming that the size of a 3D image to be reconstructed is Nx×Ny×Nz voxels (Nx=Ny=Nz=n) and the number of detectors used for imaging is M (M=n×n), the computational complexity is O(n5) for Eq. (43) and only O(n2) for Eq. (46). Furthermore, if the detection geometry is spherical, dΩ in Eq. (46) becomes a constant (dΩ=dσ/|rd|2), indicating that the computational complexity is O(1).

The approximate FBP algorithms are based on the condition of far-field approximation. This condition, however, may not be fully met in practice considering that photoacoustic signals may attenuate with propagation distance and that signal detection in space may be restricted. To solve this problem, in 2004, Finch et al. derived an exact FBP formula for image reconstruction in spherical detection geometries with odd dimensions (n≥3)^[23]. In 2005, Xu and Wang presented a universal FBP formula for image reconstruction in infinite planar, infinite cylindrical, and closed spherical detection geometries^[24]. The reconstruction formula is given as p0(rs)=∫Ωb(rd,t)δ(t−|rs−rd|v0)dΩΩ,(47)where the back-projection term b(rd,t)=2[p(rd,t)−t∂p(rd,t)∂t].(48)Ω is the solid angle of the detection surface with respect to the reconstruction point and equals 2π for the infinite planar geometry and 4π for the spherical and cylindrical geometries. FBP is one of the most commonly used algorithms in PACT.

By comparing Eqs. (47) and (41), we find that the exact FBP algorithm is very close to the approximate FBP algorithm in the formulas. They both consist of three reconstruction steps: filtering, back projection (the δ function), and summation. Their major difference is in the back-projection terms. The back-projection term b(rd,t) in the exact FBP algorithm has an extra term p(rd,t), which is not present in the approximate FBP algorithm. The reasons are as follows. According to the forward solution of the photoacoustic wave equation [Eq. (18)], the amplitude of the photoacoustic signals received by a detector is proportional to the size of the source and inversely proportional to the detection distance, i.e., p(rd,t)∝d/|rd−rs|. Under the condition of the far-field approximation, i.e., |rd−rs|≫d, the amplitude of p(rd,t) is small enough compared with the derivative term −t∂p(rd,t)/∂t in Eq. (48) to be ignored without significant loss of reconstruction accuracy. Therefore, the back-projection term b(rd,t) in the approximate FBP algorithm does not involve the term p(rd,t).

Figure 13 presents an example demonstrating the principle of FBP-based image reconstruction. In this example, the photoacoustic source is a 5-mm-diameter uniform spherical absorber, and the photoacoustic signals generated from the source are recorded by a 40-mm-diameter circular detector array [Fig. 13(a)]. Figures 13(b) and 13(c) show the representative photoacoustic signal p(rd,t) measured by a detector and the corresponding back-projection signal b(rd,t), respectively. Figure 13(d) shows the projection images of detectors at different positions. Figures 13(e)–13(g) are the reconstruction results obtained by summing the projection images of four, 16, and 256 detectors, respectively. The results demonstrate that the FBP algorithm can effectively reconstruct the structure and amplitude information of the photoacoustic source. Note that the artifacts in the reconstructed images are not caused by the FBP algorithm but by the limited view angle of the circular detector array^[137].

To evaluate the performance of the FBP algorithm for different detection geometries, a group of image reconstruction simulations was conducted. In the simulations, a multi-sphere phantom is used as the photoacoustic source. The phantom contains nine spherical absorbers with unit intensity, among which eight with diameters uniformly varying from 1 mm to 2 mm are evenly distributed on a circle with a diameter of 10 mm, while the ninth one, with a diameter of 1.6 mm, is seated at the origin. Figures 14(a)–14(c) show the relative positions of the multi-sphere phantom and a finite planar, finite cylindrical, and closed spherical detection surface. The planar detection surface is located 12 mm below the x–y plane, and the cylindrical and spherical detection surfaces are centered at the origin. The three detection surfaces have the same number (32768) of evenly distributed point-like detectors and approximately equal detection areas (plane: 150mm×150mm; cylinder: 70 mm base diameter, 100 mm height; sphere: 84 mm diameter). The reconstructed images of the x–z and x–y cross sections of the phantom are displayed in Figs. 14(d)–14(f) and 14(g)–14(i), respectively. The results show that the FBP algorithm achieves stable image reconstruction for all three detection geometries. However, the image reconstructed in the closed spherical detection surface has much fewer artifacts than those reconstructed in the planar and cylindrical detection surfaces. This is because the spherical detection surface is closer to an ideal detection geometry than the finite planar and finite cylindrical surfaces.

To extend the application scenarios of FBP, many in-depth studies have also been carried out^[138]. For example, in 2007, Kunyansky proposed a set of FBP-type formulas for image reconstruction in spherical detection geometries with arbitrary dimensions (n≥2)^[139]. In the same year, Finch et al. also developed a set of FBP-type inverse formulas for spherical detection geometries with even dimensions^[140]. In 2009, Nguyen derived a family of inverse formulas for thermoacoustic tomography, in which many previously known FBP inverse formulas can be obtained as special cases^[141]. In addition, the exact FBP algorithms mentioned above are primarily used for planar, cylindrical, and spherical detection geometries with point-like detectors. In 2007, Burgholzer et al. developed a two-step FBP algorithm for integrating line detectors^[142]. In 2012, Natterer proposed a novel FBP algorithm for elliptical detection geometries^[143], which was further developed by Palamodov in 2012^[144], Haltmeier in 2014^[145,146], and Salman in 2014^[147].

Moreover, the FBP algorithm is well suited for parallel computing. GPUs have been used to achieve real-time FBP image reconstruction in both 2D and 3D imaging^{[42,148–153]}. A field-programmable gate array (FPGA) was also employed to accelerate image reconstruction in low-cost PACT systems^[154].

Table 4 summarizes representative work on FBP-based image reconstruction in PACT.

The basic principle of series expansion (SE) algorithms is to approximate the image to be reconstructed using mathematical series. Compared with other reconstruction algorithms in PACT, SE algorithms are simple and fast for specific detection geometries such as planar geometries because they can be implemented using FFT.

Norton et al. proposed exact SE methods for the reconstruction of acoustic reflectivity in circular geometries in 1980^[156] and in planar, cylindrical, and spherical geometries in 1981^[157]. Several groups have studied similar ideas for image reconstruction in PACT. Köstli et al. in 2001 proposed an exact SE algorithm for image reconstruction in planar detection geometries^[31], and Xu et al. in 2002 reported a similar algorithm^[158].

Mathematically, SE-based PACT image reconstruction for planar detection geometries can be formulated as^[31,83]P0(kx,ky,ω)=v02kz2ωFx,y,t{p(x,y,t)},(49)P0(kx,ky,ω)→ω2=v02(kx2+ky2+kz2)P0(kx,ky,kz),(50)p0(x,y,z)=Fx,y,z−1{P0(kx,ky,kz)},(51)where p(x,y,t) is the photoacoustic signal measured by a planar surface at position (x,y) and time t; kx, ky, and kz are the spatial wavenumber components in each Cartesian direction; ω is the temporal frequency; v0 is the SOS; → represents the interpolation operation between the temporal and spatial frequencies ω and kz; p0(x,y,z) and P0(kx,ky,kz) are the initial photoacoustic pressure and its spatial Fourier transform; and F and F−1 represent the forward and inverse Fourier transforms, respectively. The reconstruction procedure above involves one interpolation and two Fourier transforms and has a computational complexity of O(n3logn) for a 3D image with a size of n×n×n voxels by use of FFT^[159]. The numerical implementation of the SE algorithm is available in the k-Wave toolbox^[83].

Figure 15 is a numerical example showing SE-based PACT image reconstruction using the k-Wave toolbox. The photoacoustic source in the xy plane contains nine spherical absorbers with unit intensity. Among them, eight absorbers with diameters uniformly varying from 1 mm to 2 mm are evenly distributed on a circle with a diameter of 10 mm, while the ninth absorber with a diameter of 1.6 mm is seated at the origin. The signal detection geometry is a planar surface 12 mm below the xy plane. The planar surface has 364×364 point detectors uniformly distributed spanning a physical size of 36mm×36mm. Figures 15(b) and 15(c) are the xz and xy cross sections of the photoacoustic source, respectively, and Figs. 15(d) and 15(e) are the corresponding reconstructed images, which show that the SE algorithm can recover major structures of the source. Note that the distortion and blurring in the reconstructed images are due to the finite size of the planar detection surface.

In 2002, the Wang group reported exact SE algorithms for PACT image reconstruction in cylindrical and spherical geometries^[135,160]. However, the image reconstruction procedures in these two cases involve complicated mathematical calculations, preventing their implementations using FFT and are thus time-consuming^[135,160]. Based on Norton’s pioneering work on circular geometries^[156], Haltmeier and Scherzer in 2007 proposed a 3D reconstruction algorithm for cylindrical detection geometries, where photoacoustic signals are recorded by linear integrating detectors^[161]. The proposed algorithm has a computational complexity of O(n4) for a 3D image with a size of n×n×n voxels, which is higher than that of the SE algorithm in planar geometries [O(n3logn)]^[159] but lower than that of FBP [O(n5)]^[24].

In 2012, Kunyansky proposed fast image reconstruction algorithms suitable for circular, cylindrical, and spherical detection geometries^[33]. In Kunyansky’s work, image reconstruction is based on the Fourier transforms of the initial photoacoustic pressure p0(rs) and the measured photoacoustic signals p(rd,t), as summarized below. Suppose that the detection geometry is a circle, as shown in Fig. 7(c). By expanding the Fourier transform P(rd,ω) of the measured photoacoustic signal p(rd,t) and the initial photoacoustic pressure p0(rs) in the Fourier series in φ and θ, we have^[33]P(rd,ω)=∑k=−∞∞Pk(ω)eikφ,rd=(rdcosφ,rdsinφ),(52)p0(rs)=∑k=−∞∞p0k(rs)eikθ,rs=(rscosθ,rssinθ),(53)where i is the imaginary unit, ω is the angular frequency, rd=|rd| is the radius of the detection circle, rs=|rs| is the distance from the reconstruction point to the center of the circular detection geometry, and Pk(ω) and p0k(rs) are expansion coefficients given by^[33]Pk(ω)=12π∫02πP(rd,ω)e−ikφdφ,(54)p0k(rs)=2π∫0∞Pk(ω)H|k|(1)(ωR)J|k|(ωrs)dω,(55)where H|k|(1) is the Hankel function of the first kind of order k and J|k| is the Bessel function of the first kind of order |k|^[162]. Equation (55) is somewhat similar to an expression in Norton’s work^[156], where the Bessel function J|k| rather than the Hankel function H|k|(1) is in the denominator and a term corresponding to the real part of Pk(ω) in the numerator.

By discretizing Eq. (55), we can obtain p0k(rs) and then the initial photoacoustic pressure p0(rs) [Eq. (53)]. However, direct discretization of Eq. (55) results in a computational complexity of O(n2) for each p0k(rs) and O(n3) for the whole reconstruction. To achieve fast reconstruction, substituting Eq. (55) into Eq. (53) yields^[33]p0(rs)=12π∫R2P0(Λ)eirs·ΛdΛ,Λ=(ωcosφ,ωsinφ),(56)where P0(Λ) is given by P0(Λ)={2π∑k=−∞∞(−i)kp0k(ω)ωH|k|(1)(ωR)eikφ,Λ≠0,∫0∞2Pk=0(ω)πωH0(1)(ωR)RJ1(ωR)dω,Λ=0,(57)where R is the radius of the circular detection geometry. Equation (56) indicates that the initial photoacoustic pressure p0(rs) can be obtained by calculating the 2D inverse Fourier transform of P0(Λ), which has a computational complexity of O(n2logn), lower than that of the direct-discretization-based reconstruction [O(n3)].

For cylindrical detection geometries with linear integrating detectors, image reconstruction can be readily realized by combining the fast image reconstruction procedure for 2D circular geometries and the 3D inverse Fourier transform^[33]. This procedure can yield fast reconstruction with a computational complexity of O(n3logn). For spherical detection geometries, the derivation of the image reconstruction procedure is similar to those for circular geometries [Eqs. (52)–(57)]. A prominent difference is that p0(rs) and P(rd,ω) are expanded into spherical harmonics. In addition, the fast spherical harmonics (FSH) transform is adopted to achieve fast reconstruction^[163], which results in a reconstruction complexity of O(n3log2n)^[33].

Wang and Anastasio in 2012 demonstrated that a mapping relationship exists between the spatial frequency components of initial photoacoustic pressure p0(rs) and the temporal frequency components of measured photoacoustic signals p(rd,t)^[164]. They thus proposed a Fourier-transform-based image reconstruction algorithm whose computational complexity is O(n2logn) for 2D image reconstruction in circular geometries and O(n3logn) for 3D image reconstruction in spherical geometries^[165]. The reconstruction formula does not involve series expansion of special functions or multi-dimensional interpolation operations in Fourier space, which are commonly used in previous work.