A Fast Reconstruction Algorithm For Electron Microscope Tomography

64K. Sandberg et al. / Journal of Structural Biology 144 (2003) 61–72specimen from its projections can be viewed as the inversion of the Radon transform. Many reconstructionalgorithms rely on this fact and solve the inversionproblem by discretizing the inverse Radon transform(Deans, 1993; Gilbert, 1972). In this section, we describethe widely used direct summation algorithm (alsoknown as the filtered or R-weighted back projection).We discretize (4) and obtainNhXgðxm ; zn Þ ¼wl ðq Rhl Þðtðxm ; zn ÞÞ;ð7Þl¼1where tðxm ; zn Þ is given by (1) and wl are weights to compensate for unequally spaced angles. For measurementsperformed over equally spaced angles, the weights wl areusually set to one. Since we have measurements only for adiscrete set of values of t, the elements ðq Rhl Þðtðxm ; zn ÞÞM 1are estimated from fðq Rhl Þ ðtk Þgk¼0 by some interpolation scheme, usually piecewise linear interpolation.Let us summarize the steps for estimating the densitygðx; zÞ from the measurements of projections.1. Filter the data to obtain ðq Rhl Þðtk Þ, k ¼ 0; 1; . . . ;M 1 by:1.1. applying the FFT along the columns of the matrix rkl defined by (3),1.2. multiplying each element of the transformed matrix by the (pre-computed) filter coefficients andthe weights wl if necessary, and1.3. applying the inverse FFT column-wise.2. Sum the result of the previous step to obtain the density by:2.1. computing tðxm ; zn Þ for each given ðxm ; zn Þ,2.2. finding ðq Rhl Þðtðxm ; zn ÞÞ by linearly interpolating ðq Rhl Þðtk Þ, and2.3. summing the result according to (7).Step 2 dominates the computational cost since we haveto sum over Nh terms N M times, for the total computational cost of OðNh MN Þ. Usually Nh , M, and N areof the same order of magnitude so the above algorithmhas the computational cost OðN 3 Þ.where!2 pxsin coswl 2pixs cosxhxhlleq vl ðxÞ ¼pxcos hlcos hlcos hl M 1XRhl ðtk Þ e2pik cosxhlð8Þ:k¼0In Section 4.2, we will describe a numerical implementation of this formula that results in a fast OðNh Mlog MÞ þ OðMN log N Þ algorithm. We note that shiftingthe coordinate system in x by a constant xsmay be necessary to account for the deformation of aspecimen during the data collection.In order to obtain (8), we will discretize the z-variablebut keep x as a continuous variable until the very end ofour derivation. Let fzn gNn¼1 be an equally spaced grid inthe z-variable. If we fix z ¼ zn while treating x as acontinuous variable, then we write gn ðxÞ ¼ gðx; zn Þ andtn ðxÞ ¼ tðx; zn Þ, whereR 1t is defined in (1). In the followingderivation f ðxÞ ¼ 1 f ðtÞ e2pitx dt denotes the Fouriertransform of a function f ðtÞ.By introducing the notation fl ðtn ðxÞÞ ¼ wl ðq Rhl Þðtðx; zn ÞÞ, we write the sum for the direct summation (7) asgn ðxÞ ¼NhXfl ðtn ðxÞÞ;ð9Þl¼1where x is a continuous variable. To derive (8), we applythe Fourier transform with respect to x to both sidesof (9) and perform the summation over angles inthe Fourier domain. The Fourier transform of (9) withrespect to x givesNh Z 1Xg n ðxÞ ¼fl ðx cos hl þ zn sin hl Þ e2pixx dxl¼1¼NhXl¼1¼NhX 1e 2pixzn tan hlcos hlZ1fl ðsÞ e2piscosxhlds 1vl ðxÞ e 2pinl ðxÞzn ;ð10Þl¼13. Inversion formula in the Fourier domainIn this section, we derive an inversion formula in theFourier domain that is equivalent to the direct summation formula in (7). We will show that ifxm ¼ xs þ m;m ¼ 1; 2; . . . ; Mf ;where xs is a shift parameter that depends on the selection of the coordinate system in the x-variable, thengðxm ; zn Þ can be written asNhXwl ðq Rhl Þðtðxm ; zn ÞÞl¼1¼Z12 12NhXl¼1!vl ðxÞ e 2pixzn tan hle 2pixxm dx;where nl ðxÞ ¼ x tan hl andZ 112pis xfl ðsÞ e cos hl ds:vl ðxÞ ¼cos hl 1ð11ÞBy definition, we havefl ðtn ðxÞÞ ¼ wl ðq Rhl Þðtðx; zn ÞÞZ 1¼ wlq ðxÞR hl ðxÞ e 2pixt dx;ð12Þ 1which combined with (11) gives us wlxx Rhlvl ðxÞ ¼q :cos hlcos hlcos hlð13ÞRecall that Rhl ðtk Þ corresponds to a discrete set ofmeasurements. In the direct summation algorithm, we

K. Sandberg et al. / Journal of Structural Biology 144 (2003) 61–72use interpolation to define Rhl ðtÞ for any t. To incorporate the linear interpolation, let us introduce the ‘‘hat’’function, or the linear spline 1 jtj; 1 t 1;ð14ÞbðtÞ ¼0;otherwise;with its Fourier transform given by 2sinðpxÞ bðxÞ ¼:pxWe express the piecewise linear interpolation of thediscrete data using bðtÞ by definingRhl ðtÞ ¼M 1Xbðt k þ xs Þrkl :ð15Þk¼0It is easily verified that the function Rhl ðtÞ is continuouswith respect to t. Using (15) we haveZ 1M 1X Rhl ðxÞ ¼rklbðt k þ xs Þ e2pitx dt 1k¼0¼M 1Xrkl e2piðk xs ÞxZ1bðsÞ e2pisx ds 1k¼0¼ e 2pixs x b ðxÞM 1Xrkl e2pikx :ð16Þk¼0¼ Fl ðxÞ rl x;cos hlvl ðxÞ ee 2pixx dx:ð21Þl¼14. Implementation4.1. DiscretizationLet us evaluate (21) at pixel locations in the finalimage gn ðxm Þ where xm is given byxm ¼ xs þ m;m ¼ 1; 2; . . . ; Mf ;and xs is a shift parameter due to possible deformationsof the specimen. In order to avoid aliasing artifacts inthe FFS algorithm, we construct an Mf N -imagewhere Mf M. The number of pixels Mf needed toavoid aliasing will be selected in Section 4.3. By discretizing x as1kxk ¼ þ;2 Mfk¼ 2k ¼ 1; 2; . . . ; Mf ;þ1Mf21 X¼MfMfk¼ ð17Þwhere nl ðxÞ ¼ x tan hl and vl ðxÞ are defined in (17).Let us consider a smooth bandlimited filter q with itsFourier transform q defined by (6). Since q is zero forjxj 12 we observe that from (17) it follows thatvl ðxÞ ¼ 0 for jxj 12. Hence (20) reduces to2ð22Þþ1wherekl¼1þ11 X 2pik Mmf;g nk eMfMf2pi xg nk ¼ e Mf sð19Þ2! NhXk 2pinl ðMk Þzn2pi k x 2pi k mfee Mf s e MfvlMfl¼1Mf2k¼ ð18ÞWe note that the factor Fl ðxÞ is independent of themeasured data once the angles hl are known.The final step is computing gn ðxÞ from g n ðxÞ. Bytaking the inverse Fourier transform of (10) we arrive at!Z 1 XNh 2pinl ðxÞzngn ðxÞ ¼vl ðxÞ eð20Þe 2pixx dx; 1! 2pinl ðxÞznThe equation (21) is equivalent to the sum (9) used in thedirect summation algorithm. We can obtain a fast reconstruction algorithm by computing the sums in (19)and (21) using the USFFT described in the Appendix A.¼and XM 1x2pik x rkl e cos hl :rl¼cos hlk¼0NhX 12where wl 2pixs cosxhxx lbFl ðxÞ ¼eq ;cos hlcos hlcos hl12we approximate (21) by using the trapezoidal rule!Mf Nh2X1 Xk 2pinl ðMk Þzn 2pi k xfe Mf megn ðxm Þ ’vlMfMfMfl¼1Combining (13) and (16) yields wl 2pixs cosxhxx lq bvl ðxÞ ¼ecos hlcos hlcos hlM 1X2pik x rkl e cos hlk¼0gn ðxÞ ¼Z65 NhXk 2pinl ðMk Þznfvl:eMfl¼1ð23ÞDue to the construction of the filter q in (6), the integrand in (21) can be considered to be smooth and periodic in x with period 1. Under this condition, thetrapezoidal rule is rapidly convergent and we canachieve any desired accuracy by choosing Mf large enough. We have found that for our applications settingMf to 1:5 2 times the number of projections M usuallysuffices. We will demonstrate the accuracy of (22) inSection 5.2.1.We note that the reconstruction formula (22) allowsthe grid in the z-variable to be shifted by a constant zs byscaling vl . Such shifts are essential in some cases whenreconstructing three-dimensional volumes as discussedin Section 5.1.

