Quantitative Photoacoustic Tomography Using Iteratively Refined Wavefield Reconstruction Inversion: A Simulation Study

The ultimate goal of photoacoustic tomography is to accurately map the absorption coefficient throughout the imaged tissue. Most studies either assume that acoustic properties of biological tissues such as speed of sound (SOS) and acoustic attenuation are homogeneous or fluence is uniform throughout the entire tissue. These assumptions reduce the accuracy of estimations of derived absorption coefficients (DeACs). Our quantitative photoacoustic tomography (qPAT) method estimates DeACs using iteratively refined wavefield reconstruction inversion (IR-WRI) which incorporates the alternating direction method of multipliers to solve the cycle skipping challenge associated with full wave inversion algorithms. Our method compensates for SOS inhomogeneity, fluence decay, and acoustic attenuation. We evaluate the performance of our method on a neonatal head digital phantom.

Generally, photoacoustic (PA) image reconstruction algorithms rely on the assumption that the imaging area is homogeneous [12], [13], [14], which implies that the speed of sound (SOS) is constant throughout the imaging region.This is an invalid assumption in most biological tissue (for example, brain tissue has many regions with different SOS) [15], [16], leading to distortion in the received PA signals, and consequently in reconstructed PA images.There have been several studies that aim to implement SOS estimation [12], [14], [17], [18], [19].Xu and Wang [12] used an iterative algorithm to compensate for SOS variation in breast tissue.Deán-Ben et al. [13] used prior information on the position of acoustic deformities, which is inferred from usual anatomical data or another imaging method to improve image quality.Zhang and Wang [19] proposed an algorithm based on the correlation of PA recorded signals to neutralize heterogeneity in acoustical properties.Treeby et al. [14] used the autofocus strategy to automatically estimate SOS in the medium by maximizing image sharpness.In Jose et al. [17], the reconstruction algorithm relies on refracted ray paths rather than direct ray paths to obtain more realistic SOS maps.Recently, Yin et al. [18] developed an algorithm that utilizes an iterative procedure to estimate an SOS map relying on the fact that the recorded signals by two neighboring ultrasonic transducers are very similar in the received waveforms.Sun and Sun [20] developed a method to estimate both the optical absorption coefficient and the SOS map simultaneously.
In contrast to methods that either use part of the recorded data, e.g., travel time or amplitude, or make unrealistic assumptions about the medium (e.g.constant speed of sound, homogeneity, etc.), full waveform inversion (FWI) exploits the full information contained in the recorded wavefields to estimate the constitutive properties of a medium without any assumption about its level of heterogeneity or requirement of a level of homogeneity [21], [22].This enriched information guarantees superior performance of wave-based tomography relative to travel time-based approaches in terms of ability to estimate several constitutive properties such as SOS, density, and acoustic attenuation.Considering emerging high-performance computing facilities and the current state of theoretical understanding of this methodology, it is reasonable to anticipate that wave-based tomography will rapidly become the preferred method to perform accurate tomography [7].Algorithms using full waveform information are capable of producing images with higher resolution (on the order of one wavelength) than algorithms using ray-based methods [23].
FWI for brain imaging is challenging due to both large variability in SOS within the brain tissue and the high SOS difference between the soft brain tissue and the surrounding skull.Classical FWI, which is based on the iterative least-squares minimization of the differences between the measured data and the numerically simulated counterparts, requires prior information about the skull to reconstruct an accurate SOS map of the brain.To relax this requirement, several approaches extend the search space of classical FWI by incorporating more convex distances between the recorded and simulated data.One such approach is a so-called adaptive waveform inversion (AWI) which has been used to reconstruct the SOS map of the human brain without resorting to prior information about the skull [24].
In this paper, we propose a qPAT method that calculates DeACs by estimating SOS, acoustic attenuation, and fluence maps.Our method uses iteratively refined wavefield reconstruction inversion (IR-WRI), which incorporates the alternating direction method of multipliers (ADMM) to solve the cycle skipping challenge associated with FWI algorithms.We evaluate the performance of our qPAT method on a digital phantom of a neonatal head.

