CN104027121A - Method for renormalization of X-ray multiple scattering simulation - Google Patents

Abstract

The invention provides a method for renormalization of X-ray multiple scattering simulation. The method comprises a first step of generating a low-difference sequence with the length of m in [0,1]<3n> space, enabling n to be scattering times and be a positive integer larger than 1, and enabling m to be the sample number; a second step of generating a scattering point sequence matched with the geometrical shape of a die body by performing isomorphic mapping transformation on the low-difference sequence; a third step of forming m photon paths for each detector grid point through the scattering point sequence; a fourth step of performing renormalization processing on each photon path among the m photon paths; a fifth step of determining the probability density of photons reaching the detector grid points through n-time scattering along each photon path after renormalization processing; and a sixth step of obtaining scattering intensity of each detector grid point according to the probability density of the photons through n-time scattering along each photon path after renormalization processing. The method can improve multiple scattering simulation efficiency and eliminate errors caused by die body edge incontinuity.

Description

The Renormalization Method of X ray Multiple Scattering simulation

[technical field]

The present invention relates to the technical field of medical image processing, relate in particular to a kind of Renormalization Method of X ray Multiple Scattering simulation.

[technical background]

Scattering numerical simulation based on geometric model is the basis of X ray scatter correction techniques.Generally carry out Monte Carlo simulation by the voxelization object model to obtaining from the CT image undergoing reconstruction at present and study the complex distributions of the scattering radiation in diagnostic roentgenology, and then reach the object of X ray being carried out to scatter correction.

In the X ray scattering analogue of existing Monte Carlo, generally adopt a series of pseudo random numbers or accurate random sequence to determine outgoing fan angle, subtended angle and the projection degree of depth of photon, and follow the tracks of every a branch of photon arrive the probability of detector lattice point and in addition integration obtain the scatter distributions of X ray.As document " An efficient Monte Carlo-based algorithm for scatter correction in keV cone-beam CT ", G.Poludniowski, P.M.Evans, V.N.Hansen and S.Webb, the CRFD method that Phys.Med.Biol.54 (2009) 3847-3864. describes.In the method, the random shooting angle producing and the projection degree of depth can not guarantee that the scattering point of each random photon drops on die body inside, and the scattering point that drops on die body outside is invalid scattering point, and it is large that these invalid scattering points can cause calculating quantitative change; And the discontinuous meeting of integrand at die body boundary causes certain error, and these errors need the point of greater number to be made up, thereby cause amount of calculation to increase.In addition, the method of stating in the use can make the probability of the effective Multiple Scattering photon producing and scattering number of times exponentially decline while carrying out the simulation of Multiple Scattering, cause invalid computation amount excessive, and then can not simulate accurately Multiple Scattering, make on the contrary amount of calculation be index and rise.

Separately there is a kind of path integral method of Monte Carlo simulation to produce a series of scattering paths by the random low diversity sequence of standard, then calculate the probability density of every paths and the probability density in all paths is carried out to integration and can obtain scattering strength.Although the method takes full advantage of the low difference character of accurate random sequence, make the variance of Monte Carlo simulation integral algorithm arrive minimum, also greatly reduced amount of calculation, can the variance of Multiple Scattering dispersed cause integration not restrained because the distance between the accurate random point of Multiple Scattering approaches very much but the method is applied in Multiple Scattering when simulation.

Therefore, the necessary Renormalization Method that a kind of X ray Multiple Scattering simulation is provided, to overcome the defect existing in prior art.

[summary of the invention]

The object of the present invention is to provide a kind of Renormalization Method of X ray Multiple Scattering simulation, can greatly improve the simulation precision of Multiple Scattering, and make integral domain be equal to the effective coverage of die body, integrand is continuous in integral domain, eliminate on die body border the discontinuous error causing, reduce amount of calculation.

For achieving the above object, the present invention is achieved through the following technical solutions: a kind of Renormalization Method of X ray Multiple Scattering simulation, comprises the steps: a) to produce [0,1] ³ⁿthe low diversity sequence that on space, length is m, wherein, n is scattering number of times and for to be greater than 1 positive integer, m is number of samples; B) by described low diversity sequence being carried out to the scattering point sequence of the geometric match of isomorphism mapping transformation generation and die body; C) by scattering point sequence, each mesh of detectors lattice point is formed to m bar photon path; D) each photon path of described m bar photon path is carried out to renormalization processing; E) determine that described photon arrives the probability density of detector trellis point through n scattering along a renormalization photon path after treatment; F) obtain scattering strength each mesh of detectors lattice point on along renormalization each photon path after treatment through the probability density of n scattering according to described photon.

Preferably, described scattering strength I _scatterto process acquisition, i.e. scattering strength I by described photon is weighted to summation along renormalization each photon path after treatment through the probability density of n scattering _scatterobtain by following formula: wherein, i article of scattering path that comprises n scattering point, photon arrives detector trellis point through n scattering probability density along i paths, the volume that V is die body, m is number of samples, n is scattering number of times.

Preferably, described scattering strength I _scatterto process acquisition by described photon is carried out to dual-integration along renormalization all photon paths after treatment through the probability density of n scattering arrival detector trellis point, that is, and scattering strength I _scattterobtain by following formula: I _scatter=∫ ∫ p (L ⁿ) d Ω, wherein, p (L _n) be the probability density of n scattering, d Ω is the integration infinitesimal of 3n dimension space.

Preferably, described scattering strength I _scatterto carry out again the acquisition of angle integration, i.e. scattering strength I by the advanced row Radial Integrals of the scattering point in 3n dimension space _scatterobtain by following formula: wherein, for the angle position of j the scattering point under the spherical coordinate system of initial point taking j-1 scattering point, for the angle integration infinitesimal of j the scattering point under the spherical coordinate system of initial point taking j-1 scattering point, a _jfor along the distance of the path of direction from j-1 scattering point to a border of die body, b _jfor along the distance of the path of direction from j-1 scattering point to another border of die body, r _jfor photon path is processed the distance of previous j-1 scattering point to j scattering point without renormalization, for the Radial Integrals that j scattering point carried out, for the angle integration that j scattering point carried out.

Preferably, described photon path is without the Radial Integrals of j scattering point of renormalization processing equal the Radial Integrals of described photon path through renormalization j scattering point after treatment ? wherein, for photon path renormalization is processed the distance of rear j-1 scattering point to j scattering point, for renormalization operator, r _jfor photon path is processed the distance of previous j-1 scattering point to j scattering point, w without renormalization _jfor renormalization weight.

Preferably, described in except the integral function of the scattering point to j scattering point in 3 (n-1) n-dimensional subspace n, wherein, be photon path while processing without renormalization angle of scattering be θ _jprobability density, P _comptonfor there is point probability of Compton scattering, the X ray attenuation quotient that μ is die body, D _jthe distance that while processing without renormalization for photon path, j section is passed in die body, Ω _{i ≠ j}3 (n-1) n-dimensional subspace n of opening for the scattering point except j scattering point, be photon path while processing without renormalization angle of scattering be θ _iprobability density, D _ithe distance that while processing without renormalization for photon path, i section is passed in die body, s _dfor the area of detector, φ _dthe incident angle of X ray on detector while processing without renormalization for photon path, L _ithe length of i section on photon path while processing without renormalization for photon path.

Preferably, described photon path renormalization is processed the distance of rear j-1 scattering point to j scattering point obtain by following formula: wherein, a _jfor along the distance of the path of direction from j-1 scattering point to a border of die body, b _jfor along the distance of the path of direction from j-1 scattering point to another border of die body, r _jfor photon path is processed the distance of previous j-1 scattering point to j scattering point without renormalization.

Preferably, described renormalization weight w _jobtain by following formula: wherein, a _jfor along the distance of the path of direction from j-1 scattering point to a border of die body, b _jfor along the distance of the path of direction from j-1 scattering point to another border of die body.

