Priority-based Mojette reconstruction from sparse noisy projections

1 Introduction

Computed Tomography (CT) technology can effectively reconstruct the internal structure of the object; thus, it has been widely used in the field of industrial detection and medical diagnosis. In the reconstruction technique, Radon transform is the theoretical basis of exact reconstruction in continuous domain, whereas the detectors in a real CT system are discrete and limited. Therefore, the image has to be sampled to acquire projection data, resulting in approximate reconstruction error. The related reconstruction algorithms have been researched for decades and also achieved high-quality results in the finite angular range [1–3]. In practice, the reconstruction of images from discrete projections is a wide field with a long history and various forms of Discrete Radon Transforms (DRT) have been constructed to take advantage of digital geometry [4–7]. Guédon et al. proposed Mojette transform [8] which is a direct discretisation of the RT and can directly implement discrete exact reconstruction. It can overcome approximate error caused by inverse Radon transform in the finite domain. As an entire discrete form of Radon transform, it provides new ideas for exact sparse-view CT reconstruction. It is inspiring to see that, the conversion problem from the acquired Radon projection to Mojette projection had obtained some achievements in [9–11].

However, the reconstruction methods based on Mojette transform are very sensitive to noise and the noise becomes even worse in a finite sampling condition. Because noise added to projection data would be gradually propagated in the iterative reconstruction process, the noise-robust Mojette inversion technique becomes a main problem in sparse-view CT reconstruction.

To address this problem, Serviéres [12] proposed the conjugate gradient method (CGM) in 2005, which obtained the optimal solution by solving the objective function in least square sense. Experimental results showed that it could effectively attenuate noise, but this method is not suitably convergent. Normand [13] used a geometrical approach to streamline the reconstruction process and found the compact reconstruction path from the starting point to destination. However, its reconstruction strategy is not beneficial in reducing the effect of noise in the slice image. Recur [14] reduced the accumulated noise with iterations at the cost of increasing the number of projections. Svalbe et al. [15] applied a discrete point-spread function (PSF) that regularised the weighting function. This approach corrected the strong aliasing on the reconstructed image from the back-projection filtration technique by direct deconvolution. However, the recovered image depended on the FFT of the PSF being well-conditioned or regularised. Kingston [16] proposed an inversion approach for Mojettte projections in the Fourier domain, which sped up the reconstruction process using exact frequency data resampling. However, when the number of projections decreased, its anti-noise ability became obviously weaker. Shekhar [17] presented a robust method based on the classical Fourier inversion technique and generated compact sets of rational projections which are directly and exactly converted to Discrete Fourier Transform (DFT) space. The method could reduce the noise via the conversion of Mojette projections to Fourier slices without interpolation. But the reconstruction of large sized images needs to be balanced against the corresponding increase in the computation time. Although it is acknowledged that Mojette reconstruction algorithms have been making some progress towards recovering a slice image from acquired noisy projections, these algorithms are not necessarily effective in the case of sparse-view noisy projections, i.e., Mojette tomographic reconstruction from sparse-view noisy projections is still a main problem.

This paper presents a robust reconstruction approach based on Mojette transform. Given that the error in the slice image is proportional to iterations, in this approach projections for an image reconstruction are taken up in order of priority according to minimising the number of iterations, then the selected prioritized projections are utilized to design a novel reconstruction path. Therefore, the proposed method can effectively suppress the accumulated noise and improve performance without increasing the number of projections, which is more suitable for sparse-view CT reconstruction.

The paper is outlined as follows: Section 2 reviews Mojette transform and its properties, and then analyzes the noise accumulation via back-projection. Section 3 proposes a priority-based Mojette inversion method for sparse-view CT reconstruction and then gives analysis of anti-noise performance. Finally, the results of a comparison experiment are showed in Section 4.

2 Mojette transform and its properties

The Dirac-Mojette transform is described over a set of angles θ _i  = arctan(q _i /p _i ), where p _i and q _i indicates the number of pixel displacements horizontally and vertically, and i is the number of projections. Note that p _i and q _i satisfy GCD (p _i , q _i ) =1. The projection direction is directly expressed as (p _i , q _i ). The Dirac-Mojette transform is defined as [18]: M p i , q i ( b ) = ∑ k = 1 P ∑ l = 1 Q f ( l , k ) δ ( b - q i · ( k - 1 ) - p i · ( l - 1 ) )(1)

where δ (·) is the discrete Kronecker function and the equi-spatial detector array includes b elements (called bins per projection). The value M_{p i ,q i}  (b) of Mojette transform in each bin is the sum of pixels centered on the parallel lines b = p _i  · (l - 1) + q _i  · (k - 1) for projection (p _i , q _i ). For example, Fig. 1 describes Dirac-Mojette transform sampling on a 3 × 3 image for the projections (1,1), (1,2) and (1,0).

OPEN IN VIEWER

Fig.1

Results of Mojette transform sampling on a 3 × 3 image for projections (1,1),(1,2),(1,0). (a) The value of bins in projection (1, 1) is {1, 6, 22, 10, 6}. (b) The value of bins in projection (1, 2) is {1, 4, 7, 9, 15, 3, 6}. (c) The value of bins in projection (1, 0) is {11, 16, 18}.

The sampling is adapted to the discrete geometry of the digital reconstructed image, i.e., the total number of bins per projection has a relationship with both the size of the reconstructed image and the projection (p _i , q _i ). Generally speaking, the total number of bins B for a P × Q image is B = ( P - 1 ) | q i | + ( Q - 1 ) | p i | + 1(2)

Its inverse transform only requires addition operation in spatial domain, and classic method is corner-based Mojette inversion algorithm (CBI). The algorithm utilizing inverse transform reconstructs the image exactly if Katz theory is satisfied [19], i.e.: ∑ i = 1 I p i ≥ P or ∑ i = 1 I q i ≥ Q(3) Where I is the total number of projections.

Fig. 2(a) shows the sampling for a 2 × 2 image with the pixels x₁ = 5, x₂ = 0, x₃ = 2, x₄ = 18 in projections (-1,1),(1,0), where the acquired projections are p₁ = x₁ + x₂ = 5, p₂ = x₃ + x₄ = 20, p₃ = x₃ = 2, p₄ = x₁ + x₄ = 23, p₅ = x₂ = 2. The result of filter back-projection (FBP) based on Radon transform is an approximate image in Fig. 2(b). By contrast, traditional CBI firstly reconstructs the unknown pixel x₂, whose corresponding line p₅ just passes through the pixel itself. Similarly, the pixel x₃ is solved from the projection line p₃. And then the reconstructed pixels x₂ and x₃ are subtracted from all the projections, and so on. The image is iteratively reconstructed pixel by pixel from corner to center. Due to the number of projections I = 2 satisfying Katz theory in Equation (3), i.e., ∑ i = 1 2 | p i | = | - 1 | + | 1 | = 2 ≥ P , the original image can be exactly reconstructed in Fig. 2(c).

OPEN IN VIEWER

Fig.2

Sampling and reconstruction for a 2 × 2 image. (a) Sampling for original image. (b) Result of filter back-projection. (c) Result of inverse Mojette transform.

To compare CBI with traditional Radon inversion techniques (Algebraic Reconstruction Technique + Total Variation and Filter Back-projection), we reconstruct the 64 × 64 Shepp-Logan phantom in Fig. 3(a) which is generated by Matlab 2012a. The acquired 5 projections are uniformly distributed in the angular range [0, π] for Algebraic Reconstruction Technique + Total Variation and Filter Back-projection algorithms (ART+TV [20, 21] and FBP), and there are 1009 bins per projection. CBI reconstructs shepp-logan phantom from Dirac-Mojette transform sampling with the five projections (15, 1) (-9, 7), (13, 3), (-14, 1), (-13, 2). Fig. 3 proves that, only CBI can exactly reconstruct the original signal for sparse-view CT.

OPEN IN VIEWER

Fig.3

Comparison of results using three kinds of reconstruction algorithms. (a) Original Image. (b) Reconstruction result by traditional CBI from five non-uniformly projections. (c) Reconstruction result by ART+TV from five uniformly projections. (d) Reconstruction result by FBP from five uniformly projections.

