Opto-acoustic imaging of relative blood oxygen saturation and total hemoglobin for breast cancer diagnosis

Abstract. Opto-acoustic imaging involves using light to produce sound waves for visualizing blood in biological tissue. By using multiple optical wavelengths, diagnostic images of blood oxygen saturation and total hemoglobin are generated using endogenous optical contrast, without injection of any external contrast agent and without using any ionizing radiation. The technology has been used in recent clinical studies for diagnosis of breast cancer to help distinguish benign from malignant lesions, potentially reducing the need for biopsy through improved diagnostic imaging accuracy. To enable this application, techniques for mapping oxygen saturation differences within tissue are necessary. Using biologically relevant opto-acoustic phantoms, we analyze the ability of an opto-acoustic imaging system to display colorized parametric maps that are generated using a statistical mapping approach. To mimic breast tissue, a material with closely matching properties for optical absorption, optical scattering, acoustic attenuation, and speed of sound is used. The phantoms include two vessels filled with whole blood at oxygen saturation levels determined using a sensor-based approach. A flow system with gas-mixer and membrane oxygenator adjusts the oxygen saturation of each vessel independently. Datasets are collected with an investigational Imagio® breast imaging system. We examine the ability to distinguish vessels as the oxygen saturation level and imaging depth are varied. At depth of 15 mm and hematocrit of 42%, a sufficient level of contrast to distinguish between two 1.6-mm diameter vessels was measured for an oxygen saturation difference of ∼4.6%. In addition, an oxygenated vessel was visible at a depth of 48 mm using an optical wavelength of 1064 nm, and a deoxygenated vessel was visible to a depth of 42 mm with 757 nm. The results provide insight toward using color mapped opto-acoustic images for diagnosing breast cancer.


Introduction
Over the past 25 years, a growing community of researchers, clinicians, and device manufacturers have been investigating and developing translational opto-acoustic (OA) imaging technologies with the aim of improving clinical care for several diseases. [1][2][3][4][5][6] For breast cancer diagnosis, one promising technique is OA imaging combined with ultrasound (OA/US). [7][8][9][10][11][12][13][14][15][16][17][18][19][20] Combining OA with ultrasound is advantageous because each modality provides information that is complementary to the other. Structural and functional information pertaining to total hemoglobin concentration (c tHb ) and blood oxygen saturation (s O 2 ) can be visualized with OA, while morphological and anatomical information can be obtained with ultrasound. In conventional ultrasound, which is commonly used to assess and characterize breast masses, there is large overlap between benign and malignant image features, which results in false-positive cases that must be confirmed with subsequent biopsy. [21][22][23] However, OA can provide functional information about the metabolism of tumors by imaging blood and vascular structure, which may further differentiate benign from malignant masses and potentially lead to improved patient care.
OA involves several physical principles that make it ideally suited for diagnostic breast imaging. Due to the composition of breast tissue, scattered light can penetrate to a depth of several centimeters, which is suitable for clinical applications. Furthermore, using light means that OA safely visualizes the optical contrast of tissue. By selecting specific optical wavelengths, OA can be optimized to endogenously detect certain molecules, which enables visualizing differences in oxygenated and deoxygenated blood. An illustration showing a handheld OA probe emitting light at two different optical wavelengths is shown in Fig. 1. The probe measures acoustic waves, which are generated by blood and other tissue structures when they are illuminated. From the measured waves, OA images are formed that display the functional composition of tissue. Typically, benign lesions have normal blood vessels and normal oxygen saturation. Malignant lesions are metabolically more active and tend to have higher vascularity, irregular blood vessels, and decreased oxygen saturation. [24][25][26][27][28] Since OA directly visualizes c tHb and s O 2 , this can potentially give clinicians the ability to perform more accurate diagnosis. [7][8][9][10][11][12][13][14][15] In OA/US, the probe also acquires conventional ultrasound images, which are coregistered and displayed simultaneously with OA.
To assess imaging performance in OA/US, there remains a need for approaches to characterize relative image maps, capable of describing how statistical color mapping relates to actual c tHb and s O 2 levels. In relative maps, instead of conveying quantitative numerical parameters, which are often difficult to determine, the colorization dynamically adjusts to display high contrast regions from a spatial distribution of values. To characterize this, different levels of relative contrast must be physically reproduced by a phantom that has properties similar to biological tissue. For OA, suitable background materials include gelatin, 45,46 gel-wax, 47 and polyvinyl chloride plastisol (PVCP). [48][49][50][51][52] Previously, a phantom comprising channels filled with bovine blood was described for quantitative OA characterization. 53,54 However, this involved a single flow-path, providing only one s O 2 level, which is insufficient for relative characterization. To deal with this limitation, in earlier work, we presented images where two s O 2 levels were used for validation purposes, 12 but performance characterization was not described.
In this article, an approach for characterizing relative OA/US color maps is presented, which relates RGB image colorization  The central nidus of the lesion is segmented with a cyan color outline. Histopathology conducted after excision revealed this grade III invasive ductal carcinoma (GR3-IDC) had extensive internal necrosis. This may account for the lack of internal OA signal within the central nidus, as vessels are not present in the necrotic region. In OA Total, strong signal (yellow) due to hemoglobin is observed in the boundary zone of the lesion. Strong relative deoxygenation (red) is observed in the boundary zone of the OA Relative map. The image features are consistent with neoangiogenesis and are interpreted here as signs of malignancy.
Journal of Biomedical Optics 121915-2 December 2019 • Vol. 24 (12) to controlled s O 2 and c tHb levels. To quantify performance, we developed a biologically relevant phantom, and a dual-path perfusion flow system, to emulate relative oxygenation differences between two vessels. The two-vessel configuration permits analyzing differences in relative colorization when both vessels are imaged simultaneously. The first objective was to evaluate relative contrast sensitivity (RCS), which is defined as the minimum resolvable difference in s O 2 between two vessels, when using statistically mapped colorization. The initial results were published in our accompanying proceedings paper. 55 Here, we also measure the impact on RCS, when varying the imaging depth and blood hematocrit level. In addition, we analyze penetration depth by simultaneously assessing the contrast-to-background ratio (CBR) of an oxygenated vessel and a deoxygenated vessel. Lastly, we discuss these results in a broader context, where relative s O 2 and c tHb maps are used to distinguish benign from malignant breast lesions.

