ADMM and Spectral Proximity Operators in Hyperspectral Broadband Phase Retrieval for Quantitative Phase Imaging

A novel formulation of the hyperspectral broadband phase retrieval is developed for the scenario where both object and modulation phase masks are spectrally varying. The proposed algorithm is based on a complex domain version of the alternating direction method of multipliers (ADMM) and Spectral Proximity Operators (SPO) derived for Gaussian and Poissonian observations. Computations for these operators are reduced to the solution of sets of cubic (for Gaussian) and quadratic (for Poissonian) algebraic equations. These proximity operators resolve two problems. Firstly, the complex domain spectral components of signals are extracted from the total intensity observations calculated as sums of the signal spectral intensities. In this way, the spectral analysis of the total intensities is achieved. Secondly, the noisy observations are filtered, compromising noisy intensity observations and their predicted counterparts. The ability to resolve the hyperspectral broadband phase retrieval problem and to find the spectrum varying object are essentially defined by the spectral properties of object and image formation operators. The simulation tests demonstrate that the phase retrieval in this formulation can be successfully resolved.


Introduction
Multispectral (MS) and hyperspectral (HS) images differ in amount and width of spectral channels. The number of spectral channels varies from a few for MS imaging and goes from tens up to hundreds for HS imaging. We do not make a difference between these two scenarios and use the word 'hyperspectral' addressing to both MS and HS ones. HS imaging is much more informative as compared with conventional RGB imaging and indispensable in many applications such as remote sensing [1] , biology [2] , medicine [3] , quality material and food characterization [4] , control of ocean [5] and earth pollution [6] , etc. A flow of publications on computational methods developed for various applications of HS imaging is growing fast (e.g., [7][8][9] ).
Conventionally, spectral information is extracted from observations of two different modes: (1) registered for each spectral channel separately, channel-by-channel, or (2) registered as the total power of entire spectral channels with subsequent spectral analysis. In this paper, we follow the latter scenario with a broadband illumination of an object of interest and a broadband sensor registering the total power of the impinging light beam.
single [25] ) intensity measurements (squared linear projections) Y t = | A t U o | 2 , t = 1 , . . . , T , where Y t ∈ R M , and A t ∈ C M×N is an image formation operator. Here and in what follows t stays for experiment number, and T is a maximum number of experiments. The latter is used also as notation for a set of experiments. The heuristic and semi-heuristic iterative algorithms with alternative projections to the object and measurement planes are well known and studied starting from the works of Gerchberg and Saxton [26] and Fienup [27] . These techniques are proven to be efficient for various optical applications.
In recent development, the setup for phase retrieval assumes a modulation of the object U o by random phase masks M t ∈ C N with the observation model where ' •' stands for the Hadamard, element-by-element, product of two vectors.
In general, the phase retrieval is a non-convex inverse imaging problem. Theoretical studies formulating restrictions on object and image formation models leading to uniqueness and convergence of the iterative algorithms are of special interest. The model (1) with the operator A t given as the Fourier transform (FT) is common for many theoretical works [28][29][30][31] . Optimization of random phase coding modulation masks M t is a hot topic in phase retrieval. It is used in order to improve the algorithm's convergence and spatial resolution of phase estimates [14,32,33] . In this paper, random phase masks are essential instruments enabling hyperspectral phase retrieval. The extended review of the recent developments in the mathematical foundation of the phase retrieval problems can be seen in [34][35][36] .
In this paper, we present a novel phase retrieval formulation that provides prospective powerful instrumentation for broadband HS complex domain (phase/amplitude) imaging. In what follows, we use a conventional vectorized representation of 2 D images as vectors. For the object of interest, it is the vector U o,k ∈ C N , N = nm , where n and m are width and height of 2 D image; integer k stays for a number of the spectral component of the hyperspectral object, k ∈ K , where K is a set of the spectral components.
For the several-bands case, straight forward approach is to expand the model (1) for the number of used wavelengths by simply getting observations separately on each wavelength (e.g., [37] ) by the analogy of placing waveband filers after a light source. However, we consider a more complex task where spectral observations are superimposed. Therefore, we introduce the HS broadband phase retrieval as a reconstruction of a complex-valued object U o,k ∈ C N , k ∈ K, from the intensity measurements: For the noisy intensity observations with additive noise ε t , Y t is replaced by Z t : Here, A t,k ∈ C M×N are linear operators modeling a propagation of 2D object images from the object plane to the sensor, and t is the index of experiment. Total intensity measurements Y t ∈ R M + are calculated over the spectrum K as the sum of spectral intensities | A t,k U o,k | 2 . It is essential in this paper that the object U o,k and the operators A t,k are spectrum dependent, varying in k . In optics, it means that the reflective and transmissive properties of the object (specimen) as well as the light propagation operators depend on the wavelength of light. There is a small number of publications relevant to the HS broadband phase retrieval from the total intensity measurements considered in this paper. We briefly review these works. The HS phase retrieval in the paper [38] is formulated as object reconstruction from the intensity observa- where the object U o is spectrally invariant, i.e., does not depends on k . This assumption differs [38] essentially from the setup considered in this paper. With the introduced observation model (2) , the interferometric HS methods also can be interpreted as phase retrieval problems. Indeed, the HS digital holography uses the observation model where R t.k stays for the harmonic reference signal varying from experiment-to-experiment [39][40][41] .
For the shearing HS digital holography, we can introduce the observations as where the reference beam R t,k is not additive to the object as in (4) but propagates through the object U o,k [42][43][44] . The harmonic modulation of the signals in (4) and (5) is a special feature in computational holography. The solution for (4) can be given by FT of the observations combined with complex domain filtering [41] . The solution for (5) is obtained by FT of the observations included in the phase retrieval iterations [44] .
Note that the measurements Y t in (2) with arbitrary A t,k are very different from those in (4) and (5) , where the harmonic modulation is used for spectral analysis of observations typical for interferometry and holography. The observations (2) -(3) do not include any spectral analysis signals or spectral analysis instruments.
The ability to resolve the HS phase retrieval problem is completely defined by the developed algorithm and the image formation operators A t,k what includes a design of encoding phase masks.
The main novelty of this paper is the formulation of the HS phase retrieval for a spectrally varying object from intensity observations (2) -(3) with spectrally varying image formation operators. The developed algorithm is derived from the variational formalization of the problem based on the alternating direction method of multipliers (ADMM) for complex-valued variables and the novel Spectral Proximity Operators obtained for Gaussian and Poissonian observations. Computations for these operators are reduced to the solution of the sets of cubic (for Gaussian observations) and quadratic (for Poissonian observations) algebraic equations. The model of the object U o,k is unknown and is not exploited in the algorithm iterations.
Pragmatically, we target on the HS phase retrieval formulation and the algorithm's development. The numerical and physical tests prove that the HS phase retrieval in this general formulation can be resolved. The proposed setup allows a simple optical implementation, possibly lensless, which is much simpler than implementations used in interferometry and holography imaging.
The rest of the paper is organized as follows. The algorithm development is a topic of Section 2 . From a simple variational formulation of the problem, one fidelity term, and one prior term, we go to ADMM considering the features concerning the formulation for complex-valued variables. The spectral proximity operators (SPO) are introduced as solutions of optimization problems with quadratic penalization. A complex-domain block matching filter is used as a regularization tool for the HS phase retrieval. The results of simulation tests are demonstrated in Section 3 . The setups and results of physical experiments are presented in Section 4 . Conclusions and discussions are given in Section 5 .

Approach
Let l({ Z t } , k ∈ K | U t,k | 2 ) be a neg-log-likelihood of the observed { Z t } , t ∈ T , and the complex-valued images at the sensor plane are Here and in what follows, the curved brackets are used as an indication of a set of variables. The neg-log-likelihood is a fidelity term measuring the misfit between the observations and the prediction of the intensities of U t,k summarized over the spectral interval. Various inverse imaging computational methods have been developed under variational optimization formulation using one fidelity term and one prior term.
We start from this tradition and introduce an unconstrained maximum likelihood optimization to reconstruct the object 3D cube { U o,k } , k ∈ K, from the criterion of the form: where a second summand f reg is an object prior (regularization function). The straightforward optimization in (6) is too complex. The alternating direction method of multipliers (ADMM) provides a valuable alternative to this sort of problem [45][46][47][48] . The following logic leads from (6) to ADMM.
Let us reformulate (6) as a constrained optimization on The problems (6) and (7) are equivalent with an advantage introduced by U t,k which can be treated as a splitting variable such that optimization can be arranged as sequential on U t,k and U o,k of the two loss functions l and f reg . The concept of splitting variables is at the root of many modern optimization methods (e.g., [46,47] ).
One of the popular ideas to resolve (7) is to replace it by the unconstrained formulation with the parameter γ > 0 : (8) Here and in what follows k,t means summation over k ∈ K and t ∈ T , and k summation over k ∈ K. The second summand is the quadratic penalty for the difference between A t,k U o,k and the split- It is recommended in these iterations to take γ varying and going to smaller values γ → 0 . The performance of the algorithm depends on this parameter.
A valuable counterpart to (8) with a decreasing sensitivity to γ is a reformulation of (7) with the Lagrange multipliers t,k and the augmented Lagrange loss function. For the complex-valued variables, this counterpart is of the form [49,50] : Here the Lagrange multipliers are complex-valued, t,k ∈ C M , and the superscript H stays for the Hermitian transpose. The complexvalued variables introduce specific features to this formulation, differing it from the conventional real-valued one.
The alternating direction algorithm of multipliers (ADMM) for the Lagrangian (10) is composed of iterations [50] : The last equation updates the Lagrange multipliers. It can be verified that As || H t,k || 2 2 does not depend on U t,k and U o,k , the iterations (11) can be rewritten for (10) as We exploit the ADMM algorithm in this form for the considered HS phase retrieval problem. The convergence of the complexvalued ADMM algorithm obtained recently in [49] can be treated as a valuable argument in favor of this type of the algorithms, despite the fact that these results are shown for the different phase retrieval problem. In the following sections, we derive the solutions for the optimizations in (12) .

