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

LOW ENERGY ACCELERATOR, RAY AND APPLICATIONS

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

Zhi-Qiang Chen，

Tiao Zhao，

Liang Li

Nuclear Science and Techniques

Vol.27, No.1

Article number 25

Published in print 20 Feb 2016

Available online 03 Mar 2016

DOI：10.1007/s41365-016-0019-4

46400

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 energyDual-energyDigital radiographyMaterial recognition

1 Introduction

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-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]. 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].

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]. 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

2.1 The principle of MeV dual-energy X-ray imaging in material recognition

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]. The corresponding attenuation coefficients, μ, have different dependences with a material atomic number, Z. We can give

\begin{array}{l} μ_{P} \propto Z^{4} < 1, \\ μ_{CS} \propto 1 1 \sim 3, \\ μ_{EP} \propto Z > 3, \end{array}

(1)

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 Em, a single substance with an atomic number of Z, an attenuation coefficient function of μ (E, Z), and a thickness t, the transparence, T is

T = \frac{I}{I_{0}} = \frac{\int_{0}^{E_{m}} N (E) e^{- μ (E, Z) t} d E}{\int_{0}^{E_{m}} N (E) d E} .

(2)

In a dual-energy situation, the boundary energy of the two X-ray sources are E₁ and E₂. We defined the logarithmic ratio of T as

R = \frac{\ln T_{1}}{\ln T_{2}} = \frac{\bar{μ} (E_{1}, t, Z)}{\bar{μ} (E_{2}, t, Z)},

(3)

where R is the ratio of the equivalent attenuation coefficient, $\bar{μ}$ . 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₁ is in the low energy range and E₂ is in the middle energy range, R can be written as

R = \frac{\ln T_{1}}{\ln T_{2}} = \frac{\ln e^{- μ (E_{1}, t, Z) t}}{\ln e^{- μ (E_{2}, t, Z) t}} = \frac{μ (E_{1}, Z)}{μ (E_{2}, Z)} \approx \frac{μ_{P} (E_{1}, Z)}{μ_{CS} E_{2}, Z} \approx c (E_{1}, E_{2}) Z^{4} .

(4)

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₁ and E₂ 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₁ to the high energy range and keep E₂ in the middle energy range, then Eq. (4) is

R = \frac{μ (E_{1}, Z)}{μ (E_{2}, Z)} \approx \frac{μ_{EP} (E_{1}, Z)}{μ_{CS} (E_{2}, Z)} \approx c (E_{1}, E_{2}) Z,

(5)

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]. 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]. Here E₁ and E₂ are usually higher than 3 MeV.

2.2 A MeV dual-energy system

We use a schematic model (See 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.

Fig. 1.

An schematic model of the cargo inspection system.

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

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).

T (x, y) = \frac{I (x, y) - I_{BK} (x, y)}{I_{0} (x, y)} \cdot L D (x) .

(6)

The coordinate x, y represents the lateral and vertical position. I is the original dual-energy image. I₀ 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₀ 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

Four kinds of classification curves, R-curve (Fig 2(a)), T₂ - T₁ curve (Fig. 2(b)), α-curve (Fig. 2(c)) [10, 11], and a separability curve which denotes the transparence difference of the two materials (Fig. 2(d)), are shown in 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 are ideal. The longitudinal coordinate of the R-curve is R, defined in Eq. (3), and the horizontal coordinate is the inverse of T₁ or T₂. Different materials have different R-curves. They are arranged in Z’s increasing order when they are in the same image [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₂ (or C), as we can see in Fig. 2.

Fig. 2.

(Color online) Four kinds of classification curves. The R curve has the best visual separability.

(x,y) is one point of the dual-energy transparence image. There are two transparence values, T₁(x,y) and T₂(x,y). They give a point, (R, 1/T₂), 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₁, T₂, and each point, (x,y), belongs to an R-curve, Rxy. 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

R_{x y} = R_{y}, y = y_{1}, y_{2}, \dots, y_{n} .

(7)

The data is

\begin{array}{l} {data}_{y} = {R_{y} [T_{1} (x_{i}, y), T_{2} (x_{i}, y)], i = 1,2, \dots, m}, \\ y = y_{1}, y_{2}, \dots, y_{n}, \end{array}

(8)

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.

\begin{array}{l} {R_{y}^{Z} | R_{y}^{Z} = fit ({data}_{y}^{Z}), y = y_{1}, y_{2}, \dots, y_{n}}, \\ C_{y_{i}} = {R_{y_{i}, Z} | Z = P b, F e, A l, C}, i = 1,2, \dots, n . \end{array}

(9)

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

\begin{array}{l} T_{1 r}^{Z} (x_{i}, Y), T_{2 r}^{Z} (x_{i}, Y), i = 1,2, \dots, n x, \\ Z = P b, F e, A l, C, \end{array}

(10)

where nx is the number of horizontal pixels. In the first step, this data is also saved as $T_{1 b}^{Z} (x_{i}, Y), T_{2 b}^{Z} (x_{i}, Y)$ .

They are the average of T on the horizontal direction, as we ignore the lateral difference. Then the revised classification curve will be

C_{x_{j} y_{i}} = {R_{x_{j} y_{j}}^{Z} | R_{x_{j} y_{j}}^{Z} = R_{y_{i}}^{Z} * F (T_{1 r}^{Z}, T_{2 r}^{Z}, T_{1 b}^{Z}, T_{2 b}^{Z}; x_{j}, Y), Z = P b, F e, A l, C}, i = 1,2, \dots, n, j = 1,2, \dots, n x .

(11)

The function F uses real-time calibrating data to estimate the difference between $R_{x v}^{Z}$ and $R_{v}^{Z}$ . The correction method in Eq. (11) relies on how function F is calculated. Considering the statistical straggling of T(xj, Y), we choose a segment of xj and use the average of the monitoring data to correct the classification curve.

2.5 Smooth of material information image

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, 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]. 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

The idea of colorization and the IHS color space was proposed in Ref. [9, 12]. We use a similar approach applying HLS color space [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,

Z < Z_{1}; Z_{1} \leq Z_{2}; \dots; Z_{p - 1} \leq Z_{p}; Z_{p} \leq Z .

(12)

There are p+1 hues H₁, H₂, …, 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 Hj.

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

Y = f (H, L, S) .

(13)

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

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]. 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, 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.

Fig. 3.

(Color online) The ’true’ curves in two states are C₁ and C₂. The revised R curve using the real-time calibration data to estimate the ’true’ curve in state 2 is C_m2.

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². 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). 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, 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) 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) 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) 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. Besides Pb, the other materials’ results show that our calibration strategy enhances classification accuracy.

Fig. 4.

(Color online) Comparison of the results with the real-time calibration (a) and without the real-time calibration (b), and (c) is the color table.

Notice that all objects in 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, 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.

Fig. 5.

(Color online) Comparison of the results with smoothing (a) and without smoothing (b), and (c) is the color table.

In 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.

Fig. 6.

(Color online) Comparison of cargo inspection result with smoothing (a) and without smoothing (b). From left to right, the irradiated objects are cigarette, salt, sugar, coffee, buckets of water and concrete. The circled regions clearly show the smoothing effect.

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

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]. Use the Cerenkov detector, which has a threshold and a good response to high energy X-rays [17]. Add the obstacle’s classification curve on the blocked material’s classification to enhance the blocked material’s recognition accuracy [18]. Use small angle scatter to realize the material recognition [19]. These methods can be applied to our system or a new imaging system based on our model and data processing flow.

References

[1]

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