1.

Introduction

Blood oxygen saturation (${\mathrm{sO}}_2$) is an important physiological indicator of tissue function and pathology. Often, the distribution of oxygen saturation values within a tissue is of clinical interest, and therefore, there is a demand for an imaging modality that can provide high-resolution images of ${\mathrm{sO}}_2$. For example, there is a known link between poor oxygenation in solid tumor cores and their resistance to chemotherapies, thus images of tumor blood oxygen saturation could be used to help stage cancers and monitor tumor therapies.^1,2 Some imaging modalities have been shown capable of providing limited information about or related to ${\mathrm{sO}}_2$ in tissue. Blood oxygenation level-dependent magnetic resonance imaging, which is sensitive to changes in both blood volume and venous deoxyhemoglobin concentration, can be used to image brain activity, but cannot respond to changes in oxygen saturation.³ Purely optical techniques, such as near-infrared spectroscopy and diffuse optical tomography, can be used to generate images of oxygen saturation.^4,5 However, because of high optical scattering in tissue, these modalities can only generate images with low spatial resolution beyond superficial depths.

Photoacoustic (PA) imaging is a hybrid modality that can be used to generate high-resolution images of vessels and tissue at greater imaging depths than purely optical modalities.⁶ PA image contrast depends on the optical absorption of the sample, so images of well-perfused tissues and vessels can, in principle, be used to generate images of ${\mathrm{sO}}_2$ with high specificity. However, unlike strictly optical techniques, information about the contrast in PA images is carried by acoustic waves that can propagate from deep within a tissue to its surface undergoing little scattering. In the ideal case of a perfect acoustic reconstruction, the amplitude of a voxel in a PA image can be described as

In some cases, it might be possible to measure an estimate of the fluence using an adjunct modality,⁹ but more commonly, attempts have been made to model the fluence. However, because the optical properties of a tissue sample are usually not known before imaging (the only reason the fluence is estimated at all is so that unknown information about the sample’s optical absorption coefficient can be recovered from the image data), it is difficult to model the fluence distribution. A variety of techniques have been developed to recover tissue absorption coefficients from PA images without total prior knowledge of the tissue’s optical properties. Progress toward solving this problem can be summarized into three key phases. In the first phase, one-dimensional analytical fluence models were used to estimate the fluence by taking advantage of assumed prior knowledge of some of the sample’s optical properties, or by extracting the optical properties of the most superficial layers from image data.^10–15 In the latter case, the effective attenuation coefficient of the most superficial tissue layer (assumed to be optically homogeneous) is usually estimated by fitting an exponential curve to the decay profile of the image amplitude above the region of interest (e.g., a blood vessel).

