CN103295188B - The path integral method of X ray Monte Carlo simulation - Google Patents

Abstract

The invention discloses a kind of path integral method of X ray Monte Carlo simulation, comprise the steps: a kind of path integral method of X ray Monte Carlo simulation, comprise the steps: a) to produce [0,1]³ⁿThe low diversity sequence that on space, length is m, wherein n is scattering number of times, m is number of samples; B) by the isomorphism mapping generation of described low diversity sequence and the scattering point sequence of the geometric match of die body; C) by those scattering point sequences, each mesh of detectors lattice point is formed to m bar photon path; D) determine the probability density that arrives detector network point via those photon paths; E), by described probability density weighted sum, obtain the scattering strength on each mesh of detectors lattice point. The path integral method of X ray Monte Carlo simulation provided by the invention, has eliminated the discontinuous error causing on border, makes the variance of sample less, under given precision, uses sample still less to complete calculating, has greatly reduced amount of calculation.

Description

The path integral method of X ray Monte Carlo simulation

Technical field

The present invention relates to Medical Image Processing, relate in particular to a kind of path integral side of X ray Monte Carlo simulationMethod.

Background technology

The basis of X-ray (X ray) scatter correction technology to the scattering numerical simulation of simple geometry model, numerical simulationEmulation is the requisite means of inspection scatter correction technology quality, and high efficiency numbered analog simulation algorithm not only can addThe R&D work of fast scatter correction technology, itself also can directly develop into a kind of scatter correction technology.

In the X ray scattering analogue of existing Monte Carlo, people adopt a series of pseudo random numbers or accurate random sequence trueDetermine outgoing fan angle, subtended angle and the transmission depth of photon, and follow the tracks of the probability of every a branch of photon arrival detector lattice point and amassedPoint. As document " AnefficientMonteCarlo-basedalgorithmforscattercorrection inkeVcone-beamCT”，G.Poludniowski，P.M.Evans，V.N.HansenandS.Webb，The CRFD method that Phys.Med.Biol.54 (2009) 3847-3864. describes.

In this method, the random shooting angle producing and transmission depth can not be guaranteed the scattering point of each random photonDrop on die body inside, these invalid scattering points cause certain additional computational overhead, and in die body boundary integrand notCause continuously certain error, can not make full use of the low difference character of accurate random sequence. It is most that these errors need to provideThe point of amount is made up, thereby causes amount of calculation to increase.

Therefore, be necessary to provide a kind of path integral method of X ray Monte Carlo simulation, integral domain is equal toThe effective coverage of die body, integrand continuously, is eliminated the discontinuous error causing on border in integral domain, reduces and calculatesAmount.

Summary of the invention

Technical problem to be solved by this invention is to provide a kind of path integral method of X ray Monte Carlo simulation, makesObtain the effective coverage that integral domain is equal to die body, integrand continuously, has been eliminated discontinuous leading on border in integral domainThe error causing, reduces amount of calculation.

The present invention solves the problems of the technologies described above the technical scheme adopting to be to provide a kind of X ray Monte Carlo simulationPath integral method, comprises the steps: a) to produce [0,1]³ⁿThe low diversity sequence that on space, length is m, wherein n is scattering timeNumber, m is number of samples; B) by described low diversity sequence being carried out to the geometric match of isomorphism mapping transformation generation and die bodyScattering point sequence; C) by those scattering point sequences, each mesh of detectors lattice point is formed to m bar photon path; D) determine via instituteState the probability density that photon path arrives mesh of detectors lattice point; E), by described probability density weighted sum, obtain each detectorScattering strength on mesh point.

Further, the computing formula of described scattering strength isWherein dΩ is three-dimensional integration infinitesimal,I article of scattering path, the volume that V is die body, m is number of samples, n is scattering timeNumber, p () is probability density.

Further, described low diversity sequence is Halton sequence, and the quasi random number scope of generation is between 0～1.

Further, the geometry of described die body is spheroid or cylindroid or tetrahedron or parallelepiped.

Further, the integer that described scattering frequency n is 0～∞, the total scatter intensity on each mesh of detectors lattice point isThe scattering strength sum obtaining after all scattering number of times calculate respectively.

