Fused Entropy Algorithm in Optical Computed Tomography

Abstract

In most applications of optical computed tomography (OpCT), limited-view problems are often encountered, which can be solved to a certain extent with typical OpCT reconstructive algorithms. The concept of entropy first emerged in information theory has been introduced into OpCT algorithms, such as maximum entropy (ME) algorithms and cross entropy (CE) algorithms, which have demonstrated their superiority over traditional OpCT algorithms, yet have their own limitations. A fused entropy (FE) algorithm, which follows an optimized criterion combining self-adaptively ME with CE, is proposed and investigated by comparisons with ME, CE and some traditional OpCT algorithms. Reconstructed results of several physical models show this FE algorithm has a good convergence and can achieve better precision than other algorithms, which verifies the feasibility of FE as an approach of optimizing computation, not only for OpCT, but also for other image processing applications.

1. Introduction

Optical computed tomography (OpCT) techniques such as interferometry tomography [1–3], light beam deflection tomography [4], emission spectral tomography [5–7], etc., are a branch of computed tomography (CT), which is mainly applied to optical testing of 3-D distributions of physical variables of a number of fluid fields [8–10]. Due to the limitations of the testing environments and devices, most OpCTs encounter limited-view problems [11–13], e.g., incomplete testing views and/or incomplete data at each view, which results in worse reconstruction precision and lower spatial resolution than that seen in medical CT. To solve this problem, many OpCT algorithms have been developed since Gordon et al. first proposed an algebraic reconstruction technique (ART) algorithm in the 1970s [14]. Gilbert introduced a simultaneous iterative reconstruction technique (SIRT) algorithm [15], in which each reconstructed pixel is revised after all the projection values have been computed in each iterative step. Aderson [16] proposed a simultaneous algebraic reconstruction technique (SART) algorithm, which combined the advantages of ART and SIRT. A natural pixel decomposition (NPD) algorithm was first proposed by Buonocore et al. [17], in which the shapes of the grids of the reconstructed plane are not rectangular but determined by the rays’ paths. Garnero et al. [18] employed the NPD to reconstruct a field of refractivity. Different from those above-mentioned row-relaxation iterative methods, a column-relaxation iterative reconstruction method was proposed by David [19]. Dean et al. [11] put forward a singular-value decomposition (SVD) algorithm that is suitable for the solution of both overdetermined and underdetermined equations. To solve the problem of the loss of projection data when reconstructing fields comprising obstacle objects, a discrete iterative reconstruction reprojection (DIRR) [20] was presented, which combines a low-pass filter with a re-projected estimation of the lost data. A Lagrange interpolation reprojection revising (LIRR) algorithm [21] adopts the pre-estimation of a Lagrange interpolation method to improve the accuracy of the re-projected estimation by DIRR, and has been demonstrated a rather improvement over DIRR.

OpCT algorithms are required to have accurate reconstruction results with incomplete data, which means optimization criteria have to be followed. Typical OpCT algorithms, such as ART, SIRT, SART, etc., mostly comply with a single optimization criterion, with which accurate reconstructed results can hardly be achieved. Although some multi-criterion OpCT algorithms [13] have been proposed, their reconstruction results are still unsatisfactory when the distribution of tested fields is relatively complex. Entropy, that first emerged in the information theory, has also been introduced into OpCT algorithms. Maximum entropy algorithms (ME) [22,23] search a most possible solution from the solution set by maximizing the entropy function of the tested target itself, which is verified as a superior approach compared to conventional OpCT algorithms when the number of views was extremely limited and tested fields are approximately symmetrical. However, the performance of ME degrades when the tested fields have poor symmetry. Cross entropy algorithms (CE) [24,25] figure out a most possible target function by minimizing the cross entropy function that measures the possible relationship between the distributions of the target function and its projections. The reconstruction precision of ME for asymmetrical targets can be improved when combined with CE [26].

In this paper, a fused entropy (FE) algorithm is proposed, which self-adaptively combines ME with CE, and hence has high reconstruction precision for both symmetrical and asymmetrical fields. The performance of FE is investigated by comparisons with ME, CE and some traditional OpCT algorithms in the reconstructions of several physical models. Results of numerical simulations show this FE has a good convergence and a better precision than other algorithms.

2. Principle of FE Algorithm

Like in medical CT, the projection data of fluid fields, which can be probed with optoelectronic sensors, are adopted to compute distributions of physical variables in the OpCT. From a mathematical point of view, the OpCT reconstruction problem can be formulated as the inverse Radon transform. As shown in Figure 1(a), the relationship between the 2-D physical function f (x, y) of a reconstructed plane of a tested fluid field and its projection q(t,θ) is given by:

$$q(t,\theta )={\int }_{-\infty }^{\infty }f(x,y)ds$$

(1)

where q(t,θ) is the integral of f (x, y) along the axis s and the Cartesian coordinate system (t, s) rotates at an angle θ with respect to the original system (x, y).

The 2-D continuous function f (x, y) is normally discretized into grids for the OpCT reconstruction; M(N) even grids are discretized in the X(Y) directions, as shown in Figure 1(b). The discrete expression of the function f (x, y) is:

$$f(x,y)=\sum _{j=1}^{MN}{f}_j{b}_j=\sum _{j=1}^{MN}{f}_{j}b(x-x_j,y-y_j)$$

(2)

where f_j is the reconstructed physical variable at the j-th grid that has a central coordinate (x_j, y_j), MN is the total number of the reconstructed variables, and b is a basis function that has the following property:

$$b(x-x_j,y-y_j)=b(\frac{x-m{l}_x}{l_x})b(\frac{y-n{l}_y}{l_y})$$

(3)

where l_x and l_y represent the span of grids in the X and Y directions respectively and m and n are the sequence number of the j-th grid along these two axes. A sinc basis function is commonly used, which can be expressed as:

$$b(k)=\frac{\text{sin}k\pi }{k\pi }$$

(4)

Then we have:

$$q_i(t,\theta )=\sum _{j=1}^{MN}{f}_j{\int }_{i}b(\frac{x-m{l}_x}{l_x})b(\frac{y-n{l}_y}{l_y})ds=\sum _{j=1}^{MN}{w}_{ij}{f}_j$$

(5)

where i is the ray’s sequence number. The number of total testing views is V, and the number of rays per view is defined as RPV. Then the number of total rays is I, the product of V and RPV. We can rewrite Equation (5) in a matrix form:

$$\mathbf{Q}=\mathbf{WF}$$

(6)

where Q is an I-dimensional measurement vector that consists of q_i(t,θ), W that consists of w_ij is a projection matrix of I rows and MN columns, and F is an image vector that comprises MN reconstructed physical variables f_j. Limited-view OpCT reconstruction is usually an underdetermined problem, since I is less than MN in most practical cases, which hence needs an optimization criterion. An optimal solution can be obtained from a number of possible solutions based on the criterion.

Conventional OpCT algorithms usually follow a single optimization criterion. For example, the ART complies with a minimum norm criterion when a suitable initial image vector is selected and the SIRT subjects to a least-squares criterion.

Entropy concepts first based on information theory have also found application in OpCT algorithms. Maximum entropy (ME) algorithms maximize the entropy function of the reconstructed physical variables f_j (the image vector F):

$${\Phi }_1(\mathbf{F})=-\sum _{j=1}^{MN}{f}_j\text{ln}{f}_j=-{\mathbf{F}}^T\text{ln}\mathbf{F}$$

(7)

to get a solution of greatest global smoothness. ME algorithms are superior to conventional OpCT algorithms in limited-view reconstructions for approximately symmetrical targets, whose performance yet degrades when targets have poor symmetry. Cross entropy (CE), unlike ME, is an information theoretic measure that quantifies the difference between two probability density functions, which has also been applied to the area of OpCT. CE function in OpCT is defined as:

$${\Phi }_2(\mathbf{F})=\sum _{i=1}^I[(\sum _{j=1}^{MN}{w}_{ij}{f}_j)\times \text{ln}(\frac{\sum _{j=1}^{MN}{w}_{ij}{f}_j}{q_i})]=\sum _{i=1}^I{\mathbf{W}}_i\mathbf{F}\times \text{ln}(\frac{{\mathbf{W}}_i\mathbf{F}}{q_i})$$

(8)

which measures the difference of information in WF and Q, and thus must be minimized to get an optimal reconstructed image F.

Here, we propose a fused entropy (FE) algorithm that self-adaptively combines ME with CE, which minimize:

$$\Phi (\mathbf{F})=-{\lambda }_1{\Phi }_1(\mathbf{F})+{\lambda }_2{\Phi }_2(\mathbf{F})$$

(9)

where λ₁ and λ₂ represent for the weight factor of ME and CE respectively. To obtain the solution of F, the matrix differential of Equation (9) must follow:

$$\frac{\text{d}\Phi (\mathbf{F})}{\mathbf{F}}={\lambda }_1(\text{ln}\mathbf{F}+\mathbf{1})+{\lambda }_2\sum _{i=1}^I\text{ln}(\frac{{\mathbf{W}}_i\mathbf{F}}{q_i}){\mathbf{W}}_i^T=\mathbf{0}$$

(10)

Then the following FE algorithm is established:

$$\begin{array}{l}{\mathbf{F}}^0=\mathbf{1}\\ {\mathbf{C}}^k=\mathbf{1}-\alpha [{\lambda }_1^k(\text{ln}{\mathbf{F}}^{\text{k}}+\mathbf{1})+{\lambda }_2^k\text{ln}(\frac{{\mathbf{W}}_i\mathbf{F}}{q_i}){\mathbf{W}}_i^T]\\ {\mathbf{F}}^{k+1}={\mathbf{C}}^k\cdot {\mathbf{F}}^k\end{array}$$

(11)

where C is an MN-dimensional vector and α is a relaxation parameter. Weight factors λ₁ and λ₂ are self-adaptively adjusted in FE algorithm based on the following scheme:

$$\begin{array}{l}{\lambda }_1^0={\lambda }_2^0=\frac{1}{2}\\ \{\begin{array}{l}{\lambda }_1^{k+1}\left|{\Phi }_1^{k+1}-{\Phi }_1^k\right|={\lambda }_2^{k+1}\left|{\Phi }_2^{k+1}-{\Phi }_2^k\right|\\ {\lambda }_1^{k+1}+{\lambda }_2^{k+1}=1\end{array}\end{array}$$

(12)

3. Numerical Simulations

3.1. Physical Models

Several physical models including a top-concaved paraboloid image (TCP), a three random peaks image (TR), a superposition image (TCPTR) of the top-concaved paraboloid and the three random peaks, and a six-peak Gaussian image (SG), which are shown in Figure 2, are chosen to investigate the performance of FE algorithm.

TCP represents a complete symmetric distribution of physical variables and can be expressed as:

$$F_1(x,y)=\{\begin{array}{rr}0,\hfill & \hfill x^2+y^2>{0.37}^2\\ 0.9[1-\frac{x^2}{{0.37}^2}-\frac{y^2}{{0.37}^2}],\hfill & \hfill {0.15}^2<x^2+y^2\le {0.37}^2\\ 0.9[1-\frac{x^2}{{0.37}^2}-\frac{y^2}{{0.37}^2}]-0.5[1-\frac{x^2}{{0.15}^2}-\frac{y^2}{{0.15}^2}],\hfill & \hfill x^2+y^2\le {0.15}^2\end{array}$$

(13)

TR denotes an asymmetric distribution of physical variables and can be formulated as:

$$\begin{array}{c}{F}_2(x,y)=\text{exp}\{-2.77[\frac{{(0.9x+0.2y+0.2)}^2}{{0.1}^2}+\frac{{(0.9\text{y}-0.2x)}^2}{{0.25}^2}]\}+\text{exp}\{-2.77[\frac{{(\frac{1}{2}x+\frac{\sqrt{3}}{2}y-0.2)}^2}{{0.11}^2}\\ +\frac{{(\frac{1}{2}y-\frac{\sqrt{3}}{2}x-0.1)}^2}{{0.3}^2}]\}+0.8\text{exp}\{-2.77[\frac{{(\frac{\sqrt{2}}{2}x+\frac{\sqrt{2}}{2}y)}^2}{{0.25}^2}+\frac{{(\frac{\sqrt{2}}{2}y-\frac{\sqrt{2}}{2}x+0.3)}^2}{{0.1}^2}]\}\end{array}$$

(14)

TCPTR is the sum of TCP and TR:

$$F_3(x,y)=F_1(x,y)+F_2(x,y)$$

(15)

SPG is a function of the rotational symmetric distribution:

$$\begin{array}{l}{F}_4(x,y)=\sum _{i=1}^6{a}_i\text{exp}\{[-\frac{4\text{ln}2}{{0.18}^2}[{(x-x_i)}^2+{(y-y_i)}^2]\}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{a}_1=a_4=0.7,\hspace{0.17em}\hspace{0.17em}{a}_2=a_5=1,\hspace{0.17em}\hspace{0.17em}{a}_3=a_6=0.5\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{x}_1=0.3,\hspace{0.17em}\hspace{0.17em}{y}_1=0;\hspace{0.17em}\hspace{0.17em}{x}_2=0.3\text{cos}(\frac{\pi }{3}),\hspace{0.17em}\hspace{0.17em}{y}_2=0.3\text{sin}(\frac{\pi }{3});\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{x}_3=0.3\text{cos}(\frac{2\pi }{3}),\hspace{0.17em}\hspace{0.17em}{y}_3=0.3\text{sin}(\frac{2\pi }{3});\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{x}_4=-0.3,\hspace{0.17em}\hspace{0.17em}{y}_4=0;\hspace{0.17em}\hspace{0.17em}{x}_5=0.3\text{cos}(\frac{4\pi }{3}),\hspace{0.17em}\hspace{0.17em}{y}_5=0.3\text{sin}(\frac{4\pi }{3});\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{x}_6=0.3\text{cos}(\frac{5\pi }{3}),\hspace{0.17em}\hspace{0.17em}{y}_6=0.3\text{sin}(\frac{5\pi }{3})\end{array}$$

(16)

Three errors are adopted to evaluate the construction quality. The first error is the average error:

$$e_1=(\sum _{j=1}^{MN}|f_j-{f_j}^{\prime }|)/(f_{j\text{max}}\times MN)$$

(17)

The second error is the maximum error that is defined as:

$$e_2={\left|f_j-{f_j}^{\prime }\right|}_{\text{max}}/f_{j\text{max}}$$

(18)

The third one is the mean square error:

$$e_3=\sqrt{\sum _{j=1}^{MN}{(f_j-{f_j}^{\prime })}^2/\sum _{j=1}^{MN}{f}_j^2}$$

(19)

Where f is the original image and f’ is the reconstruction image.

3.2. Reconstruction Results

The performance of FE is investigated by comparisons with ME, CE and two traditional OpCT algorithms, i.e., the algebraic reconstruction technique (ART) algorithm and the simultaneous iterative reconstruction technique (SIRT) algorithm, in the reconstructions of the four physical models, where M = N = 256, RPV = 256, and V = 6 (evenly distributed views over the range of 180 degrees). The reconstruction errors of the five algorithms are shown in

Table 1. Reconstruction results of FE compared with other algorithms for four physical models.

Physical Models	Errors (%)	Algorithms
		FE	CE	ME	ART	SIRT
F1	e1	0.50	0.56	0.51	0.59	1.20
	e2	4.63	4.67	4.74	6.02	8.29
	e3	2.23	2.41	2.55	2.79	5.23
F2	e1	1.51	1.56	1.27	2.17	2.99
	e2	11.87	12.17	14.09	19.02	27.52
	e3	11.84	13.66	13.32	19.27	25.75
F3	e1	1.23	1.29	1.42	2.12	2.73
	e2	9.33	9.67	17.79	16.13	23.72
	e3	8.76	9.03	11.41	15.23	19.40
F4	e1	0.70	0.72	0.93	1.33	1.71
	e2	3.29	3.33	4.87	8.09	11.17
	e3	3.70	3.80	4.87	7.33	9.79

, where numbers in bold italic are the best results.

The four reconstructed images with the FE algorithm are shown in Figure 3, where the relaxation parameter α of FE is 0.3 for all the four physical models. Furthermore, the convergence of FE has been studied. Figure 4 shows the convergence properties of FE for the reconstructions of TCP and TR.

4. Conclusions

Numerical simulations show that the FE proposed in this paper is superior to the other four algorithms tested in the reconstructions of four physical models (refer to Table 1; note: only the average error is little greater than that of ME), and images reconstructed with FE are of similar distributions as the original physical models (refer to Figure 2 and Figure 3). Besides, FE also has a good convergence (refer to Figure 4). These studies have testified the feasibility of FE as an approach of optimizing computation, which can not only be utilized for OpCT reconstructions, but also be promoted to find an optimal solution for other imaging processing problems if an FE optimization function is set up [27]. However, owing to the complexity, such as turbulences and impulses in the fields of real OpCT applications, much deeper research needs to be conducted to verify the applicability of FE.