Cite this Article
Liu Hua-Wei, Zheng Shu, Zhou Huai-Chun†. An iterative virtual projection method to improve the reconstruction performance for ill-posed emission tomographic problems. Chinese Physics B, 24(10): 104204  
Permissions
An iterative virtual projection method to improve the reconstruction performance for ill-posed emission tomographic problems
Liu Hua-Wei, Zheng Shu, Zhou Huai-Chun†
Key Laboratory for Thermal Science and Power Engineering of Ministry of Education, Department of Thermal Engineering, Tsinghua University, Beijing 100084, China

†Corresponding author. E-mail: hczh@mail.tsinghua.edu.cn

*Project supported by the China National Funds for Distinguished Young Scientists of National Natural Science Foundation of China (Grant No. 51025622), the National Natural Science Foundation of China (Grant No. 51406095), and the 100 Top Talents Program of Tsinghua University, Beijing, China (2011).

Abstract

In order to improve the reconstruction performance for ill-posed emission tomographic problems with limited projections, a generalized interpolation method is proposed in this paper, in which the virtual lines of projection are fabricated from, but not linearly dependent on, the measured projections. The method is called the virtual projection (VP) method. Also, an iterative correction method for the integral lengths is proposed to reduce the error brought about by the virtual lines of projection. The combination of the two methods is called the iterative virtual projection (IVP) method. Based on a scheme of equilateral triangle plane meshes and a six asymmetrically arranged detection system, numerical simulations and experimental verification are conducted. Simulation results obtained by using a non-negative linear least squares method, without any other constraints or regularization, demonstrate that the VP method can gradually reduce the reconstruction error and converges to the desired one by fabricating additional effective projections. When the mean square deviation of normal error superimposed on the simulated measured projections is smaller than 0.03, i.e., the signal-to-noise ratio (SNR) for the measured projections is higher than 30.4, the IVP method can further reduce the reconstruction error reached by the VP method apparently. In addition, as the regularization matrix in the Tikhonov regularization method is updated by an iterative correction process similar to the IVP method presented in this paper, or the Tikhonov regularization method is used in the IVP method, good improvement is achieved.

PACS: 42.30.Wb; 43.28.+h
Keyword: ill-posed problem; emission tomography; limited projections; Tikhonov regularization
1. Introduction

Source distribution inside an object directly depends on relevant parameter distribution. For example, the radiative source inside a high temperature enclosure such as a combustor is determined by the distribution of the temperature and radiative properties, such as spectral absorption and scattering coefficients and/or scattering phase function, etc.[1] In many cases, unknown parameters are essentially calculated from reconstructed source distributions, like radiation temperature measurement, [1] multi-wavelength thermometry, [2] and temperature pattern reconstruction by using temperature dependence of permittivity.[3] Computerized tomography (CT) imaging is used to reconstruct the source distribution in the interior of an object after collection of externally acquired projection data in different directions.[4, 5] Nevertheless, emission tomography is always a highly limited data set problem, due to the limited availability of projection angles, the coarse sampling, [6] and other practical constraints brought about by imaging hardware, scanning geometry or ionizing radiation exposure; [7] or because of limited projection data obtained from a few-view CT system to meet the requirement for high temporal resolution.[8] The challenge is how to reconstruct the source distribution by using under-sampled projection data, which is obviously mathematically ill-posed and additional knowledge is needed to help achieve reasonable reconstruction results.

Generally speaking, this kind of ill-posed problem is mainly caused by the contradiction between the lack of effective measurement information and the requirement of high resolution. In order to solve this contradiction, researches have been carried out. Huang used a simple high-resolution stereoscopic system to receive sufficient projection data to reconstruct a small number of two-dimensional (2D) unknowns, [9] thus obtaining as fine as possible reconstruction results under this experimental condition. Ishino and Ohiwa wanted to obtain three-dimensional (3D) chemiluminescence distribution with high spatial resolution by using computerized tomographic reconstruction techniques. So they used projection data taken from forty directions with forty small lenses, thus making reconstruction results with high spatial resolution feasible.[10]

Furthermore, advanced image reconstruction algorithms are also effective tools to deal with such ill-posed problems, including Tikhonov regularization (TR), [6, 11] minimum cross-entropy (MCE), [8] filtered back projection (FBP), [12] maximum entropy (ME), [13– 15] total variation (TV), [7, 16– 18] pre-iteration Landweber, [19] and many others.[20, 21] The TR method is to seek for a solution defined as the minimization of the weighted combination of the residual norm of the matrix equations and a two-norm parameter which utilizes different intrinsic smoothing principles and controls the smoothness of the regularized solution.[6] The MCE method is to maximize the likelihood function, which is equivalent to minimizing the cross-entropy between two non-negative vectors. Using regularization to deal with instability in inverse problems and a prior knowledge about the image to be recovered to reduce the ill-posedness, the problems can be solved.[8] The FBP method consists of a filtering process in which projection data are expanded by an orthogonal function system and a backprojection process. By minimizing an evaluation function for smoothness through using an evaluation matrix in the FBP calculation process, regularized solution data will be obtained.[12] In general, the reconstruction algorithms mentioned above always utilize additional prior knowledge, of which smoothness treatment is a common kind.[22] Smoothness treatment is based on the continuity within one object, constraining the difference between adjacent source terms and providing additional information.

