Scattered photons cause blurring and distortions in flash radiography, reducing the accuracy of image reconstruction significantly. The effect of the scattered photons is taken into account and an iterative deduction of the scattered photons is proposed to amend the scattering effect for image restoration. In order to deduct the scattering contribution, the flux of scattered photons is estimated as the sum of two components. The single scattered component is calculated accurately together with the uncollided flux along the characteristic ray, while the multiple scattered component is evaluated using correction coefficients preobtained from Monte Carlo simulations.The arbitrary geometry pretreatment and ray tracing are carried out based on the customization of AutoCAD. With the above model, an
I
terative
P
rocedure for image rest
OR
ation code, IPOR, is developed. Numerical results demonstrate that the IPOR code is much more accurate than the direct reconstruction solution without scattering correction and it has a very high computational efficiency.
1. INTRODUCTION
In flash radiography, the linear attenuation coefficients or the interface locations of the investigated object can be reconstructed using the information from the radiography image. However, the scattered photons in the image signal will blur the radiography or even overwhelm the useful information if there are a lot of strong scattering materials in the object. Hence, it is of great significance to estimate the scattered photons to deduct the scattering contribution.
There are many numerical methods associated with the scattering correction including the point scattered functions
[1
,
2]
, MonteCarlo scattering simulations
[3]
, deterministic approaches for transport equation
[4
,
5]
, and others
[6

11]
. The point scattered functions, employed by N. Kardjilov
[1]
and R. Hassanein
[2]
can improve the reconstruction accuracy for known simple geometries, but require a parameters library for point scattered functions based on different samples, different sampledetector distances, and other influencing factors. Thus, a large number of numerical calculations have to be done to obtain good results making it very timeconsuming especially when the geometry and/or the material of the sample are changed. Another method is based on the MonteCarlo simulations and has a high accuracy for arbitrary geometry. This method solves the transport processes by tracing the Xray track, But the computing time is tremendously long due to the timeexpensive stochastic Monte Carlo transport calculation method. Alternatively, the deterministic transport calculation method can also be employed to calculate the photon scattering flux by solving the Boltzmann photon transport equation. It is very difficult to find a method with both the capability of processing complex geometry and high computational efficiency.
In this paper, both the stochastic method and the deterministic method are adopted. The uncollided ray flux and the single scattered ray flux can be precisely calculated along the characteristic lines by use of the Method Of Characteristic (MOC) while the multiple scattered ray flux is corrected using numerical correction coefficients obtained through preliminary Monte Carlo simulations. Furthermore, the geometry pretreatment and ray tracing are expanded to 3dimensional arbitrary geometry. On one hand, this coupled method overcomes the difficulties of geometric expansion experienced when using the deterministic method. On the other hand it reduces the calculation time of the multiple scattered photons and thus has a very high computational efficiency.
The remainder of this paper is organized as follows: Section 2 shows the methodology of this paper including the theoretical models of image reconstruction and scattering flux calculation, and the outline of the scattering correction during the iterative reconstruction. Section 3 analyzes the numerical results and demonstrates the reliability and capability of the method. Finally some conclusions are drawn and a few suggestions are made in Section 4.
2. METHODOLOGY
In radiographic inspection, photons arriving at the detectors can be divided into the collided and uncollided components as shown in Eq.(1).
where
represents the uncollided photons,
represents the collided photons which can also be called scattered photons, and
is the spatial coordinates vector (
x
,
y
,
z
). In image reconstruction, the uncollided photons are used to evaluate the mass thickness of the sample which satisfies the exponential law as shown in Eq.(2).
where
is the photon source,
E
is photon energy,
is the direction vector,
μ
and
ω
are angular variables,
is the macroscopic cross section which is dependent on space and energy, and
R
denotes the photons path length along the direction vector. Treating the right side of Eq.(2) by numerical integration, it can be written as Eq.(3), and after rearranging, Eq.(3) may be rewritten as Eq. (4).
The density of the object can be obtained by calculating.
where Σ
_{k}
is the average linear attenuation coefficients of the kdistrict material of the investigated object,
d_{k}
is the corresponding path length, and
r
is the distance from the photon source to the detector.
By solving equations (4) and (5), the spatial density distribution of the inspected object can be reconstructed. But in fact the reconstructed accuracy would be reduced or even distorted because the term
on the right side of Eq.(4) is disturbed by
. Especially when there are a lot of strong scatting materials in the investigated object, the scattering noise would be much bigger than the useful information. Hence, it would be very important to estimate the scattered photons to correct the scattering disturbance.
 2.1 Theoretical Model of the Scattered Flux Calculation
The scattered photons can be described by the integral transport equation as shown in Eq.(6).
where
is the angular photon source,
R
is the photon path length from the photon source to the detector, and
E
^{'}
is the photon energy after collision.
The scattered photons can be divided into single scattered ray flux and multiple scattered ray flux as shown in Eq.(7). In this paper, the first scattered ray flux of the detectors are precisely calculated along the characteristic lines and the multiple scattered ray flux is corrected by numerical coefficients
[10]
.
where
is the first scattered ray flux, and
is the multiple scattered ray flux.
 2.1.1. Theoretical Model of the First Scattered Ray Flux
The first term on the righthand side of Eq.(7) is the contribution of the first scattered photons, which can be calculated using a numerical expression as in Eq.(8).
where
is the first collided photon angular scattering source.
Eq.(8) can be solved using a twostep algorithm. The first step is the spatial discretization and the computation of the scattering source for each space interval. The second step is the calculation of the scattered photons from the scattering source to the radiography image (detectors).
 1) Spatial discretization
The spatial discretization is provided by ANSYS packages which can discrete the object into arbitrary geometry meshes as shown in
Fig. 1
. The scattering source is assumed to be continuous in each mesh and then Eq.(8) can be rewritten as
where
and
NM
is the total number of the spatial meshes,
p
(
θ
) represents the probability of the scattering source pointing to the Δ
θ
angular region of angular
θ
. After numerical integration, Eq.(8) can be rewritten as the following
Spatial Discretization
where Σ
_{t}
(
E
^{'}
,
j
) is the total macroscopic crosssection dependent on the jdistrict material of the sample and the scattered ray energy
E
,
d_{j}
is the photon path length of the jdistrict material.
 2) Angular scheme
To solve Eq.(11), the probability of the scattering source pointing to the Δ
θ
angular region of angular
θ
must be obtained. There are several different types of scattering for those photons going through the matter, and different scattering types need different angular schemes. Photoelectric effect, Compton scattering and pair production are the main scattering types. Among them the particle energy after the photoelectric effect is not very high and thus can be assumed to deposit in situ
[8]
. In this paper,
P_{c}
(
θ
) is used to represent the angular probability of Compton scattering, and
P_{spp}
(
θ
) for pair production.
 (1) Angular scheme for Compton scattering

Compton scattering is a kind of strong anisotropic scattering, especially for highenergy initial photons. Most of the photons, after Compton scattering, fly forward.Fig. 2is a schematic diagram of the Compton scattering of the ray from the photon source to the detector. The probability of the scattering source pointing to the Δθangular region of angularθcan be expressed by Eq.(12).
where
θ
is the scattering angle,
denotes the differential cross section of Compton scattering,
d
Ω
_{d}
is the solid angle of the detector, and
σ_{c}
represents Compton scattering crosssection for an target electron. The differential crosssection of Compton scattering can be calculated by the wellknown KleinNishina expression, Eq.(13).
Schematic Diagram for Compton Scattering
where
r_{e}
is the classical electron radius,
E_{γ}
is the initial photon energy, and
E_{γ}
^{'}
is the photon energy after scattering. Suppose the coordinates of the ray source are
s
(
x
_{0}
,
y
_{0}
,
z
_{0}
), the discrete mesh coordinates are
p
(
x_{p}
,
y_{p}
,
z_{p}
), and the coordinates of the detector are
d
(
x_{t}
,
y_{t}
,
z_{t}
), it can be found that
θ
is the angle between the vectors
which can be determined by Eq.(14).
Integrating
over the 4
π
angle domain leads to
σ_{c}
.
where
 (2) Angular scheme for pair production

Considering pair production, when a photon interacts with a nucleus, an electron and a positron may be created. These particles will lose their energy during their travels in the matter. Eventually, the positron will annihilate along with other electrons by producing two photons. The energy of these two photons are both 0.511MeV and their flight vectors are exactly opposite. Based on this characteristics, pair production obeys the isotropic distribution, so the angle ratio after collision of pair production can be written as Eq.(16).
where
r
^{'}
is the distance between the scattering source and the detector.
 3) Source scheme
Eq. (17) provides the expression for the first scattered photon source
in Eq.(11).
where the first scattered photon flux for each discretized mesh obeys the exponential attenuation rule from the ray source to the mesh as shown in Eq.(18).
where
is the space,
is the dependent photon flux,
n
indicates the number of photons released by each reaction, Σ
_{ti}
(
E
) is the total macroscopic crosssection for the kdistrict material,
d_{i}
is the corresponding path length,
dV
represents the volume of each discrete mesh, and
r
is the distance from the ray source to the detector.
Because only one photon is released by the Compton effect, but two photons are released by pair production, the scattering source for Compton scattering and pair production can respectively computed by Eq.(19) and Eq.(20).
where
is the scattering source of Compton scattering,
is the scattering source of pair production, Σ
_{sc}
(
E
) is the macroscopic crosssection of Compton scattering, and Σ
_{spp}
(
E
) is macroscopic crosssection of the pair production.
Integrating the flux contribution from all the discrete meshes to the detector, we can obtain the first scattered photon flux.
where
is the first scattered photon flux produced by Compton scattering,
is the first scattered photon flux produced by pair production. Adding these together yields the first scattered photon flux of the detector, as represented in Eq.(23).
 4) Ray tracing
Based on the above theoretical model, we need to find the photon path length of each material region along the ray transport track. After geometry description and spatial discretization, ray tracing can be done. The intersections of the ray and the solid model surfaces can be searched and sorted to calculate the segment lengths in each region as shown in
Fig. 3
. Thus, information, including the material region number passed through, and the path length in every region, etc., can be obtained.
Intersections Calculation and Sorting during in Ray Tracing
For an object in a known simple geometry such as a slab or spherical shell, etc., the ray tracing can be achieved accurately and efficiently, but for arbitrary geometry the geometry description, especially the ray tracing, becomes much more complicated. In the method of this paper, the ray tracing for arbitrary geometry is established using the customization of AutoCAD
[12
,
13]
to take advantage of its powerful geometric description and editing capability. This method helps to breakdown the geometry limitation into 3dimensional arbitrary geometry, which greatly improves the capability of geometry treatment.
 2.1.2 Theoretical Model of the Multiple Scattered Ray Flux
The second term on the righthand side of Eq.(7) is the contribution of multiple scattered photons. Although its percentage compared to the total scattered flux is always small, the corresponding computation time is much longer than that of the first scattered flux. As Eq. (7) shows, there is an option to cut the multiple scattered photons off. For a sample with a long optical path, such as longer than 2 mean free paths (mfp), the ratio of the multiple scattered photons to the first scattered ones may reach up to 20% which is cannot be ignored. However, the sample contains lots of heavy media making its optical length longer than 2 mfps. Thus there is no choice but to take account of the multiple scattered flux for the sample with plenty of heavy material.
In this paper, the characteristics of the multiple scattered photons are studied in detail in terms of different parameters including geometry, material and its mass density by Monte Carlo simulations.
Figure 4
shows a sample that has a spherical shell structure and contains lots of strong scattering heavy materials. By keeping the total mass constant, the geometry and the density were compressed as shown in
Table 1
. It can be found that the percentage of the multiple scattered photons only varies slightly with heavy changes
Schematic Diagram for a Spherical Shell Geometry Object
Different Density Cases of the Object
Different Density Cases of the Object
Proportion of the Multiple Scattered Flux to the Total Scattered Flux with Different Density
of material mass density as shown in
Fig.5
. Thus it is clear that the percentage of the multiple scattered photons approximately remains the same as long as the geometry of the object does not change dramatically.
Consequently, a numerical coefficient of multiple scattered flux can be precalculated for the inspected sample using the Monte Carlo simulation. Then the coefficient can be used to make the correction to take account of the multiple scattered flux during the transport computation. Obviously, this can save a large amount of computational time.
 2.1.3 Interpolation Method for Scattered Flux