In the next phase, sample optical properties were recovered using iterative error minimization approaches.^16–19 With these techniques, knowledge of the underlying physics is used to formulate a model of image generation. The set of model parameters (which might include the concentrations of deoxyhemoglobin and oxyhemoglobin in each voxel) that minimizes the error between the images generated by the model and the experimentally acquired images are treated as estimates of the same parameters in the real images. This technique is only effective when the model of image generation is able to generate a set of simulated images very similar to the real set of images when the correct values for the chromophore concentrations are estimated. This is only possible when the image generation model is able to accurately model image acquisition in the real system. In practice, accurate models of image generation are challenging to formulate as not all aspects of the data acquisition pathway are fully characterized. Therefore, this technique has not yet been shown to be a consistently accurate method for imaging ${\mathrm{sO}}_2$ in tissue. Both iterative error-minimization and analytical techniques may require significant a priori knowledge of sample properties, such as all the different constituent chromophore types. This information is not always available when imaging tissues in vivo, and thus this requirement further reduces their viability as techniques for estimating ${\mathrm{sO}}_2$ in realistic imaging scenarios. The recent emergence of a third phase has introduced data-driven approaches for solving the problem.^20–28 With these approaches, generic models are trained to output images of ${\mathrm{sO}}_2$ or optical properties by processing a set of examples.²⁹ These data-driven models find solutions without significant a priori knowledge of sample properties and do not require the formulation of an image generation model using assumed prior knowledge of all the aspects related to image acquisition. Techniques based on data-driven models, such as deep learning, have been used to estimate ${\mathrm{sO}}_2$ from two-dimensional (2-D) PA images of simulated phantoms and tissue models.^{20–22,24,25,28} Fully connected feedforward neural networks have been trained to estimate the ${\mathrm{sO}}_2$ in individual image pixels given their PA amplitude at multiple wavelengths.²⁰ Because the fluence depends on the three-dimensional (3-D) distribution of absorbers and scatters, a pixelwise approach does not use all of the information available in an image. Encoder–decoder type networks, capable of utilizing spatial as well as spectral information, have been trained to process whole multiwavelength 2-D images of 2-D tissue models,^22,24,25,28 or 2-D images sliced from more realistic 3-D tissue models featuring reconstruction artifacts,²¹ and output a corresponding 2-D image of the ${\mathrm{sO}}_2$/optical absorption coefficient distribution. Although 2-D convolutional neural networks can take advantage of spatial information to improve estimates of ${\mathrm{sO}}_2$, networks trained on 2-D images sliced from 3-D images are missing information contained in other image slices that might improve their ability to learn a fluence correction. 3-D networks are often better at learning tasks requiring 3-D context.^30–32 Therefore, it is important to show that networks can take advantage of all four dimensions of information from a multiwavelength PA image dataset to estimate ${\mathrm{sO}}_2$. In addition to supervised learning, an unsupervised learning approach has been used to identify regions containing specific chromophores (such as oxyhemoglobin and deoxyhemoglobin) in 2-D simulated images.²³ The technique has not yet been used to estimate ${\mathrm{sO}}_2$ and has only been tested on a single simulated phantom lacking a complex distribution of absorbers and scatterers that would normally be found in in vivo imaging scenarios.

As we aim toward developing a technique for estimating 3-D ${\mathrm{sO}}_2$ distributions from in vivo image data, a more robust demonstration of a data-driven technique’s ability to acquire accurate ${\mathrm{sO}}_2$ estimates by processing whole 3-D images of realistic tissue models is desired. We trained two encoder–decoder type networks with skip connections to (1) output a 3-D image of vascular ${\mathrm{sO}}_2$ and (2) output an image of vessel locations from multiwavelength (784, 796, 808, and 820 nm) images of realistic vascular architectures immersed in three-layer skin models, featuring noise and reconstruction artifacts.

Ideally, networks would be trained on in vivo data to demonstrate their ability to cope with all the confounding effects present in real images. However, because there is no reliable technique to acquire ground truth ${\mathrm{sO}}_2$ data in vivo, generating such a dataset is very difficult. Blood flow phantoms can be used to generate images with accompanying information about the ground truth ${\mathrm{sO}}_2$.^33,34 However, these phantoms are usually much simpler than real tissue (e.g., optically homogeneous tissue backgrounds, tube-shaped vessels) and thus are not ideal for assessing whether networks can produce accurate estimates in more realistic cases. To overcome this, simulated images of realistic tissue models with known ground truths were used instead. The drawback with this approach is that simulations cannot capture every aspect of a real measurement, e.g., the noise and sensor characteristics may not be well known. Nevertheless, using training data that has been simulated in 3-D with limited-view artifacts, a gold-standard light model, realistic optical properties, and noise levels provides a good indication of the network’s ability to cope with measured data. Furthermore, given the difficulty of obtaining measured data with a ground truth, pretraining with realistic simulation data could be a very useful step prior to transfer training with a limited amount of measured data. Details about how the simulated images were generated are described in Sec. 2. Section 3 describes the network architecture and details about the training process. Section 4 describes the results.

2.

Generating Simulated Images

Ideally, a network trained to estimate ${\mathrm{sO}}_2$ from in vivo images would be capable of generating accurate estimates from a wide range of tissue samples with varying optical properties and distributions of vessels. In addition, the network should be able to do this despite the presence of reconstruction artifacts and noise. This section describes each step involved in the generation of the simulated images used in this study.

2.1.

Tissue Models

A set of several hundred tissue models, each featuring a unique vascular architecture and distribution of optical properties, were generated by immersing 3-D vessel models acquired from computed tomography (CT) images of human lung vessels into 3-D, three-layer skin models (some examples are shown in Fig. 1).^35,36 Each skin model contained three skin layers (an epidermis, dermis, and hypodermis). The thickness of each skin layer (epidermis: 0.1 to 0.3 mm, dermis: 1.3 mm to 2.9 mm, hypodermis: 0.8 mm to 2.6 mm), and the optical absorption properties of the epidermis and dermis layers were varied for each tissue model. A unique tissue model was generated for each vascular model. The equations used to calculate the optical properties of each skin layer and the vessels at each excitation wavelength (784, 796, 808, and 820 nm) are presented in Table 1 in Appendix B. These wavelengths were chosen as they fell within the near-infrared (NIR), and data were available for all skin layers at these wavelengths. The absorption properties of the epidermis layer of each tissue model were determined by choosing a random value for the melanosome volume fraction that was within expected the physiological range. The absorption properties of the dermis layer were determined by choosing random values for the blood volume fraction and dermis blood ${\mathrm{sO}}_2$ within the expected physiological range. For each tissue model, each independent vascular body was randomly assigned one of three randomly generated ${\mathrm{sO}}_2$ values between 0% and 100%. The PA efficiency throughout the tissue was set to one with no loss of generality.

Fig. 1

(a) Example of a 3-D vessel model (acquired from CT images of human lungs) used to construct 3-D tissue models. (b) Schematic of three-layer skin model used to construct tissue models.

3.

Network Architecture and Training Parameters

A convolutional encoder–decoder type network with skip connections (EDS) (shown in Fig. 2 and denoted as network $A$) was trained to output an image of the ${\mathrm{sO}}_2$ distribution in each tissue model from 3-D image data acquired at four wavelengths. Another network, network $B$, was assigned an identical architecture to network $A$ and was trained to output an image of vessel locations from the image sets (thereby segmenting the vessels). Figure 3 shows an example of the networks’ inputs and outputs. An EDS architecture was chosen for each task, as they have been shown to perform well at image-to-image regression tasks (i.e., tasks where the input data are a set of images and the output is an image).⁴⁴ The architecture takes reconstructed 3-D images of a tissue model acquired at each excitation wavelength as an input. The multiscale nature of the network allows it to capture information about features at various resolutions and use image context at multiple scales.^45–47 The network’s skip connections improve the stability of training and help retain information at finer resolutions. Finally, the network outputs a single 3-D feature map of the ${\mathrm{sO}}_2$ distribution or the vessel segmentation map.

Fig. 2

EDS network architecture. Blocks represent feature maps, where the number of feature maps generated by a convolutional layer is written above each block. Blue arrows denote convolutional layers, red arrows denote maxpooling layers, green arrows denote transposed convolutional layers, and dashed lines denote skip connections.

Fig. 3

(a) 2-D slices of 3-D $p_0$ images simulated at four wavelengths from a single tissue model used as an input for networks $A$ and $B$. (b) The corresponding 2-D slices of the 3-D outputs of the networks and the ground truths for this example.

An EDS network was trained to segment vessel locations because, although prior knowledge of the locations of vessels in the images was available for this in silico study, this information will not always be available when imaging tissues in vivo. Some technique for segmenting vessel positions from images is needed to enable the estimation of mean vessel ${\mathrm{sO}}_2$ values (the mean of the estimated ${\mathrm{sO}}_2$ values from all the voxels in a vessel) from the output map. Therefore, a vessel segmentation network was trained to show that neural networks can be used to acquire accurate mean vessel ${\mathrm{sO}}_2$ estimates without prior knowledge of vessel positions. The outputs of network $B$ were only used to enable the calculation of mean vessel ${\mathrm{sO}}_2$ values without assumed prior knowledge of vessel positions and were not used to aid the training of network $A$. As will be discussed in Sec. 4, the output of the segmentation network also provides some information about where estimates in the output ${\mathrm{sO}}_2$ map may be more uncertain. This information can be used to improve mean vessel ${\mathrm{sO}}_2$ estimates by disregarding values from these regions. Two separate networks were trained for each task to limit additional bias in the learned features that would arise from training a single network to learn both tasks simultaneously. In Ref. 22, two different loss functions were used in a single network trained to produce both an image of the vascular ${\mathrm{sO}}_2$ distribution and an image of vessel locations. A different loss function was used for each task/branch of the network, where each function was arbitrarily assigned equal weights. Training two separate networks has the benefit that it removes the need to assign arbitrary weights to multiple loss functions that may be used to train a single network.

3.1.

Training Parameters

Networks $A$ and $B$ were trained with 500 sets of images, corresponding to 500 different tissue models. An image of the true ${\mathrm{sO}}_2$ distribution of the vessels was used as the ground truth for network $A$. A binary image of true vessel locations was used as the ground truth for network $B$. Network $A$ was trained for 98 epochs (loss curve shown in Fig. 4), while network $B$ was trained for 84. Both networks were trained with a batch size of five image sets, a learning rate of ${10}^{-4}$, and with Adam as the optimizer. Training was terminated with an early stopping approach using a validation set of five examples. The networks were trained with the following error functionals, $ϵ({\theta }_A)$ and $ϵ({\theta }_B)$, (the norm of the squared difference between the network outputs and the ground truth images)

Eq. (3)

$${ϵ}_{{\theta }_A}={\Vert A_{{\theta }_A}[p_0(x,\lambda )]-{\mathrm{sO}}_2^{\text{true}}(x)\Vert }_2^2,$$

and

Eq. (4)

$${ϵ}_{{\theta }_B}=\Vert B_{{\theta }_B}[p_0(x,\lambda )]-{\mathrm{seg}}^{\text{true}}(x){\Vert }_2^2,$$

where $p_0(x,\lambda )$ are the multiwavelength images of each tissue model, ${\mathrm{sO}}_2^{\text{true}}(x)$ and ${\mathrm{seg}}^{\text{true}}(x)$ are the ground truth ${\mathrm{sO}}_2$ and vessel segmentation images, and ${\theta }_A$ and ${\theta }_B$ are the network parameters. Once trained, the networks $A_{{\theta }_A}$ and $B_{{\theta }_B}$ were evaluated on 40 test examples.

Fig. 4

Relative loss curves (${\epsilon }_A$) for the ${\mathrm{sO}}_2$ estimating network.

4.

Results and Discussion

The 3-D image outputs of both the ${\mathrm{sO}}_2$-estimating and segmentation networks were processed in order to calculate the mean ${\mathrm{sO}}_2$ in each major vessel body using only the ${\mathrm{sO}}_2$ estimates from voxels that the segmentation network was confident contained vessels (the reason for using the segmentation output was so that the mean vessel ${\mathrm{sO}}_2$ could be calculated without a priori knowledge of vessel locations that might not be available in an in vivo scenario). More details about how this process was performed are provided in Sec. 3.2. The mean of the absolute difference between the true mean vessel ${\mathrm{sO}}_2$ and the output mean vessel ${\mathrm{sO}}_2$ over all 40 sets of images was 4.4%, and the standard deviation of the absolute difference between the true mean vessel ${\mathrm{sO}}_2$ and the output mean vessel ${\mathrm{sO}}_2$ was 4.5% (some 2-D image slices taken from the networks’ 3-D outputs are shown in Fig. 5, and a plot of all the estimates is provided in Fig. 6). Therefore, on average, the predicted mean vessel ${\mathrm{sO}}_2$ was within 5% of the true value. The mean difference between the true mean vessel ${\mathrm{sO}}_2$ and the output mean vessel ${\mathrm{sO}}_2$ was $-0.3\%$ with a standard deviation of 6.3%. The typical error for a mean vessel ${\mathrm{sO}}_2$ estimate was thus between $-6.6\%$ and 6.0%.

Fig. 5

2-D slices of 3-D network outputs and corresponding ground truth ${\mathrm{sO}}_2$ and vessel segmentation images for two different tissue models (labeled a and b).

Fig. 6

(a) Plot of the output mean vessel ${\mathrm{sO}}_2$ versus the true values for all the vessels in 40 tissue models not used for training, calculated with the voxels belonging to each vessel as determined by the segmentation network output. (b) Plot of the mean ${\mathrm{sO}}_2$ values for the same 40 tissue models calculated using the voxels known to belong to each vessel as determined by the ground truth vessel positions. These plots show that using the output of the segmentation network in combination with the output of the ${\mathrm{sO}}_2$-estimating network significantly improves the accuracy of the estimates.

To assess the effect that using the output of the segmentation network may have had on the accuracy of the ${\mathrm{sO}}_2$ estimates, the mean ${\mathrm{sO}}_2$ of each vascular body in the network output was estimated using the voxels known to belong to each vessel, as opposed to the voxels assigned to each body by the segmentation network output. Curiously, the accuracy of the estimates decreased when the ground truth vessel voxels were used for calculating mean ${\mathrm{sO}}_2$ values. Figure 6 shows a plot of the results over 40 tissue models. The mean of the absolute value of the offset between the true value and the network output was 16.6%, the mean offset was 16.2%, and the standard deviation of the offset was 11.5%. This suggests that regions where the segmentation network confidently classified as belonging to a vessel corresponded to regions where the ${\mathrm{sO}}_2$ network was more accurate. Furthermore, it is clear that the accuracy of the ${\mathrm{sO}}_2$ estimates calculated using the ground truth vessels positions decreases with depth. We do not observe this in the estimates calculated using the output of the segmentation network, suggesting that the use of the segmentation network corrects for the depth dependence of the accuracy of network $A$.

Even though both networks are trained separately, they share the same input data, network architecture, and are trained with the same loss function where only the distribution of values in the corresponding ground-truth varies (continuous versus binary). As such, it is not surprising that the learned mapping properties are similar and complement each other. The L2 loss was used for training network $B$ (as opposed to a binary classification loss function that would normally be used for a segmentation task) to ensure that the network outputs would retain more information about the uncertainty of estimates.

Both networks $A$ and $B$ reflect the limited-view nature of the data in their outputs, hence the positions of the vessels in both differ similarly from the ground truth. The accuracy of the output of both networks decreased with the depth of the vessels, i.e., the distance from the detector array, as can be seen in Fig. 5. There are a couple of reasons for why this might be the case. First, image SNR decreases with depth. The image SNR decreases with depth both because the fluence decays with depth and because of the depth-dependence of the limited-view reconstruction artifacts. Second, these artifacts become more spread in out in space with depth, introducing greater uncertainty as to the shape and location of the vessels. The output of both networks is least accurate in the deepest corners of each image, where the artifacts are the most significant.

Filters in a convolutional layer are the same wherever they are applied in the image, they are spatially invariant, and are, therefore, most suited to detecting features that are also spatially invariant. However, the limited-view artifacts are not; they are small close to the center of the sensor array and become more significant further away. The multiscale nature increases the receptive field and hence locality can be learned by the network. Nevertheless, we decided to limit the receptive field using a slightly smaller network architecture than the classic U-Net. In this way, we retain uncertainty in the deeper tissue layers instead of introducing a learned bias.

The 3-D results shown in Figs. 5 and 6 are of comparable accuracy to results from other groups obtained by training 2-D convolutional neural networks to process 2-D images (lacking the presence of reconstruction artifacts) of simpler 2-D tissue models.^22,24,25 The technique presented here was not only able to handle more complex tissue models (the tissue models presented here feature more realistic vascular architectures and multiple skin layers with varying thicknesses and optical properties), but also took as the input data 3-D images featuring noise and reconstruction artifacts. Unlike networks trained on 2-D images sliced from 3-D images of tissue models (such as those used in Ref. 21), the 3-D networks were able to use information from entire 3-D image volumes to generate estimates. Because the fluence distribution and limited-view artifacts are 3-D in nature, learning 3-D features is more efficient than trying to learn to represent 2-D sections/slices through 3-D objects with 2-D feature maps. This likely increases their ability to produce accurate estimates in more complex tissue models. Despite being more sophisticated than other tissue models used to date, the tissue models used here were nevertheless created with some simplifying assumptions. Each skin layer was assigned a planar geometry, where the value of the optical properties associated with each layer at each wavelength remained constant within each layer (e.g., the scattering coefficient of the epidermis was constant within the epidermis layer). Although the absorption coefficient of each layer was varied for each tissue model, the scattering coefficient of each skin type remained constant (but did vary with wavelength). Other experimental factors that can affect image amplitude, such as the directivity of the acoustic sensors, were not incorporated into the simulation pipeline. It remains to be seen the extent to which these assumptions will hold true when this network is applied to in vivo data. To ensure that networks initially trained on simplified simulated images can output accurate estimates when provided real images, networks may have to be modified with transfer training, taking advantage of datasets of real images.^36,48,49 Looking beyond the complexity of the tissue models, there are other more fundamental challenges that will make the application to living tissue nontrivial. In order to train a network using a supervised learning approach with in vivo data (or even to validate any technique for estimating ${\mathrm{sO}}_2$ in vivo), the corresponding ground truth ${\mathrm{sO}}_2$ distribution must be available. It is unclear as to how this information might be acquired, and this poses a significant challenge that must be overcome to realize or validate the application of the technique. As an intermediate step toward generating in vivo datasets, blood flow phantoms with tuneable ${\mathrm{sO}}_2$ could be used to generate labeled data in conditions mimicking realistic imaging scenarios.^33,34 Although it is important to show that a network can cope with all the confounding effects present in real images of tissue, it is still interesting and important to know that the technique can cope with at least some of the challenges faced in such scenarios. This work provides an essential demonstration of the technique’s ability to generate accurate 3-D estimates from 3-D image data despite the presence of some confounding experimental effects that distort image amplitude, and despite some variation in the distribution of tissue types and the distribution of vessels for each tissue model.

5.

Conclusions

Data-driven approaches have been shown capable of recovering sample optical properties and maps of ${\mathrm{sO}}_2$ from 2-D PA images of fairly simple tissue models. However, because the fluence distribution and limited-view artifacts are 3-D, 2-D networks are at a disadvantage as they must learn to represent 2-D sections/slices through 3-D objects with 2-D feature maps. Networks that can process whole 3-D images with 3-D filters are more efficient as they can detect 3-D features, and this likely increases their ability to produce accurate estimates in more complex tissue models. There may be cases where accurate ${\mathrm{sO}}_2$ maps may only be generated with 3-D network architectures. Therefore, to assess whether data-driven techniques have the potential to provide accurate estimates in realistic imaging scenarios, it is essential to demonstrate a neural network’s ability to process 3-D image data to generate ${\mathrm{sO}}_2$ estimates. The capability of an EDS to generate accurate maps of vessel ${\mathrm{sO}}_2$ and vessel locations from multiwavelength simulated images (containing noise and limited view artifacts) of tissue models featuring optically heterogeneous backgrounds (with varying absorption properties) and realistic vessel architectures was demonstrated. Regions where the segmentation output was confident in its predictions of vessel locations corresponded to more accurate regions in the ${\mathrm{sO}}_2$-estimating network output. As a consequence, the accuracy of the network’s mean vessel ${\mathrm{sO}}_2$ estimates improved when the output of the segmentation network was used to determine vessel locations as opposed to the ground truth. In contrast to both analytical and iterative error-minimization techniques, the networks were able to generate these estimates without total knowledge of each tissues’ constituent chromophores, or an accurate image generation model—both of which would not normally be available in a typical in vivo imaging scenario. This work shows that fully convolutional neural networks can process whole 3-D images of tissues to generate accurate 3-D images of vascular ${\mathrm{sO}}_2$ distributions, and that accurate estimates can be generated despite some degree of variation in the distribution of tissue types, vessels, and the presence of noise and reconstruction artifacts in the data.

6.

Appendix A: Noise Test

Noise was incorporated into the simulated images by adding it to the simulated pressure time series before the reconstruction step. Noise was added to each datapoint in the simulated pressure time series by adding a random number sampled from a Gaussian distribution with a standard deviation of 1% of the maximum value over all time series data generated from the same image, resulting in realistic SNRs of 20.9, 21.3, 21.4, and 21.4 dB for a set of images of a single tissue model simulated at 784, 796, 808, and 820 nm, respectively. The details of this measurement are described in the following section.

Because the SNR depends on the optical properties of objects in the sample domain, the SNR will vary depending on the tissue model used for the test. Here, we only use a single tissue model with a single set of tissue properties to obtain some approximate idea of how much noise features in the simulated images.

7.

Appendix B: Optical Properties of Skin Layers

The refractive index, anisotropy factor, optical absorption coefficient, and optical scattering coefficient of each tissue or chromophore are required to construct a tissue model for a MCXLAB simulation. Here, we tabulate expressions for computing the relevant quantities or list the values of certain quantities for various wavelengths in Table 1. These values/resources were chosen as they featured data in the wavelength range for our simulations.

Table 1

Skin optical properties (λ is given in nm).