2.1.

Experimental Setup

A polyvinyl chloride-plastisol (PVCP) phantom was fabricated with hollow chambers, each with a diameter of 15 mm and a depth of 55 mm. To minimize the variability of ultrasound scattering among chambers due to different air bubble distributions, a single chamber was filled with either a $100\text{-}\mu \mathrm{M}$ aqueous solution of MB from Fisher Scientific (Waltham, Massachusetts), blood (Hb) or a combination of MB and Hb. A 1-mm-diameter optical fiber was inserted in the filled chamber, and the fiber tip was positioned $\sim 20\mathrm{mm}$ below the top surface. The optical fiber was connected to a Phocus Mobile laser (Opotek Inc., Carlsbad, California), transmitting laser light with wavelengths of 710 and 870 nm and with a laser energy of 3 mJ. The resulting photoacoustic signals were received by an Alpinion L3-8 linear array ultrasound probe (Alpinion Medical Systems, Seoul, South Korea) positioned on the lateral wall of the phantom, $\sim 20\mathrm{mm}$ from the hollow chamber cross section. This experimental setup, which is shown in Fig. 1(a), has been described in previous publications.^46,65

Fig. 1

(a) Experimental setup and (b) framework of the dual-wavelength atlas method for mixture estimation of MB and Hb.

Two types of Hb samples were used in this study. First, 23 vials of fresh human Hb were obtained up to two days after blood draw and storage, mixed in a single container, and exposed to air for $\sim 1\mathrm{h}$ to minimize the variability of Hb oxygenation levels. Second, experimental whole porcine blood (Innovative Research, Novi, Michigan) was used on the fourth day of its reported 24-day lifetime.

A total of 11 concentration levels were prepared using different 2-mL mixtures of Hb and MB. These concentrations ranged from 0% to 100% in 10% increments of MB volume percentage. For example, for the 2-mL MB-Hb mixture, a concentration of 60% MB consisted of 1.2 mL of $100\text{-}\mu \mathrm{M}\mathrm{MB}$ and 0.8 ml of Hb. Each concentration mixture was manually stirred with a syringe in a separate container. Unless otherwise stated, the % concentrations reported in this manuscript refer to the % MB concentration. Five trials were conducted for each concentration level by fixing the light-delivering optical fiber at 0 deg, with 20 photoacoustic frames (10 per wavelength) acquired per trial using a frame rate of 5 Hz. It is worth noting that two blocks of 64-channel aperture data were required to obtain an image created from 128 receive elements, which reduced the possible frame rate from 10 Hz (based on the 10-Hz laser pulse repetition frequency) to 5 Hz.

2.2.

Dual-Wavelength Atlas Method for Mixture Estimation

The framework for mixture estimation is shown in Fig. 1(b), based on the dual-wavelength atlas method previously proposed for the identification of individual chromophores.⁴⁶ With this setup, frequency domain information is expected to be useful because the optical fiber is placed inside the PVCP chamber, thus generating fluence maps that diminish radially from the tip of the optical fiber. This fluence distribution generates unequal acoustic frequency responses, as volume regions of different sizes are being excited. Thus, the proposed algorithm is anticipated to leverage frequency domain information to differentiate photoacoustic responses from various concentrations of chromophores.

To implement the proposed approach, conventional delay-and-sum (DAS) beamforming was first employed to create photoacoustic images for each wavelength emission, concentration, and trial. In contrast to preceding work,⁶⁶ we used $M$-weighted short-lag spatial coherence (SLSC)⁶⁷ with a cumulative lag $M=20$ instead of conventional SLSC with $M=5$ to generate the binary masks that segmented signals of interest and provided ground-truth labels. This specific $M$-weighted SLSC cumulative lag value (i.e., $M=20$) was chosen based on the comparison of SLSC⁶⁸ images with $M$-weighted SLSC. First, a representative SLSC⁶⁸ photoacoustic image from the blood dataset was generated with $M=5$. The same photoacoustic frame was then processed with M-weighted SLSC with $M$ values ranging from 10 to 30. Binary masks generated with a $-3\text{-}\mathrm{dB}$ threshold were then created for each image, and the similarity was defined with the Dice coefficient:⁶⁹

The modification from SLSC in our preceding work⁶⁶ to $M$-weighted SLSC in this work was implemented to increase the influence of lower lags over higher lags, thus creating less disjointed binary masks. A single binary mask was obtained per trial, resulting from the logical inclusive “OR” operation of 20 masks (obtained from 10 frames per two laser wavelengths).

