Reconstruction of Photoacoustic Tomography Inside a Scattering Layer Using a Matrix Filtering Method

: Photoacoustic (PA) tomography (PAT) has potential for use in brain imaging due to its rich optical contrast, high acoustic resolution in deep tissue, and good biosafety. However, the skull often poses challenges for transcranial brain imaging. The skull can cause severe distortion and attenuation of the phase and amplitude of PA waves, which leads to poor resolution, low contrast, and strong noise in the images. In this study, we propose an image reconstruction method to recover the PA image insider a skull-like scattering layer. This method reduces the scattering artifacts by combining a correlation matrix ﬁlter and a time reversal operator. Both numerical simulations and PA imaging experiments demonstrate that the proposed approach e ﬀ ectively improves the image quality with less speckle noise and better signal-to-noise ratio. The proposed method may improve the quality of PAT in a complex acoustic scattering environment, such as transcranial brain imaging.


Introduction
Photoacoustic (PA) tomography (PAT) is a noninvasive and nonionizing imaging technique that combines rich optical contrast and high ultrasonic resolution in deep tissue [1][2][3][4]. During the PAT imaging process, the biological tissue is irradiated by a laser pulse, and then optical absorbers in the tissue absorb electromagnetic energy and generate ultrasound due to rapid thermoelastic expansion. The ultrasound signal is then detected by the ultrasound transducer array arranged outside the tissue. Finally, images with optical absorption contrasts are extracted from the detected PA signal using a reconstruction algorithm. PAT breaks through the optical diffraction limit and obtains images of tissue at depths of 5-6 cm with good acoustic spatial resolution. Therefore, PAT has application potential for high-quality in vivo vascular structure imaging [5], investigations of breast tumor tissue [6][7][8][9], vasculature visualization [10,11], small animal whole-body imaging [12,13], osteoarthritis assessment [14,15], drug delivery monitoring [16], and hemodynamic functional imaging [17].
One promising application of PAT is high-resolution in vivo brain imaging, which is important in neurophysiology, neuropathology, and neurotherapy. In comparison with the often-used brain imaging modalities, including X-ray computerized tomography (CT) and magnetic resonance imaging (MRI), PAT is inexpensive and has good spatial resolution and temporal resolution [1,2]. X-ray CT involves exposure to ionizing radiation, but PAT is nonionizing and safer. PAT can provide structural images as well as functional images of the brain.
The classical image reconstruction algorithms used in PAT are usually based on the principle of coherent beamforming, which relies on an assumption of homogeneous media. In a scene with strong scattering, the coherence in signals is destroyed by the randomness of scattering, leading to a markedly 2 of 8 decreased penetration depth and significant degradation in image quality [18]. This limitation hinders the application of PAT in a wider range of scenarios, such as transcranial PA brain imaging within the skull bone. The skull bone has much higher acoustic impedance than soft tissue. A severe acoustic impedance mismatch could cause amplitude and phase distortion, and decrease the resolution and contrast of the transcranial imaging, significantly reducing the brain image quality [19][20][21][22]. Traditional approaches to overcome the limitation of acoustic scattering have mostly relied on a priori knowledge of the properties of tissue inhomogeneity or need to involve additional acoustic measurements, evidenced in the statistical reconstruction method [23], coherence factor optimization [24], and the interferometry method [25]. For a more general situation, a universal reconstruction algorithm with better performance is still required for acquiring high-quality images in scattering media.
In this study, we created an image reconstruction method that combines a correlation matrix filter and a time reversal operator using an annular transducer array. After dividing the circular array into multiple sections, each sector array can be approximated as a linear array via phase and amplitude compensation. Then, using the correlation characteristics of the direct waves, a correlation full-matrix filter can be applied to separate the direct waves from the scattering waves [26][27][28][29]. Finally, the imaging results of each section based on the time reversal method are superimposed to further improve the image quality. Both numerical simulation and phantom experiments were employed to validate the proposed method and examine its potential for extracting PA images inside a skull-like scattering layer.

