CN104331641B - A kind of fluorescence Monte-Carlo Simulation Method accelerated based on concentrating type GPU - Google Patents
A kind of fluorescence Monte-Carlo Simulation Method accelerated based on concentrating type GPU Download PDFInfo
- Publication number
- CN104331641B CN104331641B CN201410534214.7A CN201410534214A CN104331641B CN 104331641 B CN104331641 B CN 104331641B CN 201410534214 A CN201410534214 A CN 201410534214A CN 104331641 B CN104331641 B CN 104331641B
- Authority
- CN
- China
- Prior art keywords
- photon
- msub
- mrow
- mtr
- mtd
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Landscapes
- Investigating, Analyzing Materials By Fluorescence Or Luminescence (AREA)
Abstract
The present invention relates to a kind of fluorescence DSMC accelerated based on concentrating type GPU, multiple sources can be simulated simultaneously, greatly save the time that simulated photons are transmitted in biological tissues, consider fluorogen in tissue and, to the absorption of exciting light, the photon transmission of exciting light and fluorescence in biological tissue has been followed the trail of respectively.This method precision is high, can obtain the optical transport information in true biological tissue, these abundant information provide foundation for the optimization of optical imaging system, be that photodynamic diagnosis and light treatment provide accurate tutorial message.
Description
Technical field
The present invention relates to it is a kind of based on concentrating type GPU accelerate fluorescence Monte-Carlo Simulation Method, belong to computer science,
Mathematical simulation and biomedical engineering field.
Background technology
In organism optical field, the target of optical transmission process modeling is the accurate of distribution of photons characteristic in any tissue of development
Quick guantification method, it includes the implication of two levels：One is setting up model as the other application-specific moulds of evaluation
" goldstandard " of type；The second is developing fast and effeciently utility model for the specified conditions in organism optical concrete application.Cover
The essential characteristic of special Caro method is that stochastic problems are emulated, and can effectively solve stochastic problems, or even to many true
The insoluble stochastic problems of qualitative method institute can be solved more easily.Monte Carlo Method regards as gold by area research person
Standard.
Existing Monte Carlo codes only considered the propagation of exciting light in biological tissues, and ignore in biological tissue
The presence of fluorogen.Exciting light can produce fluorescence after being absorbed by fluorogen, can continue to propagate in biological tissues.Due to biological group
The different compositions knitted can send the fluorescence of different characteristic, by following the trail of these fluorescence and obtaining its information, be conducive to biology
The diagnosis and treatment of medical science.Also, with the photon number increase of simulation, the number of simulation source is continuously increased, and institutional framework is got over
Come more complicated, the time of computer simulation, existing code can not meet the demand time-consuming to simulation also at double in increase.
A kind of optimization program structure of invention, the generation of the fluorescence Monte Carlo simulation of parallel high-speed computing is speeded up to using cluster and GPU
Code can meet this demand.
The content of the invention
Present invention aims at there is provided a kind of fluorescence Monte-Carlo Simulation Method accelerated based on concentrating type GPU, the skill
Art can simulate multiple sources simultaneously, greatly save the time that simulated photons are transmitted in biological tissues, and followed the trail of life respectively
The photon transmission of exciting light and fluorescence in thing tissue.
A kind of fluorescence Monte-Carlo Simulation Method accelerated based on concentrating type GPU, it is characterised in that comprise the following steps：
(1) space structure of target biological tissue is divided into a three-dimensional voxel model, sets one and three-dimensional voxel
The numerical value of each element is the organization type of mark in model size identical 3-dimensional digital matrix, matrix；Sets target is biological
The various organization optical property parameter of tissue, optical properties of tissue is the absorption coefficient of exciting light, the scattering system of exciting light
Number, the absorption coefficient of fluorescence, the scattering coefficient of fluorescence, the absorption coefficient of fluorogen, the refraction coefficient of fluorogen and fluorogen
Anisotropy factor；
(2) by this message-passing communication agreements of MPI, host node obtains the quantity of the GPU equipment of each child node, need to
The quantity in the source to be calculated is averagely allocated to the GPU equipment of each child node, and each child node determines grid and block in GPU equipment
Dimension and size, and the subprocess of respective number is opened up according to the number of node GPU equipment, give each by distribution of computation tasks
Subprocess, starts multiple programming and calculating platform CUDA, each child node CPU storage allocations space, video memory space will be calculated
Data are copied on video memory from internal memory；
(3) incident light source is characterized as setting to the set of the photon of number, incident light source position and incident light direction are assigned
Initial position and direction to each photon；
(4) transmitting procedure of each exciting light photon is followed the trail of in target biological tissue, the big of time length will be wherein expended
The step of scale data is parallel, height calculates density is arranged on GPU parallel execution, and remaining procedure is serially performed on CPU；
CPU storage allocations, for depositing GPU output datas, by the data duplication on video memory after calculating to internal memory, each child node is returned
Operation information is to host node；
(5) such as exciting light photon is absorbed by the fluorogen of target biological tissue and produces fluorescence, follows the trail of each fluorescent photon
Transmitting procedure；Held being arranged in the step of the large-scale data for wherein expending time length is parallel, height calculates density on GPU parallel
OK, remaining procedure is serially performed on CPU；CPU storage allocations, for depositing GPU output datas, by video memory after calculating
Data duplication on internal memory, each child node returns to operation information to host node；
(6) the excitation light sub-information and fluorescent photon information that all effusions are exported after all photons, releasing memory are followed the trail of
With video memory space and exit CUDA.
Step (4) specifically includes following steps progress：
(4.1) exciting light photon is delivered；
(4.2) exciting light photon moves a step-length determined by random number, judges that exciting light photon is in the process
It is no to run into extraneous border, if running into border, step (4.3) is gone to, step (4.4) is otherwise gone to；
(4.3) if exciting light photon reflects, the step-length and direction cosines information of exciting light photon are updated, if
Exciting light photon is refracted to outside tissue, then sets photon death, stops following the trail of photon；
(4.4) judge that exciting light photon is to be absorbed or scattered by setting random number, if exciting light photon is absorbed,
Step (4.5) is then gone to, if exciting light photon is scattered, photon direction cosines information is updated, goes to step (4.2)；
(4.5) if exciting light photon, which is absorbed, does not produce fluorescence, setting photon is dead, stops following the trail of photon；
Step (5) specifically includes following steps progress：
(5.1) if exciting light photon is absorbed by fluorogen and produces fluorescence, the inceptive direction of fluorescence is dissipated by isotropism
The deflection angle penetrated and azimuth are determined, continue to follow the trail of the movement of fluorescent photon next step；
(5.2) if fluorescent photon is scattered, photon direction cosines information is updated, step (5.3) is gone to；If fluorescence light
Son is absorbed, then sets photon death, stops following the trail of photon；
(5.3) fluorescent photon moves a step-length determined by random number, judges whether fluorescent photon meets in the process
To extraneous border, if running into border, step (5.4) is gone to, step (5.2) is otherwise gone to；
(5.4) if fluorescent photon reflects, the step-length and direction cosines information of fluorescent photon is updated, step is gone to
(5.2)；If fluorescent photon is refracted to outside tissue, setting photon is dead, stops following the trail of photon.
Based on the output in step (6), statistics calculating and be converted into following required transmission characteristics obtained one are carried out
It is individual or any several：
By in same incident direction, the exciting light photon energy of all effusions is added, you can exciting light effusion energy with
Diffuse the relation curve between angle；
By in same incident direction, the fluorescence light photon energy of all effusions is added, you can obtain fluorescence effusion energy with overflowing
Relation curve between firing angle degree.
The present invention establishes a kind of superfast fluorescence Monte-Carlo Simulation Method, and the exciting light of effusion can be obtained simultaneously
Photon information and fluorescent photon information.These abundant information provide foundation for the optimization of optical imaging system, are photodynamic diagnosis
Accurate tutorial message is provided with light treatment.
Brief description of the drawings
Fig. 1 is basic flow sheet of the invention.
Embodiment
With reference to accompanying drawing, the invention will be further described.
As shown in figure 1, the implementation steps of the present invention are as follows：
(1) space structure of target biological tissue is described as a 3-dimensional digital matrix, i.e. tissue model.In matrix
The voxel of element correspondence target biological tissue, the numerical value of each voxel is the numeral of mark organization type；Voxel is smaller, description group
The tissue model for knitting structure more approaches true mechanics of biological tissue, simulates obtained photon transmission feature accuracy higher；
(2) the various organization characterisitic parameter of sets target biological tissue, optical properties of tissue is the absorption of exciting light
Coefficient, the scattering coefficient of exciting light, the absorption coefficient of fluorescence, the scattering coefficient of fluorescence, the absorption coefficient of fluorogen, fluorogen
The anisotropy factor of refraction coefficient and fluorogen；The precision of these parameter values can also influence simulation precision, can consult corresponding text
Offer to obtain these parameters.Being worth noting is, selected requirement and the target light source consistent wavelength of these parameters；
(3) by this message-passing communication agreements of MPI, host node obtains the quantity of the GPU equipment of each child node, reads in
Data, it would be desirable to which the quantity in the source of calculating is averagely allocated to the GPU equipment of each child node；
(4) each child node calculates the available resources of each GPU equipment of the node, determines the dimension and size of grid and block；
(5) each child node CPU host process opens up the subprocess of respective number according to the number of node GPU equipment, will
Distribution of computation tasks gives each subprocess；
(6) CUDA is started, each child node CPU storage allocations space, video memory space, the data that will be calculated are replicated from internal memory
Onto video memory；
(7) incident light source is characterized as setting to the set of number photon, incident light source position and incident light direction are assigned to
The initial position of each photon and direction；
(8) transmitting procedure of each exciting light photon is followed the trail of；By the large-scale data for wherein expending time length is parallel, high meter
The step of calculating density is arranged on GPU parallel execution, and remaining procedure is serially performed on CPU, and CPU storage allocations are used for
GPU output datas are deposited, by the data duplication on video memory after calculating to internal memory, each child node returns to operation information to main section
Point；
The transmitting procedure for following the trail of each exciting light photon is concretely comprised the following steps：
(8.1) exciting light photon is delivered, if photon is on tissue surface, the position is set as that photon launches position
Put；If photon is in outside tissue surface, photon is shifted into tissue surface automatically, its concrete operations uses alternative manner；Such as
Fruit photon is in organization internal, then will directly launch photon, i.e. photon according to this position and interaction with surface does not occur, and
Directly start transmission of the photon in organization internal；
(8.2) the remaining step-length sleft of photon is calculated:Sleft=-In ξ；
(8.3) calculate photon and move the step-length s moved a step
S=min (- In (ξ)/μt,min(dx,dy,dz))
Here, ξ is the pseudo random number for being evenly distributed on (0,1) produced by computer pseudorandom number generator, dx, dy,
Dz is followed successively by the length of each voxel in tissue model；When photon is the excitation light period of the day from 11 p.m. to 1 a.m：μt=μa+μs+μafp；Work as photon
During for fluorescent photon：μt=μaf+μsf。μaAnd μsIt is the absorption coefficient and scattering coefficient of exciting light respectively；μafAnd μsfIt is glimmering respectively
The absorption coefficient and scattering coefficient of light；μafpIt is the absorption coefficient of fluorogen；
(8.4) light during being somebody's turn to do is judged when the step-length calculated in front direction and (8.3) moves the step of photon one according to photon
Son whether the interface through different tissues；If it is, into step (8.5)；If not, photon is moved into a step-length,
Photon current location is updated, subsequently into step (8.8)；
(8.5) rum point of the photon on interface is determined, method is：The minterm in following expression is found out first：
Here, x, y, z represent photon current location；If Section 1 is minimum, then application point table of the photon on interface
It is up to formula：
Here, x0, y0, z0 are application point of the photon at interface, and other situations are analogized；
(8.6) direction of the photon after (8.5) middle rum point determined interacts with interface, specific method are calculated
For：Interface normal vector (a, b, c) is calculated first, and its expression formula is：
Then, incidence angle θ and refraction angle θ that photon is had an effect with interface are calculated according to the vector expressiont, it is calculated
Formula is：
Cos θ=aux+buy+cuz
cosθt=(1- (1-cos2θ)·n1 2/n2 2)1/2
Here, kx,ky,kzRepresent that photon runs into the direction cosines value before border, n1And n2Expression runs into institute before and after interface
In the refraction coefficient of tissue；
Four kinds of situations calculating photons are finally divided to run into the direction after interface：
If cos θ=0, direction is constant after setting photon runs into border；
If ε ＞ (n are worked as in cos θ=12-n1)2/(n2+n1)2, set photon run into border after direction it is constant, when ε≤
(n2-n1)2/(n2+n1)2, set photon to reflect, its direction is changed into：
If θ ＞ sin-1(n2/n1), set photon to be totally reflected, its direction is changed into：
In other cases, the direction after setting photon is had an effect with interface is：
Wherein, r expression formula is：
(8.7) the step-length s that has moved is updated to have an effect with interface into step (8.6) rum point for photon current location
Distance, then by photon current location update arrive above-mentioned rum point.
(8.8) if photon is refracted in air in (8.6), as photon effusion tissue, setting photon death stops
Photon is followed the trail of and sleft=0 is updated, and records the photon current attribute information, including position, direction, step is gone to
(8.14)；Sleft=Sleft-s* μ in the case of othert.The μ if the photon is exciting light photont=μa+μs+μafp；Work as light
Son is fluorescent photon then μt=μaf+μsf.If Sleft ＞ 0, into step (8.3)；Otherwise step (8.9) is entered.
(8.9) judge that photon is to be absorbed or scattered by setting random number, gone to if the photon is fluorescent photon
Step (8.10), otherwise then goes to step (8.11).
(8.10) as ξ≤μaf/(μaf+μafp) when, then fluorescent photon is absorbed, and setting photon is dead, and stopping is chased after to photon
Track, goes to step (8.14).If ξ ＞ μaf/(μaf+μafp) when, then fluorescent photon is scattered, and goes to step (8.13).
(8.11) as ξ≤(μa+μafp)/(μa+μafp+μs) when, then exciting light photon is absorbed, into step (8.12).When
ξ ＞ (μa+μafp)/(μa+μafp+μs) when, then exciting light photon equilibrium state, goes to step (8.13).
(8.12) judge whether photon is absorbed by fluorogen in biological tissue and produces fluorescence by setting random number, fluorescence
Inceptive direction determined by the deflection angle and azimuth of isotropic scatterning.As ξ≤μafp/(μa+μafp) when, then fluorescence is produced,
Sleft=0 is updated, step (8.2) is gone to.As ξ ＞ μafp/(μa+μafp) when, then photon death is set, stops following the trail of photon,
Go to step (8.14).
(8.13) photon is scattered, photon direction (μ ' after scatteringx,μ'y,μ'z) be updated to：
Here, (μx,μy,μz) represent the direction before photon equilibrium state, θ ' and ψ represent respectively the deflection angle before photon equilibrium state and
Azimuth, its expression formula is：
Scatter after the renewal of direction, go to step (8.2)；
(8.14) judge photon whether as last photon；If it is, continuing next step；If it is not, then on repeating
State process；
(9) CPU storage allocations, for depositing GPU output datas；By in the data duplication on video memory after calculating to internal memory,
Each child node returns to operation information to host node；
(10) the excitation light sub-information and fluorescent photon information that all effusions are exported after all photons are followed the trail of；
(11) releasing memory and video memory space and CUDA is exited；
(12) output in step (10) is based on, statistical calculation is carried out, other required transmission characteristics obtained are converted to.
Claims (3)
1. a kind of fluorescence Monte-Carlo Simulation Method accelerated based on concentrating type GPU, it is characterised in that comprise the following steps：
(1) space structure of target biological tissue is divided into a three-dimensional voxel model, setting one and three-dimensional voxel model
The numerical value of each element is the organization type of mark in size identical 3-dimensional digital matrix, matrix；Sets target biological tissue
Various organization optical property parameter, optical properties of tissue is the absorption coefficient of exciting light, the scattering coefficient of exciting light, glimmering
The absorption coefficient of light, the scattering coefficient of fluorescence, the absorption coefficient of fluorogen, the refraction coefficient of fluorogen and fluorogen it is each to different
Sex factor；
(2) by this message-passing communication agreements of MPI, host node obtains the quantity of the GPU equipment of each child node, it would be desirable to count
The quantity in the source of calculation is averagely allocated to the GPU equipment of each child node, and each child node determines the dimension of grid and block in GPU equipment
And size, and the subprocess of respective number is opened up according to the number of node GPU equipment, distribution of computation tasks is entered to each height
Journey, starts multiple programming and calculating platform CUDA, each child node CPU storage allocations space, video memory space, the data that will be calculated
Copied to from internal memory on video memory；
(3) incident light source is characterized as setting to the set of the photon of number, incident light source position and incident light direction are assigned to often
The initial position of individual photon and direction；
(4) transmitting procedure of each exciting light photon is followed the trail of in target biological tissue, the extensive of time length will be wherein expended
The step of data parallel, high calculating density, is arranged on GPU parallel execution, and remaining procedure is serially performed on CPU；CPU
Storage allocation, for depositing GPU output datas, by the data duplication on video memory after calculating to internal memory, each child node returns to fortune
Row information is to host node；Step (4) specifically includes following steps：
(8.1) exciting light photon is delivered, if photon is on tissue surface, the position is set as that photon launches position；
If photon is in outside tissue surface, photon is shifted into tissue surface automatically, its concrete operations uses alternative manner；If light
Son is in organization internal, then will directly launch photon, i.e. photon according to this position and interaction with surface does not occur, and direct
Start transmission of the photon in organization internal；
(8.2) the remaining step-length sleft of photon is calculated:Sleft=-In ξ；
(8.3) calculate photon and move the step-length s moved a step
S=min (- In (ξ)/μt,min(dx,dy,dz))
Here ξ is the pseudo random number for being evenly distributed on (0,1) produced by computer pseudorandom number generator, dx, dy, and dz is successively
For the length of each voxel in tissue model；When photon is the excitation light period of the day from 11 p.m. to 1 a.m：μt=μa+μs+μafp；When photon is fluorescence
During photon：μt=μaf+μsf；μaAnd μsIt is the absorption coefficient and scattering coefficient of exciting light respectively；μafAnd μsfIt is the suction of fluorescence respectively
Receive coefficient and scattering coefficient；μafpIt is the absorption coefficient of fluorogen；
(8.4) judge that photon is during being somebody's turn to do when the step-length calculated in front direction and (8.3) moves the step of photon one according to photon
The no interface through different tissues；If it is, into step (8.5)；If not, photon is moved into a step-length, update
Photon current location, subsequently into step (8.8)；
(8.5) rum point of the photon on interface is determined, method is：The minterm in following expression is found out first：
<mfenced open = "{" close = "">
<mtable>
<mtr>
<mtd>
<mi>x</mi>
<mo>_</mo>
<mi>t</mi>
<mo>=</mo>
<mo>(</mo>
<mn>1</mn>
<mo>+</mo>
<mo>&lsqb;</mo>
<mi>x</mi>
<mo>/</mo>
<mi>d</mi>
<mi>x</mi>
<mo>&rsqb;</mo>
<mo>-</mo>
<mi>x</mi>
<mo>/</mo>
<mi>d</mi>
<mi>x</mi>
<mo>)</mo>
<mo>/</mo>
<mi>u</mi>
<mi>x</mi>
</mtd>
</mtr>
<mtr>
<mtd>
<mi>y</mi>
<mo>_</mo>
<mi>t</mi>
<mo>=</mo>
<mo>(</mo>
<mn>1</mn>
<mo>+</mo>
<mo>&lsqb;</mo>
<mi>y</mi>
<mo>/</mo>
<mi>d</mi>
<mi>y</mi>
<mo>&rsqb;</mo>
<mo>-</mo>
<mi>y</mi>
<mo>/</mo>
<mi>d</mi>
<mi>y</mi>
<mo>)</mo>
<mo>/</mo>
<mi>u</mi>
<mi>y</mi>
</mtd>
</mtr>
<mtr>
<mtd>
<mi>z</mi>
<mo>_</mo>
<mi>t</mi>
<mo>=</mo>
<mo>(</mo>
<mn>1</mn>
<mo>+</mo>
<mo>&lsqb;</mo>
<mi>z</mi>
<mo>/</mo>
<mi>d</mi>
<mi>z</mi>
<mo>&rsqb;</mo>
<mo>-</mo>
<mi>z</mi>
<mo>/</mo>
<mi>d</mi>
<mi>z</mi>
<mo>)</mo>
<mo>/</mo>
<mi>u</mi>
<mi>z</mi>
</mtd>
</mtr>
</mtable>
</mfenced>
Here, x, y, z represent photon current location；If Section 1 is minimum, then application point expression formula of the photon on interface
For：
<mfenced open = "{" close = "">
<mtable>
<mtr>
<mtd>
<mrow>
<mi>x</mi>
<mn>0</mn>
<mo>=</mo>
<mi>x</mi>
<mo>+</mo>
<mi>u</mi>
<mi>x</mi>
<mo>*</mo>
<mi>x</mi>
<mo>_</mo>
<mi>t</mi>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mi>y</mi>
<mn>0</mn>
<mo>=</mo>
<mi>y</mi>
<mo>+</mo>
<mi>u</mi>
<mi>y</mi>
<mo>*</mo>
<mi>y</mi>
<mo>_</mo>
<mi>t</mi>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mi>z</mi>
<mn>0</mn>
<mo>=</mo>
<mi>z</mi>
<mo>+</mo>
<mi>u</mi>
<mi>z</mi>
<mo>*</mo>
<mi>z</mi>
<mo>_</mo>
<mi>t</mi>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
Here, x0, y0, z0 are application point of the photon at interface；
(8.6) direction of the photon after (8.5) middle rum point determined interacts with interface is calculated, specific method is：
Interface normal vector (a, b, c) is calculated first, and its expression formula is：
<mrow>
<mo>(</mo>
<mi>a</mi>
<mo>,</mo>
<mi>b</mi>
<mo>,</mo>
<mi>c</mi>
<mo>)</mo>
<mo>=</mo>
<mfrac>
<mrow>
<mo>(</mo>
<mo>&lsqb;</mo>
<mi>x</mi>
<mo>&rsqb;</mo>
<mo>-</mo>
<mo>&lsqb;</mo>
<mi>x</mi>
<mn>0</mn>
<mo>&rsqb;</mo>
<mo>,</mo>
<mo>&lsqb;</mo>
<mi>y</mi>
<mo>&rsqb;</mo>
<mo>-</mo>
<mo>&lsqb;</mo>
<mi>y</mi>
<mn>0</mn>
<mo>&rsqb;</mo>
<mo>,</mo>
<mo>&lsqb;</mo>
<mi>z</mi>
<mo>&rsqb;</mo>
<mo>-</mo>
<mo>&lsqb;</mo>
<mi>z</mi>
<mn>0</mn>
<mo>&rsqb;</mo>
<mo>)</mo>
</mrow>
<msup>
<mrow>
<mo>(</mo>
<msup>
<mrow>
<mo>(</mo>
<mo>&lsqb;</mo>
<mi>x</mi>
<mo>&rsqb;</mo>
<mo>-</mo>
<mo>&lsqb;</mo>
<mi>x</mi>
<mn>0</mn>
<mo>&rsqb;</mo>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mo>+</mo>
<msup>
<mrow>
<mo>(</mo>
<mo>&lsqb;</mo>
<mi>y</mi>
<mo>&rsqb;</mo>
<mo>-</mo>
<mo>&lsqb;</mo>
<mi>y</mi>
<mn>0</mn>
<mo>&rsqb;</mo>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mo>+</mo>
<msup>
<mrow>
<mo>(</mo>
<mo>&lsqb;</mo>
<mi>z</mi>
<mo>&rsqb;</mo>
<mo>-</mo>
<mo>&lsqb;</mo>
<mi>z</mi>
<mn>0</mn>
<mo>&rsqb;</mo>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mo>)</mo>
</mrow>
<mrow>
<mn>1</mn>
<mo>/</mo>
<mn>2</mn>
</mrow>
</msup>
</mfrac>
</mrow>
Then, incidence angle θ and refraction angle θ that photon is had an effect with interface are calculated according to the vector expressiont, its calculating formula is：
Cos θ=aux+buy+cuz
cosθt=(1- (1-cos2θ)·n1 2/n2 2)1/2
Here, kx,ky,kzRepresent that photon runs into the direction cosines value before border, n1And n2Expression runs into place group before and after interface
The refraction coefficient knitted；
Four kinds of situations calculating photons are finally divided to run into the direction after interface：
If cos θ=0, direction is constant after setting photon runs into border；
If cos θ=1, as ε ＞ (n2-n1)2/(n2+n1)2, direction is constant after setting photon runs into border, as ε≤(n2-n1
)2/(n2+n1)2, set photon to reflect, its direction is changed into：
<mfenced open = "{" close = "">
<mtable>
<mtr>
<mtd>
<msub>
<mi>k</mi>
<mrow>
<mi>x</mi>
<mi>i</mi>
</mrow>
</msub>
<mo>=</mo>
<msub>
<mi>k</mi>
<mi>x</mi>
</msub>
<mo>-</mo>
<mn>2</mn>
<mi>c</mi>
<mi>o</mi>
<mi>s</mi>
<mi>&theta;</mi>
<mo>*</mo>
<mi>a</mi>
</mtd>
</mtr>
<mtr>
<mtd>
<msub>
<mi>k</mi>
<mrow>
<mi>y</mi>
<mi>i</mi>
</mrow>
</msub>
<mo>=</mo>
<msub>
<mi>k</mi>
<mi>y</mi>
</msub>
<mo>-</mo>
<mn>2</mn>
<mi>c</mi>
<mi>o</mi>
<mi>s</mi>
<mi>&theta;</mi>
<mo>*</mo>
<mi>b</mi>
</mtd>
</mtr>
<mtr>
<mtd>
<msub>
<mi>k</mi>
<mrow>
<mi>z</mi>
<mi>i</mi>
</mrow>
</msub>
<mo>=</mo>
<msub>
<mi>k</mi>
<mi>z</mi>
</msub>
<mo>-</mo>
<mn>2</mn>
<mi>cos</mi>
<mi>&theta;</mi>
<mo>*</mo>
<mi>c</mi>
</mtd>
</mtr>
</mtable>
</mfenced>
If θ ＞ sin-1(n2/n1), set photon to be totally reflected, its direction is changed into：
<mfenced open = "{" close = "">
<mtable>
<mtr>
<mtd>
<msub>
<mi>k</mi>
<mrow>
<mi>x</mi>
<mi>i</mi>
</mrow>
</msub>
<mo>=</mo>
<msub>
<mi>k</mi>
<mi>x</mi>
</msub>
<mo>-</mo>
<mn>2</mn>
<mi>c</mi>
<mi>o</mi>
<mi>s</mi>
<mi>&theta;</mi>
<mo>*</mo>
<mi>a</mi>
</mtd>
</mtr>
<mtr>
<mtd>
<msub>
<mi>k</mi>
<mrow>
<mi>y</mi>
<mi>i</mi>
</mrow>
</msub>
<mo>=</mo>
<msub>
<mi>k</mi>
<mi>y</mi>
</msub>
<mo>-</mo>
<mn>2</mn>
<mi>c</mi>
<mi>o</mi>
<mi>s</mi>
<mi>&theta;</mi>
<mo>*</mo>
<mi>b</mi>
</mtd>
</mtr>
<mtr>
<mtd>
<msub>
<mi>k</mi>
<mrow>
<mi>z</mi>
<mi>i</mi>
</mrow>
</msub>
<mo>=</mo>
<msub>
<mi>k</mi>
<mi>z</mi>
</msub>
<mo>-</mo>
<mn>2</mn>
<mi>c</mi>
<mi>o</mi>
<mi>s</mi>
<mi>&theta;</mi>
<mo>*</mo>
<mi>c</mi>
</mtd>
</mtr>
</mtable>
</mfenced>
In other cases, the direction after setting photon is had an effect with interface is：
Wherein, r expression formula is：
<mrow>
<mi>r</mi>
<mo>=</mo>
<mfrac>
<mn>1</mn>
<mn>2</mn>
</mfrac>
<mo>&lsqb;</mo>
<msup>
<mrow>
<mo>(</mo>
<mfrac>
<mrow>
<msub>
<mi>n</mi>
<mn>1</mn>
</msub>
<mo>&CenterDot;</mo>
<msub>
<mi>cos&theta;</mi>
<mi>t</mi>
</msub>
<mo>-</mo>
<msub>
<mi>n</mi>
<mn>2</mn>
</msub>
<mo>&CenterDot;</mo>
<mi>cos</mi>
<mi>&theta;</mi>
</mrow>
<mrow>
<msub>
<mi>n</mi>
<mn>1</mn>
</msub>
<mo>&CenterDot;</mo>
<msub>
<mi>cos&theta;</mi>
<mi>t</mi>
</msub>
<mo>+</mo>
<msub>
<mi>n</mi>
<mn>2</mn>
</msub>
<mo>&CenterDot;</mo>
<mi>cos</mi>
<mi>&theta;</mi>
</mrow>
</mfrac>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mo>+</mo>
<msup>
<mrow>
<mo>(</mo>
<mfrac>
<mrow>
<msub>
<mi>n</mi>
<mn>1</mn>
</msub>
<mo>&CenterDot;</mo>
<mi>cos</mi>
<mi>&theta;</mi>
<mo>-</mo>
<msub>
<mi>n</mi>
<mn>2</mn>
</msub>
<mo>&CenterDot;</mo>
<msub>
<mi>cos&theta;</mi>
<mi>t</mi>
</msub>
</mrow>
<mrow>
<msub>
<mi>n</mi>
<mn>1</mn>
</msub>
<mo>&CenterDot;</mo>
<mi>cos</mi>
<mi>&theta;</mi>
<mo>+</mo>
<msub>
<mi>n</mi>
<mn>2</mn>
</msub>
<mo>&CenterDot;</mo>
<msub>
<mi>cos&theta;</mi>
<mi>t</mi>
</msub>
</mrow>
</mfrac>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mo>&rsqb;</mo>
</mrow>
(8.7) update the step-length s that has moved for photon current location had an effect into step (8.6) with interface rum point away from
From then by photon current location renewal to above-mentioned rum point；
(8.8) if photon is refracted in air in (8.6), as photon effusion tissue, setting photon death stops to light
Son is followed the trail of and updates sleft=0, and records the photon current attribute information, including position, direction, goes to step (8.14)；Its
Sleft=Sleft-s* μ in the case of himt；The μ if the photon is exciting light photont=μa+μs+μafp；When photon is fluorescence light
Sub then μt=μaf+μsf；If Sleft ＞ 0, into step (8.3)；Otherwise step (8.9) is entered；
(8.9) judge that photon is to be absorbed or scattered by setting random number, step is gone to if the photon is fluorescent photon
(8.10) step (8.11), is otherwise then gone to；
(8.10) as ξ≤μaf/(μaf+μafp) when, then fluorescent photon is absorbed, and setting photon is dead, is stopped following the trail of photon, is turned
To step (8.14)；If ξ ＞ μaf/(μaf+μafp) when, then fluorescent photon is scattered, and goes to step (8.13)；
(8.11) as ξ≤(μa+μafp)/(μa+μafp+μs) when, then exciting light photon is absorbed, into step (8.12)；As ξ ＞
(μa+μafp)/(μa+μafp+μs) when, then exciting light photon equilibrium state, goes to step (8.13)；
(8.12) judge whether photon is absorbed by fluorogen in biological tissue and produces fluorescence by setting random number, at the beginning of fluorescence
Beginning, direction determined by the deflection angle and azimuth of isotropic scatterning；As ξ≤μafp/(μa+μafp) when, then fluorescence is produced, is updated
Sleft=0, goes to step (8.2)；As ξ ＞ μafp/(μa+μafp) when, then photon death is set, stops following the trail of photon, goes to
Step (8.14)；
(8.13) photon is scattered, photon direction (μ ' after scatteringx,μ'y,μ'z) be updated to：
<mfenced open = "{" close = "">
<mtable>
<mtr>
<mtd>
<mrow>
<msubsup>
<mi>&mu;</mi>
<mi>x</mi>
<mo>&prime;</mo>
</msubsup>
<mo>=</mo>
<mfrac>
<mrow>
<msup>
<mi>sin&theta;</mi>
<mo>&prime;</mo>
</msup>
</mrow>
<msqrt>
<mrow>
<mn>1</mn>
<mo>-</mo>
<msubsup>
<mi>&mu;</mi>
<mi>z</mi>
<mn>2</mn>
</msubsup>
</mrow>
</msqrt>
</mfrac>
<mrow>
<mo>(</mo>
<msub>
<mi>&mu;</mi>
<mi>x</mi>
</msub>
<msub>
<mi>&mu;</mi>
<mi>z</mi>
</msub>
<mi>c</mi>
<mi>o</mi>
<mi>s</mi>
<mi>&psi;</mi>
<mo>-</mo>
<msub>
<mi>&mu;</mi>
<mi>y</mi>
</msub>
<mi>sin</mi>
<mi>&psi;</mi>
<mo>)</mo>
</mrow>
<mo>+</mo>
<msub>
<mi>&mu;</mi>
<mi>x</mi>
</msub>
<msup>
<mi>cos&theta;</mi>
<mo>&prime;</mo>
</msup>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<msubsup>
<mi>&mu;</mi>
<mi>y</mi>
<mo>&prime;</mo>
</msubsup>
<mo>=</mo>
<mfrac>
<mrow>
<msup>
<mi>sin&theta;</mi>
<mo>&prime;</mo>
</msup>
</mrow>
<msqrt>
<mrow>
<mn>1</mn>
<mo>-</mo>
<msubsup>
<mi>&mu;</mi>
<mi>z</mi>
<mn>2</mn>
</msubsup>
</mrow>
</msqrt>
</mfrac>
<mrow>
<mo>(</mo>
<msub>
<mi>&mu;</mi>
<mi>y</mi>
</msub>
<msub>
<mi>&mu;</mi>
<mi>z</mi>
</msub>
<mi>c</mi>
<mi>o</mi>
<mi>s</mi>
<mi>&psi;</mi>
<mo>-</mo>
<msub>
<mi>&mu;</mi>
<mi>x</mi>
</msub>
<mi>sin</mi>
<mi>&psi;</mi>
<mo>)</mo>
</mrow>
<mo>+</mo>
<msub>
<mi>&mu;</mi>
<mi>y</mi>
</msub>
<msup>
<mi>cos&theta;</mi>
<mo>&prime;</mo>
</msup>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<msubsup>
<mi>&mu;</mi>
<mi>z</mi>
<mo>&prime;</mo>
</msubsup>
<mo>=</mo>
<mo>-</mo>
<msup>
<mi>sin&theta;</mi>
<mo>&prime;</mo>
</msup>
<mi>cos</mi>
<mi>&psi;</mi>
<msqrt>
<mrow>
<mn>1</mn>
<mo>-</mo>
<msubsup>
<mi>&mu;</mi>
<mi>z</mi>
<mn>2</mn>
</msubsup>
</mrow>
</msqrt>
<mo>+</mo>
<msub>
<mi>&mu;</mi>
<mi>z</mi>
</msub>
<msup>
<mi>cos&theta;</mi>
<mo>&prime;</mo>
</msup>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
Here, (μx,μy,μz) direction before photon equilibrium state is represented, θ ' and ψ represent the deflection angle before photon equilibrium state and orientation respectively
Angle, its expression formula is：
Scatter after the renewal of direction, go to step (8.2)；
(8.14) judge photon whether as last photon；If it is, continuing next step；If it is not, then repeating above-mentioned mistake
Journey；
(5) such as exciting light photon is absorbed by the fluorogen of target biological tissue and produces fluorescence, follows the trail of the biography of each fluorescent photon
Defeated process；The step of the large-scale data for time length wherein expending is parallel, height calculates density will be arranged on GPU parallel execution,
Remaining procedure is serially performed on CPU；CPU storage allocations, for depositing GPU output datas, by video memory after calculating
Data duplication is on internal memory, and each child node returns to operation information to host node；
(6) the excitation light sub-information and fluorescent photon information that all effusions are exported after all photons are followed the trail of, releasing memory and aobvious
Deposit space and exit CUDA.
2. the fluorescence Monte-Carlo Simulation Method according to claim 1 accelerated based on concentrating type GPU, it is characterised in that step
Suddenly (5) specifically include following steps：
(3.1) if exciting light photon is absorbed by fluorogen and produces fluorescence, the inceptive direction of fluorescence is by isotropic scatterning
Deflection angle and azimuth are determined, continue to follow the trail of the movement of fluorescent photon next step；
(3.2) if fluorescent photon is scattered, photon direction cosines information is updated, step (3.3) is gone to；If fluorescent photon quilt
Absorb, then set photon death, stop following the trail of photon；
(3.3) fluorescent photon moves a step-length determined by random number, judges whether fluorescent photon runs into the process outer
Boundary border, if running into border, goes to step (3.4), otherwise goes to step (3.2)；
(3.4) if fluorescent photon reflects, the step-length and direction cosines information of fluorescent photon is updated, step is gone to
(3.2)；If fluorescent photon is refracted to outside tissue, setting photon is dead, stops following the trail of photon.
3. the fluorescence Monte-Carlo Simulation Method according to claim 1 accelerated based on concentrating type GPU, it is characterised in that：
Based on the output in step (6), statistics calculating and one or any be converted into following required transmission characteristics obtained are carried out
It is several：
By in same incident direction, the exciting light photon energy of all effusions is added, you can obtain exciting light effusion energy and diffusion
Relation curve between angle；
By in same incident direction, the fluorescence light photon energy of all effusions is added, you can obtain fluorescence effusion energy and scattering angle
Relation curve between degree.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201410534214.7A CN104331641B (en) | 2014-10-11 | 2014-10-11 | A kind of fluorescence Monte-Carlo Simulation Method accelerated based on concentrating type GPU |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201410534214.7A CN104331641B (en) | 2014-10-11 | 2014-10-11 | A kind of fluorescence Monte-Carlo Simulation Method accelerated based on concentrating type GPU |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN104331641A CN104331641A (en) | 2015-02-04 |
| CN104331641B true CN104331641B (en) | 2017-09-22 |
Family
ID=52406364
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201410534214.7A Active CN104331641B (en) | 2014-10-11 | 2014-10-11 | A kind of fluorescence Monte-Carlo Simulation Method accelerated based on concentrating type GPU |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN104331641B (en) |
Families Citing this family (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105574296A (en) * | 2016-02-03 | 2016-05-11 | 西安电子科技大学 | Electromagnetic scattering simulation method for ablated aircraft surface |
| CN107239352B (en) * | 2017-05-31 | 2019-11-29 | 北京科技大学 | The communication optimization method and its system of a kind of dynamics Monte Carlo Parallel Simulation |
| CN107944131A (en) * | 2017-11-22 | 2018-04-20 | 山东农业大学 | A kind of printing net-point coverage rate analog measurement method |
| CN110134508A (en) * | 2018-02-08 | 2019-08-16 | 北京连心医疗科技有限公司 | A kind of cloud Monte Carlo state machine system and framework method |
| CN109522106B (en) * | 2018-10-22 | 2023-01-17 | 广东工业大学 | Risk value simulation dynamic task scheduling method based on cooperative computing |
| CN109856092B (en) * | 2018-11-22 | 2021-08-31 | 上海无线电设备研究所 | Semi-analytic Monte Carlo simulation method for ocean fluorescence detection |
| CN113361080B (en) * | 2021-05-20 | 2022-05-17 | 厦门大学 | Multilayer water photon transmission semi-analytic Monte Carlo simulation method based on GPU |
| CN117272687B (en) * | 2023-11-20 | 2024-01-26 | 中国海洋大学 | Underwater optical imaging Monte Carlo vectorization isomerism parallel optimization method |
Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101268935A (en) * | 2008-04-25 | 2008-09-24 | 清华大学 | In vivo fluorescence numerator imaging modelling approach capable of calling multiple imaging algorithms |
-
2014
- 2014-10-11 CN CN201410534214.7A patent/CN104331641B/en active Active
Patent Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101268935A (en) * | 2008-04-25 | 2008-09-24 | 清华大学 | In vivo fluorescence numerator imaging modelling approach capable of calling multiple imaging algorithms |
Non-Patent Citations (3)
| Title |
|---|
| GPU集群加速的多层组织光子传输模型的蒙特卡罗仿真;蒋超 等;《中国科技论文在线》;20110228;第1-12页 * |
| Monte Carlo–based fluorescence molecular tomography reconstruction method accelerated by a cluster of graphic processing units;Guotao Quan et al.;《Journal of Biomedical Optics》;20110228;第16卷(第2期);第026018-1至026018-8页 * |
| 反射模式时域荧光分子层析原理的模拟和实验研究;和慧园;《中国优秀硕士学位论文全文数据库 医药卫生科技辑》;20090815(第08期);第24，30-33，36页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN104331641A (en) | 2015-02-04 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN104331641B (en) | A kind of fluorescence Monte-Carlo Simulation Method accelerated based on concentrating type GPU | |
| Zhou et al. | Real-time kd-tree construction on graphics hardware | |
| CN105138802B (en) | A kind of gun barrel intelligent design system and design method | |
| CN102628861B (en) | Method for simulating temperature cracking value of mass concrete | |
| CN104317655A (en) | Cluster GPU acceleration-based multi-source full path Monte-Carlo simulation method | |
| CN102282591A (en) | Ray tracing system architectures and methods | |
| CN105512366B (en) | The tree-shaped random seam net description method of compact reservoir volume fracturing containing intrinsic fracture | |
| CN110210178A (en) | A kind of construction method based on Python regeneration concrete three-dimensional random spherical shape aggregate model | |
| CN106339351A (en) | SGD (Stochastic Gradient Descent) algorithm optimization system and method | |
| CN106508012B (en) | Service-oriented group behavior parallel simulation method | |
| CN110442974A (en) | Horse shoe flame regenerator chamber of glass kiln performance optimization method and device | |
| Zhang et al. | Chaotic bean optimization algorithm | |
| CN107392446A (en) | A kind of step power station scheduling scheme evaluation method based on sensitivity analysis | |
| Yang et al. | Modeling urban design with energy performance | |
| CN103593504B (en) | A kind of based on the netting Reliability of Microprocessor emulation mode improving quality amplifying technique | |
| CN107515966A (en) | A kind of radar simulator system layering construction method based on DDS | |
| De Blas et al. | Assessment of sample size calculations used in aquaculture by simulation techniques | |
| Waibel et al. | Physics meets machine learning: Coupling FFD with regression models for wind pressure prediction on high-rise facades | |
| Passafaro et al. | Potential preferences for alternative forms of sustainable tourism: The case of rural and intergenerational tourism. | |
| Muehlbauer et al. | Automated shape design by grammatical evolution | |
| Bertels et al. | Why asynchronous parallel evolution is the future of hyper-heuristics: A cdcl sat solver case study | |
| Llabres et al. | Relational urban models: Parameters, values and tacit forms of algorithms | |
| Ye | Optimization of resource scheduling based on genetic algorithm in cloud computing environment | |
| Venetillo et al. | GPU-based particle simulation with inter-collisions | |
| Davoli et al. | Innovative methods for a sustainable retrofit of the existing building stock. A cross-path from social housing to the listed heritage |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |