Device and method for multispectral optoacoustic imaging

Abstract

A device and an according method for multispectral optoacoustic imaging of an object is provided comprising: an irradiation unit configured to irradiate the object with electromagnetic radiation at two or more different irradiation wavelengths the electromagnetic radiation having a time-varying intensity; a detection unit configured to detect acoustic waves generated in the object upon irradiating the object with the electromagnetic radiation at the different irradiation wavelengths; and a processing unit configured to reconstruct images of the object based on the detected acoustic waves generated in the object at each of the irradiation wavelengths and to determine a spatial distribution of at least one first concentration value, which relates to a concentration of at least one electromagnetic radiation absorber in the object.

Claims

1. A device for multispectral optoacoustic imaging of an object, the device comprising: a) an irradiation unit comprising at least one laser source configured to emit electromagnetic radiation including visible and/or near-infrared light, the irradiation unit configured to irradiate the object with the electromagnetic radiation at two or more different irradiation wavelengths, said electromagnetic radiation having a time-varying intensity; b) a detection unit comprising one or more ultrasound detectors, the detection unit configured to detect acoustic waves generated in the object upon irradiating the object with the electromagnetic radiation at the different irradiation wavelengths, wherein at least one electromagnetic radiation absorber in the object absorbs the electromagnetic radiation; and c) a processor configured to: reconstruct images of the object based on the detected acoustic waves generated in the object at each of the different irradiation wavelengths; and determine a spatial distribution of at least one first concentration value, which relates to a concentration of the at least one electromagnetic radiation absorber in the object, wherein the determination of the spatial distribution of the at least one first concentration value is based on the reconstructed images at the different irradiation wavelengths, at least one wavelength-dependent extinction coefficient of the at least one electromagnetic radiation absorber in the object, and a linear combination of at least two model spectra, the model spectra representing wavelength-dependent basis functions of a radiation fluence or wavelength-dependent basis functions of a normalized radiation fluence.

2. The device according to claim 1, wherein the processor is configured to determine spatial distributions of at least two first concentration values relating to concentrations of at least two electromagnetic radiation absorbers in the object, wherein the determination of the spatial distributions of the at least two first concentration values is based on the reconstructed images at the different irradiation wavelengths, the wavelength-dependent extinction coefficients of the at least two electromagnetic radiation absorbers in the object, and the linear combination of the at least two model spectra, and to derive a spatial distribution of a second concentration value from the determined spatial distributions of the at least two first concentration values.

3. The device according to claim 2, wherein the second concentration value relates to blood oxygenation saturation in a biological tissue.

4. The device according to claim 1, wherein the at least two model spectra have been determined by determining a set of spectral patterns for different depths in a medium and for different concentrations of the at least one electromagnetic radiation absorber in the medium, and applying a statistical procedure on the set of spectral patterns.

5. The device according to claim 4, the set of spectral patterns having been obtained from: analytical solutions or numerical solutions of a diffusion equation or a radiative transfer equation for the different irradiation wavelengths, Monte Carlo simulations of light propagation for the different irradiation wavelengths, and/or a set of experimental measurements for the different irradiation wavelengths.

6. The device according to claim 4, the statistical procedure applied on the set of spectral patterns comprising: a principal component analysis; a kernel principal component analysis; and/or a linear discriminant analysis.

7. The device according to claim 1, the linear combination of the model spectra corresponding to an addition of a mean model spectrum and a linear combination of a first number of further model spectra, wherein each of the further model spectra is multiplied by a model parameter.

8. The device according to claim 7, the irradiation unit being configured to irradiate the object with the electromagnetic radiation at a second number of the different irradiation wavelengths, wherein the second number is larger than or equal to the sum of the first number of further model spectra and a number of the first concentration values to be determined.

9. The device according to claim 1, further comprising a second electromagnetic radiation absorber in the object for absorbing the electromagnetic radiation such that there are two different electromagnetic radiation absorbers in the object and further comprising another first concentration value such that there are two first concentration values, wherein the determination of the spatial distribution of the two first concentration values, which relate to concentrations of the two different electromagnetic radiation absorbers in the object, comprises an inversion of a system of non-linear equations, wherein the inversion of the system of non-linear equations is performed simultaneously for a selected set of different positions distributed in an image domain.

10. The device according to claim 1, further comprising a second electromagnetic radiation absorber in the object for absorbing the electromagnetic radiation such that there are two different electromagnetic radiation absorbers in the object another first concentration value such that there are two first concentration values; and a second wavelength-dependent extinction coefficient such that there are two wavelength-dependent extinction coefficients, wherein the determination of the spatial distribution of the two first concentration values, which relate to concentrations of the two different electromagnetic radiation absorbers in the object, comprises an inversion of a system of non-linear equations given by
P(r, .sub.j) =(.sub.M(.sub.j)+.sup.P.sub.i=1m.sub.i.sup.r.sub.i(.sub.j)) (C.sub.1(r).sub.1(.sub.j)+c.sub.2(r).sub.2(.sub.j)), with j=1 . . . w, P(r, .sub.j) denoting the intensity of a reconstructed image of the object at a position r and an irradiation wavelength .sub.j, .sub.M(.sub.j)+.sup.P.sub.i=1 m.sub.i.sup.r.sub.i (.sub.j) denoting a linear combination of the wavelength-dependent model spectra .sub.M() and .sub.j(), with i=1. . . p, at an irradiation wavelength .sub.j, wherein m.sub.i.sup.r denotes model parameters at the position r, and c.sub.1(r) .sub.1(.sub.j)+c.sub.2(r) .sub.2(.sub.j) denotes a linear combination of the wavelength-dependent extinction coefficients .sub.1(.sub.j) and .sub.2(.sub.j) of the two electromagnetic radiation absorbers at an irradiation wavelength .sub.j, wherein c.sub.l(r) and c.sub.2(r) denote the two first concentration values at the position r.