Preferably, described photon along a renormalization photon path after treatment through the probability density p of n scattering (L ⁿ) obtain by following formula: $p\left(L^n\right)=\left({\mathrm{\&Pi;}}_j^n{P}_{v_0}\left({\mathrm{\&theta;}}_j^{\&prime;}\right)P_{\mathrm{Compton}}\mathrm{\&mu;}\right)\&CenterDot;e^{-{\mathrm{\&Sigma;}}_j^{n+1}\mathrm{\&mu;}{D}_j^{\&prime;}}{s}_d\cos {\mathrm{\&phi;}}_d^{\&prime;}{\mathrm{\&Pi;}}_j^{n+1}{w}_j/\left(4\mathrm{\&pi;}\right),$ Wherein, for photon path angle of scattering after renormalization is processed is θ ' _jprobability density, P _comptonfor there is point probability of Compton scattering, the X ray attenuation quotient that μ is die body, D ' _jfor the photon path distance that on photon path, j section is passed in die body after renormalization is processed, s _dfor the effective area of detector, φ ' _dfor photon path X ray angle of incidence on detector after renormalization is processed, w _jfor renormalization weight, n is scattering number of times.

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

The accurate random low diversity sequence of Renormalization Method utilization of X ray Multiple Scattering simulation provided by the invention replaces pseudo random number to produce a series of scattering series and carries out renormalization processing by the photon path that scattering series is formed, can make photon can amass continuously in integral domain in the probability density that arrives detector trellis point through n scattering through renormalization all photon paths after treatment, eliminate on die body border the discontinuous error causing, and make photon arrive the probability density of detector trellis point through n scattering along i paths variance convergence, thereby reduce amount of calculation, greatly improve the simulation precision of Multiple Scattering.

[brief description of the drawings]

Fig. 1 is the schematic flow sheet of the Renormalization Method of X ray Multiple Scattering simulation of the present invention.

Fig. 2 is the schematic diagram of the photon path of the X ray Multiple Scattering of the present invention simulation scatter distributions situation under spherical coordinate system.

Fig. 3 be a photon path in X ray Multiple Scattering of the present invention simulation in die body without renormalization process with through renormalization comparison diagram after treatment, what wherein dotted line represented is the photon path without renormalization processing in die body, and what solid line represented is through renormalization photon path after treatment in die body.

Fig. 4 is that the photon path that is made up of Halton sequence adopts method of the present invention to carry out renormalization to process with the photon path being made up of Halton sequence and do not carry out the diversity ratio of renormalization processing in relative accuracy compared with schematic diagram.

[detailed description of the invention]

Below in conjunction with the drawings and specific embodiments, the Renormalization Method of X ray Multiple Scattering simulation of the present invention is described in further detail.According to the following describes and claims, advantages and features of the invention will be clearer.It should be noted that, accompanying drawing all adopts very the form of simplifying and all uses non-ratio accurately, only for convenient, the object of the aid illustration embodiment of the present invention lucidly.

The accurate random low diversity sequence of Renormalization Method utilization that the invention provides a kind of X ray Multiple Scattering simulation replaces pseudo random number to produce a series of scattering series and carries out renormalization processing by the photon path that scattering series is formed, can make photon can amass continuously in integral domain in the probability density that arrives detector trellis point through n scattering through renormalization all photon paths after treatment, eliminate on die body border the discontinuous error causing, reduce amount of calculation, greatly improve the simulation precision of Multiple Scattering.

Fig. 1 is the schematic flow sheet of the Renormalization Method of X ray Multiple Scattering simulation of the present invention.The Renormalization Method of this X ray Multiple Scattering simulation comprises the following steps:

S11, generation [0,1] ³ⁿthe low diversity sequence that on space, length is m, wherein, n is scattering number of times and is to be greater than 1 positive integer, m is number of samples, described low diversity sequence is preferably Halton sequence, the quasi random number scope generating is between 0～1, and described low diversity sequence can produce by the coordinate position of the coordinate position of the CT value of the position of the geometry of input parameter die body, die body, die body, x-ray source, mesh of detectors lattice point, scattering frequency n, number of samples m;

S12, described low diversity sequence is carried out to the scattering point sequence of the geometric match of isomorphism mapping transformation generation and die body, and scattering point sequence is evenly distributed in die body as far as possible;

S13, by described scattering point sequence, each mesh of detectors lattice point is formed to m bar photon path;

S14, each photon path of described m bar photon path is carried out to renormalization processing;

S15, definite probability density that is arrived detector trellis point by described photon along renormalization every photon path after treatment through n scattering;

