Digital Processing Techniques for Fringe Analysis

Abstract

Digital image processing techniques are needed in order to recover the object information encoded in fringe patterns generated in a determined interferometric setup. Main fringe analysis techniques are reviewed in order to give the reader the most fundamental insights for the interpretation of interferograms. Phase shifted, open fringe, lateral shear and other types of interferograms make use of specific procedures to correctly retrieving the searched phase. Here, algorithms are described and tested in numerical simulations and real data.

1. Introduction

Interferometric techniques are widely used for measuring a wide range of physical variables including refraction index, deformations, temperature gradients, etc. Typically, an interferometer is used to generate one or several interferometric fringe patterns that contain the information of the physical variable that is being measured. Those images must be interpreted in order to recover the parameters that are encoded in the fringe patterns generated by the interferometric setup. Thus, fringe analysis methods deal with the problem of a three-dimensional reconstruction (the object information) from a two-dimensional intensity patterns (interferograms) acquired by a CCD camera. Digital interferometry became extensively used since the development of lasers and CCD devices. In those years, however, the resulting interferograms had to be interpreted visually, and only qualitative results were often achieved. Visual interpretation of an interferogram with only straight or circular fringes is not difficult, but things become more complicated for a fringe pattern that combines several regions with circular, straight and crossed fringes of varying density. Rapidly, it was recognized the need of automatic methods for fringe analysis. The first great advance arises with the development of the phase-shifting techniques. With those procedures, a set of interferograms is acquired with a phase shift among them. The phase shifts are usually introduced by a piezoelectric transducer moving the reference mirror in such way that the phase difference between two consecutive interferograms is a constant term. With phase-shifting techniques, it is possible to isolate the sine and cosine of the phase allowing the calculation of the wrapped phase distribution and consequently the continuous phase with an unwrapping algorithm. Another great success came with the method proposed by Takeda (also referred as the Fourier method) performing a band-pass filtering in the Fourier domain. The method of Takeda works only with interferograms that contain open fringes (patterns that consist in nearly straight fringes). In order to generate such interferograms, the reference beam (e.g., in a two arm interferometer) is tilted introducing a large carrier function to the phase. The Fourier transform of these interferograms is composed of three lobules, one at the center that corresponds to the background term and two lobules located symmetrically respect to the origin. One of this lobules and the one that is located at the origin are filtered out. The remaining spectrum is transformed back to the spatial domain from which the so-called wrapped phase can be calculated. A final step is to apply a phase unwrapping technique to recover the continuous phase. Interferometric measurements and fringe analysis techniques are a growing and fast-changing field of research. Through this chapter, we will review the most known procedures.

2. Interferogram acquisition

The wave nature of light can be studied theoretically by a homogeneous partial differential equation of second order, which satisfies the superposition principle:

$\nabla U\left(x,y,z,t\right)=\frac{1}{c^2}\frac{\partial U\left(x,y,z,t\right)}{\partial t}.$ E1

If two waves of the same frequency are superimposed on a point in space, they excite oscillations in the same direction:

$U_1(t)=A_1\text{cos}(\omega t+{\phi }_1)\text{ ,}$ E2

and

$U_2(t)=A_2\text{cos}(\omega t+{\phi }_2).$ E3

In the preceding equations and the subsequent ones in this section, we will drop the spatial dependence for displaying purposes. The amplitude of the resulting oscillation at that point is determined by the equation:

$A^2=A_1{}^2+A_2{}^2+2A_1{A}_2\text{cos}(\delta ),$ E4

where

$\text{δ}={\text{φ}}_1-{\text{φ}}_2.$ E5

If the phase difference, δ, of the oscillations excited by waves remains constant with time, these waves are coherent. In the case of noncoherent waves, δ varies continuosly, taking any values with equal probability, so the average value of cos(δ) is zero. Therefore,

$A^2=A_1^2+A_2^2.$ E6

So, we can conclude that the intensity observed in the superposed point by noncoherent waves equals the sum of the intensities, which create each separately. However, if the difference δ = φ ₁ − φ ₂ is constant, the cos(δ) will also have a constant value over time, but own for each point in space, so that:

$I=I_1+I_2+2\sqrt{I_1{I}_2}\text{cos}(\delta )$ E7

where I ₁ = A ₁ ² and I ₂ = A ₂ ². This superposition at a point in space results in a different sum of the intensities of the separate components intensity [1]. This phenomenon is known as interference. The light and dark areas that observed on screen placed in the region of interference are called interference fringes, and the fringes are intensity which changes from minimum to maximum, which together form a pattern commonly called interference pattern or interferogram, as can be seen in Figure 1 .

Figure 1.

Interference pattern or interferogram showing circular fringes.

The interference of two or more electromagnetic waves can be usually achieved in two ways: by division of the wave front and by division of the amplitude. A mechanism used for the division of the wave front is, for example, the Young’s experiment or double slit that is showed in Figure 2 .

