High spatial resolution computed tomography of chemiluminescence with densely sampled parallel projections

Jia Wang; Zhenyan Guo; Liang Nie; Shenjiang Wu

doi:10.1364/OE.27.021050

1. Introduction

Chemiluminescence of combustions has been recognized as essential indicators for combustion diagnostics, monitoring and control [1,2]. Comparing with other costly and cumbersome optical techniques [3–8], the optically active emitting nature, simplicity of implementation, and capability in measuring multiple important species make it attractive for applications in practical industrial environments. The projection of chemiluminescence captured by detectors is the line-of-sight integral result; therefore, by utilizing the computed tomography (CT) theory, computed tomography of chemiluminescence (CTC) can provide instantaneous three-dimensional (3D) information on flame geometry and excited species.

Recent researches about CTC mainly focused on the optimization of algorithms to improve the reconstruction accuracy [9–11] and the establishment of multi-directional projection acquisition system [12–16]. To the knowledge of authors, in all of the existing CTC systems, ordinary commodity lenses were utilized as the imaging element to corporate with CCD or fiber bundles for its convenience. Since the imaging characteristics of general lens is different from the parallel projection in classical CT, several cone beam projection models [17] and camera imaging projection models [13,18,19] have been developed to improve the precision of reconstruction. However, the calculation of weight matrices considering the imaging effect of lens is extremely complicated, which limits the computing efficiency of CTC. Besides, spatial resolution is a critical parameter for CTC system, which is affected by many factors such as the sampling rate of projections, the accuracy of the imaging model, the noise level in the measurements, the performance of reconstruction algorithms [13,20]. The sampling rate of camera equipped with general lens is relatively low as the small amplification ratio is designed to guarantee the large field of view. It results in low spatial resolution of the reconstruction and the detailed structure of the combustion fields cannot be well observed.

In the applications of classical CT, parallel projection based on 2D Radon transform is the simplest and widely used model, and the corresponding reconstruction algorithms are mature. Fortunately, in optical imaging, double telecentric lens or telescope imaging system can solve the parallax influence of ordinary industrial lens with advantages of high resolution, ultra-wide depth of field, and low distortion. It can realize the invariability of the image magnification with the object distance within a particular object distance range, and make the image be equivalent to parallel projection with high sampling rate. Therefore, in this work, double telecentric lens is used to replace the ordinary lens to implement the parallel projection CTC.

Although the theories for parallel projection CT are relatively mature, there are still a few theoretical and technological drawbacks need to be solved when utilizing double telecentric lens in CTC. Firstly, the mathematical foundation of parallel projection tomography is 2D discrete Radon transform and the weight matrix is generally determined by calculating the intersection length of the projection rays to the discrete grids in the reconstruction region [12,14]. Our research found that its accuracy decreases with the increasing of the sampling rate of projections. In order to solve the problem, this paper proposes a method based on bilinear interpolation to improve the precision of weight matrix. On the other hand, a critical issue for CTC is the establishment of the multi-directional projection acquisition system, which involves the determination of the parameters of the cameras, the number of cameras and their spatial arrangement. Although the influence of the number and spatial orientation of the projections on the accuracy and spatial resolution of tomographic reconstruction has been studied [13,20], they have not been theoretically analyzed. In this paper, an exact reconstruction condition will be studied based on the discrete Radon transform theory to describe the dependence of reconstruction accuracy on spatial resolution, the sampling rate and number of projections. Besides, the camera calibration technique for double telecentric lens and the reconstruction of combustion flame will also be studied.

The remainder of this work is organized as follows: Section 2 describes the imaging principle of double telecentric lens, the tomographic model of parallel projection and the mathematical formulation for CTC inversion; Section 3 introduces the new bilinear interpolation method to compute the weight matrix of CTC; the relationship between reconstruction accuracy and the number, sampling rate, spatial distribution of projections and spatial resolution of the reconstructed region is analyzed in Section 4; the camera calibration method is provided in Section 5 and the 3D reconstruction of a combustion flame is conducted to demonstrate the theories; the final section summarizes this work.

2. Imaging model and tomography theory

2.1 Imaging of double telecentric lens

Telecentric lens can eliminate the variation of magnification caused by the difference of distance between the lens to the measured object or detector in the ordinary camera. The double telecentric lens combines the functions of object telecentric and image telecentric lenses, which can make the image magnification constant within a certain region and realize the detection of parallel light. The imaging optical path of double telecentric lens is shown in Fig. 1.

Fig. 1 The imaging diagram of double telecentric lens.

Download Full Size | PDF

In the double telecentric lens, an aperture is placed on the focus plane of the optical system, the primary rays on both object side and image side are parallel to the optical axis, and they pass through the center point of the aperture. For an object $A B$ located on the wording distance (the designed focus plane of the lens) with height $y$ , its image $A' B'$ coincides with the detector plane and with an image height $y'$ . For an object $A_{1} B_{1}$ with the same height beyond the focus plane, the image $A_{1}' B_{1}'$ is out of the detector plane and projects a blurry circle on the detector plane. However, the center image height of the blurry circle is equal to $y'$ . The filtering aperture on the focal plane makes the parallel rays pass through and the accurate non-blurred image can be received on the detector plane, which ensures a sufficiently large depth of field. All objects within the depth of field can form a clear image on the CCD with the same magnification; therefore, the parallel projections can be acquired. The critical parameters of double telecentric lens are the amplification ratio $K = y' / y$ , the field of view and the working distance.

