CN104706320B - Fluorescent diffusion optical cross-sectional image reestablishing method based on dfMC model - Google Patents

Abstract

The invention discloses a fluorescent diffusion optical cross-sectional image reestablishing method based on a decoupling fluorescence Monte Carlo (dfMC) model, and belongs to the technical field of biomedical engineering. The method includes the steps of firstly, determining a detecting area, selecting a plurality of scanning points in the detecting area, and obtaining fluorescent intensity distribution on a detector; secondly, establishing a three-dimensional digital model for depicting a tissue optical parameter space structure, conducting forward-direction white Monte Carlo simulation of stimulating photons according to the scanning positions and directions of a light source, tracking the stimulating photons, and recording the corresponding physical quantities of the photons reaching the detector on a path; thirdly, calculating the weight of fluorescent photons through a dfMC method, and calculating a fluorescent Jacobi matrix; fourthly, calculating the positions and absorbing coefficients of fluorophores in tissue through iterative reconstruction of GPU clusters. The method has the advantage of providing an accurate and rapid reestablishing method for a three-dimensional fluorescence tomography system through the high-precision dfMC model on the basis of the accelerated iterative reconstruction process of the GPU clusters.

Description

A kind of fluorescence Diffuse Optical Tomography image rebuilding method based on dfMC model

Technical field

The invention belongs to biomedical engineering technology field, more particularly to one kind are based on decoupling fluorescence Monte Carlo The fluorescence Diffuse Optical Tomography image rebuilding method of dfMC model.

Background technology

In recent years, fluorescence Diffuse Optical Tomography imaging technique has developed into for a kind of important imaging tool^[1], and It is widely used in cancer diagnosis, medicament research and development and the visual research of gene expression.The accurate reconstruction of fluorescent target thing will Depend on and set up an accurate photon transport model (i.e. so-called forward problem) and the inverse model of photon propagation model is (inverse To problem).Forward problem is mainly based upon the diffusion approximation model under the Boltzmann radiation transmission framework of preservation of energy^[2]With Monte Carlo (MC) model^[3].Because diffusion approximation model is low for image reconstruction times cost, it is widely adopted；But it is existing Some diffusion models are suitable only for high scattering material, in absorptance and scattering coefficient quite, or source visit distance with average When transmission free path is suitable, model is just no longer applicable.Monte-Carlo model adapts to the biological tissue of any optical parametric, and appoints The complicated geometry object of meaning, but because calculation cost is very huge, computational efficiency is low." gold generally as forward model Standard "^[3].

The existing forward direction fluorescence MC model being suitable for Reverse Problem solution mainly has adjoint method model (adjoint MC )^[4]With perturbation method model (perturbation MC)^[5].Adjoint method supposes that source and detector are tradable, that is, it is right to have Etc. property, this equity is to set up for optical fiber type detection system, but the situation detecting for free space, light source and The complex characteristics of CCD make this equity hardly result in guarantee.Simultaneously for the reconstruction effect detecting situation based on free space Rate is low-down, and the detector of substantial amounts of CCD pixel size makes Monte Carlo simulation number of times increase severely.And perturbation method master If by fluorescence intensity is received to obtain on in-house MC path disturbance is come calculating detector to exciting light.Perturbation method it not Limited by imaging system, be there is the advantage of speed and precision.Become the focus of research.But perturbation method ignores fluorogen pair The impact of exciting light, have ignored the difference of exciting light and fluorescence optical parametric in tissue, and this many hypothesis makes perturbation method It is that one kind has inclined MC simulation, its fluorescence statistic has deviation with actual true fluorescent value, and this deviation can be with optics The change of parameter and number of photons and change, in most cases, error can be very big, and because photon path information is huge, in body Will make under prime model situation to rebuild inefficiency, whole process takes very big.

List of references

[1] V.Ntziachristos, " Going deeper than microscopy：the optical imaging Frontier in biology, " Nat.Meth.7,603-614 (2010).

[2] X.Song, D.Wang, N.Chen, J.Bai, and H.Wang, " Reconstruction for free- space fluorescence tomography using a novel hybrid adaptive finite element Algorithm, " Opt.Express 15,18300-18317 (2007).

