Regularized non-convex image reconstruction in digital holographic microscopy

: Inverse problem approaches for image reconstruction can improve resolution recovery over spatial ﬁltering methods while reducing interference artifacts in digital o ﬀ -axis holography. Prior works implemented explicit regularization operators in the image space and were only able to match intensity measurements approximatively. As a consequence, convergence to a strictly compatible solution was not possible. In this paper, we replace the non-convex image reconstruction problem for a sequence of surrogate convex problems. An iterative numerical solver is designed using a simple projection operator in the data domain and a Nesterov acceleration of the simultaneous Kaczmarz method. For regularization, the complex-valued object waveﬁeld image is represented in the multiresolution CDF 9 / 7 wavelet domain and an energy-weighted preconditioning promotes minimum-norm solutions. Experiments demonstrate improved resolution recovery and reduced spurious artifacts in reconstructed images. Furthermore, the method is resilient to additive Gaussian noise and subsampling of intensity measurements.


Introduction
Digital holography is a technique to capture simultaneously the amplitude and the phase of an incoming wavefront using a digital image sensor.Various configurations for holographic setups exist [1], such as Gabor holography, Fourier holography, phase-shifting holography and off-axis holography.In off-axis holography, the object wavefront is retrieved by using a high-frequency off-axis carrier wave to separate the desired first order from the unwanted zero order and twin image terms in the Fourier domain.
The off-axis geometry allows users to capture the wavefront with a single recording, unlike e.g. in phase-shifting holography, in which stability is crucial for a realization of at least three exposures.Single exposure holography is mandatory when the recorded object is dynamic.However, the drawback is that diffraction orders must be well-separated in the Fourier domain for subsequent filtering.Stringent constraints have to be put on the placement and bandwidth of the first order term.This will cause much of the digital bandwidth to be underused and thereby, discarding information about the object.Advanced reconstruction methods are therefore essential for resolution recovery.
Classical approaches for extracting a complex wavefield from intensity measurements relies on direct manipulations in the Fourier domain [2].Assuming that the object under investigation is band-limited, the Fourier filtering technique may become optimal using clever non-linear filters that are implemented using a log-transform on the intensity measurements [3,4].Unfortunately, strict bandlimitation is not a realistic hypothesis in lens-free setups.
Inverse problem approaches have been successfully applied for resolution recovery in a wide range of modalities.The seminal works in digital holography by Sotthivirat and Fessler [5][6][7] propose to reconstruct the complex wavefront using a penalized likelihood objective function with an edge-preserving filter for regularizing the solution image.They solved a surrogate problem using the optimization transfer principle and demonstrated, on horizontally-modulated artificial data, the feasibility of reconstructing complex images from arbitrary reference beams.
Subsequently, the work of Bourquard et al. [8] devised a practical method for the reconstruction of acquired data.They represented the solution complex image with a dense array of pixels and used a dual-loop method where data similarity is ensured up to a given precision, then the solution is regularized in the image space, using a total variation minimization prior [9].
In principle, this method is only able to match the measurements approximately and the interleaving of two separate objective functions in the data and solution spaces leads to slow convergence.Indeed, an exhaustive line search was used for solving the two optimization problems.Nonetheless, they demonstrated encouraging results as the method reduced greatly border artifacts and slightly increased resolution recovery compared to classical direct Fourier filtering.Furthermore, their inverse approach was robust to data undersampling and measurement noise could be mitigated by tuning the tolerance bound in the data fidelity objective.
Reconstruction from undersampled data has been initially investigated by Riverson et al. [10], in the case where complex measurements are recorded, for example through multiple acquisitions with a controlled global phase shift.In this case, the reconstruction problem is convex but underdetermined if fewer measurements than image pixels are recorded.This draws a link with the compressive sensing framework [11,12] that relies on sparse regularization.Marim et al. [13] have used compressive sensing in the hologram space for explicitly recovering missing values.
In this work, we took a different stance on the inverse problem approach for reconstructing a regularized image from off-axis intensity measurements.We propose a numerical method that is able to converge to a strictly compatible solution to the non-convex image reconstruction problem using a simple consistent projection operator in data domain.We use as solver a fast Nesterov acceleration of the simultaneous Kaczmarz iterative method.No explicit regularization is implemented in image space, but instead we select the minimum-norm representation of the solution in a multiresolution wavelet domain.Experiments conducted on a synthetic resolution chart demonstrate drastic resolution recovery improvements in comparison to the best compromise we can achieve with direct Fourier reconstruction.Moreover, robustness to noise on the intensity measurements can be obtained with a modified projection operator.
We observed that if the detector would be able to measure the complete complex values, i.e. both the amplitude and phase information of the incident light, then the inverse reconstruction problem would be convex and could be solved using well established numerical methods.We create a sequence of intermediate surrogate convex problems using inferred complete complex data.This inference step uses simple data-consistency conditions and shifts the complexity away from the numerical solver.Our method is also robust to data undersampling and implicitly interpolates the missing intensity data using both wavelet decompositions and an accurate physically-based forward propagation model using the angular spectrum operator.
For regularization, the real and imaginary parts of the complex sample image to reconstruct are represented in the Cohen-Daubechies-Feauveau (CDF) 9/7 wavelet space [14].Our prior studies [15] demonstrated that combining CDF wavelets with the angular spectrum light transport operator was a more appropriate choice over Fresnelets [16].This indirect minimum norm objective in coefficient space yields global multiresolution image regularization that is very different than direct processing in image space.
The remainder of the paper is structured as follows.We first present the geometry of off-axis holography and the image formation model in Section 2. Our regularized inverse reconstruction method is derived in Section 3. Experiments on noisy and undersampled synthetic data are shown in Section 4, as well as results from resolution charts and biological samples.Finally, we conclude with an outreach on future works in Section 5.

Off-axis holography
Following the illustration in Fig. 1, an off-axis intensity hologram is measured on the digital detector as an interference pattern I formed by the superimposition of an object beam O containing the object information and a tilted incident carrier planar wave, called the reference beam R. For more details, the fundamentals of holography have been studied extensively [17,18].Let I ∈ R N + be a rectangular recorded intensity image defined with N = N x × N y pixels that are vectorized in lexicographical order.In off-axis holography, we express intensities as energy measurements of the interference between a reference beam R ∈ C N and the sought object wavefield O ∈ C N illuminating the sample under investigation.Although the waveform is complex-valued, we acquire only an intensity image (1) where * denotes the complex conjugate and • is the Hadamard product.

Image formation
The distance from the detector to the sample is relatively short and therefore in this work, we rely on the angular spectrum method [18] for the light transport operator from the object wavefield O d ∈ C N at the sample plane towards the detector.The detector plane is located at the origin and the sample is displaced by a distance d along the z-axis.The angular spectrum transforms point sources to apodized chirp functions on the detector plane.The following forward and adjoint backward relations hold for small propagation distances d: The multiplication before inverse Fourier transform F −1 is equivalent to a convolution in the spatial domain.Specifically, the Hadamard product of K d and the spectrum of O d , followed by the inverse Fourier transformation, is equivalent to the convolution of O d with the point spread function F −1 (K d ), with the optical transfer function  where the distance field image D ∈ R N + defined by with the wavelength λ of the coherent illumination source, the distance d from the sample plane to the detector, the pixel pitch p and the integer coordinates x ∈ [1 . . .N x ] and y ∈ [1 . . .N y ].
In off-axis geometry, the reference beam R is a slightly slanted planar wave such that the time delay of arrival on the detector draws a modulating carrier signal that shifts the spectrum of O to obtain O • R * .Prior to acquiring I, the reference wavefield R is determined from recording a calibration intensity hologram I R ∈ R N + of the interference of R and the blank object wavefield O R ∈ C N that is is obtained without sample in the field of view.Therefore, both the object beam amplitude attenuations and phase shifts at the sample plane are known to be zero.
The first image in Fig. 2 shows a typical intensity hologram and its characteristic Fourier power spectrum.The vertical and horizontal off-axis angles of incidence have been calibrated such that the modulated object wavefield fills one full quadrant of the Fourier spectrum.This configuration leads to an optimal use of the available image bandwidth.

Direct Fourier reconstruction
From the recorded intensity hologram I, an approximation of the focused complex object wavefield O d can be reconstructed through simple manipulations in the Fourier domain.Intensities from I are decomposed into the -1, 0 and +1 diffraction orders: the first two terms in (2) correspond to the zero order, the third term is the so-called twin image and the last term contains the sought object beam O modulated by the complex conjugate of the reference beam R.
Since R is known, and assuming noise-free measurements, the possible solutions for R + O in every pixel lie on a circle in the complex plane.The right subfigure shows the corona of acceptable solutions, within an interval of confidence.The width of that interval depends on the assumed noise variance.As Fourier reconstructions do not allow to model noise, we consider in this section, solely the noise-free case.
First, we remove the non-modulated intensity information from the known reference beam: then, we demodulate the intensity signal by dividing with the complex conjugate of R: and finally, we crop the Fourier spectrum with a smoothed circular cropping disc for explicitly removing the two first demodulated terms and thereby isolating O.This filtering step is crucial for removing interference artifacts from the remaining zero order term |O| 2 /R * and the unwanted twin image depending on O * .We define a circularly-symmetric function with its center at the origin of the frequency domain as follow.First, we generate a smoothed circular mask , starting from a disc indicator function with linear ramp on the edge: with the quadrant cropping radius r = N x /4 − 1, the width w of the linear ramp at the window border and the pixel coordinates x and y defined above.Finally, we remove discontinuities in the first derivative by applying the Hermite polynomial "smoothstep" function: After masking the frequency response for selecting the object wavefield O in the Fourier domain, the result is refocused by reversing the effect of light transport from the detector plane to the sample position.All in one, the final propagated object wavefield is where K −d is a diagonal matrix which applies a phase delay based on the propagation distance d.

Regularized inverse reconstruction
From the recorded intensity hologram I, an improved approximation of the propagated complex object wavefield O d can be reconstructed through iterative refinement of a regularized complex solution in the CDF 9/7 wavelet domain.One important property of the sought propagated object beam in sample plane O d is that not all valid wavefronts O on the detector plane satisfying Eq. 2 are statistically equally likely: neighboring pixels will generally be highly correlated in the focused image O d .These correlations can be reduced significantly by using a sparsifying multiresolution transform [19].By doing so, the signal will be well-approximated with only a few coefficients that are selected according to their energy in the decomposition.
In the forward model shown in Fig. 3, we chose a discrete Wavelet transform (DWT) for the image transform.The DWT is a multiscale transform that allows for the simultaneous spatial and frequency analysis of signals [20].It separates a discrete signal into a low-pass and a high-pass signal by means of a dyadic transform, followed by a downsampling.This operation is applied recursively on the ensuing low-pass signal, resulting in a multiscale representation.For this application, we chose the biorthogonal CDF 9/7 wavelet transform.This wavelet was chosen for its symmetry, superior energy compaction and compression performance [21].Since the CDF wavelets are not orthogonal, the transpose operator will not be equal to the inverse one.For more details on the efficient implementation of the transpose, we refer the reader to Appendix 5.
Reconstructing the two channels of the complex sample wavefield O d from intensity measurements I only is ill-posed in two ways.First, the problem is physically ill-posed since many different sets of complex values for O d may be compatible with the intensity measurements of I, regardless of the optics and detector definition, i.e., the number of intensity samples.Secondly, the problem may be algebraically ill-posed if we record M < N samples.In that case, the number of measurements is smaller than the number of unknowns leading to an underdetemined inverse problem with more solution elements than data constrains.
For the aforementioned reasons, regularization is needed to select a unique solution among all compatible realizations.More importantly, the solution should be plausible as well.We propose in this work to use a minimum-norm regularization prior in the CDF 9/7 wavelet domain for selecting a unique solution among all combinations of amplitude and phase pairs that are compatible with the M recorded samples from the full intensity hologram I.
More formally, let x ∈ C N x × N y be the unique minimum-norm solution representing the object wavefield O d = W −1 x where the matrix W ∈ C N x × N y is a complete dictionary of complex wavelet basis functions such that the vector-matrix product can be implemented with a fast transform.One could aim first at solving the following non-convex problem: in which the data compatibility constraint only applies to squared magnitudes from complex values and therefore the phase information of the solution is implicitly lost.There are no globally convergent numerical algorithms for solving this problem and some sort of heuristic search method should be used instead.Since the problem is large, we will pursue along another path.Let's assume for a moment that the detector could measure both the real and imaginary parts of the complex object wavefield instead of the squared magnitudes only.From these complex-valued measurements of O, we could solve the simpler convex problem argmin in which the data compatibility constraint applies to both the amplitude and phase information and a convergent convex numerical method may be used to find this unique regularized solution.

Projection operators
In this work, we define two projection operators that infer new complex measurements O using a pilot estimate of wavelet coefficients x.At each step of a convex iterative solver, a projection will be applied such that data consistency conditions are respected.This procedure is similar in spirit to the Gerchberg-Saxton method for phase retrieval [22].This method alternately modifies the amplitude of the solution in the spatial and Fourier domain, while we update both channels simultaneously in a non-alternating scheme.
Starting from a zero vector of coefficients x 0 = 0, we first infer, at each iteration t > 0, the complex-valued measurements from the current object wavefield estimate with then, we use an orthogonal projection for recovering the phase and amplitude information from data consistency constraints.
In off-axis holography, we have intensity measurements that set the squared magnitude of complete complex values.When the measurements are assumed to be noise-free and unbiased, we can plug-in an approximation of the phase information into the current complex data estimate O using the following geometric non-linear data projection operator This exact projection operator is illustrated in the left panel in Fig. 4. The current solution will subsequently be updated using the minimum-norm difference vector such that data consistency is ensured and therefore, we have When intensity measurements are corrupted by Gaussian noise or some mismatch with the forward model remains, then a relaxed data consistency condition is more practical.Assuming a tolerance margin with respect to the expected noise variance 2 , we use ) A relaxed version of the projection operator is updating the current solution using the minimumnorm difference vector such that only approximate data consistency is ensured and thus, complex values are restricted to lie inside the corona of radius 2 around √ I: This relaxed approximate projection operator is illustrated in the right panel in Fig. 4.
For initialization, we assume a blank object wavefield O = ∅ and therefore we only use the reference beam inside the projection.This gives us the initial complex measurements The resulting initial wavefield does obviously not correspond to the ground truth image.A particularity is that fringes will appear both in the amplitude and the phase images.Therefore, the wavelet decomposition of this complex image will require lots of high-frequency coefficients to encode the fringes.The energy evaluated as the sum of squared wavelet coefficients will be considerably high.An energy-minimization prior will tend to suppress fringes from the solution, as this information is already present in the reference beam R.

Iterative solver
From the projected data-compatible complex measurements y = P(O), we seek for the minimumnorm regularized solution x by using the Moore-Penrose pseudo-inverse system matrix with the linear system matrix A = Ψ d W −1 .In practice, the square Gram matrix ( A * A) is not explicitly computed when inverted, but the minimum-norm solution x is found by an iterative method, given a pair of functions computing the forward operator y = Ax and its Hermitian adjoint x = A * y.The inverse wavelet transform and the conjugate transpose of the inverse wavelet transform are computed with fast lifting schemes that are detailed in Appendix 5.

Simultanous Kaczmarz method
At first, the simultaneous algebraic reconstruction technique (SART) [23] is adapted to our problem.SART is a parallel acceleration of the classical sequential row-action Kaczmarz updates [24] that can be seen as an averaged stochastic gradient algorithm [25,26].The convergence is monotonous as the norm of the solution increases from iteration to iteration while the mean squared data residuals decrease accordingly.
Starting from an empty initial image at t = 0, the norm of the solution is incrementally increased at each iteration with a simple incremental update rule: and the following simultaneous incremental update where the sum runs over any subset of measurements 1 ≤ M ≤ N.This update is a simultaneous row-action, as only the i th row [A] i is needed for computing the corresponding term in the sum.The structure of the update contributions inside the sum is simple and interesting: The difference term on the numerator is the discrepancy between the complex measurement [P(O)] i that is inferred from the current solution and the footprint of the current solution [A] i x, after propagation towards the detector plane.The differences are normalized by the energies [A] i 2 in data space before backprojection in the wavelet space using the adjoint operator A * .
The energy normalization values do not depend on the solution nor on data and are precomputed before reconstruction.Note that detector sampling and rows of the system matrix are not explicitly stored; rather, results of multiplications with all matrix row vectors [A] i and [A * ] i are implemented using fast transforms.The Kaczmarz method has been recently used for the closely related phase retrieval problem [27].However, convergence is extremely slow and accelerations are needed in practice, without resorting to hazardous numerical over-relaxation.

Nesterov accelerations
The 1983 vintage of Nesterov accelerations [28] uses the difference with wavelet coefficients at the previous iteration for increasing convergence speed.This variant is proven to be convergent, even though the norm of the solution may decrease after subsequent iterations.This first accelerated gradient update proposed by Nesterov introduces an initial relaxation value λ 0 = 1 and a temporary solution vector initialized as y 0 = x 0 .This quantity is incrementally updated as in the Kaczmarz method using then, the solution is updated using a supplemental additive term as This second step is akin to updating the solution with a finite difference approximation of the gradient.This update is therefore related to the Newton method and inherits from the quadratic O(1/t 2 ) convergence speed instead of the much slower linear O(1/t) convergence of the Kaczmarcz method.A remarkable property of the Nesterov acceleration, is that no extra calculations are needed beside a simple difference with the previous solution.At almost the same computational cost, we get a stable method with a guaranteed quadratic convergence regime.
The 2005 version of Nesterov's methods [29] further improves convergence speed by using an accumulated history of all previous incremental updates since the beginning of the iterative process.In addition to the above procedure, we keep track of the accumulated weighted updates: with u 0 = x 0 .Then, the final update is replaced by the linear combination and we iterate until the convergence is satisfactory.Note that the recent work of Kim and Fessler [30] has shown that the increment in ( 26) may be doubled while preserving convergence.We therefore used that variant in our numerical experiments.Figure 5 compares the evolution of the norm of the solution in function of the iteration number.Nesterov accelerations are essential for reaching reasonable computing performances.

Experiments
We implemented the described reconstruction methods in MATLAB.The source code is available upon request to the authors.We conducted experiments on both simulated and optically acquired data of a USAF-1951 resolution chart and two biological samples.The simulated resolution chart used the very same geometry for matching closely the experimental conditions of the acquisition setup.For the biological samples, we started from a real amplitude and phase image of a blood sample acquired on an iLine F holographic microscope from Ovizio Imaging Systems.We then generated artificial intensity holograms using the complex image as the ground truth.We used this realistic phantom for evaluating the impact of noise and subsampling on the image quality.We also reconstructed an acquisition of pollen seeds from the hologram used in the work of Seelamantula et al. [4].The optics used for that acquisition used bandlimiting filters, and in this case theoretically exact reconstruction is possible using their specific direct Fourier filtering technique.In practice, slight crosstalk between the zero order terms and the modulated complex wavefield is still present in the data whereas regularization alleviates artifacts appearing in the direct Fourier filtering image reconstructions.
Figure 7 shows the two simulated and acquired intensity holograms and their Fourier power spectra.Additional rows in this array show the amplitude images of the reference beam used for inverse reconstruction.the best result we could obtain using the direct Fourier reconstruction approach is shown as well.For simulations, we implemented a slight Gaussian fall-off intensity profile for modeling the effect of laser beam expansions to an approximate planar wave.

Data simulations and acquisitions
In our simulations and laboratory experiments on the USAF-1951 test target we built a digital holographic setup in transmission mode as shown in Fig. 1.The sample is illuminated using a JDSU Uniphase 1135P He-Ne laser of 20 mW operating at a wavelength of λ = 632.8nm.The planar sample was placed at a distance of 125 mm from the CCD camera: The Ximea MD120MU-SY, which is a 12-bit monochrome camera with a resolution of 4242 × 2830 pixels and a pixel pitch of δx = 3.1 µm.The exposure time was optimized for using the full dynamic range of the camera without signal clipping from saturation.
We recorded both the amplitude of the reference beam R and the off-axis hologram formed by the interference pattern.For the phase of the reference beam, a square crop was taken from the camera, followed by a four-fold downsampling to 256 × 256 pixels.We chose to modulate R as to have the carrier frequency of the first order to lie in the center of a quadrant of the Fourier domain.This configuration optimally uses the detector bandwidth in the off-axis geometry.The exact same geometry and procedure were used for the simulated hematology sample.

Results
We first evaluated the effective monotonous convergence of the iterative inverse method.As the norm of the solution increases, the complexity of the solution increases as well.From iteration to iteration, we observe a progressive sharpening until convergence.The convergence at early iterations is faster, suggesting that relaxation strategies could possibly be introduced for later Reference reconstruction Fig. 9. Progressive quality improvement with increasing amount of intensity measurements used for constraining the solution.In comparison to a full-data reference reconstruction, a high quality phase image is recovered from one quarter of all detector pixels only.
iterations. Figure 6 shows the progressive improvement of the solution on the phase image channel.The amplitude image is jointly sharpened in a similar way.
In terms of running time, the process converges to a high-quality solution after about 200 iterations.The required number of iterations depends on the complexity of the image content.For each iteration, the forward model and its adjoint are evaluated.Each of these operations is roughly equivalent to the computing cost of a direct Fourier extraction plus a wavelet transform, hence the required computing budget is about 600 fold in comparison to simple refocusing in Fourier space.A key advantage of the inverse problem approach is its resolution gain with flexible regularization priors, its ability to model noise and missing measurements.The phase shifting pattern also may be arbitrary instead of being restricted to a planar modulation wave.When implemented on GPU, the whole iterative image reconstruction task requires less than a second of computation.
For assessing the robustness to measurement noise, we corrupted measurements with 10% Gaussian noise.A 10% standard deviation from the maximum peak intensity is much higher than what is accustomed in current holographic exposures.The possibility to model noise allows however for new ultra-fast acquisitions protocols, using less sensitive detectors and/or low-light.
Figure 8 shows the impact of the tolerance radius in the relaxed data projection operator in the left part of Fig. 4. We obtained a pleasant visual result for a confidence interval of 95%, corresponding to a radius of 1.96 times the standard deviation of the assumed Gaussian noise realization.Quantitative evaluations in terms of PSNR confirmed the 95% rule for trading-off optimally between blur and variance in the final image.
Our inverse method may use subsets of intensity samples in detector space for updating the complex solution wavefront.We selected quasi-randomly detector pixels using a mask retaining 10%, 25% and 50% of all available data.The low-discrepancy image sampling technique was designed for avoiding regularities that could interfere with the regular fringe pattern and introduce aliasing artifacts in the reconstruction [31].
Figure 9 shows the progressive image improvement as more data are collected.As expected, the data loss translates into loss of resolution that manifests as blur artifacts instead of high frequency artifacts that could be dangerously interpreted as small image features.The norm of the solution is systematically lower as less data are processed.
For illustrating experiments on the simulated and acquired USAF test targets, we extracted crops in the image region showing the boundary of resolution capabilities.Figure 10 shows  the side-by-side comparison between the best results we could obtain using direct Fourier reconstruction and the regularized inverse reconstruction method.On simulation, the effective resolution gain was more prominent than on the optical acquisition.However, we observed a slight discrepancy between the acquisition of the reference beam intensity image and the hologram intensities.This issue has been encountered as well in the work of Bourquard et al. [8] and we used a first order gain correction factor as described in Appendix 5.The discrepancy may be explained by detector non-linearities and possible optical aberrations.We expect further gain with a more precise modeling of the detector response.
Finally, we reconstructed the pollen seeds sample images that were acquired in a setup with a bandlimited numerical aperture and magnification lenses.Figure 11 shows the whole image result with the iterative approach and makes the comparison with the direct Fourier method in regions of interest.The effect of non-exact lens rectification can be observed by the wave patterns at the border of the image domain.We do not observe any discrepancy between the acquisition of the blank scan for the reference beam intensity and the actual hologram, therefore no adaptive gain correction was used for this reconstruction.
The resolution improvement is more noticeable in the amplitude image, suggesting that the surface of pollens is relatively smooth.Note that the optical bandlimitation may also be responsible for limited resolution gain in the phase image.Moreover, not any parameter tweaking was needed for using the reconstruction method on the two different acquisition setups.The sole regularization through norm minimization is generic enough for selecting a plausible single complex solution image among all possible object beams that are compatible with measurements.

Conclusion
With transmissive off-axis holographic acquisition setups, it is possible to capture both the amplitude and phase information of light that is encoded implicitly in fringe interferences of a known modulating reference beam and the sought object wavefield.This work presented a convergent iterative method for image reconstruction of the complex wavefront from intensity holograms only.The reconstruction problem is ill-posed and a minimum-norm regularization prior selects a low-energy unique solution.Unfortunately, this objective function is non-convex and we solve instead a sequence of surrogate convex linear problems using a data projection operator at each iteration.The projection steps infer complex measurements that are consistent with recorded intensities.The method progressively sharpens the reconstructed image and does not require any tunable parameter except for the number of iterations.Arbitrary reference beam may be used beside off-axis modulation.Furthermore complete measurement data are not required and sparse sampling may be applied for reducing input data.Experiments demonstrate improvements in terms of resolution recovery, compared to direct Fourier reconstruction.

Appendix A. Lifting schemes for the CDF 9/7 wavelet transforms and their adjoints
The DWT may be implemented as a pair of convolution filters, followed by downsampling operations.However, it is more practical to use a lifting scheme representation [32].Lifting schemes are generic methods for constructing "second-generation" wavelets, and can be seen as a special case of a filter bank.They have more flexibility, as they allow to define wavelet bases on intervals, irregularly sampled grids and even non-linear transforms.Additionally, lifting steps have reduced computational complexity and memory constraints w.r.t. the convolutional approach.The generic lifting scheme is summarized on Fig. 12.
Lifting methods can be explained by polyphase matrix representations.The z-transform of a signal is defined as Fig. 12. Lifting scheme block diagram where an input signal x is split into even and odd samples (x e and x o ).Then a series of convolution-accumulate operations is applied alternately on the two divided signals, using prediction (s i ) and update (t i ) filters.Finally, the channels are scaled with constants K i , resulting in the final approximation signal λ and detail signal γ.
thus, the z-transform of a finite impulse response filter h = {h k b , ..., h k e } where h k b , h k e 0 is Filtering a signal x by a filter h is simply described in the z-domain as an ordinary multiplication: y(z) = h(z)x(z).Using this notation we can decompose a signal x(z) into even (x e (z)) and odd (x o (z)) components as x(z) = x e (z 2 ) + z −1 x o (z 2 ), where x e (z 2 ) = x(z) + x(−z) After splitting the signal into odd and even parts, the signal can be transformed by a filter pair (h, g) represented by the polyphase matrix P(z): λ(z) γ(z) = P(z) x e (z) x o (z) .
Any polyphase matrix P(z) representing a wavelet transform with finite impulse response filters can be factored into a product of a diagonal matrix and multiple unit upper and lower triangular matrices [33], matching the primal and dual lifting steps: Furthermore, we can directly compute the inverse transform as For efficiently solving the gradient descent problem using Nesterov's method, we not only need the inverse transform but its transpose as well.The inverse transpose transform will only be equal to the forward transform when the DWT is orthogonal.Using the formalism of the lifting scheme we can construct the transpose transform of any non-orthogonal wavelet, or more generally any linear lifting-based filter bank, as follows: In particular, we have implemented our system using the CDF 9/7 wavelet transform.The specific prediction, update and scaling constants are given in Table 1.

Fig. 1 .
Fig.1.Geometry of an off-axis holographic acquisition setup.The intensity of the interference between a slightly tilted reference wavefield and the object wavefield is recorded by the digital CCD sensor.The phase information of the observed sample is then encoded in slight local displacements of fringe patterns appearing on the detector plane.

Fig. 2 .
Fig. 2. Direct image reconstruction in the Fourier domain.The complex-valued object wavefield image O is recovered from the intensity hologram I = |O + R| 2 by demodulation and band-pass filtering of higher frequencies for removing the fringe modulation pattern of the reference beam R. Finally, the wavefield is refocused on the sample plane.

Fig. 3 .
Fig. 3. Forward image formation model used in inverse image reconstruction.CDF 9/7 wavelet coefficients represent the amplitude and phase components of the sample wavefield O d .Angular spectrum light transport and the interference between the propagated wavefield O and the near-uniform off-axis reference beam R generates a fringe pattern in I.

Fig. 4 .
Fig.4.Exact (left) and relaxed (right) projection operators infering complex-valued measurements from the current solution O.The closest value from O is selected such that addition to the reference beam R intersects the circular region of all compatible intensity values.

Fig. 5 .
Fig. 5. Comparison of the overall increasing norm as a function of the iteration number for the reference Kaczmarz method and the two Nesterov accelerated methods discussed in this work.Nesterov methods reach quadratic convergence speed instead of the much slower linear convergence of the Kaczmarz method plotted in light gray.

Fig. 6 .
Fig.6.Progressive sharpening of the solution until convergence.The regularization prior in the inverse reconstruction method selects the minimum-norm solution that is compatible with the intensity measurements.Stopping early the iterative process yields a blurry intermediate solution that may already be of sufficient quality for the operator.

Fig. 7 .
Fig. 7.The four simulated and optically acquired intensity hologram data used in experiments (top row) and their frequency power spectrum (second row).The blank scan for recording the reference beam amplitude are shown as well (third row) along with the amplitude (fourth row) and phase shift (bottom row) reference results with direct Fourier reconstruction.

Fig. 8 .
Fig. 8. Impact of the relaxed data projection operator for avoiding the introduction of noise artifacts into the solution.A bias-variance trade-off is driven by the radius of the intensity tolerance region, expressed as a fraction of the standard deviation of the noise realization.

Fig. 10 .Fig. 11 .
Fig. 10.Side-by-side comparison of resolution recovery in a similar region of interest in both the simulated and acquired USAF-1951 test targets.Parallel bars are crisper with the inverse reconstruction technique, while spurious low-frequency background noise is reduced.