Superpixel spectral unmixing framework for the volumetric assessment of tissue chromophores: A photoacoustic data-driven approach.

The assessment of tissue chromophores at a volumetric scale is vital for an improved diagnosis and treatment of a large number of diseases. Spectral photoacoustic imaging (sPAI) co-registered with high-resolution ultrasound (US) is an innovative technology that has a great potential for clinical translation as it can assess the volumetric distribution of the tissue components. Conventionally, to detect and separate the chromophores from sPAI, an input of the expected tissue absorption spectra is required. However, in pathological conditions, the prediction of the absorption spectra is difficult as it can change with respect to the physiological state. Besides, this conventional approach can also be hampered due to spectral coloring, which is a prominent distortion effect that induces spectral changes at depth. Here, we are proposing a novel data-driven framework that can overcome all these limitations and provide an improved assessment of the tissue chromophores. We have developed a superpixel spectral unmixing (SPAX) approach that can detect the most and less prominent absorber spectra and their volumetric distribution without any user interactions. Within the SPAX framework, we have also implemented an advanced spectral coloring compensation approach by utilizing US image segmentation and Monte Carlo simulations, based on a predefined library of optical properties. The framework has been tested on tissue-mimicking phantoms and also on healthy animals. The obtained results show enhanced specificity and sensitivity for the detection of tissue chromophores. To our knowledge, this is a unique framework that accounts for the spectral coloring and provides automated detection of tissue spectral signatures at a volumetric scale, which can open many possibilities for translational research.

Introduction

The non-invasive assessment of molecular tissue components plays an essential role to facilitate improved medical diagnosis and treatment monitoring. In the past decade, there has been a clear emphasis to monitor physiological and pathological changes in a volumetric vision instead of specific 2-dimensional (2D) planar views of the respective organs [1]. Thus, most of the recent developments in pre-clinical imaging are evolving in the direction of whole-body visualization. This enables the systemic follow-up of disease progression and also the visualization of multiple organs, by promoting the reduction and refinement of laboratory animals [2]. Many established tomographic imaging modalities facilitate whole-body imaging in pre-clinical studies. Although these imaging techniques are widely used, they also have some drawbacks. For example, preclinical MRI suffers from long acquisition time and expensive hardware to attain a homogeneous magnetic field to guarantee high spatial resolution [3], [4]. Other alternatives, such as X-ray, CT [5], [6], and PET [7] use ionizing radiation that can lead to undesired side effects in longitudinal studies [8].

Optical imaging modalities, such as fluorescence [9] and bioluminescence [10], represent a radiation-free option that enables whole-body imaging [11]. Further optical imaging modalities, such as optical coherence tomography [12], [13], diffuse optical tomography [14], [15], and hyperspectral imaging [16], [17] have been investigated in pre-clinical research and also translated towards clinical applications [18], [19], [20]. Although optical imaging technologies are promising in monitoring molecular information due to their high sensitivity, most of these approaches suffer from lower spatial resolution and imaging depth [21].

The hybrid imaging techniques that combine complementary modalities, represent the ideal solution to mitigate these limitations of imaging depth and spatial resolution. Recently, different combinations such as optical-CT, PET-CT, PET-MRI, and Ultrasound-Photoacoustic (US-PA), have been integrated to facilitate multimodal information of the tissues [22]. The combination of US-PA is a hybrid technique as it can simultaneously provide the anatomical, functional, and molecular information of the tissue of interest.

Photoacoustic (PA) is a non-ionizing technology that is based on laser-induced ultrasound waves [23], [24], [25]. In PA, near-infrared (NIR) pulsed laser light is used for tissue excitation. The optical energy absorption induces tissue thermoelastic expansion that leads to the generation of acoustic signals, which can be detected just as conventional ultrasound. Since ultrasound scattering in biological tissues is almost three orders lower than optical scattering, US-PA guarantees higher spatial resolution and imaging depth (up to several ) as compared to pure optical techniques [26].

In photoacoustic the signal generation is proportional to the optical absorption contrast. Thus, illuminating the tissues at different wavelengths enables the detection of molecular tissue constituents. Particularly, since tissue chromophores have distinct spectral absorption signatures, spectroscopic PA imaging is intrinsically sensitive to detect biomarkers, endogenous, and exogenous contrast agents [27]. Recently, spectral photoacoustic imaging (sPAI) has shown potential for a wide range of applications such as tumor theranostics, oxygen saturation imaging, atherosclerotic plaque detection, and the imaging of sentinel lymph nodes [28], [29], [30], [31], [32], [33].

Linear spectral unmixing is a commonly used method to differentiate specific tissue components from sPAI [34]. This is a fitting-based method to unmix the spectra, therefore it requires the user interaction to input the expected tissue chromophores. In the case of disease conditions, these types of supervised spectral unmixing can be challenging, as the actual absorption spectra of the chromophores could differ from the theoretical curves [35].

A highly sensitive spectral unmixing method that can accurately distinguish all the molecular tissue components, without any user interactions, is auspicious. Unsupervised Machine Learning (ML) approaches such as Principal Component Analysis (PCA), Independent Component Analysis (ICA), Vertex Component Analysis, and Non-negative Matrix Factorization (NNMF) are label-free methods that extract characteristic components in a data-driven way [36], [37]. On the other hand, various neural network architectures have been also proposed to automatically unmix oxy- and deoxy- hemoglobin spectra, by utilizing large training datasets of sPAI [38]. Generally, these algorithms follow an iterative optimization procedure that minimizes the cost function. Besides, the strong PA signal from the oxy/deoxy hemoglobin can obscure the presence of less prominent components, such as lipids and collagen, leading to incomplete identification of tissue constituents. Thus, if the majority of the observed pixels contain high absorbing tissue components, there is a significant probability that the less prominent absorbers are discarded in the iterative optimization, as the spectra of these can be overlapped. Hence, these approaches mainly detect the most prominent absorbers, such as melanin and oxy/deoxy hemoglobin [39].

Spectral coloring effect is a known artifact, that plays a crucial role in the sensitivity and specificity of tissue biomarkers detection and thus limits the current translational research applications. This is a corruption effect that alters the spectrum of the incident light along the tissue heterogeneity at depth and consequently leads to misinterpretations of the absorption spectrum shape of the tissue components [40]. Specifically, due to light fluence attenuation along the depth, it is difficult to retrieve the actual absorption fingerprint of the tissue components from sPAI. For example, the same absorber at different depths might show a diverse absorption spectrum as the absorbed light is different. To compensate for the light fluence variations along the depth, pre-defined simulations can be used, but this approach is also complex as it requires detailed information about the tissue to be imaged [41]. Recently, Tzoumas et al. [42] have proposed an eigenspectra-based fluence correction approach and automatically estimated the blood oxygen saturation in deep tissue. Their proposed approach is based on fluence modeling through the eigenspectra analysis which is mainly used to compensate for . Hence, it is still challenging to automatically and accurately extract multiple molecular components, such as prominent and weaker tissue absorbers, from sPAI.

In this work, we developed a novel data-driven superpixel PA unmixing (SPAX) framework to enable the differentiation of molecular tissue components, without any user interactions. The framework is also extended to compensate for the spectral coloring artifact. To compensate for the spectral coloring we are proposing an automated US image segmentation and spectral Monte Carlo (MC) light fluence simulations based on a predefined library of tissue optical properties. The approach is also optimized for the volumetric assessments of tissue composition. Besides, it includes an optimized visualization of the molecular components’ distribution with unprecedented details. This is a unique data processing procedure, that differentiates the SPAX framework from the other spectral unmixing algorithms. We benchmarked the sensitivity and accuracy of the proposed approach on tissue-mimicking phantoms, volumetric tissue composition, and also on whole-body animal imaging.

Material and methods

Theoretical background

In photoacoustic (PA), nanosecond pulsed laser light is used as an excitation source, causing thermoelastic expansion of the tissues. The corresponding initial pressure distribution leads to the generation of acoustic waves that propagate toward the ultrasound (US) detectors. Therefore, the reconstructed PA image is a spatial representation of the initial pressure distribution defined as follows:

where represents the spatial coordinates of the voxel, is the wavelength, is the optical absorption coefficient, is the scattering coefficient, is the optical anisotropy factor, is the Grneisen parameter, and is the light fluence. As established by Eq. (1) the initial pressure distribution depends on the wavelength of excitation and the fluence . During spectroscopic photoacoustic imaging (sPAI), the intensity of each pixel at different wavelengths is proportional to the absorption spectrum of a combination of tissue components , where and is the total number of wavelengths of the excitation light. This is a low-rank approximation of the sPAI. Hence, the pixels intensity of the sPAI and consequently the underlying initial pressure distribution can be represented as:

where is the number of distinct dominant absorbers, represent the concentrations of the different components, and are the molar extinction coefficients of the components, , where .

In a Linear Mixture Model (LMM), the sPAI can be formalized as:

where is the mixture matrix of observations and variables. contains per each column the vectored PA image at a specific wavelength. The matrices and contain the abundance maps and the absorption spectra of source components, respectively.

As described by Li et al. [43] the conventional approach to solve Eq. (3) is to input the source spectra as a priori information. Besides, blind source separation algorithms [44] can also be used to iteratively solve Eq. (3), where both and are unknown variables.

Superpixel Photoacoustic Unmixing (SPAX) framework

Fig. 1 shows the flowchart of the SPAX framework. Each procedure within the framework has been implemented and refined in order to minimize the user interaction. A detailed description of the SPAX framework procedures is provided below.

Fig. 1

Flowchart of the Superpixel Photoacoustic Unmixing (SPAX) Framework. The numbers represent the sequence of procedures implemented within the SPAX framework. Symbols legend: : light wavelength; : absorption coefficient; : scattering coefficient; : anisotropy coefficient; : refractive index; : initial pressure; : spatial coordinates; : light fluence; : Grneisen parameter; : mixture matrix; : number of observations; : number of variables; : abundance maps;: source spectra; : left-eigenvectors; : eigenvalues; : right-eigenvectors; : number of endmembers; : cross-correlation; : mixed subset; : number of mixed subsets; : source spectra of .

Input data

The high-resolution ultrasound (US) and multi-spectral photoacoustic images (sPAI), within the wavelength range of , represent the input of the proposed SPAX framework.

Automated segmentation

An automated segmentation approach that differentiates various tissue structures has been implemented to guide the light fluence compensation. Within this approach by utilizing the high-resolution US images, the skin line, tissue structures, and coupling medium can be detected. Specifically, the skin-line boundary was automatically obtained by determining a threshold from the ditch of the gradient histogram of the US image. After thresholding the US image, the Sobel filter is applied as edge detection method for finding the tissue boundaries. Thus, the detected skin line is then used as a watershed to distinguish the tissue structures and the background. Besides, in the case of phantom imaging, circular inclusions can be recognized by using the Hough Transform, which results in fully-automated segmentation.

Once segmented, the images are labeled and tagged with optical properties gathered from literature in terms of absorption and scattering [45], [46], [47], [48]. Thus, each cluster is linked with the spectral optical properties collected in a predefined library. This predefined collection of spectral optical properties of the tissue components is provided in detail in the Supplementary Materials (see Supplementary Materials; SupMatTable). The tissue optical properties values have been retrieved from literature at each wavelength of the sPAI acquisition. For exogenous agents, the absorption characteristics have been measured via spectrophotometer and the scattering properties were evaluated according to the concentration of the pigments by following the Mie scattering theory [49]:

where denotes the scattering coefficient in , is the wavelength of the light in , is the concentration in , is the reduced scattering coefficient, and is the anisotropy factor approximated by Eq. (4).

Fluence compensation