Theoretical Background
To provide context for our approach, we describe how OA signal is generated from biological tissue and how this affects image formation. In tissue, blood is the dominant optical absorber of near-infrared light. It contains both oxygenated hemoglobin (HbO 2 ) and deoxygenated hemoglobin (HbR), which are molecules that have different wavelength-specific optical absorption. For example, light at 757 nm is absorbed more strongly by HbR compared to HbO 2 , but at 1064 nm this is reversed. When molecules absorb light, they heat up, which causes thermal expansion. The optical absorption coefficient μ a ðλ; xÞ of a material determines the amount of light converted to heat. This is a function of optical wavelength λ, as well as position x ¼ ðx; y; zÞ, since molecular composition varies spatially. A material property called the Grüneisen parameter ΓðxÞ quantifies how the added heat is converted to pressure. The radiant fluence Φðλ; xÞ is the amount of light, at wavelength λ, that reaches position x. Following a rapid heating, an initial excess pressure p 0 ðλ; xÞ is generated, which is equal to E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 6 3 ; 3 4 4 p 0 ðλ; xÞ ¼ μ a ðλ; xÞ × ΓðxÞ × Φðλ; xÞ: (1) The generated pressure propagates acoustically, from each position x to transducers on the tissue surface, delayed according to the speed of sound. The received signal amplitude is proportional to p 0 ðλ; xÞ. Since blood has a much stronger μ a compared to background tissue, the comparatively strong signal from blood leads to high endogenous image contrast. In addition to dependence on λ, the optical absorption of blood depends on its s O 2 level, which represents the percentage of hemoglobin that is oxygenated. 1,56,57 The s O 2 can be determined from the hemoglobin mass concentrations c HbO 2 and c HbR , according to E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 2 a ; 6 3 ; 1 9 2s O 2 ¼ The total hemoglobin concentration is E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 2 b ; 6 3 ; 1 3 8c tHb ¼ c HbO 2 þ c HbR : The optical absorption from biological tissue is described in Appendix A, Sec. 7.1. By considering only HbO 2 and HbR, this can be written in terms of s O 2 and c tHb as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 3 2 6 ; 7 5 2 which is obtained by substituting Eqs. (2a) and (2b) into Eq. (7). The optical extinction coefficients ε HbO 2 ðλÞ and ε HbR ðλÞ describe absorption per unit concentration. Equation (3) illustrates that μ a is proportional to c tHb . Since c tHb depends on the hematocrit level α hct of blood, 55,58,59 the OA signal strength is a function of α hct . The OA signal strength is also a function of depth. In Eq. (1), Φ depends on the system illumination geometry, as well as the wavelength-specific optical properties of the medium, including its optical absorption coefficient μ a ðλ; xÞ, optical scattering coefficient μ s ðλ; xÞ, and anisotropy parameter gðλ; xÞ. [59][60][61] In tissue, Φðλ; xÞ can be approximated based on a wide-beam diffusion approximation using an effective optical attenuation coefficient μ eff ðλÞ. 60 In this case, at depth y and with surface exposure q 0 ðλÞ, the fluence, at x ¼ ðx; y; zÞ, is E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 3 2 6 ; 5 6 7 Φðλ; xÞ ¼ q 0 ðλÞ × e −μ eff ðλÞy : Since Φðλ; xÞ influences the received signal, it can impact the OA image quality as well as the penetration depth. Also, in tissue, it is commonly assumed that ΓðxÞ in Eq. (1) does not vary significantly, compared to Φðλ; xÞ and μ a ðλ; xÞ, so it can be treated as a constant. For OA imaging, estimates of s O 2 and c tHb are computed. Typically, this involves using data recorded by transducers to reconstruct a map of the initial pressure p 0 ðλ; xÞ to solve unknown molecular concentrations. For c HbO 2 and c HbR , this can be computed as described in Appendix A, Sec. 7.2. To compute s O 2 and c tHb , the computed c HbO 2 and c HbR are substituted into Eqs. (2a) and (2b).
Our approach to display relative parametric maps based on s O 2 and c tHb , which uses a statistical mapping procedure, is presented in Sec. 3.2.1. The characterization of these relative maps is described in the next section.

