Estimation of blood oxygenation with learned spectral decoloring for quantitative photoacoustic imaging (LSD-qPAI)

One of the main applications of photoacoustic (PA) imaging is the recovery of functional tissue properties, such as blood oxygenation (sO2). This is typically achieved by linear spectral unmixing of relevant chromophores from multispectral photoacoustic images. Despite the progress that has been made towards quantitative PA imaging (qPAI), most sO2 estimation methods yield poor results in realistic settings. In this work, we tackle the challenge by employing learned spectral decoloring for quantitative photoacoustic imaging (LSD-qPAI) to obtain quantitative estimates for blood oxygenation. LSD-qPAI computes sO2 directly from pixel-wise initial pressure spectra Sp0, which are vectors comprised of the initial pressure at the same spatial location over all recorded wavelengths. Initial results suggest that LSD-qPAI is able to obtain accurate sO2 estimates directly from multispectral photoacoustic measurements in silico and plausible estimates in vivo.


. Introduction
Photoacoustic (PA) imaging is a medical imaging modality that o ers an optical signal contrast up to several centimeters deep inside tissue. In the last two decades, a lot of progress has been made towards its translation into clinical routine. However, the accurate and robust quanti cation of optical tissue properties or derived functional tissue properties still remains a major challenge ( , ). In order to obtain functional tissue properties such as blood oxygen saturation (sO2), spectral unmixing algorithms are used to decompose a multispectral signal and to determine the quantity of speci c chromophores that contributed to the signal spectrum over the wavelengths. For this process, the core assumption is that the acquired PA signal image I -which is ideally an approximate reconstruction from the initial pressure distribution p 0 -is proportional to the optical absorption coe cient (I ≈ p 0 ∝ µa) ( ). However, this assumption does not always hold because p 0 is mainly proportional to the product of optical absorption µa and the light uence φ (I ≈ p 0 ∝ µa · φ). Hence, the uence has an inuence on the recorded spectra and potentially leads to large errors in sO2 quanti cation.
To overcome this issue, part of the research in the eld of quantitative PA imaging (qPAI) aims at compensating for these uence e ects (cf. e.g. ( -)). For a long time, the eld has focused on model-based approaches to extract quantitative information on optical tissue properties ( , ) which typically su er from long computation times. More recently, also machine learning-based approaches to qPAI have been published ( -). These methods are often substantially faster than model-based algorithms but have not yet been demonstrated to work accurately and robustly in realistic in vitro settings or in vivo. Due to the lack of fast and accurate uence compensation algorithms, most researchers and applications default to the simple linear unmixing algorithms in order to provide qualitative rather than quantitative values for relevant functional tissue parameters.
This work tackles quanti cation of the functional tissue property sO2 by introducing learned spectral decoloring for quantitative photoacoustic imaging (LSD-qPAI). LSD-qPAI is based on the assumption that there is a substantial bene t in considering uence e ects when quantifying sO2 ( , , ), compared to spectrally unmixing hemoglobin (Hb) and oxyhemoglobin (HbO2) with commonly used linear methods ( ) that neglect the aforementioned uence e ects. Previous deep learning approaches to qPAI try to estimate optical absorption from PA measurements and then derive functional tissue properties from these estimations, potentially leading to error propagation ( ). In contrast to this, we propose to directly estimate the functional tissue parameter sO2 from pixel-wise p 0 spectra Sp 0 , which is a vector of the initial pressure at the same spatial location over all recorded wavelengths. For this, we use multispectral in silico p 0 training data, as illustrated in Figure . This way, we force the machine learning algorithm to account for uence e ects in the p 0 spectra during sO2 estimation. According to initial results, our in vivo sO2 estimates are physiologically more plausible when compared to linear spectral unmixing techniques.

