Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
IEEE Trans Nucl Sci. Author manuscript; available in PMC 2017 September 20.
Published in final edited form as:
PMCID: PMC5606164

On Few-View Tomography and Staircase Artifacts


This paper investigates some potential methods for few-view tomography, and investigates the cause and remedy of the staircase artifacts. This paper also discusses whether there is any benefit to use edge-preserving filters in emission tomography. We formulate a general Green’s one-step-late algorithm, so that it can incorporate any linear or non-linear filters. We argue that the derivative of the penalty function can be “artificially” created, not naturally derived from a penalty function. We have gained more insight into constrained image reconstruction especially with edge-preserving constraints. Our numerical results show that the bilateral method can have better performance than the TV method, and higher order methods may not be necessary for non-piecewise-constant objects.

Index Terms: Artifacts, edge-preserving filters, image reconstruction, iterative algorithm, tomography


SINCE the early 1990’s, the TV (total variation) minimization methods have become popular [1]–[4]. In x-ray CT and MRI, the TV minimization methods have wide applications, especially when data are under-sampled. The TV minimization methods are fairly successful and can produce almost artifact-free images in few-view low-noise tomography [5]–[12]. As an alternative methed, a model-based solution is suggested in [13].

The TV method may give a staircase appearance when the true object is smooth. The TV is calculated with only the first-order derivatives. When the TV-norm is minimized, it is believed that piecewise constant regions are encouraged. Based on this consideration, methods using higher-order derivatives were suggested to battle the staircase artifacts or for non-piecewise-constant objects [14]–[17].

Some initial attempts have been made in applying the TV method to emission tomography [18]–[20]. This paper will focus on emission tomography, in which the images are not piecewise constant and the noise magnitude can be larger than the edge contrast. In this paper, we make the following arguments: (i) Green’s one-step-late EM (expectation maximization) algorithm can be used to incorporate most penalty functions and non-linear filters [21], [22]. (ii) The TV-norm minimization is not the only way to reconstruct few-view tomographic images. Any edge-preserving filter can be used in few-view tomography. (iii) The few-view tomography can work for non-piecewise images. (iv) The higher-order derivative methods can cause staircase artifacts too. (v) The staircase artifacts can be reduced by using proper control parameters.


A. Green’s One-Step-Late Algorithm

Developed for emission tomography with the Poisson noise model, Green’s one-step-late algorithm is an iterative algorithm given as [21]


where xi,j(n) is the value of image pixel (i, j) at the nth iteration, pk is the kth projection, a(i, j) k represents the contribution from the image pixel (i, j) to the kth projection, Ui,j(n) is the derivative of the penalty function with respect to the image pixel (i, j) at the nth iteration, and β is a control parameter. In (1), the backprojection of constant 1 should be normalized to one; otherwise the controlling parameter β will depend on the value of the backprojection of constant 1 [23]. In fact, in the ML-EM algorithm xi,j(n+1)=xi,j(n)ka(i,j)kka(i,j)kpkî,ĵa(î,ĵ)kxî,ĵ(n), one can always scale the coefficients a(i, j) k by any positive constant C and scale the image xi,j(n) by 1/C: (1Cxi,j(n+1))=(1Cxi,j(n))kCa(i,j)kkCa(i,j)kpkî,ĵCa(î,ĵ)k(1Cxî,ĵ(n)) without changing the algorithm nor affecting its convergence. Normalization ka(i,j)k=1 makes the algorithm more consistent with its derivation, in which the coefficient a(i, j) k is assumed to be the probability. This normalization may make the selection of β a little easier in Green’s algorithm. However, in common implementation the coefficient a(i, j) k is the line-length within the pixel and scaled by its attenuation factor.

Laplacian, and bilateral filters

The penalty function P associated with Ui,j(n) can be designed and selected by the user. If the penalty function is chosen as the TV-norm of the image at the nth iteration, (1) becomes the TV regulated EM algorithm that we developed in 1999 [22]. For a discrete image, the TV-norm penalty function can be symbolically expressed as


The penalty function (2) is not quadratic; it is not differentiable when the quantity under the square root is zero.

When the quantity under the square root is positive, the partial derivative of P with respect to pixel (i, j) is readily calculated as [22]


When the quantity under the square root is zero, the quantity already reaches its minimum and the penalty function is no longer needed. In this situation, it is safe to set the derivative to zero. One easy way to handle this situation in an actual implementation of (3) is to add a very small positive constant ε, say, 10−8, under the square roots in all denominators. Thus (3) becomes


It is realized that the small parameter ε is artificially introduced. Strictly speaking, Ui,j in (4) is not the derivative of the TV penalty function any more.

Most penalty functions are differentiable almost everywhere. Some penalty functions may or may not be differentiable at their minimum, where we can always define the derivative to be zero. With this in mind, Green’s algorithm (1) is versatile; it does not have strict requirements for the penalty function P. Surrogate functions, which are required for some sophisticated algorithms, are not necessary here.

In fact, one can get more insight of how the TV-norm performs from (3) or (4) than from (2). The algorithm drives Ui,j to zero, because the algorithm is designed to minimize the penalty function, which forces its derivative Ui,j to zero.

The three numerators in Eq. (4) try to satisfy xi,j = (xi,j+1 + xi+1,j)/2, xi,j = xi,j−1, and xi,j = xi−1,j, respectively. In other words, the modified TV-norm (with parameter ε) encourages the image pixel value to be close to those of its neighbors. Thus, (4) encourages a smooth reconstruction, and the modified TV algorithm should perform reasonably well if the object is not piecewise constant. The magnitude of the gradient forms the denominators. The gradient magnitude is the key factor for the edge-preserving effect. A large gradient is assumed to be associated with an edge, and a small gradient is assumed to be associated with the noise or artifacts. When the gradient is large, Ui,j is small and not effective. When the gradient is small, Ui,j is large and effective. As a result, noise and artifacts are reduced and the edges are kept.

Similar to (2), a Laplacian penalty function P can be set up by replacing the first order finite differences by the second order ones, for example, [see (5) at the bottom of the next page]. Taking the partial derivative of P with respect to xi,j, we obtain [see (6) at the bottom of the next page].



A quick observation of (6) is that the numerators encourage smoothness of the image, which is desirable; the denominators are the magnitude of the curvature and do not contain edge information, which is not desirable.

We argue that the function Ui,j can be “artificially” created, not naturally derived from a penalty function. This novel idea opens a door for us to obtain many new algorithms, in which the penalty function is not well defined.

The drawback of (6) can be fixed artificially. We can only use the first term of (6) and replace its denominator so that an edge can be detected. One possible option is [see (7) at the bottom of the next page]. In this paper, we refer to (7) as a modified Laplacian penalty gradient, even though it is not quite the gradient of the Laplacian. This example shows the advantage of designing the gradient function Ui,j directly, without knowing the original penalty function P, which may not exist or its derivative may not exists.


There are many effective non-linear filters that can suppress the noise and preserve edges such as the median filters and the bilateral filters [24][25]. We do not have a penalty function to characterize these filters. Also, these non-linear filters cannot be implemented as a convolution.

We claim that Green’s algorithm (1) is able to incorporate these non-linear filters readily. We suggest the following general gradient formula


where x¯i,j(n) is pixel (i, j) of the filtered image, and the filter’s input is the image obtained at the nth iteration of (1). The filter can be any linear or non-linear filter. We will consider a bilateral filter in this paper.

A bilateral filter has many forms. One of them has the following input-output relationship [25]:


where Ω(i, j) is a sub-region of the image centered at pixel (i, j), and δ is a positive parameter, controlling the shape of the Gaussian function. This Gaussian function contains the edge information. The parameter δ defines whether an edge or noise is encountered. If it is an edge, exp(−δ·(xm,nxi,j)2) should be small (close to 0); otherwise, exp(−δ · (xm,nxi,j)2) should be large (close to 1). In this paper, Ω(i, j) is chosen as a 3 × 3 region.

In the general form (8) the x¯i,j(n) is a filtered value which may or may not have a closed-form expression. In 2002, Alenius and Ruotsalainen proposed a generalized median filter to define x¯i,j(n), which does not have a closed-form expression yet yields some good results with proper parameter setup [29].

B. Computer simulations and brain phantom experiment

The focus of this paper is on emission tomography, and the objects are not piecewise constant. Green’s one-step-late algorithm is especially developed for emission tomography with Poisson noise, even though it can be extended to transmission tomography [26].

The computer generated phantom consists of a large disc and a smaller disc inside. The discs have curved (not flat) intensity functions. The imaging geometry is parallel beam, with 128 detector bins on the detector. The image array is 128 × 128. The image reconstructed using 120 views over 360° and noiseless data is considered to be the standard for other images to compare with.

The purpose of the noiseless computer simulation studies is to investigate the performance of the reconstruction algorithms in few-view tomography, where only 20 view angles were used over 360°. Green’s algorithm with three different penalty gradient functions was used for image reconstruction: modified TV, modified Laplacian, and bilateral. Different control parameters were used to demonstrate the parameter effects.

In the next set of the studies, the projections were made noisy by using the Poisson noise model. We also used 20 views. The purpose is to show that noisy data can result in staircase artifacts for all 3 edge-preserving methods.

Finally, a Hoffman brain phantom was scanned using a Siemens Signature SPECT system with a set of parallel-beam collimators, and projections of 120 views were acquired over 360°. The phantom was filled with Tc-99 m, using a clinical dose. In SPECT, data acquisition with 120 views is considered as full data acquisition, and there should be no few-view tomography artifacts. Due to noise, staircase artifacts can be produced if the algorithm’s controlling parameters are not set properly. In all image reconstruction tasks in this paper, 50 iterations were used to produce the final image.


Fig. 1 is used as a standard for other computer simulation reconstructions to compare with. This image was reconstructed by the conventional MLEM algorithm [27] with 120 views over 360° and the projections were noiseless. Both the large and small discs in the phantom had curved (instead of flat) intensity values, as shown by the central horizontal profile.

Fig. 1
Computer generated phantom reconstructed with the MLEM algorithm using 50 iterations. The number of views is 120 and the projections are noiseless. The central horizontal profile of the image is shown below the image.

When the number of views was decreased to 20, the typical star-like few-view tomography artifacts appeared, as seen in Case 3 of Fig. 2. When Green’s one-step-late algorithm was applied with the three suggested penalty gradients: modified TV, modified Laplacian, and bilateral, the star-like artifacts were reduced, as shown in Case 2 of Fig. 2.

Fig. 2
Computer generated phantom reconstructed with Green’s one-step-late algorithm using 50 iterations. The number of views is 20 and the projections are noiseless. The central horizontal profile of each image is shown below the image. Severe star-like ...

In modified TV and modified Laplacian methods, the controlling parameter is β. When β = 0, no penalty is used. When β is set to an appropriate range of values, the algorithm is able to suppress the star-like artifacts, whose magnitude is smaller than the edge contrast. However, if the β value is increased beyond this range, the artifacts come back.

In a sense, β is the slope of the penalty function; the slope corresponds to the magnitude of the gradient, which sets up a range for variation to be suppressed. When this range contains the artifact magnitude, the artifacts can be reduced, as shown in Case 2 of Fig. 2. A larger β gives a larger range of image variation to be suppressed. If the β value is selected too large, both the edge contrast and the artifact magnitude will be reduced simultaneously, which is not an effective way to battle artifacts, as shown in Case 1 of Fig. 2.

It is not easy to select the control parameter(s) in a MAP (maximum a posteriori) algorithm, because the results are very sensitive to the parameter(s). Usually, there is a narrow range for each parameter to give a desirable result. The shape of the penalty function is controlled by the parameter(s). For example, we could use the Lp-norm of the image as the penalty function. For some p values, this penalty function encourages a smooth image. However, for some other p values, it encourages a piecewise-constant image. It is not clearly known to us what functional shape encourages what image features.

The parameter β not only controls the target to be suppressed but also controls the convergence behavior of the algorithm. We do not completely understand how the parameter β affects the algorithm. We select β by trial-and-error. For a larger parameter β, the algorithm becomes more unstable or even diverges. It is an open question how to design a penalty function (and its associated parameters) to encourage or suppress a targeted image feature.

For the bilateral method, there are two controlling parameters: β and δ, and both of them can be used to select a range of variation magnitude to be suppressed. There is no known proof of the convergence of Green’s one-step-late algorithm. Experiences teach us that the convergence may be achieved by using a small parameter β. In Green’s algorithm, parameter β cannot be chosen arbitrarily. If β is too large, βUi,j(n) can be a negative value with a large magnitude, which may make ka(i,j)k+βUi,j(n) in (1) close to zero and make the algorithm diverge. The advantage of using the bilateral method is that we have more control, using δ to adjust the filtering range and using β to control the convergence. If the artifacts are within the filtering range and the edge contrasts are not, the artifacts can be effectively reduced, as shown in Case 2 of Fig. 2. If the range is too large—large enough to include the edge magnitude, the artifact reduction is not effective while the edge is not preserved.

Fig. 2 also tells us that the modified TV method is not the only method for few-view tomography, and it may not be the optimal one.

Generally speaking, emission tomography images are much noisier than transmission tomography images due to the dramatic differences in photon counts. The edge-preserving filters may not work well for noisy emission images, where the noise magnitude is most likely larger than the edge contrast. The edge-preserving filters can cause staircase artifacts as demonstrated in Case 2 of Fig. 3, where many faulty edges are created by noise. One way to reduce the staircase artifacts is to diminish the edge-preserving property of the filters, see Case 1 of Fig. 3.

Fig. 3
Computer generated phantom reconstructed with Green’s one-step-late algorithm using 50 iterations. The number of views is 20 and the projections are added with Poisson noise. The central horizontal profile of each image is shown below the image. ...

Another side effect of edge-preserving filters is the inability to remove the salt-and-pepper noise. When the magnitude of the salt-and-pepper noise is large, the noise will be mis-classified as “edges” and kept. An effective filter to remove the salt-and-pepper noise is the median filter [24].

A median filter can be used in Green’s algorithm to replace the modified TV, modified Laplacian, or bilateral [29]. However, we chose to implement a 3-point median filter in the projection domain as a pre-filter. One reason for this choice is that it is computationally more efficient to implement a 3-point median filter in the projection domain than in the image domain. Another reason is that a pre-filter only performs once, independent from the iteration number, for any image reconstruction task. The 3-point median filter is implemented as follows:



where pi is the input of the filter and [p with macron]i is the output of the filter. The effectiveness of this 3-point median filter is shown is Fig. 4. Image reconstruction methods are the same in Fig. 3 and Fig. 4; the only difference is that the noisy projections are pre-filtered by the 3-point median filter for the images in Fig. 4. The median filter effectively reduces the salt-and-pepper noise.

Fig. 4
Computer generated phantom reconstructed with Green’s one-step-late algorithm using 50 iterations. The number of views is 20 and the projections are added with Poisson noise. The central horizontal profile of each image is shown below the image. ...

Fig. 5 shows the results for the real Hoffman data acquired by a Siemens Signature SPECT system. The SPECT projections are pre-filtered by (10). Results in Case 2 of Fig. 5 have severe staircase artifacts. Better results are obtained by using parameters that down-play the edge-preserving role and are shown in Case 1 of Fig. 5. This real-data example shows that the edge-preserving property is not helpful in very noisy emission tomography.

Fig. 5
The real Hoffman phantom reconstructed with Green’s one-step-late algorithm using 50 iterations. The number of views is 120 and the projections are acquired with a Siemens SPECT system. The central horizontal profile of each image is shown below ...


The main goal of this paper is to tell the reader that we have more freedom to develop a MAP (maximum a posteriori) algorithm by designing the gradient function U directly, without worrying too much about whether the original penalty function P is differentiable. The original penalty function P may not exist or may not have a closed-form. The knowledge of P is somehow used to design U. We can see this in two different angles. View point #1: We start with a penalty function P, and take the derivative of P to get U. When U is not well defined everywhere, some modifications may be made. For example, an artificial parameter ε can be introduced. The revised U is no longer the exact derivative of the original penalty function. View point #2: A symbolic penalty function P is set up as P = 0.5 · (x[x with macron])2, where [x with macron] is a filtered version of x. The filter can be linear or non-linear. We do not know whether the “symbolic” penalty P = 0.5 · (x[x with macron])2 is convex and has a unique minimum. Therefore, we can always start with the definition of U = x[x with macron], and not to worry about whether P is differentiable, because [x with macron] may not have a closed-form mathematical expression.

This paper argues the following points:

  1. It is known that if the object is not piecewise-constant, some higher order methods are effective [14]–[17]. This paper argues that higher order methods may not be necessary. The modified gradients U for TV, Laplacian, and bilateral methods are all effective, thank to the artificially introduced parameter(s), as demonstrated by our noiseless computer simulations with a phantom that has a curved intensity function.
  2. Few-view tomography is possible only when the magnitude of the star-like artifacts is smaller than the edge contrast and the parameters are properly chosen.
  3. Green’s one-step-late algorithm is versatile, and almost any penalty gradients or filters can be readily incorporated according to (8). They include the popular edge-preserving filters such as TV, Laplacian, bilateral, median [29], Huber [28], and so on. As a matter of fact, any gradient-type algorithm shares this versatility. We prefer Green’s algorithm because it models the Poisson noise and produces non-negative images.
  4. The TV method is not the only method that can be used in few-view tomography. It may not be the optimal method. The bilateral method is our preferred method, because the bilateral method has one more control parameter. Having one more control parameter may give us more chance to find a better solution; it also makes it more difficult to find it.
  5. The penalty function may not exist; the penalty function may not be differentiable. It is possible that one can design the penalty gradient directly without knowing the original penalty function, as illustrated by (7). The desired image is when the gradient is zero. This opens a research area of designing the gradient function directly, which is able to suppress some targeted image features.
  6. The staircase artifacts are the side product of any edge-preserving filters, when the data are noisy. In order to reduce staircase artifacts, one can make the filters less edge-preserving. Edge-preserving filters may not be desirable for emission imaging, where the noise amplitude may be larger than the edge contrast. This point is illustrated by the real Hoffman brain phantom study.
  7. Using a small kernel median filter is effective in reducing the salt-and-pepper noise. This approach was suggested by You et al. [30]. Alenius and Ruotsalainen suggested another method that applies the median filter in the image domain [29].


1. Rudin L, Osher S, Fatemi E. Nonlinear total variation based noise removal algorithms. Phys. D. 1992;60:259–268.
2. Guichard F, Malgouyres F. Total variation based interpolation. Proc. Eusipco. 1998;3:1741–1744.
3. Alter F, Durand S, Froment J. Adapted total variation for artifact free decompression of JPEG images. J. Math. Imag. Vis. 2005;33(2):199–211.
4. Candes E, Romberg J, Tao T. Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information. IEEE Trans. Inf. Theory. 2006;52:489–509.
5. Bian J, Siewerdsen JH, Han X, Sidky EY, Prince JL, Pelizzari CA, Pan X. Evaluation of sparse-view reconstruction from flat-panel-detector cone-beam CT. Phys. Med. Biol. 2010;55:6575–6599. [PMC free article] [PubMed]
6. Han X, Bian J, Eaker DR, Kline TL, Sidky EY, Ritman EL, Pan X. Algorithm-enabled low-dose micro-CT imaging. IEEE Trans. Med. Imag. 2011;30:606–620. [PMC free article] [PubMed]
7. Han X, Bian J, Ritman EL, Sidky EY, Pan X. Optimization-based reconstruction of sparse images from few-view projections. Phys. Med. Biol. 2012;57:5245–5273. [PMC free article] [PubMed]
8. Ritschl L, Bergner F, Fleischmann C, Kachelrieb M. Improved total variation-based CT image reconstruction applied to clinical data. Phys. Med. Biol. 2011;56:1545–1561. [PubMed]
9. Sidky EY, Kao C-M, Pan X. Accurate image reconstruction from few-views and limited-angle data in divergent-beam CT. J. X-Ray Sci. Technol. 2006;14:119–139.
10. Sidky EY, Pan X. Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization. Phys. Med. Biol. 2008;53:4777–4807. [PMC free article] [PubMed]
11. Tang J, Nett BE, Chen G-H. Performance comparison between total variation (TV)-based compressed sensing and statistical iterative reconstruction algorithms. Phys. Med. Biol. 2009;54:5781–5804. [PMC free article] [PubMed]
12. Chen G-H, Thériault-Lauzier P, Tang J, Nett B, Leng S, Zambelli J, Qi Z, Bevins N, Raval A, Reeder S, Rowley H. Time-resolved interventional cardiac C-arm cone-beam CT: An application of the PICCS algorithm. IEEE Trans. Med. Imag. 2012;31(4):907–923. 10.1109/TMI.2011.2172951. [PMC free article] [PubMed]
13. Liu Y, Liang Z, Ma J, Lu H, Wang K, Zhang H, Moore W. Total Variation-Stokes Strategy for Sparse-View X-ray CT Image Reconstruction. IEEE Trans. Med. Imag. 2014;33(3):749–763. [PMC free article] [PubMed]
14. Buades A, Coll B, Morel J-M. The staircasing effect in neighborhood filters and its solution. IEEE Trans. Image Process. 2006;15(6):1499–1505. [PubMed]
15. Chan T, Marquina A, Mulet P. High-order total variation based image restoration. SIAM J. Sci. Comput. 2000;22(2):503–516.
16. Yang J, Yu H, Jiang M, Wang G. High order total variation minimization for interior tomography. Inverse Problems. 2010;26(3) 10.1088/0266-5611/26/3/035013, 035013. [PMC free article] [PubMed]
17. Lee M, Ward JP, Unser M, Ye JC. Proc. 3rd Int. Conf. Image Formation in X-Ray Computed Tomography. Salt Lake City, UT, USA: 2014. Multiscale interior tomography using 1D generalized total variation; pp. 347–350.
18. Zhang J, Krol A, Li S, Wu Z, Vogelsang L, Shen L, Schmidtlein C, Lipson E, Xu Y, Feiglin D. Comparison of SPECT lesion detectability between preconditioned alternating projection algorithm (PAPA) and expectation-maximization algorithm with total variation regularizer (EM-TV) J. Nucl. Med. 2014;55(Supplement 1):543.
19. Wolf PA, Jorgensen JS, Schmidt TG, Sidky EY. Few-view single photon emission computed tomography (SPECT) reconstruction based on a blurred piecewise constant object model. Phys. Med. Biol. 2013;58:5629. 10.1088/0031-9155/58/16/5629. [PMC free article] [PubMed]
20. Sawatzky A, Brune C, Wubbeling F, Kosters T, Schafers K, Burger M. Accurate EM-TV algorithm in PET with low SNR; Proc. Nuclear Science Symp. Conf. Rec; 2008. pp. 5133–5137. 10.1109/NSSMIC.2008.4774392.
21. Green PJ. Bayesian reconstruction from emission tomography data using a modified EM algorithm. IEEE Trans. Med. Imag. 1990;9:84–93. [PubMed]
22. Panin VY, Zeng GL, Gullberg GT. Total variation regulated EM algorithm. IEEE Trans. Nucl. Sci. 1999;46:2202–2210.
23. Liang Z, Turkington T, Gilland D, Jaszczak R, Coleman E. Simultaneous compensation for attenuation, scatter, and detector response for SPECT reconstruction in three dimensions. Phys. Med. Biol. 1992;36(3):587–603. [PubMed]
24. Huang T, Yang G, Tang G. A fast two-dimensional median filtering algorithm. IEEE Trans. Acoust., Speech, Signal Process. 1979;27(1):13–18.
25. Tomasi C, Manduchi R. Bilateral filtering for gray and color images; Proc. IEEE 6th Int. Conf. Computer Vision; 1998. pp. 839–846.
26. Zeng GL, Li Y, Zamyatin A. Iterative total-variation reconstruction vs. weighted filtered-backprojection reconstruction with edge-preserving filtering. Phys. Med. Biol. 2013;58:3413–3431. [PMC free article] [PubMed]
27. Langer K, Carson R. EM reconstruction algorithms for Emission and Transmission tomography. J. Comp. Assist. Tomogr. 1984;8:302–316.
28. Huber PJ. Robust Estimation of a Location Parameter. Ann. Stat. 1964;53:73–101.
29. Alenius S, Ruotsalainen U. Generalization of median root prior reconstruction. IEEE Trans. Med. Imag. 2002;21(11):1413–1420. [PubMed]
30. You J, Zeng GL, Liang Z. FBP algorithms for attenuated fan-beam projections. Inverse Problems. 2005;21:1179–1192.