However, traditional Mojette inversion reconstruction is sensitive to noise. Each Mojette back-projection process is defined as an iteration. In the initial iteration, pixels in the corner region of an image are directly reconstructed by the corresponding noise projection just passing through the pixels. In the next iteration, the noise will spread from the pixels to bins when updating the projection data. Then the noise is propagated from bins to other pixels in the following iteration and so on. From CBI reconstruction process, it is known that the same image can have different reconstruction paths, and the different paths will cause different accumulation noise. For example, Fig. 4 shows that different sets of projections corresponding to different reconstruction paths in CBI reconstruction process [13]. Using three projections (1, 1), (-1, 1), (1, 2), Fig. 4(a) recovers the image from pixel 1 and pixel 2 to black pixel along bold lines simultaneously. Using four projections (-1, 2), (-2, 1), (1, 2), (2, 1), Fig. 4(b) recovers the image from pixels 1, 2, 3 and 4 to the black pixel along black paths simultaneously. Therefore, different reconstruction paths correspond to different accumulation noise. To address the problem of CBI noise sensitivity, a new scheme of designing novel paths and selecting the corresponding projection subsets are proposed in following section.

OPEN IN VIEWER

Fig.4

Diagram of different reconstruction paths. (a) Reconstruction path with projections (1,1),(1,2),(-1,1) and noise is propagated from pixel 1 and pixel 2 to black pixel along the bold lines. (b) Reconstruction path with projections (2,1),(-2,1),(1,2),(-1,2) and noise is propagated from pixels 1, 2, 3 and 4 to black pixel along the black paths.

In this section, Priority-based Mojette inversion algorithm (PBI) with the minimum accumulation noise is proposed. Then, a simple example is shown to suppress noise accumulation by minimizing the number of iterations, followed by analysis of anti-noise performance.

3.1 PBI algorithm with the minimum accumulation noise

Traditional CBI is accurate and reliable in the absence of noise, but as we showed in Section 2, the same image can have different reconstruction paths and the different paths will cause different accumulation noise. To address the problem, this paper actively chooses a set of projections to design a reconstruction path with the minimal accumulation noise. As mentioned above, the error in the slice image is directly related with the number of iterations, i.e., a lower number of iterations results in less error in the reconstructed image. Aiming to minimize the total number of iterations for an image to be completely reconstructed, we establish a kind of priority list for projections where the selected projection in each iteration can recover the maximum of pixels. The proposed PBI algorithm, as an improved CBI, can approximately get the minimum number of projections to reconstruct an image, its reconstruction is summarized in following procedure:

(1) Setting up index matrix for each projection: The image is iteratively reconstructed pixel by pixel from corner to center. To determine which pixels to be reconstructed in each iteration, we associate the pixel coordinates with the parallel rays per projection. The index matrix Proj_P×Q represents the correspondence between the coordinate (l, k) of the image pixels and the parallel rays b with the slope q _i /p _i , i.e., b = p _i  · (l - 1) + q _i  · (k - 1). The dimension of index matrix Proj_P×Q is the same as that of the reconstructed image f (l, k). The value of each element Proj_P×Q (k, l) is equal to b, where b denotes that the pixel coordinate (l, k) is centered on the ray b = p _i  · (l - 1) + q _i  · (k - 1). Note that if the ray b only passes through the center of one pixel (called single bin-pixel correspondence), the bin value on the ray would be directly mapped to the pixel.

Fig. 5 shows index matrices for projections (1, 1), (1, 2), (1, 0). Fig. 5(a) is a 3 × 3 original image. Fig. 5(b) shows that in projection (1, 1), the pixel coordinate (1, 1) center on the ray b = p · (l - 1) + q · (k - 1) =1 · (1 - 1) +1 · (1 - 1) =0, the pixel coordinate (1, 2) on the ray b = p · (l - 1) + q · (k - 1) =1 · (1 - 1) +1 · (2 - 1) =1, the pixel coordinate (2, 1) on the ray b = p · (l - 1) + q · (k - 1) =1 · (2 - 1) +1 · (1 - 1) =1 and so on. There is index matrix [012 ; 123 ; 234] for projection (1,1). The number of single bin-pixel correspondence is two, i.e., the ray b = 0 only passes through the center of pixel 1 and the ray b = 4 only passes through the center of pixel 6. The bin values bin₀ = 1, bin₄ = 6 on the rays are directly mapped to the associated pixels. Similarly, Fig. 5(c) corresponds to index matrix [024 ; 135 ; 246] for projection (1, 2). The number of single bin-pixel correspondence is five, i.e., the ray b = 0 only across pixel 1, the ray b = 1 only across pixel 4, the ray b = 3 only across pixel 9, the ray b = 5 only across pixel 3, the ray b = 6 only across pixel 6. The bin values bin₀ = 1, bin₁ = 4, bin₃ = 9, bin₅ = 3, bin₆ = 6 on the rays are directly mapped to the associated pixels. Likewise, Fig. 5(d) corresponds to index matrix [000 ; 111 ; 222] for projection (1,0). There is no single bin-pixel correspondence.

