Automated Fourier space region-recognition filtering for off-axis digital holographic microscopy

Automated label-free quantitative imaging of biological samples can greatly benefit high throughput diseases diagnosis. Digital holographic microscopy (DHM) is a powerful quantitative label-free imaging tool that retrieves structural details of cellular samples non-invasively. In off-axis DHM, a proper spatial filtering window in Fourier space is crucial to the quality of reconstructed phase image. Here we describe a region-recognition approach that combines shape recognition with an iterative thresholding to extracts the optimal shape of frequency components. The region recognition technique offers fully automated adaptive filtering that can operate with a variety of samples and imaging conditions. When imaging through optically scattering biological hydrogel matrix, the technique surpasses previous histogram thresholding techniques without requiring any manual intervention. Finally, we automate the extraction of the statistical difference of optical height between malaria parasite infected and uninfected red blood cells. The method described here pave way to greater autonomy in automated DHM imaging for imaging live cell in thick cell cultures.


Introduction
Phase contrast microscopy [1] enables the transformation of phase into intensity that allows clear observation of cellular boundaries, but remains a qualitative imaging tool. The emergence of quantitative phase imaging aims to extract information of the optical height of biological samples. Off-axis holography [2] operates by projecting two equal copies of light into two separate paths, one undisturbed and the other perturbed (phase delay) by a sample or an object. The two paths combine to reveal the amount of phase delay that is recorded on a photosensitive film. Phase delay represents the optical path length difference between the two copies of the beams, which is then used to retrieve the optical height of the object. Combining holography with microscopy and digital imaging devices (charge coupled or complementary metal-oxide), Digital holographic microscope [3,4] records interference pattern (hologram) that emerges from superposition of a reference wavefront with that emerging from the microscopic sample. Both optical height and refractive index fluctuations are then used to quantify differences between biological samples [5][6][7][8][9]. DHM has been utilized in many biological applications such as determination of cellular motility [5,10], morphology[11] and biomechanics [12][13][14].
DHM heavily relies on numerical calculations to accurately reconstruct the phase of the biological sample. The numerical steps used in off-axis DHM technique include: discrete 2D Fourier transform, spatial frequency filtering [15], numerical propagation, phase unwrapping and aberration correction [16]. To achieve high throughput diseases diagnosis, these numerical calculations should be programmed to be automated and adaptive to different sample and imaging conditions. The complexity lies in the spatial frequency filtering process. Since the distribution of spatial frequencies of different samples varies significantly, the filter window needs to change correspondingly for the precise extraction of the desired frequency orders (first orders). One straightforward filtering method for precise extraction is manually selecting the first order [16] in spatial frequency domain. Obviously, this approach requires manual input for each hologram, which is therefore not suitable for automation. While there are some histogram analysis techniques [17,18] aiming to provide adaptive fitering that could possibly provide automated phase retrieval. Unfortuately, the histogram analysis process still require manual intervention (e.g. initial windowing parameter) to calculate the appropriate threshold. Hence, these techniques is still unable to provide sufficient flexibility or simplicity for automated selection of the appropriate spatial filter for a wide variety of samples especially confluent live cell cultures with highly scattering backgrounds.
Here we propose an adaptive filtering process based on iterative thresholding and regionbased selection. This combination gradually selects the optimal frequency component boundary and using shape recognition to extact the optimal frequency component for different holograms. We find that this approach firstly avoids fixing a single thresholding level and removes the need for any manual intervention and any initial input parameter. In addition, the automated adaptive thresholding approach operates successfully with a range of samples including scattering samples when compared to previous techniques[17, 18] (histogrambased). The region recognition approach also demonstrates to be robust for all standard experimental conditions i.e. field of view, sample density, orientation and spacing of interference fringes. The automated DHM system was showed to extract quantitative phase measurements of biological red blood cells (RBCs) infected and uninfected with malaria parasite. All the processes are conducted in a single loop, which allows the continuous imaging and controlled use of the user interface. In the following sections, we explain the steps taken to implement the region recognition and iterative thresholding process as well as the automated sample classification of RBCs.

Method
Interferometry measures the phase of cross-correlations between the sample and reference wave, shown as following and R(r , t) are wave functions of object wave and reference wave,  indicates the angular frequency of both waves and kr  indicates the dot product of wave vector k ( two waves interfere with each other with same incident angle) and optical length r .  is the phase delay produced by the sample. The phase change derived from the sample (O) that is recorded as a shift in the intensity pattern of fringes I, Based on Eq. (4), the phase shift is recorded in two spatial frequency components  that are spatially overlapping in an inline holography. In off-axis system, the sample and reference beams interfere at a slightly different angle  creating two separate wave vector , ro kk in two waves: Hence the interference pattern is rewritten as: The phase shifts are spatially separated ()   5)] is collimated and then expanded by lenses L1 (f = 30mm) and L2 (f = 100mm) and finally recombined with the sample beam by a non-polarizing beam splitter (BS, THORLABS CM-BS013) onto a charge coupled device (CCD) camera (Retiga-2000R Fast 1394 Color Cooled) by tube lens L8 built in microscope frame. Figure  1(b) shows experimental setup where the sample is placed in an inverted microscope system. The incident directions of reference and sample beams at the CCD plane are indicated by r k and o k in Fig. 1(a). o k is added at a slight tilted angle  which results in off-axis holograms digitally recorded as shown in Fig. 1(a) beside the CCD. After fast Fourier transform (FFT) of the recorded hologram [ Fig. 2(a)], the phase shift is recorded in the spatial frequency on two symmetrical areas (+1 and -1 terms) in the frequency domain (FFT hologram). The shapes and sizes of these two symmetrical areas can vary according to different imaging conditions. Figure 4 shows some examples that different samples and imaging conditions give different symmetrical frequency components outlines. This implies that the frequency component needs to be carefully selected for accurate reconstruction. Here we propose an adaptive filtering process that is robust to all standard experimental conditions. The method is built upon an iterative thresholding and region-based selection. The procedures is broken down into 3 steps:

Region-recognition of spatial frequency
(1) Apply global threshold level (GTL)[21] to the intensity of FFT hologram to get the binary image and then implement the region recognition process of it (regionprops MATLAB is used to provide three properties i.e. size of objects, centroid of objects in binary image and the box boundaries; graythresh MATLAB is used to get the GTL). (2) Increase the threshold level by one percent of GTL and repeat the first step until the number of regions reaches three (one for background intensity, one for the desired frequency component and one for the mirrored frequency component).
(3) Use box boundary data from regionprops function and the binary image form first step to get the right frequency component boundary and use it as a filtering window (additional Gaussian function is applied to smooth the edge of final filtering window). Figure 2 shows an example of the adaptive filtering process. Figure 2(a) shows the recorded hologram of RBCs and Fig. 2(b) is the spatial frequency after FFT. Figures 2(c), (d) and (e) show the iterative threshold process with incremental threshold level. The processing time of this method is shown in Fig. 2(i). Once the appropriate spatial frequency is selected and centered [ Fig. 2(g)], an inverse FFT (iFFT) is performed to retrieve the complex amplitude of the sample, which can be seen from Fig. 2(h). Figure 2(j) shows the flow chart of iterative thresholding and region-recognition program. Since the phase map acquired from the recorded hologram is limited from −π to π, so called wrapped, the true continuous phase value corresponding to the optical height of cells needs to be unwrapped [22]. Then we apply Zernike polynomials to compensate the aberration [23] in the system.
There have been previous techniques [17, 18] aimed at providing adaptive filtering methods by using histogram analysis to settle the appropriate threshold. However, these methods entail a series of steps that require manual settings or prior knowledge of the frequency distribution. For example, the first procedure in pervious technique [18] is to apply a negative Laplacian [24] to suppress zero order in FFT hologram, which requires distance between first and zero order frequency components to be known prior or well-defined in the FFT hologram. Intuitively, the approach may not operate well when positions of the spatial frequencies change (i.e. different samples and interference fringes on the CCD camera or that the samples has large optical scattering).

FIG. 3. Comparison of histogram analysis filtering technique [18] with region-recognition method
under different sample conditions. (a), (i) and (ii) are the holograms of a 6 µm diameter microsphere through water and two 50 µm diameter microspheres through opaque peptide hydrogel. (iii) and (iv) are the FFT holograms of (i) and (ii) respectively. In (iii) and (iv), the red outlines represent the best filters from region recognition method and blue circles represent manual pre-cropping area used for histogram analysis. Reconstructed phase image with histogram analysis are shown in (b) and region recognition method (c) .Color bar h1 = 6.7 µm Color bar h2 = 50 µm Color bar h3 = 6.7 µm, Color bar h4 = 50 µm The caveat of previous histogram analysis technique [18] is the requirement to first manually pre-cropping the area containing the first order frequency. The histogram of the area would then be analysed accordingly (shown in more detail in Appendix A1). The drawback of using such manual pre-cropping process is that it influences the quality of the final reconstruction from sample to sample. We illustrate this point by directly applying the same pre-cropping and histogram analysis to retrieve phase images of polymer microspheres in a transparent medium into a scattering medium. The reconstructed phase is then compared with the iterative threshold technique.
The purpose of this study is to prepare an automated DHM system that is suitable to study living cells develops into a living tissue. In thick samples, such as thick tissue culture, the effect of scattering is more pronounced and hence the reconstruction becomes more complex due to higher frequency noises. We therefore test the thresholding process on their performance sample with thick (~ 100 µm) optically scattering biological-relevant matrix (peptide-hydrogel). The scattering medium (Appendix Fig. A1) is prepared with 50 µm polymer microspheres (Duke Scientific) deposited onto a thin film of biological hydrogel matrix peptide [25] hydrogel (Fmoc-FRGDF). The peptide hydrogel has been used in thick cell growth culture for neuronal growth. The hologram of a polymer microsphere in water, Fig. 3.(a)(i), is significantly different from the hologram of hydrogel matrix [ Fig. 3.(a)(ii)], as the latter displays a large background scattering with significant reduction in opacity. The scattering sample generates additional spatial frequency content that "blurs" the boundaries and "extends" the dimensions of three spatial frequencies components [FIG. 3.(a)(iv)]. That adds complexity to precise selection of first order frequency content for phase retrieval process. The results of the reconstructed phase maps from the histogram analysis technique, Fig 3. (b), and our approach, Fig. 3(c), show that the histogram analysis approach works in a standard transparent medium Fig. 3(b) (i). but clearly fails to extract the appropriate spatial frequency in the scattering medium as shown in Fig. 3 (b) (ii). The iterative methodworks for both samples as shown in Fig. 3(c) (i) and (ii). The comparison indicates that the filtering process is complicated by the broad distribution of spatial frequency. Previous histogram analysis techinique failed to obtain the desired reconstruction. While our approach, iterative thresholding and region recognition, successfully retrieved optimal shape of the filtering windows [blue curves in Fig. 3(a)(iii) and (iv)] and obtained accurate reconstruction of the microsphere through a diffusive medium [ Fig. 3(c)]. A more detailed breakdown of the steps of both the histogram analysis technique and threshold increment determination are shown in Appendix A1. We also include another histogram analysis [17] in Appendix A2 that also fails to retrieve the spatial frequency of scattering medium. After unwrapping and compensation, the remaining operation is to calculate the optical height of the sample based on wavelength 0  and the refractive index difference ∆n between the sample (polymer sphere: 1.59 or cell: 1.36) and the media (air: 1, water: 1.33 or oil: 1.51). The height differences within a sample can be plotted as a map using phase value φ(x,y).
xy h x y n     (8) To demonstrate the adaptability of our system to a range of experimental conditions we tested it using murine primary cortical neurons (cells that are adhesive, polymorphic, and relatively large), human RBCs (cells in suspension, biconcave and small) and microspheres (polystyrene, uniform structures of intermediate size) . Figures 4(a), (b)  These results show that the iterative thresholding and region-recognition method can process a diverse range of shapes of spatial frequencies.
The entire numerical processing is implemented in a MATLAB program with a graphical user interface (GUI) as shown in Fig. 5(a). In the GUI, the hologram will be displayed at the right side of control panel with variable settings on gain and exposure time (based on camera settings). A region of interest is selected, cropped and displayed. The phase reconstruction is shown with two views: perspectives and plan.

Region-recognition for classification of infected RBCs
Analysis of cellular characteristics are often restricted to one of two scenarios: either averages of cell populations are examined (missing any heterogeneity that might be present within the population) or individual cells are studied (with the disadvantage of limited numbers and restricted in-vivo compatibility). Furthermore, most analytical imaging methods involve the addition of labels (e.g. fluorescent tags, antibodies, chemical dyes), which potentially perturb the examination. Our system facilitates the classification of cellular populations without the aid of labels rapidly by retrieving phase values for individual cells. This is critical for various biological applications. For example, RBCs are the most abundant cell type in the human body [26]. Mutations in RBCs are very common and some have serious implications for the wellbeing of the carrier (e.g. sickle cell anemia, beta thalassemia, ovaloycytosis to name only a few [27]). Moreover, infections (e.g. malaria, babesiosis [28]) and drugs (e.g. iodinated contrast media [29]) can also influence the characteristics of RBCs and the analysis of the characteristics is important for diagnostics, disease outcome and research.
Here the region-recognition program is expanded toward characterization of RBCs. Fig.6 illustrates the key steps (thresholding, perimeter location) used to extract individual RBCs. The first step is to read image and apply a Gaussian filter for thresholding (Otsu's method [30]). Next, holes are filled to remove small objects and extract perimeter pixels. The third step is to get the center regions of individual cells and do the same processing of second step with these center regions. The fourth step is to inverse image intensity. The fifth step is to separate objects using watershed transform and measure the pixel areas of each object in the combined binary image. The final step is to label cell regions and remove false cell detections based on the pixel areas of objects. The phase value of each identified cell is saved into arrays for further analysis.  N.A 1.25). The color bar beside shows the phase values by unit of radians. Seen from the figures, the phase maps present the consistency with the SEM pictures, where infected RBC has anomalous shape while uninfected RBC shows normal biconcave shape.
In order to pick out the morphological differences between infected and uninfected RBCs, the phase value over the cell surfaces is selected out to plot the histograms. The black dots in Figs. 8(a), (b) are the histograms of normalized phase values over RBCs shown as thumbnails respectively. The histogram indicates how height of a cell distributes over the 2D plane. The solid curves in Figs. 8(a), (b) are Gaussian fittings of the histograms. Quantification of the phase is inferred from the full width at half-maximum (FWHM) of the fitted curve of the phase distribution provides an indication of the spatial distribution of the phase (height of a cell) variation over an entire cell as shown in Fig. 8 (a), (b) shaded region. Figure 8(c) shows the distribution of FWHMs of the spatial distribution of the phases on 12 infected and 12 uninfected RBCs. The trends indicate a marked difference between infected and uninfected RBCs. The infected RBCs is most likely have lower FWHMs (the histogram is more centralized), which means the differences in cell height are not major. Therefore, the blue line in Fig. 8(c) is the histogram of FWHMs of infected RBCs. It shows most of infected RBCs have FWHMs at below 0.1, while uninfected cells (red line) tend to have higher FWHMs (0.25 to 0.3). This shift indicates that the height of uninfected RBCs changes largely spatially. In other words, the uninfected RBCs have broader height (phase) distribution. It can be explained by the fact that the uninfected RBCs have thinner central height. Therefore, the height difference between the centre and the periphery is larger than that of infected RBCs where the volume of the centre part is occupied by the parasite.

Conclusion
DHM relies heavily on accurate execution of numerical processes to reconstruct the phase of the biological sample. While off-axis DHM technique is perhaps the fastest holographic imaging technique, the accurate filtering of the spatial frequency can affect the final reconstructed image. We demonstrated that the combination of region-recognition and iterative thresholding conducted in a single loop can aid in selecting the optimal spatial frequency components with marginal delay (~0.1 to 0.3s). The effectiveness of our process is illustrated in the comparison study where the iterative approach clearly identifies reconstructed microspheres embedded in opaque hydrogel matrix (~ 100 µm) has addressed the use of the technique in thick samples. More importantly, as compared with past techniques, we showed that this technique can operate well in both transparent and turbid system, thus providing potential flexibility for imaging live cell activites [32]. In the next step, we will aim to conduct live cell imaging in confluent cell cultures that will aid in the study of cellular process in engineered tissue.