Figure 2.

Young’s experiment, where a double slit produces two wave fronts that interfere on a screen.

A mechanism used for dividing the amplitude of the wave is, for example, the Michelson interferometer that is shown in Figure 3 .

Figure 3.

Michelson interferometer.

Both mechanisms produce a pattern of light and dark intensities in the plane of interference fringes. The image resulting from the interference is known as interferogram. When the interferogram is captured by a recording medium, that is, a photographic film or a CCD camera, the process commonly involves some optical system, which introduces imperfections with respect to the ideal image. Such imperfections are known as aberrations. Aberrations can be classified as chromatic and monochromatic. Chromatic aberrations are present to illuminate the object with white light or polychromatic light, that is, light with different wavelengths. These aberrations are the only ones that can be predicted by the theory of the first order, which states that an optical system consisting of lenses has different focal lengths for different wavelengths. These variations are related to the change of refractive index with respect to wavelength causing that both the position and the image size are different for each wavelength. Monochromatic aberrations occur when the object is illuminated with monochromatic light (i.e., light of a single wavelength), and the reflected or transmitted light is registered by a recording medium. This type of aberration causes that the captured image of a punctual object is no longer a point, but a blurred point. Monochromatic aberrations can be calculated roughly in the third-order theory, using the first two terms of the expansion in power series of the sinθ and cosθ. Another alternative to calculate more accurately is to make the exact trigonometric trace rays through the system, where the deviation of the rays is calculated. These aberrations were studied in detail first by Ludwig von Seidel, hence often called Seidel aberrations, and are sphericity, coma, astigmatism, field curvature and distortion. These aberrations have a distinctive fringe distribution that appears frequently in the testing of a variety of objects with an interferometer. As the interferogram acquisition is only the capture of the interference phenomenon or interference pattern in a recording medium, it is convenient to study correcting aberrations problems. From these patterns, it can be known as the aberration coefficients using Zernike polynomials. Zernike polynomials have been successfully used in the recognition of patterns and image processing. Additionally, these have been used in astronomy to describe wavefront aberrations due to atmospheric turbulence and to describe wavefront aberrations in the human eye. This is because of the Seidel aberrations are related to the Zernike polynomials. Due to this relationship, polynomials are used to describe wavefront aberrations in order to calculate the aberration coefficients of a set of interferograms generated by an adaptive lens. These coefficients enable to describe the behavior of the aberrations present in the lens.

3. Interferogram filtering

In order to enhance the interferogram images and reduce the noisy caused by external factors at the interference phenomenon in an interferometer (interference produced by two or more controlled wave fronts), low-pass techniques are used as preprocessing filtering step. Interferogram smoothing and denoising are the principal purposes for the application of low-pass filters in interferometry. Filtering techniques can be grouped in spatial or frequency domains, where the spatial filtering is directly applied to the interferogram image, pixel to pixel, while the frequency filtering is usually performed in the Fourier domain. Band-pass and band-stop are some filtering techniques that can be applied in the Fourier domain, and these filters are used to attenuate some frequencies of specific noise.

3.1. Spatial filtering

The mathematical entity applied in spatial filtering is the convolution operation, also known as windowing [2], written as follows:

$f\left(x,y\right)\text{*}g\left(x,y\right)={\displaystyle \sum }_{m=-\frac{M+1}{2}}^{\frac{M+1}{2}}{\displaystyle \sum }_{n=-\frac{N+1}{2}}^{\frac{N+1}{2}}f\left(x-m,y-n\right)g(m,n),$ E8

where f, g, (x, y), (m, n) and M×N are the image to be filtered, the convolution mask, the original image coordinates, the coordinates where the convolution is performed and the size of the convolution mask, respectively. The kind of convolution filter is determined by the chosen convolution filter (averaging, Gaussian, quadratic, triangular, trigonometric, etc.). These convolution functions are presented in a mathematical form as follows: [3]

$g\left(m,n\right)=\frac{1}{{{\displaystyle \sum }}_{m=1}^M{{\displaystyle \sum }}_{n=1}^N{\omega }_{m,n}}\left[\begin{array}{ccc}{\omega }_{1,1}& \cdots & {\omega }_{1,N}\\ ⋮& \ddots & ⋮\\ {\omega }_{M,1}& \cdots & {\omega }_{M,N}\end{array}\right],$ E9

where

${\omega }_{m,n}=\{\begin{array}{cc}\begin{array}{c}\begin{array}{c}A{e}^{\left(-B\left[m^2+n^2\right]\right)},\\ A-B\left(m^2+n^2\right),\end{array}\\ \begin{array}{c}\frac{A}{2}+\frac{A}{4}\left[\cos \left(Bm\right)+\text{cos}(Bn)\right],\\ A,\end{array}\end{array}& \begin{array}{c}\begin{array}{c}\text{Gaussian}\\ \text{quadratic}\end{array}\\ \begin{array}{c}\text{trigonometric}\\ \text{average}\end{array}\end{array}\end{array}.$ E10

Here A, B and ω are the amplitude, the width factor and the weight of the spatial filter function, respectively. The study about the magnitude spectrums of Eq. (10) (low-pass masks) was reported in Ref. [2], where the Gaussian and quadratic masks delivered the best results. Low frequencies were conserved by these filters, while the high frequencies were attenuated. However, the filtering results may vary depending on the interferogram under process, as can be seen in Figure 4 , where the results of the filter on a ronchigram (a particular kind of interferogram) can be seen in Figure 4a . The four kinds of masks presented in Eq. (10) were employed to generate the filtered fringe patterns shown in Figure 4b . The parameters used by the filters were A = 1 and B = 0.1, with mask sizes of 3×3, 5×5 and 7×7 pixels.

Figure 4.

Spatial low-pass filters: (a) original ronchigram, (b) convolution image filtering with Gaussian, quadratic, trigonometric and average masks, as well as 3 × 3, 5 × 5 and 7 × 7 convolution mask sizes.

3.2. Frequency filtering

Frequency filtering is usually performed in the Fourier domain. The Fourier transform represents the change from spatial to frequency domain. Eq. (11) and Eq. (12) represent a pair of discrete Fourier transforms in two dimensions [2]

$F\left(u,v\right)=\mathcal{F}\left\{f\left(x,y\right)\right\}={\displaystyle \sum }_{x=1}^U{\displaystyle \sum }_{y=1}^{V}f\left(x,y\right)e^{-2\pi i\left(\frac{ux}{U}+\frac{vy}{V}\right)},$ E11

and

$f\left(x,y\right)={\mathcal{F}}^{-1}\left\{F\left(u,v\right)\right\}={\displaystyle \sum }_{u=1}^U{\displaystyle \sum }_{v=1}^{V}F\left(u,v\right)e^{2\pi i\left(\frac{ux}{U}+\frac{vy}{V}\right)}$ E12

where (u, v), U×V, $\mathcal{F}$ and $\mathcal{F}$ ⁻¹ are the frequency coordinates, the image size in pixels, the Fourier transform and the inverse Fourier transform operators, respectively. A significant attribute obtained from the Fourier transform is that it gives the frequencies content of the image. Due to this property, frequency filter design is a very straight forward task for image processing. Low frequencies are located into the Fourier domain around the central coordinates; as frequencies gradually increase, they spread out from the center in a radial form. This characteristic is ideal for frequency filtering (low-pass, high-pass, band-pass and band-stop) [3]. The frequency filtering development consists in the multiplication between the Fourier transform with some kind convolution function. The kind of convolution mask will determine the class of filtering to be performed. The following is a summary of the most usual filtering masks: a white centered circle for low-pass; a black centered circle for high-pass; a white centered ring for band-pass and a black centered ring for band-stop.

In Figure 5 , it is showed the masks described above along with the results of the filtering process. Low-pass filtering masks and results are presented in Figure 5b , band-pass filters and filtered images are seen in Figure 5c and finally, band-stop filters and filtering results can be appreciated in Figure 5d . The kind of filter or the size of geometrical mask depends directly in the image and in its noise content. There is no ideal filter; the kind of applied filter to process an image is dependent of the characteristics that are pretended to enhance or eliminate.

Figure 5.

Frequency filtering: (a) Fourier transform of the original ronchigram, (b) low-pass filtering with different radius of binary circle mask, (c) band-pass filtering with two size of binary ring mask, and (d) band-stop filtering with two size of binary ring mask.

4. Phase-shifting interferometry

Interference fringe patterns obtained by means of interferometric techniques can be evaluated by using digital image processing techniques for the estimation of phase map distributions. The intensity of an image in an interference fringe pattern, according with Eq. (7), can be represented by Eq. (13)

$I\left(x,y\right)=a\left(x,y\right)+b\left(x,y\right)cos\left[\varphi \left(x,y\right)\right]$ E13

where ϕ(x, y) is the phase difference between the reference and testing wave fronts that interfere, a(x, y) = I ₁(x, y) + I ₂(x, y) is the mean intensity and b(x, y) = 2(I₁ I₂ )^1/2 is the intensity modulation registered by each pixel (x, y) of a CCD camera. When analyzing a single image for the estimation of the relative phase difference in applied metrology, the sign of the phase cannot be assessed because of the argument in the sinusoidal function of the modulation term in Eq. (13) results in an identical interferogram for ±ϕ(x, y) values. In order to overcome this difficulty, multiple interferograms can be sequentially registered introducing a known small phase change amount which is linear in time; this method is well known as phase shifting [4]. In practical phase shifting, a piezoelectric device PZT is included in the interferometric system to produce the phase-shifted interferograms. In general, when a voltage signal is used to polarize a PZT actuator, this electrical signal is converted directly into linear displacement motion. Then, phase-displaced intensity images can be represented by

$I\left(x,y\right)=a\left(x,y\right)+b\left(x,y\right)cos\left[\varphi \left(x,y\right)+\theta \right].$ E14

Since there are three unknown terms in the representation of the interference intensity equation, then the measurement of at least three interferograms at known phase shifts is needed to determine the relative phase difference. However, one of the simplest modes to determinate the phase considers the use of four interferograms equally spaced by θ = π/2, obtaining:

$I_1=a+b\text{cos}(\varphi ),$ E15

$I_2=a-bsin(\varphi ),$ E16

$I_3=a-bcos(\varphi ),$ E17

and

$I_4=a+bsin(\varphi ).$ E18

Rearranging the simultaneous equation system from the above formulas, in which the spatial dependence (x, y) was not included, the phase distribution from intensity images can be calculated with:

$tan\left(\varphi \right)=\frac{I_4-I_2}{I_1-I_3}.$ E19

The phase difference estimated from the four equations is determined in the range between −π and π when using the arctangent function, hence producing a wrapped phase distribution. In order to analyze an interference pattern by the phase-shifting techniques, a Twyman-Green interferometer is considered for the testing of optical elements; the basic arrangement is shown in Figure 6 . A laser diode is used in the interferometer as illumination source, and then the emerging laser beam is collimated by a lens to obtain a plane wave front that propagates throughout a 50:50 beam splitter to produce two wave fronts of same amplitude. One beam is deviated to the reference mirror M1, and the second beam, to the mirror M2 under test. Next, the beams are reflected back toward the beam splitter, and part of the intensity overlaps on the observation plane, where a CCD camera sensor is placed to register an image of the resulting interferogram that is then stored in a PC for subsequent processing.

Figure 6.

Twyman-Green interferometer for the testing of optical elements.

In Figure 7 , a set of four phase-shifted interferograms with phase shifts of π/2 is shown. A first interferogram with θ = 0 is seen in Figure 7a , and subsequent interferograms with θ = π/2, θ = π and θ = 3π/2 are observed in Figure 7b – d , respectively. The wrapped phase is showed in Figure 7e , and the unwrapped phase related with the shape of the optical element being tested is shown in Figure 7f . A three-dimensional representation of the unwrapped phase seen in Figure 7f is presented in Figure 8 .

Figure 7.

Four-step phase-shifted interferograms from (a) to (d), the wrapped phase (e) and the unwrapped phase (f).

Figure 8.

The unwrapped continuous phase obtained from a set of four phase-shifted fringe patterns and an unwrapping method.

In applied phase-shifting interferometry, there are concerns about the presence of errors that may affect the accurate phase extraction from phase-shifted interferograms. A typical systematical source of error introduced PZT arises when there is miscalibration of the phase-shifting actuator, causing detuning in the phase extraction process. The inclusion of phase-shifting algorithms with more than three or four interferograms can be implemented in order to reduce this systematic error [5].

5. The Fourier method

The method of Takeda or the Fourier method was developed in 1982 [6]. Unlike phase-shifting methods, see Ref. [7], a single interferogram with open fringes can be analyzed. This is useful when the object under study changes dynamically or environmental disturbances (vibrations or air turbulence) do not allow the use of phase-shifting methods unless special, and often, very expensive hardware is used to acquire several images simultaneously. However, in general, the accuracy and the dynamical range of the phase that can be measured are reduced. The Fourier method makes use of the fast Fourier transform technique to separate, in the frequency domain, the background and phase terms of the interferogram. Employing complex notation, an interferogram can be written as follows:

$I\left(x,y\right)=a\left(x,y\right)+\frac{1}{2}b\left(x,y\right)e^{i\left[\varphi \left(x,y\right)+\phi (x,y)\right]}+\frac{1}{2}b\left(x,y\right)e^{-i\left[\varphi (x,y)+\phi (x,y)\right]},$ E20

where a(x, y) and b(x, y) are the background intensity and the modulation term, respectively. The phase term is denoted by ϕ(x, y) and finally, the symbol φ(x, y) denotes the lineal carrier function or tilt that is introduced usually by tilting the reference mirror in a two arm interferometer. The Fourier transform of the above expression can be written as:

$\mathcal{F}\left\{I(x,y)\right\}=\tilde{I}\left(u,v\right)=\delta \left(u,v\right)+E\left(u+\alpha ,v+\beta \right)+E^{\text{*}}(u-\alpha ,v-\beta ),$ E21

where (u, v) is the coordinates in the frequency domain, δ(u, v) is a delta function and E(u + α, v + β) and E*(u-α, v-β) are complex conjugate functions that correspond to the transforms of the second and third terms of Eq. (20), respectively. The introduction of the linear carrier function, φ(x, y) = αx + βy, shifts the terms E(u, v) and E*(u, v) in opposite directions in the frequency spectrum as can be seen noting that:

$E\left(u,v\right)=\mathcal{F}\left\{\frac{1}{2}b\left(x,y\right)e^{i\varphi (x,y)}\right\},$ E22

and

$E\left(u+\alpha ,v+\beta \right)=\mathcal{F}\left\{\frac{1}{2}b\left(x,y\right)e^{i\varphi (x,y)}{e}^{i\phi (x,y)}\right\},$ E23

similarly

$E^{\text{*}}\left(u,v\right)=\mathcal{F}\left\{\frac{1}{2}b\left(x,y\right)e^{-i\varphi (x,y)}\right\},$ E24

and

$E^{\text{*}}\left(u-\alpha ,v-\beta \right)=\mathcal{F}\left\{\frac{1}{2}b\left(x,y\right)e^{-i\varphi (x,y)}{e}^{-i\phi (x,y)}\right\}.$ E25

The separation of the terms in the Fourier spectrum due to the introduction of a linear carrier can be observed in Figure 9 . An open-fringe interferogram and its Fourier spectrum can be seen in Figure 9a and b . It can be noted that the three terms of Eq. (21) are clearly separated. The central peak corresponds to the delta function δ(u, v) while the adjacent lobules correspond to E(u + α, v + β) and E*(u - α, v - β). The interferogram was constructed as follows:

Figure 9.

Interferogram with open fringes (a) and its Fourier spectrum (logarithm of the module) (b).

$a\left(x,y\right)=120e^{-(x^2+y^2)}$ E26

$b\left(x,y\right)=124e^{-(x^2+y^2)},$ E27

$\varphi \left(x,y\right)=2\pi \left\{4e^{-4\left[{(x-0.2)}^2+{(y-0.3)}^2\right]}-5e^{-6\left[{(x+0.2)}^2+{(y+0.3)}^2\right]}\right\},$ E28

and

$\phi \left(x,y\right)=2\pi \left(16x+20y\right)$ E29

where x and y vary from 1 to −1 along the vertical and horizontal directions and the interferogram was multiplied by a circular function of radius one.

In order to recover the phase, we need to isolate one of the lateral lobules of the Fourier domain. To this end, we employ a band-pass filter. The filtered spectrum is then transformed back to the spatial domain to obtain I_f (x, y) and the wrapped phase is found with the arctangent function of the ratio of the imaginary and real parts of I_f (x, y) as shown in Figure 10 . Figure 10a shows the band-pass filter H(u, v), the filtered spectrum can be observed in Figure 10b and finally the wrapped phase is shown in Figure 10c .

Figure 10.

Band-pass filtering process of the Fourier spectrum seen in Figure 9 . Pass-band filter (a), filtered spectrum (b) and recovered wrapped phase (c).

The band-pass filter has the following form:

$H\left(u,v\right)=e^{-3000{\left[{(u-0.2)}^2+{(v-0.2)}^2\right]}^2},$ E30

where u and v vary from 1 to −1 in vertical and horizontal directions, respectively. The filtered spectrum becomes,

${\tilde{I}}_f\left(u,v\right)=\tilde{I}\left(u,v\right)H(u,v).$ E31

Transforming back to spatial domain we obtain:

$I_f\left(x,y\right)={\mathcal{F}}^{-1}\left\{{\tilde{I}}_f(u,v)\right\}.$ E32

The wrapped phase is found by:

${\psi }_w\left(x,y\right)=atan2\left\{imag\left[I_f(x,y)\right],real\left[I_f(x,y)\right]\right\}$ E33

where the atan2() function accepts two arguments corresponding to the sine and cosine and returns the results modulo 2π. This wrapped phase, however, is not the desired one because it contains the introduced tilt that is not part of the object information. The wrapped desired phase is found with:

${\varphi }_w\left(x,y\right)=atan2\left\{sin\left[{\psi }_w\left(x,y\right)-\phi (x,y)\right],cos\left[{\psi }_w\left(x,y\right)-\phi (x,y)\right]\right\}\text{.}$ E34

The final step is to apply an unwrapping method to obtain the continuous phase related with the object under study. This last procedures are shown in Figure 11a and b where the wrapped phase φ_w (x, y) and the unwrapped reconstructed phase φ_r (x, y) can be appreciated, respectively.

Figure 11.

Phase reconstruction results. Wrapped phase (a) and unwrapped continuous phase (b).

The Fourier method is not the unique procedure for phase retrieval from one interferogram with open fringes. Besides the Fourier approach there are other procedures in the spatial domain including phase locked loop [8] and spatial carrier phase-shifting methods [9, 10].

6. Phase unwrapping

Phase unwrapping is a common step to finally find a continuous phase for several fringe analysis techniques such as phase shifting, Fourier, the phase synchronous and others methods that use the arctangent function of the ratio of the sine and cosine of the phase to obtain a wrapped phase. In its simplest form, phase unwrapping consists in adding or subtracting 2π terms to the pixel being unwrapped if a difference greater that π is found with a previous pixel already unwrapped [11]. The phase unwrapping problem in one dimension can be observed in Figure 12 . The wrapped phase found with the arctangent function is seen in Figure 12a and the unwrapped continuous phase is showed in Figure 12b .

Figure 12.

Wrapped (a) and unwrapped phase (b).

The described procedure works well only for wrapped phases with no inconsistencies and low noise levels, however, delivers wrong results when dealing with noisy wrapped phases or those obtained from interferograms with broken or unconnected fringes. A more consistent approach is achieved with the least squares phase unwrapping method [12]. The least square technique integrates the discretized laplacian of the phase. To this end, the laplacian of the phase is calculated as follows:

$L_{i,j}={\varphi }_{i+1,j}^x-{\varphi }_{i,j}^x+{\varphi }_{i,j+1}^y-{\varphi }_{i,j}^y,$ E35

where

${\varphi }_{i,j}^x=atan2\left[\sin \left({\varphi }_{i,j}^w-{\varphi }_{i-1,j}^w\right),co\text{s}({\varphi }_{i,j}^w-{\varphi }_{i-1,j}^w)\right]p_{i,j}{p}_{i-1,j},$ E36

and

${\varphi }_{i,j}^y=atan2\left[\sin \left({\varphi }_{i,j}^w-{\varphi }_{i,j-1}^w\right),co\text{s}({\varphi }_{i,j}^w-{\varphi }_{i,j-1}^w)\right]p_{i,j}{p}_{i,j-1}.$ E37

In the last equations, we have used pixel subscript notation in order to limit the extension of the equations. The pupil function p_{i, j} is defined as equal to one where we have valid data and zero otherwise. One may note that if p_{i, j} = p_{i+1, j} = p_{i-1, j} = p_{i, j+1} = p_{i, j-1} = 1, then

$L_{i,j}=-4{\varphi }_{i,j}+{\varphi }_{i+1,j}+{\varphi }_{i-1,j}+{\varphi }_{i,j+1}+{\varphi }_{i,j-1}.$ E38

Solving for the phase in the above equation, we obtain that the unwrapping problem under the least square approach consist in the resolution of a linear system of equations, as:

${\varphi }_{i,j}=\frac{{\varphi }_{i+1,j}+{\varphi }_{i-1,j}+{\varphi }_{i,j+1}+{\varphi }_{i,j-1}-L_{i,j}}{4}.$ E39

An iterative technique that solves the above system of linear equations is the overrelaxation method in which the following equation is iterated until convergence:

${\varphi }_{i,j}^{k+1}={\varphi }_{i,j}^k-\frac{\left(d{\varphi }_{i,j}^k-{\varphi }_{i-1,j}^{k+1}-{\varphi }_{i+1,j}^k-{\varphi }_{i,j-1}^{k+1}-{\varphi }_{i,j+1}^k+L_{i,j}\right)w}{d},$ E40

where

$d=p_{i-1,j}+p_{i+1,j}+p_{i,j-1}+p_{i,j+1}.$ E41

In the last equations, k is the iteration number and w is a parameter of the overrelaxation method that must be set between 1 and 2. Results of the least square technique are showed in Figure 13 . A two-dimensional wrapped phase is seen Figure 13a , the phases differences are shown in Figure 13b and c , respectively. The laplacian of the phase is showed in Figure 13d , the reconstructed phase in a two-dimensional view is observed in Figure 13e and finally, for comparison purposes the reconstructed rewrapped phase is also showed in Figure 13f . One may observe that the original wrapped phase and the reconstructed wrapped phase are slightly different, this is because the least square method recovers the phase with an arbitrary constant term, however this term is usually neglected since doesn’t carry any information of the object being measured.

Figure 13.

Least square results to unwrap the phase. Wrapped phase (a), phase differences in the x (b) and y (c) directions, laplacian of the phase (d), unwrapped phase (e) and rewrapped phase (f).

If desired, the constant term in the retrieved phase may be corrected easily in the following form:

${\varphi }_c={\varphi }_r-atan2\left[\sin \left({\varphi }_r-{\varphi }_w\right),\text{cos}({\varphi }_r-{\varphi }_w)\right].$ E42

The wrapped phase seen in Figure 13a was constructed as follows:

${\varphi }_w=atan2\left[\sin \left(\varphi \right),\text{cos}(\varphi )\right],$ E43

where

${\varphi }_{i,j}=2\pi \{\begin{array}{c}16{\left[{\left(x_i\right)}^2+{\left(y_j\right)}^2\right]}^2-18\left[{\left(x_i\right)}^2+{\left(y_j\right)}^2\right]+\\ 5y_j{\left(x_i\right)}^2-4x_i{\left(y_j\right)}^2+3y_j{x}_i+2x_i-2y_j+0.56\end{array}\}.$ E44

In the above equations x_i and y_j are range variables that vary from −1 to 1, the wrapped phase was multiplied by an annular pupil function with an exterior radius of 198 pixels, while the interior radius was 60 pixels for an image size of 400 × 400. The convergence of the reconstructed phase seen in Figure 2e was reached after 700 iterations; the overrelaxation parameter w was set to 1.99, which is usual for large images. The reconstructed phase with the constant term corrected and the phase error, ε = ϕ = ϕ _c , can be appreciated in Figure 14a and b . The maximum error was about 0.00005 radians.

Figure 14.

Unwrapped phase and phase error. Three-dimensional view of the reconstructed and corrected phase (a) and phase error (b).

Finally, results on a noisy wrapped phase are presented. Random noise with uniform distribution in the range of −π/4 to π/4 radians were added to the phase to obtain the wrapped phase seen in Figure 15a , the unwrapped phase is shown in Figure 15b , and the rewrapped reconstructed phase is observed in Figure 15c . As can be noticed the least square method is a very reliable technique that works with any pupil configuration and stands noisy measurements.

Figure 15.

Results on noisy measurements. Wrapped phase with noise (a), retrieved phase (b) and wrapped retrieved phase (c).

7. Phase recovery from lateral shearing interferograms

Lateral shearing interferometry is a very important field in experimental optical measurements, in which, the test beam interferes with a laterally displaced version of itself instead of a reference beam. The resulting fringe patterns are thus related with the object wavefront derivative in a given direction. This is very useful when the object information of interest is related with the derivative as in strain analysis or when the dynamical range of the object wave front is too high that cannot be measured with direct interferometry. Let us consider a laterally shear interferogram with a beam displacement in the x direction a quantity Δx, we obtain:

$I^x\left(x,y\right)=B\left(x,y\right)+C\left(x,y\right)cos\left[{\psi }^x(x,y)\right],$ E45

where

${\psi }^x\left(x,y\right)=\varphi \left(x,y\right)-\varphi (x-\text{Δ}x,y).$ E46

In Eq. (45) B(x, y) is the background intensity and C(x, y) is the modulation term. The objective is to retrieve the undisplaced phase φ(x, y). To this end, we need, at least, another laterally shear interferogram with a beam displacement in the y direction another quantity Δy, obtaining:

$I^y\left(x,y\right)=B\left(x,y\right)+C\left(x,y\right)cos\left[{\psi }^y(x,y)\right],$ E47

where

${\psi }^y\left(x,y\right)=\varphi \left(x,y\right)-\varphi (x,y-\text{Δ}y).$ E48

Let us considered that we have retrieved the phase differences ψ^x (x, y) and ψ^y (x, y) by means of a phase-shifting technique as is depicted in Figure 16 and Figure 17 for the x and y directions, respectively. A set of four phase-shifted interferograms can be seen in Figure 16a . The wrapped phase differences ψ^wx (x, y) and the unwrapped phase ψ^x (x, y) are observed in Figure 16e and f , respectively. The modulation term and a sheared pupil function p^x (x, y) are shown in Figure 16g and h , respectively. The sheared pupil p^x (x, y) is found after normalization from zero to one and thresholding the normalized modulation term. A second set of phase-shifted interferograms are seen Figure 17a and d . The wrapped phase differences ψ^wy (x, y) and the unwrapped phase ψ^y (x, y) are observed in Figure 17e and f , respectively. The modulation term and a sheared pupil function p^y (x, y) are shown in Figure 4g and h , respectively. In a similar way to the first set of phase-shifted interferograms, the pupil p^y (x, y) is found after normalization from zero to one and thresholding the normalized modulation term. The sheared pupils p^x (x, y) and p^y (x, y) are useful to find the undisplaced pupil p(x, y), where the original wave front φ(x, y) have valid data, since p^x (x, y) = p(x, y)p(x − Δx, y) and p^y (x, y) = p(x, y)p((x, y) − Δy)

Figure 16.

Recovery of the phase differences in the x direction. Set of four sheared interferograms acquired under a phase-shifting technique (a) to (d), wrapped phase differences (e), unwrapped phase differences (f), modulation (g) and recovered sheared pupil (h).

Figure 17.

Recovery of the phase differences in the y direction. Set of four sheared interferograms acquired under a phase-shifting technique (a) to (d), wrapped phase differences (e), unwrapped phase differences (f), modulation (g) and recovered sheared pupil (h).

Once the phase differences θ^x (x, y) and θ^y (x, y) are known, we can use the next procedure to find the searched phase φ(x, y) observing that:

${\psi }_{i,j}^x{p}_{i-a,j}-{\psi }_{i+a,j}^x{p}_{i+a,j}={\varphi }_{i,j}\left(p_{i-a,j}-p_{i+a,j}\right)-{\varphi }_{i-a,j}{p}_{i-a,j}-{\varphi }_{i+a,j}{p}_{i+a,j},$ E49

and

${\psi }_{i,j}^y{p}_{i,j-b}-{\psi }_{i,j+b}^y{p}_{i,j+b}={\varphi }_{i,j}\left(p_{i,j-b}-p_{i,j+b}\right)-{\varphi }_{i,j-b}{p}_{i,j-b}-{\varphi }_{i,j+b}{p}_{i,j+b}.$ E50

In the above equations we have changed the (x, y) dependence by pixel subscript notation. As described before, p_{i, j} is the undisplaced pupil and the displacement quantities Δx and Δy are rounded to the nearest integer in pixel dimensions obtaining a and b, respectively. Adding Eq. (49) and Eq. (50) and solving for the phase, we obtain:

${\varphi }_{i,j}=\frac{F+G}{H}\text{ ,}$ E51

where

$H=p_{i+a,j}+p_{i-a,j}+p_{i,j+b}+p_{i,j-b},$ E52

$F={\psi }_{i,j}^x{p}_{i-a,j}-{\psi }_{i+a,j}^x{p}_{i+a,j}+{\psi }_{i,j}^y{p}_{i,j-b}-{\psi }_{i,j+b}^y{p}_{i,j+b},$ E53

and

$G={\varphi }_{i-a,j}{p}_{i-a,j}+{\varphi }_{i+a,j}{p}_{i+a,j}+{\varphi }_{i,j-b}{p}_{i,j-b}+{\varphi }_{i,j+b}{p}_{i,j+b}.$ E54

Eq. (51) represents a linear system of equations; however, it is ill posed because there are more unknowns than equations due to the effects of the sheared pupils. Nevertheless, a regularization term may be aggregated to overcome this problem [13–15]. The regularization term is in the form of discrete Laplacians of the phase among adjacent pixels. The following equation that incorporates the regularization term is iterated until convergence:

${\varphi }_{i,j}^{k+1}={\varphi }_{i,j}^k-\frac{H{\varphi }_{i,j}^k-\left(F+G\right)+\alpha (L_{i+1,j}^x-2L_{i,j}^x+L_{i-1,j}^x+L_{i,j+1}^y-2L_{i,j}^y+L_{i,j-1}^y)}{H},$ E55

where

$L_{i,j}^x=\left({\varphi }_{i+1,j}-2{\varphi }_{i,j}+{\varphi }_{i-1,j}\right)p_{i+1,j}{p}_{i-1,j},$ E56

and

$L_{i,j}^y=\left({\varphi }_{i,j+1}-2{\varphi }_{i,j}+{\varphi }_{i,j-1}\right)p_{i,j+1}{p}_{i,j-1}.$ E57

Here, α is a parameter that controls the effects of the regularization term. The phase reconstruction is seen in Figure 18 . A two-dimensional view of the retrieved phase is observed in Figure 18a , the same phase but in a three-dimensional perspective is shown in Figure 18b and the phase error (the actual phase minus the reconstructed one) can be appreciated in Figure 18c . This reconstruction was achieved after 800 iterations with a parameter α = 0.1 obtaining a maximum error of about 0.0004 radians.

Figure 18.

Phase reconstruction from lateral sheared interferograms. Reconstructed phase (a), three-dimensional view of the reconstructed phase (b) and the phase error (c).

The actual phase was constructed as follows:

${\varphi }_{i,j}=2\pi \left\{\begin{array}{c}20{\left[{\left(x_i\right)}^2+{\left(y_i\right)}^2\right]}^2-24\left[{\left(x_i\right)}^2+{\left(y_i\right)}^2\right]-6{\left(x_i\right)}^2{y}_j\\ +8x_i{\left(y_j\right)}^2+5x_i{y}_j+4x_i-5y_j+0.45\end{array}\right\}\text{ ,}$ E58

where x_i and y_j are range variables that vary from −1 to 1 in both vertical and horizontal directions. The pupil function is a circular one with a radius of 170 pixels. The shear distances were set to Δx = a = 12 pixels and Δy = b = 12 pixels.

8. Conclusions

Digital processing techniques applied to interferometric measurements allow to obtain the phase from fringe patterns. The fringe analysis methods described here can be used to recover the phase that is associated with the physical variable under study. Under controlled conditions, phase-shifting techniques are the most used methods to retrieve the wrapped phase. If experimental conditions suffer from vibrations, air turbulences or the object changes dynamically, among other factors, then a Fourier method may be preferable to analyze an open-fringe interferogram. Those procedures deliver a wrapped phase. Then, an unwrapping algorithm is needed to reconstruct a continuous phase related with the object being studied. The aim of this chapter is to present, to the reader, the fundamentals of principal fringe analysis techniques. Numerical simulations are provided, in such way that the reader can reproduce them by its own. The extension of the chapter is insufficient to introduce many important techniques. However, the methods presented here were described as clearly and briefly as we could. We hope that the reader finds this information useful in the interpretation of interferograms obtained in the study of some object or phenomena by using an interferometric setup.