Recently, the flatpanel detector imaging system has been employed more and more widely in 3dimensional conebeam flash radiography systems. For highresolution imaging systems, there are thousands or even hundreds of thousands of pixels on each single image. Hence, it takes up a huge amount of computing time and memory if the flux of every pixel is simulated in detail. Taking the method of Monte Carlo simulations as an example, largescale workstations or even supercomputers have to be used as alternatives to a microcomputer or a PC machine to simulate the photon transporting to the flatpanel detector imaging.
In this paper, the methods of least square approximation and interpolation are used for the scattering flux calculation. Firstly, typical detectors along the radius vector in the image are selected and their detected flux is calculated. Secondly, by using the method of least square approximation, the scattered flux distribution over the image can be inputted into a function to calculate the distance between the pixel to the center of the image. Thirdly, based on the above function, the scattered flux for any pixels can been obtained by the interpolation method. This method will save a lot of computing time and memory obviously.
 2.2 Scattering Correction Method
Taking the scattering effect into account as described in the previous sections, the following scattering correction and evaluation algorithm is proposed.

(1) Substituting the detected flux image into Eqs.(4) and (5) as an initial guess of the uncollided photon component yields an initial density distribution of the investigated object.

(2) Compute the scattering flux of the image, including both first and multiple scattered components, according to the information reconstructed in step (1) and the method described previously.

(3) Update the uncollided component of the flux image by subtracting the scattering contribution obtained in step (2) from the detected total photon flux.

(4) By substituting the updated uncollided photons information from step (3) into equations. (4) and (5) provides an updated density distribution can be derived.
When steps 24 are performed iteratively until the density distribution converges, the information of the image reconstruction will become more and more accurate.
3. NUMERICAL RESULTS
Based on the above model, an
I
terative
P
rocedure for image rest
OR
ation code (IPOR) for 3D arbitrary geometry was developed. In this section, two test problems are presented to validate the reliability and capability of the scattering correction method. Each of the two test cases includes two parts. The first part estimates the accuracy and efficiency of the scattering flux calculation method, while the second part reconstructs the density distribution of the sample. The correction effect is evaluated by comparing the results before and after scattering correction. All calculations are performed on a personal computer with an 3.2GHz Intel(R) Core™ i5 CPU. The convergence criterion for relative density error between two iterations is 10
^{4}
. Typically, the number of iterations is around 15 for those cases with heavy metal materials and less than 5 for those cases without heavy metal materials.
 3.1 Case 1: Irregular 3D Cylindrical Geometry Problem
This is an irregular 3D cylindrical geometry problem as shown in
Fig.6
and the height of the sample is 2 cm in
Schematic Diagram of the Irregular 3D Cylindrical Geometry Problem
the Z direction. It contains two types of materials as shown in
Table 2
. The ray source is located at (0.0, 20.0,0.0) and emits photons with energy of 20MeV. The detectors are located along the plane of
y
=20.0. It is notable that the unit of length in this paper is cm.
 3.1.1 The Accuracy and Efficiency of the Scattering Flux Calculation Method
