1. Introduction

At microscopic scale, unstained biological tissues are often transparent, so that the main alteration that light suffers when it passes through them are in its phase. Given the numerous biological processes which are best investigated by phase imaging, this minimally invasive technique have been the object of constant research in Optics. In recent years, systems have been proposed that can detect not only phase changes − as the traditional methods of phase contrast do − but also measure the phase variations quantitatively, allowing the detailed measurement of variation in the refraction index in the samples.

Among the many techniques proposed for quantitative phase imaging, a number of methods are based on digital holography (DH) [1,2]. A common characteristic of DH [3] is the use of an image sensor to register an interferogram, formed by bringing together a reference beam and an object beam. From the acquired image, after applying the diffraction theory, the object − a complex valued matrix − is reconstructed in the computer by calculating the field that will propagate from the hologram assuming as illumination a specific function. The phase of the object − the complex argument − is computed using the inverse tangent function. If the value of the phase is greater that $2\pi$, the true phase is not obtained, only a map of the wrapped phase, i.e., its principal value. As is common to other interferometric techniques, to obtain the real phase, algorithms are needed to process these maps by adding to each pixel the correct multiple of $2\pi$radians. This is not a trivial task given that the algorithms must be able of deal with noise, under-sampling, phase dislocations and discontinuities to produce a path-independent phase map. Over the years, numerous phase unwrapping strategies have been proposed but the problem is still challenging [4,5] and has led to special DH systems being proposed to unwrap the phase optically, rather than numerically [6].

However, it is interesting to notice that the phase is not intrinsically wrapped in DH: wrapping only arises when the complex argument is computed from the matrix representing the object. Might it be possible to obtain the phase without computing the principal value? A method to do so is proposed in this paper.

Given that the numerical reconstruction of the digital hologram provides a complex representation of the field in a given plane, diffraction theory can be used again to propagate the field to another plane or through any digitally modeled optical system. Unfortunately, there are no optical systems able to directly provide the phase from the light field but there are systems that can produce measurements of the phase gradient that can be easily modeled in the computer. Then, once the gradient has been obtained, error-minimization integration algorithms can be used to provide an estimate of the unwrapped phase. Two advantages should be emphasized in this approach; first, the inverse tangent function is never used so phase jumps never appear, and secondly, the integration algorithms can adequately manage the noise transferred to the gradient.

Iterative methods to obtain information of the phase gradient from interference holography has been proposed before [7]. Also, the advantages of using the gradient of the wrapped phase instead of the wrapped phase for pixel-by-pixel unwrapping has been suggested and applied to process interferometric measurements using iterative algorithms [4] but, as mentioned above, in DH the situation is different because it is possible to obtain directly the gradient unaffected by the wrapping.

Of the different possibilities [8] that can be considered for implementing a computational version of an optical sensor to obtain directly gradient information from the field, an especially convenient one is the pyramid sensor [9–11]. As we will show, it is easy to implement and can produce high-resolution gradient maps for integration. A brief description of the sensor basis and its computational implementation is presented in the next section.

2. A computational implementation of the pyramid sensor for digital holography

The pyramid sensor is closely related with the Foucault knife test, which exploits the relation between the gradient of the phase in a plane in front of a lens and the amplitude distribution in the Fourier plane that can be found at a focal distance from the lens. This relation permits to make an approximate assessment of the phase gradient by blocking (or selecting) part of the focal plane, observing the emerging light with the naked eye or an additional lens and a projection screen. In a similar way, the pyramid sensor selects parts of the field in the Fourier plane to provide, in this case, not an estimate but the complete phase gradient.

A detailed diffraction optics description of the pyramid sensor can be found elsewhere [12]; we offer below a brief intuitive explanation of the sensor basis. Lets us first consider an input field consisting of a constant phasor, $\exp \left(imx\right)$, i.e., a phase ramp of slope . Its Fourier transform is $\delta \left(m/2\pi +\nu \right)$, a delta function placed at a distance from the Fourier plane origin. So, it is the same to measure the phase slope as to measure the displacement of the delta function from the origin. A convenient method to measure the displacement from the center of a spot in a plane is to use a quad-cell sensor, which obtains this information by comparing the relative energy falling on each one of four photosensors arranged as a squared array, in our case, placed on the Fourier plane. The problem with the quad-cell in this particular case is that the spot is a function without extension and so can only fall on one of the four photosensors; one solution is to modulate the signal by moving the quad-cell in a symmetric path – whose size will determine the measurement gain and range − around the axis, averaging the signals generated by each photosensor. These quad-cell signals can also be obtained as amplitudes in the input-plane. To do that, the field in the Fourier plane can be spliced into four parts, each of them corresponding to a quadrant of the plane, then back-propagate the fields to the input plane by using an inverse Fourier transform and to compute the modulus. In the example we are using here − a constant phase − only one of the four amplitude distributions will be constant and different from zero. If the modulation is introduced − by moving the origin of the quadrant selection around the origin of the Fourier plane − the value of the constant amplitude will be proportional to the slope.

A more complex case is a field in the input plane consisting of a constant amplitude, A, and a phase function. In a small area located at a given position in the plane, the field can be approximated by the phasor, $A\exp \left(i\left(mx+\eta \right)\right)$, where $\eta$ is a constant, represented in the Fourier plane by the expression$A\exp \left(i\eta \right)\delta \left(m/2\pi +\nu \right)$, where the slope m can be, as before, obtained by the process of split plus back-propagation, but now selecting as values for the quad-cell signals the amplitudes of the fields in the position of the selected area. Briefly: the division process incorporates the displacement information of the corresponding delta function in the field without losing the phase needed for back-propagation.

Let us now consider a complex matrix $u_H$ of size $N\times N$ describing the field in a given input plane, which can represent the object reconstructed from a digital hologram.

The procedure begins by computing a Fast Fourier Transform to obtain a discrete version of the angular spectrum of the field

Among the different symmetric paths around the origin of the Fourier plane, the origin of the aforementioned quadrant selection can follow, a square one with side $\delta$was chosen. Then, as Fig. 1(a) depicts, the splicing of $U$begins by inserting it in a larger null-matrix, $U^{\prime }$, of dimensions $M\times M$, where$M=N+\delta$, at the pixel positions $\left(s,t\right)\in \left[0,\delta \right]$ belonging to the squared path shown as a dotted red square in the figure. From $U^{\prime }$, four $M/2\times M/2$ sub-matrices are extracted

 