11. A method for multispectral optoacoustic imaging of an object, comprising the following steps: a) irradiating the object with electromagnetic radiation at two or more different irradiation wavelengths, said electromagnetic radiation having a time-varying intensity; b) detecting acoustic waves generated in the object upon irradiating the object with the electromagnetic radiation at the different irradiation wavelengths; and c) absorbing the electromagnetic radiation by at least one electromagnetic radiation absorber in the object; d) reconstructing images of the object based on the detected acoustic waves generated in the object at each of the irradiation wavelengths; and e) providing instructions for execution by a processor that, when executed by the processor, cause the processor to determine a spatial distribution of at least one first concentration value, which relates to a concentration of the at least one electromagnetic radiation absorber in the object, wherein the determination of the spatial distribution of the at least one first concentration value is based on the reconstructed images at the different irradiation wavelengths, at least one wavelength-dependent extinction coefficient of the at least one electromagnetic radiation absorber in the object, and a linear combination of at least two model spectra, the model spectra representing wavelength-dependent basis functions of a radiation fluence or wavelength-dependent basis functions of a normalized radiation fluence.

Description

(1) The above and other elements, features, characteristics and advantages of the present invention will become more apparent from the following detailed description of preferred embodiments with reference to the attached figures showing:

(2) FIG. 1 an example of a device for multispectral optoacoustic imaging;

(3) FIG. 2 an example of an Eigenspectrum model of wavelength-dependent optical fluence: (a-d) Model spectra .sub.M(), .sub.1(), .sub.2() and .sub.3() at utilized wavelengths . (e) Mean Square Error (MSE) of the model on a training dataset for different model dimensionalities. (f-h) Values of parameters m.sub.1, m.sub.2 and m.sub.3 on the training dataset as a function of depth (y axis) and tissue oxygenation (x axis);

(4) FIG. 3 an example of an incorporation of constraints on the spatial characteristics of fluence parameters m.sub.1, m.sub.2 and m.sub.3 in the case of simultaneous inversion on a grid of points. (a) A directed graph on the grid of points enforces a constraint on the values of m.sub.2 model parameter with depth. (b) A non-directed graph on the grid of points penalizes great variations of the fluence parameters between neighbor points;

(5) FIG. 4 an example of an sO.sub.2 estimation error for all grid points using the Eigenspectrum method with simultaneous inversion on a grid of points and a non-negatively constrained linear unmixing approach.

(6) FIG. 5 (a) an MSOT image (one wavelength presented) of a tumor-bearing nude mouse, wherein the points in the image correspond to the initial grid selected for inversion. (b) Automatic refinement of the initial grid for maximizing the SNR in the pixels used for inversion. (c) Map of fluence parameter m.sub.2 in the convex hull of the inversion grid. (d) Estimation of tissue blood sO.sub.2 after fluence correction; and

(7) FIG. 6 a schematic overview on an example of processes and algorithms of a method for multispectral optoacoustic imaging.

(8) FIG. 1 shows an example of a device for multispectral optoacoustic imaging comprising an irradiation unit 10 which is configured to irradiate an object 1 with electromagnetic radiation at two or more different irradiation wavelengths .sub.j and a time-varying intensity.

(9) For example, the irradiation unit 10 may comprise at least one laser source which is configured to emit visible and/or near-infrared light, for example between 690 nm and 900 nm, preferably at 5 to 20 excitation wavelengths. Preferably, the emitted light is pulsed and/or amplitude- or frequency-modulated.

(10) By irradiating the object 1 with transient, in particular pulsed and/or modulated, electromagnetic radiation, acoustic waves are generated, which are detected by detection unit 20, which comprises one or more ultrasound detectors, e.g. one or more piezoelectric detectors. Preferably, the detection unit 20 comprises a two-dimensional or three-dimensional array of ultrasound detector elements by which the acoustic waves generated in the object 1 can be detected and spatially resolved.

(11) The object 1 can preferably be a biological tissue, like a part of a human or animal body. Depending on the particularities of the desired application, pre-clinical or clinical applications, a holder (not shown) can be provided for holding and positioning the object 1, e.g. a tissue sample or a small animal, relative to the irradiation unit 10 and the detection unit 20. Alternatively or additionally, the irradiation unit 10 and detection unit 20 can be integrated in a handheld device which can be manually placed onto the object 1 and/or moved relative to the object 1.

(12) In the detector unit 20 the detected acoustic waves are converted into corresponding detector signals which are forwarded to processing unit 30 which is configured to reconstruct images of the object 1 based on the detector signals that are obtained at each of the irradiation wavelengths .sub.j.

(13) The processing unit 30 is further configured to determine a spatial distribution of a relative concentration of at least one radiation absorber contained in the object 1 by considering the reconstructed images of the object 1 at different irradiation wavelengths .sub.j and the radiation fluence or a function of the radiation fluence in the object 1, e.g. the normalized radiation fluence, wherein the radiation fluence or the function of the radiation fluence is based on a linear model of the radiation fluence. According to the linear radiation fluence model, at least two model spectra representing wavelength-dependent basis functions of the radiation fluence or a function thereof are linearly combined.

(14) Preferably, a display unit 40, e.g. a computer monitor, is provided for displaying the reconstructed images and/or the determined spatial distribution of the relative concentration of the at least one radiation absorber contained in the object 1.

(15) In the following, preferred aspects of determining a spatial distribution of the relative concentration of the at least one radiation absorber contained in the object 1 will be elucidated in detail.

(16) An algorithmic solution for obtaining quantitative measures of blood oxygenation from multispectral optoacoustic images based on a new formulation of the quantification problem as a per-pixel non-linear spectral unmixing problem will be elucidated in more detail in section PROBLEM FORMULATION below.

(17) Preferably, a simple linear model of wavelength-dependent basis functions is created for capturing the spectral perturbations of the optical fluence in tissue using training data, which will be elucidated in more detail in section THE EIGENSPECTRUM MODEL OF OPTICAL FLUENCE below.

(18) Preferably, based on this model, a non-linear inversion scheme is formulated for obtaining quantitative blood sO.sub.2 estimates from MSOT images, which will be elucidated in section INVERSION below. The proposed inversion scheme further exploits prior knowledge on structural and spectral characteristics of the optical fluence for enhancing inversion accuracy.

(19) Finally the application of the novel algorithmic solution on in vivo multispectral optoacoustic images will be demonstrated in section APPLICATION ON TISSUE IMAGES below.

(20) Problem Formulation

(21) Optical Fluence and Acquiring Accurate sO.sub.2 Estimates

(22) With oxygenated (HbO.sub.2) and deoxygenated (HHb) hemoglobin being the main tissue absorbers, optoacoustic images of tissues P(r,) can be related to the concentrations of HbO.sub.2 and HHb through a non-linear relation, i.e.

(23) $\begin{array}{cc}P\left(r,\right)=C\left(r\right){.Math.\left(r\right).Math.}_2\frac{\left(r,\right)}{{.Math.\left(r\right).Math.}_2}.Math.\left(c_{\mathrm{HbO}2}\left(r\right).Math.{\mathrm{.Math.}}_{\mathrm{HbO}2}\left(\right)+c_{\mathrm{HHb}}\left(r\right).Math.{\mathrm{.Math.}}_{\mathrm{HHb}}\left(\right)\right).& \left(1\right)\end{array}$

(24) In Eq. (1), C(r) is a space dependent scaling factor that includes the effect of the Grneisen coefficient ((r)) and possibly additional physical phenomena (e.g. ultrasound attenuation) and system effects (e.g. the space and frequency dependent sensitivity of the ultrasound sensors). .sub.HbO2() and .sub.HHb() are the wavelength dependent molar extinction coefficients of hemoglobin, while c.sub.HbO2(r) and c.sub.HHb(r) are the associated concentrations at a position r. (r).sub.2 is the norm of the optical fluence across all wavelengths at a position r, while (r,)/(r).sub.2 is the normalized wavelength dependence of optical fluence at a specific position, a term that will also be referred to as (r,) in connection with the present invention.

(25) The space-only dependent factors C(r) and (r).sub.2 do not affect the estimation of sO.sub.2 (Eq. (2)) which is calculated as a ratio once the relative concentrations of HbO.sub.2 and HHb are known (c.sub.HbO2(r)=Const(r).Math.c.sub.HbO2(r) and c.sub.HHb(r)=Const(r).Math.c.sub.HHb(r), respectively).

(26) $\begin{array}{cc}{\mathrm{sO}}_2\left(r\right)=\frac{c_{\mathrm{HbO}2}^{}\left(r\right)}{c_{\mathrm{HbO}2}^{}\left(r\right)+c_{\mathrm{HHb}}^{}\left(r\right)}.& \left(2\right)\end{array}$

(27) An estimate of the wavelength dependence of the optical fluence (r,) is required for accurately extracting the relative values of c.sub.HbO2(r) and c.sub.HHb(r). However, (r,) is in general an unknown function.

(28) According to a preferred aspect of the invention, the wavelength dependence of optical fluence (r,) is not an arbitrary function of wavelength, but it is assumed to only lie within a limited subspace of spectral patterns within tissue. This aspect is based on the notion that the spectrum of optical fluence is a result of the cumulative light absorption by certain tissue absorbers, so the spectrum of optical fluence will always be related to the spectra of these absorbers through a complex relation.

(29) According to another preferred aspect, it is assumed that this subspace of spectral patterns can be accurately described using a low-dimensional linear model. In the section THE EIGENSPECTRUM MODEL OF OPTICAL FLUENCE the creation of such a model that consists of a linear combination of a small number of characteristic tissue optical fluence spectra, also referred to as Eigenspectra, will be described in detail.

(30) Non-linear Spectral Unmixing Formulation

(31) According to yet another aspect, it is assumed that it is possible to define a linear model that is based on a superposition of a mean spectrum .sub.M() and a number of optical fluence Eigenspectra .sub.i() to accurately characterize (r,) at each arbitrary position r within tissue, then the non-linear mixture model for MSOT images at a position r can be written as:

(32) $\begin{array}{cc}P\left(r,\right)=\left({}_M\left(\right)+\underset{i=1}{\overset{P}{.Math.}}{m}_i^r{}_i\left(\right)\right).\left(c_{\mathrm{HbO}2}^{}\left(r\right).Math.{\mathrm{.Math.}}_{\mathrm{HbO}2}\left(\right)+c_{\mathrm{HHb}}^{}\left(r\right).Math.{\mathrm{.Math.}}_{\mathrm{HHb}}\left(\right)\right).& \left(3\right)\end{array}$

(33) The solution for the P+2 unknowns (namely the P model parameters, m.sub.1 . . . P.sup.r and the relative blood concentrations c.sub.HbO2(r) and c.sub.HHb(r)) can then be given using an appropriate non-linear inversion approach and at least P+2 excitation wavelengths. As stated before, the relative blood concentrations c.sub.HbO2(r) and c.sub.HHb(r) are proportional to the originals c.sub.HbO2(r) and c.sub.HHb(r) with regard to a scaling factor, however this fact does not affect the computation of sO.sub.2, which is a ratiometric (relative) measurement.

(34) In order to achieve the aforementioned non-linear spectral un-mixing formulation, a linear model for the optical fluence, in the form of (.sub.M(), .sub.1(), . . . , .sub.P()) is established. The model is able to accurately describe the wavelength dependence of the optical fluence (r,), at any arbitrary position within tissue and independently of the surrounding optical properties and sO.sub.2 values.

(35) The Eigenspectrum Model of Optical Fluence

(36) The wavelength-dependent optical fluence in tissue may take a number of different spectral patterns. The variation of optical fluence per wavelength is not random but follows a pattern typically determined by tissue depth, tissue scattering, blood sO.sub.2 values and overall blood concentration. Preferably, by using Principal Component Analysis (PCA) on a set of different optical fluence spectral patterns, also referred to as training set, a linear model (.sub.M(), .sub.1(), . . . , .sub.P()) consisting of a small number of fluence Eigenspectra .sub.i() is determined. PCA is used for offering a minimum square error property in capturing the spectral variability of the optical fluence spectral patterns of the training set in a linear manner.

(37) By testing the Eigenspectrum model in a set of optical fluence spectral patterns produced in completely heterogeneous media with greatly varying and randomly distributed optical properties and oxygenation values, it was found that the Eigenspectrum model can accurately capture all such optical fluence spectral patterns with a very small mean square error (<2%), and it was confirmed that, independent of such conditions (i.e. optical property variation), the spectral patterns of optical fluence lie in a confined spectral subspace that can be accurately modeled by the Eigenspectrum model.

(38) In the following subsections, the light propagation models that are hereby used as well as an analytical description of the procedure for creating the Eigenspectrum model from training data are elucidated in detail.

(39) Light Propagation Models

(40) Preferably, a light propagation model based on the diffusion equation is employed, wherein the wavelength dependent optical fluence model (.sub.M(), .sub.1(), . . . , .sub.P()) is created using the solution of the diffusion equation for infinite homogeneous media, i.e.,
.sub.z=.sub.0exp(.sub.effz),(4)
where z denotes a distance from the source (tissue depth) and .sub.eff is the effective attenuation coefficient: .sub.eff={square root over (3.sub.(.sub.+.sub.s))}, .sub.s being the reduced scattering coefficient and .sub. the absorption coefficient.

(41) Alternatively, a finite element solution of the diffusion approximation for arbitrary structure of optical properties may be used.

(42) The Eigenspectrum Model

(43) Preferably, for producing a training dataset to generate the optical fluence Eigenspectrum model, simulations based on the analytical solution of the diffusion equation (Eq. (4)) are performed.

(44) For example, a number of optical fluence spectral patterns .sub.z,ox() are computed for wavelengths =700 nm to 890 nm with a step size of 10 nm, for depths from z=0 cm to z=1 cm with a step size of 0.143 mm, and for uniform hemoglobin distribution at a constant oxygenation, varying from ox=0% to ox=100% with a step of 5%. Constant optical properties of .sub.=0.3 cm.sup.1 (at the isosbestic point of hemoglobin at 800 nm) and .sub.s=10 cm.sup.1 (constant for all wavelengths) are used. The so computed 7021 optical fluence functions .sub.z,ox() were normalized to disregard any scaling .sub.z,ox()=.sub.z,ox()/.sub.z,ox.sub.2. The tissue depth of 1 cm was selected for relating the application to small animal optoacoustic imaging.

(45) The computed vectors .sub.z,ox () are used in the following as training data in the context of PCA to derive a linear model of Eigenspectra (a mean spectrum .sub.M() and three Eigenspectra .sub.1(), .sub.2(), .sub.3()) as exemplarily shown in FIG. 2(a)-(d).

(46) PCA captures the dominant variations of all spectral patterns used in the training dataset in the few first coefficients with a minimum mean square error. A model dimension of 3 Eigenspectra was selected for providing a rather simple model with a small associated error (see FIG. 2(e)).

(47) Regarding the parameters m.sub.1, m.sub.2 and m.sub.3 of the model on the training data set (see FIG. 2(f)-(h), respectively), it was surprisingly found that the values of m.sub.2 can be greatly associated with tissue depth, while the values of m.sub.1 can be greatly associated with the average levels of background tissue oxygenation. This indicates that the second Eigenspectrum .sub.2() is mainly associated with the modifications of optical fluence due to depth and the average optical properties of the surrounding tissue, while the first Eigenspectrum .sub.1() is associated with the spectral shape of optical fluence that relates to the average oxygenation of the surrounding tissue. Moreover, the highest and lowest values of these parameters on the training dataset can offer an indication of their limits [lim.sup.min.sub.i, lim.sup.max.sub.i] expected in tissue.

(48) A validation of the accuracy of the Eigenspectrum Model (.sub.M(), .sub.1(), . . . , .sub.P()) over simulation-based optical fluence spectra created in arbitrary tissues yielded a very small model error, supporting the inventive approach according to which a simple linear model with only a few, e.g. three, Eigenspectra can capture the spectral variability of (r,) in complex tissue structures, independently of the optical properties, depth and blood oxygenation.

(49) Similarly, the validity and accuracy of the Eigenspectrum model to capture optical fluence spectral patterns present in tissue was experimentally investigated by means of measurements from small animals in vivo and ex vivo, wherein the optical fluence in tissue was measured by inserting a reference chromophore with well characterized spectrum in deep tissue. The result of this investigation yielded that the Eigenspectrum model can capture realistically the spectral variability of optical fluence in tissue, even in the case of noisy experimental data obtained in vivo.

(50) Inversion

(51) Given the Eigenspectrum model for the optical fluence, the quantitative estimation of blood sO.sub.2 using MSOT reduces to solving Eq. (3) using a non-linear optimization approach. Inversion can be performed independently for a single pixel in the image or simultaneously for a grid of pixels distributed appropriately in the image domain. In the second case additional constraints on the values of the Eigenspectrum model parameters can be imposed, based on the spatial and spectral characteristics of optical fluence.

(52) In the following, three different inversion approaches are presented, namely (1) a per-pixel inversion using local or global optimization, (2) a constrained optimization on a grid of points, where the spatial characteristics of the optical fluence are further modelled, and (3) another constrained optimization, wherein both spatial and spectral characteristics of the optical fluence are taken into account for enhancing inversion stability.

(53) (1) Single Pixel Inversion

(54) In the single pixel inversion case, the optimization problem for a specific point r formulates as the minimization of the function in Eq. (5) over the values of m.sub.1.sup.r, m.sub.2.sup.r, m.sub.3.sup.r, c.sub.HbO2(r) and c.sub.HHb(r), constrained according to Eq. (6).

(55) $\begin{array}{cc}f={.Math.\begin{array}{c}P\left(r,\right)-\left({}_M\left(\right)+\underset{i=1}{\overset{P}{.Math.}}{m}_i^r{}_i\left(\right)\right).Math.\\ \left(c_{\mathrm{HbO}2}^{}\left(r\right).Math.{\mathrm{.Math.}}_{\mathrm{HbO}2}\left(\right)+c_{\mathrm{HHb}}^{}\left(r\right).Math.{\mathrm{.Math.}}_{\mathrm{HHb}}\left(\right)\right)\end{array}.Math.}_2& \left(5\right)\\ {\lim }_i^{\min }<m_i^r<{\lim }^{\max },c_{\mathrm{HbO}2}^{}\left(x\right)0,c_{\mathrm{HHb}}^{}\left(x\right)0.& \left(6\right)\end{array}$

(56) Using a constrained local optimization algorithm to solve the non-linear optimization problem of Eq. (5) and (6) in the case of a simulated data set it can be observed that the problem is non-convex and thus the inversion converges accurately for only a small percentage of pixels. Accordingly, the optimization problem is a so-called ill-posed problem. Preferably, in order to further improve the accuracy of the inversion and to achieve a unique solution, further constraints are imposed on the values of the fluence parameters m.sub.i.sup.r.

(57) (2) Inversion on a Grid of Points Incorporating Spatial Fluence Characteristics

(58) Preferably, the spatial characteristics of optical fluence are further exploited for overcoming the ill-posed nature of the optimization problem. In contrary to tissue absorption which can vary arbitrarily, the optical fluence is bound to vary smoothly in space due to the nature of diffuse light propagation.

(59) According to a preferred aspect of the Eigenspectrum model inversion, such a priori information is incorporated by attempting simultaneous inversion on a grid of points in the image domain, and penalizing large variation of the fluence model parameters m.sub.1.sup.r, m.sub.2.sup.r and m.sub.3.sup.r between neighbor pixels. This is preferably achieved through incorporation of appropriate Lagrange multipliers .sub.i to the non-linear inversion function for constraining the variation of one or all model parameters. The values of the Lagrange multipliers can be optimized using cross-validation on simulated data sets. Additionally, since the second Eigenspectrum is strongly associated with depth, the m.sub.2 model parameter can be constrained to take lower values as tissue depth increases (in accordance with FIG. 2(g)).

(60) FIG. 3 illustrates the incorporation of constraints on the spatial characteristics of fluence parameters m.sub.1.sup.r, m.sub.2.sup.r and m.sub.3.sup.r, in the case of simultaneous inversion on a grid of points, wherein a directed graph on the grid of points enforces a constraint on the values of m.sub.2 model parameter with depth (see FIG. 3(a)) and a non-directed graph on the grid of points penalizes great variations of the fluence parameters between neighbor points (see FIG. 3(b)).

(61) Assuming a circular grid of P arcs and L radial lines (see FIG. 3) with a total of PL points r.sub.p,l, and let the vector m.sub.i=[m.sub.i.sup.r1,1, m.sub.i.sup.r1,2, . . . , m.sub.i.sup.r1,L, m.sub.i.sup.r2,1, . . . , m.sub.i.sup.rp,l, . . . , m.sub.i.sup.rP,L] correspond to the values of the optical fluence parameter i (i=1 . . . 3) over all such points, the new inverse problem is preferably defined as the minimization of function f.sub.grid defined in Eq. (7) under the constraints defined in Eq. (8).

(62) $\begin{array}{cc}{f}_{\mathrm{grid}}=\underset{i}{.Math.}f\left(r_i\right)+{}_1{.Math.{\mathrm{Wm}}_1.Math.}_2+{}_2{.Math.{\mathrm{Wm}}_2.Math.}_2+{}_3{.Math.{\mathrm{Wm}}_3.Math.}_2& \left(7\right)\\ {\lim }_i^{\min }<m_i^{r_k}<{\lim }_i^{\max },i,k,m_2^{r_{p+1,1}}<m_2^{r_{p,l}},m_2^{r_{p+1,1+1}}<m_2^{r_{p,l}},m_2^{r_{p+1,1-1}}<m_2^{r_{p,l}}{c}_{\mathrm{HbO}2}\left(r_k\right)0,k,c_{\mathrm{HHb}}\left(r_k\right)0,k.& \left(8\right)\end{array}$

(63) In Eq. (7), W is the weighted connectivity matrix corresponding to grid of points assumed. Each matrix element corresponds to a pair of grid points r.sub.p1,l1 r.sub.p2,l2 and is zero if the points are not directly connected or inverse proportional to their distance (w(r.sub.p1,l1 r.sub.p2,l2)1/r.sub.p1,l1r.sub.p2,l2) if the points are connected.

(64) Simultaneous inversion in a grid of points has surprisingly turned out to significantly enhance inversion stability.

(65) (3) Inversion on a Grid of Points Incorporating Spatio-Spectral Fluence Characteristics

(66) For further enhancing the inversion stability, additional constraints can be imposed to the values of fluence parameters based on prior experimental measurements of the optical fluence in similar conditions or based on alternative prior knowledge that can be derived directly from the reconstructed multispectral optoacoustic images. Two possible examples of enforcing such constraints in the context of the grid inversion are described in the following:

(67) Global Constraints on m.sub.1

(68) For instance, in certain cases of in vivo imaging, it can be preferably assumed that the average tissue oxygenation is higher than 50%, and thus the model parameter m.sub.1 will obtain values lower or equal to zero as indicated by FIG. 2(f). This preferred assumption has been verified by multiple controlled in vivo experiments (where the animal is breathing 100% O.sub.2 gas and imaging is performed in the brain or in the abdomen), where the measured fluence parameter m.sub.1 was repeatedly found negative. Under this condition, and by constraining the values of m.sub.1.sup.r to be negative in all grid points excellent convergence rate can be achieved also in the case of high levels superimposed noise. It is noted that no assumptions are made on the specific oxygenation values of the grid points, but only on the average tissue oxygenation.

(69) Local Constraints on the Model Parameters

(70) Alternatively, if no prior fluence measurements are available some constraints for the model parameter m.sub.1 can be extracted from the available data themselves. By performing a first linear unmixing approach on the data a first estimation map of sO.sub.2 levels can be obtained. It is noted that this sO.sub.2 map is incrementally erroneous with tissue depth, however it can serve as a first approximation.

(71) Using this sO.sub.2 map and average tissue optical properties (e.g. .sub.=0.2 cm.sup.1 at 800 nm and .sub.s=10 cm.sup.1) the optical fluence within tissue can be simulated and a first estimate on the values of one or all model parameters m.sub.1.sup.r, m.sub.2.sup.r and m.sub.3.sup.r (i.e. {circumflex over (m)}.sub.1.sup.r, {circumflex over (m)}.sub.2.sup.r, {circumflex over (m)}.sub.3.sup.r) can be extracted, where r corresponds to the positions of the grid points. In the following, the optimization problem of Eq. (7) to (8) is solved, with the only difference that the values of one or all of the model parameters m.sub.1.sup.r, m.sub.2.sup.r and m.sub.3.sup.r are constrained to lie in a region close to the initial estimates {circumflex over (m)}.sub.1.sup.r, {circumflex over (m)}.sub.2.sup.r, {circumflex over (m)}.sub.3.sup.r. It is noted that these regions need to be incrementally larger with depth since in deep tissue the original estimates deviate usually from the true values. Such an optimization approach incorporating prior information for the model parameters corresponding to each grid pixel can also be formulated in the context of Bayesian inversion. By constraining the problem in such a way the inversion stability is also greatly enhanced also in cases of high random noise.