In our previous work, for a reconstruction system with six asymmetric projection directions and equilateral triangle plane meshes, a random interpolation method implemented between adjacent projections is used to fabricate virtual projection lines lineally independent of the real projection lines, which is another manifestation mechanism of the smoothness of the interior of an object rather than those embedded in the reconstruction algorithms such as TR mentioned above. It is qualitatively demonstrated that this method with using virtual projection lines can effectively reduce the degree of ill-posedness and improve the solvability of the inverse problem.[23] If a similar method can form virtual projections with better reconstruction ability, a new way to solve this kind of ill-posed problem will be found.

In this paper, with a six asymmetric projection detection system and equilateral triangle plane meshes similar to our previous work, a more generalized interpolation method to fabricate independent, virtual projection lines for every single measured projection line is proposed for much less errors and better reconstruction results. In addition, an iterative correction method for integral lengths according to the reconstruction results from the last iteration is proposed to gradually reduce the errors brought about by virtual projections. The reconstruction problem, interpolation method, and iterative correction method will be described first. Then, numerical simulations and experimental verification will be conducted to test the reconstruction ability with a non-negative linear least squares method and the Tikhonov regularization method. Finally, the corresponding conclusions will be drawn.

2. Methodology
2.1. Object studied

The projection data of a cross section are one-dimensional (1D), whereas the parameter distribution in a cross section is 2D. If we want to reconstruct source distribution with high spatial resolution when the size of the 2D meshes equals that of projection data, projection data from several directions will not be sufficient, obviously. Like our previous work, [23] as shown in Fig.  1, the square area of study is divided into equilateral triangle plane meshes whose identical edge lengths equal twice the initial inherent width of the line of sight in the projection. Therefore, measurements in each projection are independent of each other and can obtain the best use to achieve the maximum spatial resolution. For an equilateral triangle, the three altitudes are equivalent, and the three median lines are equivalent, too. As illustrated in Fig.  1, corresponding to the six characteristic lines, six asymmetric projection directions are arranged to capture projection data, ensuring that for each projection line, the integral path through a mesh element is one of the six characteristic lines and integral path lengths through the relevant mesh elements are identical for each projection, thus reducing reconstruction difficulty.

Fig.  1. Detection system with six asymmetric projections. The square area of study is divided into equilateral triangle plane meshes whose identical edge lengths equal twice the initial inherent width of the line of sight in the projections. Six asymmetric projection directions are arranged to capture projection data.

In this paper, the square area of study is 10.9  mm × 10.9  mm. The first three projections pass through altitudes of corresponding mesh elements, whereas the last three projections pass through median lines of corresponding mesh elements. The width of one measurement line in the projections is set to be 1/4.7  mm as illustrated in Fig.  1, so the initial number of measurement lines for the 10.9-mm projection is 52. The mesh number along the boundary of the square area is equal to that of the projection lines parallel to this boundary, so the number of vertical meshes in Fig.  1 is also 52. In consequence, the numbers of horizontal projection lines and meshes are both 30, due to the characteristics of equilateral triangles. Therefore, the total number of discrete meshes is 52× 30, i.e., 1560 meshes are to be reconstructed, whereas the total number of effective measurement lines in the six projections is only 246, which is significantly smaller than the total number of discrete meshes.

2.2. Forward problem

For a path S shown in Fig.  1, the projection intensity I of a line of sight detected from the projection direction 1 can be expressed in a matrix form:

or

For situations with noise, the noise level can be evaluated by SNR, which is defined as:

where σ ξ j represents the normal random error superimposed on the projections; a detailed description of the projection expression and noise level can be found in Refs.  [9] and [23].

In the present paper, five noise levels are considered, i.e., values of normal random noise σ equal 0.005, 0.01, 0.02, 0.03, and 0.05, and values of SNR are 45.54, 40.5, 34.24, 30.4, and 25.58, respectively.

2.3. Non-negative linear least squares method for reconstruction

Non-negativity constraints are of great importance for practical models. Unlike the algorithms mentioned in the introduction section, the non-negative linear least squares method imposes no additional constraints except that the solution should be non-negative, where the basic requirement is to form a physically meaningful solution. In this paper, the non-negative linear least squares method is used to solve Eq.  (1) and test the reconstruction ability of the method proposed here. For a system of linear equations as expressed in Eq.  (2), the solution obtained by the non-negative linear least squares method is expressed as described by Garda and Galias, [24] Chiao et al., [25] Bro and Jong:[26]

In our previous work, [23] the performance of the virtual lines of sight and its effects on the solution of the problem mentioned above have been examined; so in this paper, it is the intention to find a suitable generalized method to solve this problem. Fabrication of virtual projection lines is another manifestation of smoothness of the interior of an object, and no prior knowledge is necessary to stabilize or smooth the solution, so a non-negative linear least squares method is selected to examine the effect of the virtual projection lines on the reconstruction performance in this paper.

2.4. Virtual projection (VP) method

In our previous work, [23] virtual projection lines fabricated by randomly selecting mesh elements from the ones in corresponding positions of two adjacent projections along the projection direction have been used to test the potential reconstruction ability of virtual projections which has been demonstrated. However, this kind of virtual projection relies on the relationship between the positions of adjacent projections and thus brings about significant error worsening reconstruction results. Therefore, in order to develop a method with less error and a wider range of applications, a generalized interpolation method, called the virtual projection (VP) method, is described in this section.

As shown in Fig.  2, for each equilateral triangle mesh, there exist three adjacent meshes with similar source values due to the spatial continuity, except for some meshes near the boundaries, where the numbers of adjacent meshes are different. Each mesh of a projection path is randomly replaced by a corresponding adjacent mesh with a certain probability, and then a virtual projection line forms as shown in Fig.  2. For a virtual projection line formed in this method, it is assumed that the projection intensity does not have too much change, so the virtual projection intensity is represented by the measured one primarily.

Fig.  2. Illustration of random interpolation mesh for the virtual projection. (a) A measured projection along an integral line passing through 8 upwards and downwards meshes. (b) The meshes forming a virtual projection from the measured projection in panel  (a) after a random replacement process: the top first and seventh meshes are replaced by the corresponding left meshes; the second, the third, the fifth, and the eighth meshes are unchanged; the fourth mesh is replaced by the right mesh; and the sixth mesh is replaced by the upper mesh. (c) Three possible adjacent meshes for every upwards and downwards meshes in panels  (a) and (b)

The virtual projections are in the following form:

where j and j′ represent the measured projection and virtual projection fabricated from measured projection j, respectively, i′ represents the corresponding mesh after the replacement, and Δ Sj, i is the integral length of the j-th line path through the i-th mesh, which is unchanged in this step and will be corrected later in a further discussion. By comparison with Eq.  (1), only meshes in the line of sight are treated in this step. Adjacent meshes not involved in the measured projection are selected randomly to replace the meshes in the line of sight along the measured projection, thus forming mathematically linearly independent virtual projections. In this paper, the probability for a mesh to be replaced by any of the three adjacent meshes is a quarter, and there is still a probability of a quarter being unchanged. That is to say, the probability for a mesh to be replaced by any adjacent mesh or remaining unchanged is identical, and it is also true for the meshes near the boundaries where the identical probabilities need to undergo corresponding change.

As shown in Fig.  2, as a demonstration, for a measured projection, the integral line passes through 8 meshes which belong to two position types along the projection direction, upwards and downwards, corresponding to leftwards and rightwards in Fig.  1, respectively. For either position type, there are three possible adjacent positions, and there are four kinds of possibilities combined with the original position. For the situation shown in Fig.  2, after the random replacement process, the top first and the seventh meshes are replaced by the corresponding left meshes; the second, the third, the fifth, and the eighth meshes are unchanged; the fourth mesh is replaced by the right mesh, and the sixth mesh is replaced by the upper mesh, thus forming a virtual projection shown in the figure. After taking similar steps for all the measured projections, the virtual lines of sight form a new matrix equation expressed as

Combining the measured and virtual projection matrix as shown in Eq.  (6b), the source distribution F can be reconstructed.

2.5. Iterative virtual projection (IVP) method

The VP method mentioned above can indeed increase the amount of “ measured” information, at the same time, randomly replacing meshes and roughly estimating the virtual projection intensity with the corresponding measured one bring some errors, since actually the source term in each mesh is not exactly equal to those in its adjacent meshes. This replacement is not favorable for the reconstruction procedure. Unlike the method in our previous work, [23] each virtual projection using the VP method mentioned above is based on a single measured projection, enabling it to undergo further correction operation. To reduce the error brought about by the VP method mentioned above, an iterative correction method for integral length based on the reconstructed source distribution is adopted in this paper. As shown in Fig.  3, the integral lengths are corrected according to the source term ratio between the initial mesh in the measured projection and the one in the corresponding position after replacement determined in the last iteration, which can be expressed as

where

Here, i represents the initial mesh in the measured projection, i′ represents the mesh in the corresponding position after the replacement, j′ is the virtual projection, k is the iteration number, and is the correction coefficient, which is used to compensate for the influence of replacement of the i-th mesh by the i′ -th mesh, according to the reconstruction results in the (k − 1)-th iteration. By comparison with Eq.  (5), the only difference is the correction coefficient which is determined by using the reconstructed distribution in the last iteration or uniform distribution for the first iteration when all coefficients should be 1. It is obvious that the closer to the real value the solution of the (k − 1)-th step, the less error the virtual projections bring in. Combined with measured projection in every iteration, the solution is shown to approach to the real distribution by iteration. The whole process to solve this ill-posed problem is shown in Fig.  4.

Fig.  3. Integral length correction. The integral lengths are corrected according to the source term ratio between the initial mesh in the measured projection and the one in the corresponding position after the replacement determined in the last iteration.

Fig.  4. Flow chart. Firstly, virtual projections will be fabricated after random replacement of meshes based on measured projections. Then the correction coefficients will be determined and used to reconstruct the 2D distribution. After each iteration, smoothing filtering is used to calculate the integral length correction coefficient for the next iteration.

As shown in Fig.  4, firstly, based on measured projections, virtual projections will be fabricated after random replacement of meshes. Then with the reconstructed distribution or uniform distribution for the first iteration, the correction coefficients will be determined and used to reduce the error brought about by virtual projections. The measured and virtual projections are combined and used to reconstruct the 2D distribution in each iteration. After each time of iteration, smoothing filtering is used to calculate the integral length correction coefficient for the next iteration. The smoothing filtering adopted in this paper is the average filtering which can be expressed as

In this equation, the plus and minus signs are dependent on the mesh orientation shown in Fig.  1, and for the leftwards mesh, the sign is selected to be plus, indicating the right mesh, and for the rightwards mesh, the sign is selected to be minus, indicating the left mesh.

After each iteration, a convergence criterion will be evaluated to judge whether the results are converged. In this paper, the convergence criterion can be expressed by using the following equation:

The equation demonstrates that the convergence criterion is based on the difference in result between the adjacent two iterations compared with the reconstruction results in the last iteration. If the difference mentioned above is smaller than the limit value ɛ or seems to be constant or fluctuates up and down around a certain value, the iteration process is thought to be converged. The method is called the iterative virtual projection (IVP) method.

Virtual projection mentioned above is based on the smoothness of the interior of an object which is inherent in an actual object and unrelated to the mesh elements and the number and angles of projections. Therefore, the VP and IVP methods do not depend on the regular mesh elements and can be generalized to irregular grids by using Eqs.  (5) and (7), while the only difference is the random replacement process, relying on the number and distributed situation of the adjacent mesh elements.

3. Simulation results and discussion
3.1. Cases studied, measured projection intensities and virtual projections

In this paper, four cases with different source distributions are considered to test the reconstruction ability of the reconstruction method proposed in this paper. As shown in Fig.  5(a), the assumed source distribution in Case A is composed of four identical squares and peaks. As shown in Fig.  5(b), the assumed source distribution in Case B is composed of three superposed peaks. The assumed source terms in Case C shown in Fig.  5(c) are distributed into rings, and in Case D as shown in Fig.  5(d), the assumed source distribution is composed of half a ring and a peak.[6] Numerical simulations are conducted with these four initial assumed source distributions. To compare reconstruction results in different cases, average relative error is calculated according to the following equation:[27]

where Rrecon, i represents the reconstructed source term for the i-th mesh, Rassum, i represents the assumed source term for the i-th mesh, and N is the total mesh number. WR with this equation can describe the total error between the reconstructed and assumed distributions.

Fig.  5. Assumed source distributions. (a) Case A: four identical squares and peaks. (b) Case B: three superposed peaks. (c) Case C: a ring. (d) Case D: a half of a ring and a peak.

With the assumed source distribution, we can calculate the 246 measured projection intensities. The measured projection intensities for the four cases by using Eq.  (1) through using a detection system in experiments are shown in Fig.  6. So the problem is to reconstruct the unknown source terms in 1560 meshes shown in Fig.  1 from projection intensities shown in Fig.  6.

Fig.  6. Measured projection intensities using Eq.  (1).

Virtual projections for the four cases with the VP method without iterative correction are shown in Fig.  7. It can be seen that the virtual projections are similar to those measured projections shown in Fig.  6. Actually, the virtual projections are just deduced from the measured projections without “ creating” any different projections. The difference between the measured and virtual projections is only the projection paths. The virtual projections are related to those “ virtual paths” which are selected randomly from among the meshes adjacent to the real paths of the measured projections.

Fig.  7. Virtual projection intensities (without length correction). The virtual projections are just deduced from the measured projections without “ creating” any different projections. The difference between the measured and virtual projections is only the projection paths.

Reconstruction errors varied with interpolation time for Case A and Case B with the VP method are shown in Fig.  8(a). The figure shows that as the interpolation time increases, the reconstruction results are improved and the error converges to a limit value near 5% for Case A. The error converges to a limit value near 2.5% for Case B. The reconstruction results for Case A and Case B without the VP method or with the VP method with 20 times virtual projections are also shown in Figs.  8(b)– 8(e). It is obvious that the VP method can improve the reconstruction results obtained by the non-negative linear least squares method. Nevertheless, considering the total reconstruction error and the calculated quality, even the reconstruction with the VP method only is satisfactory, there is still a lot of room for improvement, and then the IVP method is examined.

Fig.  8. Results with different interpolation times. (a) Reconstruction errors varying with interpolation time without iterative correction method. (b)– (e) The reconstruction results for Case A and Case B without VP method or with 20 times interpolation projections without iterative correction.

3.2. Effects of virtual projections with iterative correction on the reconstructed results

The problem is to reconstruct 1560 unknowns with the intensities in 246 projection lines. Each time virtual projection makes 246 additional projections, making the reconstruction problem realizable. With the IVP method, for the cases mentioned above, the source distributions are reconstructed. Without any de-noising techniques to reduce noise and increase the reliability, [28] the five noisy levels mentioned above and a noise-free situation are taken into account to test the effects of virtual projections on the reconstructed results by this method.

Using an ordinary personal computer with an Intel(R) Core(TM) i5-2320 CPU and 6  GB RAM, it takes tens of seconds to fabricate and correct 15 times, i.e., 3690 virtual projections in each iteration, demonstrating that fabricating virtual projections will not require much extra running time.

The reconstruction results for Case A when the mean square deviation, σ , equals 0 or 0.01 without any iteration or with 30 iterations are shown in Fig.  9. It is demonstrated that the iteration process can reduce the reconstruction errors apparently. The situation is the same when dealing with Cases B, C, and D shown in Figs.  10, 11, and 12, respectively.

Fig.  9. The comparison among reconstruction results (15 times virtual lines, Case A). (a) Without iteration process (σ = 0), (b) 30 iterations (σ = 0), (c) without iteration process (σ = 0.01), and (d) 30 iterations (σ = 0.01).

Fig.  10. Comparison among reconstruction results (15 times virtual lines, Case B). (a) Without iteration process (σ = 0), (b) 30 iterations (σ = 0), (c) without iteration process (σ = 0.01), and (d) 30 iterations (σ = 0.01).

Fig.  11. Comparison among reconstruction results (15 times virtual lines, Case C). (a) Without iteration process (σ = 0), (b) 30 iterations (σ = 0), (c) without iteration process (σ = 0.01), and (d) 30 iterations (σ = 0.01).

Fig.  12. Comparison among reconstruction results (15 times virtual lines, Case D). (a) Without iteration process (σ = 0), (b) 30 iterations (σ = 0), (c) without iteration process (σ = 0.01), (d) 30 iterations (σ = 0.01).

Figure  13 shows the reconstruction results for the first two cases with 15 times virtual projections after 30 iterations. As shown in Fig.  13(a), the tendency of the reconstructed distribution for Case A is in good agreement with the assumed one except that there are slight fluctuations in areas near the boundaries and between two adjacent peaks. Considering the limited measurement information, the reconstruction results are acceptable. Figure  13(b) shows the relative error distribution between the reconstructed and assumed distributions. Large errors occur in areas near the boundaries and between two adjacent peaks because of the slight fluctuations and the smaller source terms in the corresponding areas, whereas in areas with a bigger source term, the corresponding reconstructed values are in good agreement with the assumed ones. Errors in the four peak areas are slightly bigger than those in the adjacent areas due to the poor smoothness in peak areas. The reconstructed source and relative error distributions for Case B are shown in Figs.  13(c) and 13(d), respectively. As with Case A, there appears the tendency towards good reconstruction except for some slight fluctuations near the boundaries, where there are bigger relative errors with smaller source terms. Generally speaking, the reconstructed results are acceptable, demonstrating that the IVP method with iterative correction for integral lengths has good reconstruction ability.

Fig.  13. Reconstruction source term and relative error distribution (15 times virtual projections, 30 iterations, σ = 0) (a) and (c) Reconstruction distribution. (b) and (d) Reconstruction error distribution. (a) Reconstructed distribution (Case A), (b) reconstructed error (Case A), (c) reconstructed distribution (Case B), (d) reconstructed error (Case B).

3.3. Variation of the errors with the time of virtual projection

Reconstruction errors using Eq.  (10) and convergence criteria using Eq.  (9) for different noise levels for the first two cases are shown in Fig.  14. As the iteration number increases, the relative error in the reconstruction results decreases, indicating an improvement in the reconstruction performance. As shown in Fig.  14(a), for Case A, as the iteration number increases with different virtual line times, the tendencies of the results are similar except for the initial values, which are obtained in the first iteration with uniform assumed distribution in which all the correction coefficients are 1. Three situations with three virtual line times have the same tendency and similar convergence values, which are between 0.02 and 0.025 and are significantly smaller than the corresponding initial values in the first iteration obtained by the VP method, calculated without iterative correction for integral lengths. The reconstruction results with different virtual projection times for Case B are shown in Fig.  14(c). Results with 10 times virtual projections show the worsening of the situation with a big iteration number. When calculated with 15 and 20 times virtual projections, the situations improve and the trend is clear, indicating that the former situation is due to a lack of interpolation time. For this case, the range of value variation ratio is wider than that for Case A, and correction coefficients in areas where peaks are superposed are susceptible to bias, making 10 times virtual projections not enough to achieve stable results available to realize accurate correction. Therefore, the times of virtual projections are chosen to be 15 in the subsequent calculation. Like Case A, the latter two situations have the same tendency and similar convergence values, which are between 0.005 and 0.01, significantly smaller than the corresponding initial values calculated by the VP method shown in Fig.  8(a). Figures  14(b) and 14(d) show the convergence limits for two cases. It is obvious that the convergence limit decreases as the reconstruction results improve. In this paper, it is intended to show the reconstruction ability with iteration process, so in the following sections, reconstruction results after a certain number of iterations will be discussed.

