Patent Document ID: 9672639
Application ID: 15108394
Patent Flag: 1

Claim One:
1. A method for bioluminescence tomography reconstruction using multitask Bayesian compressed sensing, the method comprising: 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 (SP N ) 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, SP 7 model includes four equations at band τ i and position r: - ∇ · 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 ) wherein S represents light function, μ a is a light absorption coefficient; μ s is a light scattering coefficient; μ ai =μ a +μ s (1−g i ), g i is anisotropic factor, φ i represents a linear combination of Legendre moments φ i of radiosity and is represented as: { φ 1 = ϕ 0 + 2 ⁢ ϕ 2 φ 2 = 3 ⁢ ϕ 2 + 4 ⁢ ϕ 4 φ 3 = 5 ⁢ φ 4 + 6 ⁢ ϕ 6 φ 4 = 7 ⁢ ϕ 6 , ⁢ , based on the SP 7 equation, setting φ 6 =φ 4 =0 and remain formula (1) and (2), and SP 3 equations are obtained: - ∇ · 1 3 ⁢ μ a ⁢ ⁢ 1 ⁡ ( r , τ i ) ⁢ ∇ φ 1 ⁡ ( r , τ i ) + μ a ⁡ ( r , τ i ) ⁢ φ 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 ) based on the boundary conditions of SP 7 , setting φ 6 =φ 4 =0, wherein the boundary conditions of SP 3 are obtained: ( 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 ) wherein {right arrow over (n)} represents the outward unit normal vector which is perpendicular to the boundary, A i B i and C i are a series of constants, related to the angle distance of boundary reflectivity; the corresponding surface emission light Φ(r, τ i ) is: Φ ⁡ ( 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 ) , wherein, J 0 ,. .. , J 3 are a series of constants; (Step two) constructing the linear relationship based on finite element method to discretize SP 3 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 SP 3 equation and related boundary conditions using finite element method: ∫ Ω ⁢ ( - ∇ · 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 ) ⁢ dr = ∫ Ω ⁢ S ⁡ ( r , τ i ) ⁢ Ψ ⁡ ( r , τ i ) ⁢ dr , ⁢ ∫ Ω ⁢ ( - ( 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 ) ⁢ dr = ∫ Ω ⁢ ( - 2 3 ⁢ S ⁡ ( r , τ i ) ) ⁢ Ψ ⁡ ( r , τ i ) ⁢ ⁢ dr wherein Ψ(r, τ i ) is test function, Incorporating two boundary conditions for the formula (8) and (9) in corresponding weak solution forms: ∫ Ω ⁢ ( ( 1 2 + A 1 ) ⁢ φ 1 ⁡ ( r , τ i ) + ( 1 + B 1 3 ⁢ μ a ⁢ ⁢ 1 ⁡ ( r , τ i ) ) ⁢ n → ⁡ ( r , τ i ) · ∇ φ 1 ⁡ ( r , τ i ) ) ⁢ Ψ ⁡ ( r , τ i ) ⁢ dr = ∫ Ω ⁢ ( ( 1 8 + C 1 ) ⁢ φ 2 ⁡ ( r , τ i ) + ( D 1 μ a ⁢ ⁢ 3 ⁡ ( r , τ i ) ) ⁢ n → ⁡ ( r , τ i ) · ∇ φ 2 ⁡ ( r , τ i ) ) ⁢ Ψ ⁡ ( r , τ i ) ⁢ ⁢ dr , ⁢ ∫ Ω ⁢ ( ( 7 24 + A 2 ) ⁢ φ 2 ⁡ ( r , τ i ) + ( 1 + B 2 7 ⁢ μ a ⁢ ⁢ 3 ⁡ ( r , τ i ) ) ⁢ n → ⁡ ( r , τ i ) · ∇ φ 2 ⁡ ( r , τ i ) ) ⁢ Ψ ⁡ ( r , τ i ) ⁢ dr = ∫ Ω ⁢ ( ( 1 8 + C 2 ) ⁢ φ 2 ⁡ ( r , τ i ) + ( D 2 μ a ⁢ ⁢ 1 ⁡ ( r , τ i ) ) ⁢ n → ⁡ ( r , τ i ) · ∇ φ 1 ⁡ ( r , τ i ) ) ⁢ Ψ ⁡ ( r , τ i ) ⁢ dr , incorporating the weak solution forms of the boundary conditions into that for the SP 3 equation using Green's formula; ∫ Ω ⁢ ( 1 3 ⁢ μ a ⁢ ⁢ 1 ⁢ ∇ φ 1 ⁡ ( r , τ i ) + μ a ⁢ φ 1 ⁡ ( r , τ i ) ⁢ ψ ⁢ ( r , τ i ) - 2 3 ⁢ μ a ⁢ φ 2 ⁡ ( r , τ i ) ⁢ ψ ⁡ ( r , τ i ) ) ⁢ dr - ∫ ∂ Ω ⁢ ( e 1 3 ⁢ μ a ⁢ ⁢ 1 ⁢ ∇ φ 1 ⁢ ( r , τ i ) + e 2 3 ⁢ μ a ⁢ ⁢ 1 ⁢ φ 2 ⁡ ( r , τ i ) ) ⁢ ψ ⁡ ( r , τ i ) ⁢ ⁢ dr = ∫ Ω ⁢ S ⁡ ( r , τ i ) ⁢ ψ ⁡ ( r , τ i ) ⁢ dr , ⁢ ∫ Ω ⁢ ( 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 ) ⁢ dr = ∫ Ω ⁢ ( - 2 3 ⁢ S ⁡ ( r , τ i ) ) ⁢ ψ ⁡ ( r , τ i ) ⁢ dr , wherein e 1 e 2 e 3 and e 4 are respectively: 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 ) , wherein in the above formulas, the form for W is as follow: W = ( 1 + B 1 ) ⁢ ( 1 + B 2 ) 21 ⁢ μ a ⁢ ⁢ 1 ⁢ μ a ⁢ ⁢ 3 - D 1 ⁢ D 2 μ a ⁢ ⁢ 1 ⁢ μ a ⁢ ⁢ 3 , 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 M ij is assembled, discretizing the SP 3 equation into the matrix equation as follow: [ M 11 ⁢ M 12 M 21 ⁢ M 22 ] ⁡ [ φ 1 φ 2 ] = [ B 0 0 B ] ⁡ [ S ⁡ ( τ i ) - 2 3 ⁢ S ⁡ ( τ i ) ] , wherein B is the matrix of size N×N, φ 1 and φ 1 are two partitioned submatrices and are represented as: 
 φ 1 =⅓ B [3 M 12 −1 +2 M 22 −1 IM 11 −1 M 12 +2 M 21 −1 M 22 ] −1 S (τ i ), 
 φ 2 =⅓ B [3 M 11 −1 +2 M 21 −1 IM 11 −1 M 12 +2 M 21 −1 M 22 ] −1 S (τ 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: 
 Φ(τ i )=β 1 φ 1 +β 2 φ 2 =⅓ B (β 1 [3 M 12 −1 +2 M 22 −1 IM 11 −1 M 12 +2 M 21 −1 M 22 ] −1 +β 2 [3 M 11 −1 +2 M 21 −1 IM 11 −1 M 12 +2 M 21 −1 M 22 ] −1 S (τ i )= A (τ i ) S (τ i ), 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=1 p ω(τ 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: 
 A mul S=Φ mul Now, A mul = [ ω ⁡ ( τ 1 ) ⁢ A ⁡ ( τ 1 ) ω ⁡ ( τ 2 ) ⁢ A ⁡ ( τ 2 ) ⋮ ω ⁡ ( τ P ) ⁢ A ⁡ ( τ P ) ] , Φ mul = [ Φ ⁡ ( τ 1 ) Φ ⁡ ( τ 2 ) ⋮ Φ ⁡ ( τ P ) ] ; (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: p ⁡ ( Φ ⁡ ( τ i ) | S ⁡ ( τ i ) , α 0 ) = ( α 0 2 ⁢ π ) M 2 ⁢ exp ⁢ { - α 0 2 ⁢  Φ ⁡ ( τ i ) - A ⁡ ( τ i ) ⁢ S ⁡ ( τ i )  2 2 } , 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 S j (τ 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=1 N N ( S j (τ i )|0,α 0 −1 α j −1 ), 
 p (α 0 |a,b )= Ga (α 0 |a,b ), wherein α=(α 1 ,. .. α j ,. .. α N ) T NY 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: L = ∑ i = 1 p ⁢ log ⁢ ⁢ p ⁡ ( Φ ⁡ ( τ i ) | α ) = ∑ i = 1 p ⁢ log ⁡ ( ∫ P ⁡ ( Φ ⁡ ( τ i ) ⁢  S ⁡ ( τ i ) , α 0 ) ⁢ P ( S ⁡ ( τ i )  ⁢ α , α 0 | a , b ) ) ⁢ dS τ ⁢ ⁢ i ⁢ d ⁢ ⁢ α 0 = - 1 2 ⁢ ∑ i = 1 p ⁢ [ ( n i + 2 ⁢ a ) ⁢ log ⁡ ( Φ ⁡ ( τ i ) T ⁢ B i - 1 ⁢ Φ ⁡ ( τ i ) + b ) + log ⁢  B i  2 ] + const wherein B i =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 ) ⁢ d ⁢ ⁢ α 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 =Σ mul A mul T Φ mul 
 Σmul=( A mul T A mul +Ξ) −1 , and wherein, in the formula, Ξ=diag({circumflex over (α)} 1 , {circumflex over (α)} 1 ,. .. , {circumflex over (α)} N ).