(72) FIG. 4 relates to an example of an inversion of Eq. (7) and (8) on a grid of 64 points in the case of a simulated dataset and shows the error in the estimation of local blood sO.sub.2 for all grid points when using the inventive method (blue curve) and when using simple linear non-negatively constrained unmixing (red curve). As apparent from the comparison, the inventive method is associated with a very small error in the estimation of blood oxygenation as compared to linear unmxing which is frequently used as an approximation.

(73) Application on Tissue Images

(74) In the following, the inventive Eigenspectrum method is exemplarily demonstrated in the case of experimental MSOT data. A nude mouse bearing an orthotopic Mammarian 4T1 tumor was imaged in the tumor area in eight excitation wavelengths from 690 nm to 900 nm.

(75) FIG. 5(a) depicts the anatomical optoacoustic image acquired at wavelength 900 nm. The points on the image correspond to the initial grid of pixels that were selected for the application of the non-linear inversion scheme described above. The grid was selected to cover the tumor area. The grid of pixels is in the following automatically adjusted to avoid using pixels of low image intensity (and thus low SNR) in the inversion.

(76) After applying the non-linear inversion scheme on the refined grid of points (see FIG. 5(b)), estimates for fluence parameters m.sub.1.sup.r, m.sub.2.sup.r and m.sub.3.sup.r are obtained for each grid point r. In the following, the fluence parameters for each pixel in the convex hull of the grid is obtained using cubic interpolation, thus resulting in fluence maps for the image sub-region.

(77) FIG. 5(c) presents such a map for fluence parameter m.sub.2. The wavelength dependence of optical fluence is computed for each pixel within these maps as in (r,)=.sub.M()+m.sub.1.sup.r.sub.1()+m.sub.2.sup.r.sub.2()+m.sub.3.sup.r.sub.3(), where .sub.i() is the ith fluence Eigenspectrum.

(78) Finally, a fluence-corrected MSOT image is obtained after diving the original image P(r,) with the normalized wavelength dependent optical fluence (r,) at each position r and wavelength : P.sup.corr(r,)=P(r,)/(r,). After fluence correction, non-negatively constrained linear spectral unmixing with the spectra of oxyand deoxy-hemoglobin can result in an estimation of tissue blood oxygenation.

(79) FIG. 5(d) presents the estimated tissue oxygenation in the area of the tumor as computed using the Eigenspectrum method. Tissue oxygenation is overlaid to the anatomical image with pseudocolor, represented in gray scale.

(80) FIG. 6 shows an overview on an example of processes and algorithms of a method for multispectral optoacoustic imaging. In present example, oxygenated (HbO.sub.2) and deoxygenated (HHb) hemoglobin being the main tissue absorbers in the object 1.

(81) In a first step, MSOT images are reconstructed at distinct irradiation wavelengths .sub.j based on detected acoustic waves which are generated in the object 1 (see FIG. 1) upon irradiating the object 1 with transient electromagnetic radiation at the irradiation wavelengths .sub.j.

(82) In a second step, the reconstructed MSOT images undergo a point selection process, wherein reliable pixels r.sub.i in the image domain are selected in a semiautomatic way or by user annotation.

(83) Optoacoustic image intensities P(r.sub.i,.sub.j) at the selected pixels r.sub.i for all irradiation wavelengths .sub.j are processed in a third step, wherein fluence parameters m.sub.k.sup.ri and concentration values, e.g. sO.sub.2 values, are computed for each pixel r.sub.i in the grid using the image intensities P(r.sub.i,.sub.j) and a pre-determined model of the fluence in the object, wherein the model is based on a linear combination of model spectra .sub.k(). This is performed by means of an inversion algorithm solving the following equation:

(84) $P\left(r_1,{}_1\right)=\left({}_M\left({}_1\right)+\underset{k=1}{\overset{K}{.Math.}}{m}_k^{r_i}{}_k\left({}_1\right)\right)\left({c^{}}_{\mathrm{HbO}2}\left(r_1\right){\mathrm{.Math.}}_{\mathrm{HbO}2}\left({}_1\right)+{c^{}}_{\mathrm{HHb}}\left(r_1\right){\mathrm{.Math.}}_{\mathrm{HHb}}\left({}_1\right)\right)$$P\left(r_1,{}_2\right)=\left({}_M\left({}_2\right)+\underset{k=1}{\overset{K}{.Math.}}{m}_k^{r_i}{}_k\left({}_2\right)\right)\left({c^{}}_{\mathrm{HbO}2}\left(r_1\right){\mathrm{.Math.}}_{\mathrm{HbO}2}\left({}_2\right)+{c^{}}_{\mathrm{HHb}}\left(r_1\right){\mathrm{.Math.}}_{\mathrm{HHb}}\left({}_2\right)\right)$$\mathrm{.Math.}$$P\left(r_1,{}_{\mathrm{Nw}}\right)=\left({}_M\left({}_{\mathrm{Nw}}\right)+\underset{k=1}{\overset{K}{.Math.}}{m}_k^{r_i}{}_k\left({}_{\mathrm{Nw}}\right)\right)\left({c^{}}_{\mathrm{HbO}2}\left(r_1\right){\mathrm{.Math.}}_{\mathrm{HbO}2}\left({}_{\mathrm{Nw}}\right)+{c^{}}_{\mathrm{HHb}}\left(r_1\right){\mathrm{.Math.}}_{\mathrm{HHb}}\left({}_{\mathrm{Nw}}\right)\right)$