2.2 Parallel projection model and reconstruction

Assuming that the detector is placed right facing to the measured object and perpendicular to the horizontal plane, the imaging of 3D object through double telecentric lens can be simplified to the geometric model of parallel projection, as shown in Fig. 2(a). The tested object locates in the world coordinate system $(x, y, z)$ , and its center point coincides with the origin. The origin of image coordinate system $(ξ, η, γ)$ is on the center pixel of the detector, and the optical axis γ has an angle $ϕ$ to axis $x$ in the horizontal plane.

Fig. 2 The parallel projection model (a) three-dimensional, (b) two-dimensional.

Download Full Size | PDF

The optical axis $γ = (\cos ϕ, \sin ϕ)$ determines the projection direction, the projections recorded by each row pixels of the detector $ξ = (- \sin ϕ, \cos ϕ)$ are 2D Radon transform of the 2D horizontal plane, as shown in Fig. 2(b). For an object with the function of $f (x, y, z)$ , the continuous projection can be mathematically expressed as

I (ξ) = \int_{- \infty}^{\infty} \int_{- \infty}^{\infty} f (x, z) δ (- \sin ϕ x + \cos ϕ z - ξ) d x d z

Detectors with discrete pixels are used to record the projection and the tomographic reconstruction is conducted in discrete domain. Assuming that the pixels record the rays striking to the center of the pixel, therefore the continuous projections are spatially discretely sampled with a period of pixel size. The tested object is divided into P × Q discrete grids with the size $Δ g$ , and assumes that the object function in each grid is constant. The sampling rate of projection is defined as the number of pixels imaged by each discrete grid. If the amplification ratio of the double telecentric lens is $K$ and the size of pixels is $Δ p$ , the image of each grid is with the size of $K Δ g$ , and it is sampled by $K Δ g / Δ p$ pixels, therefore, the sampling rate is $R = K Δ g / Δ p$ .

The discretized projection can be expressed as

I_{q p} = \sum_{i = 1}^{P} \sum_{j = 1}^{Q} w (i, j, p, q) f (i, j)

In Eq. (2),

w (i, j, p, q)

represents the contribution of the grid

(i, j)

to the

p th

projection of the

q th

view. The contributions are defined as weight matrix of the tomographic system, which depends on the projection geometry.

Although many algorithms have been studied for tomographic reconstruction, the algebraic reconstruction technique (ART) algorithm is still the most basic and widely used one [12,18,19]. The mathematical description of the reconstruction process can be expressed as

f^{(h + 1)} = f^{(h)} - λ \frac{w_{p q} \cdot f^{(h)} - I_{p q}}{w_{p q} \cdot w_{p q}} w_{p q}

An initial estimation

f^{(0)}

is made to start the iteration. Within each iteration, the projections are used one by one to update the solution, where

f^{(h)}

is the solution after

h th

iteration. The relaxation factor

λ \in (0, 2)

improves convergence of the process. The reconstruction is converged once the relative change in the 2-norm of the solutions between two consecutive iterations is smaller than the threshold value

Δ c

.

3. The calculation of weight matrix

3.1 The method based on bilinear interpolation

In the existing literatures, the weight matrix was generally considered to be the intersection length of the projection rays and the grids (called intersection-length method in the following content), which can be mathematically expressed as [21]

