Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group

Weighted reconstruction of three-dimensional refractive index in interferometric tomography

Open Access Open Access

Abstract

Interferometric tomography can reconstruct 3D refractive-index distributions through phase-shift measurements for different beam angles. To reconstruct a complex refractive-index distribution, many projections along different directions are required. For the purpose of increasing the number of the projections, we earlier proposed a beam-angle-controllable interferometer with mechanical stages; however, the quality of reconstructed distribution by conventional algorithms was poor because the background fringes cannot be precisely controlled. To improve the quality, we propose a weighted reconstruction algorithm that can consider projection errors. We demonstrate the validity of the weighted reconstruction through simulations and a reconstruction from experimental data for three candle flames.

© 2017 Optical Society of America

1. INTRODUCTION

Computed tomography (CT) is a widely used imaging modality in medical imaging and diagnostics. In CT, 2D projections of the absorption of tissues of the body from different incident beam angles are measured; subsequently, a 3D distribution of the attenuation coefficient is reconstructed by computations. Upon replacing the projection from the absorption with phase-shift information, we can obtain a 3D distribution of the refractive index using the same reconstruction algorithm used in medical CT such as the filtered backprojection method [1], iterative methods [210], and methods reviewed in Refs. [11,12]. According to classical electromagnetism, the refractive index of an object or medium depends on its molecular density and the number of electrons in a molecule. Therefore, the distribution of electron density in plasma, gas temperature, or concentration of a solute can be indirectly obtained through the reconstructed refractive index. The use of interferometry to measure gas temperature has been extensively studied, as reviewed in Refs. [13,14]. In this paper, we will subsequently focus solely on gas temperature measurements based on interferometry techniques. In most of these studies, the object of interest has a simple distribution such as an axisymmetric distribution, in which projections from multiple angles are not required because the Abel transform is applicable in such cases. In regards to a single asymmetric object, several projections have been measured by rotating the object [1517] or by using several fixed views for a non-rotating object [18]. In this context, Yan and Cha employed four projections with fixed angles for three heat sources consisting of electrical resistance wires [19].

The quality of the reconstructed refractive-index profile depends on the number of projections. There are two approaches to increase the number of projections. One approach involves the rotation of the object of interest; however, this approach is restricted to stable objects that are not affected by air. The second approach involves the rotation of the interferometric system itself. Because rotating the entire interferometer is difficult, we have previously proposed a system that consists of two interferometers: one is an incident-angle-controllable interferometer with a limited angular range; the other is a conventional interferometer that can acquire an additional projection, the angle of which is selected from the angular range not covered by the first interferometer [20]. By using the proposed system, it has been demonstrated that the quality of the reconstructed profile is improved by increasing the weight for the additional projection. In that study, the weight was determined from the importance of information from the additional projection. In this paper, we also employ the weight. However, the concept of the weight is extended such that the weight represents the reliability or quality of the projections.

The controllable interferometer mentioned above has two major drawbacks. Because the mirrors are aligned by moving or rotating mechanical stages, the mirrors vibrate during measurements; therefore, the resulting interferogram is blurred. To avoid this blurring, the exposure time should be short; however, the measured interferogram in this case becomes noisy. The second drawback is the uncontrollability of the background fringes in the interferogram because the mirrors on the moving stages cannot be controlled with wavelength precision. These problems affect the accuracy of the 2D phase-shift distribution. In other words, the error in the phase shift is different between projections with different angles.

When the projection error is known, a weighted least-squares method is applicable. Fessler [8] and Anderson et al. [9] applied the weighted least-squares reconstruction to emission tomography in medical diagnostics. Particularly for emission tomography, Shepp and Vardi proposed an iterative algorithm [4] termed the “maximum likelihood expectation maximization” (ML-EM). In emission tomography, radiation sources distribute in the domain to be reconstructed, and the radiation is detected as counts by several detectors. Because the radiation probability obeys the Poisson distribution, the error can be estimated from the observed data as the square root of the count for each detector. However, in the case of projection data obtained by interferometry, the error cannot be directly evaluated as in emission tomography. Against this backdrop, to evaluate the error, we employ singular points around which the phase-shift distribution does not satisfy the sampling theorem.

The purpose of this study is to demonstrate an improvement in the quality of refractive-index reconstruction by the proposed weighted reconstruction method through the example of temperature profile measurements around candle flames. The remainder of this paper is organized as follows. The weighted reconstruction algorithm is detailed in Section 2 with numerical simulations. In Section 3, experimental data acquisition, data processing to obtain the phase-shift distributions, and the reconstructed temperature distribution are presented. Section 4 presents our discussions and conclusions.

2. WEIGHTED RECONSTRUCTION

A. Original ML-EM Algorithm

The proposed weighted reconstruction algorithm shown below is based on the ML-EM method. However, the proposed approach is applicable to other iterative reconstruction algorithms because the weighting process is simply implemented by multiplications of the weights with quantities related to the projections, as described in Section 2.B. ML-EM was originally proposed for emission tomography [4]; however, there is the same implementation for transmission tomography called a “multiplicative type of simultaneous iterative reconstruction” (multiplicative SIRT) [3]. In this case, the internal distribution, f, must always be f0 or always f0.

The ML-EM method searches an internal distribution of the media constant, while applying both a backward projection and a forward projection iteratively such that the ratio between the measured projection and the forward-projected one approaches unity. When we consider the coordinate system shown in Fig. 1, the following equations hold:

pkl(n)=Fkl{fij(n)},
fij(n+1)=Bij{pklMpkl(n)}fij(n),
where f expresses the internal distribution (called density hereafter), which corresponds to the wavenumber being proportional to the refractive index in the case of optical tomography, and p expresses the projection data set called the “sinogram,” which corresponds to the phase shift; the superscripts (n) and M denote the iteration number and measured signals, respectively; the subscripts i, j, k, and l denote the indices of the discrete coordinates of (xi,yj) and (ξk,θl), and Fkl{·} and Bij{·} represent the forward-projection and backward-projection operators, respectively.

 figure: Fig. 1.

Fig. 1. Coordinate system and contribution from internal pixel to projection at a slice corresponding to z=const.

Download Full Size | PDF

The projection data, pkl, which is detected at the kth sensor with width Δs for the lth incident beam angle, can be expressed as an average over the ξ axis of the integral along the η axis of f:

pkl=1ΔsξkΔs/2ξk+Δs/2ηlf(x(ξ,η),y(ξ,η))dηdξ.
Let us assume that the grid size in (x,y) space is also equal to Δs. By denoting the contribution from the pixel at (xi,yj) to the sensor at (ξk,θl) as cijkl, the forward-projection operator in the discrete form can be written as
Fkl{fij}=i,jcijklfij.
The contribution cijkl represents the overlapped area of the (i,j)th pixel and the (k,l)th sensor. The following equations depict the relations between cijkl and some geometrical constants:
A=kCkl,Ckl=i,jcijkl,Lkl=CklΔs,Δs2=kcijkl,
where A represents the entire area of the domain to be reconstructed, Ckl the area of the beam with width Δs incident on the (k,l)th sensor, and Lkl the path length.

The backward projection is evaluated as an angular average of pkl/Lkl, where pkl/Lkl is an averaged density along the η axis:

Bij{pkl}=1NlΔs2k,lcijklpklLkl,
where Nl denotes the number of the beam angles. For a uniform averaged density, pkl/Lkl=f¯, the summation on the right-hand side becomes NlΔs2f¯ from Eq. (5). The factor 1/(NlΔs2) is determined such that the right-hand side becomes f¯.

In Eq. (2), the backward projection is applied to the ratio of the measured projection, pklM, to the current estimated projection, pkl(n); subsequently, the new density, fij(n+1), is evaluated as the product of the result of the backward projection and the current density, fij(n). If the current projection is overestimated, i.e., pkl(n)>pklM, the ratio becomes smaller than unity, and the new density, fij(n+1), is smaller than the current density, fij(n).

Because Bij is applied to the ratio, the iterative procedure becomes unstable if pkl(n)pklM. In our implementation, the density is evaluated by the following modified backward-projection operator, Bij{·}, which can limit the upper bound of the ratio to avoid unstable iterations:

Bij{p^kl}=Bij{min(p^kl,p^lim)},
fij(n+1)=Bij{pklMpkl(n)}fij(n),
where p^lim represents the upper limit of the ratio, which is taken as 2 in our implementation. In regards to the initial estimation of the density, fij(0), in Eq. (1), we use fij(0)=Bij{pklM}, which represents the backward projection of the measured projection data and not the ratio.

B. Weighted Reconstruction Algorithm

With our experimental apparatus (Section 3.A), projection data from certain incident angles cannot be measured due to experimental restrictions; furthermore, the projection data set includes errors due to noise in the measurement process or due to errors in the processing required to obtain the projection data from the measured data. In these cases, the original ML-EM has certain drawbacks. Even when the error is not included, the convergence speed of iterations is slow. Moreover, when noise is included, the iteration does not converge to the true solution. To overcome these drawbacks, we modify the backward-projection procedure to

fij(n+1)=M{Bij{WlpklMpkl(n)}fij(n)},
where the M{·} denotes a moving average operator, and Wl denotes a weight for θl that satisfies
lWl=1.
The purpose of utilizing the moving average is to avoid convergence to incorrect solutions. Meanwhile, the purpose of the weight, Wl, is both rapid convergence and a reduction of the contribution to the reconstructed density from uncertain projection data arising from errors or noise in the measurement process. The idea of a weight to consider the uncertainty of the projection is simple. Figure 2 shows a simulation for the case wherein one of eight projection data sets in a sinogram is incorrect. When we use the incorrect projection data, the reconstructed result is blurred. In contrast, when we do not employ the incorrect data, the result is close to the true density. The latter case is equivalent to the situation that the weight of the incorrect projection data is zero. More detailed characteristics of the modified method are discussed below.

 figure: Fig. 2.

Fig. 2. Improvement of the reconstructed density via exclusion of an incorrect projection. (a) Projections for a square-shaped object along eight directions, where arrows denote the incident directions and the acquired projection data are indicated by the rows of pixels at the end of the arrows. The object for θ4 is replaced by a rectangle instead of the square, i.e., pk,4M is incorrect. (b) Sinogram, pk,lM, and reconstructed density, fi,j(n), at n=100 when pk,4M denoted by the arrow is included. (c) Same as (b) but with pk,4M excluded.

Download Full Size | PDF

1. Simulation for Limited View

As previously mentioned, in the case where noise or error is not included and the density function is a smooth function, the ML-EM converges around the true solution even when some projections are absent. However, the convergence speed varies. Figure 3 compares simulations of reconstruction for different limited sinograms. In the case shown in Fig. 3(b), the reconstructed density is spread because the data of the projections are not complete; i.e., there are sparse ranges of the incident angles. This spread is suppressed by merely adding a projection whose angle is selected from the sparse range, as shown in Fig. 3(C). In the case where the weight for the sparse incident direction is increased, rapid convergence is observed, as shown in Fig. 3(d).

 figure: Fig. 3.

Fig. 3. Simulation for limited projections: (a) True density with three Gaussian functions. (b)–(d) Conditions of the incident beams (top radial chart), and reconstructed densities at n=10 (middle), and those at n=50 (bottom). The projection beam angles, θl, in (b) correspond to θl[30°,30°] with 2° intervals. In both (c) and (d), a horizontal projection (90°) is added to the projections in (b). The weights, Wl, for (b) and (c) are uniform, while Wl(90°)=5Wl(90°) in (d).

Download Full Size | PDF

2. Simulation for Sinogram with Speckle Noise

There are two types of errors in the sinogram pkl, obtained with experimental data processing. The first is a random noise corresponding to each pixel, and the second is a line error that appears as a discontinuity for different l.

The random noise is caused by speckle noise or electrical noise. The sinogram with noise is represented as

pklM=pklTrue+Nkl,pklTrue=Fkl{fijTrue},
where Nkl denotes the random noise. Figure 4 demonstrates a simulation for speckle noise, in which the performances of the moving average, M, are compared. We employed a simple moving average in which the density is convoluted with a uniform kernel with a 3×3 square window. The definition of the relative error plotted in Fig. 4(c) is
Δ^{gpq(n),ref}=||gpq(n)gpqref||2||gpqref||2,||gpq||2=gpq2pq,
where ·pq denotes the average with respect to p and q, and pkl appearing in Fig. 4(c) is pkl=Wlpkl. From Fig. 4(c), it can be observed that the error profiles of the projection, Δ^{pkl(n),M}, decrease monotonically with increasing n, regardless of whether or not the moving average is applied, and they converge to a saturated level when n10. However, the error profiles of the density, Δ^{fij(n),True}, are different. When the moving average is applied, Δ^{fij(n),True} decreases in a manner similar to Δ^{pkl(n),M}, and the results corresponding to fij(10) and fij(50) shown in the right-hand panels of Fig. 4(b)are similar. In contrast, when the moving average is not applied, Δ^{fij(n),True} reaches a minimum value at n=6, and it increases with fij(n) becoming noisy with increasing n. These characteristics can be understood from the principle of ML-EM reconstruction. The original ML-EM minimizes pklMFkl{fij(n)} with increasing n. Let us consider the case where pklM can be expressed as the sum of the true projection, pklTrue, and the error caused by noise, pklδ:
pklM=pklTrue+pklδ.

Although the ML-EM procedure is not an exact linear system, we assume here that the ML-EM is a linear system and the norm of the difference at n= is zero for the true solution. Under this assumption, the solution can be represented as

fij()=Fij1{pklTrue}+Fij1{pklδ}=fijTrue+Fij1{pklδ},
where Fij1{·} denotes the inverse operator of Fkl{·}. This equation indicates that the reconstructed density includes the reconstruction of the noise, Fij1{pklδ}, as well as the true density. When the noise in the projection data is random, Fij1{pklδ} may exhibit a random distribution with higher-frequency components. Because the moving average can reduce the higher-frequency components [21], the contribution of Fij1{pklδ} is reduced. It should be noted that the higher-frequency component in fijTrue is also reduced by the moving average; therefore, the spatial resolution of the reconstructed density is degraded.

 figure: Fig. 4.

Fig. 4. Comparison of moving average for projection data with speckle noise. The true density and the condition of the projection beams are the same as those corresponding to Figs. 3(a) and 3(d), respectively. The speckle noise included in pklM obeys the normal distribution, with the standard deviation being 10% of the maximum of the true projection data.

Download Full Size | PDF

3. Simulation for Sinogram with Line Error

The line error is introduced in the procedures for each θl value because of the following reason. The sinogram, pkl, is a 2D data “slice” on the (ξ,θ) plane represented as p(ξ,θ), which is obtained by selecting a certain slice on z=const from the 3D projection data set in the (ξ,z,θ) coordinate system. The 3D data are obtained by acquiring 2D projection data on the (ξ,z) plane as p(ξ,z) for different θ. Because p(ξ,z) is processed from a measured image, I(ξ,z), for each θl, the error of the sinogram, pkl, does not have any correlation between different l. This means that the error of p(ξ,θ) becomes a discontinuous function even when the error of p(ξ,z) is a continuous function. The sinogram with line error for a simulation is generated by using random numbers, αl, as

pklM=(1+αl)pklTrue.
Figure 5 demonstrates the validity of weighted projection for the line error. The weight for each l in the weighted reconstruction is expressed as
Wlexp((pklMpklTrue)2k/σW2),
where σW denotes a given control parameter. Similar to the case of speckle noise, the convergence of the relative error shown in Fig. 5(c) indicates that Δ^{fij(n),True} for the weighted reconstruction decreases with increasing n, and Δ^{fij(n),True} for the unweighted reconstruction has a local minimum. In the case of the unweighted reconstruction, the maxima of fij(5) and fij(50) are 30% and 70% larger than fijTrue, respectively. Among the reconstructed densities in Fig. 5(b), fij(50) for the weighted reconstruction is the best solution, whose maximum is 10% larger than fijTrue.

 figure: Fig. 5.

Fig. 5. Comparison of unweighted and weighted reconstructions for projection data with line error. The true density and the conditions of the projection beams are identical to the corresponding ones of Figs. 3(a) and 3(d), respectively. The random number, αl, to determine pklM in Eq. (15), obeys the normal distribution with a standard deviation of 1. Parameter σW in Eq. (16) is determined as max{pklTrue}/σW=0.1.

Download Full Size | PDF

To evaluate the effect of the moving average operator for this line error, we simulated the case where the moving average is applied in addition to the weighted reconstruction. Although the detailed results are not shown because of a lack of space, the convergence profile was similar to that shown in Fig. 5. By adding the moving average operator, the relative density error was reduced to Δ^{fij(50),True}=0.219, whereas Δ^{fij(50),True}=0.224 by only applying the weighted reconstruction. This means that, in the case of the sinogram including the line error, the improvement by the moving average is smaller than that by the weighted reconstruction.

From these results, we note that weighted reconstruction can improve the reconstructed density. However, there is a problem in the determination of the weights. In these simulations, because fijTrue is given as a known function to determine pklM, the true projection, pklTrue, required to determine Wl, is also given. In actual applications, however, pklTrue is unknown. Therefore, we must consider another measure to determine the weight. One solution is presented in Section 3.C.

3. EXPERIMENT AND RESULT

A. Experimental Setup

Figure 6 shows our experimental setup, which consists of two Mach–Zehnder interferometers. In the first interferometer, incident angles to an object can be controlled within θ=[45°,15°] by four mirrors (FM1, FM2, HM1, and HM2) on rotatable mechanical stages and by two linearly movable mechanical stages, which are controlled by a computer so that the interferogram can be acquired at each angle. In contrast, the angle in the second interferometer cannot be controlled; the angle is fixed to θ=90°. The object of interest is three flames of different sizes on a candle, as shown in Fig. 8(a). To reduce swaying of flames, small-sized flames were prepared; the diameters of the light-emitting portions were 2mm, 3 mm, and 4 mm.

 figure: Fig. 6.

Fig. 6. Experimental arrangement. FM1-5, full-mirrors; HM1-5: half-mirrors; X1 and X2, linear stages; LS, He–Ne laser (wavelength: 633 nm); L1 and L2, lenses (expanded beam diameter: 50 mm); OBJ, object; S, half-transparent screen. FM1, FM2, HM1, and HM2 are mounted on rotatable stages, which in turn are mounted on X1 and X2.

Download Full Size | PDF

 figure: Fig. 7.

Fig. 7. Process of phase evaluation from an interferogram. (a) Observed interferogram (θ=20.0°). (b) Modified interferogram obtained by both normalization and the Wiener filter. (c)–(f) Spectral intensities in the spectral domain limit k[kmax/2,+kmax/2], where kmax denotes the Nyquist frequency; I^twin(k)=I^1ϕ(kk)+[I^1ϕ(kκ)]* in (e). (g) Wrapped phase with range [π,+π]. (h) Grouped tree of singular points to determine compensated local area for phase unwrapping. (i) Unwrapped phase with contour lines with interval 2π. (j) Wrapped difference between rewrapped phase of (i) and wrapped phase in (g). (k) Additional phase shift to satisfy equi-phase condition on both sides of the interferogram. (l) Corrected unwrapped phase. Image sizes of (a) and (b) are 128×128 pixels, and those of (g)–(l) are 512×512 pixels.

Download Full Size | PDF

 figure: Fig. 8.

Fig. 8. Reconstructed density from experimental data with different operators for a slice (z=10mm). (a) Object with three different-sized flames on a candle. (b) Frequency histogram of SDl[εkl] and weight function, Wl, for weighted reconstruction. (c) Sinogram, pkl; the horizontal and vertical axes represent k and l, respectively. (d) Weight, Wl; the horizontal axis represents Wl within [0,1] and the vertical axis is the same as that in (c). (e) Reconstructed densities obtained with application of different operators in the backward-projection procedure; WGT denotes weighted projection and M.A. moving average. Intervals between contour lines are 2×105 for Δn.

Download Full Size | PDF

The incident angles were controlled for one round-trip in 30 s, and 154 interferograms were acquired. During the measurements, because the mirrors were moving in the angle-controllable interferometer, the mirrors were subject to vibration, and the interferogram was blurred for long exposure times. When we set short exposure times to reduce the blurring, the brightness of the interferogram reduced. To obtain brighter interferograms, we employed a half-transparent screen made of polyethylene film, and a camera with increased gain was placed behind the screen. With increase in the gain, the signal-to-noise ratio decreased. For reducing the effects of the laser profile in the interferogram analysis, we also measured the laser profile. The background fringes, which can be considered as interferograms without the object, were not measured because the mechanical stages cannot be controlled with wavelength-order precision. Instead of background measurements, the carrier frequency, which is the spatial frequency of the background fringe, was obtained from the interferograms, including the object information via choosing regions wherein the phase was not modulated by the object.

B. Interferogram Processing

1. Outline of Processing

For each interferogram, several steps were applied to obtain a 2D phase shift, which is the line integral of the 3D refractive-index distribution of the object.

The intensity of the interferogram, I(r) where r=(ξ,z), on the screen is represented as

I(r)=I0(r)+2I1(r)cos(Δϕ(r)),
Δϕ(r)=ϕobj(r)ϕref(r),
ϕobj(r)=ηk0objnobj(r)·dr,
ϕref(r)=ηk0refnref·dr,
where I0(r) and I1(r) represent intensities related to incident laser profile; n, k0, and ϕ denote the refractive index, wave vector, and phase, respectively; and the superscripts “obj” and “ref” correspond to the object wave and the reference wave, respectively. The wave vectors k0obj and k0ref are directed in different directions, but |k0obj|=|k0ref|=k0, where k0 denotes the wavenumber in vacuum. From these equations, Δϕ(r) can be rewritten as
Δϕ(r)=ϕ(r)+κ·r,
ϕ(r)=ηk0objΔn(r)·dr=k0ηΔn(r)dη,
Δn(r)=nobj(r)nref,
κ=(k0objk0ref)nref.
In the case without the object, because Δn=0 and Δϕ(r)=κ·r, a parallel fringe pattern with carrier frequency, κ, appears in the interferogram. In the case with the object, the carrier pattern is modulated by ϕ(r), which is the 2D phase shift distribution by the object.

The phase shift, ϕ(r), can be estimated as follows. The Fourier transform of I(r) in Eq. (17) distributes to three peaks:

I^(k)=I^0(k)+I^1ϕ(kκ)+[I^1ϕ(kκ)]*,
I^1ϕ(k)=F{I1(r)exp(jϕ(r))},
where F{·} denotes a 2D Fourier transform operator with the integral kernel exp(jk·r), where j denotes the imaginary unit and k the spectral domain; “^” refers to the variables in k space and “*” denotes the complex conjugate. The first term on the right-hand side of Eq. (25) is a dc peak located around the origin in the spectral domain. The second term, corresponding to the carrier peak, and the last term, corresponding to the adjoint peak, distribute around ±κ, respectively. If the peaks are split, the first and last terms can be made zero by using an appropriate spectral filter. After these eliminations, shifting the origin of the spectral domain to κ as I^1ϕ(kκ)I^1ϕ(k), and applying the inverse Fourier transform, we obtain I1(r)exp(jϕ(r)) as a complex-valued image. Because I1(r) denotes a real-valued image, we can obtain a principal value of phase of ϕ(r) by taking the complex argument
ϕw(r)=W{ϕ(r)}=arg(I1(r)exp(jϕ(r))),
where W{·} denotes a wrapping operator to represent the principal value of the angle, which returns a value within [π,+π]. Because the wrapped phase has an ambiguity of multiples of 2π, to solve the ambiguity, a phase unwrapping operator, U{·}, is applied:
ϕ(r)=U{ϕw(r)}.

2. Noise Reduction

An example of the phase evaluation process is shown in Fig. 7. The observed interferogram, Iobs(r), is shown in (a). We note that the observed interferogram, Iobs(r), has a spatial profile with a low frequency, and the brightness at either end is low. The spectrum intensity, |I^obs(k)|, shown in (c) has three bright narrow peaks. The widths of the peaks are mainly determined by the profile I0(r) and I1(r) in Eq. (17). The information of the phase-shift ϕ is spread in a tail around the carrier peak which is one of the bright peaks; however, it is overshadowed by the bright peaks of the profile.

To reduce the effect of the profile, Iobs(r) was normalized by the measured laser source profile, and a Wiener filter was subsequently applied. The result, I(r), and its spectral intensity, |I^(k)|, are shown in Figs. 7(b) and 7(d), respectively.

3. Extraction of Carrier Peak

The carrier frequency, κ, was obtained by using a Fourier transform of the cropped image, not including the object, and a maximum search algorithm with an interpolation using a sinc function in the k domain.

The tail of the dc peak partially overlaps with the tails of the carrier and adjoint peaks. A mask-based filtering of the dc peak corrupts the tail of the carrier peaks due to discontinuous truncation. To reduce the corruption due to the filtering, a filter based on a spectral shifting technique [22] was applied to remove the dc peak. To eliminate the adjoint peak, a half-plane filter in the spectral domain [23] was applied. The dc-removed spectrum, which equals the sum of the carrier and the adjoint peaks, I^twin(r), is shown in Fig. 7(e), and the shifted carrier peak, I^1ϕ(k), is shown in Fig. 7(f).

Because the experimental system could not be used to control the carrier frequency, κ, certain interferograms included a dense fringe region, and unnatural patterns were observed in the wrapped phase. To solve this dense-fringe problem, an interpolation in the real domain was applied by padding zero spectral data in the spectrum domain, by which the sampling interval in the real domain was reduced to one-fourth the original value along each axis. The wrapped phase, ϕw(r), is shown in Fig. 7(g).

4. Phase Unwrapping

An ill-conditioned wrapped phase includes many singular points (SPs) at which the sampling theorem is not satisfied, i.e., the phase difference between adjacent pixels exceeds π. For the noisy wrapped phase, the phase unwrapping method using a localized compensator [24], called LC hereafter, yields a smaller error than other methods such as path-following methods or least-square methods [2431], which are compared in Ref. [32] because the LC method can limit spread of the singularities to local areas. Figure 7(h) shows the SP trees in which the SPs in each local area are connected. We note that the regions including the SPs are the same as those with dense fringes in 7(b). The unwrapped phase, ϕ(r), obtained with the use of the LC method is shown in 7(i).

When the object is located around the central region of the interferogram, the phase around both ends of the interferogram must be identical; this condition is called the “equi-phase condition” hereafter. However, ϕ(r) in Fig. 7(i) does not satisfy this equi-phase condition. This is due to the error in the determination of the local area because the SPs are aligned across the whole region. To fix the “breaking” of the equi-phase condition, we added an additional phase, Δϕ(r), as shown in Fig. 7(k) , which is expressed as the sum of piece-wise linear functions. The factors in the functions were manually determined from the images shown in Figs. 7(g)7(i). The corrected unwrapped phase, ϕ(r)=ϕ(r)+Δϕ(r), satisfying the equi-phase condition is shown in Fig. 7(l).

C. Evaluation of Error in Unwrapped Phase

The SPs in the ill-conditioned, wrapped phase introduce uncertainty of the unwrapped phase, and the additional phases are required to satisfy the equi-phase condition. In the LC method, the rewrapped phase of the unwrapped phase, W{ϕ(r)}, is not identical to the original wrapped phase, ϕw(r), as shown in Fig. 7(j), in which ε(r) is defined as

ε(r)=W{ϕ(r)ϕw(r)}.

We employed the standard deviation of ε(r) as a measure of the quality of the unwrapped phase. Although Δϕ(r) is also a candidate of the measure to evaluate the quality, the dependency between Δϕ(r) and the quality of ϕ(r) was small if the factors to represent Δϕ(r) were determined correctly.

D. Reconstructed Flames

Among the 154 observed interferograms, 89 interferograms were analyzed; others could not be analyzed because they were blurred due to mirror vibration. Upon acquiring the unwrapped phase, ϕ(ξ,z), for all θl, we obtained a 3D projection data set for (ξ,z,θ). From the 3D data set, a 2D sinogram, pkl at (ξk,θl), was extracted by selecting a certain slice at z=zs; an example of pkl is shown in Fig. 8(c). As a measure for weighting, the wrapped difference, εkl at (ξk,θl) was similarly selected. For the weighted reconstruction, the weights for incident angles were given by “non-normalized” weights, Wl:

Wl(θ90°)=exp(SDl[εkl]2σW2),
Wl(θ=90°)=WFix,
Wl=Wl/lWl,
where SDl[εkl] denotes the standard deviation corresponding to samples εkl for a certain l, and WFix is either 1 or 10. An example of a frequency histogram of SDl[εkl] is shown in Fig. 8(b). Parameter σW in Eq. (30) was determined from this histogram as σW=0.4. The weight function is also plotted in Fig. 8(b), and the weight for each projection is shown in Fig. 8(d).

Figure 8(e) compares the reconstructed densities, fij, between different operators during the procedure of backward projection, as expressed by Eq. (9). When only the weight is considered (top row), i.e., when the moving average is not applied, fij exhibits a jagged appearance. When only the moving average operator is applied (middle row), penetrations of the tails are observed at the bottom, and the estimated peak intensity is smaller than in the other cases. In the case when both the operators applied for WFix=1 (bottom left), fij exhibits an inconsistency such that the order of the flame sizes is increasing in the order of upper left, lower, and upper right in Fig. 8(a); however, that of fij does not follow this trend. In contrast, the density for WFix=10 (bottom right) does not show such an inconsistency. Along with the simulations discussed in Section 2.B, these results demonstrate the validity of the moving average and the weighted reconstruction for noisy, limited-projection sinograms.

The internal refractive-index change, Δnij, can be estimated from the reconstructed density, fij, from a comparison between Eqs. (3) and (22) as

Δnij=fijk0Δs.
According to the Gladstone–Dale relation, ignoring the concentrations of CO2 and H2O, the gas temperature, T, can be expressed as a nonlinear function of the refractive index, n, as
T=n01n1T0,n=n0+Δn<n0,
where n0=1.0002765 under standard conditions (T0=290K). Figure 9 shows the reconstructed 3D temperature when both the moving average and the weight expressed in Eq. (30) and WFix=10 are applied. From the views in the xy plane, it can be understood that the temperature profile is reasonably accurate from comparisons of the flame sizes in Fig. 8(a). In the vertical distributions shown in the xz and zy profiles, the temperature field is tilted. This is due to the tilted positioning of the wick of the large flame. The width of the distribution decreases with increasing z, which can be explained by the Coandă effect [33]. In addition, we find that the temperature around the candle wick at z=0 is lower than that of the flames. From these points, we can conclude that the reconstructed results demonstrate qualitative correctness. However, there are certain problems. First, some artifacts still remain around the bottom half of the midsized flame in the xy views; these are unnaturally steep gradients or tails. These artifacts are due to the larger weight assigned for θ=90° over the others. Second, the reaction zone that should distribute in the thin hot layer cannot be observed due to insufficient spatial resolution. We suppose that this low resolution is due to the phase-unwrapping error for ill-conditioned interferograms. Furthermore, the quantitative temperature is different between the reconstructed one and the one measured with a thermocouple. The reconstructed maximum temperatures in the regions above the wicks of the small, middle, and large flames, were 824°C, 825°C, and 1245°C, respectively. Although the measurements with the thermocouple were inaccurate, they were still nearly 650°C, 740°C, and 810°C. In particular, the error between reconstruction and measurement is large for high temperatures. One of the reasons for this error is the nonlinearity present in the evaluation of the temperature. Equation (34) has Δn in the denominator, and T= at Δn=1n0. For the case where |Δn|<|1n0|, T is a monotonic function of Δn, and |dT/d(Δn)| increases with increasing |Δn|. Therefore, the error of T caused by the error of Δn grows with increasing T. This also can be understood from the scale bar in Fig. 9.

 figure: Fig. 9.

Fig. 9. Reconstructed 3D temperature distribution of gas around candle flames.

Download Full Size | PDF

4. CONCLUSION

The reconstruction algorithm proposed in this paper is based on two functions. The first is the weighted reconstruction in which the weight for each projection is determined from the error or uncertainty included in the projection. The second is the moving average operator utilized to reduce the jagged reconstructed distribution. These functions can be implemented with a few modifications of the backward-projection procedure in any iterative reconstruction algorithm, although we use the ML-EM algorithm as the base method. The validity of this reconstruction technique has been demonstrated through simulations.

From the viewpoint of gas temperature measurements around candle flames as an optical tomography application, the quality of the reconstructed 3D distribution is improved, and the qualitative result is reasonably accurate. However, errors still remain. To increase the reconstruction accuracy, we need to focus on improving the measurement system, which includes the use of other interferometric methods and longer wavelengths, improvement of accuracy in the phase-unwrapping process, and tuning of the weights. These will form the topic of our future tasks.

The weighted reconstruction can be applied for tomography in other fields. In this study on interferometric tomography, the error or the uncertainty of each projection was estimated from singular points in the wrapped phase as a measure to determine the weight. If an appropriate measure of the error is found, it can be applied to medical diagnostics in which an object to be measured is moving during measurement of projections.

Funding

Japan Society for the Promotion of Science (JSPS) KAKENHI (JP15K06097).

REFERENCES

1. L. A. Shepp and B. F. Logan, “The fourier reconstruction of a head section,” IEEE Trans. Nucl. Sci. 21, 21–43 (1974). [CrossRef]  

2. R. Gordon, R. Bender, and G. T. Herman, “Algebraic reconstruction techniques (art) for three-dimensional electron microscopy and x-ray photography,” J. Theoret. Biol. 29, 471–481 (1970). [CrossRef]  

3. P. Gilbert, “Iterative methods for the three-dimensional reconstruction of an object from projections,” J. Theoret. Biol. 36, 105–117 (1972). [CrossRef]  

4. L. A. Shepp and Y. Vardi, “Maximum likelihood reconstruction for emission tomography,” IEEE Trans. Med. Imaging 1, 113–122 (1982). [CrossRef]  

5. K. Lange and R. Carson, “Em reconstruction algorithms for emission and transmission tomography,” J. Comput. Assist. Tomogr. 8, 306–316 (1984).

6. A. H. Andersen and A. C. Kak, “Simultaneous algebraic reconstruction technique (sart): a superior implementation of the art algorithm,” Ultrason. Imag. 6, 81–94 (1984). [CrossRef]  

7. H. M. Hudson and R. S. Larkin, “Accelerated image reconstruction using ordered subsets of projection data,” IEEE Trans. Med. Imaging 13, 601–609 (1994). [CrossRef]  

8. J. A. Fessler, “Penalized weighted least-squares image reconstruction for positron emission tomography,” IEEE Trans. Med. Imaging 13, 290–300 (1994). [CrossRef]  

9. J. M. Anderson, B. A. Mair, M. Rao, and C.-H. Wu, “Weighted least-squares reconstruction methods for positron emission tomography,” IEEE Trans. Med. Imaging 16, 159–165 (1997). [CrossRef]  

10. T. Gerzen and D. Minkwitz, “Simultaneous multiplicative column normalized method (smart) for the 3d ionosphere tomography in comparison with other algebraic methods,” in Annales Geophysicae (European Geosciences Union, 2016), Vol. 34, pp. 97–115.

11. M. Beister, D. Kolditz, and W. A. Kalender, “Iterative reconstruction methods in x-ray ct,” Physica Medica 28, 94–108 (2012). [CrossRef]  

12. A. Iriarte, R. Marabini, S. Matej, C. O. S. Sorzano, and R. M. Lewitt, “System models for pet statistical iterative reconstruction: a review,” Comput. Med. Imaging Graph. 48, 30–48 (2016). [CrossRef]  

13. C. Shakher and A. K. Nirala, “A review on refractive index and temperature profile measurements using laser-based interferometric techniques,” Opt. Lasers Eng. 31, 455–491 (1999). [CrossRef]  

14. D. Naylor, “Recent developments in the measurement of convective heat transfer rates by laser interferometry,” Int. J. Heat Fluid Flow 24, 345–355 (2003). [CrossRef]  

15. G. W. Faris and R. L. Byer, “Beam-deflection optical tomography of a flame,” Opt. Lett. 12, 155–157 (1987). [CrossRef]  

16. H. El-Ghandoor, “Tomographic investigation of the refractive index profiling using speckle photography technique,” Opt. Commun. 133, 33–38 (1997). [CrossRef]  

17. H. M. Hertz, “Experimental determination of 2-d flame temperature fields by interferometric tomography,” Opt. Commun. 54, 131–136 (1985). [CrossRef]  

18. G. N. Blinkov, N. A. Fomin, M. N. Rolin, R. I. Soloukhin, D. E. Vitkin, and N. L. Yadrevskaya, “Speckle tomography of a gas flame,” Exp. Fluids 8, 72–76 (1989). [CrossRef]  

19. D. Yan and S. S. Cha, “Computational and interferometric system for real-time limited-view tomography of flow fields,” Appl. Opt. 37, 1159–1164 (1998). [CrossRef]  

20. S. Tomioka, S. Nishiyama, S. Heshmat, Y. Hashimoto, and K. Kurita, “Three-dimensional gas temperature measurements by computed tomography with incident angle variable interferometer,” Proc. SPIE 9401, 94010J (2015). [CrossRef]  

21. A. Savitzky and M. J. Golay, “Smoothing and differentiation of data by simplified least squares procedures,” Anal. Chem. 36, 1627–1639 (1964). [CrossRef]  

22. S. Tomioka, S. Nishiyama, and S. Heshmat, “Carrier peak isolation from single interferogram using spectrum shift technique,” Appl. Opt. 53, 5620–5631 (2014). [CrossRef]  

23. D. J. Bone, H.-A. Bachor, and R. J. Sandeman, “Fringe-pattern analysis using a 2-d fourier transform,” Appl. Opt. 25, 1653–1660 (1986). [CrossRef]  

24. S. Tomioka and S. Nishiyama, “Phase unwrapping for noisy phase map using localized compensator,” Appl. Opt. 51, 4984–4994 (2012). [CrossRef]  

25. R. M. Goldstein, H. A. Zebker, and C. L. Werner, “Satellite radar interferometry: two-dimensional phase unwrapping,” Radio Sci. 23, 713–720 (1988). [CrossRef]  

26. T. J. Flynn, “Two-dimensional phase unwrapping with minimum weighted discontinuity,” J. Opt. Soc. Am. A 14, 2692–2701 (1997). [CrossRef]  

27. M. A. Herráez, D. R. Burton, M. J. Lalor, and M. A. Gdeisat, “Fast two-dimensional phase-unwrapping algorithm based on sorting by reliability following a noncontinuous path,” Appl. Opt. 41, 7437–7444 (2002). [CrossRef]  

28. D. C. Ghiglia and L. A. Romero, “Robust two-dimensional weighted and unweighted phase unwrapping that uses fast transforms and iterative methods,” J. Opt. Soc. Am. A 11, 107–117 (1994). [CrossRef]  

29. R. Yamaki and A. Hirose, “Singularity-spreading phase unwrapping,” IEEE Trans. Geosci. Remote Sens. 45, 3240–3251 (2007). [CrossRef]  

30. S. Tomioka, S. Heshmat, N. Miyamoto, and S. Nishiyama, “Phase unwrapping for noisy phase maps using rotational compensator with virtual singular points,” Appl. Opt. 49, 4735–4745 (2010). [CrossRef]  

31. S. Heshmat, S. Tomioka, and S. Nishiyama, “Reliable phase unwrapping algorithm based on rotational and direct compensators,” Appl. Opt. 50, 6225–6233 (2011). [CrossRef]  

32. S. Heshmat, S. Tomioka, and S. Nishiyama, “Performance evaluation of phase unwrapping algorithms for noisy phase measurements,” Int. J. Optomechatron. 8, 260–274 (2014). [CrossRef]  

33. L. Pera and B. Gebhart, “Laminar plume interactions,” J. Fluid Mech. 68, 259–271 (1975). [CrossRef]  

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (9)

Fig. 1.
Fig. 1. Coordinate system and contribution from internal pixel to projection at a slice corresponding to z = const .
Fig. 2.
Fig. 2. Improvement of the reconstructed density via exclusion of an incorrect projection. (a) Projections for a square-shaped object along eight directions, where arrows denote the incident directions and the acquired projection data are indicated by the rows of pixels at the end of the arrows. The object for θ 4 is replaced by a rectangle instead of the square, i.e., p k , 4 M is incorrect. (b) Sinogram, p k , l M , and reconstructed density, f i , j ( n ) , at n = 100 when p k , 4 M denoted by the arrow is included. (c) Same as (b) but with p k , 4 M excluded.
Fig. 3.
Fig. 3. Simulation for limited projections: (a) True density with three Gaussian functions. (b)–(d) Conditions of the incident beams (top radial chart), and reconstructed densities at n = 10 (middle), and those at n = 50 (bottom). The projection beam angles, θ l , in (b) correspond to θ l [ 30 ° , 30 ° ] with 2° intervals. In both (c) and (d), a horizontal projection (90°) is added to the projections in (b). The weights, W l , for (b) and (c) are uniform, while W l ( 90 ° ) = 5 W l ( 90 ° ) in (d).
Fig. 4.
Fig. 4. Comparison of moving average for projection data with speckle noise. The true density and the condition of the projection beams are the same as those corresponding to Figs. 3(a) and 3(d), respectively. The speckle noise included in p k l M obeys the normal distribution, with the standard deviation being 10% of the maximum of the true projection data.
Fig. 5.
Fig. 5. Comparison of unweighted and weighted reconstructions for projection data with line error. The true density and the conditions of the projection beams are identical to the corresponding ones of Figs. 3(a) and 3(d), respectively. The random number, α l , to determine p k l M in Eq. (15), obeys the normal distribution with a standard deviation of 1. Parameter σ W in Eq. (16) is determined as max { p k l True } / σ W = 0.1 .
Fig. 6.
Fig. 6. Experimental arrangement. FM1-5, full-mirrors; HM1-5: half-mirrors; X1 and X2, linear stages; LS, He–Ne laser (wavelength: 633 nm); L1 and L2, lenses (expanded beam diameter: 50 mm); OBJ, object; S , half-transparent screen. FM1, FM2, HM1, and HM2 are mounted on rotatable stages, which in turn are mounted on X1 and X2.
Fig. 7.
Fig. 7. Process of phase evaluation from an interferogram. (a) Observed interferogram ( θ = 20.0 ° ). (b) Modified interferogram obtained by both normalization and the Wiener filter. (c)–(f) Spectral intensities in the spectral domain limit k [ k max / 2 , + k max / 2 ] , where k max denotes the Nyquist frequency; I ^ twin ( k ) = I ^ 1 ϕ ( k k ) + [ I ^ 1 ϕ ( k κ ) ] * in (e). (g) Wrapped phase with range [ π , + π ] . (h) Grouped tree of singular points to determine compensated local area for phase unwrapping. (i) Unwrapped phase with contour lines with interval 2 π . (j) Wrapped difference between rewrapped phase of (i) and wrapped phase in (g). (k) Additional phase shift to satisfy equi-phase condition on both sides of the interferogram. (l) Corrected unwrapped phase. Image sizes of (a) and (b) are 128 × 128 pixels, and those of (g)–(l) are 512 × 512 pixels.
Fig. 8.
Fig. 8. Reconstructed density from experimental data with different operators for a slice ( z = 10 mm ). (a) Object with three different-sized flames on a candle. (b) Frequency histogram of SD l [ ε k l ] and weight function, W l , for weighted reconstruction. (c) Sinogram, p k l ; the horizontal and vertical axes represent k and l , respectively. (d) Weight, W l ; the horizontal axis represents W l within [0,1] and the vertical axis is the same as that in (c). (e) Reconstructed densities obtained with application of different operators in the backward-projection procedure; WGT denotes weighted projection and M.A. moving average. Intervals between contour lines are 2 × 10 5 for Δ n .
Fig. 9.
Fig. 9. Reconstructed 3D temperature distribution of gas around candle flames.

Equations (34)

Equations on this page are rendered with MathJax. Learn more.

p k l ( n ) = F k l { f i j ( n ) } ,
f i j ( n + 1 ) = B i j { p k l M p k l ( n ) } f i j ( n ) ,
p k l = 1 Δ s ξ k Δ s / 2 ξ k + Δ s / 2 η l f ( x ( ξ , η ) , y ( ξ , η ) ) d η d ξ .
F k l { f i j } = i , j c i j k l f i j .
A = k C k l , C k l = i , j c i j k l , L k l = C k l Δ s , Δ s 2 = k c i j k l ,
B i j { p k l } = 1 N l Δ s 2 k , l c i j k l p k l L k l ,
B i j { p ^ k l } = B i j { min ( p ^ k l , p ^ lim ) } ,
f i j ( n + 1 ) = B i j { p k l M p k l ( n ) } f i j ( n ) ,
f i j ( n + 1 ) = M { B i j { W l p k l M p k l ( n ) } f i j ( n ) } ,
l W l = 1 .
p k l M = p k l True + N k l , p k l True = F k l { f i j True } ,
Δ ^ { g p q ( n ) , ref } = | | g p q ( n ) g p q ref | | 2 | | g p q ref | | 2 , | | g p q | | 2 = g p q 2 p q ,
p k l M = p k l True + p k l δ .
f i j ( ) = F i j 1 { p k l True } + F i j 1 { p k l δ } = f i j True + F i j 1 { p k l δ } ,
p k l M = ( 1 + α l ) p k l True .
W l exp ( ( p k l M p k l True ) 2 k / σ W 2 ) ,
I ( r ) = I 0 ( r ) + 2 I 1 ( r ) cos ( Δ ϕ ( r ) ) ,
Δ ϕ ( r ) = ϕ obj ( r ) ϕ ref ( r ) ,
ϕ obj ( r ) = η k 0 obj n obj ( r ) · d r ,
ϕ ref ( r ) = η k 0 ref n ref · d r ,
Δ ϕ ( r ) = ϕ ( r ) + κ · r ,
ϕ ( r ) = η k 0 obj Δ n ( r ) · d r = k 0 η Δ n ( r ) d η ,
Δ n ( r ) = n obj ( r ) n ref ,
κ = ( k 0 obj k 0 ref ) n ref .
I ^ ( k ) = I ^ 0 ( k ) + I ^ 1 ϕ ( k κ ) + [ I ^ 1 ϕ ( k κ ) ] * ,
I ^ 1 ϕ ( k ) = F { I 1 ( r ) exp ( j ϕ ( r ) ) } ,
ϕ w ( r ) = W { ϕ ( r ) } = arg ( I 1 ( r ) exp ( j ϕ ( r ) ) ) ,
ϕ ( r ) = U { ϕ w ( r ) } .
ε ( r ) = W { ϕ ( r ) ϕ w ( r ) } .
W l ( θ 90 ° ) = exp ( SD l [ ε k l ] 2 σ W 2 ) ,
W l ( θ = 90 ° ) = W Fix ,
W l = W l / l W l ,
Δ n i j = f i j k 0 Δ s .
T = n 0 1 n 1 T 0 , n = n 0 + Δ n < n 0 ,
Select as filters


Select Topics Cancel
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.