Gaussian observations
The loss function in the first row of (12) is of the form For minimization min U t,k J, we calculate the derivatives ∂ J/∂ U * t,k and consider the necessary minimum conditions ∂ J/∂ U * t,k (r) = 0 .
These calculations lead to a set of complex-valued cubic equations with respect to U t,k (r) : Here, r stays for the spacial coordinates instead of (x, y ) in 2 D image representations. Note that these equations are separated on r as well as with respect to t.
The following manipulations resolve the set (14) . Calculating the squared absolute values for both sides of these equations and producing summation on l ∈ K, we arrive to With the notations: we rewrite (15) in a compact form as a cubic Cardano's equation with respect to x : The Cardano's formulas give the solutions for (17) . The coefficients of this equation are real-valued. It has three roots, which, depending on the discriminant D , are real-valued for D ≤ 0 and one root is real, and two others are complex-valued for D > 0 . We are looking for real-valued nonnegative solutions x ≥ 0 . If there are several such solutions, we select the one which provides the smallest value to the criterion (13) . If such ˆ x is selected, the solution for U t,k (r) according to (14) is as follows: Equation (17) should be solved for each r and t separately. It follows that the criterion for the selection of the best ˆ x t (r) and ˆ U t,k (r) should be pixel-wise: where ˆ U t,k (r) is given by (18) .
Thus, the solution for ˆ U t,k (r) is obtained by the following three stage procedure: (1) Solution of (17) ; (2) Analysis of the roots of the Cardano's equation selecting the one minimizing (19) ; The solution ˆ U t,k (r) is a crucial point of the developed algorithm as it allows extracting spectral information from the total observations summarized over entire spectral components. This solution can be treated as a synchronous detector with a modulating where ˆ x t (r) is an estimate of Z t (r) . The spectral resolution of this detector is restricted by differences in the spectral properties of in the spectral properties of the image formation operator A t,k and in the spectral properties of the object U o,k (r) .

Poissonian observations
The observations take random integer values, which can be interpreted as a counted number of photons detected by the sensor.
This discrete distribution has a single parameter μ and is defined Here p(Z t (r) = n ) is the probability that a random Z t (r) takes value n , where n ≥ 0 is an integer. The parameter μ is the intensity flow of Poissonian random events. The parameter μ is different for different experiment t and r with values μ = Y t (r) . The probabilistic Poissonian observation model is given by the formula where χ > 0 is introduced as a scaling parameter for photon flow.
According to the properties of the Poisson distribution, we have for the mean value and the variance of the observed Z t (r) , In these formulas, χ defines the level of the random component in the signal and can be interpreted as an exposure time. A larger χ means a larger exposure time, and a larger number of the photons with the intensity The necessary minimum condition on U t,k can be calculated as To resolve the problem, we use the methodology applied to the Gaussian data. We calculate the sums over l for the squared magnitude values of the left and right sides of (22) and using the notations (16) obtain the set of the quadratic equations for x : As b 2 − 4 ac = q 2 + 4(1 + γ χ ) γ Z t (r) ≥ 0 , the quadratic equations have two real roots. The estimate for U k,t (r) corresponding to the solution ˆ x according to (22) is as follows As there are two real-valued roots, to select the proper root, which is nonnegative and the best minimizer for J, we use pixelated summands of this criterion with fixed t and r : Both solutions for min U k,t J (Gaussian and Poissonian) can be interpreted as the proximity operators being obtained minimizing the likelihood items (first item in (8) ) regularized by the quadratic penalty (second item in (8) ). With the standard compact notation for the proximity operators [38,51] , we denote these solutions as: where f stays for the neg-log-likelihood part of J and γ > 0 is a parameter.
These proximity operators resolve two problems. Firstly, the complex domain spectral components U t,k (r) are extracted from the observations where all of them are measured jointly as the total power of the signal. Thus, we obtained the spectral analysis of the observed signals. Secondly, the noisy observations are filtered with the power controlled by the parameter γ compromising the noisy intensity observations Z t and the power of the predicted signal A t,k U o,k at the sensor plane. We introduce a term 'Spectral Proximity Operator' (SPO) for (27) . The proximity operators have already been introduced for the phase retrieval in its monochromatic version (single wavelength, k = 1 ) in [52,53] . A generalization of these proximity operators to a sum of intensity measurements is given in [52] , where it was exploited for phase retrieval from sub-sampled data. However, this generalization is not sufficient to deal with the considered setup of the HS phase retrieval.
The proposed Spectral Proximity Operator is new. The cubic and quadratic equations have appeared for Gaussian and Poissonian observations, respectively, but they are scalar equations for the monochromatic phase retrieval, while for the considered HS problem, we need to resolve the sets of cubic and quadratic equations. The spectral resolution properties of these proximity operators define their unique peculiarity.
In our MATLAB implementation of SPO the calculations for all pixels r are done in parallel using analytical representations for solutions of the quadratic and cubic Equations (17) and (24) .

