Propagation phasor approach for holographic image reconstruction

To achieve high-resolution and wide field-of-view, digital holographic imaging techniques need to tackle two major challenges: phase recovery and spatial undersampling. Previously, these challenges were separately addressed using phase retrieval and pixel super-resolution algorithms, which utilize the diversity of different imaging parameters. Although existing holographic imaging methods can achieve large space-bandwidth-products by performing pixel super-resolution and phase retrieval sequentially, they require large amounts of data, which might be a limitation in high-speed or cost-effective imaging applications. Here we report a propagation phasor approach, which for the first time combines phase retrieval and pixel super-resolution into a unified mathematical framework and enables the synthesis of new holographic image reconstruction methods with significantly improved data efficiency. In this approach, twin image and spatial aliasing signals, along with other digital artifacts, are interpreted as noise terms that are modulated by phasors that analytically depend on the lateral displacement between hologram and sensor planes, sample-to-sensor distance, wavelength, and the illumination angle. Compared to previous holographic reconstruction techniques, this new framework results in five- to seven-fold reduced number of raw measurements, while still achieving a competitive resolution and space-bandwidth-product. We also demonstrated the success of this approach by imaging biological specimens including Papanicolaou and blood smears.

illumination angles 33 , or illumination wavelengths [45][46][47][48] . Each one of these holograms essentially serve as independent physical constraints on the amplitude of the optical field, which enables the use of an iterative algorithm to force the complex object field to be consistent with all these measurements 36,37,44,49,50 . Although this sequential implementation of pixel super-resolution followed by phase retrieval has enabled digital holographic microscopy to deliver high-resolution and wide-field reconstructions with giga-pixel level throughput, they currently require large amounts of holographic data. For instance, in a multi-height configuration (i.e., using multiple sample-to-sensor distances) 25,36,37,44 , if 4 × 4 pixel super-resolution is implemented at eight different heights, the total number of raw holograms to be captured becomes 128, which could be a limitation for e.g., high-speed imaging applications.
Here, we present a new computational method, termed as propagation phasor approach, which for the first time, combines pixel super-resolution and phase retrieval techniques into a unified mathematical framework, and enables new holographic image reconstruction methods with significantly improved data efficiency, i.e., using much less number of raw measurements to obtain high-resolution and wide-field reconstructions of the specimen. Based on our analytical derivations, the twin image noise and spatial aliasing signals, along with other digital holographic artifacts, can be interpreted as noise terms modulated by digital phasors, which are all analytical functions of the imaging parameters including e.g., the lateral displacement between the hologram and the sensor array planes, sample-to-sensor distance, illumination wavelength, and the angle of incidence. Based on this new propagation phasor approach, we devised a two-stage holographic image reconstruction algorithm that merges phase retrieval and pixel super-resolution into the same unified framework. Compared to previous holographic reconstruction algorithms, our new method reduces the number of raw measurements by five to seven fold, while at the same time achieving a competitive spatial resolution across a large field-of-view.
Based on the same propagation phasor framework, we also created two new digital methods to achieve pixel super-resolution using (1) the diversity of the sample-to-sensor distance (i.e., multi-height based pixel super-resolution), and (2) the diversity of the illumination angle (i.e., multi-angle based pixel super-resolution). We demonstrated the success of these methods by imaging biological specimens such as Papanicolaou and blood smears. We believe that with its significantly improved data efficiency, this new propagation phasor based approach could be broadly applicable to increase the space-bandwidth-product of various digital holographic microscopy systems.

Methods
Optical setup for holographic imaging. Figure 1a depicts our configuration of an in-line holographic imaging system: the coherent or partially coherent incident light first impinges on the specimen, the directly-transmitted light and the scattered light then interfere at an image sensor chip, which samples and digitizes the intensity of this interference pattern. To demonstrate our propagation phasor approach for holographic image reconstruction, we selected to implement it using lensfree holographic microscopy although it is broadly Figure 1. (a) Configuration of an in-line holographic imaging system. Some of the controllable parameters of the imaging system are marked in blue color, including the illumination angle (θ k and ϕ k ), wavelength λ k , sample-to-sensor distance z k , and the lateral displacements between the hologram and the image sensor planes (x shift,k and y shift,k ). (b) Schematic of the optical setup of a lensfree on-chip holographic microscope. applicable to other holographic microscopy platforms. As depicted in Fig. 1b, our lensfree holographic microscope includes three parts: a fiber-coupled wavelength-tunable light source (WhiteLase-Micro, model VIS, Fianium Ltd, Southampton, UK), an image sensor chip (IU081, Sony Corporation, Japan), and a thin specimen mounted above the sensor chip. The optical fiber's outlet is placed at e.g. ~10 cm away from the sample whereas the sample-to-sensor distance is typically 0.1-1 mm and thus the illumination at the object plane can be considered as a plane wave. By bringing sample close (sub-mm) to an image sensor chip, lensfree on-chip holography allows the utilization of the image sensor active area as the object field-of-view, creating a unit magnification in-line holographic imaging system, where the spatial resolution and field-of-view can be independently controlled and adjusted by the pixel design and the number of pixels, respectively 38 . The fiber optic cable is mounted on a rotational arm (PRM1Z8, Thorlabs, New Jersey, USA) that can move across a dome above the specimen so that the incidence light can also be adjusted to an arbitrary angle. The rotational arm is loaded on a mechanical linear stage that moves in lateral directions to introduce sub-pixel displacements between the hologram and the image sensor-array. The specimen is held by a piezo-driven positioning stage (MAX606, Thorlabs, New Jersey, USA), which can move vertically to change the distance between the sample and the image sensor chip. During the holographic data acquisition, the tunable source, the mechanical stages, and the image sensor chip are all automated and coordinated by a PC running a custom-written LabVIEW program (Version 2011, National Instruments, Texas, USA).
Sample preparation. Besides a standard 1951 USAF resolution test target, we also demonstrated the success of our propagation phasor approach by imaging biological samples, including unstained Papanicolaou (Pap) smear slides and blood smears. For this purpose, we used existing and anonymous specimen, where any subject related information cannot be retrieved. Pap smears are prepared using ThinPrep ® method (Hologic, Massachusetts, USA). The blood smear samples are prepared using EDTA (ethylenediaminetetraacetic acid) anticoagulated human blood and stained with Wright's Stain 51 .
Mathematical formalism of propagation phasor approach in digital holography. In this sub-section we present the concept of propagation phasors by deriving the analytical expressions that contain not only the holographic information of the specimen, but also the twin image noise, spatial aliasing signal, and upsampling related spatial artifacts. In this manuscript, we use lower case letters to represent the functions in spatial domain, and the upper case letters for functions in spatial frequency domain. Throughout our analysis, we assume a plane wave illumination as also supported by our imaging set-up, Fig. 1. The transfer function of the optical system between the specimen and the image sensor plane can be written as h k (x, y, z k , λ k , θ k , ϕ k ), where x and y are the lateral coordinates at the sensor plane, z k is the vertical sample-to-sensor distance, λ k is the illumination wavelength, and (θ k , ϕ k ) defines the angle of incidence. The subscript k denotes different imaging configurations, achieved by e.g., vertically moving the specimen or sensor chip to record the holograms at different sample-to-sensor distances z k , changing the illumination wavelength λ k , or tilting the illumination beam to change the angle of incidence, θ k and ϕ k . One additional pair of variables in our imaging configuration is the lateral displacements between the image sensor and the object planes, i.e., x shift,k and y shift,k , see Fig. 1a. Such sub-pixel displacements are utilized as one way of mitigating the spatial undersampling at the image sensor chip due to a large pixel size.
Under these different imaging configurations, each labeled with index k, the transmission properties of a two-dimensional (2D) specimen can be generally expressed as o k (x, y) = 1 + s k (x, y), where s k refers to the scattered object field that interferes with the background unscattered light. The frequency spectrum O k (f x , f y ) of o k (x, y) can be written as: Similarly, we can write the 2D spatial frequency spectrum of the transfer function h k (x, y, z k , λ k , θ k ) as: where FT refers to the Fourier Transform operation. From now on, we will simplify the expressions of all the frequency spectra in our equations by hiding the spatial frequency variables f x , and f y . The frequency spectrum of the field intensity i k (x, y) on the image sensor plane can then be expressed as: where '·' represents the multiplication operation, the superscript '− ' represents using variable set (− f x , − f y ) instead of (f x , f y ) and the asterisk stands for complex conjugate operation. SS k represents the self-interference terms, which can be written as SS k = Γ fx,fy {H k · S k }, where Γ fx,fy refers to the autocorrelation operation. T k is determined by the transfer function H k , i.e.,: where f x,k = n k · sin θ k ·cos ϕ k /λ k , f y,k = n k · sin θ k · sin ϕ k /λ k , and n k is the refractive index of the medium, which is assumed to be a function of only the illumination wavelength. It is important to notice that H k is a complex function with a unit magnitude, defining a phasor 52 . Based on Eq. (4), as a product of ⁎ H f f ( , ) k x k y k , , and H k , the function T k is also a phasor, and we term T k as a propagation phasor, the function of which in our reconstruction framework will be more clear later on. When any intensity distribution i k (x, y) is sampled by an image sensor-array with a pixel pitch of Δx and Δy in lateral directions, the discrete Fourier transform (DFT) of the sensor's output can be expressed as: , 0, 1, 2, In Eq. (5) u and v are integers representing the aliasing orders, and (u, v) = (0, 0) denotes the non-aliased target signal of the object. P k (f x , f y ) is the 2D FT of the pixel function that defines the responsivity distribution within each pixel of the image sensor chip 30 . Originally, f x , and f y in Eq. (5) are discrete frequency values confined within the Nyquist window. Based on the periodic nature of DFT, Eq. (5) and all of our further derivations can be numerically extended to a broader frequency domain by simply upsampling the raw measurements. Therefore, without change of notations, I sampled,k refers to the DFT of the upsampled version of our raw measurements. Now we will incorporate the lateral displacements between the holograms and the image sensor chip into Eq. (5). If we add lateral shifts (x shift,k , y shift,k ) to each hologram, then Eq. (5) can be re-written as: where we simplify the expression of spatial aliasing order by using the subscript uv, and φ shift,uv,k represents the phase change caused by a lateral shift: On the right side of Eq. (8), we can see that, for each aliasing order (i.e., each combination of u and v, including the target signal: u = 0, v = 0), there are four items inside the square brackets. The first item, δ uv , represents the background light, the second item, T uv,k · S uv,k , represents the real image, the third item, , , represents the twin image; and the last item, SS uv,k , is the self-interference term.
In the next sub-sections, we will present a generic, two-stage holographic reconstruction algorithm using propagation phasors, which aims to recover the object term δ 00 + S 00,k from a series of measured holograms.