Further, the total scatter intensity on described each mesh of detectors lattice point is for stating scattering frequency n value 1,2,3 respectivelyThe scattering strength sum obtaining after calculating.

Further, described scattering frequency n value is that the scattering strength of 2 o'clock is that the definite photon path of two scattering points is generalRate is at R⁶Integration on space.

Further, described photon arrives the probability density of mesh of detectors lattice point through primary scatteringBy following formulaCalculate:

Wherein,The probability density of angle of scattering while being θ, P_ComptonPoint probability that Compton scattering occurs, L₁BeX-ray source is to the distance of scattering point, L₂The distance of scattering point to detector lattice point, D₁That X ray arrives scattering point before at mouldThe distance of passing in body, D₂The distance of X ray through passing in die body after scattering, φ_dThat X ray is on detectorIncident angle, μ is the X ray attenuation coefficient of die body, s_dIt is the effective area of detector.

Further, described photon arrives the probability density of detector network point through n scatteringBy following formulaCalculate:

Wherein, n is scattering number of times, and n is greater than 1 positive integer, θ_iThe angle of scattering of the i time scattering,It is angle of scatteringFor θ_iTime probability density, P_ComptonPoint probability that Compton scattering occurs, L_jThe length of j section broken line on photon pathDegree, D_jThe length of j section broken line process in die body on photon path, φ_dThe incident angle of X ray on detector,μ is the X ray attenuation coefficient of die body, s_dIt is the effective area of detector.

The present invention contrasts prior art following beneficial effect: the road of X ray Monte Carlo simulation provided by the inventionFootpath integration method, utilizes accurate random low diversity sequence to replace pseudo random number to produce a series of scattering paths, and shines upon by isomorphismMake the scattering point on these paths evenly drop on die body inside, there is no extra invalid scattering point, thereby eliminated on borderThe discontinuous error causing, makes the variance of sample less, under given precision, uses sample still less to complete calculating, greatlyReduce amount of calculation.

Brief description of the drawings

Fig. 1 is the flow chart of steps of the path integral method of X ray Monte Carlo simulation of the present invention;

Fig. 2 is the bidimensional primary scattering schematic diagram in the photon scattering path of X ray Monte Carlo simulation of the present invention;

Fig. 3 is the three-dimensional primary scattering schematic diagram in the photon scattering path of X ray Monte Carlo simulation of the present invention;

Fig. 4 is the bidimensional double scattering schematic diagram in the photon scattering path of X ray Monte Carlo simulation of the present invention;

Fig. 5 is Halton sequence and the pseudo-random sequence that X ray Monte Carlo simulation path integral method of the present invention is usedThe diversity ratio of integration in relative accuracy is compared with schematic diagram;

Fig. 6 is contrast effect figure before and after isomorphism mapping in the path integral of X ray Monte Carlo simulation of the present invention;

Fig. 7 is number of samples and precision pair before and after isomorphism mapping in the path integral of X ray Monte Carlo simulation of the present inventionCompare design sketch.

Detailed description of the invention

Below in conjunction with drawings and Examples, the invention will be further described.

Fig. 1 is the flow chart of steps of the path integral method of X ray Monte Carlo simulation of the present invention.

Please refer to Fig. 1, the path integral method of X ray Monte Carlo simulation provided by the invention, comprises the steps:

Step S101, produces [0,1]³ⁿThe low diversity sequence that on space, length is m, wherein n is scattering number of times, m is sampleNumber; Low diversity sequence is preferably Halton sequence, and the quasi random number scope of generation, can be by the following ginseng of input between 0～1Number produces: the geometry of die body, position, CT value; The coordinate position of x-ray source and detector trellis; Scattering frequency n; SampleNumber m.

Step S102, by the low diversity sequence in step S101 by generating after isomorphism mapping transformation and the geometric form of die bodyThe scattering point sequence of shape coupling, and shine upon scattering point sequence is evenly distributed in die body as far as possible by isomorphism.

Step S103, forms m bar photon path by these scattering point sequences to each mesh of detectors lattice point.

Fig. 2 is the bidimensional primary scattering schematic diagram in the photon scattering path of X ray Monte Carlo simulation of the present invention; Fig. 3For the three-dimensional primary scattering schematic diagram in the photon scattering path of X ray Monte Carlo simulation of the present invention.