PVCP background material
To characterize the ability to display relative OA/US images, a PVCP-based OA phantom was constructed to mimic breast tissue, as shown in Fig. 3. Breast consists of parenchymal and adipose tissue, which have different optical and acoustical properties that are affected by several factors, including age, menopausal state, and the presence of lesions. 50,62,63 The biologically relevant ranges for breast tissue properties are shown in Table 1. This assumes the range 0.5 to 0.3 cm −1 for optical absorption of near-infrared light; 50 1430 to 1520 m∕s for acoustic speed of sound; 50 and 1 to 25 dB∕cm for acoustic attenuation (at broadband frequency range). 50 The volume fraction of blood f v;blood is ∼1.0% for parenchymal breast tissue, and 4% to 10% for invasive breast cancers. 64 The reduced optical scattering coefficient follows a nonlinear fit of μ 0 59,65 Thus, for optical wavelengths of 757 and 1064 nm, we obtain μ 0 s ð757 nmÞ ¼ 9.4 cm −1 and μ 0 s ð1064 nmÞ ¼ 6.8 cm −1 , with variability of around AE2.5, using a ¼ 6700, which is consistent with other published values. In our PVCP phantom, the optical and acoustical properties were tuned, using a similar formulation to Vogt et al., 50 which included TiO 2 and carbon black. The optical absorption coefficient μ a , anisotropy parameter g, and the reduced scattering coefficient μ 0 s of the phantom were measured using inverse adding and doubling 66 as shown in Table 2.

Vessel configuration
Two vessels of inner diameter 1.6 mm were embedded in the phantom, as shown in Fig. 3(b). The vessels consisted of Teflon tubing and were sloped at a 17-deg incline to permit varying of depth during imaging.