Fig. 7
presents the scattered fluxes of the detectors provided by IPOR and MCNP. Taking the results from MCNP as a reference, the largest relative error of IPOR is 6.7%, where the statistical estimated relative error of MCNP is in the range of 1.0 × 10
^{4}
~ 2.0 × 10
^{4}
. The computational time of MCNP is 153.35 minutes, but that of
Material Information Data for Case 1
Material Information Data for Case 1
The Scattered Fluxes Calculated by IPOR and MCNP for Case 1
IPOR is only 7.99 minutes. It is demonstrated that the computational efficiency of IPOR is much higher than that of MCNP.
 3.1.2 The Effect of the Scattering Correction
The density of each material region within the sample are rebuilt using the image signal of the detectors (3.0, 20.0, 1.0) and (3.0, 20.0, 1.0) and supposing that the initial densities are all 0.0g/cm
^{3}
.
Table 3
presents a comparison of the reconstructed results before and after the scattering correction. It can be seen that the scattering correction considerably reduces the error introduced during the determination of the sample densities. The largest relative error after the correction is 0.46% , which is much smaller than the error of 7.58% before the correction. It means that the application of this correction can improve the rebuilding accuracy significantly. The computing time of this case is 141.65s which is pretty short compared to the direct transport simulations.
 3.2 Case 2: Spherical Shell Geometry Problem
Case two considers a spherical shell geometry problem as shown in
Fig.4
and as listed in
Table 4
. The objective contains lots of heavy materials. In this case, the detector is a flatpanel detector with a 28 × 28 pixel
The Comparison of the Reconstruction Results with and without Scattering Correction
The Comparison of the Reconstruction Results with and without Scattering Correction
Schematic Diagram of the Flatpanel Detector
array. The distance between the detector and the center of the object is 100 cm. Because the sample is rotationally symmetric, oneeighth of the image is selected to compute and study as the shaded part shown in
Fig.8
.
 3.2.1 The Accuracy and Efficiency of the Scattering Flux Calculation Method
The reference solution is also given by MCNP.
Fig. 9
shows the relative error distribution of the scattered flux in the image with the largest value of 2.9%. The statistical estimated relative error of MCNP is in the range of 1.0 × 10
^{4}
~ 7.5 × 10
^{3}
. It can been seen that IPOR provides satisfying accuracy. The computational time of MCNP is 3692.15 minutes, while that of IPOR is only 0.142 minutes. In this scenario, IPOR is superior to MCNP due to the fact that it reduces the computing time by a factor of four
Material and Geometry Information Data for Case 2
Material and Geometry Information Data for Case 2
The Relative Error Distributions of the Scattered Fluxes in the Shaded Image by Comparing IPOR with MCNP for Case 2
Radiography of Case 2
degrees of magnitude. The superiority of the computational efficiency seems much more significant in this case than in case one because in this case the detector is a flatpanel detector with 28 × 28 pixels array. All of the pixels are calculated in MCNP, but in IPOR only the selected typical pixels are calculated and the scattering distribution is obtained using least square fitting and interpolating. Thus the computational efficiency will be improved more significantly if there are more pixels in the image.
 3.2.2 The Effect of the Scattering Correction
The density distribution of the sample is rebuilt based on the image signal.
Fig.10
shows the radiography of this
The Ratio of Uncollided Photons to Collided Photons for Case 2
case. It can be seen obviously that the scattered photons blur the radiographies seriously. In the region near to the center of the image, the noise overwhelms the useful information. The ratio of the uncollided photons to collided photons is presented in
Fig.11
. It can be found that in some regions the values are even less than 10
^{3}
. It is very difficult to reconstruct the density distribution accurately in such strong scattering situations.
The density reconstructions in terms of different scattered photons percentages are studied and the rebuilt results before and after correction are presented in
Fig.12
. It can be seen that the reconstruction without correction is approximately
The Comparison of the Reconstructed Densities before and after Correction
The Average Relative Error of the Reconstructed Density for the Middle Region Material
The Average Relative Error of the Reconstructed Density for the Middle Region Material
distorted for the strong scattering sample.
Table 5
shows the average relative error of the reconstructed density for the middle region material and it can been seen that the rebuilding accuracy is improved by one degrees of magnitude after correction. The computing time for this case is only 160.16s.So it can been concluded that the scattering correction method not only improves the rebuilding accuracy significantly but also has a high computational efficiency.
4. CONCLUSIONS
In photon radiography, the application of scattering correction can reduce or even eliminate the disruption of the collided photons in the image, thus improving the reconstruction accuracy. The correction is based on the iterative reconstruction and scattering correction.

