Ting Su

Valérie Kaftandjian
email: valerie.kaftandjian@insa-lyon.fr

Philippe Duvauchelle

Yuemin Zhu

Ting Author

Su

Classification-based material decomposition method for photon counting-based spectral radiography: application to plastic sorting

Keywords: X-ray, radiography, photon-counting detector, plastic sorting

Dual-energy X-ray radiography has been commonly used for materials separation. However, its performance is limited, especially for the separation of close materials or the identification of multiple materials. To cope with this problem, we propose to investigate the material decomposition ability of spectral radiography based on photon-counting detector. The latter having energy resolving capability can provide spectral information of several energy bins and thus enables selective imaging of multiple materials. In this framework, a classification-based patchwise regularized decomposition method was proposed to gain better differentiation between materials. It consists of performing several decompositions with reduced number of materials in the basis and classifying these decompositions using their cost function values. The results on simulations showed that, in the presence of Poisson noise, the method without classification can separate acrylonitrile-butadiene-styrene (ABS) from three kinds of flame retardants (FRs: brominated FR, chlorinated FR and phosphorus FR), but that the type of FR cannot be identified. With the classification technique, ABS and three kinds of FRs can be both separated and identified at the same time when the

Introduction

X-ray radiography is an imaging technique that allows to view the internal structure of an object. Since the transmitted image is a 2-D projection of the 3-D object, useful information can be hidden because of physical overlapping.

Dual-energy radiography using either different kVps or dual-layer detector is able to obtain two radiographs or images, one at higher energy and the other at lower energy. By making weighted subtraction of these two images [START_REF] Brody | Dual-energy projection radiography: initial clinical experience[END_REF][START_REF] Vock | Dual energy subtraction: principles and clinical applications[END_REF] or by performing a two-material decomposition [START_REF] Brody | A method for selective tissue and bone visualization using dual energy scanned projection radiography[END_REF], elimination of overlaying materials and enhancement of the selective target material can be achieved.

Dual-energy imaging has broad applications, such as rib suppression in chest radiography in medical domain [START_REF] Manji | Comparison of dual energy subtraction chest radiography and traditional chest x-rays in the detection of pulmonary nodules[END_REF][START_REF] Zhou | Visualization of coronary artery calcium in dual energy chest radiography using automatic rib suppression[END_REF], threat detection in security inspections [START_REF] Singh | Explosives detection systems (eds) for aviation security[END_REF][START_REF] Pashby | Radiation detection and dualenergy x-ray imaging for port security[END_REF][START_REF] Lee | Efficient material decomposition method for dual-energy x-ray cargo inspection system[END_REF] and waste sorting in industrial applications [START_REF] Montagner | Dual energy radioscopy applied to waste sorting[END_REF][START_REF] Gundupalli | A review on automated sorting of source-separated municipal solid waste for recycling[END_REF][START_REF] Duvillier | Inline multi-material identification via dual energy radiographic measurements[END_REF]. However, its material discrimination ability is limited and materials with close atomic number (Z) can barely be separated.

Thanks to the recent development of X-ray detector technology, a variety of pixelated photon-counting detectors (PCD) with imaging capability have been developped [START_REF] Ogawa | Development of an energy-binned photoncounting detector for x-ray and gamma-ray imaging[END_REF][START_REF] Babichev | Photon counting detector for the personal radiography inspection system "sibscan[END_REF][START_REF] Gimenez | Development of a schottky cdte medipix3rx hybrid photon counting detector with spatial and energy resolving capabilities[END_REF][START_REF] Zambon | Spectral response characterization of cdte sensors of different pixel size with the ibex asic[END_REF]. By setting different thresholds (energy bins), each pixel of PCD is capable to either discriminate the transmitted photons corresponding to the selected energy bins, or firstly count the photons above the thresholds and then make subtraction to obtain the counts of each energy bin.

With such spectral information, it is possible to set apart multiple components of the object with a single scanning and obtain improved performance in the aforementioned areas [START_REF] Beldjoudi | An optimised method for material identification using a photon counting detector[END_REF][START_REF] Fredenberg | Spectral and dual-energy x-ray imaging for medical applications[END_REF][START_REF] Xu | Systematic implementation of spectral ct with a photon counting detector for liquid security inspection[END_REF][START_REF] Maruhashi | Evaluation of a novel photon-counting ct system using a 16-channel mppc array for multicolor 3-d imaging[END_REF][START_REF] Si-Mohamed | Review of an initial experience with an experimental spectral photon-counting computed tomography system[END_REF]. Most researches about the material decomposition ability of PCD concentrate on computed tomography (CT), called spectral CT. In the work of Schlomka et al. [START_REF] Roessl | K-edge imaging in x-ray computed tomography using multi-bin photon counting detectors[END_REF][START_REF] Schlomka | Experimental feasibility of multi-energy photon-counting k-edge imaging in pre-clinical computed tomography[END_REF], the measured attenua-tion coefficients in spectral CT were decomposed into photoelectric absorption, Compton scattering and K-edge components, showing that dual contrasts can be well separated and quantified. Such K-edge technique is confined to the discrimination of high-Z materials with their K-edge inside the detection energy range. Other works [START_REF] Le | Least squares parameter estimation methods for material decomposition with energy discriminating detectors[END_REF][START_REF] Alessio | Quantitative material characterization from multi-energy photon counting ct[END_REF][START_REF] Zeng | Penalized weighted least-squares approach for multienergy computed tomography image reconstruction via structure tensor total variation regularization[END_REF][START_REF] Xie | Material decomposition in x-ray spectral ct using multiple constraints in image domain[END_REF] considered the attenuation as a combination of several basis materials, whose distribution can therefore be quantified through decomposition process. Such approach is particularly suitable for the discrimination of lighter materials that have much lower K-edge energy beyond the detection energy range. In our previous work [START_REF] Su | A spectral x-ray ct simulation study for quantitative determination of iron[END_REF][START_REF] Su | Quantitative material decomposition methods for x-ray spectral ct[END_REF], a patchwise regularized decomposition method based on basis materials combination was proposed for iron determination; the method allows iron to be quantitatively separated from calcium, potassium and water. However, potassium and calcium cannot be discriminated from each other since they have close Z numbers and hence close attenuation properties. The above material decomposition methods initially proposed for CT can be applied to performing radiography, thus resulting in the so-called PCD-based spectral radiography named as spectral radiography in the following for simplicity. In the present work, we propose a new material decomposition method for spectral radiography and investigate its feasibility for plastic sorting.

Waste electrical and electronic equipment (WEEE) has been increasing rapidly due to the development of electronic industry. In European Union, the amount of WEEE generated in 2005 is 9 million tons and this number is supposed to grow to 12 million by 2020 [START_REF] Beigbeder | Study of the physico-chemical properties of recycled polymers from waste electrical and electronic equipment (weee) sorted by high resolution near infrared devices[END_REF]. WEEE covers a wide variety of products such as lamps, hairdryers, computers, TV-sets, fridges and cell phones. This kind of waste contains various compositions among which polymers play an important part [START_REF] Cui | Mechanical recycling of waste electric and electronic equipment: a review[END_REF]. Acrylonitrile-butadiene-styrene (ABS) is a major component among all the polymers. Due to the presence of electronic power of EEE, the plastic materials should meet high fire safety standards. However, this is not possible to realize with pure polymers, therefore flame retardant (FR) are added to change flammability of plastics and increase the fire resistance [START_REF] Peeters | Closed loop recycling of plastics containing flame retardants[END_REF]. Halogenated organic compounds (e.g., brominated aromatic compounds, chlorinated paraffins and alicyclic compounds) and phosphorus compounds (e.g. organophosphates, halophosphates, phosphine oxides and red phosphorus) are two important FRs [START_REF] Fromme | Brominated flame retardantsexposure and risk assessment for the general population[END_REF]. Recycling of plastics from WEEE is challenging because of the existence of FRs (especially the brominated and chlorinated FRs), which might result in serious environmental pollution. Therefore, sorting technique is important for subsequent treatment and recycling of WEEE to avoid major environmental and health problems. As mentioned before, some researches using dual-energy radiography for waste management have been reported. For example, Montagner et al. [START_REF] Montagner | Dual energy radioscopy applied to waste sorting[END_REF] showed that heavy FR containing bromine (Br) is easier to be detected, but that chlorinated (Cl) FR can be sorted in the same family as Br, while the lighter one containing phosphor (P) remains undetectable when material thickness is less than 10 mm. To overcome these limitations, PCDbased spectral radiography appears an interesting solution due to its ability of decomposing multiple components of an object with a single scanning.

The purpose of the present work is to investigate the method of material decomposition for spectral radiography to discriminate plastic material from several FRs and further identify the type of FR. Based on our preliminary results reported in the ECNDT 2018 conference [START_REF] Su | Plastic sorting by x-ray radioscopy with photon counting detector[END_REF], the present paper aims to propose a new material decomposition method, which consists of performing several decompositions with reduced number of materials in the basis and classifying these decompositions using their cost function values. Both simulation and experimental results will be demonstrated. This paper reuses a substantial part from one of the co-author's own thesis [START_REF] Su | Quantitative material decomposition methods for x-ray spectral ct[END_REF] with permission. More precisely, the reused contents appear mainly in sections I to IV and in Appendix.

Classification-based decomposition method

In our previous work [START_REF] Su | A spectral x-ray ct simulation study for quantitative determination of iron[END_REF], a patchwise regularized decomposition (PRD) method has been proposed, which divides the projection images into small patches and performs decomposition on each patch. The limitation of this method is to distinguish materials with close properties, and the performance of decomposition decreases when increasing the number of materials in the basis.

To cope with the difficulty, the idea here is to perform several decompositions with reduced number of materials in the basis and classify these decompositions using their cost function values, thus leading to the so-called classification-based PRD (CPRD) method. CPRD method also works in the way of patch-by-patch decomposition, but for several groups of materials. Figure 1 illustrates the procedure of this method for one patch. The principle is to firstly choose several groups of basis materials ( Group 1, Group 2, ... Group g, ..., Group G), in which the same material can belong to several groups. Then, for each patch, we perform G times decomposition by carrying out PRD on each of basis material groups. For a Group g, the output is a set of decomposed density integrals at each pixel of the patch and a unique cost function value. After this step, we select the basis material group having the smallest cost function and we adopt the set of concentration values obtained from the decomposition based on this basis material group, while the concentration of other materials inside this patch are set to 0. Finally, by repeating the above procedure for all the patches, we obtain the final image of each material.

More precisely, we consider the attenuation µ( -→ x , E) as a combination of M basis materials, where -→ x represents location and E represents energy. The basis materials can be chosen according to the prior knowledge of the scanned objects. We then have:

µ( -→ x , E) = M ∑ α=1 ρ α ( -→ x )µ mα (E), (1)
where ρ α ( -→ x ) denotes the density of material α at point -→ x and µ mα (E) the mass attenuation coefficient of material α at energy E. According to Beer-Lambert law, the expected number of photons λ i , in energy bin B i (i = 1, 2...N, N is the total number of energy bins) can be expressed as:

λ i = E f (i) ∑ E=Es(i) D(E)N 0 (E)exp[- ∫ µ( -→ x , E) ds] = E f (i) ∑ E=Es(i) D(E)N 0 (E)exp[- ∫ M ∑ α=1 ρ α ( -→ x )µ mα (E) ds] = E f (i) ∑ E=Es(i) D(E)N 0 (E)exp[- M ∑ α=1 µ mα (E)P α (sx, sy)], (2)
with

P α (sx, sy) = ∫ ρ α ( -→ x ) ds, (3)
where E s (i) and E f (i) denote respectively the start and final energies of bin B i , D(E) is the detector absorption efficiency, N 0 (E) is the number of photons in the initial spectrum at energy E, and P α (sx, sy) is the density integral that is line integral of object density ρ α ( -→ x ) along the measured projection path, with (sx, sy) representing the index of detector pixels.

To estimate P α (sx, sy), we build an objective function that combines the log-least squares [START_REF] Su | Quantitative material decomposition method for spectral ct imaging[END_REF][START_REF] Su | Material decomposition for spectral ct: application to calcium and iodine identification[END_REF], and a regularization term R(P C α ) for each small patch C on the attenuation image to reduce the effect of noise and enforce smoothness:

P C α (sx, sy) = arg min P C α (sx,sy) { ∑ (sx,sy)∈C N ∑ i=1 [ln(λ i (P C α )) -ln(m C i )] 2 + rR(P C α )}, (4)
where m C i is the measured number of photons in energy bin B i and within patch C. r denotes the relaxation parameter. R(P C α ) is the sum of the L2 regularizations of gradient images of P C α :

R(P C α ) = M ∑ α=1 || ▽ P C α || 2 2 = M ∑ α=1 ∑ (sx,sy)∈C (sx-1,sy)∈C (sx,sy-1)∈C {[P C α (sx, sy) -P C α (sx -1, sy)] 2 +[P C α (sx, sy) -P C α (sx, sy -1)] 2 }. ( 5
)

Description of the spectral radiography simulation

We use the Virtual X-ray Imaging (VXI) software [START_REF] Duvauchelle | A computer code to simulate x-ray imaging techniques[END_REF], which was developed 125 in INSA Lyon, to simulate spectral radiography imaging. This software allows to define multiple parameters of the object to scan, including X-ray source, detector and system geometry.

(α) = ρ(mixture) × ω(α).
For example, ρ eff (Br) = ρ(mixture) × ω(Br) = 1060 × 8.82% = 93.5 mg/cm 3 ; 1.

ρ eff (ABS) = ρ(mixture) × ω(ABS) = 1060 × (1 -15%) = 901 mg/cm 3 .

ABS-FRs phantom

We simulate a phantom containing multiple polymer+FR mixtures. As men-130 tioned above, ABS is a major component among the polymers used in electrical and electronic equipments. Therefore, we selected it as the polymer material for investigation. We also chose three kinds of commonly used flame retardants, including the brominated, chlorinated and phosphorus FRs. They are respectively tetrabromobisphenol A (TBBPA), dechlorane plus (DDC-CO) and resorcinol bis(diphenyl phosphate) (RDP). Three ABS-FR materials were obtained by mixing each FR with ABS at mass percentage of 15%. Detailed information of these materials is summarized in Table 1.

The simulated phantom is composed of multiple cubes with height of 10 mm, width of 10 mm and different thicknesses. As shown in Figure 2, each column of cubes are of the same material (denoted on the top) and each row of cubes are of the same thickness (denoted on the left). Figure 3 plots the mass attenuation coefficients of the components contained in the phantom: ABS, Br, Cl and P, which will be used to distinguish different ABS-FR materials. From the decomposition method described in the above, the finally obtained P α for material α is the line integral of material density ρ α (Equation 3). In this case, with known thickness d of each cube, theoretical values for P α at position

P α (sx, sy) = d(sx, sy) × ρ α (sx, sy), (6)
where ρ α equals the effective density calculated in Table 1. The simulated system used a point X-ray source of tungsten target material and the target angle was 17 • . According to reference [START_REF] Montagner | Nouvelles méthodes de tri des déchets par rayons x[END_REF], tube voltage was set to be 100 kVp, tube current 15 mA and exposure time 1 s. The X-ray spectrum was computed based on Birch & Marshall model [START_REF] Birch | Computation of bremsstrahlung x-ray spectra and comparison with spectra measured with a ge (li) detector[END_REF] without filtration. A 90×112 CdTe detector array with pixel size of 0.5 mm × 0.5 mm and thickness of 3 mm was simulated. Six energy bins were set to be evenly distributed from 30 keV to 90 keV. We assumed that the detector has ideal energy resolution, detector absorption efficiency was simulated as the fraction of photons absorbed by the detector material. The distance from X-ray source to phantom center was 2000 mm and the distance of detector to phantom center was 3 mm. Figure 4 is a schematic view of the simulated spectral radiography system. Photonic noises were simulated by simple Poisson distribution with means equaling the expected numbers of photons acquired for each energy bin [START_REF] Whiting | Properties of preprocessed sinogram data in xray computed tomography[END_REF][START_REF] Ma | Variance analysis of x-ray ct sinograms in the presence of electronic noise background[END_REF]. This process was implemented in Matlab using command poissrnd.

Simulation Results

We have simulated two sets of spectral radiography acquisition data: with and without photon noise. To better evaluate the performance of the proposed CPRD method, it was also compared with PRD method without classification.

Decomposition results of PRD method

Regarding the components of the phantom, it is normal to consider a four basis material decomposition: ABS, Br, Cl and P. However, our experiment results show that P and Cl can barely be separated from each other since they have too close atomic number and hence close attenuation properties. Therefore, we decompose the radiographic images into three basis images corresponding to ABS, Br and Cl and expect that P will also be present in the Cl basis image.

The patch size is set to be 2×2 pixels.

For simulated acquisition without noise, the relaxation parameter r is set to 0 since it is not necessary to enforce smoothness within patches without noise.

The decomposition results are shown in Figure 5(a). It is observed that: , the obtained images are less noisy than those in the former case, however, the separation of different materials is not enhanced. Another method is thus necessary.

Decomposition results with the new CPRD method

Instead of the 3-material (ABS, Br and Cl) decomposition using PRD method, we perform two independent decompositions with two groups of basis materials, group 1 is ABS + Br, and group 2 ABS + Cl. We give up a third group of ABS + P because it generates comparable cost function values with group 2, resulting in severe cross-talk between Cl and P. With the selected two groups of decomposition, by comparing the cost function values (fval1 and fval2) at the end of each decomposition, one of the two decompositions is chosen for a given patch. For example, if fval1<fval2, we choose the results of the decomposition that decomposes the data into ABS and Br, consequently, the values of Cl in this patch will be set to 0. After all patches being considered, full images of ABS, Br and Cl are obtained. Figure 6 shows the decomposition results of CPRD method and a reference method in the same condition as that in the above method without classification. The reference method employs a widely used least squares error objective function to solve the decomposition problem pixel by pixel [START_REF] Lee | Quantitative material decomposition using spectral computed tomography with an energy-resolved photon-counting detector[END_REF]. 2. The dash lines on images of (a) mark the position where 1-D profiles are plotted in later analysis.

In the noise-free condition (Figure 6 In contrast, the reference method fails to separate Br and Cl. As can be seen in 6(d), the cubes containing Br also appear in the Cl basis image, no matter what thickness the sample has. Therefore, the CPRD method outperforms the reference method for material decomposition.

A more quantitative comparison is done in Figure 7, where the profiles of computed P α with PRD and CPRD method are plotted together with the theoretical values (along dash lines in both Figure 5 and6). It is observed that the CPRD method yields closer results to true values with respect to the PRD method especially for the Cl basis curve, where ABS-TBBPA cubes are mistaken as Cl containing FR by PRD method. It is noteworthy that the theoretical values of P α for the ABS-RDP cube in the Cl basis image are set to 1.72mg/cm 2 (corresponding to Table 2, the P α of P in ABS-RDP with thickness of 1 mm) instead of 0, because we expect that P will appear in the Cl basis image.

So, the CPRD method has better decomposition ability than PRD for the sorting application, no matter when Poisson noise is or is not present. ABS and Br containing FR can be identified in their corresponding basis images.

Furthermore, when the cube thickness is increased to 2 mm or 4 mm, the decomposition results of Poisson noise condition become much closer to those of noise-free condition, where ABS and the FRs containing Br, Cl and P can be identified simultaneously. 2)for perfect decomposition. Blue curves with triangle marker represent the calculated values using PRD method and the magenta curves with cross marker using CPRD method.

Preliminary experimental results

In order to evaluate the performance of our proposed method on physical data, we carried out radiography measurements using a CdTe photon-counting Because of the linear structure of the detector, we have acquired single line data of 256 pixels. Therefore, the patchwise regularization, which requires 2 dimensional gradient, can not be directly applied. We hence adjusted our CPRD method by removing the regularization term and using patch size of 1x1 for experimental validation. 3 groups of basis materials were selected: group 1 "ABS + Br", group 2 "ABS + Cl", group 3 "ABS + Ca". However, a mistaken peak of PP with 20% CaCO 3 also appears, indicating that Ca probably cannot be separated from Cl. Also, the sample containing 5% Cl remains undetected.

In short, the proposed method achieves to identify ABS, Br as well as Cl with relatively higher concentration, but materials with close atom numbers, such as Ca and Cl remain unseparated.

Discussions

We have proposed a classification-based material decomposition method, which achieves to identify several different plastic materials in both simulation and experimental validations. There are several parameters and factors that may influence the decomposition performance as discussed below.

Choice of basis materials. During simulation, we have found that using basis materials that have close chemical formulas to those actually present in the object can improve the decomposition performance. For example, using the same simulation data with Poisson noise, besides the above presented results, we have also performed a 3-group decomposition with basis materials like, Group 1 (ABS + TBBPA), Group 2 (ABS + DDC-CO) and Group 3 (ABS + RDP). The results showed that with cube thickness of 2 mm and 4 mm, ABS and three FRs can be identified according to their corresponding decomposition image. That is to say, DDC-CO and RDP can be distinguished from each other directly by decomposing them into their corresponding basis, instead of by their density difference in the same basis (Cl) image, as described in the above results of the CPRD method. We have not chosen these groups of basis materials with better performance because they are too specific. However, for the application where the components of object are well known, appropriate choice of basis materials can yield better decomposition performance.

Patch size.

In the present study, we have used a patch size of 2 × 2 pixels.

To investigate the influence of patch size on decomposition performance, we plot in Figure 10 the curve of quantification error versus patch size, where quantification error represents root-mean-square-error (RMSE) of decomposed Br basis image. As observed, 2 × 2 patch size gives the smallest RMSE. Quantification error increases when patch size is larger than 2. This could be caused by the fact that large patch size results in over-smoothed image [START_REF] Zhang | Low-dose ct reconstruction via l1 dictionary learning regularization using iteratively reweighted least-squares[END_REF]. As a result, since there exist sharp grid-like edges throughout our entire image, over-smoothing of these edges leads to larger quantification errors. Relaxation parameter r. Although the separation between materials is not enhanced (for PRD method) or slightly enhanced (for CPRD method) by choosing r according to L-curve method, it is still meaningful that noise (indicated by the variation of pixel values) of decomposition image can be decreased, which is obvious in Figure 5. We have applied L-curve method to certain region of the radiographic images, and used the same r for all patches. However, r should perhaps be changed for different patches if the signal-to-noise ratio varies. Therefore, further improvement could be achieved by tuning the relaxation parameter in each patch.

Placement of energy bins. In the present work, we have used 6 uniformly distributed bins within detecting energy range to evaluate the performance of our proposed method in a general way. However, according to some researches [START_REF] Leng | Noise reduction in spectral ct: Reducing dose and breaking the trade-off between image noise and energy bin selection[END_REF][START_REF] Wang | Optimal energy thresholds and weights for separating materials using photon counting x-ray detectors with energy discriminating capabilities[END_REF], appropriate selection of energy bins may have large influence on decomposition performance, especially when the total number of bins are small. In Detector response. In the present work, simulation validation only considers the absorption efficiency of the detector. However, in reality, photon counting detector's response function can be modeled by mixture of Gaussian functions allowing taking into account several artifacts such as electronic noise, charge sharing, pulse pile-up, K-escape and Compton scattering [START_REF] Schmidt | Ct energy weighting in the presence of scatter and limited energy resolution[END_REF]. Electronic noise is an important factor that limits the energy resolving capability of detector, it results in spectral resolution degradation and potential false counts especially at low energy. [START_REF] Fredenberg | Energy resolution of a photon-counting silicon strip detector[END_REF]. Charge sharing effect occurs when the charge cloud generated from X-ray interaction is shared by neighboring pixels and results in two or more counts in several pixels and with lower energies. This effect is especially serious when pixel size is less than 0.5 mm [START_REF] Ren | Tutorial on x-ray photon counting detector characterization[END_REF]. Pulse pileup takes place at high count rate; two nearly simultaneously incident photons are counted as one, leading to loss of count and wrong recorded energy. This effect can be reduced by controlling incident count rate [START_REF] Fredenberg | Energy resolution of a photon-counting silicon strip detector[END_REF]. K-escape effect is due to the characteristic X-rays induced by photoelectric effect. It can lead to different consequences: the recorded energy is lowered by K-shell energy, or two counts are generated, or still pile-up effect is induced [START_REF] Taguchi | Vision 20/20: Single photon counting x-ray detectors in medical imaging[END_REF].Compton scattering refers to Compton effect-induced scattered photons that are detected by an adjacent pixel or leave the PCD completely, resulting in long tail in lower energy. All these effects lead to distorted detector response and reduction of energy resolution. The degraded performance of material decomposition caused by discrepancy between theoretical and measured data can be addressed by two approaches. The first one is to characterize detector response by synchrotron radiation, radioactive isotopes or x-ray fluorescence [START_REF] Ding | Characterization of energy response for photon-counting detectors using x-ray fluorescence[END_REF], and then incorporate this information into material decomposition [START_REF] Liu | Material reconstruction for spectral computed tomography with detector response function[END_REF]. The second one is to calibrate material's linear attenuation coefficient by scanning samples with known components and thickness, so as to take into consideration the influence of detector response. The second approach is more flexible, easy to accomplish and efficient according to our preliminary test on experimental data. Though this strategy was not used for the demonstrated results here, we can still observe pretty good decomposition performance, indicating that our proposed method is able to process real experimental data.

There are certain limitations of the proposed method. Firstly, the CPRD method is not enough robust to noise when the object thickness is small. Secondly, the decomposition of the two FRs with similar properties (DDC-CO and RDP), even in the noise-free case, depends on their density difference in the Cl image. Therefore, in case where the thickness information of objects is unknown, or the two materials overlap in the X-ray direction, the decomposed density will be largely influenced, yielding the materials indistinguishable. Although an optional solution has been mentioned in the above, which consists in using directly DDC-CO and RDP as basis materials, but in a more generalized

Conclusions

To investigate the ability of PCD-based spectral radiography for plastic sorting, we have proposed a classification-based decomposition method and evaluated its performance through both simulation and experimental validations.

The simulation results showed that the proposed CPRD method is able to distinguish ABS and the FRs containing Br, Cl and P at the same time if the cube thickness is as large as 2 mm or 4 mm. Moreover, the quantified material densities agree well with theoretical values. Experimental validations further confirm that ABS, Br and Cl can be separated from each other, but that materials with close atom numbers, i.e. Ca and Cl remain unseparated. Despite some limitations, our results suggest the possibility that PCD-based spectral radiography can be used to sort different types of plastic materials.

In addition to plastic sorting presented in the present work, the proposed method can be readily applied to other applications such as industrial and medical ones. For example, in explosive detection, the explosive substances have close X-ray attenuation properties with some common materials such as sugar and polyethylene, a spectral radiography scan followed by the CPRD decomposition may improve their separation. Furthermore, the proposed method can be adapted to spectral photon-counting CT by adding a reconstruction step to transform the projection density integral information into cross-sectional density distributions. Thus spectral CT inspections could benefit from the proposed method too, such as abdominal imaging and atherosclerosis imaging, where material-specific imaging is of great importance for characterizing lesions or high risk plaques.

Figure 1 :

 1 Figure 1: Flowchart of the decomposition method with classification for one patch. fvalg represents the cost function value obtained by the gth decomposition with basis materials of Group g.

Figure 2 :

 2 Figure 2: Illustration of the ABS-FR phantom used for spectral radiography imaging. Materials in the cubes are given in Table1.

Figure 3 :

 3 Figure 3: Mass attenuation coefficients (µm) of components (ABS, Br, Cl and P) within the detection energy range from 30 keV to 90 keV. µm of ABS is calculated according to formula µm = ∑ α ωαµmα, with ωα and µmα representing the weight fraction and mass attenuation coefficient of each element α of the compound. Data taken from [37].

Figure 4 :

 4 Figure 4: Scheme of spectral radiography system geometry. This figure is only a simplified schematic view of the simulation system where the dimensional information (D SC , D DC , size of phantom and detector) is not in real scale.

Figure 5 :

 5 Figure 5: Decomposition results (Pα) under different conditions using PRD method. (a): simulation without noise, r = 0; (b): simulation with Poisson noise, r = 0; (c): simulation with Poisson noise, r = 2.51 • 10 4 . The first column represents ABS basis, second column Br basis, and third column Cl basis image. The colorbars indicate decomposed Pα values in unit of mg/cm 2 , whose true values are given in Table2. The dash lines on images of (a) mark the position where 1-D profiles are plotted in later analysis.

Figure 6 :

 6 Figure 6: Decomposition results under different conditions using CPRD method and reference method: (a) CPRD method, simulation without noise, r = 0; (b) CPRD method, simulation with Poisson noise, r = 0; (c) CPRD method, simulation with Poisson noise, r = 2.51 • 10 4 ; (d) reference method, simulation with Poisson noise. The first column represents ABS basis, second column Br basis, and third column Cl basis image. Each column shares the same colorbar which indicates decomposed Pα values in unit of mg/cm 2 , whose values are given in Table2. The dash lines on images of (a) mark the position where 1-D profiles are plotted in

4 (

 4 (a)), FRs containing ABS and Br are identified respectively by CPRD method in the ABS and Br basis images as expected. Meanwhile, only the cubes of FR materials containing Cl and P appear in the Cl basis image, but they have significant density difference for cubes with the same thickness, therefore, they can be easily distinguished by the observer even if the concentration of FRs or thickness of cubes changes more or less.In the presence of Poisson noise, Figure6(b) shows the decomposition results of CPRD method when r = 0 (without regularization). Noise can be observed in decomposition images, however, material separation is better than with PRD method (compared with Figure5): only few pixels containing Cl are misidentified in the Br basis image and vice versa. The cross-talk situation improves with the increase of thickness. In Figure6(c), when r = 2.51 • 10 with regularization), the material separation is further enhance to some extent compared with Figure6(b): less pixels of ABS-DDC-CO cubes appear in the Br basis image, and only the ABS-DDC-CO and ABS-RDP cubes are highlighted in the Cl image when d = 2 mm or 4 mm, like the results in noise-free case.

Figure 7 :

 7 Figure 7: Performance comparison of PRD and CPRD methods in noise-free condition: 1-D profiles along the dash lines in Figures 5 and 6. ABS basis (top), Br basis (middle) and Cl basis (bottom). Black curves with point marker represent the theoretical density integrals of basis materials (Pα for thickness of 1 mm in Table2)for perfect decomposition. Blue curves

Figure 8 :

 8 Figure 8: Illustration of experimental specimen (left) and system geometry (right).

Figure 9 :

 9 Figure 9: Decomposition results of ABS (a), Br (b), Cl (c) and Ca (d). The green check mark indicates a good decomposition while the red cross mark indicates a false decomposition.

Figure 10 :

 10 Figure 10: RMSE of quantified Br densities versus patch size.