II. NOTATION
The mathematical symbols adopted in this paper are as follows.We use italics for scalar quantities, boldface lowercase letters for vectors, and boldface capital letters for matrices and tensors.We use the superscript T to denote the adjoint of an operator.The ith component of the column vector x is shown by x i and its absolute value is returned by |x i |.
For the real-valued n-length column vectors x and y, the dot product is defined by ⟨x, y⟩ = x T y = n i=1 x i y i and their Hadamard product, denoted by x • y, is another vector made up of their component-wise products, i.e. (x • y) i =x i y i .The ℓ2-and ℓ1-norms of x are, respectively, defined by ||x||

A. Digital Phantom of Neonate Head and Simulation Setup
We created a digital phantom of a 30-week-old neonatal head by averaging 324 MRI atlases made by Brigadoi et al. [58].The model is composed of extra-cranial tissue (ECT), cerebrospinal fluid (CSF), gray matter (GM), and white matter (WM).Assigned optical properties for all layers were taken from the literature [59], [60], [61].ECT includes scalp and skull; the averaged optical properties of scalp and skull were assigned to this layer.To fully mimic the neonatal head, we included complex vasculature as the primary imaging target of the phantom.The size of the computational grid for this head model is 1000 × 1000 pixels, in which pixel size is 120 microns.In this work, 64 optical fibers with size of 0.6 mm are equally interspersed between 128 US transducers (active aperture = 1 mm) to form a ring (diameter of 120 mm) around the digital phantom.Fig. 2 shows the regions in the digital phantom of the neonatal head as well as the locations and arrangement of the optical fibers and US transducers.

B. SOS and Acoustic Attenuation Estimation
We implemented FWI in the frequency domain as opposed to the more conventional time domain formulation.The frequency domain decreases the nonlinearity of the problem, which helps to bypass the local minima (cycle-skipping) problem.It is also computationally efficient when the inversion can be limited to a few discrete frequencies.This is the case when the acquisition device (in our case, a ring array) provides a redundancy of information due to abundant illumination of the target in both the US and PA setups.The frequency domain FWI can be formulated as the following multivariate constrained optimization problem [62] minimize R(m) where R(m) is a regularization function and the constraints are the wave equation A(m, ω)U(ω) = B(ω) and the observation equation PU(ω) = D(ω).Moreover, m ∈ C N×1 denotes complex-valued model parameters, constructed by a mapping between SOS and acoustic attenuation.The mathematical model describing the attenuating medium is given in Appendix A. In addition, N is the number of discretized points of the medium, ω is the angular frequency, n f is the number of frequencies, U(ω) = [u 1 (ω), u 2 (ω), . . ., u n s (ω)] ∈ C N ×n s is a matrix, whose columns represent the wavefields triggered by each source, A(m, ω) = + ω 2 diag(m)∈ C N ×N is the impedance matrix in the acoustic environment, is the discretized Laplace matrix, and diag (•) is a diagonal matrix with • in its diagonal.In addition, ×n s is the matrix of recorded data, n r is the number of receivers and P ∈ R n r ×N is the linear observation (sampling) operator that samples the wavefields U at receiver positions.Finally, M is a convex set defined according to our prior knowledge of m.For example, if we know the lower and upper bounds on m then ( The level of non-linearity of FWI in the optimization problem ( 1) is moderate [63], which means that local optimization methods, like Newton-family solvers [64] even with an inaccurate starting model, converge to a reasonable solution.
To solve the optimization problem (1) in the full parameter space, we should deal simultaneously with all the optimization variables (U(ω), m).This approach which is called all-at-once [65] is not feasible for practical applications of FWI because of massive computational burden and memory requirements.The classical formulation of FWI was implemented on the reduced parameter space to keep the computational burden and memory requirements reasonable [22], [66].In classical FWI, the wave-equation constraint is solved for U(ω), i.e.U(ω) = A(m, ω) −1 B(ω), and implemented in the second constraint, i.e.PA(m, ω) −1 B(ω) = D(ω), which implies that the optimization variables reduce down to the constitutive parameters of the medium.In this case, the inverse problem is highly nonlinear due to dependency on A(m, ω) −1 , meaning that local optimization methods require a good starting model to converge to the global minimum.Without a good starting model, the problem is trapped in a local minimum, a process which is called cycle skipping [67].In order to solve the cycle skipping challenge of FWI, we recently proposed IR-WRI [68], [69] which uses an all-at-once approach for (1) in which the optimization variables are updated in an alternating mode based on ADMM [70] instead of simultaneous update of all optimization variables.IR-WRI mitigates cycle skipping by keeping the wavefield as an independent optimization variable.In this method, first the wavefields that jointly satisfy the observation equation and the wave equation are computed in a least-squares sense.Then, the constitutive parameters of the medium are updated by minimizing the wave-equation errors.Then, the Lagrange multipliers, which are required for solving a constraint optimization, are updated.We solve (1) with an augmented Lagrangian method [64] whose objective function is given by: Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
+ λ 1 2 where V(ω) ∈ C N ×n s and W(ω) ∈ C n r ×n s are Lagrange multipliers and the scalars λ 1 , λ 2 > 0 are the penalty parameters assigned to the wave equation and the observation equation constraints, respectively.The primal variables m and U(ω) and the dual variables (i.e., the Lagrange multipliers) V(ω) and W(ω) can be updated in alternating mode within the framework of ADMM to break the full multivariate problem into several manageable subproblems [68].The stationary point of ( 3) is determined iteratively as where ω = ω 1 , ω 2 , . . ., ω n f .The wavefield U(ω) k+1 in (4a) can be obtained in a closed-form expression as given in [68].
Once the wavefields are reconstructed, the model parameters are updated by minimizing the wave-equation errors keeping the wavefields fixed (4b).The ill-posedness of any FWI method, including IR-WRI, also requires regularization to remove the undesired features in the reconstructed model that may be due to inhomogeneous illumination, sparse acquisition, or parameter cross-talk during multiparameter reconstruction [71].
A suitable regularization renders the solution unique, increases its stability, and prevents data over-fitting that may cause coherent noise in the reconstructed model.The regularization should be able to (mathematically) describe the solution while being easy to implement [72].Conventional regularizations such as those promoting smoothness, sparseness, and blockiness are derived according to some prior assumptions about the targeted model and are independent of the data.In contrast, plug-and-play regularizations (such as those based on nonlocal means filters and block matching 3D filters (BM3D) [73]) are more flexible for regularizing complex models.Another rationale for using BM3D is that proximal methods, such as ADMM [70], provide a versatile and computationally efficient approach to implement non-differentiable and hybrid regularization functions (such as BM3D) through variable splitting, as reviewed in [71].Here, we use BM3D regularization for SOS and acoustic attenuation estimation.SOS and acoustic attenuation can be easily represented by one complex-valued frequency-dependent term (see Appendix A).In [69], Aghamiry et al. presents a method to process this complex-valued parameter as the optimization variable in IR-WRI, instead of processing SOS and attenuation as two

C. Derived Absorption Coefficient Calculations
1) Image Segmentation: In order to calculate the fluence map and the DeACs inside a neonate head, we start by identifying the topological structure including the boundaries of the ECT, CSF, GM, WM, and the vessel regions.These regions are specified by applying a segmentation method, the nearest neighbor algorithm described in [74], on the IR-WRI derived SOS map.As published in the literature [75] and data available on the website (https://itis.swiss), the following reference SOS values are used: 2300 m/s, 1600 m/s, 1650 m/s, 1700 m/s, 1510 m/s for ECT, CSF, GM, WM, and vessel regions, respectively.With the segmentation algorithm, each pixel is assigned to the closest reference SOS.This results in identification of five segmented regions.As an example, if one pixel of the estimated SOS was 1655m/s, the reference value of 1650m/s is the closest value and therefore GM is assigned to that pixel.IR-WRI could theoretically estimate a DeAC for each of the segments described above.However, in reality, the tissue is optically heterogeneous.Determining only five DeACs, one for each of the five segmented head regions, is a rough estimation.To better estimate a real-life scenario of varying DeAC in the brain tissue, we further divide the main regions (i.e., ECT, CSF, GM, WM, and vessel regions) into subregions (n g = total number of subregions or segments) by overlaying a grid structure and using regional divisions to subdivide each gridded square (see Fig. 5(b)).We then determine DeAC for each subregion separately.
2) Fluence Calculation: Fluence refers to the amount of optical energy delivered to tissue per unit area.MCX, also known as Monte Carlo eXtreme, is a software specifically developed for simulating the transport of photons in turbid media using the Monte Carlo method.It specializes in time-resolved photon transport simulations [76].In our study, we used MCX (with the input parameters in Table I) to calculate fluence.Iteratively, fluence was calculated based on scattering coefficients (µ s ) and anisotropy values (g) extracted from literature (referenced as [77] and listed in Table II), and DeACs (generated by IR-WRI in each iteration) for various subregions.The simulation involved 640 million photon packets and utilized a grid size of 120 microns.
3) DeAC Estimation: Generally, the imaging problem is very sensitive to variations in the source term b, which is the element-wise product of fluence and DeAC as defined in (5).However, there is a significant trade-off between b and the properties of the medium when both quantities need to be estimated.This ill-posedness prevents assignment of a degree of freedom to each sample of the medium for DeAC.Instead, a parsimonious parametrization of DeAC is preferred.Thus, in order to reduce the dimensionality of search space, we use  [77] the segmented image (section III-C.1)based on the SOS map obtained from IR-WRI (section III-B).The PA source or initial pressure b ∈ C N ×1 can be decomposed as where F ∈ R N×N is the diagonal form of the fluence that was estimated in section III-C.2, and the columns of M ∈ R N×ng contain the location of the segments according to the information extracted in section III-C.1.Also, x ∈ R ng×1 is a vector of attenuation coefficient values for different segments.
We direct interested readers to [78], [79], and [80] for further information on the discretized Helmholtz equation ( 1) and source model (5).Next, we demonstrate the IR-WRI method to extract values of DeACs in this parsimonious parametrization.
Frequency-domain FWI with known model parameters m 0 and unknown PA source b can be written as: minimize R(x) where R(x) is an appropriate regularization function of x, u(ω) ∈ C N ×1 and d(ω) ∈ C N ×1 denote the wavefield and the recorded data for frequency ω, respectively, and ω1 , ω2 , . . ., ωn f indicate the frequency content of the PA signal.Also, A(ω) ≡ A(m 0 , ω), and the dependency to the passive SOS model is removed (for the sake of simplicity).Determination of the optimum solution of ( 6) is extremely difficult and ill-posed, requiring sophisticated regularizations.In this study, we have used the total variation (TV) regularization [81] for x, i.e.R(x) = ∥ ∇x ∥ 1 where ∇ is the first-order difference operator.The augmented Lagrangian function associated with the problem in ( 6) is given by where the scalars γ 1 , γ 2 > 0 are the penalty parameters assigned to the wave equation and the observation equation constraints, respectively, and v(ω) ∈ C N ×1 and w(ω) ∈ C nr ×1 are the Lagrange multipliers.Beginning with the model m 0 obtained in section III-B, v 0 (ω) = 0, w 0 (ω) = 0 ∀ω and initial value of absorption coefficients x 0 , ADMM solves the multivariate optimization problem (7) iteratively as follows: where ω = ω1 , ω2 , . . ., ωn f .Problem (8a) is another wavefield reconstruction problem (a problem that jointly satisfies the wave equation and the observation equation in a least-squares sense), whose solution can be obtained in a closed form similar to the one of (4a).The optimization problem (8b) reads By adding and subtracting the term ||w(ω) k || 2 2 to (9), we have where we have ignored the −||w(ω) k || 2 2 -term as it does not impact the optimization result.Equation ( 10) is a one dimension TV regularization and thus can be solved efficiently [82].TV regularization has been used to control the large null space of the inversion.As previously mentioned, this algorithm requires an initial guess of DeACs (x 0 ).It is estimated jointly with a set of monochromatic wavefields u 0 (ω).To this end, we may neglect the regularization term in (7) and solve the optimization problem jointly for u(ω) and xfor a few frequencies.The closed-form expression of this system for r frequencies, i.e., ω1 , . . ., ωr can be written as: Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
After solving (11) for x 0 , it, along with the extracted SOS and acoustic attenuation from section III-B can be used as inputs to subproblems (8) to extract accurate DeACs.In summary, our method involves the following steps: • IR-WRI [68] is employed to find the SOS and acoustic attenuation maps based on US data.
• Pixels of the estimated SOS map are segmented into five major regions of the head.
• Finally, IR-WRI is employed to estimate DeAC from the MCX-updated fluence map and the acquired PA data in which the SOS and acoustic attenuation maps.Fig. 3 shows the block diagram explaining the methodology we have used.