(1) In this paper, the uncollided ray flux and the single scattered ray flux are precisely calculated along the characteristics and the multiple scattered ray flux is corrected using numerical coefficients that have been precalculated. Compared with the Monte Carlo simulation, this method has both a very high computational efficiency and a satisfying accuracy.

(2) Based on the customization of AutoCAD, the geometry description and ray tracing for arbitrary geometry is carried out to greatly expand the capability of handling complex geometries.

(3) The numerical results indicate that the scattering correction can improve the reconstruction accuracy and the computational efficiency.
This scattering correction method has a high accuracy, high efficiency and expanded treatment capability for 3dimensional arbitrary geometry. So it provides a good scattering correction method for 3D conebeam flash radiography.
Acknowledgements
This research was carried out under the financial support of the National Science Foundation of ChinaNSAF (Grant No.10976021) and a grant from the National Magnetic Confinement Fusion Science Program (Approved number 2010GB111007).
Kardjilov N
,
Beer FD
,
Hassanei R.
(2005)
“Scattering corrections in neutron radiography using point scattered functions”
Nuclear Instruments and Methods in Physics Research A
542
336 
341
DOI : 10.1016/j.nima.2005.01.159
Hassanein R
,
Lehmann E
,
Vontobel P.
(2005)
“Methods of scattering corrections for quantitative neutron radiography”
Nuclear Instruments and Methods in Physics Research A
542
353 
360
DOI : 10.1016/j.nima.2005.01.161
Boone JM
,
Seibert JA.
(1988)
“Monte Carlo simulation of the scattered radiation distribution in diagnostic radiology”
Medical Physics
15
(5)
713 
720
DOI : 10.1118/1.596185
Akinao S.
(2002)
“Calculation of gammaray buildup factors up to depths of 100MFP by the method of invariant embedding”
Journal of Nuclear Science and Technology
39
(5)
477 
486
DOI : 10.1080/18811248.2002.9715225
Li BY
,
Shi JJ
,
Liu J
(2005)
“Numerical simulation of distribution of scattered exposure and reduced scatter in flash radiographic system”
High Power Laser and Particle Beams
17
(5)
788 
792
Shi JJ
,
Li BY
,
Liu J
(2006)
“Analytic determination of scatter exposure for radiography”
High Power Laser and Particle Beams
18
(7)
1211 
1214
Peng GH
,
Cai XH
,
Qiao NS
(2007)
“Scattering correction for Xray TICT in testing composites workpiece”
Non Destructive Testing
29
(4)
181 
184
Qian L.
(2002)
“Compensation for Compton scatter effects in TI201 spectrum”
Chinese Journal of Computational Physics
(in Chinese)
19
(1)
38 
42
Wang MQ
,
Cao LZ
,
Wu HC
“Study on scattering corrections in flash radiography”.
Atomic Energy Science and Technology
(in Chinese), to be published.
Cheng YX
,
Cao LZ
,
Wu HC
(2010)
“Direct solutions for inverse transport problems”
Nuclear Power Engineering
(in Chinese)
31
(S2)
38 
42
Chen QC
,
Wu HC
,
Cao LZ.
(2008)
“Auto MOCA 2D neutron transport code for arbitrary geometry based on the method of characteristics and customization of AutoCAD”
Nuclear Engineering and Design
238
(10)
2828 
2833
Liu ZY
,
Wu HC
,
Cao LZ
(2011)
“A new threedimensional method of characteristics for neutron transport calculation”
Annals of Nuclear Energy
38
447 
454
DOI : 10.1016/j.anucene.2010.09.021