|Home | About | Journals | Submit | Contact Us | Français|
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.
SINCE the early 1990’s, the TV (total variation) minimization methods have become popular –. 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 –. As an alternative methed, a model-based solution is suggested in .
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 –.
Some initial attempts have been made in applying the TV method to emission tomography –. 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 , . (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.
Developed for emission tomography with the Poisson noise model, Green’s one-step-late algorithm is an iterative algorithm given as 
where 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, 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 . In fact, in the ML-EM algorithm , one can always scale the coefficients a(i, j) k by any positive constant C and scale the image by 1/C: without changing the algorithm nor affecting its convergence. Normalization 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.
The penalty function P associated with 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 . 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 
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 . 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 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 :
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,n − xi,j)2) should be small (close to 0); otherwise, exp(−δ · (xm,n − xi,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 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 , which does not have a closed-form expression yet yields some good results with proper parameter setup .
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 .
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  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.
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.
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, can be a negative value with a large magnitude, which may make 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.
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 .
A median filter can be used in Green’s algorithm to replace the modified TV, modified Laplacian, or bilateral . 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 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. 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.
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 − )2, where 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 − )2 is convex and has a unique minimum. Therefore, we can always start with the definition of U = x − , and not to worry about whether P is differentiable, because may not have a closed-form mathematical expression.
This paper argues the following points: