A curve-based material recognition method in MeV dual-energy X-ray imaging system

A curve-based material recognition method in MeV dual-energy X-ray imaging system 1 2 first-author 1 2 1 2 corresp Corresponding author, lliang@mail.tsinghua.edu.cn Corresponding author, lliang@mail.tsinghua.edu.cn lliang@mail.tsinghua.edu.cn Department of Engineering Physics, Tsinghua University, Beijing 100084, China Department of Engineering Physics, Tsinghua University, Beijing 100084, China Key Laboratory of Particle & Radiation Imaging (Tsinghua University), Ministry of Education, Beijing 100084, China Key Laboratory of Particle & Radiation Imaging (Tsinghua University), Ministry of Education, Beijing 100084, China High energy dual-energy X-ray Digital Radiography (DR) imaging is mainly used in the material recognition of cargo inspection. We introduce the development history and principle of the technology and describe the data process flow of our system. The system corrects original data to get a dual-energy transparence image. Material categories of all points in the image are identified by the classification curve, which is related to the X-ray energy spectrum. For the calibration of classification curve, our strategy involves a basic curve calibration and a real-time correction devoted to enhancing the classification accuracy. Image segmentation and denoising methods are applied to smooth the image. The image contains more information after colorization. Some results show that our methods achieve the desired effect. High energy Dual-energy Digital radiography Material recognition 1 Introduction 1 The X-ray imaging technique has become one of the most important tools in customs inspection. Presently, there are mainly two X-ray imaging modalities: radiography and computed tomography (CT). Although CT can provide 3-D structures and an accurate attenuation map of the cargo, its complexity and high price limit its application [ 1 1 - 3 3 1-3 ]. X-ray radiography, including single-energy and dual-energy, is still the mainstream technology. The development of X-ray radiography undergoes three stages: X-ray film photography, Computed Radiography (CR) and Digital Radiography (DR). The single-energy X-ray DR image merely gives the cumulative density information of the irradiated objects in one direction. It is used in preliminary medical diagnosis and simple security inspection. Since single energy X-ray DR provides limited information, the dual-energy method was developed. Low energy dual-energy X-ray DR imaging has been widely used in current security inspection equipment, which can detect and distinguish contraband by determining material atomic number Z . The X-ray’s energy here is usually lower than 1 MeV. This technology is inapplicable in high Z material recognition or cargo inspection, as the energy of the X-ray which can penetrate the object in these situations needs to be a few MeV. The British company Cambridge Imaging first proposed the idea of high energy dual-energy X-ray imaging. There were some disputes about the validity of the high energy dual-energy method in material recognition. The Russian Efremov research institute proved the feasibility of this method with their experimental prototype [ 4 4 ]. The Germany company Heimann and the American company EG & G applied X-ray hardening technology to this field and proposed the filter method. The department of engineering physics at Tsinghua University and its cooperative enterprise, Nuctech, established a platform and made some achievements on material recognition and related studies. The theory of high energy dual-energy X-ray DR imaging and material recognition has been deeply studied and the corresponding experiment results further validated the feasibility of dual-energy imaging material recognition [ 5 5 ]. In this paper, we construct an imaging system model and a whole data processing flow. For the best visual effects of the final results, we used the image smoothing strategy and image colorization processing. Some realization details are also given. The R -curve material recognition method is a typical high energy dual-energy X-ray DR material recognition method [ 6 6 ]. We developed a real-time R -curve calibration method. It deals with the differences of the R -curves of different energy spectra caused by the system status fluctuation and inconsistency. In Sec. 2, we introduced the technology principle and elaborated on the methods of a MeV dual-energy imaging model. We focused on the calibration strategy we designed. In Sec. 3, we gave and discussed some experimental results. We concluded and envisioned future work in Sec. 4. 2 THEORY AND METHOD 1 2.1 The principle of MeV dual-energy X-ray imaging in material recognition 2 The three main interactions between a photon and matter are the photoelectric effect, Compton scatter effect, and the electron pair effect. They respectively dominate the low ( < 1 MeV), middle (1 MeV∼3 MeV) and high ( > 3 MeV) energy range [ 7 7 ]. The corresponding attenuation coefficients, μ , have different dependences with a material atomic number, Z . We can give <math display="block" id="M1"><mrow><mtable columnalign="left"><mtr columnalign="left"><mtd columnalign="left"><mrow/></mtd><mtd columnalign="left"><mrow><msub><mi>μ</mi><mtext>P</mtext></msub><mo>∝</mo><msup><mi>Z</mi><mn>4</mn></msup><mtext> </mtext><mo><</mo><mn>1,</mn></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow/></mtd><mtd columnalign="left"><mrow><msub><mi>μ</mi><mrow><mtext>CS</mtext></mrow></msub><mo>∝</mo><mn>1</mn><mtext> </mtext><mn>1</mn><mo>∼</mo><mn>3,</mn></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow/></mtd><mtd columnalign="left"><mrow><msub><mi>μ</mi><mrow><mtext>EP</mtext></mrow></msub><mo>∝</mo><mi>Z</mi><mtext> </mtext><mo>></mo><mn>3,</mn></mrow></mtd></mtr></mtable></mrow></math> where P, CS, and EP are the abbreviations of the three effects. Consider an X-ray source whose energy spectrum is N ( E ) and the highest energy is E m , a single substance with an atomic number of Z , an attenuation coefficient function of μ ( E , Z ), and a thickness t , the transparence, T is <math display="block" id="M2"><mrow><mi>T</mi><mo>=</mo><mfrac><mi>I</mi><mrow><msub><mi>I</mi><mn>0</mn></msub></mrow></mfrac><mo>=</mo><mfrac><mrow><mstyle displaystyle="true"><mrow><msubsup><mo>∫</mo><mn>0</mn><mrow><msub><mi>E</mi><mi>m</mi></msub></mrow></msubsup></mrow></mstyle><mi>N</mi><mo stretchy="false">(</mo><mi>E</mi><mo stretchy="false">)</mo><msup><mi>e</mi><mrow><mo>−</mo><mi>μ</mi><mo stretchy="false">(</mo><mi>E</mi><mn>,</mn><mi>Z</mi><mo stretchy="false">)</mo><mi>t</mi></mrow></msup><mtext>d</mtext><mi>E</mi></mrow><mrow><mstyle displaystyle="true"><mrow><msubsup><mo>∫</mo><mn>0</mn><mrow><msub><mi>E</mi><mi>m</mi></msub></mrow></msubsup></mrow></mstyle><mi>N</mi><mo stretchy="false">(</mo><mi>E</mi><mo stretchy="false">)</mo><mtext>d</mtext><mi>E</mi></mrow></mfrac><mn>.</mn></mrow></math> In a dual-energy situation, the boundary energy of the two X-ray sources are E 1 and E 2 . We defined the logarithmic ratio of T as <math display="block" id="M3"><mrow><mi>R</mi><mo>=</mo><mfrac><mrow><mi>ln</mi><msub><mi>T</mi><mn>1</mn></msub></mrow><mrow><mi>ln</mi><msub><mi>T</mi><mn>2</mn></msub></mrow></mfrac><mo>=</mo><mfrac><mrow><mover accent="true"><mi>μ</mi><mo>¯</mo></mover><mo stretchy="false">(</mo><msub><mi>E</mi><mn>1</mn></msub><mn>,</mn><mi>t</mi><mn>,</mn><mi>Z</mi><mo stretchy="false">)</mo></mrow><mrow><mover accent="true"><mi>μ</mi><mo>¯</mo></mover><mo stretchy="false">(</mo><msub><mi>E</mi><mn>2</mn></msub><mn>,</mn><mi>t</mi><mn>,</mn><mi>Z</mi><mo stretchy="false">)</mo></mrow></mfrac><mn>,</mn></mrow></math> where R is the ratio of the equivalent attenuation coefficient, <math id="M4"><mover accent="true"><mi>μ</mi><mo>¯</mo></mover></math> . When the X-ray source is monochromatic, which means the energy spectrum N ( E ) is a single line, Eq. (2) and Eq. (3) can be simplified. Suppose E 1 is in the low energy range and E 2 is in the middle energy range, R can be written as <math display="block" id="M5"><mrow><mi>R</mi><mo>=</mo><mfrac><mrow><mi>ln</mi><msub><mi>T</mi><mn>1</mn></msub></mrow><mrow><mi>ln</mi><msub><mi>T</mi><mn>2</mn></msub></mrow></mfrac><mo>=</mo><mfrac><mrow><mi>ln</mi><msup><mi>e</mi><mrow><mo>−</mo><mi>μ</mi><mo stretchy="false">(</mo><msub><mi>E</mi><mn>1</mn></msub><mn>,</mn><mi>t</mi><mn>,</mn><mi>Z</mi><mo stretchy="false">)</mo><mi>t</mi></mrow></msup></mrow><mrow><mi>ln</mi><msup><mi>e</mi><mrow><mo>−</mo><mi>μ</mi><mo stretchy="false">(</mo><msub><mi>E</mi><mn>2</mn></msub><mn>,</mn><mi>t</mi><mn>,</mn><mi>Z</mi><mo stretchy="false">)</mo><mi>t</mi></mrow></msup></mrow></mfrac><mo>=</mo><mfrac><mrow><mi>μ</mi><mo stretchy="false">(</mo><msub><mi>E</mi><mn>1</mn></msub><mn>,</mn><mi>Z</mi><mo stretchy="false">)</mo></mrow><mrow><mi>μ</mi><mo stretchy="false">(</mo><msub><mi>E</mi><mn>2</mn></msub><mn>,</mn><mi>Z</mi><mo stretchy="false">)</mo></mrow></mfrac><mo>≈</mo><mfrac><mrow><msub><mi>μ</mi><mtext>P</mtext></msub><mo stretchy="false">(</mo><msub><mi>E</mi><mn>1</mn></msub><mn>,</mn><mi>Z</mi><mo stretchy="false">)</mo></mrow><mrow><msub><mi>μ</mi><mrow><mtext>CS</mtext></mrow></msub><msub><mi>E</mi><mn>2</mn></msub><mn>,</mn><mi>Z</mi></mrow></mfrac><mo>≈</mo><mi>c</mi><mo stretchy="false">(</mo><msub><mi>E</mi><mn>1</mn></msub><mn>,</mn><msub><mi>E</mi><mn>2</mn></msub><mo stretchy="false">)</mo><msup><mi>Z</mi><mn>4</mn></msup><mn>.</mn></mrow></math> From Eq. (3), R is easily computable. From Eq. (4), R is a clear indication of Z . Besides, when the X-ray is polychromatic, a great dependence between R and Z still exists. Based on these facts, low energy dual-energy X-ray DR imaging technology has been widely used in small security inspection devices for material discrimination. Low energy dual-energy means that E 1 and E 2 are usually lower than 1 MeV. When Z is high and the irradiated object is thick, low energy X-ray imaging becomes useless. If we change E 1 to the high energy range and keep E 2 in the middle energy range, then Eq. (4) is <math display="block" id="M6"><mrow><mi>R</mi><mo>=</mo><mfrac><mrow><mi>μ</mi><mo stretchy="false">(</mo><msub><mi>E</mi><mn>1</mn></msub><mn>,</mn><mi>Z</mi><mo stretchy="false">)</mo></mrow><mrow><mi>μ</mi><mo stretchy="false">(</mo><msub><mi>E</mi><mn>2</mn></msub><mn>,</mn><mi>Z</mi><mo stretchy="false">)</mo></mrow></mfrac><mo>≈</mo><mfrac><mrow><msub><mi>μ</mi><mrow><mtext>EP</mtext></mrow></msub><mo stretchy="false">(</mo><msub><mi>E</mi><mn>1</mn></msub><mn>,</mn><mi>Z</mi><mo stretchy="false">)</mo></mrow><mrow><msub><mi>μ</mi><mrow><mtext>CS</mtext></mrow></msub><mo stretchy="false">(</mo><msub><mi>E</mi><mn>2</mn></msub><mn>,</mn><mi>Z</mi><mo stretchy="false">)</mo></mrow></mfrac><mo>≈</mo><mi>c</mi><mo stretchy="false">(</mo><msub><mi>E</mi><mn>1</mn></msub><mn>,</mn><msub><mi>E</mi><mn>2</mn></msub><mo stretchy="false">)</mo><mi>Z</mi><mn>,</mn></mrow></math> and R is dependent on Z to a certain extent, so it can be still used to classify material. This ideal conclusion without consideration of the subordinated interaction of photons with matter is based on the assumption of the single-line X-ray energy spectrum. In fact, a MeV X-ray DR system uses the linear accelerator as the X-ray source, which generates X-rays with a broad energy spectrum [ 8 8 ]. Most of the photons with distributed in the middle energy range, where μ and Z have no correlations. The effectiveness of R in material recognition is not obvious. It was found that, although R changes when thickness t changes, R is still dependent on Z [ 9 9 ]. Here E 1 and E 2 are usually higher than 3 MeV. 2.2 A MeV dual-energy system 2 We use a schematic model (See Fig. 1 Fig. 1 ), including an accelerator, which emits a vertical fan shaped X-ray beam, a scanning track, which is perpendicular to the X-ray’s main beam direction, an L-shape detector, and a data processing unit. The cargo moves along the scanning track while the L-shape detector receives the photons passing through the cargo and form the dual-energy X-ray images. The data processing unit consists of four steps. First, correct the acquired original dual-energy X-ray images and calculate the dual-energy transparence images. Second, use the classification curve on the dual-energy transparence images to form the material information image. Third, to improve the image quality, a smoothing process is implemented. The final step is the colorization of the gray image. In the next four sections, we will introduce these four steps of the data processing unit and mainly concentrate on the calibration of the classification curve. We propose a real-time R -curve calibration method. In Sec. 3, we can see that our method enhances the classification accuracy and gives a satisfactory visual result. 2.3 Data acquisition and pre-processing 2 We assume that the accelerator produces the dual-energy X-ray simultaneously. Accordingly, the detector is able to distinguish high and low energy X-rays and form the dual-energy X-ray images separately. Other aspects of our model are basically the same as reality. The fan shaped X-ray has an angular distribution. Its main beam is located near the middle of the fan and has a maximum intensity decreasing towards both sides. It causes different vertical positions of the X-ray image which have a different X-ray intensity and energy spectrum. The accelerator state fluctuation in the scanning process causes different lateral position of the X-ray image which have a different X-ray intensity and energy spectrum. The detector background and response inconsistencies also exist. In data processing unit, the first step is to correct the original data and get the dual-energy transparence image. We give Eq. (6) based on Eq. (2). <math display="block" id="M7"><mrow><mi>T</mi><mo stretchy="false">(</mo><mi>x</mi><mn>,</mn><mi>y</mi><mo stretchy="false">)</mo><mo>=</mo><mfrac><mrow><mi>I</mi><mo stretchy="false">(</mo><mi>x</mi><mn>,</mn><mi>y</mi><mo stretchy="false">)</mo><mo>−</mo><msub><mi>I</mi><mrow><mtext>BK</mtext></mrow></msub><mo stretchy="false">(</mo><mi>x</mi><mn>,</mn><mi>y</mi><mo stretchy="false">)</mo></mrow><mrow><msub><mi>I</mi><mn>0</mn></msub><mo stretchy="false">(</mo><mi>x</mi><mn>,</mn><mi>y</mi><mo stretchy="false">)</mo></mrow></mfrac><mo>⋅</mo><mi>L</mi><mi>D</mi><mo stretchy="false">(</mo><mi>x</mi><mo stretchy="false">)</mo><mn>.</mn></mrow></math> The coordinate x , y represents the lateral and vertical position. I is the original dual-energy image. I 0 is the air image. I BK is the detector background image. The correction factor, LD ( x ), obtained by monitoring the accelerator state fluctuation in the scanning process is a function of the lateral position, x . In Eq. (6), the division by I 0 corrects the intensity angular distribution and detector response inconsistencies. Rest corrections are also done by Eq. (6). After the point-by-point correction and simple denosing, we can get the dual-energy transparence image, T . To take advantages of two transparence images, we use a point-wise weighted fusion of them. In the thick places of the irradiated object, the high energy transparence image gives more information than the low energy transparence image. The conclusion is opposite in the thin places. The gray value, representing one point’s thickness, determines the weight. Using the dual-energy transparence image and the classification curve, which we will elaborate on next, we can get the material information image. The colorization assembling the fused transparence image and the material information image gives a final result. 2.4 Curve-based material recognition method 2 Four kinds of classification curves, R -curve ( Fig 2(a) Fig 2(a) ), T 2 - T 1 curve ( Fig. 2(b) Fig. 2(b) ), α-curve ( Fig. 2(c) Fig. 2(c) ) [ 10 10 , 11 11 ], and a separability curve which denotes the transparence difference of the two materials ( Fig. 2(d) Fig. 2(d) ), are shown in Fig. 2 Fig. 2 . Among them, the R -curve has the best visual separability, so we use the R -curve to show the calibration output and result comparison. The data in Fig. 2 Fig. 2 are ideal. The longitudinal coordinate of the R -curve is R , defined in Eq. (3), and the horizontal coordinate is the inverse of T 1 or T 2 . Different materials have different R -curves. They are arranged in Z ’s increasing order when they are in the same image [ 6 6 ]. The classification curve here has several single R -curves related to the same number of typical objects. We assume that there are four typical objects, including Pb, Fe, Al, and CH 2 (or C), as we can see in Fig. 2 Fig. 2 . ( x , y ) is one point of the dual-energy transparence image. There are two transparence values, T 1 ( x , y ) and T 2 ( x , y ). They give a point, ( R , 1/ T 2 ), on the classification curve, C . Because R -curves of different materials are arranged in Z ’s increasing order, the point ( x , y ) is obtained by interpolating between two adjacent R -curves in C . By repeating this procedure on each point of the dual-energy transparence image, the material information image Z ( x , y ) is formed. Equations (2) and (3) show that the R -curve correlates to the energy spectrum. We already know that each point of the dual-energy transparence image has different X-ray intensity and energy spectrum because of the angular distribution and the accelerator state fluctuation. This fact causes a problem that different points of the dual-energy transparence image need different classification curves. However, it is impossible to calibrate all of classification curves. We designed a new calibration method to obtain all of the points’ approximate R -curves. The calibration strategy we employ takes two steps. First, we get the basic classification curve before the system scans cargo or after the system state significantly changes. This step requires scanning several typical materials. In one scanning process, the system only scans one typical material with a different mass thickness. Like the common scan, we get two transparence images T 1 , T 2 , and each point, ( x , y ), belongs to an R -curve, R xy . Assuming a steady system state in the process, we ignore the differences of energy spectrum in different lateral positions and just use several vertical points to represent the whole angular distribution. So we have <math display="block" id="M8"><mrow><msub><mi>R</mi><mrow><mi>x</mi><mi>y</mi></mrow></msub><mo>=</mo><msub><mi>R</mi><mi>y</mi></msub><mn>,</mn><mtext> </mtext><mi>y</mi><mo>=</mo><msub><mi>y</mi><mn>1</mn></msub><mn>,</mn><msub><mi>y</mi><mn>2</mn></msub><mn>,</mn><mo>…</mo><mn>,</mn><msub><mi>y</mi><mi>n</mi></msub><mn>.</mn></mrow></math> The data is <math display="block" id="M9"><mtable columnalign="left"><mtr><mtd><msub><mtext>data</mtext><mi>y</mi></msub><mo>=</mo><mo>{</mo><msub><mi>R</mi><mi>y</mi></msub><mo stretchy="false">[</mo><msub><mi>T</mi><mn>1</mn></msub><mo stretchy="false">(</mo><msub><mi>x</mi><mi>i</mi></msub><mn>,</mn><mi>y</mi><mo stretchy="false">)</mo><mn>,</mn><msub><mi>T</mi><mn>2</mn></msub><mo stretchy="false">(</mo><msub><mi>x</mi><mi>i</mi></msub><mn>,</mn><mi>y</mi><mo stretchy="false">)</mo><mo stretchy="false">]</mo><mn>,</mn><mi>i</mi><mo>=</mo><mn>1,2,</mn><mo>…</mo><mn>,</mn><mi>m</mi><mo>}</mo><mn>,</mn></mtd></mtr><mtr><mtd><mi>y</mi><mo>=</mo><msub><mi>y</mi><mn>1</mn></msub><mn>,</mn><msub><mi>y</mi><mn>2</mn></msub><mn>,</mn><mo>…</mo><mn>,</mn><msub><mi>y</mi><mi>n</mi></msub><mn>,</mn></mtd></mtr></mtable></math> where m is the number of different mass thicknesses and n is the number of vertical points. The R -curve is fitted with this data and the classification curve is formed with these fitted curves of several typical objects. Here, we use Pb, Fe, Al, and C. <math display="block" id="M10"><mrow><mtable columnalign="left"><mtr columnalign="left"><mtd columnalign="left"><mrow/></mtd><mtd columnalign="left"><mrow><mo>{</mo><msubsup><mi>R</mi><mi>y</mi><mi>Z</mi></msubsup><mo>|</mo><msubsup><mi>R</mi><mi>y</mi><mi>Z</mi></msubsup><mo>=</mo><mtext>fit</mtext><mo stretchy="false">(</mo><msubsup><mrow><mtext>data</mtext></mrow><mi>y</mi><mi>Z</mi></msubsup><mo stretchy="false">)</mo><mn>,</mn><mi>y</mi><mo>=</mo><msub><mi>y</mi><mn>1</mn></msub><mn>,</mn><msub><mi>y</mi><mn>2</mn></msub><mn>,</mn><mo>…</mo><mn>,</mn><msub><mi>y</mi><mi>n</mi></msub><mo>}</mo><mn>,</mn></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow/></mtd><mtd columnalign="left"><mrow><msub><mi>C</mi><mrow><msub><mi>y</mi><mi>i</mi></msub></mrow></msub><mo>=</mo><mo>{</mo><msub><mi>R</mi><mrow><msub><mi>y</mi><mi>i</mi></msub><mn>,</mn><mi>Z</mi></mrow></msub><mo>|</mo><mi>Z</mi><mo>=</mo><mi>P</mi><mi>b</mi><mn>,</mn><mi>F</mi><mi>e</mi><mn>,</mn><mi>A</mi><mi>l</mi><mn>,</mn><mi>C</mi><mo>}</mo><mn>,</mn><mtext> </mtext><mi>i</mi><mo>=</mo><mn>1,2,</mn><mo>…</mo><mn>,</mn><mi>n</mi><mn>.</mn></mrow></mtd></mtr></mtable></mrow></math> In the cargo scan, the system monitors the state variation to complete the real-time calibration of the classification curve. A small device consisted of the typical materials with the single thickness set in a certain vertical position, Y , and scanned synchronously with the cargo. The data is <math display="block" id="M11"><mrow><mtable columnalign="left"><mtr columnalign="left"><mtd columnalign="left"><mrow/></mtd><mtd columnalign="left"><mrow><msubsup><mi>T</mi><mrow><mn>1</mn><mi>r</mi></mrow><mi>Z</mi></msubsup><mo stretchy="false">(</mo><msub><mi>x</mi><mi>i</mi></msub><mn>,</mn><mi>Y</mi><mo stretchy="false">)</mo><mn>,</mn><msubsup><mi>T</mi><mrow><mn>2</mn><mi>r</mi></mrow><mi>Z</mi></msubsup><mo stretchy="false">(</mo><msub><mi>x</mi><mi>i</mi></msub><mn>,</mn><mi>Y</mi><mo stretchy="false">)</mo><mn>,</mn><mi>i</mi><mo>=</mo><mn>1,2,</mn><mo>…</mo><mn>,</mn><mi>n</mi><mi>x</mi><mn>,</mn></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow/></mtd><mtd columnalign="left"><mrow><mi>Z</mi><mo>=</mo><mi>P</mi><mi>b</mi><mn>,</mn><mi>F</mi><mi>e</mi><mn>,</mn><mi>A</mi><mi>l</mi><mn>,</mn><mi>C</mi><mn>,</mn></mrow></mtd></mtr></mtable></mrow></math> where nx is the number of horizontal pixels. In the first step, this data is also saved as <math id="M12"><mrow><msubsup><mi>T</mi><mrow><mn>1</mn><mi>b</mi></mrow><mi>Z</mi></msubsup><mo stretchy="false">(</mo><msub><mi>x</mi><mi>i</mi></msub><mn>,</mn><mi>Y</mi><mo stretchy="false">)</mo><mn>,</mn><msubsup><mi>T</mi><mrow><mn>2</mn><mi>b</mi></mrow><mi>Z</mi></msubsup><mo stretchy="false">(</mo><msub><mi>x</mi><mi>i</mi></msub><mn>,</mn><mi>Y</mi><mo stretchy="false">)</mo></mrow></math> . They are the average of T on the horizontal direction, as we ignore the lateral difference. Then the revised classification curve will be <math display="block" id="M13"><mrow><msub><mi>C</mi><mrow><msub><mi>x</mi><mi>j</mi></msub><msub><mi>y</mi><mi>i</mi></msub></mrow></msub><mo>=</mo><mo>{</mo><msubsup><mi>R</mi><mrow><msub><mi>x</mi><mi>j</mi></msub><msub><mi>y</mi><mi>j</mi></msub></mrow><mi>Z</mi></msubsup><mo>|</mo><msubsup><mi>R</mi><mrow><msub><mi>x</mi><mi>j</mi></msub><msub><mi>y</mi><mi>j</mi></msub></mrow><mi>Z</mi></msubsup><mo>=</mo><msubsup><mi>R</mi><mrow><msub><mi>y</mi><mi>i</mi></msub></mrow><mi>Z</mi></msubsup><mo>∗</mo><mi>F</mi><mo stretchy="false">(</mo><msubsup><mi>T</mi><mrow><mn>1</mn><mi>r</mi></mrow><mi>Z</mi></msubsup><mn>,</mn><msubsup><mi>T</mi><mrow><mn>2</mn><mi>r</mi></mrow><mi>Z</mi></msubsup><mn>,</mn><msubsup><mi>T</mi><mrow><mn>1</mn><mi>b</mi></mrow><mi>Z</mi></msubsup><mn>,</mn><msubsup><mi>T</mi><mrow><mn>2</mn><mi>b</mi></mrow><mi>Z</mi></msubsup><mo>;</mo><msub><mi>x</mi><mi>j</mi></msub><mn>,</mn><mi>Y</mi><mo stretchy="false">)</mo><mn>,</mn><mi>Z</mi><mo>=</mo><mi>P</mi><mi>b</mi><mn>,</mn><mi>F</mi><mi>e</mi><mn>,</mn><mi>A</mi><mi>l</mi><mn>,</mn><mi>C</mi><mo>}</mo><mn>,</mn><mi>i</mi><mo>=</mo><mn>1,2,</mn><mo>…</mo><mn>,</mn><mi>n</mi><mn>,</mn><mi>j</mi><mo>=</mo><mn>1,2,</mn><mo>…</mo><mn>,</mn><mi>n</mi><mi>x</mi><mn>.</mn></mrow></math> The function F uses real-time calibrating data to estimate the difference between <math id="M14"><mrow><msubsup><mi>R</mi><mrow><mi>x</mi><mi>v</mi></mrow><mi>Z</mi></msubsup></mrow></math> and <math id="M15"><mrow><msubsup><mi>R</mi><mi>v</mi><mi>Z</mi></msubsup></mrow></math> . The correction method in Eq. (11) relies on how function F is calculated. Considering the statistical straggling of T ( x j , Y ), we choose a segment of x j and use the average of the monitoring data to correct the classification curve. 2.5 Smooth of material information image 2 The detector signal noise inevitably exists and it is even amplified in the data processing. It has an impact on material recognition accuracy and makes the material information image rough. The visual effect of final result after colorization is not good enough. Material information image quality promotion is necessary. Some literature proposes image segmentation smooth strategy [ 9 9 , 12 12 ]. We apply the same idea. The first step of our smooth process is the segmentation of the fused transparence image, which is obtained previously. The image is segmented into regions, which keep the continuity of the irradiated objects interior as much as possible and discriminate different irradiated objects as clearly as possible. The average of all the Z values in a corresponding region in the material information image is assigned to all points in this region. The general image segmentation algorithm, like the single-pass spilt-merge algorithm, or data clustering algorithm, like the Leader algorithm, can be used here with some adjustment [ 13 13 ]. The irradiated objects may be mixed and disorderly. So the segmentation or cluster result may have too many small areas with only several pixels. This over-segmentation can be solved by merging the small area into the nearest large area. The ’nearest’ means not only the distance but also the similarity between them. Using the average of the Z values in one segmentation region to replace all points in this region brings some loss of the original information. The denoising approaches can give a better material information image and the majority of the image remains the same. It is a challenge to find a better method, which can smooth the image while maximizing the retention of the original information. 2.6 Colorization 2 The idea of colorization and the IHS color space was proposed in Ref. [ 9 9 , 12 12 ]. We use a similar approach applying HLS color space [ 14 14 ]. In the colorization, different colors represent different materials. We use three color spaces, including RGB, HLS, and YUV. If all three values in a color space are known, a color is determined. We divide the range of Z into several parts as follow, <math display="block" id="M16"><mrow><mi>Z</mi><mo><</mo><msub><mi>Z</mi><mn>1</mn></msub><mo>;</mo><msub><mi>Z</mi><mn>1</mn></msub><mo>≤</mo><msub><mi>Z</mi><mn>2</mn></msub><mo>;</mo><mo>…</mo><mo>;</mo><msub><mi>Z</mi><mrow><mi>p</mi><mo>−</mo><mn>1</mn></mrow></msub><mo>≤</mo><msub><mi>Z</mi><mi>p</mi></msub><mo>;</mo><msub><mi>Z</mi><mi>p</mi></msub><mo>≤</mo><mi>Z</mi><mn>.</mn></mrow></math> There are p +1 hues H 1 , H 2 , … , H p +1 . When the Z value of one point in the material information image falls into the j th part, the H value in the HLS color space equals H j . The sensitivity of the human eye to different color is different. There are red, green, and blue with the same L value in the HLS color space. The brightness felt by eyes is different. Green is the brightest and red is brighter than blue. If the L value equals the gray value, the points with same gray value will give a different brightness when we look at the final result. It is inappropriate. In the YUV color space, colors with same Y value give the closest brightness feeling. Let the Y value equal the gray value of the fused transparence image. The YUV color space and HLS color space can be converted to each other. So we have <math display="block" id="M17"><mrow><mi>Y</mi><mo>=</mo><mi>f</mi><mo stretchy="false">(</mo><mi>H</mi><mn>,</mn><mi>L</mi><mn>,</mn><mi>S</mi><mo stretchy="false">)</mo><mn>.</mn></mrow></math> As Y , H , f is known and S is given, L is the solution of Eq. (13). Then all three values in the HLS color space are already set, and a color is determined. Repeat this procedure on every point of the material information image to get the final result. The table of hue and the saturation value, S , are changeable and directly influences of the image is visual effect. The mapping relationship between Y and the gray value can be optimized by some transformation, like logarithm stretching. 3 EXPERIMENTAL RESULT 1 The data is provided by a 6/3 MeV X-ray DR imaging cargo inspection system, which applies our basic design of a data processing unit. The accelerator alternatively emits the high and low energy X-rays. The emission frequencies are both 40 Hz. We use a single column of the CWO detector and a scanning speed is 0.2 m/s. The basic classification curve is formed by scanning three single materials, C, Al, and Fe. A stair-step object with a single material component is scanned to form one R -curve [ 15 15 ]. All R -curves together form a classification curve. For clarity of the results, we use one R -curve representing a classification curve. There are two different system states. In Fig. 3 Fig. 3 , the dotted curve is the classification curve in state 1 and the dashed curve represents the classification curve in state 2. The solid curve is the revised classification curve based on the classification curve in state 1. The revision uses the difference between the real-time calibrating data of the two states. We think that the dashed curve is the ’true’ classification curve in state 2 and the solid curve is an estimation of the ’true’ one. Their closeness shows the effectiveness and rationality of our calibration strategy. We arrange eight objects in the order of Pb, Fe, Al, C, Al, C, Pb, and Fe and divide them into two groups according to size. Their mass thicknesses are all 40 g/cm 2 . The scanning of these eight objects is in state 2. Suppose we do not know the classification curve in state 2 (dashed curve in Fig. 3 Fig. 3 ). Our calibration strategy means that we can get the classification curve in state 2 if we know the classification curve in state 1 and the real-time calibration data of the two states. In Fig. 4 Fig. 4 , the image a) on the left is the final result without the use of the real-time calibrating data. For the four larger objects Al, C, Pb, and Fe, the specific Z values obtained using the classification curve in state 1 (dotted curve in Fig. 3 Fig. 3 ) are 46, 27, 62, and 54. The image b) in the middle is the final result with the use of the real-time calibrating data. For the four larger objects Al, C, Pb, and Fe, the specific Z values obtained by using the calibrated curve (solid curve in Fig. 3 Fig. 3 ) are 18, 9, 54, and 45. The image c) on the right is the color table. The system colorization settings give that the hues of four typical objects, C, Al, Fe, and Pb, are orange, green, blue and purple. We also use the classification curve in state 2 (dashed curve in Fig. 3 Fig. 3 ) and get the Z values 16, 5, 54, and 45. Note that there is no R -curve of Pb, the Z values of the Pb object from three sets of results are 62, 54, and 54 and far form 82. The consequent color of the Pb object is deviated from the righteous color. Although the Z value of the Fe object has an obvious deviation due to the lack of the R -curve of Pb, it is in the righteous region and, thus, the color is also righteous. The data and figure comparison clearly show the effectiveness of the real-time calibration and also matches the comparison of the curves in Fig. 3 Fig. 3 . Besides Pb, the other materials’ results show that our calibration strategy enhances classification accuracy. Notice that all objects in Fig. 4 Fig. 4 are made up of the uniform single material and the two classification curves under two states have distinct differences. The conditions are ideal and, accordingly, the results are good. In the real scan, the scanned objects are always complex. We may be unable to figure out true Z values of all the points, thus, we cannot verify the accuracy of the recognition results. What we can do is to make the classification curve, whether directly calibrated or real-time revised, as accurate as possible. There are two points to note. First, the calibration of the basic classification curve, which is used as the basis for a steady system state. Second, the real-time monitoring data has statistical fluctuation even exceeding the state fluctuation of the accelerator or angular distribution. Using the average of a piece of data to reflect the variation is better than using every single point. In Fig. 5 Fig. 5 , the comparison is the final color image with and without the smooth process. The larger orange object in image a) looks not uniform, although it should be the same color. We can see that the smoothing improves the image quality. However, the effect is obsolete because of the monotony and uniformity of the irradiated objects. When the objects are complex, the smoothing process will influence the classification accuracy because the rearrangement of the Z values may change the well classified points, which are segmented into a region with a fairly different Z value. In Fig. 6 Fig. 6 , we give a cargo inspection result. The emission frequencies of the dual-energy X-rays are both 33 Hz. The scanning speed is 0.2 m/s. The irradiated objects from left to right are a cigarette, salt, sugar, coffee, buckets of water, and concrete. According to the continuity of the irradiated objects, we can assure that the spots and stripes in the object regions of the image on the top are noise and need to be removed. The smoothing process eliminates the noise spots and nonuniformity and significantly promotes the image quality. In the red circled region of the image on the top, the bottom margin of the bucket is overwhelmed by the noise and almost disappears. After the smoothing process, the margin is recovered. Our smoothing method may strengthen the details of simple object regions. In the cigarette region, the thickness of the cigarette is small. When the irradiated object is thin, the calibration curve separability is worse. With the data fluctuation, the final color image will be full of stripes and spots, as we can see in the amplified region of the image on the top. The cigarette belongs to the orange category, so the blue pixels are noise. The smoothing process take the average of in the cigarette region of the material information image. So the final color will be the middle color between orange and blue, according to the color table. The color changes to green, as we can see in the amplified region of the image on the bottom. 4 DISCUSSION AND CONCLUSION 1 We described a simple MeV dual-energy X-ray DR imaging cargo inspection system with a detailed description of the data processing unit. The preliminary treatment converts the original data into images for subsequent processing needs. The calibration strategy of the classification curve enhances the classification performance. The smoothing of the material information image is to enhance the image quality. Segmentation is the key part of our smoothing process. Better segmentation methods lead to better image quality. Color imaging can carry more information and give a better visual effect. Colorization can be regulated in different application environments. This system design has a certain guiding significance to the engineering. Our calibration method is devoted to give the correct classification curve and enhance the classification accuracy. There are several different directions to achieve this goal or push forward the further development of dual-energy X-ray imaging technology. Add a low energy detector to promote recognition ability of thin objects [ 16 16 ]. Use the Cerenkov detector, which has a threshold and a good response to high energy X-rays [ 17 17 ]. Add the obstacle’s classification curve on the blocked material’s classification to enhance the blocked material’s recognition accuracy [ 18 18 ]. Use small angle scatter to realize the material recognition [ 19 19 ]. These methods can be applied to our system or a new imaging system based on our model and data processing flow. References L. Li , Z.Q. Chen , Y.X. Xing , et al . A general exact method for synthesizing parallel-beam projections from cone-beam projections via filtered backprojection . Phys Med Biol , 2006 , 51 : 5643 - 5654 . DOI: 10.1088/0031-9155/51/21/017 http://doi.org/10.1088/0031-9155/51/21/017 A general exact method for synthesizing parallel-beam projections from cone-beam projections via filtered backprojection G.W. Zhang , J.P. Cheng , L. Zhang , et al . A practical reconstruction method for dual energy computed tomography . J X-ray Sci Technol , 2008 , 16 : 67 - 88 . A practical reconstruction method for dual energy computed tomography L. Li , K.J. Kang , Z.Q. Chen , et al . A general region-of-interest image reconstruction approach with truncated Hilbert transform . J X-ray Sci Technol , 2009 , 17 : 135 - 152 . DOI: 10.3233/XST-2009-0218 http://doi.org/10.3233/XST-2009-0218 A general region-of-interest image reconstruction approach with truncated Hilbert transform M.F. Vorogushin , Experiments on material recognition for 8 MeV customs inspection system for trucks and large-scale container . Proceeding of the 20th Int’l Linear Accelerator Conference , Monterey, CA, USA , 2000 . Experiments on material recognition for 8 MeV customs inspection system for trucks and large-scale container X.W. Wang , J.M. Li , C.X. Tang , et al . Material discrimination by high-energy X-Ray dual-energy imaging . High Energy Physics and Nuclear Physics , 2007 , 31 : 1076 - 1081 . Material discrimination by high-energy X-Ray dual-energy imaging Z.Q. Chen , T. Zhao , L. Li , The review of high energy dual-energy X-ray DR imaging and material recognition technology . CT Theor Appl , 2014 , 23 : 731 - 742 . (in Chinese) The review of high energy dual-energy X-ray DR imaging and material recognition technology X.W. Wang , Research on the theory and experiment of material discrimination by dual-energy method in high-energy X-ray imaging. Ph.D. Thesis , Tsinghua University , 2005 . (in Chinese) Research on the theory and experiment of material discrimination by dual-energy method in high-energy X-ray imaging. Ph.D. Thesis Z.Q. Chen , L. Li , J.C. Feng , New development of high energy industrial computed tomography . CT Theor Appl , 2005 , 14 : 1 - 4 . (in Chinese) DOI: 10.3969/j.issn.1004-4140.2005.04.001 http://doi.org/10.3969/j.issn.1004-4140.2005.04.001 New development of high energy industrial computed tomography S. Ogorodnikov , V. Petrunin , Processing of interlaced images in 4 MeV-10 MeV dual energy customs system for material recognition . Phys Rev Spec Top-AC , 2002 , 5 : 104701 . DOI: 10.1103/PhysRevSTAB.5.104701 http://doi.org/10.1103/PhysRevSTAB.5.104701 Processing of interlaced images in 4 MeV-10 MeV dual energy customs system for material recognition S.A. Ogorodnikov , V.I. Petrunin , M.F. Vorogushin , Application of high-penetrating introscopy systems for recognition of materials . Proceeding of the 7th European Particle Accelerator Conference , Vienna, Austria , 2000 . Application of high-penetrating introscopy systems for recognition of materials V.L. Novikov , S.A. Ogorodnikov , V.I. Petrunin , Dual energy method of material recognition in high energy introscopy systems . Proceeding of the 16th Int’l. Workshop on Charged Particle Linear accelerators , Alushta, Crimea, Ukraine , 1999 . (in Chinese) Dual energy method of material recognition in high energy introscopy systems S. Ogorodnikov , V. Petrunin , M. Vorogushin , Radioscopic discrimination of materials in 1/10 MeV range for customs applications . Proceeding of the 8th European Particle Accelerator Conference , Paris,France , 2002 . Radioscopic discrimination of materials in 1/10 MeV range for customs applications K. Fu , Performance enhancement approaches for a dual energy X-ray. Ph.D. Thesis , University of California , 2010 . Performance enhancement approaches for a dual energy X-ray. Ph.D. Thesis X.P. Wu , Research and application on dual energy X-ray method of material recognition. Ph.D. Thesis , Tsinghua University , 2004 . (in Chinese) Research and application on dual energy X-ray method of material recognition. Ph.D. Thesis X.P. Wu , Z.Q. Chen , X.W. Wang , Application of LCIS for material discrimination with dual energy method . Nucl Electron Detect Technol , 2005 , 25 : 782 - 784 . (in Chinese) DOI: 10.3969/j.issn.0258-0934.2005.06.055 http://doi.org/10.3969/j.issn.0258-0934.2005.06.055 Application of LCIS for material discrimination with dual energy method Q.H. Li , X.W. Wang , J.M. Li , et al . Research on multi-spectrum detector in High-energy dual-energy X-ray imaging system . Nucl Electron Detect Technol , 2008 , 28 : 821 - 825 . (in Chinese) DOI: 10.3969/j.issn.0258-0934.2008.04.037 http://doi.org/10.3969/j.issn.0258-0934.2008.04.037 Research on multi-spectrum detector in High-energy dual-energy X-ray imaging system S.W. Li , K.J. Kang , Y. Wang , et al . Employing cerenkov detectors in material effective material atomic number detection of dual-energy X-ray beams . Nucl Electron Detect Technol , 2010 , 30 : 1012 - 1015 . (in Chinese) DOI: 10.3969/j.issn.0258-0934.2010.08.003 http://doi.org/10.3969/j.issn.0258-0934.2010.08.003 Employing cerenkov detectors in material effective material atomic number detection of dual-energy X-ray beams Q.H. Li , J.M. Li , X.W. Wang , et al . Overlapped objects discrimination using dual-energy, high energy X-ray imaging . Journal of Tsinghua University: Science and Technology , 2008 , 48 : 1256 - 1259 . (in Chinese) DOI: 10.3321/j.issn:1000-0054.2008.08.008 http://doi.org/10.3321/j.issn:1000-0054.2008.08.008 Overlapped objects discrimination using dual-energy, high energy X-ray imaging X.W. Wang , H.Q. Zhong , Q.H. Li , et al . The pilot study in the small angle forward scattering used in high energy dual energy imaging method . Prog Report China Nucl Sci Technol , 2009 , 1 : 230 - 235 . (in Chinese) The pilot study in the small angle forward scattering used in high energy dual energy imaging method 10.1007/s41365-016-0019-4 Supported by National Natural Science Foundation of China (Nos. 11235007 and 10905030) 2015-01-05 2015-05-09 2015-05-14 2016-03-03 2016-02-20 2021-04-11 © China Science Publishing & Media Ltd. (Science Press), Shanghai Institute of Applied Physics, the Chinese Academy of Sciences, ChineseNuclear Society and Springer Nature Singapore Pte Ltd. 2016 This is a post-peer-review, pre-copyedit version of an article published in Nuclear Science and Techniques. The final authenticated version is available online at: http://dx.doi.org/10.1007/s41365-016-0019-4 http://dx.doi.org/10.1007/s41365-016-0019-4 . 2016 Nuclear Science and Techniques 01 27