Introduction

Computed tomography (CT) is widely used in industrial applications as a nondestructive testing technique [1-4]. However, when imaging plate-like objects such as fossils, paintings, composite panels in the aerospace industry, and printed circuit boards (PCB), it is difficult to obtain high-precision three-dimensional (3D) images with the commonly-used circular cone-beam CT (CBCT) owing to the limitations of the imaging space and radiation source energy [5-7]. At the same time, computed laminography (CL) only requires rays to pass through an object in the thickness direction, and thus has great potential for imaging plate-like objects [8]. Initially, CL could only record images of the focal plane of an object. With the development of computers, digital detectors, and CL reconstruction algorithms, CL can now obtain 3D images of objects such as in CT [9-11].

According to the difference in the scan trajectory, CL can be divided into translational CL [12, 13], rotational CL [14], and swing CL [15], among which rotational CL is widely used because of its strong adaptability, rich projection information, and the same resolution in the xy direction [16]. In terms of the scanning geometry, rotational CL is analogous to CBCT: the detector and X-ray source rotate 360° around the rotation axis (i.e., z axis) to collect projection information. However, the angles between the central ray and the rotation axis (i.e., tilt angle α in Fig. 1(a)) in both systems are different. In CBCT, the central ray is perpendicular to the rotation axis (α=90°), whereas the tilt angle in rotational CL is less than 90° (α<90°), as shown in Fig. 1(b). This characteristic enables only the X-rays to pass through a plate-like object in the thickness direction during the 360° scanning process [17].

Fig. 1Full-screen View

(Color online) Imaging geometric diagram of CT and CL: (a) CBCT, (b) rotational CL. O is the center of the object, S is the source, D is the center of the detector, and v, u are the vertical and horizontal axes of the detector local coordinate system, respectively. The ray passing through S, O, and D in sequence is called the central ray

In both CBCT and CL, flat-panel detectors are widely used. However, the detector in CBCT is set vertically and faces the source during rotation, making full use of the detector [18, 19]. However, in contrast to CT, there are various settings for flat-panel detectors in rotational CL. As shown in Fig. 2(a), in the first setting, the detector is parallel to the rotation axis, which is similar to CBCT. The second and third settings are commonly used. In the second setting, the detector is placed perpendicular to the central ray during rotation, as shown in Fig. 2(b). In the third setting, shown in Fig. 2(c), the detector is perpendicular to the rotation axis and with in-plane rotation such that its v-axis always points to the rotation axis. The fourth setting is our proposed method. Similar to the third setting, the detector is set perpendicular to the rotation axis. However, the detector only exhibits transitional motion, and the orientation of the detector remains unchanged during rotation. Different detector settings imply different scanning geometries, which have a direct impact on the image reconstruction [20].

Fig. 2Full-screen View

(Color online) System configuration with different detector settings: (a) detector parallel to the rotation axis, (b) detector perpendicular to the central ray, (c) detector perpendicular to the rotation axis and with in-plane rotation, (d) detector perpendicular to the rotation axis and with transition movement only. β is the angle between the projection line of the central ray on the xy plane and the positive y-axis direction

Image reconstruction is an important aspect of CL imaging [6, 21]. The existing CL reconstruction methods can be divided into three categories: analytical [16, 22], iterative [23, 24], and deep learning methods [2, 25]. Although some studies [17, 18, 26-28] have shown that deep learning methods have excellent performance in terms of computational efficiency and accuracy, there are still many challenges (e.g., lack of training data), and further optimization is needed before it can be extensively accepted. Analytical and iterative methods are widely used for practical applications. Iterative methods exhibit good noise resistance and the ability to process incomplete projection data. However, they need large calculations; thus, it is difficult to achieve real-time reconstruction. In contrast, analytical algorithms have less computational complexity and no parameters are required. Hence, they are widely used in commercial applications. However, analytical methods are specifically bound to imaging geometries, and different geometries require different analytical algorithms [29].

Different reconstruction methods have different application scenarios [30-33]. Although the reconstructed images of analytical methods have worse artifacts compared with those of iterative methods, they are suitable for scenarios requiring efficiency, e.g., online detection of circuit board defects. In an analytical algorithm of rotational CL, Yang et al. [34] proposed a filtering backprojection reconstruction formula suitable for rotational CL. However, this method focused only on the backprojection process and did not consider the filtering process. Sun et al. [16] proposed a reconstruction algorithm based on projection transformation (PT-FDK). In this method, the CL scanning data and parameters were converted into those of CT that conform to the FDK conditions [35]. Then, the filtering backprojection operation was carried out on the converted CL data. In this manner, the CL projection data was reconstructed using the standard FDK algorithm. Compared with Yang’s work, this method converts projection data to standard geometry and adopts standard FDK. Thus, it has high applicability. However, it requires a large amount of computation, and the interpolation error can be sufficiently large to degrade the image reconstruction.

In this study, for fast and high-precision imaging of circuit boards, we first propose a rotational CL detector setting and compare its field of view (FOV) with other detector settings. Subsequently, an FDK-type analytical reconstruction algorithm for the proposed detector setting is derived and verified through numerical experiments. Finally, the proposed rotational CL scheme is validated using a real system for PCB inspection.

FOV analysis with different detector settings

As shown in Fig. 3, during rotational CL imaging, the imaging range under the projection angle β is the quadrangular pyramid region formed by the X-ray source S and the four vertices of the detector. The intersection of the imaging ranges for all projection angles is the FOV of the CL imaging system. The projection information of voxel points within the FOV can be recorded by the detector at all projection angles. To ensure the quality of reconstruction in CL imaging, all voxels of interest must be located within the FOV. Therefore, a larger FOV enables larger objects to be scanned.

Fig. 3Full-screen View

(Color online) Schematic diagram of the field-of-view (FOV) at projection angle β

Under a projection angle β, let R₁, R₂, R₃, and R₄ be the intersection points of the rays SP₁, SP₂, SP₃, SP₄, and the z=z₀ plane, respectively. The quadrilateral region R₁R₂R₃R₄ is the imaging range of CL on the z=z₀ plane. Correspondingly, the intersection of the quadrilateral regions R₁R₂R₃R₄ of the CL at all projection angles is the FOV of the CL on the z=z₀ plane. In general, the FOV of the CL on the z=z₀ plane varies with the coordinates z₀. However, as the circuit board is small in the thickness direction (i.e., the z-direction), it is important to evaluate the FOV of the CL system during imaging on the circuit board by directly analyzing its FOV on the z=0 plane.

According to the derivation of the equation (see the Appendix), as shown in Fig. 4 with the first three settings, the shape of the imaging region R₁R₂R₃R₄ of the CL on the z=0 plane does not change with the projection angle and only rigidly rotates around the origin O during imaging. Therefore, their FOV shapes on the z=0 plane are circles, and the radii of these circles can be determined by finding the minimum distance from the origin O to the four sides (i.e., line R₁R₂, R₂R₃, R₃R₄, and R₄R₁) of the quadrangle region R₁R₂R₃R₄. In the fourth setting, the quadrangle R₁R₂R₃R₄ not only has a constant shape, but also does not rotate around the origin O. Therefore, the FOV is the quadrangle R₁R₂R₃R₄, which is a rectangle.

Fig. 4Full-screen View

(Color online) FOVs in z=0 plane: (a) first, (b) second, (c) third, and (d) the fourth settings. The red circle/box is the corresponding FOV

Let , , , and be the distance in the ith (i=1, 2, 3, 4) setting from the origin O to the line R₁R₂, R₂R₃, R₃R₄, and R₄R₁, respectively. As , , (the explanation is provided in appendix A), the circle radius of the first three settings can be expressed as:(1)In practical applications, the tilt angle of the CL is less than 60°, i.e., 0°<α≤60°. At this point, we obtain(2)Therefore, . The areas of the circular FOVs are(3)As the FOV shape in the fourth setting is rectangular, its area can be calculated as(4)To summarize, when 0°<α≤60°, S⁽¹⁾<S⁽²⁾<S⁽³⁾<S⁽⁴⁾. Hence, under the same imaging conditions, the fourth setting has the largest FOV, followed by the third, second, and first settings.

To more intuitively compare the FOVs, we compared the FOVs for different detector settings using a numerical test. In the numerical tests, four rotational CL systems with different detector settings (Fig. 2) were simulated using the ASTRA toolbox [36]. These systems had the same imaging parameters, except for the detector settings.

Table 1 lists these imaging parameters. In the simulation, if the projections of a reconstruction point are located inside the detector at all projection angles, then this point belongs to the FOV. For greater number of points, the FOV is larger.

To illustrate this more intuitively, Fig. 5 shows two mutually perpendicular sections and the areas of the four FOVs. Fig. 5(a) shows the coronal plane (i.e., yz cross-section, x=0 voxel), and Fig. 5(b) shows the transverse plane (i.e., xy cross section, z=0 voxel). The volumes of the FOVs are also provided. It is evident that in the FOV distribution, the shapes of the four FOVs are irregular on the coronal planes. Meanwhile, as shown by the theoretical equations, the xy cross sections of the first three FOVs are circular, whereas that of the fourth is a special rectangle: square because the detector has the same size in two directions. The volume of the FOV in the first setting is the smallest, followed by the third, second, and fourth settings. Although the volumes of FOVs in the second and fourth directions are similar, the second direction is more slender along the z-direction and is not suitable for imaging plate-like objects.

Fig. 5Full-screen View

(Color online) Comparison of FOVs in four settings: (a) coronal plane (x=0 voxel), (b) transverse plane (z=0 voxel)

From the above analysis, it can be concluded that the proposed setting has the largest FOV under the same imaging parameters. In addition, the shape of the xy cross-section is the largest and most rectangular, which is beneficial for imaging plate-like objects using CL such as circuit boards. Most of these objects are rectangular. Finally, the detector is horizontal and requires smaller installation space. As this setting has a rectangular FOV shape, it is named ‘rectangular cross-section FOV rotational CL (RC-CL) in this study.

Analytical reconstruction algorithm for RC-CL

As the imaging geometry in RC-CL is different from CBCT, the classical FDK algorithm cannot be used directly. It is possible to transfer the projection data of RC-CL to fit into the CBCT geometry by 2D interpolation such that FDK can be applied to reconstruction similar to PT-FDK. However, two unfavorable factors need to be considered: 1) transferring projections of RC-CL to CBCT may require a significantly larger virtual detector because RC-CL projections correspond to a large cone angle in CBCT. 2) 2D interpolation error in this situation can significantly reduce the image quality of the reconstruction. Therefore, an analytical reconstruction method specifically for RC-CL is necessary for efficient and effective reconstruction.

A 3D schematic of the RC-CL system is shown in Fig. 6(a), where a global coordinate system O-xyz is defined with z being the rotation axis, and the origin O is the intersection of axis z and the center-ray connecting the source (S) and the center of the detector (D). The zenith angle α is referred to as the CL tilt angle. The plane E includes the detector and O’ is the intersection of plane E and the rotation axis. S’ is the projection of S onto plane E. Fig. 6(b) shows a 2D schematic of the top view of plane E. A native coordinate system D–uv exists in the detector. During rotation, the directions of the axes u and v in D–uv remain parallel to the axes x and y, respectively.

Fig. 6Full-screen View

(Color online) RC-CL imaging: (a) 3D schematic diagram, (b) schematic diagram on plane E, and are the locations of the source/detector at two different projection angles, respectively

3.1

Formulation of analytical reconstruction on a virtual 2D problem

To derive the reconstruction formula for RC-CL, we followed the idea of the FDK method and started with a 2D filtered back projection (FBP) reconstruction on plane E. First, a rotational coordinate system D-u’v’ is configured on the detector with axis v’ pointing to the rotation axis during rotation to form a 2D virtual problem. For convenience, we define the angle between the axes v’ and v as the projection angle β. The relation between D-uv and D-u’v’ at the projection angle β is:(5)As shown in Fig. 7, on the detector plane E, if we regard point S’ as the X-ray source in the 2D problem, the projection data along v’ at a certain provides a standard view of the fan beam CT. Hence, we can apply an FBP reconstruction algorithm in this situation.(6)where (x,y) are the coordinates of the reconstruction point R. is a point on the v’ axis, where . is the projected distance of on ; and are the projection positions of point R on the detector; and are the distances from the source S’ to the detector center D and origin O’, respectively; represents projection data under the coordinate system D-u’v’, and is a ramp filter.

Fig. 7Full-screen View

(Color online) Geometry in the 2D situation with D, O’, and S’ representing the locations of the detector center, isocenter, and source, respectively, and P is the projection of R on the detector

In Eq. 6, FBP filtering is performed along the axis. However, the projection data in RC-CL is recorded along the u and v axes. Therefore, a reconstruction formula for RC-CL is derived based on Eq. 6.

To simplify the derivation, let’s define(7)This can be further expressed as:(8)Substituting Eq. 5 into Eq. 8 yields(9)where is the projection data on the physical detector grids recorded in the coordinate system D–uv; u^* and u^* are the corresponding coordinates in the D–uv coordinate system of and .

In Eq. 9, based on the scaling property of the Dirac delta function, we obtain(10)By substituting Eq. 10 into Eq. 9, we obtain(11)From ,(12)By substituting Eq. 12 into Eq. 11,(13)According to the Fourier transform property, , we obtain(14)By substituting Eq. 14 into Eq. 13,(15)By combining Eq. 15 and Eq. 6, we obtain an FBP-type reconstruction formula for RC-CL:(16)

3.2

Extension to a 3D scenario

In a 3D situation, the z-dimension must be considered during reconstruction. Similar to the derivation of the FDK algorithm for CBCT, when extending the FBP algorithm from 2D to 3D in RC-CL, two aspects in FBP need to be modified.

The first is the weighting factor before the filtering operation (recorded as η₁). According to Eq. 16, the expression for η₁ in FBP is(17)Physically, in the 2D case, η₁ represents the cosine of the fan angle (i.e., γ in Fig. 7) of the reconstruction point R. In the 3D case, as shown in Fig. 8, the fan angle of the reconstruction point R is . The influence of the cone angle (i.e., in Fig. 8) needs to be considered. Therefore, the expression for η₁ is(18)According to the cosine theorem, Eq. 18 can be written as(19)In the global coordinate system O-xyz, the coordinates of point Os, S, P, and Ps are , , , , respectively. Eq. 19 can be expressed as(20)By substituting Eq. 12 into Eq. 20, we obtain(21)The second is the weighting factor for backprojection (recorded as η₂). According to Eq. 16, the expression for η₂ in FBP is(22)Physically, η₂ is determined by the source-to-detector distance and the projected distance between the source and reconstruction point on the central ray. Therefore, as illustrated in Fig. 8, the expressions for the 3D case is(23)According to the triangle similarity theorem,(24)where z is the coordinate of the reconstruction point R.

Therefore, Eq. 23 can be written as(25)Replacing η₁ and η₂ with Eq. 21 and Eq. 25, the FDK-type reconstruction formula for RC-CL can be obtained.(26)The implementation of the proposed algorithm can be summarized in three steps.

Fig. 8Full-screen View

(Color online) Geometry in the 3D situation, the subscript s represents the projection of each point on the plane where S is located, Q is the point on the v’ axis with , R is the reconstruction point and its projection on detector is P, R₁ is the projection of point R on line SQ and R₂ is the projection of R₁ on the S-located plane

1. Preweighting: Multiply the two-dimensional projection data by a weighting factor computed by(27)

2. Filtration: In numerical implementation, to lower the discretization error and avoid cos β = 0 at or , we divide [0, 2π) into four parts, , , , to perform filtration, as shown in Fig. 9.

Fig. 9Full-screen View

(Color online) The diagram of filtration in RC-CL. The brown area represents filtration along the u axis, the blue area represents filtration along the v axis. The gray value at point P can be obtained through 1D linear interpolation of the values of point P₁ and point P₂ on the grid

Specially, when , . (i.e., Eq. 12) is adopted to replace v in Eq. 11. The derivations of Eq. 13–Eq. 26 are based on this scenario. Accordingly, the u coordinates can be specified manually, and filtration is performed along the u-axis.

Meanwhile, when , . can be used to replace u in Eq. 11. For conciseness, we do not provide a detailed derivation here, but the reader can easily derive it according to Eq. 13–Eq. 26. Now, the v coordinates can be specified manually, and filtration is performed along the v-axis.

During filtration, although the u and v coordinates can be assigned as integers when filtering along the u and v axes, respectively, the corresponding v and u coordinates need to calculated as and , which are usually not integers. Therefore, interpolation is required to obtain the projection values at these positions. However, in contrast to 2D interpolation in the projection data transfer algorithm, we only need 1D interpolation, which is easy to perform as shown in Fig. 9.

3. Weighted back projection: The 3D back projection weighted by is similar to those of other FDK type reconstructions.

Experimental analysis

4.1

Simulation study

To verify the proposed reconstruction method (referred to as CL-FDK), we simulated an RC-CL system. In the simulation, the radiation source was regarded as a point source. Ray- and voxel-driven models were chosen as the forward and back-projectors, respectively. The printed PCB phantom shown in Fig. 10 was used. The phantom contained three interconnected copper circuit layers. The mass attenuation coefficients were from the table of X-ray mass attenuation coefficients of the National Institute of Standards and Technology (NIST) [37, 38], and the range was [0.05, 0.46]. The detailed imaging parameters were as follows: the title angle α was 45°, and the distances from the source to the origin and detector center were 45.79 mm and 194.58 mm, respectively. The detector was simulated using a 768 pixel×768 pixel array and 0.17 mm×0.17 mm pixel size, and 256 projection images were acquired. The reconstructed image grids were 300 voxel×300 voxel×80 voxel with a 0.07 mm×0.07 mm×0.07 mm voxel size.

Fig. 10Full-screen View

(Color online) Numerical PCB phantom for the simulation study

To quantitatively evaluate the quality of the reconstructed PCB images, three metrics were used to measure the similarity between the reconstructed and reference images: root mean square error (RMSE), mean structural similarity index (MSSIM) and peak signal-to-noise ratio (PSNR). For a smaller RMSE, the reconstruction quality was better. In contrast, the reconstruction quality was better for a larger MSSIM and PSNR.

For comparison, the reconstruction results obtained using the PT-FDK method [16] and the simultaneous iterative reconstruction technique (SIRT) are also presented with the SIRT iteration count set to 200. The reconstructed results of the three algorithms are shown in Fig. 11 and Fig. 12. It is evident that all three methods can reconstruct the main structural features in the phantom. The reconstructed 2D cross-sectional images of slice #40 (i.e., z=40 voxel) shown in Fig. 12 clearly show that artifacts are unavoidable because of incomplete CL scan data. All reconstructed images are darker than the original images. The difference between the reference and reconstructed slice images shows that the SIRT result has the least artifacts and the PT-FDK result has the most significant error. The horizontal profiles along lines a and b in Fig. 12 are shown in Figs. 13, which confirm that the difference between the intensity of SIRT reconstruction and that of the phantom is the smallest, and the error in the CL-FDK result is smaller than that in PT-FDK.

Fig. 11Full-screen View

(Color online) Volumetric rendering of reconstruction results from the three algorithms: (a) SIRT, (b) PT-FDK, (c) CL-FDK

Fig. 12Full-screen View

(Color online) Reconstruction result at slice #40: (a) reconstructed image, the display window is [-0.24, 0.57], (b) difference between the reconstructed image and reference, the display window is [-0.43, 0.31]

Fig. 13Full-screen View

(Color online) Horizontal profiles: (a) along line a, (b) along line b

Figure 14 shows the values of the three metrics for different reconstruction methods. It can be observed that the SIRT reconstruction method is the best, followed by CL-FDK and PT-FDK. As a filtering backprojection algorithm, the PT-FDK algorithm has lower accuracy than CL-FDK mainly because PT-FDK requires out-of-plane interpolation of the projection image during the transformation process, and the additional interpolation operation not only increases the computational load, but also introduces interpolation errors.

Fig. 14Full-screen View

(Color online) Comparison of evaluation indicators: (a) RMSE, (b) MSSIM

4.2

Influence of tilt angle α

The tilt angle α is an important parameter for CL imaging. To study its influence on RC-CL, we experimented by setting the tilt angle α to 25°, 35°, 45°, 55°, and 65°, while keeping the other parameters constant. Figure 15 presents the reconstruction results for slice #30 and their differences from the reference. As shown, the CL-FDK can reconstruct the main internal features of the phantom at different tilt angles. However, more severe artifacts are evident for a smaller α. The superimposed structure from the other layers become weaker when the tilt angle is increased.

Fig. 15Full-screen View

(Color online) Reconstruction results under different tilt angles:(a) reconstructed image, the display window is [-0.45, 1.18], (b) difference between the reconstructed image and reference, the display window is [-0.68, 0.93]

Figure 16 shows the variations in the quantitative metrics with respect to the tilt angle. The RMSE decreases with an increase in the tilt angle. Its value at 25° is 1.69 times larger than that at 65°. However, the RMSE of MSSIM increases with an increase in the tilt angle.

Fig. 16Full-screen View

(Color online) Comparison of evaluation indicators: (a) RMSE, (b) MSSIM, (c) PSNR

4.3

Real experimental study

In this section, we present the experimental CL scan of a PCB sample. The experiment was performed using an RC-CL system, as shown in Fig. 17(a). A microfocus X-ray source was used as the radiation source in the experiment. Its main characteristics included: X-ray tube voltage operational range 60 kVp to 110 kVp, X-ray tube current operational range 10 μA to 800 μA, and X-ray focal spot size (nominal value) of 4 μm. In this study, the X-ray tube was set at 80 kVp and 20 μA. The scanned PCB sample was an L-shaped computer motherboard (Fig. 17(b)). As the bottom of the PCB sample was not flat, it was placed on an aluminum base during imaging. In the experiment, the tilt angle was set to 45° and 512 projections uniformly distributed over 2π were acquired. Each projection had 2048×2048 detector bins, and the bin size was 0.14 mm × 0.14 mm. The distance between the source and detector was 263.101 mm, and the distance between the source and origin was 28.681 mm.

Fig. 17Full-screen View

(Color online) Real experiment: (a) CL system, (b) PCB sample

Owing to the large size of the PCB, we selected only a few representative areas of the sample for imaging during the experiment. Figure 18(a) shows the reconstructed results for ball grid array (BGA) solder joints. Multiple bubble defects featured by black holes can be observed in the reconstructed images, for example, the ones indicated by red arrows in the figure. Figure 18(b) shows the reconstructed results for the quad-flat no-lead (QFN) package solder joints. The square area in the image represents the QFN solder joints and the irregular circular area inside represents the internal bubble defects. Figure 18(b) also shows the gray-value profile along the yellow line in the three QFN, and the change pattern of the gray values is highly correlated with the location of defects. Based on these results, it can be concluded that the proposed CL-FDK algorithm can accurately reconstruct the main internal features of the tested objects and can be applied to real systems.

Fig. 18Full-screen View

(Color online) Reconstruction result of a PCB: (a) two BGA images with multiple bubble defects, (b) three QFN images with multiple bubble defects

Conclusion

This study proposes a new rotational CL imaging system with a horizontal and fixed-orientation detector, and develops a suitable analytical reconstruction algorithm. The results showed that the proposed imaging system had the largest field-of-view (FOV) under the same conditions. However, it exhibited superior performance to the commonly-used projection resorting reconstruction algorithm. On this basis, the influence of the tilt angle on the reconstruction result was analyzed and a larger tilt angle was suggested for better performance. Finally, the proposed imaging system and its reconstruction algorithm were validated on a system to image circuit boards for defect detection.

Although the proposed method provides a new method for plate-type object 3D imaging, owing to the intrinsic shortcomings of rotational CL (i.e., a lack of projection information under the same angles), the reconstructed image contains interlayer aliasing artifacts, which are difficult to eliminate through traditional methods. In this situation, deep-learning methods can be a good choice. We intend to conduct further studies on this topic in the future.

Meanwhile, it should be noted that although motion and scattering artifacts are two common artifacts in CT/CL imaging [39], they were not the main errors in this study. For motion artifacts, the step-and-shoot mode was used to record the projection images, and the motions of the detector and source were well-controlled. Therefore, the motion artifact in the reconstructed image was negligible. However, if continuous mode is used to record projection images, for example, in the online detection of circuit boards, the influence of motion artifacts cannot be ignored. However, in this study, scattering artifacts were not considered as the research object was a circuit board with a high contrast ratio and the effect of scattering on the reconstruction was small. However, we also found that the CT value of the air region was not zero during CL image reconstruction, indicating the presence of scattering artifacts in CL imaging. Hence, we conclude that when using CL for objects with a low contrast ratio, a detailed analysis of scattering artifacts is essential.