CN108618796A - A kind of Monte Carlo scattered photon analogy method of task based access control driving - Google Patents
A kind of Monte Carlo scattered photon analogy method of task based access control driving Download PDFInfo
- Publication number
- CN108618796A CN108618796A CN201810133436.6A CN201810133436A CN108618796A CN 108618796 A CN108618796 A CN 108618796A CN 201810133436 A CN201810133436 A CN 201810133436A CN 108618796 A CN108618796 A CN 108618796A
- Authority
- CN
- China
- Prior art keywords
- photon
- path
- scattering
- probability
- point
- 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.)
- Pending
Links
- 238000000034 method Methods 0.000 title claims abstract description 69
- 238000005070 sampling Methods 0.000 claims abstract description 66
- 238000004088 simulation Methods 0.000 claims abstract description 33
- 230000008021 deposition Effects 0.000 claims abstract description 9
- 238000005295 random walk Methods 0.000 claims abstract description 6
- 239000006185 dispersion Substances 0.000 claims description 7
- 238000009826 distribution Methods 0.000 claims description 7
- 230000008569 process Effects 0.000 claims description 7
- 239000007787 solid Substances 0.000 claims description 4
- 230000003993 interaction Effects 0.000 claims description 3
- 230000007704 transition Effects 0.000 claims description 3
- 238000007408 cone-beam computed tomography Methods 0.000 description 8
- 238000012937 correction Methods 0.000 description 7
- 238000010586 diagram Methods 0.000 description 5
- 241000272443 Penelope Species 0.000 description 3
- 239000002245 particle Substances 0.000 description 3
- 238000012545 processing Methods 0.000 description 3
- 230000007812 deficiency Effects 0.000 description 2
- 238000013461 design Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 210000005075 mammary gland Anatomy 0.000 description 2
- 239000011159 matrix material Substances 0.000 description 2
- 238000001959 radiotherapy Methods 0.000 description 2
- 239000000243 solution Substances 0.000 description 2
- 206010008111 Cerebral haemorrhage Diseases 0.000 description 1
- 238000000342 Monte Carlo simulation Methods 0.000 description 1
- 206010028980 Neoplasm Diseases 0.000 description 1
- RTAQQCXQSZGOHL-UHFFFAOYSA-N Titanium Chemical compound [Ti] RTAQQCXQSZGOHL-UHFFFAOYSA-N 0.000 description 1
- 230000003044 adaptive effect Effects 0.000 description 1
- 210000000988 bone and bone Anatomy 0.000 description 1
- 238000000205 computational method Methods 0.000 description 1
- 238000002591 computed tomography Methods 0.000 description 1
- 238000012790 confirmation Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 239000003814 drug Substances 0.000 description 1
- 230000008030 elimination Effects 0.000 description 1
- 238000003379 elimination reaction Methods 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 230000036541 health Effects 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 239000004615 ingredient Substances 0.000 description 1
- 238000002347 injection Methods 0.000 description 1
- 239000007924 injection Substances 0.000 description 1
- 230000005433 particle physics related processes and functions Effects 0.000 description 1
- 230000035515 penetration Effects 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 230000011218 segmentation Effects 0.000 description 1
- 238000012546 transfer Methods 0.000 description 1
Classifications
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
- A61B6/02—Arrangements for diagnosis sequentially in different planes; Stereoscopic radiation diagnosis
- A61B6/03—Computed tomography [CT]
- A61B6/032—Transmission computed tomography [CT]
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
- A61B6/40—Arrangements for generating radiation specially adapted for radiation diagnosis
- A61B6/4064—Arrangements for generating radiation specially adapted for radiation diagnosis specially adapted for producing a particular type of beam
- A61B6/4085—Cone-beams
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
- A61B6/52—Devices using data or image processing specially adapted for radiation diagnosis
- A61B6/5258—Devices using data or image processing specially adapted for radiation diagnosis involving detection or reduction of artifacts or noise
- A61B6/5282—Devices using data or image processing specially adapted for radiation diagnosis involving detection or reduction of artifacts or noise due to scatter
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
- A61B6/58—Testing, adjusting or calibrating thereof
- A61B6/582—Calibration
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T11/00—2D [Two Dimensional] image generation
- G06T11/003—Reconstruction from projections, e.g. tomography
- G06T11/008—Specific post-processing after tomographic reconstruction, e.g. voxelisation, metal artifact correction
Landscapes
- Engineering & Computer Science (AREA)
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Physics & Mathematics (AREA)
- Medical Informatics (AREA)
- Animal Behavior & Ethology (AREA)
- Molecular Biology (AREA)
- Biophysics (AREA)
- High Energy & Nuclear Physics (AREA)
- Theoretical Computer Science (AREA)
- Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
- Optics & Photonics (AREA)
- Pathology (AREA)
- Radiology & Medical Imaging (AREA)
- Biomedical Technology (AREA)
- Heart & Thoracic Surgery (AREA)
- Veterinary Medicine (AREA)
- Surgery (AREA)
- General Physics & Mathematics (AREA)
- General Health & Medical Sciences (AREA)
- Public Health (AREA)
- Data Mining & Analysis (AREA)
- Mathematical Analysis (AREA)
- Mathematical Physics (AREA)
- Pure & Applied Mathematics (AREA)
- Mathematical Optimization (AREA)
- Computational Mathematics (AREA)
- Pulmonology (AREA)
- Evolutionary Biology (AREA)
- Bioinformatics & Computational Biology (AREA)
- Operations Research (AREA)
- Probability & Statistics with Applications (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Algebra (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- General Engineering & Computer Science (AREA)
- Measurement Of Radiation (AREA)
Abstract
A kind of Monte Carlo scattered photon analogy method of task based access control driving, it is as follows, photon equilibrium state parameter is preset and initialized, generate a primary photon path, position sampling is carried out to scattering point using uniform sampling algorithm and random walk sampling algorithm, generate a simulated photons path, compare the probability in primary photon path and simulated photons path, it is sampled automatically to simulated photons path using path sampling algorithm, and judge whether the photon path quantity of deposition meets preset photon path quantity, if, then end operation, otherwise second step is returned.This method promotes photon equilibrium state simulation precision from the level of Sampling, actively simulation mission requirements are included in controllable path variation space, it is sampled into row variation to photon path by automatic sampling model, then the importance of relatively more front and back sample path, select the path relatively important to simulation task as present energy deposition path, to realize the automatic importance sampling of photon path.
Description
Technical field
The present invention relates to scattered photon analogue technique fields, more particularly to a kind of Monte Carlo of task based access control driving
Scattered photon analogy method.
Background technology
Pencil-beam computed tomography (Cone-beam Computed Tomography, CBCT) system is due to its body
Product is small, and light-weight, sweep speed is widely used in the numerous clinical departments of hospital soon, especially in dentistry CT, mammary gland CT and
In termed image-guided radiotherapy.The detector used due to CBCT systems can not be blocked for flat panel detector using collimator technology
Scattered rays, so there are serious scatter artefacts in the CBCT images rebuild.
To make CBCT images deeply be applied in clinic, domestic and international researcher has carried out greatly the correction of scatter artefacts
Quantifier elimination.Boone J M are summarized as two kinds of scatter correction techniques types based on hardware and software.It is wherein based on hard
The scatter correction method of part can bring other various related problems and defect, finally lose more than gain.Another kind of method is to be based on
The scatter correction method of software realizes scatter correction by the estimation of scatter distributions is realized in calculating.Wherein to design scattering
The method that kernel function carries out deconvolution is representative.This method calculating speed is fast, has good clinical value, but scatter
Nuclear design difficulty causes scattering estimation inaccurate, can only carry out the confirmation of pendulum position and school using bone tissue form and profile at present
Positive (Sun M, Star-Lack J:Improved scatter correction using adaptive scatter
kernel superposition.Physics in medicine and biology 2010,55(22):6695-6720.)。
It is by Monte Carlo (Monte Carlo, MC) SIMULATED SCATTERING light in addition there are another kind of more accurate method
Muon physics effect obtains scattering component and realizes scatter correction.It is that tumour radiotherapy is controlled based on the computational methods that MC simulation particles transport
Goldstandard of the treatment field in relation to particle physics transport issues.MCNP (the Brown of conventional MC analogy methods such as Brown FB exploitations
FB:MCNP–A General Monte Carlo N-Particle Transport Code,Version 5.Los Alamos
National Laboratory, Oak Ridge, TN 2003.) and Xun Jia exploitation gMCDRR (Jia X, Yan H, L,Folkerts M,Jiang SB:A GPU tool for efficient,accurate,and
Realistic simulation ofcone beam CT projections.Medical Physics 2012), compared to
PENELOPE(Schmidtke J,Engel W.PENELOPE:An algorithm for Monte Carlo simulation
of the penetration and energy loss of electrons and positrons in matter[J]
.Nuclear Instruments&Methods in Physics Research,1995,100(1):31-46.) these are classical
MC is simulated, and has promotion on simulation precision and accuracy.But either PENELOPE, MCNP or gMCDRR, they
Simulation is all passive type and inactive mode.Passive type simulation refers in photon MC simulations, successively to the injection side of photon
To position, type and angle, position, type and the scattering angle etc. that second scattering point occurs occur for first scattering point
Independent sampling finally judges whether photon reaches detector to decide whether to deposit the path until photon leaves scanning object
The energy of photon.Passive type photon equilibrium state independent sampling strategy thinking relatively nature and code structure is simple, but the disadvantage is that can only
After simulating the photon transport of all scattering exponent numbers and ingredient, is passively picked out from result and deposit to detector region
Scattering component.And the photon for largely failing to reach detector is then dropped, and is caused to spend a large amount of time and is simulated to final
Photon of the scattered signal without contribution.So routine MC analogy methods are due to can not calculate the probability in photon entirety path,
Also the significance distribution for just ignoring photon entirety path, causes simulation precision too low.
The bottleneck problem that convergence rate to solve MC simulations is slow and efficiency is too low, common way are to introduce various subtract
Skill (Haghighat A, the Wagner JC of variance:Monte Carlo variance reduction with
deterministic importance functions. Progress in Nuclear Energy 2003,42(1):25-
53.) and use graphics processing unit (Graphics Processing Unit, GPU) parallel processing technique.But likewise, this
Root cause problems of the method for two kinds of improvement without solving passive type simulation.And, the skill for subtracting variance only has user to provide conjunction again
Suitable particle position significance distribution, is just improved the possibility of computational efficiency, otherwise result meeting altering error, and provides properly
Significance distribution be typically the most difficult.Exactly because the root cause problems of passive type simulation are not resolved, to such as
Mammary gland CBCT and cerebral hemorrhage CBCT etc. are impossible to the scattering higher scattering analogue task of required precision.
Therefore, in view of the shortcomings of the prior art, providing a kind of Monte Carlo scattered photon analogy method of task based access control driving
It is very necessary to solve prior art deficiency.
Invention content
A kind of Meng Teka of task based access control driving is provided it is an object of the invention to avoid the deficiencies in the prior art place
Lip river scattered photon analogy method, the Monte Carlo scattered photon analogy method of task based access control driving, from the level of Sampling
Photon equilibrium state simulation precision is promoted, actively simulation mission requirements are included in controllable path variation space, by automatic
Sampling model samples to photon path into row variation, then the importance of relatively more front and back sample path, and selection is to simulating task phase
To important path as present energy deposition path, to realize the automatic importance sampling of photon path.Task-driven
Path sampling method has abandoned the independent sampling strategy of photon, and the active regulation and control scattering task of energy so that simulation precision is big
It is big to improve.
The above-mentioned purpose of the present invention is realized by following technological means.
A kind of Monte Carlo scattered photon analogy method of task based access control driving is provided, is as follows：
The first step：Photon equilibrium state parameter is preset and initialized, a primary photon path is generated；
Second step：Position sampling is carried out to scattering point using uniform sampling algorithm and random walk sampling algorithm, generates one
Simulated photons path；
Third walks：Compare the probability in primary photon path and simulated photons path；
4th step：It is sampled automatically to simulated photons path using path sampling algorithm；
5th step：Simulated photons energy is deposited in corresponding detector, and judges that the photon path quantity of deposition is
It is no to meet preset photon path quantity, if so, end operation, otherwise returns to second step.
Specifically, the parameter of photon equilibrium state includes photon equilibrium state energy, type, exponent number, path probability and photon road
Diameter quantity, and the product that the total probability in photon equilibrium state path is each segment probability in all segments.
Preferably, the algorithm of third step photon path probability is specific as follows：
If photon is in from source point to dA detector pixels, after A1、…、Ai、…、ANN number of scattering point scattering process,
By N+1 segment, each segment is denoted as lkThe length of (k=1,2 ..., N+1), homologous segment are sk, scattering point AiIt is corresponding
Any one scattering type T in Rayleigh scattering and Compton scatteringj(j=1 or 2), T1For Rayleigh scattering, T2For Compton
Scattering, scattering point AiIn xiCompton scattering occurs for position or the linear attenuation coefficient of Rayleigh scattering isScattering point Ai
In xiIt is μ (x that Compton scattering and the linear attenuation coefficient of Rayleigh scattering, which occur, for positioni)；
As k=1, photon passes through initial segment l1Probability function P1Method for solving it is as follows：
Wherein, s1For segment l1Length, the Direction Probability function F (x) of distribution of photons, according to Beer-Lambert Law, light
Son reaches x1Point probability function bePhoton is in x1T occurs for pointjThe probability of the scattering process of type
Function isIt is tired that multiply three kinds of probability functions be initial segment l1Probability function, as shown in formula (1)；
As 1 ＜ k≤N, photon passes through intermediate segment l2~lNIn any one segment probability function it is as follows：
Wherein,For Ai、Ai+1The differenta1 scattering cross-section coefficient in 2 points of be in line directions, and Tm∈Tj, photon from
Scattering point AiMove to Ai+1And T occursjType interaction probability be
As k=N+1, photon is by ending segment lN+1Probability function PN+1Method for solving it is as follows：
Wherein, α lN+1The normal direction angle for transporting direction and detector of segment, photon is in ANIt is inclined that scattering occurs for point
The differential scattering coefficient for going to dA detector members isDetector pixel dA identical element corresponding AsNThe solid angle of point
ForPhoton is from ANThe probability for reaching detector dA positions is
The total probability function of photon fullpath is as follows：
Specifically, the concrete operations of the 4th step are as follows：
S1, the probability function P that primary photon path is calculated separately according to formula (4)totalWith the probability letter in simulated photons path
Number P 'total；
S2, acceptance probability P is judged according to formula (5)acceptValue, if Paccept>=1, then probability function P 'totalPlace path
For final path, and enter the 5th step；
If Paccept＜ 1, then enter step S3；
S3, sample a random number ξ from [0,1], if ξ ＜ Paccept, then probability function P 'totalPlace path is final
Path simultaneously enters the 5th step；
Otherwise, probability function P ' is definedtotalPlace path is primary photon path, and enters the 5th step.
Further, probability function PacceptCalculating steps are as follows：
Wherein, q (x) and q (x) ' correspond to P respectivelyTotalAnd P 'totalThe path transition function of place photon path.
Specifically, uniform probability density function or Gaussian probability-density function are symmetric function, so formula (5) can
It is reduced to：
Further, the 5th step records the photon path simulation number of Compton scattering and Rayleigh scattering signal respectively.
Preferably, in second step, corresponding sampling algorithm is chosen according to scattering point position difference and carries out random sampling behaviour
Make, it is specific as follows：
P1, the position for judging scattering point；
If single order scattering point, then randomly selected using uniform sampling algorithm；
Otherwise, P2 is entered step；
P2, if multistage scattering point, then judge scattering point occur scattering type；
If happens is that Compton scattering, is randomly selected using uniform sampling algorithm；
If happens is that Rayleigh scattering, is randomly selected using RWS algorithms.
Further, in multistage scattering point, the method using uniform sampling algorithm is specific as follows：
A scatteringangleθ is randomly selected within the scope of 360 degreeC, in scatteringangleθCOpening direction on choose some conducts
Compton second order dispersion point.
Specifically, in multistage scattering point, the method using RWS algorithms is specific as follows：
K1, judge probability function P 'totalPlace path is final path；
If probability function P 'totalPlace path is final path, then by carrying out height by a small margin in current angle of scattering
This offset, obtains new angle of scattering；
Otherwise, K2 is entered step；
K2, a random scatteringangleθR, with σ2Gauss sampling is carried out for variance, obtains new scatteringangleθR', in θR' side
It uniformly chooses a little as new scattering point upwards.
The present invention promotes photon equilibrium state simulation precision from the level of Sampling, and actively simulation mission requirements are included in can
In the path variation space of control, sampled into row variation to photon path by automatic sampling model, then relatively more front and back sampling
The importance in path selects the path relatively important to simulation task as present energy deposition path, to realize photon
The automatic importance sampling in path.The path sampling method of task-driven has abandoned the independent sampling strategy of photon, and can be actively
Formula regulates and controls scattering task so that simulation precision greatly improves.
Description of the drawings
Using attached drawing, the present invention is further illustrated, but the content in attached drawing does not constitute any limit to the present invention
System.
Using attached drawing, the present invention is further illustrated, but the content in attached drawing does not constitute any limit to the present invention
System.
Fig. 1 is a kind of system block diagram of the Monte Carlo scattered photon analogy method of task based access control driving of the present invention.
Fig. 2 is that circular orbit CBCT scans schematic diagram.
Fig. 3 is the photon path figure of segmentation display.
Fig. 4 is the schematic diagram of three photon paths generated after path makes a variation.
Fig. 5 is the partial code schematic diagram of automated path sampling algorithm.
Fig. 6 is the photon that the photon path simulated using the method for the present invention in embodiment 2 is simulated with gMCDRR methods
The effect contrast figure simulated using the method for the present invention and gMCDRR methods in the comparison diagram and embodiment 2 in path.
Specific implementation mode
The invention will be further described with the following Examples.
Embodiment 1.
As shown in Figs. 1-5, a kind of Monte Carlo scattered photon analogy method of task based access control driving, is as follows：
The first step：Photon equilibrium state parameter is preset and initialized, a primary photon path is generated.
Wherein, the parameter of photon equilibrium state includes photon equilibrium state energy, type, exponent number, path probability and photon path number
Amount, and the product that the total probability in photon equilibrium state path is each segment probability in all segments.
Second step：Position sampling is carried out to scattering point using uniform sampling algorithm and random walk sampling algorithm, generates one
Simulated photons path.
In second step, corresponding sampling algorithm is chosen according to scattering point position difference and carries out random sampling operation, specifically such as
Under：
P1, the position for judging scattering point.
If single order scattering point, then randomly selected using uniform sampling algorithm.
Otherwise, P2 is entered step.
P2, if multistage scattering point, then judge scattering point occur scattering type.
If happens is that Compton scattering, is randomly selected using uniform sampling algorithm.
For in multistage scattering point, the method using uniform sampling algorithm is specific as follows：
A scatteringangleθ is randomly selected within the scope of 360 degreeC, in scatteringangleθCOpening direction on choose some conducts
Compton second order dispersion point.
If happens is that Rayleigh scattering, is randomly selected using RWS algorithms.
For in multistage scattering point, the method using RWS algorithms is specific as follows：
K1, judge probability function P 'totalPlace path is final path.
If probability function P 'totalPlace path is final path, then by carrying out height by a small margin in current angle of scattering
This offset, obtains new angle of scattering.
Otherwise, K2 is entered step.
K2, a random scatteringangleθR, with σ2Gauss sampling is carried out for variance, obtains new scatteringangleθR', in θR' side
It uniformly chooses a little as new scattering point upwards.
Third walks：Compare the probability in primary photon path and simulated photons path.
The algorithm that third walks photon path probability is specific as follows：
If photon is in from source point to dA detector pixels, after A1、…、Ai、…、ANN number of scattering point scattering process,
By N+1 segment, each segment is denoted as lkThe length of (k=1,2 ..., N+1), homologous segment are sk, scattering point AiIt is corresponding
Any one scattering type T in Rayleigh scattering and Compton scatteringj(j=1 or 2), T1For Rayleigh scattering, T2For Compton
Scattering, scattering point AiIn xiCompton scattering occurs for position or the linear attenuation coefficient of Rayleigh scattering isScattering point Ai
In xiIt is μ (x that Compton scattering and the linear attenuation coefficient of Rayleigh scattering, which occur, for positioni)。
As k=1, photon passes through initial segment l1Probability function P1Method for solving it is as follows：
Wherein, s1For segment l1Length, the Direction Probability function F (x) of distribution of photons, according to Beer-Lambert Law, light
Son reaches x1Point probability function bePhoton is in x1T occurs for pointjThe probability of the scattering process of type
Function isIt is tired that multiply three kinds of probability functions be initial segment l1Probability function, as shown in formula (1).
As 1 ＜ k≤N, photon passes through intermediate segment l2~lNIn any one segment probability function it is as follows：
Wherein,For Ai、Ai+1The differenta1 scattering cross-section coefficient in 2 points of be in line directions, and Tm∈Tj, photon from
Scattering point AiMove to Ai+1And T occursjType interaction probability be
As k=N+1, photon is by ending segment lN+1Probability function PN+1Method for solving it is as follows：
Wherein, α lN+1The normal direction angle for transporting direction and detector of segment, photon is in ANIt is inclined that scattering occurs for point
The differential scattering coefficient for going to dA detector members isDetector pixel dA identical element corresponding AsNThe solid angle of point
ForPhoton is from ANThe probability for reaching detector dA positions is
The total probability function of photon fullpath is as follows：
4th step：It is sampled automatically to simulated photons path using path sampling algorithm.
The concrete operations of 4th step are as follows：
S1, the probability function P that primary photon path is calculated separately according to formula (4)totalWith the probability letter in simulated photons path
Number P 'total。
S2, acceptance probability P is judged according to formula (5)acceptValue, if Paccept>=1, then probability function P 'totalPlace path
For final path, and enter the 5th step.
If Paccept＜ 1, then enter step S3.
S3, sample a random number ξ from [0,1], if ξ ＜ Paccept, then probability function P 'totalPlace path is final
Path simultaneously enters the 5th step.
Otherwise, probability function P ' is definedtotalPlace path is primary photon path, and enters the 5th step.
Wherein, probability function PacceptCalculating steps are as follows：
Wherein, q (x) and q (x) ' correspond to P respectivelyTotalAnd P 'totalThe path transition function of place photon path.
Uniform probability density function or Gaussian probability-density function are symmetric function, so formula (5) can be reduced to：
Specific algorithm is as shown in figure (5).
5th step：Simulated photons energy is deposited in corresponding detector, and judges that the photon path quantity of deposition is
It is no to meet preset photon path quantity, if so, end operation, otherwise returns to second step.
5th step records the photon path simulation number of Compton scattering and Rayleigh scattering signal respectively.
Photon equilibrium state simulation precision is promoted from the level of Sampling, actively simulation mission requirements are included in controllable
Path makes a variation in space, is sampled into row variation to photon path by automatic sampling model, then compares front and back sample path
Importance selects the path relatively important to simulation task as present energy deposition path, to realize photon path
Automatic importance sampling.The path sampling method of task-driven has abandoned the independent sampling strategy of photon, and can active regulation and control
Scattering task so that simulation precision greatly improves.
Embodiment 2.
A kind of Monte Carlo scattered photon analogy method of task based access control driving, other feature is same as Example 1, no
It is with place：As shown in fig. 6, the embodiment operates for realistic simulation, data that are used below and obtaining are practical behaviour
Make gained, there is authenticity.
The place platform of the simulated operation is ubuntu-12.04.4 system GPU architectures, and video card equipment is NVIDIA
GeForce GTX TITAN Z, x-ray source energy are 60kVp, and used homogeneous phantom size is 10 × 10 × 2.8cm3,
Matrix size is 250 × 250 × 280, and the size of detector is 40cm × 30cm, and matrix size is 512 × 384, and radiographic source arrives
Rotation center distance and the distance of rotation center to detector are respectively 15.59cm and 49.41cm, below specifically to simulate
Step：
The first step carries out photon equilibrium state energy, type, exponent number, the default and initialization of the projects such as path probability, generation
Original path X.
Photon energy is set as single energy 60kVp in the debugging stage, and scattering type is mainly Compton scattering and Rayleigh scattering
Two kinds, a Type Control threshold value CONTROLTYPE=0.5 is set.
Then random number ζ is uniformly extracted in [0,1], if ζ<0.5, then it is assumed that Compton scattering occurs for the point, otherwise recognizes
For Rayleigh scattering occurs.
It is interfered with each other in order to prevent in the present invention, difference scattering exponent number is separately simulated.
Second step：Using uniform sampling algorithm and random walk sampling algorithm (Random-Walk Sampling, RWS)
It is scattered a position sampling, generates new photon path Y.
For single order scattering point, either Compton scattering or Rayleigh scattering, it is all made of uniform sampling strategy.
A position is randomly selected in homogeneous phantom as scattering process point.For second order dispersion point, Compton
Scattering and Rayleigh scattering take the different methods of samplings respectively.
Uniform sampling algorithm is used to Compton scattering, randomly selects a scatteringangleθ within the scope of 360 degree firstC,
Then it is uniformly chosen on this angle direction again and is a little used as Compton second order dispersion point.
RWS is taken to Rayleigh scattering, i.e., if receiving current path, height by a small margin is carried out under current angle of scattering
This offset, to obtain new angle of scattering.
Specifically, one scatteringangleθ of random initializtionR, with this θRFor mean value, with σ2Gauss sampling (σ is carried out for variance2's
Selection is arbitrary, but generally takes smaller value), obtain new scatteringangleθR', and in θR' some conducts are uniformly chosen on direction
Rayleigh second order dispersion point.
Third walks：Calculate length l of each photon segment in scanning objectk, linear attenuation coefficient mu (x) is inquired in terms of
Calculate photon path integralInquiry differential scattering data simultaneously calculate solid angleAccording to public affairs
Formula (1), (2), (3) calculate P1、PkAnd PN+1.The probability function P in primary photon path is finally calculated according to formula (4)totalAnd mould
The probability function P ' of quasi- photon pathtotal。
4th step：It samples automatically photon path using path sampling algorithm, since the path transfer chosen in the present invention is general
Rate function is uniform probability density function or Gaussian probability-density function, these probability density functions are symmetrical, so q (x)
=q (x) ', formula (5) are as follows after simplifying：
Judge acceptance probability PacceptValue, if Paccept>=1, then probability function P 'totalPlace path is final path,
And enter the 5th step.
If Paccept＜ 1, then enter step S3.
S3, sample a random number ξ from [0,1], if ξ ＜ Paccept, then probability function P 'totalPlace path is final
Path simultaneously enters the 5th step.
Otherwise, probability function P ' is definedtotalPlace path is primary photon path, and enters the 5th step.
5th step：Photon energy is to corresponding detector pixel point in deposition path, records Compton scattering and auspicious respectively
Sharp scattered signal, and judge whether to complete required path simulation number.
If so, end simulation, otherwise returns to second step.
Simulation result is as follows：
Figure (a)-(d) in Fig. 6 is respectively the single order Compton scattering that this method is simulated, single order Rayleigh scattering, and two
Rank Compton scattering and second order Rayleigh scattering signal.
The single order Compton scattering that figure (e)-(h) in Fig. 6 is simulated by gMCDRR, single order Rayleigh scattering, second order health
Pu Dun is scattered and second order Rayleigh scattering signal.
Figure (i)-(l) in Fig. 6 corresponds to the line drawing of yellow reference line in figure (a)-(h) figures in Fig. 6.
Contrast experiment does not know to compare under horizontal δ (uncertainty level) same, whereinN is detector pixel number, stdiSignal standards for ith pixel point is poor.
As it can be seen that the result of this method simulation and the result that gMCDRR is simulated are almost the same, also from figure (a)-figure (i)
It is to say that this method can reach dummy level identical with gMCDRR.
Secondly, the simulated time of two methods has also been counted, as shown in table (1)：
Table (1)
As seen from the table, this method analog result under the same conditions, the speed of simulation is obviously better than the side gMCDRR
Method.
For single order scattering, method accelerates about 208 times；For second order dispersion, about 12 are accelerated
Times.
This illustrates that it is feasible to improve MC scattering analogue efficiency from Sampling level.
Finally it should be noted that the above embodiments are merely illustrative of the technical solutions of the present invention rather than is protected to the present invention
The limitation of range is protected, although being explained in detail to the present invention with reference to preferred embodiment, those skilled in the art should
Understand, technical scheme of the present invention can be modified or replaced equivalently, without departing from the reality of technical solution of the present invention
Matter and range.
Claims (10)
1. a kind of Monte Carlo scattered photon analogy method of task based access control driving, it is characterised in that：It is as follows：
The first step：Photon equilibrium state parameter is preset and initialized, a primary photon path is generated；
Second step：Position sampling is carried out to scattering point using uniform sampling algorithm and random walk sampling algorithm, generates a mould
Quasi- photon path；
Third walks：Compare the probability in primary photon path and simulated photons path；
4th step：It is sampled automatically to simulated photons path using path sampling algorithm；
5th step：Simulated photons energy is deposited in corresponding detector, and judges whether the photon path quantity of deposition is full
The preset photon path quantity of foot, if so, end operation, otherwise returns to second step.
2. a kind of Monte Carlo scattered photon analogy method of task based access control driving according to claim 1, feature exist
In：The parameter of photon equilibrium state includes photon equilibrium state energy, type, exponent number, path probability and photon path quantity, and photon equilibrium state
The total probability in path is the product of each segment probability in all segments.
3. a kind of Monte Carlo scattered photon analogy method of task based access control driving according to claim 2, feature exist
In：The algorithm that third walks photon path probability is specific as follows：
If photon is in from source point to dA detector pixels, after A1、…、Ai、…、ANN number of scattering point scattering process, by N
+ 1 segment, each segment are denoted as lkThe length of (k=1,2 ..., N+1), homologous segment are sk, scattering point AiCorresponding Rayleigh dissipates
It penetrates and any one scattering type T in Compton scatteringj(j=1 or 2), T1For Rayleigh scattering, T2For Compton scattering, scattering
Point AiIn xiCompton scattering occurs for position or the linear attenuation coefficient of Rayleigh scattering isScattering point AiIn xiPosition occurs
The linear attenuation coefficient of Compton scattering and Rayleigh scattering is μ (xi)；
As k=1, photon passes through initial segment l1Probability function P1Method for solving it is as follows：
Wherein, s1For segment l1Length, the Direction Probability function F (x) of distribution of photons, according to Beer-Lambert Law, photon reaches
x1Point probability function bePhoton is in x1T occurs for pointjThe probability function of the scattering process of type isIt is tired that multiply three kinds of probability functions be initial segment l1Probability function, as shown in formula (1)；
As 1 ＜ k≤N, photon passes through intermediate segment l2~lNIn any one segment probability function it is as follows：
Wherein,For Ai、Ai+1The differenta1 scattering cross-section coefficient in 2 points of be in line directions, and Tm∈Tj, photon is from scattering
Point AiMove to Ai+1And T occursjType interaction probability be
As k=N+1, photon is by ending segment lN+1Probability function PN+1Method for solving it is as follows：
Wherein, α lN+1The normal direction angle for transporting direction and detector of segment, photon is in ANPoint occurs scattering and deflects into dA
The differential scattering coefficient of detector member isDetector pixel dA identical element corresponding AsNPoint solid angle bePhoton is from ANThe probability for reaching detector dA positions is
The total probability function of photon fullpath is as follows：
4. a kind of Monte Carlo scattered photon analogy method of task based access control driving according to claim 3, feature exist
In：The concrete operations of 4th step are as follows：
S1, the probability function P that primary photon path is calculated separately according to formula (4)totalWith the probability function in simulated photons path
P’total；
S2, acceptance probability P is judged according to formula (5)acceptValue, if Paccept>=1, then probability function P 'totalPlace path is most
Whole path, and enter the 5th step；
If Paccept＜ 1, then enter step S3；
S3, sample a random number ξ from [0,1], if ξ ＜ Paccept, then probability function P 'totalPlace path is final path
And enter the 5th step；
Otherwise, probability function P ' is definedtotalPlace path is primary photon path, and enters the 5th step.
5. a kind of Monte Carlo scattered photon analogy method of task based access control driving according to claim 4, feature exist
In：Probability function PacceptCalculating steps are as follows：
Wherein, q (x) and q (x) ' correspond to P respectivelyTotalAnd P 'totalThe path transition function of place photon path.
6. a kind of Monte Carlo scattered photon analogy method of task based access control driving according to claim 5, feature exist
In：Uniform probability density function or Gaussian probability-density function are symmetric function, so formula (5) can be reduced to：
7. a kind of Monte Carlo scattered photon analogy method of task based access control driving according to claim 6, feature exist
In：5th step records the photon path simulation number of Compton scattering and Rayleigh scattering signal respectively.
8. a kind of Monte Carlo scattered photon analogy method of task based access control driving according to claim 7, feature exist
In：In second step, corresponding sampling algorithm is chosen according to scattering point position difference and carries out random sampling operation, it is specific as follows：
P1, the position for judging scattering point；
If single order scattering point, then randomly selected using uniform sampling algorithm；
Otherwise, P2 is entered step；
P2, if multistage scattering point, then judge scattering point occur scattering type；
If happens is that Compton scattering, is randomly selected using uniform sampling algorithm；
If happens is that Rayleigh scattering, is randomly selected using RWS algorithms.
9. a kind of Monte Carlo scattered photon analogy method of task based access control driving according to claim 8, feature exist
In：For in multistage scattering point, the method using uniform sampling algorithm is specific as follows：
A scatteringangleθ is randomly selected within the scope of 360 degreeC, in scatteringangleθCOpening direction on choose a little be used as Compton
Second order dispersion point.
10. a kind of Monte Carlo scattered photon analogy method of task based access control driving according to claim 9, feature exist
In：For in multistage scattering point, the method using RWS algorithms is specific as follows：
K1, judge probability function P 'totalPlace path is final path；
If probability function P 'totalPlace path is final path, then inclined by carrying out Gauss by a small margin in current angle of scattering
It moves, obtains new angle of scattering；
Otherwise, K2 is entered step；
K2, a random scatteringangleθR, with σ2Gauss sampling is carried out for variance, obtains new scatteringangleθR', in θR' on direction
It is even to choose a little as new scattering point.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810133436.6A CN108618796A (en) | 2018-02-09 | 2018-02-09 | A kind of Monte Carlo scattered photon analogy method of task based access control driving |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810133436.6A CN108618796A (en) | 2018-02-09 | 2018-02-09 | A kind of Monte Carlo scattered photon analogy method of task based access control driving |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN108618796A true CN108618796A (en) | 2018-10-09 |
Family
ID=63706024
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810133436.6A Pending CN108618796A (en) | 2018-02-09 | 2018-02-09 | A kind of Monte Carlo scattered photon analogy method of task based access control driving |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN108618796A (en) |
Cited By (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110009707A (en) * | 2019-04-09 | 2019-07-12 | 上海联影医疗科技有限公司 | Image dispersion bearing calibration, device, computer equipment and storage medium |
| CN112763519A (en) * | 2020-12-22 | 2021-05-07 | 清华大学 | Method for simulating photon scattering by using quasi-Monte Carlo method |
| CN113361080A (en) * | 2021-05-20 | 2021-09-07 | 厦门大学 | Multilayer water photon transmission semi-analytic Monte Carlo simulation method based on GPU |
| CN113804709A (en) * | 2021-09-22 | 2021-12-17 | 清华大学 | Method for correcting CT scattering signal based on quasi-Monte Carlo and forced detection |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20100179792A1 (en) * | 2009-01-09 | 2010-07-15 | Kabushiki Kaisha Toshiba | Monte carlo simulation method, simulation apparatus, and medium storing simulation program |
| CN104027121A (en) * | 2013-03-05 | 2014-09-10 | 上海联影医疗科技有限公司 | Method for renormalization of X-ray multiple scattering simulation |
| US8983887B2 (en) * | 2011-04-11 | 2015-03-17 | Xerox Corporation | Probabilistic sampling using search trees constrained by heuristic bounds |
| CN105740573A (en) * | 2016-03-02 | 2016-07-06 | 苏州网颢信息科技有限公司 | Double-step Monte Carlo simulation method applied to radioactive ray dose computation |
-
2018
- 2018-02-09 CN CN201810133436.6A patent/CN108618796A/en active Pending
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20100179792A1 (en) * | 2009-01-09 | 2010-07-15 | Kabushiki Kaisha Toshiba | Monte carlo simulation method, simulation apparatus, and medium storing simulation program |
| US8983887B2 (en) * | 2011-04-11 | 2015-03-17 | Xerox Corporation | Probabilistic sampling using search trees constrained by heuristic bounds |
| CN104027121A (en) * | 2013-03-05 | 2014-09-10 | 上海联影医疗科技有限公司 | Method for renormalization of X-ray multiple scattering simulation |
| CN105740573A (en) * | 2016-03-02 | 2016-07-06 | 苏州网颢信息科技有限公司 | Double-step Monte Carlo simulation method applied to radioactive ray dose computation |
Non-Patent Citations (4)
| Title |
|---|
| 希莫•萨日伽: "《贝叶斯滤波与平滑》", 31 August 2015 * |
| 李婷: "《生物医学光子传输》", 30 April 2015 * |
| 许淑艳: "《蒙特卡罗方法在实验核物理中的应用》", 31 August 2006 * |
| 陈宇思: "A Metropolis Monte Carlo Simulation Scheme for Fast Scattering Calculation in Cone Beam CT", 《FULLY3D 2017 PROCEEDINGS》 * |
Cited By (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110009707A (en) * | 2019-04-09 | 2019-07-12 | 上海联影医疗科技有限公司 | Image dispersion bearing calibration, device, computer equipment and storage medium |
| CN112763519A (en) * | 2020-12-22 | 2021-05-07 | 清华大学 | Method for simulating photon scattering by using quasi-Monte Carlo method |
| CN112763519B (en) * | 2020-12-22 | 2022-03-01 | 清华大学 | Method for simulating photon scattering by using quasi-Monte Carlo method |
| CN113361080A (en) * | 2021-05-20 | 2021-09-07 | 厦门大学 | Multilayer water photon transmission semi-analytic Monte Carlo simulation method based on GPU |
| CN113361080B (en) * | 2021-05-20 | 2022-05-17 | 厦门大学 | Multilayer water photon transmission semi-analytic Monte Carlo simulation method based on GPU |
| CN113804709A (en) * | 2021-09-22 | 2021-12-17 | 清华大学 | Method for correcting CT scattering signal based on quasi-Monte Carlo and forced detection |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Schuemann et al. | Site-specific range uncertainties caused by dose calculation algorithms for proton therapy | |
| Park et al. | Proton dose calculation on scatter‐corrected CBCT image: feasibility study for adaptive proton therapy | |
| Doolan et al. | Patient-specific stopping power calibration for proton therapy planning based on single-detector proton radiography | |
| De Jaeger et al. | Incorporating an improved dose-calculation algorithm in conformal radiotherapy of lung cancer: re-evaluation of dose in normal lung tissue | |
| Papagiannis et al. | Current state of the art brachytherapy treatment planning dosimetry algorithms | |
| Giantsoudi et al. | Validation of a GPU-based Monte Carlo code (gPMC) for proton radiation therapy: clinical cases study | |
| Bauer et al. | Integration and evaluation of automated Monte Carlo simulations in the clinical practice of scanned proton and carbon ion beam therapy | |
| CN108618796A (en) | A kind of Monte Carlo scattered photon analogy method of task based access control driving | |
| Plautz et al. | An evaluation of spatial resolution of a prototype proton CT scanner | |
| Botta et al. | Use of the FLUKA Monte Carlo code for 3D patient-specific dosimetry on PET-CT and SPECT-CT images | |
| Meyer et al. | Dosimetric accuracy and radiobiological implications of ion computed tomography for proton therapy treatment planning | |
| McCowan et al. | An in vivo dose verification method for SBRT–VMAT delivery using the EPID | |
| Cheng et al. | Improved dose–volume histogram estimates for radiopharmaceutical therapy by optimizing quantitative SPECT reconstruction parameters | |
| Penfold et al. | A more accurate reconstruction system matrix for quantitative proton computed tomography | |
| Dickmann et al. | Prediction of image noise contributions in proton computed tomography and comparison to measurements | |
| Penfold | Image reconstruction and Monte Carlo simulations in the development of proton computed tomography for applications in proton radiation therapy | |
| Schultze et al. | Particle-tracking proton computed tomography—data acquisition, preprocessing, and preconditioning | |
| Berndt et al. | Application of single-and dual-energy CT brain tissue segmentation to PET monitoring of proton therapy | |
| Yepes et al. | Monte Carlo fast dose calculator for proton radiotherapy: application to a voxelized geometry representing a patient with prostate cancer | |
| Schmid et al. | Monte Carlo study on the sensitivity of prompt gamma imaging to proton range variations due to interfractional changes in prostate cancer patients | |
| Maneval et al. | pGPUMCD: an efficient GPU-based Monte Carlo code for accurate proton dose calculations | |
| Cassetta et al. | Accuracy of low‐dose proton CT image registration for pretreatment alignment verification in reference to planning proton CT | |
| Collins Fekete et al. | Quantifying the effect of seed orientation in postplanning dosimetry of low‐dose‐rate prostate brachytherapy | |
| Pastor-Serrano et al. | How should we model and evaluate breathing interplay effects in IMPT? | |
| Van der Heyden et al. | Artificial intelligence supported single detector multi-energy proton radiography system |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| RJ01 | Rejection of invention patent application after publication |
Application publication date: 20181009 |
|
| RJ01 | Rejection of invention patent application after publication |