In this paper, a simple isoclinic phase map unwrapping method is proposed to retrieve map with ambiguities at photoelastic isotropic points. Regional phase unwrapping method is also utilized to enhance the retrieving efficiency after all phase inconsistencies have been fully detected and branch cutting works have been properly done to ensure blockings of all the paths which could cause incorrect integrations while involuntarily crossing them. The correctly retrieved isoclinic data are then fed into isochromatic formulation, and as a consequence an inconsistency free isochromatic phase map will be obtained. This map can be unwrapped by any simple and fast unwrapping algorithm accurately and effectively. Circular disk and ring under diametric compression samples are both applied for the verification of the proposed algorithm. The experimental results show the proposed algorithm can successfully solve the annoying problems occurred at photoelastic isotropic points with a processing time of roughly 2 seconds for a 420 x 420 pixels map by a general personal computer.
©2010 Optical Society of America
In the field of stress analysis, photoelasticity plays an important role because it is a non-contact whole field optical method and provides isoclinic (principal stress direction) and isochromatic (principal stress difference) data as well, which serve as the two most significant parameters in the field. Applications are allowed in either 2D or 3D forms. Through epoxy resin modeling structure under test  the three dimensional stress distributions in the interior of a component can be studied. Or by photoelastic sheet form plastic adhesive component under test  the surface stresses in a component can also be investigated. The recent developments of powerful computer techniques and fast growing digital cameras demand have greatly enhanced the usefulness of automated fringe analysis skills and promoted the enhancement of related automated technologies [3–8]. But, unfortunately, there is a coupling problem existing between the isoclinic and isochromatic parameters [9–17] and this phenomenon results in phase ambiguity problems of the isochromatic phase map unless the isoclinic data have been correctly procured before being substituted into the isochromatic formulation for its calculation. Among them, the temporal phase unwrapping method is an effective way of circumventing isoclinic-isochromatic coupling problem either with load-stepping [14,16] or wavelength-stepping  techniques. However, the load-stepping technique is not applicable for stress frozen samples and in addition the birefringence error of optical wave plates for different wavelength also has to be calibrated to minimize its error. Ramesh  has proposed a method which uses a dark field image for ambiguity zone identification and its correction based on a priori knowledge of the stress distribution but this method can be difficult to achieve for a complex distribution of stress. Prasad et al.  have proposed an algorithm for the identification of the boundary pixels of the ambiguous zone and then applies the reversal of intensity levels to these boundary pixels. With this pixel wise detection and correction technique, phase ambiguities can be removed but it can be failed at points with big practical error. Wang and Patterson  proposed the use of fuzzy set classification to classify periods of isoclinic and relative retardation data. Using this algorithm, the user has to define the fringe order at two arbitrary points that allows the algorithm to identify which map corresponds to (σ1 –σ2). Siegmann et al.  proposed an algorithm based on processing area rather than individual pixels to demodulate the data and employing quality assessments of the phase derivative data to guide the unwrapping procedure. In general, quality map guided unwrapping can be greatly affected by the phase derivative whose data could be influenced easily by experimental noise. Pinit and Umezaki  proposed a pixel wise unwrapping algorithm with the retrieval work done from the detected largest valid region then expanding point by point outward till the whole sample region are fully restored. The identification of isotropic point is very critical and crucial in the refinement of the largest valid region. Any error identification can lead to error of the result.
The motivation of this work was to produce an algorithm to resolve the isoclinic ambiguity directly. Once the isoclinics have been correctly retrieved the following isochromatic unwrapping will be an ambiguity free retrieval work and can be effectively done by any 2D phase unwrapping algorithm .
In this paper, a novel isoclinic phase unwrapping algorithm  is utilized for the isoclinic and isochromatic parameters calculation. The newly developed method can successfully focus on the ambiguity problems of isoclinic map. Through phase inconsistencies identifying, branch cuts building, and regional phase unwrapping the proposed algorithm can easily restore photoelastic map with isotropic points  like the circular ring case under diametric compression model.
2. Description of the algorithm
The aim of the algorithm is generating a correct map of isoclinic distribution and then producing a continuous map of isochromatic fringe order that is independent of the isoclinic angle. Two well known phase shifting algorithms [9,13] are considered for the calculation of the isoclinic parameter – the Wang-Patterson algorithm and the Pinit-Umezaki algorithm which are formulated by circular and plane polariscope respectively and with their isoclinic calculations listed as Eq. (1) and (2):Eq. (1) and (2), respectively) simultaneously existing in the numerator and denominator of the above equations and that could cause undefined region when their values approach zero. Pinit and Umezeki  suggested a multi-wavelength compensation method which will be adopted in the following isoclinic calculation to circumvent the undefined problems mentioned above.
Figure 1 shows the schematic diagram of the plane polariscope used in this work for the phase stepping isoclinic calculation. Dark field setup is used and its four phase shifting frames are set as β angle equals 0, π/8, π/4, and 3π/8, respectively and the four phase shifting light intensities are formulated as Eqs. (3)a) ~(3d):
I1,i ~I4,i are the four phase stepping frames and the second sub index i stands for R (red), G (green) and B (blue) images of the color digital camera under white light illumination. Ib and Ia are respectively the background light intensity and light amplitude of the plane polariscope. φ and Δ are of the same physical meanings as they have appeared in Eqs. (1) and (2). The modified isoclinic formulation is as follows:
N<> represents a normalization calculation and will convert all image data into an interval ranging between 0 and 1. Since isoclinic calculation should be independent of light wavelength it used, thus if there exists any undefined term resulted from null sine square term in either one of the wavelengths of RGB, the summation operation of which in Eq. (4) will automatically average out the undefined result. In addition, the sin2Δ/2 term multiplier won’t change the signs of both of the numerator and denominator, therefore, the arctangent operation can be expanded into a 2π full module and the calculated wrapped isoclinic is thus spanning over a π/2 interval.
A 1D simulated isoclinic distribution (wrapped and unwrapped) is provided for illustration (see Fig. 2 .) Basically, the simulated true isoclinic data should be continuously distributed except passing through isotropic point which causes abrupt phase discontinuity thereat. In Fig. 2, the red solid line is to simulate a linear phase distribution increasing gradually from left to right with an abrupt phase jumping at the isotropic point. Its wrapped data should be spanned over a range of π/2 (-π/4 ~ + π/4 in this work) as depicted by the blue dash line. While checking these two lines (i.e., the wrapped and unwrapped data), a preliminary operation should be done first to shift all the negative wrapped data (i.e., data ranging between -π/4 ~0) by a translation of π/2 into the interval of ( + π/4 ~ + π/2). This operation can be formulated as Eq. (8) and will convert the red data into the green one (see Fig. 2). Further checking between the green data and the red one, not much difficulty is found that to convert the green wrapped data into the final isoclinic can be done simply by shifting each successively neighboring areas, which are set apart by the phase discontinuities, by π/2 phase or not alternatively. To implement the work, the regional phase unwrapping algorithm  is the best choice because of its great efficiency in restoring ensemble data as a whole.
According to the discussions above, the phase unwrapping work is divided into two parts – the retrievals of isoclinic and isochromatic data, respectively. The proposed algorithm is described as follows.
2.1 Unwrapping of isoclinic data
A 2D simulated circular ring under diametrical compression is used for demonstration. Following are the proposed isoclinic unwrapping steps:(a) Use Eq. (4) to achieve the isoclinic wrapped map (see Fig. 3(a) ).(b) Add π/2 to all negative isoclinic data to convert calculated wrapped data of the previous step into a range of 0 ~π/2 (see Fig. 3(b)).(c) Find all the phase inconsistencies points (see appendix) of step (b).(d) Build branch cuts on the phase inconsistencies to avoid any non-conservative phase integration.(e) Detect all the π/2 phase discontinuities boundaries (with branch cut repair) which divide the wrapped map into different regions.(f) Designate two labeled numbers (e.g., 1 and 2) alternatively to the neighboring regions of the previous step. For example, set the biggest area with label number 1, then its neighboring areas with label number 2. Pick out the largest one from the neighboring regions with latest labeled number 2 and find its neighboring areas then set label number 1 to these unlabeled neighboring areas. Repeat the rule on and on till all the regions are labeled (see Fig. 3 (c).)(g) Subtract π/2 to all the isoclinic data on regions with label number 1 to yield the first associated unwrapped result. Similarly, subtract π/2 to all the isoclinic data on regions with label number 2 to yield the other associated unwrapped result.(h) Find the area in Fig. 3(a) of step (a) with isoclinic data ranging in the region of -π/8 to + π/8.(i) Check each individual of step (g) respectively with the result of step (h) to determine the total number of pixel on which both data of (g) and (h) are correspondingly the same.(j) The two results of step (g) are the two in-plane principal stress directional maps of the tested model with the one in higher matching number of pixel as the first principal stress direction (see Fig. 3 (d)).
2.2 Unwrapping of isochromatic data
Provided the isoclinics have been correctly retrieved the following isochromatic calculation will be very straight forward. The Wang-Patterson algorithm [9,12] is used for the isochromatic calculation. A big difference between the present study and the previous works is that the isoclinic substitution of this work is the restored data instead of the wrapped ones. By doing so, no extra phase ambiguity will be brought into the isochromatic calculation. Consequently, the following isochromatic unwrapping jobs can be successfully done by any 1D or 2D efficient algorithm .
Two experiments are conducted to prove the effectiveness of the proposed algorithm. The first case - a stress frozen sample of circular ring under diametric compression is tested. We use the schematic diagram of Fig. 1 with white light as its light source and take pictures by a colored digital camera modeled Fujifilm FinePix 6900z. The three color components (RGB) of the grabbed images are substituted into Eq. (4) to minimize any possible undefined effect due to any single wavelength then the isoclinic phase map can be obtained and shown as Fig. 4(a) . From the experimental result one can see that comparing to the simulation result of Fig. 3(a), certain degree of noise do exist even after color filter suppression. Then, π/2 is added to all negative isoclinic data (ranging in the range of -π/4 ~0) and yields the result of Fig. 4(b). After inconsistencies detection, branch cuts pairing and regions partition, two alternatively labeling regions are formed and depicted as Fig. 4(c). In addition, the -π/8 ~ + π/8 regions of Fig. 4(a) are found and shown as Fig. 4(d). The two labeled regions of Fig. (c) are subtracted respectively by π/2 to yield the result of Fig. 4(e) and 4(f). The data of these two results are respectively checked with the data of Fig. 4(d) to figure out all the pixels whose data are the same in the comparison. An easy operation is to subtract the data of Fig. 4(e) and Fig. 4(f) respectively with that of Fig. 4(d) and to check whether the subtraction result is zero or not. (Near) zero result means the comparison data are the same. All the pixels whose data are the same in the comparison are counted. For this experimental case, 80% and 20% pixels of Fig. 4(d) whose data are the same as the corresponding pixels’ data of Fig. 4(e) and Fig. 4(f), respectively. Therefore, we can conclude that Fig. 4(e) and 4(f) represent the first and second principal stress direction of the tested sample, respectively. The isochromatic wrapped map Fig. 4(g) can be obtained under substituting the data of Fig. 4(e) into the isochromatic formulation of Wang-Patterson algorithm [9,12]. This isochromatic map is free from isoclinic affection and can be unwrapped easily by most of the 2D unwrapping algorithm.
A second example – a stress frozen sample of circular disk under diametric compression is also conducted here. The wrapped isoclinic data is shown in Fig. 5(a) and the corresponding 0 ~π/2 distribution is shown as Fig. 5(b). Figure 5(c) shows the alternative region labeling and its first principal direction is shown as Fig. 5(d). The two tested samples are very popular in the photoelastic stress analysis field and the experimental results show that the proposed algorithm is effective in restoring them. The two practical results are also compared with the simulated ones (see Figs. 6 and 7 ) with the red and green representing the experimental one and simulated one, respectively. Two horizontal lines (the upper and lower ones are in 0.9 and 0.75 diameter height from the button, respectively) are drawn in each sample for verification. The root mean square value of difference between the experiment and simulation for the ring case is 0.037 and 0.046 rad, respectively. For the disk case, their values are 0.024 and 0.057 rad, respectively.
4. Discussion and conclusion
A simple and effective isoclinic phase inconsistencies removal technique is proposed and has been proven to be effective for the phase unwrapping of digital photoelasticity. A simulated case study of a circular ring under diametric compression is detailed to illustrate each consecutive steps of the proposed algorithm. The demonstrated case in this article has isotropic points as well as singular points on it, which makes the isoclinic unwrapping process difficult and complicated. However, the isoclinic unwrapping problems around these points are successfully solved by the proposed algorithm in present study.
Two experimental cases, a disk and a ring under diametric compression, are conducted also. Without having to find the isotropic and singular points on the tested sample, as well as the largest valid region as the initial for the following isoclinic unwrapping work, the experimental results prove the effectiveness and robustness of the proposed algorithm. Thus, the annoying isoclinic-isochromatic interaction problem can be successfully solved. In addition, the proposed algorithm is very simple and can be implemented easily. Its regional retrieving characteristic also makes the unwrapping work very time effective. For the two practical cases in the work it takes about two seconds for each isoclinic map retrieval with image size roughly 420 x 420 pixels.
Phase inconsistency detection is applied to the isoclinic wrapped map to find if there exists any phase inconsistency on the phase map. The phase inconsistency detection rule  is formulated as follows:
The “int” operator is an operator that rounds the value of the bracket to its nearest integer. A 2 x 2 inconsistency check is performed for every point (m, n) on the phase map. Provided the result of N (m, n) equals zero, the point (m, n) is an inconsistency-free point. A positive or negative inconsistency is existed depending on the result of N (m, n) is +1 or −1.
The authors would like to thank the National Science Council of the Republic of China for financially supporting this research under the contract Numbers NSC 97-2221-E-005-010 and NSC98-2221-E-005-009.
References and links
1. S. K. Mangal and K. Ramesh, “Determination of characteristic parameters in integrated photoelasticity by phase-shifting technique,” Opt. Lasers Eng. 31(4), 263–278 (1999). [CrossRef]
2. D. R. Petersen, R. E. Link, M. N. Pacey, S. J. Haake, and E. A. Patterson, “A novel instrumentation for automated principal strain separation in reflection photoelasticity,” J. Test. Eval. 28(4), 229–235 (2000). [CrossRef]
3. G. Petrucci, “Full-field automatic evaluation of an isoclinic parameter in white light,” Exp. Mech. 37(4), 420–426 (1997). [CrossRef]
4. M. J. Ekman and A. D. Nurse, “Completely automated determination of two-dimensional photoelastic parameters using load stepping,” Opt. Eng. 37(6), 1845–1851 (1998). [CrossRef]
5. T. Y. Chen and C. H. Lin, “Whole-field digital measurement of principal stress directions in photoelasticity,” Opt. Lasers Eng. 30(6), 527–537 (1998). [CrossRef]
6. A. D. Nurse, “Automated photoelasticity: weighted least-squares determination of field stresses,” Opt. Lasers Eng. 31(5), 353–370 (1999). [CrossRef]
7. T. Liu, A. Asundi, and C. G. Boay, “Full field automated photoelasticity using two-load-step method,” Opt. Eng. 40(8), 1629–1635 (2001). [CrossRef]
8. D. Zhang, Y. Han, B. Zhang, and D. Arola, “Automated determination of parameters in photoelasticity,” Opt. Lasers Eng. 45(8), 860–867 (2007). [CrossRef]
9. Z. F. Wang and E. A. Patterson, “Use of phase-stepping with demodulation and fuzzy sets for birefringence measurement,” Opt. Lasers Eng. 22(2), 91–104 (1995). [CrossRef]
10. V. Prasad, K. Madhu, and K. Ramesh, “Towards effective phase unwrapping in digital photoelasticity,” Opt. Lasers Eng. 42(4), 421–436 (2004). [CrossRef]
11. P. Siegmann, D. Backman, and E. A. Patterson, “A robust approach to demodulating and unwrapping phase-stepped photoelastic data,” Exp. Mech. 45(3), 278–289 (2005). [CrossRef]
12. K. Ashokan and K. Ramesh, “A novel approach for ambiguity removal in isochromatic phasemap in digital phtoelasticity,” Meas. Sci. Technol. 17(11), 2891–2896 (2006). [CrossRef]
13. P. Pinit and E. Umezaki, “Digitally whole-field analysis of isoclinic parameter in photoelasticity for four-step color phase-shifting technique,” Opt. Lasers Eng. 45(7), 795–807 (2007). [CrossRef]
14. A. Baldi, F. Bertolino, and F. Ginesu, “A temporal phase unwrapping algorithm for photoelastic stress analysis,” Opt. Lasers Eng. 45(5), 612–617 (2007). [CrossRef]
15. J. Villa, J. A. Quiroga, and E. Pascual, “Determination of isoclinics in phtoelasticity with a fast regularized estimator,” Opt. Lasers Eng. 46(3), 236–242 (2008). [CrossRef]
16. T. Y. Chen, “A simple method for the determination of the photoelastic fringe order,” Exp. Mech. 40(3), 256–260 (2000). [CrossRef]
17. T. Y. Chen, “Digital determination of phtoelastic birefringence using two wavelengths,” Exp. Mech. 37(3), 232–236 (1997). [CrossRef]
18. D. C. Ghiglia, and M. D. Pritt, Two-dimensional phase unwrapping: theory, algorithms, and software (John Wiley & Sons, Inc., New York, 1998).
19. M. J. Huang, and B. C. Song, “Ambiguity elimination of photoelastic phase map,” presented at the 2009 SEM Annual Conference and Exposition of Experimental and Applied Mechanics, Albuquerque, New Mexico, USA, 1–4 Jun. 2009.
20. J. J. Gierloff, “Phase unwrapping by regions,” Proc. SPIE 818, 2–9 (1987).
21. D. C. Ghiglia, G. A. Mastin, and L. A. Romero, “Cellur-automata method for phase unwrapping,” J. Opt. Soc. Am. A 4(1), 267–280 (1987). [CrossRef]