Phase distribution analysis of tissues based on the off-axis digital holographic hybrid reconstruction algorithm

: Off-axis digital holography (DH) has great potential in histopathology for its high efficiency and precision. Phase distribution, usually extracted by the angular spectrum (AS) algorithm from a digital hologram, reflects important structural information of biological tissues. However, the complex structure of tissues introduces spectrum aliasing of the hologram, making the AS algorithm hard to realize and accurate phase analysis difficult to conduct. Here, we present a hybrid reconstruction algorithm, combining Fresnel reconstruction in spatial domain with the AS algorithm in frequency domain, to solve aliasing by spatial filtering. Through simulation, we demonstrate the feasibility and superiority of the hybrid algorithm and verified the precision (10 − 3 rad) of the hybrid algorithm with spectrum aliasing. We extract phase distributions from normal urothelial and bladder cancer tissues by the hybrid algorithm and make quantitative analysis through histogram and standard deviation. The result shows the hybrid algorithm in off-axis DH has great advantage for the high-precision phase extraction of tissues and supplies significant information for cancer diagnosis.


Introduction
The gold standard of cancer diagnosis is mainly based on manual investigation of hematoxylin and eosin (H&E) stained tissue slides and mainly relies on the experience of pathologists [1,2], making it subject to investigator bias and resulting in low throughput.
Recently, some studies showed there was structural discrepancy between normal and cancer tissues. From biological view, the process of carcinogenesis usually comes with a change in cellular organelles, the extracellular matrix, blood vessels and cancer associate fibroblast [3]. It has been confirmed that the bidirectional interaction between cancer cells and surrounding stromal components is an important contributor to cancer progression and metastasis [4]. Such differences in both cellular and tissue level between cancer and benign indicate the differences in their refractive index (RI) distribution. When light illuminates through the tissue, the RI difference leads to the change of phase distribution. In this way, through phase extraction, we may have an insight of the structure and morphology of tissues and then make cancer diagnosis according to phase analysis [2][3][4][5].
Many optical techniques are investigated to extract the phase distribution, such as optical coherence tomography (OCT) [6][7][8], phase-shifting interferometry [9][10][11][12], digital holography (DH) [13][14][15][16][17][18] and so on. Most of them are based on interferometry and have the ability to obtain phase information with high sensitivity and without contact and stain. However, the OCT needs mechanical scanning, which is time-consuming and precision is low. The phaseshifting interferometry usually depends on phase shifting operation, which gets the phase images quickly via simple mathematical operations but is demanding to the experimental environment and the stability of optical system. DH attracts large attention for its high precise measurement. It was widely used in cells research [19], eardrum tissue study [20,21], microfluidic channel technique [22], cancer diagnosis [23,24], etc. Moreover, the spatial light interference microscopy (SLIM) improved the resolution of phase distribution greatly [25,26], and got primary achievement in breast cancer diagnosis [27]. Most of the studies were based on in-line DH, which has wide field and fully uses recording sensor. However, the in-line DH exists unavoidable image overlapping problem [15,16], which always needs phase shifting to extract object phase distribution. Off-axis DH comes into sight for its high-efficiency. It solves the overlapping problem in in-line DH through introducing an incident angle between object and reference waves and enable to reconstruct object wave from single hologram [15][16][17]28,29].
In off-axis DH, the precision of reconstruction algorithm is important, due to the limited resolution of CCD [29]. Fresnel and AS algorithms are both common reconstruction algorithms in off-axis DH [30,31]. The Fresnel algorithm is usually used to reconstruct object image, while the AS algorithm, an analytical solution of wave equation, has higher precision for phase extraction and thus are widely-used algorithm in biological analysis. The AS algorithm needs spectral selection firstly [15][16][17][18]. Thus, the accurate and complete selection of spectrum has a direct bearing on extracting high-precision phase distribution. The object spectrum is usually selected by band-pass filtering [17,18,32,33]. However, for biological tissues, which have complex structure, the wide spectrum of the tissues will cause serious spectrum aliasing with zero-order diffraction spectrum, which reduces the precision of AS algorithm and make the biological phase analysis not accurate enough. Several methods have been proposed to suppressing zero-order and thus eliminate the spectrum aliasing, such as nonlinear reconstruction technique [34], iterative approach [35-37] and image plane filtering method based on spherical reconstructed wave [38,39]. These methods can solve the spectrum aliasing, but the nonlinear technique may have complex processing, the iterative approach is comparatively time-consuming and the spherical wave introduced in image plane filtering is not suitable for phase extraction because of the introduction of other phase information.
In this paper, to extract high-precision phase distribution of biological tissues, we present a hybrid algorithm to solve spectrum aliasing in off-axis digital hologram. The algorithm takes advantage of separable images of Fresnel amplitude reconstruction, achieves complete object spectrum through spatial filtering and back calculation, and finally extracts highprecision phase distribution by AS algorithm. Through simulation, we demonstrate the feasibility and superiority of hybrid algorithm and verify that the precision of hybrid algorithm with spectrum aliasing is 10 −3 rad, which is equal to the precision of AS algorithm without spectrum aliasing. Based on the algorithm, we extract high-precision phase distributions from normal urothelial and bladder cancer tissues and make quantitative analysis through phase histogram and standard deviation. The result shows that the hybrid algorithm in off-axis DH has advantage in phase extraction of biological tissues and supplies significant information for cancer diagnosis.

