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 [2–10], 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 [15–17] 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, , must always be or always .
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:
where 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 expresses the projection data set called the “sinogram,” which corresponds to the phase shift; the superscripts and denote the iteration number and measured signals, respectively; the subscripts , , , and denote the indices of the discrete coordinates of and , and and represent the forward-projection and backward-projection operators, respectively.The projection data, , which is detected at the th sensor with width for the th incident beam angle, can be expressed as an average over the axis of the integral along the axis of :
Let us assume that the grid size in space is also equal to . By denoting the contribution from the pixel at to the sensor at as , the forward-projection operator in the discrete form can be written as The contribution represents the overlapped area of the th pixel and the th sensor. The following equations depict the relations between and some geometrical constants: where represents the entire area of the domain to be reconstructed, the area of the beam with width incident on the th sensor, and the path length.The backward projection is evaluated as an angular average of , where is an averaged density along the axis:
where denotes the number of the beam angles. For a uniform averaged density, , the summation on the right-hand side becomes from Eq. (5). The factor is determined such that the right-hand side becomes .In Eq. (2), the backward projection is applied to the ratio of the measured projection, , to the current estimated projection, ; subsequently, the new density, , is evaluated as the product of the result of the backward projection and the current density, . If the current projection is overestimated, i.e., , the ratio becomes smaller than unity, and the new density, , is smaller than the current density, .
Because is applied to the ratio, the iterative procedure becomes unstable if . In our implementation, the density is evaluated by the following modified backward-projection operator, , which can limit the upper bound of the ratio to avoid unstable iterations:
where represents the upper limit of the ratio, which is taken as 2 in our implementation. In regards to the initial estimation of the density, , in Eq. (1), we use , 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
where the denotes a moving average operator, and denotes a weight for that satisfies The purpose of utilizing the moving average is to avoid convergence to incorrect solutions. Meanwhile, the purpose of the weight, , 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.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).
2. Simulation for Sinogram with Speckle Noise
There are two types of errors in the sinogram , 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 .
The random noise is caused by speckle noise or electrical noise. The sinogram with noise is represented as
where denotes the random noise. Figure 4 demonstrates a simulation for speckle noise, in which the performances of the moving average, , are compared. We employed a simple moving average in which the density is convoluted with a uniform kernel with a square window. The definition of the relative error plotted in Fig. 4(c) is where denotes the average with respect to and , and appearing in Fig. 4(c) is . From Fig. 4(c), it can be observed that the error profiles of the projection, , decrease monotonically with increasing , regardless of whether or not the moving average is applied, and they converge to a saturated level when . However, the error profiles of the density, , are different. When the moving average is applied, decreases in a manner similar to , and the results corresponding to and shown in the right-hand panels of Fig. 4(b)are similar. In contrast, when the moving average is not applied, reaches a minimum value at , and it increases with becoming noisy with increasing . These characteristics can be understood from the principle of ML-EM reconstruction. The original ML-EM minimizes with increasing . Let us consider the case where can be expressed as the sum of the true projection, , and the error caused by noise, :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 is zero for the true solution. Under this assumption, the solution can be represented as
where denotes the inverse operator of . This equation indicates that the reconstructed density includes the reconstruction of the noise, , as well as the true density. When the noise in the projection data is random, may exhibit a random distribution with higher-frequency components. Because the moving average can reduce the higher-frequency components [21], the contribution of is reduced. It should be noted that the higher-frequency component in is also reduced by the moving average; therefore, the spatial resolution of the reconstructed density is degraded.3. Simulation for Sinogram with Line Error
The line error is introduced in the procedures for each value because of the following reason. The sinogram, , is a 2D data “slice” on the plane represented as , which is obtained by selecting a certain slice on from the 3D projection data set in the coordinate system. The 3D data are obtained by acquiring 2D projection data on the plane as for different . Because is processed from a measured image, , for each , the error of the sinogram, , does not have any correlation between different . This means that the error of becomes a discontinuous function even when the error of is a continuous function. The sinogram with line error for a simulation is generated by using random numbers, , as
Figure 5 demonstrates the validity of weighted projection for the line error. The weight for each in the weighted reconstruction is expressed as where 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 for the weighted reconstruction decreases with increasing , and for the unweighted reconstruction has a local minimum. In the case of the unweighted reconstruction, the maxima of and are 30% and 70% larger than , respectively. Among the reconstructed densities in Fig. 5(b), for the weighted reconstruction is the best solution, whose maximum is 10% larger than .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 , whereas 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 is given as a known function to determine , the true projection, , required to determine , is also given. In actual applications, however, 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 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 . 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 , 3 mm, and 4 mm.
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, where , on the screen is represented as
where and represent intensities related to incident laser profile; , , 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 and are directed in different directions, but , where denotes the wavenumber in vacuum. From these equations, can be rewritten as In the case without the object, because and , a parallel fringe pattern with carrier frequency, , appears in the interferogram. In the case with the object, the carrier pattern is modulated by , which is the 2D phase shift distribution by the object.The phase shift, , can be estimated as follows. The Fourier transform of in Eq. (17) distributes to three peaks:
where denotes a 2D Fourier transform operator with the integral kernel , where denotes the imaginary unit and the spectral domain; “^” refers to the variables in 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 , and applying the inverse Fourier transform, we obtain as a complex-valued image. Because denotes a real-valued image, we can obtain a principal value of phase of by taking the complex argument where 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 , to solve the ambiguity, a phase unwrapping operator, , is applied:2. Noise Reduction
An example of the phase evaluation process is shown in Fig. 7. The observed interferogram, , is shown in (a). We note that the observed interferogram, , has a spatial profile with a low frequency, and the brightness at either end is low. The spectrum intensity, , shown in (c) has three bright narrow peaks. The widths of the peaks are mainly determined by the profile and 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, was normalized by the measured laser source profile, and a Wiener filter was subsequently applied. The result, , and its spectral intensity, , 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 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, , is shown in Fig. 7(e), and the shifted carrier peak, , 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, , 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 [24–31], 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, , 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, 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, , 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, , 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, , is not identical to the original wrapped phase, , as shown in Fig. 7(j), in which is defined as
We employed the standard deviation of as a measure of the quality of the unwrapped phase. Although is also a candidate of the measure to evaluate the quality, the dependency between and the quality of was small if the factors to represent 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, , for all , we obtained a 3D projection data set for . From the 3D data set, a 2D sinogram, at , was extracted by selecting a certain slice at ; an example of is shown in Fig. 8(c). As a measure for weighting, the wrapped difference, at was similarly selected. For the weighted reconstruction, the weights for incident angles were given by “non-normalized” weights, :
where denotes the standard deviation corresponding to samples for a certain , and is either 1 or 10. An example of a frequency histogram of is shown in Fig. 8(b). Parameter in Eq. (30) was determined from this histogram as . 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, , 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, 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 (bottom left), 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 does not follow this trend. In contrast, the density for (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, , can be estimated from the reconstructed density, , from a comparison between Eqs. (3) and (22) as
According to the Gladstone–Dale relation, ignoring the concentrations of and , the gas temperature, , can be expressed as a nonlinear function of the refractive index, , as where under standard conditions (). Figure 9 shows the reconstructed 3D temperature when both the moving average and the weight expressed in Eq. (30) and are applied. From the views in the 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 and 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 , which can be explained by the Coandă effect [33]. In addition, we find that the temperature around the candle wick at 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 views; these are unnaturally steep gradients or tails. These artifacts are due to the larger weight assigned for 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 in the denominator, and at . For the case where , is a monotonic function of , and increases with increasing . Therefore, the error of caused by the error of grows with increasing . This also can be understood from the scale bar in Fig. 9.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]