Fig. 1 (a) Scheme of the data organization; (b) Response,$m\text{'}$, to a constant phase ramp with slope m.

Download Full Size | PDF

It is important to note the limitation in the maximum gradient that can be determined from a sampled complex field. Because of the use of the discrete Fourier transform, the slope must conform to $\left|m\right|\le \pi$ so as not to exceed the Nyquist limit. This is not a limitation of the procedure, but an intrinsic characteristic of the sampling: a digitized phase with slope greater that $\pi$ is always under-sampled.

To evaluate the response of the algorithm, matrices (200x200 px) representing the complex field $u_H$ of constant phase ramps with different slopes were generated. These matrices were inserted into a 1024x1024 px null-matrices prior to compute Eq. (1). The values of $\left(s,t\right)$ correspond to the coordinates of a square path of $\delta =1020\text{px}$side travelled in 10 px steps; the value of $\delta$ was chosen to extend the measurement range to the maximum admissible by the Nyquist limit, while also taking into account the step size. Equation (6) generates two images of 1022x1022 px for the two gradient components. Only the central region of these images (200x200 px) contains useful information, but, in order to avoid border effects, slightly smaller images (190x190 px) are selected to form the matrices of $g_x\left(v_1\right)$ and $g_y\left(v_1\right)$that, in this particular case, the first is constant and the second null. Figure 1(b) shows the value,$m^{\prime }$, of the median of this selection of $g_x\left(v_1\right)$ in the abscises, for the corresponding tilts in one direction in the range $m=\left[0,2.9\right]=\left[0^{\circ },{71.1}^{\circ }\right]$ represented on the ordinate axis. It can be observed that the response is almost linear, only departing slightly from the line at slopes close to the Nyquist limit where the external lobes of the Fourier transform (a shifted delta function convolved with a sinc function) begin to generate aliasing.

3. Obtaining the unwrapped phase in digital holographic microscopy

To demonstrate the proposed phase unwrapping procedure, we used off-axis digital transmission holograms of samples obtained with the system described in references [13] and [14]. Figures 2(a) and (b) show the images (1024x1024 px) of the modulus and phase, respectively, of the object reconstructed from a digital hologram of a sample of human red blood cells (RBC). The reconstruction was done using the angular spectrum approach [3], filtering out the zero diffraction order and discarding the twin image [15]. A quadratic phase term − caused by the non-telecentric (non-afocal) optical system − can be clearly seen.

 

Fig. 2 (a) The modulus and (b) the phase of the reconstructed field of a digital hologram of a sample of human RBCs. (c) The modulus of the Fourier transform corresponding to the region marked I in (a) and the corresponding gradient x-component (upper inset) and y-component (bottom inset). (d) Normalized representation of the phase, corresponding to the equally labeled ROIs in (a). (e) The modulus of the reconstruction of the digital hologram of a slice of a Drosophila head and the corresponding phase, (f), and the integrated phase (g) for the marked ROIs. The color scale is valid for all the contour plots.

Download Full Size | PDF

The unwrapped phase was obtained for different regions of interest (ROIs) of 200x200 px. Each of them was inserted in the centre of a larger null-matrix (1024x1024 px) prior to computing Eq. (1) to increase the sampling of the complex field in the Fourier domain.

Figure 2(c) shows the central region (100x100 px) for the ROI labelled as I in Fig. 2(a) of the matrix corresponding to$\left|U\left(m_2\right)\right|$. It can be seen that the maximum is decentred by the global quadratic phase that can be approximated in the ROI region by an added phase tilt, which was removed by re-centering$U\left(m_2\right)$. Equation (6) was computed using the described algorithm with the same parameters as before; the result is shown in the upper and bottom panels on the right of Fig. 2(c). These data were used to obtain the phase estimation, which, given that the hologram is registered in transmission, corresponds to a map of the optical path difference that the light has suffered traversing the sample. There are several approaches to obtain a surface from the digitized version of its gradient [16,17]. We chose the least squares solution to the Poisson equation with Neumann boundary conditions, which we used previously with good results [10]. The normalized phase map − obtained in 8 milliseconds on a 2.8 GHz Intel Core i7 − is shown in Fig. 2(d) together with the phase distributions for the ROIs labelled as II and III in Fig. 2(a). These results are in good agreement with the expected oval biconcave disk shape of the RBCs.

The holograms in the previous example can be considered regular and the generated wrapped phase an easy target for a conventional unwrapping algorithm. Where our method shows its potentiality is with fields of higher complexity. As an example, we applied it to the object reconstructed from a hologram of a sample consisting of a slice of a drosophila head obtained with the same optical system as above. The modulus is shown in Fig. 2(e). Two ROIs where chosen and independently processed. As can be seen, despite the numerous phase discontinuities in the reconstruction (see Fig. 2(f)), the method provides reliable continuous integrated unwrapped phase maps, which, as Fig. 2(g) shows, can be stitched to display the phase map of a larger region.

6. Conclusions

A method is proposed to obtain the phase from a digital hologram, circumventing the need for unwrapping and so avoiding the problems associated with this difficult process in the presence of noise.

The method is based on computing the gradient of the phase from the angular spectrum of the reconstructed complex field. The procedure allows the propagation of the noise in the hologram data to the noise in the gradient, which, the integration algorithm inherently assumes.

The proposed procedure to determine the gradient is simple, and involves the computation of one direct Fourier transform and several independent inverse Fourier transforms, whose number and dimensions can be optimized and computation parallelized. Also, taking in to account that the reconstruction of the complex field of the object from the digital hologram involves the computation of two Fourier transforms, and that the integration algorithm is fast, the proposed method, once optimized, can provide unwrapped holographic quantitative phase maps in real time.