. Materials and Methods
A. Concept overview. With LSD-qPAI we approximate a function f to estimate sO2 from initial pressure p 0 spectra Sp 0 (f : Sp 0 ∈ R n → sO2 ∈ R). It is a data-driven method in which a neural network learns to compensate for di erent extents of spectral coloring within a given p 0 spectrum, where the term spectral coloring refers to changes in the spectrum at a given spatial location due to wavelength-dependent absorption in the surrounding tissue ( ). To account for the lack of real data comprising reliable reference or even ground truth sO2 values given p 0 recordings, we create a dataset of pixel-wise p 0 spectra with various degrees of spectral coloring obtained from Monte Carlo simulated in silico multispectral p 0 images (cf. Figure ). The variety of possibilities for spectral coloring is simulated by extracting p 0 spectra from di erent spatial locations within the image (e.g. shallow or deep background structures and super cial of deep vascular structures). With the generation of a representative dataset, the network learns to account for spectral coloring when inferring sO2 from input p 0 spectra. During training, the algorithm is presented tuples (Sp 0 , sO2), with Sp 0 ∈ R n and sO2 ∈ R. Here, each spectrum is normalized such that n i p 0λ i = 1. When estimating oxygenation for recorded in vivo spectra these are normalized as wellsacri cing p 0 amplitude information to eliminate the need to calibrate the in silico p 0 training data to the target domain.
B. Prototype implementation. In the following paragraphs, we provide detailed descriptions of the prototype implementation of our approach, namely the simulation pipeline, the used deep learning model, and a linear spectral unmixing method.
Simulation pipeline. We create simple homogeneous tissue volumes comprising tubular vessel structures that run orthogonal to the imaging plane. We simulated the light transport in this medium with a Monte Carlo method and simulated each of the in silico vessel phantoms with wavelengths equidistant from nm to nm in nm steps. For the Monte Carlo simulation, we used a multithreaded adaptation of the Monte Carlo framework mcxyz ( ) with 10 7 photons for each simulation on a . mm grid.
Deep learning model. For reconstruction of sO2 from Sp 0 , we use a simple fully connected neural network architecture (cf. Figure ). We implement the network with eight hidden layers with four times the size of the input layer using pytorch ( ) and perform the experiments with the trixi framework ( ). We use an L loss function, a learning rate of −4 , a batch size of , with batches per epoch, and train the network for epochs. The size of the input layers corresponds to the number of simulated or measured wavelengths ( in this study) and we have a one dimensional output layer, which corresponds to the sO2 estimation as the target parameter. Linear spectral unmixing. We use a linear spectral unmixing algorithm as a reference for the in vivo imaging results. Specifically, we use the C Eigen ( ) implementation of the QR householding matrix decomposition to unmix the recorded spectra for Hb and HbO2 and calculate sO2 from the ratio of these.

. Experiments and Results
In the following section we outline the experiments we conducted for validation of the presented approach. The purpose of our experiments was two-fold: we rst show the feasibility of our method with an in silico evaluation on a held-out set from the same distribution as the training data. Secondly, we conduct feasibility tests on two di erent in vivo datasets, of which the rst contains PA images of an open porcine brain and the second dataset comprises PA images of the human forearm. These datasets are chosen because the compositions of relevant chromophores in tissue are very di erent in both cases. For both of these in vivo settings we also compute sO2 with a linear spectral unmixing algorithm to provide a reference.

Simulation parameters.
The key to the general applicability of the LSD-qPAI method is the simulation of a representative dataset of p 0 spectra. For this initial study, we created a total of in silico phantoms containing -tubular vesselmimicking structures in a homogeneous background at random locations orthogonal to the imaging plane. For tissue mimicking a priori conditions, we assumed blood vessels to have a hemoglobin concentration of g/dl ( ), generic background tissue to have a blood volume fraction of , and also consider an average of fat and water in this background tissue. The tubular vessel structures have a radius randomly drawn from a uniform distribution between . mm and mm. Each vessel phantom was assigned a distinct constant blood oxygenation in vessel structures and the same constant blood oxygenation in the background. The tissue absorption coe cient was calculated for each voxel based on these assumptions and we set the reduced scattering coe cient to be constant over all wavelengths with µ s = cm −1 .
In silico validation. We use of the in silico dataset for training, for supervision of the training process and hyperparameter optimization, and the remaining as a test set to report the in silico results. We report the relative sO2 estimation error e sO 2 = (ŝ O2 − sO2)/sO2, whereŝ O2 represents the estimated blood oxygen saturation and sO2 represents the ground truth simulated blood oxygenation. All in silico estimation results for the test data are shown in Figure . Additionally, violin plots present the error distribution within each interval of ten percentage points ( -, -, etc). The median relative estimation error e sO 2 was . with an interquartile range of ( . , . ). Imaging of porcine brain. As our rst in vivo experiment, we applied our method to images of a porcine brain during open surgery. These images were recorded at the same wavelengths as in the training dataset. We evaluate a single series of images which were normalized by the recorded laser energy and reconstructed with the delay-and-sum algorithm using a hamming window. For the calculation of the mean oxygenation within the ROI, only those point estimates were taken into consideration, where the signal at the isosbestic point of nm was greater than a noise equivalent level. In the image, the red bounding box marks an area where blood with a high oxygenation was present (cf. Figure ). Our method estimated a blood oxygenation of about , whereas linear spectral unmixing using a QR decomposition yielded an estimate of in the same area.
Imaging of human forearm. We also applied our method to in vivo images of a the forearm of a healthy human volunteer also with the same wavelengths as used in the simulations. This imaging scenario was speci cally chosen as an out-ofdistribution test for the method. This is, because no well matching spectra should be contained in our training set, as we did not consider the presence of melanin in our simulation. We normalized each image by the laser energy and reconstructed the images with the delay-and-sum algorithm using a hamming window. To decrease the inter frame variability, we averaged over consecutive frames of the same wavelength. Prior to averaging we registered the images with an optical ow-based method to account for motion artifacts ( ).
In the oxygenation images shown in Figure , the upper signal originating from the radial artery is chosen as the region of interest (ROI) and marked by a red bounding box. In this ROI, only those point estimates were taken into consideration for the averaged sO2 result, where the signal at nm was greater than a noise equivalent level. The arterial blood oxygenation as determined by our algorithm was about , whereas spectral unmixing using a QR linear unmixing algorithm yielded a value of about .

