Iterative image reconstruction in transcranial photoacoustic tomography based on the elastic wave equation

Photoacoustic computed tomography (PACT) is an emerging computed imaging modality that exploits optical contrast and ultrasonic detection principles to form images of the photoacoustically induced initial pressure distribution within tissue. The PACT reconstruction problem corresponds to a time-domain inverse source problem, where the initial pressure distribution is recovered from the measurements recorded on an aperture outside the support of the source. A major challenge in transcranial PACT of the brain is to compensate for aberrations and attenuation in the measured data due to the propagation of the photoacoustic wavefields through the skull. To properly account for these effects, a wave equation-based inversion method can be employed that can model the heterogeneous elastic properties of the medium. In this study, an optimization-based image reconstruction method for 3D transcranial PACT is developed based on the elastic wave equation. To accomplish this, a forward-adjoint operator pair based on a finite-difference time-domain discretization of the 3D elastic wave equation is utilized to compute penalized least squares estimates of the initial pressure distribution. Computer-simulation and experimental studies are conducted to investigate the robustness of the reconstruction method to model mismatch and its ability to effectively resolve cortical and superficial brain structures.


Abstract
Photoacoustic computed tomography (PACT) is an emerging computed imaging modality that exploits optical contrast and ultrasonic detection principles to form images of the photoacoustically induced initial pressure distribution within tissue. The PACT reconstruction problem corresponds to a time-domain inverse source problem, where the initial pressure distribution is recovered from the measurements recorded on an aperture outside the support of the source. A major challenge in transcranial PACT of the brain is to compensate for aberrations and attenuation in the measured data due to the propagation of the photoacoustic wavefields through the skull. To properly account for these effects, a wave equation-based inversion method can be employed that can model the heterogeneous elastic properties of the medium. In this study, an optimization-based image reconstruction method for 3D transcranial PACT is developed based on the elastic wave equation. To accomplish this, a forward-adjoint operator pair based on a finite-difference time-domain discretization of the 3D elastic wave equation is utilized to compute penalized least squares estimates of the initial pressure distribution. Computer-simulation and experimental studies are conducted to investigate the robustness of the reconstruction method to model mismatch and its ability to effectively resolve cortical and superficial brain structures. such as DOT, suffers from inherently low spatial resolution in transcranial imaging (Culver et al 2003). On the other hand, PACT provides the same functional contrasts as optical imaging but has a superior spatial resolution due to the principles of ultrasonic detection. Hence, PACT is well-suited for anatomical and functional brain imaging applications. As PACT exploits haemoglobin-based endogenous optical contrast, it can image anatomical structures such as blood vessels as well as measure functional parameters, such as oxygen saturation (sO 2 ), and the metabolic rate of oxygen consumption (MRO 2 ) (Zhang et al 2006, Yao et al 2011.
In vivo transcranial PACT studies have revealed structure and haemodynamic responses in small animals (Wang et al 2003, Li et al 2010, Xu et al 2011. Because the skulls in these studies were relatively thin (∼1 mm), they did not significantly aberrate the photoacoustic (PA) wavefields. Hence, conventional PACT image reconstruction algorithms such as backprojection (BP) (Finch andRakesh 2004, Kunyansky 2007), that assume that the entire medium is acoustically homogeneous were successfully employed. This assumption is violated in transcranial PACT applications involving primate skulls. As a result, images produced by conventional reconstruction methods in the presence of a primate skull can contain significant distortions and degraded spatial resolution (Xu and Wang 2006b, Jin et al 2008, Yang and Wang 2008, Nie et al 2011, Huang et al 2012. Thus, to render transcranial PACT an effective imaging modality for use in humans, it is necessary to develop image reconstruction methodologies that can accurately compensate for the skull-induced aberrations and attenuation of the recorded PA signals.
Numerous image reconstruction methods have been proposed to compensate for aberrations and attenuation of the measured PA wavefields induced by an acoustically heterogeneous fluid medium (Jin et al 2008, Huang et al 2012, 2013. However, a limitation of such works is that they assume a simplified wave propagation model in which longitudinal-to-shear-wave mode conversion within the skull (Fry 1978, White et al 2006, Schoonover and Anastasio 2011 was neglected. As a result of the simplified models employed, only modest improvements in image quality were observed as compared to the use of a standard BP-based reconstruction algorithm. Therefore, there remains an important need for the development of improved reconstruction methodologies for transcranial PACT that are based on more accurate models of the imaging physics. To circumvent limitations of the aforementioned approaches, a numerical framework for image reconstruction in transcranial PACT based on an elastic wave equation that describes an isotropic, lossy and heterogeneous medium was proposed and validated (Mitsuhashi et al 2017). In that work, a discrete forward operator (i.e. imaging model) and an associated adjoint operator for transcranial PACT were implemented based on the 3D elastic wave equation. The adjoint operator was validated and employed as a reconstruction operator to compensate for aberrations introduced by the skull. Even though the adjoint operator may perform better than existing filtered BP methods in compensating for the acoustic and elastic properties of the skull, the use of the adjoint operator as a reconstruction operator possesses limitations. For instance, for a lossy medium and/or a sparsely sampled measurement aperture, the adjoint operator is not a good approximation of the inverse operator.
Furthermore, the use of an adjoint operator as a reconstruction operator does not provide a general framework for incorporating regularization that can help mitigate the effects of measurement noise and incomplete measurement data. Hence, the use of the adjoint operator to directly reconstruct images can result in suboptimal contrast and suboptimal tradeoffs between image variance and spatial resolution. The use of the optimizationbased image reconstruction methods proposed in this work can help circumvent these shortcomings. In other works, the adjoint of the continuous wave operator, which describes the propagation of PA waves in linear isotropic viscoelastic media has been proposed (Javaherian and Holman 2018). Such studies have not explored image reconstruction problems with realistically sized 3D volumes or utilized experimentally acquired data.
In transcranial PACT, the head is illuminated from the outside with a near-infrared optical pulse resulting in the generation of acoustic pressure signals via the PA effect. As the scalp vessels can strongly absorb the illuminated light, only a small fraction of the illuminated light or optical fluence makes its way to the cortex resulting in the generation of initial pressure distribution in the cortex. Hence, there is a large difference in amplitude between the photoacoustically induced initial pressure distribution in the cortex and the photoacoustically induced initial pressure distribution in the scalp. Thus, spatially resolving the weak cortical structures from the strong scalp or superficial structures poses a significant challenge for image reconstruction . A challenge for any transcranial PACT reconstruction algorithm is to be able to accurately reconstruct the cortical structures in the presence of large differences in amplitude between the initial pressure distribution generated in the cortex and the initial pressure distribution generated in the scalp. Furthermore, in most transcranial PACT image reconstruction algorithms it is assumed that the acoustic properties of the skull are known accurately beforehand. However, such information may be difficult to obtain in practice. Any inaccuracies in modeling the acoustic properties of the skull will be a source of model mismatch error in conventional image reconstruction algorithms. Thus, an additional challenge for a transcranial PACT reconstruction algorithm is to be able to robustly model errors arising from inaccurate modeling of the acoustic characteristics of the skull.
Given the previously developed numerical framework for implementing the forward-adjoint operator pair (Mitsuhashi et al 2017), a natural topic of investigation is to explore the use of these operators in studies of optimization-based image reconstruction. Hence, in this work a modern optimization-based iterative image reconstruction algorithm for transcranial PACT is proposed and its performance in overcoming the aforementioned challenges is evaluated through computer-simulation and experimental studies. The computer-simulation studies explore the robustness of the optimization-based image reconstruction algorithm to noise as well as the model mismatch errors. In the presence of large differences in amplitude between the initial pressure distribution generated in the scalp and the initial pressure distribution generated in the cortex, the performance of the optimization-based reconstruction method (OBRM) in accurately spatially resolving the cortical and superficial structures is also investigated through computer-simulation and experimental studies.
The article is organized as follows. In section 2, the physics of elastic wave propagation in elastic media is reviewed along with the modern optimization-based formulation of the image reconstruction problem. A description of the computer-simulation studies and a discussion of the corresponding results is provided in section 3. In section 4, the experimental phantom studies are described and the article concludes with a summary in section 5.

Photoacoustic wavefield propagation: continuous and discrete formulation
Let the photoacoustically-induced stress tensor at location r ∈ R 3 and time t 0 be defined as where σ ij (r, t) represents the stress in the ith direction acting on a plane perpendicular to the j th direction. Additionally, let p 0 (r) denote the photoacoustically-induced initial pressure distribution within the object, referred to as the object function, and let u(r, t) ≡ (u 1 (r, t),u 2 (r, t),u 3 (r, t)) represent the vector-valued acoustic particle velocity. Here, ρ(r), λ(r) , and µ(r) represent the spatial distribution of the medium's density and the spatial distribution of the Lamé parameters, respectively. The compressional and shear wave propagation speed are given by c l (r) = λ(r)+2µ(r) ρ (r) and c s (r) = µ(r) ρ(r) , respectively. The object function p 0 (r) and all functions that describe the elastic isotropic medium are assumed to be bounded with compact supports. Here, we model acoustic absorption of the skull by a diffusive absorption model (Moczo et al 2007). The diffusive absorption model assumes that the wavefield absorption is independent of temporal frequency; this model does not accurately describe the power law based frequency-dependent absorption characteristics of soft tissue and bone (Szabo 1994, Treeby andCox 2014). This model, however, is sufficiently accurate in cases where the bandwidth of the PA signals is limited. Moreover, the model also assumes that the shear absorption to compressional absorption ratio is given by the compressional velocity to shear velocity ratio. This is approximately true in bone, because slower shear waves are, in fact, more attenuated than faster compressive waves (Pinton et al 2011). The model also assumes that for fluid media (soft tissue), the diffusive absorption value only describes the absorption of compressional waves as shear waves are not supported in fluid media.
In a heterogenous isotropic elastic medium with a diffusive acoustic absorption distribution α(r), the propagation of u(r, t) and σ(r, t) can be modeled by the following two coupled equations (Alterman and Karal 1968, Boore 1972, Virieux 1986, Madariaga et al 1998: subject to the initial conditions Here, tr (·) is the operator that calculates the trace of a matrix and I ∈ R 3×3 is the identity matrix. In equation (2c), the object function p 0 (r) is assumed to be compactly supported in a fluid medium where the shear modulus µ(r) = 0. Hence, the formulation of the initial value condition assumes that there are no optical absorbers inside the elastic solid. Consider that p(r 0 , t) = tr(σ(r 0 , t)) is recorded outside the support of the object for r 0 ∈ ∂Ω and t ∈ [0, T], where ∂Ω ⊂ R 3 denote a measurement aperture. In this case, the imaging model can be described by a continuous-to-continuous (C-C) mapping as: where r ∈ Ω and H : is a linear operator that denotes the action of the wave equation given in equation (2). Moreover, p(r 0 , t) ∈ L 2 (∂Ω × [0, T]) denotes the measured data function, and the operator M is the restriction In practice, the detected pressure wavefield is discretized temporally and spatially at specific transducer locations. Let p ∈ R NrL denote the discretized pressure signal, where N r represents the number of transducers and L represents the total number of discrete temporal samples. In the imaging model described in equation (3), the effects of the acousto-electrical impulse response (EIR) as well as the spatial impulse response (SIR) of the transducer have been ignored. However, the transducer responses can be incorporated into the imaging model readily (Huang et al 2013). A continuous-to-discrete (C-D) PACT imaging model, can be generally expressed as Here, ∆ t is the temporal sampling interval and the vectors r k 0 ∈ R 3 , k = 0, 1, 2, . . . , N r − 1, represent the position vectors of the N r receivers on the aperture ∂Ω.
To obtain a discrete-to-discrete (D-D) imaging model for use with an OBRM, finite-dimensional approximate representations of the object function p 0 (r) and the material parameters ρ(r), α(r), λ(r), and µ(r) need to be introduced. In this study, the numerical method employed to solve the wave equation in equation (2) was a 10th-order staggered grid finite difference time domain (FDTD) scheme. Note that for the staggered finitedifference (FD) scheme, the material properties, stress, and particle velocity functions are sampled at different points of a staggered FDTD cell as shown in figure 1. Let ρ, α, λ, µ, and p 0 ∈ R N be the finite dimensional representations of ρ(r), α(r), λ(r), µ(r) and p 0 (r), respectively, where N is the total number of grid points on the 3D grid. For a given ρ, α, λ, µ, and p 0 , the FDTD method for solving the elastic wave equation can be described in operator form as where H ∈ R NL×N is the discrete approximation of the wave operator H that solves the initial value problem defined in equation (2) and M ∈ R NrL×NL is a sampling matrix that maps pressure data sampled on the computational grid onto the transducer locations. In this study, the transducers are assumed to be point-like and thus when the receiver and grid point locations do not coincide, an interpolation method is required. In this study, the elements of M are chosen such that trilinear interpolation is performed. The goal of image reconstruction in a discrete setting is to determine an estimate of p 0 given the measured data p and the forward model in equation (5). In our previous work, the use of the adjoint operator as a reconstruction operator was investigated (Mitsuhashi et al 2017). The reconstructed estimate of p 0 obtained by use of the adjoint operator is given bŷ

Optimization-based image reconstruction method (OBRM)
By use of the proposed D-D imaging model defined in equation (5), a variety of OBRMs can be employed for determining estimates of p 0 . In this work, we utilize an iterative reconstruction algorithm that seeks to compute penalized least squares estimates (PLS) by solving where γ is a regularization parameter and R(p 0 ) is a regularizing penalty term. In the studies below, the total variance (TV) semi-norm penalty, given by was employed. Here, [p 0 ] n denotes the nth grid node, and are the neighboring nodes before the nth node along the first, second and third dimension, respectively. The fast iterative shrinkage/thresholding algorithm (FISTA) with backtracking linesearch and adaptive restart was employed to solve the optimization problem in equation (7)

Computer-simulation studies
The performance and robustness of the proposed OBRM was initially evaluated by use of computer-simulation studies as described below. Additionally, computer-simulation studies that assess the performance of the OBRM in the presence of measurement noise and the effects of TV regularization are summarized in appendix A.

Resolution of superficial and cortical structures
As a consequence of the large differences in amplitude between the initial pressure distribution generated in the cortex and the initial pressure distribution generated in the scalp, spatially resolving the weak cortical structures from the strong scalp or superficial structures poses a significant challenge for transcranial image reconstruction . A computer-simulation study was performed to evaluate the performance of the proposed OBRM in accurately reconstructing cortical structures in the presence of strong absorbers in the scalp.

Methods
To study the performance of the OBRM in resolving scalp and cortical structures, the numerical blood vessel phantom shown in figure 2 was employed. The superficial blood vessels located in the scalp and the cortical vessels are visible in the maximum intensity projection (MIP) images. For display purposes, the dynamic range was adjusted so that the cortical vessels are visible; as a result, the superficial vessels appear saturated. The optical fluence ratio between the superficial and cortical vessel was assumed to be 50 . Hence, the ratio of initial pressure distribution between superficial and cortical blood vessels in the phantom shown in figure 2 was set to 50. A 3D computational volume of 408.0 mm × 408.0 mm × 190.2 mm was employed.
A 3D approximately hemispherical measurement system, as shown in figure 3(a), was emulated. The measurement system consisted of a total 38 400 transducers distributed evenly across 64 rings. The measurement system was similar to the experimental system used for acquisition of pressure data as described in section 4. The position of the measurement system with respect to the skull is shown in figure 3(b).
A 3D x-ray CT image of an intact human skull was utilized to generate the isotropic, elastic medium employed for the simulation studies. The CT images of the skull were employed to infer the thickness and contour of the skull and place it within the 3D simulation volume. For this set of simulation studies, the skull was assumed to be an acoustically homogeneous. In reality, however the skull bone is acoustically heterogeneous and consists of three relatively homogeneous layers: the outer table, inner table and the central diploe layer (Sadleir andArgibay 2007, Marieb 1997).
In order to extract the contour and location of the skull from CT images, a segmentation algorithm was employed. The segmentation algorithm generated a binary mask specifying the location of the skull within the 3D volume. The value of the medium parameters in the 3D volume were assigned such that the skull acoustic parameters (ρ = 1850 kg m −3 , c l = 3.0 mm µs −1 , c s = 1.48 mm µs −1 and α = 0.75 µs −1 ) were set at all grid positions where the mask was equal to one and the background acoustic parameters (ρ = 1000 kg m −3 , c l = 1.5 mm µs −1 , c s = 0.0 mm µs −1 and α = 0.0 µs −1 ) were set all grid positions where the mask was equal to zero. At the material interface between the skull and the background fluid medium, the density and the absorption values were averaged to avoid any instability with the FDTD wave equation solver (Moczo et al 2007).
Different discretization strategies were employed to generate the simulated pressure data and to reconstruct the initial pressure distribution, thereby avoiding inverse crime (Colton and Kress 2013). The simulated pressure data were generated by use of a spatial grid size of ∆x = 0.225 mm and a temporal sampling rate of 50 MHz. The simulated pressure data were corrupted white Gaussian noise with a standard deviation 10 percent of the maximum amplitude of cortical signals.

Results
The MIP images of the 3D volumes reconstructed by use of the adjoint operator and the OBRM are shown in figures 4(a)-(f), respectively. Each column of images represents the MIPs of the reconstructed initial pressure distribution along the x-, y -and z-axes, respectively. All images were reconstructed with a uniform spatial grid sampling of ∆x = 0.3 mm. From figures 4(a)-(c), it is observed that the images produced by use of the adjoint operator do not clearly reveal the cortical vessels. This is because the bleed down effect from the superficial structures overwhelms the signal from the cortical structures. However this is not the case in the images reconstructed by use of the OBRM, as shown in figures 4(d)-(f). In these images, the superficial vessels are accurately resolved because the strong signals from the superficial vessels are spatially separated and do not overwhelm the cortical vessels. The OBRM was successful in improving the resolution of the superficial structures, thereby reducing the amount of contamination that affects the cortical structures. Hence, compared to the images reconstructed by use of the adjoint operator, the cortical vessels are prominent and visible in the images reconstructed by use of the OBRM.

Model mismatch studies
The proposed OBRM requires that the acoustic properties of the skull are known. However, such information may be difficult to obtain in practice. Any inaccuracies in modeling the acoustic properties of the skull will be a source of model mismatch error in equation (5). Model mismatch errors can result in artifacts in the reconstructed images. Hence, in this section, computer-simulation studies were conducted to evaluate the robustness of the proposed OBRM to model mismatch.
There are at least two major types of skull model errors/mismatches that are challenging to account for in OBRM. The first source of model mismatch resides in the inaccuracies in modeling the frequency-dependent absorption characteristics of the skull (Fry 1978, White et al 2006. In transcranial PACT imaging applications, the acoustic absorption is not negligible. In the initial value problem described in equation (2), the acoustic absorption within the skull was described by a diffusive absorption model (Aubry et al 2003). The model ignores the fact that the wavefield absorption within the skull is dependent on the temporal frequency of the PA signals (White et al 2006). This approximation, however, may be only reasonable for cases where the bandwidth of the PA signals are limited (Aubry et al 2003). Computer-simulation studies were performed to explore the feasibility of the proposed OBRM when there was significant model mismatch with respect to accurately modeling the acoustic absorption characteristics of the skull. The second major source of model mismatch that we commonly encounter originates from the failure to model the presence of acoustic heterogeneities within the skull. Various adjunct imaging data based methods as well as joint reconstruction methods can be employed to account for acoustic variations in the skull in the proposed optimization-based image reconstruction framework (Huang et al 2016. A study of these techniques is beyond the scope of this study.

Methods
To study the impact of model mismatch associated with acoustic absorption, the cortical blood vessel phantom shown in figure 5 was employed. The forward pressure data were generated employing a skull model that had linear frequency-dependent absorption characteristics. The linear dependence of the compressional wave absorption value with respect to frequency is shown in figure 6(a). The slope of the linear frequency-dependent was 16 dBcm −1 MHz −1 . The shear wave absorption value is dependent on the ratio between the compressional speed and shear speed. Hence, for the configuration described in figure 6(a), the shear wave absorption characteristics will be linear with respect to frequency and will have a slope of 32 dBcm −1 MHz −1 .
In order to generate the simulated pressure data that modeled the frequency-dependent absorption characteristics of the skull, multiple runs of the imaging model with different diffusive absorption values had to be performed. This is due to the fact that the imaging operator in equation (5) assumes a diffusive absorption model where the absorption coefficient value is independent of frequency of PA signals. In order to account for the frequency-dependent absorption characteristics of the skull, multiple runs of the imaging operator and filtering operations need to be performed. A schematic describing the generation of the forward pressure data from a skull model with frequency-dependent absorption characteristics is shown in figure 6(b). The set of discrete diffusive absorption coefficient values {α i } N bin −1 i=0 used to generate the forward data are given by α i = α 0 f i , where α 0 is the linear slope value of the frequency-dependent absorption model. Here, N bin are the number of frequency bins and f i represents the central frequency of the each bin. Let ∆f be the size of each frequency bin. To generate the pressure data, the frequency spectrum was discretized into different bins and the corresponding diffusive absorption value is assigned based on the central frequency of the bin. Hence, each frequency band was assigned its own constant diffusive absorption coefficient value that was calculated based on the frequency-dependent absorption characteristics shown in figure 6(a). Subsequently, pressure data corresponding to each of the diffusive absorption value α i were generated using the imaging model H(α i ). Here, H(α i ) represents the D-D imaging model when the diffusive acoustic absorption value of the skull was set to α i . The pressure data corresponding to each f i are then subsequently filtered and summed to generate the final pressure data. Hence, the final pressure data contains contribution from all the frequency bands that have been appropriately weighted based on the action of the imaging operator H(α i ) and the filtering operation. For the simulation studies, we assume N bin = 10 and ∆f = 0.2 MHz. Figure 4. The MIP of the reconstructed 3D initial pressure distribution. The first row of images corresponds to images reconstructed using the matched adjoint operator. The second row of images corresponds to images reconstructed using the OBRM. Each column of images corresponds to the MIP of the 3D initial pressure distribution along x-, y -and z-axes, respectively.
The same 3D computational volume and the measurement geometry employed in section 3.1 were utilized. Furthermore, different discretization strategies were employed to generate the simulated pressure data and to reconstruct the initial pressure distribution. The simulated pressure data were generated using a spatial grid size of ∆x = 0.225 mm and a temporal sampling rate of 50 MHz. The simulated pressure data were corrupted with white Gaussian noise with a standard deviation 5 percent of the maximum amplitude of cortical signals before the application of the reconstruction method. The OBRM employed a diffusive absorption model with a fixed skull absorption coefficient for image reconstruction. All other acoustic properties of the skull other than the absorption coefficient value were assumed to be known.

Results
The MIPs of the 3D images reconstructed by use of the OBRM are shown in figure 7. Each column of images represent the maximum intensity projections of the reconstructed initial pressure distribution along the x-, yand z-axes, respectively. Each row of figure 7 displays images reconstructed with a different diffusive absorption value. All the reconstruction was performed with a spatial grid sampling of ∆x = 0.3 mm.
From figure 7, one can observe that the minimum root mean square error (RMSE) was obtained for the initial pressure distribution reconstructed with a diffusive absorption value of α = 0.75 µs −1 . The optimal value of the diffusive absorption value of α = 0.75 is associated with a central frequency of 0.5 MHz and the frequency bin of [0.4 MHz, 0.6 MHz]. The optimal diffusive absorption value is defined as the diffusive absorption value associated with the OBRM that produced images with the lowest RMSE. The optimal absorption value was associated with the frequency band that contains the most energy in the measured pressure data. From the results, one can observe that the model mismatch due to inaccuracy in modeling for frequency-dependent characteristics of the skull can be accounted for by using a diffusive absorption model, as long as the diffusive absorption value is representative of the dominant frequency band in the measured pressure data. The line profiles through the reconstructed 3D images, shown in figure 9(a), also emphasize these observations. The MIPs of the reconstructed 3D images reconstructed by use of the adjoint operator are shown in figures 8(a)-(c). The adjoint operation was performed with the optimal diffusive absorption value of α = 0.75 µs −1 . From figures 7 and 8(a)-(c), the contrast of the images obtained using the adjoint as the reconstruction operator is observed to be poor compared to the images obtained using the OBRM. The sharp difference in contrast can be also be observed by comparing the line profiles through the reconstructed images shown in figure 9(b). In summary, one can observe that the impact of the failure to model the frequency-dependent absorption characteristics of the skull results in blurring of the reconstructed images. However, the impact of model mismatch errors on the accuracy of the reconstruction can be minimized by employing diffusive absorption value associated with the dominant frequency band in the measured pressure data. This conclusion was confirmed by calculating at the RMSE of the reconstructed pressure distribution with different diffusive absorption values.

Experimental studies
In this section, studies that utilized experimental PACT data produced by use of a physical phantom are presented. Two sets of experimental PACT data were employed to validate the proposed OBRM. The experimental phantoms and the setup employed to acquire the measurement data are described below.

Phantoms
The design of the physical phantom employed for the experimental studies was motivated by transcranial PACT. The phantom was comprised of a spherical acrylic globe of thickness of 2.5 mm and an inner radius of 76.2 mm, placed within a 3D volume filled with water. The acrylic globe is an elastic solid that possesses acoustic property values that are representative of a human skull. For the experimental studies, two configurations of optical absorbing vessel-like structures were employed with the acrylic globe. In the first configuration, referred to as phantom #1, optically absorbing vessel-like structures were painted on the inner surface of the acrylic globe, as shown in figure 10. These vessels were intended to mimic cortical vessels that reside near the top surface of a brain. Some of the cortical vessels painted on the acrylic globe have been labeled in figure 10 and will be referred to in the subsequent discussions. In the second configuration, also referred to as phantom #2, strips of superficial absorbers were placed outside the acrylic globe along with cortical vessels of phantom #1 as shown in figure 11. Hence, the experimental data obtained from illuminating phantom #2 mimicked the more realistic scenario described in section 3.1 whereby the PACT data were generated from cortical as well as superficial structures.

System
A 64 element transducer arc was scanned over 600 evenly distributed locations in the azimuthal direction in an approximately hemispherical configuration as shown in figure 12 to acquire the experimental data. The location of the acquisition system with respect to the acrylic globe is shown in figure 12(a). A short 6 ns laser pulse (Nd-YAG Quantel Brilliant B laser) with a repetition rate of 10 Hz at a wavelength of 1064 nm was used to irradiate a sample located in the center of the measurement system as shown in figure 12(a). The generated acoustic signals were detected by unfocused transducers, with a center frequency of 1 MHz and a −6 dB fractional bandwidth of 80%. The electrical signals recorded by the transducers were sampled at a temporal sampling rate of 25 MHz. This system was employed to image the physical phantoms.

Data preprocessing
Prior to image reconstruction, the measured data were preprocessed. The preprocessing involved deconvolving the acquired EIR of the transducer. The Wiener deconvolution method was employed to estimate the deconvolved PA signals from the raw electrical transducer measurements. After deconvolution, the data were filtered with a Hann-window low-pass filter with a cutoff frequency of 2 MHz. The filtered data were also upsampled by a factor of 2.5, with the goal of circumventing numerical stability issues with the wave equation solver.

Results
When employing either the adjoint operator or the OBRM for image reconstruction, the acoustic parameters of the acrylic globe were set to be ρ = 1200 kg m −3 , c l = 2.8 mm µs −1 , c s = 1.4 mm µs −1 , and α = 0.1 µs −1 . Similar to the computer-simulation study, at the material interface between the globe and the background fluid medium, the density and the absorption values were arithmetically averaged to avoid any instability issues with the FDTD wave equation solver. The images reconstructed from experimental data are shown in figures 13 and 14. The first set of results shown in figure 13 corresponds to an experiment where the pressure data were generated by illuminating phantom #1. Each column of images shown in figure 13, corresponds to 2D maximum intensity projections of the reconstructed 3D volume along three mutually perpendicular directions. The first row of   figure 13. The contrast was evaluated by dividing the difference in mean signal intensity of the vessel structure and the background in the selected 3D volume and dividing it by the variance of the background signal. It was found that the contrast of the image reconstructed by use of the the OBRM was 12.6 compared to a contrast of 3.2 for the image reconstructed by use of the adjoint operator. Hence, an improvement of contrast by a factor of approximately 4 was observed when the image was reconstructed by use of the OBRM as opposed to reconstructed by use of the matched adjoint operator. Furthermore, one can observe vessel like structures labelled v1, v2, v3 towards the outer periphery of the images are distinctly visible in the images reconstructed by use of the OBRM while the same structures are not in the images reconstructed by use of the adjoint operator. The results in figure 13 demonstrate that the use of OBRM to reconstruct initial pressure distribution within an elastic media can lead to significant improvement in image quality compared to the images generated using the adjoint operator. The second set of experimental results corresponding to pressure data generated by illuminating phantom #2 is shown in figure 14. Each column of images shown in figure 14 have been generated by performing MIP of the reconstructed 3D volume along three mutually perpendicular directions. In addition, the first row of images shown in figure 14 corresponds to images reconstructed by use of the adjoint operator while the second row of images corresponds to images reconstructed by use of the OBRM with no TV regularization. In the images reconstructed by use of the adjoint operator as shown in figures 14(a)-(c), the signals from the superficial structures spread into the cortical region to considerably degrade the image quality of the reconstructed cortical structures. In the computer-simulation results obtained using the adjoint as a reconstruction operator, the cortical structures were not visible as they were completely overwhelmed by signals from the superficial structures. In the images reconstructed by use of the adjoint operator from the experimental data, the signal from the superficial structures do not overwhelm the cortical structures to the same extent as shown in the images reconstructed with the use of the adjoint operator in the computer-simulation studies. This is due to the fact that the acrylic globe is far less absorbing than a human skull. Hence, the fluence mismatch in the experimental setup using a acrylic skull is not high enough for the signals from the superficial structures to completely overwhelm the cortical vessels in the images reconstructed by use of an adjoint operator.
However, the signals from the superficial structures that contribute to obscuring the cortical structures, is significantly mitigated in the images reconstructed with the aid of the OBRM, as shown in figures 14(d)-(f). The vessel like cortical structures are much more clearly visible in the images reconstructed by use of the OBRM as compared to the images reconstructed by use of the adjoint operator. Hence, the results in figure 14 demonstrate that the OBRM is effective in accurately resolving the superficial structures and is successful in mitigating the signals from superficial structures that contribute to obscuring the cortical structures. This leads to significant improvements in image quality of the reconstructed cortical structures when employing the OBRM as compared to the adjoint operator.

Conclusion
In this work, an OBRM for use in transcranial PACT was proposed and investigated. The proposed method was based on the 3D elastic wave equation that accurately modeled the physics of acoustic wave propagation in linear, isotropic elastic lossy, heterogeneous media. Such modern reconstruction methods allow for accurate modeling of the imaging physics and provide a general framework to incorporate regularization, which can help mitigate the effects data incompleteness, model mismatch and noise. The proposed OBRM was validated and investigated in computer-simulation and experimental phantom studies whose designs were motivated by transcranial PACT applications. The computer-simulation and experimental studies compared and contrasted the performance of the OBRM and the adjoint as a reconstruction operator. The results from the studies demonstrated that the use of adjoint as a reconstruction operator yielded reconstructed images that possessed suboptimal contrast. The images reconstructed by use of the OBRM had improved contrast and also allowed for flexibility in reconstructing images with improved tradeoff between image variance and spatial resolution. Furthermore, the robustness of the OBRM to model mismatch was also studied. The model mismatch associated with failure to model accurately the acoustic absorption capabilities of the skull was explored through computer-simulation studies. The results demonstrated that when employing the OBRM, the failure to model the frequency-dependent characteristics did not have a significant effect on the reconstructed image quality when an appropriate diffusive absorption coefficient value was employed, that corresponded to the dominant frequency band in the measured pressure data.
A challenge for any transcranial PACT reconstruction algorithm is to be able to accurately reconstruct the cortical structures in the presence of large differences in amplitude between the initial pressure distribution generated in the cortex and the initial pressure distribution generated in the scalp. Hence, computer-simulation studies and experimental phantom studies were conducted to explore the feasibility of the proposed OBRM in resolving superficial and cortical structures. The results from both computer-simulation and experimental phantom studies demonstrated that the OBRM is effective in accurately resolving the superficial structures and cortical structures. The significant improvement in the image quality of the cortical structures when employing the OBRM was more pronounced when the difference in optical fluence between the scalp and the cortex was high. There remain several important topics for further investigation. One of the major sources of model mismatch that we commonly encounter in transcranial PACT originates from the failure to model the presence of acoustic heterogeneities within the skull. Therefore, future studies should work towards formulation of joint reconstruction problem for transcranial PACT applications. In such joint reconstruction formulations, the optimizationbased reconstruction procedure can be employed to reconstruct both the initial pressure distribution and the spatial distribution of the skull acoustic parameters concurrently.
For brevity, the following notation is employed The algorithm described in algorithm 1 has a number of advantages over alternative first-order optimization methods. First-order optimization methods refer to the set of methods for which the computation of the descent direction involves a linear or first order taylor approximation of the cost function. Hence, information about the gradient of the function is employed to compute the descent direction. The FISTA method belongs to a subset of the first-order methods called the proximal-gradient methods that permits the use of non-smooth regularization terms, such as the TV semi-norm. Furthermore, it also belongs to a family of optimization methods that for weakly convex optimization problems achieves the optimal asymptotic convergence rate (Beck and Teboulle 2009a, 2009b, Beck 2014. Namely, for a weakly convex function C(p 0 ), FISTA obtains the convergence rate where k is the iteration number and p (k) 0 is the estimate of the sought-after quantity for the kth iteration. The problem given by equation (7) represents a weakly convex optimization problem whenever H is not full rank and R(θ) is convex.