As noted in Fig. 1(b) and in our previous work,⁶⁶ for each image, only those pixels included in the coherence masks were used for feature extraction, training, and classification. In-phase and quadrature (IQ) demodulation was implemented with 2.75 MHz and 85% bandwidth, and the resulting analytic spectra from a sliding window of axial kernels was stacked along the frequency axis. For the last step of feature extraction, principal component analysis (PCA) was applied to the power of the stacked spectra to reduce the complexity of the feature space.

The training dataset consisted of a random trial of 0% MB (i.e., 100% Hb) and 100% MB, and the remaining trials and concentration levels were included in the testing datasets. After the dual-wavelength atlas was constructed [i.e., ${\mathrm{PCA}}_{\mathrm{MB}}$ and ${\mathrm{PCA}}_{\mathrm{Hb}}$ in Fig. 1(b)], a mixture estimation algorithm was designed to obtain the concentration level of MB from a testing sample. First, $N=1000$ samples were randomly selected from the MB and Hb atlas. Then, assuming that mixtures are linear combinations of pure MB and Hb concentrations, a concentration distribution $C^{\prime }\in {\mathbb{R}}^{N\times 1}$ was obtained using the following equation:

Because Eq. (2) is applied on frequency data, the three most influential parameters of the dual-wavelength atlas method for the estimation of $C$ are the axial kernel size, threshold of spectral log compression, and number of principal components. An increase in axial kernel size improves the accuracy of the estimation of spectral amplitudes in the Fourier space. The threshold of spectral log compression adds a tolerance to the discrimination between spectral peaks, which affects the variance that is computed in the PCA. The number of principal components controls the amount of information that is reduced to the feature space.

To provide complementary information regarding the location of photoacoustic targets in the imaging plane, the displayed concentration maps were overlayed with an ultrasound image of the plastisol chamber processed with locally weighted SLSC (LW-SLSC),⁷⁰ using a regularization factor $\alpha =1$ and a $1.25\mathrm{mm}\times 1.2\mathrm{mm}$ kernel. The LW-SLSC parameters were obtained from our previous optimization study.⁶⁶ A factor of $\alpha =1$ corresponded to an equal weight between cost and penalty function for calculating the optimal coefficients in LW-SLSC. An axial kernel size of 1.25 corresponded to approximately $4.5\lambda$, where $\lambda$ is the wavelength associated with the center frequency of the L3-8 ultrasound probe.

2.3.

Quantitative Metrics

To evaluate the overall accuracy of the dual-wavelength atlas method for mixture estimation, the estimated concentration levels $C$ were first ordered as a function of the ground-truth labels $k$. The function $f_1:C(k)$ was then compared with $f_2:C=k$, representing a 1:1 relationship between estimated and true concentrations. Finally, the coefficient of determination $R^2$ was used to quantify the performance of the dual-wavelength atlas method:

Similarly, monotonicity was evaluated in two cases: (1) between photoacoustic spectral amplitudes of a specific frequency and ground-truth concentrations and (2) between estimated and ground-truth concentrations. The monotonic trend of these data pairs were quantified with Spearman’s rank correlation coefficient $\rho$, described in the following equation:⁷¹

The mean absolute error (MAE) was used to assess the accuracy of the generated concentration maps:

Finally, processing times were computed for each stage of the dual-wavelength atlas method applied on the human Hb dataset. These stages were (1) generation of coherence masks, (2) generation of acoustic spectra and spectra stacking, (3) PCA, and (4) estimation of MB concentrations, with the last stage including projection to the feature space using the principal component coefficient obtained from the atlas. These processing times were evaluated individually when training our method (i.e., constructing the atlas) and testing our method on a single frame (i.e., estimating a concentration map). The processing times for stages (3) and (4) were only calculated for training and testing, respectively. For training, robustness in the estimation of computation times was achieved by averaging the processing times from 5 different trials of 0% concentration. Similarly, robustness in the estimation of computation times for testing was achieved by averaging the processing times from every frame of the testing dataset detailed in Sec. 2.2. Computation times were measured using the MATLAB (Natick, Massachusetts) environment executed on an Intel Core i5-6600K processor with 32 GB RAM and a TITAN Xp graphical processing unit (GPU). To enhance the computation speed of stage (1) for all results reported in this manuscript, $M$-weighted SLSC was implemented on the GPU described above, following the architecture of preceding work.⁶⁸ In addition, this GPU-$M$-weighted SLSC approach and the previously-reported GPU-SLSC approach⁶⁸ were both implemented to obtain results for Eq. (1) in Sec. 2.2.

3.1.

Concentration Estimations from MB and Human Hb

Figure 2 shows the results from the photoacoustic spectra obtained with mixtures of MB and human Hb. Stacked photoacoustic spectra are shown in Fig. 2(a) with the $y$ axis denoting mixture concentration and the $x$ axis denoting acoustic frequency for optical wavelengths 710 (left) and 870 nm (right). Each concentration block separated by the dashed white lines in Fig. 2(a) shows 20,000 spectra samples that were randomly selected from the total number of training trials, frames, and kernels. Although a subtle frequency shift is observed in the spectra obtained with the 710-nm wavelength when increasing the MB concentration level, the spectral shift across the mixture for the 710-nm response was not strongly monotonic ($\rho =-0.63$). Similarly, no apparent spectral shift was observed for the 870-nm response ($\rho =-0.27$). One possible cause of the frequency shift in the spectra obtained with the 710-nm wavelength is the photoacoustic interaction with residual particles of MB in the PVCP chamber. These particles were unable to be removed as they stained the chamber walls, and they likely added frequency components to the overall spectral response of different chromophore concentrations.

Fig. 2

Evaluation of monotonicity and linearity of photoacoustic spectra obtained from mixtures of MB and Hb. (a) Stacked photoacoustic spectra of several mixture concentration ($y$ axis) when using 710- and 870-nm laser wavelengths ($x$ axis). Each spectrum is normalized and log compressed to a 60-dB dynamic range. (b) Spearman’s $\rho$ coefficient and $R^2$ values of linear regression of concentration levels versus acoustic frequency while varying the log compression dynamic range of the acoustic spectra. (c) Example of an acoustic frequency with a spectral amplitude of concentration levels that follow a monotonic trend ($\rho =-0.88$) and a linear trend ($R^2=0.76$).

Figure 2(b) shows the linearity and monotonicity evaluation results as a function of the acoustic frequency for multiple image dynamic ranges. The $y$ axes denote the $\rho$ (left) and $R^2$ (right) for and linearity and monotonicity evaluations, respectively. The $x$ axis of each plot is divided into acoustic frequencies for optical wavelengths of 710 nm (left) and 870 nm (right). In the acoustic frequency range of 0.5 to 5 MHz for the 870-nm wavelength, as the dynamic range decreases from 40 to 10 dB, the mean ± one standard deviation of improvements in Spearman’s $\rho$ and $R^2$ were $-0.47\pm 0.15$ and $0.45\pm 0.14$, respectively. For dynamic ranges $>40\mathrm{dB}$, the results primarily overlap, and there is less improvement. This overlap likely occurs because the log compression reaches the resolution limit of the acquired photoacoustic amplitudes. However, for the 710-nm wavelength, only the acoustic frequency range from 2.4-5 MHz shows $\rho$ and $R^2$ changes when varying the log compression from 10 to 40 dB. In the acoustic frequency range of 0.5 to 2.4 MHz, results primarily overlap regardless of the chosen dynamic range. This overlap likely occurs because of the strong spectral peak observed from 1.5 to 2.5 MHz for nearly every MB concentration in Fig. 2(a) (i.e., the spectral shift pattern described in the preceding paragraph). This strong spectral peak is robust against log compression, which generates minimal changes to Spearman’s $\rho$ and variation of $R^2$.

Focusing on the $R^2$ results, multiple strong linear trends were observed when results were obtained with 870-nm laser wavelength, yielding a mean $\pm$ one standard deviation $R^2$ of $0.74\pm 0.14$ for frequencies ranging from 0.5 to 5 MHz and thresholds from 10 to 100 dB. In contrast, no strong linearity was observed for the 710-nm region, yielding an $R^2=0.32\pm 0.21$ for the same range of frequencies and thresholds as the 870-nm region. The strong spectral peak observed in Fig. 2 at the 710 nm wavelength yielded a maximum $R^2$ of 0.55 with a dynamic range of 10 dB, whereas the maximum peak observed at 870 nm was 0.81 with a dynamic range of 50 dB. These results suggest that Hb has a stronger influence on the linear relationship of the acoustic frequency versus concentration compared with MB.

Figure 2(c) shows the linear fit for the acoustic frequency of 2.3 MHz and the dynamic range of 60 dB [i.e., parameters that produced the maximum $R^2$ in Fig. 2(b)]. The $y$ axis denotes the log-compressed spectral amplitudes, and the $x$ axis denotes the ground-truth concentrations. For each violin plot, the shape, horizontal line, white circle, and gray box represent the kernel density estimate of the data points, mean, median, and interquartile range, respectively. $R^2=0.78$, suggesting a strong linearity and monotonicity between the spectral amplitudes of photoacoustic signals and ground-truth concentrations. Spearman’s $\rho$ was 0.89, which indicates that there is a high level of monotonicity.