OPEN IN VIEWER

Fig.5

Diagram of index matrices for projections (1, 1), (1, 2), (1, 0). (a) A 3 × 3 original image. (b) Index matrix [012 ; 123 ; 234] for projection (1, 1). (c) Index matrix [024 ; 135 ; 246] for projection (1, 2). (d) Index matrix [000 ; 111 ; 222] for projection (1,0).

(2) Projection selection for initial iteration: Given a P × Q image and the total number of bins B on projection direction, projections (p _i , q _i ) must satisfy the acquisition requirement for Equation (4), ( Q - 1 ) | p i | + ( P - 1 ) | q i | + 1 ≤ B(4) Based on Mojette transform theory, noise is propagated from bins to pixels in the reconstruction process, the first reconstructed pixels only include the error from the first projection. Aiming to avoid the dispersion of noise from the start, the projection sampling the maximal number of single bin-pixel correspondence in initial iteration is prioritized to picked out. Because of corner-based Mojette inversion algorithm (CBI), the maximal number of reconstructed pixels in one corner of the image is |p| · |q|.

Suppose an image corner can be reconstructed first from one of the projections (5,4), (4,3) and (5,1). Fig. 6 shows the reconstructed results via single bin-pixel correspondence. Fig. 6(a) shows the rays passing through the discrete object in projection (5,4). Based on Mojette transform theory, the centers of pixels in the 5 × 4 shadow region are only crossed by one ray (single bin-pixel correspondence), i.e., the ray only passes through the center of one pixel in the image corner, e.g., the rays 1, 2, 3 pass through the center of pixels 1,2,3 respectively. These pixels can be directly reconstructed from projections and the total number is |p| · |q|. Likewise, Fig. 6(b) and 6(c) show 4 × 3 and 5 × 1 single bin-pixel correspondences in projections (4,3) and (5,1) respectively. Therefore, we give priority to the projection (5,4) with the largest |p| · |q| for the initial iteration, which ensures the largest number of reconstructed pixels.

OPEN IN VIEWER

Fig.6

Diagram of single bin-pixel correspondence for projections (5,4), (4,3), (5,1). (a) The centers of pixels in the (5,4) shadow region are only crossed by one ray and 5 × 4 single bin-pixel correspondence in the projection. (b) 4 × 3 single bin-pixel correspondence in projection (4,3). (c) 5 × 1 single bin-pixel correspondence in projection (5,1).

(3) Index matrix and projection data update: The matrix element Proj_P×Q (k, l) that corresponds to the coordinate (k, l) of the reconstructed pixel is removed from all the index matrices Proj_P×Q. Meanwhile, the reconstructed pixel is subtracted from the bins that include this pixel on projections.

(4) Projection selection for following iterations: We need to select index matrix Proj _k with the maximum number of single bin-pixel correspondence among current index matrices Proj_P×Q. The projection corresponding to Proj _k is prioritized to be picked out for the next iteration. Therefore, we can maximize the number of single bin-pixel correspondences from the prioritized selected projection, and then the maximum number of reconstructed pixels can be obtained for the following iterations.

(5) Loop until the reconstruction done: Finally, we repeat steps 3 to 5 until all the pixels of slice image are reconstructed.

Given that all the pixels have been completely solved in the above reconstruction process, the priority-based subset must satisfy Katz theorem, i.e., ∑ i = 1 I * p i ≥ P or ∑ i = 1 I * q i ≥ Q , where I^* is the total number of all the selected projections. When ∑ i = 1 I * p i = P or ∑ i = 1 I * q i = Q is satisfied, PBI can reconstruct the image from the minimum set of projections. This is a kind of ideal situation in the case of the least number of projections to achieve the best quality of images. When ∑ i = 1 I * p i > P or ∑ i = 1 I * q i > Q is satisfied, the number of projections may not be the least and include more redundancy and alternative projections to reconstruct the image. The proposed method chooses the rational projection which can obtain the largest number of reconstructed pixels in each iteration. Therefore, it requires the minimum number of iterations for an image reconstruction and correspondingly reduces the noise accumulation. In a word, in the above two cases, the proposed PBI has achieved the best quality of images with the least number of projections.

3.2 An example of accumulation noise analysed between different subsets of projections

In this section, an example is given to compare the priority-based reconstruction path with random reconstruction paths. The times of error propagation from bins to pixels, also called as a noise multiplication factor, is calculated according to reconstruction paths.

Firstly, we take Fig. 7 as an example to design a priority-based reconstruction path and random reconstruction paths. The detector array includes 136 equi-spaced elements and the size of reconstructed image is 16 × 16. According to Equation (4), all the projections satisfying the inequalities are gathered together and expressed as follows: {(0, 1), (1, 1), (1, 2), (1, 3), (1, 4), (1, 5), (2, 1), (2, 3), (2, 5), (3, 1), (3, 2), (3, 4), (3, 5), (4, 1), (4, 3), (4, 5), (-1, 1), (-1, 2), (-1, 3), (-1, 4), (-1, 5), (-2, 1), (-2, 3), (-2, 5), (-3, 1), (-3, 2), (-3, 4), (-3, 5), (-4, 1), (-4, 3), (-4, 5)}. Based on step 2 in Section 3.1, the largest number of single bin-pixel correspondence exists on projections (5,4) or (-5,4) for the first iteration. Without loss of generality, here projection (5,4) in initial iteration is selected to reconstruct the maximum number of pixels shown in upper left corner of Fig. 7(a). Next, the proposed method purposefully maps the bins on selected projections to the corresponding pixels by Mojette inversion technique. Then, we update index matrix and projection data. For the next iteration, the projection (-5,4) recovering the maximum number of pixels would be picked out shown in upper right corner of Fig. 7(a), and so on. Finally, the projections (5, 4), (-5, 4), (5, 1), (5, 4), (-5, 4), (1, 5), (5, 4) are achieved, where the projection (5,4) is repeated three times in iterations. To compare with the priority-based subset of projections, the same number of randomly selected projections (5, 1), (-4, 1), (3, 1), (4, 1) makes up the other subset of random projections.

OPEN IN VIEWER

Fig.7

An example for reconstruction paths with the selected projections. (a) Reconstruction path with the purposefully selected projections corresponds to seven iterations. The error accumulated on pixel 7 reaches seventeen times. (b) Reconstruction path with the randomly selected projections corresponds to one hundred and nine iterations. The error accumulated on pixel 109 is far greater than one hundred and nine times.

Then, we calculate a noise multiplication factor according to reconstruction paths in Fig. 7. The lines in Fig. 7(a) and 7(b) represent different reconstruction paths with different subsets of projections. The solid and bold dashed lines in Fig. 7(a) are the priority-based reconstruction paths. Number in each gridding unit is used to record the number of iterations when the pixels are reconstructed and also represents the reconstructed sequences of pixels, e.g., pixel 1 in this gridding image, labeled by number 1, can be reconstructed in the first iteration. The section takes the solid line as an example, which passes through pixels 1, 3, 4, 5, 6 and 7. Thus, up to seven iterations are required to reconstruct pixel 7. Specifically, the solution along reconstruction paths is roughly described as follows: firstly, projections (5,4) and (-5,4) are used to solve pixels 1 and 2 in the image corners. And these pixels are subtracted from projection (5,1) to get pixel 3, as the lines on projection (5,1) only pass through the center of pixels 1, 3 and 2. Therefore, the error on pixel 3 is accumulated three times when pixel 3 is solved by using projections (5, 4), (-5, 4), (5, 1).