S16, obtain scattering strength each mesh of detectors lattice point on along renormalization each photon path after treatment through the probability density of n scattering according to described photon.

Fig. 2 shows a distribution situation in die body in the scattering point sequence described in step S12, step S13, taking the scattering point of the j-1 time scattering as initial point, and the j time scattering point distribution situation on assigned direction.What Fig. 2 a showed the is scattering point of the j-1 time scattering is in the situation of die body inside, Fig. 2 b demonstration be that the scattering point of the j-1 time scattering is in the situation of die body outside.

Particularly, the scattering strength I in S16 _scatterto process acquisition, i.e. scattering strength I by described photon is weighted to summation along renormalization each photon path after treatment through the probability density of n scattering _scatterobtain by following formula: wherein, i article of scattering path that comprises n scattering point, photon arrives detector trellis point through n scattering probability density along i paths, the volume that V is die body, m is number of samples, n is scattering number of times.

Described scattering strength I _scatteralso can carry out dual-integration and process acquisition by described photon is arrived to the probability density of detector trellis point along renormalization all photon paths after treatment through n scattering, that is, and scattering strength I _scatteralso can obtain by following formula: I _scatter=∫ ∫ p (L _n) d Ω, wherein, p (L ⁿ) be the probability density of n scattering, d Ω is the integration infinitesimal of 3n dimension space.

Particularly, it is exactly that the advanced row Radial Integrals of the scattering point in 3n dimension space is carried out to angle integration again that described photon is carried out to dual-integration along renormalization all photon paths after treatment through the probability density of n scattering arrival detector trellis point, wherein, for the angle position of j the scattering point under the spherical coordinate system of initial point taking j-1 scattering point, for the angle integration infinitesimal of j the scattering point under the spherical coordinate system of initial point taking j-1 scattering point, a _jfor along the distance of the path of direction from j-1 scattering point to a border of die body, b _jfor along the distance of the path of direction from j-1 scattering point to another border of die body, r _jfor photon path is processed the distance of previous j-1 scattering point to j scattering point without renormalization, for the Radial Integrals that j scattering point carried out, for the angle integration that j scattering point carried out.Wherein, described photon path is without the Radial Integrals of j scattering point of renormalization processing equal the Radial Integrals of described photon path through renormalization j scattering point after treatment ? wherein, for photon path renormalization is processed the distance of rear j-1 scattering point to j scattering point, for renormalization operator, r _jfor photon path is processed the distance of previous j-1 scattering point to j scattering point, w without renormalization _jfor renormalization weight, shown in Fig. 3.Fig. 3 illustrates the path renormalization schematic diagram of photon die body scattering 3 times.Described except the integral function of the scattering point to j scattering point in 3 (n-1) n-dimensional subspace n, wherein, be photon path while processing without renormalization angle of scattering be θ _jprobability density, P _comptonfor there is point probability of Compton scattering, the X ray attenuation quotient that μ is die body, D _jthe distance that while processing without renormalization for photon path, j section is passed in die body, Ω _{i ≠ j}3 (n-1) n-dimensional subspace n of opening for the scattering point except j scattering point, be photon path while processing without renormalization angle of scattering be θ _iprobability density, D _ithe distance that while processing without renormalization for photon path, i section is passed in die body, s _dfor the area of detector, φ _dthe incident angle of X ray on detector while processing without renormalization for photon path, L _ithe length of i section on photon path while processing without renormalization for photon path.

Described photon path renormalization is processed the distance of rear j-1 scattering point to j scattering point obtain by following formula: wherein, a _jfor along the distance of the path of direction from j-1 scattering point to a border of die body, b _jfor along the distance of the path of direction from j-1 scattering point to another border of die body, r _jfor photon path is processed the distance of previous j-1 scattering point to j scattering point without renormalization.

Described renormalization weight w _jobtain by following formula: wherein, a _jfor along the distance of the path of direction from j-1 scattering point to a border of die body, b _jfor along the distance of the path of direction from j-1 scattering point to another border of die body.