Stage I of Propagation Phasor based Holographic Reconstruction: Generation of an Initial
Guess. As depicted in Fig. 2, the first stage of the reconstruction is to generate a high-resolution initial guess of the specimen, and this Stage I is composed of three steps (i.e., Steps 1-3 in Fig. 2).
Step 1: Upsampling of each raw measurement serves as the first step in our holographic reconstruction algorithm. This upsampling factor, although does not introduce any new information, should be large enough to expand the expression of I sampled,k to cover the entire passband of the optical system. Since the computation cost of the reconstruction increases quadratically with the upsampling factor, it should also be limited to avoid unnecessary computational burden/time. For our lensfree microscopy platform reported here, we typically set an upsampling factor of ≤ 7.
Step 2: The second step of the holographic reconstruction is to offset the lateral displacements x shift,k , and y shift,k , and then perform back-propagation on the upsampled raw measurements. To do so, we multiply both sides of Eq.
,00, and reorganize the terms to extract the true object signal, i.e., the target signal: On the left side of Eq. (9), we have kept the pixel function P 00,k multiplied with δ 00 + S 00,k ; note, however, that it can be later removed using deconvolution techniques as the last step of the holographic reconstruction 30 . The right side of Eq. (9) shows that in order to extract(δ 00 + S 00,k · P 00,k , there are five terms that need to be eliminated from the back-propagated intensity (i.e., ). The first term, 00, 00, 00, , represents the twin image noise; the second and third terms which contain S uv,k or −⁎ S uv k , (u ≠ 0, v ≠ 0) represent the spatial aliasing signals for real and twin images, respectively; the fourth term with δ uv (u ≠ 0, v ≠ 0) is the high frequency artifacts generated during the upsampling process. The last term with SS uv,k is the self-interference signal.
Scientific RepoRts | 6:22738 | DOI: 10.1038/srep22738 Step 3: Summation of all the upsampled and back-propagated holograms to generate an initial guess. This initial summation can greatly suppress the twin image noise, aliasing signal and other artifact In each subfigure, we scan one of the imaging parameters while keeping all the others constant, and plot the values of the phasors in color-coded points on the unit circle. The first row shows the twin image phasor values as a function of (b) the lateral shifts between the hologram and the image sensor-array, x shift,k and y shift,k , (c) the sample-to-sensor-distance, z k , (d) the illumination wavelength, λ k , and (e) the illumination angle θ k . Similarly, the second row shows the spatial aliasing phasor values as a function of (f) x shift,k and y shift,k , (g) z k , (h) λ k , and (i) θ k .
Scientific RepoRts | 6:22738 | DOI: 10.1038/srep22738 terms outlined above in Step 2. To better explain the impact of this summation step, we can simplify the expression of the phasor terms in Eq. (9) as: are also phasors with unit amplitudes, and their phases change as a function of all the imaging parameters (i.e., z k , λ k , θ k , ϕ k , x shift,k , and y shift,k ), see e.g., Figs 3 and 4.
Also notice that except the illumination wavelength λ k , the changes of the imaging parameters z k , θ k , ϕ k , x shift,k , and y shift,k do not affect the transmission properties of the 2D specimen. During the imaging process, we confine the illumination wavelengths within a narrow spectral range, typically less than 10 nm, so that the transmission properties of the specimen and the image sensor's pixel function can be approximately considered identical when generating an initial guess of the object, i.e., S uv,k ≈ S uv , and P uv,k ≈ P uv . If we list Eq. (10) for all the possible K imaging conditions (e.g., as a function of various illumination wavelengths, sub-pixel shifts, etc.), and then sum them up with a set of weighting factors, {c k }, we can have:  By finding a set of weighting factors {c k } that satisfy , we can have "complete elimination" of the twin image noise, aliasing signals and upsampling related spatial artifacts, while still maintaining the target object function, (δ 00 + S 00 )·P 00 . However, considering the fact that Q uv k twin , , Q uv k real , and Q uv k artifact , are also functions of spatial frequencies (f x , f y ), it is computationally expensive to obtain a set of ideal {c k } values. Therefore we adopt an alternative strategy as shown in Fig. 2   In fact, as illustrated in Fig. 3, with proper selection of the imaging configuration, the summations of these phasors can be significantly smaller than K. This implies that, by simply summing up Eq. (11) for all K imaging configurations, the twin image noise −⁎ S ( ) 00 , aliasing signals (S uv,k and −⁎ S uv k , , u ≠ 0,v ≠ 0) and upsampling related can be significantly suppressed in comparison with the target signal (δ 00 + S 00 ) · P 00 . Therefore, we consider a simple summation as a good initial guess of the specimen at this Stage I of our propagation phasor based holographic reconstruction approach, i.e., This initial guess is then used as the input to an iterative algorithm (Stage II) to reconstruct and refine the object function/image, which will be detailed in the next sub-section.

Stage II of Propagation Phasor based Holographic Reconstruction: Iterative Image
Reconstruction. Using the initial guess defined by Eq. (13), we next implement an iterative process as the second stage of our propagation phasor based holographic reconstruction algorithm to eliminate the remaining twin image noise, aliasing signal, and the upsampling related artifacts. Each iteration of Stage II is comprised of four steps (i.e., Steps 4 through 7-see Fig. 2): Step 4: Based on the parameters of each illumination condition, (i.e., z k , λ k , θ k , ϕ k ), we apply a phase modulation on the initial guess of the specimen, defined by Eq. (13), and propagate the field from the object plane to the image sensor using the angular spectrum approach 52 . For this wave propagation, we use the free space transfer function: We term the wave propagation from the object plane to the image sensor as forward-propagation, and denote the spatial form of the forward-propagated field as g forward,k (x, y). Note that the Fresnel transform based digital wave propagation can also be used at this step, although for high-resolution imaging applications the angular spectrum approach is more suitable without any low NA approximations.
Step 5: On the image sensor plane, we use the raw measurements (i.e., the low-resolution, undersampled holograms) to update the amplitude of the high-resolution, forward-propagated field g forward,k (x, y). To do so, we first convolve the intensity of the field, | g forward,k (x, y)| 2 , with the pixel function of the image sensor 30 , and shift the convolved intensity by an amount of (x shift,k , y shift,k ) to compensate the corresponding lateral displacement. Next, this shifted intensity is downsampled to the same resolution as the raw measurement, and the difference between this downsampled intensity and the raw measurement is considered as a low-resolution correction map. In order to apply this low-resolution correction map to each shifted intensity, we upsample this correction map by taking its Kronecker product with the pixel function, and add the upsampled correction map to the shifted intensity with a relaxation factor (typically ~0.5). Then this 'corrected' intensity is deconvolved with the pixel function using Wiener deconvolution, and shifted back in place by the amount of (− x shift,k , − y shift,k ). The Wiener filter takes into account the measured noise level of the image sensor to avoid over-amplification of noise during each iteration. We then use the square root of the deconvolved and shifted intensity to replace the amplitude of g forward,k (x, y), while keeping its phase unaltered.
Step 6: Back-propagate the amplitude-updated, high-resolution field to the object plane, and remove the phase modulation caused by the illumination angle.
Step 7: The back-propagated field is then used to update the transmitted field on the object plane. Different from Step 6, this update on the object plane is carried out in the spatial frequency domain. The spatial frequency region for this update is a circular area centered at f x,k = n k · sin θ k ·cos ϕ k /λ k , f y,k = n k · sin θ k ·sin ϕ k /λ k , and we choose the radius of the circle so that all the spatial frequencies within it experience less than 3dB amplitude attenuation during wave propagation. This update in the spatial frequency domain is also smoothened using a relaxation factor of ~0.5. In other words, the updated frequency region is the weighted sum of the old transmitted field and the back-propagated field, and the weighting factor (i.e., relaxation factor) for the back-propagated field Scientific RepoRts | 6:22738 | DOI: 10.1038/srep22738 is ~0.5. After this update, we convert the phase of the field into an optical path length map of the object, and the amplitude of the field gives us the object's final transmission image, i.e., reconstruction. Note that for relatively thick specimen, phase unwrapping needs to be performed before converting the reconstructed phase into an optical path length 35 .
These above outlined steps (Steps 4 to 7) are performed for every imaging configuration. It is considered as one iteration cycle when all the K raw measurements are used for once. Similar to the convergence condition defined in Ref. 53, we determine the convergence of our iterations and the reconstruction when the sum-squared error (SSE avg itr ) between the raw measurement and the downsampled intensity map satisfies the following criterion: where 'itr' is the index of the iteration cycle, and ε is the convergence constant, empirically defined as ~0.2% of SSE avg itr .

Computation platform used for propagation phasor based holographic reconstructions. For
proof-of-concept implementation, our propagation phasor approach based reconstruction algorithm has been implemented using MATLAB (Version R2012a, MathWorks, Massachusetts, USA) on a desktop computer with 3.60-GHz central processing unit (Intel Xeon E5-1620) and 16 GB random-access memory. Using an upsampling factor of seven, the computation time of one iteration in reconstruction Stage II (detailed in the previous sub-section) is ~1.2 seconds for a region-of-interest of ~1 × 1 mm 2 . As for the total computation time including Stages I and II, assuming that the number of intensity distribution updates is ~8-10 per iteration (see e.g. Fig. 5b,d and Fig. 6b,d), and that the convergence can be reached within ~6-7 iteration cycles, the total image reconstruction time ranges between ~1-1.5 minutes per 1 mm 2 . More than 85% of this computation time is spent on wave propagation between the sample and the image sensor planes, which heavily relies on Fast Fourier Transforms (FFTs). Therefore, the adoption of graphic processing units (GPUs) or other parallel computing architectures could significantly reduce the total computation time 21 .

Results and Discussion
The main challenges of wide field-of-view, high-resolution holographic imaging include: (1) phase retrieval, and (2) mitigating the undersampling caused by an image sensor chip. The propagation phasor approach of this manuscript relies on the fact that in the digital hologram of a specimen, the twin image noise and spatial aliasing signals vary under different imaging configurations. Such variations enable us to eliminate these unwanted noise terms (twin image noise and aliasing signal) and obtain phase-retrieved and high-resolution (i.e., super-resolved) reconstructions of the object. The imaging configuration in a holographic microscope can in general be changed by varying different parameters: (1) the lateral displacements between the holograms and the sensor-array (i.e., lateral relative shifts x shift,k and y shift,k ), (2) the sample-to-sensor distance (z k ), (3) the illumination wavelength (λ k ), and (4) the angle of incidence (θ k , ϕ k ). In this section, to better illustrate the inner workings of our propagation phasor approach, we will first demonstrate the dependencies of the twin image noise and the aliasing signal on these controllable imaging parameters and then explore and summarize the combinations of these imaging parameters that can create phase-retrieved and high-resolution reconstructions while also improving the data efficiency of holographic imaging. ). From the perspective of our propagation phasor approach, we desire that the phasors that modulate these unwanted noise terms or artifacts exhibit sufficient variations across [0, 2π ], so that they can be significantly suppressed during the initial summation in the reconstruction Stage I (detailed in the Methods Section). In this manuscript, we focus our discussion on twin image phasor Q k twin 00, and aliasing related phasors Q uv k real , , Q uv k twin , , (u ≠ 0, v ≠ 0), where the conclusions would be broadly applicable to a wide range of holographic imaging systems (lens-based or lensfree). Meanwhile, the self-interference patterns/artifacts are much weaker in signal strength compared to the holographic interference terms and can be easily suppressed by the iterative reconstruction algorithm (Stage II) that is detailed in the Methods Section.