Methods
As shown in Figure 1, an imaging target is placed in the region of interest (ROI), which is surrounded by a skull-like scattering layer. Under the illumination of the pulse laser, the optical absorber will radiate ultrasound waves to the surrounding medium as a result of the PA effect. After the acoustic waves propagate through the random scattering layer, a passive ultrasound annular array with 576 elements detects the ultrasound signal. The annular array is then divided into N S sectors for data processing and imaging. Each sector has N = 45 elements and can be approximated as a linear array so that the filtered imaging method proposed in our previous studies can be applied for the image reconstruction [26,27]. For each sector array, the recorded N-channel signals are arranged in a vector H(t) = [H 1 (t),H 2 (t), . . . ,H n (t), . . . ,H N (t)] where n = 1,2, . . . ,N. markedly decreased penetration depth and significant degradation in image quality [18]. This limitation hinders the application of PAT in a wider range of scenarios, such as transcranial PA brain imaging within the skull bone. The skull bone has much higher acoustic impedance than soft tissue. A severe acoustic impedance mismatch could cause amplitude and phase distortion, and decrease the resolution and contrast of the transcranial imaging, significantly reducing the brain image quality [19][20][21][22]. Traditional approaches to overcome the limitation of acoustic scattering have mostly relied on a priori knowledge of the properties of tissue inhomogeneity or need to involve additional acoustic measurements, evidenced in the statistical reconstruction method [23], coherence factor optimization [24], and the interferometry method [25]. For a more general situation, a universal reconstruction algorithm with better performance is still required for acquiring high-quality images in scattering media.
In this study, we created an image reconstruction method that combines a correlation matrix filter and a time reversal operator using an annular transducer array. After dividing the circular array into multiple sections, each sector array can be approximated as a linear array via phase and amplitude compensation. Then, using the correlation characteristics of the direct waves, a correlation full-matrix filter can be applied to separate the direct waves from the scattering waves [26][27][28][29]. Finally, the imaging results of each section based on the time reversal method are superimposed to further improve the image quality. Both numerical simulation and phantom experiments were employed to validate the proposed method and examine its potential for extracting PA images inside a skull-like scattering layer.

Methods
As shown in Figure 1, an imaging target is placed in the region of interest (ROI), which is surrounded by a skull-like scattering layer. Under the illumination of the pulse laser, the optical absorber will radiate ultrasound waves to the surrounding medium as a result of the PA effect. After the acoustic waves propagate through the random scattering layer, a passive ultrasound annular array with 576 elements detects the ultrasound signal. The annular array is then divided into NS sectors for data processing and imaging. Each sector has N = 45 elements and can be approximated as a linear array so that the filtered imaging method proposed in our previous studies can be applied for the image reconstruction [26,27].    The sector array can be approximated to a linear array by calculating the corresponding delay. For any imaging target, its propagation distances to the nth transducers of the two arrays are l 1 and l 2 , respectively. The corresponding delay is d n = (l 1 − l 2 )F s /c, where c is the sound speed and F s is the sampling frequency. The time domain signals for the linear array are then The signals can be written in their frequency domain by applying a Fourier transform to P(T,f ) = [P 1 (T,f ),P 2 (T,f ), . . . ,P n (T,f ), . . . ,P N (T,f )], where P n (T,f ) is the short-time Fourier transform of the nth channel signal in the time window [T − ∆t/2, T + ∆t/2], T is the time of flight, and ∆t is the width of the window. P(T,f ) contains both the direct wave P D (T,f ) and the scattering wave P S (T,f ). When the concentration of the scatterers is sufficiently high, the scattering wave occupies a dominant position and the imaging quality is greatly degraded due to the randomness of the scattering process.
To reduce the influence of the randomly scattering, a matrix filtering method can be applied to extract the direct wave from the received signals of a linear array [26]. Firstly, a response matrix can be constructed as: Depending on whether or not contain the term P S , the matrix K can be divided into two parts: K R and K C , where K R is a random term due to the random nature of the scatterer distribution and the scattering paths and K C is independent of the scattering paths and shows a deterministic relation of the phase between its elements along the antidiagonal: where k is the wave number, w is the pitch size, and Z is the depth of the time window. The phase relation along the antidiagonal only depends on the channel position qw, indicating that the direct waves from the optical absorbers manifest themselves as a particular coherence on the antidiagonals of the matrix. Using the coherent property, we can extract the direct wave by the following steps. Firstly, the original matrix K is rotated 45 • counterclockwise so that the coherent property occurs on the columns. The elements are divided into two parts based on their symmetry characteristic and are expanded to two new matrices: A 1 and A 2 [26]. The matrices A 1 and A 2 are called indifferently A as they are filtered in the same way. A filtering matrix F = CC † is then applied to extract the coherent part from matrix A as: is the characteristic spacing of direct waves. The coherent part remains unchanged during the filtering process, as FA C = A C and the random part decreases by a factor (N + 1)/2. Finally, a filtered matrix K F with the original coordinate can be obtained by applying the inverse process of the first step to matrices A 1 F and A 2 F .
Using a time reversal operator, we can recover the image of the optical absorbers from the matrix K F corresponding to the single sector area [28,29]. The singular value decomposition is applied to the matrix K F = U F Λ F V F † , and each optical absorber is mainly associated with one non-zero singular λ i F contained in diagonal matrix Λ [30]. Therefore, by numerically back-propagating the singular vector, the image of the ith singular valued at a depth of Z = cT can be achieved by the time reversal operator I i = λ i F |G × V|, where G with the component g m = exp( jkr m )/ √ r m is the propagating operator between the array and the focal plane in a homogeneous media, and r m is the distance between the mth element in the array and the point in the focal plane. By superimposing the imaging results of multiple sectors, we can further improve the imaging quality and reduce the influence of the inhomogeneous aberrating effect of the scattering layer. For the wanted imaging point A, we calculated the image value I m A for all the linear arrays with m = 1, 2, . . . , N s . Taking the first sector area as an example, the image value for points corresponding to the array can be obtained by applying the time reversal operator. If point A is within the imaging area, we can find two imaging points, B and C, closest to point A in the focal plane and then perform a weighted evaluation based on the distance from point A to points B and C: where d AC and d AB represent the distances between A, B and A, C, respectively. If point A is outside the imaging area, the imaging value is then recorded as 0. When we calculate the imaging value of point A for all sector areas, its final value is the result of taking the summation and average to all the non-zero values: where N Z is the number of zero values. As such, we can fully use the valid information received by all the transducers of the annular array for image reconstruction of each point in the ROI.

Numerical Simulation
Numerical simulations were used to validate the proposed method. The simulation is depicted in Figure 2. The radius of the annular array was set to R a = 40 mm, and there were N = 576 elements in the array. The acoustic scattering layer consisted of 120 acoustic scatterers with a diameter of d = 0.8 mm. The inner and outside radius of the scattering layer were R 1 = 10 mm and R 2 = 15 mm, respectively. The sound velocity and the density of these scatterers were 5200 m/s and 7,870 kg/m 3 , respectively. The sound velocity of the surrounding medium was 1500 m/s, the density of which was 1,000 kg/m 3 , which is close to that of soft tissues or water. Therefore, the scatterers had a severe impedance mismatch with the background media, which resulted in strong scattering. The imaging targets were five optical absorbers with a diameter of d = 0.8 mm in the ROI within the scattering layer. The coordinates of these targets were (−2, image value for points corresponding to the array can be obtained by applying the time reversal operator. If point A is within the imaging area, we can find two imaging points, B and C, closest to point A in the focal plane and then perform a weighted evaluation based on the distance from point A to points B and C: where dAC and dAB represent the distances between A, B and A, C, respectively. If point A is outside the imaging area, the imaging value is then recorded as 0. When we calculate the imaging value of point A for all sector areas, its final value is the result of taking the summation and average to all the non-zero values: where NZ is the number of zero values. As such, we can fully use the valid information received by all the transducers of the annular array for image reconstruction of each point in the ROI.

Numerical Simulation
Numerical simulations were used to validate the proposed method. The simulation is depicted in Figure 2. The radius of the annular array was set to Ra = 40 mm, and there were N = 576 elements in the array. The acoustic scattering layer consisted of 120 acoustic scatterers with a diameter of d = 0.8 mm. The inner and outside radius of the scattering layer were R1 = 10 mm and R2 = 15 mm, respectively. The sound velocity and the density of these scatterers were 5200 m/s and 7,870 kg/m 3 , respectively. The sound velocity of the surrounding medium was 1500 m/s, the density of which was 1,000 kg/m 3 , which is close to that of soft tissues or water. Therefore, the scatterers had a severe impedance mismatch with the background media, which resulted in strong scattering. The imaging targets were five optical absorbers with a diameter of d = 0.8 mm in the ROI within the scattering layer. The coordinates of these targets were (-2,-2), (-2,2), (2,-2), (2,2), and (0,0). Under the illumination of the laser pulse, these optical absorbers in the scattering layer emitted ultrasonic pulses with a central frequency of 2.0 MHz and a -3 dB bandwidth of 1.26-2.68 MHz to the surrounding medium. The PA signals were subsequently detected by the annular array enclosing the region of interest. The PA signals awee sampled with a frequency of 40 Mhz. The imaging targets had the same sound velocity and the same density as the scatterers.   Figure 3a illustrates the typical PA signals detected by one transducer. The measured broadband signal was partly overwhelmed by the noise and the amplitude was unstable. The PA signals detected by all transducers are shown in Figure 3b to visualize the effect of the scattering layer on the received signals, where the gray dotted line corresponds to the results in Figure 3a. The coherence of the signals was broken by the randomness of the scattering. The presence of a scattering layer inevitably causes a drop in the image quality.For the sake of comparison, Figure 3c illustrates the imaging results reconstructed by the classic delay-and-sum (DAS) beamforming method. As shown, for the DAS beamforming method, the image reconstructed in the case with a scattering layer (Figure 3c) suffered from low resolution, and the background showed strong noise. As a result, the contrast between the targets and the background decreased distinctly. Figure 3d depicts the image recovered by the time reversal method combined with a correlation matrix filter. The image shows good contrast, as the background noise has been significantly reduced after filtering.  Figure 3a illustrates the typical PA signals detected by one transducer. The measured broadband signal was partly overwhelmed by the noise and the amplitude was unstable. The PA signals detected by all transducers are shown in Figure 3b to visualize the effect of the scattering layer on the received signals, where the gray dotted line corresponds to the results in Figure 3a. The coherence of the signals was broken by the randomness of the scattering. The presence of a scattering layer inevitably causes a drop in the image quality.For the sake of comparison, Figure 3c illustrates the imaging results reconstructed by the classic delay-and-sum (DAS) beamforming method. As shown, for the DAS beamforming method, the image reconstructed in the case with a scattering layer (Figure 3c) suffered from low resolution, and the background showed strong noise. As a result, the contrast between the targets and the background decreased distinctly. Figure 3d depicts the image recovered by the time reversal method combined with a correlation matrix filter. The image shows good contrast, as the background noise has been significantly reduced after filtering. We used the signal-to-noise ratio (SNR) and the full width at half maximum (FWHM) of the targets' image to quantify the image quality. The SNR is defined as the ratio of the mean signal intensity in the signal area and the mean noise intensity in the background area. We defined the area inside the circles with a diameter of 0.8 mm (equal to the actual size of the targets) at the locations of the five targets as the signal area and the area outside as the background area. The averaged FWHM along the X and Z axes of all targets was calculated to quantify the degree of the distortion. For the results in Figure 3c, the SNR of the image was 17.2 dB, and the averaged FWHM along the X and Z axes were 1.42 mm and 2.6 mm, respectively. The difference between the FWHM along the X and Z axes were caused by the non-uniform distribution of the scatterers in the scattering layer. The SNR of the image in Figure 3d increased to 21.0 dB due to the matrix filter. In addition, the FWHM along the X and Z axes were 0.87 mm and 1.63 mm, respectively, which are both closer to the actual size of the targets. The results demonstrate that the matrix filtering method can more accurate depict the shape and position of the targets in the ROI.

Experimental Study
We also performed an experiment on a vessel shape source to validate the performance of the proposed method. A schematic diagram of the experimental system is provided in Figure 4a. A vessel-shaped source placed in the ROI was irradiated by a Q-switched Nd: YAG pulse laser with a wavelength of 532 nm. The pulse width was about 8 ns. The ultrasonic waves were detected by a 5 We used the signal-to-noise ratio (SNR) and the full width at half maximum (FWHM) of the targets' image to quantify the image quality. The SNR is defined as the ratio of the mean signal intensity in the signal area and the mean noise intensity in the background area. We defined the area inside the circles with a diameter of 0.8 mm (equal to the actual size of the targets) at the locations of the five targets as the signal area and the area outside as the background area. The averaged FWHM along the X and Z axes of all targets was calculated to quantify the degree of the distortion. For the results in Figure 3c, the SNR of the image was 17.2 dB, and the averaged FWHM along the X and Z axes were 1.42 mm and 2.6 mm, respectively. The difference between the FWHM along the X and Z axes were caused by the non-uniform distribution of the scatterers in the scattering layer. The SNR of the image in Figure 3d increased to 21.0 dB due to the matrix filter. In addition, the FWHM along the X and Z axes were 0.87 mm and 1.63 mm, respectively, which are both closer to the actual size of the targets. The results demonstrate that the matrix filtering method can more accurate depict the shape and position of the targets in the ROI.

Experimental Study
We also performed an experiment on a vessel shape source to validate the performance of the proposed method. A schematic diagram of the experimental system is provided in Figure 4a   To mimic the effect of the skull on brain imaging, we placed the sample in a glass beaker. The glass has a thickness of about 1.0 mm, and the corresponding speed of sound and density were 2,730 m/s and 2,500 kg/m 3 , respectively. The acoustic impedance of the glass beaker was much higher than that of the agar and the surrounding water. It caused strong acoustic scattering of the PA waves, as would the skull. Serious acoustic impedance mismatch between the glass beaker and the surrounding media could result in a significant drop in image quality.
For comparison, the images recovered by the DAS method and the proposed method in the cases with and without the beaker are shown in Figure 5. Figures 5a,b show the imaging results of the DAS method and the proposed method without the beaker. In both cases, the pattern is clear enough that the shape of the source is distinguishable.   Figure 4b is a picture of the imaging sample. A vessel-shaped source made of hair was embedded in agar for fixing, with an agar to water ratio of about 1.2%, and the sound velocity was measured as 1,460 m/s at 25°C. To mimic the effect of the skull on brain imaging, we placed the sample in a glass beaker. The glass has a thickness of about 1.0 mm, and the corresponding speed of sound and density were 2,730 m/s and 2,500 kg/m 3 , respectively. The acoustic impedance of the glass beaker was much higher than that of the agar and the surrounding water. It caused strong acoustic scattering of the PA waves, as would the skull. Serious acoustic impedance mismatch between the glass beaker and the surrounding media could result in a significant drop in image quality.
For comparison, the images recovered by the DAS method and the proposed method in the cases with and without the beaker are shown in Figure 5. Figure 5a,b show the imaging results of the DAS method and the proposed method without the beaker. In both cases, the pattern is clear enough that the shape of the source is distinguishable.   Figure 4b is a picture of the imaging sample. A vessel-shaped source made of hair was embedded in agar for fixing, with an agar to water ratio of about 1.2%, and the sound velocity was measured as 1,460 m/s at 25 ℃. To mimic the effect of the skull on brain imaging, we placed the sample in a glass beaker. The glass has a thickness of about 1.0 mm, and the corresponding speed of sound and density were 2,730 m/s and 2,500 kg/m 3 , respectively. The acoustic impedance of the glass beaker was much higher than that of the agar and the surrounding water. It caused strong acoustic scattering of the PA waves, as would the skull. Serious acoustic impedance mismatch between the glass beaker and the surrounding media could result in a significant drop in image quality.
For comparison, the images recovered by the DAS method and the proposed method in the cases with and without the beaker are shown in Figure 5. Figures 5a,b show the imaging results of the DAS method and the proposed method without the beaker. In both cases, the pattern is clear enough that the shape of the source is distinguishable.  However, for the sample in a beaker, the image recovered by the DAS method (Figure 5c) shows low contrast and strong background noise. The shape of the source is hard to distinguish. In the imaging result produced by our proposed method, we obtained a clear pattern of the vessel shape source with high contrast, which is similar to the result in Figure 5b. We conclude that the proposed method can effectively reduce the degradation of imaging quality caused by the beaker, and the experiments validated the performance of the proposed method.

Discussion and Conclusions
In summary, we developed an image reconstruction method for PAT. This method combines a correlation matrix filter and a time reversal operator. We applied this method to reconstruct PA images inside a scattering layer. Both numerical simulation and phantom experiments were used to examine the performance of the proposed method. The simulation quantitively validated that the proposed imaging method provides higher imaging quality in complex scattering media compared with the classic beamforming method. The advantages are reflected in better a SNR and FWHM.
Though demonstrating potential in preclinical and clinical applications, PAT is still in its early age of development. Among the several limitations preventing PAT's application, the main one is the amplitude and phase distortion of the received signals caused by the randomly distributed scatterers. In this research, we proposed an image reconstruction method to reduce the negative effect of the scattering waves by separating the direct waves from the received signals. With the filtered signals, PA images can be optimized, especially the resolution and the contrast between the targets and the background. This work might be valuable for understanding the physics of the interaction of sound in the complex media. The proposed method may potentially be used for acquiring high-quality images of the brain through the skull.
In this study, the experiments were performed by a point-by-point scanning with a single ultrasound transducer. Point-by-point scanning is not fast enough to meet the imaging requirements in pre-clinical or clinical scenarios. This problem can be solved by using a circular ultrasound array, which is what we are preparing to study in our next work. Additionally, as future work, we will evaluate the proposed method in an in vivo study on a small animal.