Particularly, described photon along a renormalization photon path after treatment through the probability density p of n scattering (L ⁿ) also can obtain by following formula: $p\left(L^n\right)=\left({\mathrm{\&Pi;}}_j^n{P}_{v_0}\left({\mathrm{\&theta;}}_j^{\&prime;}\right)P_{\mathrm{Compton}}\mathrm{\&mu;}\right)\&CenterDot;e^{-{\mathrm{\&Sigma;}}_j^{n+1}\mathrm{\&mu;}{D}_j^{\&prime;}}{s}_d\cos {\mathrm{\&phi;}}_d^{\&prime;}{\mathrm{\&Pi;}}_j^{n+1}{w}_j/\left(4\mathrm{\&pi;}\right),$ Wherein, for photon path angle of scattering after renormalization is processed is θ ' _jprobability density, P _comptonfor there is point probability of Compton scattering, the X ray attenuation quotient that μ is die body, D ' _jfor the photon path distance that on photon path, j section is passed in die body after renormalization is processed, s _dfor the effective area of detector, φ ' _dfor photon path X ray angle of incidence on detector after renormalization is processed, w _jfor renormalization weight, n is scattering number of times.

Fig. 4 shows the renormalization of the photon path employing the inventive method being made up of Halton sequence and processes rear and do not adopt the renormalization of the inventive method to process the difference comparison diagram in relative accuracy.Fig. 4 a is the schematic diagram that is related to of the photon path that is made up of the Halton sequence error of not carrying out the number of samples that adopts when renormalization is processed and generation.The error of Fig. 4 b number of samples that method renormalization of the present invention adopts after processing for the photon path being made up of Halton sequence adopts and generation be related to schematic diagram.Comparison diagram 4a and Fig. 4 b can clearly draw, be greater than 10 at number of samples m ⁴time, will be much smaller than the error without renormalization processing through renormalization error after treatment in the same number of number of samples situation of employing, the number of samples adopting in the situation that reaching corresponding precision after renormalization is processed at photon path can tail off greatly, also just can greatly reduce amount of calculation, improve the efficiency of Multiple Scattering simulation.

To sum up, these are only preferred embodiment of the present invention, should not limit the scope of the invention with this, i.e. every simple equivalence of doing according to claims of the present invention and description of the present invention changes and modifies, and all should still remain within the scope of the patent.

Claims (10)

1. a Renormalization Method for X ray Multiple Scattering 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 and for to be greater than 1 positive integer, m is number of samples;

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

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

D) each photon path of described m bar photon path is carried out to renormalization processing;

E) determine that described photon arrives the probability density of detector trellis point through n scattering along a renormalization photon path after treatment;

F) obtain scattering strength each mesh of detectors lattice point on along renormalization each photon path after treatment through the probability density of n scattering according to described photon.

2. the Renormalization Method of X ray Multiple Scattering simulation as claimed in claim 1, is characterized in that described scattering strength I _scatterto process acquisition, i.e. scattering strength I by described photon is weighted to summation along renormalization each photon path after treatment through the probability density of n scattering _scatterobtain by following formula:

$I_{\mathrm{scatter}}=\underset{m\&RightArrow;\&infin;}{\lim }\underset{i}{\overset{m}{\mathrm{\&Sigma;}}}p\left(L_i^n\right)\frac{V^n}{m}$

Wherein, i article of scattering path that comprises n scattering point, photon arrives detector trellis point through n scattering probability density along i paths, the volume that V is die body, m is number of samples, n is scattering number of times.

3. the Renormalization Method of X ray Multiple Scattering simulation as claimed in claim 1, is characterized in that described scattering strength I _scatterto process acquisition by described photon is carried out to dual-integration along renormalization all photon paths after treatment through the probability density of n scattering arrival detector trellis point, that is, and scattering strength I _scatterobtain by following formula: I _scatter=∫ ∫ p (L ⁿ) d Ω, wherein, p (L ⁿ) be the probability density of n scattering, d Ω is the integration infinitesimal of 3n dimension space.

4. the Renormalization Method of X ray Multiple Scattering simulation as claimed in claim 3, is characterized in that described scattering strength I _scatterto carry out again the acquisition of angle integration, i.e. scattering strength I by the advanced row Radial Integrals of the scattering point in 3n dimension space _scatterobtain by following formula: wherein, for the angle position of j the scattering point under the spherical coordinate system of initial point taking j-1 scattering point, for the angle integration infinitesimal of j the scattering point under the spherical coordinate system of initial point taking j-1 scattering point, a _jfor along the distance of the path of direction from j-1 scattering point to a border of die body, b _jfor along the distance of the path of direction from j-1 scattering point to another border of die body, r _jfor photon path is processed the distance of previous j-1 scattering point to j scattering point without renormalization, for the Radial Integrals that j scattering point carried out, for the angle integration that j scattering point carried out.

5. the Renormalization Method of X ray Multiple Scattering simulation as claimed in claim 4, is characterized in that, described photon path is without the Radial Integrals of j scattering point of renormalization processing equal the Radial Integrals of described photon path through renormalization j scattering point after treatment ? wherein, for photon path renormalization is processed the distance of rear j-1 scattering point to j scattering point, for renormalization operator, r _jfor photon path is processed the distance of previous j-1 scattering point to j scattering point, w without renormalization _jfor renormalization weight.

6. the Renormalization Method of X ray Multiple Scattering simulation as claimed in claim 5, is characterized in that, described in except the integral function of the scattering point to j scattering point in 3 (n-1) n-dimensional subspace n, wherein, be photon path while processing without renormalization angle of scattering be θ _jprobability density, P _comptonfor there is point probability of Compton scattering, the X ray attenuation quotient that μ is die body, D _jthe distance that while processing without renormalization for photon path, j section is passed in die body, Ω _{i ≠ j}3 (n-1) n-dimensional subspace n of opening for the scattering point except j scattering point, be photon path while processing without renormalization angle of scattering be θ _iprobability density, D _ithe distance that while processing without renormalization for photon path, i section is passed in die body, s _dfor the area of detector, φ _dthe incident angle of X ray on detector while processing without renormalization for photon path, L _ithe length of i section on photon path while processing without renormalization for photon path.

7. the Renormalization Method of X ray Multiple Scattering simulation as claimed in claim 5, is characterized in that, described photon path renormalization is processed the distance of rear j-1 scattering point to j scattering point obtain by following formula: wherein, a _jfor along the distance of the path of direction from j-1 scattering point to a border of die body, b _jfor along the distance of the path of direction from j-1 scattering point to another border of die body, r _jfor photon path is processed the distance of previous j-1 scattering point to j scattering point without renormalization.

8. the Renormalization Method of X ray Multiple Scattering simulation as claimed in claim 7, is characterized in that described renormalization weight w _jobtain by following formula: wherein, a _jfor along the distance of the path of direction from j-1 scattering point to a border of die body, b _jfor along the distance of the path of direction from j-1 scattering point to another border of die body.

9. the Renormalization Method of X ray Multiple Scattering as claimed in claim 8 simulation, is characterized in that, described photon along a renormalization photon path after treatment through the probability density p of n scattering (L ⁿ) obtain by following formula $:p\left(L^n\right)=\left({\mathrm{\&Pi;}}_j^n{P}_{v_0}\left({\mathrm{\&theta;}}_j^{\&prime;}\right)P_{\mathrm{Compton}}\mathrm{\&mu;}\right)\&CenterDot;e^{-{\mathrm{\&Sigma;}}_j^{n+1}\mathrm{\&mu;}{D}_j^{\&prime;}}{s}_d\cos {\mathrm{\&phi;}}_d^{\&prime;}{\mathrm{\&Pi;}}_j^{n+1}{w}_j/\left(4\mathrm{\&pi;}\right),$ Wherein for photon path angle of scattering after renormalization is processed is θ ' _jprobability density, P _comptonfor there is point probability of Compton scattering, the X ray attenuation quotient that μ is die body, D ' _jfor the photon path distance that on photon path, j section is passed in die body after renormalization is processed, the effective area that sd is detector, φ ' _dfor photon path X ray angle of incidence on detector after renormalization is processed, w _jfor renormalization weight, n is scattering number of times.

10. the Renormalization Method that X ray Multiple Scattering is simulated as in one of claimed in any of claims 1 to 9, is characterized in that, described low diversity sequence is Halton sequence, and the quasi random number scope of generation is between 0～1.