Skull acoustic aberration correction in photoacoustic microscopy using a vector space similarity model: a proof-of-concept simulation study.

Skull bone represents a highly acoustical impedance mismatch and a dispersive barrier for the propagation of acoustic waves. Skull distorts the amplitude and phase information of the received waves at different frequencies in a transcranial brain imaging. We study a novel algorithm based on vector space similarity model for the compensation of the skull-induced distortions in transcranial photoacoustic microscopy. The results of the algorithm tested on a simplified numerical skull phantom, demonstrate a fully recovered vasculature with the recovery rate of 91.9%.

The findings about the effect of skull on acoustic pressure waves are as follows. Skull bone represents a highly acoustical impedance mismatch and dispersive barrier for the propagation of acoustic waves [23][24][25], that distorts the amplitude and phase of the received acoustic waves [20,21]; the acoustic attenuation occurs due to the absorption and scattering of the skull tissue and affects the magnitude of the acoustic waves [19,26,27]; the acoustic dispersion is the frequency dependency of the speed of sound in the skull and it distorts the phase of the acoustic wave [26]; the degree of attenuation and dispersion are defined by the density, porosity, and thickness of the skull [28]; frequency-dependent reduction of acoustic wave amplitude contributes to the broadening of the received acoustic signal at the transducer [29]; the significantly higher speed of sound in the bone ( 2900 m/s [16]) as compared to brain's soft tissue ( 1500 m/s [30]) makes the acoustic waves travel faster through the skull and be detected earlier, leading to a different time shift for individual frequency components, and contributes to broadening of the received acoustic signal at the transducer [9,[19][20][21]; skull-originated reverberations occur due to the reflection of acoustic waves from the skull-tissue interface [9,20,21]; longitudinal to shear wave mode conversion occurs when a wave encounters an interface between materials of different acoustic impedances with the incident angle not being normal to the interface [22].
Several analytical and numerical compensation algorithms have been designed for correcting the skull-induced aberrations in transcranial ultrasound wave propagation. These methods can be broadly classified to ray-based [19,31] and wave-based [31][32][33][34][35][36][37]. Ray-based methods calculate the correcting phases through ray tracing and numerical simulation. In the wave-based methods such as time-reversal, the pressure waveform over time is recorded, reversed, and re-emitted to focus on the location of the desired target, providing corrections for both phase and amplitude. A review of phase correction methods can be found in [38]. In addition to the above methods, recently, a deconvolution-based algorithm, was proposed by Estrada et al. for skull-induced distortions correction in transcranial optoacoustic microscopy (OAM) [20].
The skull's aberration correction methods implemented so far, are either fast but not accurate, time consuming or require an axillary imaging modality to acquire the structural information of the imaging target. On the other hand, among modern computational methods, although used for improved image reconstruction [39][40][41][42][43][44][45][46][47], there has not yet been any study on machine learning (ML) for skull's aberration correction. Here, we introduce a novel algorithm based on the vector space similarity (VSS) model, for the first time, in conjunction with a ray-tracing simulation to correct for the skull-induced distortions for the images generated by a TsPAM. VSS has been used, for the first time, in Cornell system for the mechanical analysis and retrieval of text (SMART) in 1960s [48]; it has widely been used in intelligent information retrieval from search engines [49] to big data platforms such as biomedical documents [50]. We employ a modern take on of the VSS model in the context of matching between the extremum information in the distorted skull-induced PA time-frequency domain signal and the reference signals generated by our recently developed ray-tracing-based simulation [21,[51][52][53]. Additional justification for the application of VSS model is based on the successful use of it in PAM post-processing for tissue vasculature classification [54], and quantification of tissue response to cancer treatment [55].

PA wave propagation model
For the purpose of the simulation of the PA initial pressure wave propagation from a point source, a semi-analytical numerical acoustic solver is used that was recently developed by us [21]. The solver is based on a deterministic ray-tracing approach in the time-frequency domain that considers a homogenized single layer model for the skull, taking into account dispersion, reflection, refraction and mode conversion between the skull surfaces. Attenuation of light due to the skull or depth is ignored and it is assumed that sufficient initial pressure is generated at the imaging target location, and the acoustic attenuation for the initial pressure traveling towards the skull surface has been considered.
The imaging target is assumed to be a combination of many point sources. The impulse rays propagated to the first fluid-solid interface, multiplied by the skull's transmission coefficient and propagated further. This procedure is repeated until all the rays reach the transducer's surface. The PA signal is produced by convolving each individual frequency component of the initial pressure with its corresponding impulse response. When the skull is not present, the ray is simply propagated through a free space and travels directly towards the transducer's surface. Figure 1 shows the 2-D illustration of the simplified numerical skull phantom that is used in our simulations. The imaging target is a vessel represented by 27007 point sources. The diameter of the vessel varies between 30 µm and 150 µm, and the entire vessel is located at the averaged depth of 5 mm. In this simulation the axial and lateral resolutions are considered 6 µm and 10 µm, respectively. We modeled every initial pressure point source as a sphere with the radius of r 0 . The spherical point sources generate broadband spherical acoustic waves. A single layer homogenized skull tissue with a thickness variation between 0.5 mm and 2.5 mm is considered. The outer-skull surface considered to be flat but the inner-skull surface is quantized as locally-flat surfaces. This is done to approximately study the effect of varying skull thickness on the resulting image. The space between the skull and the target is modeled with the same acoustic properties as the brain soft tissue [22,30]. A flat ultrasound transducer with an element diameter of 2r d , 30 MHz center frequency, and 100% bandwidth is used. The transducer is in contact with the outer-skull surface through ultrasound gel. The acoustic properties of the skull, brain soft tissue, and the ultrasound gel that are used in our simulations are listed in Table 1. The initial pressure point source modeled as a sphere with the radius of r 0 at a depth of d from the outer-skull surface. A single layer homogenized skull tissue with a thickness variation between 0.5 mm and 2.5 mm is considered. The solid acceptance angle of the transducer is indicated by the dashed lines.

Vector space similarity model
Our proposed method is based on VSS model, in which no direct signal amplification or shift will be performed, instead the compensated signal is reconstructed according to the similarity between the skull-affected signal and the signals in the training dataset [49,57]. To describe the VSS model, let's suppose we have a numeric database, q = (f q1 , f q2 , . . . , f qn ) , and the goal is to find the most similar data from the training dataset, i.e., . . , f in ), to the desired query with the defined similarity features. Where, d i is the i th data in the training dataset, q is the desired query, and f ij and f qj are the j th feature vectors of d i and q, respectively. The dot product of each query feature vector in all corresponding feature vectors of the training dataset is calculated, then the cosine similarity measure is used to find the minimum angle between query and training dataset as indicated in Eqs. (1) and (2); "." is the notation of the intersection or dot product and " " is the notation of the norm of vector.

Aberration correction algorithm
The aberration correction algorithm uses the PA signal extremum information in the timefrequency domain as feature vector. The algorithm is as follow. The input to the algorithm is "with skull" signals. The signals are initially decomposed to their time-frequency components using the short-time Fourier transform (STFT) [58]. For the implementation of the STFT in this study, the signals are discretized at a rate of 250 Megasamples per second and 32 frequencies are modeled. Then for each frequency given, vectors of features, namely TimeVector and AmpVector, are extracted as follow.

PA signal vector space (For real and imaginary parts individually)
For each frequency f i in the signal: TimeVector: Temporal delay Time points at which minimum occurs (mnpt 1 ,mnpt 2 ,. . . ,mnpt k ) Time points at which maximum occurs (mxpt 1 ,mxpt 2 ,. . . ,mxpt k ) AmpVector: Minimum peak amplitudes (mnpa 1 ,mnpa 2 ,. . . mnpa k ) Maximum peak amplitudes (mxpa 1 ,mxpa 2 ,. . . ,mxpa k ) Then, the dot product of the obtained feature vectors in each of the corresponding vectors of the training dataset (d i ) are calculated and divided by the norms of the vectors to yield the cosine of the angle between them. The similarity is then calculated via the arccos of the obtained value. For each frequency, the overall similarity of each reference signal from the training dataset to the query data (q) (the skull affected signal) is calculated as the mean of its feature vector similarity (see Eq. (3)). The mean similarity is then calculated between similarities in all frequencies for each two pairs and created a similarity metric vector.
Where, T i and T q are the TimeVectors of d i and q, respectively, and A i and A q are the AmpVectors of d i and q, respectively.
Finally, the minimum of the similarity metric vector is used to select the best matched reference signal. For each frequency, according to the temporal delay and maximum amplitude difference between the input signal and the corresponding matched reference signal, the appropriate amplification coefficients and time shifts are performed on the "without skull" signal (see Algorithm (1)).

Input
:ws -photoacoustic signal with skull. ts -training dataset of with/without tuples Output :cwo -compensated signal Data: ws, cwo are time traces of photoacoustic signal. ts is a set of with-without tuples decomposed by frequencies and represented by feature vectors namely TimeVector (signal temporal delay and time points at which signal extremum occurs) and AmpVector (signal extremum amplitudes We created a large training dataset of "without skull" and "with skull" signals (with different skull thicknesses (from 0.3 mm to 2.5 mm with 0.1 mm step size), and different imaging depths (from 0.1 mm to 25 mm below the skull surface with 0.5 mm step size)), and then extracted the abovementioned features from the "with skull" signals, creating 1150 training dataset.

Skull-induced aberration compensation: investigation on time traces
Initially, we evaluated the compensation algorithm in correcting the induced aberration of the PA signal produced from a 0.1 mm absorbing sphere and passes through a skull tissue. The tested scenarios are as follows. (i) Skull thickness of 0.5, 1, and 2 mm, when the absorbing sphere is located at the depth of 5 mm, (ii) skull thickness of 1 mm, when the absorbing sphere is located at the depths of 2, 8, and 20 mm. The unaberrated (i.e., without skull), aberrated (i.e., with skull) and compensated signals are shown in Fig. 2 and Fig. 3.  We extracted the PA background noise, comprised of a combination of electronic noise and system thermal noise from an experimental setup, and added that to our simulated signals to evaluate the tolerance of our compensation algorithm to noise. The experimental setup was as follow. An Nd:YAG laser (PhocusMobil, Opotek Inc., CA, USA) with a repetition rate of 10 Hz and a pulse width of 8 ns was used. For light delivery, we used a custom fiber bundle (Newport Corporation, Irvine, CA, USA). For data acquisition, a Verasonics Vantage 128 system was used. For PA signal detection, L22-14v linear array (Verasonics Inc., USA) ultrasound probe with 128-elements and 18.5 MHz central frequency and 65% bandwidth was used. On the imaging end, the transducer was placed and held perpendicularly to the sample. A 2 mm diameter carbon lead phantom in water was imaged at 690 nm. The noise was extracted following the deconvolution algorithm explained in [59]. We used 100 frames of data and modeled the noise distribution. The noise signal was normalized, and two levels of noise were formed: 10% and 20%. The noise signals were then added to some of the distorted signals in Fig. 2 and Fig. 3. The compensation algorithm were applied to the noisy PA signals. The results of this experiment are shown in Fig. 4.
In order to improve the accuracy of the compensation algorithm, we considered a pre-processing step before compensation, where we thresholded low amplitude samples (< 5% of the signal peak) to zero.
The results in Fig. 2 and 3 show that the signal distortion due to both time shift and amplitude distortion have well been recovered. The gradient of the recovered signal and undistorted signals are almost identical which suggest that both the depth of target and skull thickness do not introduce phase distortions to the recovered signal. Figure 4 shows that with a noisy distorted signal, the compensation algorithm still recovers the undistorted signal. The phase information however is not affected unless there is a steep rise or fall in the signal; such distortion is translated in displacement of the components of the imaging target in axial direction and slight speckle-looking artifact in the image. Comparing the results in Figs. 4 a, b, c, and d, one can conclude that the amplitude recovery of our proposed compensation algorithm is more sensitive than time shift recovery.
To quantify the performance of the compensation algorithm, we defined a quantitative measure, we called it: "recovery percentage", calculated as: Where, l 2 norm is calculated as the square root of the sum of the squared signal sample values. Table 2, shows the recovery percentage for the results showed in Figures (2-4).

Skull-induced aberration compensation: investigation on synthetic TsPAM images
The compensation algorithm was then applied to aberrated signals collected from a synthetic PAM experiment (see the setup and the imaging target in Fig. 1). The results are shown in Fig. 5.
A representative depth profile of the unaberrated, aberrated and compensated images (indicated Fig. 4. Skull-induced aberration compensation of PA signals produced from a 0.1 mm absorbing sphere passing through a skull tissue. (a) Skull thickness is 2 mm and the target is located at 5 mm depth (PA signal is contaminated with 10% background noise), (b) skull thickness is 2 mm and the target is located at 5 mm depth (PA signal is contaminated with 20% background noise), (c) skull thickness is 1 mm and the target is located at 20 mm depth (PA signal is contaminated with 10% background noise), and (d) skull thickness is 1 mm and the target is located at 20 mm depth (PA signal is contaminated with 20% background noise). (i) Signal amplitudes and (ii) signal gradients. In this simulation, there is a 5 mm layer of ultrasound gel between the ultrasound transducer and the skull.
with green dotted lines in the images in Fig. 5(a)) are plotted in Fig. 5(b). As can be seen in Fig. 5(b), the axial profile has been almost perfectly recovered, in terms of both amplitude and phase. The recovery percentage calculated for the compensated image was 91.9%. Using the same method explained in the results section (a), we created two noisy aberrated TsPAM images, one with 10% and one with 20% noise (see Fig. 6(ii)). We then compensated the images, the results are shown in Fig. 6(iii). The recovery percentage calculated for these compensated images are 90.60% and 87.95%, respectively.

Execution time analysis
Both the signal simulator [21] and the compensation algorithm were implemented in Java openJDK 13. All the signal processing were conducted in MATLAB R2016a. Utilizing a hexa-core Intel Core i7 CPU with 6 cores, 32 GB of RAM, and 3.60 GHz, the compensation algorithm, took 24 ± 2 ms for one task of compensation; with the current code, the compensation is done offline. By implementing the compensation algorithm in graphical processing unit (GPU) or in field programmable gate array (FPGA), hundreds or thousands fold speed-up is achievable, that could make the real-time compensation of the signal in the data acquisition line feasible.

Non-flat skull aberration compensation: preliminary results
We evaluated our compensation algorithm to correct aberrated PA signal produced from a 0.1 mm spherical absorber and passed through an angled skull tissue when the absorbing sphere was located at the depth of 5 mm from the outer-skull surface. The pressure wave travels a longer path when passes the angled skull; the larger the angle, the longer the travel time. Table 3, shows the preliminary results for skull thickness variation (i.e., ∆h), the percentage of the signal transmitted, and the recovery percentage of the distorted signal after compensation at angles θ = 5 • to 30 • relative to the transducer axis. Since the critical angle for the fluid-solid interface is about 31 • [22], we reported simulation results for angles up to 30 • only.

Discussion
We study a novel algorithm based on VSS model for quasi real-time compensation of the skull-induced distortions in transcranial photoacoustic microscopy. Although VSS is an effective and efficient similarity measurement algorithm, it is not the only similarity algorithm that is explored in the literature. There are other similarity metric such as Pearson Correlation Coefficient (PCC) [60] that works based on similar principle but with a heavier computational complexity. Both VSS and PCC have trouble in distinguishing different importance of features. To deal with this problem, many variations of similarity measurement including weighting approaches, combination measures, and rating normalization methods have been developed [61]. The skull-induced aberration compensation algorithm described here is designed for TsPAM imaging. The proposed compensation method is based on the longitudinal PA waves generated from a single point source which is placed in the axis of the point detector and the effect of other adjacent sources on the PA signal is neglected. Therefore, the transducer receives waves with only small incident angle with almost no shear wave component. The consideration of the point source can be accurately assumed if the lateral resolution of the image is governed by the diffraction-limited size of the focused light beam. The simulated numerical phantom results (Figs. 5 and 6), confirmed the ability of the proposed algorithm to accurately compensate for the four previously explored skull-induced acoustic distortions [20], including the signal amplitude attenuation, time shift, signal broadening, and multiple reflections (Figs. 2 and 3); due to the high attenuating effect of the skull, the multiple reflections cannot be visualized in simulated aberrated signals. In Fig. 3(a), only one peak can be visualized after the main peak which is because of the strength of the first reflection at the skull surface when the thickness of the skull was 0.5 mm. With thicker skulls, the attenuation induced by the skull makes the reflected signals too weak such that they cannot be visualized. The noise tolerance of the compensation algorithm was also tested. It was shown that the recovery rate is only slightly affected with the presence of 10% and 20% noise; the error appears as a mild speckle-looking artifact in the image that can easily be removed by median and mean filters. We also observed that the amplitude recovery of our proposed compensation algorithm is more sensitive than the time shift recovery.
The main feature of the algorithm is its simplicity and fast execution time that are translated in its light computation. Introducing more parallelism to the implementation of the algorithm is possible in various ways. One possible way is to convert the code from currently used CPU-based multi thread execution to GPU accelerated execution. There are two independent anchor points within the compensator code that can be utilized for this purpose: (i) using 32 cores in parallel to calculate similarity of a given signal against the training dataset through evaluation of feature-vector angles in which 32 is the number of frequency channels used, and (ii) employing a GPU/many-core system with each core responsible to check similarity with one or a predefined subset of the dataset. Obviously, use of the first method would decrease the computation time for a single A-line by a factor of 32 while the latter can decrease the entire end-to-end compensation time with a factor equal to the number of cores; this can be up to 7000 for modern GPUs. The fact that our method can easily be segregated amongst parallel threads argues for a FPGA-based hardware implementation to be plausible which can be integrated into the transducers, outputting aberration compensated image in real-time.
The proposed method is independent from the skull anatomy. This is a valuable feature of our proposed algorithm because in a TsPAM experiment, the anatomy of the skull as well as the spatial characteristics of the skull are not available. In this preliminary work, we assumed that the skull is flat and perpendicular to the transducer axis.
Although the compensation algorithm is based on preparing a training set, it is not considered a machine learning algorithm, mainly because it does not have a layered kernel to yield the compensated signal from the input aberrated signal.
The aberration compensation proposed algorithm has several limitations: (i) the skull is considered as a homogenized single layer bone with smooth surfaces and no curvature; (ii) the brain is considered as a homogenized soft tissue with constant acoustic properties; (iii) attenuation of light due to the skull or depth is ignored and it is assumed that sufficient initial pressure is generated at the imaging target location; (iv) PA signal is assumed to be generated from a single point spherical source and the effect of other adjacent sources on the PA signal is neglected; (v) the simulation framework and wave propagation model could be extended to simulate a line-shaped absorbing source to account for a realistic shape of the optical absorption source in the tissue.
In a real-world application after we acquire raw data from a TsPAM system, we will use our proposed algorithm as explained in Section 2.3, with only one change which is, adding more data to the VSS training dataset. So far we have trained the VSS algorithm only with flat skulls. In the future, by using finite-element method the skull surface will be segmented into very small regions (that are comparable to the acoustic wavelength). These small regions can each be approximated as a layer with a flat surface but angled versus the axis that connects the transducer (defined as a point detector) and the absorbing target (defined as a point source). Our simulator [21] will then generate data with different angles of the flat skull (see Table 3 representing preliminary data related to angled skull transmission) to determine how much of the incident signal is diffracted and how much of it received by the transducer; such data will be added to the VSS training dataset.

Conclusion
We developed a skull-induced aberration compensation algorithm based on vector space similarity model and ray-tracing-based simulations. The main feature of the algorithm is its simplicity and fast execution time. We demonstrated the effectiveness of the algorithm tested on numerical phantoms with a recovery percentage of 91.9%; i.e., 91.9% of the distorted signal due to the amplitude attenuation, time shift, and signal broadening were retrieved. By adding noise to the aberrated signals, the noise tolerance of the algorithm was evaluated; the recovery percentage was decreased to 90.60% (adding 10% noise) and 87.95% (adding 20% noise). Using GPU and FPGA for parallel implementation of the code, considering more sophisticated skull and tissue models, and taking into account the effect of the fluence, are the future plans of the current study.