Ao 57 8 1944
Ao 57 8 1944
Ao 57 8 1944
Received 27 October 2017; revised 6 February 2018; accepted 9 February 2018; posted 9 February 2018 (Doc. ID 310218);
published 9 March 2018
This paper presents a method for extracting multiple phases from a single moire fringe pattern in digital holo-
graphic interferometry. The method relies on component separation using singular value decomposition and
an extended Kalman filter for demodulating the moire fringes. The Kalman filter is applied by modeling the
interference field locally as a multi-component polynomial phase signal and extracting the associated multiple
polynomial coefficients using the state space approach. In addition to phase, the corresponding multiple phase
derivatives can be simultaneously extracted using the proposed method. The applicability of the proposed method
is demonstrated using simulation and experimental results. © 2018 Optical Society of America
OCIS codes: (120.3180) Interferometry; (100.2650) Fringe analysis; (090.2880) Holographic interferometry; (100.5070) Phase
retrieval.
https://doi.org/10.1364/AO.57.001944
filtering, or complex two-dimensional phase unwrapping oper- a matrix enhancement technique [18,22] where the given block
ations. In addition to multiple phases, the mathematical formu- B is converted into a larger block by partitioning and stacking
lation in the proposed method also allows to simultaneously its rows, which is shown in Eq. (3). The enhanced matrix, Be
extract the multiple phase derivatives. The theory of the pro- with size equal to JK × M b − J 1N b − K 1, encloses
posed method is explained in the next section. Simulation and the maximum amount of noise that can be filtered out using
experimental results are illustrated in Section 3. Subsequently, singular value decomposition (SVD).
discussion about the proposed method and conclusions are The enhanced matrix is given by
outlined. 2 3
B0 B1 BM b −J
6 B 7
6 1 B2 BM b −J1 7
2. THEORY 6 7
Be 6 . .. . . .. 7; (3)
6 .. 7
The experimental setup for recording multi-dimensional dis- 4 . . . 5
placements using DHI is shown in Fig. 1. In the dual-wave BJ−1 BJ BM b −1
setup, the object is placed along the x–y plane, and the holo-
grams are recorded at the image plane using a camera placed where
normal to z axis. The object is illuminated symmetrically by 2 3
two object beams, and their scattered beams are made to inter- bn; 0 bn; 1 bn; N b − K
6 7
fere with a reference beam at the image plane. The reference 6 bn; 1 bn; 2 bn; N b − K 1 7
6 7
beam and the object beams are derived from the same laser Bn 6 .. .. .. .. 7:
6 7
source. The two-dimensional numerically reconstructed com- 4 . . . . 5
plex field can be represented as [6] bn; K − 1 bn; K bn; N b − 1
Γm;n a1 m;ne jϕ1 m;n a2 m;ne jϕ2 m;n ηm;n; (1) (4)
where a1 and a2 are the amplitudes of two signal components, The values of J and K must be selected, such that [22]
ϕ1 and ϕ2 are their corresponding interference phases, m ∈
0; M − 1 and n ∈ 0; N − 1 are the indices along rows and 1
M b 1 ≥ J ≥ N s 1
columns, respectively, and η represents complex additive white 2
Gaussian noise. The real part of Γm; n represents a moire 1
N 1 ≥ K ≥ N s 1; (5)
fringe pattern that encodes information about the multiple 2 b
interference phases ϕ1 and ϕ2 .
with N s denoting the number of signal components present in
In our method, Γm; n is divided into a number of small
the noisy input. In our case, N s is equal to two. The enhanced
blocks of size M b × N b , over which the field is approximated
matrix, B e , is decomposed using SVD into singular values and
as a sum of complex exponential signals. Let a given block be
singular vectors as
represented by
2 3 Be USV H ; (6)
b0; 0 b0; 1 b0; N b − 1
6 7 where H represents the Hermitian transpose operation, U and
6 b1; 0 b1; 1 b1; N b − 1 7
6 7
B6 .. .. .. .. 7; V contain left and right singular vectors, respectively, and S is
6 7
4 . . . . 5 a diagonal matrix containing singular values arranged in
bM b − 1; 0 bM b − 1; 1 bM b − 1; N b − 1 descending order.
The main advantage of applying SVD here is that the largest
(2) singular value always corresponds to the signal component
where the individual element bi; j Γi; j. In order to make having the largest amplitude. Further, the singular values cor-
the signal extraction process more robust against noise, we use responding to noise are always less than the values correspond-
ing to signal components. In other words, in Eq. (1), if we
assume that the amplitude a1 is greater than a2 , the component
corresponding to a1 can be separated by retaining only the
highest singular value and making all the remaining elements
in S equal to zero. Experimentally, the amplitude inequality or
discrimination criterion can be achieved by placing a neutral
density filter in one of the arms of the dual-wave digital holo-
graphic moire setup to attenuate the corresponding compo-
nent’s strength or amplitude. The new enhanced matrix is
converted back to a block matrix of size M b × N b using
Eqs. (3) and (4), which now contains only the first signal com-
ponent. The separated first component is subtracted from origi-
nal block B, and the SVD procedure is applied for retrieving the
Fig. 1. Schematic diagram of digital holographic interferometry second component. Subsequently, the process is repeated for all
setup to record moire fringe patterns. blocks. If the two signal components have significant spectral
1946 Vol. 57, No. 8 / 10 March 2018 / Applied Optics Research Article
overlap, the aforesaid component separation procedure may be RefΓc ng
repeated for a few iterations to improve accuracy. yn ; (14)
ImfΓc ng
Subsequent to the separation of the individual components,
we apply the state space formulation to extract the multiple an cosϕn x 0 cosx 1
phases and their derivatives. We divide each row of the sepa- hXn ; (15)
an sinϕn x 0 sinx 1
rated component into L number of windows and model the
phase in each window as a polynomial of degree P. For a and
given window, the complex field of the component can be Refηc ng
expressed as wn ; (16)
Imfηc ng
Γc n ac ne jϕc n ηc n ∀ n ∈ 0; N L − 1; (7)
with Re and Im denoting real and imaginary parts.
where the subscript c indicates c-th component, N L is the num- If we define the coefficient matrix as,
ber of elements in each window, and ηc is complex additive Θ a α0 α1 αP T ; (17)
white Gaussian noise. For the sake of mathematical simplicity,
the subscript c will be dropped in further calculations without where a indicates the assumed constant amplitude, then the
loss of generality. The phase, modeled as a local polynomial of state-vector can be written as,
degree P such that P < N L , can be represented by Xn AnΘ; (18)
X
P where
ϕn αp np ; (8) 2 3
p0 1 0 0 0 0
6 7
with αp denoting the polynomial coefficients. Using Taylor’s 60 1 n n2 nP 7
6 7
series expansion, ϕn 1 can be expressed as 6 PnP−1 7
An 6 0 0 1 2n 7: (19)
6. .. .. .. .. .. 7
X
P 6. 7
1 4. . . . . . 5
ϕn 1 ϕp n; (9)
p0
p! 0 0 0 0 P!
with Thus, when the state vector is known, the coefficients
matrix is calculated as
X
P
1
ϕq n 1 ϕp n; (10) Θ A −1 nXn; (20)
pq
p − q!
from using which the phase and its derivative can be calculated.
where the superscript “q” represents the derivative of order q. We applied an extended Kalman filter (EKF) [24] as our
All the unknown parameters to be estimated are formulated as a state estimator, mainly because of its efficiency and ease of im-
state vector [23] given by plementation. Estimation using EKF is carried out in two steps,
namely, a predict step and an update step, each of which are
Xn an ϕn ϕ1 n ϕPn T ; (11)
represented by a set of matrix operations as given below.
where T denotes vector transpose. The dimension of the state Predict step:
vector is P 2 × 1.
X̂ − n FX̂n
We also assume that the amplitude variations are much
weaker than the variations in phase within the segment, such P− n FPnFT : (21)
that the amplitude can be approximated to be reasonably Update step:
constant with an 1 ≈ an ≈ a. Using Eqs. (9)–(11), we
can write Kn P− nHT nHnP− nHT n R−1
2 3 2 32 3
an 1 1 0 0 0 an ŷn hX̂ − n
6 7 6 1 76 7
6 ϕn 1 7 6 0 1 1!1 P! 7 6 7 X̂n X̂ − n Knyn − ŷn
6 7 6 76 ϕn 7
6 1 7 60 0 1 1 76 1 7
6 ϕ n 1 7 6 P−1! 76 ϕ n 7: Pn I − KnHnP− n: (22)
6 7 6 76 7
6 .. 7 6. . . . .. 7 6 7
6 . 7 6. . . . 76 ... 7 Here, X̂n is the estimated state vector, Pn is the covariance
4 5 4. . . . . 54 5 error matrix, and Kn is the Kalman gain matrix after “n” num-
ϕP n 1 0 0 0 1 ϕP
n ber of measurements. Further, X̂ − n and P− n denote a priori
The same can be written as state estimate and a priori error covariance matrix, respectively.
H is a Jacobian matrix calculated by linearizing the non-linear
Xn 1 FXn; (12) equation hXn around the predicted state vector X̂ − n. Hence,
where F is called the state transition matrix. Similarly, from ∂hXn
Eq. (7), we can write the measurement equation as Hn ; (23)
∂Xn XnX̂ − n
yn hXn wn; (13)
with ∂ denoting partial derivative. R contains variance of the
where measurement error. In practical situations, precise calculation
Research Article Vol. 57, No. 8 / 10 March 2018 / Applied Optics 1947
of measurement error is very difficult. We initialized R with window. This process is repeated for all windows to get com-
large variance by over-estimating the measurement error. plete estimate of the phase and its derivative. To reduce the
Substituting Eqs. (11) and (15) in Eq. (23), we get small discontinuities along the edges of the windows, we used
overlapping segments and a phase-offsetting operation [23],
cosx̂ 1 −x̂ 0 sinx̂ 1 0 0 where the last pixel of a segment overlaps with the first pixel
Hn ; (24)
sinx̂ 1 x̂ 0 cosx̂ 1 0 0 of the next segment, and a suitable phase offset is added to
obtain a smooth phase map.
with x̂ 0 , x̂ 1 denoting the first two elements of X̂ − n. The proposed method does not require complex two-
The Kalman filter is initialized using X̂ − 0 dimensional phase unwrapping operations, as direct phase es-
jλ0j ∠λ0 0 0 and P− 0 diag 1 1 1 timates are obtained. Note that the proposed method relies on
as the prior estimate and prior error covariance matrix, respec- one-dimensional processing of the moire signal by using row-
tively. The predict and update steps are applied for each or column-based segmenting, and applying polynomial phase
measurement in a given window. After considering all mea- approximation within the segments. This approach reduces the
surements, the polynomial coefficients are calculated using overall computational overhead as compared to scanning two-
Eq. (20). Subsequently, the phase and derivative are obtained dimensional window-based moire fringe processing. Further,
from the estimated polynomial coefficients for the given it needs to be emphasized that due to the local modeling of
the phase as a polynomial, the derivative information can also
be easily ascertained once the polynomial coefficients are
computed. The feasibility of simultaneous estimation of multi-
ple phase derivatives is an added advantage of the proposed
method.
3. RESULTS
We simulated a complex interference moire field of size 256 ×
(a) 256 pixels using Eq. (1). The amplitudes of the signal compo-
nents were chosen as, a1 2 and a2 1. A complex zero-
mean white Gaussian noise was added with intensity-based
signal to noise ratio (SNR) equal to 15 dB. The SNR is com-
puted using the ratio of signal power and noise power. The
simulated moire fringe is shown in Fig. 2(a). Component ex-
traction was carried out using multiple blocks of size 32 × 32.
We used J K 10 and 10 iterations for component sepa-
ration. EKF was applied by dividing each row of separated
(b) (c) components into L 2 windows. The phase in each window
was modeled as a second-order polynomial.
(d) (e)
(a) (b)
(f) (g)
Fig. 2. (a) Simulated moire fringe pattern with two signal compo-
(c) (d)
nents and additive white Gaussian noise. Estimated phases ϕ1 m; n
and ϕ2 m; n in radians using the proposed method in (b) and (c). Fig. 3. Estimated phase derivatives ∂ϕ1∂n
m;n
and ∂ϕ2∂n
m;n
in radians/
Corresponding estimation error in radians in (d) and (e). Mesh plots pixel in (a) and (b). Corresponding errors in radians/pixel in (c)
for estimated phases in (f ) and (g). and (d).
1948 Vol. 57, No. 8 / 10 March 2018 / Applied Optics Research Article
Table 1. RMSE of Estimated Phases and Their We also tested the proposed method for demodulating the
Derivatives for Different SNRs experimental moire fringe pattern obtained in DHI. We used a
RMSE in Estimation of:
dual-beam setup similar to the one shown in Fig. 1. For experi-
ment, a circularly clamped object was deformed by applying
∂ϕ1 m;n ∂ϕ2 m;n
SNR ϕ1 m;n ϕ2 m;n ∂n ∂n external load, and holograms were recorded prior to and after
(dB) rad rad rad∕pixel rad∕pixel deformation. For the setup, we used a Coherent Verdi laser
0 0.7384 1.4559 0.0228 0.0658 with wavelength of 532 nm as the illumination source and a
5 0.6387 0.6674 0.0172 0.0261 Sony XCL-U1000 CCD digital camera for recording holo-
10 0.5911 0.6243 0.0162 0.0312 grams. The hologram reconstruction and computation of moire
15 0.4638 0.5769 0.0157 0.0265 fringes were performed according to the procedure outlined in
20 0.4546 0.4827 0.0156 0.0262 our previous work [6]. The experimental moire fringe pattern
is shown in Fig. 4(a). For component extraction, we chose
M b N b 30 and J K 10. To extract the phase maps,
The estimates of phases ϕ1 m; n and ϕ2 m; n in radians
we applied a local second-order polynomial phase approxima-
are shown in Figs. 2(b) and 2(c). The corresponding estimation
tion with a window size of N L 128. The estimated first and
errors in radians are shown in Figs. 2(d) and 2(e). For the sake
second phase maps in radians using the proposed method are
of visualization, the mesh plots for the two estimated phases are
shown in Figs. 4(b) and 4(c), respectively. For better visualiza-
also shown in Figs. 2(f )–2(h). The estimated phase derivatives
tion, the cosine fringe patterns cos ϕ1 and cos ϕ2 correspond-
and their corresponding errors are shown in Fig. 3. The overall
ing to individual estimated phase maps are also shown in
computational time required by the proposed method to esti-
Figs. 4(d)–4(e).
mate the two phases and their derivatives from the moire fringe
pattern was approximately 28 s.
To validate the performance of the proposed method, we 4. DISCUSSION
computed the root mean square errors (RMSE) for phase and The proposed method offers an elegant approach for extracting
phase derivative estimation with varying noise levels or SNR. multiple phase maps associated with a digital holographic moire
These results are shown in Table 1. It is evident that the pro- fringe pattern. Using the SVD-based component separation
posed method shows decent performance even for SNR as and subsequent local polynomial phase-modeling-based EKF
low as 0 dB. operation, the method eliminates the need of applying spectral
filtering and complex two-dimensional phase unwrapping algo-
rithms and the requirement of prior phase information for
initializing the filter. In addition, the proposed method offers
the possibility of simultaneously estimating the multiple phase
derivatives in the given framework without extra effort.
However, the proposed method also has some limitations.
Foremost, amplitude inequality for the multiple components
has to be ensured for reliable working of the proposed tech-
nique. In addition, as the offset operation is carried at the
(a) end of each segment to minimize the phase discontinuity be-
tween adjacent segments, errors may arise due to improper off-
sets. Hence, even small offset errors might propagate as more
segments are sequentially joined. This can be remedied by us-
ing segments with more overlapping regions to minimize offset-
based errors, albeit at the cost of increased complexity and
higher computational overhead. The accuracy of the proposed
technique would also be affected by non-uniform illumination,
spatially varying direct current (dc) term, and presence of severe
(b) (c) noise and faulty regions such as shadows in the holographic
moire fringe patterns. Further, the choice of the length of seg-
ment also plays a role in estimation accuracy. Normally, the
segment length should be large enough so that state space iter-
ations can converge to provide reliable estimates for the poly-
nomial coefficients. However, a very large value of segment
length would lead to an increase in computational cost and also
induce local polynomial phase approximation errors. Regarding
future directions, the temporal performance of the proposed
(d) (e) technique has to be quantified by recording and processing sev-
Fig. 4. (a) Experimental digital holographic moire fringe pattern. eral holograms temporally. This would provide insights about
Estimated first and second phases in (b) and (c). Corresponding cosine the temporal evolution of errors for the technique. In addition,
fringes in (d) and (e). the impact of spatially varying dc terms in the holograms needs
Research Article Vol. 57, No. 8 / 10 March 2018 / Applied Optics 1949
to be investigated. In the current technique, we applied the 10. M. Takeda, H. Ina, and S. Kobayashi, “Fourier-transform method of
fringe-pattern analysis for computer-based topography and interfer-
mean subtraction method [1] for dc removal, which would be
ometry,” J. Opt. Soc. Am. 72, 156–160 (1982).
inadequate in the case of dc term with high spatial variations. 11. Q. Kemao, “Two-dimensional windowed Fourier transform for fringe
pattern analysis: principles, applications and implementations,” Opt.
5. CONCLUSION Lasers Eng. 45, 304–317 (2007).
12. L. Watkins, S. Tan, and T. Barnes, “Determination of interferometer
The paper proposed a method for demodulating moire fringe phase distributions by use of wavelets,” Opt. Lett. 24, 905–907 (1999).
patterns obtained in digital holographic interferometry. The 13. S. S. Gorthi, G. Rajshekhar, and P. Rastogi, “Strain estimation in dig-
authors believe that the proposed method has good application ital holographic interferometry using piecewise polynomial phase
approximation based method,” Opt. Express 18, 560–565 (2010).
potential for multi-dimensional deformation analysis in preci-
14. A. M. S. Maallo and P. F. Almoro, “Numerical correction of optical vor-
sion metrology and non-destructive testing. tex using a wrapped phase map analysis algorithm,” Opt. Lett. 36,
1251–1253 (2011).
Funding. Science and Engineering Research Board (SERB) 15. F. C. I. Catalan, A. M. S. Maallo, and P. F. Almoro, “Fringe analysis
India (ECR/2015/000180). and enhanced characterization of sub-surface defects using fringe-
shifted shearograms,” Opt. Commun. 285, 4223–4226 (2012).
16. M. Servin, J. Marroquin, and F. Cuevas, “Demodulation of a single
REFERENCES interferogram by use of a two-dimensional regularized phase-tracking
1. U. Schnars and W. P. Jüptner, “Digital recording and numerical technique,” Appl. Opt. 36, 4540–4548 (1997).
reconstruction of holograms,” Meas. Sci. Technol. 13, R85–R101 17. G. Rajshekhar, S. S. Gorthi, and P. Rastogi, “Estimation of multiple
(2002). phases from a single fringe pattern in digital holographic interferom-
2. G. Pedrini, Y.-L. Zou, and H. Tiziani, “Simultaneous quantitative etry,” Opt. Express 20, 1281–1291 (2012).
evaluation of in-plane and out-of-plane deformations by use of a 18. R. Kulkarni and P. Rastogi, “Iterative signal separation based multiple
multidirectional spatial carrier,” Appl. Opt. 36, 786–792 (1997). phase estimation in digital holographic interferometry,” Opt. Express
3. C. Kohler, M. R. Viotti, and G. A. Albertazzi, Jr., “Measurement of 23, 26842–26852 (2015).
three-dimensional deformations using digital holography with radial 19. R. Kulkarni and P. Rastogi, “Digital holographic moire for the direct
sensitivity,” Appl. Opt. 49, 4004–4009 (2010). and simultaneous estimation of strain and slope fields,” Opt.
4. G. Pedrini and H. Tiziani, “Quantitative evaluation of two-dimensional Express 22, 23192–23201 (2014).
dynamic deformations using digital holography,” Opt. Laser Technol. 20. R. G. Waghmare, D. Mishra, and G. S. Subrahmanyam, “UKF based
29, 249–256 (1997). multi-component phase estimation in digital holographic Moiré,”
5. P. Picart, E. Moisson, and D. Mounier, “Twin-sensitivity measurement in Fourth National Conference on Computer Vision, Pattern
by spatial multiplexing of digitally recorded holograms,” Appl. Opt. 42, Recognition, Image Processing and Graphics (NCVPRIPG) (IEEE,
1947–1957 (2003). 2013), pp. 1–4.
6. G. Rajshekhar, S. S. Gorthi, and P. Rastogi, “Simultaneous multi- 21. R. Kulkarni and P. Rastogi, “Simultaneous estimation of multiple
dimensional deformation measurements using digital holographic phases in digital holographic interferometry using state space analy-
moire,” Appl. Opt. 50, 4189–4197 (2011). sis,” Opt. Lasers Eng. (2017), in press.
7. Y. Morimoto, T. Matui, M. Fujigaki, and A. Matsui, “Three-dimensional 22. Y. Hua, “Estimating two-dimensional frequencies by matrix enhance-
displacement analysis by windowed phase-shifting digital holographic ment and matrix pencil,” IEEE Trans. Signal Process. 40, 2267–2280
interferometry,” Strain 44, 49–56 (2008). (1992).
8. E. Kolenovic, W. Osten, R. Klattenhoff, S. Lai, C. von Kopylow, and 23. G. Rajshekhar and P. Rastogi, “Phase estimation using a state-space
W. Jüptner, “Miniaturized digital holography sensor for distal three- approach based method,” Opt. Lasers Eng. 51, 1004–1007 (2013).
dimensional endoscopy,” Appl. Opt. 42, 5167–5172 (2003). 24. R. G. Brown and P. Y. Hwang, Introduction to Random Signals and
9. G. Rajshekhar and P. Rastogi, “Fringe analysis: premise and perspec- Applied Kalman Filtering: With MATLAB Exercises and Solutions
tives,” Opt. Lasers Eng. 50, iii–x (2012). (Wiley, 1997).