Direct integration of the inverse Radon equation for X-ray computed tomography

1 Introduction

The image reconstruction method based on the inverse Radon equation is the most popular one in X-ray tomography [1–5]. It is used widely on software and hardware levels of numerical implementation. Radon’s inversion equation, which gives the explicit expression for the phantom density, is usually done by using Fourier transformation. But Radon [1] derived his inversion formula otherwise without Fourier transformation. In this article, we propose a direct integration method for calculating Radon’s inversion equation, without using Fourier transformation. As shall be seen in the examples they we consider, the proposed direct integration method provides more accurate results than standard methods such as filtered back-projection (which uses Fourier transformation).

First, let us consider the derivation of Radon’s inversion equation. Radon’s equation has a clear geometrical sense. Let’s first reproduce the derivation of this equation on base of the original Radons article [1]. Consider a double integral on a plane: J ( q , x , y ) = ∫ ( x - α ) 2 + ( y - β ) 2 > q 2 μ ( α , β ) d α d β ( x - α ) 2 + ( y - β ) 2 - q 2 .(1)

This integral is calculated over all plane with the exception of the circle of radius q and with center (x, y). If we change the Cartesian coordinates of the integral to ray coordinates, which are related by α = x - s sin φ + q cos φ and β = y - s cos φ + q sin φ, then we get J ( q , x , y ) = ∫ 0 2 π d φ ∫ s > 0 μ ( x - s sin φ + q cos φ , y - s cos φ + q sin φ ) ds .(2)

The meaning of Equation (2) is that the value of J (q, x, y) is the sum of all ray-sums for rays tangent to the circle of radius q. This value can always be calculated (using the sinogram), and this allows us to consider the dependance on J (q, x, y) known. The problem is now reduced to the determination of the density μ (α, β) from the integral Equation (1). To do this, it is necessary to make another change of variables, introducing polar coordinates (r, ψ), related by α = x + r cos ψ and β = y + r sin ψ, and then we get J ( q , x , y ) = ∫ r = q ∞ r dr r 2 - q 2 ∫ 0 2 π μ ( x + r cos ψ , y + r sin ψ ) d ψ .(3)

The internal integral in Equation (3) corresponds to the mean optical density of the phantom on circles with radius r > q, i.e. μ ˆ ( r , x , y ) = 1 2 π ∫ 0 2 π μ ( x + r cos ψ , y + r sin ψ ) d ψ .(4)

Thereby, the expression in Equation (3) is an integral equation for the unknown value of the mean optical density of the phantom, and we have ∫ r = q ∞ μ ˆ ( r , x , y ) r dr r 2 - q 2 = J ( q , x , y ) 2 π .(5)

Equation (5) is Abel’s equation, which has analytical solution of the form μ ˆ ( r , x , y ) = - 1 π 2 ∫ q = r ∞ ∂ J ∂ q q 2 - r 2 dq ,(6) where now r < q.

The Equation (6) is the result of Radon. There q is the distance from the point (x, y) to the respective beam in the phantom plane, and r is the radius of the circle with same center (x, y). As it was noted, the value J (q, x, y) is the sum of all ray-sums for rays tangent to the circle with center (x, y) and radius q. But in tomography, each ray is characterized by the number of its detector, p, and by the direction of the radiation source, φ. The two-dimensional array of points (p, φ) forms a space of rays passing through the detectors, and the measured ray-sums are stored in the sinogram S (p, φ). To write Equation (6) in terms of beams (p, φ), it is necessary to return to formula (2) and write its right-hand side by J ( q , x , y ) = ∫ 0 2 π S ( p , φ ) d φ .(7)

Substitution of Equation (7) into Equation (6) yields μ ˆ ( r , x , y ) = - 1 π 2 ∫ 0 2 π d φ ∫ q = r ∞ ∂ S ∂ p q 2 - r 2 dp .(8)

The formula (8) calculates the mean density of the phantom on circles of radius r and center (x, y) via the ray-sums contained in sinogram S (p, φ). Suppose that the relation between distance q and values p, φ, x, y is known, that is, q = q ( p , φ , x , y ) .(9)

This relation is defined from the type of phantom’s scanning. For r → 0, the mean density is naturally converted to the phantom’s density in point (x, y), and we get μ ( x , y ) = - 1 π 2 ∫ 0 2 π d φ ∫ - L L ∂ S ( p , φ ) ∂ p dp q ( p , φ , x , y ) ,(10) where q (p, φ, x, y) has to be defined for each scanning geometry (i.e., parallel rays, fan-beam, etc.) separately. The inner integral in Equation (10), namely ∫ - L L ∂ S ( p , φ ) ∂ p dp q ( p , φ , x , y ) displaymath corresponds to the convolution operation. There, q (p, φ, x, y) is the distance from the point (x, y) to the beam (p, φ). Equation (10) is Radon’s inversion equation; we have the following properties:

For a fixed value of the angle φ = φ _k , we get an image of the k-th projection μ _k  (x, y).

The rest of the article is organized as follows. In the next section, we provide a complete derivation of adapting Equation (10) to the special case of fan-beam scanning. Then, in the third section, we describe a direct integration method for calculating Radon’s inversion equation. To show the results obtained through the proposed direct integration method, we provide two examples in the fourth section. In the last section, we provide conclusions and ideas for using the proposed method in practice.

For the fan-beam scanning each ray passes through the source of the radiation z₀ and the point z₁ lying on the line of detectors (Fig. 1). It is necessary to obtain a mathematical expression for the distance q from an arbitrary point z to the beam z₀z₁. The dashed line denotes the area of the phantom, and point C is the center of this area. Let R be the distance from the radiation source to point C, d the distance from line of detectors to the center C. Then the positions of the points z₀ and z₁ in the complex plane are defined by the expressions: z 0 = Re i φ , z 1 = e i φ ( - d + ip )(11)

A normal n to the ray z₀ - z₁ equals z 0 - z 1 | z 0 - z 1 | i = p cos φ - ( R + d ) sin φ p 2 + ( R + d ) 2 + i p sin φ + ( R + d ) cos φ p 2 + ( R + d ) 2 .(12)

The distance q is the scalar product of the vector z₀ - z on the normal n. Performing the calculations, we find q ( p , φ , x , y ) = R - x cos φ - y sin φ p 2 + ( R + d ) 2 [ p - ( R + d ) y cos φ - x sin φ R - x cos φ - y sin φ ] .(13)

Radon’s inversion Equation (10) for the case of fan-beam scanning has the following form: μ ( x , y ) = - 1 π 2 ∫ 0 2 π d φ ∫ - L L ∂ S ∂ p p 2 + ( R + d ) 2 R - x cos φ - y sin φ dp p - ( R + d ) y cos φ - x sin φ R - x cos φ - y sin φ .(14)

The case of parallel-beam scanning is obtained from (15) when R→ ∞: μ ( x , y ) = - 1 π 2 ∫ 0 2 π d φ ∫ - L L ∂ S ∂ p dp p - ( y cos φ - x sin φ ) .(15)

The formula (16) is the classic one, and is often mentioned in the literature in connection with its application in tomography.

Integral expression (15) or (16) is a specific algorithm for image restoration from the sinogram. It is only necessary to pay attention to the correct calculation of the integral with respect to the variable p included in the formula. We write it in the form: S 1 ( p 0 ) = ∫ - L L ∂ S ∂ p dp p - p 0 .(16)

Such a singular integral with difference kernel, called convolution, needs to be understood in the sense of the principal value. Such integrals are difficult to calculate by usual numerical methods, and therefore for their calculation are used Fourier series and the associated spectral theorem. This approach does not work quite reasonable for discrete data and the accompanying known Gibbs phenomenon can lead to significant errors. Therefore, for the singular integrals of convolution type with finite limits, it is desirable to consider an alternative method of calculating them without the use of spectroscopic methods.

Suppose we want to calculate the following integral, considered as a function of the variable p₀: S ( p 0 ) = ∫ - 1 1 F ( p ) dp p - p 0 .(17)

Such integrals are common in mathematics and physics, and they are usually written as a sum of two integrals: S ( p 0 ) = ∫ - 1 1 F ( p ) - F ( p 0 ) p - p 0 dp + F ( p 0 ) ∫ - 1 1 dp p - p 0 .(18)

The first integral in this formula can always be calculated suitable numerical method, and the second one appears in an explicit form, i.e. S ( p 0 ) = ∫ - 1 1 F ( p ) - F ( p 0 ) p - p 0 dp + F ( p 0 ) ln | 1 - p 0 1 + p 0 | .(19)

From (19) it follows that S (p₀) has a logarithmic singularity, if p₀ lies on the boundary of the integration interval. Further, we note that the integrand of the integral function (19) is defined for all values of the argument p, except for the values p = p₀, where there is uncertainty, as the numerator and denominator in it vanish simultaneously. Revealing this uncertainty, we get: F ( p ) - F ( p 0 ) p - p 0 | p = p 0 = F ′ ( p 0 ) .(20)

In the numerical implementation of the formula (19), it is convenient to calculate convolution not for a single fixed value of the argument p₀, but immediately for all its discrete values that coincide with the variable of integration p = p _k , for all k. To calculate the integral in (19), we apply the following sequence of operations:

The values of F (p), for all p = p _k , are determined by its derivative F′ (p), using difference approximations.

After numerical integration of matrix C (by the trapezoid method) along columns, we obtain the convolution S (p₀) represented as a column vector.

The proposed direct integration method for calculating integral (20) is very fast. For a cross-section image of size about 250 by 250 pixels, the method usually takes around the same time as the filtered back-projection method (around 1 second if implemented, for example, in the C language on a common CPU). For showing the results of the proposed direct integration method, we take the following two example.

Consider the cross-section image shown in Fig. 2 of size 250×250 pixels; it is a cross-section of a pipe welding process. As shown in Figs. 3 and 4, the direct integration method (60 projection angles, 359 detectors) gives more details compared to the filtered back-projection algorithm (also, 60 projection angles, 359 detectors). In fact, the direct integration method shows details which are (almost) not visible in the reconstruction obtained by the filtered back-projection algorithm. Both the contrast and the signal-to-noise ratio are improved by the proposed direct integration method. For the region of interest shown in Fig. 5 (S meaning the region of interest of the signal that we have randomly chosen, and N meaning the region of interest of the noise that we have randomly chosen), the contrast is 0.216 for the direct integration method and 0.196 for the filtered back-projection algorithm; the signal-to-noise ratio (SNR) is 3.315 for the direct integration method and 2.655 for the filtered-back projection algorithm. As a second example, consider the cross-section image shown in Fig. 6. In practice, this could be for example the cross-section of a metal tube inspected longitudinally. As shown in Fig. 7 the filtered back-projection method gives poor results, while the results obtained through the proposed direct integration method (Fig. 8) are significantly closer to the real image.

In both these examples we have considered the case of 60 projections, but better results are obtained also for more projections, up to 360 projections in both examples. However, for 60 projections, the difference between the proposed direct integration method and filtered back-projection is more clearly visible, this being the reason for choosing this case as illustration. More improvements in the reconstruction, in general, might be possibly obtained if, for example, the trapezoid method for approximation is replaced by other method, more suitable for the specific characteristics of the cross-section examined.

5 Discussion

We have proposed a method of tomographic reconstruction based on the direct integration of Radon’s inversion equation. We have compared the proposed method and shown that it provides for better reconstruction results, not only in terms of details, but also in terms of contrast and signal-to-noise ratio. In practice, the direct integration method, as described, is useful for example:

instead of the Filtered Back-Projection method, for any number of projections up to, say, 360 projections;