[3] J.Swartling, A.Pifferi, A.M.K.Enejder, and S.Andersson-Engels, “Accelerated Monte Carlo models to simulate fluorescence spectra from layered Tissues, " J.Opt.Soc.Am.A 20,714-727 (2003).

[4] G.Ma, J.F.Delorme, P.Gallant, and D.A.Boas, " Comparison of simplified Monte Carlo simulation and diffusion approximation for the fluorescence Signal from phantoms with typical mouse tissue optical properties, " Appl.Opt.46,1686-1692 (2007).

[5] Sassaroli, " Fast perturbation Monte Carlo method for photon Migration in heterogeneous turbid media, " Opt.Lett.36,2095-2097 (2011).

Content of the invention

In order to overcome the above-mentioned deficiencies of the prior art, the invention provides a kind of be based on decoupling fluorescence Monte Carlo (dfMC) the fluorescence Diffuse Optical Tomography image rebuilding method of model, the purpose of the present invention is to keep MC simulation " goldstandard " Precision, improves the time consuming nature of MC simulation simultaneously, realizes high accuracy rapid fluorescence tomographic image reconstructing.

The technical solution adopted in the present invention is：

A kind of fluorescence Diffuse Optical Tomography image rebuilding method based on dfMC model is it is characterised in that include following step Suddenly：

(1) determine search coverage, select light source scanning point, and the fluorescence obtaining under different scanning position on detector overflows Penetrate light distribution；

(2) digitized object to be detected, sets up voxel model, carries out white Monte Carlo according to scanning light source position and direction Simulation, records the weight of the excitation photon collected by detectorAnd the routing information in each voxel and corresponding thing Reason amount；

(3) according to dfMC model, calculate fluorescent photon weight；

(4) fluorescence Jacobian matrix fluorescence Jacobi is set up by excitation photon routing information and corresponding physical quantity；

(5) weight of the fluorescent photon according to output and Fluorescence Fluorescence Jacobian matrix, iterative approximation calculates in tissue Fluorogen position and its absorptance.

In step (2), entirely white Monte Carlo process is accelerated using GPU cluster, record reaches the excitation photon of detector Related physical quantity on path, the path related physical quantity L in voxel v_1v、L_2vIt is respectively：

In formula：N is the scattering imaging in voxel v,It is the unit scattering direction vector of i ＆ lt scattering in voxel v, P_IIt is Isotropism directivity function, P_AIt is anisotropic orientation function, l_jFor the jth in voxel v time scattering path length.

In step (2), GPU cluster accelerates entirely white Monte Carlo process specific as follows：

2.1 on cluster, according to the number of times that MC simulation is carried out on number of light sources and each node of clustered node number assignment, And corresponding light source position and directional information are transferred on node；

2.2 on each node, according to photon number storage allocation to be simulated, can be assigned to tissue model information Texture depositor in GPU, optical parametric information is assigned to literal register to improve reading speed, and the technology according to GPU Information determines the number of parallel block and thread in optimum CUDA operation；

2.3 on each thread, simulates a light attached bag, and follows the tracks of the mistake in scatter within tissue and absorption after photon is launched Journey, records introductory path information in each voxel for the photon and physical quantity, and photon terminates photonic analogy after overflowing tissue.

Step (3), calculates fluorescence weight w^emP () computing formula is as follows：

In formula：P is the six-vector in position and unit direction,For unit direction vector, k is photon on in-house road The voxel number that footpath is passed through, m is photon in in-house scattering imaging, μ_afvIt is fluorescence coefficient in voxel v,It is voxel v Interior excites light scattering coefficient, and η is quantum yield, w^exFor excitation photon weight, g_v ^emFor fluorescence in voxel v anisotropy because Son, g_v ^exFor the anisotropy factor of exciting light in voxel v, μ_av ^emFor the absorptance of fluorescence in voxel v, μ_av ^exFor in voxel v The absorptance of exciting light, and

In step (4), the computing formula of fluorescence Jacobian matrix is as follows：

In formula：P is the six-vector in position and unit direction,For unit direction vector, m is that photon dissipates in-house Penetrate number of times, μ_afvIt is fluorescence coefficient in voxel v,It is the light scattering coefficient that excites in voxel v, η is quantum yield, w^exFor exciting Photon weight, g_v ^emFor the anisotropy factor of fluorescence in voxel v, g_v ^exFor the anisotropy factor of exciting light in voxel v, μ_av ^em For the absorptance of fluorescence in voxel v, μ_av ^exFor the absorptance of exciting light in voxel v, andJ is the light source numbering in node i, and d numbers for detector, and v is the model body needing to rebuild Element numbering, for photon in in-house scattering imaging, J is fluorescence Jacobian matrix to m,For corresponding voxel v, light source j, detection The fluorescence Jacobean matrix array element of device d.

The whole fluorescence process of reconstruction of step (5) is all accelerated by GPU cluster, and on each node, calculated fluorescence is refined can Ratio matrix-block, computing formula is as follows：

D_i=[D₁₁…D_1t；…；D_s1…D_st]^T

In formula：S is the number of light sources in node i, and t is detector sum, and n is the model number of voxels mesh needing to rebuild, D_i The difference of the fluorescence distribution obtaining with dfMC simulation for the fluorescence intensity distribution that the detector under light sources all in node i detects Value, D_stThe fluorescence intensity distribution detecting for t-th probe unit in the detector under s light source in node i with corresponding The difference of the fluorescence distribution that dfMC simulation obtains.

Step (5), intermediate quantity J_iJ_i ^TX, J_iD_iIt is transferred to and carry out conjugate gradient iterative fluorescence distribution again on main section, These intermediate quantities are fewer than the memory space that fluorescence Jacobian matrix block takies, and iterative calculation formula is as follows：

In formula：I is node serial number, and N is total interstitial content, and X is fluorogen absorptance distribution matrix J_i ^TFor matrix J_i's Transposed matrix, f (X) is object function.

Compared with prior art, beneficial effects of the present invention are：It is capable of the three-dimensional of fluorogen in biological tissue is divided Cloth and the reconstruction of fluorescence coefficient, have reconstruction precision height, fireballing advantage.

Brief description

Fig. 1 is white Monte Carlo simulation flow chart under the parallel acceleration of GPU of the present invention；

Fig. 2 is that the lower dfMC three-dimensional fluorescence of the parallel acceleration of GPU cluster of the present invention rebuilds flow chart；

Fig. 3 a, Fig. 3 b are respectively the present invention and rebuild body mould experimental result picture.

Specific embodiment

The present invention is further described below in conjunction with the accompanying drawings.

The implementation steps of the present invention are as follows：

1st, determine search coverage, select light source scanning point, and the fluorescence obtaining under different scanning position on detector overflows Reflection intensity distribution D, so that fluorescence inverting solves；

2nd, digitized detects tissue, foundation group organization model, and foundation scanning light source position and direction carry out GPU and accelerate lower white illiteracy Special Monte Carlo Simulation of Ions Inside, obtains routing information in each voxel for the excitation photon collected by detector and related physical quantity L_1vWith L_2v, the white Monte Carlo process under GPU cluster acceleration is as shown in figure 1, step is as follows：

2.1 on cluster, according to the number of times that MC simulation is carried out on number of light sources and each node of clustered node number assignment, And corresponding light source position and directional information are transferred on node；

2.2 on each node, according to photon number storage allocation to be simulated, can be assigned to tissue model information Texture depositor in GPU, optical parametric information is assigned to literal register to improve reading speed, and the technology according to GPU Information determines the number of parallel block and thread in optimum CUDA operation；

2.3 on each thread, simulates a light attached bag, and follows the tracks of the mistake in scatter within tissue and absorption after photon is launched Journey, records introductory path information in each voxel for the photon and physical quantity, and photon terminates photonic analogy after overflowing tissue.

3rd, carry out under GPU cluster, based on the fluorescence intensity that detects under each light source be distributed that D and MC simulation obtains a large amount of Physical message along photon path record carries out fluorescence distribution inverting solution, and concrete steps are as shown in Figure 2：

3.1 on cluster, carries out the cluster in process of reconstruction according to cluster distribution information during white Monte Carlo simulation and divides Join so that the corresponding physical quantity of the excitation photon of some light source along paths record is corresponded on each node, if altogether There is n light source, then the information with path-dependent of 1/n light source MC simulation is had on each node；

3.2 on each node, is loaded into the along-path information of corresponding excitation photon, and calculates phase according to dfMC method The Fluorescence Fluorescence Jacobian matrix block J of corresponding node i_i, computing formula is as follows：

In formula：P is the six-vector in position and unit direction,For unit direction vector, m is that photon dissipates in-house Penetrate number of times, μ_afvIt is fluorescence coefficient in voxel v,It is the light scattering coefficient that excites in voxel v, η is quantum yield, w^exFor exciting Photon weight, g_v ^emFor the anisotropy factor of fluorescence in voxel v, g_v ^exFor the anisotropy factor of exciting light in voxel v, μ_av ^em For the absorptance of fluorescence in voxel v, μ_av ^exFor the absorptance of exciting light in voxel v, andJ is the light source numbering in node i, and d numbers for detector, and v is the model body needing to rebuild Element numbering, for photon in in-house scattering imaging, J is fluorescence Jacobian matrix to m,For corresponding voxel v, light source j, detection The fluorescence Jacobean matrix array element of device d.

3.3 can obtain calculated fluorescence Jacobian matrix on node according to the light source number in node and detector number Block, its computing formula is as follows：

D_i=[D₁₁…D_1t；…；D_s1…D_st]^T

In formula：S is the number of light sources in node i, and t is detector sum, and n is the model number of voxels mesh needing to rebuild, D_i The difference of the fluorescence distribution obtaining with pfMC simulation for the fluorescence intensity distribution that the detector under light sources all in node i detects Value, D_stThe fluorescence intensity distribution detecting for t-th probe unit in the detector under s light source in node i with corresponding The difference of the fluorescence distribution that dfMC simulation obtains.

3.4 according to the weight of fluorescent photons of output and fluorescence Jacobian matrix block, iterate to calculate out light blob position and its Absorptance, its computing formula is as follows：

In formula：I is node serial number, and N is total interstitial content, and X is fluorogen absorptance distribution matrix J_i ^TFor matrix J_i's Transposed matrix, f (X) is object function.

Tested with a simple fluorophor mould below and provided in an embodiment of the present invention covered based on decoupling fluorescence to verify The feasibility of the fluorescence Diffuse Optical Tomography imaging method of special Caro model, described below：

Make one and filled with and be inserted with two in fat emulsion solution glass (internal diameter 30mm, high 40mm) and fill DIA-BOA Carefully excellent (internal diameter 2mm) the body mould of fluorescent dye glass, by Lamda 950 measure Intralipid solution optical parametric and The fluorescence coefficient of glo-stick, fatty aqueous solution μ_s=214.7cm^-1And μ_a=0.01cm^-1, fluorescence system in two glass glo-sticks Number is respectively μ_af=0.6,0.8cm^-1.With dfMC method, fluorescence reconstruction is carried out to this model, result as shown in Figure 3 a, 3 b, The glo-stick position rebuilding in Fig. 3 b is consistent with glo-stick position true in Fig. 3 a.The experiment of body mould shows the method high precision, with When pass through GPU cluster energy accelerated reconstruction.

Claims (7)

1. a kind of fluorescence Diffuse Optical Tomography image rebuilding method based on dfMC model is it is characterised in that comprise the following steps：

(1) determine search coverage, select light source scanning point, and obtain the fluorescence diffused light on detector under different scanning position Strong distribution；

(2) digitized object to be detected, sets up voxel model, carries out white Monte Carlo mould according to scanning light source position and direction Intend, record the weight of the excitation photon collected by detectorAnd the routing information in each voxel and corresponding physics Amount, ex represents exciting light, and for photon in in-house scattering imaging, p represents the state of photon to m；

(3) according to dfMC model, calculate fluorescent photon weight；

(4) fluorescence Jacobian matrix is set up by excitation photon routing information and corresponding physical quantity；

(5) weight of the fluorescent photon according to output and fluorescence Jacobian matrix, iterative approximation calculates in-house fluorogen Position and its absorptance.

2. the fluorescence Diffuse Optical Tomography image rebuilding method based on dfMC model according to claim 1, its feature exists In step (2), entirely white Monte Carlo simulation process is accelerated using GPU cluster, record reaches the excitation photon road of detector Related physical quantity on footpath, the path related physical quantity L in voxel v_1v、L_2vIt is respectively：

$L_{1v}=\underset{i=1}{\overset{n}{\Σ}}{P}_I({\hat{s}}_i\·{\hat{s}}_{i+1})/P_A({\hat{s}}_i\·{\hat{s}}_{i+1})$

$L_{2v}=\underset{i=1}{\overset{n}{\Σ}}\underset{j=1}{\overset{i}{\Σ}}{l}_j{P}_I({\hat{s}}_i\·{\hat{s}}_{i+1})/P_A({\hat{s}}_i\·{\hat{s}}_{i+1})$

In formula：N is the scattering imaging in voxel v,It is the unit scattering direction vector of i ＆ lt scattering in voxel v, P_IIt is each to together Property directivity function, P_AIt is anisotropic orientation function, l_jJth time scattering path length for photon in voxel v.

3. the fluorescence Diffuse Optical Tomography image rebuilding method based on dfMC model according to claim 1, its feature exists In：In step (2), GPU cluster accelerates entirely white Monte Carlo simulation process specific as follows：

2.1 on cluster, according to the number of times that MC simulation is carried out on number of light sources and each node of clustered node number assignment, and handle Corresponding light source position and directional information are transferred on node；

2.2 on each node, according to photon number storage allocation to be simulated, tissue model information can be assigned to GPU Interior texture depositor, optical parametric information is assigned to literal register to improve reading speed, and the technical information according to GPU Determine the number of parallel block and thread in optimum CUDA operation；

2.3 on each thread, simulates a light attached bag, and follows the tracks of the process in scatter within tissue and absorption after photon is launched, Record introductory path information in each voxel for the photon and physical quantity, photon terminates photonic analogy after overflowing tissue.

4. the fluorescence Diffuse Optical Tomography image rebuilding method based on dfMC model according to claim 1, its feature exists In step (3), calculate fluorescent photon weight w^emP () computing formula is as follows：

$w_m^{em}\left(p\right)=w_m^{ex}\left(p\right)\underset{v=1}{\overset{k}{\Σ}}\frac{{\mathrm{\η\μ}}_{afv}}{{\mathrm{\μ}}_{sv}^{ex}}(L_{1v}-{\mathrm{\μ}}_{afv}{L}_{2v})\underset{j=i}{\overset{m}{\Π}}\left(\frac{P_A({\hat{s}}_j\·{\hat{s}}_{j+1},g_v^{em})}{P_A({\hat{s}}_j\·{\hat{s}}_{j+1},g_v^{ex})}\right)\mathrm{\Γ}(p_i\→p_e,{\mathrm{\μ}}_{av}^{em}-{\mathrm{\μ}}_{av}^{ex})$

In formula：P is the six-vector in position and unit direction,For unit direction vector, k is for photon in in-house path institute Through voxel number, m be photon in in-house scattering imaging, μ_afvIt is fluorescence coefficient in voxel v,It is swashing in voxel v Luminous scattering coefficient, η is quantum yield, w^exFor exciting light photon weight, g_v ^emFor the anisotropy factor of fluorescence in voxel v, g_v ^exFor the anisotropy factor of exciting light in voxel v, μ_av ^emFor the absorptance of fluorescence in voxel v, μ_av ^exExcite in voxel v The absorptance of light, andEm is expressed as fluorescence, L_1vAnd L_2vRelated for the path in voxel v Physical quantity, P_AIt is anisotropic orientation function, p_i→p_eIt is expressed as from photon states p_iTo photon states p_e,Represent photon shape State p_iTo photon states p_eBetween distance.

5. the fluorescence Diffuse Optical Tomography image rebuilding method based on dfMC model according to claim 1, its feature exists In step (4), the computing formula of fluorescence Jacobian matrix is as follows：

$\begin{array}{c}{J}_{jd}^v=\underset{m=0}{\overset{\mathrm{\∞}}{\Σ}}\frac{\∂w_m^{em}\left(p\right){|}_{\det ectord}}{\∂{\mathrm{\μ}}_{afv}}\\ =\underset{m=0}{\overset{\mathrm{\∞}}{\Σ}}{w}_m^{em}\left(p\right){|}_{\det ectord}\left(\frac{\mathrm{\η}}{{\mathrm{\μ}}_{sv}^{ex}}{I}_{Iv}-\frac{{\mathrm{\η\μ}}_{afv}}{{\mathrm{\μ}}_{sv}^{ex}}{L}_{2v}\right)\underset{j=i}{\overset{m}{\Π}}\left(\frac{P_A({\hat{s}}_j\·{\hat{s}}_{j+1},g_v^{em})}{P_A({\hat{s}}_j\·{\hat{s}}_{j+1},g_v^{ex})}\right)\mathrm{\Γ}(p_i\→p_e,{\mathrm{\μ}}_{av}^{em}-{\mathrm{\μ}}_{av}^{ex})\end{array}$

In formula：P is the six-vector in position and unit direction,For unit direction vector, m is photon in in-house scattering time Number, μ_afvIt is fluorescence coefficient in voxel v,It is the light scattering coefficient that excites in voxel v, η is quantum yield, w^exFor excitation light Sub- weight, g_v ^emFor the anisotropy factor of fluorescence in voxel v, g_v ^exFor the anisotropy factor of exciting light in voxel v, μ_av ^emFor The absorptance of fluorescence, μ in voxel v_av ^exFor the absorptance of exciting light in voxel v, andj For the light source numbering in node i, d is detector numbering, and v is the model voxel numbering needing to rebuild, and m is photon in-house Scattering imaging, J is fluorescence Jacobian matrix,For corresponding to the fluorescence Jacobean matrix array element of voxel v, light source j, detector d, Detector d represents detector d, L_1vAnd L_2vFor the path related physical quantity in voxel v, P_AIt is anisotropic orientation function, p_i →p_eIt is expressed as from photon states p_iTo photon states p_e,Represent photon states p_iTo photon states p_eBetween distance.

6. the fluorescence Diffuse Optical Tomography image rebuilding method based on dfMC model according to claim 1, its feature exists All accelerated by GPU cluster in the whole fluorescence process of reconstruction of step (5), calculated fluorescence Jacobian matrix on each node Block, computing formula is as follows：

D_i=[D₁₁…D_1t；…；D_s1…D_st]^T

In formula：S is the number of light sources in node i, and t is detector sum, and n is the model number of voxels mesh needing to rebuild, D_stFor section The fluorescence intensity distribution that in detector under the upper s light source of point i, t-th probe unit detects is simulated with corresponding dfMC The difference of the fluorescence distribution arriving.

7. the fluorescence Diffuse Optical Tomography image rebuilding method based on dfMC model according to claim 6, its feature exists In step (5), intermediate quantity J_iJ_i ^TX, J_iD_iIt is transferred to and carry out conjugate gradient iterative fluorescence distribution again, in these on main section The area of a room is fewer than the memory space that fluorescence Jacobian matrix block takies, and iterative calculation formula is as follows：

$\underset{X}{\min }f\left(X\right)=\left|\right|J^{T}JX-J^{T}D\left|\right|=\left|\right|\underset{i=1}{\overset{N}{\Σ}}\left(J_i{J}_i^{T}X-J_i{D}_i\right)\left|\right|$

D_i=[D₁₁…D_1t；…；D_s1…D_st]^T

In formula：I is node serial number, and N is total interstitial content, and X is fluorogen absorptance distribution matrix Ji T is turning of matrix J i Put matrix, f (X) is object function, Di is the fluorescence intensity distribution and dfMC that the detector under all light sources in node i detects Simulate the difference of fluorescence distribution obtaining, Dst is that in detector under s light source in node i, t-th probe unit detects The difference of the fluorescence distribution that fluorescence intensity distribution is obtained with corresponding dfMC simulation.