Quantitative phase imaging using a combination of flat fielding and windowed Fourier filtering demodulated by a graph cuts algorithm for screening opaque and transparent objects

Dahi Ghareab Abdelsalam Ibrahim

doi:10.1364/OPTCON.448128

1. Introduction

Objects viewed in highly coherent light acquire a peculiar granular appearance [1]. The detailed structure of this peculiar granular appears chaotic and unordered. When the coherent light is diffusely reflected from a rough surface with roughness on the scale of an optical wavelength, the obtained field consists of many wavelets. Interference of the coherent wavelets results in the granular intensity pattern is defined as speckles. Using mean and median filters [2] are good examples for speckle noise suppression. However, they do not work well for waving structures in fringe patterns. A good exception is the sine or cosine filter [3], but it fails to deal with fringe patterns with fast waving structures [4]. Ibrahim used convolution of windowed Fourier transform with Chebyshev type 2 and elliptic filters to improve the intensity-contrast image of a noisy digital hologram [5]. However, enhancement of the phase-contrast image was not considered. Due to the nature of optical interferometry, the extracted phase may present some ambiguities [6,7]. Such ambiguities are phase-order ambiguity and sign ambiguity. The phase-order ambiguity can be removed by phase unwrapping, while the sign ambiguity can be mitigated by forcing the derivatives of local frequency. The high-quality pixels are processed with no error, however, the low quality pixels may cause an error. In this paper, we proposed a method to increase the efficiency of low quality pixels to avoid the error in phase-contrast image resulting in the sign ambiguity. The method is a combination of flat fielding and windowed Fourier filtering demodulated by graph cuts algorithm. The flat fielding removes dust particles attached with the interferogram as well as adjusts inhomogeneity of its intensity. The corrected interferogram is then convolved by windowed Fourier filtering (WFF) to produce an interferogram free from speckle noise. The obtained interferogram is then reconstructed by the Fourier transform (FT) method, or wavelet transform (WT) method, or angular spectrum method (ASM) to extract the wrapping phase map of the object being screened. A modified graph cuts algorithm has been used to unwrap the wrapping phase map to remove the 2π ambiguity [8]. The algorithm segmented the image into background and foreground images and the procedure continuous for multiple regions simultaneously until no further segment is found. The proposed method has been applied to quantify the phase of opaque and transparent objects precisely. For the opaque objects, off-axis interferograms of large scale step height of 200µm and small scale step height of 1.50µm were captured and corrected by the proposed method then reconstructed by the FT method and the WT method. For the transparent objects, inline interferograms were captured and corrected by the proposed method then reconstructed by the ASM. The purpose for using inline interferograms in case of transparent objects is to make full use of the space-bandwidth product (SBP) of the camera and hence can resolve finer spatial details of the biological cells. Since the conventional reconstructed inline interferogram is always blurred by the transmitted zeroth-order light and twin image, the proposed method mitigates the zero-order light and twin image of the inline interferogram significantly. Experimental results show that a systematic error [9] due to phase ripples, which it depends on the random speckle phase and characterizes the accuracy of the system [10,11] have been minimized to the least measure. Such ripples are inevitable due to the side lobes effect when the FT method is used [12,13]. Quantitatively, over one hundred measurements, the mean standard deviation of the measured error due to phase ripples has been estimated to be in the range of 0.6rad, that is path length in the range of 30nm. Since the proposed method minimizes such noise to the least measure, so we claim that quantitative phase of the object being screened is measured precisely. Moreover, the lateral resolution, which is represented by slope of the roll-off is slightly improved without applications of digital filters. Based on the merits of the proposed method, we applied it to extract the phase map of the shapes of the echinocytes of the cancerous blood cells which have the sharpest spatial features. To the best of our knowledge, this is the first time that a combination of flat fielding and WFF demodulated with a graph cuts algorithm is used to precisely quantify the phase-contrast images.

2. Single shot reconstruction methods

In this paper, we used three reconstruction methods to extract the phase map of the object being screened. The three methods are the FT, the WT, and the ASM. Unlike the phase shifting (PS) method [14–16], these three methods have the ability to demodulate a single interferogram. The FT method is useful for analysis of stationary signals as seen in Fig. 1(a).

Fig. 1. (a) An example of a stationary signal. (b) An example of a non-stationary signal.

Download Full Size | PDF

The FT has poor ability to localize the frequency components of a signal with respect to position [17]. Windowed-Fourier Transform (WFT) [18,19] has been modified for demodulation to provide time-frequency representation with better signal localization. WFT is based on division of the non-stationary signal into small portions, which are assumed to be stationary. This is done using a window function of a selected width, which is shifted and multiplied with the signal to obtain the small stationary signals. The Fourier transform is then applied to each of these portions to obtain the required transform of the signal. This WFT method uses a fixed window length and therefore gives a fixed resolution at all times. The WL method [20] improves the fixed resolution in the WFT method to some extent. Moreover, it is a suitable method for demodulation of non-stationary signals as seen in Fig. 1(b) and works well with noisy interferograms with no necessity for choice of window in the transform plane. The ASM is valid for a shorter distance, so we used it to reconstruct the modulated interferograms of the transparent objects (biological blood cells). For the opaque objects, we used the FT method to reconstruct the modulated interferograms. Since the proposed method based FT may shift the roll-off of the object to the interior edge (Fig. 5(a)), therefore, we applied the WL method to spatially relocalize the roll-off of the object to the exterior edge.

2.1. Fourier transform method