After the segmentation, the light fluence correction step has been incorporated into the SPAX framework. The light fluence distribution was simulated with MCXLAB [50] at each wavelength of light used during the sPAI acquisition. Fluence simulations were run with photons, all launched isotropically at each wavelength, and voxel linear length of . Although in the experiments we are using a double-side illumination as a light source, the two beams are converging, and thus guarantee a uniform illumination at the surface. Since during the US-PA image acquisition, we are normally positioning the surface of the sample/animal close to this converging point, we could approximate the light source as a truncated Gaussian beam in the simulations. Thus, the excitation light source was set as a truncated Gaussian beam with a waist radius of voxels, with its center on the top of the skin line. The fluence map is obtained by integrating the flux output from MCXLAB over time, using time steps of for a total range of . The light fluence distribution has been simulated in 2D per each slice of the volumetric domain. Finally, since the same energy has been used at each wavelength to attain the fluence in each voxel, this is normalized to the total energy from the light source. Thus, the obtained fluence maps were normalized and the results have been interpolated to render and match the volumetric imaging domain.

Automatic estimation of the hyperparameter

Singular Value Decomposition (SVD) is implemented to estimate the number of significant components, named hyperparameter , by inspecting the mixed sPAI data along orthogonal directions. The automatic tune of the rank hyperparameter is crucial to have an unbiased characterization of the tissue composition. In particular, the SVD method enables to determine and sort the directions along which the variables exhibit higher variation, which represents the low-rank approximation of the original data [51]. To estimate the number of significant components, the so-called elbow method has been applied [52]. Thus, the mixed sPAI matrix , which is also described as a linear mixture in Eq. (3), is expressed as , where and are matrices of the left- and right-eigenvectors respectively. is the diagonal matrix of eigenvalues in decreasing order. The gradient operation has been applied to the eigenvalues to determine the “elbow”. By interpolating the eigenvalues after the elbow with linear fitting, the noise level is defined. We selected the SNR as a watershed to automatically distinguish the relevant information above the defined noise level. Consequently, all the eigenvalues that have SNR higher than are selected and thus it defines the actual number of endmembers .

Unsupervised unmixing of prominent absorbers

The unsupervised unmixing of the prominent absorbers is based on a modified and optimized version of Non-negative Matrix Factorization (NNMF). This approach recognizes an initial set of underlying highly absorbing molecular components, such as oxy/deoxy hemoglobin. Thus, source spectra () and the respective abundance maps () are obtained. The details of the optimized NNMF algorithm are described in our previous work [53].

Superpixel subsampling

After the unmixing of the prominent absorbers, to differentiate the less prominent molecular components, a novel spectral superpixel subsampling (SS) approach has been implemented and integrated within the SPAX framework. In this procedure, we further analyze the spectra and the respective distribution pixels , where there is an overlapping of the prominent absorbers. A detailed workflow of the superpixel subsampling procedure is provided and described below.

SS-Step1: Cross-correlation of the source components. To select the overlapping wavelengths and pixels a cross-correlation is performed between all the prominent detected spectra , obtained from 2.2.5. The cross-correlation matrix , where , is symmetric. The correlation coefficients can range from to , with representing no correlation, and representing a direct negative or positive correlation respectively.

SS-Step 2: Subsampling. Each matching pair of spectra () with a positive correlation () is analyzed in detail and the wavelengths in which the components have similar intensity are selected for further analysis.

In addition to the spectral analysis, the abundance maps () obtained from 2.2.5 of the positively correlated components were also further analyzed. A pixel-to-pixel intensity comparison is performed between the respective abundance maps. A kernel of dimension three around each pixel is used to determine the pixels’ average used as a threshold. The pixels with similar intensities, which are above the threshold have been selected. In this way, subgroups of meaningful wavelengths and pixels are obtained for further analysis.

In particular, by using the subgroup of wavelengths and pixels, a reduced subset () is created from the original raw dataset (). Where and represents the number of positively correlated components.

SS-Step 3: Subsets Unmixing. To detect less prominent absorbers, per each subset , the unmixing procedures 2.2.4 and 2.2.5 are reiterated till the process converges. Thus, per each subset , after SVD and NNMF a subgroup of spectral components is obtained. This procedure is repeated till all the subsets have been analyzed and thus the condition is verified.

SS-Step 4: Interpolation & Clustering. The spectral components , obtained from each reduced subset , are linearly interpolated to have the original number of variables/wavelengths. Finally, k-means clustering approach was used to classify all the detected source components and resized , obtained from the original mixture and all the subsets respectively.

Further details of the superpixel subsampling data processing, with an example, are provided in Section 3.2.

Output data

Finally, the tissue spectral components are extracted in an unsupervised way from sPAI. A comparison with the reference absorption spectra known from the literature [47] has been performed. Thus, some of the endogenous tissue chromophores, such as oxy/deoxy hemoglobin, and fat can be recognized. Especially in the case of disease models, there is also a good chance to detect additional spectra that are unknown and that cannot be compared with spectra from the literature used as a reference. The increased number of automatically extracted spectra demonstrates the higher sensitivity of the proposed approach to reveal any spectral feature.

Besides the spectral components, the respective distribution maps have been also obtained and refined to enhance the vessel structures by using the Frangi filter [54].

Data acquisition

High-resolution Ultrasound (US) and spectral Photoacoustic imaging (sPAI) have been acquired by using the platform Vevo LAZR-X (FUJIFILM VisualSonics, Inc., Toronto, ON, Canada). The imaging setup includes a high-frequency US system (Vevo 3100), an optical parametric oscillator (OPO) integrated Nd:YAG nanosecond pulsed laser, and the animal imaging platform. The US system is equipped with a linear US transducer array (MX201) that consists of elements at a nominal center frequency of and bandwidth of . This transducer guarantees an axial and lateral resolution of and respectively. Light from the laser is delivered to the tissue through optical fibers, mounted on either side of the transducer. To obtain the homogenous light illumination, the sample to be imaged is placed on the converging area of the two light beams. The repetition rate of the laser is , and the spectral photoacoustic images can be acquired by tuning the laser wavelengths within the range of . The system allows the selection of any wavelength within this NIR range. In particular, for the experiments, we have used a step size of between the wavelengths to acquire the spectral PA images.

During volumetric US-PA acquisitions, a stepper motor is used for the linear translation of the US transducer and optical fibers along the sample. The linear stepper motor moves in steps of a minimum of while capturing 2-D parallel images, for a maximum 3D range distance of . The voxel size used for the reconstructed US and PA imaging was . A modified animal setup (see Supplementary Materials; SupMatFig) has been used for the whole-body sPAI acquisition. A small water container with an optically and acoustically transparent thin polyurethane membrane was used to guarantee the optimal acoustic coupling during the acquisitions. The water temperature was kept at . During the whole-body acquisitions, the stepper motor translated along the animal body from cranial to caudal in steps of . The configuration in epi-illumination, the ECG monitoring, and the respiratory gating have been maintained during all the in vivo whole-body experiments.

Experimental validation

The SPAX framework has been tested and validated by using tissue-mimicking phantoms and in vivo experiments.

Agarose phantom

The spectral coloring compensation approach included within the SPAX framework has been tested on a phantom made from agarose. The schematic illustration of the phantom is depicted in Fig. 2 (A). This includes two absorbers secured into a 3D printed chamber, fixed at and from the transducer. The inclusions have been fabricated by mixing native gel wax (FF1 003, Mindsets Online, Waltham Cross, UK) with of black oil color (Winsor & Newton, London, UK) as described elsewhere [55]. The bulk material of the phantom is made of Agar (Alfa Aesar, Heysham, Lancaster), mixed with Intralipid (IL) (20%, Sigma-Aldrich, Canada) and Verdye (Diagnostic Green GmbH, Aschheim, Germany). A dilution of Verdye and IL have been used to achieve the absorption as and the scattering as of the bulk material (Agar-IL-Verdye).

Fig. 2

Schematic of the agarose phantom (A); Absorption spectra of Verdye and Black Ink measured by spectrophotometer (B).

The absorption spectra of Verdye and Black ink were optically characterized with the spectrophotometer and the results were shown in Fig. 2 (B). The black ink used for the inclusions shows a broad absorption spectrum while the Verdye included in the bulk material has a characteristic spectrum with an absorption peak at . The spectral optical properties of the phantom absorber materials are reported in Supplementary Materials (see Supplementary Materials; SupMatTable) and these were used to set the Monte Carlo fluence simulations. In addition, water has been used as a coupling material to fill depth between the solid phantom and the US transducer.

Tissue phantom

A tissue phantom is used to evaluate the detection of multi-spectral components by SPAX framework. The phantom combines endogenous muscle-fat tissue layers and exogenous inclusions with different absorption characteristics. The advantage of this approach is that this tissue phantom offers good approximations of realistic samples. A schematic illustration of the phantom is depicted in Fig. 3 (A). The phantom was prepared by using a muscle layer from veal meat of thickness. The tissue was selected to contain pure and untreated muscle tissue and not blood. Three cylindrical absorbers have been embedded within the muscle tissue. The cylindrical inclusions have been prepared by mixing gel wax (FF1 003, Mindsets Online, Waltham Cross, UK) with three different oil colors (Winsor & Newton, London, UK). Specifically, of black, green, and blue oil colors have been mixed with melted gel wax respectively. Finally, these have been poured into 3D printed cylindrical curved molds, used to shape the three inclusions (see Fig. 3 (A)). The absorption spectra of the oil colors have been measured via spectrophotometer and reported in Fig. 3 (B). Besides, a thick fat layer, from porcine lard, has been embedded on top of the muscle layer. Thus, the tissue phantom is composed of multiple components: the cylindrical inclusions that are prominent absorbers and the fat layer that has less prominent absorption characteristics. Water has been used as a coupling medium during the US-PA imaging acquisition to fill depth between the phantom and the US transducer.

Fig. 3

Schematic of the muscle-fat tissue phantom (A); Absorption spectra of Blue, Green, and Black oil colors measured by spectrophotometer (B).

In vivo experiments

The in vivo experiments were conducted at FUJIFILM Sonosite/VisualSonics facility in Amsterdam, The Netherlands. All the experiments involving animals were in full compliance with the protocol (AV D2450020173644) evaluated and approved by the Animal Use and Ethics Committee (CEUA) of The Netherlands. These were in accordance with FELASA guidelines and the National Law for Laboratory Animal Experimentation (Law No.18.611). For the hindlimb and whole-body imaging, a female CD-1 mouse model (Envigo, Horst, the Netherlands) and athymic nude-Foxn1nu mice (Envigo, Horst, the Netherlands) were used respectively. Mice were anesthetized with isoflurane ( by volume with gas flow). The US and sPAI acquisition of the hindlimb region was performed with the animal in its supine position, while the whole-body imaging have been acquired from supine and prone positions with the transducer aligned perpendicularly to image the region of interest.

Results

This section presents the results of the SPAX framework on the tissue-mimicking phantoms and also on the volumetric animal imaging.

Spectral coloring compensation

The agarose phantom has been used to validate the spectral coloring compensation procedure included within the SPAX framework. Besides, the unmixing results obtained with and without the spectral fluence compensation have been evaluated. Fig. 4 (A) shows the US-PA image of the agarose phantom in cross-section at , where the inclusions at two depths are visible. The sPAI of the agarose phantom has been acquired within the range of . This wavelength range has been selected in accordance with the known spectral characteristics of the absorbers included within the agarose phantom, thus reducing the acquisition time.

Fig. 4

US-PA image of the agarose phantom cross-section at , including the manually drawn ROIs around the inclusions at two depths ( and ) (A); Averaged PA intensity of the pixels within the ROIs of and (B).

The two inclusions with black ink absorber (named and ) show higher PA signal intensity than the surrounding. Since some fragments of black oil color may remain suspended in the gel wax, the US-PA signal within the inclusions could appear not homogenous. Fig. 4 (B) shows the PA intensity profiles from the absorbers in the wavelength range of . By following the conventional approach used to evaluate the spectral PA intensity, the shape of the inclusions and the respective spectra were obtained by manually drawing the ROI around the inclusions. The spectral measurements have been performed over multiple slices of the phantom. To highlight the differences in the spectral shapes of the two inclusions, due to spectral coloring along depth, the spectra were normalized to visually enhance this effect. Specifically, the spectral data are shown in the graph as the normalized average of the pixels’ intensity within the ROIs and the standard deviation (shown as the lighter color band along the plotline).

