CN112763519B - Method for simulating photon scattering by using quasi-Monte Carlo method - Google Patents

Method for simulating photon scattering by using quasi-Monte Carlo method Download PDF

Info

Publication number
CN112763519B
CN112763519B CN202011526102.9A CN202011526102A CN112763519B CN 112763519 B CN112763519 B CN 112763519B CN 202011526102 A CN202011526102 A CN 202011526102A CN 112763519 B CN112763519 B CN 112763519B
Authority
CN
China
Prior art keywords
scattering
photon
photons
probability
detector pixel
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
Application number
CN202011526102.9A
Other languages
Chinese (zh)
Other versions
CN112763519A (en
Inventor
林桂元
邓世沃
王小群
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Tsinghua University
Capital Normal University
Original Assignee
Tsinghua University
Capital Normal University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Tsinghua University, Capital Normal University filed Critical Tsinghua University
Priority to CN202011526102.9A priority Critical patent/CN112763519B/en
Publication of CN112763519A publication Critical patent/CN112763519A/en
Application granted granted Critical
Publication of CN112763519B publication Critical patent/CN112763519B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N23/00—Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00
    • G01N23/02—Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by transmitting the radiation through the material
    • G01N23/04—Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by transmitting the radiation through the material and forming images of the material
    • G01N23/046—Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by transmitting the radiation through the material and forming images of the material using tomography, e.g. 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/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

Landscapes

  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Pathology (AREA)
  • General Health & Medical Sciences (AREA)
  • Medical Informatics (AREA)
  • Physics & Mathematics (AREA)
  • Radiology & Medical Imaging (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Biophysics (AREA)
  • Theoretical Computer Science (AREA)
  • Analytical Chemistry (AREA)
  • General Physics & Mathematics (AREA)
  • Immunology (AREA)
  • Chemical & Material Sciences (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Pulmonology (AREA)
  • High Energy & Nuclear Physics (AREA)
  • Biochemistry (AREA)
  • Optics & Photonics (AREA)
  • Biomedical Technology (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Molecular Biology (AREA)
  • Surgery (AREA)
  • Animal Behavior & Ethology (AREA)
  • Public Health (AREA)
  • Veterinary Medicine (AREA)
  • Nuclear Medicine (AREA)

Abstract

The invention discloses a method for simulating photon scattering by using a quasi-Monte Carlo method, which is applied to a system comprising a light source S, a die body M and a detector D, wherein preset photons are emitted from the light source S, and are preset at an intersection point A0Initially into the phantom M at the interaction point AiAnd (i-1.. n.) the i-order scattering occurs, the method tracks the distribution of photons from the light source S to the measured phantom M, simulates the scattering path of the photons in the phantom M, and forcibly selects the detector pixel corresponding to the interaction point, thereby realizing the simulation of the photon scattering process.

Description

Method for simulating photon scattering by using quasi-Monte Carlo method
Technical Field
The invention relates to the field of scattered photon simulation, in particular to a method for simulating photon scattering by using a quasi-Monte Carlo method.
Background
Cone-beam Computed Tomography (CBCT) has been widely used in the fields of clinical medicine, industry, security, and the like. However, photons passing through the phantom will scatter, and both scattered and unscattered signals (useful signals) will be collected by the detector, which will affect the useful signal, especially when using a flat panel detector and collecting under high-energy X-rays. This is one of the major challenges in obtaining high quality CBCT images. If the photon scattering problem is not considered, the contrast resolution of the reconstructed image is reduced, cup-shaped artifacts, shadows, stripes and inhomogeneities are generated, and the accuracy of the CBCT image is greatly influenced. Therefore, it is an important topic to study how to eliminate photon scattering.
Monte Carlo (MC) simulation methods proposed in the 40 s of the 20 th century were used to solve numerical problems related to the diffusion of random neutrons in fissile materials in atomic bomb design. In the last two thirty years, the MC method has played an important role in solving the problem of the application of the radiation medicine physics, with its extremely high flexibility and powerful functions. Many MC tools have been developed to study photon scattering phenomena, such as PENELOPE, MC-GPU, gMCDRR, gMMC, and the like. Currently, photon scattering simulation methods commonly used in China include a Monte Carlo-graphics processing unit (MC-GPU) method, a GPU-based sink MC (gMCDRR) method and a GPU-based Metropolis MC (gMMC) method. However, these photon scattering simulation methods have a problem that the scattered photon calculation time is too long.
Accordingly, it is desirable to provide a method for simulating photon scattering using a monte carlo-like method.
Disclosure of Invention
It is an object of the present invention to overcome or at least mitigate the above-mentioned disadvantages of the prior art by providing a method of simulating photon scattering using a quasi-monte carlo method.
In order to achieve the above object, the present invention provides a method for simulating photon scattering by using a quasi-Monte Carlo method, which is applied to a system comprising a light source S, a phantom M and a detector D, wherein a preset photon is emitted from the light source S, and is positioned at an intersection point A0Initially into the phantom M at the interaction point AiI-order scattering occurs, (i 1.. multidot.n) and the intersection point a is selected0Corresponding detector pixel
Figure GDA0003307176460000021
And are provided with
Figure GDA0003307176460000022
Selecting k-1 and A in a preset domain range as the centeriCorresponding detector pixel
Figure GDA0003307176460000023
j is 1, … k-1, and the photons reach the detector pixel after i-order scattering
Figure GDA0003307176460000024
The path of j-0, …, k-1 is noted
Figure GDA0003307176460000025
Wherein B isiIndicating the edge of the photon after i-order scattering
Figure GDA0003307176460000026
The intersection point of the direction running out of the motif M and the boundary of M comprises the following steps:
step 1, giving a specific numerical value of a scattering order n to generate a (6n +1) -dimensional Sobol point u1,…,u6n+1；
Step 2, using the 1 st component u of the (6n +1) -dimensional Sobol sequence1Sampling the initial energy E of incident photon from a preset energy spectrum phi by an inverse function sampling method0(ii) a Using the 2 nd and 3 rd components u of the (6n +1) -dimensional Sobol sequence2、u3Solid angle region omega for incident light0The incident direction and the intersection point A of the photons are obtained by uniform sampling0From A) in accordance with the following equation 1) to calculate that the photons are not scattered and0edge of
Figure GDA0003307176460000027
Probability p of direction passing through motif0；
Figure GDA0003307176460000028
Wherein the content of the first and second substances,
Figure GDA0003307176460000029
denotes the initial incident direction of the photon, μt(x，E0) Representing photon energy as E0A linear attenuation coefficient space coordinate function, wherein the total attenuation coefficient is the sum of a scattering coefficient and a photoelectric absorption coefficient;
and is provided with W1＝W0*(1-p0)，W1Indicating the initial direction of incidence of the photons
Figure GDA00033071764600000210
No scatter stays at the weight of the phantom M; w0＝fl(u2,u3) Wherein f islIs the photon initial incident direction probability density;
step 3, random number u is used according to the following formula 2)4Sampling first order interaction point A1：
Figure GDA00033071764600000211
Wherein u is4Is the 4 th component of the (6n +1) -dimensional Sobol sequence and is represented by a random number u6,u7Selection of and A1Corresponding detector pixel
Figure GDA00033071764600000212
And is arranged at
Figure GDA00033071764600000213
Within a preset domain of k-1 pixels
Figure GDA00033071764600000214
j is 1, …, k-1 and A1Correspondingly, the photon path is calculated according to
Figure GDA00033071764600000215
j ∈ {0, …, k-1} reaches the detector pixel
Figure GDA00033071764600000216
The weight of (c):
Figure GDA00033071764600000217
j＝0,…k-1
wherein
Figure GDA00033071764600000218
Step 4, if 1 is less than n, the photons continue to generate second-order scattering; let the photon be at the interaction point Ai-1I is 2, …, the sample probability of rayleigh scattering for n +1 is γ, and the (6(i-1) +5) th component u of the (6n +1) -dimensional Sobol sequence is used6(i-1)+5If u is6(i-1)+5If < gamma, it is judged as being Ai-1Rayleigh scattering occurs at the position; otherwise, it is determined as Ai-1Is treated with ComptonScattering; is provided with
Figure GDA0003307176460000031
Figure GDA0003307176460000032
Wherein p isc(Ai-1)、pr(Ai-1) And pa(Ai-1) Respectively representing the photons at the next interaction point Ai-1The proportion of compton scattering, rayleigh scattering and photoelectric absorption occurring;
step 5, sampling scattering angle direction: the (6(i-1) +2) th component u of the (6n +1) dimensional Sobol sequence is used in the RITA algorithm6(i-1)+2To sample the cosine of the scatter polar angle cos (theta)i-1) And using the (6(i-1) +3) th component u of the (6n +1) dimensional Sobol sequence6(i-1)+3Obtaining an azimuthal scattering angle phii＝2π*u6(i-1)+3Then calculating new scattering direction
Figure GDA0003307176460000033
Then calculates the (i-1) order scattering of the photon and then calculates the (i-1) order scattering of the photon at the interaction point Ai-1Edge of
Figure GDA0003307176460000034
Probability of direction running out of motif M
Figure GDA0003307176460000035
And is provided with Wi＝Wi-1*(1-pi-1)；
Step 6) in the direction according to the following formula 3)
Figure GDA0003307176460000036
Upper random number u6(i-1)+4Sampling the next order interaction point Ai,i＝2,…,n：
Figure GDA0003307176460000037
Wherein u is6(i-1)+4Is (6n +1) dimensionThe (6(i-1) +4) th component of the Sobol sequence;
step 7, using the 6i and 6i +1 components u of the (6n +1) dimensional Sobol sequence6iAnd u6i+1Selecting and interacting with a point AiN, i-2, …
Figure GDA0003307176460000038
Then
Figure GDA0003307176460000039
Selecting k-1 interaction points A within the preset field rangeiCorresponding detector pixel
Figure GDA00033071764600000310
j is 1, … k-1, and the photon is calculated at the interaction point AiProbability density of scattering angle at which Compton scattering occurs
Figure GDA00033071764600000311
Probability density of scattering angle at which Rayleigh scattering occurs
Figure GDA00033071764600000312
And solid angle
Figure GDA00033071764600000313
The probability density of the scattering angle is the normalized scattering differential cross-section DCS, so
Figure GDA00033071764600000314
Wherein muc(Ai,Ei) And mup(Ai,Ei) Respectively representing photons at the interaction point AiAt an energy of EiCompton and rayleigh differential scattering cross sections;
Figure GDA00033071764600000315
j is 0, …, k-1; calculating the photon from A byiStarting along the direction
Figure GDA00033071764600000316
Probability of running out of motif M:
Figure GDA00033071764600000317
calculating photon path by
Figure GDA0003307176460000041
j ∈ {0, …, k-1} reaches the detector pixel
Figure GDA0003307176460000042
The weight of (c):
Figure GDA0003307176460000043
i＝2,…n,j＝0,…k-1；
step 8, if i is less than n, continuing to scatter the photons for the next time, and returning to the step 4 if i is equal to i + 1;
step 9, if i is equal to n, returning
Figure GDA0003307176460000044
A value of (d);
step 10, calculating the detected scattered signal
Figure GDA0003307176460000045
i-1, …, n, j-1, …, k-1, wherein
Figure GDA0003307176460000046
For the first photon along the path
Figure GDA0003307176460000047
Reach the detector pixel
Figure GDA0003307176460000048
The weight of (a) is determined,
Figure GDA0003307176460000049
is a function of the detector response, which, in relation to the path,
Figure GDA00033071764600000410
and N is the number of analog photons.
Preferably, the method further comprises:
and displaying the result of simulating photon scattering according to the detected scattering signal.
Preferably, after the photons enter the die body M, the motion paths of the photons are distributed
Figure GDA00033071764600000411
Preferably, before step 1, the method further comprises: calculating the probability of a path of a photon scattered in a measured phantom M and the probability of the photon reaching a detector pixel from an interaction point after scattering by using a high-dimensional integration mode, wherein the probability comprises the following steps:
calculating the photon from A according to0From the beginning to the edge
Figure GDA00033071764600000412
Directional motion, with a probability of staying at the phantom M of
Figure GDA00033071764600000413
Wherein omega0Representing the solid angle of incident light, s representing a point in the solid angle of incident light, flIs the probability density of the initial direction of incidence of the photon,
Figure GDA00033071764600000414
is the value space of the energy spectrum, and phi (E) is the spectral probability density function.
Preferably, the calculating the probability of the path of the photon scattering in the phantom M by using a high-dimensional integration method includes:
calculating the scattering of the photon from A after the photon generates i-1, i epsilon {2, …, n } order according to the following formulai-1Starting in the direction
Figure GDA00033071764600000415
Motion, probability of staying within phantom M:
Figure GDA00033071764600000416
wherein
Figure GDA0003307176460000051
Is ai-1Scattering polar angle theta at apex angleδThe solid angle of (a) is,
Figure GDA0003307176460000052
i-2, …, n denotes the photon interaction point ai-1The energy after compton scattering has occurred,
Figure GDA0003307176460000053
i-2, …, n denotes the photon interaction point ai-1The energy after rayleigh scattering.
Preferably, the i-order scattering of the photons is calculated by high-dimensional integration according to the following formula, and then the photons are transmitted from the interaction point AiI e { 1.. multidot.n } runs out of the phantom M to reach the detector pixel
Figure GDA00033071764600000518
Probability of j ═ 0, …, k-1:
Figure GDA0003307176460000054
wherein the content of the first and second substances,
Figure GDA0003307176460000055
is the scattering point AiAnd corresponding detector pixel
Figure GDA0003307176460000056
The solid angle of (1) is given by
Figure GDA0003307176460000057
And is
Figure GDA0003307176460000058
Wherein
Figure GDA0003307176460000059
Is a detector pixel
Figure GDA00033071764600000510
The length of (a) of (b),
Figure GDA00033071764600000511
to represent
Figure GDA00033071764600000512
And a detector pixel
Figure GDA00033071764600000513
The angle in the normal direction.
Preferably, the photon along path is calculated by high-dimensional integration according to the following formula
Figure GDA00033071764600000514
Reaches a detector pixel after i-order scattering
Figure GDA00033071764600000515
Probability of (c):
Figure GDA00033071764600000516
due to the adoption of the technical scheme, the invention has the following advantages:
the method for simulating photon scattering by using the quasi-Monte Carlo method is provided, the distribution of photons reaching a measured die body M from a light source S is tracked, the scattering path of the photons in the die body M is simulated, and detector pixels corresponding to interaction points are forcibly selected, so that the simulation of the photon scattering process is realized.
Drawings
FIG. 1 is a diagram of a model for simulating scattered photons in an embodiment of the present invention.
Fig. 2 is a schematic diagram of a scattering path of a photon in the model established in fig. 1 according to an embodiment of the present invention.
FIG. 3 is a diagram illustrating a path obtained by high-dimensional integration in a method for simulating photon scattering by a Monte Carlo simulation method according to an embodiment of the present invention
Figure GDA00033071764600000517
Schematic diagram of the probability integral expression of (1).
Fig. 4 is a schematic flow chart of a method for simulating photon scattering by using a quasi-monte carlo method according to an embodiment of the present invention.
Fig. 5 is a schematic diagram of a preset energy spectrum provided by an embodiment of the invention.
FIG. 6 shows a schematic diagram for generating azimuthal scattering angles provided by an embodiment of the present invention.
FIG. 7(a) shows DCS of Rayleigh scattering of bone material in an energy range of 5-150 keV.
Fig. 7(b) and (c) are interpolation errors between PDF by RITA interpolation and the Generalized Inverse Transformation Method (GITM) and the original rayleigh hdcs.
FIG. 8(a) shows DCS for Compton scattering of bone material in the energy range of 5-150 keV.
Fig. 8(b) and (c) are interpolation errors between PDF by RITA interpolation and the Generalized Inverse Transformation Method (GITM) and the original rayleigh hdcs.
FIGS. 9(a) and (a1) show the Shepp-Logan motif and its geometric schematic.
FIGS. 9(b), (c) and (d) show transverse, coronal and sagittal sections, respectively, of the Shepp-Logan motif.
FIG. 10 shows the scatter signals of 5-120keV cone beam X-rays through the Shepp-Logan phantom.
Detailed Description
The invention is described in detail below with reference to the figures and examples.
The embodiment of the invention provides a method for simulating photon scattering by using a quasi-Monte Carlo method, which is applied to a scattered photon simulation model shown in figure 1.
Fig. 1 shows a scattered photon simulation model comprising a light source S, a phantom M and a detector D. FIG. 1 shows a geometric model of a photon imaging system including a light source S, modeA volume M and a detector D. The geometrical illustration of the third-order scattering is also given exemplarily in fig. 1. Photons passing through the phantom M undergo third-order scattering, A1,A2,A3Respectively representing 1,2,3 order interaction points. Photons arrive at each interaction AiAfter i is 1,2,3, the direction changes and with a certain probability passes through the phantom M and reaches the detector D, and a1,A2,A3The pixels of the detectors corresponding to each other are D1,D2,D3. Also illustratively shown in FIG. 1 is a three-dimensional coordinate system established within the phantom.
Fig. 2 shows an exemplary illustration of the photon scattering path in the model established in fig. 1. Wherein the initial energy is considered to be E0The photons of (a) pass through the path l of the three-dimensional irregular phantom M. As shown in FIG. 2, photons emerge from the light source S at an intersection A0Initially into the phantom M. The light source is an X-ray source. After entering the mold body M, photons are scattered or absorbed. After some scattering, the photons reach the detector or escape with a certain probability. In this example, the point at which the photons are scattered is called the interaction point, denoted Ai(i 1.., n), where i represents the order in which the photons are scattered. First, forced random selection of a detector pixel corresponding to each interaction point, and AiThe corresponding detector pixel is noted
Figure GDA0003307176460000071
Then
Figure GDA0003307176460000072
Select k-1 detector pixels within a predetermined range of fields
Figure GDA0003307176460000073
j is 1, …, k-1 and AiAnd (7) corresponding. The photons reach the detector pixel after i-order scattering
Figure GDA0003307176460000074
Is marked as
Figure GDA0003307176460000075
Wherein
Figure GDA0003307176460000076
Indicating the edge of the photon after i-order scattering
Figure GDA0003307176460000077
The direction runs out of the motif M and intersects the boundary of M. C in FIG. 2i-1And (i 1.. multidot.n) represents the trailing edge of the photon after i-order scattering
Figure GDA0003307176460000078
The intersection of the directional motion with the border of the phantom M. For the convenience of description, herein, will be described
Figure GDA0003307176460000079
The direction is called the probing direction, and
Figure GDA00033071764600000710
the direction is called the scattering direction.
The embodiment of the invention provides a method for simulating photon scattering by using a quasi-Monte Carlo method, which is applied to the systems shown in figures 1 and 2. In the method provided by the invention, the abstract particle transport path problem can be firstly converted into a high-dimensional integral problem, and the path is given
Figure GDA00033071764600000711
The probability integral expression of (2) and a simulation method based on a quasi-Monte Carlo method are provided.
As shown in fig. 3, for the path
Figure GDA00033071764600000712
The calculation of the integral expression of (a) may include the following steps S310 to S350.
S310, the probability of the photon initialization energy and the initial incidence direction and the probability of escaping from the phantom without scattering are calculated.
S320, sampling and calculating the position A of the 1 st order interaction point in the initial incidence direction1。
S330, sampling photons at an interaction point Ai-1I 2, …, scattering type of n +1, scattering angle distribution, in
Figure GDA00033071764600000713
Directionally sampling the next order interaction point AiAnd calculating photon passing interaction point Ai-1To AiThe probability of (c).
S340, forced random selection and AiCorresponding detector pixel
Figure GDA00033071764600000714
And is arranged at
Figure GDA00033071764600000715
Within a preset field range, to select k-1 detector pixels
Figure GDA00033071764600000716
j is 1, …, k-1 and AiAnd (7) corresponding.
S350, calculating the photon from AiTo the detector pixel
Figure GDA00033071764600000717
j is 0, …, k-1, and accumulates; the scattering matrix is recorded as is (D),
Figure GDA00033071764600000718
wherein
Figure GDA00033071764600000719
For the first photon along the path
Figure GDA00033071764600000720
Reach the detector pixel
Figure GDA00033071764600000721
The probability of (a) of (b) being,
Figure GDA00033071764600000722
is a function of the detector response, which, in relation to the path,
Figure GDA00033071764600000723
and N is the number of analog photons.
Specific illustrative examples of the above steps are as follows.
In S310, after the photons enter the die body M, the motion paths of the photons obey distribution
Figure GDA0003307176460000081
Wherein mut(x,E0) A function representing the linear spatial coordinates of the attenuation coefficient of the photon, E0As an initial value of photon energy, A1Is the incident direction of photons
Figure GDA0003307176460000082
Random point on. So that photons are from A0From the beginning to the edge
Figure GDA0003307176460000083
The direction of movement, the probability of staying at the phantom M being (i.e. the direction of incidence of the photons from the source S)
Figure GDA0003307176460000084
Move to A1Probability of (1)
Figure GDA0003307176460000085
Wherein omega0Is the solid angle of incident light, flIs the photon initial incident direction probability density, s represents the point in the solid angle of incident light,
Figure GDA0003307176460000086
indicating the direction of incidence in which the photons are not scattered,
Figure GDA0003307176460000087
representing photons from the point of incidence A0Move to interaction point a1，μt(x,E0) Representing photon energy as E0Linear space of time attenuation coefficientTarget function, E0Is the initial energy of the incident light.
So that the photons are not scattered from A0From the beginning to the edge
Figure GDA0003307176460000088
Directional motion with a probability of escaping from the phantom M of
Figure GDA0003307176460000089
Note the book
Figure GDA00033071764600000810
Representing photons from A0From the beginning to the edge
Figure GDA00033071764600000811
Probability of direction escaping the motif M.
In S320, in the initial incidence direction, the 1 st order interaction point position A is sampled and calculated1。
Firstly according to the formula
Figure GDA00033071764600000812
First order interaction point position A of sampled photons1Wherein u is4Is the 4 th component of the (6n +1) -dimensional Sobol sequence. After the photons generate i-1, i-2, …, n +1 order scattering, the photons are sampled at the interaction point AiThe direction of movement and the length of movement to determine the next interaction point location.
In S330, the sampled photon is at interaction point Ai-1I 2, …, n +1, scattering type, scattering angle distribution, and according to the formula
Figure GDA00033071764600000813
In that
Figure GDA0003307176460000091
Directionally up-sampling next order interactionUsing point AiWherein u is6(i-1)+4Is the 6(i-1) +4 th component of the (6n +1) -dimensional Sobol sequence. Finally calculating photon interaction point Ai-1To AiThe probability of (c).
Photon is in A1In the direction of motion of photons at A1Type of interaction of dots Tδδ ∈ {0,1} and the scattering angle distribution. Interaction type Tδδ ∈ {0,1} includes T0And T1Compton scattering and rayleigh scattering, respectively. The probability density of the scattering angle is a normalized DCS (differential scattering cross section) function.
Figure GDA0003307176460000092
Respectively representing photons at the interaction point A1Probability of compton scattering and probability density of scattering angle,
Figure GDA0003307176460000093
respectively representing photons at the interaction point A1Probability of occurrence of rayleigh scattering and probability density of scattering angle.
After the first-order scattering of the photons occurs, the motion path of the photons in the die body M follows the distribution:
Figure GDA0003307176460000094
wherein E1Is the energy, mu, remaining after 1 st order scattering of a photont(x,E1) Representing photon energy as E1Function of linear spatial coordinates of the attenuation coefficient of time, A2Is that
Figure GDA0003307176460000095
Random point on. If Rayleigh scattering occurs to the photon, the energy does not changei＝Ei-1. If the photon is Compton scattered, the energy E after each scatteringi(i ═ 1.., n) is changed,
Figure GDA0003307176460000096
wherein m isec2Representing electron mass energy, EiN represents the number i of scattered photons,at interaction point AiEnergy of thetaiIs shown at interaction point AiThe scattering angle of (c).
So that photons are from A1From the beginning to the edge
Figure GDA0003307176460000097
The direction of movement, the probability of staying in the motif M being complementary (i.e. photons at A)1After scattering occurs, the edge
Figure GDA0003307176460000098
Direction of movement to A2Probability of (1)
Figure GDA0003307176460000099
Wherein
Figure GDA00033071764600000910
Is a1Scattering angle theta at apex angleδThe solid angle of (1).
Figure GDA00033071764600000911
Respectively representing photons at the interaction point A1Probability of compton scattering and probability density of scattering angle,
Figure GDA00033071764600000912
respectively representing photons at the interaction point A1The probability of occurrence of rayleigh scattering and the probability density of scattering angles,
Figure GDA00033071764600000913
i is 1, …, n represents photon interaction point AiThe energy after compton scattering has occurred,
Figure GDA00033071764600000914
i is 1, …, n represents photon interaction point AiThe energy after rayleigh scattering. For example,
Figure GDA00033071764600000915
indicates the photon interaction point A1The energy after compton scattering has occurred,
Figure GDA0003307176460000101
indicates the photon interaction point A1The energy after rayleigh scattering.
In the same way, after i-1, i-E { 2.,. n } order scattering occurs to the photon, the photon is scattered from Ai-1Starting in the direction
Figure GDA0003307176460000102
The probability of motion, staying within the motif M, is (i.e. the photon is at A)i-1After the i-1 th scattering occurs, the edge
Figure GDA0003307176460000103
Direction of movement to AiProbability of (1)
Figure GDA0003307176460000104
Wherein
Figure GDA0003307176460000105
Is ai-1Theta at the apex angleδThe solid angle of (a) is,
Figure GDA0003307176460000106
i-2, …, n denotes the photon interaction point ai-1The energy after compton scattering has occurred,
Figure GDA0003307176460000107
i-2, …, n denotes the photon interaction point ai-1The energy after rayleigh scattering.
In S340, force random selection and AiCorresponding detector pixel
Figure GDA0003307176460000108
And is arranged at
Figure GDA0003307176460000109
In the preset field of the image processing system, k-1 detector pixels are selected
Figure GDA00033071764600001010
j is 1, …, k-1 and AiAnd (7) corresponding.
In S350, the photon from A is calculatediTo the detector pixel
Figure GDA00033071764600001011
j equals 0, …, k-1, and accumulates.
Wherein photons are emitted from the interaction point AiI e {1, …, n } runs out of the phantom M to the detector pixel
Figure GDA00033071764600001012
I.e. the passing path
Figure GDA00033071764600001013
Has a probability of
Figure GDA00033071764600001014
Wherein
Figure GDA00033071764600001015
Is the scattering point AiAnd corresponding detector pixel
Figure GDA00033071764600001016
The solid angle of (1) is given by
Figure GDA00033071764600001017
And is
Figure GDA00033071764600001018
Wherein
Figure GDA00033071764600001019
To represent
Figure GDA00033071764600001020
And a detector pixel
Figure GDA00033071764600001021
The angle of the normal direction is such that,
Figure GDA00033071764600001022
is a detector pixel
Figure GDA00033071764600001023
The length of (a) of (b),
Figure GDA00033071764600001024
is an interaction point AiTo the corresponding detector pixel
Figure GDA00033071764600001025
The distance of (c).
So that photons follow the path
Figure GDA00033071764600001026
The probability of i ∈ {1, …, n }, j ∈ { 0., k-1} is
Figure GDA00033071764600001027
Wherein A is0,A1,...,Ai(i 1.. n.) is a random point, related to the initial energy and initial incident direction, and aiAnd Ai-1In connection with this, the present invention is,
Figure GDA0003307176460000111
for forced selection of andia corresponding detector pixel.
The steps shown in fig. 3 above convert the abstract particle transport path problem into a high-dimensional integration problem, thereby providing a new scattered photon simulation concept. However, it is often difficult to directly compute the high-dimensional integral, and a new method based on quasi-monte carlo simulation, gQMC, is proposed herein to compute the high-dimensional integral problem. The method can simultaneously calculate the path probability of different scattering orders in one simulation process, and the traditional simulation method can only calculate the path probability of a certain order of scattering in one simulation process.
Fig. 4 is a schematic flow chart of a method for simulating photon scattering by using a quasi-monte carlo method according to an embodiment of the present invention.
As shown in FIG. 4, steps S410-S4100 are included. Wherein the photon path is related to the photon initial energy, the scatter angle distribution, the scatter order, the location of the interaction, the scatter interaction type, and the location of the detector. And randomly selecting a detector corresponding to the interaction point in advance in the simulation process.
Step S410, giving a specific numerical value of the scattering order n to generate a (6n +1) -dimensional Sobol point u1,...,u6n+1。
The specific value of n is determined by the user according to actual needs, and the specific value is not limited in this document.
The generation of the Sobol point is a known technology, and the meaning thereof is not described in detail herein.
Step S420, using random number u from preset energy spectrum1Sampling photon energy E0Using u2、u3Solid angle region omega for incident light0Uniformly sampling to obtain the intersection point A of the incident direction of photons and the entering die body0And calculating that no scattering occurs in the photons according to the following formula 8), from A0From the beginning to the edge
Figure GDA0003307176460000112
Weight p of direction through motif0，
Figure GDA0003307176460000113
Wherein the content of the first and second substances,
Figure GDA0003307176460000114
indicating the direction of incidence of the photons without scattering; mu.st(x,E0) Representing photon energy as E0The linear space coordinate function of the time attenuation coefficient, the value of which is obtained by looking up the table; in this context makeThe PENELOPE physical database in Geant4(GEometryANd Tracking) was used;
and is provided with W1＝W0*(1-p0)，W1Indicating the initial direction of incidence of the photons
Figure GDA0003307176460000115
No scatter stays at the weight of the phantom M; w0＝fl(u2,u3) Wherein f islIs the photon initial incidence direction probability density.
Fig. 5 shows a schematic diagram of a preset energy spectrum. Wherein u may be used1From the spectral probability density function distribution phi (E)0) Middle sampling E0。
Step S430, using random number u4In the initial incident direction
Figure GDA0003307176460000121
Sampling first order interaction point A1And force random selection of interaction point A1Corresponding detector pixels and calculating the weight of the photons reaching the detector.
Since the photons are in the motif M, from A0→A1Obeying the distribution:
Figure GDA0003307176460000122
the interaction point A may be sampled by an inverse function sampling method1. Therefore A1The specific position of (c) can be given by the following equation:
Figure GDA0003307176460000123
wherein u is4Is the 4 th component of the (6n +1) -dimensional Sobol sequence,
Figure GDA0003307176460000124
the photons are not scattered and are incident along the direction
Figure GDA0003307176460000125
Probability of escaping motifs M.
And using random number u6,u7Forced random selection with A1Corresponding detector pixel
Figure GDA0003307176460000126
Then
Figure GDA0003307176460000127
In the preset field, forcibly prepare k-1 pixels
Figure GDA0003307176460000128
j is 1, …, k-1 and A1And (7) corresponding. So that photons follow the path
Figure GDA0003307176460000129
j ∈ {0, …, k-1} reaches the detector pixel
Figure GDA00033071764600001210
The weight of (A) is:
Figure GDA00033071764600001211
wherein
Figure GDA00033071764600001212
Step S440 determines the weight occupied by the scattering type and the corresponding type.
If 1 is less than n, the photon continues to generate second-order scattering; let the photon be at the interaction point Ai-1I is 2, …, the sample probability of rayleigh scattering for n +1 is γ, and the (6(i-1) +5) th component u of the (6n +1) -dimensional Sobol sequence is used6(i-1)+5If u is6(i-1)+5If < gamma, it is judged as being Ai-1Rayleigh scattering occurs at the position; otherwise, it is determined as Ai-1Compton scattering occurs; is provided with
Figure GDA00033071764600001213
Figure GDA00033071764600001214
Wherein p isc(Ai-1)，pr(Ai-1) A partial sum of pa(Ai-1) Respectively representing the photons at the next interaction point Ai-1The proportion of compton scattering, rayleigh scattering and photoelectric absorption that occurs.
In step S450, the scattering angle direction is sampled.
After (i-1) -order scattering of a photon, i belongs to { 2., n +1}, and the probability density function p (cos (theta) of the scattering angle of Rayleigh scatteringi-1) Is given by the following equation:
Figure GDA0003307176460000131
wherein A is a normalization constant, Ei-2Is the energy of the photon after (i-2) -order scattering, c is the constant of the speed of light, and F is the atomic structural factor.
Furthermore, the scattering angle probability density distribution of compton scattering is described herein using Klein-Nishina differential cross-sections:
Figure GDA0003307176460000132
wherein B is a normalization constant, Ei-1Is the energy remaining after (i-1) -order scattering of a photon, and G is the atomic scattering function.
Rayleigh scattering of photons with constant energy Ei＝Ei-1. Compton scattering of photons, energy E after each scatteringi(i ═ 1.., n) is changed,
Figure GDA0003307176460000133
wherein m isec2Representing electron mass energy, EiWhere i is 1, …, n indicates that the photon has been scattered the ith time and then at the interaction point AiThe energy of (a). ThetaiIs shown at interaction point AiThe scattering angle of (c).
Photon generation (i-1) order powderAfter irradiation, Polar Scattering Angle θ of Rayleigh Scattering and Compton Scattering was achieved using a Rational Inverse Transform with Aliasing (RITA) algorithmi-1Is sampled.
Using the (6(i-1) +2) th component u of the (6n +1) -dimensional Sobol sequence6(i-1)+2Extracting cos (theta) from DCSi-1) And obtaining the interval probability density function of the cosine of the scattering angle. From the uniform distribution in (0,2 π), an Azimuthal Scattering Angle (Azimuthal Scattering Angle) φ is generatedi-1Therefore, it is
Figure GDA0003307176460000134
As shown in fig. 6.
In an implementation, the DCS pre-computes the data used in RITA. At each energy, the interpolation error between the interpolated PDF and the original DCS is minimized using the appropriate nodes, here 64 nodes.
For Rayleigh scattering, in determining phii-1Thereafter, the DCS data can be directly calculated by equation (10).
For compton scattering, the DCS data cannot be obtained directly from equation (11) because the specific structure of the atomic scattering function G is not known. Sampling Compton scattering by utilizing Geant4 software, sampling each specific material in a given energy range to obtain a scattering angle, carrying out statistics to obtain an angle histogram, and obtaining the Compton scattering DCS data based on the angle histogram.
To illustrate the accuracy of the RITA sampling method, at the kth grid xkIntroducing interpolation error, which is defined as
Figure GDA0003307176460000141
Table 1 shows the rayleigh scatter angular distribution sampling error of the bone.
Figure GDA0003307176460000142
TABLE 1
As can be seen from Table 1, for Rayleigh scattering angles, the generalized quasi-transform method (GIT) is applied at an energy of 30keVM) method using 4096 grid points, the average error is 1.15 × 10-3While the RITA method only needs to use 64 grid points, the average error is 3.18 × 10-4. This indicates that RITA has higher accuracy and speed than GITM. Furthermore, when single precision floating point numbers are used, the GITM effect is poor, and the interpolation error increases as the grid point increases. In this section, when the rayleigh scattering angle is sampled by GITM, a double precision floating point number will be used.
FIG. 7(a) shows DCS of Rayleigh scattering of bone material in an energy range of 5 to 150 keV. Fig. 7(b) and (c) are interpolation errors between PDF by RITA interpolation and Generalized Inverse Transformation Method (GITM) and original Rayleigh DCS sampled from Rayleigh DCS with energies of 30keV and 60keV, respectively, and the interpolation grid point number of GITM is 2048. The above data show that using RITA sampling scatter angles, the simulation speed and accuracy can be optimized.
Table 2 shows the compton scattering angle distribution sampling error of the bone.
Figure GDA0003307176460000151
TABLE 2
As can be seen from Table 2, for the Compton scattering angle, the average error is 1.08X 10 when 4096 grid points are used for the GITM method at an energy of 60keV-4While the RITA method only needs to use 64 grid points, with an average error of 1.57 × 10-4. This indicates that RITA has higher accuracy and speed than GITM. FIG. 8(a) shows DCS for Compton scattering of bone material in the energy range of 5-150 keV. Fig. 8(b) and (c) are interpolation errors between PDF by RITA interpolation and Generalized Inverse Transformation Method (GITM) and the original Rayleigh DCS sampled from compton DCS with energies of 30keV and 60keV, respectively, and the interpolation grid point number of GITM is 2048. These above show that using RITA sample scatter angles, the simulation speed and accuracy can be optimized.
After sampling the angular distribution of the scatter, the photon generation (i-1) is calculated, i-2, …, n +1 order scatter, at the interaction point ai-1Edge of
Figure GDA0003307176460000152
Probability of direction running out of motif M
Figure GDA0003307176460000153
And is provided with Wi＝Wi-1*(1-pi-1)。
In step S460, the next-order scattering angle position is determined by sampling.
By using a random number u6(i-1)+4Sampling interaction point Ai,i＝2,…,n。
Since the photons are in the motif M, from Ai-1→AiN obeys the distribution:
Figure GDA0003307176460000161
the interaction point A may be sampled by an inverse function sampling methodi. Therefore AiThe specific position of (c) can be given by equation (3):
Figure GDA0003307176460000162
wherein u is6(i-1)+4Is the (6(i-1) +4) th component of the (6n +1) -dimensional Sobol sequence,
Figure GDA0003307176460000163
is a photon at Ai-1After (i-1) -order scattering occurs, along the scattering direction
Figure GDA0003307176460000164
Probability of escaping motifs M.
In implementation, we optimize the algorithm by reducing the number of computations. The method comprises the following specific steps: in equation (3) the calculation is required
Figure GDA0003307176460000165
Is integrated, but
Figure GDA0003307176460000166
Is a path integral
Figure GDA0003307176460000167
A part of (a). Can be found in
Figure GDA0003307176460000168
Reducing paths by pre-calculating integration results at different positions
Figure GDA0003307176460000169
Integral calculation of (a) wherein
Figure GDA00033071764600001610
The different locations being generally paths
Figure GDA00033071764600001611
The bisector point of (a).
Step S470, forcing a random selection of detectors and calculating the probability of photons reaching the detectors, includes: with the 6i and 6i +1 components u of the (6n +1) -dimensional Sobol sequence6iAnd u6i+1Forced random selection of interaction point AiN, i-2, …
Figure GDA00033071764600001612
Then
Figure GDA00033071764600001613
Forcibly selecting k-1 interaction points A in the small neighborhoodiCorresponding detector pixel
Figure GDA00033071764600001614
j is 1, … k-1, and the photon is calculated at the interaction point AiProbability density of scattering angle at which Compton scattering occurs
Figure GDA00033071764600001615
Probability density of scattering angle at which Rayleigh scattering occurs
Figure GDA00033071764600001616
And solid angle
Figure GDA00033071764600001617
The probability density of the scattering angle is the normalized scattering Differential Cross Section (DCS), so
Figure GDA00033071764600001618
Wherein muc(Ai,Ei) And mup(Ai,Ei) Respectively representing photons at the interaction point AiAt an energy of EiThe compton and rayleigh differential scattering cross-sections are also known as compton scattering and rayleigh scattering coefficients.
Figure GDA00033071764600001619
0, …, k-1; calculating photon from AiStarting in a forced direction
Figure GDA00033071764600001620
Probability of escaping a motif M
Figure GDA00033071764600001621
So that photons follow the path
Figure GDA00033071764600001622
i ∈ { 2., n }, j ∈ {0, …, k-1} reaches detector pixel DijThe weight of (c):
Figure GDA0003307176460000171
i＝2,…n,j＝0,…k-1；
steps S480 and S490, if i is less than n, the photon continues to scatter next time, i is set to i +1, and the procedure returns to step S440;
step S4100, if i is n, returns
Figure GDA0003307176460000172
The value of (c). Calculating the detected scatter signal is
Figure GDA0003307176460000173
i＝1,…,n,j＝1,…,k-1, wherein
Figure GDA0003307176460000174
For the first photon along the path
Figure GDA0003307176460000175
Reach the detector pixel
Figure GDA0003307176460000176
The weight of (a) is determined,
Figure GDA0003307176460000177
is a function of the detector response, which, in relation to the path,
Figure GDA0003307176460000178
and N is the number of analog photons.
After step S4100, the method may further include: and displaying the result of simulating photon scattering according to the value of the detected scattering signal.
In the embodiment of the invention, the abstract particle transport path problem is converted into a high-dimensional integration problem, and a new method-gQMC based on Monte Carlo simulation is provided for calculating the high-dimensional integration problem. Is different from the existing MC method in principle and sampling implementation. The method for simulating photon scattering by using the Monte Carlo simulation method provided by the embodiment of the invention can greatly shorten the calculation time of scattered photons while achieving the same precision.
The experimental results of the method for simulating photon scattering by using the quasi-Monte Carlo method provided by the embodiment of the invention are as follows:
FIGS. 9(a) and (a1) show the Shepp-Logan motif and its geometric schematic. FIGS. 9(b), (c) and (d) show transverse, coronal and sagittal sections, respectively, of the Shepp-Logan motif. As shown in fig. 9, each voxel size is 0.5 × 0.5mm when composed of 320 × 400 × 360 voxels3The Shepp-Logan motif of (a) verifies the validity of the method. Wherein, the distance from the X-ray source S to the die body and the distance from the die body to the detector are both 500mm, and the cone angle of the light source S is 9.45 degrees. The detector has a resolution of 512 pixels by 512 pixels with a pixel size of 0.8 pixels by 0.8mm2. As computing hardware, a GeForce RTX 2080 Ti graphics card equipped with 4352 processors and 11.0GB global memory was used.
FIG. 10 shows the scatter signals of a 5-120keV cone beam X-ray through the Shepp-Logan phantom, where (a) and (b) are the total scatter signals (sum of first and second order scatterings) calculated by the gMCDRR and gQMC methods, respectively. FIG. 10 shows a gMCDRR simulation 243Photonic and gQMC analog 228The photon, k, is taken as the scattered signal obtained at 16. These images visually demonstrate good agreement between the gQMC results, the gMCDRR results and the gMMC results, indicating the accuracy of the gQMC code. Wherein gMCDRR simulation 243One photon takes about 549.8 minutes, gQMC simulation 228When k is 16, about 2.4 minutes is required, and the relative error between the two is 0.82%.
Finally, it should be pointed out that: the above examples are only for illustrating the technical solutions of the present invention, and are not limited thereto. Those of ordinary skill in the art will understand that: modifications can be made to the technical solutions described in the foregoing embodiments, or some technical features may be equivalently replaced; such modifications or substitutions do not depart from the spirit and scope of the corresponding technical solutions of the embodiments of the present invention.

Claims (7)

1. A method for simulating photon scattering by using a quasi-Monte Carlo method is applied to a system comprising a light source S, a die body M and a detector D, and is characterized in that preset photons are emitted from the light source S and are positioned at an intersection point A0Initially into the phantom M at the interaction point AiI is 1, …, n is i-order scattered, and the intersection A is selected0Corresponding detector pixel
Figure FDA0003307176450000011
And are provided with
Figure FDA0003307176450000012
Selecting k-1 and A in a preset domain range as the centeriCorresponding detector pixel
Figure FDA0003307176450000013
The photons reach the detector pixel after i-order scattering
Figure FDA0003307176450000014
Is marked as
Figure FDA0003307176450000015
Wherein B isiIndicating the edge of the photon after i-order scattering
Figure FDA0003307176450000016
The intersection point of the direction running out of the motif M and the boundary of M comprises the following steps:
step 1, giving a specific numerical value of a scattering order n to generate a (6n +1) -dimensional Sobol point u1,…,u6n+1；
Step 2, using the 1 st component u of the (6n +1) -dimensional Sobol sequence1Sampling the initial energy E of incident photon from a preset energy spectrum phi by an inverse function sampling method0(ii) a Using the 2 nd and 3 rd components u of the (6n +1) -dimensional Sobol sequence2、u3Solid angle region omega for incident light0The incident direction and the intersection point A of the photons are obtained by uniform sampling0From A) in accordance with the following equation 1) to calculate that the photons are not scattered and0edge of
Figure FDA0003307176450000017
Probability p of direction passing through motif0；
Figure FDA0003307176450000018
Wherein the content of the first and second substances,
Figure FDA0003307176450000019
denotes the initial incident direction of the photon, μt(x，E0) Representing photon energy as E0Linear attenuation coefficient of time space coordinate function, total attenuation coefficient being scattering coefficientAnd the sum of the photoelectric absorption coefficients;
and is provided with W1＝W0*(1-p0)，W1Indicating the initial direction of incidence of the photons
Figure FDA00033071764500000110
No scatter stays at the weight of the phantom M; w0＝fl(u2,u3) Wherein f islIs the photon initial incident direction probability density;
step 3, random number u is used according to the following formula 2)4Sampling first order interaction point A1：
Figure FDA00033071764500000111
Wherein u is4Is the 4 th component of the (6n +1) -dimensional Sobol sequence and is represented by a random number u6,u7Selection of and A1Corresponding detector pixel
Figure FDA00033071764500000112
And is arranged at
Figure FDA00033071764500000113
Within a preset domain of k-1 pixels
Figure FDA00033071764500000114
And A1Correspondingly, the photon path is calculated according to
Figure FDA00033071764500000115
Reach the detector pixel
Figure FDA0003307176450000021
The weight of (c):
Figure FDA0003307176450000022
wherein
Figure FDA0003307176450000023
Step 4, if 1 is less than n, the photons continue to generate second-order scattering; let the photon be at the interaction point Ai-1I is 2, …, the sample probability of rayleigh scattering for n +1 is γ, and the (6(i-1) +5) th component u of the (6n +1) -dimensional Sobol sequence is used6(i-1)+5If u is6(i-1)+5If < gamma, it is judged as being Ai-1Rayleigh scattering occurs at the position; otherwise, it is determined as Ai-1Compton scattering occurs; is provided with
Figure FDA0003307176450000024
δ＝0，1，
Figure FDA00033071764500000217
Wherein p isc(Ai-1)、pr(Ai-1) And pa(Ai-1) Respectively representing the photons at the next interaction point Ai-1The proportion of compton scattering, rayleigh scattering and photoelectric absorption occurring;
step 5, sampling scattering angle direction: the (6(i-1) +2) th component u of the (6n +1) dimensional Sobol sequence is used in the RITA algorithm6(i-1)+2To sample the cosine of the scatter polar angle cos (theta)i-1) And using the (6(i-1) +3) th component u of the (6n +1) dimensional Sobol sequence6(i-1)+3Obtaining an azimuthal scattering angle phii＝2π*u6(i-1)+3Then calculating new scattering direction
Figure FDA0003307176450000025
Then calculates the (i-1) order scattering of the photon and then calculates the (i-1) order scattering of the photon at the interaction point Ai-1Edge of
Figure FDA0003307176450000026
Probability of direction running out of motif M
Figure FDA0003307176450000027
And is provided with Wi＝Wi-1*(1-pi-1)；
Step 6) in the direction according to the following formula 3)
Figure FDA0003307176450000028
Upper random number u6(i-1)+4Sampling the next order interaction point Ai,i＝2,…,n：
Figure FDA0003307176450000029
Wherein u is6(i-1)+4Is the (6(i-1) +4) th component of the (6n +1) dimensional Sobol sequence;
step 7, using the 6i and 6i +1 components u of the (6n +1) dimensional Sobol sequence6iAnd u6i+1Selecting and interacting with a point AiN, i-2, …
Figure FDA00033071764500000210
Then
Figure FDA00033071764500000211
Selecting k-1 interaction points A within the preset field rangeiCorresponding detector pixel
Figure FDA00033071764500000212
Calculating photon at interaction point AiProbability density of scattering angle at which Compton scattering occurs
Figure FDA00033071764500000213
Probability density of scattering angle at which Rayleigh scattering occurs
Figure FDA00033071764500000214
And solid angle
Figure FDA00033071764500000215
The probability density of the scattering angle is the normalized scattering differential cross-section DCS, so
Figure FDA00033071764500000216
Wherein muc(Ai,Ei) And mup(Ai,Ei) Respectively representing photons at the interaction point AiAt an energy of EiCompton and rayleigh differential scattering cross sections;
Figure FDA0003307176450000031
calculating the photon from A byiStarting along the direction
Figure FDA0003307176450000032
Probability of running out of motif M:
Figure FDA0003307176450000033
calculating photon path by
Figure FDA0003307176450000034
Reach the detector pixel
Figure FDA0003307176450000035
The weight of (c):
Figure FDA0003307176450000036
step 8, if i is less than n, continuing to scatter the photons for the next time, and returning to the step 4 if i is equal to i + 1;
step 9, if i is equal to n, returning
Figure FDA0003307176450000037
A value of (d);
step 10, calculating the detected scattered signal
Figure FDA0003307176450000038
Wherein
Figure FDA0003307176450000039
For the first photon along the path
Figure FDA00033071764500000310
Reach the detector pixel
Figure FDA00033071764500000311
The weight of (a) is determined,
Figure FDA00033071764500000312
is a function of the detector response, which, in relation to the path,
Figure FDA00033071764500000313
and N is the number of analog photons.
2. The method of claim 1, further comprising:
and displaying the result of simulating photon scattering according to the detected scattering signal.
3. The method of claim 1, wherein the photon motion path follows a distribution after the photons enter the phantom M
Figure FDA00033071764500000314
4. The method of claim 1, further comprising, prior to step 1: calculating the probability of a path of a photon scattered in a measured phantom M and the probability of the photon reaching a detector pixel from an interaction point after scattering by using a high-dimensional integration mode, wherein the probability comprises the following steps:
calculating the photon from A according to0From the beginning to the edge
Figure FDA00033071764500000315
Directional motion, with a probability of staying at the phantom M of
Figure FDA00033071764500000316
Wherein omega0Representing the solid angle of incident light, s representing a point in the solid angle of incident light, flIs the probability density of the initial direction of incidence of the photon,
Figure FDA0003307176450000041
is the value space of the energy spectrum, and phi (E) is the spectral probability density function.
5. The method of claim 4, wherein calculating the probability of the path of the photon scattering within the phantom M using high-dimensional integration comprises:
calculating the scattering of photons from A after the photons generate i-1, i-E {2i-1Starting in the direction
Figure FDA0003307176450000042
Motion, probability of staying within phantom M:
Figure FDA0003307176450000043
wherein
Figure FDA0003307176450000044
Is ai-1Scattering polar angle theta at apex angleδThe solid angle of (a) is,
Figure FDA0003307176450000045
indicates the photon interaction point Ai-1The energy after compton scattering has occurred,
Figure FDA0003307176450000046
indicates the photon interaction point Ai-1The energy after rayleigh scattering.
6. According to the claimsThe method of claim 5, wherein the i-order scattering of the photons is calculated by high-dimensional integration according to the following formula, and the photons are transmitted from the interaction point AiI e { 1.. multidot.n } runs out of the phantom M to reach the detector pixel
Figure FDA0003307176450000047
Probability of (c):
Figure FDA0003307176450000048
wherein the content of the first and second substances,
Figure FDA0003307176450000049
is the scattering point AiAnd corresponding detector pixel
Figure FDA00033071764500000410
The solid angle of (1) is given by
Figure FDA00033071764500000411
And is
Figure FDA00033071764500000412
Wherein
Figure FDA00033071764500000413
Is a detector pixel
Figure FDA00033071764500000414
The length of (a) of (b),
Figure FDA00033071764500000415
to represent
Figure FDA00033071764500000416
And a detector pixel
Figure FDA00033071764500000417
The angle in the normal direction.
7. The method of claim 6, wherein the photon along path is computed by high-dimensional integration according to
Figure FDA00033071764500000418
Reaches a detector pixel after i-order scattering
Figure FDA00033071764500000419
Probability of (c):
Figure FDA00033071764500000420
CN202011526102.9A 2020-12-22 2020-12-22 Method for simulating photon scattering by using quasi-Monte Carlo method Active CN112763519B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011526102.9A CN112763519B (en) 2020-12-22 2020-12-22 Method for simulating photon scattering by using quasi-Monte Carlo method

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011526102.9A CN112763519B (en) 2020-12-22 2020-12-22 Method for simulating photon scattering by using quasi-Monte Carlo method

Publications (2)

Publication Number Publication Date
CN112763519A CN112763519A (en) 2021-05-07
CN112763519B true CN112763519B (en) 2022-03-01

Family

ID=75694492

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011526102.9A Active CN112763519B (en) 2020-12-22 2020-12-22 Method for simulating photon scattering by using quasi-Monte Carlo method

Country Status (1)

Country Link
CN (1) CN112763519B (en)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113804709B (en) * 2021-09-22 2022-08-12 清华大学 Method for correcting CT scattering signal based on quasi-Monte Carlo and forced detection
CN115840243B (en) * 2022-12-26 2024-04-05 成都理工大学 Numerical calculation method for passive efficiency scale of X/gamma ray detector

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2013034709A1 (en) * 2011-09-07 2013-03-14 Katholieke Universiteit Leuven Non-invasive in-situ radiation dosimetry
CN103295188A (en) * 2012-02-28 2013-09-11 上海联影医疗科技有限公司 Path integral method for X-ray Monte Carlo simulation
CN104700448A (en) * 2015-03-23 2015-06-10 山东大学 Self adaption photon mapping optimization algorithm based on gradient
CN108618796A (en) * 2018-02-09 2018-10-09 南方医科大学 A kind of Monte Carlo scattered photon analogy method of task based access control driving
CN109342367A (en) * 2018-09-30 2019-02-15 华中科技大学 A kind of diffusion optical imaging method and system based on control Monte Carlo method

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7952583B2 (en) * 2000-06-19 2011-05-31 Mental Images Gmbh Quasi-monte carlo light transport simulation by efficient ray tracing
JP5426177B2 (en) * 2009-01-09 2014-02-26 株式会社東芝 Simulation method, simulation apparatus, and simulation program
US8983162B2 (en) * 2011-05-11 2015-03-17 Korea Advanced Institute Of Science And Technology Method and apparatus for estimating monte-carlo simulation gamma-ray scattering in positron emission tomography using graphics processing unit

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2013034709A1 (en) * 2011-09-07 2013-03-14 Katholieke Universiteit Leuven Non-invasive in-situ radiation dosimetry
CN103295188A (en) * 2012-02-28 2013-09-11 上海联影医疗科技有限公司 Path integral method for X-ray Monte Carlo simulation
CN104700448A (en) * 2015-03-23 2015-06-10 山东大学 Self adaption photon mapping optimization algorithm based on gradient
CN108618796A (en) * 2018-02-09 2018-10-09 南方医科大学 A kind of Monte Carlo scattered photon analogy method of task based access control driving
CN109342367A (en) * 2018-09-30 2019-02-15 华中科技大学 A kind of diffusion optical imaging method and system based on control Monte Carlo method

Also Published As

Publication number Publication date
CN112763519A (en) 2021-05-07

Similar Documents

Publication Publication Date Title
US10034646B2 (en) Material decomposition of multi-spectral X-ray projections using neural networks
Jia et al. A GPU tool for efficient, accurate, and realistic simulation of cone beam CT projections
CN103384498B (en) Sniffer
US9330458B2 (en) Methods and systems for estimating scatter
CN112763519B (en) Method for simulating photon scattering by using quasi-Monte Carlo method
JP2008520257A5 (en)
CN104254786B (en) Computerized axial tomography imaging method and system
Vernekohl et al. Feasibility study of Compton cameras for x-ray fluorescence computed tomography with humans
JPWO2007145154A1 (en) Compton camera device
Zhao et al. Iterative beam hardening correction for multi-material objects
Wu et al. Xcist—an open access x-ray/ct simulation toolkit
Banjak X-ray computed tomography reconstruction on non-standard trajectories for robotized inspection
de Jong et al. Rapid SPECT simulation of downscatter in non-uniform media
US7272205B2 (en) Methods, apparatus, and software to facilitate computing the elements of a forward projection matrix
Staelens et al. Fast hybrid SPECT simulation including efficient septal penetration modelling (SP-PSF)
Xie et al. Scatter correction for cone-beam computed tomography using self-adaptive scatter kernel superposition
Bowen et al. Design and performance evaluation of a 20-aperture multipinhole collimator for myocardial perfusion imaging applications
Jiang et al. A feasibility study of enhanced prompt gamma imaging for range verification in proton therapy using deep learning
Bouwens et al. Image-correction techniques in SPECT
Peterson et al. Monte Carlo-based quantitative pinhole SPECT reconstruction using a ray-tracing back-projector
KR101493683B1 (en) Super-resolution Apparatus and Method using LOR reconstruction based cone-beam in PET image
KR20110124685A (en) Method and apparatus for gamma-ray scattering estimation in positron emission tomography using graphic processing unit
Zhuo et al. Scatter correction for cone-beam CT via scatter kernel superposition-inspired convolutional neural network
Fin et al. A practical way to improve contrast‐to‐noise ratio and quantitation for statistical‐based iterative reconstruction in whole‐body PET imaging
Mao et al. A tailored ML-EM algorithm for reconstruction of truncated projection data using few view angles

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
GR01 Patent grant
GR01 Patent grant