Many techniques allow the measurement of physical properties based on the retrieval of phase information encoded in an interference pattern. Techniques such as profilometric1,2 and interferometric3,4 methods measure mechanical properties (e.g., train or deformation) of materials. In the current paper, we focus on digital holography interferometry techniques.5 Digital holography is a technology of obtaining and processing holographic measurement data, typically via a charge-coupled diode (CCD) camera or a similar device. However, a major problem with digital holography and other interferometric techniques,6 which recover phase information, is that the reconstructed phase is mathematically limited to the interval . In general, the true phase may range over an interval larger than , which implies that the obtained phase may contain discontinuities. Therefore, unwrapping these discontinuities in the phase is required. Phase unwrapping is a technique of finding the true phase values on the assumption that the unwrapped phase map is a continuous surface. However, the process of phase unwrapping cannot avoid complexity due to a number of serious issues. For example, the presence of noise or singular points (SP), e.g., rising from regions of low fringe contrast in the hologram, is a primary concern for effective unwrapping of the reconstructed phase. Noise impedes the ability to judge the existence of a phase jump, while SPs mean that the unwrapped phase will be path-dependent.
Although various methods to estimate a correct unwrapped phase map have been proposed, they can be divided into three categories. The first category contains algorithms based on following the paths.188.8.131.52.12.13.–14 These methods search for SPs, then pair these SPs by placing branch cuts. By examining the branch cuts and determining if any appear to be placed poorly or if any isolate a region, it can be determined whether or not the paths can be followed to retrieve the phase maps, as well as whether these methods succeed or fail. The second category includes methods which use the least-squares approach.15184.108.40.206.20.–21 These algorithms use the same idea of minimization of discrete gradients difference squares as used in the leased-squares approach. These differences are taken between the wrapped phase gradients and supposed unwrapped phase gradients. Algorithms based on spreading singularity2223.24.–25 are also classified into the same category of least-square methods since they spread the singularity like the least-square methods do. In terms of accuracy, the method using localized compensator (LC)25 is superior to the other methods. However, it has a major disadvantage of computational cost since this method has a high time cost to produce its unwrapped results. The last category is denoising-unwrapping methods.2627.28.–29 This type of methods performs phase map denoising to remove noise from wrapped phase by using a filtering process. The filtering process is sometimes applied before the unwrapping process such as windowed Fourier transform method,26 or it is applied simultaneously with the unwrapping process such as dynamic filter method.27,28 These methods can reduce the noise within the original spatial resolution. However, they have a minimum signal-to-noise (S/N) ratio for the wrapped phase data in which these methods start to fail to obtain accurate unwrapped phase result. In addition, the unwrapped results of these methods are highly depending on the relaxation parameter, which can control the cut-off frequency of the filter.
Based on the above discussion, existing phase unwrapping methods suffer from various problems that can affect the time cost and the accuracy of unwrapped results. Therefore, investigation into and improvement upon existing unwrapping methods are needed. In this paper, an effective and fast phase unwrapping algorithm is presented to be a powerful tool in digital holography measurements. It was found that the distribution of SP dipole distances shows that there are many dipole pairs a short distance apart.10 According to this finding, we have introduced a new technique to compensate the singularities at adjoining SP pairs, which is a direct compensator (DC).24 The method uses DCs for adjoining SP pairs, and rotational compensators (RC)23 for other SPs. By applying our algorithm, the phase inconsistency is canceled, and then we can unwrap the phase data simply by summing the phase differences between neighboring pixels. Furthermore, in the fringe analysis of digital holographic measurement, we need to remove the effect of background phase modulation. There are two ways to eliminate background fringe from the measured data. In this paper, we present these ways to extract the phase information about the object in hologram maps.
The paper is organized as follows. Section 2 presents the basic principles about experimental setup system of digital holography and phase measurement. Moreover, Sec. 2 introduces phase singularity and unwrapping. In addition, methods to extract phase information about the object from hologram maps are also shown. Section 3 reviews a description of the developed phase unwrapping algorithm based on singularity compensation for the holographic measurements. Meanwhile, Sec. 4 gives a performance comparison of the studied methods for the measurement of hologram data and the object extraction methods. Lastly, a conclusion of this paper is drawn in Sec. 5.
Phase Measurement and Object Extraction for Digital Holographic Data
Basic Principles of a Digital Holography System
Holography is a method of recording the phase modulation from the light reflected or transmitted from a projected object to screen in the form of an interference pattern. The optical setup of a holographic system is shown in Fig 1. For recording a digital hologram, light from a laser is divided by a beam splitter into beams, one for object illumination and another for reference. The object beam illuminates the object, and the illuminating light is scattered back by the object toward the detector, where it forms an image of the object on a CCD camera. A hologram is a fringe pattern formed on the CCD as a result of the interference between the reference beam, , and the object beam, . The intensity recorded on the CCD is given by
The last two terms of Eq. (1) contain information about the amplitude and the phase of the object. To retrieve the phase information, two techniques are partly used in digital holography. One is the phase shift interferometry,3031.–32 in which several fringe patterns are recorded by varying the known phase shift introduced to one of the beams in the interferometer system. The other one is spatial filtering for the hologram using the Fourier transform method.33,34 The Fourier transform method requires only one fringe pattern; however, the high frequency part of spectrum in the Fourier domain cannot be used. Fortunately, in cases where the object is heated gas, the fringe patterns do not have high frequency components.35 For this reason, we apply the Fourier transform method to extract phase data in our experiment.23 A schematic diagram describing the use of the Fourier transform method is shown in Fig. 2. In the figure, the Fourier spectrum, , shows a symmetrical distribution of the origin. The is the filtered spectrum in which the zero-frequency component and one of the symmetrical distributions are eliminated. The center of gravity of the filtered spectrum can be also obtained. There is a problem in the filtering of the spectrum, namely that we normally cannot determine which nonzero spectrum should be eliminated. However, even when the eliminated term is incorrect, the sign of the phase modulation is only inverted. Therefore, if we know the sign of the phase modulation, we can determine it after the phase unwrapping. After filtering and inverse Fourier transformation, the complex amplitude of the wave front is obtained. From the complex digitalized amplitude, the phase of the wave front is calculated by the relation
SPs and Phase Unwrapping
A major problem with digital holography and other interferometric techniques that recover phase information is that the reconstructed phase is mathematically limited to the interval . Therefore, an unwrapping process is needed to apply. Phase unwrapping is the process of retrieving the absolute phase data from its wrapped phase. In general, phase unwrapping is an ill-posed problem. Hence, the existing noise during measurements or phase discontinuities resulting in SPs makes the process to recover unambiguous phase data more complicated. However, certain assumptions of underlying process can solve the phase unwrapping problem. The most common assumption is that the true unwrapped phase data changes slowly enough to make the phase difference values of neighboring points within one half cycle ( radian), which means . The relationship between the wrapped phase, , and the unwrapped phase, , can be stated as
The two-dimensional (2D) phase map consists of a sequence of horizontal and vertical segments joined at their adjacent points. SPs are defined as the local inconsistencies in the phase map, and they mark the beginning and the end of discontinuities. To calculate the SPs, a closed path was applied which had four segments () starting in every point defined by the corners of a square. These SPs are identified when the summation of the wrapped phase gradient, , along the closed path in the clockwise direction have nonzero value, as follows:4) is or . The SP is called a positive residue when is ; otherwise, it is called a negative residue when is . Basically, SPs appear as a dipole, which are pairs of SPs with the opposite sign. However, some isolated SPs or SPs of dipoles with long distances, which are called monopoles, may appear near the boundaries of the phase map because the measurement domain is finite. This can lead to that the numbers of positive SPs and negative SPs in the measurement area are different. Monopoles spread errors throughout the entire measurement area.23 Thus, the difficult problems in phase unwrapping process are, how to deal with the SPs, and to process their singularities in the phase data, in Sec. 3, we try to solve these problems.
Phase Extraction Methods of Object Information from Holograms
In the holographic measurement system, two fringe patterns should be measured to produce information on an object by using a similar setup to the experiment shown in Fig. 1. First, a fringe pattern was obtained for the existence of an object. This fringe is referred as an object fringe and it is a superposed result of the object light passing through the object upon the reference light. The other fringe pattern is measured for background, which is the result of the same system but without an object. To compute the information about the object, the background fringe from the measured data needs to be eliminated. There are two methods to extract the phase shift caused by the object from the measured data, and these methods depend on the time of excluding the background. Figure 3 explains schematic diagrams to compute the phase shift of an object from experimental data. The first method was the pre-rejection of background data by subtracting the unwrapped phase of background data obtained without the existence of the object, from the wrapped phase data obtained with the existence of the object. Then the phase difference was unwrapped to get the unwrapped phase shift result, as shown in Fig. 3(a). Meanwhile the other method, the post-rejection of background data was carried out, as illustrated in Fig. 3(b). In this method, the phase difference is computed by excluding the unwrapped background phase data from the unwrapped phase of the object. Section 4 examines the effect of the extraction methods on the unwrapped results of the phase shift caused by the object for digital hologram maps.
Phase Unwrapping Based on Rotational and Direct Compensation
Basic Concept of Compensator
Phase unwrapping is an essential process of removing inconsistencies by local neighborhood tests and corrections to produce continuous phase maps. In order to solve inconsistencies caused by SPs, we have presented our proposed algorithm based on compensating the phase singularities to cancel their effect.24 The idea of singularities compensation to remove their effect in the phase map was proposed before, in the singular-spreading phase unwrapping (SSPU) method22 and in methods using a RC23 or a LC.25 However, these methods compensate the singularities in different ways. The SSPU method requires an iteration process to compute the compensators. On the other hand, the RC method can compute the compensator through superposing the effect of each SP by adding an integral of isotropic singular function along any loop. Meanwhile, the LC method regularizes the inconsistencies only in a local areas, which are clusters, around the SPs by integrating the solution of Poisson’s equation for each cluster to compute the compensators. In terms of accuracy, the method using LC25 is superior to the other methods. Despite this, LC method has the major disadvantage of computational cost since this method requires a long time to generate cluster groups and to compute the compensators. The algorithm similarly compensates the inconsistencies and confines the effect of each one in a local region. The main issues which determine the behavior of the algorithm are similar to RC method.23 However, the way of computing the compensators for adjoining SPs pairs in the method is different from in RC.
The RC method uses local phase information to compensate the singularity parts of phase map caused by existence of SPs. The RC method is based on three techniques: an RC, unconstrained singular point (USP) positioning and virtual singular points (VSP). The RC technique acts to compensate the singularity of each SP for all unwrapping paths. The USP provides freedom to adjust the SP positioning in order to improve the accuracy of compensation. Since it can make some dipoles of SPs that have shorter distances than the pixel size, the undesired longer effect of the compensator is suppressed.23 The VSPs with an opposite sign of the isolated SPs are put outside the measurement area, so that these VSPs and the isolated SPs make dipoles. Thus, the error caused by the existence of isolated SPs may be reduced, and the effect of the compensator is confined in local narrow regions around SPs. However, the RC method has a drawback of an undesired phase error because the RC should be applied to the regular region with no SPs as well as to the singular region. To reduce this error, we have applied a new method for a SP dipole pair that has a short distance.
Every SP has a residue of , and a pair of two SPs with different polarities is considered as a dipole. It was found that the distribution of SP dipole distances shows that there are many dipole pairs with short distances.10 Figure 4 shows a distribution of dipole distances. It can be seen that many dipoles are distributed around one pixel distance. Based on this finding, we have proposed the method. This algorithm reduces the drawbacks of the RC method, which are high computational time cost and undesired phase error due to its effect on the regular region. The proposed method compensates the singularities of adjoining SP pairs by adding a DC. Thus, the effect of each SP is confined in a closer local region. If the distance between two SPs with opposite polarity is one pixel, the DC can be applied. The sum of the DC along the smallest path surrounding one of the two SPs, which consists of four segments, equals one cycle ( radian). Furthermore, the sum of the DC along the path surrounding both SPs must vanish. A solution satisfying these conditions is obtained by setting the branch cut. The branch cut is placed between the two SPs. When the segment to unwrap crosses the branch cut, the sum of DC for the two SPs is defined as one cycle. No DC is applied for the other segments.
Description of the RC + DC Method
The method is coupling RC and DC methods to compute the compensators depending on the SPs locations. In other words, it uses DC as compensators for the adjoining SP pairs, and uses RC to compute the compensators of nonadjoining pairs. The adjoining pair is a dipole which consists of two SPs with opposite signs, separated by one pixel horizontally or vertically. Figure 5 shows the configuration of the branch cuts placed between adjoining SPs in the phase map and the concept of the direct compensation. Figure 5(a) shows a case in which the branch cut is placed horizontally between a pair of adjoining SPs, so that the DC will be added to the vertical segment that crosses the short branch cut. In contrast, Fig. 5(b) shows the case in which the branch cut is placed vertically between the adjoining SPs and the DC is added horizontally. The compensator value of the segment is divided into two, and distributed through the two adjacent loops, which contain adjoining SP pairs, as illustrated in Fig. 5. The direction of the DC for the segment is based on the position of each segment with respect to the location of the tested SP. Therefore, the confinement of the DC effect in a closer region around the SPs leads to the improvement in the accuracy of the unwrapped phase results. When the segments are far from SPs, their DCs have zero value so that the computation of DC is not needed. In contrast, the RC requires computation for all segments. For this reason, the computation time requirements of the algorithm for computing total compensators will be reduced and the accuracy of the unwrapped phase will be improved, as discussed in Sec. 4. After computing the total compensators for each segment, the true unwrapped phase values can be retrieved by summing the phase differences between the adjoining pixels and the total compensators, as follows:
The steps of the algorithm can be summarized as follows:
1. Calculate SPs in the wrapped phase map, as defined in Eq. (4).
2. Append VSPs for the monopoles SPs, then analyze the SP pairs. After that, define the adjoining SP pair positions.
3. For each SP ,
a. If the ’th SP is one of the adjoining pair SPs: The ’th segment position related to this SP is determined, then DC of this SP with respect to the ’th segment, , will be used, and then it will be added to the total compensator of this segment, . The value of is and its sign is dependent on the position of the segment with respect to the location of the tested SP. Meanwhile, the other computations of compensators for this ’th SP respect to other segments are zero, therefore, these computations will be skipped.
b. If the ’th SP is not related to any adjoining pairs: RC, , will be computed for this SP with respect to each segment in the phase map. Then these values of compensators of this tested SP will be added to the total compensators of each segment, .
4. Finally, the unwrapped phase data can be retrieved by adding the compensators, to wrapped phase differences by using Eq. (5) where is the total compensator by DC and RC, which regularizes the singularity of .
This description of direct compensation for adjoining SPs pairs makes the algorithm simple and easy to implement. It provides a fast and efficient way to unwrap the phase map.
Results and Discussion
In this section, two examples of noisy wrapped phase maps are presented. One is a simulated phase map where the true phase is known to evaluate the accuracy of the compared methods. The other is the experimental data obtained with interferometer to demonstrate the performance of the method for holographic data.
Computer Simulation Results
In order to demonstrate applicability of the algorithm, a simulated noisy phase map which has Gaussian distribution shape is generated. This phase data has image size of and the standard deviation of noise is 0.15 cycle. The original and wrapped phase data are shown in Fig. 6(a) and 6(b), respectively. The wrapped phase data has 675 SPs, and the numbers of positive and negative SPs are 338 and 337, respectively. Their distribution map is shown in Fig. 6(c). In addition, the unwrapped phase results obtained by Goldstein et al.’s path-following method,7 Flynn’s minimum weighted discontinuity algorithm,11 the least-square method with discrete cosine transform (LS-DCT),21 and the algorithm are shown in Fig. 6(d), 6(e), 6(f), and 6(g), respectively. Also, the rewrapped phase results of these algorithms are illustrated in Fig. 6(h), 6(i), 6(j), and 6(k). It is clear that Goldstein’s algorithm gives poor accuracy in the unwrapped result. It can be also observed that the rewrapped result of LS-DCT method is not identical to the wrapped data. Compared to that, the algorithm and Flynn’s method give satisfactory results. In addition to this data, Table 1 presents a quantitative comparison for the unwrapped results, in terms of whether or not there are phase jumps in the unwrapped results, and in terms of the standard deviation () of the difference between the true phase and the unwrapped phase obtained by each algorithm. It can be observed that there are phase jumps in the unwrapped results obtained by Goldstein et al.’s method, while the unwrapped results of the other methods do not have phase jumps. In addition, Table 1 shows that the smallest error in terms of is found for Flynn’s result and is followed by the result of the algorithm. The computation time required for each algorithm to obtain the unwrapped results is also shown in Table 1. The computational time for each phase unwrapping algorithm is measured using a PC with Intel Core 2 DUO CPU installed, with 2.13 GHz clock in a single CPU operation mode. The computing language used to implement the compared phase unwrapping algorithms is C language. From Table 1, the highest time cost can be found for Flynn’s method. The results in Fig. 6 and Table 1 explain that the method gives an unwrapped result with acceptable quality in both accuracy and low computational time cost.
Comparison of the accuracy of the unwrapped phase results corresponding to the execution time cost for each algorithm.
|Algorithm Name||Phase Jump||σ (cycle/pixel)||Required Time (s)|
|Goldstein et al.||Exists||0.33||0.08|
|Flynn||Does not exist||0.08||12.00|
|LS-DCT||Does not exist||0.35||0.15|
|RC+DC||Does not exist||0.14||0.48|
The algorithm has also been tested experimentally on a 2D wrapped phase map that resulted from the analysis of real fringe pattern. The object of this experiment is the temperature measurement of the heated gas (air) around a candle flame through measuring the phase shift caused by the flame using a Mach–Zehnder interferometer.23 The setup of the experiment is similar to the one shown in Fig. 1. The fringe pattern obtained in existence of candle is shown in Fig. 7(a); it is referred as object fringe. The phase data has an image size of and number of SPs is 2532 (1267 positive SPs and 1265 negative SPs). The wrapped phase data for the object and its corresponding SPs distribution map are shown in Fig. 7(b) and 7(c), respectively. In this measurement, the exposure time cannot be set long enough, because the flame varies in time by convection flow around the flame itself. For this reason, the exposure time is set to 1 ms. This setting cause two problems: firstly, the fringe has low S/N; secondly, we cannot apply the phase shift techniques3031.–32 that use several fringes with different reference lights to obtain a wrapped phase. Spatial filtering for the hologram using Fourier transform method is hence applied to extract the phase information, as shown in Fig. 2. In addition, the background phase map in this experiment is obtained by fitting a planar function by using information from the wrapped phase data extracted from the object fringe pattern. This information is taken from the area where the object light did not pass through the flame in the wrapped data.
To obtain the phase shift caused by the candle flame from the measured data, two methods to eliminate background can be applied, as shown in Fig. 3. Figures 8 and 9 show the unwrapped and rewrapped results of the phase shift of the candle flame depending on the extracting way of the object. Figure 8 presents a comparison of the accuracy of the examined phase unwrapping algorithms’ results of the candle flame for pre-rejection of the background way to extract the object. Meanwhile, Fig. 9 shows the compared results for the post-rejection way. From the figures, it can be found that the unwrapped result of the Goldstein et al. method causes phase jumps; however, the other three methods have no phase jumps. Although there is no phase jump in the unwrapped results obtained by the LS-DCT method for both ways (pre-rejection and post-rejection) of background, its rewrapped results produced by both ways are different, as shown in Figs. 8(c) and 9(c). This implies that the accuracy of LS-DCT method remains in doubt. On the other hand, the rewrapped phase results in both ways for object extraction, which are pre-rejection and post-rejection, are quite similar for both the Flynn method or the algorithm, as shown in Figs. 8 and 9. Therefore, it can be said that Goldstein et al.’s method and the LS-DCT method provide inaccurate phase results. Meanwhile, the Flynn method and the algorithm produce accurate results. However, the unwrapped results of the examined algorithms are affected by the manner of object extraction. It is understood that in the way of post-rejection for the background, Goldstein et al., Flynn and LS-DCT methods give better results than the pre-rejection way does. The reason is that the number of SPs from the wrapped phase data in post-rejection way, which is 2532 for the studied unwrapping algorithms, is smaller than its number of SPs in the pre-rejection way, which is 3046 for Goldstein et al’s method, 3208 for the Flynn algorithm and 2690 for the LS-DCT method. In contrast, the unwrapped phase result obtained by the method in the way of pre-rejection for the background is better than its unwrapped result obtained in the post-rejection way. This is due to that the ratio of adjoining pair of SPs for the wrapped phase data in pre-rejection way, which is 70.41%, is larger than its ratio of 60.54% in post-rejection way. This implies the benefit of the algorithm which uses DC to compensate the singularities of adjoining pair of SPs to reduce the unwrapping error.
In addition, the execution time required for each studied algorithm to obtain the unwrapped results is also evaluated. It is found that the highest time cost for producing the unwrapped results in the both ways of object extraction are belongs to the Flynn method, which is 736.60 s in the pre-rejection way and 450.90 s in the post-rejection way. In the meantime, the execution time of the algorithm to obtain its unwrapped result for both ways of the object extraction showed better performance than the Flynn method, being 8.84 s in pre-rejection way and 8.85 s in post-rejection way. Hence, it can be concluded from the above discussion that the method gives results of acceptable quality with low computational time costs.
We have presented a phase unwrapping algorithm for digital hologram measurements. The method is coupling the RC method and the DC method. DC compensates the singularity of adjoining pair SPs that are connected by branch cut. The compensator along the segment that crosses the branch cut is just . The performances of our developed phase unwrapping algorithm and of other existing phase unwrapping methods for digital holographic data are compared. In addition, the methods to extract phase information about the object from hologram maps are also investigated. The results show that the algorithm gives results with acceptable quality and with low computational time costs compared to the existing methods.
This research was supported by Japan Society for the Promotion of Scientific Research (c), 24560216, 2013.
Samia Heshmat received her MEng degree in quantum science engineering from the Graduate School of Engineering, Hokkaido University, in 2010. She joined the Faculty of Engineering, at Aswan University in 2005. At Aswan University, she worked as a teaching assistant for undergraduate students and conducted research on image processing. Currently, she is a doctoral student in the Division of Quantum Science and Engineering, Graduate School of Engineering, Hokkaido University. Her current research interests include image processing, phase unwrapping, interferometry data.
Satoshi Tomioka is an associate professor at Hokkaido University. He received his BEng, MEng, and DrEng degrees in the field of nuclear engineering from Hokkaido University, Sapporo, Japan, in 1986, 1988, and 1996, respectively. His current research interests include plasma diagnostics, boundary element method for electro-magnetic field, genetic algorithm, three-dimensional measurement of refractive index distribution.
Shusuke Nishiyama received his BEng degree in nuclear engineering from the Department of Engineering, Hokkaido University, Sapporo, Japan, in 1994, and his MEng degree in nuclear engineering and his DrEng degree in quantum energy engineering from the Graduate School of Engineering, Hokkaido University, in 1996 and 2000, respectively. He is currently an assistant professor at the Division of Quantum Science and Engineering, Graduate school of Engineering, Hokkaido University. His research interests include plasma diagnostics and numerical analysis of electromagnetic field.