Figure 3 shows the optimization of the dual-wavelength atlas method using the human Hb dataset by maximizing the $R^2$ fit of a 1:1 slope of estimated versus ground-true concentrations. The variation of $R^2$ values as a function of the three most influential parameters of the dual-wavelength atlas method is shown in Fig. 3(a). The z axis denotes the axial kernel size, and the orthogonal axes denotes the number of principal components used and the threshold for spectra log compression, respectively. Although increasing the axial kernel size improved the accuracy of the spectral amplitude generated with fast Fourier transforms, no clear improvement was observed when measuring the overall $R^2$ values. Similarly, parameter sets using two principal components yielded $R^2\ge 0.7$, suggesting a good estimation performance when using more than one principal component. However, no improvement was observed for $\mathrm{PC}\ge 3$. Finally, increasing the threshold from 20 to 60 dB improved the $R^2$ values when other parameters remained fixed, and no additional improvement was observed for dynamic range thresholds $\ge 60\mathrm{dB}$.

Fig. 3

Performance optimization of the dual-wavelength atlas algorithm using the human Hb dataset. (a) $R^2$ values as a function of axial kernel size ($z$ axis), threshold for spectra log compression ($x$ axis), and number of principal components used ($y$ axis). (b) Examples of mixture estimations ($y$ axis) versus ground truth concentrations ($x$ axis) when using the parameter set that yielded the highest $R^2$ value (i.e., 0.80), represented by the dotted square in (a).

Figure 3(b) shows an example of concentration estimation versus ground-truth concentration when using the optimal set of parameters that yielded the greatest $R^2$ value, denoted by the dashed black box in Fig. 3(a). For this optimized estimation result, the $R^2$ value and Spearman’s $\rho$ were 0.80 and 0.89, respectively, indicating a high degree of linearity and monotonicity, respectively.

3.2.

Concentration Estimations from MB and Porcine Hb

Figure 4 shows the estimated concentrations when using porcine Hb. The sensitivity of estimated $R^2$ values while varying two parameters of the dual-wavelength atlas method is shown in Fig. 4(a) with the $y$ axis denoting the threshold for spectra log compression and the $x$ axis denoting the number of principal components used. These $R^2$ values were computed by fixing the axial kernel size to 1.8 mm [i.e., optimal parameter obtained in Fig. 3(a)]. Overall, the effect of varying the number of principal components and the dB threshold implemented prior to estimating the $R^2$ values showed similar trends to those observed in Fig. 3(a) (i.e., decreasing estimation performance with increasing the number of principal components). However, the optimal set of parameters included a threshold of 40 dB instead of the 60 dB threshold obtained in the experiment using human Hb [i.e., Fig. 3(a)], with $R^2=0.86$ and $\rho =0.93$.

Fig. 4

Estimated concentrations from porcine Hb. (a) $R^2$ value of a 1:1 linear fit between estimated and true concentration levels while varying the threshold for spectra log compression ($y$ axis) and the number of principal components used ($x$ axis). (b) Mixture estimation ($y$ axis) versus ground truth concentration ($x$ axis) when using a 40 dB threshold and 1 principal component.

Figure 4(b) shows detailed estimation results when using the optimal set of parameters from Fig. 4(a). The overall performance of the dual-wavelength atlas method when using the porcine Hb was measured by computing the MAE between the vector constructed with the mean of each violin plot and the ground-truth labels, yielding a value of 6.53%. For comparison, the MAE obtained from the optimal parameter set using the human Hb [i.e., Fig. 3(b)] was 10.49%. Thus, the dual-wavelength atlas method generated more accurate estimates of MB concentration when using porcine Hb instead of human Hb, which can be attributed to the uniformity in chemical composition and oxygenation levels of the porcine Hb samples.

Figure 5 shows example DAS photoacoustic images, coherence masks, and an estimated concentration map for a ground-truth concentration of 60% of a single trial. These DAS images were normalized to the maximum amplitude value obtained from both 710- and 870-nm laser wavelength responses and displayed with 35-dB dynamic range to offer a direct comparison across wavelengths. The contours within each image represent the $-3\mathrm{dB}$ boundary of the coherence mask with an area of $5.32{\mathrm{mm}}^2$. It is worth noting that these masks include the low-amplitude regions of the photoacoustic images that appear if minimal signal is present, which highlights the benefit of the coherence mask that is used to determine the presence of coherent signals, regardless of amplitude. The MAE between the concentration map and the ground truth is 9.68%.