Figure 11 (Figure 11 :

 1111 Figure 11(b), we demonstrate the decomposition results using a different strategy to place 6 energy bins. More bins were placed in lower energy on purpose since lower energy carries more material-differentiation information. We can observe performance improvement compared to Figure 11(a). From the ROI marked by the red circles, Br was better recognized with less cross-talk in Cl basis. Moreover, the calculated RMSE of Br basis image decreased from 1.23 to 0.88, indicating less quantification errors.

  situation, close materials are still hard to be distinguished. Thirdly, the proposed CPRD method performs several decompositions into different groups of basis materials and then uses a classification technique to decide the right type of material. Although a decreased number of basis materials in each group re-duces the running time of each decomposition, compared with a normal one with multiple basis materials, but on the other hand, the computation time can be increased because of extra decomposition procedures. For instance, in the present study, the 2-group (two basis materials in each group) CPRD decomposition has about 36% longer running time than the 3-material PRD decomposition.

Figure A2 :

 A2 Figure A2: The L-curve as a function of residual norm for different values of the relaxation parameter r (a), and a magnified view of the blue circle region, where r varies from 10 4 to 10 5 (b).

Table 1 :

 1 Components of three ABS-FR materials used for the phantom. Chemical formula of FR C 15 H 12 Br 4 O 2 C 18 H 12 Cl 12 C 30 H 24 O 8 P 2 Cl and P, respectively 93.5 mg/cm 3 103.5 mg/cm 3 17.2 mg/cm 3

	Material(ABS-FR)	ABS-TBBPA	ABS-DDC-CO	ABS-RDP
	Density of material ρ(mixture)	1060 mg/cm 3 1060 mg/cm 3 1060 mg/cm 3
	Mass % (ω)of FR	15%	15%	15%
	Mass % (ω) of Br, Cl and P, respectively	8.82%	9.76%	1.62%
	ρ eff

*

of Br, * : ρ eff is the effective density calculated by ρ eff

Table 2 :

 2 Theoretical values of Pα for each material. Material ABS has two columns: the first is ABS pure, which corresponds to the material of the 4th column of cubes in the phantom shown in Figure 2, and ABS in mixture corresponds to the ABS compound in the other cubes of Figure 2.

	Thickness (mm)	ABS pure (mg/cm 2 )	ABS in mixtures (mg/cm 2 )	Br in ABS-TBBPA (mg/cm 2 )	Cl in ABS-DDC-CO (mg/cm 2 )	P in ABS-RDP (mg/cm 2 )
	0.3	31.80	27.03	2.81	3.11	0.52
	0.5	53.00	45.05	4.68	5.18	0.86
	1	106.00	90.10	9.35	10.35	1.72
	2	212.00	180.20	18.70	20.70	3.44
	4	424.00	360.40	37.40	41.40	6.88

(sx, sy) of the transmitted image is given by

Table 2

 2

	DDC=3 mm
	DSC=2000 mm	PCD
	X-ray source	Planar array

lists the theoretical values of P α for materials ABS, Br, Cl and P at various thicknesses.

150

3.2. Acquisition parameters and system geometry

Appendix A. Determination of r with L-curve method

The L-curve is a log-log plot of the norm of a regularized solution versus the norm of the corresponding residual, where the "corner" of the curve corresponds to the selected relaxation parameter [START_REF] Hansen | The l-curve and its use in the numerical treatment of inverse problems[END_REF]. The solution and residual norms correspond to the objective function value and R(P C α ) value in Equation 4. We choose a region of interest on the radiographic image for the determination of r, which is shown in Figure A1. The ROI contains 20×20 pixels. Given that the patch size is 2×2, there are 10 ×10=100 patches within the ROI. For a given r, we define the solution norm and residual norm as the average values calculated over the 100 patches.

We plot the L-curve with r varying from 10 -10 to 10 10 , as shown in Figure A2 (a). It is observed that the "corner" exists between r = 10 4 and r = 10 as the relaxation parameter for all patches.


N D Le

B Davier

P Dollfus

M Pala

A Bournel

J Saint-Martin

Full Band Monte Carlo simulation of phonon transfer at interfaces

Our home made Full-Band particle Monte Carlo is used to investigate the thermal interface conductance at Silicon/Germanium heterojunctions.

Introduction

As self-heating effect could degrade electronic device performances [1], understanding the nature of heat transport at the nanoscale is important for the CMOS industry [2]. In particular, the contribution of out-ofequilibrium phonons in heat transfer at interfaces remains an important issue.

Molecular Dynamics (MD) [3,4] and Non Equilibrium Green's function (NEGF) [5,6] simulations have been used to investigate this topics. The method of MD, provided important results to the study of heat transport as the complex structures of interfaces can be described at the atomic level in terms of roughness, strain and anharmonic phonon/phonon scattering mechanisms [7][8][9]. However, the size of the system that can be studied is limited to few thousands of atoms. Whereas MD only takes into account the classical nature of phonon transport, NEGF, by contrast, considers its quantum nature as interference and Bose-Einstein statistics can be captured. It is a very efficient method to calculate the transmission coefficient through interface, but is preferred for small and low-dimensional systems rather than micrometer long systems [10][11][12][13].

The particle Monte Carlo (MC) algorithm [14][15][16] is a powerful and complementary tool to investigate the heat transfer as the simulated system size can vary from the nanoscale to the microscale. MC solves the Boltzmann Transport Equation (BTE) for phonons without any assumption on the distribution function, thus the out-of-equilibrium transport regime is naturally taken into account. Besides, with convenient modeling, particle scattering mechanisms can be accurately described in the full phase space especially at the interfaces.

In this work, our Full-Band (FB) particle Monte Carlo [17] approach has been extended to treat the phonon transport at semi-transparent interfaces and used to extract and analyze the thermal interface conductance at Silicon/Germanium heterojunctions.

Full Band Monte Carlo

In the particle Monte Carlo method dedicated to phonon transport, the trajectories of a large number of phonons are randomly selected according to relevant scattering rates and phonon dispersion that is in our case a Full-Band description. "Full Band Monte Carlo" means that all the phonon states in the first 3D Brillouin zone are considered as detailed in Ref. [17]. Here, the phonon dispersion, phonon velocity in cubic Si and cubic Ge have been calculated by using the Adiabatic Bond Charge Model (ABCM) [18] and employed as an input for the Monte Carlo simulations. The first Brillouin zone was discretized into a sampling of 29791 (31×31×31) wave vectors.

For Silicon, the phonon scattering rates have been derived from DFT as presented in Ref. [19]. For Germanium which is a heavier atom, DFT approach is more complex.

Consequently a set of Holland parameters specific to a FB approach [20] has been adjusted to fit the bulk thermal conductivity in the full temperature range. The considered relaxation times are those of the three-phonon scattering and phononimpurity scattering processes.

According to the definition used in Ref. [21], three out-of-equilibrium temperatures T, T + and T -are defined. The pseudo temperature T is the temperature of an equilibrium distribution of phonons having the same energy density considering the total population of phonon. T + and T -are called "hemispherical temperatures" and are defined by considering only the phonons with positive or only those with negative velocity, respectively.

Semi-transparent interface modeling

When a particle collides with a semi-transparent interface, both elastic and diffusive mechanisms are considered. Besides, the probability of transmission of a phonon colliding the interface is obtained within the Full-Band Diffusive Mismatch Model (DMM) [22]. In this model, the interface is assumed to have many defects so that each phonon colliding the interface loses totally the memory of its initial state. The wave vector and the mode of the phonon in the final state are independent of those in the initial state. The only conserved quantity is the phonon energy. Thus, the transmission coefficient depends only on the phonon energy, and neither on the wave vector, nor on the state of the phonon.

As the reflection and transmission processes are assumed to be diffusive and elastic, the angular distribution of the wave vector of reflected or transmitted phonon is chosen randomly over the states belonging to the same isoenergy surface to satisfy the energy conservation. Moreover, the appearance of a net flux along the interface must be prohibited and a Full-Band version of the Lambert's emission law [START_REF] Smith | The Design of Optical Systems: General[END_REF] must be applied to avoid a non-physical accumulation of phonons near the interface. Thus, among the states belonging to the same isoenergy surface, the final state j and thus the final wave vector 𝑞 ⃗ 𝑗 (and mode) are selected according to this following statistic weight:

𝑃 j ∝ 𝑣 𝑗 ⃗⃗ ⃗⃗. 𝑛 ⃗⃗ (1)
where 𝑣 𝑗 ⃗⃗ ⃗⃗ is the group velocity of the state j and 𝑛 ⃗⃗ is the vector normal to the interface. Fig. 1 shows the simulated Si/Ge heterostructure of length L. One device end is in contact with a hot thermostat at temperature TH and the other end is in contact with a cold thermostat at temperature TC. The Si/Ge junction is located at the middle of the device (x = L/2). The interface between Si and Ge is assumed to be abrupt. The width w and height h of the heterostructure (along y and z directions, respectively) are equal to 100 nm. As the cross-plane transport configuration is simulated, the structures are supposed to be infinite along the y and z directions. In that purpose, external interfaces that are normal to y and z directions follow periodic boundary conditions. (TH+TC)/2=300K.

Simulated Devices

The system is primarily uniformly meshed along the x-axis by using 20 cells. Then, this mesh could be refined at 5 nm around the heterojunction and near the thermostats to keep the cell length equal to or lower than 1 nm. It should be noticed that the size of a cell should be significantly higher than the atomic radii of Ge (125 pm) and Si (110 pm) [START_REF] Slater | [END_REF] i.e. it should remain in the order of magnitude of nanometer.

Results

Fig. 2 shows the simulated heat flux density Q as a function of the temperature difference between the two thermostats for a 20 nm long Si/Ge heterostructure. For all studied temperature biases, the cold and hot thermostats have been swapped. The evolution of the flux is clearly a linear function of the temperature difference. The linear coefficient gives a total thermal conductance of the heterostructure (for both thermostat configurations) of Gtot = 183 MW/m 2 /K.

Fig. 3. Heat flux density in a Si/Ge heterostructure as a function of the device length L. TH=302 K and TC=298 K.

However for TH -TC = 0 K, a non-zero flux of δQ = 45.3 MW/m 2 is extracted from the simulation. This value is assumed to be a good estimate of the MC simulation resolution. The statistical error for the calculation of a physical quantity f by using the Monte Carlo method is [START_REF] Volz | Microscale and Nanoscale Heat Transfer[END_REF]:

|𝑓 -〈𝑓〉| < 𝛼√ 〈𝑓 2 〉-〈𝑓〉 2 𝑁 , ( 2
)
where N is the number of random observations of the quantity f. Assuming that f follows the normal distribution, the interval of 95% confidence corresponds to α = 1.96. By increasing N, the statistical error can be reduced as needed. Nevertheless, the resolution of our MC simulation is the lower limit of the error related to the heat flux density estimation and will be used in the following to define the error bars. In Fig. 3, the heat flux density Q is plotted as a function of the device length L for a temperature bias TH -TC of 4.0 K. In ultra-short devices (L<100 nm), the fluxes tend to saturate and to become independent of the length as the phonon transport becomes ballistic. In long devices (L > 100 nm) in which the transport is diffusive, the heat flux density tends to decrease linearly when L increases. Moreover, the thermostat swap induces a shift of the flux smaller than the error bar of the MC simulation, thus no significant thermal rectification can be observed.

Fig. 4. Temperature profiles in a Si/Ge heterostructure of length L=20 nm in contact with TH = 302 K (298 K), TC = 298 K (302 K) in continuous lines (in dashed lines, respectively).

Fig. 5. Thermal interface conductance in a Si/Ge heterostructure as a function of the device length L at 300 K and with a temperature bias of TH-TC=4 K

Fig. 4 presents the temperature profiles: T, T + and T - in a 20 nm long device for both positive and negative temperature biases of 4 K. In the first configuration, the Ge bar is in contact with the hot thermostat and the Si bar is in contact with the cold thermostat (solid line) and, as the flux is positive along the x axis, T + is the highest temperature. By swapping the heterostructure bias (dashed line), the thermal flux is reversed and the phonons having positive x-velocity, i.e. propagating from the cold to the hot thermostat, are minority. Thus, in this configuration T + is the lowest temperature.

The temperature drop at the interface is sharp while the temperature gradient in the Si and Ge bars are weak. Besides, all the temperatures gradients (of T, T + and T -) are higher in Ge than in Si as Ge exhibits lower thermal conductivity. Moreover, the of T is slightly higher than those of T + and T -. It should be noted that, the net thermal flux is governed by the gradient of T within the Fourier's formalism or the difference between T + and T -at a given position within the formalism of Ref. [21].

The interface thermal conductance is finally computed from dividing the heat flux density Q by the difference of the hemispherical temperatures T + and T - on each side of the interface [21]:

𝐺 𝑖𝑛𝑡 = 𝑄 𝑇 + (𝐿 2 ⁄ -𝛿𝑥)-𝑇 -(𝐿 2+𝛿𝑥 ⁄ ) (3)
where δx stands for an infinitely small distance.

From the previous results of the resolution of the heat flux density and the total thermal conductance, the resolution of the temperature bias is estimated by using the following ratio:

𝛿(∆𝑇) = 𝛿𝑄 𝐺 𝑡𝑜𝑡 = 0.25 𝐾 (2)
For the device of length L, the lower bound for the statistical error of the interface thermal conductance 𝛿𝐺 𝑖𝑛𝑡 (𝐿) is deduced from the resolutions of the heat flux density and the temperature bias according to:

𝛿𝐺 𝑖𝑛𝑡 (𝐿) 𝐺 𝑖𝑛𝑡 (𝐿) = 𝛿𝑄(𝐿) 𝑄(𝐿) + 𝛿(∆𝑇) 𝑇 𝐻 -𝑇 𝐶 (3)
The resulting Gint are plotted in Fig. 5 for several device lengths. The width of error bar increases when the length L increases as the temperature difference (T H -T C ) near the interface also decreases leading to higher uncertainty in terms of extracted temperature.

However, the variation of Gint as a function of the device length remains lower than the error bars. Thus, Gint appears to be equal to 243.3±41.9 MW/m 2 independently of the device length. This is consistent with the DMM theory which predicts a unique value for the Si/Ge interface conductance of 225 MW/m 2 /K [22]. The modification of the transport regime of phonon in the studied length range does not affect Gint.

It should be noticed that the standard definition of the interface thermal conductance G int 𝑇 = Q/[𝑇(𝐿 2 ⁄ -𝛿𝑥) -𝑇(L 2 + δx ⁄ )] using only the pseudo temperature T leads to different results. Indeed, the values of G int 𝑇 that are 360 MW/m 2 for L=1nm and 350 MW/m 2 for L = 100 nm. are significantly higher than the DMM results.

Ensuring by assumptions spectral energy balance, the DMM does not provide by itself any physical ingredient that could induced thermal rectification. However, the presence of out of equilibrium phonons could have generated such asymmetry in the phonon transport. In Fig. 5, the difference between the values of Gint when heat flux is reversed is still smaller than the error bars and no thermal rectification due to "hot" phonons is observed.

Conclusion

In this work, the heat transfer through Si/Ge interfaces was investigated by using Full Band Monte Carlo simulation for phonon using DMM for modeling the phonon transmission at semi-transparent interfaces. Our MC simulator is able to calculate all the relevant temperature profiles, the heat flux density and finally the interface thermal conductance at heterojunctions. It can simulate short systems working in ballistic regime as well as long devices working in diffuse transport regime.

A linear relationship was obtained between the heat flux density and the applied temperature difference across a short heterostructure of 20 nm. Thus, even in the case of a strongly ballistic phonon transport, the concept of thermal conductance remains relevant. Moreover, the value of the interface thermal conductance obtained from this MC approach are consistent with the semi-analytical DMM model by using the concept of hemispherical temperature T + and T -. Besides, no thermal rectification is observed in our model at the Si/Ge interface.

Finally, this work demonstrates the relevance of this MC simulator to investigate the heat transport at the nanoscale as it provides very detailed insight into thermal transfer at interfaces.

Fig. 1 .

 1 Fig. 1. Si/Ge heterojunction. Orange and violet cells are made of Si and Ge, respectively. Hot and cold thermostats are red and blue, respectively. DMM interface is in Yellow. White faces stand for periodic boundaries.

Fig. 2 .

 2 Fig. 2. Heat flux density in a Si/Ge heterostructure of length L=20 nm as a function of the applied temperature bias TH-TC.