We repeat the above process for the other pixels in this region where is crossed by the solid lines, until the corresponding pixels are completely reconstructed. The noise is constantly accumulated when the pixels are reconstructed from the corner of the image to the center. The larger values in this gridding image, the bigger error accumulated on the reconstructed pixels. For example, the error accumulated on pixel 7 reaches seventeen times in Fig. 7(a). The bold dashed lines are fully symmetrical with the solid lines. Similar to the reconstruction path in the solid line, the error on pixel 7 is also propagated seventeen times in the blue line.

The lines in Fig. 7(b) describe random reconstruction path with the same number of projections as that of the priority-based reconstruction path in Fig. 7(a). Numbers in Fig. 7(b) are obtained from reconstructions with the randomly selected projections. Fig. 7(b) shows that the largest pixel in gridding image is one hundred and nine. Based on the above calculation method, the noise multiplication factor on pixel 109 is far greater than one hundred and nine times. Clearly, the reconstructed sequences of pixels in Fig. 7(b) are larger than those of Fig. 7(a).

In conclusion, the number of iterations based on the priority-based reconstruction path is far fewer than that of random reconstruction paths when recovering the image. In general, noise is uniformly distributed. Thus, smaller error on the reconstructed image will be accumulated by using the priority-based reconstruction path. Therefore, the proposed method using a minimum number of iterations effectively reduces the noise influence and obtains more accurate reconstruction results. In following section, analysis of anti-noise performance based on the reconstruction path isgiven.

4 Experiments and results

To demonstrate the effectiveness of accumulation noise reduction using the proposed method, experiments for the original hand and head image are implemented from different subsets of projections with a random uniform noise. An equi-spatial detector array is set that includes 1009 elements. The 64 × 64 images are simulated by Matlab 2012a. Based on discrete image geometry, the acquisition is performed along parallel rays to obtain exact and noise-free Mojette projections M. The prioritized sparse projections are non-uniformly distributed over [0, π]. The reconstruction process of PBI algorithm is described as follows.

Without loss of generality, we establish a kind of priority list for the reconstructed pixels in each iteration. The priority-based subset of projections with minimum accumulation noise is subjectively selected from all projections using the proposed method. The noise added to projection data comes from a random uniform distribution in interval [-0.1 % M_max, 0.1 % M_max], where M_max = max{ M_p,q (b), (b, p, q) ∈ M } [14]. An image is completely reconstructed from these noisy projections in 38 iterations. Specifically, the location of the reconstructed pixels in each iteration is distributed in Fig. 8, where first image shows the pixels in first iteration and second image expresses the pixels in second iteration, and so on. An image is iteratively reconstructed pixels by pixels. Similarly, intermediate results from the priority-based subset of projections in each iteration are shown in Fig. 9.

OPEN IN VIEWER

Fig.8

The location of the reconstructed pixels in each iteration.

OPEN IN VIEWER

Fig.9

The reconstructed results in each iteration.

To compare with traditional CBI, the other two subsets of projections are randomly selected from all projections. We reconstruct the original hand phantom from the first random subset and second random subset. Table 1 shows specific projections from three different subsets and the number of iterations. Three different subsets have the same number of projections but the main differences are the number of iterations. When reconstructions for the original hand image are done, there are 38 iterations for the priority-based subset, 137 iterations for the first random subset, and 1061 iterations for the second random subset. The computed times for the priority and random projection schemes are also shown in Table 1.

Fig. 10 shows the results reconstructed by PBI and the comparison results using traditional CBI, and the difference images between the reconstructed results and original image so the reconstruction noise can be more easily seen. Subfigure (a) is the original hand image. Subfigure (b) is the reconstructed result from the priority-based projection subset with the noise [-0.1 % M_max, 0.1 % M_max], which indicates that the PBI reconstruction is with good quality in the case of small noise, and especially the intensity of reconstructed pixels in the image corner is almost equal to the original pixels. Increasing the noise range to [-0.5 % M_max, 0.5 % M_max], the reconstruction result shown in subfigure (c) includes more errors, i.e., more accumulated noise is introduced into pixels of the center region. Nevertheless, traditional CBI fails to reconstruct an image from the random subsets even in the [-0.1 % M_max, 0.1 % M_max] noise environment, as shown in subfigure (d) and (e). Thus, subfigures (b) and (c) obtained from the proposed PBI algorithm have higher quality than the subfigures (d) and (e) from traditional CBI algorithm. Similarly, we have extended the proposed PBI algorithm to a complex head phantom with more grey values, its reconstructed results for the priority and random projection schemes are shown in Fig. 11. Image quality metrics between a distorted image and a reference image are computed by the Structural SIMilarity (SSIM) criterion [22]. The closer the SSIM value is to 1, the better the reconstruction quality is. Moreover, Peak Signal to Noise Ratio (PSNR) is used as an objective criterion to denoise. The SSIM and PSNR scores of the reconstructed results show that, PBI algorithm can yield higher noise reduction than traditional CBI in sparse-view imaging, even if in the larger [-0.5 % M_max, 0.5 % M_max] noise environment.