Fig. 5

Example of photoacoustic DAS images, coherence masks used for segmentation, and estimated concentration map for a ground-truth concentration of 60% (i.e., 60% MB and 40% experimental porcine Hb). The purple contours represent the masks obtained by thresholding the coherence images at –3 dB and merging the results from 710 and 870 nm.

Table 1 reports the average processing times for each stage of the dual-wavelength atlas method using the porcine Hb dataset. Without considering acquisition times and transferring times of raw photoacoustic data, the construction of the dual-wavelength atlas was completed in $<6\min$, and the mean ± one standard deviation processing time of a single concentration map was $9.2\pm 0.5\mathrm{s}$. This $<10\min$ processing time is relatively fast in comparison with the total procedure times of cardiovascular interventions (e.g., the average procedure time of perfusion coronary interventions is $76\pm 31\min$⁷²), and it can be performed prior to the initiation of the procedure.

Table 1

Processing times (in seconds) for each training and testing stage of the dual-wavelength atlas method.

Figure 6 shows examples of concentration maps obtained from one trial per ground-truth concentration results. The 30% concentration map produced the greatest MAE of 14.80%, whereas the 100% concentration maps produced the lowest MAE of 0.72%. These deviations agree with the overall errors shown in Fig. 4(b), where 30% and 100% ground-truth concentrations reported MAEs of 18.29% and 0.75%, respectively. Overall, the dual-wavelength atlas method achieved a greater accuracy when estimating concentration levels that primarily include MB.

Fig. 6

Examples of concentration maps of MB and porcine Hb generated by the dual-wavelength atlas method in a phantom experiment. The concentration maps (colored) are overlaid on ultrasound images generated with LW-SLSC.

Our initial dual-wavelength atlas method used LW-SLSC beamforming to extract meaningful photoacoustic data.⁴⁶ The method herein uses M-Weighted SLSC beamforming, rather than LW-SLSC beamforming to reduce the extensive processing times required to generate coherence masks for a total of 2200 frames (i.e., 11 concentrations $\times 2$ wavelengths $\times 10$ frames $\times 5$ trials $\times 2$ datasets) of photoacoustic data. However, considering that only a single ultrasound image is needed to show the structural detail surrounding the concentration maps in Figs. 5 and 6, LW-SLSC was employed to beamform the background ultrasound image and improve the contrast, edges, and interpretation of the structural details surrounding the photoacoustic signals.

3.3.

Effect of Mask Size on Concentration Estimation Performance

To illustrate the impact of mask size on estimation performance, coherence masks were generated from M-weighted SLSC images with coherence thresholds ranging from 0.3 to 0.9 in increments of 0.02. These masks were used to segment photoacoustic signals from the human blood dataset and generated concentration maps of different sizes. Then, for each coherence threshold, the MAE was calculated for the all 10 frames per 11 concentration levels.

Figure 7 shows the performance of the dual-wavelength atlas method when using segmented masks generated with varying coherence thresholds. Overall, the MAE decreases when using higher coherence thresholds to generate the segmentation masks. Therefore, we selected a coherence threshold of 0.7 to generate all coherence masks in the experiment results.

Fig. 7

MAE of the estimated chromophore concentration obtained with the dual-wavelength atlas method as a function of $M$-weighted SLSC coherence thresholds.

Figure 8 shows examples of concentration maps with a mixture of 30% MB (top) and 90% MB (bottom) when using coherence thresholds of 0.32, 0.48, and 0.70 (from left to right, respectively). The measured MAE are reported in the lower right corner of each image. For the 30% MB mixture, when increasing the coherence threshold from 0.32 to 0.70, a MAE decrease from 12.06% to 6.23% (i.e., 5.83% decrease) is observed. The reduced error is attributed to the absence of inaccurate estimates near the bottom of the 0.70 concentration map. In contrast, for the 90% MB mixture, when increasing the coherence threshold from 0.32 to 0.70, a MAE decrease of from 5.81% to 5.06% (i.e., 0.75% decrease) is observed, which can be attributed to a more uniform distribution of the mixture across the chamber. Overall, decreasing the mask size of the region of interest with more appropriate thresholding, increases the performance of the dual-wavelength atlas method, as observed in Figs. 7 and 8.

Fig. 8

Example concentration maps generated with different coherence threshold values.