Dependency of Twin Image Noise and Aliasing
To illustrate the dependencies of the twin image noise and the aliasing signal on the holographic imaging parameters, we choose the twin image phasor  1, v = 1), as examples and visualize them as a function of the imaging parameters (x shift,k , y shift,k , z k , λ k , θ k , and ϕ k ) as shown in Fig. 3. In each sub-figure of Fig. 3, we only change one of the imaging parameters while keeping all the others constant. For instance, in Fig. 3b that shows e jφ twin , we only change the lateral shift x shift,k from 0 μm to 1.12 μm (i.e., the pixel pitch of the image sensor chip used in our experiments) with a step size of ~0.11 μm, while the other parameters are fixed at z k = 150 μm, λ k = 500 nm, θ k = 0°, and ϕ k = 0°. Similarly, Fig. 3c through Fig. 3e depict e jφ twin as a function of z k , λ k , and θ k separately, while Fig. 3g through Fig. 3i show e jφ alias as a function of x shift,k , z k , λ k , and θ k , respectively.
From Figs 3b-i we can see that, except the twin image phasor's insensitivity to lateral shifts, the diversity of all the other imaging parameters can cause both the twin image phasor and the aliasing phasors to be modulated. To better illustrate these phasors' sensitivities to various imaging parameters, we calculated in Fig. 4  derivatives of φ twin and φ alias with respect to x shift,k , y shift,k , z k , λ k , θ k and ϕ k . Next we will analyze the values of these partial derivatives along the f x axis (i.e., f y = 0), and summarize each imaging parameter's effect on φ twin and φ alias (see Fig. 4a-h).

Lateral shifts (x shift,k , y shift,k ). Since the twin image phasor
00, 00, twin (see Eq. 4) does not contain variables x shift,k or y shift,k , the absolute value of its partial derivatives with respect to x shift,k and y shift,k is zero, i.e., |∂φ twin /∂ x shift,k | = 0 and |∂φ twin /∂ y shift,k | = 0 (Fig. 4a). In other words, lateral shifts do not introduce any variations in the twin image noise term as a result of which they are not directly useful for twin image elimination or phase retrieval. On the other hand, as illustrated in Fig. 4e, when spatial aliasing exists in either x or y direction (i.e., u ≠ 0, v ≠ 0), we then have |∂φ alias /∂ x shift,k | > 0 and |∂φ alias /∂ y shift,k | > 0 , which suggests that x shift,k and y shift,k introduce linear phase modulations (see Eq. 7) in the spatial aliasing phasor term. This linear relationship between φ alias and (x shift,k , y shift,k ) makes the lateral shifts ideal choice for aliasing signal elimination. As shown in the Supplementary Materials, if we set the lateral shifts to be evenly distributed within one pixel pitch, where x shift,k ∈ {m/(M · Δx)|m = 1,2,… M} and y shift,k ∈ {n/(N · Δy)|n = 1,2,… N} summing up the upsampled and back-propagated holograms (i.e., Stage I of the reconstruction algorithm detailed in the Methods Section) can lead to complete elimination of the aliasing signals. This summation is mathematically equivalent to back-propagating the pixel super-resolved holograms 16,22,25,26,30-33,40-43,54, . To conclude, the diversity of the lateral shifts can only contribute to the aliasing signal elimination, i.e., pixel super-resolution.
Sample-to-sensor distance (z k ). Using the diversity of the sample-to-sensor distance (z k ) to eliminate the twin image noise has been one of the most widely-used phase retrieval techniques in holographic image reconstruction 13,25,27,32,36,37,49,50 . For completeness of our discussion, here we analyze the effect of z k on the twin image noise from the perspective of the propagation phasor approach. As shown in Fig. 4b, |∂φ twin /∂z k |rises as spatial frequency f x increases. Except at very low spatial frequencies (e.g., |f x | < 0.1 μm −1 ), φ twin exhibits strong sensitivity to z k . For example, even at |f x | ≈ 0.1 μm −1 , changing the sample-to-sensor distance by ~100 μm can make the twin image phasor e jφ twin reverse its polarity. This sensitivity makes z k a very useful variable for twin image noise elimination. For aliasing signal elimination, as depicted in Fig. 4f, we can see that φ alias also shows a good sensitivity to z k , i.e. |∂φ alias /∂z k | ≥ 0.01π except for a very limited number of spatial frequency points. Therefore, besides twin image elimination, the diversity of z k can also be used for aliasing signal elimination.
Wavelength (λ k ). The diversity of illumination wavelength can be used for twin image elimination (i.e., phase retrieval) 46,55 . We have previously reported that it can also be used for eliminating the spatial aliasing signals 35 . As shown in Fig. 4c,g, one important property of |∂φ twin /∂λ k | and |∂φ alias /∂λ k | is that they show strong dependencies on the illumination wavelength only when the sample-to-sensor distance z k is large enough (e.g., z k > ~100 μm). Stated differently, by changing the illumination wavelength λ k , the holographic interference patterns at the sensor-array will surely vary, but such variations become more pronounced and useful at larger distances, z k . Therefore, in a point-to-point focused imaging system (using e.g., a lens-based imaging set-up), the diversity of  wavelength is of no use for phase retrieval or resolution enhancement unless a slight defocus (i.e., z k ) is introduced in the imaging system.
Angle of incidence (θ k , ϕ k ). We have previously reported the use of the diversity of illumination angles (θ k and ϕ k ) for phase retrieval 16,22,33,38 as well as for expanding/improving the frequency bandwidth, i.e., the spatial resolution through a synthetic aperture approach in lensfree on-chip microscopy 33 . As shown in Fig. 4d,h, similar to the case of wavelength diversity, to make use of the illumination angle for phase retrieval and elimination of aliasing signal, sufficient sample-to-sensor distance (e.g., z k > 100 μm) is needed. Fig. 4d also suggests that, for phase retrieval, relatively large angular variations (e.g., Δθ > 10°) are preferred since |∂φ alias /∂θ k | > 0.1π·degree −1 . Another important observation from Fig. 4h is that at different illumination angles θ k , |∂φ alias /∂θ k | remains non-zero in most of the spatial frequencies, which is similar in behavior to |∂φ alias /∂x shift,k |as shown in Fig. 4e. Intuitively, this implies that slight perturbations on the illumination angle will introduce lateral shifts of the interference patterns on the image sensor plane, which can be considered as one method of generating x shift,k , and y shift,k . In fact, shifting the light source by small amounts has been proven as an effective way of performing lateral shift-based pixel super-resolution in lensfree holography 16,21,22 .
Regarding the parameter ϕ k , although not depicted in Fig. 4, it is important to emphasize that |∂φ twin /∂ϕ k | = 0 and |∂φ alias /∂ϕ k | = 0 when θ k = 0, and that the sensitivity of both φ twin and φ alias to ϕ k increases with θ k . Therefore, both θ k and ϕ k can be used for the elimination of twin image noise and spatial aliasing signal.
The above-described contributions of various imaging parameters to eliminate twin image noise and spatial aliasing signal terms are summarized in Table 1. From Table 1 we can see that the propagation phasor approach of this manuscript: (1) provides a unique mathematical formalism that combines/merges various existing phase retrieval and pixel super-resolution techniques used in digital holography into the same unified framework, and (2) creates two new techniques to eliminate the aliasing signal in digital holography, namely using the diversity of the sample-to-sensor distance, and the diversity of the illumination angle. For consistency with the previous used terminology, we name these two new methods as multi-height based pixel super-resolution and multi-angle based pixel super-resolution, respectively. Next, we will experimentally demonstrate the imaging results and the advantages of these two new methods.