Tissue preparation
Tissue samples were obtained from 10 bladder cancer patients undergoing radical cystectomy at Ren Ji Hospital, Shanghai Jiao Tong University School of Medicine, following an institutional review board approval. Two samples were taken from each patient, one of which was normal urothelium, and the other was pure cancer tissue. The tissue used in this procedure was embedded in paraffin and cut on a Leica RM2255 microtome at 5 μm thickness. The sections were placed on a water bath at 38 °C. Each section was collected on a separate glass slide and heat fixed in a 70 °C slide drying oven for 15 min. After experiment, the samples are stained by H&E for pathological confirmation.

Experimental system
The modified Mach-Zehnder interferometer system used in experiment is shown in Fig. 1. The light source is MSL-FN-532nm single longitudinal mode laser. The resolution of CCD is 2452 × 2054 (M × N) pixels. The pixel pitch (Δx 0 × Δy 0 ) of CCD is 3.45 μm × 3.45 μm. In optical system, the output of laser is attenuated to an appropriate intensity through attenuator 1 and split into two beams of light by beam splitter 1(BS1). One beam passes through objects, utilized as object wave, carrying structural information of the sample. The other beam acts as the plane reference wave. It is attenuated to an appropriate intensity by attenuator 2 and then reflected by mirror 1 (M1) and M2. The object wave and reference wave are amplified by 40 × (0.65 NA) microscope objectives 1 (MO1) and MO2 separately. The recording distance d is from the output of MO1 to the CCD recording plane. Through adjusting BS2, we get appropriate incident angle θ between the object and reference wave, which interfere on the CCD sensor interface and generate the off-axis digital hologram, saved in the computer for further processing. In the experiment, we extract phase distributions of holograms with and without object information and eliminate additional phase noise through phase subtracting. So the phase result we get is the phase difference between the tissue and paraffin in fact [17, 18].

Digital holography
Holography is mainly based on optical interference between object and reference wave, turning the phase information of object wave into intensity one. The digital hologram I = I(mΔx 0 , nΔy 0 ) is a discrete interference figure recorded by CCD: where O is object wave, R is reference wave, Δx 0 and Δy 0 denote the sampling intervals in recording plane, and M and N are equal to the horizontal and vertical pixel number of CCD. Then, we generate an illuminating wave C, the same as R in this paper, and restore the wavefront U in the CCD recording plane through multiplying the hologram by C: The equation shows that the digital hologram contains three parts: zero-order diffraction |O| 2 × R + |R| 2 × O, object image RR*O and conjugation image RRO*. The zero-order diffraction part contains intensity information of object and reference wave. This part is nothing to do with phase extraction and thus is the main noise of digital hologram. RR*O contains useful object information we need, and RRO* is the conjugation of object image.

Off-axis DH reconstruction algorithm
AS algorithm, an analytical solution of wave equation, simulates the process of diffraction in frequency domain. Spectral selection is the first step of the algorithm to get complete object spectrum. Corresponding to the three parts of digital hologram in Eq. (2), the spectrum of hologram also has three parts, whose location are separated in off-axis DH theory [13] and related to the incidence angle θ between object and reference wave and wavelength λ. To satisfy the Nyquist sampling theorem and to avoid spectrum aliasing between object and zeroorder, the value of θ must be limited within the following range [15,16].
where B is the bandwidth of object spectrum, mainly related to the structure of object. The limitation of θ is usually within four degrees in general, which is difficult to control in experiments. Furthermore, when the object has complex structure, making the bandwidth B wide and the value of θ min close to or even more than the value of θ max , the spectrum aliasing is unavoidable, leading to serious precision reduction of AS algorithm. Fresnel algorithm is also a common method for reconstruction in spatial domain. It efficiently simulates the process of diffraction by only once Fourier transform through paraxial approximation. However, because of the approximation, there is scale variation problem in the algorithm. Equation (4) shows the relation between the sampling intervals Δx and Δy in Fresnel reconstruction plane and the sampling intervals Δx 0 and Δy 0 in CCD plane: where z is reconstructed distance, which is the same as recording distance d. The character changes the size of Fresnel reconstructed image and thus decreases the pixels to record the object image, making Fresnel algorithm not suitable to phase extraction. However, the Fresnel reconstructed image contains complete information of object wave and the image size can be controlled by z, which has nothing to do with θ. If we choose a suitable z to control the size of Fresnel reconstructed image, separating the three parts exactly and only selecting the object image we need in spatial domain, we can get complete object spectrum based on the selected object image. In this way, we solve spectrum aliasing by filtering in spatial domain and thus give full play to the high-precision advantage of AS algorithm.

Hybrid reconstruction algorithm
According to the idea said above, we propose a hybrid algorithm, combining Fresnel amplitude reconstruction in spatial domain with AS algorithm in frequency domain, to solve the spectrum aliasing of digital hologram and thus improve the precision of biological tissues phase extraction in off-axis DH.
The basic procedure contains Fresnel amplitude reconstruction, spatial filtering, back calculation and AS algorithm, as shown in Fig. 2. Based on off-axis digital hologram, both Fresnel amplitude reconstruction and spatial filtering are performed in spatial domain, as framed in the figure, to achieve complete object image. After that, back calculation turns the spatial domain into frequency domain, getting complete object spectrum, and then AS algorithm is used to extract phase distributions of samples. The following is the mathematical expression and analysis of hybrid algorithm. In spatial domain, reconstructed image is needed rather than its final phase distribution, so we only use a part of Fresnel algorithm, called Fresnel amplitude reconstruction, to get spatial reconstructed image U 0 (kΔx, lΔy), as shown in Eq.
where FFT{} is fast Fourier transform. According to the separable images in U 0 (kΔx, lΔy), we select object image of it and define the filtered image as U 0f (kΔx, lΔy). As this selection process is done in spatial domain, we define the operation as spatial filtering. After that, the back calculation is carried out, turning the spatial domain into frequency domain and getting complete object spectrum V(pΔξ, qΔη) as following: where Δξ and Δη denote the sampling intervals in frequency domain. Then, the traditional AS algorithm is used to get object wave U(mΔx 0 , nΔy 0 ): ( ) Equation (8) shows the phase distribution Δφ(mΔx 0 , nΔy 0 ), which is between (-π, π]. After least-square unwrapping, we can get the real phase distribution of object wave:

Simulation
In this section, numerical simulation is performed to verify the feasibility and superiority of hybrid algorithm in biological tissues phase extraction. First, we generate phase and amplitude distributions based on a biological tissue to simulate practical situation, as shown in Fig. 3(a) and (b). The incidence angle θ is set to four degrees. The recording distance d is set to 1.5 m, as well as the reconstructed distance z to separate the reconstructed image exactly in spatial domain. Other parameters are assumed according to the real parameters, as described in experiment section. After simulating interference, we get the simulated off-axis hologram, shown in Fig. 3(c) and the image in top right corner is partial enlarge image framed in the figure.

Phase extraction without spectrum aliasing
To verify the feasibility of hybrid algorithm, we consider to get a digital hologram without spectrum aliasing and compare the phase distributions extracted by hybrid and AS algorithm. Figure 4(a) shows the spectrum of original simulated hologram. Obviously, the three parts of the spectrum, marked in the figure, exists serious aliasing. Figure 4(b) is the spectrum of zeroorder diffraction, which can be separated from hologram in simulation. We remove the zeroorder diffraction spectrum and thus get the new hologram shown in Fig. 4(c). Figure 4(d) and Fig. 4(e) is the spectrum and Fresnel reconstructed image of the new hologram respectively. After spatial filtering (framed in Fig. 4(e)) and back calculation, we get the complete object spectrum, as shown in Fig. 4(f). We choose object spectrum (framed in Fig. 4(d)) of the new hologram and extract its twodimension (2D) phase distribution by AS algorithm, as shown in Fig. 5(a). Meanwhile, based on the spectrum in Fig. 4(f), we extract the phase distribution by hybrid algorithm, as shown in Fig. 5(b). Then, the phase distributions of these two algorithms are compared with simulated phase distribution. In Fig. 5(c), we choose phase curves in the 800th column of the phase distributions (drawn in Fig. 5(a), (b)) and compare them with the same column of original phase distribution in Fig. 3(a). In addition, we subtract phase distribution extracted by hybrid algorithm from original one in Fig. 3(a) to get the hybrid algorithm error value distribution, as shown in Fig. 6(a), and the same operation is done to get the AS algorithm error value distribution shown in Fig. 6(b). Figure 6(c) shows the histogram associated with the error value distributions of two algorithms. In the histogram, the error values mainly concentrate between −5 × 10 −3 and 5 × 10 −3 rad, which means the extracted phase distributions are consistent with original one and the precisions of both algorithms achieve 10 −3 rad. Thus, to off-axis digital hologram without spectrum aliasing, the precision of hybrid algorithm is similar to the one of traditional AS algorithm and thus the hybrid algorithm is a feasible method for phase extraction.

Phase extraction with spectrum aliasing
In this part, we compare hybrid and AS algorithm in the condition of spectrum aliasing. Because of serious spectrum aliasing, it is hard for AS algorithm to select the spectrum precisely. However, to make comparison, we choose different spectral ranges framed by yellow and red box in Fig. 7(a). The yellow one tries to select complete object spectrum, while the red one avoid the zero-order as possible. Then, the corresponding phase distributions are extracted by AS algorithm, as shown in Fig. 8(a) and (b). In hybrid algorithm, we get complete object spectrum, as shown in Fig. 7(c), through spatial filtering (framed in Fig. 7(b)) and back calculation. Figure 8(c) is the phase distribution extracted by hybrid algorithm. Figure 8(d) is the phase curve in the 800th column of 2D phase distributions (drawn in Fig. 8(a), (b) and (c)) and compare them with the same column of the phase distribution in Fig. 3(a). Obviously, the phase distribution extracted by hybrid algorithm is still agreement with simulated one, but there are serious noises in AS algorithm 1 and much information loss in AS algorithm 2.   Figure 9(d) is the histogram of error values. According to the histogram, the precision of AS algorithm reduces seriously, while the precision of hybrid algorithm still keeps within 10 −3 rad. The result indicates that the hybrid algorithm has superiority in extracting phase information from the hologram with spectrum aliasing. The phase distribution extracted by hybrid algorithm contains complete information and shows more details of biological tissues, making the phase analysis based on off-axis DH more accurate.

Experiments and results
In this section, we extract high-precision phase distributions from normal urothelial and bladder cancer tissues by hybrid algorithm and make phase analysis based on them.  Figure 10(b) is the hologram spectrum with spectrum aliasing. Figure 10(c) is Fresnel amplitude reconstructed image of hybrid algorithm. After spatial filtering and back calculation, we get complete spectrum in Fig. 10(d). Similarly, we choose different ranges (framed by yellow and red box in Fig. 10(b)) of the hologram spectrum, and extract their phase distributions respectively by AS algorithm, as shown in Fig.  11(a) and (b). Figure 11(c) is the phase distribution extracted by hybrid algorithm based on the spectrum in Fig. 10(d). Figure 11(d)-(f) are the 3D display maps of Fig. 11(a)-(c).  In our optical system, single digital hologram shows 100 × 100 μm 2 region of one sample. To achieve wider field of view (2 × 2 mm 2 ), we get 400 holograms by scanning and mosaic them together. Figure 12(a) and (d) show H&E stained images of normal urothelial and bladder cancer tissues from one patient. Figure 12(b) and (e) are the phase distribution of them with the same areas. Figure 12(d) and (f) are local zooms of Fig. 12(b) and (e). Obviously, the phase distribution of bladder cancer tissue is much more dispersive without well-defined structure than that of normal one and has lower phase value in general than the normal one. Regarding the pathological features, urothelial cancers differ from normal one by having a loss of cell polarity, increased prominent nucleoli and clumping of chromatin [3], which may cause these discrepancies of their phase distributions. In order to quantitatively analyze the information contained in phase distributions for normal urothelial and bladder cancer tissues, we present histograms (in Fig. 13(a) and (b)) associated with phase distributions of Fig. 12(b) and (e). In histograms, we divide each phase distribution into 10 regions in transverse direction and thus each curve represents statistical data of one region. The horizontal ordinate shows the phase distribution, while the vertical ordinate shows the statistic pixels of each phase value. Obviously, the phase distributions are consistent within the normal (cancer) tissues but different between normal urothelial and bladder cancer tissues. The distributions of normal urothelial tissue is wider than that of cancer tissue. Based on phase distributions, we calculate standard deviation of 20 samples from 10 patients. The standard deviation shows the degree of difference between the phase distribution and its average value. Green and red columns in each group represents the data of normal urothelial and bladder cancer tissues from one patient, as shown in Fig. 13(c). It is shown that the standard deviation values of cancer tissues are generally lower than that of normal urothelial tissues. Based on existing samples, the threshold of standard deviation between cancer and normal tissues is around 1.5.

Conclusion
Phase distribution has an insight to the structural and morphological information of biological tissues and thus becomes research hotspot in histopathology. The accuracy of phase analysis in histopathology is directly linked to the precision of phase distribution. In this paper, to extract high-precision phase distributions from biological tissues, we present a hybrid reconstruction algorithm in off-axis DH. The algorithm combines Fresnel amplitude reconstruction in spatial domain with AS algorithm in frequency domain to solve spectrum aliasing in off-axis digital hologram. Through simulation, we demonstrate the feasibility and superiority of hybrid algorithm in biological tissues phase extraction. The result shows that the precision of hybrid algorithm with spectrum aliasing is 10 −3 rad, which is equal to the precision of AS algorithm without spectrum aliasing. We extract high-precision phase distribution of normal and tumor bladder tissues by hybrid algorithm, make quantitative analysis through histogram and standard deviation and find the phase discrepancy of normal and tumor tissues.
The nature of modern pathology could be concluded as that the biological information of tissue samples is turned into optical information by staining, and then captured and analyzed by human. In our study, this procedure is conducted by a modified machine system, by using which the benign tissues can be distinguished from malignancies automatically in computer. With the improvement of phase extraction precision, the off-axis DH using hybrid algorithm is expected to be an effective method for gaining more detailed and quantitative data from biological tissues. In the near future, validation study in a larger sample population will be made, which has a certain potential to perform or assist digital pathological diagnosis.