Step 104, determines the probability density that arrives mesh of detectors lattice point via these photon paths;

Described photon arrives the probability density of mesh of detectors lattice point through primary scatteringCan calculate by following formula (1):

Refer to Fig. 2 and Fig. 3, wherein,The probability density of angle of scattering while being θ, P_ComptonThat Compton occursPoint probability of (Compton) scattering, L₁The distance of x-ray source to scattering point, L₂The distance of scattering point to detector lattice point, D₁That x-ray source arrives the distance of passing in die body before scattering point, D₂That X ray passes in die body through after scatteringDistance, φ_dBe the incident angle of X ray on detector, μ is the X ray attenuation coefficient of die body, and the value of μ is by the CT value of die bodyDetermine s_dBe the effective area of detector, h is the dimension of the little voxel of scattering point, Δ φ_sThe subtended angle of little voxel to x-ray source,Δφ_dThe subtended angle of mesh of detectors lattice point to little voxel,It is the normal direction of detector.

About the source of formula (1), please refer to Fig. 3, photon arrives probability and the Δ φ of detector along free routing_s、Δφ_dBe directly proportional respectively, and decay toDoubly, wherein Δ φ_sWithBe inversely proportional to, Δ φ_dWithBe inversely proportional to, Δ φ_dWithTime again with s_dWith cos φ_dBe directly proportional. There is probability and μ and the P of Compton scattering at this corpusculum element along this path in photon_ComptonBe directly proportional respectively. The probability density that photon, along this path, the scattering that angle of scattering is θ occurs isThese independent eventsThe product of probability or probability density is exactly photon arrives mesh of detectors lattice point after 1 scattering probability density along this pathThereby obtain formula (1).

Fig. 4 is the photon of X ray Monte Carlo simulation of the present invention bidimensional double scattering schematic diagram in scattering path. ThreeDimension Multiple Scattering by that analogy.

Please refer to Fig. 4, under three-dimensional case, in the time of photon path process Multiple Scattering, by following formula (2) calculating pathProbability density

Wherein, n is scattering number of times, and is to be greater than 1 positive integer, θ_iThe angle of scattering of the i time scattering,It is scatteringAngle is θ_iTime probability density, P_ComptonPoint probability that Compton scattering occurs, L_jIt is j section broken line on photon pathLength, D_jThe length of j section broken line process in die body on photon path. φ_dThe incidence angle of X ray on detectorDegree, μ is the X ray attenuation coefficient of die body, s_dIt is the effective area of detector.

Step 105, the probability density weighted sum that step S104 is calculated, obtains on each mesh of detectors lattice pointScattering strength. Scattering strength is calculated by following formula (3):

Wherein, d Ω is three-dimensional integration infinitesimal,I article of scattering path, the volume that V is die body, m is sampleNumber, n is scattering number of times. The equation right side of formula (3) represents to approach integrated value with the discrete random point summation that is uniformly distributed,Monte carlo method.

Scattering frequency n value is that the scattering strength of 2 o'clock is at R by the definite photon path probability of 2 scattering points⁶On spaceIntegration, equally can be by R⁶On space, monte carlo method approaches. , scattering strength is loose by n for Multiple ScatteringThe definite photon path probability of exit point is at R³ⁿIntegration on space provides.

Above-mentioned steps flow process is that scattering frequency n can be 0～∞'s about calculating the scattering strength of photon through n scatteringArbitrary integer, if need to calculate the overall strength of multiple different number of times scatterings, the total scatter intensity on each mesh of detectors lattice pointIt is the scattering strength sum obtaining after all scattering number of times calculate respectively. Generally, need to calculate 1～3 scatteringOverall strength, respectively taking n=1|2|3 as input parameter, the scattering strength summation calculating.

More general, n=0 means 0 scattering, i.e. the value of direct projection. In theory, if faling apart n=0～∞Penetrate all summations of intensity and obtain whole intensity that detector receives, photon arrives certain detector from light sourceProbability.

In step S102, for the die body of simple geometric shape, as spheroid, cylindroid, tetrahedron, parallelepipedDeng all easily distributing completely uniformly by isomorphism Mapping implementation. Concrete methods of realizing is as follows:

The geometric parameter of S102.1 die body represents by following four vectors:

$\left[\begin{array}{ccc}{x}_0& y_0& z_0\\ x_1& y_1& z_1\\ x_2& y_2& z_2\\ x_3& y_3& z_3\end{array}\right]$

A. for spheroid, [x₀y₀z₀] be the coordinate position at ellipsoid center, $\left[\begin{array}{ccc}{x}_1& y_1& z_1\\ x_2& y_2& z_2\\ x_3& y_3& z_3\end{array}\right]$ Be respectively three of ellipsoidThe axial vector of individual axle, these three vectors should orthogonal or full rank.

B. for cylindroid, [x₀y₀z₀] be the coordinate position at the center of cylindroid bottom surface, $\left[\begin{array}{ccc}{x}_1& y_1& z_1\\ x_2& y_2& z_2\end{array}\right]$ Be respectivelyThe axial vector of two axles of cylindroid bottom surface, [x₃y₃z₃] be the axial vector of cylindroid central shaft.

C. for parallelepiped, [x₀y₀z₀] be the origin position of parallelepiped, $\left[\begin{array}{ccc}{x}_1& y_1& z_1\\ x_2& y_2& z_2\\ x_3& y_3& z_3\end{array}\right]$ RespectivelyFor from [x₀y₀z₀] the rib vector of three ribs setting out.

D. for tetrahedron, [x₀y₀z₀] be tetrahedral origin position, $\left[\begin{array}{ccc}{x}_1& y_1& z_1\\ x_2& y_2& z_2\\ x_3& y_3& z_3\end{array}\right]$ Be respectively from [x₀y₀z₀] the rib vector of three ribs setting out.