w = {\begin{matrix} 1, | x | \leq \frac{Δ g (| | \cos θ | - | \sin θ | |)}{2} \\ \frac{- 2 | x | + (| \cos θ | + | \sin θ |) Δ g}{2 Δ g \min (| \cos θ |, | \sin θ |)}, \frac{Δ g (| | \cos θ | - | \sin θ | |)}{2} < | x | \leq \frac{Δ g (| | \cos θ | + | \sin θ | |)}{2} \\ 0, | x | > \frac{Δ g (| | \cos θ | + | \sin θ | |)}{2} \end{matrix}

where

θ

is the projection angle,

x

is the distance between the projection pixel and the center of the grid along the axis

ξ

, and

Δ g

is the size of discrete grid.

However, when the projection rays intersect the grid at the corners, the length is approximately 0, as shown in Fig. 2(b), which can introduce discrete errors in the projections. To solve the problem, a weight matrix calculation method based on bilinear interpolation (called interpolation method in the following) is studied.

As shown in Fig. 3(a), seeing from the projection direction, the parallel projection can be expressed as

I (ξ) = \int_{γ_{\min}}^{γ_{\max}} f (ξ, γ) d ξ

In Eq. (5),

γ_{\min}

and

γ_{\max}

are the lower and upper intersection boundaries of the rays with the reconstructed region. The integral interval

[γ_{\min}, γ_{\max}]

is evenly divided into G infinitesimal regions (i.e. the step number of integral is G), and the step length is

Fig. 3 Schematic diagram of the interpolation method.

Download Full Size | PDF

Δ γ = \frac{γ_{\max} - γ_{\min}}{G}

Assuming that the object function in each infinitesimal interval is constant, therefore

I (ξ) = \sum_{m = 0}^{G} f (ξ, γ_{m}) Δ γ = \sum_{m = 0}^{G} f (ξ, γ_{\min} + m Δ γ) Δ γ

Based on the coordinate relationship of

(x, z)

and

(ξ, γ)

, the position of the infinitesimal in

(x, z)

coordinate system can be obtained as

\begin{array}{l} x_{m} = \cos ϕ (γ_{\min} + m Δ γ) - \sin ϕ ξ \\ z_{m} = \sin ϕ (γ_{\min} + m Δ γ) + \cos ϕ ξ \end{array}

Thus,

I (ξ) = \sum_{m = 0}^{G} f (x_{m}, z_{m}) Δ γ

According to the bilinear interpolation principle, the value $f (x_{m}, z_{m})$ can be determined by the adjacent four discrete grids. The contribution or weight of the adjacent discrete grid to the infinitesimal is inversely proportional to their distances, shown as Fig. 3(b). The index of the grid $(i_{m}, j_{m})$ where the infinitesimal $(x_{m}, z_{m})$ located can be determined by

\begin{array}{l} i_{m} = f l o o r (\frac{P + 1}{2} - \frac{z_{m}}{Δ g}) \\ j_{m} = f l o o r (\frac{x_{m}}{Δ g} + \frac{Q + 1}{2}) \end{array}

And

\begin{array}{l} f (x_{m}, z_{m}) = \frac{1}{Δ g^{2}} [(Δ g - Δ x_{m}) (Δ g - Δ z_{m}) f (i_{m}, j_{m}) + Δ x_{m} (Δ g - Δ z_{m}) f (i_{m}, j_{m} + 1) \\ + Δ z_{m} (Δ g - Δ x_{m}) f (i_{m} + 1, j_{m}) + Δ x_{m} Δ z_{m} f (i_{m} + 1, j_{m} + 1)] \end{array}

Therefore, the projection can be written in the form of

\begin{array}{l} I (ξ) = \sum_{m = 0}^{G} \frac{Δ γ}{Δ g^{2}} [(Δ g - Δ x_{m}) (Δ g - Δ z_{m}) f (i_{m}, j_{m}) + Δ x_{m} (Δ g - Δ z_{m}) f (i_{m}, j_{m} + 1) \\ + Δ z_{m} (Δ g - Δ x_{m}) f (i_{m} + 1, j_{m}) + Δ x_{m} Δ z_{m} f (i_{m} + 1, j_{m} + 1)] \end{array}

The weight factor of each discrete gird to the projection is the sum of the contributions of all infinitesimal in the grid.

3.2 The accuracy of the interpolation method

The accuracy of the proposed method is verified by calculating the projections of a square on different angles. The square is divided into 64 × 64 discrete grids with the size $Δ g = 1$ , and the object function values in the square are all equal to 1.

The projections calculated by the interpolation method and traditional intersection-length method are compared with the theoretical projection. The theoretical projections of the square are equal to the intersection length of the projection rays with the square, which can be computed by Eq. (4) when $x$ is the distance of the projection rays to the square center and $Δ g$ is the whole length of the square region.

Figure 4 shows the projections of the square at 72° under different sampling rate. It can be seen that, when the sampling rate $R = 1$ , the projections calculated by both methods are in good agreement with the theoretical values. However, as the increase of sampling rate, the errors of the intersection-length method increase, while the interpolation method maintains high precision.

Fig. 4 Comparison of the interpolation method and the intersection-length method.

Download Full Size | PDF

In addition, in the interpolation method, the step number impacts the accuracy of weight matrix. In theory, the smaller the step length, the closer it is to the actual projection. However, small step length will reduce the computational efficiency, which requires that the appropriate step number should be selected. The root means square error (RMSE) of the projections calculated with different step numbers are shown in Fig. 5. The results indicate that, when the step number is five times greater than the grid number of the reconstruction region, the errors are small enough, and the variation of errors tends to be gentle. Therefore, the step number in the actual calculation is selected to be 5 times the number of grids.

Fig. 5 The effect of step number to the projection errors.

Download Full Size | PDF

4. Exact reconstruction condition

In this section, the dependence of the reconstruction accuracy on the sampling rate of projections, the number of the view of projections and their spatial arrangement, and the spatial resolution of the reconstructed region will be studied based on the discrete Radon transform theory.

In tomography theory, Mojette transform is a special form of discrete Radon transform [22–24]. A set of co-prime integer pairs ${(p_{i}, q_{i}), p_{i} \in Z, q_{i} \in Z^{+}}$ determines the geometric directions of Mojette projections. The corresponding projection angles are $ϕ_{i} = \tan^{- 1} q_{i} / p_{i}$ with $ϕ_{i} \in (0, π)$ . The projection values of Mojette transform are the sum of all pixels whose centers are intersected with the projection rays. As Radon transform contains all the Mojette transform information when the projection conditions are the same, the theories of Mojette transform has valuable guidance for Radon transform [21,25].

According to Katz criterion for exact reconstruction [22,23], a $P \times Q$ tested image can be accurately recovered if the projection set satisfy the condition

\sum_{i = 1}^{M} | p_{i} | \geq P o r \sum_{i = 1}^{M} | q_{i} | \geq Q

It manifests that the number of projections for exact reconstruction is related to the choice of projection vectors

(p_{i}, q_{i})

and the number of discrete grids in the reconstruction region. When the larger projection vectors

| p_{i} |

or

| q_{i} |

are chosen, or the number of discrete grids

P

or

Q

is smaller, exact reconstruction can be conducted with fewer projections.

The resolution of Mojette projections varies with projection vector $(p_{i}, q_{i})$ . The space interval between each ray $h_{i}$ and the number of projection rays $B$ on the record plane under the vector $(p_{i}, q_{i})$ are

\begin{array}{l} h_{i} = \frac{1}{\sqrt{p_{i}^{2} + q_{i}^{2}}} \\ B (P, Q, p_{i}, q_{i}) = (Q - 1) | p_{i} | + (P - 1) | q_{i} | + 1 \end{array}

Regardless of the limitation of pixels number, for a detector with the constant sampling rate $R$ for a single grid, i.e. the ray spacing $h_{i} = 1 / R$ , the projection vector should satisfy the relation

| p_{i} | + | q_{i} | \geq \sqrt{p_{i}^{2} + q_{i}^{2}} = R

The projection set

(p_{i}, q_{i})

satisfying the condition can be automatically computed by Farey series.

For a tomographic system with $N$ equally spaced projections, the view angle of each projection is $ϕ_{i} = 180 (i - 1) / N$ , where $i = 1, 2, \dots, N$ .The vector closest to the angle in the projection set is selected as the characteristic of each angle; therefore the maximum summation of the $N$ vectors $\max (\sum_{i = 1}^{N} | p_{i} |, \sum_{i = 1}^{N} | q_{i} |)$ can be determined. When the sampling rate is 10, the summation values for equally spaced projections with different numbers are plotted Fig. 6(a). As the vectors vary irregularly with the projection number; the curve of maximum summation of the vectors is noisy. It clearly shows that, in the overall trend, the values increase with the number of projections, but the values of even numbers are smaller than which of the adjacent odd ones. This is because that, when the number of projections is even, the angle $N / 2 + 1$ is 90°, and the corresponding projection vector (0, 1) has the minimum value in projection set. Its existence decreases the sum of the projection vectors.

Fig. 6 Relationship between the number of projections and (a) the sum of projection vectors, (b) reconstruction accuracy.

Download Full Size | PDF

In the following, the Shepp-Logan phantom is used as the tested object in simulation experiments to numerically study the factors affecting reconstruction accuracy. In the forward process, the projections are numerically computed by the interpolation method in Section 3. In the backward process, the computed projections are used as inputs, and the reconstructions are carried out by using ART algorithm, the relaxation factor $λ = 1.0$ and the threshold value $Δ c = 10^{- 4}$ .

A. Effects of number of projections and number of discrete grids

Firstly, the effect of the projection number on the reconstruction accuracy of objects with different size is studied. In the simulation, the size of the discrete grids in the reconstruction area is constant as 1, while the number of discrete grids adjusts the length of the region. Reconstructions for three cases with grid number of 40, 50, and 64 are conducted and shown in Fig. 7, in which the sampling rate is 10 and the projections are evenly distributed around the measured object. Also, Fig. 6(b) shows the reconstruction errors for three cases with different number of projections. From the results, two conclusions can be obtained: (1) for the same object, the reconstruction error decrease as the increase of the number of projections in the overall trend, but the errors of even numbers are larger than which of the adjacent odd ones. Comparing Fig. 6(a) and (b), it can be found that the reconstruction error is inversely proportional to the sum of the projection vectors. Therefore, the sum of the projection vectors is an essential indicator for reconstruction quality; (2) with the same projections, the fewer the number of discrete grids in the reconstruction region, the higher the reconstruction quality.

Fig. 7 Reconstructions for three objects with different number of projections.

Download Full Size | PDF

The results reveal that satisfactory reconstruction quality can be obtained when the numbers of projections and discrete grids meet the following condition

\max (\sum_{i = 1}^{N} | p_{i} |, \sum_{i = 1}^{N} | q_{i} |) \geq P / 2

It is called exact reconstruction condition in the following. The minimum projection numbers satisfying the condition for three cases are 7, 9, and 9 respectively, with which the phantoms are well reconstructed. With the same number of projections, when the number of discrete grids is small, the exact reconstruction condition is easier to satisfy, so the reconstruction quality is better.

B. Effect of projections’ arrangement

The projections’ spatial arrangement significantly impacts the reconstruction quality. In the existing studies, projections can be placed in either uniform or non-uniform manner [26]. This section numerically studies the effect of the spatial positions of the projections on the reconstruction accuracy. In the simulation, the sampling rate is 10, the girds number is 50, and the number of projections is 9.

First, the projections are non-uniformly sampled in a relatively small angle range. The vectors with large $| p_{i} |$ values are selected from the projection set satisfying condition $| p_{i} | + | q_{i} | = 10$ , and (p, q) = [(9, 1), (−9, 1), (−8, 1), (8, 1), (−7, 3), (7, 3), (−7, 2), (7, 1), (−7, 1)]. The maximum absolute sum of vector is 55. The corresponding angles are [6.3, 173.7, 172.9, 7.1, 156.8, 23.2, 164.1, 8.1, and 171.9], which are limited in the range of 23.2°. The reconstruction result is shown in Fig. 8(a) and RMSE is 0.1067.

Fig. 8 Reconstructions with different projections’ arrangements (a) in the angle range of 23.2°, (b) in the angle range of 63.4°, (c) evenly spaced in 180°, and (d) randomly spaced in 180°.

Download Full Size | PDF

Second, extend the projection angles to a large range and take the projection vectors as (p, q) = [(9, 1), (−9, 1), (7, 3), (−7, 3), (5, 4), (−5, 4), (1, 1), (−1, 1), (1, 2)], the maximum absolute sum is 51. The corresponding angles are [6.3, 173.7, 23.2, 156.8, 38.7, 141.3, 45, 135, and 63.4], which are limited in the range of 63.4°. The reconstruction result is shown in Fig. 8(b) and the root means square error is 0.0896.

Third, when the 9 projections are uniformly spaced in 180°, the projection vectors are (p, q) = [(1, 0), (3, 1), (5, 4), (3, 5), (1, 6), (−1, 6), (−3, 5), (−5, 4), (−3,1)], the sum of vectors is 32. The reconstruction result is shown in Fig. 8(c) and the root means square error is 0.0629.

Finally, nine vectors are randomly selected from the projection set, and (p, q) = [(3, 4), (−3, 1), (4, 3), (−4, 5), (−5, 1), (−1, 8), (−8, 1), (4, 1), (1,5)]. The maximum absolute sum of vector is 33 and the angles are [53.1, 161.6, 36.9, 128.7, 168.7, 97.1, 172.9, 14.0, and 78.7]. The reconstruction result is shown in Fig. 8(d), and the root mean square error is 0.0702.

From the results, it can be concluded that, for the same number of projections, the best reconstruction quality can be obtained by uniform arrangement. When the projections are limited in small angular range, the reconstruction quality is reduced even if the absolute sum values of the projection vectors are large. When the projections are sampled in the range of 180°, the uniform distribution is more sensible than the random distribution. The exact reconstruction condition in Eq. (13) is only satisfied when the projections are evenly placed in 180°.

C. Effect of projection sampling rate

According to Eq. (13) and (14), if the sampling rate of projection is changed, the number of projections satisfying exact reconstruction condition also changes. For the phantom with the grid number of 50, when the projection sampling rate is 20, the minimum number that satisfies Eq. (14) is 5 (the absolute sum of projection vectors is 25), and the reconstruction results are as shown in Fig. 9(a). When the sampling rate is reduced to 5, the minimum number satisfying Eq. (13) is 13 (the absolute sum of projection vectors is 26), the results are shown in Fig. 9(b). Comparing it with the second row in Fig. 7, it can be found that when the sampling rate is increased, better reconstruction results can be obtained with fewer angles.

Fig. 9 Reconstructions with different sampling rates.

Download Full Size | PDF

It should be specially noted that the exact reconstruction condition is applicative only when the sampling rate is high enough. When the sampling rate is low, the number of vectors in projection set is few, so the vector to each projection determined by the proposed theory has a large error. Practical simulations have shown that the absolute sum of all vectors in the projection set determined by $| p_{i} | + | q_{i} | = R$ must be higher than the half number of the grids; otherwise the theory cannot be used. For example, when $| p_{i} | + | q_{i} | = 5$ , the absolute sum of vectors is 37, therefore the condition is suitable for the reconstruction of phantom with 50 discrete grids. While when $| p_{i} | + | q_{i} | = 4$ , the vector corresponding to each uniformly distributed projection cannot be accurately determined, the reconstruction results no longer follow the exact reconstruction condition.

D. Effect of the resolution of the reconstructed region

To study the effect of the resolution of the reconstruction region, in the simulations, the region is divided into finer grid while the total length is kept constant, therefore the grid number increase.

When the grid number is 50, the grid size is 1, and the sampling rate of projection is 10 (the pixel size is 0.1), the minimum number of projections required is 9, and the reconstruction result is shown in Fig. 10(a).

Fig. 10 Reconstructions for phantom with different resolutions.

Download Full Size | PDF

The resolution of the reconstructed area is increased by dividing the area into 100 discrete grids with the grid size of 0.5. Firstly, the parameters of the projections are not changed, that is, the pixel size of projections remains at 0.1. In this case, the sampling rate for each grid is 5. The reconstruction result by 9 projections is shown in Fig. 10(b), and it can be seen that serious artifact arises. This is because that the projection vectors $(p_{i}, q_{i})$ take small values and the grid number N increases, the exact reconstruction condition cannot be satisfied by 9 projections. Secondly, keep the sampling rate 10 constant, that is, the pixel size of projections becomes 0.05, and the reconstruction result by 9 projections is shown in Fig. 10(c). It can be seen that the reconstruction quality is improved, but artifacts still exist. Comparing with Fig. 10(a), $(p_{i}, q_{i})$ are unchanged while the value of N is increased, the condition also cannot be satisfied.

In the case when the sampling rate is 10, the minimum number of projections satisfying Eq. (14) is 15 with the sum of vectors is 50. As shown in Fig. 10(d), the phantom is well reconstructed.

It can be concluded that, regardless of the resolution or the grids number of the reconstructed area, as long as the exact reconstruction condition in Eq. (14) is satisfied, the phantom can be well reconstructed.

5. Experiment and flame reconstruction

To expediently illustrate the theories detailed in the previous section, an example experiment is performed with an axisymmetric flame. Only one camera composed by a CCD and a double telecentric lens is used for projection recording in this work. Any number of projections can be easily obtained by making copies of its image. The CCD used is Basler scA780-54gm, with a pixel resolution of 782 × 528 and a pixel size of 8.3μm × 8.3μm. A double telecentric lens GCO-23 is used to collect the parallel projection, with the amplification ratio of 0.16, the field of view of 40mm × 30mm, and the working distance of 146.3mm. The camera is placed facing to the flame and make its principal point located in the symmetrical axis of the flame, and the detector plane is perpendicular to the horizontal plane.

5.1 Camera calibration and projection correction

As inevitable errors exist in the actual installation of cameras, there are certain angular deviations and planar translations to the pre-set positions of projections. Camera calibration is one of the key techniques in CTC to exactly determine the true orientation of each projection. In the existing CTC with ordinary lens, many standard camera calibration methods based on pinhole imaging can be used [18,19,27]. However, the imaging of double telecentric lens follows the affine camera model.

In the affine camera model, if the distortion is ignored, the fixed world coordinates system $(x, y, z)$ and the image coordinates system $(ξ, η, γ)$ of a camera arbitrarily placed in 3D space with Euler angles $(ψ, θ, ϕ)$ (Yaw, Pitch, and Roll) and translations $T = {(T x, T y, T z)}^{T}$ satisfy the relationship:

[\begin{matrix} ξ \\ η \\ γ \\ α \end{matrix}] = α [\begin{matrix} R & T \\ 0 & 1 \end{matrix}] [\begin{matrix} x \\ y \\ z \\ 1 \end{matrix}]

where the rotation matrix

R = [\begin{matrix} \cos ϕ \cos ψ & \cos ψ \sin ϕ & - \sin ψ \\ \cos ϕ \sin ψ \sin θ - \cos θ \sin ϕ & \sin ϕ \sin ψ \sin θ + \cos ϕ \cos θ & \sin θ \cos ψ \\ \cos ϕ \sin ψ \cos θ + \sin ϕ \sin θ & \sin ϕ \sin ψ \cos θ - \sin θ \cos ϕ & \cos ψ \cos θ \end{matrix}]

The rotation matrix and translations uniquely determine the position and orientation of the camera, which are called extrinsic parameters of the camera. The magnification ratio

α

is the intrinsic parameter.

If the world coordinates and imaging coordinates of N (N>4) feature points are known, linear equations can be obtained that

B_{2 N \times 8} y_{8 \times 1} = b_{2 N \times 1}

where

B = [\begin{matrix} x_{1} & y_{1} & z_{1} & 1 \\ 0 & 0 & 0 & 0 \\ ⋮ \\ \begin{matrix} x_{N} \\ 0 \end{matrix} & \begin{matrix} y_{N} \\ 0 \end{matrix} & \begin{matrix} z_{N} \\ 0 \end{matrix} & \begin{matrix} 1 \\ 0 \end{matrix} \end{matrix} \begin{matrix} 0 & 0 & 0 & 0 \\ x_{1} & y_{1} & z_{1} & 1 \\ \begin{matrix} 0 \\ x_{N} \end{matrix} & \begin{matrix} 0 \\ y_{N} \end{matrix} & \begin{matrix} 0 \\ z_{N} \end{matrix} & \begin{matrix} 0 \\ 1 \end{matrix} \end{matrix}]

y = {[\begin{matrix} α r_{1} & α r_{2} & α r_{3} & α T x \end{matrix} α \begin{matrix} r_{4} & α r_{5} & α r_{6} & α T y \end{matrix}]}^{T}

b = {[\begin{matrix} ξ_{1} & η_{1} & \dots & ξ_{N} & η_{N} \end{matrix}]}^{T}

Solving the equations, the vector

y

containing the unknown parameters can be obtained. As the parameters in rotation matrix satisfy the relation

\sqrt{{(r_{1} + r_{5})}^{2} + {(r_{2} - r_{4})}^{2}} + \sqrt{{(r_{1} - r_{5})}^{2} + {(r_{2} + r_{4})}^{2}} = 2

the magnification ratio can be computed by

α = \frac{\sqrt{{(y (1) + y (6))}^{2} + {(y (2) - y (5))}^{2}} + \sqrt{{(y (1) - y (6))}^{2} + {(y (2) + y (5))}^{2}}}{2}

And

r_{1} = y (1) / α

,

r_{2} = y (2) / α

,

r_{3} = y (3) / α

,

r_{4} = y (5) / α

,

r_{5} = y (6) / α

,

r_{6} = y (7) / α

,

T x = y (4) / α

, and

T y = y (8) / α

. Then

r_{7}

,

r_{8}

and

r_{9}

can be determined by the orthonormal and right-handed property. Since the parallel imaging does not be affected by

T z

, it is not calibrated. The Euler angles can be determined by

\begin{array}{l} ψ = arc \sin (- r_{3}), ψ \in [- \frac{π}{2}, \frac{π}{2}] \\ θ = arc \sin (\frac{r_{6}}{\sqrt{1 - r_{3}^{2}}}), θ \in [- \frac{π}{2}, \frac{π}{2}] \\ ϕ = arc \sin (\frac{r_{2}}{\sqrt{1 - r_{3}^{2}}}), ϕ \in [- \frac{π}{2}, \frac{π}{2}] \end{array}

Equation (16) requires that the feature points used for calibration cannot be lying in the same plane; otherwise the linear equations cannot be solved. In the calibration process, a 3D dot-matrix board with 2mm × 2mm spaces shown in Fig. 11 was placed in the center position where the flame was supposed to be. The image coordinates of the dots are extracted from the recorded images by computing the barycentric coordinates of each dot. By using the coordinates of the dots in the both world and image coordinate system as inputs, the extrinsic and intrinsic parameters can be solved by above calibration theory and are listed in Table 1.

Fig. 11 The 3D dot-matrix board for camera calibration.

Download Full Size | PDF

Table 1. Camera calibration results

View Table

The precision of camera calibration is checked by re-projecting the feature points with the same world coordinate to the image coordinate system. By using the calibration result, the maximum absolute error of the re-projected image coordinates and the extracted image coordinates is 0.47 pixels, proving that the camera calibration is accurate.

Replace the calibration board by the tested flame and make sure that their symmetric axes coincide with each other. Record the luminescence projection of the flame by the camera. Since the camera is rotated −1.41° in the detection plane and translated 1.29mm and −4.26mm on both horizontal and vertical directions, the projection of the flame is askew and asymmetric, which will seriously affect the reconstruction quality. Therefore, before the tomographic reconstruction, the projection must be corrected by rotating 1.41° and translating 1.29mm and −4.26mm in the horizontal and vertical direction respectively. The corrected projection is shown in Fig. 12.

Fig. 12 The corrected flame projection.

Download Full Size | PDF

5.2 Flame reconstruction

The effective height of the flame in the projection is 430 pixels size, i.e. 430 × 8.3μm. Therefore, the measured flame is divided into 430 cross-section slices along the vertical direction, and the slices are reconstructed separately. The reconstruction region in each slice is 24mm × 24mm. First, the reconstruction region is discretized into 40 × 40 grids, resulting in a spatial resolution of 0.6mm. In this case, the sampling rate of each grid is 12. According to the exact reconstruction condition, 6 evenly spaced projections are required for reconstruction. The 3D structure of the flame is obtained by performing the tomographic reconstruction using the ART algorithm, and the cross sections and longitudinal section are shown in Fig. 13(a). It can be seen that, although certain radial artifacts are existing in the reconstructed cross section, the 3D structure of the flame is well reconstructed.

Fig. 13 Reconstructions of the flame with different spatial resolution and projection number.

Download Full Size | PDF

If the reconstructed region is divided into 60 × 60 grids with higher spatial resolution of 0.4mm, the sampling rate is about 8 and at least 9 evenly spaced projections are required for reconstruction. The reconstruction result is shown in Fig. 13(d).

In addition, comparative reconstructions are performed. When 9 projections are used for the reconstruction of the resolution 0.6mm, as the number of projections is greater than the minimum number determined by the exact reconstruction condition, satisfactory result is obtained. However, the reconstruction result for spatial resolution 0.4mm with 6 projections is poor for the insufficiency of projection.

The experiment verifies the reliability and accuracy of high spatial resolution CTC with telecentric lens, and further demonstrates the exact reconstruction condition.

6. Conclusion

In summary, this work studies the theories for high spatial resolution CTC with densely sampled parallel projection by using double telecentric lens. It has been shown that the proposed bilinear interpolation for the calculation of weight matrix has higher accuracy over the existing intersection-length method, especially when the sampling rate is high. In addition, the exact reconstruction condition developed based on discrete Radon transform theory describes the dependence of the reconstruction accuracy on the sampling rate, the number of projections, the spatial resolution and the number of discrete grids in the tested region. When the condition is satisfied, the object can be well reconstructed. Moreover, the spatial resolution of tomographic reconstruction can be improved by increasing the sampling rate and number of projections. The camera with double telecentric lens can be calibrated by affine calibration method, and the projections can be corrected to the preset positions. The tomographic reconstruction of an axisymmetric flame proves the feasibility and accuracy of the proposed theories. This work provides valuable theoretical and technical guidance for high spatial resolution tomographic reconstruction of CTC and other tomography techniques.

Funding

National Natural Science Foundation of China (NSFC) (61701385, 61701239); Joint Technical and Field Fund for equipment Pre-Research of China (61406190501); Xi'an Key Laboratory of Intelligent Detection and Perception (201805061ZD12CG45).

References

1. J. Kojima, Y. Ikeda, and T. Nakajima, “Basic aspects of OH(A), CH(A), and C2(d) chemiluminescence in the reaction zone of laminar methane-air premixed flame,” Combust. Flame 140(1-2), 34–45 (2005). [CrossRef]

2. Y. K. Jeong, C. H. Jeon, and Y. J. Chang, “Evaluation of the equivalence ratio of the reacting mixture using intensity ratio of chemiluminescence in laminar partially premixed CH4-air flames,” Exp. Therm. Fluid Sci. 30(7), 663–673 (2006). [CrossRef]

3. Y. Song, J. Wang, Y. Jin, Z. Guo, Y. Ji, A. He, and Z. Li, “Implementation of multidirectional moiré computerized tomography: multidirectional affine calibration,” J. Opt. Soc. Am. A 33(12), 2385–2395 (2016). [CrossRef] [PubMed]

4. W. Cai and C. F. Kaminski, “A tomographic technique for the simultaneous imaging of temperature, chemical species, and pressure in reactive flows using absorption spectroscopy with frequency-agile lasers,” Appl. Phys. Lett. 104(3), 034101 (2014). [CrossRef]

5. C. Atkinson, S. Coudert, J. M. Foucaut, M. Stanislas, and J. Soria, “The accuracy of tomographic particle image velocimetry for measurements of a turbulent boundary layer,” Exp. Fluids 50(4), 1031–1056 (2011). [CrossRef]

6. A. Bauknecht, B. Ewers, C. Wolf, F. Leopold, J. Yin, and M. Raffel, “Three-dimensional reconstruction of helicopter blade-tip vortices using a multi-camera BOS system,” Exp. Fluids 56(1), 1–13 (2015). [CrossRef]

7. S. Sharma, G. Sheoran, and C. Shakher, “Investigation of temperature and temperature profile in axi-symmetric flame of butane torch burner using digital holographic interferometry,” Opt. Lasers Eng. 50(10), 1436–1444 (2012). [CrossRef]

8. M. Kumar and C. Shakher, “Measurement of temperature and temperature distribution in gaseous flames by digital speckle pattern shearing interferometry using holographic optical element,” Opt. Lasers Eng. 73, 33–39 (2015). [CrossRef]

9. T. Yu and W. Cai, “Benchmark evaluation of inversion algorithms for tomographic absorption spectroscopy,” Appl. Opt. 56(8), 2183–2194 (2017). [CrossRef] [PubMed]

10. K. J. Daun, S. J. Grauer, and P. J. Hadwin, “Chemical species tomography of turbulent flows: discrete ill-posed and rank deficient problems and the use of prior information,” J. Quant. Spectrosc. Radiat. Transf. 172, 58–74 (2016). [CrossRef]

11. Y. Jin, Y. Song, X. Qu, Z. Li, Y. Ji, and A. He, “Hybrid algorithm for three-dimensional flame chemiluminescence tomography based on imaging overexposure compensation,” Appl. Opt. 55(22), 5917–5923 (2016). [CrossRef] [PubMed]

12. J. Floyd, P. Geipel, and A. M. Kempf, “Computed tomography of chemiluminescence instantaneous 3D measurements and phantom studies of a turbulent opposed jet flame,” Combust. Flame 158(2), 376–391 (2011). [CrossRef]

13. X. Li and L. Ma, “Volumetric imaging of turbulent reactive flows at kHz based on computed tomography,” Opt. Express 22(4), 4768–4778 (2014). [CrossRef] [PubMed]

14. Y. Ishino and N. Ohiwa, “Three-dimensional computerized tomographic reconstruction of instantaneous distribution of chemiluminescence of a turbulent premixed flame,” JSME Int. J. Ser. B Fluids Therm. Eng. 48(1), 34–40 (2005). [CrossRef]

15. G. Gilabert, G. Lu, and Y. Yan, “Three-dimensional tomographic reconstruction of the luminosity distribution of a combustion flame,” IEEE Trans. Instrum. Meas. 56(4), 1300–1306 (2007). [CrossRef]

16. Y. Jin, Y. Song, X. Qu, Z. Li, Y. Ji, and A. He, “Three-dimensional dynamic measurements of CH* and C2* concentrations in flame using simultaneous chemiluminescence tomography,” Opt. Express 25(5), 4640–4654 (2017). [CrossRef] [PubMed]

17. K. T. Walsh, J. Fielding, and M. B. Long, “Effect of light-collection geometry on reconstruction errors in Abel inversions,” Opt. Lett. 25(7), 457–459 (2000). [CrossRef] [PubMed]

18. J. Wang, Y. Song, Z. H. Li, A. Kempf, and A. Z. He, “Multi-directional 3D flame chemiluminescence tomography based on lens imaging,” Opt. Lett. 40(7), 1231–1234 (2015). [CrossRef] [PubMed]

19. T. Yu, H. Liu, and W. Cai, “On the quantification of spatial resolution for three-dimensional computed tomography of chemiluminescence,” Opt. Express 25(20), 24093–24108 (2017). [CrossRef] [PubMed]

20. W. Cai, X. Li, F. Li, and L. Ma, “Numerical and experimental validation of a three-dimensional combustion diagnostic based on tomographic chemiluminescence,” Opt. Express 21(6), 7050–7064 (2013). [CrossRef] [PubMed]

21. B. Recur, P. Desbarats, and J. Domenger, “Radon and Mojette projections’ equivalence for tomographic reconstruction using linear systems,” International Conferences in Central Europe on Computer Graphics, Visualization and Computer Vision, Plzen, Czech Republic, 191–198 (2008).

22. M. B. Katz and R. Gordon, “Questions of uniqueness and resolution in reconstruction from projections,” Phys. Today 32(12), 52–56 (2008). [CrossRef]

23. J. Guédon, “The Mojette transform: theory and applications,” Wiley-ISTE (2009).

24. M. Servieres, N. Normand, J. P. Guédon, and Y. Bizais, “The Mojette transform: discrete angles for tomography,” Electron. Notes Discrete Math. 20, 587–606 (2005). [CrossRef]

25. B. Recur, H. Sarkissian, and M. Servieres, “Validation of Mojette reconstruction from Radon acquisitions,” 20th IEEE International Conference on Image Processing, Melbourne, 1041–1045 (2013). [CrossRef]

26. H. Liu, T. Yu, M. Zhang, and W. Cai, “Demonstration of 3D computed tomography of chemiluminescence with a restricted field of view,” Appl. Opt. 56(25), 7107–7115 (2017). [CrossRef] [PubMed]

27. J. Wang, W. Zhang, Y. Zhang, and X. Yu, “Camera calibration for multidirectional flame chemiluminescence tomography,” Opt. Eng. 56(4), 041307 (2016). [CrossRef]

High spatial resolution computed tomography of chemiluminescence with densely sampled parallel projections

Abstract

1. Introduction

2. Imaging model and tomography theory

2.1 Imaging of double telecentric lens

2.2 Parallel projection model and reconstruction

3. The calculation of weight matrix

3.1 The method based on bilinear interpolation

3.2 The accuracy of the interpolation method

4. Exact reconstruction condition

A. Effects of number of projections and number of discrete grids

B. Effect of projections’ arrangement

C. Effect of projection sampling rate

D. Effect of the resolution of the reconstructed region

5. Experiment and flame reconstruction

5.1 Camera calibration and projection correction

5.2 Flame reconstruction

6. Conclusion

Funding

References

Cited By

Figures (13)

Tables (1)

Equations (25)

Optics Express