. Discussion and Conclusion
In this work, we investigate the feasibility of inferring sO2 from pixel-wise initial pressure spectra Sp 0 using a deep learning model trained on simulated data. While a photoacoustic signal is proportional to the product of optical absorption µa and the light uence φ, standard unmixing algorithms ignore the dependence on φ, which can potentially result in large errors, especially deep in tissue. To overcome this issue, we propose the LSD-qPAI algorithm that learns a mapping function from pixel-wise multispectral initial pressure data to blood oxygenation. Our initial in silico results look promising, yielding a median sO2 estimation error of in silico. The LSD-qPAI method is potentially relevant especially for clinical imaging with handheld PA scanners, as sO2 estimations with high accuracy and robustness can potentially be computed in real time. This has been demonstrated in similar di use re ectance multispectral imaging applications ( ).
In prior work, we achieved a comparable accuracy in silico when estimating the optical absorption coe cient ( , ). However, in contrast to these methods, we now are also able to obtain plausible results in vivo when imaging an open pig brain. Here, the estimations of our method are closer to the expected arterial blood oxygenation values of healthy subjects (near ( )). Especially in comparison to the results of classical linear unmixing, LSD-qPAI seems to provide physiologically more plausible estimations.
Our method utilizes a normalization of the simulated and recorded spectra, sacri cing p 0 amplitude information to eliminate the need to calibrate the in silico p 0 training data to the target domain. However, this also means that we discard one dimension of the feature vector, restricting the minimum amount of spectra needed to be able to reliably infer optical properties.
For the in vivo experiments it was to be expected that the recorded in vivo spectra of the brain images match our training distribution more closely compared to the forearm images. This is because light needs to penetrate through the skin in order to obtain images of the forearm, which contains melanin, a chromophore that was not included in the simulation framework. It can be seen that neither the proposed method nor linear unmixing can handle skin tissue well, estimating implausible and even impossible sO2 results. Projecting the in vivo spectra to the rst two principal components of the training data seems to con rm this hypothesis (cf. Fig. ). In the gure it is apparent that most of the porcine spectra (colored in black) are contained within the support of the rst two principal components of our simulation space, whereas this does not appear to be the case for the forearm spectra. In this context it should be noted that the rst two principal components (computed on the training data) account for .
of the variation in our dataset. The projection image illustrates that the distribution of the training data match the distribution of the test (in vivo) data more closely when imaging pig brain instead of human forearm.
When considering the large number of plausible tissue geometry and oxygenation con gurations, the number of training samples used for this study was low. Also, no noise model was applied to the simulated spectra in addition to the simulation noise inherent to the Monte Carlo procedure. Hence, a thorough investigation of the in uence of a realistic noise model on the spectra would be very interesting. Our a priori assumptions for the tissue parameters used in the simulation pipeline led to very di cult situations for the unmixing algorithm to resolve. For example, due to the blood volume fraction in the background, light uence was the dominating factor even in shallow point absorbers -a behavior not usually observed in realistic scenarios, where less than would be a more realistic assumption ( ). Also, the simulation assumption that blood oxygenation is constant throughout the tissue is not generally correct and more realistic variations on these assumptions should be considered in future work.
Overall, the presented initial results in vivo are very promising. However, a thorough and well-designed in vivo validation needs to be performed to deduce meaningful conclusions regarding the general applicability of the presented method. Such future studies should cover a much broader range of possible sO2 values and include more realistic scenarios when assessing the unmixing accuracy. Future work should also include a comprehensive comparison of this method to state-of-theart methods, such as eigenspectra multispectral optoacoustic tomography (eMSOT) ( ), the end-to-end qPAI method presented by Cai et al. ( ), as well as other linear and nonlinear spectral unmixing techniques ( , ). -StG-. The authors would like to thank E. Santos, M. Herrera, and A. Hernández-Aguilera for provisioning of PA brain data, K. Dreher, N. Holzwarth, A. Klein, and S. Onogur for proofreading the manuscript and the ITCF of the DKFZ for enabling extensive use of their computing cluster for data simulation. . Wang, L. V. and Hu, S., "Photoacoustic tomography: in vivo imaging from organelles to organs," science ( ), -( ). . Li, M., Tang, Y., and Yao, J., "Photoacoustic tomography of blood oxygenation: