Applying Mojette discrete Radon transforms to classical tomographic data H Fayad 1,2, JP Guédon 2, I Svalbe 3, Y Bizais 1 , N Normand 2.. 1. LATIM, INSERM U650, Brest, Bât. 2bis (I3S) CHU Morvan - 5, Av. Foch, 29609 Brest - France 2. IRCCyN CNRS 6597, Ecole polytechnique University of Nantes, La chantrerie, BP 50609, F44306 Nantes France. 3. School of Physics, Clayton Campus, Monash University, 3800 Melbourne, Australia. ABSTRACT Tomographic acquisition uses projection angles evenly distributed around 2π. The Mojette transform and the discrete Finite Radon Transform (FRT) both use discrete geometry to overcome the ill-posedeness of the inverse Radon transform. This paper focuses on the transformation of acquired tomographic projections into suitable discrete projection forms. Discrete Mojette and FRT algorithms can then be used for image reconstruction. The impact of physical acquisition parameters (which produce uncertainties in the detected projection data) is also analysed to determine the possible useful interpolations according to the choice of angle acquisitions and the null space of the transform. The mean square error (MSE) reconstruction results obtained for data from analytical phantoms consistently shows the superiority of these discrete approaches when compared to the classical “continuous space” FBP reconstruction. Keywords: projection acquisition, tomographic reconstruction, Mojette transform, Finite Radon Transform
1. INTRODUCTION Classical tomography starts with a set of projections of the object under reconstruction. In Radon’s initial paper [Radon 1917], the definition of the transform was entirely continuous, both for the 2D function under projection and reconstruction, and for the projections themselves. These projections are generally equally spread around π or 2π (2π is not needed mathematically but because of the physics of acquisition), and each projection is sampled with a fixed number (Nbin) of bins. For an NxN image grid, the tomographic process inverts the linear system p=Rf, where p represents the set of projection values, f the image grid and R stands for the Radon transform. Classically, the number of projections is chosen to be close to N, in order to get the same range of unknowns and data. However, this unitary transform (compared, for example, to Fourier, Hadamard, wavelet, Haar) can not be inverted in a simple manner because it belongs to the class of ill-posed problems. The ill posedness comes from the sampling of the transform. Figure 1 exemplifies that situation. The Algebraic Reconstruction Technique of Herman [Herman 79] focused on the discrete characteristics of the resultant image. It gives a remarkable iterative algorithm to find a consistent solution for examples such as the configuration shown in figure 1.a (even where there are not enough data for perfect reconstruction). The operator R*R (projecting – backprojecting data) ensures a consistent solution is obtained after each iteration. The focus of Katz’ work [Katz 78] was to associate the discrete grid of the reconstructed image with the projections of the object, in order to obtain a sufficient number of samples that are determined by the set of discrete angles selected for projection. Discrete tomography arose from the field of discrete geometry, albeit that Myron Katz’ pioneering work in [6] was only focused on ill-posedness of the Radon transform. Initially, the focus was to reconstruct a binary object from a pair of (orthogonal) projections using various algorithms, or to decide, from the available data, if there were none, one or many
.
Further author information. Send correspondence to
[email protected], -
[email protected] [email protected]
Medical Imaging 2008: Physics of Medical Imaging, edited by Jiang Hsieh, Ehsan Samei, Proc. of SPIE Vol. 6913, 69132S, (2008) · 1605-7422/08/$18 · doi: 10.1117/12.770478
Proc. of SPIE Vol. 6913 69132S-1 2008 SPIE Digital Library -- Subscriber Archive Copy
possible solutions [4,5]. The Mojette transform [7,8] and the Discrete Radon Transform (FRT) also provide solutions to exactly this kind of problem, but for more general image values and more flexible projection constraints.
a b+d c
a
b
a+b
a
b
a+b
d
c
c+d
d
c
c+d
Figure 1 : A 2x2 image with two set of projections, 1.a : the constant density of sampling on each projection does not give a unique solution; 1.b the adaptive angle sampling of the Mojette transform gives a unique solution through use of a simple linear matrix inversion.
There exists a natural antagonism between the null space and the classical continuous Radon transform, compared to the same situation that pertains for discrete spaces. The mechanical and physical constraints for the acquisition of real projection data also may compromise the use of the “exact” angles that are suited to discrete transforms. Servières, Normand and Guédon [1,2] have shown how to use the Mojette transform for tomographic reconstruction from data already in the Mojette space. In order to do that, discrete angles, denoted (p, q) where p, q ∈ Z and gcd(p, q) = 1, are taken, such that the angle tan(θ) = q/p. There are discrete geometry conditions (the Katz criterion) that permit exact reconstruction of the original image from such projections for noise-free data. There are also different algorithms, based either on Filtered BackProjection (FBP) or on iterative algorithms, such the conjugate gradient method, that have been derived to used Mojette approach for noisy projections, even in the case where not enough data is available. During the same period of time, Svalbe and Kingston have derived efficient algorithms to work with the Finite Radon Transform, which uses the same discrete angles, albeit with a different bin sampling and pixel filling methods [10, 11]. The next section is devoted to describing the essence of the Finite Radon Transform and the Mojette transform, in terms of projections, sampling and main properties of both transforms. Section 3 presents the method of transforming classical projections into discrete projections. The results shown in Section 4 compare an original image, its FBP reconstruction from discrete projections, and its FBP reconstruction using the approximated projections as computed using some of the methods presented in section 3.
2. DISCRETE ANGLES FOR THE FINITE RADON AND MOJETTE TRANSFORMS 2.1 Discrete angles A classical tomograph acquires, say, M projections, equally distributed over 2π, with a fixed number Nbin of bins onto each projection. Then, for standard direct or iterative algorithms, the linear problem is inverted using many different kinds of physical constraints. A first rebinning step is mandatory to convert the acquired projection into a Mojette or DRT one. Figure 2 presents these two kinds of projections. Here we present simple method to generate discrete projections from classical ones. A Farey series of order N, denoted FN, is the set of ordered ratios q/p with 0 ≤ p ≤ q ≤ N, where p, q ∈ Z and gcd(p, q) = 1, such that the angle θ = atan(q/p) belongs to [0, π/4]. For instance, F3 = {0/0, 1/3, 1/2, 2/3, 1/1}. However, in this paper, we extend the set of angles by symmetry to cover the full [0, π] range. An classical tomograph indexes data by angles theta(i) = 2πi / I, i=0 to I-1; and onto the projection (or ray index) by another index k = 1 to Nbin (typically to obtain around 64 or 128 samples for SPECT and TEP, or 256 or 512 samples
Proc. of SPIE Vol. 6913 69132S-2
for MRI) to cover the whole (disk) field of acquisition. The method presented here starts with the choice of the image size N×N for the reconstruction grid, the set of discrete angles (p, q) and an angle theta(i) that is closest to (p, q). Notice that, in accord with the physical constraints explained in section 3, the closest angle is not always the only (or best) choice for performing the interpolation. This can be easily seen in figure 2 where it is evident that special effects occur at different angles, especially around 0 and 45°.
21
19
17
15
23
21
19
17
15
13
9
11
7
5
0 3
5 6 7 8 9 10 11 12 13 14 15 16 17
1
1 2 3 4
13
10
0
9
20
10
11
30
20
7
30
8 7 6 5 4 3 2 1 0 5
40
3
50
40
1
50
Figure 2 : Comparison for (0, 45°) range between regular angles increments (step of 2,8125 °) and Farey series (N=8) increments (2.c shows the differences between two successive angles)
2.2 Finite Radon Transform The Finite Radon Transform (FRT), as defined by Matus and Flusser [9], is built from components well suited to discrete projection geometry: integer (whole pixel) translates for projection rays and rational fractions (from the Farey sequence) to define the discrete projection angles. The problem of inverting discrete FRT projections is tackled up front by restricting (for the 2D case) the reconstruction image array not just to be square, N×N, but requiring N to be prime (eg 257×257, rather than the traditional 256×256). The other major assumption made is to impose periodic boundary conditions. A discrete FRT projection ray samples whole pixel values located at separation intervals p pixels across and q pixels down. When the ray reaches any border of the image space, the sampling ray wraps to the opposite side of the image and continues sampling values at the same p, q steps but with an offset reset by the modulus N operation. The prime N array size requirement means all discrete projection rays cross each other once and only once (as the gradient of any ray is never a multiple of the gradient of any other ray). Similarly, by insisting on taking exactly N integer translates for the parallel rays that make up each projection (creating one ray through each pixel on each row, or column), all pixels in the image space are sampled once and only once by each projection. To formalise the above, the discrete FRT projection, R(t, m), of an N×N discrete image f(k, l), is given as: R(t, m) = Σf(t + mi, i)
for 0 ≤ i ≤ N and 0 ≤ t < N, and where N is prime.
(1)
The integer variable m represents the gradient of a discrete FRT projection. There is a one-to-one relationship between the value of m and the discrete projection angle, (p, q). The assignment of m to some angle (p, q) is fixed by finding the minimum distance between the pixel samples, as measured along each ray direction [10]. For example, m = 1 always corresponds to the (1, 1) or 45 degree projection. We assign m = 0 and N to be projections along the row and column directions. The m = N projection keeps the transform symmetric (it can be formulated using translates either along rows or columns). The FRT is an optimal projection system, as it stores N translates for each of (N + 1) projections to make N2 + N projection samples to represent N×N independent image pixel values. As each projection samples all image pixels once, the sum of each projection is equal to the sum of all image values (denoted S in the following). Then only N-1 translate values are needed for each projection provided the image sum S is known at least once, yielding (N - 1)(N + 1) + 1 = N2 independent projection samples. The FRT then automatically meets Katz’ criterion and a direct and exact reconstructed image can be retrieved from any complete set of FRT projections.
Proc. of SPIE Vol. 6913 69132S-3
The ease of inversion of the FRT is its best attribute (as it was achieved by deliberate design). Because each projection ray samples all pixels just once, there is no need to compensate for variable density sampling (as in traditional back projection) nor are there any problems of unpacking the mixed content of projection rays, as each ray is known to cross each other only once. Inversion can be performed either by doing FRT projection sums on the FRT projection space (the FRT forward and inverse is very close to symmetric) or by simple back projection of the FRT data into the image space. The latter method is very efficient, as each whole row from R(t, m) can be copied intact into the reconstructed version of f(k, l), with only periodic translation shifts (by m pixels) are required as the data from the mth projection is copied to each row.
t=3
Figure 3. Selection of the image pixels sampled for FRT projection R(t, m) at m =2, t = 3 for N = 7.
The apparent disadvantage of the FRT relates not to need to enforce use of a prime array size, but to its assumptions of image periodicity and the wrapping of projections. However, as shown in [10], the Mojette projection data can be summed into the periodic FRT form, without requiring any interpolation or rebinning of the data. The FRT is then a special case of the more general Mojette approach, as outlined in the next section. The other main objection to the use of the FRT (and Mojette) is that the gap between pixels sampled along a ray (a separation of √(p2 + q2) for ray angle (p, q)), is unphysical for real projection data. This is a serious point that needs to be addressed by the detailed pixel model used to represent the projected and reconstructed data. For the FRT, the largest gap distance between pixel samples can be shown to scale as √N, so that there is a relative convergence to the continuous case as N becomes large (this justifies the mapping of m to the shortest distance between neighbouring (p, q) pixels). 2.3 Mojette Transform The Mojette transform is a discrete, exact version of the Radon transform, just as is the FRT. From a digital image f(k, l), a set of discrete projections is computed, but the orientation and the number of bins onto each projection varies: the sampling path onto projection (p, q) is given by: t2 (p2+q2) =b2 and the number of bins for a discrete projection is given by: NB(p, q, N) = (N - 1)(IpI + IqI) + 1.
(2)
As emphasised in the previous subsection, this sampling leads (for N = 512) to NB(1, 1, 512) = 1023 and NB(30, 17, 512) = 24018, compared to the typical classical values for Nbin in tomography of around 512. The simplest Dirac Mojette definition corresponding to the FRT (in terms of the grid sampling modeled by a 2D Dirac function inside a pixel) is:
Proc. of SPIE Vol. 6913 69132S-4
proj(b,p,q) = M f(k,l) = Σ
Σ f(k,l) ∆(b + qk – pl) = Mf(k, l), with ∆ (b) = 1 if b=0 and 0 elsewhere.
(3)
This discrete geometry can be used efficiently for tomographic reconstruction, as reported in work by Servières et al [7, 8]. Different algorithms can be used, even in the presence of noise, to invert the Mojette transform. However, this work assumes we start with data acquired exactly in the Mojette domain. The topic of this paper is to discuss the optimal mapping from classical tomographic data, as acquired from real x-ray projections, into a form compatible with the Mojette and fast Radon domains. 2.4 Relationship between Mojette and Finite Radon Transform The Mojette transform can be inverted using a variety of methods [1], and even via the same algorithm as for the Finite Radon transform. It has been independently described in the Ph.D Theses of M. Servières (respectively A. Kingston) with the Mojette operator (resp. with the FRT) as reported in [1,10]. In both cases, the entire Farey series of projection angles is directly back projected onto image g(k, l) (one bin per projection per pixel), from which the sum of bins S is subtracted, and the result divided by the number of projections minus one: (I – 1) f(k,l) = g(k,l) – S.
(4)
The backprojected image g(k, l) is obtained directly either from the relationship M*M f =g, using all of the relevant Farey angles, or by inverting the finite Radon data after also using the backprojection operator.
3. METHOD 3.1 Bin values The method presented here is simple, since it reduces the problem to how best to reformat available data, from some given projection angle, into a form suitable for one at a predefined Farey angle. In this case, the projection data is only resampled from the Nbin original bins to the correct (and generally much bigger) number of the nearest Mojette projection bins, using the generalized sampling theorem of Unser and AlDroubi. The discrete tj values are first converted into a continuous signal using a cardinal interpolator kernel ck using the classical formulae: proj(s) = Σ proj(tj) ck(tj - s) .
(5)
This continuous signal can now be sampled at any rate. Of course, as written in (5), the mathematical properties of the cardinal interpolator are transmitted to the proj(s) data. For these experiments, a simple linear interpolator was chosen.
Proc. of SPIE Vol. 6913 69132S-5
tj b
(p,q)
teta
Figure 3: Classical and Mojette projections correspondence
An equivalent method performs the re-sampling into the FRT domain, where each interpolated projection, in Mojette form, is re-summed from Nbin to match the number of translates (of prime size) that are selected for the FRT data [11]. The measure of discrepancy of the initial data interpolation can only be measured in the image original space because of the differences in the projections definitions. Thus, the results presented here are only done with the two sets of projections: original and Mojette. Reconstructions that use the FRT algorithm (after Mojette re-sampling of the projection data as used here) have yet to be completed and will be reported elsewhere. 3.2 Choice of the near discrete angle A Farey series order N has to be chosen. If N is small, many real angles will map to the same corresponding Farey discrete angle. This is shown by the number of non-similar discrete angles corresponding to 360 real angles for different Farey series: for N = 5, 40 discrete angles spread the [0, π] range, but N = 20 corresponds to 311 angles (see also the work in [3]). On the other hand, if a large N is chosen, a huge number of bins must be computed. Our algorithm for choosing a discrete angle starts with a small Farey order (N = 10) and, on finding a discrete angle already allocated to a previous theta candidate, the order of the Farey series is increased by one.
4. RESULTS AND DISCUSSION 4.1 Results This method allows for the reconstruction of image data using any discrete (FRT or Mojette) method. In order to perform a fair comparison, we use almost the same number of projections M. To permit the computation of the mean square error (MSE) distance with respect to an original image, analytical phantom data was used. The first phantom is a simple square inside a 129 x 129 image having four intensities: a background of 0, a central square of 0.5, corners of 0.25 and the inside square values of 1. The classical Shepp-Logan phantom was also used. The reconstruction method was the Filtered Back Projection, in order to perform a fair comparison. The classical FBP and its ideal Mojette counterpart were implemented as described in [8]. The proposed method was to use only those classical FBP angles that were remapped onto discrete angles (p, q) from the relevant Farey series, before performing the Mojette FBP. Thus we have computed an error measure taken with respect to the original image for three cases.
Proc. of SPIE Vol. 6913 69132S-6
In order to compare the Ideal method (using exact discrete angles) and the Proposed method (approximate discrete angles before the reconstruction), the distances between projections generated by each method were calculated. Here smaller error distances indicate a more faithful reconstruction. The results are shown in the following table: Number of projections 30
90
180
360
Total difference
689.6
2091.3
5454.3
13997.85
number of Bins
41758
119514
297524
825704
Difference per bin
6.28E-4
3.82E-4
2.48E-4
1.43E-4
Table 1 : Distance onto the projections
N13O
N19O
N118O
N136O
MSE: 5.
]\ISE: U. U U42
T1SE: 1.
9SE— 4
MSE: 1.
USE— 4
T1SE: 1.
U7E— 4
MSE: 2.
42E—5
MSE: 7.
92E—6
MSE: 2. 24E—6
MSE: 1.
14E— 4
MSE 2. 97E—5
MSE: 1.
U7E—5
]\ISE: S.
62E—5
7—6
Figure 4. MSE Comparison for square reconstruction. First line: classical FBP reconstructions. Second line Discrete Mojette reconstruction algorithm with exact Mojette projections. Third line: same method as the second line, but after rebinning classical projections onto discrete ones.
Proc. of SPIE Vol. 6913 69132S-7
LJ
bJ
P C N
H
P C
H LJ
bJ
co
P C
C
P C
bJ bJ
H
bJ
C
Q
C
P
LI
a
C
The first (unsurprising) conclusion is the same for the three methods: the error decreases in proportion to the increasing number of projections. For the Square phantom reconstruction, the proposed method gives results better than the FBP (MSE difference with a factor of 10) and is only a little worse than the ideal method. The reconstruction error images are compiled in Figure 4, together with the MSE, as taken with respect to the original image.
Figure 5. MSE Comparison for Shepp-Logan reconstruction: First line, classical FBP reconstructions. Second line, Discrete Mojette transform algorithm with exact Mojette projections. Third line, proposed discrete method after rebinning classical projections into discrete ones.
This factor of 10 (which differs only a little according to the number of projections) between the mean square error in the proposed method and the mean square error of the FBP, is retained for the Shepp-Logan phantom reconstructions, as shown on Figure 5. In addition, the difference between the mean square error in the Proposed method and in the Ideal method still remains very small.
Proc. of SPIE Vol. 6913 69132S-8
4.2 Comparison between Mojette FBP and the completely discrete algorithm The exact discrete tomographic back-projection can also be implemented as given by (4). However, this reconstruction will have a drawback in terms of computational efficiency (because of the large number of angles) and in sensitivity to noise, which degrades the reconstructions. This is shown by using the same set of angles generated by our algorithm, and then using (4) to perform the reconstruction.
Figure 7. Reconstructions using (4) and using the angle interpolation part of the Proposed method applied to reconstruct: (a) Square phantom. (b) Shifted square phantom. (c) Shepp-Logan phantom.
Proc. of SPIE Vol. 6913 69132S-9
As a result, we can make a first comparison to the case where only some fraction of the angles within a given Farey series are used. The results are shown in Figure 7. As exemplified by the poor results, a reduced angle interpolation scheme is not suitable when using this method. 4.3 Discussion The results presented here show that the use of discrete angles, allied with a suitable discrete reconstruction algorithm, always performs better than FBP reconstruction for a classical, evenly-spaced angle set of acquired projections. It also demonstrates that the approach proposed here always reduces the MSE distance compared to the classical FBP. Moreover, the distance between original discrete angle projections and the rebinned projections can always be kept small. There are hence two major and novel outcomes resulting from this study: 1) Any set of real, acquired tomographic data can be rebinned into a compatible Mojette or DRT projection space, without any loss of reconstruction power. 2) The distance between an original phantom and the reconstructed image after classical FBP from projections is always greater than the distance between original and discrete reconstructed images using the same data set. It should also be stressed that the interpolation scheme was crude and could have been refined by using the projection of a 2D interpolator onto the projection instead of a simple 1D interpolator (not connected to the image properties). This paper improves the prospect for the increased use of discrete reconstruction techniques to be applied to real projection data, and flags that this approach may be more than competitive by providing results superior to those obtained from classical iterative or backprojection algorithms.
5. CONCLUSION This paper demonstrates a practical interest in revisiting the use of classical FBP via the application of discrete Radon transforms (such as the Mojette or DRT) for tomographic image reconstruction for real data sets. The results presented here (at least as assessed just through the error distance measure MSE) clearly show the advantage of using these discrete reconstruction methods, and opens the door to further comparison with iterative classical methods. The next step will be to simulate the effect that blur from a real, finite detector device will have on the projection values. It will then be interesting to see if the advantages, as demonstrated clearly above, will persist with noisy data.
Proc. of SPIE Vol. 6913 69132S-10
REFERENCES [1] M. Servières, N. Normand, Jp. Guédon, Y Bizais «The mojette transform : discrete angles for tomography » Electronic Notes in Discrete Mathematics, Elsevier, pp. 587 – 606, July 2005. [2] M. Servières, N Normand, P. Subirats, JP Guédon, «Some links between continuous and discrete Radon transform » SPIE Medical Imaging 2004, Feb.2004. [3] I. Svalbe, S. Chandra, A. Kingston and JP Guédon «Quantised Angular Momentum Vectors and Projection Angle Distributions for Discrete Radon Transforms» , LNCS 4245, Discrete Geometry and Computer Imagery, Kuba ed., Oct. 2006. [4] R. Gordon, R. Bender, G. Herman, «Algebraic reconstruction techniques (ART) for three-dimensional electron microscopy and x-ray photography. » Journal of Theoretical Biology 29:471-482, 1970. [5] R. Gardner, Geometric tomography, Cambridge University Press, 1995. [6] M. Katz, Questions of Uniqueness and Resolution in Reconstruction from projections, Lecture Notes in Biomathematics, Springer-Verlag, Vol. 26, 1979. [7] M Servières, J. Idier, N. Normand, and JP Guédon. « Conjugate gradient Mojette reconstruction. » In SPIE Medical Imaging 2005 : Image Processing, volume 5747, pages 2067–2074. J. Michael Fitzpatrick and Joseph M. Reinhardt, 2005. [8] M. Servières, N. Normand, Y. Bizais, and JP Guédon. Noise behavior for spline mojette FBP reconstruction. In SPIE Medical Imaging 2005 : Image Processing, volume 5747, pages 2100–2109. J. Michael Fitzpatrick and Joseph M. Reinhardt, 2005. [9] F. Matus and J. Flusser, Image Representation via a Finite Radon Transform, IEEE Trans. PAMI, 15(10), 1993, 9961006. [10] A. Kingston and I. Svalbe, Projective Transforms on Periodic Discrete Image Arrays, Advances in Imaging and Electron Physics, 139, 2006, 75-177 [11] I. Svalbe and D. van der Spek, Reconstruction of Tomographic Images Using Analog Projections and the Digital Radon Transform, Linear Algebra and its Applications, 339, 2001, 125-145.
Proc. of SPIE Vol. 6913 69132S-11 View publication stats