(85) The calculated fluence parameters m.sub.k.sup.ri (i=1 to Np, Np being the number of selected pixels, k=1 to K, K being the number of model spectra) can be used to infer the wavelength dependence of the optical fluence in the whole area corresponding to the convex hull of the grid by means of interpolation, i.e. the m.sub.k values are interpolated from m.sub.k.sup.ri for all intermediate pixels in the convex hull of the grid of points.

(86) Subsequently, the wavelength dependence of the optical fluence is computed for each pixel within the convex hull as in (r,)=.sub.M()+m.sub.1.sup.r.sub.1()+m.sub.2.sup.r .sub.2()+m.sub.3.sup.r.sub.3()+ . . . .

(87) A fluence-corrected MSOT image is obtained after dividing the originally reconstructed MSOT image P(r,) with the normalized wavelength-dependent optical fluence (r,) at each position r and wavelength : P.sup.corr(r,)=P(r,)/(r,).

(88) After fluence correction, non-negatively constrained spectral unmixing within the spectra of oxy- and deoxy-hemoglobin result in the determination of a particularly accurate sO.sub.2 map.

(89) Preferably, the model of the fluence, in particular the model spectra .sub.k(), is established independently, in particular in advance, of the current measurements by means of a model creation algorithm, wherein base functions of the wavelength dependence of the optical fluence in the object are created for the excitation wavelengths .sub.j employed. Preferably, in distinction to previous fluence modeling methods, only the wavelength dependence of optical fluence is modeled but not the actual fluence in absolute terms. Moreover, modeling is performed by applying a statistical procedure on training data and not via a spatial light propagation model that is fitted to the MSOT images.

SUMMARY

(90) The invention relates to a novel solution for the problem of quantitative multispectral optoacoustic imaging of tissue oxygenation. According to preferred aspects of the invention, a simple linear model for successfully capturing the spectral variability of optical fluence with wavelength in tissue was created and validated. Using this model, a non-linear inversion scheme for correcting for the optical fluence and obtaining accurate quantitative estimates of tissue blood oxygenation was formulated. Formulating the quantification problem of MSOT as a per-pixel non-linear unmixing problem instead of a model-based optical property estimation approach allows for greater flexibility and easier applicability in experimental images of tissue. The inventive methodology does not assume full and accurate spatial information, and does not rely on absolute image values that are prone to error due a multitude of phenomena that are often not modelled in image reconstruction (e.g. spatial sensitivity field of the ultrasound arrays, ultrasound attenuation and scattering, Grneisen coefficient). In this respect the inventive method is readily applicable in in vivo MSOT images. Furthermore, since the presented inverse problem is of low dimensionality (typically 300-500 unknown parameters for the case of grid inversion), the computational cost is rather low. Specifically, the required computational time for the inversion of a grid of 64 points using an Intel Core i7-4820K can be less than 1 minute.

(91) Based on a validation of the inventive approach, a rather simple 3-dimensional model has turned out sufficient for accurately capturing the spectral variability of optical fluence, independently of the optical properties and the physiology of the surrounding tissue.

(92) Basically, the absorption of water, lipids and melanin also present in tissue may be ignored under an appropriate wavelength selection, assuming that the influence of hemoglobin is much more prevalent in the near-infrared region. Moreover, the wavelength dependence of scattering may be ignored. However, according to preferred embodiments of the invention, the influence of such parameters could also be incorporated similarly for optimizing the Eigenspectrum model.

(93) Despite the aforementioned simplifying assumptions, a model fit in the case of in vivo and ex vivo experimental optoacoustic images of mice indicated a very promising match between theory and experimental reality.

(94) Because the inverse problem, in its simplest form (per pixel inversion), is non-convex and ill-posed, in a per pixel inversion case different combinations of optical fluence and local oxygenation may produce an almost identical result after multiplication. Thus, the nonlinear inversion may converge to a wrong solution unless additional restrictions are applied.

(95) Therefore, it is particularly preferred to incorporate spatio-spectral fluence characteristics using available prior knowledge into a more complex inversion scheme. In this way, inversion stability was achieved even in cases where the spectra are corrupted by random noise. The preferred inversion scheme provides accurate estimates of tissue oxygenation with far lower estimation error than the one of linear spectral un-mixing, that is typically used today in experimental studies as an approximation.