|
1.IntroductionBiomedical photoacoustic imaging has been a rapidly growing field of research in the last two decades, with significant development in the instrumentation, reconstruction algorithms, and exogenous contrast agents. This has enabled a range of clinical applications of photoacoustic imaging, such as breast imaging,1,2 cardiovascular imaging,3–5 and skin imaging,6,7 as well as preclinical studies using endogenous contrast, such as hemoglobin, for imaging the vasculature in the mouse brain8–10 and in tumor models11–13 or exogenous contrast agents14,15 and reporter genes16,17 for molecular or cellular imaging. To acquire a photoacoustic image, the tissue is illuminated with nanosecond laser pulses, and the optical energy is absorbed by the chromophores in the tissue. This leads to a small temperature and pressure rise, generating broadband acoustic waves that propagate to the surface of the tissue where they are detected and used to reconstruct images of the pressure rise. By acquiring multiple images using excitation light with different wavelengths, photoacoustic imaging can be extended to a spectroscopic technique that can in principle quantify the concentration of the chromophores and hence reveal functional information about the subject. However, the image contrast depends on the absorbed optical energy, which is nonlinearly related to the chromophore concentrations, because it is a product of the local light fluence and the absorption coefficient. The light fluence is not only spatially and spectrally nonuniform but also depends on the unknown chromophore concentrations. Therefore, it presents the key challenge for quantitative photoacoustic imaging. The most straightforward method of accounting for the nonlinear fluence distortion is to divide the multiwavelength photoacoustic images by an approximation of the fluence at each wavelength, which can be obtained by for example assuming homogeneous optical properties.18,19 In this way, the effect of the fluence is approximately “cancelled out,” such that the images represent a linear sum of chromophore concentrations that can be unmixed using linear spectroscopic inversion (SI), as in conventional optical spectroscopy, which uses the pseudoinverse of the matrix of the known specific absorption spectra. Linear methods with approximate fluence correction are easy to implement and fast to compute; hence, they are attractive for in vivo studies. However, the accuracy of SI is limited for nonsuperficial imaging depths,20 because it depends on the accuracy of the fluence correction, which is only an approximation of the true fluence distribution in practice. Improved quantification accuracy in some scenarios has been demonstrated using more sophisticated inversion schemes, including methods involving light source modulation,21 sparse signal decomposition,22 fluence modeling using light transport equations,23–25 or generating spectral libraries.26 However, these methods are more challenging to implement for in vivo images due to issues such as model mismatch, increased level of complexity, and/or high computational demand. One way to enable quantification in a wider range of applications would be to retain the simplicity of a linear model but increase the robustness to errors in the fluence correction. This may, in some cases, be achieved using independent component analysis (ICA),27–29 which is a linear source separation method that decomposes the multiwavelength image data into components representing the individual chromophores by assuming that they are mutually statistically independent. ICA has the advantages of being fast, easy to implement, and computationally inexpensive. It was first proposed for unmixing photoacoustic images by Glatz et al.,28 who showed that it can provide spectral unmixing with less cross talk error than SI. It has also been shown in some cases that applying ICA to the logarithm of the photoacoustic images can result in more accurate estimations of the relative concentrations compared to using ICA on the images without taking the logarithm.30 However, since the components of interest are represented as a ratio to the background chromophore concentration, this approach requires that the chromophore concentration in the background tissue is known. ICA has been used in an in vivo multiwavelength photoacoustic imaging study of near-infrared fluorescent protein-expressing glioma cells,31 as well as in studies using exogenous probes that bind to injured cells32 or cells undergoing apoptosis.33 In the above-mentioned studies, the main purpose of applying ICA was to aid the visualization of the probe by distinguishing it from the background tissue, rather achieving quantitatively accurate representation of the relative chromophore concentrations. As a result, despite having been applied in several in vivo studies, the quantitative accuracy of ICA has been not rigorously assessed for photoacoustic imaging. Therefore, the aim of this paper is to investigate the conditions required for ICA to provide quantitatively accurate unmixing, which required an approximate fluence correction step. Simulated multiwavelength photoacoustic images are used to assess the robustness of ICA against inaccuracies in the fluence correction and demonstrate the effect of retaining different numbers of dimensions as a preprocessing step to ICA. The results of ICA are also compared to SI for experimentally acquired images of a tissue-mimicking phantom. Furthermore, the performances of ICA and SI are analyzed for inversions using the experimental images with varying spectral ranges. We demonstrate both the advantages of ICA and its limitations, by including examples of cases where it provides accurate quantification, as well as cases where it breaks down. 2.Principles of Independent Component AnalysisAssuming that the acoustic reconstruction of the initial pressure rise due to the absorption of optical energy is perfect, so there are no image artifacts, a set of photoacoustic images with M voxels acquired at N wavelengths, whose contrast originates from K chromophores, can be described using matrix notation as where ∘ denotes the Hadamard product, P is a (N×M) matrix with each row corresponding to the reconstruction of the initial pressures at one of the wavelengths, ϕ is the (N×M) fluence matrix with each row representing the spatially varying fluence at one of the wavelengths, A is the (N×K) mixing matrix with each column representing the wavelength-dependent specific absorption coefficient of one of the chromophores, C is the (K×M) concentrations matrix with each row representing the concentration distribution of one of the chromophores, and Γ is a (K×M) matrix, where all rows are identical and equal to the spatially varying Grüneisen parameter.Both ICA and SI are unmixing methods based on a linear model. The initial photoacoustic pressure is, however, nonlinearly related to the concentrations, due to the spectrally and spatially varying fluence matrix in Eq. (1). An approximation to linearity can be achieved by dividing P element-wise by an approximate estimate of the fluence, ˘ϕ, such that In SI, the mixing matrix with the specific absorption coefficients is known. Hence, the chromophore concentrations can be estimated using the pseudoinverse of the mixing matrix, A† SI cannot separate Γ from the concentrations in general, as Γ does not vary with wavelength and scales with the concentrations. In in vivo imaging studies, Γ is often assumed to be a spatially constant scaling factor, even though it may vary significantly between different tissue types.34 Unlike in SI, the mixing matrix is not known in ICA. Instead, ICA aims to decompose the multiwavelength photoacoustic images into the source components, which represent the individual chromophores, based on the assumption that they are statistically mutually independent of each other. Statistical independence is defined as follows: two events, X and Y, are independent of each other if the probability of X occurring does not in any way influence the probability of Y occurring. In quantitative photoacoustic imaging, this requires that the knowledge of the concentration of a chromophore at a location is not affected by the information provided about another chromophore at the same location. In other words, if the distributions of certain probes are unrelated to other chromophores, such as blood and lipids, then the concentration of the probe does not reveal any information about the amount of blood in the same pixel; hence, they are independent from each other. Example of such independent probes can be found in the imaging of reporter genes, such as genetically encoded fluorescent proteins18 or tyrosinase-expressing cells,35 and exogenous contrast agents that accumulate in regions unrelated to the blood distribution, such as in lymphatic vessels and lymph nodes.36 Examples of chromophores that are unlikely to be independent include oxy- and deoxyhemoglobin, as their distributions are typically highly correlated. Therefore, ICA is not suitable for unmixing oxy- and deoxyhemoglobin for the estimation of blood oxygenation. In ICA, the independent components are found from the multiwavelength images by searching for a mixing matrix, W, whose inverse can be multiplied by the matrix with the mixed signals to give the matrix of the source components, S where the rows of S have maximal statistical independence. The unmixed component in each row of S corresponds to one of the chromophores, such that S≈Γ∘C, provided that the true spatial distributions of the chromophore concentrations are independent.In the widely used ICA algorithm FastICA,27 the search for the most independent components is based on the central limit theorem, which states that the probability distribution of the sum of multiple independent variables will always be more Gaussian than the distribution of one of the independent variables alone. The Gaussianity of the probability distribution of the output components is measured using negentropy, which is a normalized version of entropy. Hence, the independent source components are found by iteratively searching for a mixing matrix that maximizes the negentropy of each row of S. Unlike SI, ICA does not recover the concentrations for different chromophores to a common scale, because the magnitude of each independent source is unknown. It is clear from the ICA model in Eq. (4) that since neither W nor S is fixed, the magnitude of each source component can be arbitrarily changed by multiplying any row in S and dividing the corresponding column in W† with the same constant. This means that while the relative concentration to other spatial locations for each chromophore can be recovered, the concentration of one chromophore cannot be compared to that of another chromophore. 3.Generating Multiwavelength Photoacoustic ImagesIf the fluence was estimated exactly, such that ˘ϕ=ϕ, the approximate model in Eq. (2) would be an exact equality, and if the chromophore distributions were independent, both ICA and SI would achieve accurate unmixing of the chromophores, such that the outputs of both inversions are exactly equal to Γ∘C. In practice, however, ˘ϕ is only an approximation of ϕ, and both methods will lead to errors. Both experimental and numerically simulated multiwavelength photoacoustic images of tissue-mimicking phantoms were generated to investigate the accuracy of ICA for various levels of absorption and heterogeneity in the tissue structure, as well as the effect of dimension reduction and wavelength selection, in order to understand the extent to which ICA can be used as a reliable quantification method. The same data sets are also unmixed with SI for comparison. The experimental and numerical phantoms were designed such that the following criteria were satisfied:
Hence, using these phantoms, it is possible to see how the inaccuracies in an approximate fluence adjustment affect the performance of the unmixing methods. 3.1.Experimental Photoacoustic Image AcquisitionEight capillary tubes (Paradigm Optics, Morcap 83) with inner diameter 590 μm and wall thickness 66.5 μm were used to construct the tissue-mimicking phantom. The tubes were arranged horizontally in two vertical lines as shown in Fig. 1, which shows a cross section of the phantom in the experimental setup. The left column of tubes was filled with four different concentrations of copper(II) chloride dihydrate (CuCl2·2H2O) dissolved in deionized water. The concentrations were in the ratios 1:2:3:4, where the uppermost tube had the lowest concentration of 5.2 gL−1 while the bottom tube had the highest concentration of 20.8 gL−1. The right column of tubes was filled with solutions of nickel(II) chloride hexahydrate (NiCl2·6H2O), also with concentration ratio of 1:2:3:4. These contrast agents simulate different absorbers in the tissue and will be referred to as CuCl2 and NiCl2. The absorption spectra of all chromophores are shown in Fig. 2. Since the specific absorption coefficient of NiCl2 is approximately one order of magnitude lower than that of CuCl2 [see Fig. 2(a)], the concentrations of NiCl2 were set to be ∼10 times higher than CuCl2 to give similar optical absorption. Therefore, the uppermost and bottom tubes in the right column had NiCl2 concentrations of 55.1 and 220.3 gL−1, respectively. The average absorption of all tubes was 0.25 mm−1 at 810 nm. Fig. 1DownloadFig. 2DownloadThe tubes were submerged at depths between 1.0 and 6.1 mm in a background solution containing India ink (951 black, Winsor & Newton) and 1% Intralipid diluted in deionized water, such that they provide an absorption of ∼0.013 mm−1 at 810 nm and a scattering coefficient of ∼1 mm−1, which are comparable to realistic values in soft tissue.38 The Grüneisen parameter of CuCl2 and NiCl2 is both known to scale with their concentrations such that where the βi coefficients are 5.80×10−3 Lg−1 and 2.25×10−3 Lg−1 for CuCl2 or NiCl2, respectively,39 ci denotes the concentrations of CuCl2 or NiCl2, and ΓH2O is the Grüneisen parameter of water, which is 0.12 Lg−1 at 22°C.40 The phantom was illuminated at the surface using pulsed laser radiation (SpitLight 600 OPO, InnoLas Laser GmbH) of 6-ns duration. The beam diameter was ∼10 mm and the pulse energy was 7 to 13 mJ at the surface of the phantom depending on wavelength. A small fraction of the excitation light was reflected into an integrating sphere with a photodiode to measure the laser pulse energy, which was used to normalize the photoacoustic signals. A planar Fabry–Perot polymer film interferometer (FPI)41 situated at the bottom of the tray containing the phantom was used for ultrasound detection. The acoustic pressure waves that propagate to the sensor modulate the thickness and hence the reflectivity of the FPI. Thus, by scanning the surface of the FPI with an interrogation laser and recording the reflected intensity, a time-varying spatial mapping of the photoacoustic signals was generated. The sensor was interrogated along a 20-mm line with 10-μm step size to reconstruct two-dimensional (2-D) images of the cross section of the tubes at 18 equally spaced wavelengths between 750 and 1090 nm.3.2.Numerically Simulated Photoacoustic ImagesA total of 37 sets of simulated 2-D multiwavelength photoacoustic images based on the same structure and chromophores as the experimental phantom were generated to evaluate the performance of ICA for different levels of absorption. The dimensions of the simulated images were 7.4×20.0 mm2 and spacing between the elements was 20 μm. The main differences to the experimental phantom are the absence of the tube walls and the fact that the India ink and Intralipid are also present within the tubes. The concentration ratios of 1:2:3:4 between the tubes in each column, as well as the 10-fold ratio between CuCl2 and NiCl2 concentrations were kept constant for all data sets, while the absolute concentrations of CuCl2 and NiCl2 were varied in different ways in three case studies:
Fig. 3DownloadThe absorption coefficient was calculated for the same 18 wavelengths as the experimental phantom for each data set based on the chromophore concentrations and their known specific absorption coefficients. The reduced scattering coefficient was kept the same as the experimental phantom for all data sets. Based on the distribution of the absorption and reduced scattering coefficients, the fluence was modeled using the diffusion approximation to the radiative transfer equation with the MATLAB software package TOAST++42 for a 10-mm wide Gaussian light source at the top of the phantom. The initial pressure was found by multiplying the fluence by the optical absorption coefficient and the Grüneisen parameter. The photoacoustic wave propagation from the initial pressure distribution was modeled using the MATLAB toolbox k-Wave43 based on a k-space pseudospectral method. The time-varying photoacoustic pressure was recorded at the bottom of the numerical phantom. 3.3.Image ProcessingBoth the experimental and simulated images were reconstructed using the time-reversal method,44 which backpropagates the recorded time series of the photoacoustic pressure into the image domain to form an image of the initial pressure. Gaussian noise with standard deviation equal to 3% of the peak intensity of data set in Case I with the same concentrations as the experimental phantom was added to all simulated images. An approximate fluence adjustment was performed by dividing the reconstructions of the initial pressure by the estimation of the fluence, ˘ϕ. The ˘ϕ was calculated based on the simplifying assumptions that the medium is a semi-infinite optically homogeneous slab illuminated by infinitely wide plane waves, such that the one-dimensional (1-D) solution to the diffuse approximation can be applied where z is the depth from the illuminated surface, ϕ0 is the fluence at the illuminated surface, and μeff is the effective absorption coefficient given by μeff=√3μa(μa+μ′s), where the μa and μ′s are assumed to be known and equal to that of the background solution. This exponential fluence adjustment is chosen in this study because it can be applied straightforwardly in practice for in vivo photoacoustic images, with μeff estimated as an average parameter for the tissue. A similar correction can also be applied for imaging systems with circular illumination, in which case the fluence can be approximated as a radially symmetric function.45 While the fluence correction step is necessary for obtaining accurate estimates of the relative chromophore concentrations, it is still an approximate method and does not fully remove the fluence distortion.The fluence adjustment overamplifies the noise in the regions away from the light source where the values of ˘ϕ are small. To avoid this, a 6.5×6.5-mm2 region of interest in the fluence adjusted experimental images was decomposed using ICA and SI. The raw experimental images, as well as the fluence estimations and the fluence adjusted images, are shown in Fig. 4. As expected based on the specific absorption coefficients of the chromophores, the raw images at 750 and 1090 nm show higher intensity for the right column of tubes containing NiCl2, while the image acquired at 810 nm shows higher intensity for the left column of tubes containing CuCl2. The streak patterns extending from the tubes and the negative values in the reconstruction are artifacts due to the limited detection aperture of the planar sensor. The bottom row of figures shows that after the fluence adjustment, the image intensity at the tubes increases with depth, which is expected as the chromophore concentrations are higher inside the deeper tubes. Since the negative values in the images lack physical meaning, they were set to zero before the unmixing. Fig. 4DownloadThe FastICA algorithm requires further preprocessing of the multiwavelength fluence adjusted images before unmixing. The images are mean subtracted and whitened using principal component analysis, and the dimensions of the image data are reduced such that only the three principal components (PCs) with the largest eigenvalues of the covariance matrix of ˘P are unmixed with FastICA. Since the ordering of the estimated independent components is arbitrary, the output components are identified to the corresponding chromophores by calculating the sum of squared differences between the normalized columns in the estimated mixing matrix and the normalized known absorption spectra of the chromophores of interest. For comparison, the fluence adjusted images were also unmixed using SI, based on the known specific absorption spectra of CuCl2, NiCl2, and the background solution. 4.Unmixing Using Independent Component Analysis and Spectroscopic Inversion4.1.Accuracy as a Function of AbsorptionThe fluence adjusted and preprocessed simulated photoacoustic images at 18 wavelengths were decomposed using ICA and SI. The unmixed components were normalized to their average value at the tube with the highest concentration of the relevant chromophore. The normalized components were compared to the expected components, which are the true concentrations of the relevant chromophores multiplied by the Grüneisen parameter. Three types of errors are defined to provide a quantitative assessment of the accuracy of the unmixing methods:
For example, in the estimated CuCl2 component, δc is the average error at the left column of tubes, δf is the average error at the right column of tubes, and δb is the average error outside the tubes. For Case III, the errors in the additional region are included in the definition of δc and δf. The three types of errors are plotted in Fig. 5 as functions of the average absorption of the tubes at 810 nm, μtubesa, for Case I, the absorption of the background solution at 810 nm, μbkga, for Case II, and the absorption of the additional region at 810 nm, μadda, for Case III. In Case I, the δc error for SI increases with μtubesa at a nearly constant rate for the whole range of absorptions between 0.05 and 0.75 mm−1, as shown with asterisks. This is expected as the estimated fluence is based on an optically homogeneous region, and therefore, the accuracy of ˘ϕ is reduced for increasing μtubesa. The δf and δb of SI errors remain below 7% for this entire range of absorption. The errors of ICA are represented with dots and all three types of errors are >10% for μtubesa<0.15 mm−1 in Case I. The large errors at low absorption may be due to the lower signal-to-noise ratio (SNR), which leads to the first three PCs containing a smaller fraction of the total variance of the data (see Sec. 4.2). When μtubesa is increased beyond 0.15 mm−1, the errors of ICA are comparable to SI for the data points with relatively low μtubesa, and as μtubesa is further increased, ICA results in smaller δc than SI. However, when μtubesa>0.55 mm−1, the performance of ICA abruptly deteriorates as the errors increase to large values beyond the scale of the plots. This is because beyond this threshold, the accuracy of the fluence estimation is sufficiently low to lead to nonlinearities so severe that the independent components found by ICA cannot be correctly identified to the relevant chromophores based on the estimated absorption spectra, hence resulting in a significant increase in errors. This absorption level is comparable to that of blood, which is 0.46 mm−1 at 810 nm38 (assuming a total hemoglobin concentration of 150 gL−1 and 100% oxygenation). The blue circles show the errors that would have been obtained if the correct components were manually selected. This manual selection may not be possible in practical applications because prior knowledge of the expected chromophore distribution may not be available. Fig. 5DownloadIn Case II, accurate unmixing with all three types of errors ≲10% was obtained for μbkga<0.06 mm−1 at 810 nm using ICA and <0.10 mm−1 using SI. When μbkga is increased further, the errors for both unmixing methods increase rapidly. However, these thresholds are significantly higher than the absorption coefficient of the common types of biological tissue, which is typically ≈0.01 mm−1 at 810 nm.38 The performance of both ICA and SI is dependent on the accuracy of the fluence adjustment, which in this study is mainly determined by the level of spatial inhomogeneity in the optical properties of the phantom. This is highlighted by the fact that both methods fail to produce accurate results in the presence of the additional region with increasing concentration of CuCl2 in Case III, where the errors increase with μadda. 4.2.Component SelectionThe thresholds above which the independent components cannot be identified to the correct chromophores can potentially be shifted toward higher absorption levels if the multiwavelength image data are decomposed into fewer independent components. This can be realized by keeping fewer PCs in the preprocessing stage to reduce the dimensions of the data. Figure 6 shows the δc, δf, and δb errors from unmixing the images in Case I with ICA using two, three, or four PCs. The results show that if only two PCs were processed with ICA, the unmixed independent components can be correctly identified as CuCl2 or NiCl2 for the full range of absorption levels investigated in Case I, and no abrupt increase in error is observed. However, using only two PCs also results in larger unmixing errors for the higher absorption levels. When four PCs are used, the error trends are similar to using three PCs, but a slight lower shift of the range of the absorption in which ICA results in accurate unmixing is observed. This shift may be explained by the fact that, when the main features in the image have lower contrast, a smaller fraction of the total variance of the data is included in the first three PCs. Fig. 6Download4.3.Experimental Images and the Condition Number of the Mixing MatrixThe experimental images were acquired at 18 wavelengths over a spectral range where the absorption spectra of the chromophores have distinct features. However, in in vivo imaging applications, some chromophores in the tissue may have flatter and/or less unique absorption spectra. To investigate how the unmixing methods deal with poorly conditioned absorption spectra, ICA and SI were applied on data sets with reduced spectral ranges. The spectral range was reduced in 15 steps by removing the images of the longest wavelengths, such that in the final inversion (the inversion with the fewest wavelengths) only the images at the wavelengths 750, 770, and 790 nm were used. The condition number of the mixing matrix is plotted as a function of the longest wavelength, λmax, used in the inversion in Fig. 7. The figure shows that the condition number increases significantly when the range of wavelengths is limited to λmax<830 nm. Fig. 7DownloadThe unmixed components were normalized and the errors δc, δf, and δb are defined in the same way as for the simulated images. The three types of errors are plotted as a function of λmax used in ICA and SI in Fig. 8. When the full range of wavelengths between 750 and 1090 nm is used, both methods result in accurate unmixing with the δc, δf, and δb errors equal to 11%, 4%, and 3% for ICA and 10%, 5%, and 3% for SI, respectively. This result is in agreement with the simulated data which showed high accuracy for lower absorption levels for both ICA and SI. The high accuracy is maintained for the inversions using λmax>930 nm. When λmax is further reduced, the accuracy of SI rapidly deteriorates, while the errors associated with ICA remain relatively low. An example of the unmixed components is shown in Fig. 9, where λmax=850 nm. Both the CuCl2 and NiCl2 components are estimated accurately with ICA. On the other hand, SI is able to unmix the NiCl2 component with relatively low errors, but its estimate of the CuCl2 component is noisy and contains large errors. Fig. 8DownloadFig. 9DownloadTwo key factors are likely to affect the unmixing results when the spectral range is reduced: first, due to the spectral variation in the absorption of the contrast agents in the tubes and the background solution, the fluence correction is more accurate for some wavelength combinations than others. This results in variations in the accuracy of SI when certain wavelengths are excluded. Second, the condition number of the known mixing matrix increases by several orders of magnitude when λmax<830 nm, as shown in Fig. 7. Therefore, the fact that the inversion is more ill-conditioned is likely to be the dominant cause of the large errors of SI for the inversions where λmax<830 nm. However, both these factors affect ICA less significantly. As shown with the simulated image data in Sec. 4.1, ICA is more robust against nonlinearities caused by poor fluence correction, and therefore, the errors of ICA remain low when wavelength range changes. Furthermore, the low errors of ICA when λmax<830 nm suggest that unmixing based on statistical independence can potentially be used to obtain accurate separation of the chromophores when SI performs poorly due to ill-conditioning of the mixing matrix. This is possible because ICA does not rely on the known spectra for unmixing and can, therefore, tolerate ill-conditioning better than SI. 5.Summary and DiscussionAccurate unmixing of chromophores from multiwavelength photoacoustic images in the absence of accurate knowledge of the fluence is a challenging task. This study has shown that a simple exponential fluence correction, which is straightforwardly applicable in practice, allows linear models—ICA and SI—to be used to unmix the chromophores to useful degree of accuracy under certain conditions. Experimental and numerical phantoms were designed for which the assumptions underlying both ICA (statistical independence) and SI (full rank molar absorption coefficient matrix) were true. Either approach would, therefore, give accurate results for a perfect fluence correction. When using the approximate, but practical, fluence correction [Eq. (6)], it was shown that ICA results in smaller unmixing errors compared to SI in two circumstances: first, when the absorption level of the contrast agents is ∼0.4 to 0.5 mm−1, ICA is more robust against nonlinearities caused by inaccurate fluence estimation than SI, because it allows the mixing matrix to vary in order to produce the most independent components. Second, when the mixing matrix with the known absorption spectra is more ill-conditioned, ICA provides significantly more accurate results as it avoids using the fixed mixing matrix by searching for the chromophore components based on their statistical independence instead. When the average absorption of the features of interest is higher than a certain threshold (which in this particular phantom is at 0.55 mm−1), leading to greater differences between the estimated and the true fluence, the errors of ICA were shown to abruptly increase as the output components corresponding to the chromophores of interest could not be identified by comparing the estimated spectra with the known. This suggests that, with sufficient SNR, ICA can provide accurate unmixing for absorptions up to a threshold. The threshold will likely be dependent on the absorption spectra, the spatial structure of the chromophores, and the accuracy of the fluence adjustment. Given the potential presence of these upper absorption thresholds under which ICA provides relatively low errors, it would be ideal to use lower concentrations of contrast agents. This would require the imaging system to have high sensitivity to obtain sufficient SNR. ICA is a suitable unmixing method that is robust to the small errors in the fluence correction that are unavoidable in practical scenarios, provided that the chromophores are known to be mutually statistically independent. This assumption is valid for the chromophores in the phantoms used in this study and in applications such as unmixing some exogenous contrast agents from the tissue in the background. However, the independence criterion is not always fulfilled for all tissue chromophores, hence limiting the range of applications of ICA. The approximate fluence correction based on spatially homogenous optical properties cannot account for highly absorbing structures in the vicinity of the regions of interest. This is demonstrated in Case III, where the additional absorbing region is shown to cause errors for both ICA and SI. Despite this being an extreme situation where the chromophores of interest are completely surrounded by the additional absorber, it illustrates that the general applicability of linear methods relying on a simple exponential fluence correction is limited in biological tissues with large highly absorbing regions. The number of components that will be estimated is fixed in SI and determined by the dimensions of the mixing matrix, which is equal to the number of chromophores. This study showed that in ICA, it is not straightforward to choose how many independent components the multiwavelength images should be decomposed into. Since larger dimensionality leads to difficulty in identifying the components as chromophores, one should ideally retain the minimum number of PCs that contain a sufficient fraction of the total variance to explain the data well. In this study, PCs representing >75% of the variance needed to be kept for further processing with ICA to provide accurate results. However, there are no general guidelines based on theoretical principles for the optimal choice of dimension reduction. As discussed in Sec. 2, the magnitude of the independent components estimated using ICA is arbitrary because both W and S are unknown. However, some prior knowledge often exists for the absorption spectra of the chromophores. These known specific absorption spectra can potentially be used to fix the magnitude of the independent components, such that the relative concentration between different chromophores is scaled correctly, and hence reducing ambiguities of ICA. One simple approximate scaling method would be to divide the estimated independent components by a scaling factor equal to the ratio between the mean of each estimated spectrum and the corresponding known spectrum. 6.ConclusionIn conclusion, ICA offers a fast and simple alternative to unmixing multiwavelength photoacoustic images into components representing individual chromophores, provided that the spatial distributions of the chromophore concentrations are statistically independent. When a first-order fluence adjustment has been applied and the absorption is within certain ranges and relatively spatially homogeneous, ICA can provide accurate quantification of the relative chromophore concentrations. The results of ICA depend on the choice of dimensions retained in the preprocessing step, as accurate results require that the components can be identified as the correct chromophores. It was shown that ICA outperforms SI when mixing matrix is ill-conditioned, and that ICA is more robust to errors in the fluence correction compared to SI. DisclosuresThe authors have no relevant financial interests in this article and no potential conflicts of interest to disclose. AcknowledgmentsThis work was funded by the University College London EPSRC Centre for Doctoral Training in Medical Imaging. ReferencesR. A. Kruger et al.,
“Photoacoustic angiography of the breast,”
Med. Phys., 37
(11), 6096
–6100
(2010). https://doi.org/10.1118/1.3497677 MPHYA6 0094-2405 Google Scholar
S. Manohar et al.,
“The Twente photoacoustic mammoscope: system overview and performance,”
Phys. Med. Biol., 50
(11), 2543
–2557
(2005). https://doi.org/10.1088/0031-9155/50/11/007 PHMBA7 0031-9155 Google Scholar
T. J. Allen et al.,
“Spectroscopic photoacoustic imaging of lipid-rich plaques in the human aorta in the 740 to 1400 nm wavelength range,”
J. Biomed. Opt., 17
(6), 061209
(2012). https://doi.org/10.1117/1.JBO.17.6.061209 JBOPFO 1083-3668 Google Scholar
B. Wang et al.,
“In vivo intravascular ultrasound-guided photoacoustic imaging of lipid in plaques using an animal model of atherosclerosis,”
Ultrasound Med. Biol., 38
(12), 2098
–2103
(2012). https://doi.org/10.1016/j.ultrasmedbio.2012.08.006 USMBA3 0301-5629 Google Scholar
K. Jansen et al.,
“Intravascular photoacoustic imaging of human coronary atherosclerosis,”
Opt. Lett., 36
(5), 597
–599
(2011). https://doi.org/10.1364/OL.36.000597 OPLEDP 0146-9592 Google Scholar
J.-T. Oh et al.,
“Three-dimensional imaging of skin melanoma in vivo by dual-wavelength photoacoustic microscopy,”
J. Biomed. Opt., 11
(3), 034032
(2006). https://doi.org/10.1117/1.2210907 JBOPFO 1083-3668 Google Scholar
C. P. Favazza, L. A. Cornelius and L. V. Wang,
“In vivo functional photoacoustic microscopy of cutaneous microvasculature in human skin,”
J. Biomed. Opt., 16
(2), 026004
(2011). https://doi.org/10.1117/1.3536522 JBOPFO 1083-3668 Google Scholar
X. Wang et al.,
“Noninvasive laser-induced photoacoustic tomography for structural and functional in vivo imaging of the brain,”
Nat. Biotechnol., 21
(7), 803
–806
(2003). https://doi.org/10.1038/nbt839 NABIF9 1087-0156 Google Scholar
J. Laufer et al.,
“Three-dimensional noninvasive imaging of the vasculature in the mouse brain using a high resolution photoacoustic scanner,”
Appl. Opt., 48 D299
–D306
(2009). https://doi.org/10.1364/AO.48.00D299 APOPAI 0003-6935 Google Scholar
E. W. Stein, K. Maslov and L. V. Wang,
“Noninvasive, in vivo imaging of blood-oxygenation dynamics within the mouse brain using photoacoustic microscopy,”
J. Biomed. Opt., 14
(2), 020502
(2009). https://doi.org/10.1117/1.3095799 JBOPFO 1083-3668 Google Scholar
G. Ku et al.,
“Imaging of tumor angiogenesis in rat brains in vivo by photoacoustic tomography,”
Appl. Opt., 44
(5), 770
–775
(2005). https://doi.org/10.1364/AO.44.000770 APOPAI 0003-6935 Google Scholar
Y. Lao et al.,
“Noninvasive photoacoustic imaging of the developing vasculature during early tumor growth,”
Phys. Med. Biol., 53
(15), 4203
–4212
(2008). https://doi.org/10.1088/0031-9155/53/15/013 PHMBA7 0031-9155 Google Scholar
J. Laufer et al.,
“In vivo preclinical photoacoustic imaging of tumor vasculature development and therapy,”
J. Biomed. Opt., 17
(5), 056016
(2012). https://doi.org/10.1117/1.JBO.17.5.056016 JBOPFO 1083-3668 Google Scholar
A. De La Zerda et al.,
“Carbon nanotubes as photoacoustic molecular imaging agents in living mice,”
Nat. Nanotechnol., 3
(9), 557
–562
(2008). https://doi.org/10.1038/nnano.2008.231 NNAABX 1748-3387 Google Scholar
P.-C. Li et al.,
“In vivo photoacoustic molecular imaging with simultaneous multiple selective targeting using antibody-conjugated gold nanorods,”
Opt. Express, 16 18605
–18615
(2008). https://doi.org/10.1364/OE.16.018605 OPEXFF 1094-4087 Google Scholar
L. Li et al.,
“Photoacoustic imaging of lacZ gene expression in vivo,”
J. Biomed. Opt., 12
(2), 020504
(2007). https://doi.org/10.1117/1.2717531 JBOPFO 1083-3668 Google Scholar
C. Qin et al.,
“Tyrosinase as a multifunctional reporter gene for photoacoustic/MRI/PET triple modality molecular imaging,”
Sci. Rep., 3 1490
(2013). https://doi.org/10.1038/srep01490 SRCEC3 2045-2322 Google Scholar
D. Razansky et al.,
“Multispectral opto-acoustic tomography of deep-seated fluorescent proteins in vivo,”
Nat. Photonics, 3
(7), 412
–417
(2009). https://doi.org/10.1038/nphoton.2009.98 NPAHBY 1749-4885 Google Scholar
S. Kim et al.,
“In vivo three-dimensional spectroscopic photoacoustic imaging for monitoring nanoparticle delivery,”
Biomed. Opt. Express, 2
(9), 2540
–2550
(2011). https://doi.org/10.1364/BOE.2.002540 BOEICL 2156-7085 Google Scholar
R. Hochuli, P. C. Beard and B. Cox,
“Accuracy of approximate inversion schemes in quantitative photoacoustic imaging,”
Proc. SPIE, 8943 89435V
(2014). https://doi.org/10.1117/12.2039825 PSISDG 0277-786X Google Scholar
S. Choi et al.,
“Wavelength-modulated differential photoacoustic spectroscopy (WM-DPAS) for noninvasive early cancer detection and tissue hypoxia monitoring,”
J. Biophotonics, 9
(4), 388
–395
(2016). https://doi.org/10.1002/jbio.v9.4 Google Scholar
A. Rosenthal, D. Razansky and V. Ntziachristos,
“Quantitative optoacoustic signal extraction using sparse signal representation,”
IEEE Trans. Med. Imaging, 28 1997
–2006
(2009). https://doi.org/10.1109/TMI.2009.2027116 ITMID4 0278-0062 Google Scholar
B. T. Cox, S. R. Arridge and P. C. Beard,
“Estimating chromophore distributions from multiwavelength photoacoustic images,”
J. Opt. Soc. Am. A, 26 443
–455
(2009). https://doi.org/10.1364/JOSAA.26.000443 JOAOD6 0740-3232 Google Scholar
G. Bal and K. Ren,
“On multi-spectral quantitative photoacoustic tomography in diffusive regime,”
Inverse Prob., 28
(2), 025010
(2012). https://doi.org/10.1088/0266-5611/28/2/025010 INPEEY 0266-5611 Google Scholar
T. Saratoon et al.,
“A gradient-based method for quantitative photoacoustic tomography using the radiative transfer equation,”
Inverse Prob., 29
(7), 075006
(2013). https://doi.org/10.1088/0266-5611/29/7/075006 INPEEY 0266-5611 Google Scholar
S. Tzoumas et al.,
“Eigenspectra optoacoustic tomography achieves quantitative blood oxygenation imaging deep in tissues,”
Nat. Commun., 7 12121
(2016). https://doi.org/10.1038/ncomms12121 NCAOBW 2041-1723 Google Scholar
A. Hyvärinen and E. Oja,
“Independent component analysis: algorithms and applications,”
Neural Network, 13 411
–430
(2000). https://doi.org/10.1016/S0893-6080(00)00026-5 Google Scholar
J. Glatz et al.,
“Blind source unmixing in multi-spectral optoacoustic tomography,”
Opt. Express, 19
(4), 3175
–3184
(2011). https://doi.org/10.1364/OE.19.003175 OPEXFF 1094-4087 Google Scholar
L. An and B. Cox,
“Independent component analysis for unmixing multi-wavelength photoacoustic images,”
Proc. SPIE, 9708 970851
(2016). https://doi.org/10.1117/12.2208137 PSISDG 0277-786X Google Scholar
X. L. Deán-Ben et al.,
“Estimation of optoacoustic contrast agent concentration with self-calibration blind logarithmic unmixing,”
Phys. Med. Biol., 59 4785
–4797
(2014). https://doi.org/10.1088/0031-9155/59/17/4785 PHMBA7 0031-9155 Google Scholar
N. C. Deliolanis et al.,
“Deep-tissue reporter-gene imaging with fluorescence and optoacoustic tomography: a performance overview,”
Mol. Imaging Biol., 16
(5), 652
–660
(2014). https://doi.org/10.1007/s11307-014-0728-1 Google Scholar
A. Taruttis et al.,
“Multispectral optoacoustic tomography of myocardial infarction,”
Photoacoustics, 1
(1), 3
–8
(2013). https://doi.org/10.1016/j.pacs.2012.11.001 Google Scholar
A. Buehler et al.,
“High resolution tumor targeting in living mice by means of multispectral optoacoustic tomography,”
EJNMMI Res., 2
(1), 14
(2012). https://doi.org/10.1186/2191-219X-2-14 Google Scholar
D.-K. Yao et al.,
“Photoacoustic measurement of the Gruneisen parameter of tissue,”
J. Biomed. Opt., 19
(1), 017007
(2014). https://doi.org/10.1117/1.JBO.19.1.017007 JBOPFO 1083-3668 Google Scholar
A. P. Jathoul et al.,
“Deep in vivo photoacoustic imaging of mammalian tissues using a tyrosinase-based genetic reporter,”
Nat. Photonics, 9
(4), 239
–246
(2015). https://doi.org/10.1038/nphoton.2015.22 NPAHBY 1749-4885 Google Scholar
C. Kim et al.,
“Sentinel lymph nodes and lymphatic vessels: noninvasive dual-modality in vivo mapping by using indocyanine green in rats—volumetric spectroscopic photoacoustic imaging and planar fluorescence imaging,”
Radiology, 255
(2), 442
–450
(2010). https://doi.org/10.1148/radiol.10090281 RADLAX 0033-8419 Google Scholar
L. Kou, D. Labrie and P. Chylek,
“Refractive indices of water and ice in the 0.65- to 2.5-μm spectral range,”
Appl. Opt., 32
(19), 3531
–3540
(1993). https://doi.org/10.1364/AO.32.003531 APOPAI 0003-6935 Google Scholar
S. L. Jacques,
“Optical properties of biological tissues: a review,”
Phys. Med. Biol., 58 R37
–R61
(2013). https://doi.org/10.1088/0031-9155/58/11/R37 PHMBA7 0031-9155 Google Scholar
J. Laufer, E. Zhang and P. Beard,
“Evaluation of absorbing chromophores used in tissue phantoms for quantitative photoacoustic spectroscopy and imaging,”
IEEE J. Sel. Top. Quantum Electron., 16
(3), 600
–607
(2010). https://doi.org/10.1109/JSTQE.2009.2032513 IJSQEN 1077-260X Google Scholar
L. V. Wang and H.-I. Wu, Biomedical Optics: Principles and Imaging, John Wiley and Sons, Inc., Hoboken, New Jersey
(2009). Google Scholar
E. Zhang, J. Laufer and P. Beard,
“Backward-mode multiwavelength photoacoustic scanner using a planar Fabry–Perot polymer film ultrasound sensor for high-resolution three-dimensional imaging of biological tissues,”
Appl. Opt., 47 561
–577
(2008). https://doi.org/10.1364/AO.47.000561 APOPAI 0003-6935 Google Scholar
M. Schweiger and S. Arridge,
“The Toast++ software suite for forward and inverse modeling in optical tomography,”
J. Biomed. Opt., 19 040801
(2014). https://doi.org/10.1117/1.JBO.19.4.040801 JBOPFO 1083-3668 Google Scholar
B. E. Treeby and B. T. Cox,
“k-Wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wave fields,”
J. Biomed. Opt., 15 021314
(2010). https://doi.org/10.1117/1.3360308 JBOPFO 1083-3668 Google Scholar
B. E. Treeby, E. Z. Zhang and B. T. Cox,
“Photoacoustic tomography in absorbing acoustic media using time reversal,”
Inverse Prob., 26
(11), 115003
(2010). https://doi.org/10.1088/0266-5611/26/11/115003 INPEEY 0266-5611 Google Scholar
D. Razansky and V. Ntziachristos,
“Hybrid photoacoustic fluorescence molecular tomography using finite-element-based inversion,”
Med. Phys., 34
(11), 4293
–4301
(2007). https://doi.org/10.1118/1.2786866 MPHYA6 0094-2405 Google Scholar
|