Blood preparation
During testing, the vessels were filled with bovine blood (Innovative Research Inc., Catalog No.: IR1-040N), which has optical properties that are similar to human blood. 67,68 The blood was treated with an anticoagulant (sodium citrate). The hematocrit of the whole blood was analyzed by filling capillary tubes (Fisher, P/N: 22-362-574) with the blood sample and sealing one end with sealing compound (Fisher, P/N: 02-678). Samples were then placed in a microhematocrit centrifuge for 5 min. Using calipers, the length of the red cell volume and the length of total volume were measured. The ratio represents the packed cell volume which corresponds to the hematocrit level. For testing that involved using blood at lowered hematocrit levels, the whole blood was diluted using sterile phosphate-buffered saline (MP Biomedicals, Catalog#: 1860454). The hematocrit of the diluted blood was then remeasured using the same procedure described above.

Dual-path perfusion flow system
A dual-path perfusion flow system was used to independently control the s O 2 of each vessel. This is schematically shown in Fig. 4. The perfusion flow system was driven by a peristaltic pump (Masterflex, model: 75728-10) that pulled blood from a reservoir to a membrane oxygenator (Permselect, model: PDMSXA-2500), to balance the dissolved gas levels of oxygen, nitrogen, and carbon dioxide. The gas supply consisted of a tank of each pure compressed gas. In each path, the tanks were connected to a single gas-mixer (MCO, model: GB 100), to control the flow and combine it into a single output connected to an oxygenator. To determine s O 2 level, the partial O 2 pressure, pH level, and temperature were monitored using sensors. Dissolved CO 2 influences the pH level of blood. This was initially adjusted to establish a pH level of 7.4 using a pH meter (Denver Instruments, model: UB-10). The partial O 2 pressure was measured by a pair of dissolved oxygen (DO) meters (Oakton, model: DO-450), which also record temperature to  AE0.1°C. The meters output DO as a percentage of the calibrated partial pressure of oxygen at atmospheric pressure. As shown in Fig. 4, in each path, the blood flowed from an oxygenator directly to the phantom and then to a DO meter before returning to the reservoir.
To determine known s O 2 levels for blood in both vessels, we established a measurement approach based on DO sensor readings. Oxygen saturation in blood varies according to known hemoglobin dissociation curves. 69 These curves model how s O 2 changes as a function of partial oxygen pressure, pH level, and temperature, and to a lesser extent on partial pressure of CO 2 and DPG (2,3-diphosphoglycerate) concentration. A mathematical relationship to determine oxygen saturation from these variables is described by Dash and Bassingthwaighte. 69 A published Matlab function 69 for computing s O 2 based on the Dash model was used in the present analysis to relate sensor readings to s O 2 values.
An investigation of the measurement accuracy of using hemoglobin dissociation curves was also performed by varying input parameters in the Dash model to compute error-bars for s O 2 measurement. It was determined that the highest accuracy occurs near full deoxygenation and full oxygenation. This is due to the lower slope of the oxygen dissociation curves at the extremes. A greater uncertainty was found for midlevel s O 2 values. This analysis of measurement accuracy is provided in the Supplemental Material.

Statistical color mapping
We performed characterization on relative s O 2 and c tHb maps of the PVCP phantom that were acquired using an investigational Imagio ® system. The maps were generated using a statistical color mapping approach, consisting of two steps. First, statistics were analyzed from a reference region of each map to automatically determine color scaling and offset parameters. Next, using these parameters, a color mapping was applied, to display OA information from areas of strong signal, relative to background tissue. Figure 5 shows how statistical mapping applies a color palette to oxygenation levels for several different scaling and offset parameters.
For determining the mapping, the reference region of each map is selected to include average background tissue, so statistical information about the background can be gathered. The choice of the reference region influences the colorization, since different regions have different statistical properties. Commonly, it is set to encompass the entire midfield tissue region (e.g., 5 to 30 mm depth). However, to analyze certain objects, the reference region can also be set to encompass the upper and lower bounds of a target (or targets) of interest (e.g., the vessels of the phantom), which puts higher weight on the statistics of the targets themselves. Color mapping parameters are obtained by computing a mathematical characteristic applied to s O 2 or c tHb intensity values of the reference region. In particular, the average signal level of reference region is used as a color midpoint, to determine the offset parameter. The color scaling parameter is computed based on the standard deviation of the reference region to maximize the contrast level displayed. As a result, for s O 2 , this qualitatively highlights regions of significance that have higher or lower levels of oxygen saturation than the reference level. In addition, a transparency mask is generated to permit only regions exceeding a specific colorization level to be displayed. Figure 6 shows the RGBA palettes that are used for s O 2 and c tHb . Mathematically, an RGBA color palette defines a colormapping function f RGBA : R → R 4 , by mapping a real value to three color channels (red, green, and blue) and a transparency channel (alpha).
To perform color mapping, let the estimated values for either s O 2 ðxÞ or c tHb ðxÞ be represented by a function fðxÞ.
The estimated values arise from using Eq. (9) with Eq. (2). Since position is represented by x ¼ ðx; y; zÞ, for 2-D imaging, we assume z ¼ 0. Let R represent a reference region with N pixels. The mapping function is offset and scaled by the mean μ and standard deviation σ of fðxÞ in R, where E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 3 2 6 ; 3 7 8 The four-channel color-mapped RGBA image is E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 3 2 6 ; 3 1 6 where Fig. 5 Statistical color mapping is used to maximize functional contrast for a region of tissue. The mapping parameters are automatically computed to display areas of high and low s O 2 , relative to a background reference level. In (a)-(c), the vertical axis illustrates color mapping corresponding to three images with different statistical properties.
(a) (b) To produce RGB output, the RGBA data of Eq. (6) are overlaid with a grayscale ultrasound image by blending of the alpha transparency channel. The function clipðuÞ is used to ensure the values are between k min and k max , which represents the domain of f RGBA . Values below k min or above k max will be saturated at the maximum or minimum colors of the palette. For the color palette in Fig. 6(a), the values are mapped from the domain −1 to 1; for Fig. 6(b), they are mapped from 0 to 1. The parameter k g represents a gain constant.

Display of OA/US maps
To display OA information, Imagio ® generates six individual image maps, arranged as a six-on-one composite image. A representative image is shown in Fig. 7 Figure 7 shows image maps from scanning the biologically relevant PVCP phantom at controlled s O 2 levels.
The appearance of vessels in the maps is described with the results of Sec. 4.1.

Sensitivity to Differences in Oxygen Saturation
To determine RCS, the biologically relevant PVCP-based phantom was imaged using an investigational Imagio ® system. The positioning of the probe relative to the phantom is shown in Fig. 3. The vessels were oriented perpendicular to the imaging plane. A representative image of the configuration is shown in Fig. 7. A reference region with a vertical size of 5 mm encompassed the vessels, as indicated by markers (magenta arrows). Several tests were performed, with the procedure repeated for varied imaging depth, oxygenation conditions, and hematocrit levels, according to Table 3. In each test, image output in six-on-one format (described in Sec. 3.2.2) was continually collected while the s O 2 levels of both vessels were individually controlled. Since the vessels in the phantom are sloped at an incline, the vessel depth was adjusted for each test by an initial positioning of the probe. A clamp was used to keep the probe in a fixed position. Sensor readings of pH, DO, and temperature for each vessel were monitored and recorded as image annotations. Each collected six-on-one image corresponds to a single measurement and was stored as a separate RGB image file (in a lossless PNG format). For each measurement, the sensor readings were used to determine the s O 2 for Vessel 1 and Vessel 2. This corresponds December 2019 • Vol. 24 (12) to the measured oxygen saturation, which was plotted as a function of measurement number. A series of sequential measurements was stored as a dataset. The datasets that were analyzed are listed in Table 3. For Dataset 1A, the results are in Sec. 4.2.
To analyze RCS, a single numerical value ΔM was computed to characterize factors that affect the ability to distinguish small differences oxygen saturation. This was performed for varied depths and hematocrit levels using Datasets 1B and 1C of Table 3. For computing ΔM, the s O 2 of one vessel was varied while holding the other constant. The s O 2 of Vessel 1 was set to full oxygenation as a fixed reference vessel. Vessel 2 s O 2 was initially set to 40% and then gradually increased to full oxygenation. When one vessel remains fully oxygenated, this reduces variability in its s O 2 and facilitates measuring the oxygenation level at which a loss of contrast is reached. Image colorization was measured from the RGB data (see Sec. 3.5). The lowest measurement number was identified where the difference in measured image colorization (as described in Sec. 3.5), between the vessels is less than a contrast distinction threshold ΔC. This was set to an empirically selected threshold value of ΔC ¼ 0.5, where visual distinction is apparent and above the noise floor. The result ΔM was computed by determining, at the measurement number where threshold occurs, the difference in the measured oxygen saturation between the vessels. The procedure was applied to Datasets 1B and 1C, with results in Sec. 4.3.

Penetration Depth
The penetration depth for the biologically relevant phantom was investigated by imaging the vessels at two different s O 2 levels (Dataset 2 from Table 3). The s O 2 of Vessel 1 was set to full oxygenation (∼99% s O 2 ), and Vessel 2 was set to very low deoxygenation (∼5% s O 2 ). The probe was mounted on a linear actuator configured to slide along its elevational axis at a constant speed. Since the vessels are sloped, the acquisition of data at multiple depths could be automated. The geometry of the phantom permitted vessel depth to be varied from ∼14 to 53 mm. While the linear actuator was active, data acquisition was performed at a constant frame rate. The reference region was set to a vertical size of 5 mm and was positioned to encompass the vessels as their depth was varied. The six-on-one composite image maps were collected and stored as individual RGB data files.
The measurements were analyzed in increments of 1 mm. At each depth, the measured image colorization was computed and recorded for the OA Combined and OA Total maps, as described in Sec. 3.5. Contrast-to-noise ratio (CNR) is defined by CNR ¼ C σ , where C represents the contrast of a target, divided by the standard deviation σ of the background. In statistical mapping using Eq. (6), the color mapping already includes division by σ. Therefore, the measured image colorization, which can be directly associated with C, is also a measure of each target's CNR, because the background standard deviation has already been normalized during the color mapping. In addition, CBR measurements were performed for each vessel using the OA Short and OA Long maps. The equation CBR ¼ S−B B was used, where S is the peak signal within the vessel, and B is the average background level. Penetration depth corresponds to the point at which CBR is less than a threshold when the vessel can no longer be visualized. The results are presented in Sec. 4.4.

Image Analysis
In each dataset, RGB image output was produced by the statistical mapping approach. To analyze this output, an inverse mapping was applied to convert RGB values into a scalar value, which corresponds to the amount of relative colorization of each pixel. Since the hue of RGB colors remains unaltered when RGBA is overlaid with grayscale, this was used to determine the inverse mapping for the f RGBA palette in Eq. (6). The palettes have a unique hue for each color, so the inverse mapping f −1 RGBA exists and is not affected by the grayscale wherever the transparency is nonzero. The hue of each pixel was used with f −1 RGBA to determine a real number between k min and k max . OA Relative and OA Combined are mapped from −1 to 1 (from red to green); OA Total is mapped from 0 to 1 (from transparent to yellow). The gray RGB pixels were set to zero, corresponding to transparent areas.
To measure colorization of the vessels in OA Combined, a rectangular window surrounding each vessel was mapped into numerical values, using the inverse mapping described above. The average of this was computed, but only for pixels with mapped values having a magnitude exceeding the median magnitude of the window, which avoided noise by skipping highly transparent pixels. This numerical result, called the measured image colorization, was recorded for each vessel and plotted as a function of measurement number. Figure 7 shows a composite image of the phantom, captured by an investigational Imagio ® system, which consists of the six simultaneously displayed maps described in Sec. 3 Fig. 7, a different intensity for each vessel is observed because the optical absorption coefficient of blood is wavelength specific and depends on oxygen saturation. Vessel 2, which has a low s O 2 (31%), has a low signal in the long wavelength; however, Vessel 2 has strong signal for the short wavelength. Vessel 1, which has high oxygen saturation (99%), has moderately weak signal at the short wavelength and strong signal in the long wavelength.

Varying s O 2 of both vessels
In Fig. 8, the RGB image colorization of the phantom is displayed when s O 2 of both vessels is varied using the dual-path perfusion system described in Sec. 3.1. For this figure, the parameters for Dataset 1A, described in Table 3, are used. The image colorization was computed from the RGB values of the OA Combined map, as described above. To graphically show the RGB colorization of the vessels from the series of measurements, the region surrounding the vessels is arranged according to measurement number [ Fig. 8(b)]. The line style colors plotted in Figs. 8(a) and 8(c) were selected to match the initial RGB measurements in Fig. 8(b). Initially, the s O 2 of Vessel 1 is lower than Vessel 2. The measured colorization of Vessel 1 begins negative (red), while Vessel 2 has a positive (green) colorization. As the s O 2 levels converge, the colorization becomes weaker, and around measurement number 45 this switches colors as Vessel 1 becomes more oxygenated than Vessel 2.

Effect of vessel s O 2 on relative maps
In Fig. 9, some of the OA Combined and OA Total image maps used to generate Fig. 8 are shown at full-size for several s O 2 levels. The oxygen saturation difference Δs O 2 (the difference of measured s O 2 between the two vessels) is indicated for each case. For OA Combined, statistical color mapping causes the strongly oxygenated vessel (relative to the background level) to appear green, and the strongly deoxygenated vessel (relative to the background level) to appear red. Signal that is midway between these levels is transparent. As described in Fig. 5, the colorization adjusts dynamically based on the image content to present regions that have significant s O 2 differences relative to background levels. In the top row of Fig. 9, it is observed that with the relative colorization, the resulting color mapped images remain virtually the same until Δs O 2 is below a background noise level. Below this level, when the vessels are no longer distinguishable, the colorization represents amplified background variability, since the s O 2 levels of both vessels have become nearly the same (Fig. 9, top right). In the bottom row  Journal of Biomedical Optics 121915-8 December 2019 • Vol. 24 (12) of Fig. 9, the corresponding OA Total maps are shown. When s O 2 levels are adjusted, the amount of hemoglobin in the blood does not change. Thus, changing the s O 2 does not have a substantial impact on the yellow colorization displayed here. Since the visibility in OA Total remains high, characterizing the ability to detect small differences applies only for relative s O 2 maps and not for c tHb .

Relative contrast sensitivity measurement
Results from Dataset 1B are shown in Fig. 10, for vessel depth of 15 mm and hematocrit of 42%. Vessel 1 was held at a constant fully oxygenated s O 2 level of 99%, while Vessel 2 was varied from low oxygenation to full oxygenation. The line style colors plotted in Figs. 10(a) and 10(c) match the RGB measurements in Fig. 10(b). Less distinction is observed between the two vessels as their oxygen saturation measurements become closer together. The measurement number at which ΔC ≤ 0.5 first occurs was identified and is indicated in the plots with a vertical blue dotted line. In Fig. 10, this corresponds with measurement 24. The displayed vessel images [ Fig. 10(b)] highlight how this corresponds to the last measurement where sufficient color distinction was detected. At this measurement number, the difference in actual measured oxygen saturation ΔM between the two vessels, which corresponds to RCS, was determined as 4.6%.

RCS measurement
The values for ΔM of Dataset 1B, where depth was varied, are recorded in Table 4(a). The observed trend is that as depth is increased (from 15 to 30 mm), ΔM is increased (from 4.6% to 9.7%). This occurs because the CBR becomes lower as depth is increased, which increases the s O 2 difference required for color distinction. The parameter ΔM is able to assist with characterizing this phenomenon. A similar trend is observed for Dataset 1C in Table 4(b), where the testing was repeated for varied levels of hematocrit. As hematocrit is decreased (from 42% to 10%), ΔM is increased (from 4.6% to 7.8%). Similar to Table 4(a), the ΔM is increased because the CBR becomes lower as hematocrit decreases.

Image colorization versus s O 2 difference
From the datasets, we can also examine the difference in image colorization as a function of Δs O 2 , the difference in oxygen saturation. This is plotted in Fig. 11. The threshold where ΔC ¼ 0.5 is indicated by a horizontal dotted line. The curves, as depth is varied for Dataset 1B, are shown in Fig. 11(a). In Fig. 11(b), the trend is plotted for the three hematocrit levels from Dataset 1C. From these curves, the value of ΔM can be derived by finding the x-coordinate where each curve intersects the dotted line. This also characterizes RCS by illustrating the relationship between Δs O 2 and ΔC and provides information about how ΔM changes with the threshold for ΔC. In addition, Fig. 11 shows that for cases where CBR is higher (15 mm depth, 42% hematocrit), the ability to distinguish the vessels based on image colorization increases. It also shows that the maximum difference in image colorization plateaus beyond a certain level of s O 2 , and this level is greater at high CBR. It is seen that below ΔC ≈ 0.5, the curves have greater variability, which illustrates why this cutoff level was relevant.
The effect of varying depth on CBR is shown visually in Fig. 12. OA Combined and OA Total maps are shown for different imaging depths, when the vessels are held at constant s O 2 levels. At 15 mm depth, the intensity of the background is low. Background intensity visually becomes greater when the vessels are imaged at 23 mm and is substantially increased when the vessels are at 30 mm. The signal received from a vessel becomes weaker as its depth is increased. However, since the reference Journal of Biomedical Optics 121915-9 December 2019 • Vol. 24 (12) region encompasses the vessels, the statistics for color mapping automatically maximize the contrast of the vessels. This means that when the color mapping is targeted to deep vessels, the colorization scales automatically to visualize these weaker signals. As a result, when the vessels are imaged at 30 mm depth, amplification of the near-field background is observed in the image; however, the background for region surrounding the vessels is transparent. This is what is desired for the analysis described above, where the RGB values for the region surrounding the vessels determine the image characterization.

Penetration Depth
The results of penetration depth testing from Dataset 2 are shown in Fig. 13. This is displayed graphically in Fig. 13(a). Figure 13(b) demonstrates that both targets are visible in OA Combined map to depths of around 48 mm. In Fig. 13(c), both targets are visible in OA Total beyond 42 mm, while Vessel 1 remained visible to around 48 mm.

38.
40.  In Figs. 13(d) and 13(e), the CBR measurement was used to determine the penetration depth of the targets at each wavelength. A CBR threshold value of 1.5 was used, which corresponds to when the targets could be distinguished from the background. Below this threshold, it was difficult to visualize the vessels clearly. In Fig. 13(d), for OA Short (757 nm), the oxygenated target maintains CBR > 1.5 up to a depth of 32 mm, while the deoxygenated target maintains CBR > 1.5 up to a depth of 42 mm. In Fig. 13(e), for OA Long (1064 nm), the oxygenated target maintains CBR > 1.5 up to a depth of 48 mm, while the deoxygenated target did not have significant contrast even at the minimum depth of 14 mm. Therefore, the penetration depth for OA Short (757 nm), which strongly detects HbR, is determined to be 42 mm. For OA Long (1064 nm), which strongly detects HbO 2 , the penetration depth was determined to be 48 mm.

Characterization of Relative s O 2 and c tHb with Dual-Path Flow System
We investigated the use of a dual-path perfusion flow system to perform device characterization for statistical color mapping of relative s O 2 and c tHb . Controlling the s O 2 level of two independent vessels permits determining the ability to resolve differences in s O 2 level. In a single-vessel configuration, which was described in previously published approaches, 50,54 the device characterization would be difficult for several reasons: (i) first, the signal intensity and oxygenation level of the background medium are not adjustable. If the relative colorization of a single vessel was compared to the background, the background would serve as a fixed reference. When the s O 2 level of the vessel reached that of the background, the vessel would become transparent according to the color map of Fig. 6. The transparency would be a problem for analyzing the amount of image colorization. (ii) The background produces a relatively weak OA signal, whereas a vessel can produce a strong signal from blood. The CBR of the background is weak (by definition it is equal to zero), whereas for a vessel it is higher. Therefore, using the phantom's background as a reference point for a quantitative s O 2 level, which has an inherently low CBR, would be problematic. (iii) The homogeneity of the phantom background causes it to predominantly emit low-frequency OA signal. This tends to be filtered out by the system transfer-function response, so the DC portion of the received signal is blocked. Therefore, the absence of vascularity makes it difficult to use the phantom's background as a quantitative reference point for s O 2 .
On the other hand, the two-vessel configuration permits the ability to characterize mappings of relative s O 2 and c tHb by examining the difference between each vessel as they are simultaneously imaged. This avoids the issues of the single-vessel configuration described above.

Differences between Phantoms and Breast Tissue
When performing characterization on a phantom, there are several potential limitations in interpreting the results due to differences between phantoms and breast tissue. However, phantom-based results can be used to make some general interpretations about relative color mapping in tissue. When performing characterization, both targets in the phantom were imaged at approximately equal depth, and due to the homogeneous material properties, the targets received nearly equal fluence and had nearly equal acoustic propagation paths for measurement. This represents an ideal scenario for performing controlled characterization, as the independent variables of the blood can be isolated based on the amount of signal generated by either target.
Conversely, biological tissue has increased heterogeneity of optical and acoustical properties compared to the phantom, which can impact signal generation. As a result, tissue may have decreased CBR compared to the phantom. In addition, benign or malignant lesions may occur at different depths and have different hematocrit levels. By definition, the CBR increases with OA signal strength, which is impacted by depth [Eq. (4)] and hemoglobin concentration [Eq. (3)]. Therefore, this situation may lead to reduced CBR during imaging for some lesions. Furthermore, while the phantom only has two targets, acoustic interference from other sources or targets in tissue can lead to decreased CBR. A decrease in CBR may reduce penetration depth. This may also impact relative colorization, since it was demonstrated in Sec. 4.3 that RCS is impacted by the CBR. When CBR is reduced, for example, at increased depth, or at decreased hematocrit level, the ΔM increases (as shown in Fig. 11). In this case, a greater difference in s O 2 between each vessel would be required to maintain distinguishable colorization.
Other factors can also impact the RCS and ability to display relative s O 2 . To illustrate this, consider a perturbation ΔΦ from a modeled fluence distribution Φ. From Eq. (1), the generated OA signal depends on Φ. Variability of optical properties or spatial nonhomogeneity, which may be present in tissue, tends to increase ΔΦ for targets being imaged. If ΔΦ is increased, then ΔM, which represents the difference in actual s O 2 values to distinguish the targets, may also increase. With the phantom characterization apparatus, at either specific wavelength, both targets receive nearly equal fluence, and perturbation from the modeled Φ is assumed to be small. However, in a scenario where targets receive a greater uncompensated ΔΦ, such as in tissue, this may correspond with greater ΔM. In addition, in tissue, the Grüneisen parameter Γ can vary spatially, leading to additional perturbation in signal generation. For tissue, it is assumed that there is small variability in Γ. Spatial variability in Γ can lead to increased ΔM. Likewise, spatial variation in hematocrit level could produce a similar effect.

Visualizing Relative s O 2 and c tHb with Statistical Color Mapping
In statistical color mapping, which performs relative colorization of s O 2 and c tHb (as described in Sec. 3.2.1), the goal is to maximize the contrast of a target or lesion compared to background tissue, rather than to measure precise numerical values of a parameter. Colorization is determined dynamically, to highlight areas that are relatively more significant compared to others. The approach is considered qualitative because numerical measurements of s O 2 and c tHb are not conveyed with the generated RGB maps. One advantage for statistical color mapping is the ability to display relative s O 2 and c tHb when quantitative parameter estimates could have offset and scaling that is difficult to determine in advance. Specifically, the procedure can help compensate for some variability in tissue OA properties and fluence distribution that impact OA signal generation, as described in Sec. 5.2. To dynamically maximize contrast, colorization is determined relative to a background reference region, by computing μ and σ in Eq. (6). This tends to cancel out perturbation to parameters that affect signal generation, when stationary over the reference region. For example, if the assumed parameters and coefficients contributing to Eq. (1) or Eq. (9) are offset or scaled from their actual values (e.g., by stronger signal from one wavelength), then μ and σ of the estimated s O 2 and c tHb quantities will also be offset and scaled, at least to a first approximation. However, Selecting a reference region to visualize deep vessels causes an increase of statistical gain parameter σ, which amplifies near-field background while contrast of targets is maximized. The region surrounding vessels is used for analysis.
Journal of Biomedical Optics 121915-12 December 2019 • Vol. 24 (12) with statistical mapping, any offset and scaling to the distribution of parameter estimates will be dynamically compensated, relative to the background. In other approaches, where quantitative parameter estimates are directly output or mapped to a fixed color, variations in fluence distribution and OA properties may result in the displayed color or output being shifted. For purposes of diagnostic imaging in tissue, the statistical reference region used to compute μ and σ (described in Sec. 3.2.1) is commonly set to encompass the entire mid-field of the image (e.g., from 5 to 30 mm). This is selected to provide and maximize statistical contrast of the image based on properties of the full image region. This is applicable when imaging tissue, which can have more than two vessels as configured in our phantom. By using the statistics of this region, the mapping detects areas of contrast that appear as outliers from background levels, which presents information about the relative spatial distributions of c tHb and s O 2 . In clinical studies, this was demonstrated to be an effective method to visualize breast lesions compared to the background tissue. 7,9,13,14

Application of Relative s O 2 and c tHb for OA/US Clinical Breast Diagnosis
Several clinical studies for breast cancer diagnosis [7][8][9][10][11][12][13][14][15] have examined the ability to distinguish benign from malignant lesions using relative OA/US maps generated with an investigational Imagio ® breast imaging system. For radiological assessment during these studies, the OA/US maps were formatted as a six-on-one composite image, similar to Fig. 7 A distinction between statistical mapping and quantitative approaches is illustrated by the RGB image characterization of Sec. 3.3, where the measured image colorization was used instead of the s O 2 level. In particular, to quantitatively determine s O 2 or c tHb from a statistically mapped RGB image, the values of σ and μ in Eq. (6) must be known (to uniquely relate the color palette to the physical parameter values). However, since σ and μ are internal variables that are not provided as output, this is underdetermined, and values were mapped between k min and k max . In this sense, only the relative spatial information is accessible from the statistically mapped RGB output.

Conclusion
We have presented a method for characterizing colorized OA image maps that display spatial information regarding relative blood oxygen saturation and hemoglobin in tissue. A biologically relevant PVCP-based phantom was constructed to mimic the optical and acoustical properties of breast tissue. The phantom comprised two independently oxygenated vessels filled with bovine blood. A statistical color mapping procedure was used to generate RGB colorization of OA signal. A technique was proposed to characterize the ability to resolve small differences in s O 2 when using maps with relative colorization. In the PVCP phantom, it was measured that distinguishing the two vessels required a s O 2 difference of 4.6%, when imaged at 15 mm depth and at 42% hematocrit. By analyzing the effects of varying vessel depth and hematocrit level, it was determined that a reduction of CBR increases the s O 2 difference required to distinguish vessels. In addition, by analyzing the color mapping approach, we examined how the spatial distributions of relative s O 2 and c tHb may potentially be used during clinical assessment to visualize differences between benign and malignant lesions for breast cancer diagnosis. Future work involves performing further device characterization and developing techniques to improve visualization of diagnostic imaging features.

Optical Absorption of Biological Tissue
In biological tissue, absorption can be modeled as a mixture of HbR, HbO 2 , melanin, and several other less significant chromophores (e.g., water and lipids). 59 At optical wavelength λ, the overall absorption coefficient μ a ðλÞ is a weighted combination of the i'th chromophore's absorption μ a;i ðλÞ times its volume fraction f v;i . For example, E Q -T A R G E T ; t e m p : i n t r a l i n k -; s e c 7 . 1 ; 3 2 6 ; 2 5 3 μ a ðλÞ ¼ f v;HbO 2 × μ a;HbO 2 ðλÞ þ f v;HbR × μ a;HbR ðλÞ þ f v;melanin × μ a;melanin ðλÞ þ f v;water × μ a;water ðλÞ þ f v;lipids × μ a;lipids ðλÞ þ : : : : Equivalently, this can be written as a summation of the mass concentration c i times the optical extinction coefficient ε i ðλÞ. Commonly, the optical absorption coefficient μ a ðλÞ is computed by considering absorption only from hemoglobin, since melanin is present mainly in the epidermal layer, and the other chromophores are assumed negligible. Therefore, at position x, the absorption can be approximated as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 7 ; 6 3 ; 7 5 2 μ a ðλ; xÞ ¼ c HbO 2 ðxÞ × ε HbO 2 ðλÞ þ c HbR ðxÞ × ε HbR ðλÞ: This expression can also be grouped in terms of s O 2 and c tHb . By substituting Eq. (2) into Eq. (7), this results in Eq. (3), which has a linear dependence of μ a ðλÞ on c tHb .

Computation of c HbO 2 and c HbR
To determine c HbO 2 and c HbR from transducer measurements, Eq. (1) must be inverted, which relies on the wavelengthspecific optical absorption of Sec. 7.1. For a tissue sample with known hemoglobin concentrations, the generated OA pressure can be obtained by combining Eq. (1) with Eq. (7). For wavelength λ, and position x, this is For each position x, this can be written more compactly as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 8 ; 6 3 ; 3 7 9 p ¼ Ac: When imaging tissue, the unknown concentrations of c HbO 2 ðxÞ and c HbR ðxÞ, which are represented byc, can be solved by inverting the 2 × 2 matrix for A in Eq. (8). Thus, E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 9 ; 6 3 ; 3 1 5c Here,p represents a map of the initial pressure p 0 ðλ; xÞ, which has been reconstructed by an imaging algorithm using data recorded from transducers, because p from Eq. (8)