Abstract
We study a new MART-AP algorithm of few-views discrete tomography. Its efficiency for high-frequency structure reproduction is investigated in a numerical experiment where we reconstruct a 2D model for the estimation of the spatial resolution limit. We estimate the modulation transfer function of the reconstruction algorithm and compare it with the modulation transfer function of projection distortions. Our results show that MART-AP weakly influences the contrast of spatial structures being reproduced and can be used for high-resolution reconstruction when only a few projections are registered.
1. Introduction
We have recently seen a growing interest to few-views computed tomography [1]. A need to reconstruct the inner structure of an object in conditions when a few projections (usually less than 10) are only registered arises in many areas of tomography applications. These are, for example, nondestructive testing in industry [2, 3], luggage inspection [3, 4], plasma emission tomography [5], tomography of explosion-compressed metal shells [6], fast detonation research [7], and other. The main problem of few-views tomography is that reconstructed images are not free from artifacts which result from strongly limited data and make the reproduction of spatial structures less accurate. It is not easy to compensate the artifacts and as a rule, the problem is resolved through the use of a priori information about the object [8]. In doing so, the leading role is given to methods of discrete tomography [9–11], which use a priori information on the discrete values of the object function to be reconstructed. What is of great importance is the way in which this information is incorporated into the image. Until recently, the Bayesian reconstruction techniques [8–12] were thought to be most effective in terms of the ease of a priori information incorporation. Now more often different modifications to the well-known algebraic reconstruction technique [1] are applied to solve this problem. As an example, we can refer to an algorithm [13] based on a multistep procedure where each iteration includes segmentation with a priori information and then correction of intensities in cells which are near the boundaries of the contrast structures. In our recent work [14], we proposed a discrete tomography algorithm based on the multiplicative algebraic reconstruction technique (MART). We call it MART-AP (an MART which uses a priori information). Our algorithm, like that one described in [13], is based on a multistep solution correction procedure. But we apply corrections multiplicatively to all cells in the reconstruction region with an account for the a priori known discrete values of the object function. We incorporate this a priori information in the tomogram with the help of an image mask which is constructed through segmentation and adjusted at each iteration. The MART-AP algorithm has demonstrated a satisfactory performance in the reconstruction of simple binary images with low-frequency structures. Here we investigate the efficiency of the MART-AP algorithm for the reconstruction of high-frequency contrast structures and evaluate the spatial resolution of few-views discrete tomography with MART-AP. For this end, we use modulation-transfer-function analysis. Restricting ourselves to numerical experiment, we, however, compare the modulation transfer function (MTF) of MART-AP with the MTF of projection distortions present in real physical experiment due to blurring. The MTF of projection distortions on tomograms is estimated with an analytical approach known in optics [15] and radiography [16], but we adapted it to tomography.
The paper is organized as follows. Section 2 provides basic relationships which we use to estimate the MTF of projection distortions for a hypothetic X-ray facility which allows objects to be imaged from different angles. As an example, we provide here estimates for the MTF of projection distortions and the spatial resolution of tomograms obtained with the hypothetic facility, but without considering the effect of the reconstruction algorithm. Section 3 briefly describes the MART-AP algorithm and a numerical experiment for the reconstruction of a test phantom; reconstruction results are also presented. Section 4 provides results obtained in MTF analysis of data from the numerical experiment. It is shown how MART-AP can influence the accuracy of reproduction of high-frequency contrast structures on tomograms. We conclude with inferences on MART-AP applicability for the reconstruction of high-frequency contrast structures.
2. MTF of Projection Distortions
Both in radiography [16] and in tomography [17, 18], blurring is caused, first of all, by a finite focal spot of the source and a finite aperture of the detector. What also controls projecting accuracy is the interval between samples on a projection due to the sampling process [15]. Blurring can be estimated with the MTF which characterizes the contrast with which spatial frequencies are reproduced on the image. We can obtain an expression for the MTF if we use analytical formulas relating the MTF to other characteristics of blurring, such as the edge spread function, the line spread function, and the point spread function (see, e.g., [15, 16]). For tomographic images, it is suitable to use the notions of ideal forward and back projecting [1, 5].
With no loss of generality, we can restrict our consideration to 2D tomography where 2D images are reconstructed from 1D projections. Consider a point mass at point as an object. Here denotes the Dirac delta function. An ideal projection of such object can easily be obtained through applying, for instance, the singular Radon transform [5] to it: Here, stands for counts on the projection, and is an angle at which the projection is registered. Thus, the result of the ideal forward projecting of a pulse signal at point is also a point mass at point on the projection. Let the point mass be uniformly distributed within an interval of size after detection with a detector aperture of size . Then (1) transforms to
Take, for an example, an Abel inversion [5] for ideal back projecting, then for the isotropic point spread function of the tomographic image we obtain where is radial distance. In order to change to the frequency domain and estimate the MTF, we need to apply a 2D Fourier transform to , which is converted in the case of axial symmetry into the Hankel transform [19]. So, for the MTF which allows for the effect of the detector aperture, we obtain
where is the radial spatial frequency, and is the zero-order Bessel function of the first kind. Equation (4) holds for parallel ray geometry. Since fan (or conic) geometry is most often used in practice, a correction which allows for the geometric magnification of an object in projecting should be introduced in (4): where is the geometrical magnification factor.
Proceeding in the same manner, we can also obtain relations of type (5) for the MTF of the focal spot, , and the MTF of sampling, . In the end, the MTF of projection distortions which considers all the above three factors takes the form where is the focal spot size and is the sampling interval. Just this expression we are using here to estimate potential blurring for a hypothetic tomographic facility. Figure 1 exemplifies an MTF for such a facility. For the detector aperture we took mm which approximately corresponds to the intrinsic resolution of a cassette registration system based on photochromic screens with storage. Such systems are rather widely used in computed radiography [20], including its applications in industry [21]. For the focal spot we took mm. Such a focal spot corresponds to the fine-focus betatron [18], which is a high energy X-ray source and allows images of large and dense objects to be taken. As for the sampling interval, it seems appropriate to take it equal to the detector aperture size: mm. For the magnification factor we chose . This value is typical of many X-ray facilities.
The limit of spatial resolution was estimated with the Foucault-Rayleigh criterion [19]. As per the criterion, the ultimately reproducible spatial frequency corresponds to a point where the modulation transfer function meets a straight line characterizing 20% contrast (Figure 1). Analysis of Figure 1 shows that for the chosen parameters of the hypothetic X-ray facility, the ultimately reproducible frequency is 1.19 cycles/mm, and the associated limit of resolution is 0.42 mm. Later on we will see how the limit of resolution changes if not only blurring, but also the MART-AP algorithm checks the accuracy of reconstruction.
3. Reconstruction of a Resolution Test with MART-AP
3.1. MART-AP Description
It is well known (see, e.g., [1]) that the reconstruction problem can be described with the system of linear algebraic equations: where is a set of projection data ( is the number of projections), is a matrix of weight coefficients (weights), and is a vector of discrete values of the sought object function ( is the number of cells in the reconstruction region). In X-ray tomography where the object function stands for the distribution of linear attenuation coefficients, the trajectory along which particles move from source to receiver is assumed to be an infinitely narrow straight line (ray). Then the weight is the length of the intersection of the th ray with the th cell.
The problem of discrete tomography assumes that in addition to the measurement data , we also know a set of discrete values such that each th value a priori corresponds to a group of cells united in the th cluster. In monochromatic tomography, each value in will unambiguously correspond to one of materials constituting the object. As a result, the reconstruction problem reduces to search for such a distribution of values over the reconstruction region that would provide a best fit to the measurement data , that satisfy (7). In polychromatic tomography, we do not know the values exactly as the object function depends on photon energy. However, we can assume that for each material, it is possible to define a region of known values , where is an average value of the region and is a limit deviation. Thus, the above-stated problem of discrete tomography can also be stated and solved in the case of a polychromatic source. The difference is that the average value is used instead of .
The main idea of MART-AP is as follows. It organizes a cycle of external iterations where the traditional MART works. When MART reconstruction is finished, we create an adjuvant image which is referred to as a mask. This mask takes into account the result of MART reconstruction and a priori knowledge of discrete values of the object function. The mask is used at the next iteration of MART approximation.
The mask is a result of a transform , where is a set of natural numbers from the interval . At each external iteration the mask’s cells take the values from the set . In the interval these values correspond to the indices of object materials, which were found at the previous iteration. Zero values of the mask’s cells mean that none of the indices has yet been found to correspond to these cells. Iterations start with zeros for all cells of the mask. As iterations progress, the mask is being filled with nonzero values, while the current approximation is taking values closer and closer to the a priori ones. When an overwhelming majority of mask values become nonzero, the reconstruction ends and its result will be maximally close to the sought distribution of values over the cells in the reconstruction region. The basic steps of the MART-AP algorithm are presented in Algorithm 1.
| ||||||||||||
Corrections into the cells at each internal iteration of MART are applied in the standard manner [1]. However, if a cell is a member of a cluster for which the mask is nonzero, the similar corrections will also be applied to the other cells of this cluster. The correction procedure at internal th iteration can be defined by where is a parameter which regulates the rate of iteration convergence, is the sum of weights of all cells in the th cluster along the th trajectory, is the index of the cell in the th cluster, and is the number of cells in this cluster. Equation (8) is the basic formula of MART-AP. Equation (9) means that in order to correct the weight and then apply (8), it is necessary to calculate weight sums along trajectories within each cluster. Finally, (10) dictates that similar corrections need to be applied to all cells in the associated cluster and this is to be done for all clusters with nonzero values of the mask.
3.2. Numerical Experiment on the Reconstruction of the Resolution Test
Efficiency of MART-AP for high-frequency structure reproduction was investigated in a numerical experiment where we reconstructed a 2D model shown in Figure 2.
The model has two circular iron shells with the centers shifted from each other. The outer and inner diameters of the outer shell are equal to 100 and 80 mm. The same diameters of the inner shell are equal to 60 and 40 mm, respectively. The shift between the shell centers is 10 mm. Thus, the minimal gap between the shells is zero. Projection modeling and image reconstruction were done using fan beam geometry with the following parameters. The projections are determined in a 180° sector at equal angular steps. The number of projections is varied and equal to 15, 9, and 5, respectively. The geometrical magnification factor is equal to 1.5. The sampling interval is equal to 0.75 mm. The registration system is a line of 240 point detectors. Reconstruction was done onto a grid. For a priori discrete values, we took and . They characterized air and iron linear attenuation coefficients, respectively. In this particular case, the transform is written as , and the mask at the -th iteration subject to is where is a parameter which defines a neighborhood of point , inside which the object function must take the value equal to . In our numerical experiment, was varied within the range . The variation of helped us control the rate of iteration convergence. The reconstruction results are given in Figure 3. The left column of images presents the results obtained with MART (the first iteration of MART-AP). And the right column demonstrates what was obtained after 100 iterations of MART-AP.
(a)
(b)
Visual analysis of Figure 3 shows that the images of the right column are very close to the original of Figure 2. Only the tomogram reconstructed from 5 projections (the right image of the bottom row) seems to be of interior quality. It has visible distortions in the region where the shells touch each other (see the rectangle in Figure 2).
4. Analysis of the Numerical Experiment
The resolution test of Figure 2 is very suitable for MTF estimation as it allows a continuous scale of spatial frequencies to be determined. Indeed, the radial profile of the image taken at the angle (see Figure 2) can be put in conformity with the gap between the shells and, therefore, with some spatial frequency. On the basis of the profile of the reconstructed image, the modulation transfer coefficient can be determined as the relative depth of a trough between two intensity peaks. The MTF is estimated as a dependence of the modulation transfer coefficient on spatial frequency. The graduated formula defining a dependence of spatial frequency on angle is derived from the geometry of Figure 2 and is written as where is the inner radius of the outer shell, is the outer radius of the inner shell, and is the shift between the shell centers.
Figure 4 demonstrates the results of MTF estimation. The MTFs of the MART-AP algorithm, which were estimated from the right images of Figure 3, are shown by dotted lines. The solid lines correspond to the total MTFs which were calculated by the multiplication of the MTF of projection distortions (see Figure 1) and the MTFs of the reconstruction algorithm. We can see that MART-AP weakly influences the contrast of the spatial structures which can be reproduced by a hypothetic tomographic facility. The quantitative analysis of the solid lines in Figure 4 shows that, with contrast no less than 20%, the spatial frequencies of no more than 1.19, 1.03, and 1.01 cycles/mm can be reproduced on tomograms reconstructed from 15, 9, and 5 projections, respectively. It means that the spatial structures with linear sizes of 0.42, 0.49, and 0.5 mm can be resolved on reconstructed images. These estimates are very close to the resolution limit calculated from the MTF of projection distortions. Thus, the MART-AP algorithm seems to be quite effective for the reconstruction of high-frequency contrast structures. However, since the right bottom tomogram of Figure 3 has visible distortions, we can assume that 5 projections are not enough to reconstruct with MART-AP the image of very good quality.
(a)
(b)
(c)
In order to estimate the limit for the number of views, we additionally reconstructed the model of Figure 2 from 6, 7, and 8 projections and calculated the image error in the region of expected distortions (see Figure 2) by the following formula: where vectors and present the model and the reconstructed image, respectively. Sum in (13) is taken over all cells of the region of expected distortions, which is shown in Figure 2 by the rectangle. Figure 5 shows the fragments of the subtraction images , which visualize the region of distortions. And Figure 6 plots the dependence of the image error on the number of projections. We can see that this dependence has a dip when the number of projections changes from 6 to 7. Then the image error varies weakly if the number of projections grows from 7 to 15. It does mean that the minimal number of projections, which can be recommended for MART-AP applications, is equal to 7.
Our codes are written in MATLAB. Calculations were done on a 1.7 GHz Pentium 4, 512-Mb RAM Intel PC. The reconstruction of an image onto a grid with MART-AP (100 iterations of the traditional MART) takes approximately 20 minutes. No optimization was used to make calculations faster.
5. Conclusion
We have studied a new iterative reconstruction algorithm for few-views discrete tomography. The algorithm, we called MART-AP, is based on the known multiplicative algebraic reconstruction technique and incorporates into the image a priori information on the discrete values of the object function to be reconstructed. To check the efficiency of MART-AP for high-frequency structure reproduction, we have conducted a numerical experiment on the reconstruction of the resolution test. We have estimated the modulation transfer function of the MART-AP algorithm for a different number of projections and compared it with the modulation transfer function of projection distortions, which was calculated for a hypothetic tomographic facility.
Our results show that the MART-AP algorithm is quite effective for the reconstruction of high-frequency contrast structures and theoretically allows submillimeter spatial resolution to be obtained. The minimal number of projections, which we can recommend for MART-AP applications, is equal to 7.
In the future, it is interesting to consider more complicated tests for several materials which differ in image contrast. Of course, it is of special interest to study the algorithm in a physical experiment.
Acknowledgment
The authors thank S. V. Kolchugin for fruitful and helpful discussions during the paper preparation.