CN102940482B - Adaptive tomographic fluorescence imaging (TFI) reconstructing method - Google Patents

Abstract

An adaptive TFI reconstructing method includes converting a diffusion equation into a linear equation through a finite element method; constructing a linear relation between unknown fluorescent light source distribution and a surface fluorescence measurement value; calculating a current regularization parameter and selecting an element with the largest absolute value in a difference correlation coefficient into a support set I; extracting all the current elements in the updated support set from corresponding lines in a matrix A, constituting a matrix, and obtaining a next-step searching direction; calculating the length of the next step and updating the support set; obtaining a next-step result according to the solved searching direction and step length iteration, and updating the regularization parameter; and determining whether stopping requirements are satisfied, if yes, finishing the reconstructing process, or otherwise, turning to the fourth step. According to the method, the regularization parameter is not needed to predict but determined adaptively during the reconstructing method, the reconstructing robustness is improved through adaptive regularization parameter selection, and the reconstructing efficiency is improved greatly.

Description

Adaptive fluorescence tomography rebuilding method

Technical field

The present invention relates to optical molecular video imaging modal excitation fluorescence fault imaging (TomographicFluorescence Imaging, TFI) technology, be related specifically to a kind of adaptive fluorescence tomography rebuilding method.

Background technology

As a kind of optical molecular video imaging mode, fluorescence excitation imaging technique has obtained development rapidly and has applied widely.By the biological tissue's mark optical molecular probe to area-of-interest, and gather the fluorescence signal of body surface, we can obtain the information of organism molecular cell level with a kind of noninvasive method.But due in visible ray and near infrared range, can produce serious scattering phenomenon when photon is propagated in biological tissue, therefore traditional two-dimentional fluorescence excitation imaging technique cannot show the accurate location of fluorescence light source.And exciting tomography fluorescence imaging (TFI) is a kind of three-dimensional reconstruction, by collection surface fluorescence data and based on specific inverse problem model, the three-dimensional that can realize fluorescence light source is accurately located.

TFI is a kind of typical ill-posed problem, this is because measurement data be only surface fluorescence distribution, and need solve be whole imaging space fluorescence light source distribution, in this case, the position reconstruction of light source does not have unique solution, and very responsive to noise.In order to obtain accurate and stable reconstructed results, usual way comprises regularization term in optimization problem, and this regularization term can be regarded as the priori of distribution of light sources.Modal regularization is Tikhonov regularization method, by adding L2 norm constraint in optimization problem, TFI problem can be made more stable.Tradition TFI method for reconstructing pre-estimates the value of regularization parameter based on experience value, and regularization parameter is remained unchanged in whole iterative reconstruction process.But the regularization parameter required by different Problems of Reconstructions is usually different, traditional reconstructing method is difficult to the suitable regularization parameter of Accurate Prediction thus completes Exact Reconstruction.Therefore suitable regularization parameter is selected to be still a challenge.

Summary of the invention

The present invention is directed to the regularization On The Choice in TFI reconstruction technique, propose a kind of fast robust TFI method for reconstructing of adaptive regularization parameter.

For achieving the above object, a kind of adaptive fluorescence tomography rebuilding method, comprising:

S1 utilizes Finite Element Method that diffusion equation is converted into linear equation;

S2 sets up the linear relationship between unknown fluorescence light source distribution and surface fluorescence measured value;

S3 calculates current regularization parameter, and the element of maximum absolute value in residual error coefficient correlation is selected into support set I;

Row corresponding in matrix A for current all elements in support set I after renewal take out by S4, form a matrix, obtain next step the direction of search;

S5 calculates next step step-length and upgrades support set;

S6 obtains next step result according to the direction of search solved and step iteration, and upgrades regularization parameter;

S7 judges whether to reach stop condition, if reach, then process of reconstruction terminates, otherwise forwards step S4 to.

The present invention without look-ahead regularization parameter, but determines regularization parameter adaptively in process of reconstruction.The present invention does not improve the robustness of reconstruction technique by means of only adaptive regularization parameter Selection Strategy, and greatly improves reconstruction efficiency.

Accompanying drawing explanation

Fig. 1 is the flow chart of the inventive method;

Fig. 2 is that multiple light courcess is imitated the tomograph of body and is positioned at z=0 plan cross-sectional view;

Fig. 3 is the single light source reconstructed results in 4 groups of measurement data situations;

Fig. 4 is the two light source reconstruction results in 4 groups of measurement data situations;

Fig. 5 is three light source reconstruction results in 4 groups of measurement data situations.

Detailed description of the invention

Describe method for reconstructing of the present invention in detail below in conjunction with accompanying drawing, be to be noted that described embodiment is only intended to be convenient to the understanding of the present invention, and any restriction effect is not play to it.

Fig. 1 is the overall procedure of method for reconstructing of the present invention:

Step 101: the present invention adopts diffusion equation as the imaging model of TFI, model comprises excitation process and emission process two diffusion equations be coupled, by using Finite Element Method, can by diffusion equation discretization, so that carry out follow-up process;

Step 102: for TFI inverse problem, surface fluorescence distribution is known, and Internal Fluorescent distribution of light sources is unknown, by carrying out matrixing to the diffusion model after discretization, the linear equation between surface measurement data and unknown distribution of light sources can be set up.Simultaneously for the situation of many group measurement data, multiple linear equation can be combined as unified equation;

Step 103: calculate current regularization parameter, and the element of maximum absolute value in residual error coefficient correlation is selected into support set;

Step 104: the direction of search being obtained next step by solving equation;

Step 105: calculate next step step-length and upgrade support set;

Step 106: obtain next step result according to the direction of search solved and step iteration, and upgrade regularization parameter;

Step 107: judge whether to reach stop condition.If reach, then process of reconstruction terminates; Otherwise forward step 104 to.

Describe in detail one by one the committed step involved by TFI method for reconstructing of the present invention below, concrete form is as described below:

Step 101: utilize Finite Element Method that diffusion equation is converted into linear equation;

The Mathematical Modeling that accurate description light transmits in nonuniformity biological tissue is radiation transfer equation, and it is a complicated differential-intergral equation, solves very difficult.Fortunately, at visible ray and near infrared light spectral coverage, there is when photon transmits in biological tissues the feature of high scattering, low absorption, in this case, diffusion equation can Approximate radiative transmission equation well, and the computation complexity of diffusion equation is lower, can solve in finite time, so be suitable as very much TFI imaging model.When utilizing continuous wave confined laser as excitaton source, the following imaging process of two diffusion equations be coupled to TFI can be used to be described:

$\left\{\begin{array}{c}-\▿\·\left[D_x\left(r\right){\▿\mathrm{\Φ}}_x\left(r\right)\right]+{\mathrm{\μ}}_{\mathrm{ax}}\left(r\right){\mathrm{\Φ}}_x\left(r\right)=\mathrm{\Θ\δ}(r-r_l)\\ -\▿\·\left[D_m\left(r\right){\▿\mathrm{\Φ}}_m\left(r\right)\right]+{\mathrm{\μ}}_{\mathrm{am}}\left(r\right){\mathrm{\Φ}}_m\left(r\right)={\mathrm{\Φ}}_x\left(r\right){\mathrm{\η\μ}}_{\mathrm{af}}\left(f\right)\end{array}(r\∈\mathrm{\Ω})\right.$

The boundary condition of above-mentioned diffusion equation can be described as:

${\mathrm{\Φ}}_{x,m}\left(r\right)+2q{D}_{x,m}\left(r\right)[\stackrel{\→}{v}\left(r\right)\·{\mathrm{\Φ}}_{x,m}\left(r\right)]=0,(r\∈\∂\mathrm{\Ω})$

Wherein, Ω represents whole imaging space, representation space border, subscript x and m represents exciting light and utilizing emitted light respectively, μ _{ax, am}absorption coefficient, D _{x, m}=1/3 (μ _{ax, am}+ μ ' _{sx, sm}) be diffusion coefficient, μ ' _{sx, sm}reduced scattering coefficient, Φ _{x, m}represent photon density, Θ δ (r-r _l) representing isotropic point-like excitation source, Θ represents the intensity of light source, be the export-oriented unit normal vector on border, q is the particular constant depending on border optical reflection factor deviation, η μ _afrepresent fluorescence quantum yield to be reconstructed.

Under finite element theory framework, first above-mentioned diffusion equation and boundary condition thereof are represented as following weak form:

$\left\{\begin{array}{c}{\∫}_{\mathrm{\Ω}}(D_x\left(r\right){\▿\mathrm{\Φ}}_x\left(r\right)\·\▿\mathrm{\Ψ}\left(r\right)+{\mathrm{\μ}}_{\mathrm{ax}}\left(r\right){\mathrm{\Φ}}_x\left(r\right)\mathrm{\Ψ}\left(r\right))\mathrm{dr}\\ +{\∫}_{\∂\mathrm{\Ω}}\frac{1}{2A(r;n,n^{\′})}{\mathrm{\Φ}}_x\left(r\right)\mathrm{\Ψ}\left(r\right)\mathrm{dr}={\∫}_{\mathrm{\Ω}}\mathrm{\Θ\δ}(r-r_l)\mathrm{\Ψ}\left(r\right)\mathrm{dr}\\ {\∫}_{\mathrm{\Ω}}(D_m\left(r\right){\▿\mathrm{\Φ}}_m\left(r\right)\·\▿\mathrm{\Ψ}\left(r\right)+{\mathrm{\μ}}_{\mathrm{am}}\left(r\right){\mathrm{\Φ}}_m\left(r\right)\mathrm{\Ψ}\left(r\right))\mathrm{dr}\\ +{\∫}_{\∂\mathrm{\Ω}}\frac{1}{2A(r;,n,n^{\′})}{\mathrm{\Φ}}_m\left(r\right)\mathrm{\Ψ}\left(r\right)\mathrm{dr}={\∫}_{\mathrm{\Ω}}{\mathrm{\Φ}}_x\left(r\right){\mathrm{\η\μ}}_{\mathrm{af}}\left(r\right)\mathrm{\Ψ}\left(r\right)\mathrm{dr}\end{array}\right.$

Wherein, Ψ (r) represents any test function.Then, utilize tetrahedron to carry out subdivision to whole imaging space, and using corresponding basic function as test function Ψ (r), thus above-mentioned weak form can be turned to two following equatioies by discrete:

K _xΦ _x＝S _x

K _mΦ _m＝FX

Wherein, K _{x, m}be sytem matrix, matrix F is that vectorial X represents the light-source quantum yield that needs are rebuild by obtaining the fluorescence light source quantum yield discretization of the unknown.

Step 102: set up the linear relationship between unknown fluorescence light source distribution and surface fluorescence measured value;

For excitation process, photon density Φ _xcan by solving K _xΦ _x=S _xand directly obtain, the Φ obtained _xusing the energy source as emission process.For TFI inverse problem, due to K _mfor symmetric positive definite matrix, therefore following matrix equation can be obtained:

${\mathrm{\Φ}}_{m,l}=K_{m,l}^{-1}\mathrm{FX}=B_{l}X$

Remove vectorial Φ _{m, l}in non-surface measurement element and matrix B _lthe row of middle correspondence, above-mentioned equation can be further converted to:

${\mathrm{\Φ}}_{m,l}^{\mathrm{meas}}=A_{l}X$

For TFI, usually use multiple point-like excitation source to obtain many group fluorescence measurement data, each group data can obtain above-mentioned matrix equation, multiple equation group is loaded and can obtain following equation:

Φ＝AX

Wherein,

$\mathrm{\Φ}=\left\{\begin{array}{c}{\mathrm{\Φ}}_{m,1}^{\mathrm{meas}}\\ {\mathrm{\Φ}}_{m,2}^{\mathrm{meas}}\\ \·\\ \·\\ \·\\ {\mathrm{\Φ}}_{m,L}^{\mathrm{meas}}\end{array}\right\},$ $A=\left[\begin{array}{c}{A}_1\\ A_2\\ \·\\ \·\\ \·\\ A_L\end{array}\right]$

Because TFI is an ill posed inverse problem, therefore need to incorporate certain regularization method to make TFI problem more stable, what the inventive method adopted is L1 norm sparsity constraints, and under this constraints, TFI is converted into the optimization problem with following form:

$\underset{X}{\min }{F}_{\mathrm{\λ}}\left(X\right)=\frac{1}{2}{\left|\right|\mathrm{AX}-\mathrm{\Φ}\left|\right|}_2^2+\mathrm{\λ}{\left|\right|X\left|\right|}_1$

λ is regularization parameter, F _λ(X) be object function to be optimized.

Step 103: for realizing minimizing object function, first to F _λ(X) differentiate obtains:

${\∂}_{X_{\mathrm{\λ}}}{F}_{\mathrm{\λ}}\left(X_{\mathrm{\λ}}\right)=A^T({\mathrm{AX}}_{\mathrm{\λ}}-\mathrm{\Φ})+{\mathrm{\λ}\∂}_{X_{\mathrm{\λ}}}{\left|\right|X_{\mathrm{\λ}}\left|\right|}_1$

Wherein be || X _λ|| ₁partial differential, and to be determined by following formula:

${\∂}_{X_{\mathrm{\λ}}}{\left|\right|X_{\mathrm{\λ}}\left|\right|}_1=\{u\∈R^n|\left.\begin{array}{cc}{u}_i=\mathrm{sgn}\left(X_{\mathrm{\λ},i}\right),& X_{\mathrm{\λ},i}\≠0\\ u_i\∈[-\mathrm{1,1}],& X_{\mathrm{\λ},i}\≠0\end{array}\right\}$

Set I={i:X _λi () ≠ 0} is used for representing unknown vector X _λsupport set.The relative coefficient c of current residue is determined by following formula:

c＝A ^T(Φ-AX _λ)

After obtaining the coefficient correlation c of residual error, then determine initial regularization parameter λ:

I＝{j：|c _k(j)|＝||c _k|| _∞＝λ}

Wherein j is the element of maximum absolute value in vectorial c, is added by jth element and supports in set I, and makes regularization parameter λ assignment be the value of the element of maximum absolute value in c.

Step 104: row corresponding in matrix A for current all elements in the support set I after renewal are taken out, forms a matrix A _i.Next, by separating following equation, next step direction of search p is obtained _k:

$A_I^T{A}_I{p}_k\left(I\right)=\mathrm{sgn}\left(c_k\left(I\right)\right)$

At acquisition p _kvalue after, by p _kbe not set to 0 in the value supporting the component elements in set I in vector, thus obtain p _k(I).

Step 105: the step-length γ calculating next step direction of search _k, and upgrade support set.

Step 106: obtain next step direction of search p in step 104 _k, and step 105 obtains the step-length γ in next step direction of search _kafter, current unknown number X can be upgraded _k:

X _k＝X _k-1+γ _kp _k

After obtaining new iteration result, next will automatically upgrade regularization parameter λ:

λ＝λ-γ _k

Step 107: this step judges whether the condition of iteration stopping reaches.Calculate the value F of object function after upgrading _λ(X _k), and before upgrading, target function value is F _λ(X _k-1).If | _fλ (X _k)-F _λ(X _k-1) |/F _λ(X _k) < T, then stop iteration, otherwise forward step 104 to and start next round iteration.Wherein the value of threshold value T is 0.0001 ~ 0.01.

Wherein, step 105 comprises following 5 sub-steps:

Step 201: material calculation

${\mathrm{\&gamma;}}_k^{+}=\underset{i\&Element;I^c}{\min }\left\{\frac{\mathrm{\&lambda;}-c_k\left(i\right)}{1-a_i^T{A}_I{p}_k\left(I\right)}\right\}$

Wherein I ^cit is the supplementary set of I.A _ithe i-th row of matrix A.Obtain while record minimum of a value occur position pos+.

Step 202: material calculation

${\mathrm{\&gamma;}}_k^{-}=\underset{i\&Element;I}{\min }\{-X_k\left(i\right)/p_k\left(i\right)\}$

Obtain while record minimum of a value occur position pos-.

Step 203: compare with size.If forward step 204 to; Otherwise forward step S205 to;

Step 204: the position pos+ according to record finds corresponding element, by it from supplementary set I ^cin move on to and support in set I, thus the number supporting set element increases, and makes

Step 205: the position pos-according to record finds corresponding element, removes to supplementary set I by it from support set I ^cin, thus the number supporting set element reduces, and make

Operation result

In order to verify the accuracy of the inventive method, we utilize the imitative body of three-dimensional nonuniformity emulation to carry out TFI and rebuild experiment.

In this experiment, it is cylindrical that the nonuniformity adopted imitates body, wherein contain 4 regions, represent muscle (M), lung (L), heart (H) and bone (B) respectively, table 1 corresponding to different tissues gives at optical coefficient.The diameter of this imitative body and be highly 20mm, Fig. 2 shows the tomograph of this imitative body and is positioned at the sectional view of z=0 plane.In this experiment, we have carried out single light source, two light source, three light source reconstruction experiments respectively, and as shown in Figure 2, all light sources are spherical, and diameter is 2mm, and be centrally located in z=0 plane, the quantum yield of light source is set as 0.5.Black round dot in the sectional view of Fig. 2 represents the position of point-like excitation source, and for each spot light, surface fluorescence distribution is that the periphery in 160 ° of visuals field is over there measured.In order to rebuild light source, this nonuniformity is imitated body and is turned to 3430 points and 17623 tetrahedrons by discrete.

We used 4 point-like excitaton sources and generate fluorescence measurement data, Fig. 3, Fig. 4 and Fig. 5 respectively illustrate the reconstructed results of single light source, two light source and three light sources.Reconstruction time is respectively 0.39 second, 0.42 second and 0.48 second.Result shows, this method for reconstructing can without look-ahead regularization parameter, but in process of reconstruction the value of adaptive decision regularization parameter, result very accurately can be obtained, and have and rebuild speed faster.

Optical coefficient (the unit: mm of zones of different in body imitated by table 1 ^-1)

The above; be only the detailed description of the invention in the present invention; but protection scope of the present invention is not limited thereto; any people being familiar with this technology is in the technical scope disclosed by the present invention; the conversion or replacement expected can be understood; all should be encompassed in and of the present inventionly comprise within scope, therefore, protection scope of the present invention should be as the criterion with the protection domain of claims.

Claims (3)

1. an adaptive fluorescence tomography rebuilding method, comprising:

S1 utilizes Finite Element Method that diffusion equation is converted into linear equation;

S2 sets up the linear relationship between fluorescence light source distribution to be solved and surface fluorescence measured value, obtains the matrix A represented the linear relationship between fluorescence light source distribution and surface fluorescence measured value, and obtains objective function F to be optimized based on described matrix A _λ(X);

S3 calculates described objective function F _λ(X) current regularization parameter λ _k, and by residual error coefficient correlation c _kthe element of middle maximum absolute value is selected into and supports set I;

Row corresponding in matrix A for current all elements in support set I after renewal take out by S4, form a matrix, obtain next step the direction of search;

S5 calculates next step step-size in search γ _kand upgrade support set;

S6 obtains next step result according to the direction of search solved and step iteration, and upgrades regularization parameter;

S7 judges whether to reach stop condition, if reach, then process of reconstruction terminates, otherwise forwards step S4 to.

2. the method for claim 1, is characterized in that described step S5 comprises:

Material calculation and record the position pos+ of minimum of a value appearance;

Material calculation and record the position pos-of minimum of a value appearance;

Relatively with size, if forward step S4 to; Otherwise forward step S5 to;

By supplementary set I ^cin be arranged in pos+ position element join and support set I, and to make

Its supplementary set I is removed to by supporting the element being positioned at pos-position in set I ^c, and make

Wherein, λ is regularization parameter, c _kthe residual error relative coefficient of kth time iteration, a _ibeing that the row that in the support set I after upgrading, current all elements is corresponding in matrix A take out the matrix formed, is the i-th row of matrix A, p _kthe described direction of search, X _kit is fluorescence light source distribution vector to be solved.

3. the method for claim 1, is characterized in that judging whether that reaching stop condition comprises: the value F calculating object function after upgrading _λ(X _k), if | F _λ(X _k)-F _λ(X _k-1) |/F _λ(X _k) < T, then stop iteration, otherwise forward step S4 to and start next round iteration, wherein the value of threshold value T is 0.0001-0.01, F _λ(X _k-1) represent the object function before upgrading.