Bioluminescence tomography reconstruction based on multitasking Bayesian compressed sensing

Info

Publication number: 20170148193
Type: Application
Filed: Jul 21, 2015
Publication Date: May 25, 2017
Inventors: Jinchao Feng (Beijing), Kebin Jia (Beijing), Huijun Wei (Beijing)
Application Number: 15/108,394

Abstract

Implementations of the present disclosure relate to methods for reconstruction for bioluminescence tomography based on a method of multitask Bayesian compressed sensing in the field of medical image processing. The method includes the following operations. Firstly the high order approximation model is used to model the law of light propagation in biological tissues, then the inner-correlation among multispectral measurements is researched based on multitask learning method and incorporated into a reconstruction algorithm of bioluminescence tomography as prior information to reduce ill-posedness of BLT reconstruction, and then on this basis, three-dimensional reconstruction of bioluminescent source is realized. Compared with other reconstruction algorithms for BLT, the correlation among multispectral measurements is incorporated into the disclosure and the ill-posedness of BLT reconstruction is reduced. The bioluminescent source can be reconstructed and located accurately using the proposed algorithm, and computational efficiency can be greatly improved.

Description

Description

CROSS REFERENCE TO RELATED APPLICATION

This application is a national stage application of International application number PCT/CN2015/084555, filed Jul. 21, 2015, titled “Bioluminescence tomography reconstruction based on multitasking Bayesian compressed sensing,” which claims the priority benefit of Chinese Patent Application No. 201510397424.0, filed on Jul. 8, 2015, which is hereby incorporated by reference in its entirety.

TECHNICAL FIELD

The disclosure belongs to the field of medical image processing, and more particularly to a reconstruction method for bioluminescence tomography based on a method of multitask Bayesian compressed sensing.

BACKGROUND

Optical molecular imaging is a rapidly developing technique of molecular imaging, which combines optical process with certain molecular properties. Optical molecular imaging analyzes and processes bioluminescence or excited fluorescence in objectives, and studies them qualitatively and quantitatively.

Representative imaging methods of optical molecular imaging technology are fluorescence imaging and bioluminescent imaging (BLI). Both of them belong to two-dimensional luminescent imaging technology. Although the technology is convenient and simple, it has limitations in application, especially for imaging depth, two-dimensional fluorescence image that may not reflect the information of light source depth and hard to quantify it. The two technologies mentioned above can only reflect projection information of fluorescent probes in organisms at a certain angle, which is generated by superposition of multiple depths signals. So, two-dimensional imaging technology has a very low resolution.

BLT developed based on two-dimensional bioluminescent imaging technology has become an important branch of optical molecular imaging because BLT can reflect the information of signal depth. Bioluminescence tomography technology works not by excitation of an external light source, but by a biochemical luminescence reaction in organisms. Photons generated in organisms propagate in a certain way in the biological tissue and constantly interact with it such as to reach body surface. Finally, distribution of the bioluminescent source in a small animal can be reconstructed using the bioluminescent images acquired on the surface of biological tissue by a high sensitivity detector, which reveals activity laws of molecules in vivo in essence.

The photons in biological tissue do not transmit along a straight line, but experience a lot of scattering process, which leads to BLT inverse problems: a highly ill-posed problem in mathematics. And small measurement disturbance from outside will bring great changes to reconstruction results. The researchers did a lot of work to reduce the ill-posedness usually by providing different prior knowledge to the problem from different angles.

The existing methods reducing ill-posedness are mostly based on multiplying spectral information and sparsity of light sources to study BLT reconstruction methods. However, the correlation among multi-spectra is not considered. In mathematics, the method of multi-spectra increases the known quantity of the equation to be solved, which can improve solving results. For BLT, in the method based on multi-spectra, the light source is reconstructed by the data obtained using multiple filters from different wavelengths, while the light source is remained unchanged. In fact, measurements among multi-spectra are not independent, but interrelated, and their common task is the reconstruction of the bioluminescent source. Therefore, the quality of reconstructed BLT images will be likely improved if the correlation among multi-spectra is considered into reconstruction. Based on this, a BLT reconstruction method based on inner-correlation of multi-spectra is proposed in the disclosure.

The disclosure is based on a multitask learning method. Firstly a high order approximation model is used to model the law of light transmission in the biological tissue, then inner-correlation among multi-spectra is researched based on the multitask learning method, and finally on this basis three-dimensional reconstruction of the bioluminescent source is realized.

SUMMARY

To solve the above problems existing in reconstruction algorithms for bioluminescence tomography, the disclosure proposes a method based on multitask learning. A high order approximation model is used to model the law of light propagation in the biological tissue, and inner-correlation among multi-spectra is researched based on multitask learning method, and finally on this basis three-dimensional reconstruction of a bioluminescent source is realized.

The technical scheme adopted in the disclosure is: firstly, a high order approximation model is used to model the law of light propagation in the biological tissue, and inner-correlation among multi-spectra is researched based on multitask learning method and incorporated into the reconstruction algorithm as prior information to reduce the ill-posedness of BLT reconstruction, finally on this basis three-dimensional reconstruction of the bioluminescent source is realized. Specifically, the following steps are provided as follow.

Step one, defining the problems and initialization—obtaining values of P spectral bands at M measurement points, and setting gamma prior distribution α₀with the shape parameter a and scale parameter b, wherein α₀⁻¹is the variance of emission light measurements Φ(τ_i) with Gaussian distribution.

Two models can be adopted to simulate the law of light propagation in a biological tissue: radiative transfer equation and diffusion equation. Radiative transfer equation is a mathematical model describing accurately the law of light propagation, but the solution of RTE is very difficult even under a very simplified condition, and the amount of calculation is very large. Diffusion equation is the low-order approximation of the radiative transfer equation, and the solution of which is simple but the precision is low. Furthermore, the application of the diffusion equation must meet conditions of biological tissues with strong scattering and low absorption. It is not suitable for the condition that the light source is located on the boundaries of biological tissues, which limits the application of diffusion equation in the biomedical fields. Based on this, simplified spherical harmonics approximation (SP_N) equations are adopted in the disclosure. The model may describe accurately the light propagation in a biological tissue with strong scattering and low absorption, at the same time, it has a modest computational complexity. Instead of diffusion approximation transfer equation, the model is used to accurately simulate the law of the light propagation in the biological tissue to improve precision of forward solution.

Given the effect of spectrum, SP₇model include four equations at band τ_iand position r:

$\begin{matrix} - \nabla \cdot \frac{1}{3 μ_{a 1} (r, τ_{i})} \nabla ϕ_{1} (r, τ_{i}) + μ_{a} (r, τ_{i}) ϕ_{1} (r, τ_{i}) = S (r, τ_{i}) + (\frac{2}{3} μ_{a} (r, τ_{i})) ϕ_{2} (r, τ_{i}) - (\frac{8}{15} μ_{a} (r, τ_{i})) ϕ_{3} (r, τ_{i}) + (\frac{16}{35} μ_{a} (r, τ_{i})) ϕ_{4} (r, τ_{i}) & (1) \\ - \nabla \cdot \frac{1}{7 μ_{a 3} (r, τ_{i})} \nabla ϕ_{2} (r, τ_{i}) + (\frac{4}{9} μ_{a} (r, τ_{i}) + \frac{5}{9} μ_{a 2} (r, τ_{i})) ϕ_{2} (r, τ_{i}) = - \frac{2}{3} S (r, τ_{i}) + (\frac{2}{3} μ_{a} (r, τ_{i})) ϕ_{1} (r, τ_{i}) + (\frac{16}{45} μ_{a} (r, τ_{i}) + \frac{4}{9} μ_{a 2} (r, τ_{i})) ϕ_{3} (r, τ_{i}) - (\frac{32}{105} μ_{a} (r, τ_{i}) + \frac{8}{21} μ_{a 2} (r, τ_{i})) ϕ_{4} (r, τ_{i}) & (2) \\ - \nabla \cdot \frac{1}{11 μ_{a 5} (r, τ_{i})} \nabla ϕ_{3} (r, τ_{i}) + (\frac{64}{225} μ_{a} (r, τ_{i}) + \frac{16}{45} ϕ_{a 2} (r, τ_{i}) + \frac{9}{25} μ_{a 4} (r, τ_{i})) ϕ_{3} (r, τ_{i}) = \frac{8}{15} S (r, τ_{i}) - (\frac{8}{15} μ_{a} (r, τ_{i})) ϕ_{1} (r, τ_{i}) + (\frac{16}{45} μ_{a} (r, τ_{i}) + \frac{4}{9} μ_{a 2} (r, τ_{i})) ϕ_{2} (r, τ_{i}) + (\frac{128}{525} μ_{a} (r, τ_{i}) + \frac{32}{105} μ_{a 2} (r, τ_{i}) + \frac{54}{175} μ_{a 4} (r, τ_{i})) ϕ_{4} (r, τ_{i}) & (3) \\ - \nabla \cdot \frac{1}{15 μ_{a 7} (r, τ_{i})} \nabla ϕ_{4} (r, τ_{i}) + (\frac{256}{225} μ_{a} (r, τ_{i}) + \frac{64}{245} μ_{a 2} (r, τ_{i}) + \frac{324}{1225} μ_{a 4} (r, τ_{i}) + \frac{13}{49} μ_{a 6} (r, τ_{i})) ϕ_{4} (r, τ_{i}) = - \frac{16}{35} S (r, τ_{i}) + (\frac{16}{35} μ_{a} (r, τ_{i})) ϕ_{1} (r, τ_{i}) - (\frac{32}{105} μ_{a} (r, τ_{i}) + \frac{8}{21} μ_{a 2} (r, τ_{i})) ϕ_{2} (r, τ_{i}) + (\frac{128}{525} μ_{a} (r, τ_{i}) + \frac{32}{105} μ_{a 2} (r, τ_{i}) + \frac{54}{175} μ_{a 4} (r, τ_{i})) ϕ_{3} (r, τ_{i}) & (4) \end{matrix}$

Wherein S represents light function; μ_ais light absorption coefficient; μ_iis light scattering coefficient; μ_ai=μ_a+μ_s(1−gⁱ), gⁱis anisotropic factor; φ_irepresents the linear combination of Legendre moments φ_iof radiosity and can be represented as:

${\begin{matrix} ϕ_{1} = φ_{0} + 2 φ_{2} \\ ϕ_{2} = 3 φ_{2} + 4 φ_{4} \\ ϕ_{3} = 5 φ_{4} + 6 φ_{6} \\ ϕ_{4} = 7 φ_{6} \end{matrix}$

Based on the SP₇equation, set φ₆=φ₄=0, and remain formula (1) and (2), and SP₃equations are obtained:

$\begin{matrix} - \nabla \cdot \frac{1}{3 μ_{a 1} (r, τ_{i})} \nabla ϕ_{1} (r, τ_{i}) +  μ_{a} (r, τ_{i}) ϕ_{1} (r, τ_{i}) - \frac{2}{3} μ_{a} (r, τ_{i}) ϕ_{2} (r, τ_{i}) = S (r, τ_{i}) & (6) \\ - (\frac{2}{3} μ_{a} (r, τ_{i})) ϕ_{1} (r, τ_{i}) - \nabla \cdot \frac{1}{7 μ_{a 3} (r, τ_{i})} \nabla ϕ_{2} (r, τ_{i}) + (\frac{4}{9} μ_{a} (r, τ_{i}) + \frac{5}{9} μ_{a 2} (r, τ_{i})) ϕ_{2} (r, τ_{i}) = - \frac{2}{3} S (r, τ_{i}) & (7) \end{matrix}$

Based on the boundary conditions of SP₇, set φ₆=φ₄=0, and the boundary conditions of SP₃are obtained:

$\begin{matrix} (\frac{1}{2} + A_{1}) ϕ_{1} (r, τ_{i}) + (\frac{1 + B_{1}}{3 μ_{a 1} (r, τ_{i})}) \vec{n} (r, τ_{i}) \cdot \nabla ϕ_{1} (r, τ_{i}) = (\frac{1}{8} + C_{1}) ϕ_{2} (r, τ_{i}) + (\frac{D_{1}}{μ_{a 3} (r, τ_{i})}) \vec{n} (r, τ_{i}) \cdot \nabla ϕ_{2} (r, τ_{i}) & (8) \\ (\frac{7}{24} + A_{2}) ϕ_{2} (r, τ_{i}) + (\frac{1 + B_{2}}{7 μ_{a 3} (r, τ_{i})}) \vec{n} (r, τ_{i}) \cdot \nabla ϕ_{2} (r, τ_{i}) = (\frac{1}{8} + C_{2}) ϕ_{1} (r, τ_{i}) + (\frac{D_{2}}{μ_{a 1} (r, τ_{i})}) \vec{n} (r, τ_{i}) \cdot \nabla ϕ_{1} (r, τ_{i}) & (9) \end{matrix}$

Wherein {right arrow over (n)} represents the outward unit normal vector perpendicular to the boundary; A_iB_iand C_iare a series of constants which are related to the angle distance of boundary reflectivity. The corresponding surface emission light Φ(r,τ_i) is:

$\begin{matrix} Φ (r, τ_{i}) = (\frac{1}{4} + J_{0}) ϕ_{1} (r, τ_{i}) - (\frac{0.5 + J_{1}}{3 μ_{a 1} ϕ_{1} (r, τ_{i})}) \vec{n} \cdot \nabla ϕ_{1} (r, τ_{i}) - (\frac{1}{16} + \frac{2 J_{0} - J_{2}}{3}) ϕ_{2} (r, τ_{i}) - (\frac{J_{3}}{7 μ_{a 1} ϕ_{2} (r, τ_{i})}) \vec{n} \cdot \nabla ϕ_{2} (r, τ_{i}) \end{matrix}$

Wherein J₀, . . . , J₃are a series of constants.

Step two, constructing the linear relationship—based on finite element method to discretize SP₃diffusion equation, and constructing the linear model between boundary measurements and unknown light source for each wavelength;

In the formula (6) and (7), S(r,τ_i) represents the distribution of bioluminescent source and also the unknown quantity to be solved. In order to solve the SP₃equation and its boundary conditions using finite element method, firstly two SP₃equations should be written in corresponding weak solution forms:

$\begin{matrix} \int_{Ω}^{} (- \nabla \cdot \frac{1}{3 μ_{a 1} (r, τ_{i})} \nabla ϕ_{1} (r, τ_{i}) + μ_{a} (r, τ_{i}) ϕ_{1} (r, τ_{i}) - \frac{2}{3} μ_{a} (r, τ_{i}) ϕ_{2} (r, τ_{i})) Ψ (r, τ_{i}) \partial r = \int_{Ω}^{} S (r, τ_{i}) Ψ (r, τ_{i}) \partial r \\ \int_{Ω}^{} (- (\frac{2}{3} μ_{a} (r, τ_{i})) ϕ_{1} (r, τ_{i}) - \nabla \cdot \frac{1}{7 μ_{a 3} (r, τ_{i})} \nabla ϕ_{2} (r, τ_{i}) + (\frac{4}{9} μ_{a} (r, τ_{i}) + \frac{5}{9} μ_{a 2} (r, τ_{i})) ϕ_{2} (r, τ_{i})) Ψ (r, τ_{i}) \partial r = \int_{Ω}^{} (- \frac{2}{3} S (r, τ_{i})) Ψ (r, τ_{i}) \partial r \end{matrix}$

Wherein, Ψ(r,τ_i) is test function, and two boundary conditions for formula (8) and (9) also should be written in corresponding weak solution forms:

$\begin{matrix} \int_{Ω}^{} ((\frac{1}{2} + A_{1}) ϕ_{1} (r, τ_{i}) + (\frac{1 + B_{1}}{3 μ_{a 1} (r, τ_{i})}) \vec{n} (r, τ_{i}) \cdot \nabla ϕ_{1} (r, τ_{i})) Ψ (r, τ_{i}) \partial r = \int_{Ω}^{} ((\frac{1}{8} + C_{1}) ϕ_{2} (r, τ_{i}) + (\frac{D_{1}}{μ_{a 3} (r, τ_{i})}) \vec{n} (r, τ_{i}) \cdot \nabla ϕ_{2} (r, τ_{i})) Ψ (r, τ_{i}) \partial r \\ \int_{Ω}^{} ((\frac{7}{24} + A_{2}) ϕ_{2} (r, τ_{i}) + (\frac{1 + B_{2}}{7 μ_{a 3} (r, τ_{i})}) \vec{n} (r, τ_{i}) \cdot \nabla ϕ_{2} (r, τ_{i})) Ψ (r, τ_{i}) \partial r = \int_{Ω}^{} ((\frac{1}{8} + C_{2}) ϕ_{1} (r, τ_{i}) + (\frac{D_{2}}{μ_{a 1} (r, τ_{i})}) \vec{n} (r, τ_{i}) \cdot \nabla ϕ_{1} (r, τ_{i})) Ψ (r, τ_{i}) \partial r \end{matrix}$

Then the weak solution forms of the boundary conditions are incorporated into that for the SP₃equation using Green's formula.

$\begin{matrix} \int_{Ω}^{} (\frac{1}{3 μ_{a 1}} \nabla ϕ_{1} (r, τ_{i}) + μ_{a} ϕ_{1} (r, τ_{i}) ψ (r, τ_{i}) - \frac{2}{3} μ_{a} ϕ_{2} (r, τ_{i}) ψ (r, τ_{i})) \partial r - \int_{\partial Ω}^{} (\frac{e_{1}}{3 μ_{a 1}} \nabla ϕ_{1} (r, τ_{i}) + \frac{e_{2}}{3 μ_{a 1}} ϕ_{2} (r, τ_{i})) ψ (r, τ_{i}) \partial r = \int_{Ω}^{} S (r, τ_{i}) ψ (r, τ_{i}) \partial r \\ \int_{Ω}^{} (\frac{1}{7 μ_{a 3} (r, τ_{i})} \nabla ϕ_{2} (r, τ_{i}) \cdot \nabla ψ (r, τ_{i}) + (\frac{4}{9} μ_{a} (r, τ_{i}) + \frac{5}{9} μ_{a 2} (r, τ_{i})) ϕ_{2} (r, τ_{i}) ψ (r, τ_{i}) - (\frac{2}{3} μ_{a} (r, τ_{i})) ϕ_{1} (r, τ_{i}) ψ (r, τ_{i}) \partial r = \int_{Ω}^{} (- \frac{2}{3} S (r, τ_{i})) ψ (r, τ_{i}) \partial r \end{matrix}$

Wherein e₁e₂e₃and e₄are respectively:

$\begin{matrix} e_{1} = \frac{1}{W} (\frac{D_{1} (1 + 8 C_{2})}{8 μ_{a 3}} - \frac{(1 + B_{2}) (1 + 2 A_{1})}{14 μ_{a 3}}) \\ e_{2} = \frac{1}{W} (\frac{(1 + B_{2}) (1 + 8 C_{2})}{56 μ_{a 3}} - \frac{D_{1} (7 + 24 A_{2})}{24 μ_{a 3}}) \\ e_{3} = \frac{1}{W} (\frac{(1 + B_{1}) (1 + 8 C_{2})}{24 μ_{a 1}} - \frac{D_{2} (7 + 2 A_{1})}{2 μ_{a 1}}) \\ e_{4} = \frac{1}{W} (\frac{D_{2} (1 + 8 C_{1})}{8 μ_{a 1}} - \frac{(1 + B_{1}) (7 + 24 A_{2})}{72 μ_{a 1}}) \end{matrix}$

In the above formulas, the form for W is as follow:

$W = \frac{(1 + B_{1}) (1 + B_{2})}{21 μ_{a 1} μ_{a 3}} - \frac{D_{1} D_{2}}{μ_{a 1} μ_{a 3}}$

For mesh generation, the distribution of light source S(r,τ_i) is written as S(τ_i) only related to the spectral band, when the corresponding basis function is selected, the whole stiffness matrix M_ijis assembled, the SP₃equations are discretized into the matrix equation as follow:

$[\begin{matrix} M_{11} M_{12} \\ M_{21} M_{22} \end{matrix}] [\begin{matrix} ϕ_{1} \\ ϕ_{2} \end{matrix}] = [\begin{matrix} B & 0 \\ 0 & B \end{matrix}] [\begin{matrix} S (τ_{i}) \\ - \frac{2}{3} S (τ_{i}) \end{matrix}]$

Wherein B is the matrix of size N×N, φ₁and φ₁are two partitioned submatrices and can be represented as:

φ₁=⅓B[3M₁₂⁻¹+2M₂₂⁻¹τM₁₁⁻¹M₁₂+2M₂₁⁻¹M₂₂]⁻¹S(τ_i)

φ₂=⅓B[3M₁₁⁻¹+2M₂₁⁻¹τM₁₁⁻¹M₁₂+2M₂₁⁻¹M₂₂]⁻¹S(τ_i)

The linear relationship between the surface emission light Φ(τ_i) and the distribution of unknown bioluminescent source S(τ_i) should be constructed. In the experiment, only the optical signals on the surface of organisms can be collected, so only the node values on the boundary of φ(τ_i) is remained and the rows which do not contain the nodes on the boundary of matrix are removed, and the equation is obtained as follow:

$Φ (τ_{i}) = β_{1} ϕ_{1} + β_{2} ϕ_{2} = \frac{1}{3} B ({β_{1} [3 M_{12}^{- 1} + 2 M_{22}^{- 1}] [M_{11}^{- 1} M_{12} + 2 M_{21}^{- 1} M_{22}]}^{- 1} + {β_{2} [3 M_{11}^{- 1} + 2 M_{21}^{- 1}] [M_{11}^{- 1} M_{12} + 2 M_{21}^{- 1} M_{22}]}^{- 1} S (τ_{i}) = A (τ_{i}) S (τ_{i})$

The above equation constructs the relationship between the distribution of unknown bioluminescent source and the boundary emission light, i.e.:

Φ(τ_i)=A(τ_i)S(τ_i)

In multispectral conditions, the energy ratio ω(τ_i) of spectral bands is measured by spectrum analysis in advance, set S represent the total energy of all spectral bands of light sources, i.e. S=Σ_i=1^pω(τ_i)S(τ_i), and p represent the number of spectral bands. All spectral bands of light sources and boundary measurements are incorporated, and the relationship between the distribution of unknown bioluminescent source and the boundary emission light in multispectral conditions is as follow:

A_mulS=Φ_mul

Now,

$A_{mul} = [\begin{matrix} ω (τ_{1}) A (τ_{1}) \\ ω (τ_{2}) A (τ_{2}) \\ ⋮ \\ ω (τ_{P}) A (τ_{P}) \end{matrix}], Φ_{mul} = [\begin{matrix} Φ (τ_{1}) \\ Φ (τ_{2}) \\ ⋮ \\ Φ (τ_{p}) \end{matrix}]$

Step three, estimating the shared prior—based on empirical Bayesian maximum likelihood function, inferring parameter α which represents the correlation among multi-wavelength;

In practical problems, the emission light measurement Φ(τ_i) usually contains noise, if the noise obeys Gaussian random distribution with mean zero and variance α₀⁻¹, the maximum likelihood function of Φ(τ_i) related to the distribution of bioluminescent source S(τ_i) and α₀is represented as:

$p (Φ (τ_{i}) | S (τ_{i}), α_{0}) = {(\frac{α_{0}}{2 π})}^{\frac{M}{2}} \exp {- \frac{α_{0}}{2} { Φ (τ_{i}) - A (τ_{i}) S (τ_{i}) }_{2}^{2}}$

In order to describe the correlation among spectral bands, a multitask idea is introduced, and the whole light source is regarded as a whole task T, wherein S(τ_i) is the ith task model in T, and S_j(τ_i) represents the jth element of S(τ_i). If S(τ_i) obeys the Gaussian prior distribution with mean zero and parameter α₀obeys Gamma prior distribution, then:

p(S(τ_i)|α,α₀)=Π_j=1^NN(S_j(τ_i)|0,α₀⁻¹α_j⁻¹)

p(α₀|a,b)=Ga(α₀|a,b)

Wherein α=(α₁, . . . , α_j, . . . α_N)^Tis hyper-parameter, representing prior information. Then, the values of α_jare identical in all spectral bands, i.e. α representing the correlation among various spectral bands, the empirical Bayesian method estimate the hyper-parameter α is used by maximizing the marginal likelihood function, then:

$L = \sum_{i = 1}^{p} \log p (Φ (τ_{i}) | α) = \sum_{i = 1}^{p} \log (\int P (Φ (τ_{i}) | S (τ_{i}), α_{0}) P (S (τ_{i}) | α, α_{0}) P (α_{0} | a, b)) \partial S_{\dot{π}} \partial α_{0} = - \frac{1}{2} \sum_{i = 1}^{p} [(n_{i} + 2 a) \log ({Φ (τ_{i})}^{T} B_{i}^{- 1} Φ (τ_{i}) + b) + \frac{\log \langle B_{i} \rangle}{2}] + const$

Wherein,

B_i=I+A(τ_i)Λ⁻¹A(τ_i)^T, Λ⁻¹=diag(α₁,α₁, . . . ,α_N).

Step four, reconstructing unknown light source—based on the estimated {circumflex over (α)} of the known hyper-parameter and boundary measurements Φ_mul, the bioluminescent source is reconstructed using maximum likelihood function.

After obtaining the estimated {circumflex over (α)}, all spectral bands of light sources and boundary measurements are incorporated, and the maximum likelihood function of the unknown light source S can be represented as:

$P (S | Φ_{mul}, \hat{α}) = \int P (S | Φ_{mul}, \hat{α}, α_{0}) P (α_{0} | a, b) \partial α_{0} = \frac{{Γ (a + \frac{N}{2}) [1 + \frac{1}{2 b} {(S - μ_{mul})}^{T} \sum_{mul}^{- 1} (S - μ_{mul})]}^{- (a + \frac{N}{2})}}{Γ (a) {(2 π b)}^{N / 2} {\langle \sum_{mul} \rangle}^{1 / 2}}$

Wherein,

μ_mul=Σ_mulA_mul^TΦ_mul

Σ_mul=(A_mul^TA_mul+Ξ)⁻¹

In the formula, Ξ=diag({circumflex over (α)}₁, {circumflex over (α)}₁, . . . , {circumflex over (α)}_N), the framework of multitask model is shown in FIG. 1.

The disclosure is based on multitask learning method. Firstly a high order approximation model is used to model the law of light propagation in the biological tissue, then inner-correlation among multi-spectra is researched based on multitask learning method, and finally on this basis three-dimensional reconstruction of the fluorescent source is realized. The reconstruction results show not only bioluminescent sources in a mouse body are accurately reconstructed and located in the disclosure, but computational efficiency has been greatly improved.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 shows a framework of a multitask model;

FIG. 2 shows a collection area of measurements of a digital mouse chest;

FIG. 3 shows reconstruction results of light sources in a lung (a) shows the reconstruction positions of the light sources in the lung, (b)-(d) show coronal slices of reconstructed results using BLT corresponding XCT. In (b)-(d), positions of real light sources marked by circles;

FIG. 4 shows reconstruction results of light sources based on L1 regularization reconstruction method (a)-(c) show coronal slices of reconstructed results using BLT corresponding to XCT, positions of the real light sources marked by circles;

FIG. 5 shows reconstruction results of the light sources in a liver (a) and (b) show reconstructed slices using micro-CT. (c)-(f) show the distribution of light sources reconstructed based on the algorithm proposed in this disclosure, (g)-(j) show the reconstruction results using L1 norm regularization method. (c) and (i) show the coronal slices when y=51 mm, (f) and (j) show the cross section when z=12.97 mm. The positions of the real light sources are marked by circles.

FIG. 6 shows reconstruction results based on two algorithms. (a) and (c) show the reconstruction results using the algorithm proposed in the disclosure, (b) and (d) show the reconstruction results using L1 norm regularization method. (a) and (b) show the coronal section when y=51.05 mm, (c) and (d) show the cross section when z=15.97. Positions of real light sources are marked by circles.

FIG. 7 shows reconstruction results of light sources in a liver (a) shows the coronal slice based on micro-CT data. (b) and (d) show the reconstruction results based on RIRM and L1 RMRM, respectively. (c) and (e) show the reconstruction results based on two algorithms on XCT. Positions of the real light sources are marked by circles.

FIG. 8 shows reconstruction results of light sources (a) shows the reconstruction results based on Bayesian compressed sensing method, (b) shows the reconstruction results based on the proposed algorithm. The positions of the real light sources are marked by circles.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

The disclosure will be described in more detailed below accompanying drawings with preferred embodiments.

In order to test effectiveness of the proposed method in the disclosure, two simulation experiments were performed using heterogeneous digital mouse phantoms. In the experiments, parameters of optical properties of biological tissues were supposed to be known, and the detailed parameters were shown in table 1. The forward boundary measurements were obtained using finite element method in this disclosure. In order to avoid “inverse crime” effectively, forward data was obtained using SP₃equation, a mesh of which contain 42342 nodes and 208966 tetrahedral units, and the mesh used in reconstruction contain 20196 nodes and 108086 tetrahedral units. Given the spectral distribution, two spectral bands adopted in the experiment are 620 nm and 700 nm, respectively. In order to simulate a collection of bioluminescent data, only the measurements of digital mouse chest (a total of 1311 boundary measurement points) were adopted, as shown in FIG. 2.

TABLE 1 Parameters of Optical Properties of Biological Tissues Light Absorption Light Absorption Biological Spectral Coefficient μ_a Coefficient μ_s Tissue Band (mm⁻¹) (mm⁻¹) Muscle 620 nm 0.086 0.429 700 nm 0.0027 0.3792 Liver 620 nm 0.345 0.678 700 nm 0.23 0.648 Lung 620 nm 0.195 2.173 700 nm 0.13 2.124 Bone 620 nm 0.06 2.495 700 nm 0.039 2.34 Heart 620 nm 0.058 0.963 700 nm 0.038 0.905

Experiment 1

First group: in the first group experiment, a spherical light source with the radius of 1.5 mm was placed in a mouse lung, and the central position of which was (14.5, 44.0, 12.8) mm and the distance from the organism surface was 4.8 mm.

The bioluminescent source was reconstructed using the reconstruction algorithm proposed in the disclosure, and the reconstruction results were analyzed and shown in FIG. 3. It can be seen that the bioluminescent source can be reconstructed accurately using the method proposed in the disclosure, the central position of which was (14.69, 44.17, 14.69) mm and the absolute distance from the real light source was 0.28 mm.

Further, in order to verify its effectiveness, the algorithm proposed in the disclosure was compared with the reconstruction algorithm based on L1 regularization method, and the reconstruction results from the latter were shown in FIG. 4. The center position of the light source reconstructed based on L1 regularization method was (14.59, 43.63, 14.59) mm, and the distance from the real light source was 0.39 mm.

In addition, the computation time spent on reconstruction is 1773.9 seconds using L1 regularization method, whereas only 0.67 seconds spent on reconstruction using the algorithm proposed in this disclosure.

Second group: in the second group experiment, the light source with same size was placed in the liver, the central position of which was (22.0, 51.0, 22.0) mm, and the distance from the organism surface was 7.06 mm.

Then, the light source was reconstructed using the algorithm proposed in the disclosure, and the results were shown in FIG. 5. It can be seen that the bioluminescent source can be reconstructed accurately using the algorithm proposed in the disclosure, the central position of which was (22.37, 50.68, 22.37) mm, and the absolute distance from the real light source was 0.5612 mm.

In order to verify the effectiveness of the algorithm, it was compared with the reconstruction algorithm based on L1 regularization method, and the reconstruction results of the latter were shown in FIG. 6. The center position of the light source reconstructed based on L1 regularization reconstruction method was (23.01, 50.06, 14.87) mm, and the distance from the real light source was 2.3239 mm.

Accordingly, it can be seen that the reconstruction based on L1 regularization method had poor effect and big errors when the reconstructed bioluminescent source was located at a deeper position of the biological tissue; whereas the fluorescent source can be located accurately using the method proposed in this disclosure. At the same time, the computation time of reconstruction was 1738.8 seconds using L1 regularization method, whereas only 0.936 seconds using the algorithm in this disclosure.

The summarizing quantitative comparison of the above two group experiments were performed, the reconstruction algorithm based on multitask compression sensing proposed in this disclosure was abbreviated to IRIRM, and the reconstruction algorithm based on L1 norm regularization was abbreviated to L1-RMRM. The comparison results were shown in table 2.

TABLE 2 The Summarizing Quantitative Comparison Between Two Algorithms Real Positions of Reconstruction the Light Reconstructed Position Reconstruction Algorithm Sources Positons Error (mm) Time (s) IRIRM (14.5, 44.0, 12.8) (14.69, 44.17, 12.9) 0.28 0.67 (22.0, 51.0, 13.0) (22.37, 50.68, 13.15) 0.5612 0.936 L1-RMRM (14.5, 44.0, 12.8) (14.59, 43.63, 12.88) 0.39 1773.9 (22.0, 51.0, 13.0) (23.01, 50.06, 14.87) 2.3239 1738.8

Experiment 2

In order to further test the effectiveness of the algorithm, the second experiment of double light sources was performed.

The first group: in the first experimental group, two spherical light sources with the radius of 1.5 mm were respectively placed in the liver and lung, the central position of which were (14, 51, 16) mm and (21, 51, 16) mm, and the distance from the organism surface were 5.56 mm and 5.82 mm, respectively. The reconstruction was respectively performed using the algorithm proposed in this disclosure and L1 regularization method, and the reconstruction results were shown in FIG. 7. The bioluminescent source can be reconstructed accurately using both algorithms, as shown in FIG. 7.

The second group: in the second experimental group, both of two light sources were placed in the liver, the light sources were reconstructed using the above two algorithms and reconstruction results were shown in FIG. 8. As shown in FIG. 8, only one light source was reconstructed based on L1 regularization method; whereas both of two bioluminescent sources can be reconstructed accurately using the algorithm proposed in the disclosure, which further verified the effectiveness of the proposed method.

The summarizing quantitative comparison of two group experiments of double light sources was performed, and the results were shown in table 3.

TABLE 3 The Summarizing Quantitative Comparisons between Two Algorithms Based on Double Light Source Biological Tissues Light Position Error (mm) Sources Reconstruction Light Light Reconstruction Located in Algorithm Source 1 Source 2 Time (s) Lung and Liver IRIRM 0.6772 0.3735 0.9828 L1-RMRM 0.6772 0.3735 1770.6 Liver IRIRM 1.0973 0.3330 0.4836 L1-RMRM NaN 0.3330 1767.7

In the disclosure, the simulation experiments were performed by setting a single light source in a lung and liver of the mouse-shape phantoms, and the algorithm proposed in the disclosure and reconstruction method based on L1 regularization were compared. Further, the experiments of double light sources were performed, the light sources were respectively reconstructed using the algorithm proposed in the disclosure and the reconstruction method based on L1 regularization, and the results were compared. Experimental results suggested the bioluminescent source can be reconstructed and located accurately using the proposed algorithm, and computational effectiveness can be greatly improved, which further testify the effectiveness of the proposed algorithm.

Claims

1. A method for bioluminescence tomography reconstruction using multitask Bayesian compressed sensing, the method comprising: - ∇ · 1 3  μ a   1  ( r, τ i )  ∇ ϕ 1  ( r, τ i ) +  μ a  ( r, τ i )  ϕ 1  ( r, τ i ) = S  ( r, τ i ) + ( 2 3  μ a  ( r, τ i ) )  ϕ 2  ( r, τ i ) - ( 8 15  μ a  ( r, τ i ) )  ϕ 3  ( r, τ i ) + ( 16 35  μ a  ( r, τ i ) )  ϕ 4  ( r, τ i ), ( 1 ) - ∇ · 1 7  μ a   3  ( r, τ i )  ∇ ϕ 2  ( r, τ i ) + ( 4 9  μ a  ( r, τ i ) + 5 9  μ a   2  ( r, τ i ) )  ϕ 2  ( r, τ i ) = - 2 3  S  ( r, τ i ) + ( 2 3  μ a  ( r, τ i ) )  ϕ 1  ( r, τ i ) + ( 16 45  μ a  ( r, τ i ) + 4 9  μ a   2  ( r, τ i ) )  ϕ 3  ( r, τ i ) - ( 32 105  μ a  ( r, τ i ) + 8 21  μ a   2  ( r, τ i ) )  ϕ 4  ( r, τ i ), ( 2 ) - ∇ · 1 11  μ a   5  ( r, τ i )  ∇ ϕ 3  ( r, τ i ) + ( 64 225  μ a  ( r, τ i ) + 16 45  μ a   2  ( r, τ i ) + 9 25  μ a   4  ( r, τ i ) )  ϕ 3  ( r, τ i ) = 8 15  S  ( r, τ i ) - ( 8 15  μ a  ( r, τ i ) )  ϕ 1  ( r, τ i ) + ( 16 45  μ a  ( r, τ i ) + 4 9  μ a   2  ( r, τ i ) )  ϕ 2  ( r, τ i ) + ( 128 525  μ a  ( r, τ i ) + 32 105  μ a   2  ( r, τ i ) + 54 175  μ a   4  ( r, τ i ) )  ϕ 4  ( r, τ i ), ( 3 ) - ∇ · 1 15  μ a   7  ( r, τ i )  ∇ ϕ 4  ( r, τ i ) + ( 256 225  μ a  ( r, τ i ) + 64 245  μ a   2  ( r, τ i ) + 324 1225  μ a   4  ( r, τ i ) + 13 49  μ a   6  ( r, τ i ) )  ϕ 4  ( r, τ i ) = - 16 35  S  ( r, τ i ) + ( 16 35  μ a  ( r, τ i ) )  ϕ 1  ( r, τ i ) - ( 32 105  μ a  ( r, τ i ) + 8 21  μ a   2  ( r, τ i ) )  ϕ 2  ( r, τ i ) + ( 128 525  μ a  ( r, τ i ) + 32 105  μ a   2  ( r, τ i ) + 54 175  μ a   4  ( r, τ i ) )  ϕ 3  ( r, τ i ), ( 4 ) { ϕ 1 = φ 0 + 2   φ 2 ϕ 2 = 3   φ 2 + 4   φ 4 ϕ 3 = 5   φ 4 + 6   φ 6 ϕ 4 = 7   φ 6, - ∇ · 1 3   μ a   1  ( r, τ i )  ∇ ϕ 1  ( r, τ 1 ) + μ a  ( r, τ 1 )  ϕ 1  ( r, τ i ) - 2 3  μ a  ( r, τ i )  ϕ 2  ( r, τ i ) = S  ( r, τ i ) ( 6 ) - ( 2 3  μ a  ( r, τ i ) )  ϕ 1  ( r, τ i ) - ∇ · 1 7   μ a   3  ( r, τ i )  ∇ ϕ 2  ( r, τ i ) + ( 4 9  μ a  ( r, τ i ) + 5 9  μ a   2  ( r, τ i ) )  ϕ 2  ( r, τ i ) = - 2 3  S  ( r, τ i ), ( 7 ) ( 1 2 + A 1 )  ϕ 1  ( r, τ i ) + ( 1 + B 1 3   μ a   1  ( r, τ i ) )  n →  ( r, τ i ) · ∇ ϕ 1  ( r, τ i ) = ( 1 8 + C 1 )  ϕ 2  ( r, τ i ) + ( D 1 μ a   3  ( r, τ i ) )  n →  ( r, τ i ) · ∇ ϕ 2  ( r, τ i ) ( 8 ) ( 7 24 + A 2 )  ϕ 2  ( r, τ i ) + ( 1 + B 2 7   μ a   3  ( r, τ i ) )  n →  ( r, τ i ) · ∇ ϕ 2  ( r, τ i ) = ( 1 8 + C 2 )  ϕ 1  ( r, τ i ) + ( D 2 μ a   1  ( r, τ i ) )  n →  ( r, τ i ) · ∇ ϕ 1  ( r, τ i ), ( 9 ) Φ  ( r, τ i ) = ( 1 4 + J 0 )  ϕ 1  ( r, τ i ) - ( 0.5 + J 1 3   μ a   1  ϕ 1  ( r, τ i ) )  n → · ∇ ϕ 1  ( r, τ i ) - ( 1 16 + 2  J 0 - J 2 3 )  ϕ 2  ( r, τ i ) - ( J 3 7   μ a   1  ϕ 2  ( r, τ i ) )  n → · ∇ ϕ 2  ( r, τ i ), ∫ Ω  ( - ∇ · 1 3   μ a   1  ( r, τ i )  ∇ ϕ 1  ( r, τ i ) + μ a  ( r, τ i )  ϕ 1  ( r, τ i ) - 2 3  μ a  ( r, τ i )  ϕ 2  ( r, τ i ) )  Ψ  ( r, τ i )   r = ∫ Ω  S  ( r, τ i )  Ψ  ( r, τ i )   r,  ∫ Ω  ( - ( 2 3  μ a  ( r, τ i ) )  ϕ 1  ( r, τ i ) - ∇ · 1 7   μ a   3  ( r, τ i )  ∇ ϕ 2  ( r, τ i ) + ( 4 9  μ a  ( r, τ i ) + 5 9  μ a   2  ( r, τ i ) )  ϕ 2  ( r, τ i ) )  Ψ  ( r, τ i )   r = ∫ Ω  ( - 2 3  S  ( r, τ i ) )  Ψ  ( r, τ i )   r, ∫ Ω  ( ( 1 2 + A 1 )  ϕ 1  ( r, τ i ) + ( 1 + B 1 3   μ a   1  ( r, τ i ) )  n →  ( r, τ i ) · ∇ ϕ 1  ( r, τ i ) )  Ψ  ( r, τ i )   r = ∫ Ω  ( ( 1 8 + C 1 )  ϕ 2  ( r, τ i ) + ( D 1 μ a   3  ( r, τ i ) )  n →  ( r, τ i ) · ∇ ϕ 2  ( r, τ i ) )  Ψ  ( r, τ i )   r,  ∫ Ω  ( ( 7 24 + A 2 )  ϕ 2  ( r, τ i ) + ( 1 + B 2 7   μ a   3  ( r, τ i ) )  n →  ( r, τ i ) · ∇ ϕ 2  ( r, τ i ) )  Ψ  ( r, τ i )   r = ∫ Ω  ( ( 1 8 + C 2 )  ϕ 1  ( r, τ i ) + ( D 2 μ a   1  ( r, τ i ) )  n →  ( r, τ i ) · ∇ ϕ 1  ( r, τ i ) )  Ψ  ( r, τ i )   r, ∫ Ω  ( 1 3   μ a   1  ∇ ϕ 1  ( r, τ i ) + μ a  ϕ 1  ( r, τ i )  ψ  ( r, τ i ) - 2 3  μ a  ϕ 2  ( r, τ i )  ψ  ( r, τ i ) )   r - ∫ ∂ Ω  ( e 1 3   μ a   1  ∇ ϕ 1  ( r, τ i ) + e 2 3   μ a   1  ϕ 2  ( r, τ i ) )  ψ  ( r, τ i )   r = ∫ Ω  S  ( r, τ i )  ψ  ( r, τ i )   r,  ∫ Ω  ( 1 7   μ a   3  ( r, τ i )  ∇ ϕ 2  ( r, τ i ) · ∇ ψ  ( r, τ i ) + ( 4 9  μ a  ( r, τ i ) + 5 9  μ a   2  ( r, τ i ) )  ϕ 2  ( r, τ i )  ψ  ( r, τ i ) - ( 2 3  μ a  ( r, τ i ) )  ϕ 1  ( r, τ i )  ψ  ( r, τ i )   r = ∫ Ω  ( - 2 3  S  ( r, τ i ) )  ψ  ( r, τ i )   r, e 1 = 1 W  ( D 1  ( 1 + 8  C 2 ) 8   μ a   3 - ( 1 + B 2 )  ( 1 + 2   A 1 ) 14   μ a   3 ),  e 2 = 1 W  ( ( 1 + B 2 )  ( 1 + 8  C 2 ) 56   μ a   3 - D 1  ( 7 + 24  A 2 ) 24   μ a   3 ),  e 3 = 1 W  ( ( 1 + B 1 )  ( 1 + 8  C 2 ) 24   μ a   1 - D 2  ( 7 + 2  A 1 ) 2   μ a   1 ),  e 4 = 1 W  ( D 2  ( 1 + 8  C 1 ) 8   μ a   1 - ( 1 + B 1 )  ( 7 + 24  A 2 ) 72   μ a   1 ), W = ( 1 + B 1 )  ( 1 + B 2 ) 21   μ a   1  μ a   3 - D 1  D 2 μ a   1  μ a   3, [ M 11 M 12 M 21 M 22 ]  [ ϕ 1 ϕ 2 ] = [ B 0 0 B ]  [ S  ( τ i ) - 2 3  S  ( τ i ) ], Φ  ( τ i ) =  β 1  ϕ 1 + β 2  ϕ 2 =  1 3  B ( β 1  [ 3  M 12 - 1 + 2  M 22 - 1 ]  [ M 11 - 1  M 12 + 2  M 21 - 1  M 22 ] - 1 +  β 2  [ 3  M 11 - 1 + 2  M 21 - 1 ]  [ M 11 - 1  M 12 + 2  M 21 - 1  M 22 ] - 1  S  ( τ i ) =  A  ( τ i )  S  ( τ i ), A mul = [ ω  ( τ 1 )  A  ( τ 1 ) ω  ( τ 2 )  A  ( τ 2 ) ⋮ ω  ( τ P )  A  ( τ P ) ], Φ mul = [ Φ  ( τ 1 ) Φ  ( τ 2 ) ⋮ Φ  ( τ p ) ]; p  ( Φ  ( τ i )  S  ( τ i ), α 0 ) = ( α 0 2   π ) M 2  exp  { - α 0 2   Φ  ( τ i ) - A  ( τ i )  S  ( τ i )  2 2 }, L =  ∑ i = 1 p  log   p  ( Φ  ( τ i )  α ) =  ∑ i = 1 p  log  ( ∫ P  ( Φ  ( τ i )  S  ( τ i ), α 0 )  P  ( S  ( τ i )  α, α 0 )  P  ( α 0  a, b ) )   S τ   i   α 0 =  - 1 2  ∑ i = 1 p  [ ( n i + 2  a )  log  ( Φ  ( τ i ) T  B i - 1  Φ  ( τ i ) + b ) + log   B i  2 ] + const

performing three-dimensional reconstruction of a bioluminescent source based on a high order approximation model associated with light propagation in a biological tissue, inner-correlation among multispectral measurements associated with a multitask learning method and incorporation of a reconstruction algorithm as prior information to reduce ill-posedness of BLT reconstruction by:

(Step one) defining problems and performing initialization by acquiring P wavelengths' boundary data at M measurement points, and setting gamma prior distribution α0 with a shape parameter a and scale parameter b, wherein α0−1 is a variance of emission light measurements Φ(τi) with Gaussian distribution, wherein each of two models is adopted to simulate the law of light propagation in the biological tissue, the two models include radiative transfer equation and diffusion equation, the radiative transfer equation is a mathematical model that describes the law of light propagation accurately, diffusion equation is low-order approximation of the radiative transfer equation, the application of the diffusion equation meets the conditions of biological tissues with strong scattering and low absorption, and the light source is not located on the boundaries of biological tissues, simplified spherical harmonics approximation (SPN) equation is adopted thereof, the high order approximation model describes the light propagation in the biological tissue with strong scattering and low absorption and simulate the law of the light propagation in biological tissues to improve the precision of forward solution,

given effect of spectrum, SP7 model includes four equations at band τi and position r:

wherein S represents light function, μa is a light absorption coefficient; μs is a light scattering coefficient; μai=μa+μs(1−gi), gi is anisotropic factor, φi represents a linear combination of Legendre moments φi of radiosity and is represented as:

based on the SP7 equation, setting φ6=φ4=0 and remain formula (1) and (2), and SP3 equations are obtained:

based on the boundary conditions of SP7, setting φ6=φ4=0, wherein the boundary conditions of SP3 are obtained:

wherein {right arrow over (n)} represents the outward unit normal vector which is perpendicular to the boundary, Ai Bi and Ci are a series of constants, related to the angle distance of boundary reflectivity; the corresponding surface emission light Φ(r, τi) is:

wherein, J0,..., J3 are a series of constants;

(Step two) constructing the linear relationship based on finite element method to discretize SP3 diffusion equation, and constructing the linear model between boundary measurements and unknown light source for each wavelength, wherein, in the formula (6) and (7), S(r,τi) represents the distribution of bioluminescent light source and also the unknown quantity to be solved,

incorporating two SP3 equations should be written in corresponding weak solution forms to solve the SP3 equation and related boundary conditions using finite element method:

wherein ψ(r, τi) is test function,

Incorporating two boundary conditions for the formula (8) and (9) in corresponding weak solution forms:

incorporating the weak solution forms of the boundary conditions into that for the SP3 equation using Green's formula;

wherein e1 e2 e3 and e4 are respectively:

wherein in the above formulas, the form for W is as follow:

wherein, for mesh generation, the distribution of light source S(r, τi) is written as S(τi) related to the spectral band, when the corresponding basis function is selected, the whole stiffness matrix Mij is assembled,

discretizing the SP3 equation into the matrix equation as follow:

wherein B is the matrix of size N×N, φ1 and φ1 are two partitioned submatrices and are represented as: φ1=⅓B[3M12−1+2M22−1τM11−1M12+2M21−1M22]−1S(τi), φ2=⅓B[3M11−1+2M21−1τM11−1M12+2M21−1M22]−1S(τi),

wherein the linear relationship between the surface emission light Φ(τi) and the distribution of unknown bioluminescent light source S(τi) are constructed, in the experiment, only the optical signals on the surface of organisms are collected such that only the node values on the boundary of Φ(τi) is remained and the rows which do not contain the nodes on the boundary of matrix are removed, and the equation is obtained as follow:

wherein the above equation constructs the relationship between the distribution of unknown bioluminescent light source and the boundary emission light, which is: Φ(τi)=A(τi)S(τi)

wherein, in multispectral conditions, the energy ratio ω(τi) of spectral bands is measured by spectrum analysis in advance,

setting S to represent the total energy of all spectral bands of light sources, which is S=Σi=1pω(τi)S(τi), and p to represent the number of spectral bands,

wherein all spectral bands of light sources and boundary measurements are incorporated, the relationship between the distribution of unknown bioluminescent light source and the boundary emission light in multispectral conditions is formulated as follow: AmulS=Φmul Now,

(Step three) estimating the shared prior based on empirical Bayesian maximum likelihood function, infer parameter α which represents the correlation among multi-wavelength;

wherein, the emission light measurement Φ(τi) contains noise, if the noise obeys Gaussian random distribution with mean zero and variance α0−1, the maximum likelihood function of Φ(τi) related to the distribution of bioluminescent light source S(τi) and α0 is represented as:

introducing a multitask idea is introduced to describe the correlation among spectral bands, wherein the whole light source is regarded as a whole task T, S(τi) is the ith task model of T, and Sj(τi) represents the jth element of S(τi), and if S(τi) obeys the Gaussian prior distribution with mean zero and the parameter α0 obeys Gamma prior distribution, then: p(S(τi)|α,α0)=Πj=1NN(Sj(τi)|0,α0−1αj−1) p(α0|a,b)=Ga(α0|a,b),

wherein α=(α1,... αj,... αN)T is hyper-parameter, representing prior information, the values of αj are identical in all spectral bands, i.e. α representing the correlation among various spectral bands, an empirical Bayesian method that estimates the hyper-parameter α is used by maximizing a marginal likelihood function, then:

wherein Bi=I+A(τi)Λ−1 A(τi)T, Λ−1=diag(α1, α1,..., αN); and

(Step four) reconstructing an unknown light source based on the estimated {circumflex over (α)} of the known hyper-parameter and boundary measurements Φmul, wherein the bioluminescent light source is reconstructed using maximum likelihood function,

after obtaining the estimated {circumflex over (α)}, incorporating all spectral bands of light sources and boundary measurements, wherein the maximum likelihood function of the unknown light source S is represented as: P  ( S  Φ mul, α ^ ) =  ∫ P  ( S  Φ mul, α ^, α 0 )  P  ( α 0  a, b )   α 0 =  Γ  ( a + N 2 )  [ 1 + 1 2  b  ( S - μ mul ) T  ∑ mul - 1  ( S - μ mul ) ] - ( a + N 2 ) Γ  ( a )  ( 2   π   b ) N / 2   ∑ mul  1 / 2, wherein, μmul=ΣmulAmulTΦmul Σmul=(AmulTAmul+Ξ)−1, and

wherein, in the formula, Ξ=diag({circumflex over (α)}1, {circumflex over (α)}1..., {circumflex over (α)}N); a framework of multitask model is shown below: