Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group

A practical inverse-problem approach to digital holographic reconstruction

Open Access Open Access

Abstract

In this paper, we propose a new technique for high-quality reconstruction from single digital holographic acquisitions. The unknown complex object field is found as the solution of a nonlinear inverse problem that consists in the minimization of an energy functional. The latter includes total-variation (TV) regularization terms that constrain the spatial amplitude and phase distributions of the reconstructed data. The algorithm that we derive tolerates downsampling, which allows to acquire substantially fewer measurements for reconstruction compared to the state of the art. We demonstrate the effectiveness of our method through several experiments on simulated and real off-axis holograms.

© 2013 Optical Society of America

1. Introduction

Given an intensity hologram that is acquired in a proper configuration, it is generally possible to reconstruct the complex field, commonly through demodulation, either in the temporal domain with the so-called phase-shifting approach [1], or in the spatial domain in an off-axis configuration, where complex information can be retrieved through Fourier filtering [2]. Several solutions have been proposed for reconstruction following these approaches [3, 4]. The corresponding algorithms are presently used in practice and incorporated in digital holographic microscopy for instance. Their common characteristic is that the spectral information that is actually used for reconstruction reduces to the first diffraction orders.

Besides direct Fourier-based reconstruction, iterative procedures that are based on more implicit formulations have been proposed for the recovery of complex fields from intensity data. This includes the well known Gerchberg-Saxton algorithm used for phase retrieval in its original form [5] as well as more evolved strategies [6]. The recent trend in the literature is to formulate the task of digital holographic reconstruction as an inverse problem. For instance, Cetin et al.[7] have applied total-variation-type regularization to reconstruct amplitude profiles based on a linear physical model. Zhang and Lam have followed a similar approach for optical scanning holography [8]. In particular, they provided reconstruction results using model-matched phantom experiments. In their work, Bardy et al.[9] and Marim et al.[10] have demonstrated the possibility of reconstructing amplitude profiles from undersampled holograms in the context of Gabor and off-axis holography, respectively. Along the same line of research, Rivenson et al.[11] have investigated distinct undersampling strategies. Note that, although the aforementioned methods yield promising results compared to standard reconstruction techniques, they either require the acquisition of two or more holographic planes or assume a linearized model for reconstruction. Meanwhile, Sotthivirat and Fessler [12] have developed an optimization-transfer method for reconstruction that involves a nonlinear model and that also handles the case of a single acquisition. Focusing on image-plane holography, their work has shown substantial improvements on simulated data compared to the conventional approaches used for that particular setting.

In this paper, we propose a new variational-reconstruction approach that can be applied to experimentally-acquired holographic data. Our technique accurately models the measurement system based on a nonlinear forward model and allows to recover the complex object field from one single off-axis intensity hologram. Compared to the state of the art, the main contribution of our reconstruction algorithm is to employ side information apart from the measured data to improve the quality of reconstruction. Specifically, we make prior assumptions on some shared object properties, and also use the knowledge of a reference correction hologram (RCH) to simplify the task of reconstruction and thus find suitable solutions. Another asset of our method is to exploit the complete image information at hand. The assumptions on the solution are used as a mean to regularize the amplitude and phase components of the complex field, which allows to perform reconstructions from spatially undersampled measurements. Note that an undersampling scheme was initially proposed in [13] for a particular off-axis protocol in low-level imaging conditions where the forward model can be linearized. However, unlike in our case, the protocol of [13] requires the complex amplitudes of the field to be first determined at the hologram plane through the combination of several intensity acquisitions. This is a step that can be avoided with our new nonlinear-model-based method presented in the sequel.

The paper is organized as follows. In Sect. 2, we describe the nonlinear model that corresponds to our own holographic setup, and introduce the relevant notations. In Sect. 3, we introduce our formulation of the reconstruction problem, and discuss the relevance of this new approach for obtaining high-quality solutions in practice. We then derive our reconstruction algorithm in Sect. 4. In Sect. 5, we conduct experiments on synthetic and real holograms, and compare our method with the state of the art from both qualitative and quantitative standpoints. We conclude our work in Sect. 6.

2. Forward model

The general holographic model that we consider is standard in the literature. It involves a coherent object field o(x;z) and a reference field r(x;z) that interfere to form the measured intensities

I(x)=|o(x;0)+r(x;0)|2.
The coordinate vector x = (x1, x2) corresponds to the spatial position parallel to the hologram plane, while the coordinate z denotes the signed distance from the hologram plane in the direction of light propagation. Note that the squared norm appearing in Eq. (1) makes the intensities I depend nontrivially on the field o. This will cause our own reconstruction problem to be non-convex, as further discussed in Sect. 3.

Following the MKS system of units, the physical parameters of the acquisition system consist in the wavelength λ, the camera pixel size Δx, the detector resolution N, the numerical aperture NA of the microscope objective, the object-space immersion medium ni, and the magnification factor M. Assuming a diffraction-limited configuration, the spatial bandwidth of the object corresponds to Δω0 = 2NAΔxN/(Mniλ) [14]. In the case of an off-axis hologram, a demodulated and non-aberrated object field can potentially be retrieved by employing an estimate r̃ of the reference wave r(x;0) employed for recording, up to a multiplicative factor α. One common way to get r̃ is to use a calibration measurement with an empty field of view which accounts for the phase terms induced by the optical system, denoted RCH [15].

We assume in this work that the holographic device measures an infinitely thin object profile. The field o is thus expected to be focused at some position z = −d, the corresponding planar profile being denoted as Ψd. While the in-focus propagation distance d is classically independent from the actual object-field reconstruction, it is relevant to our method because we shall introduce some prior knowledge on Ψd in Sect. 3. Accordingly, we consider Ψd as our unknown, and reformulate Eq. (1) as

I(x)=|(𝒫dΨd)(x)+r(x;0)|2.
The linear operators 𝒫d and denote digital propagation with distance d and spatial-frequency cut-off associated to the microscope objective (MO), respectively. These operators behave as a function of the parameters of the holographic setup, and their mathematical expressions are provided in Appendix A.

Finally, according to the pixel size Δx of the sensor array, the intensities I(x) are sampled to

I[m]=|(PdΨd)(x)+r(x;0)|2|x=mΔx,
and a binary mask w predefines what samples are actually measured according to some down-sampling factor D. Each sample I[m] is thus counted as a measurement if and only if the value w[m] is unity for the same m. The mask follows the pseudorandom structure of [13], and consists in binary values that are independent, identically distributed, and nonzero with probability 1/D. As shown in our experimental section, this generalized formulation successfully handles cases where the amount of data that is available for reconstruction is substantially reduced. Note that downsampled acquisitions potentially allow for faster imaging when using cameras where pixel values are addressed in random access [16].

3. Reconstruction problem

According to Eq. (3), the goal of holographic reconstruction is to produce the most precise estimate of the object wave Ψd at reconstruction distance d, given noisy and quantized versions of the intensity samples I. In this section, we propose an estimate Ψ̃d of Ψd that involves all the available measured data. Since the complex field Ψd corresponds to an unknown quantity for our reconstruction algorithm, it is treated as a variable in the sequel.

3.1. Consistent formulation

Following a variational framework, we define the estimated complex field Ψ̃d as the solution of an inverse problem. This solution must be matched to the available measurements, and has to be regularized so as to make the problem well posed.

According to Eq. (3), we define the class of suitable solutions as those yielding measurements that are compatible with the image data after reintroduction into our forward model. Accordingly, given the samples of I acquired by the sensor array, one can express a so-called consistency constraint on Ψd as

I(x)=|(𝒫dΨd)(x)+αr¯(x)|2,xmΔxwithw[m]=1.
The function r̄ corresponds to a normalized version of r̃ whose phases are left untouched and whose amplitudes are replaced by their mean pointwise. This allows the actual reference amplitudes to be fully determined by the factor α obtained through the reconstruction process. Due to the presence of noise, we relax the strict data-fidelity constraint of Eq. (4) and rewrite it as
𝒟(Ψd,α)=m2w[m](|(𝒫dΨd)(x)+αr¯(x)|2I(x))2|x=mΔx<ε.
Besides Ψd, the multiplicative factor α influences the value of 𝒟 in Eq. (5); both quantities shall be alternately optimized in our algorithm, as described in Sect. 4. Finally, the positive scalar ε determines to what extent the reconstructed solution can depart from the available measurements, depending on the level of noise that affect them. Based on the sampled intensities, this constant can be deduced from the signal-to-noise ratio (SNR) of the acquisition device εSNR through the relation
ε=exp(εSNRln(10)/10)m2w[m]I2(x)|x=mΔx,
as shown in Appendix B. Under the sole data-fidelity constraint of Eq. (5), the solution Ψ̃d is underdetermined. In order to make the problem well posed, we regularize Ψd by minimizing the total variation (TV) [17] of its phase and amplitude components separately. From a qualitative standpoint, our regularization approach is derived from the assumption that the unwrapped-phase and amplitude maps of the focused hologram are both well approximated by piecewise-constant functions as they are in-focus. This assumption shall prove fruitful for improving the reconstruction quality, as shown in the experimental section. Note also that our reconstruction algorithm proposed in Sect. 4 does not require any phase unwrapping to regularize the solution. Denoting the phase and amplitude TV values by TVϕd) and TVρd), respectively, we write
TVϕ(Ψd)=2(arguΨd)(x)dx,
TVρ(Ψd)=2|Ψd|(x)dx,
where argu Ψd denotes the unwrapped phase of Ψd. The above thus corresponds to the constrained minimization problem
Ψ˜d=argminΨd[TVϕ(Ψd)+γTVρ(Ψd)]s.t.𝒟˜(Ψd)<ε,
where γ is a positive scalar weight that balances the respective influence of phase and amplitude regularization, and where
𝒟˜(Ψd)=argminα𝒟(Ψd,α).
Given the structure of Eq. (1), the expression of 𝒟d, α) in Eq. (5) contains two imbricated squaring operations. This implies that 𝒟̃(Ψd) is a non-convex functional of the field Ψd. The problem of Eq. (9) can thus have several local optima; its solution is non-unique in general.

3.2. Discretization

Using vector notation, we discretize the unknowns and the aforementioned operators according to the pixel size Δx. Each vector contains the components of the corresponding discretized map in lexicographic order. Accordingly, the discrete solution is defined as

Ψ˜d=argminΨd[TVϕ(Ψd)+γTVρ(Ψd)]s.t.𝒟˜(Ψd)<ε.
The corresponding data term 𝒟 takes the form
𝒟(Ψd,α)=iwi(|(AΨd)i+αr¯i|2yi)2,
where the vectors w, Ψd, r̄, and y correspond to the sampled versions of w, Ψd, r̄, and I, respectively. The matrix A corresponds to the successive application of the linear operators and 𝒫d in discretized form. The discretized regularization terms read
TVϕ(Ψd)=i[(R1(arguΨd))i2+(R2(arguΨd))i2+ν]1/2,
TVρ(Ψd)=i[(R1|{Ψd}|)i2+(R2|{Ψd}|)i2+ν]1/2,
where the R1,2 are convolution matrices implementing the discretized first derivatives in the horizontal and vertical directions, respectively, and where ν is an additive parameter that ensures the differentiability of both TV functionals [18]. The horizontal and vertical filters associated with R1,2 are each one-dimensional, and correspond to ( z1,211) in the 𝒵 -transform domain. The symbol |{·}| denotes the absolute value applied to a vector pointwise. Similarly, argu takes the unwrapped argument of the vector components pointwise.

4. Reconstruction algorithm

We solve the problem of Eq. (11) based on a gradient-descent approach that we combine with a specifically devised line-search strategy. Specifically, we propose to apply the following iterative procedure:

  1. Initialize the solution as Ψd = 0.
  2. Find the optimal reference factor α̃ = argminα𝒟(Ψd, α) to determine 𝒟̃(Ψd).
  3. If 𝒟̃(Ψd) > ε:
    1. Compute the data-term gradient ∇𝒟̃ = ∂𝒟̃ (Ψd)/Ψd.
    2. Compute the step τ̃ ∈ ℝ+ s.t. 𝒟̃(Ψdτ̃∇𝒟̃) = ε if it exists; optimize τ̃ = argminτ𝒟̃(Ψdτ𝒟̃) otherwise.
    3. Update the solution as ΨdΨdτ̃∇𝒟̃.
  4. If 𝒟̃(Ψd) ≤ ε:
    1. Compute the phase-regularization gradient ∇TVϕ = TVϕ(Ψd)/Ψd.
    2. Compute the amplitude-regularization gradient ∇TVρ = TVρ(Ψd)/Ψd.
    3. Update the solution as ΨdΨdτreg (∇TVϕ + γ∇TVρ).
  5. Return to Step 2 until V iterations are reached.

The positive scalar τreg is a step parameter that determines the strength of the regularization flow. The corresponding regularization gradients ∇TVϕ, ∇TVρ used in Step 4 and the data-term gradient ∇𝒟̃ used in Step 3 are complex-valued with ∇· = · / Re(Ψd) + i( · /Im(Ψd)). The explicit forms of the update terms are then obtainable from Eqs. (12), (13), and (14) using differential calculus. Our regularization technique does not require any explicit phase unwrapping, as explained in Appendix C. Finally, determining either of the two scalars α̃ and τ̃ in Steps 2 and 3(b) amounts to solving an elementary line-search problem. This issue is addressed in more detail in Appendix D.

Since the problem of Eq. (11) is non-convex, the solution obtained with our algorithm potentially depends on the initial estimate. In practice, however, no significant changes were observed, hence our generic choice of zero initialization in Step 1. While rigorous convergence analysis is nontrivial in non-convex settings [19], experimental evidence suggests that our algorithm converges to a fixed point. Indeed, the solution and the corresponding value of the regularization functional in Eq. (11) have been observed to stabilize as the iterations proceed. Each iteration of our algorithm is comparable to the standard technique in terms of computational cost.

5. Experiments

Based on the above developments, we conduct experiments in simulated and practical off-axis configurations. The holographic reconstructions are performed with our method as well as with the conventional approach consisting in demodulation followed by back-propagation. Besides the standard reconstruction paradigm, we investigate cases where the amount of available measurements is reduced through random subsampling. In the sequel, a first set of experiments on simulated holograms evaluates the algorithms quantitatively in each case, using the oracle complex fields for comparison. The experiments that follow allow to determine the relevance and the competitiveness of our approach when real data are involved.

The standard and proposed algorithms are implemented in MATLAB, using the angular-spectrum approximation to model light propagation [20]. The reconstructed objects are shown with the same resolution and field of view as those of the measured holograms. In practical configurations, an additive linear-phase field can appear in the reconstruction when tilt coming from sample positioning occurs between acquisition of the RCH and the one of the object. This additive field is removed to simplify the visual interpretation of the results. Note that our algorithm is set to regularize the reconstructions on a slightly larger field of view than the one displayed in order to minimize the effects of border conditions during digital propagation. As shown in the sequel, this improves the quality of reconstruction at the boundaries compared to the standard technique that uses apodization.

5.1. Synthetic data

In this set of experiments, the original complex fields are available, their amplitudes and phases at the focus plane being defined from the two-dimensional spatial maps of Table 1. The Pentagon picture originates from the BM3D database at http://www.cs.tut.fi/ foi/GCF-BM3D/ and the others from the USC-SIPI database at http://sipi.usc.edu/database/. After normalization of its spatial-domain amplitude and phase values to [0, 1] and [0, π], respectively, each complex field is used to generate a distinct 512 × 512 hologram according to our forward model. These holograms were first obtained from larger objects and then restricted to a limited field of view—centered and with 512 × 512 samples—in analogy with a physical setup. The acquisition parameters that are used for hologram generation correspond to an off-axis configuration, choosing the wavelength as corresponding to the He-Ne laser line, i.e., λ = 633 · 10−9, and choosing Δx = 6.45 · 10−6, N = 1,024, NA = 0.25, ni = 1, M = 10, d = 4 · 10−2, and r ̃(x) = exp(2πj(k1 + k2)/5Δx). The oracle maps mentioned in Table 1 are shown along with the corresponding holograms in Fig. 1. Note that the scales at which our object maps, holograms, and reconstructions are displayed are normalized in each figure for convenience.

 figure: Fig. 1

Fig. 1 Full-sized oracle maps and corresponding hologram acquisitions used in the synthetic experiments. From (A) to (D): Pentagon, Man, Airplane, and Airport maps defining the objects of Table 1. From (E) to (I): intensity holograms of size 512 × 512 associated with the objects no. 1, 2, 3, 4, and 5.

Download Full Size | PDF

Tables Icon

Table 1. Reconstruction Quality in Synthetic Experiments

Given the synthetic holograms and their corresponding parameters, our goal is to determine how accurately the oracle complex fields can be reconstructed by the standard and proposed methods. For each hologram, we consider the classical setting where all samples are available for reconstruction, as well as the setting D = 2 where only 50% of the data are kept through random downsampling. As reconstruction parameters, our algorithm uses ν = 10−10, τreg = 1.5 · 10−2, εSNR = 35, and V = 2,000 in all synthetic experiments. The regularization weight is set to γ = 0.5. In order to neglect the possible influence of boundary artifacts, the SNR values are evaluated on fields of view that are centered and whose size is reduced by 200 pixels in each dimension compared to those of the corresponding holograms.

The numerical results that are obtained are reported in Tables 1 and 2. As an example, the reconstruction is shown for the phase-only hologram #5 in Fig. 2; the phase values [0, π] are mapped to the grayscale range [black, white]. In the classical paradigm where all samples are available, the values of Table 1 indicate an overall improvement in terms of reconstruction quality when using our approach. Our visual results also demonstrate that less artifacts are produced by the proposed method at the object boundaries. Furthermore, besides backward compatibility with the standard approach, our reconstruction algorithm is able to successfully recover object profiles when the amount of sampled data is reduced. In particular, the corresponding SNR values of Table 1 remain stable as compared to the fully sampled configurations. By contrast, the SNR values that correspond to the standard reconstruction approach decrease dramatically when downsampling is used. The strong aliasing effects that are due to downsampling can indeed not be handled properly in that case since no regularization is used.

 figure: Fig. 2

Fig. 2 Reconstruction of the phase-only Man object from the synthetic hologram #5 without downsampling (D = 1) and with downsampling (D = 2).

Download Full Size | PDF

Tables Icon

Table 2. Effective Field of View in Synthetic Experiments

In Table 2, the proposed field-of-view (FOV) ratio measures the relative area of the solution where the reconstruction quality remains highest with respect to the oracle, considering the downsampled case only for our method. This area is determined through local MSE computations performed at every distance from the boundaries of the given object. It is obtained as the centered region where this local MSE does not exceed twice the global MSE. The FOV results further document the improvement of our method over the standard one. Our algorithm yields data-dependent values because its regularization strategy is itself data-adaptive.

At this stage, these observations confirm the relevance of our approach for holographic reconstruction in both classical and downsampled configurations. They remain to be corroborated in the real-data experiments that follow.

5.2. Real data

In this second experimental part, we consider holograms that are acquired practically from distinct physical objects. The 1,024 × 1,024 hologram Neuron and the 512 × 512 hologram Epithelial consist in acquisitions of neural and epithelial cell samples, respectively. The 512 × 512 holograms USAF 5-4 and USAF 9-8 correspond to USAF targets with mixed phase and amplitude information. As in the synthetic case of Sect. 5.1, the acquisition settings that have been used in the optical setup are off-axis. The parameters that are common to all acquisitions are Δx = 6.45 · 10−6 and ni = 1. The other scalar parameters are reported in Table 3, while the approach followed to generate the reference waves is described in [14].

Tables Icon

Table 3. Object-Dependent Parameters Used for Optical Acquisition

The proposed algorithm is parameterized as above, except τreg = 3 · 10−4, γ = 1 for the hologram USAF 9-8, and except εSNR. In this practical setting, the latter quantity depends nontrivially on the acquisition conditions and is determined experimentally. The values of εSNR specific to each hologram are thus reported in Table 3 along with the acquisition parameters. Moreover, since the amplitude of the reference wave tends to vary slowly in space due to the non-ideality of the holographic device, we propose to determine α in Eq. (5) in a spatially adaptive manner. Specifically, we perform one distinct optimization for each 8 × 8-pixel block of the solution. Note that each optimized value of α is still estimated based on the principles of Appendix D, although the sums involved in the computations are now restricted to the corresponding block.

The reconstructions from Neuron, Epithelial, and USAF 5-4 are shown in Figs. 3, 4, 5, and 6. While the standard technique only handles the fully sampled scenario, distinct downsampling factors D = 1,2,4 are considered with our algorithm. When all samples are available, the visual results show similar advantages of our algorithm over the classical method as in the synthetic case. For instance, the field of view of the reconstructed objects is extended when using our iterative approach, unlike in the standard approach where the signal must be attenuated at borders in order to properly handle the boundary conditions. These results thus demonstrate that the boundary issues arising in digital propagation can be suitably addressed when using an implicit formulation of the reconstruction problem such as Eq. (11). In this set of real experiments, our method is also observed to reduce the amount of noise in the recovered phase and amplitude profiles. In particular, the continuity of their background is improved. When increasing D, the reconstruction quality remains acceptable even though some fine details are inevitably lost due to the associated decrease in resolution. This decrease is associated to small patterns of “pointillism” type appearing in the holograms.

 figure: Fig. 3

Fig. 3 Reconstruction from the phase-only hologram Neuron.

Download Full Size | PDF

 figure: Fig. 4

Fig. 4 Reconstruction from the phase-only hologram Epithelial.

Download Full Size | PDF

 figure: Fig. 5

Fig. 5 Reconstructed amplitudes from the hologram USAF 5-4.

Download Full Size | PDF

 figure: Fig. 6

Fig. 6 Reconstructed phases from the hologram USAF 5-4.

Download Full Size | PDF

The contrast and the sharpness of object features are also well preserved when using our algorithm. In particular, the results on USAF 5-4 shown in Figs. 5 and 6 demonstrate that sharp transitions and large regular structures are reconstructed accurately. The resolution performance of the standard and proposed methods can be compared in Figure 7 where the reconstruction from USAF 9-8 is shown in the fully sampled setting. In this example, the resolution capability of both methods is similar. While the use of TV regularization produces a slight smearing effect in some places, the level of noise is strongly reduced in the reconstruction obtained with our method. This set of practical experiments thus corroborates the observations made in the synthetic case, and validates our approach as applied to real off-axis holograms subject to noise.

 figure: Fig. 7

Fig. 7 Reconstructed amplitudes and phases from the fully sampled hologram USAF 9-8.

Download Full Size | PDF

6. Conclusions

We have devised an algorithm for off-axis holographic reconstruction that is based on a consistent problem formulation. Based on suitable regularity assumptions, our technique has allowed to reconstruct complex-valued object profiles satisfactorily, including in the case where the amount of measured samples is decreased. Compared to the standard technique, the proposed method avoids the presence of boundary artifacts in the solutions, and reduces the perceived level of noise in practical configurations. From a general perspective, the obtained results further illustrate the interest of inverse-problem approaches in holography.

A. Definition of and 𝒫d

The operator corresponds to an amplitude filter defined in the Fourier domain as

H^(ω)={1,ωΔωo/20,ω>Δωo/2.
We remind that this frequency limitation fundamentally comes from the NA of the objective. Given the distance d and the wavenumber k = 2π/λ, 𝒫d corresponds to a phase filter in the Fourier domain that is defined as
P^d(ω)=exp(ikd(1λ2ω2)1/2).
This propagation model is based on the angular-spectrum method [20].

B. Data-fidelity constant

As mentioned in Sect. 3, we constrain suitable estimates of the unknown complex field to be consistent with the measured intensity samples. According to our forward model, a given solution Ψd should thus yield intensities of the form

I˜(x)=|(𝒫dΨd)(x)+αr¯(x)|2
that are close to I(x), ∀xmΔx with w[m] = 1. Assuming that the noise level affecting the samples of I corresponds to a SNR of εSNR in dB, we constrain the similarity between Ĩ and I to match this value as a lower limit. Taking only the non-masked samples into account, we obtain
εSNR<10log10(m2w[m]I(x)2|x=mΔxm2w[m](I˜(x)I(x))2|x=mΔx),
which expands as
m2w[m](I˜(x)I(x))2|x=mΔx<exp(εSNRln(10)/10)m2w[m]I2(x)|x=mΔx.
This expression is equivalent to Eq. (5). The left-hand side of Eq. (19) corresponds to the data term 𝒟d, α), while its right-hand side is the definition of the data-fidelity constant ε in Eq. (6). Note that this derived result is specific to the quadratic form of Eq. (5) that is matched to a Gaussian noise model. Another noise model (e.g., Poisson) would correspond to a distinct formula for ε.

C. Spatial derivatives of unwrapped phase

The phase derivative can be estimated independently of wrapping [21]. More specifically, since our finite-difference filters only consider two adjacent phases at the same time, only the relative phase differences are relevant. Assuming that the phase jumps cannot exceed π between two pixels, and defining vi = ((R1,2(argΨd))i mod 2π), we have the relation

(R1,2(arguΨd))i={vi,vi<πvi2π,viπ.
The above assumption on maximal phase jumps corresponds to the required properties for proper phase sampling [21]. In settings such as [22] where this assumption is violated, our method may still be applied, provided that only the solution amplitudes are regularized. While preliminary investigations indicate satisfactory results on objects whose phases are random and uniformly distributed in [−π, π[, this issue remains to be addressed in further research.

D. Line search

Each of the optimization problems that arise in our algorithm when determining α̃ and τ̃ can be recast as the resolution of a simple polynomial equation. According to Eq. (5), both data-term quantities appearing in Steps 2 and 3(b) are written in the generic form

𝒞(χ)=m2(|a[m]+χb[m]|2c[m])2,
where a and b, and c are known constant sequences—the latter being real-valued and nonnegative—and where χ is the scalar variable of interest to optimize. Expanding the terms of Eq. (21), and factoring out χ, we obtain the simplified expression
𝒞(χ)=Aχ4+Bχ3+Cχ2+Dχ+E,
where
A=m2|b[m]|4,B=m24Re{a[m]b*[m]}|b[m]|2,C=m22|b[m]|2(|a[m]|2c[m])+4(Re{a[m]b*[m]})2,D=m24Re{a[m]b*[m]}(|a[m]|2c[m]),E=m2(|a[m]|2c[m])2.
Given their definition, the coefficients A, B, C, D, and E are real scalars that only depend on a, b, and c. Depending on the optimization problem, 𝒞 has to be either equalized to ε, or minimized with respect to the scalar parameter. In the first case, the solution—if it exists—corresponds to the smallest real and nonnegative root χi of the fourth-degree polynomial
Aχ4+Bχ3+Cχ2+Dχ+(Eε).
In the second case, the argument minimizing 𝒞 must cancel the derivative of 𝒞 with respect to χ. Accordingly, the solution is the real and nonnegative root of the polynomial
4Aχ3+3Bχ2+2Cχ+D
that yields the smallest value of 𝒞.

Acknowledgments

This work was supported by the European Commission under Grant ERC-2010-AdG 267439-FUN-SP.

References and links

1. E. Cuche, P. Marquet, and C. Depeursinge, “Simultaneous amplitude-contrast and quantitative phase-contrast microscopy by numerical reconstruction of Fresnel off-axis holograms,” Appl. Opt. 38(34), 6994–7001 (1999). [CrossRef]  

2. I. Yamaguchi and T. Zhang, “Phase-shifting digital holography,” Opt. Lett. 22(16), 1268–1270 (1997). [CrossRef]   [PubMed]  

3. G. Popescu, T. Ikeda, R. Dasari, and M. Feld, “Diffraction phase microscopy for quantifying cell structure and dynamics,” Opt. Lett. 31(6), 775–777 (2006). [CrossRef]   [PubMed]  

4. Y. Awatsuji, T. Tahara, A. Kaneko, T. Koyama, K. Nishio, S. Ura, T. Kubota, and O. Matoba, “Parallel two-step phase-shifting digital holography,” Appl. Opt. 47(19), D183–D189 (2008). [CrossRef]  

5. R. W. Gerchberg and W. O. Saxton, “A practical algorithm for the determination of phase from image and diffraction plane pictures,” Optik 35(2), 227–246 (1972).

6. V. Katkovnik, A. Migukin, and J. Astola, “3D wave field reconstruction from intensity-only data: variational inverse imaging techniques,” in 9th Euro-American Workshop on Information Optics (Helsinki, Finland, 2010).

7. M. Cetin, W. Karl, and A. Willsky, “Edge-preserving image reconstruction for coherent imaging applications,” in Proceedings of ICIP 2002 International Conference on Image Processing (Rochester, NY, USA, 2002).

8. X. Zhang and E. Lam, “Edge-preserving sectional image reconstruction in optical scanning holography,” J. Opt. Soc. Am. A 27(7), 1630–1637 (2010). [CrossRef]  

9. D. Brady, K. Choi, D. Marks, R. Horisaki, and S. Lim, “Compressive holography,” Opt. Express 17(15), 13,040–13,049 (2009).

10. M. Marim, M. Atlan, E. Angelini, and J. Olivo-Marin, “Compressed sensing with off-axis frequency-shifting holography,” Opt. Lett. 35(6), 871–873 (2010). [CrossRef]   [PubMed]  

11. Y. Rivenson, A. Stern, and B. Javidi, “Compressive fresnel holography,” J. Disp. Technol. 6(10), 506–509 (2010). [CrossRef]  

12. S. Sotthivirat and J. Fessler, “Penalized-likelihood image reconstruction for digital holography,” J. Opt. Soc. Am. A 21(5), 737–750 (2004). [CrossRef]  

13. M. Marim, E. Angelini, J.-C. Olivo-Marin, and M. Atlan, “Off-axis compressed holographic microscopy in low-light conditions,” Opt. Lett. 36(1), 79–81 (2011). [CrossRef]   [PubMed]  

14. N. Pavillon, C. Seelamantula, J. Kuhn, M. Unser, and C. Depeursinge, “Suppression of the zero-order term in off-axis digital holography through nonlinear filtering,” Appl. Opt. 48(34), H186–H195 (2009). [CrossRef]  

15. T. Colomb, J. Kühn, F. Charrière, C. Depeursinge, P. Marquet, and N. Aspert, “Total aberrations compensation in digital holographic microscopy with a reference conjugated hologram,” Opt. Express 14(10), 4300–4306 (2006). [CrossRef]   [PubMed]  

16. M. Bigas, E. Cabruja, J. Forest, and J. Salvi, “Review of CMOS image sensors,” Microelectron. J. 37(5), 433–451 (2006). [CrossRef]  

17. L. Rudin, S. Osher, and E. Fatemi, “Nonlinear total variation based noise removal algorithms,” Phys. D 60, 259–268 (1992). [CrossRef]  

18. T. Chan and P. Mulet, “On the convergence of the lagged diffusivity fixed point method in total variation image restoration,” SIAM J. Numer. Anal. 36(2), 354–367 (1999). [CrossRef]  

19. O. Scherzer, ed., Energy Minimization Methods (Springer, 2011).

20. J. W. Goodman, Introduction to Fourier Optics, 2nd ed. (McGraw Hills Companies, Inc., 1996).

21. M. Guizar-Sicairos, A. Diaz, M. Holler, M. Lucas, A. Menzel, R. Wepf, and O. Bunk, “Phase tomography from x-ray coherent diffractive imaging projections,” Opt. Express 19(22), 345–,357 (2011). [CrossRef]  

22. K. Choi, R. Horisaki, J. Hahn, S. Lim, D. L. Marks, T. J. Schulz, and D. J. Brady, “Compressive holography of diffuse objects,” Appl. Opt. 49(34), H1–H10 (2010). [CrossRef]  

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (7)

Fig. 1
Fig. 1 Full-sized oracle maps and corresponding hologram acquisitions used in the synthetic experiments. From (A) to (D): Pentagon, Man, Airplane, and Airport maps defining the objects of Table 1. From (E) to (I): intensity holograms of size 512 × 512 associated with the objects no. 1, 2, 3, 4, and 5.
Fig. 2
Fig. 2 Reconstruction of the phase-only Man object from the synthetic hologram #5 without downsampling (D = 1) and with downsampling (D = 2).
Fig. 3
Fig. 3 Reconstruction from the phase-only hologram Neuron.
Fig. 4
Fig. 4 Reconstruction from the phase-only hologram Epithelial.
Fig. 5
Fig. 5 Reconstructed amplitudes from the hologram USAF 5-4.
Fig. 6
Fig. 6 Reconstructed phases from the hologram USAF 5-4.
Fig. 7
Fig. 7 Reconstructed amplitudes and phases from the fully sampled hologram USAF 9-8.

Tables (3)

Tables Icon

Table 1 Reconstruction Quality in Synthetic Experiments

Tables Icon

Table 2 Effective Field of View in Synthetic Experiments

Tables Icon

Table 3 Object-Dependent Parameters Used for Optical Acquisition

Equations (25)

Equations on this page are rendered with MathJax. Learn more.

I ( x ) = | o ( x ; 0 ) + r ( x ; 0 ) | 2 .
I ( x ) = | ( 𝒫 d Ψ d ) ( x ) + r ( x ; 0 ) | 2 .
I [ m ] = | ( P d Ψ d ) ( x ) + r ( x ; 0 ) | 2 | x = m Δ x ,
I ( x ) = | ( 𝒫 d Ψ d ) ( x ) + α r ¯ ( x ) | 2 , x m Δ x with w [ m ] = 1.
𝒟 ( Ψ d , α ) = m 2 w [ m ] ( | ( 𝒫 d Ψ d ) ( x ) + α r ¯ ( x ) | 2 I ( x ) ) 2 | x = m Δ x < ε .
ε = exp ( ε SNR ln ( 10 ) / 10 ) m 2 w [ m ] I 2 ( x ) | x = m Δ x ,
TV ϕ ( Ψ d ) = 2 ( arg u Ψ d ) ( x ) d x ,
TV ρ ( Ψ d ) = 2 | Ψ d | ( x ) d x ,
Ψ ˜ d = arg min Ψ d [ TV ϕ ( Ψ d ) + γ TV ρ ( Ψ d ) ] s . t . 𝒟 ˜ ( Ψ d ) < ε ,
𝒟 ˜ ( Ψ d ) = arg min α 𝒟 ( Ψ d , α ) .
Ψ ˜ d = arg min Ψ d [ TV ϕ ( Ψ d ) + γ TV ρ ( Ψ d ) ] s . t . 𝒟 ˜ ( Ψ d ) < ε .
𝒟 ( Ψ d , α ) = i w i ( | ( A Ψ d ) i + α r ¯ i | 2 y i ) 2 ,
TV ϕ ( Ψ d ) = i [ ( R 1 ( arg u Ψ d ) ) i 2 + ( R 2 ( arg u Ψ d ) ) i 2 + ν ] 1 / 2 ,
TV ρ ( Ψ d ) = i [ ( R 1 | { Ψ d } | ) i 2 + ( R 2 | { Ψ d } | ) i 2 + ν ] 1 / 2 ,
H ^ ( ω ) = { 1 , ω Δ ω o / 2 0 , ω > Δ ω o / 2.
P ^ d ( ω ) = exp ( i k d ( 1 λ 2 ω 2 ) 1 / 2 ) .
I ˜ ( x ) = | ( 𝒫 d Ψ d ) ( x ) + α r ¯ ( x ) | 2
ε S N R < 10 log 10 ( m 2 w [ m ] I ( x ) 2 | x = m Δ x m 2 w [ m ] ( I ˜ ( x ) I ( x ) ) 2 | x = m Δ x ) ,
m 2 w [ m ] ( I ˜ ( x ) I ( x ) ) 2 | x = m Δ x < exp ( ε SNR ln ( 10 ) / 10 ) m 2 w [ m ] I 2 ( x ) | x = m Δ x .
( R 1 , 2 ( arg u Ψ d ) ) i = { v i , v i < π v i 2 π , v i π .
𝒞 ( χ ) = m 2 ( | a [ m ] + χ b [ m ] | 2 c [ m ] ) 2 ,
𝒞 ( χ ) = A χ 4 + B χ 3 + C χ 2 + D χ + E ,
A = m 2 | b [ m ] | 4 , B = m 2 4 Re { a [ m ] b * [ m ] } | b [ m ] | 2 , C = m 2 2 | b [ m ] | 2 ( | a [ m ] | 2 c [ m ] ) + 4 ( Re { a [ m ] b * [ m ] } ) 2 , D = m 2 4 Re { a [ m ] b * [ m ] } ( | a [ m ] | 2 c [ m ] ) , E = m 2 ( | a [ m ] | 2 c [ m ] ) 2 .
A χ 4 + B χ 3 + C χ 2 + D χ + ( E ε ) .
4 A χ 3 + 3 B χ 2 + 2 C χ + D
Select as filters


Select Topics Cancel
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.