66K. Sandberg et al. / Journal of Structural Biology 144 (2003) 61–724.2. Numerical algorithmOur first goal in designing the FFS algorithm is tomatch the results with those of the direct summationalgorithm. We do it for two reasons. First, since thedirect summation algorithm has been used for a longtime in TEM and significant experience has been accumulated for interpretation of the images, we avoid theissue of acceptance. Second, we demonstrate the flexibility of the FFS algorithm. As it turns out by changingparameters we can achieve a higher order interpolationin the input data in comparison with the linear interpolation used within the direct summation.In order to match the direct summation algorithm,we consider linear interpolation and use a bandlimitedversion of the filter q ðxÞ ¼ jxj. The discretization of theradial weighting filter is given by8jxk j 6 xc ; jxk j; ðx xc Þ2k 2q ðxk Þ ¼ x e 2xs ; x jx j 6 1 2;ck : c0;jxk j 1 2;where xc is a user specified cutoff frequency and xs is auser specified parameter for the Gaussian roundoff. Thisfilter matches the filtering routinely used for direct summation, but other choices of filters are also possible. By anappropriate choice of the Gaussian roundoff, the filter is(approximately) continuous with respect to x. The discretization of the linear interpolation filter is given by 2sinðpxk Þb ðxk Þ ¼:pxkWe discuss other choices of this filter in Section 4.4below.Next we summarize the main steps of the FFS algorithm. To estimate the number of operations, let usconsider Nh projection angles with M samples each toreconstruct an image with M N pixels. We consider athree-dimensional volume consisting of Ny slices. In whatfollows Mf (to be estimated in Section 4.3 below) is thenumber of spatial frequency modes in the x-direction.The steps are illustrated in Fig. 2. Using the symmetry ofthe Fourier transform f ðxÞ of the real data f ðtÞ, wedouble the speed of the algorithm by using f ð xÞ ¼f ðxÞ.The algorithm assumes a limited range of projectionangles. For the applications in this paper the angles aretypically between 70 . We note that the algorithmcan be modified for the full angular range by firstcomputing the density gðx; zÞ based on the projectionangles satisfying jhl j 6 45 , then interchanging the roleof x and z, and computing g based on the projectionangles satisfying jhl j 45 . Finally, we add the tworeconstructions. A similar procedure was used by Pottsand Steidl (2001).In most applications, the angles hl and filters are thesame for all slices, which means that the precomputationstep in the algorithm only needs to be done once whilesteps 2.1–2.3 are repeated for each slice.4.3. OversamplingApplying the backprojection operator, we note thatthe collected data will produce non-zero intensity notonly within but also outside the shaded portion of thespecimen in Fig. 1. We refer to the shaded portion inFig. 1 as ‘‘the region of interest.’’ The extension of thesupport outside the region of interest does not cause anyproblems in the direct summation algorithm. However,by using (22) we see that since gðx þ Mf ; zÞ ¼ gðx; zÞ, wereconstruct a periodic image. If Mf is not large enough,the reconstructed area outside the region of interest willwrap around and overlap with the region of interestcausing undesirable artifacts. By knowing the full extentof the reconstruction, Mext in Fig. 1, we can choose thenumber of frequencies Mf large enough. Since the size ofMf affects the computational speed of the algorithm, it isimportant to choose Mf as small as possible.From Fig. 1 we observe thatDM ¼Mext M¼ N tan hmax ;2NAlgorithm1. Precomputation: For each angle hl and eachfrequency xm , compute Fl ðxm Þ defined by (18).2. For each image, evaluate (22):2.1. For each angle hl and each frequency xm , compute the sum (19) using the USFFT and multiplythe result by Fl ðxm Þ to obtain vl ðxm Þ in (17). Seethe Appendix A for details. Computational cost:OðNh Mf Þ þ OðNh M log MÞ.2.2. For each frequency xm , compute the sum in (23)using the USFFT. See the Appendix A for details.Computational cost: OðMf Nh Þ þ OðMf N log N Þ.2.3. Compute the sum in (22) using the FFT. Computational cost: OðNMf log Mf Þ.h. This is illustrated in Fig. 3,where hmax ¼ maxfjhl jgl¼1where there is no wrap-around between the left and rightsides of the reconstruction.If Mf is smaller than the spatial support Mext ¼M þ 2DM, we will observe aliasing but as long as theoverlapping region is outside the region of interest, noharm is done to the reconstruction. Hence, in order toavoid aliasing artifacts, we must chooseMf P M þ N tan hmax :ð24ÞThis is illustrated in Fig. 4, where the wrap-around doesnot overlap the region of interest. For our applications,we note that choosing Mf according to (24) also is sufficient to achieve the desired accuracy in discretizing theintegral in (21).

K. Sandberg et al. / Journal of Structural Biology 144 (2003) 61–7267Fig. 2. Interpretation of the FFS in the Fourier domain. The measurement given by the data acquisition is represented by the filled dots. Step 2.1 in thealgorithm amounts to interpolating the data to the positions indicated by the squares. Step 2.2 of the algorithm computes the matrix g n ðxk Þ columnwise by adding data along the summation lines. The final step of the algorithm computes the image by applying the inverse FFT row-wise on g n ðxk Þ.Fig. 3. For this reconstruction M ¼ 572, N ¼ 140, and hmax ¼ 73:31 . For this reconstruction, we chose Mf ¼ 1512. We observe that for this choice ofMf , the entire support of the image fits in the reconstructed image.Choosing Mf larger than M can be thought of asoversampling the image. The oversampling factor is given by the ratio Mf M and it is desirable to make thisfactor as close to one as possible. For typical data sets,we have found that this oversampling factor rangesbetween 1.5 and 2. In (24), the size of Mf depends on

68K. Sandberg et al. / Journal of Structural Biology 144 (2003) 61–72Fig. 4. We construct the same data set as in Fig. 3, but here we chose Mf ¼ 1080 which satisfies (24). Here, we do see some wrap-around in the left andright part of the upper image, but it does not overlap the region of interest cropped out in the lower image.N tan hmax . We note that for physical reasons, large hmaxrequires thin specimens, which implies a small N . Hence,the oversampling factor is bounded in most applicationsand we have found it is typically sufficient to chooseMf 6 2M. Therefore, the total computational cost of thefull three-dimensional reconstruction algorithm is givenby OðNy Nh M log MÞ þ OðNy Nh Mf log Mf Þ. This shouldbe compared to OðNy MNNh Þ for the direct summation.Actual speed comparisons are given in Section 5.2.2.the order as a parameter in the algorithm. Such choiceshould be guided by practical considerations and experiences. Since using high-order interpolation has abandlimiting effect, it also implies that we need fewerfrequencies for the reconstruction which, in turn, provides a marginal speed improvement.4.4. Interpolation5.1. Incorporation of the FFS algorithm into IMODIn Sections 3 and 4, we matched the FFS to reproduce the result of the direct summation algorithm thatuses piecewise linear interpolation of the data. In directsummation, the interpolation is applied in step 2.2 of thealgorithm in Section 2.3. We now show that higher order interpolation schemes can be easily incorporatedinto the FFS without additional computational costs.As in the case of linear interpolation, the piecewiseinterpolation of higher order can be described in the spacedomain by using the B-splines, (Schoenberg, 1973). Forthe linear interpolation, the linear B-spline is given by thehat function in (14). Although higher order interpolationschemes may be cumbersome to implement in the spacedomain, in Fourier based algorithms the interpolation isimplemented in the Fourier domain where it is simply amultiplication by an appropriate filter. From the definition of the B-splines as a repeated convolution of thecharacteristic function, it follows that the Fourier transform of the B-spline of odd order k is given by kþ1sinðpxÞb ðkÞ ðxÞ ¼; k ¼ 1; 3; 5; . . .ð25ÞpxThe FFS algorithm was incorporated into the Tiltprogram of the IMOD package (http://www.bio3d.colorado.edu/imod; Kremer et al., 1996). The basicframework of this program was originally written byMike Lawrence while at the Medical Research Councilin Cambridge. Before the speed comparisons reportedhere, the direct summation procedures in Tilt were optimized in two ways. First, all boundary checks weremoved outside of the inner summation loops. Second,each line of input data was stretched by the cosine of theThus, all we need to use higher order interpolation inour reconstruction algorithm, is to apply (25) whencomputing the factor Fl ðxÞ given by (18). This is done instep 1 in the algorithm in Section 4.2.Since the higher order interpolation filter (25) decaysfaster with order k, increasing the order of interpolationeffectively low-passes the data. We leave the choice of5. ResultszαyFig. 5. x-axis tilting. Side view of section to be reconstructed where thesection is tilted by the angle a around the x-axis (pointing out from thepaper). Dashed lines are slices computed by the algorithm and the solidline is the output slice interpolated from vertical slices.

K. Sandberg et al. / Journal of Structural Biology 144 (2003) 61–72tilt angle, with the stretched data oversampled by afactor of two to minimize the filtering effect of interpolating twice. The result is that a line of stretched inputdata can be added into a line of the tomographic slice bystepping through the input line at a fixed interval andwith fixed interpolation factors. The advantages of FFSwould have been considerably higher than describedhere without these improvements to the original code.Further developments in the Tilt program and in theFFS algorithm were spurred by the desire to correct forthe plane of the specimen not being parallel to the tilt69axis. When the specimen is tilted in the plane of the reconstruction, the thickness of the reconstruction must beincreased to contain the material of interest, sometimesby as much as 50% for large data sets. To avoid this, theTilt program can apply a tilt around the x-axis whencomputing the reconstruction by direct summation.However, this means that one output slice is based onmultiple lines of input data, making fast backprojectionunusable for this computation. To solve this problem,the FFS algorithm was modified to shift the outputslice vertically (in the z dimension) by a chosen amount.Fig. 6. Image reconstructions. (A) A test data set computed with the FFS algorithm. (B) The difference between the image in (A) and the corresponding slice reconstructed by direct summation with the same contrast as in (A). (C) The same dataset as in (B), but with the contrast amplified 29times. (D) Test data set for the Fourier ring correlation test. (E) Portion of reconstruction using the FFS of the data set in (D) without noise added.(F) Portion of reconstruction using the FFS of the data set in (D) with noise added.

70K. Sandberg et al. / Journal of Structural Biology 144 (2003) 61–

Gilbert (1972) via direct summation (for a review, see Frank, 1992). The well-known Fourier slice theorem relates pro-jection data to the Fourier transform of the image. The one-dimensional Fourier transform of the collected projection data corresponds to samples on a polar grid in the Fourier domain where, in our case, the polar