CN104679946A - Voxel-based perturbing fluorescent light Monte Carlo simulation method - Google Patents

Abstract

The invention relates to a voxel-based perturbing fluorescent light Monte Carlo simulation method. The voxel-based perturbing fluorescent light Monte Carlo simulation method includes the steps of representing an incident light source in an assembly of a set number of photons, and determining the initial position and the direction of the light source and the position of a detector; putting in an exciting light photon, tracking the exciting light photon transmitted in a biological tissue, computing the percentage of the exciting light photon being converted into a fluorescent light photon, computing the weight of fluorescent light photon received by the detector along the exciting light path, and saving photon path information; when the absorbance of a fluorophore changes slightly, computing the weight of the fluorescent light photon on the detector directly according to the saved photon path information. The voxel-based perturbing fluorescent light Monte Carlo simulation method has the advantages that transmitting time of a simulative photon in the biological tissue is saved greatly, high computing efficiency is achieved, and computational efficiency of fluorescent tomography reconstruction can be improved greatly.

Description

A kind of perturbation fluorescence Monte-Carlo Simulation Method based on voxel

Technical field

The invention belongs to mathematical simulation and biomedical engineering field, relate to a kind of perturbation fluorescence Monte-Carlo Simulation Method based on voxel.

Background technology

Most of biological tissue is the three-dimensional turbid media of high scattering, sets up a high precision and high efficiency computing method and has great significance for the quantitative accuracy of fluorescence fault imaging.Monte Carlo is a kind of discrete statistical methods based on random sampling procedure.Be compared to other method, Monte Carlo method can simulate random geometry, the photon transport process under boundary condition and optical parametric.Due to its applicability widely, it transports the most direct of actual physics process as a simulated photons, the most effectively and the most believable method.Thus, it becomes the goldstandard evaluating other application-specific method.

A.J.Welch describes the rule that fluorescence excites in biological tissues and propagates the earliest, and proposes standard fluorescence Monte Carlo method [1].In the stratiform turbid media that semiinfinite is large, the simulation result of the method has been proved to be accurately [2].Liebert proposes the perturbation fluorescence Monte Carlo method [3] applying to stratiform turbid media.In fluorescence tomography rebuilding, when fluorophore absorption coefficient changes, often need to carry out once complete standard fluorescence Monte-Carlo Simulation to obtain result, calculated amount is often quite large; And the mechanics of biological tissue information provided by imaging means such as CT and MRI is all based on voxel.The method that A.J.Welch proposes can cause larger computational burden to computing machine, and the method that Liebert proposes can only be applicable to stratified model, and above method all cannot meet the demand of fluorescence tomography rebuilding completely.Therefore invent the Monte-Carlo Simulation Method based on voxel that a kind of utilization once calculates fluorescence weight on the direct calculating detector of the photon path information of preservation, greatly can save time, and meet this demand.

[1]Welch A J，Gardner C，Richards-Kortum R，et al.Propagation of fluorescent light[J].Lasers in surgery and medicine，1997，21(2)：166-178.

[2]Vishwanath K，Pogue B，Mycek M A.Quantitative fluorescence lifetime spectroscopy in turbid media：comparison of theoretical，experimental and computational methods[J].Physics in Medicine and Biology，2002，47(18)：3387.

[3]Liebert A，Wabnitz H，Zolek N，et al.Monte Carlo algorithm for efficient simulation of time-resolved fluorescence in layered turbid media[J].Optics express，2008，16(17)：13188-13202.

Summary of the invention

The object of the invention there are provided a kind of when fluorophore absorption coefficient subtle change, calculates the Monte-Carlo Simulation Method of the upper fluorescence weight of detection fast.The method is by once calculating the photon path information of preservation, and secondary calculating directly utilizes the photon path information of preservation to calculate fluorescence, greatly saves the time that simulated photons is transmitted in biological tissues.

Based on a perturbation fluorescence Monte-Carlo Simulation Method for voxel, it is characterized in that comprising the following steps:

(1) target biological tissue is determined, three-dimensional dividing is carried out to target biological tissue, build a three-dimensional voxel model, set up a 3-dimensional digital matrix, voxel one_to_one corresponding in each element and three-dimensional voxel model in 3-dimensional digital matrix, a kind of biological tissue of numerical identity of each element, arranges the optical property parameter of biological tissue: absorption coefficient, scattering coefficient, fluorophore absorption coefficient, refractive index and anisotropy factor;

(2) incident light source is characterized by the set of setting number photon, determine the initial position of light source and the position of direction and detector, throw in exciting light photon, follow the trail of exciting light photon to transmit in biological tissues, calculate exciting light photon and be converted into fluorescent photon ratio, fluorescent photon weight received on the calculating detector of exciting light path, preserves the routing information of photon;

(3) when the fluorophore absorption coefficient in biological tissue changes, the photon path information of preserving is utilized, fluorescent photon weight received on calculating detector.

In step (2), exciting light photon is converted in fluorescent photon process, assuming that fluorescence is anisotropic scattering, fluorescent scattering direction is consistent with exciting light scattering direction.

Assuming that the path of exciting light photon and fluorescent photon random walk is in biological tissues identical in step (2).

Assuming that exciting light goes out fluorescence in fluorescence area continuous agitation in step (2).

Step (2) is specifically carried out according to the following steps:

(2.1) follow the trail of the transmission in biological tissues of exciting light photon, along exciting light path, utilize formulae discovery exciting light photon weight w (r _s, r);

$w(r_s,r)=w_0\exp (-\underset{j=1}{\overset{p_i}{\mathrm{\Σ}}}({\mathrm{\μ}}_a^{\mathrm{ex}}\left(r_j\right)+{\mathrm{\μ}}_{\mathrm{af}}\left(r_j\right))l\left(r_j\right))$

In formula: r _sfor light source position, r are fluorescence excitation position;

for absorption coefficient, the μ of exciting light _affor fluorophore absorption coefficient;

R _j(j=1 ... p _i) for exciting light photon is from r _spath between the jth-1 time and jth time scattering events of r;

W ₀for the initial weight of exciting light photon, l are at r _jthe path that microcell experiences;

(2.2) if excitation photon fluorescence area is in biological tissues absorbed and excites generation fluorescence, calculate according to fluorophore coefficient and quantum efficiency the Probability p (r) that exciting light is converted into fluorescence:;

p(r)＝η(1-exp(-μ _afl(r)))

In formula: η is quantum efficiency;

(2.3) along exciting light path, formulae discovery fluorescent photon weight w (r, r is utilized _d);

$w(r,r_d)={w_0}^{\′}\exp (-\underset{j=p_i+1}{\overset{q_i}{\mathrm{\Σ}}}({\mathrm{\μ}}_a^{\mathrm{em}}\left(r_j\right)+{\mathrm{\μ}}_{\mathrm{af}}\left(r_j\right))l\left(r_j\right))$

In formula: r _dfor detector position;

for the absorption coefficient of fluorescence;

R _j(j=p _i+ 1 ... .q _i) for fluorescent photon is from r to r _djth-1 time and jth time scattering events between path;

W ₀' for the initial weight of fluorescent photon, l be at r _jthe path that microcell experiences;

(2.4) follow the trail of all fluorescent photons, until fluorescent photon effusion tissue or dead, the fluorescent photon weight on calculating detector, preserves the routing information from source to detector photon.Source position is at r _splace, excite position at r place, detector is at r _dplace receives all fluorescent photon weights and W (r _s, r _d, r) be:

$\begin{array}{c}W(r_s,r_d,r)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\exp (-\underset{j=1}{\overset{p_i}{\mathrm{\Σ}}}({{\mathrm{\μ}}_a}^{\mathrm{ex}}\left(r_j\right)+{\mathrm{\μ}}_{\mathrm{af}}\left(r_j\right))l_i\left(r_j\right))\×(1-\exp \left({-\mathrm{\μ}}_{\mathrm{af}}\left(r\right)l_i\left(r\right)\right)\×\\ \exp (-\underset{j=p_i+1}{\overset{q_i}{\mathrm{\Σ}}}({{\mathrm{\μ}}_a}^{\mathrm{em}}\left(r_j\right)+{\mathrm{\μ}}_{\mathrm{af}}\left(r_j\right))l_i\left(r_j\right))\mathrm{\η}\end{array}$

At non-fluorescence region μ _af(r _j) be 0;

Step (3) is specifically carried out according to the following steps:

(3.1) determine the index value of the voxel of fluorophore absorption coefficient change in biological tissue, and extract the path that in this voxel, fluorescent photon is walked;

(3.2) utilize the fluorescent photon routing information preserved, substitute into the fluorescent photon weight on the direct calculating detector of following formula.

$\begin{array}{c}W(r_s,r_d,r)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\exp (-\underset{j=1}{\overset{p_i}{\mathrm{\Σ}}}({{\mathrm{\μ}}_a}^{\mathrm{ex}}\left(r_j\right)+{\mathrm{\μ}}_{\mathrm{af}}\left(r_j\right))l_i\left(r_j\right))\×(1-\exp \left({-\mathrm{\μ}}_{\mathrm{af}}\left(r\right)l_i\left(r\right)\right)\×\\ \exp (-\underset{j=p_i+1}{\overset{q_i}{\mathrm{\Σ}}}({{\mathrm{\μ}}_a}^{\mathrm{em}}\left(r_j\right)+{\mathrm{\μ}}_{\mathrm{af}}\left(r_j\right))l_i\left(r_j\right))\mathrm{\η}\end{array}$

The inventive method counting yield is high, greatly can improve the counting yield of fluorescence tomography rebuilding.

Accompanying drawing explanation

Fig. 1 is basic flow sheet of the present invention.

Fig. 2 is three-dimensional model diagram.

Fig. 3 a is normalization fluorescence intensity profile on the detector that obtains in standard fluorescence monte carlo modelling.

Fig. 3 b is normalization fluorescence intensity profile on the detector that obtains in the perturbation fluorescence monte carlo modelling based on voxel.

Fig. 4 is standard and based on normalization fluorescence intensity level line distribution plan on the detector obtained in the perturbation fluorescence monte carlo modelling of voxel.

Embodiment

The invention will be further described by reference to the accompanying drawings.

As shown in Figure 1, implementation step of the present invention is as follows:

(1) target biological tissue is carried out three-dimensional dividing, build a three-dimensional voxel model, set up a 3-dimensional digital matrix, voxel one_to_one corresponding in each element and three-dimensional voxel model in 3-dimensional digital matrix, a kind of biological tissue of numerical identity of each element, arranges the optical property parameter absorption coefficient of biological tissue, scattering coefficient, fluorophore absorption coefficient, refractive index and anisotropy factor;

(2) incident light source is characterized by the set of setting number photon, determines the initial position of light source and the position of direction and detector;

(3) exciting light photon is thrown in, using incident light source position and incident light direction as the initial position of each photon and direction, follow the trail of exciting light photon to transmit in biological tissues, calculate exciting light photon and be converted into fluorescent photon ratio, fluorescent photon weight received on the calculating detector of exciting light path, preserves the routing information of photon;

Following the trail of photon transmitting procedure concrete steps is in biological tissues:

(3.1) exciting light photon is thrown in, exciting light photon initial weight is 1, if exciting light photon is outside biological tissue, then by iterative algorithm along the exciting light direction of propagation, photon is moved on tissue surface, if exciting light photon in biological tissues, then this position is set as photo emissions position;

(3.2) scattering length, the distance between double scattering event, is sampled by the scattering coefficient of current location, arranges step-length Sleft=-In ξ/μ that photon bag often walks _s;

(3.3) photon moves along scattering path and moves a step; If move after moving a step, photon enters into other voxel, and so photon moves to the termination of voxel interphase place, determines photon rum point and the new direction of propagation on voxel interphase.After entering next voxel, photon continues to have moved residue step-length;

(3.4) along exciting light path, utilize Bill-youth Bai Dingli, calculate exciting light photon weight; Often make a move, if exciting light photon does not also enter fluorescence area, the attenuation ratio of photon weight is if exciting light photon enters fluorescence area, the attenuation ratio of photon weight is $\exp \left((-{\mathrm{\μ}}_a^{\mathrm{ex}}+{\mathrm{\μ}}_{\mathrm{af}})l\right);$

(3.5) if exciting light photon fluorescence area is in biological tissues absorbed and excites generation fluorescence.The ratio that exciting light photon is converted into fluorescent photon is η (1-exp (-μ _afl));

(3.6) after fluorescent photon is excited, along exciting light path, utilize Bill-youth Bai Dingli, calculate fluorescent photon weight; Often make a move, at fluorescence area fluorescent photon weight attenuation ratio be in non-fluorescence region, the attenuation ratio of fluorescent photon weight is

(3.7) in scattering position, new scattering direction vector is calculated according to Henyey-Greenstein function;

(3.8) step (3.3) is repeated, until all photons are dead or effusion medium;

(3.9) source position is at r _splace, excite position at r place, detector is at r _dplace receives all fluorescent photon weights and W (r _s, r _d, r) be:

$\begin{array}{c}W(r_s,r_d,r)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\exp (-\underset{j=1}{\overset{p_i}{\mathrm{\Σ}}}({{\mathrm{\μ}}_a}^{\mathrm{ex}}\left(r_j\right)+{\mathrm{\μ}}_{\mathrm{af}}\left(r_j\right))l_i\left(r_j\right))\×(1-\exp \left({-\mathrm{\μ}}_{\mathrm{af}}\left(r\right)l_i\left(r\right)\right)\×\\ \exp (-\underset{j=p_i+1}{\overset{q_i}{\mathrm{\Σ}}}({{\mathrm{\μ}}_a}^{\mathrm{em}}\left(r_j\right)+{\mathrm{\μ}}_{\mathrm{af}}\left(r_j\right))l_i\left(r_j\right))\mathrm{\η}\end{array}$

At non-fluorescence region μ _af(r _j) be 0;

(3.10) after photon is followed the trail of and is terminated, record each photon the path that photon is walked in the voxel index value of process and this voxel.

(4) when fluorescence index variation in biological tissue, the photon path information of preserving is utilized, fluorescent photon weight received on quick calculating detector;

(4.1) determine the voxel index value of fluorescence index variation in biological tissue, and extract the path that in this voxel, photon is walked.

(4.2) the photon path information will extracted, substitutes into equation

$\begin{array}{c}W(r_s,r_d,r)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\exp (-\underset{j=1}{\overset{p_i}{\mathrm{\Σ}}}({{\mathrm{\μ}}_a}^{\mathrm{ex}}\left(r_j\right)+{\mathrm{\μ}}_{\mathrm{af}}\left(r_j\right))l_i\left(r_j\right))\×(1-\exp \left({-\mathrm{\μ}}_{\mathrm{af}}\left(r\right)l_i\left(r\right)\right)\×\\ \exp (-\underset{j=p_i+1}{\overset{q_i}{\mathrm{\Σ}}}({{\mathrm{\μ}}_a}^{\mathrm{em}}\left(r_j\right)+{\mathrm{\μ}}_{\mathrm{af}}\left(r_j\right))l_i\left(r_j\right))\mathrm{\η}\end{array}$

Thus the fluorescent photon weight on direct calculating detector;

(5) after following the trail of all photons, the fluorescent photon weight that output detector receives.

The present invention is set forth further below by example.

Embodiment:

With the model of three-dimensional voxel shown in Fig. 2.This model contains 4 kinds of dissimilar tissues, and be muscle, bone, kidney, heart respectively, fluorophore is positioned at kidney portion; Fluorophore absorption coefficient is set to 1cm ^-1, quantum efficiency is 1; Altogether containing 301401 voxels in voxel model, the size of each voxel is 0.05cm, and media size is 3.05cm × 3.05cm × 4.05cm.The position in source is taken at Fig. 2 marker location, and detector is taken in the region on 180 degree, opposite, source, and it is evenly distributed on 80 layers, every layer of 180 detector; We list the optical parameter value of each tissue in Table 1, and all the other parameters such as g is set to 0.9, and refractive index is set to 1.37, are the representative value of biological tissue under near infrared spectrum.The photon number of simulation is 10 ⁹.

The optical parameter value that table 1 is respectively organized

Fig. 3 a is normalization fluorescence intensity profile on the detector that obtains in the simulation of standard fluorescence Monte Carlo (sfMC) method, Fig. 3 b is normalization fluorescence intensity profile on the detector that obtains in perturbation fluorescence Monte Carlo (pfMC) the method simulation based on voxel, and Fig. 4 is that on the detector that obtains in standard and the monte carlo modelling of perturbation fluorescence, the level line of normalization fluorescence intensity profile picture compares.From detector, normalization fluorescence intensity profile and contour map can find out that two kinds of methods and resultses meet very well.

Claims (6)

1., based on a perturbation fluorescence Monte-Carlo Simulation Method for voxel, it is characterized in that comprising the following steps:

(1) target biological tissue is determined, three-dimensional dividing is carried out to target biological tissue, build a three-dimensional voxel model, set up a 3-dimensional digital matrix, voxel one_to_one corresponding in each element and three-dimensional voxel model in 3-dimensional digital matrix, a kind of biological tissue of numerical identity of each element, arranges the optical property parameter of biological tissue: absorption coefficient, scattering coefficient, fluorophore absorption coefficient, refractive index and anisotropy factor;

(2) incident light source is characterized by the set of setting number photon, determine the initial position of light source and the position of direction and detector, throw in exciting light photon, follow the trail of exciting light photon to transmit in biological tissues, calculate exciting light photon and be converted into fluorescent photon ratio, fluorescent photon weight received on the calculating detector of exciting light path, preserves the routing information of photon;

(3) when the fluorophore absorption coefficient in biological tissue changes, the photon path information of preserving is utilized, fluorescent photon weight received on calculating detector.

2. the perturbation fluorescence Monte-Carlo Simulation Method based on voxel according to claim 1, it is characterized in that: in step (2), exciting light photon is converted in fluorescent photon process, assuming that fluorescence is anisotropic scattering, fluorescent scattering direction is consistent with exciting light scattering direction.

3. the perturbation fluorescence Monte-Carlo Simulation Method based on voxel according to claim 1, is characterized in that: assuming that the path of exciting light photon and fluorescent photon random walk is in biological tissues identical in step (2).

4. the perturbation fluorescence Monte-Carlo Simulation Method based on voxel according to claim 1, is characterized in that: assuming that exciting light goes out fluorescence in fluorescence area continuous agitation in step (2).

5. the perturbation fluorescence Monte-Carlo Simulation Method based on voxel according to claim 1, is characterized in that: step (2) is specifically carried out according to the following steps:

(2.1) follow the trail of the transmission in biological tissues of exciting light photon, along exciting light path, utilize formulae discovery exciting light photon weight w (r _s, r);

$w(r_s,r)=w_0\exp (-\underset{j=1}{\overset{p_i}{\mathrm{\Σ}}}({\mathrm{\μ}}_a^{\mathrm{ex}}\left(r_j\right)+{\mathrm{\μ}}_{\mathrm{af}}\left(r_j\right))l\left(r_j\right))$

In formula: r _sfor light source position, r are fluorescence excitation position;

for absorption coefficient, the μ of exciting light _affor fluorophore absorption coefficient;

R _j(j=1 ... p _i) for exciting light photon is from r _spath between the jth-1 time and jth time scattering events of r;

W ₀for the initial weight of exciting light photon, l are at r _jthe path that microcell experiences;

(2.2) if excitation photon fluorescence area is in biological tissues absorbed and excites generation fluorescence, calculate according to fluorophore coefficient and quantum efficiency the Probability p (r) that exciting light is converted into fluorescence;

p(r)＝η(1-exp(-μ _afl(r)))

In formula: η is quantum efficiency;

(2.3) along exciting light path, formulae discovery fluorescent photon weight w (r, r is utilized _d);

$w(r,r_d)={w_0}^{\′}\exp (-\underset{j=p_i+1}{\overset{q_i}{\mathrm{\Σ}}}({\mathrm{\μ}}_a^{\mathrm{em}}\left(r_j\right)+{\mathrm{\μ}}_{\mathrm{af}}\left(r_j\right))l\left(r_j\right))$

In formula: r _dfor detector position;

for the absorption coefficient of fluorescence;

R _j(j=p _i+ 1 ... .q _i) for fluorescent photon is from r to r _djth-1 time and jth time scattering events between path;

W ₀' for the initial weight of fluorescent photon, l be at r _jthe path that microcell experiences;

(2.4) follow the trail of all fluorescent photons, until fluorescent photon effusion tissue or dead, the fluorescent photon weight on calculating detector, preserves the routing information from source to detector photon.Source position is at r _splace, excite position at r place, detector is at r _dplace receives all fluorescent photon weights and W (r _s, r _d, r) be:

$\begin{array}{c}W(r_s,r_d,r)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\exp (-\underset{j=1}{\overset{p_i}{\mathrm{\Σ}}}({{\mathrm{\μ}}_a}^{\mathrm{ex}}\left(r_j\right)+{\mathrm{\μ}}_{\mathrm{af}}\left(r_j\right))l_i\left(r_j\right))\×(1-\exp (-{\mathrm{\μ}}_{\mathrm{af}}\left(r\right)l_i\left(r\right))\×\\ \exp (-\underset{j=p_i+1}{\overset{q_i}{\mathrm{\Σ}}}({{\mathrm{\μ}}_a}^{\mathrm{em}}\left(r_j\right)+{\mathrm{\μ}}_{\mathrm{af}}\left(r_j\right))l_i\left(r_j\right))\mathrm{\η}\end{array}$

At non-fluorescence region μ _af(r _j) be 0;

6. the perturbation fluorescence Monte-Carlo Simulation Method based on voxel according to claim 1, is characterized in that: step (3) is specifically carried out according to the following steps:

(3.1) determine the index value of the voxel of fluorophore absorption coefficient change in biological tissue, and extract the path that in this voxel, fluorescent photon is walked;

(3.2) utilize the fluorescent photon routing information preserved, substitute into the fluorescent photon weight on the direct calculating detector of following formula,

$\begin{array}{c}W(r_s,r_d,r)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\exp (-\underset{j=1}{\overset{p_i}{\mathrm{\Σ}}}({{\mathrm{\μ}}_a}^{\mathrm{ex}}\left(r_j\right)+{\mathrm{\μ}}_{\mathrm{af}}\left(r_j\right))l_i\left(r_j\right))\×(1-\exp (-{\mathrm{\μ}}_{\mathrm{af}}\left(r\right)l_i\left(r\right))\×\\ \exp (-\underset{j=p_i+1}{\overset{q_i}{\mathrm{\Σ}}}({{\mathrm{\μ}}_a}^{\mathrm{em}}\left(r_j\right)+{\mathrm{\μ}}_{\mathrm{af}}\left(r_j\right))l_i\left(r_j\right))\mathrm{\η}\end{array}.$