Minimization on U o,k
Minimizing (12) on U o,k we drop the regularization term and replace its effects by joint filtering of { U o,k , k ∈ K} . This approach is in line with a recent tendency to use efficient filters for the regularization of ill-posed solutions without formalization of the regularization penalty in variational setups of imaging problems (e.g., [54] ). As a replacement for the prior, we use the Complex Cube Filter (CCF) [41] . Then, the solution for U o,k is of the form where the regularization parameter β reg > 0 is included if t A H t,k A t,k is singular or ill-conditioned.

HS complex domain phase retrieval (HSPhR) algorithm
Let us present the iterations of the algorithm developed based on the above solutions.
(5) Backward propagation and preliminary object estimation: The complex domain initialization (Step 1) is required for the considered spectral domain k ∈ K. In our experiments, we assume 2D random white-noise Gaussian distribution for phase and a uniform 2D positive distribution on (0,1] for amplitude, which are independent for each k . The algorithm robustness to the randomness of initialization was tested on 100 different initializations which resulted in successful reconstructions with a standard deviation of only 0.001 in relative error value. The Lagrange multipliers are initialized by zero values, k,t = 0 . The forward propagation is produced for all k ∈ K and t ∈ T (Step 2). The update of the wavefront at the sensor plane (Step 3) is produced by the proximal operators. For the Gaussian observations, this operator is defined by (17) - (18) and for the Poissonian observations by (24) - (25) . It requires solving the polynomial equations, cubic or quadratic for the Gaussian and Poissonian cases, respectively.
The initial iterations are produced with zeroed Lagrange variables. In Step 4, the Lagrange variables are updated starting from s > s 1 , where s 1 is given. This delay in involving the Lagrange variables can be used to shorten the calculation time and in order to improve the convergence of the algorithm. A similar delay with similar reasons is used in the application of the CCF filter in Step 6.
The backward propagation of the wavefront from the sensor plane to the object plane is combined with an update of the spectral object estimate in Step 5. The sparsity-based regularization by Complex Cube Filter (CCF) is relaxed by the weight-parameter 0 < β s < 1 at Step 6. The iteration number n is fixed (Step 6) in this algorithm implementation. However, the stop condition might be estimated by the stagnation of the reconstruction by calculating relative error for consecutive reconstructions of U s o,k and U s −1 o,k and stop the iterations when it becomes lower than a threshold.
The CCF algorithm is designed to deal with 3D complex-valued cube data and is introduced in detail in [41] . Examples of its application and features of this algorithm can be seen in [44,55,56] .
A few notes on this algorithm. The CCF algorithm is based on SVD analysis of the HS complex-valued cube. It identifies an optimal subspace for the HS image representation, including both the dimension of the eigenspace and eigenimages in this space. The Complex-Domain Block-Matching 3D (CDBM3D) algorithm [57,58] filters these eigenimages. Going from the eigenimage space back to the original image space, we obtain the reconstruction of the object cube. The CDBM3D is developed as an extension to the complex domain of the popular BM3D algorithm [59] . CDBM3D is implemented in MATLAB as Complex Domain Image Denoising (CDID) Toolbox [58] . A sparsity modeling in CDBM3D is based on patch-wise 3D/4D Block-Matching grouping and 3D/4D High-Order Singular Decomposition (HOSVD) is exploited for block-wise spectrum design, analysis, and filtering (thresholding and Wiener filtering). Three types of sparsity are developed in CDID, corresponding to the following representations of complex-valued variables: I. Complex variables (3D grouping); II. Real and imaginary parts of complex variables (4D grouping); III. Phase and amplitude (4D grouping). In the forthcoming simulation tests, we use the joint real/imaginary sparsity (Type II).

Numerical experiments
We present the results of numerical experiments in order to evaluate the empirical performance of the developed algorithms. Phase and amplitude imaging by the lensless computational microscopy system is demonstrated. The optical setup illustrating the image formation is shown in Fig. 1 . A broadband coherent laser beam, wavelength range = [400 : 700] nm, of the uniform intensity distribution impinges on a transparent thin object O . The Fig. 1. Schematic optical setup corresponding to our tests and data formation model: 'Laser' is a broadband coherent light source; the object 'O' and the phase mask (implemented by SLM) are flat located in the same plane; 'Sensor' is a registration camera, and d is a distance from the object/mask plane to the sensor. phase coding modulation mask M t is attached to the object. The phase encoding by the mask is implemented by a spatial light modulator (SLM), which is a pixel-wise programmable device allowing to change a mask phase-delay M t (phase coding) from experiment to experiment. The laser beam propagates through the object, goes through the phase mask M t , and freely propagates in the air to the sensor. After that, the intensity of the beam is registered by the sensor as a coded broadband diffraction pattern. The free space propagation linking U o,k and U t,k is modeled by the angular spectrum operator calculated using the forward F{·} and inverse F −1 {·} Fourier transforms: where the angular spectrum transfer function in the 2 D Fourier domain ( f x , f y ) is defined according to [60] as: The λ k in Eq. (29) are the wavelengths corresponding to the spectral components from K, and d is a distance from the object/mask plane to the sensor plane (see Fig. 1 ), d = 2 mm in our tests. The Eqs. (29) and (30) are written for 2 D variables, thus the coding phase mask M t,k and the object U o,k are 2 D complex-valued functions.
Comparing (29) with (2) , the linear operator A t,k as function of k is defined by the spectrum varying AS λ k and the spectrum varying phase-mask M t,k . The variations of the operator A t,k from experiment-to-experiment (index t) is governed by the masks M t,k .
The developed HS phase retrieval algorithm can be treated as a nonlinear comb-filter. The independent random phase masks M t,k are crucial instruments to enable both spatial and spectral resolution of this comb-filter. A larger number K of the spectral channels requires a larger number of experiments T and corresponding T different phase masks (see tests in Section 3.1 ). It can be shown that the matrices t A H t,k A t,k in Step 5 of HSPhR can be well conditioned with a dominated main diagonal provided large T and properly selected randomness for M t . In this case, a high spatial and spectral resolution is guaranteed for HSPhR.
The HS wavefronts and measurements are modeled by propagating the broadband laser beam through the object and the phase modulation masks. Both of them are modeled by the complexvalued transfer function [60] : where h is the thickness of the object/mask, n λ is a refractive index depending on λ.
The phase delay of the wavefront propagating through g is defined by the argument of the exponential function in (31) . It is assumed in our tests that the amplitude a and the thickness h are spatially varying on (x, y ) but invariant with respect to the spectral (wavelength) arguments k and λ. The model (31) shows that the properties of the object and the masks are spectrally varying due to λ in the argument and to dependence on λ of the refractive index n λ .
The set of the wavenumbers λ k , k ∈ K, used in the algorithm defines a spectral sampling of the object spectrum and the wavelengths for which HS phase retrieval is produced. The spatial sampling corresponds to the camera pixel size equal to 3 . 45 μm.
We assume that the amplitude and phase of the object are 64 × 64 images: ' peppers ' for amplitude and ' cameraman ' for phase. We take these images quite different to test how far phase and amplitude can be separated by the proposed algorithm. The phase is scaled in such a way that the object phase delay for all λ ∈ would be in the range [0, π ] . The amplitude a is scaled to the interval (0,1]. For the coding mask, a = 1 , and h is pixel-wise random with equal probabilities taking one of the following five val- For each of the T experiments, the masks are generated independently. Thus, overall we have T different masks. The forward propagation (image formation) operator A t,k in (2) is defined by the propagation model (29) and the modulation masks. It is clear from the model (29) that the propagation model is wavelength varying. The same is true for the modulation mask, as it is fixed for each experiment, but its spectral properties vary with λ according to (31) . As the phase mask is included in the operator A t,k , the complex-valued image at the sensor plane is defined as U t,k = A t,k U o,k , where the object is also spectrally varying according to (31) . The measured intensities are Y t = k ∈ K | U t,k | 2 .
In our simulation, we assume that the refractive index in (31) as a function of λ is known and calculated according to Cauchy's equation [61] with parameters taken for the glass BK7 [62] . The input laser beam is uniform in both phase and amplitude, the amplitude is equal to 1, and the phase is equal to 0. We formulate the phase retrieval as the reconstruction of U o,k .
The transfer function (31) is used for the calculation of the observations Y t and the phase delay in the masks but is not used in the algorithm iterations.
We evaluate the accuracy of the complex-valued reconstruction by the relative error criterion introduced in [28] : where x and ˆ x are the true signal and its estimate respectively. The noise level in observations is characterized by the signal-tonoise ratio (SNR) in dB. For the Gaussian noise case, SNR is defined as where Y t is noiseless and Z t is noisy observations. For the Poissonian obserations: SNR Poisson = 10 log 10 where χ following (20)   The efficiency of the developed algorithm is demonstrated in simulation tests with Gaussian and Poissonian observations. In these experiments, we present the results achieved by the developed algorithms using both the Lagrange multipliers (Step 4) and the CCF filtering (Step 6) as it is in the HSPhR algorithm, Subsection 2.4 . In order to evaluate the value of these two key components of the algorithm, we also show results obtained when these components are disabled.

Accuracy as a function of K and N
A dependence of the algorithm accuracy on the number of wavelengths K and the number of observations T is of special interest. We calculate the accuracy criterion ERROR rel for K = [2 , 4 , 8 , 12] and T = [2 , 6 , 12 , 24 , 36] . The wavelengths for the varying K (spectral channels) are defined as uniformly covering the interval [400 : 700] nm. The number of iterations is fixed to n = 300 .
The relative errors ERROR rel obtained in these experiments are shown in Table 1 .
A larger value of T for a fixed K leads to a more accurate reconstruction with smaller ERROR rel . We found that visually a good phase imaging is achieved provided ERROR rel < 0 . 1 . This threshold for K = 2 is overcome for T = 6 . For K = 4 , 8 , 12 , it happens for T = 6 , 12 , 24 , respectively. We may conclude that T ≥ 3 K is sufficient for good accuracy HS phase imaging. The results in Table 1 are given for nearly noiseless data with SNR = 54 dB, but the conclusion that the inequality T ≥ 3 K is sufficient for good accuracy phase imaging holds for noisy data also. This conclusion is one of the reasons to exploit in the forthcoming tests T = 18 for K = 6 .

Gaussian observations
The relative error maps for the HSPhR algorithm are shown in Fig. 2 for SNR taking values [14 : 55] dB. The left image Fig. 2 (a) demonstrates the accuracy as a function of SNR and λ, while the right image Fig. 2 (b) provides the accuracy as a function of SNR and a number of iterations. In this latter image, the relative errors are averaged over λ, K = 6 . We may conclude from the left image that the acceptable quality imaging, ERROR rel < 0 . 1 , is achieved for all wavelengths provided that SNR is larger than 28 dB. The accuracy for smaller wavelengths is higher than the accuracy for larger values of wavelengths.
The convergence rate can be evaluated from the right image. After 200 iterations and for SNR > 28 dB the accuracy becomes acceptable, ERROR rel < 0 . 1 . The bend in the accuracy map in Fig. 2 (b) well seen at iteration 50 is of special interest. The algorithm's iterations on CCF and the Lagrange variables are disabled in this experiment up to this 50th iteration. Thus, we can observe the dramatic improvement in the convergence rate due to CCF and the Lagrange variables.
In Fig. 3 , we show the relative error maps provided that both the Lagrange variables and CCF are disabled in the HSPhR algorithm completely. The degradation of the algorithm in values of ERROR rel is clear in this case. In particular, it is demonstrated by comparing the error maps in Fig. 2 and Fig. 3 . The relative errors in Fig. 3 are always larger than 0.1, i.e., the imaging accuracy is unacceptable.
Both the Lagrange variables and CCF can be activated from the first iterations, which may improve the convergence rate but at the price of more expensive calculations. The accuracy of the ADMM algorithm as a function of s 1 is demonstrated in Fig. 4 . This curve shows that the best accuracy for 100 iterations is achieved for s 1 = 50 . This result supports our choice for s 1 = 50 in our experiments.

Poissonian observations
For the scenario identical to the considered for the case of Gaussian noise, we introduce results obtained for Poissonian observations. The level of the noise is controlled by the parameter χ, which takes values corresponding to SNR in the interval [14 : 44] dB. The maps of the relative errors are shown in Fig 5 . The left image Fig. 5 (a) demonstrates the accuracy as a function Fig. 3. Relative error maps for the HSPhR algorithm with the disabled Lagrange variables and the CCF filter, Gaussian noise: a) ERROR rel as a function of SNR and wavelength; b) ERROR rel as a function of SNR and iteration number. In b) ERROR rel is averaged over the wavelengths. The achieved accuracy is much worse than that for the HSPhR algorithm with enabled iterations on the Lagrange variables and the CCF filter shown in Fig. 2 . of SNR and λ, while the right image Fig. 5 (b) provides the accuracy as a function of SNR and a number of iterations. In this latter image, the relative errors are averaged over λ, K = 6 . We may conclude from the left image that the acceptable quality imaging, ERROR rel < 0 . 1 , is achieved for all wavelengths provided that SNR is larger than 32 dB. The convergence rate is well seen from the right image. After 150 iterations and for SNR > 32 dB the accuracy becomes acceptable. A visual bend in the accuracy map in Fig. 5 (b) happened at the 50th iteration demonstrates the effects of CCF and the Lagrange variables, which are disabled up to the 50th iteration. We can observe the dramatic improvement in the convergence rate due to CCF and the Lagrange variables.
In Fig. 6 we demonstrate the results obtained by the algorithm where the Lagrange iterations and CCF filtering are disabled completely. Similar to what was discussed for the Gaussian case, we may note a serious degradation in the algorithm performance with ERROR rel > 0 . 2 , which confirms the essential role of both the Lagrange iterations and the CCF filtering for the algorithm performance.

Imaging
In this subsection, we evaluate the visual quality of reconstructions. The algorithm provides the HS broadband imaging, i.e., reconstruction of 3D complex-valued cubes. In Fig. 7 , we show 2D images of amplitude and phase for the middle wavelength of the interval , λ = 640 nm. It is a nearly noiseless case as SNR = 54 dB.      nian data. In all cases, the images of the amplitude and phase are of high visual quality. Relative errors are low, equal to 0.019 and 0.018 for Gaussian and Poissonian data, respectively.
The images are shown in square frames in order to emphasize that the size and location of the object support are assumed to be unknown and reconstructed automatically by the HSPhR algorithm. The true image's support is used only for computation of observations produced for the zero-padded object and is not exploited in the algorithm's iterations. For the amplitude images, the estimates in these frames have nearly ideal zero values, while the phase estimates take random values in the frame areas as the phase cannot be defined for the amplitude equal to zero. Practically, these variations of the phase do not influence the calculation of ERROR rel as the amplitude is estimated quite accurately in these areas and close to zero. Nevertheless, ERROR rel shown in the images are calculated for the central parts of images corresponding to the true location of the image support. Fig. 8 shows the results for the noisy cases of SNR = 34 dB. With this level of SNR, the noise effects are essential. The relative errors are much higher with values 0.058 and 0.08 for Gaussian and Poissonian data, respectively.
In Fig. 9 we show the images obtained by the HSPhR algorithm for the same noisy scenario, SNR = 34 dB, with the Lagrange iterations and CCF filtering disabled. The obtained amplitude/phase images are very noisy, and the visual quality of imaging is low, as compared with the results in Fig. 8 . The relative errors take higher values equal to 0.39 and 0.36 for Gaussian and Poissonian observations, respectively. These experiments provide a visual and numer-ical confirmation of the discussed above role of the Lagrange iterations and CCF filtering as the key components of the HSPhR algorithm. The complexity of the algorithm is illustrated by the computational time per iteration required for the test images shown, K = 6 , T = 18 . This time is equal to 0.8 sec. for calculations without CCF and equal to 21.5 sec. for calculations with CCF. All calculations are done in MATLAB R2019b on a computer with 32 Gb RAM and CPU with a 3.40 GHz Intel® Core TM i7-3770 processor. the MATLAB demo code of the developed algorithm is publicly available Supplementary Code 1.

Optical setup of hyperspectral QPI
In this section, we demonstrate the performance of the developed algorithms in a real-life scenario. The phase object and modulation by the masks M t are realized by a spatial light modulator (SLM). The optical setup implemented in the physical experiments is shown in Fig. 10 . The SLM is a GAEA-2 Holoeye with a resolution of 4160 × 2464 pixels and a pixel size of 3 . 74 μm. The super-continuum laser, = 550 − 2400 nm (YSL photonics CS-5), is an illumination source. To work in the spectral range of the sensor, we limit the laser's spectral bandwidth to = 550 − 10 0 0 nm by a bandpass filter. The camera is monochrome Blackfly S, model BFS-U3-200S6M, FLIR, with a pixel size of 2 . 4 μm.
In Fig. 10 , the illumination wavefront expanded by lenses ' L 1 ' and ' L 2 ' propagates to SLM through the beamsplitter (BS), where SLM changes the wavefront phase distribution according to the object and modulation mask phases. After, by the achromatic doublet lenses ' L 3 ' and ' L 4 ' (with diameter 12.7 mm and focal length 50 mm), the wavefront is transferred to the 'Object' plane. Further, it propagates freely to the registration camera 'CMOS. ' The free propagation distance is d = 2 mm. The normalized laser's spectrum multiplied by the camera's sensor spectral response is shown as an insert in Fig. 10 .

Experiment organization
The object to be imaged is a phase object, i.e. the amplitude is invariant, and the phase is spatially varying. For the distribution of the object phase variation, we selected the 'cameraman ' test image. The phase variation is scaled to the interval [0 , π ] for λ = 520 nm, corresponding to our simulation tests. In a similar way, encoding phase mask M t is random with pixel-wise phase values taking with equal probabilities the following five values [0 , 1 , −1 , 1 / 2 , −1 / 2] · π / 4 , for λ = 520 nm. In each measurement t, the phase shift implemented by SLM is a sum of the object phase and the encoding mask phase. Thus, both the object and the encoding masks are implemented by SLM.
In the shown results, the number of measurements T = 243 (intensity observations Z t ) with a reconstruction of K = 81 spectral channels. It is the ratio T /K = 3 , i.e., the number of spectral channels is three times less than the number of observations. The physical distance d between the 'object' plane and camera was fixed to 2 mm as it is in the simulation tests. In order to improve the focusing sharpness of the back-propagation in AS calculations, we tested various values of d. It was found that the best results are achieved with spectrally varying distances of d calculated as d k = 2 . 3 − 0 . 5 80 (k − 1) mm, k = 1 : K. This modification of the algorithm compensates the supercontinuum laser's phase curvature for separate wavelengths.
To take into consideration the irregularity of the wavefront from the laser, following to [65] , as the preliminary tests, we perform experiments without the object and use these observations to reconstruct the corresponding HS wavefront amplitude and phase distributions at the object plane. Let us denote these complexvalued distributions calculated for all spectral bands as U k (n +1) . We use these distributions for corrections of the results obtained for the observations with the object. We subtract these phase distributions U k (n +1) from the results U (n +1) o,k obtained for the measurements with the object. The following HS phase retrieval results are obtained using these corrections.

Noise in intensity observations
Special tests are produced in order to evaluate noise characteristics in measured intensity observations. The results are obtained from 30 0 0 experiments performed with a fixed random encoding phase mask programmed on SLM. Fig. 11 (a) shows the image of the pixel-wise mean value of 30 0 0 intensity images. The histogram of these mean values is given in Fig. 11 (c). The pixel-wise image of standard deviation (STD) for the intensity observations is shown in Fig. 11 (b). The histogram of this STD image is given in Fig. 11 (d). We may conclude from the latter histogram, that STD is spatially nearly invariant with its distribution highly concentrated at the value equal to 3 a.u. It follows that the additive Gaussian model for the noise in intensity observation with invariant STD is quite relevant model to our experiments. Based on this analysis, we use the Gaussian version of the algorithm HSPhR for data processing in the physical tests.

Reconstruction results
The imaging results for object phase and amplitude are shown in Fig. 12 . The top row corresponds to the true amplitude and phase images as they are programmed on SLM, and the bottom row is for reconstructed images. Fig. 12 (a,b,c) are obtained for λ = 60 0 , 745 , 90 0 nm, respectively. See all spectral reconstructions in Supplementary Video 1. For each wavelength, the reconstructed phase images correspond to the true phases with great detail. We would like to mention here that in the HSPhR algorithm, we do not use the support constraint, but we can note that the support is reconstructed automatically with high accuracy. The reconstructed amplitudes are not spatially invariant as it is in the true images. These amplitude variations are defined by leakage of phase information to amplitude, which can be attributed mainly to non-ideality of SLM performance as well as to the misalignment effects in the optical setup.
Despite the point that the quality of imaging is much lower in physical experiments than it is in simulation tests, the amazing result is that spectral analysis is produced for 81 spectral channels,   and it is done using the numerical algorithm without traditional optical spectrometry. Remind that the ability to resolve the spectral channel is based on the spectral properties of the algorithm.
The blue curve in Fig. 13 shows a normalized spectral distribution of the total image intensity (intensities summarized in the spatial domain) as a function of wavelength. This intensity is calculated for the images reconstructed by the HSPhR algorithm. For comparison, we show in this figure the normalized spectral intensity of the laser beam (red curve) multiplied by the camera sensor response function, as it is in Fig. 10 . This curve is obtained by intensity measurements using a commercial spectrometer. The red and blue curves are close to each other. It shows that the devel-oped algorithm properly preserves the spectral intensity distribution of the illumination source.

Conclusion
A novel formulation of the HS broadband phase retrieval problem is proposed, where the object and image formation operators vary spatially and spectrally. The proposed algorithm is based on the complex domain version of the ADMM technique. The derived spectrum proximity operators and CCF noise suppression are important elements of this algorithm. The spectral proximity operators are defined for Gaussian and Poissonian observations and calculated by solving the sets of cubic (for Gaussian) and quadratic (for Poissonian) algebraic equations. The ability to resolve the HS phase retrieval problem and to find the spectral varying object components U o,k , k ∈ K, thoroughly depends on spectral properties of the spectral channels A t,k . The object model is not used in the algorithm, and the modulation phase masks are calculated in advance and fixed in the algorithm iterations. The simulation and experimental tests demonstrate that the HS phase retrieval in this formulation can be resolved with a quite general modeling of object and image formation operators. The MATLAB demo codes of the developed algorithm are made publicly available, Supplementary Code 1.

Declaration of Competing Interest
Authors declare no competing interests.

Data availability
Data will be made available on request.

Supplementary material
Supplementary material associated with this article can be found, in the online version, at 10.1016/j.sigpro.2023.109095