Propagation Phasor Approach Using Multi-height and Multi-angle Holographic Data.
Using this new propagation phasor based reconstruction framework, the diversities of sample-to-sensor distance or illumination angle can enable not only twin image elimination, but also resolution enhancement, i.e., super-resolution. To demonstrate the resolution enhancement brought by the diversity of z k (i.e., multi-height based pixel super-resolution - Table 1), we captured the holograms of a standard resolution test target at eight different heights, where the values of z k are evenly distributed between 200 μm and 305 μm with a spacing of ~ 15 μm. For comparison, we first reconstructed the specimen using a previous technique: multi-height based phase retrieval algorithm 25,27,32 (see Fig. 5a). For the same set of raw data, compared to this previous technique our propagation phasor based reconstruction delivers a half-pitch resolution improvement from ~0.87 μm to 0.69 μm, corresponding to a numerical aperture (NA) improvement from 0.3 to 0.4 (wavelength: 530 nm), see Fig. 5b.
In addition to multi-height based pixel super-resolution, a similar resolution enhancement can also be achieved using the diversity of illumination angles (i.e., multi-angle based pixel super-resolution - Table 1). As shown in Fig. 5c,d, we demonstrated multi-angle pixel super-resolution using the data captured from 9 different illumination angles, where one of them is vertical (0°), and rest of the angles are placed at ± 15° and ± 30° along two axes above the specimen (see Fig. 1b). The half-pitch resolution improvement brought by the diversity of illumination angle is also similar: from ~0.87 μm down to 0.69 μm.
In the next sub-section we will demonstrate that much higher resolution images can be reconstructed using our propagation phasor approach by simply adding lateral shift based pixel super resolution to only one of the measurement heights or angles, which is used as an initial guess at Stage I of our reconstruction algorithm detailed in the Methods Section. As will be presented next, this approach is also quite efficient in terms its data requirement compared to existing approaches.