Fig.  14. Iteration reconstruction results (σ = 0), showing (a) and (c) reconstruction errors using Eq.  (10). (b) and (d) Convergence limits using Eq.  (9). (a) Reconstructed distribution (Case A), (b) convergence criterion (Case A), (c) reconstructed distribution (Case B), (d) convergence criterion (Case B).

3.4. Comparison and analysis of noise immunity

To evaluate the reconstruction ability of the IVP method proposed in this paper, a comparison is conducted. The reconstructions were in three situations where no virtual projections are fabricated, 15 times virtual projections are fabricated with the VP method, and 15 times virtual projections are fabricated with the IVP method with 30 times of iterative correction for integral lengths. Reconstruction errors calculated using Eq.  (10) are shown in Table  1.

Table 1. Reconstruction errors for different situations with VP method and IVP method.

As shown in Table  1, the VP method proposed in this paper can increase the amount of effective information and the iterative virtual projection method can indeed improve the quality of the additional information. As the noise level increases, the improvement effect decreases, because the iterative correction method of integral lengths can reduce the error brought about by virtual projections, but cannot reduce the error existing in the measured projections. On the contrary, the error of the measured projections can reduce the reconstruction quality in each iteration, thus making the correction coefficient inaccurate and reducing the auxiliary force of the iteration process for integral lengths. When σ , the mean square deviation of noise is smaller than 0.03, i.e., the signal-to-noise ratio is bigger than 30.4, the iterative virtual projection method has a significant effect. When σ increases to 0.03, the large fluctuations of the reconstruction results make the iteration process invalid, demonstrating the applicable scope for the iteration process proposed in this paper. In that situation, the VP method is still able to generate acceptable additional information to help to finish the reconstruction process. Besides, reconstruction errors with 15 times virtual projections with the method in our previous work[23] are also given in this table, and the significant improvement in situations with low noise demonstrates that the iterative virtual projection method proposed in this paper has made great progress in terms of accuracy.

4. Effect on TIKHONOV regularization method
4.1. Updating the regularization matrix by iterative correction method

At first, the iterative correction method can be integrated into the Tikhonov regularization (TR) method through updating its regularization matrix D.[22] In each iteration, the regularization matrix is updated according to reconstruction results from the last iteration, which can be expressed as

where k is the iteration number, f is the solution of the last iteration, and D is the regularization matrix expressed as

where n is the number of adjacent meshes for mesh i, which is 3 for most meshes except some meshes near the boundary, and j1– jn represent the adjacent meshes, respectively, N is the total number of meshes, which is 1560 in this paper. Except elements defined in the equation, the other elements of matrix D are all 0. Reconstruction errors obtained by using Eq.  (10) for Case D vary with iteration number when σ equals 0 or 0.05 as shown in Fig.  15. This figure shows that the iterative correction method for the regularization matrix can gradually improve reconstruction results obtained by the TR method, even under big noise levels.

Fig.  15. Variations of reconstruction error for the TR method with iterative number (α = 0.1, Case D). In each iteration, the regularization matrix is updated according to reconstruction results from the last iteration as expressed in Eq.  (11).

4.2. TR method with iterative virtual projections

In addition, results in three situations are compared. In the first situation, the Tikhonov regularization method is used with a best regularization coefficient α with one significant digit, which makes the reconstruction error determined by Eq.  (10) minimal. In the second situation, the TR method is combined with 15 times virtual projections without iterative correction (TR + VP), and in the last situation, the TR method is combined with 15 times virtual projections after 10 times iterative correction (TR + IVP). All the results are shown in Fig.  16. Besides, the regularization coefficients are also shown in the figure.

Fig.  16. Comparison among reconstruction errors. In each figure, in the first situation, the Tikhonov regularization method is used with a best regularization coefficient α with one significant digit which makes the reconstruction error determined by Eq.  (10) minimal. In the second situation, the TR method is combined with 15 times virtual projections without iterative correction (TR + VP), and in the last situation, the TR method is combined with 15 times virtual projections after 10 times iterative correction (TR + IVP). (a) Case A, (b) Case B, (c) Case C, and (d) Case D.

Using the ordinary personal computer mentioned above, it takes tens of seconds to fabricate and correct 3690 virtual projections in each iteration and a few seconds to solve Eq.  (6b) with the TR method. The total time of only a few minutes has shown good real-time performance in practical application.

