Fast unmixing of multispectral optoacoustic data with vertexcomponent analysis

X. Luís Deán-Ben a, Nikolaos C. Deliolanis b, Vasilis Ntziachristos a, Daniel Razansky a,n

a Institute for Biological and Medical Imaging, Technische Universität München and Helmholtz Zentrum München, Ingolstädter Landstraße 1,85764 Neuherberg, Germanyb Fraunhofer Institut für Produktionstechnik und Automatisierung IPA, Theodor-Kutzer-Ufer 1-3, 68167 Mannheim, Germany

a r t i c l e i n f o

Article history:Received 1 September 2013Received in revised form18 January 2014Accepted 30 January 2014

Keywords:Optoacoustic imagingPhotoacoustic imagingHyperspectral unmixingMolecular imaging

a b s t r a c t

Multispectral optoacoustic tomography enhances the performance of single-wavelength imaging interms of sensitivity and selectivity in the measurement of the biodistribution of specific chromophores,thus enabling functional and molecular imaging applications. Spectral unmixing algorithms are used todecompose multi-spectral optoacoustic data into a set of images representing distribution of eachindividual chromophoric component while the particular algorithm employed determines the sensitivityand speed of data visualization. Here we suggest using vertex component analysis (VCA), a method withdemonstrated good performance in hyperspectral imaging, as a fast blind unmixing algorithm formultispectral optoacoustic tomography. The performance of the method is subsequently compared witha previously reported blind unmixing procedure in optoacoustic tomography based on a combination ofprincipal component analysis (PCA) and independent component analysis (ICA). As in most practicalcases the absorption spectrum of the imaged chromophores and contrast agents are known or can bedetermined using e.g. a spectrophotometer, we further investigate the so-called semi-blind approach, inwhich the a priori known spectral profiles are included in a modified version of the algorithm termedconstrained VCA. The performance of this approach is also analysed in numerical simulations andexperimental measurements. It has been determined that, while the standard version of the VCAalgorithm can attain similar sensitivity to the PCA–ICA approach and have a robust and fasterperformance, using the a priori measured spectral information within the constrained VCA does notgenerally render improvements in detection sensitivity in experimental optoacoustic measurements.

& 2014 Elsevier Ltd. All rights reserved.

1. Introduction

Multispectral optoacoustic tomography (MSOT) extends the capa-cities of monochromatic (single-wavelength) imaging by making useof the spectral signatures of the various intrinsic tissue chromo-phores and extrinsically administered contrast agents present in theimaged tissue [1,2]. Thereby, whereas single-wavelength excitationprovides high-resolution mapping of the background absorption inbiological tissues for depths up to a few centimeters [3,4], a bettermapping of the specific absorbers present is attained by recordingimages at multiple wavelengths. In this way, it is possible to image e.g. the oxygenation level of blood or the distribution of externallyadministered photo-absorbing agents, which facilitates emergence ofa powerful set of functional and molecular optoacoustic imagingapplications [5–8].

Improvements in the quantification capacity and sensitivityof MSOT are the two main challenges presented by this techniquefor its applicability in molecular imaging. Quantitative MSOTof photoabsorbing probes requires the accurate determinationof their spatial distribution within the imaged tissue in relativeamounts or actual units of concentration, where imaging inarbitrary units can be sufficient if the system is properly calibrated.Quantification errors are originated in the reconstructed single-wavelength images mainly due to effects of the ultrasonic trans-ducer [9,10], acoustic attenuation [11,12] and optical attenuation[13,14]. If the single-wavelength images are properly corrected toavoid quantification errors, the multispectral processing proceduredetermines the sensitivity of MSOT, i.e., its capacity to imagelow concentrations of the probe. The optical absorption withina biological tissue is a linear combination of the absorptioncontributed by the different tissue chromophores and externalcontrast agents, whereas each absorbing substance has a char-acteristic absorption spectrum. In the simplest form, the distribu-tion of a given substance within a biological tissue can in principle

Contents lists available at ScienceDirect

journal homepage: www.elsevier.com/locate/optlaseng

Optics and Lasers in Engineering

http://dx.doi.org/10.1016/j.optlaseng.2014.01.0270143-8166 & 2014 Elsevier Ltd. All rights reserved.

n Corresponding author.E-mail address: dr@tum.de (D. Razansky).

Optics and Lasers in Engineering 58 (2014) 119–125

be imaged by exciting with an optical wavelength close to its peakabsorption [15]. However, this procedure is limited to those casesin which the absorption of the substance of interest is considerablyhigher than the background absorption. If the optical absorption ofthe different chromophores is comparable, the spectral profile ofeach must be taken into account in order to efficiently distinguishthe spatial distribution of the substance of interest using imagestaken at different optical wavelengths. This procedure is furtherreferred to as spectral unmixing [16].

The simplest spectral unmixing method consists in calculatingthe difference of the images taken at two wavelengths correspond-ing to high and low absorption of the substance of interest [17].This method is especially convenient for unmixing fluorescentdyes having a sharp drop in their optical absorption in a relativelynarrow spectral window. However, the detection sensitivity of thevarious substances can be substantially increased by consideringimages at more wavelengths. If the spectral signature of thedifferent components is known, which is usually the case as thebackground absorption is mainly due to blood (oxygenated anddeoxygenated hemoglobin), this information can subsequently beused in direct spectral inversion procedures such as least-squarefitting unmixing [6,18]. A different approach, the so-called blindunmixing, does not include any a priori spectral information butthe spectra are instead recovered as a part of the inversionprocedure. Blind unmixing in MSOT can be achieved by e.g.combining principal component and independent componentanalysis methods (PCA–ICA) [19], which was shown to performbetter in terms of sensitivity compared to the spectral fittingapproach. The reasons for this enhanced performance are mainlyrelated to errors in the reconstructed spectral profiles of theimages and can be multifarious. For instance, errors in thereconstruction due to wrong assumptions, wrong normalizationof the laser energy incident upon the tissue or the spectraldependence of the light attenuation might influence the spectralsignature of the optoacoustic tomographic reconstructions. Thecomputational time of the unmixing procedure has also becomean important factor with the emergence of fast (real-time)optoacoustic imaging systems [15,20–23], so that fast unmixingof multispectral data is essential for the visualization of spectrally-distinctive chromophores.

In this paper we suggest a different blind unmixing approachtermed vertex component analysis (VCA) [24] for unmixing ofmultispectral optoacoustic data. Its performance in terms ofsensitivity is compared with the previously reported PCA–ICAmethod. Since in this approach the spectral components aredetermined individually (separately), we further investigate theeffect of including some of the a priori information on the spectralsignatures of absorbers in the spectral inversion process, which isherein termed the constrained VCA or semi-blind unmixingapproach.

2. Theory

2.1. Spectral fitting

For each given laser wavelength λi, the reconstructed optoa-coustic image represents the spatial distribution of the opticalabsorbed energy Hðr; λiÞ. If m different substances contribute tothe absorption, Hðr; λiÞ in arbitrary units is given by

Hðr; λiÞ ¼ Uðr; λiÞϵ1ðλiÞC1ðrÞþUðr; λiÞϵ2ðλiÞC2ðrÞþ⋯þUðr; λiÞϵmðλiÞCmðrÞ; ð1Þ

with i¼ 1;2;…;n and n being the total number of wavelengthsemployed. Uðr; λiÞ represents the space and wavelength dependentlight fluence, ϵjðλiÞ and CjðrÞ are the wavelength-dependent molar

extinction coefficient for the j-th absorber and its spatial distribu-tion (concentration) respectively. Here, a homogeneous distribu-tion of the Grueneisen parameter was assumed. For each pixel k ofthe image, Eq. (1) can be expressed in a matrix form as

Hk ¼ ϵkCk; ð2Þwhere Hk is a column vector with elements HkðiÞ ¼Hðrk; λiÞ, Ck is acolumn vector with elements CkðjÞ ¼ CjðrkÞ and ϵk is a matrix withelements ϵkði; jÞ ¼ Uðrk; λiÞϵjðλiÞ. Generally, the spectral variation ofthe light fluence is also spatially dependent. Even if the illumina-tion profile on the surface of the object is kept constant fordifferent wavelengths, the changes in the optical attenuationcreate a wavelength dependence of the light fluence for deeperlocations. Thus, estimation of absorber concentration must bedone on a per-pixel basis by means of e.g. spectral fitting, whichconsists of a least square minimization of the form

Cksol ¼ argmin

Ck

‖Hkexp�ϵkCk‖22; ð3Þ

where Hkexp corresponds to the experimental measurement of the

column vector Hk.As a first order approximation, one can assume that the

spectral dependence of the light fluence is less significant thanthe spectral variation of the molar extinction coefficients of thedifferent substances, so that Eq. (2) can be approximated as

Hk �U ðrkÞϵCk; ð4Þwhere U ðrkÞ is the mean value of Uðrk; λiÞ for all the wavelengthsand ϵ is a matrix with elements ϵði; jÞ ¼ ϵjðλiÞ. In this case, Eq. (3) istransformed in arbitrary units to

U ðrkÞCksol ¼ argmin

Ck

‖Hkexp�ϵCk‖22: ð5Þ

The minimization in Eq. (5) can be done with the pseudoin-verse ϵþ of matrix ϵ for all the pixels of the image, i.e.,

CU;sol ¼ ϵþHexp; ð6Þbeing CU;sol and Hexp two matrices whose columns are U ðrkÞCk

soland Hk

exp respectively.As indicated in Eq. (5), the resulting images obtained by

spectral fitting are affected by light fluence variations, so that theymust be corrected in order to obtain a quantitative distribution ofthe absorber concentration. Even though several methods can beused to correct for the light fluence variations [13,25,26], this workis only concerned with the procedure to retrieve CU;sol.

2.2. Vertex component analysis

In the previous section, it was shown that the spectral depen-dence of the light fluence may induce errors in the unmixed imagesif the spectral variations in the images are directly fitted to the molarextinction coefficients of the absorbing substances, especially forabsorbers located deeper in tissue where spectrally-dependentlight attenuation may become significant. In practice, the spectralvariations of the reconstructed images may also be affected by otherfactors. For example, even if the signals for each wavelength arenormalized to the respective laser pulse energy, the illuminationprofile may also change with the laser wavelength due to unstable orotherwise low laser beam quality. Furthermore, errors in the tomo-graphic reconstructions may also influence the spectral variationsin the images. The blind unmixing methods, including VCA, areexpected to show improved sensitivity in those cases.

The VCA algorithm consists of a two step procedure [24]. First,the spectral signatures of the different components are calculatedwhereas distribution of the different components is subsequentlyobtained by spectral fitting. The spectral signatures correspond to

X. Luís Deán-Ben et al. / Optics and Lasers in Engineering 58 (2014) 119–125120

the columns of the matrix ϵVCA given by

ϵVCA ¼

ϵVCA;1ðλ1Þ ϵVCA;2ðλ1Þ ⋯ ϵVCA;mðλ1ÞϵVCA;1ðλ2Þ ϵVCA;2ðλ2Þ ⋯ ϵVCA;mðλ2Þ

⋮ ⋮ ⋱ ⋮ϵVCA;1ðλnÞ ϵVCA;2ðλnÞ ⋯ ϵVCA;mðλnÞ

0BBBB@

1CCCCA: ð7Þ

Once the matrix ϵVCA is determined, distribution of the differentcomponents is obtained as described in the previous section, i.e.,

CU;sol ¼ ϵþVCAHexp: ð8Þ

Two hypotheses are made in deriving the algorithm. The firstone assumes that there is at least one pixel of the image contain-ing a pure component, i.e., for every component j, there is aposition k in the reconstructed image that fulfills

ϵVCA;j ¼Hkexp; ð9Þ

where ϵVCA;j is the j-th column of the matrix ϵVCA. A maximumof n components can be obtained with the VCA algorithm, wherethe second hypothesis assumes that a linear combination of theabsorption at different wavelengths is constant for every pixelin the image, i.e., for every pixel k of the image, the followingcondition is fulfilled:

∑iaiHðrk; λiÞ ¼ 1; ð10Þ

where ai are wavelength-dependent coefficients. Fig. 1 illustratesthe principle of the VCA algorithm. Here it is considered that twoabsorbing components (with absorption spectra corresponding tovectors a! and b

!) are present and the sample is imaged at two

different wavelengths. The black dots correspond to a representa-tion of the absorbed energy at the imaged wavelengths for eachpixel of the image. We further assumed that the first approxima-tion above is verified, i.e., there are two pixels A and B whoseabsorption spectra correspond to the vectors a! and b

!respec-

tively. If the second assumption mentioned above is furtheraccurate, the pixels lie in a hyperplane defined by Eq. (10),represented as a dashed line in Fig. 1a. In general, this hypothesismay turn inaccurate in cases where the light fluence variations arestrong.

The procedure to obtain the spectral signatures works asfollows. First, a random direction is considered (m!1 in Fig. 1a).The vectors representing the absorption of the different pixels(black dots) are projected onto this direction while the pixelcorresponding to the maximum value of the projection is con-sidered (point A in Fig. 1a). The direction defined by this point(vector a! in Fig. 1a) is assigned to the spectral signature of thefirst component. Subsequently, a random direction normal to the

spectra of the already determined components is taken, so thatthe new component is obtained by projecting the absorption ofeach pixel onto the new direction as described above. In theexample shown in Fig. 1a, the spectral signature of the secondcomponent (vector b

!) is retrieved by taking a direction m!2

normal to a!. In Fig. 1a, the correct spectra of the compo-nents are retrieved. However, wrong spectra may be rendered insome cases. For example, the maximum value of the projectionalong direction m!1 in Fig. 1b corresponds to vector c!, which doesnot correspond to the spectrum of an actual component. This kindof errors may appear in case that two or more components with ahigh concentration are located in the same region or in case thatareas with a high light fluence (hot spots) are present.

2.3. Constrained vertex component analysis

In many practical applications, the purpose of the unmixing inMSOT is to determine the distribution of a specific substance withknown absorption spectrum. In particular, in molecular imagingapplications one often seeks after a specific bio-marker targetinga certain molecular pathway. In such a case, the spectral signatureof the absorption of the substance of interest can be a priorimeasured (for example with an absorption spectrometer). Thisinformation can be then included in the VCA algorithm by takingadvantage of the fact that the spectral signatures are determinedseparately. Thus, one can fix the first l components and use theVCA algorithm to only determine the spectral signatures of otherl�m components, so that the matrix ϵCVCA is given by

ϵCVCA ¼

ϵ1ðλ1Þ ⋯ ϵlðλ1Þ ϵVCA;lþ1ðλ1Þ ⋯ ϵVCA;mðλ1Þϵ1ðλ2Þ ⋯ ϵlðλ2Þ ϵVCA;lþ1ðλ2Þ ⋯ ϵVCA;mðλ2Þ

⋮ ⋮ ⋮ ⋮ ⋱ ⋮ϵ1ðλnÞ ⋯ ϵlðλnÞ ϵVCA;lþ1ðλnÞ ⋯ ϵVCA;mðλnÞ

0BBBB@

1CCCCA: ð11Þ

This variation of the suggested algorithm, which is hereintermed constrained VCA, may help to improve the outcomein cases where errors are produced due to the two simplifiedassumptions mentioned above.

3. Materials and methods

Numerical simulations and experimental measurements weredone in order to test the performance of the VCA and theconstrained VCA algorithms in MSOT imaging.

A numerical phantom was utilized in order to simulate theore-tical improvement of constrained VCAwith respect to the standardversion of the algorithm. The phantom was constructed using

Fig. 1. Principle of the vertex component analysis algorithm.

X. Luís Deán-Ben et al. / Optics and Lasers in Engineering 58 (2014) 119–125 121

cryosection image of a mouse in the kidney region [27], whichmimics realistic tissue morphology to be imaged by optoacoustictomography. For this, it was considered that the blood concentra-tion is given by the image in Fig. 2a while a random map wastaken to represent blood oxygenation levels (Fig. 2b). Also, twoinsertions containing AlexaFluor 750 (AF750) fluorochrome wereincluded. In this way, it was assumed that absorption of light arisessolely from the linear combination of optical absorption byoxygenated hemoglobin, deoxygenated hemoglobin and AF750.A uniform light fluence was considered throughout the specimen,so that optoacoustic image is directly represented by the distribu-tion of the optical absorption coefficient. Thereby, no wavelength-dependent light fluence was taken into account, which wouldprobably generate further distortion to the spectral variations ofthe images. Random noise was added to the images. The peakabsorption of AF750 was taken equivalent to the absorption ofblood at 797 nm, so that the optical absorptions of the threechromophores are in the same order. Two different scenarios wereconsidered. In the first scenario, the AF750 insertions were locatedin two regions with a relatively high blood concentration but

different oxygenation levels (diagonal crosses in Fig. 2b). In thesecond scenario, the insertions were in regions with low back-ground absorption (crosses Fig. 2b). Since the first assumption ofthe VCA algorithm is only valid when the background bloodconcentration at the location of the insertion is low, the perfor-mance of the algorithm is expected to be better for the firstscenario.

Experimental measurements on a post-mortem mouse were alsoperformed with two polyethylene tubings inserted subcutaneously.The tubings contained different chromophoric (fluorescent) sub-stances, namely AF750 and indocyanine green (ICG) dyes. The tubingmaterial is transparent and does not introduce multispectral noise.Furthermore, the inner and outer diameters are 0.6 mm and 1 mmrespectively, so that no significant distortion in the images isproduced due to acoustic mismatch. The concentration of the twoabsorbers was varied so that the peak absorption coefficient rangedbetween μa ¼ 2 cm�1 and μa ¼ 0:2 cm�1 as was also verified withan absorption spectrophotometer. The overall purpose was to experi-mentally compare performance of the VCA versus PCA–ICA algo-rithms in terms of sensitivity. Also, the absorption spectra of the twoabsorbers, as measured with the spectrophotometer, were used inthe constrained VCA algorithm in order to analyze the potentialimprovement in experimental measurements.

A whole-body optoacoustic scanner was used to perform theexperimental measurements [15]. The scanner consists of a tunableshort-pulsed (nanosecond) laser, in which a fiber bundle was usedin order to create a ring-type illumination around the imaged cross-section of the mouse. The laser wavelength was varied between700 nm and 900 nm in a 10 nm step to acquire multispectral data.The generated ultrasonic waves were detected with a phased-arraytransducer consisting of 64 elements covering an angle of 1721around the imaged cross-section. Each element of the array iscylindrically focused with a focal length of 40 mm and a centralfrequency of 5 MHz. The signals acquired for each wavelength wereaveraged 20 times and a band-pass filter with cut-off frequencies0.1 MHz and 7 MHz was applied to the recorded signals. Thesubsequent image reconstruction was performed with a model-based algorithm described in [28] and the images normalized withthe measured per-pulse energy of the laser.

4. Results

The simulation results for the two different scenarios presentedin Fig. 2b are shown in Fig. 2c,e and d,f. The images and spectrarendered with the VCA algorithm by considering 3 componentsare displayed in Fig. 2c and d. It is shown that a distorted spectrumis retrieved if the background blood absorption is high (firstcolumn). The unmixed image in that case is also affected bycross-talk noise. On the other hand, a more faithful spectrumand overall a better quality image are rendered when the back-ground blood concentration is low. The corresponding resultsobtained with the constrained VCA algorithm by including thespectral signature of AF750 are displayed in Fig. 2e and f.The cross-talk noise present in Fig. 2c is reduced, which indicatesa better performance of the constrained version of the VCA. It isimportant to take into account that the results obtained with bothalgorithms depend on directions along which the projections aredone, which are generally random. As a result, the representativecomponent may or may not be obtained, yet the constrainedVCA algorithm generally provides better results in numericalsimulations.

The single-wavelength tomographic reconstructions for themouse experiment described in Section 3 are displayed in Fig. 3overlaid with the full spectra of the two probes employed. Thelocations of the tubings are indicated with white arrows in the

0

1

).U.

A()

(

700(nm)

0

1

).U.

A()

(

700(nm)

0

1

).U.

A()

(

700(nm)

0

1

).U.

A()

(

700

800 900 800 900

800 900 800 900(nm)

SO2

0

1

Fig. 2. (a) The simulated numerical phantom. (b) Random distribution of bloodoxygenation levels as applied to the simulated data. The crosses and diagonalcrosses in (b) show the two scenarios for placement of the AlexaFluor 750absorbers. (c) and (d) Unmixing results attained with the standard version of thevertex component analysis algorithm. (e) and (f) Unmixing results for theconstrained version of the algorithm.

X. Luís Deán-Ben et al. / Optics and Lasers in Engineering 58 (2014) 119–125122

figure. The optoacoustic images were acquired at wavelengthscorresponding to the peak absorption of the probes, for which theabsorption coefficient was μa ¼ 2 cm�1 for both probes, as mea-sured with a spectrometer. Fig. 4 shows comparison of unmixingresults obtained with the combined PCA–ICA (first and second

column) versus the VCA algorithm (third and fourth column).In this case, four different concentrations of the probe were usedso that each row corresponds to a different value of opticalabsorption in the tubings. The peak absorptions measured with aspectrometer are indicated within the figure. Both algorithmsdepend on random initializations and on the number of compo-nents assumed, so that representative components cannot alwaysbe efficiently retrieved. The results shown in Fig. 4 correspond tothe most representative solutions obtained in both cases. For eachunmixed image, the signal-to-noise ratio (SNR) is estimated as themaximum value in the region marked with a white square in thefirst row of Fig. 4 divided by the standard deviation outside thisregion. For the first column the resulting values of the SNR are 11.5,6.0, 3.9 and 3.5, corresponding to absorption coefficients values of2 cm�1, 1 cm�1, 0.4 cm�1 and 0.2 cm�1 respectively. The equiva-lent values are 8.8, 9.2, 11.1 and 9.3 for the second column, 9.2, 7.9,5.2 and 3.8 for the third column and 9.1, 10.6, 9.0 and 9.7 for thefourth column. Generally, both algorithms present a similar perfor-mance in terms of sensitivity, whereas the computational effort forthe VCA algorithm is significantly lower. Specifically, MATLABimplementations of the VCA and the PCA–ICA algorithms for 21images with 200�200 pixels take approximately 0.2 and 1 srespectively on a Intel Core2 Duo CPU E8400 computer operating

4mm

0

1

700

).U.

A()

(

(nm) 4mm700800 900 800 900

(nm)

0

1

).U.

A()

(

Fig. 3. Optoacoustic tomographic images at 750 nm (a) and 800 nm (b). Theabsorption spectra of AF750 and ICG measured with a spectrometer are alsooverlaid. The locations of the tubings containing AF750 and ICG are indicated withwhite arrows in (a) and (b) respectively.

Fig. 4. Unmixing results for different concentrations of AF750 and ICG with PCA–ICA (columns one and two) and with VCA (columns three and four). Each row correspondsto a peak absorption of the probe indicated in the figure.

X. Luís Deán-Ben et al. / Optics and Lasers in Engineering 58 (2014) 119–125 123

at 3.00 GHz with 8 GB of RAM. It is also important to emphasizethat the VCA algorithm can be readily parallelized so that a moreefficient implementation of the code can further substantiallyaccelerate the unmixing procedure.

Finally, performance of the constrained VCA algorithm wastested experimentally. As a representative example, Fig. 5 showsthe results obtained with constrained VCA by assuming 6 compo-nents with the two first spectra (AF750 and ICG) being a prioriknown and used in the unmixing procedure. The absorption peakof the probes measured with a spectrometer is μa ¼ 2 cm�1 in thiscase. It is expected that the constrained VCA algorithm retrievesthe distribution of AF750 and ICG for the first two components,where the spectra are a priori fixed. However, the distribution ofthese two absorbers is retrieved for the sixth and fourth compo-nents respectively, i.e., the spectral variation of the reconstructedabsorbed energy fits better to the blindly-determined spectra thanto the theoretical ones and the constrained approach failsto retrieve the correct results. This limits the capability of theconstrained VCA algorithm to unmix probes with low concentra-tions in practical cases.

5. Discussion and conclusions

Derivation of efficient unmixing procedures for sensitive detec-tion of tissue bio-markers in the presence of strongly absorbingbackground is a challenging yet an important task for the devel-opment of efficient molecular imaging applications by means ofmultispectral optoacoustic tomography.

In this work, applicability of the vertex component analysisalgorithm in unmixing multispectral optoacoustic reconstructionshas been evaluated. Specifically, the performance of vertex com-ponent analysis in terms of sensitivity was compared to anotherblind unmixing procedure based on a combination of principalcomponent analysis and independent component analysis. Thepossibility of using a priori known spectra during the unmixingprocedure was also suggested in a variation of the algorithmtermed constrained vertex component analysis. The performanceof this suggested semi-blind approach was also evaluated innumerical simulations and experimental measurements.

The experimental results indicate that the overall sensitivityof the VCA method in optoacoustic unmixing is similar to that ofthe combination of PCA–ICA, whereas its computational time issubstantially lower, an important advantage when dealing with

large amounts of data (e.g. in three-dimensional imaging) or tryingto perform unmixing in real time [21,29]. A more accuratecomparison between the sensitivity of the methods would involveconsidering different regions in mice, different probes and differ-ent relative concentrations of the probe with respect to bloodbackground, as well as a statistical analysis on the interpretationof the results by independent users. We have also shownthat, whereas constrained VCA performs better in numericalsimulations, it turned less sensitive in experimental measure-ments where it was found that the spectral profile of thereconstructed absorbed energy fits better to the spectra of thecomponents which were determined blindly than to the absorp-tion spectrum measured with the spectrophotometer. This resultmight be attributed to the fact that the wavelength-dependentbeam profile of the laser and errors in the reconstructions havealso affected spectral signatures of probes located close to thesurface of the object. Another possible reason is the fact that somechromophores can present an intrinsic optoacoustic spectrumdiffering substantially from the absorption spectrum measuredwith a spectrophotometer [30]. Moreover, the spectral depen-dence of the light attenuation further affects the spectral signa-tures of absorbers located deeper in tissue while fluctuations inthe laser energy or otherwise motion artefacts may also contributeto the spectral variations in the reconstructions. Accurate model-ing of the wavelength-dependent optical attenuation is then animportant issue to address in future work. As was previouslyshown in [19], the PCA–ICA approach can generally achieve highersensitivity compared to spectral fitting, and the results presentedhere confirmed the better performance of blind unmixing. Recon-structions obtained with the constrained VCA are neverthelessalways easier to interpret and reproduce than those yielded withthe standard version of the algorithm, for which the correctcomponent may not be retrieved correctly depending on randominitial settings. The constrained approach may also be beneficialwhen the probe is mixed with other tissue chromophores e.g.haemoglobin, which is common in blood-pool agents. However,since the spectrum measured with the spectrophotometer maynot be reliable, an important next step is the accurate estimationof the intrinsic optoacoustic spectra of photoabsorbing probeswith a dedicated optoacoustic spectrometer. The determination ofthe intrinsic optoacoustic spectra may also improve the sensitivityof the constrained VCA algorithm.

In conclusion, the vertex component analysis algorithm can beemployed in multispectral optoacoustics to successfully unmix thedistribution of relatively low concentrations of spectrally-distinctbio-markers and optical probes. The fast computational time of thealgorithm anticipates its applicability for real-time unmixing andhandling of large multi-spectral datasets.

Acknowledgments

D.R. acknowledges support from the European Research Coun-cil under Starting Independent Researcher Grant ERC-SG-LS7(Dynamit). V.N. acknowledges support from European Unionproject FAMOS (FP7 ICT, Contract 317744). The authors thankSarah Glasl, Uwe Klemm and Florian Jurgeleit for excellenttechnical assistance.

References

[1] Wang LV, Hu S. Photoacoustic tomography: in vivo imaging from organelles toorgans. Science 2012;335(6075):1458–62.

[2] Beard P. Biomedical photoacoustic imaging. Interface Focus 2011;1(4):602–31.[3] Brecht HP, Su R, Fronheiser M, Ermilov SA, Conjusteau A, Oraevsky AA. Whole-

body three-dimensional optoacoustic tomography system for small animals.J Biomed Opt 2009;14(6):064007.

0

1

700 800 900

).U.

A()(

(nm)700 800 900

(nm)

0

1

).U.

A()(

700 800 900

(nm)

0

1

).U.

A()(

700 800 900(nm)

0

1

).U.

A()(

0

1

).U.

A()(

700 800 900(nm)

0

1

).U.

A()(

700 800 900(nm)

Fig. 5. Unmixing results obtained with the constrained VCA by considering6 components in which the first 2 spectral signatures (AF750 in (a) and ICG in(b)) are a priori fixed.

X. Luís Deán-Ben et al. / Optics and Lasers in Engineering 58 (2014) 119–125124

[4] Ma R, Distel M, Deán-Ben XL, Ntziachristos V, Razansky D. Non-invasivewhole-body imaging of adult zebrafish with optoacoustic tomography. PhysMed Biol 2012;57(22):7227–37.

[5] Zhang HF, Maslov K, Stoica G, Wang LHV. Functional photoacoustic microscopyfor high-resolution and noninvasive in vivo imaging. Nat Biotechnol 2006;24(7):848–51.

[6] Razansky D, Distel M, Vinegoni C, Ma R, Perrimon N, Koster RW, et al.Multispectral opto-acoustic tomography of deep-seated fluorescent proteinsin vivo. Nat Photon 2009;3(7):412–7.

[7] Wang L, Maslov K, Wang LV. Single-cell label-free photoacoustic flowoxigra-phy in vivo. Proc Natl Acad Sci 2013;110(15):5759–64.

[8] Deán-Ben XL, Razansky D. Adding fifth dimension to optoacoustic imaging:volumetric time-resolved spectrally enriched tomography. Light: Science &Applications 2014;3:e137 http://dx.doi.org/10.1038/lsa.2014.18.

[9] Wang K, Su R, Oraevsky AA, Anastasio MA. Investigation of iterative imagereconstruction in three-dimensional optoacoustic tomography. Phys Med Biol2012;57(17):5399–423.

[10] Queiros D, Deán-Ben XL, Buehler A, Razansky D, Rosenthal A, Ntziachristos V.Modeling the shape of cylindrically focused transducers in three-dimensionaloptoacoustic tomography. J Biomed Opt 2013;18(7):076014.

[11] Deán-Ben XL, Razansky D, Ntziachristos V. The effects of acoustic attenuationin optoacoustic signals. Phys Med Biol 2011;56(18):6129–48.

[12] Treeby BE. Acoustic attenuation compensation in photoacoustic tomographyusing time-variant filtering. J Biomed Opt 2013;18(3):036008.

[13] Jetzfellner T, Rosenthal A, Buehler A, Dima A, Englmeier KH, Ntziachristos V,et al. Optoacoustic tomography with varying illumination and non-uniformdetection patterns. J Opt Soc Am A 2010;27(11):2488–95.

[14] Cox B, Laufer JG, Arridge SR, Beard PC. Quantitative spectroscopic photo-acoustic imaging: a review. J Biomed Opt 2012;17(6):061202.

[15] Buehler A, Herzog E, Razansky D, Ntziachristos V. Video rate optoacoustictomography of mouse kidney perfusion. Opt Lett 2010;35(14):2475–7.

[16] Farkas DL, Du CW, Fisher GW, Lau C, Niu W, Wachman ES, et al. Non-invasiveimage acquisition and advanced processing in optical bioimaging. ComputMed Imaging Graphics 1998;22(2):89–102.

[17] Razansky D, Vinegoni C, Ntziachristos V. Multispectral photoacoustic imagingof fluorochromes in small animals. Opt Lett 2007;32(19):2891–3.

[18] Keshava N. A survey of spectral unmixing algorithms. Linc Lab J 2003;14(1):55–78.

[19] Glatz J, Deliolanis NC, Buehler A, Razansky D, Ntziachristos V. Blind sourceunmixing in multi-spectral optoacoustic tomography. Opt Express 2011;19(4):3175–84.

[20] Gamelin J, Maurudis A, Aguirre A, Huang F, Guo P, Wang LV, et al. A real-timephotoacoustic tomography system for small animals. Opt Express 2009;17(13):10489–98.

[21] Buehler A, Deán-Ben XL, Claussen J, Ntziachristos V, Razansky D. Three-dimensional optoacoustic tomography at video rate. Opt Express 2012;20(20):22712–9.

[22] Xiang L, Wang B, Ji L, Jiang H. 4-d photoacoustic tomography. Sci Rep2013;3(1):1113.

[23] Deán-Ben XL, Razansky D. Portable spherical array probe for volumetric real-timeoptoacoustic imaging at centimeter-scale depths. Opt Express 2013;21(23):28062–71.

[24] Nascimiento JMP, Dias JMB. Vertex component analysis: a fast algorithmto unmix hyperspectral data. IEEE Trans Geosci Remote Sens 2005;43(4):898–910.

[25] Shao P, Cox B, Zemp RJ. Estimating optical absorption, scattering andGrueneisen distributions with multiple-illumination photoacoustic tomogra-phy. Appl Opt 2011;50(19):3145–54.

[26] Ripoll J, Ntziachristos V. Quantitative point source photoacoustic inversionformulas for scattering and absorbing media. Phys Rev E 2005;71(3):031912.

[27] Sarantopoulos A, Themelis G, Ntziachristos V. Imaging the bio-distribution offluorescent probes using multispectral epi-illumination cryoslicing imaging.Mol Imaging Biol 2011;13(5):874–85.

[28] Deán-Ben XL, Ntziachristos V, Razansky D. Acceleration of optoacousticmodel-based reconstruction using angular image discretization. IEEE TransMed Imaging 2012;31(5):1154–62.

[29] Deán-Ben XL, Ozbek A, Razansky D. Volumetric real-time tracking of periph-eral human vasculature with GPU-accelerated three-dimensional optoacoustictomography. IEEE Trans Med Imaging 2013;32(11):2050–5.

[30] Laufer J, Jathoul A, Pule M, Beard P. In vitro characterization of geneticallyexpressed absorbing proteins using photoacoustic spectroscopy. Biomed OptExpress 2013;4(11):2477–90.

X. Luís Deán-Ben et al. / Optics and Lasers in Engineering 58 (2014) 119–125 125