Improving the Data Efficiency in High-resolution Holographic Reconstructions Using the
Propagation Phasor Approach. Using the multi-height imaging configuration outlined earlier, we performed 4 × 4 lateral shift-based pixel super-resolution at only one sample-to-sensor distance (i.e., ~190 μm), which added 15 extra raw measurements/holograms to the original data set that is composed of measurements at 8 heights. In our propagation phasor based reconstruction, we directly used the back-propagation of this super-resolved hologram at this height (190 μm) as our initial guess (Stage I of our algorithm -see the Methods Section). The resolution improvement that we have got by using these additional 15 raw measurements in our propagation phasor approach is significant: we achieved a half-pitch resolution of ~0.55 μm (corresponding to an NA of ~0.48 at 530 nm illumination), which is the same level of resolution that is achieved by performing lateral shift-based super-resolution at every height (see Fig. 6a,b). In other words, to achieve the same resolution level, the propagation phasor approach utilized 5.5-fold less number of raw measurements (i.e., 23 vs. 128) compared to the conventional lateral shift-based multi-height method 25,27,32 . A similar level of improvement in data efficiency of our propagation phasor approach is also observed in the multi-angle imaging configuration: by simply performing 6 × 6 pixel super-resolution at only the vertical illumination, the propagation phasor based reconstruction can achieve a half-pitch resolution of ~0.49 μm (corresponding to an NA of ~0.53 at 530 nm illumination). As a comparison, the synthetic aperture approach 33 achieves a half-pitch resolution of ~0.44 μm; however it uses 6 × 6 pixel super-resolution at every illumination angle (Fig. 6c), and therefore our propagation phasor approach (Fig. 6d) has 7-fold improvement in its data efficiency (i.e., 44 vs. 324 raw measurements). This improvement and significant reduction in the number of raw measurements/holograms are especially important to make wide-field, high-resolution holographic imaging suitable for high speed applications.

Imaging Biological Samples Using the Propagation Phasor Approach.
To demonstrate the success of our propagation phasor approach in imaging biological specimen, we imaged unstained Papanicolaou (Pap) smears (see Fig. 7a-d) and stained blood smears (see Fig. 7e-h). For Pap smear imaging, we captured the holograms of the specimen at multiple sample-to-sensor distances, and at each z k , only one raw measurement is recorded. For comparison, we first reconstructed the Pap smear using a previously reported multi-height phase retrieval algorithm 25,27,32 (Fig. 7a). Using the same holographic data set and raw measurements, the reconstructions created by our propagation phasor approach (Fig. 7b) show resolution improvements compared to the previously reported method. To further improve the resolution without significantly increasing the burden of data acquisition, we added eight extra raw measurements for shift-based pixel super-resolution (with a super-resolution factor of 3 × 3) at only one of the heights, which is used as an initial guess (in Stage I) of our reconstruction algorithm. As shown in Fig. 7c, our propagation phasor approach based reconstruction shows a good agreement with the images captured using a conventional phase contrast microscope (40 × objective lens, NA = 0.6). For imaging of stained blood smears, we captured the lensfree holograms at multiple illumination angles. The comparison between Fig. 7e and Fig. 7f also confirms the resolution improvement brought by our propagation phasor based reconstruction algorithm. By adding lateral shift-based pixel super-resolution (with a super-resolution factor of 3 × 3) at only the vertical illumination angle (i.e., θ k = 0), we further improved the resolution of our reconstructed image (Fig. 7g), which shows comparable performance against a bright-field microscope with a 40 × objective lens (NA = 0.6), Fig. 7h.  33 . At each of the 13 illumination angles (− 30° : 10° : 30°, two axes), only one raw measurement is used. (f) Propagation phasor approach-based reconstruction using the same data set used in (e). (g) Propagation phasor approach-based reconstruction, where at the normal illumination, 9 raw measurements are used for pixel super-resolution, while at each one of the other angles, only one raw measurement is used. (h) Conventional bright-field microscope image of the same blood smear using a 20× objective lens (NA = 0.45).
Based on these results, we confirm that our propagation phasor approach would greatly increase the speed of high-resolution and wide-field holographic microscopy tools. In previously reported holographic imaging modalities, multiple laterally shifted images are captured to achieve pixel super-resolution at every one of the sample-to-sensor distances 25,27,32 or illumination angles 33 . As demonstrated in Figs 6 and 7, the propagation phasor approach can reduce the number of required raw holograms by five to seven fold while also achieving a competitive resolution. This reduction in raw data also lowers the need for data transmission and storage, which could further improve the cost-effectiveness of holographic imaging modalities such as handheld lensfree microscopy tools 22,27,34 for telemedicine applications.
Although our experimental demonstrations in this manuscript utilized a lensfree on-chip imaging set-up, we would like to once again emphasize that this propagation phasor approach is broadly applicable to a wide range of holographic imaging modalities, including lens-based holographic microscopy techniques. For instance, in a lens-based undersampled holographic imaging system, multi-height pixel super-resolution can simply be achieved by capturing a series of defocused images at different heights. Considering the fact that the depth focusing operation is naturally required and performed every time a sample is loaded onto a lens-based traditional microscope, this propagation phasor approach provides a unique method to enlarge the space-bandwidth-product of the final image without compromising the image acquisition time.

Conclusions
In this manuscript, we demonstrated a propagation phasor approach for high-resolution, wide-field holographic imaging with significantly improved data efficiency. Different from previous holographic reconstruction methods, our propagation phasor approach merges phase retrieval and pixel super-resolution techniques into a unified mathematical framework, where the twin image noise, spatial aliasing signals and other digital artifacts are all interpreted as noise terms that are modulated by phasors. These propagation phasors analytically depend on and can be controlled by various imaging parameters such as the lateral displacement between the hologram and the sensor-array, sample-to-sensor distance, illumination wavelength, and the angle of incidence. We systematically investigated and summarized the sensitivities of both the twin image noise and the aliasing signal to these imaging parameters, which enabled us to establish two new super-resolution methods that utilize the diversity of the sample-to-sensor distance and the diversity of the illumination angle. Compared to previous reconstruction algorithms, this propagation phasor framework can deliver phase-retrieved reconstructions with a competitive resolution using five-to seven-fold reduced number of raw measurements/holograms, which makes it especially appealing for high speed and cost effective microscopy applications. We further confirmed the success of this approach by imaging biological samples including unstained Papanicolaou smears and stained blood smears.