S102.2 is when the low diversity sequence { [uvw] existing on three dimensions_i}，[uvw]_i∈[0，1]³Time, correspondingHow much die bodys on low diversity sequence { [xyz]_iMapping relations as follows,

A. spheroid:

$\begin{array}{c}\left[\begin{array}{ccc}x& y& z\end{array}\right]\\ =\[\begin{array}{ccc}\sqrt[3]{u}\·2\sqrt{(1-v)v}\·\cos \left(2\mathrm{\π}w\right)& \sqrt[3]{u}\·2\sqrt{(1-v)v}\·\left(2\mathrm{\π}w\right)& \sqrt[3]{u}\end{array}\·(1-2v)\]\\ \·\left[\begin{array}{ccc}{x}_1& y_1& z_1\\ x_2& y_2& z_2\\ x_3& y_3& z_3\end{array}\right]+\left[\begin{array}{ccc}{x}_0& y_0& z_0\end{array}\right]\end{array}.$

B. cylindroid:

$\begin{array}{c}\left[\begin{array}{ccc}x& y& z\end{array}\right]\\ =\left[\begin{array}{ccc}\sqrt{u}\·\cos \left(2\mathrm{\π}v\right)& \sqrt{u}\·\sin \left(2\mathrm{\π}v\right)& w\end{array}\right]\·\left[\begin{array}{ccc}{x}_1& y_1& z_1\\ x_2& y_2& z_2\\ x_3& y_3& z_3\end{array}\right]\\ +\left[\begin{array}{ccc}{x}_0& y_0& z_0\end{array}\right]\end{array}.$

C. parallelepiped:

$\left[\begin{array}{ccc}x& y& z\end{array}\right]=\left[\begin{array}{ccc}u& v& w\end{array}\right]\·\left[\begin{array}{ccc}{x}_1& y_1& z_1\\ x_2& y_2& z_2\\ 3& y_3& z_3\end{array}\right]+\left[\begin{array}{ccc}{x}_0& y_0& z_0\end{array}\right].$

D. tetrahedron:

$\begin{array}{c}\left[\begin{array}{ccc}x& y& z\end{array}\right]\\ =\left[\begin{array}{ccc}u\·\sqrt{1-v}\·\sqrt[3]{1-w}& (1-\sqrt{1-v})\·\sqrt[3]{1-w}& 1-\sqrt[3]{1-w}\end{array}\right]\\ \·\left[\begin{array}{ccc}{x}_1& y_1& z_1\\ x_2& y_2& z_2\\ x_3& y_3& z_3\end{array}\right]+\left[\begin{array}{ccc}{x}_0& y_0& z_0\end{array}\right]\end{array}.$

S102.3, to 2 scatterings, should adopt the low diversity sequence on 6 dimension spaces. When having low the diversity sequence { [u of 6 dimension₁v₁w₁u₂v₂w₂]_i, according to the mapping relations of S102.2 respectively by [u₁v₁w₁] generation [x₁y₁z₁], by [u₂v₂w₂] generation [x₂y₂z₂]. For n scattering, adopt 3n to tie up low diversity sequence, mapping relations are by that analogy. Loose for 2 timesPenetrate, the low diversity sequence of 6 dimension can not be reduced to the 3 dimension low diversity sequences replacement separate with 2, in fact 2 the 3 low differences of dimensionSequence is possibility " separate " not. For n scattering, 3n ties up low diversity sequence can not tie up low diversity sequence with n individual 3 equallyReplace.

When there are multiple objects in S102.4, and each object is while having above-mentioned simple geometric shape, adopts in die bodyThe method of stating is to each object configurations mapping relations separately. If any k object, generate the low difference of k group 3 dimension for 1 scatteringSequence; For 2 scatterings, need to generate k²The low diversity sequence of group 6 dimension; For n scattering, need kⁿGroup 3n ties up low differenceSequence.

Below in conjunction with contrast effect figure, the effect of X ray Monte Carlo simulation path integral method of the present invention is verifiedAnd explanation.

Fig. 5 is that X ray Monte Carlo simulation integration method of the present invention is used Halton sequence and pseudo-random sequence integration to existDiversity ratio in relative accuracy is compared with schematic diagram.

As long as make random scatter point enough drop on uniformly die body inside, just can not produce invalid scattering point, and subtractSmall sample variance. Fig. 5 has compared in two-dimensional elliptic shape die body, adopt Halton sequence (one of low diversity sequence) with pseudo-withThe difference of machine sequence integration in relative accuracy; The solid line that indicates " Random " and " Halton " in Fig. 5 represents respectively pseudorandomThe relative error of sequence and Halton sequence is with the relation of sample number N, and dotted line represents error estimate separately. Can see,Reach 10^-4Precision adopts Halton sequence only to need 1000 left and right samples, and pseudo-random sequence reaches same accuracy far awayBe greater than this value. As can be seen here, utilizing low diversity sequence to carry out isomorphism mapping can make scattering point sequence be uniformly distributed as much as possibleIn die body, for simple geometry die body, as spheroid, cylindroid, tetrahedron, parallelepiped etc. all easily reflect by isomorphismPenetrate realization.

Fig. 6 is contrast effect figure before and after isomorphism mapping in the path integral of X ray Monte Carlo simulation of the present invention.

In the prior art, determine the position of scattering point by the angle of emergence and transmission depth, even if adopt low difference orderRow can not ensure that scattering point evenly drops on die body inside how much, and can cause because border is discontinuous precise decreasing. Fig. 6 isThe distribution of triaxial ellipsoid volume scattering intensity detector position, it is long-pending by the low diversity sequence of three groups of different parameters that left and right is respectivelyDivide the result obtaining. Left figure shines upon through isomorphism the even scattering point obtaining, and right figure does not do isomorphism mapping in prior artDirectly the point beyond die body is removed to (cut-off), the number that all sequences drop on the scattering point in die body is 10,000Individual. As can see from Figure 6, the curved surface that in left figure, 3 groups of different sequences obtain is almost consistent, and the different sequences of right Fig. 3 groupThe curved surface obtaining is some difference, illustrates in the path integral method of X ray Monte Carlo simulation of the present invention and adopt isomorphism to reflectThe mode precision of penetrating is higher.

Fig. 7 is number of samples and precision pair before and after isomorphism mapping in the path integral of X ray Monte Carlo simulation of the present inventionCompare design sketch.

As can be seen from Figure 7, the precision that isomorphism is shone upon the even scattering point obtaining is significantly better than the scattering point of cut-offPrecision.

Therefore, the path integral method of X ray Monte Carlo simulation provided by the invention, by the random low diversity sequence of standardProduce a series of scattering paths, and calculate the probability density of every paths, integration (summation) obtains scattering strength in addition. Pass throughSuitable isomorphism mapping, the scattering point on these paths evenly drops on die body inside, there is no extra invalid scattering point. Simultaneously long-pendingSubregion is equal to the effective coverage of die body completely, and integrand continuously, has been eliminated discontinuous leading on border in integral domainThe error causing. In addition, under same sample size (being the number of scattering point), the precision of Monte Carlo integral algorithm depends onThe variance of sample, utilizes the low difference character of accurate random sequence, can make the variance of Monte Carlo integral algorithm reach minimum, therebyUnder given precision, use sample still less to complete calculating, greatly reduce amount of calculation.

Although the present invention discloses as above with preferred embodiment, so it is not in order to limit the present invention, any this area skillArt personnel, without departing from the spirit and scope of the present invention, when doing a little amendment and perfect, therefore protection model of the present inventionEnclose and work as with being as the criterion that claims were defined.

Claims (9)

1. a path integral method for X ray Monte Carlo simulation, is characterized in that, comprises the steps:

A) produce [0,1]³ⁿThe low diversity sequence that on space, length is m, wherein n is scattering number of times, m is number of samples;

B) by described low diversity sequence being carried out to the scattering point order of the geometric match of isomorphism mapping transformation generation and die bodyRow;