OPEN IN VIEWER

Fig.10

Comparison results. (a) Original hand image. (b) Reconstructed result by PBI in [-0.1 % M_max, 0.1 % M_max] noise environment. (c) Reconstructed result by PBI in [-0.5 % M_max, 0.5 % M_max] noise environment. (d) and (e) Reconstructed results using traditional CBI with two different random schemes. (f), (g), (h) and (i) Difference images between original hand image and their reconstructed results.

OPEN IN VIEWER

Fig.11

Comparison results. (a) Original head image. (b) Reconstructed result by PBI in [-0.1 % M_max, 0.1 % M_max] noise environment. (c) Reconstructed result by PBI in [-0.5 % M_max, 0.5 % M_max] noise environment. (d) and (e) Reconstructed results using traditional CBI with two different random schemes. (f), (g), (h) and (i) Difference images between original head image and their reconstructed results.

The differences from each original image are shown in subfigures (f), (g), (h) and (i) of Figs. 10 and Fig. 11, respectively. The accumulated error in reconstruction results is at different levels in different cases. The minimum noise accumulations in subfigures (f) and (g) are from the priority-based subset of projections because that image reconstruction has the minimal number of iterations. When the number of iterations increased from 38 to 137, much more error in subfigure (h) is accumulated in the center region of the slice image when using traditional CBI. With 1061 iterations from the second random subset of projections, the accumulation noise in subfigure (i) is very large and its reconstruction result (e) can hardly be identified. In conclusion, because of the minimal number of iterations from the priority projection scheme, the reconstructed images with the minimum noise accumulation are obtained.

To see the evolution of the accumulated noises for the priority and random projection schemes, take the hand image with the noise in the interval [-0.1 % M_max, 0.1 % M_max] as an example, thereconstruction noises per iteration in subfigures (b), (d) and (e) are plotted as the line (b), the line (d) and the line (e) in Fig. 12. The graph clearly shows that the line (b) stops at the 38th iteration which is the maximum number of steps for image reconstruction with the minimum noise accumulation. The line (d) ends at the 137th iteration, and its accumulation noise becomes larger than for the line (b). In addition, at 1061 iterations along the line (e), the accumulation noise increases too much to reconstruct the image. Therefore, the noise accumulation follows different paths, and the two different random projection schemes by CBI cannot get lower noise and stable results even through earlier termination.

OPEN IN VIEWER

Fig.12

The noise accumulation per iteration for the priority and random projection schemes in the [-0.1 % M_max, 0.1 % M_max] noise environment. The noise accumulations on the hand reconstruction results are plotted as the line (b), the line (d) and the line (e).

Experimental results demonstrate that noises gradually accumulate from edge to center in Mojette reconstruction process. The number of iterations directly determines the quality of reconstructed images. Therefore, PBI with the minimum number of iterations effectively suppresses the accumulated noise without increasing the number of projections. To clearly see the differences between the PBI and CBI, we did not apply any methods to eliminate the noise on projections. PBI substantially improves the robustness of exact CBI through qualitative and quantitative analysis. Based on PBI, we will combine the statistical method to further reduce noise and apply Mojette transform in practice as soon as possible.

Traditional CBI algorithm can guarantee the accuracy of Mojette reconstruction in the absence of noise. However, its reconstruction strategy has poor performance with noisy projections. By contrast, the proposed PBI algorithm reduces its noise sensitivity and can reconstruct an image from sparse noisy projections. In our future work, we will combine the statistical method with the PBI algorithm to further reduce noise. Finally, we will apply the Mojette transform in practice as soon as possible.