Simulation of X-ray Computed Tomography

An application of the Radon transform

Sébastien Gilbert
Towards Data Science

--

Figure 1: Animation showing the Radon transform. Top: input image. Bottom: corresponding Radon transform. Generated by this code. Image by the author.

In the electronics manufacturing industry, it is common practice to inspect the interconnections between two discrete components with X-ray computed tomography. A ball grid array (BGA), for example, is a dense grid of small solder balls connecting a printed circuit board (PCB) and a chip. We want to make sure that each ball has soldered correctly to the metal pads of the chip and the PCB. If we took a perpendicular X-ray of the assembled package, the same way your dentist grabs an image of your teeth, the resulting image quality would suffer from the multiple layers of metal structures in the chip and the PCB. Each metal structure acts as an X-ray attenuator, creating occlusions that might prevent the detection of subtle defects in the BGA plane (typically isolated solder satellites, bridges, and abnormal balls).

Figure 2: The Ball Grid Array (BGA) being between the PCB and the chip, it is not entirely visible. Image by the author.

The BGA inspection scenario

You are the computer vision engineer in charge of inspecting the BGA’s on the assembly line, right after the components have been soldered in a furnace. You considered taking pictures from the side of the array, but you quickly realized that you can’t see any further than the first row of connections. You are searching for defects that can be found anywhere under the chip, so you’ve ruled out this option.

After some extensive googling, you come across a solution based on X-ray computed tomography. Being an outstanding maker, you build a one-dimensional X-ray imaging system that captures a slice of the BGA absorption at multiple angles (Cf. Figure 3).

Warning, kids: Don’t do this at home, unless under the supervision of an adult holding a Ph.D. in high energy radiation.

After thinking twice about it, don’t do this, period.

Figure 3: The X-ray imaging system. Image by the author.

Your imaging system produces a slice of the BGA integrated density along the length of a linear array of 363 X-ray sensors¹. A rotation mechanism can select the angle θ of the projection. You can find the code that simulates this system here.

After days of alignment tweaks, you finally manage to obtain an image built by horizontally stacking the projections of a BGA, at 256 discrete angles, from -π/2 to π/2:

Figure 4: The BGA plane sinogram. Image by the author.

We call such an image a sinogram because a point in the original density image creates a sinusoid in the image of stacked slices.

What do you do now? How does the sinogram of Figure 4 reveal the density map on the BGA plane?

The Radon Transform

The operation of creating a linear profile by integrating a 2D density map over parallel rays at a given angle is a Radon Transform.

Figure 5: The value of the Radon transform in (ρ, θ) is the integral of f(x, y) along the line defined by (ρ, θ). Image by the author.

The Radon transform of f(x, y) is defined by:

Our setup of rotating source and sensor shown in Figure 3 performed the Radon transform of the function corresponding to the X-ray local absorption coefficients of the BGA plane. Since we are interested in the X-ray absorption map, we need to compute the inverse Radon transform.

The Inverse Radon Transform

We can get an approximation of the inverse Radon transform by making use of the Fourier slice theorem²:

For a given θ=θ₀, the 1D Fourier transform of Rf(ρ, θ₀) (integrating over ρ) is the slice view of F(u, v) (the 2D Fourier transform of f(x, y)) passing through the origin and along θ₀ + π/2.

Each column in Figure 4 generates a 1D spectrum that gets pasted in a 2D spectrum F(u, v). Having a large number of projections, we can build a “radially dense” 2D spectrum (Cf. Fig. 6).

Figure 6: The Fourier slice theorem. Image by the author.

Do you mean we must paste 1D spectrums at arbitrary angles on a 2D spectrum defined on a rectangular grid?

Yes. As you can guess, the inverse Radon transform is an interpolation festival. Fortunately, the skimage library has done this work for us, and computing the inverse Radon transform is as simple as:

from skimage.transform import iradon
reconstruction = iradon(sinogram, theta=thetas, filter_name='ramp')

The BGA density map

Applying the inverse Radon transform to the sinogram of Figure 4, we obtain the following density map:

Figure 7: The inverse Radon transform of Figure 4. Image by the author.

The reconstruction is a bit noisy because the inverse Radon transform is not exact. The quality is affected by the number of projection angles and the interpolation method used to match the rectangular lattice of the 2D spectrum with the polar coordinates of the 1D slices.

Nevertheless, we can see the 16 x 16 individual metal connections of the BGA. There is a small dense object that is not supposed to be there, to the left of the image center. We found a microbead, one of the defects affecting a BGA, that we could only see with computed tomography.

Conclusion

We simulated the inspection of a BGA with computed tomography. That led us to introduce the Radon transform and its inverse operation.

We built the sinogram of a BGA by integrating projections at multiple angles. The inverse Radon transform of the sinogram revealed the object density map.

Bonus game

In a sinogram, a small object appears as a sinusoid. The microbead in our BGA breaks the symmetry. Can you perceive the faint sine wave that looks out of phase, in Figure 4?

[1] We can model the total absorption of X-ray along the path as the product of exponential absorptions over a series of short propagation lengths Δx where the local absorption coefficient is αi:

Figure 8: The X-ray intensity decreases as the beam propagates through the object under test. Image by the author.

The ratio of the detected intensity Id over the source intensity I₀ can be written as:

Equation (4) tells us that -ln(Id/I₀) is proportional to the sum (the integration in the limit case) of the local absorption coefficients αi encountered by the X-ray along its path across the material. That is why we can consider the X-ray projections as integrations of the density map in the object plane.

[2]

Proof of the Fourier Slice Theorem

Figure 9: The density map f(x, y) and its 2D Fourier transform F(u, v). Image by the author.

We need to show:

Inserting (1) into the definition of a 1D Fourier transform:

Since the change of variables from (z, ρ) to (x, y) corresponds to a rotation of the axes without scale change, the integration over the whole (z, ρ) plane is the same as the integration over the whole (x, y) plane.

--

--