Mean Gradient Descent: An optimization approach for single-shot interferogram analysis
Sunaina Mansi Butola Kedar Khare




 

0.2



A New Littlest Seesaw Model


    
Ping-Tao Chen^1[E-mail: chenpt@mail.ustc.edu.cn],  Gui-Jun Ding^1[E-mail: dinggj@ustc.edu.cn],  Stephen F. King^2[E-mail: king@soton.ac.uk], Cai-Chang Li^1,3[E-mail: lcc0915@mail.ustc.edu.cn]  





^1Interdisciplinary Center for Theoretical Study and  Department of Modern Physics,

University of Science and Technology of China, Hefei, Anhui 230026, China

^2Physics and Astronomy,
University of Southampton,
Southampton, SO17 1BJ, U.K.

^3School of Physics, 
Northwest University, 
Xi'an 710127, China






    
===================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================




Complex object wave recovery from single-shot interference pattern is an important practical problem in interferometry and digital holography. The most popular single-shot interferogram analysis method involves Fourier filtering of cross-term but this method suffers from poor resolution. For obtaining full pixel resolution, it is necessary to model the object wave recovery as an optimization problem. The optimization approach typically involves minimizing a cost function consisting of a data consistency term and one or more constraint terms. Despite its potential performance advantages, this method is not used widely due to several tedious and difficult tasks such as empirical tuning of free parameters. We introduce a new optimization approach “Mean gradient descent (MGD)” for single-shot interferogram analysis that is simple to implement, robust and does not require any free parameters. The MGD iteration does not try to achieve minimization of any cost function but instead aims to reach a solution point where the data consistency and the constraint terms balance each other. This is achieved by iteratively progressing the solution in the direction that bisects the descent directions associated with the error and constraint terms. Numerical illustrations are shown for recovery of a step phase object from its corresponding off-axis as well as on-axis interferograms simulated with multiple noise levels. Our results show full pixel resolution as evident from the recovery of phase step and excellent rms phase accuracy relative to the ground truth phase map. The concept of MGD as presented here can potentially find applications to wider class of optimization problems. 









§ INTRODUCTION

Demodulation of interference fringes is of importance to a wide ranging applications in optics <cit.> including optical metrology, digital holography for live cell imaging, surface inspection, particle-field holography, astronomical imaging etc. Currently the interference fringes in all these applications are predominantly recorded on CCD or CMOS array sensors that are readily available. The digitally recorded fringe pattern is then numerically processed for recovering the complex-valued object wave of interest. Traditionally there are two main methods that are used for interferogram analysis. For off-axis interference pattern the dc and the cross-terms in the interference pattern can be nominally separated in Fourier space and thus the cross term can be obtained by Fourier space filtering <cit.>. This method requires single interference pattern but as explained later, the Fourier space filtering operation inherently causes loss of resolution. Phase shifting interferometry <cit.> on the other hand can achieve full pixel resolution using a multi-shot interferogram recording approach and requires stringent vibration isolation. The CCD/CMOS arrays available today have good sensitivity so that just a milli-second of exposure (during which the recording is not affected much by ambient vibrations) is often sufficient to record good quality interferograms with nominal tabletop interferometer setups. In view of this it is highly desirable to have a practical algorithm that can demodulate a single-shot interferogram without compromising on resolution and accuracy.

In this work we are primarily interested in complex-valued object recovery from a single-shot image plane digital hologram/interferogram represented as:

    H = |R|^2 + |O|^2 + RO^* + R^*O,

Here the hologram H, the reference beam R = |R|e^iϕ_R and the object beam O = |O|e^iϕ_O are two-dimensional functions of the pixel coordinates (x,y) at the sensor plane. Further we will assume that the reference beam R (both its magnitude |R| and phase ϕ_R) is known and the problem is to recover the object wave function O from the single-shot hologram. We observe that since we are trying to fit two functions |O(x,y)| and ϕ_O (x,y) to a single hologram frame H(x,y), the problem does not have a unique solution and additional constraint is required in order to obtain a desirable solution. For example in the present case of image plane holography, it is expected that the object of interest to be imaged will typically have sharp edges and hence constraints like minimal Total Variation (TV) may be applied to the object field O(x,y). Such constraints can be included effectively if the object wave recovery is modeled as an optimization problem.  

The first optimization framework for interferogram analysis may be considered to be the regularized phase tracking (RPT) method<cit.>. In this approach the amplitude and phase of an unknown object are separately recovered pixel by pixel by fitting a local polynomial model for the unknown object wave after the removal of low frequency intensity terms from the interferogram data. In another work on complex wave retrieval algorithm <cit.> the phase and amplitude of the object are separately recovered by local least-square fitting with a variable sized window after changing the non-linear variable equation of hologram to a set of linear equations. A nominal constrained optimization approach using the complex (or Wirtinger) derivatives <cit.> for functional gradients was proposed in <cit.>. This approach demonstrated recovery of object wave from a single-shot off-axis hologram even when the dc and the cross terms showed significant overlap in Fourier space. This method was further used to demonstrate highly accurate phase recovery from low light level interference data <cit.>. An alternating amplitude and phase optimization strategy was developed in <cit.> which is an interesting approach that again demonstrated resolution improvement over the Fourier transform method. The optimization methodology in <cit.> was modified further in <cit.> where an adaptive optimization approach was proposed. In the optimization framework the problem of reconstruction of the object wave O typically modelled as minimization of a two-term cost function:

    C(O,O^*)    = ||H - |O+R|^2||^2_2 + α  TV(O,O^*) 
       = C_err(O,O^*) + α  C_TV(O,O^*),

Here C, C_err, C_TV are functions of the hologram H, the reference beam R and the unknown complex object wave O. The first term C_err  refers to the L2-norm squared error or data-consistency term and C_TV  refers to the TV penalty. The positive number α determines the weight between two terms of the cost function. The definition of TV we use for the present work is<cit.>:

    TV(O,O^*) = ∑_i= all pixels [   |∇_x O_i|   + |∇_y O_i|   ]  ,

where, ∇_x and ∇_y are the x and y partial derivative operators respectively. In <cit.>, it was shown that starting with any initial guess O^(0) for the solution, if the cost function is iteratively reduced via a gradient descent scheme, the quality of the resultant solution depends on the numerical value of α. The parameter α therefore needs to be determined empirically which is often not desirable. 

In the present work we propose a novel approach - Mean Gradient Descent (MGD) - that does not require any weight parameter α. The aim of MGD is not to achieve minimum of any cost function as in Eq. (<ref>) but to instead achieve balance between the two terms C_err and C_TV. Our methodology is inspired by a successful image reconstruction algorithm ASD-POCS <cit.> for X-ray computed tomography. In an earlier work <cit.>, a methodology very similar to ASD-POCS  was employed for the single-shot interferogram analysis. The main idea in ASD-POCS is to reach a solution point where the directions - ∇_O^∗ C_err and - ∇_O^∗ C_TV corresponding to descent directions for the reduction in C_err and C_TV make an obtuse angle. The descent directions corresponding to the two terms of the cost function can thus be thought to achieve an equilibrium. ASD-POCS achieves this balance by alternating minimization of C_err and C_TV. Starting with O = O^(n) the error term C_err is reduced leading to an intermediate solution O^(n)_int. The C_TV associated with this intermediate solution is then recursively reduced to obtain the next guess O^(n+1) such that the distances d_1 = ||O^(n) - O^(n)_int||_2  and d_2 = ||O^(n+1) - O^(n)_int||_2 are approximately matched in every iteration. Once the solution reaches near the optimal point, ensuring that d_1 ≈ d_2 in further iterations causes negligible change in the solution as the two descent directions oppose each other. 

We show here that reaching the balance point as in ASD-POCS is possible without employing an alternating minimization scheme but by iteratively progressing the solution in a direction that bisects that of the two functional gradients. The MGD iteration is computationally much simpler as compared to the alternating minimization approach. While MGD is discussed here in the context of single-shot interferogram analysis, we believe that it may be useful to two term optimization problems in general. In the context of interferogram analysis we show that procedure is uniformly applicable to various interferometric configurations (off-axis as well as on-axis) and thus amenable to be used with multiple digital holography/ interferometry system configurations.

The paper is organized as follows. In section 2 we provide a detailed description of MGD algorithm with the help of numerical illustrations for different noise levels present in the interferogram data corresponding to off-axis plane reference beam and step phase object. In section 3, we demonstrate the performance of MGD iteration for on-axis and near on-axis spherical reference beam configurations of the interferogram setup. Finally,in section 4, we summarize our observations and provide our insights over the results obtained by MGD algorithm.



§ DESCRIPTION OF MGD ITERATION

For simplicity we describe the MGD iteration with an illustration of single-shot off-axis interferogram.  For object wave in the interferogram plane, we use a square-shaped binary phase object on a 500 × 500 pixel grid. We have taken unit amplitude across the entire object and a step phase of 2π/3 radians inside the square area of 250 × 250 as shown in Fig. 1(a). A step phase object is used here as we are interested in image plane holography where the object wave O may contain sharp edges. Although in a realistic case the edge sharpness will be limited by numerical aperture of the imaging system, in this case we have assumed a phase object with ideal phase step without any band-limit which is generally considered to be a hard problem for traditional Fourier filtering approach. In the methods based on optimization approach the step phase reconstruction is still a hard problem due to involvement of empirical parameters <cit.>. The reference beam is an off-axis plane wave given by R = exp[i 2π (f_1 x + f_2 y) with f_1 = f_2 = 0.04 /pixel. The corresponding interferogram is shown in Fig. 1(b). The interferogram has been simulated with Poisson noise equivalent to an average  light level 10^4 photon/pixel. Since the dc and cross-term peaks of this interferogram are separated in Fourier space, the cross term may be filtered as shown in Fig. 1(c) where a filter of radius 0.5 times the distance between the dc and cross terms peaks in Fourier space has been used for illustration. A Hamming window is also applied to the filtered cross term prior to computing the inverse Fourier transform in order to mitigate ringing artifacts. The resultant phase map is shown in Fig. 1 (d) has poor resolution compared to the step phase object in Fig. 1(a). Recently a Hilbert transform processing methodology <cit.> has been demonstrated which provides superior resolution than what is obtained in Fig. 1(d), however, this procedure still requires a band-limited object wave.


We will now proceed with a description of MGD iteration for recovering the complex object wave O from a single-shot interferogram H. In general the goal of MGD is to find a solution O such that costs associated with both data consistency C_err and the Total Variation C_TV as in Eq. (<ref>) simultaneously achieve minimal numerical values and that additionally the two functional gradients associated with these costs balance each other. For a given guess solution O, we start by defining two unit vectors:

    û_1 =∇_O^* C_err(O,O^*)/|| ∇_O^*C_err(O,O^*)||_2

and 

    û_2 =∇_O^* C_TV(O,O^*)/|| ∇_O^*C_TV(O,O^*)||_2.

Here, the two functional gradients (or Wirtinger derivatives) are defined as:

    ∇_O^* C_err(O,O^*) = -2(H-|O + R|^2)(O+R),

and

    ∇_O^* C_TV(O,O^*) = - ∇.[∇_x O/|∇_x O|x̂+∇_y O/|∇_y O|ŷ].

Next we introduce a vector u that is along the mean direction that bisects û_1 and û_2:

    u = û_1 + û_2/2,

and consider an iteration of the form:

    O^(n+1) = O^(n)  - t ||O^(n)||_2  [ u]_O=O^(n).

Here the parameter t denotes the step size in units of the norm ||O^(n)||_2 of the guess solution after n iterations. Note that since û_1 and û_2 are unit vectors, for any arbitrary value of t, the changes in the solution 

    D_1,n   = ||  t/2 ||O^(n)||_2  [û_1]_O=O^(n)  ||_2, 
    
       D_2,n   = ||  t/2 ||O^(n)||_2  [û_2]_O=O^(n)  ||_2

due to the progression in the error and TV reducing directions are always guaranteed to be equal. The iteration is much easier computationally as compared to the alternating minimization scheme required for ASD-POCS type algorithms. The scheme for selecting t will be explained later. In order to understand the progression of the solution by MGD, we examine the behaviour of C_err, C_TV and the angle θ between the directions of û_1 and û_2 as the iterations progress. Since the two vectors û_1 and û_2 are complex-valued, for the purpose of calculating angle between them, we form new vectors by concatenating their real and imaginary parts:

    v_1 = [real(û_1j) , imag(û_1j)] 
     v_2 = [real(û_2j) , imag(û_2j)].

Here the index j runs over computational window size (j = 1, 2, ..., (500)^2). The angle between u_1 and u_2 is then defined as:

    θ = arccos[ v_1 · v_2 /|| v_1 ||_2 || v_2 ||_2].


As seen in illustrations below, we observe that following this scheme leads to successive solutions where
both C_err and C_TV nominally decrease and eventually, the angle between û_1 and û_2 becomes obtuse implying that the resultant solution balances the two terms C_err and C_TV. We term this scheme as “Mean Gradient Descent” in view of the definition Eq. (<ref>) of vector u and that the two costs are seen to nominally progress to their minimal possible values. 


For a typical illustration with Poisson noise corresponding to the average light level of 10^4 photons/pixel in the hologram data frame, we initiate the MGD iteration with a random complex valued function with amplitude and phase distributed uniformly in [0,1] and [0,2π] respectively. As per Eq. (<ref>), in each iteration the solution is simply pushed in the direction - u = -(û_1 + û_2)/2. In the following discussion we provide our thought process for selecting the step size t which is kept constant for initial iterations and then reduced slowly.

As already explained before, the problem of determining amplitude |O(x,y)| and phase ϕ_O (x,y) from a single hologram data frame H(x,y) does not have a unique solution even when the reference beam R(x,y) is known exactly. There can be multiple combinations of the amplitude and phase functions that may satisfy the hologram data. Starting with any random guess solution, if we reduce C_err alone by a gradient descent scheme, we observe that we reach a solution representing a local minimum of C_err that does not simultaneously have a low value for C_TV. The progression in the mean direction  however leads to moving away from such undesirable solutions. We note that since û_1 and û_2 are unit vectors, the maximum magnitude of u is equal to 1. We therefore initiate t with a nominal trial value of 0.1, suggesting that the solution can maximally change in norm by 10 % in a single iteration. When this initial value of t is held constant for the first few hundred iterations, we find that the solution reaches close to the desirable solution in the solution landscape. The phase and amplitude maps corresponding to the resultant solution after 200 iterations with a fixed t value are shown in Fig. 2 (a), (b) respectively. We note that the solution has the expected features of a step phase object but additionally contains undesirable fringe-like artifacts. The blue curves in Fig. 3(a), (b), (c) show the plots for C_err, C_TV and θ as a function of iteration number respectively. Fig. 3 shows plots of these quantities for three different noise levels as we will discuss later. For now we will concentrate on the blue curves in these plots that correspond to the Poisson noise corresponding to the average light level of 10^4 photons/pixel. At the end of 200 iterations we observe that the three curves for C_err, C_TV and θ have nearly flattened. (The blue, red and green curves in Fig. 3 representing different noise levels are nearly overlapping in this region). This behaviour suggests that the solution has essentially stagnated. While the solution appears to be close to what we want (with some additional artifacts) the fixed value of t is too large at this point and the solution is unable to approach the desired minimum in C_err or C_TV. 

In the further iterations, we check the value of error term C_err. If the numerical value of C_err has increased compared to that in the previous iteration, the step size t is decreased slowly by a constant factor 0.99. The plots in Fig. 3 (a), (b), (c) start showing interesting behaviour at this point - the numerical values of C_err and C_TV are seen to start nominally decreasing while the angle θ starts increasing as we reduce t. The oscillations in these curves right after t starts reducing represents the fact that t is still too large than what is desired and is being slowly adjusted to a lower value. It is important to note that every iteration simply involves a fixed straightforward computation of û_1 and û_2 followed by progressing the solution in the mean direction, thus the computational cost per iteration is minimal. 

In Fig. 4 (a),(c) and Fig. 4 (b),(d) we show the amplitude and phase maps of the resultant solutions after 500 and 2000 MGD iterations for the data with Poisson noise corresponding to the average light level of 10^4 photons/pixel. The computational time for a MATLAB implementation on a 3.5 GHz processor and 16 GB RAM was observed to be approximately 0.07 seconds per iteration. The profile of the phase function after 2000 iterations along the dotted line in Fig. 4 (d) is plotted in Fig. 4 (e) and shows excellent recovery of the sharp edge in phase function. 

In order to understand the sensitivity of MGD approach to noise we considered two additional Poisson noise realizations of the interferogram in Fig. 1(b) corresponding to an average light level of 10^3 and 10^5 photons/pixel. The behaviour of C_err, C_TV and θ for these cases as a function of iterations is also shown in Fig. 3 (a), (b), (c) respectively (red and green curves). We observe that with increasing light level, the resultant numerical values of C_err and C_TV are consistently lower and the numerical value of C_TV increasingly approaches the ground truth value of TV of the simulated step-phase object (shown in magenta). We note that the rise in angle θ after 200 iterations when we start reducing step size t is fastest for the data with the highest relative noise (10^3 photons/pixel). This is expected since the balancing of the two terms of the cost function should start happening at higher value of C_err for data with higher noise.

Table 1 shows the RMS phase error between the ground truth phase map and the reconstructed phase map for the three noise realizations and the RMS error performance is best for the data with lowest relative noise as expected. The RMS error performance for all the three cases is excellent and in fact superior to the expected performance purely based on shot noise considerations<cit.>. Superior performance compared to single-pixel based shot noise level is expected due to the sparsity of the object wave as already demonstrated in <cit.>. A more detailed analysis of performance of MGD with respect to the light level and the sparsity of the object wave will be taken up in future. For random initial guess we observe that the initial direction between û_1 and û_2 was observed to be close to 90^∘ suggesting that the two directions were independent. As the iterations progressed, the angle θ was initially acute but eventually became obtuse and close to 180^∘ as the number of iterations was made very large. This behaviour of angle θ confirms our main motivation for using the MGD approach. 



§ PERFORMANCE OF MGD FOR ON-AXIS AND NEAR ON-AXIS INTERFEROGRAMS WITH SPHERICAL REFERENCE BEAM

From the previous section it is clear that the working of MGD is independent of the various noise levels in the hologram data. It should also be noted that the MGD approach never utilized the fact that we were analyzing an off-axis interferogram. As long as the form of reference beam is known, MGD should be able to handle the hologram data in the same manner as the off-axis case as we illustrate in this section. We now test the evolution of MGD solution for two interferogram recording configurations shown in Fig. 5 (a), (b) where the reference beam is in the form of on-axis and near on-axis spherical beams. The on-axis spherical wave is taken of the form R = exp(i 2 π p (x^2 +y^2)) and the near on-axis spherical wave has form R = exp(i 2 π q [(x-51)^2 +(y-51)^2]) with p = q = 0.0004/pixel^2. The interferograms in Fig. 5(a), (b) are generated with Poisson random noise corresponding to average light level of 10^4 photons/pixel. Note that for both on-axis and the near on-axis spherical reference beam configurations, the dc and the cross terms in the interferograms substantially overlap in the Fourier domain as clearly visible in Fig. 5(c),(d). As a result there is no effective Fourier filtering strategy (as in Fig. 1) that can separate out the object wave even in an approximate sense. The MGD iterations is independent of such considerations and high quality object wave reconstructions are obtained as shown in Fig. 5 (e), (f). The edge profiles corresponding to resultant solutions in  Fig. 5 (e), (f) are plotted in Fig. 5(g),(h) along the dotted lines respectively. The phase profiles clearly show the excellent step phase recovery for both the cases. The behaviour of C_err, C_TV and angle θ for the two configurations as the iterations progress is observed to be similar to that of the previously illustrated off-axis case in Fig. 3(a)-(c) respectively. The RMS phase error, with respect to the ground truth solution, calculated after 2500 MGD iterations for the reconstructed phase solutions in Fig. 5(e), (f) are 0.0040rad and 0.0048rad respectively which is similar to the numerical value in Table 1. The computational time per iteration for both the cases illustrated above is identical to the off-axis case as the steps involved in MGD are independent of the nature of interference pattern. From these results, MGD appears to be a robust methodology that uniformly works with multiple interferometric configurations and noise levels. The algorithm is therefore expected to have applicability for interferometric systems in various geometrical configurations. To the best of our knowledge, the MGD iteration as presented here has not been explored in the prior literature. At present we are reporting results of MGD with illustrative examples for lack of a formal proof for its convergence properties. Such proofs if possible will have to be worked out in future. Our numerical trials in this and the previous section however suggest that MGD works uniformly in a robust manner and provides excellent complex object wave recoveries from single-shot interferogram data. We believe that MGD as a concept may be useful beyond interferometry in a multitude of optimization problems that are similar in nature to the present problem as described in Eq. (<ref>).




§ DISCUSSION AND FUTURE OUTLOOK

In summary we have presented a new optimization approach that we call as Mean Gradient Descent (MGD) for single-shot interferogram analysis. Unlike the usual optimization approaches which aim to minimize a cost function, we aim to reach a solution point where the data consistency and constraint penalty functions balance each other. This is achieved by iteratively progressing the solution in a direction that bisects the descent directions for the data consistency (or error) and the penalty terms. This approach does not require any free parameters. MGD scheme works uniformly for varying noise levels in the data as well as data representing different interferometry configurations. MGD involves straightforward computation per iteration which is very simple to implement compared to alternating minimization schemes. As illustrated in our work MGD showed excellent object wave recoveries when interferograms in different configurations (on axis or off axis) were used. For the step phase object used, MGD showed excellent phase step recovery indicating full pixel resolution performance. Since MGD effectively utilized the expected object wave sparsity for phase reconstruction, the rms phase error performance better than the usual definition of single-pixel based shot-noise level was observed. A detailed study of this aspect will be carried out in future. We believe that a robust approach MGD can lead to widespread employment of optimization based methodologies in interferometry and digital holographic imaging applications. The associated devices can thus operate in single-shot mode with full pixel resolution performance as well as superior accuracy. MGD as a concept can potentially work for a number of optimization problems as we will explore in future. 
1
Hariharan
P. Hariharan, Optical Holography, (Cambridge U. P.,1984).
Takeda
 M. Takeda, H. Ina, and S. Kobayashi, “Fourier-transform method of fringe-pattern analysis for computer based topography and interferometry,” J. Opt. Soc. Am. A 72, 156–160 (1982).
 Kreis
 T. Kreis, Thomas,“Digital holographic interference-phase measurement using the Fourier-transform method,” J. Opt. Soc. Am. A 3, 847-855 (1986).
Creath
K. Creath, “Phase shifting interferometry techniques,” in Progress in Optics, E. Wolf, ed. (Elsevier, 1988), Vol. 26, pp. 357–373
Yamaguchi
I. Yamaguchi and T. Zhang, “ Phase Shifting Digital Holography,” Opt. Lett. 22, 1268-1270 (1997).
Servin
M. Servin, J. L. Marroquin and F. J. Cuevas, “Demodulation of a single interferogram by use of a two-dimensional regularized phase-tracking technique,” Appl. Opt. 36(19), 4540-4548 (1997).
Servin2014
M. Servin, J. A. Quiroga, and J. M. Padilla, Fringe pattern analysis for optical metrology: theory, algorithms, and applications, (John Wiley & Sons, 2014).
Liebling
M. Liebling, T. Blu, and M. Unser, “Complex-wave retrieval from a single off-axis hologram,” J. Opt. Soc. Am. A 21(3), 367-377 (2004).
Brandwood
D. H. Brandwood, “A complex gradient operator and its application in adaptive array theory,” Proc. IEE 130, 11–16 (1983).
Kedar
K. Khare, P.T.S. Ali, and J. Joseph, “Single shot high resolution digital holography,” Opt. Exp. 21,  2581-2591 (2013).
Mandeep2015
M. Singh, K. Khare, A.K. Jha, S. Prabhakar, and R.P. Singh, “Accurate multipixel phase measurement with classical-light interferometry,” Phys. Rev. A 91, 021802 (2015).
Bouman
D. J. Lee, C. A. Bouman, and A. M. Weiner, “Single shot Digital Holography Using Iterative Reconstruction with Alternating updates of Amplitude and Phase,” Electronic imaging 19, 1-6 (2016).
Mandeep2017  
M. Singh and K. Khare, “Single-shot interferogram analysis for accurate reconstruction of step phase objects,” J. Opt. Soc. Am. A 34, 349–355 (2017).
Shi
Y. Shi, and Q. Chang, “Efficient algorithm for isotropic and anisotropic total variation deblurring and denoising,” J. Appl. Math. 2013, (2013).
Pan
E. Y. Sidky. and X. Pan, “Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization,” Phys. Med. Biol. 53(17), 4777 (2008).
Legarda
R. Legarda-Saenz, C. Brito-Loeza, and A. Espinosa-Romero, “Total variation regularization cost function for demodulating phase discontinuities,” Appl. Opt. 53, 2297–2301 (2014).
Galvan
C. Galvan, and M. Rivera, “Second-order robust regularization cost function for detecting and reconstructing phase discontinuities,” Appl. Opt. 45, 353-359 (2006).
Baek
Y. Baek, K. Lee, S. Shin, and Y. Park,“ Kramers-Kronig holographic imaging for high-space-bandwidth product,” Optica 6, 45-51 (2019).
Walls
D. F. Walls and G. J. Milburn, Quantum optics, 2nd ed. (Springer Science & Business Media, 2007).