The FT method is a popular method for interferogram demodulation [21]. To extract the quantitative phase map of the object, the Fourier transform algorithm, which includes a digital two-dimensional (2D) Fourier transform, filtering one of the cross-correlation terms, and an inverse 2D Fourier transform, was used. The wrapping phase map is obtained and then unwrapped to remove the 2π ambiguity. In brief, the intensity of the modulated interferogram parallel to the y-axis can be expressed as:

(1)$$g(x,y) = \alpha (x,y) + \beta (x,y)\cos ({2\pi {f_0}x + \varphi (x,y)} )$$

where α(x,y), β(x,y) are the background illumination and the amplitude modulation of the fringes, respectively, f_o is the spatial carrier frequency, ϕ(x,y) is the phase modulation of the interferogram and x, y are axes in x and y, respectively. The expression of Eq. 1 can be re-writing as:

(2)$$g(x,y) = \alpha (x,y) + c(x,y){e^{i2\pi {f_0}x}} + {c^\ast }(x,y){e^{ - i2\pi {f_0}x}},$$

where, $c(x,y) = \beta (x,y){e^{i\varphi (x,y)}}/2,{c^\ast }(x,y) = \beta (x,y){e^{ - i\varphi (x,y)}}/2$.

The spectrum of the interferogram can be expressed as;

(3)$$G(f,y) = A(f,y) + C(f - {f_0},y) + {C^\ast }(f + {f_0},y).$$

Where C(f - f_o, y) and C*(f + f_o, y) are the carrier signals and A(f, y) is DC term. Figure 2(a) shows the one-dimensional (1D) profile along the three spectra in frequency domain, and Fig. 2(b) shows the selected spectrum in Fig. 2(a) after centering.

A filter window is used to select the carrier signals and shifting that component towards the origin. The inverse Fourier transform is then applied to the spectrum component C(f, y) to obtain c(x, y). The real and imaginary components of c(x, y) are expressed as:

(4)$$\textrm{Re} (c(x,y)) = \beta (x,y)\cos (2\pi {f_0}x + \varphi (x,y)),$$

(5)$${\mathop{\rm Im}\nolimits} (c(x,y)) = \beta (x,y)\sin (2\pi {f_0}x + \varphi (x,y)).$$

Fig. 2. (a) Spectra of a fringe pattern in frequency domain. (b) Selected frequency component at the center.

Download Full Size | PDF

The phase is retrieved by dividing Eq. (5) by Eq. (4) as:

(6)$$\begin{array}{l} \varphi (x,y) = {\tan ^{ - 1}}\left[ {\frac{{{\mathop{\rm Im}\nolimits} (c(x,y))}}{{Re (c(x,y))}}} \right] = \\ {\tan ^{ - 1}}\left[ {\frac{{\beta (x,y)\sin (2\pi {f_0}x + \varphi (x,y))}}{{\beta (x,y)\cos (2\pi {f_0}x + \varphi (x,y))}}} \right] = \\ {\tan ^{ - 1}}[{\tan ((2\pi {f_0}x + \varphi (x,y))} ]= 2\pi {f_0}x + \varphi (x,y). \end{array}$$

The obtained phase is limited between -π/2 and π/2 due to the arctangent function. The obtained wrapped phase values are then unwrapped by the graph cuts algorithm.

2.2. Wavelet transform method

The main concept of the wavelet transform [20] is to divide a modulated signal into its various scaled and shifted versions of a mother wavelet. Here, we used the Morlet wavelet transform method to demodulate the off-axis interferogram row by row. The Morlet wavelet provides better localization in both the frequency and spatial domains. A row of the modulated interferogram is project onto the daughter wavelets, which are obtained by both translating and dilating a single basis function, called the mother wavelet. The resulting transformation, which is a 2D complex array, contains the required phase information of that deformed interferogram’s row, which can be retrieved by determining the ridge of the transform. A modulus array and phase array can be calculated from the resulting two-dimensional complex array. The direct maximum method can be used to extract the ridge from the modulus array. The Morlet wavelet is a plane wave modulated by a Gaussian function, and it is expressed as:

(7)$$\psi (x) = {\pi ^{1/4}}{e^{i\alpha x}}{e^{ - {x^2}/2}},$$

where i = $\sqrt { - 1} $, α is a fixed spatial frequency in the range of 5 to 6 to satisfy the admissibility condition. Daughter wavelets ψ_b,a are built by translation on the x-axis by b (shifting) as a location parameter and by a as a dilation parameter (scaling) of the mother wavelet ψ(x) as given by

(8)$${\psi _{b,a}}(x) = \frac{1}{a}\psi \left( {\frac{{x - b}}{a}} \right).$$

The modulated signal ƒ(x) multiplied by scaled, shifted versions of the mother wavelet. The result are wavelet coefficients that are a function of the scale a and the shifting b and can be expressed as;

(9)$$D(a,b) = \frac{1}{a}\int\limits_{ - \infty }^\infty {{\psi ^ \ast }\left( {\frac{{x - b}}{a}} \right)f(x)dx} ,$$

where ƒ(x) in our case represents a row of a fringe pattern, * denotes complex conjugation. These coefficients refer to the closeness of the signal ƒ(x) to the wavelet at a particular scale a and translation b. The modulus ∣D(a,b)∣ and the phase ϕ(a,b) arrays are expressed as

(10)$$|{D(a,b)} |= \sqrt {{{[{\textrm{Re} (D(a,b))} ]}^2} + {{[{{\mathop{\rm Im}\nolimits} (D(a,b))} ]}^2},}$$

(11)$$\varphi (a,b) = {\tan ^{ - 1}}\left( {\frac{{{\mathop{\rm Im}\nolimits} (D(a,b))}}{{\textrm{Re} (D(a,b))}}} \right).$$

The obtained wrapped phase values are then unwrapped by the graph cuts algorithm.

2.3. Angular spectrum method

In scalar diffraction theory, the diffracted field Ψ_p (x, y; z) can be obtained from the optical field of the hologram plane Ψ_p0 (x, y) as:

(12)$${\Psi _p}(x,y;z) = {\Im ^{ - 1}}\{{\Im \{{{\Psi _{p0}}(x,y)} \}T({{f_x},{f_y};z} )} \},$$

where the reference beam is assumed to be a plane wave and T (f_x, f_y; z) is the spatial frequency transfer function of propagation through a distance z expressed as [22–24]:

(13)$$T({f_x},{f_y};z) = \textrm{exp} \left( {j\frac{{2\pi }}{\lambda }z\sqrt {1 - {{(\lambda {f_x})}^2} - {{(\lambda {f_y})}^2}} } \right),$$

where j = $\sqrt { - 1} $, f_x and f_y are called spatial frequencies in x and y, λ = 635 nm is the illumination wavelength, and z is the propagation phase factor. The amplitude and phase of the object have been calculated directly using MATLAB as (abs (Ψ_p)) and (angle (Ψ_p)), respectively. The obtained wrapped phase values are then unwrapped by the graph cuts algorithm.

3. Experiment

Depending on the sample being investigated, two setup configurations were used. A reflection configuration type was used for opaque objects (e.g. metallic surfaces), and a transmission configuration type was used for transparent samples (e.g. biological cells). The experimental arrangement for reflection and transmission configurations are depicted in Figs. 3(a) and (b), respectively.

Fig. 3. Mach-Zehnder configurations at (a) reflection and (b) transmission; O, object wave, R, reference wave; NPBS, non-polarized beam splitter; M₁ and M₂, mirrors; OC, objective condenser lens; MO, microscope objective lens; CCD, charge-coupled device.

Download Full Size | PDF

The configurations consist in a Mach-Zehnder interferometer setup combined with a microscope objective (MO) and a black-and-white high-resolution USB CCD camera from Thorlabs (1024 × 768 pixels with 4.65 µm × 4.65 µm pixel size). The light beam emitted by the laser source (laser diode 635nm from Thorlabs) is collimated by an expander resulting in a collimated beam of size around 25.4mm. In case of reflection configuration (Fig. 3(a)), the collimated beam is directed to a non-polarizing beam splitter NPBS₁ (400–700nm) where the beam is divided into two copies: the reference beam (R) and the object beam (O). The object beam O reflected from a flat mirror M₁ (size of 1 inch and flatness of λ/20) to the object being investigated via a non-polarizing beam splitter NPBS₂ (400–700nm). The reference beam R reflected from a flat mirror M₂ (size of 1 inch and flatness of λ/20) to NPBS₂. The object beam passes through NPBS₂ to the object being investigated and then reflected to interfere with the reference beam at NPBS₂. The object is a sample mounted on an xyz-translation stage placed at a distance z from the hologram plane, whose magnified image is projected on the CCD camera, as well as the reference beam. In case of transmission configuration (Fig. 3(b)), the collimated beam is divided into two parts: the reference beam (R) and the object beam (O). After having interacted with the sample, the object beam is collected by the MO (20x, 0.4NA), recombined with the reference wave, and the interference is recorded on the CCD camera. Here, the reference lens matches the curvature induced on the object wave by the MO. The lateral magnification, M = - f_OC/ f_MO, where f_OC is the focal length of the objective condenser lens, and f_MO is the focal length of the microscope objective. The total magnification can be calculated by the distance between object, MO and the CCD, or using a standard resolution test version. As seen in both configurations, no moving parts are presents in the optical design resulting in simplified operations, improved accuracy, high robustness and real-time imaging with a single acquisition. When the object wave with intensity I_O and a reference wave with intensity I_R interfere, the intensity of the interference pattern recorded on the CCD may be expressed as

(14)$${I_{{\mathop{\rm int}} }} = {I_O} + {I_R} + \sqrt {{I_O}{I_R}} \cos (\varphi ),$$

where ϕ is the phase difference between the two waves. This Eq. (14) is exactly the same of Eq. (1) with α (x, y) I_O + I_R, β (x, y) (I_O I_R)^0.5, and f_o = 0 (on-axis interferogram). The modulated interferogram I_int (x, y) is then corrected by the flat fielding and the obtained interferogram is I_flat (x, y) expressed as [25–27]:

(15)$${I_{flat}} = \frac{{W({{I_{{\mathop{\rm int}} }} - {I_{th}}} )}}{{({{I_f} - {I_{th}}} )}},$$

where W is the average pixel value of the corrected flat field interferogram, I_int is the raw interferogram, I_th is the thermal noise image (Fig. 4(b)), and I_f is the flat frame image (Fig. 4(a)). The I_flat is then convolved with windowed Fourier transform (WFT) expressed as:

(16)$$WFT({f_x},{f_y}) = \int\limits_{ - \infty }^{ + \infty } {\int\limits_{ - \infty }^{ + \infty } {{I_{flat}}(x,y){r^ \ast }(x,y)dxdy} ,}$$

where r(x,y) exp[-x²/2$\sigma _x^2 - $y²${/}2\sigma _y^2]$ is the Gabor transform, which is a Gaussian function and the symbol * denotes the complex conjugate operation. Here the ${\sigma _x}$ and ${\sigma _y}$ are the standard deviations of the Gaussian function in x and y directions, and selected to be 10, respectively. The Gaussian window is selected because it has higher emphasis on the data near the window center and lower emphasis on the data further away from the window center. By shifting the central position of the Gaussian window point by point and computing the WFT of the signal, the fundamental spectral component of each local signal can be extracted. Then, by integrating all these spectral components, a universal spectral component can be retrieved as:

(17)$$\int\limits_{ - \infty }^{ + \infty } {\int\limits_{ - \infty }^{ + \infty } {WFT({f_0},{f_0})dxdy} = F({f_0},{f_0})} .$$

Finally, applying the inverse Fourier transform to the universal spectral component, a complex image is obtained.

The real values of the complex image have been taken to provide a speckle free interferogram. The obtained speckle free interferogram is then demodulated here by the FT method or by the WL method or by the ASM to produce wrapped phase map of the object being screened. The wrapped phase map is then unwrapped by the graph cuts algorithm. In preparation for measurement, all precautions have been taken into account to minimize the systematic error [10].

Fig. 4. (a) Flat frame image captured when the object arm of the interferometer is blocked. (b) Thermal noise image captured when there is no illumination. The standard deviation of (a) is calculated to be σ = 11.96. The σ of (b) is calculated to be σ = 1.09. The color bar unit is (a.u.).

Download Full Size | PDF

Fig. 5. (a) Original off-axis interferogram of step height of 200µm. (b) Correction of (a) with WFF. (c) Correction of (a) with flat fielding. (d) Correction of (a) with both WFF and flat fielding. The interferogram is an image of 512 × 512 pixels.

Download Full Size | PDF

4. Performance of the proposed method

To introduce the performance of the proposed method in screening the objects, we investigated two opaque objects and two transparent objects. The two opaque objects are large scale step height of 200µm and small scale step height of 1.50µm investigated by (Fig. 3(a)) and reconstructed by the FT method after correction with the proposed method. The two transparent objects are normal blood cells from healthy human donor and cancerous blood cells from unhealthy (with leukemia) human donor investigated by (Fig. 3(b)) and reconstructed by the ASM after correction with the proposed method. Here, we used a single wavelength, however, the proposed method can be applied to the synthetic interferogram of different wavelengths.

4.1. Screening opaque objects

The opaque object of step height of 200µm was positioned gently in the interferometer and then adjusted carefully to see the perfect off-axis interferogram. Figure 5(a) shows average of 100 off-axis interferograms with standard deviation σ = 32.59. Here, we focused on using σ, since there is a direct relationship between σ and signal to noise ratio (SNR) as SNR = 10log10 (0.5/σ²). The σ and SNR of the interferogram can be calculated directly by using the function std2 and (SNR), respectively in MATLAB. The function std2 computes the standard deviation of the array image using std (image (:)). Figure 5(c) shows correction of Fig. 5(a) with flat fielding using Eq. (15). The σ of Fig. 5(c) is calculated to be σ = 25.11. Figure 5(b) shows correction of Fig. 5(a) with WFF using Eqs. (16) and (17) only without applying flat fielding. The σ of Fig. 5(b) is calculated to be σ = 17.77. Figure 5(d) shows correction of Fig. 5(a) with combination of both flat fielding and WFF using Eqs. (15), (16) and (17), respectively. The σ of Fig. 5(d) is calculated to be σ = 11.22. As seen in Fig. 5(d), the proposed method clears the interferogram from the coherent noise and increases the SNR of the interferogram by around 30%. To extract the phase map of the off-axis interferogram, we reconstructed Fig. 5(c) as an example by the FT method in the absence and in the presence of the object and the difference is then extracted.

Figure 6(a) is the off-axis interferogram in the absence of the object and Fig. 6(b) is its wrapped phase map. Figure 6(c) is the wrapped phase map of Fig. 5(a). Figure 6(d) is the unwrapped phase map of Fig. 6(b), while Fig. 6(e) is the unwrapped phase map of Fig. 6(c). Figure 6(f) is the difference between Fig. 6(d) and Fig. 6(e).

Fig. 6. (a) Off-axis interferogram in the absence of the object. (b) Wrapped phase map (a). (c) Wrapped phase map in the presence of the object. (d) Unwrapped phase map (b). (e) Unwrapped phase map of (c). (f) Difference between (d) and (c). The color bar unit is in radians.

Download Full Size | PDF

To see the performance of the proposed method, we reconstructed the three off-axis interferograms of Figs. 5(b), (c), and (d) by the FT method following the same procedure of Fig. 6. Figures 7(a), (c), and (e) show the three-dimensional (3D) pseudocolor of the reconstructed phase maps of Figs. (a), (b), and (d), respectively. Figures 7(b), (d), (f) are zoomed in the black rectangles in Figs. 7(a), (c), (e), respectively. As seen in the black circles of Figs. 7(b) and (d), some spikes are appeared. Such spikes are disappeared completely in Fig. 7(f). This confirms that the proposed method, which is a combination of flat fielding and WFF removes the spikes completely. An interesting observation is that the lateral resolution, which is represented by slope of the roll-off is slightly improved without applications of digital filters [5].

Fig. 7. (a), (c), (e) are 3D pseudocolor of the reconstructed phase maps of Figs. 5(a), (b), (d), respectively. (b), (d), (f) are zoomed in the black rectangles in (a), (c), (e), respectively. The color bar unit is in radians.

Download Full Size | PDF

To see such improvement in lateral resolution, three phase profiles along the red, blue, and black lines in Figs. 7(b), (d), and (f) were extracted and plotted together as seen in Fig. 8(a). As seen in Fig. 8(a), the black profile is steeper than the blue profile, which confirms the performance of the proposed method in enhancing the lateral resolution, which is represented by slope of the roll-off without applications of digital filters.

Fig. 8. (a)1D phase profiles along the red, blue, and black lines of Figs. 7(b), (d), and (f) extracted by FT method. (b) 1D phase profiles along the red line of Fig. 7(b) extracted by the FT method and the WL method.

Download Full Size | PDF

As seen in Fig. 8(a), the application of flat fielding and WFF makes the roll-off of the object shifts to the interior edge. This shift (around 20 pixels) is due to the existence of two edges (interior and exterior) as shown in Fig. 5(a). Without application (W/O) of flat fielding and WFF the roll-off moves to the exterior edge (red profile of Fig. 8(a)). When we applied the Morlet wavelet method to check the right location of the profile, we found that the roll-off of the object moves to the exterior edge. Figure 8(b) shows the 1D phase profiles of Fig. 7(b) without applications of flat fielding and WFF reconstructed by the FT method (red profile) and by the Morlet WL method (purple profile). The Morlet wavelet transform method can spatially relocalize the roll-off of the object to the exterior edge. It is worth mentioning that the lateral resolution is limited by coherent diffraction limit, which is defined by NA of the system and the wavelength [28]. The opaque object of step height of 1.50µm was positioned gently in the interferometer and then adjusted carefully to capture the perfect parabolic interferogram. Since we used single wavelength, so chromatic aberration is nearly negligible [29]. Figure 9(a) shows average of 100 parabolic interferograms with standard deviation σ = 12.8. The corrected interferogram of Fig. 9(a) by the proposed method is shown in Fig. 9(b) with σ = 5.3. The parabolic interferograms of Fig. 9(a) and Fig. 9(b) were reconstructed by the FT method [10]. If the original parabolic interferogram is denoted as I_p(x,y) and the corrected interferogram by the proposed method is denoted as I_pf (x,y). The I_pf (x,y) is digitized after 2D sampling and the digitized interferogram I_p($\aleph$,l) is reconstructed by computing the Fresnel integral as:

(18)$$\begin{array}{l} \psi \textrm{(m,n)} = A\textrm{exp}\left[ {\frac{{\textrm{i}\pi }}{{\lambda \textrm{d}}}({m^2}\Delta {\xi^2} + {n^2}\Delta {\eta^2})} \right]\\ \quad \quad \quad \times FFT{\left\{ {{R_D}(\aleph ,l){I_p}(\aleph ,l) \times \textrm{exp} \left[ {\frac{{\textrm{i}\pi }}{{\lambda \textrm{d}}}({\aleph^2}\Delta {x^2} + {l^2}\Delta {y^2})} \right]} \right\}_{\textrm{m,n}}}\textrm{,} \end{array}$$

where ψ is the reconstructed wave front, Δx and Δy are the sampling intervals in the interferogram plane, m, n, $\aleph$, l are integers, λ is the wavelength of the laser light, A is a complex constant given as A = exp (i2π/λ)/(iπλ), FFT denotes the fast Fourier transform, and Δξ and Δη are the sampling intervals in the observation plane defined as Δξ = Δη = λd/L. The digital reference wave R_D ($\aleph$, l) is defined as A exp [(i2π/λ) (k_x$\aleph$Δx + k_y lΔy)]. The wave vectors k_x and k_y are adjusted such that the propagation direction of R_D($\aleph$,l) agrees as close as possible with that of the experimental reference wave. The phase of the simulated digital lens is expressed as φ(m,n) exp [(-iπ/λD) (m²Δξ² + n²Δη²⁾], where D is expressed as 1/D = 1/d (1 + d₀/d), where d₀ is the distance between the microscope objective and the object, d is the distance between the object and the CCD camera of size (L). The array of the complex numbers of the conjugate digital lens is multiplied by the reconstructed wave front to compensate for the wave front curvature. The result of multiplication is an array of complex numbers denoted as Ψ. The phase contrast image is obtained using MATLAB as arctan {Re(Ψ)/Im(Ψ)}. The obtained phase is then unwrapped by the graph cuts algorithm. Figures 9(c) and (d) show the reconstructed phase-contrast images of Fig. 9(a) and Fig. 9(b), respectively. To see the performance of the proposed in removing the ripples, phase profiles along the red and blue lines of Figs. 9(c) and (d) were extracted and plotted together as shown in Fig. 9(e). As seen in Fig. 9(d), the maximum ripple in the red phase profile is in the range of 0.6 rad, corresponds to a path length of (0.6*635)/ (4*3.14) ≅ 30nm. Such ripples, which represented by the blue profile in Fig. 9(e), are minimized to the least measure by the proposed method.

Fig. 9. (a) Original parabolic interferogram of step height of 1.50µm. (b) Correction of (a) with flat fielding and WFF. (c) Reconstructed phase-contrast image of (a). (d) Reconstructed phase-contrast image of (b). 1D phase profiles along the red and blue lines of (c) and (d).

Download Full Size | PDF

4.2. Screening transparent objects

A cell is considered typically a complex 3D object consisting of many subcellular components and organelles. It has been discovered that the refractive index values vary across these minute biological structures. The value of the refractive index quantifies the speed of light in media which is related to the regional concentration of protein (which has a higher value of the refractive index). The phase map ϕ(x, y) at pixels (x, y) contains information about the cell thickness h(x,y) and the difference of refractive index between the cells n and the surrounding medium n₀ and can allow quantification of cellular morphology under experimental conditions [30,31]. The refractive index n of the red blood cell (RBC) reflects the cellular internal composition, which is mostly defined by the hemoglobin concentration, while the cell thickness h(x,y) is strongly impacted by its morphology [31]. The thickness of the cell at the (x, y) position can be calculated as

(19)$$h(x,y) = \lambda ({\varphi (x,y)/2\pi } )/({n - {n_0}} ),$$

where λ is the wavelength of illumination light. Generally, the physical thickness of cells and different refractive index distributions will cause different phase shifts of light waves, so this kind of technology can be utilized to observe the living cells [32,33]. The refractive index of the RBCs equals nearly 1.42 [31], so the physical thickness of the cell at the (x, y) position can be calculated easily once the quantitative phase map ϕ(x, y) at pixels (x, y) is extracted. To extract the quantitative phase map of the object, ASM was used to reconstruct the inline interferogram. Unlike the Fresnel diffraction method, the ASM is valid for a shorter distance. For the preparation of the blood samples, blood from healthy human donor and unhealthy human (with leukemia) donor were obtained from the blood bank of the University Hospital of the Cairo University, Cairo, Egypt. The blood cells were displaced on glass substrates for measurements. The blood samples were investigated by the interferometer whose configuration is shown in Fig. 3(b). For the two samples, two inline interferograms were captured by the CCD camera and corrected by the proposed method then reconstructed by the ASM. Figure 10(a) shows off-axis interferogram (752 × 480) pixels of the blood healthy human donor sample captured when there is a slightly small angle between the reference and object beams.

Figure 10(b) shows the corrected inline interferogram of (100 × 93) pixels of the red rectangle of (a) when there is no angle between the reference and object beams. Figure 10(c) shows the reconstructed phase map of Fig. 10(b) with the ASM. The image plane of the MO is located approximately 50mm behind the camera. Figure 10(d) shows 3D pseudo color of an individual cell in the white square of Fig. 10(c) with size of (12 × 12) pixels. As seen in Fig. 10(d), the normal blood cell has a flattened biconcave disc of about 8µm in diameter with no deformation in its membrane. The phase profile along the black line of Fig. 10(d) is shown in Fig. 12(a). The maximum height of the cell at its center is in the range of 2.7 radians. Figure 11(a) shows the corrected inline interferogram (400 × 420 pixels) of the blood unhealthy human donor captured when there is no angle between the reference and object beams. Reconstructed amplitude-contrast and phase-contrast images of Fig. 11(a) are shown in Fig. 11(b) and Fig. 11(c), respectively. The image plane of the MO is located approximately 22mm behind the camera. Figure 11(d) shows 3D pseudocolor of an individual cell in the red rectangle of Fig. 11(c). As seen in Fig. 11(d), the membrane of the cell is clearly deformed. The 1D phase profile of Fig. 11(d) along the blue line of (d) is shown in Fig. 12(b). The maximum height is in the range of 2.7 radians. As seen in Figs. 12(a) and (b), the maximum heights of both the healthy and cancerous cells are congruent. The main difference between the two cells appears in the deformation of the membrane. The membrane of the healthy cell appears regular, however, the membrane of the cancerous cell appears deformed. As seen in Figs. 12(a) and (b), the maximum thickness at its center were estimated to be in the range of 650nm (using Eq. (19) when the refractive index of the RBCs equals 1.42 [31]. We claim that the deformation in the cancerous cell of Fig. 11(d) and its profile of Fig. 12(b) is due to changes in the packing properties of the cell [34,35].

Fig. 10. (a) Off-axis interferogram of size (752 × 480) pixels. (b) Corrected inline interferogram of the red rectangle of (a) with size of (100 × 93). (c) Reconstructed phase map of (b) by the ASM with a propagation phase factor z = −50 mm. (d) 3D pseudocolor of the white rectangle of (c) with size of (12 × 12) pixels. The phases in (c) and (d) are in radians.

Download Full Size | PDF

Fig. 11. (a) Inline interferogram of size (400 × 420) pixels. (b) Amplitude-contrast image of (a) reconstructed by the ASM. (c) Phase-contrast image of (a) reconstructed by the ASM with a propagation phase factor z = −22 mm. (d) 3D pseudocolor of the selected rectangle of (c) with size of (100 × 93) pixels. The phases in (c) and (d) are in radians.

Download Full Size | PDF

Fig. 12. 1D phase profiles along the (a) black line of Fig. 10(d) and (b) the blue line of Fig. 11(d). Blue profile in (b) is the original extracted phase, while the red profile is its quadratic fitting.

Download Full Size | PDF

Deformation in the cancerous cell membrane may be attributed to changes in the packing properties of the phospholipid macromolecules forming the RBCs cell wall leading to decreased cell membrane elasticity. An interesting observation is viewing bulges in the cancerous cell of Fig. 11(d).

Since the general shape of the normal blood cells is a discocyte, this gives an evidence that the screened cell of Fig. 11(d) is cancerous cell and it seems to be echinocyte type I with 9 bulges (E1-9). Echinocytes in general are characterized by one or more spicules. Echinocyte type I with 9 bulges (E1-9) is characterized by a 9-fold rotation axis plus up-down symmetry, reflecting the nine identical bulges that develop on the rim. Types of echinocytes with different bulges are reported in [36]. The agents of echinocytes shapes may due to cholesterol addition, high salt (hypertonic saline), high pH, intracellular adenosine triphosphate depletion, or proximity to glass.

5. Conclusion

Phase-contrast images associated with coherent phase noise rendering it unsuitable for quantitative phase-contrast imaging with high precision. To achieve precise phase measurement, phase noise should be minimized to the least measure. In this paper, we proposed a method to suppress the coherent phase noise, especially the ripples to the least measure. The method is based on a combination of flat fielding and windowed Fourier filtering demodulated by the graph cuts algorithm. The flat fielding removes dust particles attached with the interferogram as well as adjusts inhomogeneity of its intensity. The corrected interferogram has been convolved by the windowed Fourier filtering to produce an image free from coherent noise. The corrected interferogram is then demodulated by the FT method, or WT method, or ASM to extract the wrapping phase map of the object being screened. The graph cuts algorithm has been used to unwrap the wrapping phase map to remove the 2π ambiguity via segmentation of the wrapping phase image into background and foreground images and the procedure continuous for multiple regions simultaneously until no further segment is found. We introduced the performance of the proposed method by screening two opaque objects and two transparent objects. The two opaque objects were screened by the FT method, while the two transparent objects were screened by the ASM. Experimental results showed that the proposed method has the ability to minimize the phase noise, especially the ripples to the least measure. Quantitatively, the measured error due to ripples has been estimated to be in the range of 30nm. In case of objects having interior and exterior edges, the proposed method can shift the roll-off of the object to the interior edge. The Morlet wavelet transform method has been used to spatially relocalize the roll-off of the object to the exterior edge. By using the proposed method with single shot inline interferogram, we clearly viewed the bulges, which have the sharpest spatial features in the echinocytes. Conventionally, the reconstructed inline interferogram is always blurred by the transmitted zeroth-order light and twin image. By the proposed method, the zero-order and twin image were mitigated drastically. Moreover, the lateral resolution, which is represented by slope of the roll-off is slightly improved without applications of digital filters.

Funding

Science, Technology & Innovation Funding Authority (STDF) (40490).

Disclosures

The author declares no conflicts of interest.

Data Availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

References

1. P. Jacquot and J. -M. Fournier (Editors), “interferometry in speckle light: theory and applications,” Proceedings of the international conference 25-28 September, Springer, Lausanne, Switzerland (2000).

2. R. C. Gonzalez and R. E. Woods, Digital Image Processing, 3rd edition, Pearson Education Inc., Upper Saddle River, pp. 344–357 (2010).

3. H. A. Aebischer and S. Waldner, “A simple and effective method for filtering speckle-interferometric phase fringe patterns,” Opt. Comm. 162(4-6), 205–210 (1999). [CrossRef]

4. Q. Kemao, L. T. H. Nam, L. Feng, and S. H. Soon, “Comparative analysis on some filters for wrapped phase maps,” Appl. Opt. 46(30), 7412 (2007). [CrossRef]

5. D. G. A. Ibrahim, “Improving the intensity-contrast image of a noisy digital hologram by convolution of Chebyshev type 2 and elliptic filters,” Appl. Opt. 60(13), 3823–3829 (2021). [CrossRef]

6. F. Charriere, T. Colomb, F. Montfort, E. Cuche, P. Marquet, and C. Depeursinge, “Shot-noise influence on the reconstructed phase image signal-to-noise ratio in digital holographic microscopy,” Appl. Opt. 45(29), 7667–7673 (2006). [CrossRef]

7. Qian Kemao, “Windowed Fringe Pattern Analysis,” Chapter 3, Spie press, Bellimgham, Washington, USA, (2013).

8. J. M. Bioucas-Dias and G. Valadão, “Phase unwrapping via graph cuts,” IEEE Trans. on Image Process. 16(3), 698–709 (2007). [CrossRef]

9. Guide to the Expression of Uncertainty in Measurement. ICS 03.120.30; 17.020 (1996).

10. D. G. A. Ibrahim, “Estimation of an uncertainty budget and performance measurement for a dual-wavelength Twyman-Green interferometer,” J. Microsc. 282(3), 224–238 (2021). [CrossRef]

11. D. G. A. Ibrahim, “Steep Large Film Thickness Measurement with Off-axis Terahertz Digital Holography Reconstructed by a Direct Fourier and Hermite polynomial,” Appl. Opt. 57(10), 2533–2538 (2018). [CrossRef]

12. D. G. A. Ibrahim, “Enhancement of the steepness measurement of a film thickness edge using wavelet transforms with fringe thinning,” OSA Continuum 3, 1928–1937 (2020). [CrossRef]

13. D. G. A. Ibrahim, “Demodulation of a parabolic interferogram in time domain for rough surface characterization,” Appl. Phys. B: Lasers Opt. 126(9), 146 (2020). [CrossRef]

14. D. G. Abdelsalam and D. Kim, “Two-wavelength In-line Phase-shifting Interferometry Based on Polarizing Separation for Accurate Surface Profiling,” Appl. Opt. 50(33), 6153–6161 (2011). [CrossRef]

15. D. G. A. Ibrahim and T. Yasui, “Multi-object investigation using two-wavelength phase-shift interferometry guided by optical frequency comb,” Appl. Phys. Lett. 112(17), 171101–171104 (2018). [CrossRef]

16. Dahi Ghareab Abdelsalam Ibrahim, “Optical Metrology with Interferometry,” Cambridge Scholars Publishing, London, (2019).

17. J. Zhong and J. Weng, “Spatial carrier-fringe pattern analysis by means of wavelet transform: wavelet transform profilometry,” Appl. Opt. 43(26), 4993–4998 (2004). [CrossRef]

18. Q. Kemao, “Windowed Fourier transform for fringe pattern analysis,” Appl. Opt. 43(13), 2695–2702 (2004). [CrossRef]

19. Q. Kemao, “Windowed Fourier transform method for demodulation of carrier fringes,” Opt. Eng. 43(7), 1472–1473 (2004). [CrossRef]

20. Qian Kemao, H. Wang, and W. Gao, “Windowed Fourier transform for fringe pattern analysis: Theoretical analyses,” Appl. Opt. 47(29), 5408–5419 (2008). [CrossRef]

21. P. Rastogi and E. Hack, “Phase estimation in optical interferometry”, CRC press, Taylor and Francis group, Boca Raton, London, New York, (2017).

22. M. Takeda, H. Ina, and S. Kobayashi, “Fourier-transform of fringe-pattern analysis for computer-based topography and interferometry,” J. Opt. Soc. Am. 72(1), 156 (1982). [CrossRef]

23. T.-C. Poon and P. P. Banerjee, Contemporary Optical Image Processing with MATLAB (Elsevier, Oxford, UK, 2001).

24. M. K. Kim, L. Yu, and C. J. Mann, “Interference techniques in digital holography,” J. Opt. A: Pure Appl. Opt. 8(7), S518–S523 (2006). [CrossRef]

25. D. G. Abdelsalam, R. Magnusson, and D. Kim, “Single-shot, Dual-wavelength Digital Holography Based on Polarizing Separation,” Appl. Opt. 50(19), 3360–3368 (2011). [CrossRef]

26. D. G. Abdelsalam and D. Kim, “Coherent noise suppression in digital holography based on flat fielding with apodized apertures,” Appl. Opt. 19(19), 17951 (2011). [CrossRef]

27. D. G. Abdelsalam and D. Kim, “Real-time dual-wavelength digital holographic microscopy based on polarizing separation,” Opt. Express 285(3), 233–237 (2011). [CrossRef]

28. Y. Cotte, F. Toy, P. Jourdain, N. Pavillon, D. Boss, P. Magistretti, P. Marquet, and C. Depeursinge, “Marker-free phase nanoscopy,” Nat. Photonics. 7(2), 113–117 (2013). [CrossRef]

29. D. G. A. Ibrahim, “Simultaneous dual-wavelength digital holographic microscopy with compensation of chromatic aberration for accurate surface characterization,” Appl. Opt. 58(23), 6388–6395 (2019). [CrossRef]

30. M. Bardyn, B. Rappaz, K. Jaferzadeh, D. Crettaz, J.-D. Tissot, I. Moon, G. Turcatti, N. Lion, and M. Prudent, “Red blood cells ageing markers: A multi-parametric analysis,” Blood Transfus. 15, 239–248 (2017). [CrossRef]

31. B. Rappaz, A. Barbul, Y. Emery, R. Korenstein, C. Depeursinge, P. J. Magistretti, and P. Marquet, “Comparative study of human erythrocytes by digital holographic microscopy, confocal microscopy, and impedance volume analyzer,” Cytometry 73A(10), 895–903 (2008). [CrossRef]

32. C. L. Curl, C. J. Bellair, P. J. Harris, B. E. Allman, A. Roberts, K. A. Nugent, and L. M. D. Delbridge, “Quantitative phase microscopy: A new tool for investigating the structure and function of unstained live cells,” Clin. Exp. Pharmacol. Physiol. 31(12), 896–901 (2004). [CrossRef]

33. S. Bhatt, A. Butola, S. R. Kanade, and A. Kumar, “High-resolution single-shot phase-shifting interference microscopy using deep neural network for quantitative phase imaging of biological samples,” J. Biophotonics 14(7), 1–9 (2021). [CrossRef]

34. D. Kabaso, E. Gongadze, S. Perutkova, C. Matschegewski, V. Karalj-Iglic, U. Beck, U. Van Rienen, and A. Iglic, “Mechanics and electrostatics of the interactions between osteoblasts and titanium surface,” Comput. Methods Biomechanics Biomed. Eng. 14(5), 469–482 (2011). [CrossRef]

35. A. G. Petrov, “Flexoelectricity of model and living membranes,” Biochim. Biophys. Acta 1561(1), 1–25 (2002). [CrossRef]

36. G. Gompper and M. Schick, “Soft Matter: Lipid Bilayers and Red Blood Cells,” WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim, (2008).

Quantitative phase imaging using a combination of flat fielding and windowed Fourier filtering demodulated by a graph cuts algorithm for screening opaque and transparent objects

Abstract

1. Introduction

2. Single shot reconstruction methods

2.1. Fourier transform method

2.2. Wavelet transform method

2.3. Angular spectrum method

3. Experiment

4. Performance of the proposed method

4.1. Screening opaque objects

4.2. Screening transparent objects

5. Conclusion

Funding

Disclosures

Data Availability

References

Data Availability

Cited By

Figures (12)

Equations (19)

Optics Continuum