For all of the numerical computations, a computer with an
Intel Xeon Platinum 8260 24C/35.75M/2.40G/165Was a CPU and 16x 64GB DDR4-2933 RDIMM as RAM was used.Also, the coefficients optimized nine-point stencil finite-difference method implemented with anti-lumped mass and perfectly matched layer (PML) absorbing boundary conditions [78] was used to perform forward modeling for all the tests.

A. Estimation of SOS and Acoustic Attenuation
The acquisition device is comprised of 128 single element transducers of size 1 mm which collect data in the frequency band from 0.3 to 1 MHz.The direction of each transducer points towards the center of the ring.We use a set of frequencies, f = [0.30.4 0.7 0.8 1] MHz, for mono-frequency IR-WRI followed by a multiscale frequency continuation strategy.We perform three passes through the frequencies, using the final model of one pass as the initial model for the next one.The starting and finishing frequencies of the three passes are [0.3,0.4], [0.3, 0.4], and [0.3, 1] MHz, respectively.The stopping criterion with respect to iterations is IR-WRI must stop after a maximum of 30, 25, and 15 iterations for the first, second, and third passes, respectively.The hyperparameters (λ 1 and λ 2 ) are tuned based on a small fraction, e.g.0.001, of the highest eigenvalue of A(m 0 , ω 1 ) −T P −T PA(m 0 , ω 1 ).The readers can refer to appendix C of [71] for more details.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
is shown in Fig. 4(h) at different iterations in blue, and the same plot for the estimated acoustic attenuation map (Acoustic Attenuation Error) is shown in red.The computational time for 185 iterations was 25 hours.As seen, IR-WRI is able to provide a highly accurate SOS map.