Although the inclusions are made of the same material, the obtained spectra show a different spectral trend. In particular, the PA spectrum from the superficial inclusion () showed a broad trend that is more in accordance with the measured absorption spectrum of the black ink (Fig. 2 (B)). The spectrum from the deeper inclusion () has shown a lower intensity within the range of , while it follows the expected trend above . These differences in spectra at depth are a clear indication of spectral coloring artifact, as the incident light attenuates along the depth and thus changes the spectral shape.

Fig. 5 (A) shows the automated segmentation obtained from the US image by using the SPAX framework. From the segmented image, it is clear that the phantom, two inclusions, and the coupling medium are accurately outlined. Based on the segmentation mask and the defined optical properties (see Supplementary Materials; SupMatTable), spectral Monte Carlo simulations have been implemented ranging from to . Fig. 5 (B) represents the fluence map at a single wavelength (), that qualitatively models the spatial distribution of light energy along the geometry.

Fig. 5

Segmented mask obtained from the high-resolution US image (A); Light fluence distribution map obtained from the MCXLAB simulation at (B).

Fig. 6 shows two different scenarios where the automatically unmixed components were obtained with and without applying the light fluence compensation.

Fig. 6

Automatically unmixed spectra without the spectral fluence compensation (A-C); Automatically unmixed spectra including the spectral fluence compensation (D-F); The reference spectra of Black Ink* and Verdye* are included as dashed lines (D-E).

The spectra obtained without applying the fluence correction are shown in Fig. 6 (A-C), where four components have been automatically detected. Fig. 6 (A) shows two spectra that are similar to the manually obtained spectra from the inclusions (Fig. 4 (B)). In addition to the spectra of the inclusions, the SPAX framework has also detected two less prominent spectra that could represent the Verdye and the background (Fig. 6 (B-C)).

Fig. 6 (D-F) show the components detected by the SPAX framework, including the segmentation and fluence compensation approach. Unlike the previous scenario, here only three spectral components have been automatically identified. Most importantly, by applying the fluence correction, both inclusions have been identified as one component with a broad absorption spectrum as shown in Fig. 6 (D). The detected spectrum from both the inclusions is in accordance with the measured spectrum of the black ink obtained by using the spectrophotometer (depicted as a dashed line). This is clear evidence of overcoming the spectral coloring artifact by using the SPAX framework. In addition, improved spectra of Verdye and background have been detected and shown in Fig. 6 (E-F).

Fig. 6 emphasizes the effect of the spectral coloring and the compensation of the SPAX framework. Besides, the correlation values between the spectra measured by the spectrophotometer and unmixed components obtained by the SPAX framework have been evaluated and reported in Table 1. Without the fluence compensation, two spectra were detected for , and these have a correlation value of , respectively, with the black ink spectrum. On contrary, when including the fluence correction, the algorithm has detected only one spectral component from both the inclusions, which has a correlation value of . This is in accordance with the expectation from the fabricated phantom, where both inclusions are made of the same component.

Table 1

Correlation Values of the Agarose Phantom Components.

Component	Without Correction	With Correction
Black Ink	0.87 (Inclusion_1); 0.22 (Inclusion_2)	0.91
Verdye	0.92	0.97

Furthermore, the correlation with the measured Verdye absorption spectrum has also been evaluated before and after the correction, and an improved value from to has been obtained. The Verdye (in Fig. 6 (E)) shows a slight shift of the absorption peak as compared with the Verdye spectrum measured via spectrophotometer (depicted as a dashed line). Since the Verdye absorption spectrum depends on the concentration of the dye, even small changes in the concentration may cause a change in the spectral shape.

Besides, the automatically obtained background spectrum could be related to the Intalipid, since the corresponding abundance map is mainly distributed in the surroundings of the inclusions.

Superpixel multi-component unmixing

The spectral PA images of the tissue phantom have been used to validate the superpixel spectral unmixing procedures of the SPAX framework. Fig. 7 (A) depicts the cross-sectional view of the muscle-fat tissue phantom at , where the PA (red scale) is overlaid on the US image (grayscale). Fig. 7 (B) shows all the observed spectra from sPAI, after US segmentation and fluence compensation, which are defined as the original mixture matrix ().

Fig. 7

US and PA image at 700 nm of the tissue phantom in cross-section (A); PA spectra of the original mixture matrix (B); Eigenvalues analysis for the automatic estimation of the hyperparameter , where 6 components resulted above the noise floor (C); Source spectra obtained from the unsupervised unmixing of prominent absorbers (D).

The SVD analysis is performed on these data. This analysis is crucial to automatically determine the hyperparameter () which represents the number of eigenvalues () that are significantly above the noise floor. Fig. 7 (C) shows the eigenvalues obtained from the tissue phantom sPAI dataset, where shown in red are automatically selected as the hyperparameter.

After the identification of the hyperparameter, the unsupervised spectral unmixing procedure is performed on the matrix , and the result is shown in Fig. 7 (D), where the source spectra are depicted. These obtained components may not be fully unmixed or exclude some weaker absorbers. Therefore, after the detection of these components, an advanced superpixel subsampling procedure is implemented to refine the unmixed spectra. Thus, the superpixel subsampling (SS) approach is applied by following the steps as described in Section 2.2.6.

SS-step1: Within the superpixel subsampling procedure the cross-correlation values between each pair of these components are evaluated. Fig. 8 (A) shows the cross-correlation matrix with the respective values obtained per each comparison. Hence, the pairs of spectra which have a similar trend will have a positive correlation, while the spectra which are easily distinguishable will have negative correlation values. To improve the specificity and also to detect the less prominent absorbers, the pairs of spectra with positive correlation values are selected for further analysis.

Fig. 8

