- Email: [email protected]
Fast phase retrieval in slightly off-axis digital holography
Optics and Lasers in Engineering 97 (2017) 9–18
Contents lists available at ScienceDirect
Optics and Lasers in Engineering journal homepage: www.elsevier.com/locate/optlaseng
Fast phase retrieval in slightly off-axis digital holography Zhi Zhong, Hongyi Bai, Mingguang Shan∗, Yabin Zhang, Lili Guo College of Information and Communication Engineering, Harbin Engineering University, Harbin, Heilongjiang 150001, PR China
a r t i c l e
i n f o
Keywords: Complex encoding Digital holography Fast phase retrieval Slightly off-axis Spatial multiplexing Spectral cropping
a b s t r a c t In this study, three efficient algorithms are proposed for fast phase retrieval in slightly off-axis digital holography using spectrum cropping, spatial multiplexing, and complex encoding. In the first algorithm, the real spectral order of the subtracted hologram is filtered and cropped, and the number of pixels is decreased in the subsequent retrieval operations. In the second algorithm, two sequential subtracted holograms are digitally phase shifted and spatial multiplexed into one synthetic hologram, and thus only one inverse Fourier transformation is then required. In the third algorithm, two sequential subtracted holograms are encoded separately into the real part and the imaginary part of a complex hologram. Two cross-correlations can be used to reconstruct the phase, thereby improving the utilization of the spectrum. The three new algorithms speed up our previously proposed retrieval method with the assistance of specimen-free holograms. Our experiments demonstrated the validity and improved time requirements of the proposed methods. © 2017 Elsevier Ltd. All rights reserved.
1. Introduction Digital holography (DH) has emerged as an important tool for achieving quantitative phase analyses of biological specimens, surface measurements, and micro-structures [1–4]. Due to the reliance of DH on interference, three types of methods in DH are usually used: on-axis, off-axis and slightly off-axis. The on-axis methods can fully utilize the space–bandwidth product of the DH system. However, to achieve high phase retrieval accuracy, phase-shifting methods must record several holograms (typically ≥3) in time sequence or in a single shot [5–8]. Unlike phase-shifting methods, the off-axis methods just need one hologram in a single shot, which make them suitable for faster imaging by fully exploiting the field of view (FOV) of the camera. However, these benefits come at the cost of the space–bandwidth product [9–13]. By integrating off-axis and phase-shifting to record two phase-shifted slightly off-axis holograms, the slightly off-axis methods provide an intermediate solution between the on-axis and off-axis methods [14–19]. After applying a subtraction operation to the two slightly off-axis holograms to eliminate the DC term, digital phase reconstruction can also be performed using a Fourier transform algorithm with the slightly off-axis DHs. However, similar to the off-axis DHs, digital processing requires large amounts of computations so it must be conducted offline. An improved phase retrieval method was proposed by our group based on the offline prior acquisition of specimen-free holograms [18,19], but it is still necessary to speed up the phase retrieval process for slightly offaxis DHs. The retrieval operations are usually the same for different pixels in one hologram. Parallel computing techniques using graphic processing ∗
Corresponding author. E-mail address: [email protected] (M. Shan).
http://dx.doi.org/10.1016/j.optlaseng.2017.05.004 Received 19 December 2016; Received in revised form 2 May 2017; Accepted 2 May 2017 0143-8166/© 2017 Elsevier Ltd. All rights reserved.
units (GPUs) are then introduced to achieve higher retrieval efficiency [20,21]. However, special graphic units and programming skills are required. Considering that Fourier transform operation is the main timeconsuming process, Girshovitz and Shaked retained only the real spectral order and omit all the other redundant data during the retrieval process [11]. They then achieved real time phase reconstruction with 1 megapixel. To further speed up this process, Sha et al. spatially multiplexed four holograms into one synthetic hologram [12], while Girshovitz and Shaked encoded four sequential holograms into one complex hologram, and then used both the real and twin image parts [13]. However, all of these methods are only suitable for off-axis DHs. In this study, we present three more efficient algorithms for phase retrieval in slightly off-axis DHs. Those algorithms combine spectrum cropping, spatial multiplexing, and complex encoding with our previously proposed method. They speed up phase retrieval in slightly offaxis DHs, as well as compensating for the disturbing phase caused by the phase shift 𝛼, carrier frequency, and background phase. Finally, we demonstrated the improved performance of our proposed approaches based on simulations, and by processing phase maps of a phase plate and water drops. 2. Reconstruction algorithms In slightly off-axis DH, two measured holograms with a phase shift 𝛼 have the forms: [( )] 𝐻1 (𝑥, 𝑦) = 𝑎(𝑥, 𝑦) + 0.5𝑏(𝑥, 𝑦) exp 𝑖 𝜑𝑡 + 𝜑𝑏𝑔 +2𝜋𝑓𝑥 𝑥 [ ( )] + 0.5𝑏(𝑥, 𝑦) exp −𝑖 𝜑𝑡 + 𝜑𝑏𝑔 +2𝜋𝑓𝑥 𝑥 ,
(1)
Z. Zhong et al.
Optics and Lasers in Engineering 97 (2017) 9–18
Fig. 1. Algorithm A: the algorithm based on an offline prior acquisition of specimen-free holograms.
[( )] 𝐻2 (𝑥, 𝑦) = 𝑎(𝑥, 𝑦) + 0.5𝑏(𝑥, 𝑦) exp 𝑖 𝜑𝑡 + 𝜑𝑏𝑔 +2𝜋𝑓𝑥 𝑥 exp (𝑖𝛼) [ ( )] + 0.5𝑏(𝑥, 𝑦) exp −𝑖 𝜑𝑡 + 𝜑𝑏𝑔 +2𝜋𝑓𝑥 𝑥 exp (−𝑖𝛼),
is subtracted out, so the maximum circle in the BPF can reach the zerofrequency point [14]. Thus, for a spectral matrix containing N × N pixels, the maximum diameter of the circle in the BPF can be selected as N/2. A4: Two-dimensional inverse FFT (2-D IFFT): Convert the result obtained in step A3 back to the image domain by using a 2-D IFFT to product an N × N complex matrix containing the phase distribution of the specimen and disturbance. The disturbing phase includes phase shift, carrier information and background phase. A5: Load the result of prior measurement: To compensate for the disturbing phase, obtain two phase shifted holograms without a specimen before the experiment and then process them according to steps A1–A4, as described above. The result containing N × N pixels is actually stored before the measurement and it is only loaded to help retrieving phase during the reconstruction process. A6: Specimen phase retrieval: Divide the result obtained from step A4 by that from step A5 to eliminate the disturbing phase. Finally, implement an arc tangent operation to obtain the phase distribution of the specimen.
(2)
where a(x, y) is the DC term, b(x, y) is the modulation term, 𝜑t is the phase distribution of specimen at time t, 𝜑bg is the background phase caused by the aberration and the noise of the system, fx is carrier frequency in the x direction and 2𝜋fx x is its corresponding phase tilt. After applying subtraction, we can eliminate the DC term and obtain: {[ 𝑆 = 𝐻1 − 𝐻2 = 0.5𝑏(𝑥, 𝑦)
] [( )] } 1 − exp (𝑖𝛼) exp 𝑖 𝜑𝑡 + 𝜑𝑏𝑔 +2𝜋𝑓𝑥 𝑥 [ ] [ ( )] . + 1 − exp (−𝑖𝛼) exp −𝑖 𝜑𝑡 + 𝜑𝑏𝑔 +2𝜋𝑓𝑥 𝑥 (3)
The phase distribution 𝜑t of a specimen can be retrieved from the subtracted hologram using Fourier transform algorithm [17]. To increase the retrieval speed, we suggested a faster retrieval algorithm based on the offline prior acquisition of specimen-free subtracted holograms [18,19]. The faster algorithm can also compensate for the disturbing phase caused by the phase shift 𝛼, carrier frequency, and background phase. In the following, we first review the algorithm, and then propose three more efficient algorithms.
2.2. Algorithm B: spectrum cropping algorithm Shaked et al. noted that when selecting the real spectral order, there is no need to use a BPF matrix with the same size as the original images (N × N). The filtering operation can be performed by decreasing the number of pixels and the 2-D FFT operation can be followed by using a cropped matrix with the same size. For slightly off-axis DH, we can crop the real spectral order using a BPF with N/2 × N/2 pixels and continue processing with a matrix that is four times smaller. Fig. 2 shows the digital process, which is described in Algorithm B. The algorithm comprises the following step. B1: Holograms subtraction: This is the same as step A1. B2: 2-D FFT: This is the same as step A2. B3: DBPF and cropping: Apply a BPF containing N/2 × N/2 pixels to the area (N/2 + 1:N, N/4 + 1:3 N/4) of the spectrum and crop the real spectral order. B4: 2-D IFFT: Convert the result obtained from step B3 back into the image domain by using a 2-D IFFT to produce an N/2 × N/2 complex matrix containing the phase distribution of the specimen and disturbance.
2.1. Algorithm A: algorithm based on the offline prior acquisition of specimen-free holograms When the phase variance of a specimen is less than 2𝜋, we previously proposed a simple retrieval algorithm assisted by a prior measurement without any specimen [18,19]. This algorithm can achieve a fast phase retrieval speed and it comprises the following steps, as shown in Fig. 1. A1: Holograms subtraction: Subtract the two phase-shifted holograms containing N × N real pixels to eliminate the DC term and obtain a subtracted hologram. A2: Two-dimensional fast Fourier transform (2-D FFT): Convert the subtracted hologram into the spatial-frequency domain using a 2-D FFT, thereby obtaining a matrix containing N × N complex pixels. A3: Digital band pass filtering (DBPF): Use a band pass filter (BPF) with a size of N × N pixels and digitally filter out the twin spectral order, but retain the real spectral order. In slightly off-axis DH, the DC order 10
Z. Zhong et al.
Optics and Lasers in Engineering 97 (2017) 9–18
Fig. 2. Algorithm B: the spectrum cropping algorithm.
B5: Load the result of prior measurement: Acquire two specimenfree phase-shifted holograms before the experiments and process them according to steps B1–B4 to obtain a matrix containing N/2 × N/2 pixels. Load the result obtained previously during the reconstruction process. B6: Specimen phase retrieval: Divide the two results obtained from steps B4 and B5, and apply an arc tangent operation to obtain the phase distribution of the specimen containing N/2 × N/2 pixels. Finally, enlarge the N/2 × N/2 phase to the size of N × N phase matrix. 2.3. Algorithm C: spatial multiplexing algorithm
(7)
( ) { { } } 𝑐2 ′ 𝑥′ , 𝑦′ = IFT FT 𝑆𝑆 ⋅ BPF2 ( )[ ] = 0.5𝑏′ 𝑥′ , 𝑦′ 1 − exp (𝑖𝛼) [( ′ )] exp 𝑖 𝜑𝑏𝑔 +2𝜋𝑓𝑥 ′ 𝑥′ + 2𝜋𝑓𝑦𝑚 ′ 𝑦′ .
(8)
Thus, the phase retrieved can be expressed as: [ ( )] 𝑐1 𝑥′ , 𝑦′ ( ) ′ ′ ′ 𝜑𝑡 𝑥 , 𝑦 = arctan . 𝑐 1 ′ (𝑥 ′ , 𝑦 ′ ) [ ( )] 𝑐2 𝑥′ , 𝑦′ ( ) 𝜑𝑡+1 ′ 𝑥′ , 𝑦′ = arctan . 𝑐 2 ′ (𝑥 ′ , 𝑦 ′ )
According to the shift theorem for the Fourier transform, a phase shift in the spatial domain can be interpreted as a frequency shift in the frequency domain. Thus, two subtracted holograms can be digitally phase shifted in the y direction and spatially multiplexed into one image [12] as follows: ( ) ( ) 𝑆𝑆 = 𝑆𝑡 exp −𝑖2𝜋𝑓𝑦𝑚 𝑦 + 𝑆𝑡+1 exp 𝑖2𝜋𝑓𝑦𝑚 𝑦 (4)
(9)
(10)
Finally, 𝜑t and 𝜑t+1 containing N × N pixels can be acquired by enlarging 𝜑t ′ and 𝜑t+1 ′ containing N/2 × N/2 pixels. A flowchart illustrating the algorithm described above is shown in Fig. 3 and its steps are as follows. C1: Holograms subtraction: This is the same as step A1 but for two sequential subtracted holograms St and St +1 . C2: Synthesis of subtracted holograms: Digitally phase shift St and St +1 , and sum them to obtain one synthetic hologram SS . C3: 2-D FFT: Convert SS into the spatial-frequency domain using a 2D FFT and obtain a matrix containing N × N complex pixels. This step can actually be replaced by two steps to yield a 1.5-D FFT algorithm where a 1-D FFT in rows is first implemented to produce N × N pixels, and a 1-D FFT in columns is calculated to obtain an N × N/2 array containing only the real images from the two subtractions [11]. To facilitate convenient comparisons with other methods, we use a 2-D FFT in this study. C4: DBPF and cropping: Apply two BPFs containing N/2 × N/2 pixels to cover the areas (N/2 + 1:N, 1:N/2) and (N/2 + 1:N, N/2 + 1: N) of the spectrum, and crop the corresponding real spectral orders of 𝜑t and 𝜑t +1 . C5: 2-D IFFT: Convert the result obtained from step C4 back into the image domain by using a 2-D IFFT to yield two N/2 × N/2 complex matrices containing the phase distributions of 𝜑t and 𝜑t +1 . C6: Load the results of prior measurements: This is the same as step B5 but for two sequential phase distributions each containing N/2 × N/2 pixels.
where St and St +1 are two sequential subtracted holograms and fym is a spatial phase shift in the y direction. For a typical slightly off-axis DH, the bandwidth B of the subtracted hologram is equal to two times that of the recorded specimen. Thus, in order to separate the spectra of St and St +1 , fym should not be less than B/4. According to Eq. (3) and Algorithm B, the cropped real spectral orders of St and St +1 can be extracted as: ( ) { { } } 𝑐1 𝑥′ , 𝑦′ = IFT FT 𝑆𝑆 ⋅ BPF1 ( )[ ] = 0.5𝑏′ 𝑥′ , 𝑦′ 1 − exp (𝑖𝛼) [( ′ )] (5) exp 𝑖 𝜑𝑡 + 𝜑𝑏𝑔 ′ +2𝜋𝑓𝑥 ′ 𝑥′ − 2𝜋𝑓𝑦𝑚 ′ 𝑦′ . ( ) { { } } 𝑐2 𝑥′ , 𝑦′ = IFT FT 𝑆𝑆 ⋅ BPF2 ( )[ ] = 0.5𝑏′ 𝑥′ , 𝑦′ 1 − exp (𝑖𝛼) [( )] exp 𝑖 𝜑𝑡+1 ′ + 𝜑𝑏𝑔 ′ +2𝜋𝑓𝑥 ′ 𝑥′ + 2𝜋𝑓𝑦𝑚 ′ 𝑦′ .
( )[ ] = 0.5𝑏′ 𝑥′ , 𝑦′ 1 − exp (𝑖𝛼) [( ′ )] exp 𝑖 𝜑𝑏𝑔 +2𝜋𝑓𝑥 ′ 𝑥′ − 2𝜋𝑓𝑦𝑚 ′ 𝑦′ .
(6)
The corresponding cropped real spectral orders of free-specimens can also be extracted as: ( ) { { } } 𝑐1 ′ 𝑥′ , 𝑦′ = IFT FT 𝑆𝑆 ⋅ BPF1 11
Z. Zhong et al.
Optics and Lasers in Engineering 97 (2017) 9–18
Fig. 3. Algorithm C: the spatial multiplexing algorithm.
C7: Specimen phase retrieval: This is the same as step B6 but for two sequential phase distributions.
( ) { { } } 𝑑2 𝑥′ , 𝑦′ = IFT FT 𝑆𝐶 ⋅ BPF2 {[ } ] [ ( )] ′ ( ) 1 − exp (−𝑖𝛼) exp −[𝑖 𝜑(𝑡 ′ + 𝜑𝑏𝑔 ′ +2𝜋𝑓𝑥 𝑥′ )] . [ ] = 0.5𝑏′ 𝑥′ , 𝑦′ + 1 − exp (−𝑖𝛼) exp −𝑖 𝜑𝑡+1 ′ + 𝜑𝑏𝑔 ′+2𝜋𝑓𝑥 ′ 𝑥′ − 𝜋2
2.4. Algorithm D: complex encoding algorithm Like off-axis holograms, the slightly off-axis holograms also have no imaginary parts, which make the Fourier transformation operation inefficient. Thus, it can employ the imaginary part to encode two sequential subtracted holograms into a complex hologram [13]. The complex hologram is encoded as follows. 𝑆𝐶 = 𝑆𝑡 + 𝑖𝑆𝑡+1
(13) After applying the conjugate operation to Eq. (13), we can obtain the following results by the summation and subtraction of Eqs. (12) and (13):
(11)
( )[ ] [( )] ( ) ( ) 𝑏′ 𝑥′ , 𝑦′ 1 − exp (𝑖𝛼) exp 𝑖 𝜑𝑡 ′ + 𝜑𝑏𝑔 ′ +2𝜋𝑓𝑥 ′ 𝑥′ = 𝑑1 𝑥′ , 𝑦′ +𝑑2∗ 𝑥′ , 𝑦′ ,
After a 2-D FFT operation, the spectrum can be obtained for SC . The real and imaginary orders containing N/2 × N/2 pixels can then be obtained by BPF and cropping, and the 2-D IFFT calculation. This process can be expressed as:
(14)
( ) { { } } 𝑑1 𝑥′ , 𝑦′ = IFT FT 𝑆𝐶 ⋅ BPF1 {[ } ] [( )] ′ ( ) 1 − exp (𝑖𝛼) exp 𝑖[𝜑(𝑡 ′ + 𝜑𝑏𝑔 ′ +2𝜋𝑓𝑥 𝑥′ )] ′ ′ ′ [ ] . = 0.5𝑏 𝑥 , 𝑦 + 1 − exp (𝑖𝛼) exp 𝑖 𝜑𝑡+1 ′ + 𝜑𝑏𝑔 ′ +2𝜋𝑓𝑥 ′ 𝑥′ + 𝜋2
( ) ( ) ( )[ ] [( )] 𝑑1 𝑥′ , 𝑦′ −𝑑2∗ 𝑥′ , 𝑦′ 𝑏′ 𝑥′ , 𝑦′ 1 − exp (𝑖𝛼) exp 𝑖 𝜑𝑡+1 ′ + 𝜑𝑏𝑔 ′+2𝜋𝑓𝑥 ′ 𝑥′ = , 𝑖 (15) where, “∗ ” denotes the conjugate operation.
(12) 12
Z. Zhong et al.
Optics and Lasers in Engineering 97 (2017) 9–18
Fig. 4. Algorithm D: the complex encoding algorithm.
Repeating the operation above, we can obtain the same specimenfree results as follows: ( )[ ] [( )] ( ) ( ) 𝑏′ 𝑥′ , 𝑦′ 1 − exp (𝑖𝛼) exp 𝑖 𝜑𝑏𝑔 ′ +2𝜋𝑓𝑥 ′ 𝑥′ = 𝑑1 ′ 𝑥′ , 𝑦′ + 𝑑2 ′∗ 𝑥′ , 𝑦′ ,
D1: Holograms subtraction: This is the same as step C1. D2: Encoding subtracted holograms: According to Eq. (11), sum the two subtracted holograms St and St +1 to obtain a complex hologram SC , where St is the real part of SC and St +1 is the imaginary part of SC . The complex hologram SC contains N × N pixels. D3: 2-D FFT: Convert SC into the spatial-frequency domain using a 2-D FFT to obtain a matrix containing N × N complex pixels. D4: DBPF and cropping: Apply two BPFs containing N/2 × N/2 pixels to cover the area (1:N/2, N/4 + 1:3 N/4) and (N/2 + 1:N, N/4 + 1:3 N/4) of the spectrum, and crop the two cross-correlations. D5: 2-D IFFT: Convert the results obtained from step D4 back into the image domain using 2-D IFFTs, thereby producing two N/2 × N/2 complex matrices as shown in Eqs. (12) and (13). D6: Load the results of prior measurements: Load the specimenfree results obtained according to steps D1–D5 following Eqs. (16) and (17) to produce two N/2 × N/2 complex matrices. D7: Specimen phase retrieval: Process the results obtained above using Eqs. (18) and (19) to produce two sequential phase distributions each containing N/2 × N/2 pixels. Finally, enlarge the phase to the N × N phase matrix.
(16) ( ) ( ) ( )[ ] [( )] 𝑑1 ′ 𝑥′ , 𝑦′ − 𝑑2 ′∗ 𝑥′ , 𝑦′ 𝑏′ 𝑥′ , 𝑦′ 1 − exp (𝑖𝛼) exp 𝑖 𝜑𝑏𝑔 ′ +2𝜋𝑓𝑥 ′ 𝑥′ = 𝑖 (17) Thus, the retrieved phases of 𝜑t ′ and 𝜑t +1 ′ containing N/2 × N/2 pixels can be expressed as: [ ( ′ ′) ( )] 𝑑1 𝑥 , 𝑦 + 𝑑2∗ 𝑥′ , 𝑦′ ′ 𝜑𝑡 = arctan . (18) 𝑑 1 ′ (𝑥 ′ , 𝑦 ′ ) + 𝑑 2 ′ ∗ (𝑥 ′ , 𝑦 ′ ) [ ( ) ( )] 𝑑1 𝑥′ , 𝑦′ − 𝑑2∗ 𝑥′ , 𝑦′ 𝜑′𝑡+1 = arctan 𝑖 ′ . (19) 𝑑1 (𝑥′ , 𝑦′ ) + 𝑑2 ′ ∗ (𝑥′ , 𝑦′ ) After an enlarging operation, we can finally obtain 𝜑t and 𝜑t +1 containing N × N pixels. A flowchart illustrating the algorithm described above is shown in Fig. 4 and its steps are as follows: 13
Z. Zhong et al.
Optics and Lasers in Engineering 97 (2017) 9–18
Fig. 5. Simulations. (a) Specimen phase; (b) specimen hologram with 𝛼 = 0; (c) specimen hologram with 𝛼 = 𝜋/2; (d) specimen-free hologram with 𝛼 = 0; and (e) specimen-free hologram with 𝛼 = 𝜋/2.
of 0° and 90°, phase shifts of 0 and 𝜋/2 were introduced into the holograms, respectively. Compared with the method employed by Ref. [19], the proposed setup could avoid the effects of the intensity ratio between different components of the object and reference beams, and as well as making full use of the FOV of the camera. In the experiments, a He–Ne laser with a wavelength of 𝜆 = 632.8 nm was used as the light source. Holograms were recorded by a CCD camera with 1600 × 1200 pixels (each pixel size is 4.4 μm × 4.4 μm). Fig. 8(a) and (b) show the holograms captured without specimens, and Fig. 8(c) and (d) show the holograms captured with a phase plate (refractive index n = 1.5168) with a step height of 580.22 nm measured by a BRUKER atomic force microscope, which could provide an optical path difference (OPD) of 2.98 rad at 𝜆 = 632.8 nm. The four holograms all comprised 640 × 640 pixels, and the spectra of H1 ′–H2 ′ and H1 –H2 are shown in Fig. 9. The results indicate that the DC term was eliminated, which increased the frequency coverage of the real spectral order. We used a BPF with 640 × 640 pixels for Algorithm A, and BPFs with 320 × 320 pixels for Algorithms B–D. However, to achieve a trade-off between resolution and spatial phase noise [26], we selected the pass area of the BPF using F = 0.6, where F = 1.0 corresponds to the diameter of the circle in Fig. 9 that touches the zero-frequency point. Fig. 10 shows the similar results retrieved using all four algorithms. To further confirm the validity of the algorithms, multiple water drops placed on part of a slide were also measured, where Fig. 11 shows the corresponding similar results retrieved using all four algorithms. To evaluate the run time of the different algorithms, holograms with different pixel sizes were retrieved with Matlab 2014a using a simple desktop (Intel i5-4430, 3 GHz, 8 GB RAM, where only a single core was used). After a large number of repetitions, the average processing times required to reconstruct one frame of phase containing 1024 × 1024 pixels using Algorithms A and B are shown in Table 2, and the times for two frames using Algorithms C and D are shown in Table 3. These results demonstrate that compared with Algorithm A, all the new algorithms could effectively speed up the retrieval process. However, it should be noted that the operations required to load prior results were very time-
Table 1 Error comparison for different algorithms.
ARE (%) RMSE (mrad)
Algorithm A
Algorithm B
Algorithm C
Algorithm D
0.0023 0.123
0.11 2.3
0.11 2.3
0.11 2.3
3. Simulations The proposed algorithms were verified based on numerical simulations. Fig. 5(a) shows a specimen phase of <2𝜋 rad with 1024 × 1024 pixels (each pixel size is 4.4 μm × 4.4 μm). For a slightly off-axis DH system [15], the carrier frequency fx was set to be 256/1024/(4.4 μm) [22]. The simulated specimen holograms that satisfied Eqs. (1) and (2) are shown in Fig. 5(b) and (c) with 𝛼 = 𝜋/2, respectively, and the corresponding specimen-free holograms are shown in Fig. 5(d) and (e). The phase distribution was then retrieved using Algorithms A–D and the results are shown in Fig. 6. It can be seen that despite spectrum cropping, the average relative error (ARE) and root mean square error (RMSE) [23] with the proposed algorithms, as shown in Table 1, were higher compared with those using Algorithm A, and all the proposed algorithms retrieved the specimen phase with high quality. In fact, ARE and RMSE results obtained using the proposed algorithms were far lower than those using the carrier estimation algorithms without the aid of specimen-free holograms [24,25]. 4. Experiment results In order to further demonstrate the efficiency of the algorithms, experiments were conducted with slightly off-axis DH microscopy, as shown in Fig. 7. Similar to the method employed in our previous work [19], the microscopy approach utilized a 45° linearly polarized object beam and a circularly polarized reference beam for interference. However, a polarizer P3 was placed before the CCD camera to replace the polarizing beam splitter used by Ref. [19]. By setting P3 with angles 14
Z. Zhong et al.
Optics and Lasers in Engineering 97 (2017) 9–18
Fig. 6. Simulation results. (a) Phase retrieved by Algorithm A; (b) phase retrieved by Algorithm B; (c) phase retrieved by Algorithm C; (d) phase retrieved by algorithm D; and (e) 1-D profiles along the dashed lines marked in Figs. 5(c) and 6(a)–(d) (all the color bars indicate rad). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 7. Experimental setup: P1, P2 and P3, polarizers; QW, quarter-wave plate; MO, microscopic objective; S, specimen; CL, collimated lens; L1 and L2, lenses with a same focal length; BS, beam-splitter; M1, reflective pinhole mirror; and RR, retro-reflector.
Table 2 Processing time for one frame of phase containing 1024 × 1024 pixels using Algorithms A and B (ms).
Algorithm A Algorithm B
2D-FFT
2D-IFFT
Load the result of prior measurement
Specimen phase retrieval
Others
Total time
24.347 24.236
24.068 7.839
94.754 24.639
16.471 11.521
13.409 7.359
173.049 75.594
15
Z. Zhong et al.
Optics and Lasers in Engineering 97 (2017) 9–18
Fig. 8. Holograms. (a, b) Holograms captured without specimens; (c, d) holograms captured with a phase plate. Table 3 Processing time for two frames of phase containing 1024 × 1024 pixels using Algorithms C and D (ms).
Algorithm C Algorithm D
Subtracted-holograms synthesizing
2D-FFT
2D-IFFT
Load the result of prior measurement
Specimen phase retrieval
Others
Total time
24.563 4.906
22.862 22.749
12.826 12.845
49.451 49.434
21.371 22.463
6.672 7.034
137.745 119.431
Table 4 Processing time for 100 frames of phase containing different pixels (s). Pixel size
Algorithm A
Algorithm B
Algorithm C
Algorithm D
512 × 512 768 × 768 1024 × 1024
1.520 3.660 7.298
1.001 2.378 4.436
0.978 2.323 4.269
0.732 1.777 3.295
consuming with all the algorithms. In fact, after we loaded the prior results at the beginning of the retrieval process, we did not need to repeat this process and the results were simply called throughout the whole process. Thus, we present the processing times for 100 frames of phase containing different pixels without considering the operations needed for loading prior results and these results are shown in Table 4. Compared with Algorithm A, Algorithm D obtained an improvement
Fig. 9. Spectra. (a) Spectrum obtained by the subtraction of H1 ′ and H2 ′ in Fig. 8(a) and (b); (b) spectrum obtained by the subtraction of H1 and H2 in Fig. 8(c) and (d).
Fig. 10. Experimental results obtained for the phase plate using different algorithms. (a) Phase retrieved by Algorithm A; (b) phase retrieved by Algorithm B; (c) phase retrieved by Algorithm C; (d) phase retrieved by Algorithm D; and (e) 1-D profiles along the dashed lines marked in (a)–(d).
16
Z. Zhong et al.
Optics and Lasers in Engineering 97 (2017) 9–18
Fig. 11. Experimental results obtained for water drops using different algorithms. (a) Phase retrieved by Algorithm A; (b) phase retrieved by Algorithm B; (c) phase retrieved by Algorithm C; (d) phase retrieved by algorithm D; and (e) 1-D profiles along the dashed lines marked in (a)–(d).
References
over two times in the terms of the overall processing time. For the holograms containing 1024 × 1024 pixels, Algorithm D achieved high frame rate processing (typically 30 fps). However, for holograms containing 768 × 768 pixels or less, all of the new algorithms can achieved better frame rate processing.
[1] Józwik M, Kozacki T, Liżewski K, Kostencka J. Digital holography with multidirectional illumination by LCoS SLM for topography measurement of high gradient reflective microstructures. Appl Opt 2015;54(9):2283–8. [2] Ma L, Rajshekhar G, Wang R, Bhaduri B, Sridharan S, Mir M, Chakraborty A, Iyer R, Prasanth S, Millet L, GilletteMU,PopescuG. Phase correlation imaging of unlabeled cell dynamics. Sci Rep 2016;6:32702. [3] Srivastava V, Anna T, Mehta DS. Full-field Hilbert phase microscopy using nearly common-path low coherence off-axis interferometry for quantitative imaging of biological cells. J Opt 2012;14:125707. [4] Bitou Y. Digital phase-shifting interferometer with an electrically addressed liquid-crystal spatial light modulator. Opt Lett 2003;28(17):1576–8. [5] Guo R, Yao B, Gao P, Min J, Zheng J, Ye T. Reflective point-diffraction microscopic interferometer with long-term stability. Chin Opt Lett 2011;9(12):120002. [6] Ghosh N, Sarka S, Bhattacharya K. Low-magnification polarization phase-shifting interference microscope for three-dimensional profilometry. Opt Eng 2012;51(8):085601. [7] Zhong Z, Hao B, Shan M, Wang Y, Diao M, Zhang Y. Two-shot common-path phase-shifting interferometer with a four-step algorithm and an unknown phase shift. Appl Opt 2014;53(10):2067–72. [8] Rodriguez-Zurita G, Toto-Arellano NI, Meneses-Fabian C, Vázquez-Castillo JF. One-shot phase-shifting interferometry: five, seven, and nine interferograms. Opt Lett 2008;33(23):2788–90. [9] Takeda M, Ina H, Kobayashi S. Fourier-transform method of fringe-pattern analysis for computer-based topography and interferometry. J Opt Soc Am 1982;72(1):156–60. [10] Pham HV, Edwards C, Goddard LL, Popescu G. Fast phase reconstruction in white light diffraction phase microscopy. Appl Opt 2013;52(1):A97–9. [11] Girshovitz P, Shaked NT. Real-time quantitative phase reconstruction in off-axis digital holography using multiplexing. Opt Lett 2014;39(8):2262–5. [12] Sha B, Liu X, Ge X, Guo C. Fast reconstruction of off-axis digital holograms based on digital spatial multiplexing. Opt Express 2014;22(19):23067–72. [13] Girshovitz P, Shaked NT. Fast phase processing in off-axis holography using multiplexing with complex encoding and live-cell fluctuation map calculation in real-time. Opt Express 2015;23(7):8773–7. [14] Das B, Yelleswarapu CS, Rao DVGLN. Parallel quadrature phase-shifting digital holographic microscopy using polarization beam splitter. Opt Commun 2012;285:4954–60.
5. Conclusion In this study, we developed three algorithms to facilitate rapid phase retrieval in slightly off-axis DH, in where we combined spectrum cropping, spatial multiplexing, and complex encoding algorithms with our previously proposed method. For phase variance in a specimen less than 2𝜋, all the new algorithms achieved faster frame rate processing with holograms containing 768 × 768 pixels or less, and the algorithm using complex encoding also improved the frame rate processing even when the holograms contained one megapixel. All of the new algorithms can eliminate the effects of the disturbing phase. These algorithms can also be implemented in parallel using a GPU to obtain further performance improvements. We hope that our approaches will be proved useful in applications to slightly-off-axis DH.
Acknowledgments The authors thank the National Natural Science Foundation of China (61377009); Major National Scientific Instrument and Equipment Development Project of China (2013YQ290489, 2011YQ040136); China Scholarship (201506685053); Heilongjiang Science Foundation of China (F201411); Harbin Science Foundation of China (2014RFQXJ030); and Fundamental Research Funds for the Central Universities. 17
Z. Zhong et al.
Optics and Lasers in Engineering 97 (2017) 9–18
[15] Shaked NT, Zhu Y, Rinehart MT, Wax A. Two-step only phase-shifting interferometry with optimized detector bandwidth for microscopy of live cells. Opt Express 2009;17(18):15585–91. [16] Gao P, Yao B, Harder I, Min J, Guo R, Zheng J,YeT. Parallel two-step phase-shifting digital holograph microscopy based on a grating pair. J Opt Soc Am A 2011;28(3):434–40. [17] Shan M, Hao B, Zhong Z, Diao M, Zhang Y. Parallel two-step spatial carrier phase-shifting common-path interferometer with a Ronchi grating outside the Fourier plane. Opt Express 2013;21(2):2126–32. [18] Hao B, Shan M, Zhong Z, Diao M, Wang Y, Zhang Y. Parallel two-step spatial carrier phase-shifting interferometric phase microscopy with fast phase retrieval. J Opt 2015;17:035602. [19] Bai H, Zhong Z, Shan M, Liu L, Guo L, Zhang Y. Interferometric phase microscopy using slightly-off-axis reflective point diffraction interferometer. Opt Lasers Eng 2017;90:155–60. [20] Pham H, Ding H, Sobh N, Do M, Patel S, Popescu G. Off-axis quantitative phase imaging processing using CUDA: toward real-time applications. Biomed Opt Express 2011;2(7):1781–93.
[21] Backoach O, Kariv S, Girshovitz P, Shaked NT. Fast phase processing in off-axis holography by CUDA including parallel phase unwrapping. Opt Express 2016;24(4):3177–88. [22] Fan Q, Yang H, Li G, Zhao J. Suppressing carrier removal error in the Fourier transform method for interferogram analysis. J Opt 2010;12:115401. [23] Zhang Q, Wu Z. A carrier removal method in Fourier transform profilometry with Zernike polynomials. Opt Lasers Eng 2013;51:253–60. [24] Du Y, Feng G, Li H, Zhou S. Accurate carrier-removal technique based on zero padding in Fourier transform method for carrier interferogram analysis. Optik 2014;125:1056–61. [25] Lan B, Feng G, Dong Z, Zhang T, Zhou S. A carrier removal method based on frequency domain self-filtering for interferogram analysis. Optik 2016;127:5961–7. [26] Shan M, Kandel ME, Majeed H, Nastata V, Popescu G. White-light diffraction phase microscopy at doubled space-bandwidth product. Opt Express 2016;24(25):29033–9.
18
Contents lists available at ScienceDirect
Optics and Lasers in Engineering journal homepage: www.elsevier.com/locate/optlaseng
Fast phase retrieval in slightly off-axis digital holography Zhi Zhong, Hongyi Bai, Mingguang Shan∗, Yabin Zhang, Lili Guo College of Information and Communication Engineering, Harbin Engineering University, Harbin, Heilongjiang 150001, PR China
a r t i c l e
i n f o
Keywords: Complex encoding Digital holography Fast phase retrieval Slightly off-axis Spatial multiplexing Spectral cropping
a b s t r a c t In this study, three efficient algorithms are proposed for fast phase retrieval in slightly off-axis digital holography using spectrum cropping, spatial multiplexing, and complex encoding. In the first algorithm, the real spectral order of the subtracted hologram is filtered and cropped, and the number of pixels is decreased in the subsequent retrieval operations. In the second algorithm, two sequential subtracted holograms are digitally phase shifted and spatial multiplexed into one synthetic hologram, and thus only one inverse Fourier transformation is then required. In the third algorithm, two sequential subtracted holograms are encoded separately into the real part and the imaginary part of a complex hologram. Two cross-correlations can be used to reconstruct the phase, thereby improving the utilization of the spectrum. The three new algorithms speed up our previously proposed retrieval method with the assistance of specimen-free holograms. Our experiments demonstrated the validity and improved time requirements of the proposed methods. © 2017 Elsevier Ltd. All rights reserved.
1. Introduction Digital holography (DH) has emerged as an important tool for achieving quantitative phase analyses of biological specimens, surface measurements, and micro-structures [1–4]. Due to the reliance of DH on interference, three types of methods in DH are usually used: on-axis, off-axis and slightly off-axis. The on-axis methods can fully utilize the space–bandwidth product of the DH system. However, to achieve high phase retrieval accuracy, phase-shifting methods must record several holograms (typically ≥3) in time sequence or in a single shot [5–8]. Unlike phase-shifting methods, the off-axis methods just need one hologram in a single shot, which make them suitable for faster imaging by fully exploiting the field of view (FOV) of the camera. However, these benefits come at the cost of the space–bandwidth product [9–13]. By integrating off-axis and phase-shifting to record two phase-shifted slightly off-axis holograms, the slightly off-axis methods provide an intermediate solution between the on-axis and off-axis methods [14–19]. After applying a subtraction operation to the two slightly off-axis holograms to eliminate the DC term, digital phase reconstruction can also be performed using a Fourier transform algorithm with the slightly off-axis DHs. However, similar to the off-axis DHs, digital processing requires large amounts of computations so it must be conducted offline. An improved phase retrieval method was proposed by our group based on the offline prior acquisition of specimen-free holograms [18,19], but it is still necessary to speed up the phase retrieval process for slightly offaxis DHs. The retrieval operations are usually the same for different pixels in one hologram. Parallel computing techniques using graphic processing ∗
Corresponding author. E-mail address: [email protected] (M. Shan).
http://dx.doi.org/10.1016/j.optlaseng.2017.05.004 Received 19 December 2016; Received in revised form 2 May 2017; Accepted 2 May 2017 0143-8166/© 2017 Elsevier Ltd. All rights reserved.
units (GPUs) are then introduced to achieve higher retrieval efficiency [20,21]. However, special graphic units and programming skills are required. Considering that Fourier transform operation is the main timeconsuming process, Girshovitz and Shaked retained only the real spectral order and omit all the other redundant data during the retrieval process [11]. They then achieved real time phase reconstruction with 1 megapixel. To further speed up this process, Sha et al. spatially multiplexed four holograms into one synthetic hologram [12], while Girshovitz and Shaked encoded four sequential holograms into one complex hologram, and then used both the real and twin image parts [13]. However, all of these methods are only suitable for off-axis DHs. In this study, we present three more efficient algorithms for phase retrieval in slightly off-axis DHs. Those algorithms combine spectrum cropping, spatial multiplexing, and complex encoding with our previously proposed method. They speed up phase retrieval in slightly offaxis DHs, as well as compensating for the disturbing phase caused by the phase shift 𝛼, carrier frequency, and background phase. Finally, we demonstrated the improved performance of our proposed approaches based on simulations, and by processing phase maps of a phase plate and water drops. 2. Reconstruction algorithms In slightly off-axis DH, two measured holograms with a phase shift 𝛼 have the forms: [( )] 𝐻1 (𝑥, 𝑦) = 𝑎(𝑥, 𝑦) + 0.5𝑏(𝑥, 𝑦) exp 𝑖 𝜑𝑡 + 𝜑𝑏𝑔 +2𝜋𝑓𝑥 𝑥 [ ( )] + 0.5𝑏(𝑥, 𝑦) exp −𝑖 𝜑𝑡 + 𝜑𝑏𝑔 +2𝜋𝑓𝑥 𝑥 ,
(1)
Z. Zhong et al.
Optics and Lasers in Engineering 97 (2017) 9–18
Fig. 1. Algorithm A: the algorithm based on an offline prior acquisition of specimen-free holograms.
[( )] 𝐻2 (𝑥, 𝑦) = 𝑎(𝑥, 𝑦) + 0.5𝑏(𝑥, 𝑦) exp 𝑖 𝜑𝑡 + 𝜑𝑏𝑔 +2𝜋𝑓𝑥 𝑥 exp (𝑖𝛼) [ ( )] + 0.5𝑏(𝑥, 𝑦) exp −𝑖 𝜑𝑡 + 𝜑𝑏𝑔 +2𝜋𝑓𝑥 𝑥 exp (−𝑖𝛼),
is subtracted out, so the maximum circle in the BPF can reach the zerofrequency point [14]. Thus, for a spectral matrix containing N × N pixels, the maximum diameter of the circle in the BPF can be selected as N/2. A4: Two-dimensional inverse FFT (2-D IFFT): Convert the result obtained in step A3 back to the image domain by using a 2-D IFFT to product an N × N complex matrix containing the phase distribution of the specimen and disturbance. The disturbing phase includes phase shift, carrier information and background phase. A5: Load the result of prior measurement: To compensate for the disturbing phase, obtain two phase shifted holograms without a specimen before the experiment and then process them according to steps A1–A4, as described above. The result containing N × N pixels is actually stored before the measurement and it is only loaded to help retrieving phase during the reconstruction process. A6: Specimen phase retrieval: Divide the result obtained from step A4 by that from step A5 to eliminate the disturbing phase. Finally, implement an arc tangent operation to obtain the phase distribution of the specimen.
(2)
where a(x, y) is the DC term, b(x, y) is the modulation term, 𝜑t is the phase distribution of specimen at time t, 𝜑bg is the background phase caused by the aberration and the noise of the system, fx is carrier frequency in the x direction and 2𝜋fx x is its corresponding phase tilt. After applying subtraction, we can eliminate the DC term and obtain: {[ 𝑆 = 𝐻1 − 𝐻2 = 0.5𝑏(𝑥, 𝑦)
] [( )] } 1 − exp (𝑖𝛼) exp 𝑖 𝜑𝑡 + 𝜑𝑏𝑔 +2𝜋𝑓𝑥 𝑥 [ ] [ ( )] . + 1 − exp (−𝑖𝛼) exp −𝑖 𝜑𝑡 + 𝜑𝑏𝑔 +2𝜋𝑓𝑥 𝑥 (3)
The phase distribution 𝜑t of a specimen can be retrieved from the subtracted hologram using Fourier transform algorithm [17]. To increase the retrieval speed, we suggested a faster retrieval algorithm based on the offline prior acquisition of specimen-free subtracted holograms [18,19]. The faster algorithm can also compensate for the disturbing phase caused by the phase shift 𝛼, carrier frequency, and background phase. In the following, we first review the algorithm, and then propose three more efficient algorithms.
2.2. Algorithm B: spectrum cropping algorithm Shaked et al. noted that when selecting the real spectral order, there is no need to use a BPF matrix with the same size as the original images (N × N). The filtering operation can be performed by decreasing the number of pixels and the 2-D FFT operation can be followed by using a cropped matrix with the same size. For slightly off-axis DH, we can crop the real spectral order using a BPF with N/2 × N/2 pixels and continue processing with a matrix that is four times smaller. Fig. 2 shows the digital process, which is described in Algorithm B. The algorithm comprises the following step. B1: Holograms subtraction: This is the same as step A1. B2: 2-D FFT: This is the same as step A2. B3: DBPF and cropping: Apply a BPF containing N/2 × N/2 pixels to the area (N/2 + 1:N, N/4 + 1:3 N/4) of the spectrum and crop the real spectral order. B4: 2-D IFFT: Convert the result obtained from step B3 back into the image domain by using a 2-D IFFT to produce an N/2 × N/2 complex matrix containing the phase distribution of the specimen and disturbance.
2.1. Algorithm A: algorithm based on the offline prior acquisition of specimen-free holograms When the phase variance of a specimen is less than 2𝜋, we previously proposed a simple retrieval algorithm assisted by a prior measurement without any specimen [18,19]. This algorithm can achieve a fast phase retrieval speed and it comprises the following steps, as shown in Fig. 1. A1: Holograms subtraction: Subtract the two phase-shifted holograms containing N × N real pixels to eliminate the DC term and obtain a subtracted hologram. A2: Two-dimensional fast Fourier transform (2-D FFT): Convert the subtracted hologram into the spatial-frequency domain using a 2-D FFT, thereby obtaining a matrix containing N × N complex pixels. A3: Digital band pass filtering (DBPF): Use a band pass filter (BPF) with a size of N × N pixels and digitally filter out the twin spectral order, but retain the real spectral order. In slightly off-axis DH, the DC order 10
Z. Zhong et al.
Optics and Lasers in Engineering 97 (2017) 9–18
Fig. 2. Algorithm B: the spectrum cropping algorithm.
B5: Load the result of prior measurement: Acquire two specimenfree phase-shifted holograms before the experiments and process them according to steps B1–B4 to obtain a matrix containing N/2 × N/2 pixels. Load the result obtained previously during the reconstruction process. B6: Specimen phase retrieval: Divide the two results obtained from steps B4 and B5, and apply an arc tangent operation to obtain the phase distribution of the specimen containing N/2 × N/2 pixels. Finally, enlarge the N/2 × N/2 phase to the size of N × N phase matrix. 2.3. Algorithm C: spatial multiplexing algorithm
(7)
( ) { { } } 𝑐2 ′ 𝑥′ , 𝑦′ = IFT FT 𝑆𝑆 ⋅ BPF2 ( )[ ] = 0.5𝑏′ 𝑥′ , 𝑦′ 1 − exp (𝑖𝛼) [( ′ )] exp 𝑖 𝜑𝑏𝑔 +2𝜋𝑓𝑥 ′ 𝑥′ + 2𝜋𝑓𝑦𝑚 ′ 𝑦′ .
(8)
Thus, the phase retrieved can be expressed as: [ ( )] 𝑐1 𝑥′ , 𝑦′ ( ) ′ ′ ′ 𝜑𝑡 𝑥 , 𝑦 = arctan . 𝑐 1 ′ (𝑥 ′ , 𝑦 ′ ) [ ( )] 𝑐2 𝑥′ , 𝑦′ ( ) 𝜑𝑡+1 ′ 𝑥′ , 𝑦′ = arctan . 𝑐 2 ′ (𝑥 ′ , 𝑦 ′ )
According to the shift theorem for the Fourier transform, a phase shift in the spatial domain can be interpreted as a frequency shift in the frequency domain. Thus, two subtracted holograms can be digitally phase shifted in the y direction and spatially multiplexed into one image [12] as follows: ( ) ( ) 𝑆𝑆 = 𝑆𝑡 exp −𝑖2𝜋𝑓𝑦𝑚 𝑦 + 𝑆𝑡+1 exp 𝑖2𝜋𝑓𝑦𝑚 𝑦 (4)
(9)
(10)
Finally, 𝜑t and 𝜑t+1 containing N × N pixels can be acquired by enlarging 𝜑t ′ and 𝜑t+1 ′ containing N/2 × N/2 pixels. A flowchart illustrating the algorithm described above is shown in Fig. 3 and its steps are as follows. C1: Holograms subtraction: This is the same as step A1 but for two sequential subtracted holograms St and St +1 . C2: Synthesis of subtracted holograms: Digitally phase shift St and St +1 , and sum them to obtain one synthetic hologram SS . C3: 2-D FFT: Convert SS into the spatial-frequency domain using a 2D FFT and obtain a matrix containing N × N complex pixels. This step can actually be replaced by two steps to yield a 1.5-D FFT algorithm where a 1-D FFT in rows is first implemented to produce N × N pixels, and a 1-D FFT in columns is calculated to obtain an N × N/2 array containing only the real images from the two subtractions [11]. To facilitate convenient comparisons with other methods, we use a 2-D FFT in this study. C4: DBPF and cropping: Apply two BPFs containing N/2 × N/2 pixels to cover the areas (N/2 + 1:N, 1:N/2) and (N/2 + 1:N, N/2 + 1: N) of the spectrum, and crop the corresponding real spectral orders of 𝜑t and 𝜑t +1 . C5: 2-D IFFT: Convert the result obtained from step C4 back into the image domain by using a 2-D IFFT to yield two N/2 × N/2 complex matrices containing the phase distributions of 𝜑t and 𝜑t +1 . C6: Load the results of prior measurements: This is the same as step B5 but for two sequential phase distributions each containing N/2 × N/2 pixels.
where St and St +1 are two sequential subtracted holograms and fym is a spatial phase shift in the y direction. For a typical slightly off-axis DH, the bandwidth B of the subtracted hologram is equal to two times that of the recorded specimen. Thus, in order to separate the spectra of St and St +1 , fym should not be less than B/4. According to Eq. (3) and Algorithm B, the cropped real spectral orders of St and St +1 can be extracted as: ( ) { { } } 𝑐1 𝑥′ , 𝑦′ = IFT FT 𝑆𝑆 ⋅ BPF1 ( )[ ] = 0.5𝑏′ 𝑥′ , 𝑦′ 1 − exp (𝑖𝛼) [( ′ )] (5) exp 𝑖 𝜑𝑡 + 𝜑𝑏𝑔 ′ +2𝜋𝑓𝑥 ′ 𝑥′ − 2𝜋𝑓𝑦𝑚 ′ 𝑦′ . ( ) { { } } 𝑐2 𝑥′ , 𝑦′ = IFT FT 𝑆𝑆 ⋅ BPF2 ( )[ ] = 0.5𝑏′ 𝑥′ , 𝑦′ 1 − exp (𝑖𝛼) [( )] exp 𝑖 𝜑𝑡+1 ′ + 𝜑𝑏𝑔 ′ +2𝜋𝑓𝑥 ′ 𝑥′ + 2𝜋𝑓𝑦𝑚 ′ 𝑦′ .
( )[ ] = 0.5𝑏′ 𝑥′ , 𝑦′ 1 − exp (𝑖𝛼) [( ′ )] exp 𝑖 𝜑𝑏𝑔 +2𝜋𝑓𝑥 ′ 𝑥′ − 2𝜋𝑓𝑦𝑚 ′ 𝑦′ .
(6)
The corresponding cropped real spectral orders of free-specimens can also be extracted as: ( ) { { } } 𝑐1 ′ 𝑥′ , 𝑦′ = IFT FT 𝑆𝑆 ⋅ BPF1 11
Z. Zhong et al.
Optics and Lasers in Engineering 97 (2017) 9–18
Fig. 3. Algorithm C: the spatial multiplexing algorithm.
C7: Specimen phase retrieval: This is the same as step B6 but for two sequential phase distributions.
( ) { { } } 𝑑2 𝑥′ , 𝑦′ = IFT FT 𝑆𝐶 ⋅ BPF2 {[ } ] [ ( )] ′ ( ) 1 − exp (−𝑖𝛼) exp −[𝑖 𝜑(𝑡 ′ + 𝜑𝑏𝑔 ′ +2𝜋𝑓𝑥 𝑥′ )] . [ ] = 0.5𝑏′ 𝑥′ , 𝑦′ + 1 − exp (−𝑖𝛼) exp −𝑖 𝜑𝑡+1 ′ + 𝜑𝑏𝑔 ′+2𝜋𝑓𝑥 ′ 𝑥′ − 𝜋2
2.4. Algorithm D: complex encoding algorithm Like off-axis holograms, the slightly off-axis holograms also have no imaginary parts, which make the Fourier transformation operation inefficient. Thus, it can employ the imaginary part to encode two sequential subtracted holograms into a complex hologram [13]. The complex hologram is encoded as follows. 𝑆𝐶 = 𝑆𝑡 + 𝑖𝑆𝑡+1
(13) After applying the conjugate operation to Eq. (13), we can obtain the following results by the summation and subtraction of Eqs. (12) and (13):
(11)
( )[ ] [( )] ( ) ( ) 𝑏′ 𝑥′ , 𝑦′ 1 − exp (𝑖𝛼) exp 𝑖 𝜑𝑡 ′ + 𝜑𝑏𝑔 ′ +2𝜋𝑓𝑥 ′ 𝑥′ = 𝑑1 𝑥′ , 𝑦′ +𝑑2∗ 𝑥′ , 𝑦′ ,
After a 2-D FFT operation, the spectrum can be obtained for SC . The real and imaginary orders containing N/2 × N/2 pixels can then be obtained by BPF and cropping, and the 2-D IFFT calculation. This process can be expressed as:
(14)
( ) { { } } 𝑑1 𝑥′ , 𝑦′ = IFT FT 𝑆𝐶 ⋅ BPF1 {[ } ] [( )] ′ ( ) 1 − exp (𝑖𝛼) exp 𝑖[𝜑(𝑡 ′ + 𝜑𝑏𝑔 ′ +2𝜋𝑓𝑥 𝑥′ )] ′ ′ ′ [ ] . = 0.5𝑏 𝑥 , 𝑦 + 1 − exp (𝑖𝛼) exp 𝑖 𝜑𝑡+1 ′ + 𝜑𝑏𝑔 ′ +2𝜋𝑓𝑥 ′ 𝑥′ + 𝜋2
( ) ( ) ( )[ ] [( )] 𝑑1 𝑥′ , 𝑦′ −𝑑2∗ 𝑥′ , 𝑦′ 𝑏′ 𝑥′ , 𝑦′ 1 − exp (𝑖𝛼) exp 𝑖 𝜑𝑡+1 ′ + 𝜑𝑏𝑔 ′+2𝜋𝑓𝑥 ′ 𝑥′ = , 𝑖 (15) where, “∗ ” denotes the conjugate operation.
(12) 12
Z. Zhong et al.
Optics and Lasers in Engineering 97 (2017) 9–18
Fig. 4. Algorithm D: the complex encoding algorithm.
Repeating the operation above, we can obtain the same specimenfree results as follows: ( )[ ] [( )] ( ) ( ) 𝑏′ 𝑥′ , 𝑦′ 1 − exp (𝑖𝛼) exp 𝑖 𝜑𝑏𝑔 ′ +2𝜋𝑓𝑥 ′ 𝑥′ = 𝑑1 ′ 𝑥′ , 𝑦′ + 𝑑2 ′∗ 𝑥′ , 𝑦′ ,
D1: Holograms subtraction: This is the same as step C1. D2: Encoding subtracted holograms: According to Eq. (11), sum the two subtracted holograms St and St +1 to obtain a complex hologram SC , where St is the real part of SC and St +1 is the imaginary part of SC . The complex hologram SC contains N × N pixels. D3: 2-D FFT: Convert SC into the spatial-frequency domain using a 2-D FFT to obtain a matrix containing N × N complex pixels. D4: DBPF and cropping: Apply two BPFs containing N/2 × N/2 pixels to cover the area (1:N/2, N/4 + 1:3 N/4) and (N/2 + 1:N, N/4 + 1:3 N/4) of the spectrum, and crop the two cross-correlations. D5: 2-D IFFT: Convert the results obtained from step D4 back into the image domain using 2-D IFFTs, thereby producing two N/2 × N/2 complex matrices as shown in Eqs. (12) and (13). D6: Load the results of prior measurements: Load the specimenfree results obtained according to steps D1–D5 following Eqs. (16) and (17) to produce two N/2 × N/2 complex matrices. D7: Specimen phase retrieval: Process the results obtained above using Eqs. (18) and (19) to produce two sequential phase distributions each containing N/2 × N/2 pixels. Finally, enlarge the phase to the N × N phase matrix.
(16) ( ) ( ) ( )[ ] [( )] 𝑑1 ′ 𝑥′ , 𝑦′ − 𝑑2 ′∗ 𝑥′ , 𝑦′ 𝑏′ 𝑥′ , 𝑦′ 1 − exp (𝑖𝛼) exp 𝑖 𝜑𝑏𝑔 ′ +2𝜋𝑓𝑥 ′ 𝑥′ = 𝑖 (17) Thus, the retrieved phases of 𝜑t ′ and 𝜑t +1 ′ containing N/2 × N/2 pixels can be expressed as: [ ( ′ ′) ( )] 𝑑1 𝑥 , 𝑦 + 𝑑2∗ 𝑥′ , 𝑦′ ′ 𝜑𝑡 = arctan . (18) 𝑑 1 ′ (𝑥 ′ , 𝑦 ′ ) + 𝑑 2 ′ ∗ (𝑥 ′ , 𝑦 ′ ) [ ( ) ( )] 𝑑1 𝑥′ , 𝑦′ − 𝑑2∗ 𝑥′ , 𝑦′ 𝜑′𝑡+1 = arctan 𝑖 ′ . (19) 𝑑1 (𝑥′ , 𝑦′ ) + 𝑑2 ′ ∗ (𝑥′ , 𝑦′ ) After an enlarging operation, we can finally obtain 𝜑t and 𝜑t +1 containing N × N pixels. A flowchart illustrating the algorithm described above is shown in Fig. 4 and its steps are as follows: 13
Z. Zhong et al.
Optics and Lasers in Engineering 97 (2017) 9–18
Fig. 5. Simulations. (a) Specimen phase; (b) specimen hologram with 𝛼 = 0; (c) specimen hologram with 𝛼 = 𝜋/2; (d) specimen-free hologram with 𝛼 = 0; and (e) specimen-free hologram with 𝛼 = 𝜋/2.
of 0° and 90°, phase shifts of 0 and 𝜋/2 were introduced into the holograms, respectively. Compared with the method employed by Ref. [19], the proposed setup could avoid the effects of the intensity ratio between different components of the object and reference beams, and as well as making full use of the FOV of the camera. In the experiments, a He–Ne laser with a wavelength of 𝜆 = 632.8 nm was used as the light source. Holograms were recorded by a CCD camera with 1600 × 1200 pixels (each pixel size is 4.4 μm × 4.4 μm). Fig. 8(a) and (b) show the holograms captured without specimens, and Fig. 8(c) and (d) show the holograms captured with a phase plate (refractive index n = 1.5168) with a step height of 580.22 nm measured by a BRUKER atomic force microscope, which could provide an optical path difference (OPD) of 2.98 rad at 𝜆 = 632.8 nm. The four holograms all comprised 640 × 640 pixels, and the spectra of H1 ′–H2 ′ and H1 –H2 are shown in Fig. 9. The results indicate that the DC term was eliminated, which increased the frequency coverage of the real spectral order. We used a BPF with 640 × 640 pixels for Algorithm A, and BPFs with 320 × 320 pixels for Algorithms B–D. However, to achieve a trade-off between resolution and spatial phase noise [26], we selected the pass area of the BPF using F = 0.6, where F = 1.0 corresponds to the diameter of the circle in Fig. 9 that touches the zero-frequency point. Fig. 10 shows the similar results retrieved using all four algorithms. To further confirm the validity of the algorithms, multiple water drops placed on part of a slide were also measured, where Fig. 11 shows the corresponding similar results retrieved using all four algorithms. To evaluate the run time of the different algorithms, holograms with different pixel sizes were retrieved with Matlab 2014a using a simple desktop (Intel i5-4430, 3 GHz, 8 GB RAM, where only a single core was used). After a large number of repetitions, the average processing times required to reconstruct one frame of phase containing 1024 × 1024 pixels using Algorithms A and B are shown in Table 2, and the times for two frames using Algorithms C and D are shown in Table 3. These results demonstrate that compared with Algorithm A, all the new algorithms could effectively speed up the retrieval process. However, it should be noted that the operations required to load prior results were very time-
Table 1 Error comparison for different algorithms.
ARE (%) RMSE (mrad)
Algorithm A
Algorithm B
Algorithm C
Algorithm D
0.0023 0.123
0.11 2.3
0.11 2.3
0.11 2.3
3. Simulations The proposed algorithms were verified based on numerical simulations. Fig. 5(a) shows a specimen phase of <2𝜋 rad with 1024 × 1024 pixels (each pixel size is 4.4 μm × 4.4 μm). For a slightly off-axis DH system [15], the carrier frequency fx was set to be 256/1024/(4.4 μm) [22]. The simulated specimen holograms that satisfied Eqs. (1) and (2) are shown in Fig. 5(b) and (c) with 𝛼 = 𝜋/2, respectively, and the corresponding specimen-free holograms are shown in Fig. 5(d) and (e). The phase distribution was then retrieved using Algorithms A–D and the results are shown in Fig. 6. It can be seen that despite spectrum cropping, the average relative error (ARE) and root mean square error (RMSE) [23] with the proposed algorithms, as shown in Table 1, were higher compared with those using Algorithm A, and all the proposed algorithms retrieved the specimen phase with high quality. In fact, ARE and RMSE results obtained using the proposed algorithms were far lower than those using the carrier estimation algorithms without the aid of specimen-free holograms [24,25]. 4. Experiment results In order to further demonstrate the efficiency of the algorithms, experiments were conducted with slightly off-axis DH microscopy, as shown in Fig. 7. Similar to the method employed in our previous work [19], the microscopy approach utilized a 45° linearly polarized object beam and a circularly polarized reference beam for interference. However, a polarizer P3 was placed before the CCD camera to replace the polarizing beam splitter used by Ref. [19]. By setting P3 with angles 14
Z. Zhong et al.
Optics and Lasers in Engineering 97 (2017) 9–18
Fig. 6. Simulation results. (a) Phase retrieved by Algorithm A; (b) phase retrieved by Algorithm B; (c) phase retrieved by Algorithm C; (d) phase retrieved by algorithm D; and (e) 1-D profiles along the dashed lines marked in Figs. 5(c) and 6(a)–(d) (all the color bars indicate rad). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 7. Experimental setup: P1, P2 and P3, polarizers; QW, quarter-wave plate; MO, microscopic objective; S, specimen; CL, collimated lens; L1 and L2, lenses with a same focal length; BS, beam-splitter; M1, reflective pinhole mirror; and RR, retro-reflector.
Table 2 Processing time for one frame of phase containing 1024 × 1024 pixels using Algorithms A and B (ms).
Algorithm A Algorithm B
2D-FFT
2D-IFFT
Load the result of prior measurement
Specimen phase retrieval
Others
Total time
24.347 24.236
24.068 7.839
94.754 24.639
16.471 11.521
13.409 7.359
173.049 75.594
15
Z. Zhong et al.
Optics and Lasers in Engineering 97 (2017) 9–18
Fig. 8. Holograms. (a, b) Holograms captured without specimens; (c, d) holograms captured with a phase plate. Table 3 Processing time for two frames of phase containing 1024 × 1024 pixels using Algorithms C and D (ms).
Algorithm C Algorithm D
Subtracted-holograms synthesizing
2D-FFT
2D-IFFT
Load the result of prior measurement
Specimen phase retrieval
Others
Total time
24.563 4.906
22.862 22.749
12.826 12.845
49.451 49.434
21.371 22.463
6.672 7.034
137.745 119.431
Table 4 Processing time for 100 frames of phase containing different pixels (s). Pixel size
Algorithm A
Algorithm B
Algorithm C
Algorithm D
512 × 512 768 × 768 1024 × 1024
1.520 3.660 7.298
1.001 2.378 4.436
0.978 2.323 4.269
0.732 1.777 3.295
consuming with all the algorithms. In fact, after we loaded the prior results at the beginning of the retrieval process, we did not need to repeat this process and the results were simply called throughout the whole process. Thus, we present the processing times for 100 frames of phase containing different pixels without considering the operations needed for loading prior results and these results are shown in Table 4. Compared with Algorithm A, Algorithm D obtained an improvement
Fig. 9. Spectra. (a) Spectrum obtained by the subtraction of H1 ′ and H2 ′ in Fig. 8(a) and (b); (b) spectrum obtained by the subtraction of H1 and H2 in Fig. 8(c) and (d).
Fig. 10. Experimental results obtained for the phase plate using different algorithms. (a) Phase retrieved by Algorithm A; (b) phase retrieved by Algorithm B; (c) phase retrieved by Algorithm C; (d) phase retrieved by Algorithm D; and (e) 1-D profiles along the dashed lines marked in (a)–(d).
16
Z. Zhong et al.
Optics and Lasers in Engineering 97 (2017) 9–18
Fig. 11. Experimental results obtained for water drops using different algorithms. (a) Phase retrieved by Algorithm A; (b) phase retrieved by Algorithm B; (c) phase retrieved by Algorithm C; (d) phase retrieved by algorithm D; and (e) 1-D profiles along the dashed lines marked in (a)–(d).
References
over two times in the terms of the overall processing time. For the holograms containing 1024 × 1024 pixels, Algorithm D achieved high frame rate processing (typically 30 fps). However, for holograms containing 768 × 768 pixels or less, all of the new algorithms can achieved better frame rate processing.
[1] Józwik M, Kozacki T, Liżewski K, Kostencka J. Digital holography with multidirectional illumination by LCoS SLM for topography measurement of high gradient reflective microstructures. Appl Opt 2015;54(9):2283–8. [2] Ma L, Rajshekhar G, Wang R, Bhaduri B, Sridharan S, Mir M, Chakraborty A, Iyer R, Prasanth S, Millet L, GilletteMU,PopescuG. Phase correlation imaging of unlabeled cell dynamics. Sci Rep 2016;6:32702. [3] Srivastava V, Anna T, Mehta DS. Full-field Hilbert phase microscopy using nearly common-path low coherence off-axis interferometry for quantitative imaging of biological cells. J Opt 2012;14:125707. [4] Bitou Y. Digital phase-shifting interferometer with an electrically addressed liquid-crystal spatial light modulator. Opt Lett 2003;28(17):1576–8. [5] Guo R, Yao B, Gao P, Min J, Zheng J, Ye T. Reflective point-diffraction microscopic interferometer with long-term stability. Chin Opt Lett 2011;9(12):120002. [6] Ghosh N, Sarka S, Bhattacharya K. Low-magnification polarization phase-shifting interference microscope for three-dimensional profilometry. Opt Eng 2012;51(8):085601. [7] Zhong Z, Hao B, Shan M, Wang Y, Diao M, Zhang Y. Two-shot common-path phase-shifting interferometer with a four-step algorithm and an unknown phase shift. Appl Opt 2014;53(10):2067–72. [8] Rodriguez-Zurita G, Toto-Arellano NI, Meneses-Fabian C, Vázquez-Castillo JF. One-shot phase-shifting interferometry: five, seven, and nine interferograms. Opt Lett 2008;33(23):2788–90. [9] Takeda M, Ina H, Kobayashi S. Fourier-transform method of fringe-pattern analysis for computer-based topography and interferometry. J Opt Soc Am 1982;72(1):156–60. [10] Pham HV, Edwards C, Goddard LL, Popescu G. Fast phase reconstruction in white light diffraction phase microscopy. Appl Opt 2013;52(1):A97–9. [11] Girshovitz P, Shaked NT. Real-time quantitative phase reconstruction in off-axis digital holography using multiplexing. Opt Lett 2014;39(8):2262–5. [12] Sha B, Liu X, Ge X, Guo C. Fast reconstruction of off-axis digital holograms based on digital spatial multiplexing. Opt Express 2014;22(19):23067–72. [13] Girshovitz P, Shaked NT. Fast phase processing in off-axis holography using multiplexing with complex encoding and live-cell fluctuation map calculation in real-time. Opt Express 2015;23(7):8773–7. [14] Das B, Yelleswarapu CS, Rao DVGLN. Parallel quadrature phase-shifting digital holographic microscopy using polarization beam splitter. Opt Commun 2012;285:4954–60.
5. Conclusion In this study, we developed three algorithms to facilitate rapid phase retrieval in slightly off-axis DH, in where we combined spectrum cropping, spatial multiplexing, and complex encoding algorithms with our previously proposed method. For phase variance in a specimen less than 2𝜋, all the new algorithms achieved faster frame rate processing with holograms containing 768 × 768 pixels or less, and the algorithm using complex encoding also improved the frame rate processing even when the holograms contained one megapixel. All of the new algorithms can eliminate the effects of the disturbing phase. These algorithms can also be implemented in parallel using a GPU to obtain further performance improvements. We hope that our approaches will be proved useful in applications to slightly-off-axis DH.
Acknowledgments The authors thank the National Natural Science Foundation of China (61377009); Major National Scientific Instrument and Equipment Development Project of China (2013YQ290489, 2011YQ040136); China Scholarship (201506685053); Heilongjiang Science Foundation of China (F201411); Harbin Science Foundation of China (2014RFQXJ030); and Fundamental Research Funds for the Central Universities. 17
Z. Zhong et al.
Optics and Lasers in Engineering 97 (2017) 9–18
[15] Shaked NT, Zhu Y, Rinehart MT, Wax A. Two-step only phase-shifting interferometry with optimized detector bandwidth for microscopy of live cells. Opt Express 2009;17(18):15585–91. [16] Gao P, Yao B, Harder I, Min J, Guo R, Zheng J,YeT. Parallel two-step phase-shifting digital holograph microscopy based on a grating pair. J Opt Soc Am A 2011;28(3):434–40. [17] Shan M, Hao B, Zhong Z, Diao M, Zhang Y. Parallel two-step spatial carrier phase-shifting common-path interferometer with a Ronchi grating outside the Fourier plane. Opt Express 2013;21(2):2126–32. [18] Hao B, Shan M, Zhong Z, Diao M, Wang Y, Zhang Y. Parallel two-step spatial carrier phase-shifting interferometric phase microscopy with fast phase retrieval. J Opt 2015;17:035602. [19] Bai H, Zhong Z, Shan M, Liu L, Guo L, Zhang Y. Interferometric phase microscopy using slightly-off-axis reflective point diffraction interferometer. Opt Lasers Eng 2017;90:155–60. [20] Pham H, Ding H, Sobh N, Do M, Patel S, Popescu G. Off-axis quantitative phase imaging processing using CUDA: toward real-time applications. Biomed Opt Express 2011;2(7):1781–93.
[21] Backoach O, Kariv S, Girshovitz P, Shaked NT. Fast phase processing in off-axis holography by CUDA including parallel phase unwrapping. Opt Express 2016;24(4):3177–88. [22] Fan Q, Yang H, Li G, Zhao J. Suppressing carrier removal error in the Fourier transform method for interferogram analysis. J Opt 2010;12:115401. [23] Zhang Q, Wu Z. A carrier removal method in Fourier transform profilometry with Zernike polynomials. Opt Lasers Eng 2013;51:253–60. [24] Du Y, Feng G, Li H, Zhou S. Accurate carrier-removal technique based on zero padding in Fourier transform method for carrier interferogram analysis. Optik 2014;125:1056–61. [25] Lan B, Feng G, Dong Z, Zhang T, Zhou S. A carrier removal method based on frequency domain self-filtering for interferogram analysis. Optik 2016;127:5961–7. [26] Shan M, Kandel ME, Majeed H, Nastata V, Popescu G. White-light diffraction phase microscopy at doubled space-bandwidth product. Opt Express 2016;24(25):29033–9.
18