B. DeAC Calculation
1) Image Segmentation: As explained in section III-C.1, a segmented image of the head was obtained using the SOS map estimated by IR-WRI.Fig. 5(a) shows the five regions (ECT, CSF, GM, WM, vessels) extracted from the SOS map.Results were used to calculate fluence maps as well as DeACs.
2) Fluence Calculation: The fluence map was calculated iteratively based on the specifications of the optical fibers and the optical properties of the neonatal head phantom at a wavelength of 800 nm (refer to Tables I and II in the Methods section) and the DeAC values were acquired from IR-WRI.Fig. 5(c) illustrates the calculated fluence map at the first iteration, while Fig. 5(e) represents the fluence map at the final iteration.It is evident from these figures that the optical scattering effect in biological tissue results in a significant decrease in the optical energy deposited in the central part of the head.To provide a clearer visualization of the fluence distribution, particularly in deeper areas, Fig. 5(d) and 5(f) display the logarithmic scale representation of the fluence map at the first and last iterations, respectively.
3) DeAC Estimation of Subregions: One hundred frequencies selected evenly over 0.3 to 1 MHz are used.To begin, the first 20 frequencies are used to solve (11) for extracting an initial DeAC (x 0 ).Then, all 100 frequencies are used simultaneously for solving subproblems (8) using the SOS and acoustic attenuation estimation maps shown in Fig. 4(c) and 4(f) and the fluence map estimation shown in Fig. 5(c).Tuning hyperparameters like γ 1 and γ 2 in (7) were challenging.We set them according to the following recipe, but of note, it is not the optimal one.The hyperparameter γ 1 is tuned as a small fraction, e.g.0.001, of the highest eigenvalue of M T F T FM.Then γ 2 is tuned such that γ = γ 1 /γ 2 will be a small fraction of the highest eigenvalue ξ of A( ω1 ) −T P T PA( ω1 ).We set γ = 10 −4 ξ for noiseless data and increase it when analyzing noisy data to prevent data overfitting.For example, we used γ = 0.01ξ for the noisy data.
The true DeAC map used in the simulation is shown in Fig. 5(g).Next, an estimated DeAC map based on the SOS, acoustic attenuation and fluence maps that undergoes updates in each iteration is computed.Since the model we used in this study had only 5 regions with different DeACs, and real-case scenarios, because of the heterogeneity of tissue, will contain more regions, we decided to initially assume the possibility of the existence of up to 200 different DeACs and test to see if this process will correctly converge on the five DeACs in the true DeAC map.In order to quantify the accuracy of the reconstructed DeACs, we define DeAC percentage error.The DeAC percentage error for the full image is 9.41% and it is shown for each iteration in fig.5(j).The DeAC percentage error for each of the five regions is ECT=9.81%,CSF=9.62%,GM=9.01%,WM=10.40%, and vessel=10.25%.The results Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.were obtained by running iterations until the rate of change in the error was < 1%: these conditions were satisfied beyond the 252 nd iteration.The results for the 1 st and 252 nd iterations are shown in Figs.5(h)-5(i).The computational time for DeAC calculation for 450 iterations was ∼49 hours.
For further assessment of the proposed method, we compare our results with the least squares (iterative) time reversal reconstruction algorithm [83], [84] which selected for comparison because it and the IR-WRI algorithms use the same wave equations and similar starting data.The time reversal algorithm was run for 100 iterations and stopped when the error did not decrease.We reconstructed the PA image and then divided it by the fluence map to generate a time-reversal (TR) DeAC map.Fig. 6(a) utilizes a homogeneous SOS = 1500 m/s and there is no adjustment for acoustic attenuation.and for each of the five regions is ECT=82.51%,CSF=86.3%,GM=80.6%,WM=70.3%, and vessel=40.2%.
Comparing the DeAC map from our proposed algorithm (Fig. 5(i)) and those extracted using the iterative time reversal algorithm (Fig. 6(b)), shows our proposed algorithm can recover the DeAC map nearly perfectly, while the iterative time reversal algorithm, using nearly the same starting information, poorly recovers the DeAC map, especially in deeper regions.
We next assess the robustness of our method against random noise in the data.We added different amounts of noise level to the raw PA data (noise level =1 to 5 (Fig. 7(d)) where noise level 5 represents similar amplitude for noise and signal, and noise levels 1-4 represent 20% to 80% the amplitude of noise level 5.We next performed IR-WRI, and DeAC error varied from 11.6% to 56.1% (Fig. 7(c)).The DeAC maps for noise level 1 and noise level 5 are shown in Fig. 7(a) and 7(b).

V. DISCUSSION
Many commonly-used ray-based PA image reconstruction algorithms (e.g., universal back-projection [85], delay and sum [86] and multiply delay and sum [9]) assume that SOS and fluence maps are homogeneous across the entire imaging area and there is no acoustic attenuation.Therefore, to estimate DeACs, insufficient information is available, leading to large error.Prior information from a different imaging modality can be used to decrease DeAC estimation error [30], [87], [88].Because of the use of US transducers in PA data collection, US-PA dual modality is a natural choice to obtain such prior information.On the other hand, in wave-based methods, the full information contained in the recorded wavefields, i.e., amplitudes and phases, is considered without any assumption about level of heterogeneity [21], [22].Therefore, these techniques can be effective for estimating accurate SOS and acoustic attenuation maps from US data.Several algorithms have been proposed for this very purpose [89], [90], [91].Some limitations of these US methods are sensitivity to the selection of the initial model [90], massive computational burden [92], and suboptimal reconstruction of the acoustic attenuation map [93].
Several studies have attempted to estimate the SOS, the DeAC, or optical absorbance density [94], [95], [96], [97], [98], [99].However, a few limitations associated with these methods (but more recently resolved) are: utilizing the Born approximation implies assuming SOS(x) = SOS 0 (x) + SOS perturbation (x), where the magnitude of SOS perturbation / SOS 0 is significantly smaller than 1 [97].However, this assumption does not hold in our case due to the presence of the skull, where the SOS is approximately 2400 m/s while SOS 0 is around 1500 m/s.Another limitation is a requirement for incorporating prior knowledge about the geometry of the SOS map [99].
To estimate DeACs, several qPAT algorithms have been proposed [38], [40], [41], [43], [53], [100], [101].Some historical limitations of these qPAT methods are: The assumption that the transducers have a point-like nature and lack directivity, with a constrained aperture size [40]; uniform distribution of acoustic properties within the tissue, which may not always be the case in practice [38]; prior knowledge of the speed of sound in the medium is required [43]; use of diffusion approximation [100] to calculate fluence, which is not a correct approximation in the head due to the presence of cerebrospinal fluid (CSF) region, which acts as a non-scattering medium and poses a challenge for this approach.In fact, some older qPAT methods eschew fluence modeling entirely, and assume that the fluence is uniform throughout the tissue [7], [27], [29].Moreover, recently, several research studies have made use of minimization-based methods [38], [39], [40], [101], [102], [103].Nonetheless, the primary drawback associated with these newer approaches is their tendency to exhibit escalating computational complexity, ultimately resulting in substantial computation time requirements [104].The time it takes for calculations depends on various factors, including the number of photon packets traced, the specifications of the hardware used, the size of the computational area, the quantity of transducers, and the complexity of the phantom being studied.For example, in [38], during the inversion process, 10 7 photon packets were utilized to compute radiance and fluence, while 5 × 10 6 photons were employed to calculate the corresponding adjoint quantities.On a highend consumer GPU such as the NVIDIA GeForce Titan X Pascal, a single iteration of the inversion process, including the adjoint model, took approximately 84 seconds.On a different note, as described in [40], computations were performed on a computer equipped with two Intel(R) Xenon(R) Gold 6136 processors, possessing a total of 24 cores and 256 GB of memory.Each reconstruction step took around 150 minutes, with approximately 20 to 30 minutes allocated for tasks other than the primary computation.
Classical FWI can also be used for DeAC estimation.In this case the wave-equation constraint in ( 6) is solved for u(ω), i.e. u(ω) = A(ω) −1 FMx, and substituted into the second constraint in (1).The resulting optimization problem will be a regularized least-squares time-reversal.This optimization problem cannot be solved easily because of the presence of A(ω) −1 , which is a very large and dense matrix that cannot be stored or directly inverted to extract x [105].Using the second IR-WRI algorithm to iteratively optimize parameters to include for fluence correction can be implemented with no initial guess or boundary assumption for the fluence, or preconceptions for the iteration process.This task would be infeasible for existing FWI-based methods because the nonlinearity of fluence compensation would set up another requirement for inverting a large and dense matrix.
In this paper we describe a qPAT framework to estimate DeACs.The requirements for this method are as follows: • US data acquired by a ring of transducers when each transducer is activated as a transmitter in turn, • Optical source locations and specifications, and • PA data acquired by all transducers when the imaging target is irradiated by a laser beam.The method operates by first calculating the SOS and acoustic attenuation maps using the acquired US data through the IR-WRI algorithm (with BM3D regularization).
It should be noted that the calculated acoustic attenuation map (error < 14%) is less accurate than the calculated SOS map (error < 5%).This is because the US recorded data is less responsive to acoustic attenuation compared to SOS.Subsequently, the IR-WRI algorithm is used again (with TV regularization) to solve a source estimation problem on the PA data (using the obtained SOS, acoustic attenuation and fluence maps).The PA problem is then solved by segmenting the SOS map to compensate for the trade-off between not knowing the pixelwise DeACs and the reduced size of the null space (11).
Some of the limitations of the proposed method are as follows.First, the segmentation method used is simplistic.Although the starting penalties used in conjunction ( 4) and ( 8) ultimately worked, these need to be optimized.The method could also be made more robust by making the DeAC estimation step less sensitive to the quality of the SOS map.Finally, computation time needs to be reduced either by more efficient coding or use of GPU.

VI. CONCLUSION
In this paper we describe a qPAT framework to estimate DeACs.A key advantage of our qPAT method is its accuracy throughout a large imaging target due to the combined use of full wave inversion with fluence decay compensation delivered by the MCX toolbox: the estimated DeAC map is within 12% of the actual absorption map on which the simulation was based.Our results demonstrate the potential use of qPAT for transcranial imaging of neonatal brains which could be developed, for example, to detect conditions such as hemorrhage and hypoxia.Our method suffers from expensive computations; we plan to overcome this to some extent by implementing it on GPU.

A. Standard Linear Solid (SLS) Attenuation Mechanism
Basically, m is a frequency-dependent non-linear mapping between SOS and attenuation factor at the reference frequency.[106] have compiled eight different models related to some empirical attenuation mechanisms.In this paper, we use SLS mechanism defined as where ω r is a reference frequency, ν is SOS, τ ε and τ σ are relaxation times related to the constants of the effective springs and dash-pot of the model [107] and defined as where α is the acoustic attenuation.Also, by having complex valued m, one can extract acoustic attenuation as Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

Fig. 1 .
Fig. 1.In photoacoustic imaging, a pulsed laser illuminates tissue, and acoustic pressure waves (PA waves) are generated by chromophores in the tissue.These PA waves propagate through the tissue and are measured by ultrasound (US) transducers.

Fig. 2 .
Fig. 2. Configuration of the simulated US/PA tomography system.US data collection proceeds as follows.Each of the 128 transducers is activated singly, with the other 127 transducers in receive mode.ECT: extra-cranial tissue, CSF: cerebrospinal fluid, GM: gray matter, WM: white matter.

Fig. 4 .
Fig. 4. SOS and acoustic attenuation estimation by IR-WRI.(a) True SOS map, (b) initial guess for the SOS map, (c) estimated SOS map, (d-f), same as (a-c) but for acoustic attenuation, (g) cross sections of the SOS maps in a, b, and c at the location of the dashed line, (h) error graph for the SOS map (blue) and acoustic attenuation map (red).

Fig. 4 (
a) shows the true SOS model.Figs.4(b) and 4(c) show the initial guess and estimated SOS maps after 185 iterations, respectively.Figs.4(d)-4(f) are the same as Figs.4(a)-4(c), but for acoustic attenuation.Cross-section profiles of the true, initial, and estimated SOS maps at the white dashed lines shown in Figs.4(a)-4(c) are plotted in Fig. 4(g).The SOS percent error (SOS Error) defined as 100 pi xel |true SOS -estimated SOS| |true SOS| ,

Fig. 5 (
b) shows the segmented model described in III-C.1.
Fig. 6(b) utilizes the SOS and acoustic attenuation maps derived previously using IR-WRI (Figs. 4(c) and 4(f)).The TR DeAC percentage error for the model shown in Fig. 6(b) is 71.98%

TABLE II OPTICAL
SCATTERING COEFFICIENT (µ s ) AND ANISOTROPY FACTOR (g) OF THE NEONATE'S HEAD AT 800 NM