Cross-correlation matrix of the detected spectral components (A); and are an example of positively correlated spectra, which are further analyzed to select a subgroup of wavelengths (B). Overlapped cross-sectional abundance maps of and used to select a subgroup of pixels (C); The scale bar size is ; The subset of spectra obtained from the superpixel subsampling of and (D).

SS-step2: Specifically, per each pair of positively correlated spectra, a subset of mixed data would be created by selecting specific wavelengths and pixels. As an example, from the cross-correlation matrix, and show a positive correlation value of . For the wavelengths selection, the corresponding spectral intensities of and at each wavelength are shown in the scatter plot in Fig. 8 (B). Thus, the coordinates of each point of the graph represent the intensities of the pair of spectra at a specific wavelength. Overall, the graph shows the positively correlated spectra at all the wavelengths, and the red lines depict the mean intensity values of each component of the considered pair. From this process, a subgroup of meaningful wavelengths, those that have comparable intensities can be selected. Generally, these would be located in the quadrants and . In the example, in quadrant both and result above their mean intensity values (red lines). In quadrant both components have intensity below their mean values. On contrary, the wavelengths that appear in quadrants and are excluded as their intensities are distinct, and thus directly distinguishable. Specifically, in quadrants and , one component results below and the other is above their mean values and vice-versa.

After the selection of wavelengths, a similar analysis is performed to obtain the subgroup of pixels that still need to be distinguished. To this end, the abundance maps that display the spatial distribution of the components are further evaluated. The maps of the positively correlated components are used for the pixel-to-pixel intensity analysis, where the pixels with comparable intensities are selected. For the example, Fig. 8 (C) shows a composite image where the abundance maps of and are overlaid. White pixels in the composite image show where the two maps have the same intensities. Magenta and green regions show where the intensities are different. Thus, this graphically shows that the white pixels are selected for further reiterations, since the two positively correlated components are spatially overlapped in there.

Overall this process of selecting the subgroups of wavelengths and pixels is called superpixel subsampling. The same analysis is repeated for all the positively correlated pairs of components and as a result, multiple subsets with a reduced number of wavelengths and pixels are obtained.

An example, that represents a subset obtained from the original sPAI raw data is shown in Fig. 8 (D). This represents one of the reduced subsets of observed spectra . In the graph, it is evident that in comparison to the original mixture matrix shown in Fig. 7 (B), this subset includes a smaller number of variables and observations. Specifically, the wavelengths around that are characteristics of the strong black absorber have been excluded. This example shows that the framework adapts the wavelengths and pixels selection to highlight the spectra that are generally obscured by the highly absorbing components.

SS-step3: Per each subset the same procedures (2.2.4 and 2.2.5) are repeated, as for the raw matrix . The SVD and unsupervised unmixing are performed till it converges to detect the meaningful components for each subset.

SS-step4: All the spectral components , obtained from each subset , are resized to the same original number of wavelengths . Finally, k-means approach is used to classify the spectral signatures obtained from all the subset and the initial mixture . Hence, the superpixel subsampling process increases the detection of less prominent tissue chromophores, those are generally obscured from the most absorbing chromophores.

Fig. 9 (A) shows the spectral components obtained as output of the SPAX framework from the tissue phantom. In this case, in total spectral components are obtained. Specifically, after superpixel subsampling chromophores have been distinguished as black, green, blue oil colors, and fat. The additional components could be related to water and background tissue. Fig. 9 (B) shows the distribution maps of the unmixed components in the 2D cross-sectional view. The fat is the component depicted in yellow that is present as a top layer. Under the fat layer, the cylindrical inclusions made of highly absorbing dyes are depicted in blue, green, and red. Fig. 9 (C) shows the distribution of the unmixed components in 3D where the three curved cylindrical inclusions and the layer of fat are distinguished.

Fig. 9

Source spectra automatically obtained as output of the SPAX framework (A); Distribution maps of black, green, blue inclusions, and the layer of fat in the 2D cross-section (B), and in 3D view (C).

In vivo Validation

For the in vivo validations, US-PA images of the mouse hindlimb have been used as input. Fig. 10 shows one of the cross-sectional views of the mouse hindlimb.

Fig. 10

The US cross-sectional image of the mouse hindlimb before segmentation (A) where the yellow dotted region indicates the popliteal fat pad which is within the intramuscular tissue; The PA image of the mouse hindlimb obtained at before segmentation and fluence correction (B); The scale bar size is .

Specifically, Fig. 10 (A) shows the US image before the image segmentation. From the US image, different anatomical structures are visible such as the muscle tissue, part of the inguinal region, and the popliteal lymph node fat pad (highlighted by the dotted yellow line). Fig. 10 (B) shows the PA image at , before the fluence compensation, where the hemoglobin is the prominent absorber.

Fig. 11 (A) represents the segmentation mask where the skin line in yellow, the standard tissue in light blue, and the background in dark blue were identified with the automated approach. Fig. 11 (C) depicts the US image after the segmentation with the removed background. The MCXLAB simulations have been performed on the segmented masks between , by utilizing the tissue optical properties summarized in the Supplementary Materials (see Supplementary Materials; SupMatTable). As an example, the result of the light fluence distribution at is shown in Fig. 11 (B). At each wavelength, the simulated light fluence map has been used to correct the corresponding PA image. The segmentation and the spectral fluence simulation have been performed per each 2D cross-sectional slice and finally extended to match the entire volume. Fig. 11 (D) shows the respective PA image at after the segmentation and fluence compensation. As a comparison, Fig. 11 (E) represents the PA intensity at , along the white dotted line depicted in Fig. 10 (B) and Fig. 11 (D). After the light fluence compensation, the PA intensity increases at depth within the tissue, leading to a higher contrast image. A smaller increase of the signal at the surface is also induced by the compensation, since the light fluence reaching the surface is already slightly diffused. While in deeper regions, the signal has a higher increase due to the spectral coloring compensation.

Fig. 11

Segmentation mask obtained from the US image (A); The light fluence distribution map at , where the dashed line indicates the light source position (B); The US cross-sectional image of the mouse hindlimb after segmentation (C) where the yellow dotted region indicates the popliteal fat pad which is within the intramuscular tissue; The PA image of the mouse hindlimb obtained at after segmentation and fluence correction (D); The scale bar size is ; The PA intensity at before and after segmentation and fluence correction (E) of the pixels along the white dotted line in (D).

After the spectral coloring compensation, the SVD analysis is performed on these corrected sPAI data. This automatically leads to determine the hyperparameter (q), which represents the number of components significantly above the noise floor. Specifically, for the in vivo validation, is the number that has been automatically selected as the hyperparameter. Once the hyperparameter has been identified, the unsupervised spectral unmixing procedure is performed, and the initial spectral components are obtained.

Finally, the advanced superpixel subsampling procedure is implemented to refine the unmixed components. Fig. 12 (A) shows the spectral components obtained as output of the SPAX framework from the hindlimb. In this case, tissue chromophores have been detected. The correlation with theoretical spectra of endogenous tissue chromophores (see Supplementary Materials; SupMatTable) is performed and Fig. 12 (B) shows the spectra that have a positive correlation. Specifically, the obtained spectral signatures are matching with the theoretical spectra of oxy-hemoglobin, deoxy-hemoglobin, and fat with correlation values of , , and respectively. In addition, the spectral has shown a correlation value of with the theoretical absorption spectrum of melanin. Although the correlation value is lower, within the abundance map of this resulted mainly distributed along the skin line (see Fig. 14 (F)). Thus, this component could be related to skin pigmentation, and to have a real comparison with melanin a further investigation would be required. For a qualitative comparison of the spectra shapes, the theoretical absorption spectra of fat and melanin have been also overlaid in Fig. 12 (B) (depicted as dashed lines).

Fig. 12

Source spectra automatically identified by the SPAX framework (A); Source components identified as oxy-deoxy hemoglobin, skin, and fat (B) after the comparison with the theoretical absorption spectra.

Fig. 14

Abundance maps automatically unmixed by the SPAX framework of oxy-hemoglobin, deoxy-hemoglobin before (A-B) and after (C-D) Frangi filter. The distribution maps of fat (E), and skin (F) automatically obtained as output of the SPAX framework; The scale bar size is .

Besides the tissue spectral curves, the SPAX framework also provides the distribution maps of the automatically extracted components. As a comparison, Fig. 13 (A-B) shows the unmixed maps of oxy/deoxy hemoglobin obtained from the linear unmixing approach of the hindlimb. These maps have been obtained by providing the theoretical spectral curves of oxy/deoxy hemoglobin as a priori information.

Fig. 13

Distribution maps of oxy-hemoglobin (A) and deoxy-hemoglobin (B) obtained from linear unmixing; The scale bar size is .

On contrary, Fig. 14 represents the abundance maps obtained as output of the SPAX framework. The maps of oxy-deoxy hemoglobin have been reported before (Fig. 14 (A-B)) and after (Fig. 14 (C-D)) image optimization by means of the Frangi filter (FF). The vessel-like structures have been enhanced after applying the Frangi filter, thus leading to improved visualization of the oxy/deoxy hemoglobin distribution. In addition, to oxy/deoxy hemoglobin, the abundance maps of fat and skin (Fig. 14 (E-F)) have been detected. Specifically, the distribution map of fat obtained as output of the SPAX analysis, as expected, is distributed subcutaneously, within the inguinal region, and in the triangular fat pad.

The maps automatically obtained by the SPAX framework show higher accuracy than the linear unmixing. This could be explained by the iterative approach that leads to converge to a more accurate solution than a fitting-based approach.

Whole-body validation

Furthermore, the data-driven SPAX framework has been extended to whole-body animal imaging. US and multi-wavelength PA volumetric images have been used as input for the SPAX framework. Specifically, two whole-body scans have been performed with the animal in its prone and supine position. Fig. 15 (A-C) depict the coronal US images of the animal in prone and supine position respectively. The US image guides through the anatomical structures, such as kidneys and spine are visible in Fig. 15 (A). As an output of the SPAX framework, multiple spectra have been automatically detected and the correlation with known tissue chromophores has been performed. In addition to the oxy/deoxy hemoglobin, the spectrum of fat has also been detected. Fig. 15 (B-D) show the volumetric distribution of the respective chromophores along the whole-body of the animal.

Fig. 15

Whole-body US image of the mouse in prone and supine position (A-C); Oxy/deoxy hemoglobin distribution depicted in red and green respectively and fat distribution map in blue, in prone and supine view (B-D).

Specifically, the whole-body abundance maps of the oxy/deoxy hemoglobin are shown in red and green color respectively. The yellow color depicts the spatial overlapping of red and green, which are representative of oxy/deoxy hemoglobin. Since the hemoglobin content is mainly distributed within the blood vessels, from prone/supine view the dorsal vascular network and the mammary arteries are primarily observed in the 3-D rendering. The fat content has also been discriminated, and its spatial distribution along the whole-body is shown in the prone and supine view. The fat content is delineated in blue color along the whole-body and mostly distributed in the neck region and in the abdominal area of the mammary glands, as qualitatively observed in the 3-D rendering. Further investigations would be required to validate the components' distributions.

Discussion

The main objective of this study was to automatically recover the unique spectral fingerprint of the tissue chromophores and their volumetric distribution. To this end, we implemented a novel superpixel PA unmixing (SPAX) framework that can automatically detect the distribution of tissue chromophores with improved sensitivity and specificity. The highlight of this framework is the data-driven spectral unmixing approach without any a priori information and user interaction to detect the tissue chromophores. The SPAX framework also includes the modeling of light fluence distribution to compensate for the spectral coloring, and thus prevent the unmixing misinterpretations. The compensation approach utilizes US image segmentation and spectral Monte Carlo simulations, which are based on a predefined library of optical properties.

Spectral coloring is one of the major limitations in sPAI, which causes spectral changes at depth due to light fluence attenuation. The conventional approach to overcome the spectral coloring artifact is either an assumption of homogeneous fluence distribution along the depth or using simulations of pre-defined geometries. This is a complex approach as it requires a lot of a priori information such as the tissue geometry and optical properties of tissue components. Thus, this has the limitation to be generalized to other tissue types. Therefore, within the SPAX framework, we have implemented the automated segmentation of the US image that is simultaneously obtained with the multi-spectral PA images. This enables the precise modeling and simulation of the light fluence distribution at each wavelength of light used during the sPAI acquisition. Overall the SPAX framework facilitates overcoming the spectral coloring artifact and improving the automated extraction of the spectral features.

For the highly sensitive unmixing, an initial SVD analysis and the NNMF are applied to blindly reveal the prominent tissue components from sPAI. Besides, the framework includes a novel spectral superpixel subsampling approach to also identify less prominent tissue chromophores. Hence, less prominent tissue chromophores, such as fat, that generally remain obscured by the most prominent absorbers (like hemoglobin), can be identified. The conventional approach to detect the fat is to use the wavelengths within the Far Infrared (FIR) range, where the absorption from the hemoglobin is negligible. Since the PA excitation sources in the FIR range are having lower energy compared to NIR, the use of this range is constrained to superficial structures and microscopic applications. As the SPAX framework is optimized in the NIR range of , where we have more availability of nanosecond pulsed laser sources, this can be used to detect the distribution of fat along the whole-body of the small animals.

The SPAX framework has been validated on tissue-mimicking phantoms and animal imaging. Specifically, in the agarose phantom, we showed the correction of the spectral coloring distortions at depth. This spectral compensation also promotes the accurate unsupervised unmixing of the molecular components. Besides, the muscle-fat tissue phantom has been used to benchmark the sensitivity of the SPAX framework to detect less and most prominent components. Hence, the proposed SPAX framework corrects the spectral coloring and enables the accurate identification of multiple tissue chromophores. The SPAX framework has also been used for whole-body animal imaging. In this case, the fat component has been detected, despite the oxy- and deoxy-hemoglobin show absorption intensities of two orders greater than fat.

Since the source components are extracted in an unsupervised way from sPAI, we could obtain more spectral characteristics in an unknown order. To facilitate the identification of known tissue chromophores, a correlation matrix could be evaluated between the automatically detected source spectra and a library of theoretical spectra. The correlation values are of significant importance to recognize known components from noise or unknown constituents. In healthy conditions, these values could enable to identify some of the automatically detected spectra as known components. While the identification in disease models requires further investigations, cause the tissue components tend to alter their spectral signatures in pathological conditions, due to changes at a molecular level. For example, the presence of methemoglobin is more pronounced during inflammation, and by using the SPAX framework we may be able to detect these early changes, as the correlation value with the respect to the oxy-hemoglobin might change. Thus, the SPAX framework may have the potential to monitor the spectral changes in disease conditions. However, some of the components could still remain unidentified and require further investigations to be interpreted.

Currently, the framework design includes a fully-automated segmentation step to cluster the skin line, the tissue structures, and the background. For complex tissue structures, the user can guide the US segmentation by selecting some regions of interest for active contours. This refined segmentation will result in a semi-automated approach that could lead to an improved fluence compensation. Although we are using a predefined library of optical tissue properties for the light fluence simulations, it is challenging to retrieve the optical properties from literature or characterize these a priori. Besides, this represents a limitation for disease models where the prediction of the tissue properties is less accurate. Thus, Brochu et al. [41] have shown a possible optimization to retrieve the tissue optical properties iteratively and speed up the computational time using finite element-based simulations, implemented on GPU cards. Very recently Gröhl et al. [56] have proposed an open-source simulation toolkit that could be feasible for these iterative approaches, by enabling realistic high-speed PA simulations. In addition to the tissue properties, to obtain the absolute quantification in PA imaging, improved measurement of the system response, and Grneisen parameter are required. In the future, we aim to address the absolute quantification, by implementing deep learning approaches based on non-explicit light fluence estimation [57].

Finally, further efforts would be required to adjust the tradeoff between acquisition time and spectral resolution. As already demonstrated by Luke et al. [58], the selection of significant wavelengths facilitates the unsupervised unmixing and minimizes the acquisition time. In addition, the current volumetric acquisition setup is based on the translation of the linear array transducer along the animal body. For preclinical studies with high-frequency transducers, this approach guarantees volumetric imaging with a better elevation focus and resolution. However, for clinical translation, the use of 2D matrix array transducers would be ideal as recently shown by Ron et al. [59].

Although the PA technology is already an established preclinical imaging modality, it has still some challenges for clinical translation. Therefore within the SPAX framework, we have included a combination of procedures to overcome the current limitations such as supervised unmixing techniques and misinterpretations caused by spectral coloring. Besides, a superpixel subsampling approach has been implemented to automatically extract multiple tissue components at a volumetric scale without any a priori information. The framework also compensates for the changes of the spectral shape along depth due to light fluence variations. Thus, the SPAX framework would be beneficial and open many possibilities to monitor molecular changes at a volumetric scale in clinical applications.

Conclusions

In summary, we developed the SPAX framework as a strategy to detect molecular tissue components and their volumetric distribution from spectral photoacoustic imaging. The fully-automated unmixing approach has been validated in phantoms and healthy animals. Our initial results have shown that the proposed data-driven algorithm overcomes the spectral coloring limitations and prevents unmixing misinterpretations. To our knowledge, this is a unique algorithm that accounts for the spectral coloring and provides automated detection of tissue spectral signatures in the whole-body of the animal. The initial results show that the SPAX approach is sensitive to spectral unmixing and it has the potential to discern any spectral change, that might occur in pathological conditions. Hence, the SPAX framework could be used as a predictive and monitoring tool to identify early pathological conditions and facilitate the therapeutic follow-up. To this end, in our future work, we also plan to test the SPAX framework in disease models. Finally, the SPAX framework has the potential to expedite the clinical translation of photoacoustic, as it can provide enhanced volumetric tissue characterizations.

Declaration of Competing Interest

The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: Jithin Jose reports a relationship with FUJIFILM VisualSonics that includes: employment. Jithin Jose has patent pending to FUJIFILM VisualSonics.