C) by those scattering point sequences, each mesh of detectors lattice point is formed to m bar photon path;

D) determine the probability density that arrives mesh of detectors lattice point via described photon path;

E), by described probability density weighted sum, obtain the scattering strength on each mesh of detectors lattice point.

2. the path integral method of X ray Monte Carlo simulation as claimed in claim 1, is characterized in that, described scattering is strongDegree is calculated as follows:

Wherein, d Ω is three-dimensional integration infinitesimal,I article of scattering path, the volume that V is die body, m is sampleNumber, n is scattering number of times, p () is probability density.

3. the path integral method of X ray Monte Carlo simulation as claimed in claim 1, is characterized in that, described low differenceSequence is Halton sequence, and the quasi random number scope of generation is between 0～1.

4. the path integral method of X ray Monte Carlo simulation as claimed in claim 1, is characterized in that, described die bodyGeometry is spheroid or cylindroid or tetrahedron or parallelepiped.

5. the path integral method of X ray Monte Carlo simulation as claimed in claim 1, is characterized in that, described scattering timeThe integer that number n is 0～∞, the total scatter intensity on each mesh of detectors lattice point is to obtain after all scattering number of times calculate respectivelyScattering strength sum.

6. the path integral method of X ray Monte Carlo simulation as claimed in claim 5, is characterized in that, each detectorTotal scatter intensity on mesh point is the scattering frequency n scattering strength sum that value 1,2,3 obtains after calculating respectively.

7. the path integral method of X ray Monte Carlo simulation as claimed in claim 2, is characterized in that, described scattering timeNumber n values are that the scattering strength of 2 o'clock is that the definite photon path probability of two scattering points is at R⁶Integration on space.

8. the path integral method of the X ray Monte Carlo simulation as described in claim 1 or 7, is characterized in that, described photonArrive the probability density of mesh of detectors lattice point through primary scatteringBe calculated as follows:

Wherein,The probability density of angle of scattering while being θ, P_comptomPoint probability that Compton scattering occurs, L₁It is X raySource is to the distance of scattering point, L₂The distance of scattering point to detector lattice point, D₁That X ray arrives scattering point before in die bodyThe distance of passing, D₂The distance of X ray through passing in die body after scattering,The incident of X ray on detectorAngle, μ is the X ray attenuation coefficient of die body, s_dIt is the effective area of detector.

9. the path integral method of the X ray Monte Carlo simulation as described in claim 1 or 7, is characterized in that, described photonArrive the probability density of detector network point through n scatteringBe calculated as follows:

Wherein, n is scattering number of times, and n is greater than 1 positive integer, θ_iThe angle of scattering of the i time scattering,That angle of scattering is θ_iTime probability density, p_comptomPoint probability that Compton scattering occurs, L_jThe length of j section broken line on photon path,D_jThe length of j section broken line process in die body on photon path,Be the incident angle of X ray on detector, μ isThe X ray attenuation coefficient of die body, s_dIt is the effective area of detector.