From Fig.  16, it can be seen that the results obtained by the TR + IVP method for the four cases are best. With a regularization coefficient with one significant digit chosen from a certain range originating from that for the TR method, which is maybe not the best in the whole scope as shown in Fig.  16, the method composed of the TR and VP has similar results to the TR method, while the method composed of TR and IVP has better results in all simulation points, especially in Cases C and D. This is because the last two distributions are not so smooth as the former two, while for the TR method the smoother the distribution, the better the results will be.

Results show that the IVP method proposed in this paper can improve reconstruction results for algorithms like the non-negative linear least squares method and the Tikhonov regularization method to deal efficiently with ill-posed emission tomographic reconstruction problems.

5. Experimental verification

Experiment is conducted to verify the effectiveness of the proposed method. Figure  17(a) shows a photograph of the experimental system composed of six CCD devices. As shown in this figure, an ethylene laminar diffusion flame is generated with 194-mL/min ethylene and 240-L/min air with a co-flow laminar diffusion flame burner reported by Snelling et al.[2] Six CCD devices are asymmetrically arranged in accordance with Fig.  1 to take radiation images.

Fig.  17. Experimental set-up and detected intensity for a cross section. (a) Photograph of experimental set-up and (b) detected intensity for a cross section (R channel).

Figure  17(b) shows the distribution of detected intensities (R channel) for a cross section with a height of 35  mm above the burner, composed of projection data collated into 52  pixels for the first three CCDs or 30  pixels for the last three CCDs. As shown in Fig.  17(b), distribution of detection intensity for each CCD is not symmetrical, indicating that the flame is not completely axisymmetric. Reconstruction results with the non-negative least squares method only, VP method, TR + VP, and TR + IVP methods are shown in Fig.  18. As shown in the figure, the non-negative linear least squares method only cannot reconstruct acceptable results, while the other three methods can. Tendencies of reconstruction distributions for the last three methods are similar. For the last three methods, reconstructed radiation source terms are distributed into rings, while the distribution of positions of peaks is discontinuous. For the ethylene laminar diffusion flame studied in this paper, input flows of ethylene and air are axisymmetric. Therefore, ideally, source term distributions should also be axisymmetric. However, detection intensity asymmetry due to variances in burner geometry, noise and disturbance when experimenting, along with inadequate meshing, i.e., 52 × 30 meshes makes the discontinuity, which can also be found in Fig.  11 for Case C, where noise causes significant fluctuation in the area near the ring. Situation for the experiment is more complex and source terms change more rapidly, thus making the discontinuity more obvious. Comparing results shown in Fig.  18, the distributions of results with the TR method shown in Figs.  18(c) and 18(d) are smoother and much closer to rings. Table  2 shows the mean square errors of results obtained with these four methods. For results obtained with the last two methods, the smoothnesses and scopes of numerical value are similar, while the mean square error for TR + IVP is 1.165× 1011, which is larger than 1.137× 1011 for TR + VP, demonstrating that the trend of results obtained with Tikhonov regularization and iterative virtual projection methods is clearer, with peak distribution much closer to a thin ring as shown in Fig.  18, and indicating that the iterative correction process works. The virtual projection method and iterative correction process proposed in this paper are general methods, which do not depend on meshes and can be generalized to meshes of other types, as long as majority meshes have at least one adjacent mesh. For the case with more projection directions and finer meshes, reconstruction results based on corresponding meshes would be better. Through the reconstruction results introduced above, the effectiveness of the proposed method is verified.

Fig.  18. Reconstruction results. (a) Non-negative linear least squares method only, (b) VP (virtual projection) method. (c) TR (Tikhonov regularization) and VP method, (d) TR and IVP (iterative virtual projection) method.

Table 2. Mean square error of reconstruction results.
6. Conclusions

In order to reconstruct the 2D source distribution with high spatial resolution, an interpolation method to fabricate virtual lines of projections (VP method), and the method combined with an iterative correction method for integral lengths (IVP method), are proposed in this paper. Based on a scheme of equilateral triangle plane meshes and a six asymmetric arranged detection system, the numerical simulations and experimental verification are conducted with the non-negative linear least squares method and the Tikhonov regularization method. Simulation results show that the VP method can gradually reduce the reconstruction error and converge to the desired one by fabricating additional effective projections. When the mean square deviation of normal error superimposed on the simulated measured projections is smaller than 0.03, i.e., SNR for the measured projections is higher than 30.4, the IVP method can further reduce the reconstruction error reached by the VP method apparently. In addition, simulation results show that the iterative correction method on regularization matrix and the IVP method can both improve the reconstruction performance of the Tikhonov regularization method.

Reference
1 Lou C and Zhou H C 2005 Combust. Flame 143 97 DOI:10.1016/j.combustflame.2005.05.005 [Cited within:2]
2 Snelling D R, Thomson K A, Smallwood G J, Gülder Ö L, Weckman E J and Fraser R A 2002 AIAA J. 40 1789 DOI:10.2514/2.1855 [Cited within:2]
3 Kimoto A, Nakatani T, Matsuoka Y I and Shida K 2005 IEEE Trans. Instrum. Meas. 54 2407 DOI:10.1109/TIM.2005.858526 [Cited within:1]
4 Yang F Q, Zhang D H, Huang K D, Wang K and Xue Z 2014 Acta Phys. Sin 63 058701(in Chinese) DOI:10.7498/aps.63.058701 [Cited within:1]
5 Wang Q and Wang H X 2011 ISA Trans. 50 256 DOI:10.1016/j.isatra.2010.11.001 [Cited within:1]
6 Craciunescu T, Bonheure G, Kiptily V, Murari A, Tiseanu I and Zoita V 2009 Nucl. Instrum. Method Phys. Res. A 605 374 DOI:10.1016/j.nima.2009.03.2242009 [Cited within:4]
7 Sidky E Y, Kao C M and Pan X C 2006 J. x-ray Sci. Technol. 14 119 [Cited within:2]
8 Wang Q, Wand H X and Yan Y 2011 Flow Meas. Instrum. 22 295 DOI:10.1016/j.flowmeasinst.2011.03.010 [Cited within:3]
9 Huang Q X, Wang F, Liu D, Ma Z Y, Yan J H, Chi Y and Cen K F 2009 Combust. Flame 156 5653 [Cited within:2]
10 Ishino Y and Ohiwa N 2005 JSME Int. J. B 48 34 DOI:10.1299/jsmeb.48.34 [Cited within:1]
11 Balima O, Favennec Y and Rousse D 2013 J. Comput. Phys. 251 461 DOI:10.1016/j.jcp.2013.04.043 [Cited within:1]
12 Zhu N, Jiang Y and Kato S 2005 Energy 30 509 DOI:10.1016/j.energy.2004.09.005 [Cited within:2]
13 Reis M L and Robery N C 1992 Inverse Probl. 8 623 DOI:10.1088/0266-5611/8/4/011 [Cited within:1]
14 Wei W D 1986 Acta Phys. Sin. 35 171(in Chinese) DOI:10.7498/aps.63.058701 [Cited within:1]
15 Yang F Q, Zhang D H, Huang K D, Wang K and Xue Z 2014 Acta Phys. Sin 63 058701(in Chinese) DOI:10.7498/aps.62.1689022013 [Cited within:1]
16 Lin M B, Zhao C X, Yu X D, Yu F S, Yang T Y and Yan K 2014 Acta Phys. Sin. 63 138701(in Chinese) DOI:10.7498/aps.63.1387012014 [Cited within:1]
17 Gu Y F, Yan B, Li L, Wei F, Han Y and Chen J 2014 Acta Phys. Sin 63 058701(in Chinese) DOI:10.7498/aps.63.058701 [Cited within:1]
18 Fang S and Guo H 2014 Chin. Phys. B 23 057401 DOI:10.1088/1674-1056/23/5/057401 [Cited within:1]
19 Zhao J, Xu Y B, Tan C and Dong F 2014 Meas. Sci. Technol. 25 085401 DOI:10.1088/0957-0233/25/8/085401 [Cited within:1]
20 Yang D W, Xing D, X, Zhao X H, Pan C N and Fang J S 2010 Chin. Phys. Lett. 27 054301 DOI:10.1088/0256-307X/27/5/054301 [Cited within:1]
21 Kim B S, Khambampati A K, Kim S and Kim K Y 2011 Meas. Sci. Technol. 22 104009 DOI:10.1088/0957-0233/22/10/104009 [Cited within:1]
22 Zhou H C, Han S D, Sheng F and Zheng C G 2002 JQSRT 72 361 DOI:10.1016/S0022-4073(01)00130-3 [Cited within:2]
23 Liu H W, Lv W and Zhou H C 2013IEEE Int. Conf. Imaging Systems and TechniquesOctober 22–23, 2013Beijing, China 321 DOI:10.1109/IST.2013.6729714 [Cited within:7]
24 Garda B and Galias Z 2012Int. Conf. Signals and Electronic SystemsSeptember 18–21, 2012Wroclaw, Poland 1 DOI:10.1109/ICSES.2012.6382220 [Cited within:1]
25 Chiao P, Fessler J A, Zasadny K R and Wahl R L 1995IEEE Nuclear Science Symposium and Medical Imaging Conference October 21–28, 1995San Francisco, CA, USA 1680 DOI:10.1109/NSSMIC.1995.501909 [Cited within:1]
26 Bro R and Jong S D 1997 J. Chemom. 11 393 DOI:10.1002/(ISSN)1099-128X [Cited within:1]
27 [Cited within:1]
28 Attivissimo F, Cavone G, Lanzolla A M L and Spadavecchia M 2005 IEEE Trans. Instrum. Meas. 591251101109/TIM. 20102040932 [Cited within:1]