CN101396262B - Fluorescent molecule tomography rebuilding method based on linear relationship - Google Patents

Abstract

A tomography reconstruction method of fluorescence molecules based on a linearity relation belongs to the technical field of the image reconstruction in biomedicine. The reconstruction method is characterized in that the method sequentially comprises the following steps: 1. the outside surface outline is obtained by picture processing and numerical fitting based on the map data of a mouse used inexperiment; 2. the outside surface outline is led into a finite element software package for meshing; 3. every excitation light source, and finite element grid nodes which are corresponding to corresponding fluorescence monitoring points are found; 4. aiming at one given excitation light source sk, a linear relation matrix Wsk between fluorescence intensity which is detected from surface and unknown fluorescent probe distribution is established; 5. aiming at all excitation light sources, a general linear relation matrix is generated; 6. after the linear relation is obtained, a algebra reconstruction method is used for iterating the fluorescent probe distribution which approaches the reality. The method is characterized in that non-uniform optical parameters can be handled and the precise linear relation between the fluorescence intensity which is detected from surface and the unknown fluorescent probe distribution is rapidly established.

Description

A kind of fluorescent molecule tomography rebuilding method based on linear relationship

Technical field

The invention belongs to the technical field of image reconstruction in biomedicine, relate to a kind of 3D image reconstruction method.

Background technology

Molecular imaging is a popular direction of biomedical sector in recent years.Molecular imaging broadly can be defined as the qualitative and quantitative measurement of the bioprocess under the condition of living organism being carried out at molecule and cellular level [R.Weissleder, Radiology, 2001].Fluorescent molecular tomography (FMT) is an importance in the molecular imaging.It is applied in biology widely, and medical research field makes the fluorescence molecule imaging become a kind of can be used in accurately, the instrument of 3 dimensional imagings from a kind of instrument that is used as cut sections for microscopic examination qualitatively.

The key of fluorescent molecular tomography is that the data for projection that obtains a plurality of angles is used for rebuilding.Different based on the equipment of high-energy ray with CT etc., (400nm～900nm) propagation in biological tissue is not to follow straightline propagation to light.This propagation path of light can be similar to diffusion equation and describe.Utilization just can reconstruct the CONCENTRATION DISTRIBUTION of fluorescent probe based on the method for reconstructing of diffusion equation model tomography from a plurality of angle data for projection.

People such as R.B.Schulz [R.B.Schulz, IEEE Transactions on Medical Imaging, 2004] proposed a method for reconstructing, can access the surface and detect fluorescence intensity with the linear relationship between the inner phosphor concentration distribution based on analytic solutions.The optical parametric of this method hypothesis imaging object is uniform.Yet the optical parametric of the experimental rat of actual imaging is heterogeneous.Thereby the inaccuracy that optical parametric is estimated can cause the very big distortion of reconstructed results.People such as S.Bjoern [S.Bjoern, Advances in Medical Engineering, 2007] cause distortion to study to the optical parametric estimation is inaccurate.

Fair and clear professor leader's medical imaging laboratory has been delivered the method for reconstructing [X.Song based on Finite Element Method in 2007, Optics Express, 2007], can handle optical parametric heterogeneous, in method, calculate simultaneously surperficial detected fluorescence intensity with the linear relationship between the unknown fluorescent probe distribution, and then rebuild.These methods have accurately, do not need to find the solution repeatedly the characteristics of gradient and two gradients in iterative process, thereby computational speed are in 5-10 minute.Yet, in the process of calculating linear relationship, the operation of inverting at extensive sparse matrix being arranged, the computational speed ground that has limited method further improves.

Summary of the invention

The objective of the invention is for the fluorescent molecular tomography system provides a kind of method for reconstructing, be characterized in the linear relationship between can set up the detected fluorescence intensity in experimental rat surface rapidly and accurately distributes with unknown fluorescent probe, realize rebuilding.In this process, method can be handled irregular imaging object and optical parametric heterogeneous.In the method:,, on the basis of cutting apart good experimental rat collection of illustrative plates that obtains, utilize the bianry image contour line to extract and the method for numerical fitting obtains the irregular outer surface profile of experimental rat as CT at other imaging pattern; On the experimental rat collection of illustrative plates, determine the intravital non-homogeneous optical parametric of experimental rat by the method for tabling look-up, comprise absorbing and scattering parameter; Right at each excitation source-fluoroscopic examination, utilize the method for a vector-sparse matrix product to obtain the fluorescence intensity at this test point place with the linear relationship between the unknown fluorescent probe distribution.

The invention is characterized in that described method realizes successively according to the following steps with a computer:

Step (1). import in the described computer with the stratified collection of illustrative plates of Mus cut apart the biotic experiment that finishes with the CT imaging device, this computer obtains outer surface profile by following image processing method, and the statement surface profile is modeled as solid:

Step (1.1). to each tomographic image I of described collection of illustrative plates, carry out binaryzation with an identical threshold value thr=0.5: get 1 greater than selected threshold value thr person, be less than or equal to this threshold value thr person and get 0, use tissue in this collection of illustrative plates and background described " 1 " and " 0 " separately successively, obtain the binary map BI of each tomographic image:

$\mathrm{BI}[\mathrm{row},\mathrm{col}]=\left\{\begin{array}{cc}1,& I[\mathrm{row},\mathrm{col}]>\mathrm{thr},\\ 0,& I[\mathrm{row},\mathrm{col}]\≤\mathrm{thr},\end{array}\right.$

Wherein, row represents the row of this image I, and col represents the row of this image I,

Step (1.2). obtain the outer contour S of described binary map BI with the method for binary map outer contour extraction:

To each the pixel BI[row on this binary map BI, col] the following operation of execution:

If: BI[row, col] be 1, and its 8 adjacent pixels have at least one to be 0, this pixel BI[row then, col] on described outer contour S,

If: BI[row, col] be 0, and its 8 adjacent pixels have at least one to be 1, this pixel BI[row then, col] on described outer contour S,

Step (1.3). the described outer contour S that step (1.2) is obtained corresponds to cartesian coordinate from pixel coordinate, and then any one point [row, col] on this outer contour S corresponds on the described cartesian coordinate [x, y] and is,

x＝(row-row ₀)*Δ，

y＝(col-col ₀)*Δ，

Wherein, point [row ₀, col ₀] be the center of this outer contour, Δ is the physical dimension of each pixel, represents with cm,

Step (1.4). this cartesian coordinate [x, y] that step (1.3) is obtained corresponds on the polar coordinate [α, R], α is an angle, R is polar radius, on the described outer contour S have a few, each the point on described radius R do fitting of a polynomial with angle [alpha] by following formula:

$R_n=\underset{c=0}{\overset{c=12}{\mathrm{\Σ}}}{p}_c{\mathrm{\α}}_n^c,$

C is the match number of times, c=0, and 1 ... 12, p _cBe the multinomial coefficient that match obtains, with { p _cExpression, and n=1,2 ... S_num, S_num are the numbers of being had a few on the described contour line S,

At [0,2 π] evenly selected α _ num=30 angle { α ₁..., α _{α _ num}, with multinomial { p _cRadius { R on these angles of match ₁..., R _{α _ num}, make described each layer outer contour S of each tomographic image be expressed as:

{{α ₁，…，α _{α_num}}，{R ₁，…，R _{α_num}}，{p ₀，…，p ₁₂}}，

With S} represents the contour line that obtains on the I on all tomographic images,

Step (1.5). the outer contour at described all tomographic image I { among the S}, according to highly evenly selecting M shell, is represented the outer contour of described all tomographic image I with this, then adjacent two layers H ₁＜H ₂Between, represent to be expressed as the coordinate of any one deck outer contour with following formula with radius R:

$R_1=P_{H_1}\left(\mathrm{\α}\right),$

$R_2=P_{H_2}\left(\mathrm{\α}\right),$

H ₁＜H ₂，

R＝R ₁*(H ₂-H)/(H ₂-H ₁)+R ₂(H-H ₁)/(H ₂-H ₁)，

Wherein, R ₁, R ₂Be respectively in described height H ₁, H ₂The pairing radius of layer angle [alpha],

Be respectively in height H ₁, H ₂Pairing 12 order polynomial of layer,

Step (1.6). the three-dimensional outer contour that step (1.4)～(1.5) are obtained is transformed into three-dimensional cartesian coordinate system [x, y, z] from polar coordinate [α, R, H], has obtained the geometric model G under the described three-dimensional cartesian coordinate system;

Step (2). a biotic experiment that carry out the living body fluorescent molecular tomographic is placed on the rotation platform with Mus, the exciting light that laser instrument sends shines described living body biological experimental rat from a side, excitation goes out fluorescence, the fluorescence that is inspired in a relative side by CCD phase machine testing, again in the testing result input computer, the branch angle is by making described rotation platform rotation, gather the fluoroscopic image that the different rotary angle is excited, be equivalent to realize that the exciting light point source when different angles excites;

Step (3). the geometric model G that step (1) is obtained imports in the finite element software bag, simultaneously the pairing point of described each exciting light sources is added on the node of described geometric model, reuse finite element software bag is done grid and is divided, be separated into a tetrahedral grid Fem_Mesh, all node automatic numberings of this tetrahedral grid are i=1,, N, N are total grid node numbers;

Step (4). determine the pairing finite element grid node on described tetrahedral grid Fem_Mesh of each exciting light point source s according to the following steps:

Step (4.1). the laser beam that impinges perpendicularly on described living body biological experimental rat is modeled as next free path ltr=1/ (3 μ ' that freely transmit of incident illumination skin _Se(r _s)) position r _sThe exciting light point source δ (r-r at place _s), this δ (r-r _s) represent at position r _sAn impulse function at place, μ ' _Se(r _s) be r _sThe scattering coefficient of place's excitation wavelength,

Step (4.2). go up the finite element grid node i for each described tetrahedral grid Fem_Mesh, if r _s=r _i, then this finite element grid node i is exactly the pairing finite element grid node of exciting light point source s;

Step (5). determine the finite element grid node i of described each exciting light point source s on the pairing described tetrahedral grid Fem_Mesh of the lip-deep test point of finite element grid on the described Fem_Mesh:

If: described surperficial node i satisfies following relation:

r _i∈ { angle_low～angle_high} and r _i∈ height_low～height_high},

Described angle_low, angle_high, height_low, height_high is setting value, then described node i belongs to the grid node of described exciting light point source s correspondence;

Step (6). after obtaining described each test point of step (5), the fluoroscopic image that utilizes described CCD phase machine testing to obtain obtains the fluoroscopic image that each test point sends;

Step (7). at described given position r _sThe place exciting light point source s, according to following steps calculate the described surperficial test point of corresponding described exciting light point source s d1, d2 ..., the last detected fluorescence intensity of dL} with the unknown fluorescence intensity between linear relationship, wherein L is the test point sum:

Step (7.1). be calculated as follows N * N dimension stiffness matrix K _eIn element K _e[i, j], this stiffness matrix K _eBe biological tissue excitation wavelength on the finite element node on the described tetrahedral grid Fem_Mesh absorption and the embodiment of scattering coefficient, this stiffness matrix K _eRow and column represent with i and j respectively, all grid node N on the corresponding described Fem_Mesh of described capable i arrange the row that forms according to number order, all grid node N on the corresponding described Fem_Mesh of described row j arrange the row that form according to number order, down together, and wherein any element K _e[i, j] is expressed as:

$K_e[i,j]={\∫}_{\mathrm{\Ω}}(D_e\left(r\right)\▿{\mathrm{\ψ}}_i\left(r\right)\·\▿{\mathrm{\ψ}}_j\left(r\right)+{\mathrm{\μ}}_{\mathrm{ae}}\left(r\right){\mathrm{\ψ}}_i\left(r\right){\mathrm{\ψ}}_j\left(r\right))\mathrm{dr}+\frac{1}{2q}{\∫}_{\∂\mathrm{\Ω}}{{\mathrm{\ψ}}_i\left(r\right)\mathrm{\ψ}}_j\left(r\right)\mathrm{dr},$

Wherein, Ω represents the inner corresponding area of space of described geometric model G,

Represent the surface-boundary of described geometric model G, D _e(r) be the diffusion coefficient that is positioned at r position excitation wavelength, described D _e(r)=1/ (3 μ ' _Se(r)), μ ' _Se(r) for being positioned at the scattering coefficient that is organized in described excitation wavelength of r position, be given value, μ _Ae(r) for being positioned at the absorptance that is organized in described excitation wavelength of r position, be given value, ψ _i(r) and ψ _j(r) be that line number is i, column number is the linear shape function of Lagrange of the finite element node correspondence of j, provide by described finite element software bag, $\▿{\mathrm{\ψ}}_i\left(r\right)=[\∂{\mathrm{\ψ}}_i\left(r\right)/\∂x,{\∂\mathrm{\ψ}}_i\left(r\right)/\∂y,{\∂\mathrm{\ψ}}_i\left(r\right)/\∂z]$ With $\▿{\mathrm{\ψ}}_j\left(r\right)=[\∂{\mathrm{\ψ}}_j\left(r\right)/\∂x,{\∂\mathrm{\ψ}}_j\left(r\right)/\∂y,{\∂\mathrm{\ψ}}_j\left(r\right)/\∂z]$ Be respectively the linear shape function ψ of described Lagrange _i(r), ψ _j(r) gradient, q=0.1511,

Step (7.2). use Matlab software kit function PCG to be calculated as follows the intensity distributions vector Φ of described exciting light point source s on described finite element grid Fem_Mesh node _e, Φ _eBe N * 1 column vector:

K _eΦ _e＝B _s，

N * 1 column vector B _sColumn number be the element B of j _s[j] represents with following formula:

The intensity distributions vector Φ of cicada on node _e, obtain the intensity substep at r place, optional position ${\mathrm{\Φ}}_e\left(r\right)=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\Φ}}_e\left[j\right]$

Step (7.3). be calculated as follows N * N matrix F _sElement F _s[i, j], this matrix F _sBe described excitating light strength Φ _e(r) to the embodiment of excitation on the finite element node on the described tetrahedral grid Fem_Mesh of fluorescence, wherein:

$F_s[i,j]={\∫}_{\mathrm{\Ω}}{\mathrm{\Φ}}_e\left(r\right){\mathrm{\ψ}}_i\left(r\right){\mathrm{\ψ}}_j\left(r\right)\mathrm{dr},$

Step (7.4). be calculated as follows N * N dimension stiffness matrix K _mIn element K _m[i, j], this stiffness matrix K _mThe expression biological tissue excitation wavelength on the finite element node on the described tetrahedral grid Fem_Mesh absorption and the embodiment of scattering coefficient, wherein:

$K_m[i,j]={\∫}_{\mathrm{\Ω}}(D_m\left(r\right)\▿{\mathrm{\ψ}}_i\left(r\right)\·\▿{\mathrm{\ψ}}_j\left(r\right)+{\mathrm{\μ}}_{\mathrm{am}}\left(r\right){\mathrm{\ψ}}_i\left(r\right){\mathrm{\ψ}}_j\left(r\right))\mathrm{dr}+\frac{1}{2q}{\∫}_{\∂\mathrm{\Ω}}{\mathrm{\ψ}}_i\left(r\right){\mathrm{\ψ}}_j\left(r\right)\mathrm{dr},$

Wherein, D _m(r) be the diffusion coefficient that is positioned at r position wavelength of fluorescence, described $D_m\left(r\right)=1/\left(3{\mathrm{\μ}}_{\mathrm{sm}}^{\′}\left(r\right)\right),{\mathrm{\μ}}_{\mathrm{sm}}^{\′}\left(r\right)$ For being positioned at the scattering coefficient that is organized in wavelength of fluorescence of r position, be given value, μ _Am(r) for being positioned at the absorptance that is organized in wavelength of fluorescence of r position, be given value,

Step (7.5). imagination is when the point source of a wavelength of fluorescence is placed on described test point dl place, dl ∈ d1, d2 ... dL}, the intensity distributions vector Φ of fluorescence field on the finite element node on all described tetrahedral grid Fem_Mesh that uses Matlab software kit function PCG calculating to be produced _Dl, vectorial Φ _DlBe N * 1 column vector, wherein:

K _mΦ _dl＝B _dl，

N * 1 column vector B _DlColumn number be the element B of j _Dl[j] represents with following formula:

Step (7.6). according to following formula generate described each test point d1, d2 ..., the linear relationship between the detected fluorescence intensity in dL} place distributes with unknown fluorescent probe is used matrix W _sExpression:

$W_s=\left[\begin{array}{c}{\mathrm{\&Phi;}}_{\mathrm{<figure-callout\; id="1"\; label="d"\; filenames="H200810225444XE00032.png"\; state="\{\{state\}\}">d</figure-callout>}1}\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ {\mathrm{\&Phi;}}_{\mathrm{dL}}\end{array}\right]F_s;$

Step (8). at the exciting light point source sk that spatial order distributes, k=1 ..., N _s, and corresponding described test point, make up total linear relationship matrix W, following Φ according to following formula _{M, meas}Be the detected fluorescence of each exciting light point source correspondence, be given value:

${\mathrm{\&Phi;}}_{m,\mathrm{meas}}=\left[\begin{array}{c}{W}_{s1}\\ .\\ .\\ .\\ W_{s{N}_s}\end{array}\right]U=\mathrm{WU},$

N * 1 column vector U is the unknown distribution of fluorescent probe on described finite element grid Fem_Mesh node, Φ _{M, meas}Be M * 1 column vector, M is total test point number of all exciting light point source correspondences, and W is M * N matrix;

Step (9). by following algebraic expression, the total linear relationship matrix W that obtains according to step (8) is found the solution the fluorescent probe distribution U of position:

$U{\left[j\right]}^{u\_\mathrm{iter}+1}=U{\left[j\right]}^{u\_\mathrm{iter}}+\mathrm{\&lambda;}\underset{{\mathrm{<figure-callout\; id="1"\; label="t"\; filenames="H200810225444XE00032.png"\; state="\{\{state\}\}">t</figure-callout>}}_1=1}{\overset{M}{\mathrm{\&Sigma;}}}\frac{V\left[t_1\right]-\underset{t_2=1}{\overset{N}{\mathrm{\&Sigma;}}}E[{\mathrm{<figure-callout\; id="1"\; label="t"\; filenames="H200810225444XE00032.png"\; state="\{\{state\}\}">t</figure-callout>}}_1,t_2]U{\left[t_2\right]}^{u\_\mathrm{inter}}}{\underset{t_2=1}{\overset{N}{\mathrm{\&Sigma;}}}W[{\mathrm{<figure-callout\; id="1"\; label="t"\; filenames="H200810225444XE00032.png"\; state="\{\{state\}\}">t</figure-callout>}}_1,t_2]W({\mathrm{<figure-callout\; id="1"\; label="t"\; filenames="H200810225444XE00032.png"\; state="\{\{state\}\}">t</figure-callout>}}_1,t_2)}W[{\mathrm{<figure-callout\; id="1"\; label="t"\; filenames="H200810225444XE00032.png"\; state="\{\{state\}\}">t</figure-callout>}}_1,j],$

Wherein:

U[j] ^U_iter+1, U[j] ^U_iterBe respectively the value of j element when u_iter+1 and u_iter iteration among the U, U[t ₂] ^U_iterBe t among the U ₂The value of individual element when the u_iter time iteration, W[t ₁, t ₂] be matrix W t ₁OK, t ₂The value of row, W[t ₁, j] and be matrix W t ₁OK, the value of j row, Φ _{M, meas}Represent V[t with V ₁] the V t that is ₁Individual element, λ=0.1, U iteration initial value is made as U=0.

Described angle_low=-48 °, angle_high=48 °,

Described height_low=-1.5cm, height_high=1.5cm.

In order to verify the reliability of this method for reconstructing, experimental selection one be enclosed within the emulated data that obtains on the experimental rat collection of illustrative plates.From reconstructed results as can be seen, the fluorophor of reconstruction has good location and quantitative effect.

Description of drawings

The mechanism of Fig. 1 fluorescence molecule imaging: 1 the representative be fluorescent optical filter.

The mechanism of Fig. 2 fluorescent molecular tomography.

The outer surface profile of Fig. 3 experimental rat collection of illustrative plates.

Fig. 4 noncontact fluorescent molecular tomography system construction drawing: 1 the representative be fluorescent optical filter.

The part that Fig. 5 experimental rat outer surface profile is divided the finite element grid that obtains shows.

Emulated data reconstructed results on Fig. 6 experimental rat spectrum model.

The FB(flow block) of Fig. 7 computer program.

The specific embodiment

The mechanism of fluorescence molecule imaging can illustrate with Fig. 1: the in-house fluorescent probe of the excitation that exciting light sources sends (albumen of labelling, tumor, fluorescent dye etc.) emitting fluorescence, the fluorescence that propagates into the surface is detected by the CCD photographing unit.In Fig. 1, be positioned at position r _sThe exciting light that sends of exciting light point source in the experimental rat body, propagate, represent with solid arrow.After exciting light runs into fluorescent probe, can be absorbed by fluorescent probe, fluorescent probe is excited then, launches the fluorescence of higher wavelength.Fluorescence is propagated in the experimental rat body, and the with dashed lines arrow is represented.The fluorescence that arrives excitation source offside experimental rat surface can detect with the CCD photographing unit.Before the CCD photographing unit fluorescent optical filter is arranged, guarantee that the fluorescence that only is excited is just detected by CCD.Because biological tissue has the strong scattering characteristic, thereby exciting light is not straightline propagation in the experimental rat body with the fluorescence that is excited.The propagation in biological tissue of the mechanism of this imaging and light can be represented with two coupled partial differential equations that are called diffusion equation:

$\left\{\begin{array}{c}-\&dtri;\&CenterDot;[D_e\left(r\right)\&dtri;{\mathrm{\&Phi;}}_e\left(r\right)]+{\mathrm{\&mu;}}_{\mathrm{ae}}\left(r\right){\mathrm{\&Phi;}}_e\left(r\right)=\mathrm{\&delta;}(r-r_s),\\ -\&dtri;\&CenterDot;[D_m\left(r\right)\&dtri;{\mathrm{\&Phi;}}_m\left(r\right)]+{\mathrm{\&mu;}}_{\mathrm{am}}\left(r\right){\mathrm{\&Phi;}}_m\left(r\right)={\mathrm{\&Phi;}}_e\left(r\right)U\left(r\right),\end{array}\right.$

Wherein: subscript e represents excitation wavelength, and subscript m is represented wavelength of fluorescence.Φ _e(r) and Φ _m(r) represent the exciting light and the intensity of fluorescence that is excited to distribute respectively.R represents the locus coordinate in the whole geometric model.In this article, (...) representative function operation, and [...] expression vector, image or matrix manipulation.For example, Φ _e(r) represent the value of excitation light high field at the r place, and Φ _e[j] then represents exciting light high field vector Φ _eJ element.The propagation of light in tissue is subjected to the influence of two factors, is called to absorb and scattering.Absorption characteristic can light absorbing part energy, and scattering properties can change the direction of propagation of light, make light in tissue can not straightline propagation.The absorption of tissue and scattering properties can be inequality at different wave length.D _e(r) and D _m(r) be respectively to be positioned at the exciting light at r place and the diffusion coefficient of wavelength of fluorescence, correspond to $D_e\left(r\right)=1/\left(3{\mathrm{\&mu;}}_{\mathrm{se}}^{\&prime;}\left(r\right)\right)$ With $D_m\left(r\right)=1/\left(3{\mathrm{\&mu;}}_{\mathrm{sm}}^{\&prime;}\left(r\right)\right),$ Wherein

Be respectively that what to be organized in exciting light and wavelength of fluorescence is scattering coefficient.μ _Ae(r) and μ _Am(r) be the absorptance that is organized in exciting light and wavelength of fluorescence respectively.In fluorescent molecular tomography, the laser beam that impinges perpendicularly on the experimental rat surface is modeled as the next transmission of incidence point skin free path $\mathrm{ltr}=1/\left(3{\mathrm{\&mu;}}_{\mathrm{se}}^{\&prime;}\right)$ The position r at place _sExciting light point source δ (r-r _s).δ (r-r _s) represent at position r _sAn impulse function at place.U (r) is the distribution that is positioned at the fluorescent probe at r place.Fluorescent molecular tomography is to utilize knownly in the detected fluorescence intensity in surface location place, rebuilds the CONCENTRATION DISTRIBUTION U (r) of unknown fluorescent probe.Its principle Fig. 2 brief description: being distributed in intravital all fluorescent probes of experimental rat all has contribution to the detected fluorescence intensity at test point d1 and d2 place, but the weight of contribution has different.In the drawings, briefly represent intravital all fluorescent probes of experimental rat to distribute with fluorescent probe 1,2,3.With test point d1 is example, fluorescent probe 1,2, and 3 pairs of d1 contribution weights are respectively U1 * w1, U2 * w2, U3 * w3.U1 wherein, U2, U3 are respectively fluorescent probes 1,2,3 concentration, and w1, w2, w3 are positioned at fluorescent probe 1, and prosposition is put the fluorescent probe of unit concentration at place to the contribution weight of test point d1, is called unit contribution weight.The contribution weight w1 of unit, w2, w3 can find the solution by fluorescence molecule cross sectional reconstruction method and obtain.Right to different excitation source-fluoroscopic examination points, can calculate corresponding unit contribution weight.Utilize all unit that calculates contribution weight and detected fluorescence intensities, just can reconstruct unknown fluorescence concentration U1, U2, U3.This unit contribution weight is called as detected fluorescence intensity in this article with the linear relationship between the unknown fluorescent probe distribution.Actual method for reconstructing calculates the unit contribution weight of the fluorescent probe of the unit fluorescence concentration on the finite element grid node to detected fluorescence intensity, can rebuild the fluorescent probe CONCENTRATION DISTRIBUTION that obtains on the finite element grid node with optimization method then.

In fluorescent molecular tomography, the exciting light point source of a plurality of angles excites in proper order, detects corresponding fluorescence data for projection simultaneously, just can go out unknown fluorescent probe CONCENTRATION DISTRIBUTION by cross sectional reconstruction from a plurality of angle data for projection.In the fluorescent molecular tomography system, obtain the anatomical information of experimental rat earlier with the CT imaging device, utilize the method for medical image segmentation to cut apart the CT image then, obtain the organ collection of illustrative plates of experimental rat.Here the method for Ti Chuing can be utilized the fluorescence data for projection of a plurality of angles that measure on the basis of the experimental rat collection of illustrative plates that obtains, and reconstructs unknown fluorescent probe CONCENTRATION DISTRIBUTION.

The key step of method comprises:

Step (1) is utilized and is cut apart the experimental rat collection of illustrative plates that finishes, and obtains outer surface profile by Flame Image Process, and outer surface profile is modeled as solid, and this method realizes successively according to the following steps:

Step (1.1). in the experimental rat collection of illustrative plates that has obtained, the gray value of organizing in the collection of illustrative plates all is more than or equal to 1, and the gray value of background is 0.Thereby selected threshold value thr=0.5 just can separate tissue and background in the collection of illustrative plates.Each tomographic image I to collection of illustrative plates utilizes threshold value thr to carry out binaryzation, greater than threshold value get 1, less than threshold value get 0, obtain binary map BI:

$\mathrm{BI}[\mathrm{row},\mathrm{col}]=\left\{\begin{array}{c}1,I[\mathrm{row},\mathrm{col}]>\mathrm{<figure-callout\; id="0"\; label="thr"\; filenames="H200810225444XE00032.png"\; state="\{\{state\}\}">thr</figure-callout>}\\ 0,I[\mathrm{row},\mathrm{col})]\&le;\mathrm{thr}\end{array}\right.,$

Wherein: the row of row representative image, and the row of col representative image.

Step (1.2). obtain the outer contour S of binary map BI with the method for binary map contour line extraction.Method is as follows: to each pixel BI[row of binary map BI, col]: (1) if BI[row, col] be 1, and its 8 adjacent pixels have at least one to be 0, this pixel is on outer contour S so; (2) if BI[row, col] be 0, and its 8 adjacent pixels have at least one to be 1, this pixel is on outer contour S so.So just obtained the outer contour of a tomographic image in the collection of illustrative plates.Layers all in the collection of illustrative plates is all carried out binaryzation, carry out contour line then and extract, just can obtain the outer contour of all layers, with { S} represents.

Step (1.3). successively outer contour S is corresponded to cartesian coordinate from pixel coordinate, give these contour lines selected center [row ₀, col ₀], the physical dimension of each pixel is that Δ cm. is so with [row ₀, col ₀] be initial point, any one point [row, col] on the contour line all corresponds to geometric coordinate [x, y]:

x＝(row-row ₀)*Δ

y＝(col-col ₀)*Δ

Then cartesian coordinate [x, y] is corresponded to polar coordinate [α, R], wherein α is the angle in the polar coordinate system, and R is the radius in the polar coordinate system.To on each bar outer contour S have a few, the radius R of each point is done 12 order polynomial matches with angle [alpha] by following formula:

$R_n=\underset{c=0}{\overset{c=12}{\mathrm{\&Sigma;}}}{p}_c{\mathrm{\&alpha;}}_n^c,$

Wherein: c is the match number of times, c=0, and 1 ... 12.p _cBe the multinomial coefficient that match obtains, with { p _cExpression.N=1,2 ... S_num, wherein S_num is the number of being had a few on the contour line S.At [0,2 π] evenly selected α _ num=30 angle { α ₁..., α _{α _ num}, with multinomial { p _cFit within the radius { R on these angles ₁..., R _{α _ num}.Like this, the pairing contour line S of each tomographic image I is expressed as { { α ₁..., α _{α _ num}, { R ₁..., R _{α _ num}, { p ₀..., p ₁₂.

Step (1.4). each tomographic image can both obtain a contour line S, with { S} represents the resulting contour line of all tomographic images.{ among the S},, representing whole outline with this M shell contour line according to highly evenly selecting M shell.Between adjacent two layers, can obtain the profile coordinate by the method for interpolation.That is to say: specify an angle [alpha] and height H arbitrarily, can both correspondingly calculate its radius R.If being in, height H highly is H ₁And H ₂Two-layer (H ₁＜H ₂) between, so:

${\mathrm{<figure-callout\; id="1"\; label="R"\; filenames="H200810225444XE00032.png"\; state="\{\{state\}\}">R</figure-callout>}}_1=P_{H_1}\left(\mathrm{\&alpha;}\right)$

$R_2=P_{H_2}\left(\mathrm{\&alpha;}\right)$

R＝R ₁*(H ₂-H)/(H ₂-H ₁)+R ₂(H-H ₁)/(H ₂-H ₁)

R wherein ₁, R ₂Be respectively at H ₁, H ₂The pairing radius of layer angle [alpha],

With

Be respectively at H ₁, H ₂Pairing 12 order polynomials.

Step (1.5). the outline that obtains is transformed into three-dimensional cartesian coordinate [x, y, z] from three-dimensional polar coordinate [α, R, H], has so just obtained the geometric model G under cartesian coordinate system, as Fig. 3.

Step (2). in imaging system shown in Figure 4, experimental rat is placed on the rotation platform.Laser shines experimental rat from a side, the excitation experiment intravital fluorescence of Mus, and the fluorescence that is inspired is detected by the CCD photographing unit at offside.By rotation platform, images acquired just can realize exciting, also with regard to corresponding a plurality of exciting light point sources in different angles then.The geometric model G that step (1) is obtained imports in the finite element software bag, simultaneously the pairing point of each exciting light point source is added on the node of this geometric model G, does grid with finite element software then and divides, and is separated into tetrahedral grid Fem_Mesh.Fig. 5 has partly shown this tetrahedral grid.All node automatic numberings of this tetrahedral grid are i=1 ..., N, N are total grid node numbers.

Step (3). determine each exciting light point source s with and the pairing finite element node of all test points.In discrete tetrahedral grid Fem_Mesh, all corresponding finite element node of each exciting light point source, the node of definite exciting light point source correspondence in the following manner: to each finite element node i, if r _s=r _i, i is exactly the pairing finite element node of exciting light point source s, wherein r so _sBe the pairing position coordinates of exciting light point source s, and r _iIt is the pairing position coordinates of finite element node i.Simultaneously, to each exciting light point source s, in the angle angle_low～angle_high of its offside correspondence and height height_low～height_high scope, detect fluorescence.In method, angular range is elected-48 °～48 ° as, and altitude range is elected as in relative exciting light point source height-1.5～1.5cm.In method, require test point to be placed on the surperficial node of finite element grid, thereby, can determine all finite element nodes of test point correspondence with the following methods: to each surperficial node i, if r _i∈ { angle_low～angle_high} and r _i{ height_low～height_high}, node i belongs to the test point of this light source correspondence to ∈ so, otherwise does not belong to.To particular excitation luminous point light source, find respective detection point after, utilize the fluorescent image photograph, obtain the fluorescence intensity that measures at these test point places.

Step (4). describe the diffusion equation group that light is propagated in the fluorescent molecular tomography and on finite element grid Fem_Mesh, can be separated into a pair of system of linear equations:

$\left\{\begin{array}{c}{K}_e{\mathrm{\&Phi;}}_e=B_s\\ K_m{\mathrm{\&Phi;}}_m=F_{s}U\end{array}\right.,$

Wherein: N * 1 column vector Φ _eAnd N * 1 column vector Φ _mBe respectively the exciting light and the intensity distributions of fluorescence on grid node that be excited; N * 1 column vector U has represented the distribution of fluorescent probe on grid node; N * 1 column vector B _sRepresented exciting light point source δ (r-r _s) in the distribution of grid node.In Finite Element Method, the distribution of cicada on discrete grid node can obtain distribution on the r at an arbitrary position by the linear shape function of the pairing Lagrange of finite element node.With the excitating light strength is example, at an arbitrary position the excitating light strength Φ on the r _e(r), can be expressed as ${\mathrm{\&Phi;}}_e\left(r\right)=\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&Phi;}}_e\left[j\right]{\mathrm{\&psi;}}_j\left(r\right),$ Φ wherein _e[j] is column vector Φ _eJ element.Φ _e[j] also is that excitating light strength is in the value that is numbered on the grid node of j.ψ _j(r) be the linear shape function of the pairing Lagrange of grid node that is numbered j.K _eAnd K _mIn Finite Element Method, be called stiffness matrix.Matrix K _eBe absorption and the embodiment of scattering properties on grid node that is organized in excitation wavelength.Matrix K _mBe absorption and the embodiment of scattering properties on grid node that is organized in wavelength of fluorescence.Matrix F _sBe excitating light strength distribution Φ _e(r) to the embodiment of excitation on grid node of fluorescent probe.These three sparse matrixes that matrix all is N * N.For matrix K _e, K _mAnd F _s, their row and column represents with i and j that respectively all the node N on the corresponding Fem_Mesh of described capable i arrange the row that forms according to number order, all the node N on the corresponding Fem_Mesh of described row j arrange the row that form according to number order.At being positioned at position r _sExcitation source s, its surperficial detection node d1, d2 ..., the linear relationship between the last detected fluorescence intensity of dL} distributes with unknown fluorescent probe can obtain according to following steps:

Step (4.1). the calculated rigidity matrix K _eK _eI capable, j column element K _e[i, j] can try to achieve with following formula,

$K_e[i,j]={\&Integral;}_{\mathrm{\&Omega;}}(D_e\left(r\right)\&dtri;{\mathrm{\&psi;}}_i\left(r\right)\&CenterDot;\&dtri;{\mathrm{\&psi;}}_j\left(r\right)+{\mathrm{\&mu;}}_{\mathrm{ae}}\left(r\right){\mathrm{\&psi;}}_i\left(r\right){\mathrm{\&psi;}}_j\left(r\right))\mathrm{dr}+\frac{1}{2q}{\&Integral;}_{\&PartialD;\mathrm{\&Omega;}}{\mathrm{\&psi;}}_i\left(r\right){\mathrm{\&psi;}}_j\left(r\right)\mathrm{dr},$

Wherein: Ω has represented the inner pairing area of space of whole geometric model G;

Then represent the surface-boundary of geometric model G; Subscript e represents excitation wavelength; R is a position coordinates, D _e(r) be the diffusion coefficient that is positioned at the r place, correspond to $D_e\left(r\right)=1/\left(3{\mathrm{\&mu;}}_{\mathrm{se}}^{\&prime;}\left(r\right)\right)$ Q is relevant with the coefficient of refraction of tissue and air, generally gets q=0.1511; ψ _i(r) and ψ _j(r) be that line number is i, column number is the linear shape function of Lagrange of the finite element node correspondence of j, is provided by described finite element software bag. $\&dtri;{\mathrm{\&psi;}}_i\left(r\right)=[\&PartialD;{\mathrm{\&psi;}}_i\left(r\right)/\&PartialD;x,\&PartialD;{\mathrm{\&psi;}}_i\left(r\right)/\&PartialD;y,\&PartialD;{\mathrm{\&psi;}}_i\left(r\right)/\&PartialD;z]$ With

$\&dtri;{\mathrm{\&psi;}}_j\left(r\right)=[\&PartialD;{\mathrm{\&psi;}}_j\left(r\right)/\&PartialD;x,\&PartialD;{\mathrm{\&psi;}}_j\left(r\right)/\&PartialD;y,\&PartialD;{\mathrm{\&psi;}}_j\left(r\right)/\&PartialD;z]$ Be respectively the linear shape function ψ of described Lagrange _i(r), ψ _j(r) gradient.Sparse matrix K _eBe organized in the pairing optical parametric μ of excitation wavelength by in the finite element software bag, setting _Ae(r) and

Interface function, utilize the finite element software bag to find the solution and obtain.Optical parametric

And μ _Ae(r) obtain like this: at first, obtain the gray value of the collection of illustrative plates of position r, determine the histoorgan that r place, position belongs to; Then, table look-up and obtain the optical parametric of this histoorgan in the excitation wavelength correspondence.Following table is exactly the optical parameter value of different tissues organ at exciting light and wavelength of fluorescence:

Wherein, the unit of optical parametric is cm ^-1

Step (4.2). calculate the intensity distributions Φ of the pairing exciting light of exciting light point source s on finite element grid Fem_Mesh node _e,

K _eΦ _e＝B _s。

Φ _eJ element Φ _e[j] represented the exciting light light intensity in the value that is numbered on the finite element node of j.Column vector B _sJ element B _s[j] represent by following formula,

Utilize the pretreated conjugate gradient method of Jacobi to come this system of linear equations of iterative, use Matlab software kit function PCG to realize input matrix K _eWith vectorial B _s, output vector Φ _e

Step (4.3). compute matrix F _sF _sI capable, j column element F _s[i, j] is expressed as:

F _s[i，j]＝∫ _ΩΦ _e(r)ψ _i(r)ψ _j(r)dr，

Φ wherein _e(r) be the value of the light intensity field of exciting light point source s, be expressed as at locus r place ${\mathrm{\&Phi;}}_e\left(r\right)=\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&Phi;}}_e\left[j\right]{\mathrm{\&psi;}}_j\left(r\right),$ This accumulation calculating can be finished by the finite element software bag.F _sBy in the finite element software bag, setting interface function, utilize the finite element software bag to find the solution and obtain.This interface function is such, input position r, the excitation light intensity intensity Φ at output r place _e(r).

Step (4.4). the calculated rigidity matrix K _mMatrix K _mI is capable, j column element K _m[i, j] is expressed as,

$K_m[i,j]={\&Integral;}_{\mathrm{\&Omega;}}(D_m\left(r\right)\&dtri;{\mathrm{\&psi;}}_i\left(r\right)\&CenterDot;\&dtri;{\mathrm{\&psi;}}_j\left(r\right)+{\mathrm{\&mu;}}_{\mathrm{am}}\left(r\right){\mathrm{\&psi;}}_i\left(r\right){\mathrm{\&psi;}}_j\left(r\right))\mathrm{dr}+\frac{1}{2q}{\&Integral;}_{\&PartialD;\mathrm{\&Omega;}}{\mathrm{\&psi;}}_i\left(r\right){\mathrm{\&psi;}}_j\left(r\right)\mathrm{dr}.$

Wherein, subscript m is represented the optical wavelength that is excited.This matrix find the solution similar this compute matrix K _e, by in the finite element software bag, setting the optical parametric μ that is organized in the optical wavelength that is excited _AmWith

Interface function utilizes finite element software to find the solution and obtains, optical parametric μ _AmWith

Obtain like this: at first, judge the histoorgan that r place, position belongs to; Then, table look-up and obtain the optical parametric of this organ in the wavelength of fluorescence correspondence.

Step (4.5). suppose the point source of a wavelength of fluorescence is placed on detection node dl, dl ∈ d1, d2 ..., the distribution vector Φ of fluorescence field on the finite element node that produces calculated at the dL} place _Dl,

K _mΦ _dl＝B _dl，

N * 1 column vector B wherein _DlJ element

This system of linear equations comes iterative by the pretreated conjugate gradient method of Jacobi, sees step (4.2)

Step (4.6). generate test point { d1, the linear relationship matrix W between d2..., the detected fluorescence in dL} place distribute with unknown fluorescent probe _s,

$W_s=\left[\begin{array}{c}{\mathrm{\&Phi;}}_{\mathrm{<figure-callout\; id="1"\; label="d"\; filenames="H200810225444XE00032.png"\; state="\{\{state\}\}">d</figure-callout>}1}\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ {\mathrm{\&Phi;}}_{\mathrm{dL}}\end{array}\right]F_s.$

Can see that right for each excitation source-fluoroscopic examination, for example s-d1 passes through a vector-sparse matrix product Φ at the detected fluorescence of test point d1 with the linear relationship between the unknown fluorescent probe distribution _DlF _sObtain.

Step (5). at the light source sk that a plurality of orders excite, k=1...N _sWith and the relevant detection point, stack generates total linear relationship matrix W.According to step (4), set up the linear relationship matrix W of each excitation source sk and its respective detection point respectively _Sk, can obtain total linear relationship matrix to these matrix stacks then, as follows:

${\mathrm{\&Phi;}}_{m,\mathrm{meas}}=\left[\begin{array}{c}{W}_{s1}\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ W_{s{N}_s}\end{array}\right]U=\mathrm{WU},$

Wherein, N * 1 column vector U is the unknown distribution of fluorescent probe on finite element grid Fem_Mesh node.M * 1 column vector Φ _{M, meas}Be the pairing fluorescence intensity levels of all test points, total test point number of all exciting light point source correspondences of M wherein.W is M * N matrix.

Step (6). rebuild unknown fluorescent probe according to the linear relationship matrix and distribute, utilize algebraic reconstruction technique (ART) to realize rebuilding, Φ _{M, meas}Represent with vectorial V, when finding the solution system of linear equations WU=V, according to the following formula iteration:

$U{\left[j\right]}^{u\_\mathrm{iter}+1}=U{\left[j\right]}^{u\_\mathrm{iter}}+\mathrm{\&lambda;}\underset{{\mathrm{<figure-callout\; id="1"\; label="t"\; filenames="H200810225444XE00032.png"\; state="\{\{state\}\}">t</figure-callout>}}_1=1}{\overset{M}{\mathrm{\&Sigma;}}}\frac{V\left[t_1\right]-\underset{t_2=1}{\overset{N}{\mathrm{\&Sigma;}}}W[{\mathrm{<figure-callout\; id="1"\; label="t"\; filenames="H200810225444XE00032.png"\; state="\{\{state\}\}">t</figure-callout>}}_1,t_2]U{\left[t_2\right]}^{u\_\mathrm{iter}}}{\underset{t_2=1}{\overset{N}{\mathrm{\&Sigma;}}}W[{\mathrm{<figure-callout\; id="1"\; label="t"\; filenames="H200810225444XE00032.png"\; state="\{\{state\}\}">t</figure-callout>}}_1,t_2]W[{\mathrm{<figure-callout\; id="1"\; label="t"\; filenames="H200810225444XE00032.png"\; state="\{\{state\}\}">t</figure-callout>}}_1,t_2]}W[{\mathrm{<figure-callout\; id="1"\; label="t"\; filenames="H200810225444XE00032.png"\; state="\{\{state\}\}">t</figure-callout>}}_1,j],$

Wherein: U[j] ^U_iter+1, U[j] ^U_iterBe the value of j element when u_iter+1 and u_iter iteration among the U.U[t ₂] ^U_iterBe t among the U ₂The value of individual element when the u_iter time iteration.W[t ₁, t ₂] be matrix W t ₁OK, t ₂The value of row.W[t ₁, j] and be matrix W t ₁OK, the value of j row.V[t ₁] be V t ₁Individual element.The λ value is that the initial value of λ=0.1. iteration is got U=0.The reconstructed results of example as shown in Figure 6.Fig. 6 is that the fluorescent probe of rebuilding is distributed on that height layer of exciting light point source, coordinate x, and the isoconcentrate that the y each point draws distributes.

Claims (2)

1. the fluorescence molecule cross sectional reconstruction method based on linear relationship is characterized in that, described method realizes successively according to the following steps with a computer:

Step (1) is imported in the described computer with the stratified collection of illustrative plates of Mus cut apart the biotic experiment that finishes with the CT imaging device, and this computer obtains outer surface profile by following image processing method, and described surface profile is modeled as solid:

Step (1.1). to each tomographic image I of described collection of illustrative plates, carry out binaryzation with an identical threshold value thr=0.5: get 1 greater than selected threshold value thr person, be less than or equal to this threshold value thr person and get 0, use tissue in this collection of illustrative plates and background described " 1 " and " 0 " separately successively, obtain the binary map BI of each tomographic image:

Wherein, row represents the row of this image I, and col represents the row of this image I,

Step (1.2). obtain the outer contour S of described binary map BI with the method for binary map outer contour extraction:

To each the pixel BI[row on this binary map BI, col] the following operation of execution:

If: BI[row, col] be 1, and its 8 adjacent pixels have at least one to be 0, this pixel BI[row then, col] on described outer contour S,

If: BI[row, col] be 0, and its 8 adjacent pixels have at least one to be 1, this pixel BI[row then, col] on described outer contour S,

Step (1.3). the described outer contour S that step (1.2) is obtained corresponds to cartesian coordinate from pixel coordinate, and then any one point [row, col] on this outer contour S corresponds on the described cartesian coordinate [x, y] and is,

x＝(row-row ₀)*Δ，

y＝(col-col ₀)*Δ，

Wherein, point [row ₀, col ₀] be the center of this outer contour, Δ is the physical dimension of each pixel, represents with cm,

Step (1.4). this cartesian coordinate [x, y] that step (1.3) is obtained corresponds on the polar coordinate [α, R], α is an angle, R is polar radius, on the described outer contour S have a few, each the point on described radius R do fitting of a polynomial with angle [alpha] by following formula:

C is the match number of times, c=0, and 1 ... 12, pc is the multinomial coefficient that match obtains, with { p _cExpression, and n=1,2 ... S_num, S_num are the numbers of being had a few on the described outer contour S,

At [0,2 π] evenly selected α _ num=30 angle { α ₁..., α _{α _ num}, with multinomial { p _cRadius { R on these angles of match ₁..., R _{α _ num}, make described each layer outer contour S of each tomographic image be expressed as:

{{α ₁，…，α _{α_num}}，{R ₁，…，R _{α_num}}，{p ₀，…，p ₁₂}}，

With S} represents the outer contour that obtains on the I on all tomographic images,

Step (1.5). the contour line at described all tomographic image I { among the S}, according to highly evenly selecting M shell, is represented the outer contour of described all tomographic image I with this, then the height H of adjacent two layers ₁＜H ₂Between, represent to be expressed as the coordinate of any one deck outer contour with following formula with radius R:

H ₁＜H ₂，

R＝R ₁*(H ₂-H)/(H ₂-H ₁)+R ₂(H-H ₁)/(H ₂-H ₁)，

Wherein, R ₁, R ₂Be respectively in described height H ₁, H ₂The pairing radius of layer angle [alpha], Be respectively in height H ₁, H ₂Pairing 12 order polynomial of layer,

Step (1.6). the three-dimensional outer contour that step (1.4)～(1.5) are obtained is transformed into three-dimensional cartesian coordinate system [x, y, z] from polar coordinate [α, R, H], has obtained the geometric model G under the described three-dimensional cartesian coordinate system;

Step (2). a biotic experiment that carry out the living body fluorescent molecular tomographic is placed on the rotation platform with Mus, the exciting light that laser instrument sends shines described living body biological experimental rat from a side, excitation goes out fluorescence, the fluorescence that is inspired in a relative side by CCD phase machine testing, again in the testing result input computer, the branch angle is by making described rotation platform rotation, gather the fluoroscopic image that the different rotary angle is excited, be equivalent to realize that the exciting light point source when different angles excites;

Step (3). the geometric model G that step (1) is obtained imports in the finite element software bag, simultaneously the pairing point of described each exciting light sources is added on the node of described geometric model, reuse finite element software bag is done grid and is divided, be separated into a tetrahedral grid Fem_Mesh, all node automatic numberings of this tetrahedral grid are i=1,, N, N are total grid node numbers;

Step (4). determine the pairing finite element grid node on described tetrahedral grid Fem_Mesh of each exciting light point source s according to the following steps:

Step (4.1). the laser beam that impinges perpendicularly on described living body biological experimental rat is modeled as next free path ltr=1/ (3 μ ' that freely transmit of incident illumination skin _Se(r _s)) position r _sThe exciting light point source δ (r-r at place _s), this δ (r-r _s) represent at position r _sAn impulse function at place, μ ' _Se(r _s) be r _sThe scattering coefficient of place's excitation wavelength,

Step (4.2). go up the finite element grid node i for each described tetrahedral grid Fem_Mesh, if r _s=r _i, r _iBe the position of finite element grid node i, then this finite element grid node i is exactly the pairing finite element grid node of exciting light point source s;

Step (5). determine the finite element grid node i of described each exciting light point source s on the pairing described tetrahedral grid Fem_Mesh of the lip-deep test point of finite element grid on the described Fem_Mesh:

If: described node i satisfies following relation:

r _i∈ { angle_low～angle_high} and r _i∈ height_low～height_high},

Described angle_low, angle_high, height_low, height_high is setting value,

Then described node i belongs to the grid node of described exciting light point source s correspondence;

Step (6). after obtaining described each test point of step (5), the fluoroscopic image that utilizes described CCD phase machine testing to obtain obtains the fluoroscopic image that each test point sends;

Step (7). at given position r _sThe exciting light point source s at place, according to following steps calculate corresponding described exciting light point source s surface test point d1, d2 ..., the last detected fluorescence intensity of dL} is with the linear relationship between the fluorescence intensity of the unknown, wherein L is the test point sum:

Step (7.1). be calculated as follows N * N dimension stiffness matrix K _eIn element K _e[i, j], this stiffness matrix K _eBe biological tissue excitation wavelength on the finite element node on the described tetrahedral grid Fem_Mesh absorption and the embodiment of scattering coefficient, this stiffness matrix K _eRow and column represent with i and j respectively, all grid nodes on the corresponding described Fem_Mesh of described capable i are arranged the row that forms according to number order, all grid nodes on the corresponding described Fem_Mesh of described row j are arranged the row that form according to number order, down together, and wherein any element K _e[i, j] is expressed as:

Wherein, Ω represents the inner corresponding area of space of described geometric model G, Represent the surface-boundary of described geometric model G, D _e(r) be the diffusion coefficient that is positioned at r position excitation wavelength, described D _e(r)=1/ (3 μ ' _Se(r)), μ ' _Se(r) for being positioned at the scattering coefficient that is organized in described excitation wavelength of r position, be given value, μ _Ae(r) for being positioned at the absorptance that is organized in described excitation wavelength of r position, be given value, ψ _i(r) and ψ _j(r) be that line number is i, column number is the linear shape function of Lagrange of the finite element node correspondence of j, provide by described finite element software bag,

With

Be respectively the linear shape function ψ of described Lagrange _i(r), ψ _j(r) gradient, q=0.1511,

Step (7.2). use Matlab software kit function PCG to be calculated as follows the intensity distributions vector Φ of described exciting light point source s on described finite element grid Fem_Mesh node _e, Φ _eBe N * 1 column vector:

K _eΦ _e＝B _s，

N * 1 column vector B _sColumn number be the element B of j _s[j] represents with following formula:

The intensity distributions vector Φ of cicada on node _e, obtain the intensity substep at r place, optional position

,

Step (7.3). be calculated as follows N * N matrix F _sElement F _s[i, j], this matrix F _sBe described excitating light strength Φ _e(r) to the embodiment of excitation on the finite element node on the described tetrahedral grid Fem_Mesh of fluorescence, wherein:

F _s[i，j]＝∫ _ΩΦ _e(r)ψ _i(r)ψ _j(r)dr，

Step (7.4). be calculated as follows N * N dimension stiffness matrix K _mIn element K _m[i, j], this stiffness matrix K _mThe expression biological tissue excitation wavelength on the finite element node on the described tetrahedral grid Fem_Mesh absorption and the embodiment of scattering coefficient, wherein:

Wherein, D _m(r) be the diffusion coefficient that is positioned at r position wavelength of fluorescence, described D _m(r)=1/ (3 μ ' _Sm(r)), μ ' _Sm(r) for being positioned at the scattering coefficient that is organized in wavelength of fluorescence of r position, be given value, μ _Am(r) for being positioned at the absorptance that is organized in wavelength of fluorescence of r position, be given value,

Step (7.5). imagination is when the point source of a wavelength of fluorescence is placed on test point dl place, dl ∈ d1, d2 ... dL}, the intensity distributions vector Φ of fluorescence field on the finite element node on all described tetrahedral grid Fem_Mesh that uses Matlab software kit function PCG calculating to be produced _Dl, vectorial Φ _DlBe N * 1 column vector, wherein:

K _mΦ _dl＝B _dl，

N * 1 column vector B _DlColumn number be the element B of j _Dl[j] represents with following formula:

Step (7.6). according to following formula generate described each test point d1, d2 ..., the linear relationship between the detected fluorescence intensity in dL} place distributes with unknown fluorescent probe is used matrix W _sExpression:

Step (8). at the exciting light point source sk that spatial order distributes, k=1 ..., N _s, and corresponding described test point, make up total linear relationship matrix W, following Φ according to following formula _{M, meas}Be the detected fluorescence of each exciting light point source correspondence, be given value:

N * 1 column vector U is the unknown distribution of fluorescent probe on described finite element grid Fem_Mesh node, Φ _{M, meas}Be M * 1 column vector, total test point number of all exciting light point source correspondences of M, W are M * N matrixes;

Step (9). by following algebraic expression, the total linear relationship matrix W that obtains according to step (8) is found the solution the fluorescent probe distribution U of position:

Wherein:

U[j] ^U_iter+1, U[j] ^U_iterBe respectively the value of j element when u_iter+1 and u_iter iteration among the U, U[t ₂] ^U_iterBe t among the U ₂The value of individual element when the u_iter time iteration, W[t ₁, t ₂] be matrix W t ₁OK, t ₂The value of row, W[t ₁, j] and be matrix W t ₁OK, the value of j row, Φ _{M, meas}Represent V[t with V ₁] the V t that is ₁Individual element, λ=0.1, U iteration initial value is made as U=0.

2. a kind of fluorescence molecule cross sectional reconstruction method based on linear relationship according to claim 1 is characterized in that,

Described angle_low=-48 °, angle_high=48 °,

Described height_low=-1.5cm, height_high=1.5cm.

Images