Lensless hyperspectral imaging by Fourier transform spectroscopy for broadband visible light: phase retrieval technique

A novel phase retrieval algorithm for broadband hyperspectral phase imaging from noisy intensity observations is proposed. It utilizes advantages of the Fourier Transform spectroscopy in the self-referencing optical setup and provides, additionally beyond spectral intensity distribution, reconstruction of the investigated object's phase. The noise amplification Fellgett's disadvantage is relaxed by the application of sparse wavefront noise filtering embedded in the proposed algorithm. The algorithm reliability is proved by simulation tests and results of physical experiments on transparent objects which demonstrate precise phase imaging and object depth (profile) reconstructions.


Introduction
Hyperspectral imaging (HSI) is a technique concerned with capturing 3D data where the first two dimensions are spatial coordinates defining 2D images (slices of 3D data cube) and the third coordinate is a spectral variable. This third dimension adds value to traditional imaging with ability to recover from spectra extra information about an investigated object, as for example in paintings investigations [1], disease detection (e.g. Alzheimer's [2], Parkinson's [3]), or for food ripeness evaluation [4].
In straightforward approaches for 3D data spectral capturing, detectors with a number of color filters are used, as for example in [5]. The hyperspectral digital holography (HSDH) is an alternative approach with a single detector and no color filters, appeared in the 90's years of the previous century starting from the work [6], where HSDH has been developed employing principles of incoherent holography and Fourier spectroscopy. The algorithm proposed by the arXiv:2002.10886v2 [eess.IV] 16 Mar 2020 authors was implemented in a self-reference optical scheme and performed only amplitude imaging with the phase information lost.
As long as HSI has been used mainly in geoscience and remote sensing [7], phase information was not necessary for investigators. Nowadays, HSI spreads in the areas of biological and medical applications [8,9], where phase becomes crucial and almost vital as it brings additional information about the thickness and refractive index of cells and their behavior (e.g. [10,11]) without dyeing and additional sample preparations. To perform phase imaging in HSI, the novel HSDH techniques were developed with various reference beams [12] .
It is possible to retrieve phase in the self-reference setup without any reference beam by modulation one of the object's beams and filtering, which keeps only zero-order diffraction and suppresses all other orders as it is done in [13]. This solution requires a quite complex system with special lenses and other optical components as well as a projection of imaging on a sensor plane.
In this paper, we propose and study a different approach to phase imaging in the HSDH. The system is lensless self-reference and much simpler in implementation as compared with [13], without lenses and projection on the sensor plane. To resolve the phase imaging for the object, the original iterative HS phase retrieval algorithm was developed for HS quantitative phase imaging from noisy intensity measurements. In this way, the complexity of the problem was shifted from the optical hardware components to the software of the algorithm. The diversity needed for HS phase retrieval is enabled by the joint processing of spectral wavefronts covering a broadband range of wavelengths.
The developed algorithm and its verification by simulation and physical experiments are the main contributions of this paper. The advantage of the proposed solution is the simplicity of the lensless optical setup which provides large field-of-view and makes results free from color aberrations which could be extremely corrupting for HS illumination.

Problem description
The self-reference optical scheme assumes that the basic light beam and its phase-shifted copy go through or reflect from the object simultaneously. The interference of these two beams is registered at the sensor plane. The phase-shift (phase-delay) is varying and known. It is usually performed by an interferometric scheme with a moving mirror in one of the two arms, with N steps of size ∆z covering the whole phase-delay distance Z = N · ∆z.
Let V (λ) be the complex-valued spectrum of a broadband light beam after free propagation through the object to the sensor plane, λ ∈ Λ, where Λ is a spectral range of illumination. Then the recorded intensity for the mirror-shift distance z in the delay line is calculated as integral over Λ: The Fourier transform calculated for J(z) over z allows to get the intensity spectrum for V (λ): where FT stays for the Fourier Transform.
As a result, we obtain a set of squared spectra |V (λ)| 2 for a set of wavelengths Λ at the sensor plane but the phase of V (λ) information is lost.
The intensities J(z) calculated over the spectral interval Λ are registered by a multi-pixel detector and the same spectral calculations are valid for each pixel of the detector.

Algorithm
We present the phase retrieval algorithm in the context of thickness (profile) measurement for a transparent object with HS illumination in the self-reference configuration of the optical system. The coherent laser beam of the wavelength λ propagates through the object. It results in the phase delay, ϕ λ (x, y), between the input and output beams. The corresponding output wavefronts A λ (x, y) = |A λ (x, y)| · exp(jϕ λ (x, y)) are different for different wavelengths, λ ∈ Λ, where Λ is a set of the HS wavelengths. These wavefronts at the sensor plane are calculated as V (λ) = P λ {A(λ)}, where P λ stays for the forward propagation operators depending on wavelength.
For the finite discrete set of the shifts ∆z we can calculate the spectral distribution |V (λ)| 2 only for the corresponding finite set of the wavelength, i.e. for Λ = [λ 1 , ..., λ L ]. Then, integration in (1)-(2) is replaced by summation and the Input: {|Ṽ (λs)|}, λs ∈ Λ Initialization: 3. Phase update: Sensor wavefront noise suppression and amplitude update: integral Fourier transform by the discrete Fourier transform. Therefore, according to (1) the measurements by the sensor are defined as This equation gives the intensity calculated over the spectrum range and registered by each pixel of the sensor.
The developed Hyperspectral Phase Retrieval (HSPR) algorithm is composed of two separate stages. The first stage is the spectrum analysis defining the intensity spectral distributions |V (λ)| 2 from the registered set of J(z). These calculations are based on Discrete and Fast Fourier Transform (FFT), for the details of this part of the algorithm we refer to [14].
The second stage is phase retrieval, i.e. a reconstruction of the object complex transfer functions A λ (x, y) from the already given spectra |V (λ)| 2 , λ ∈ Λ. A flowchart of the proposed phase retrieval algorithm is presented in Fig. 1. The algorithm is iterative using typical forward/backward propagation between the object and sensor planes. The following points define the originality of this algorithm: (1) At step 2, the complex domain CDBM3D filter is used for denoising of object transfer functions for each wavelength and in each iteration; (2) At step 5, special filtering is produced for amplitudes of the wavefronts at the sensor plane.
In the case of HS illumination for a transparent object, it is common to rely on the assumption that object wavefronts produced by neighboring wavelengths are similar [15]. It is quite a natural assumption since neighboring wavelengths are close to each other and phase properties of the object are smooth functions of wavelengths.We utilize this assumption and (5) to provide a connection between neighboring wavelength phases by introducing the coefficient such that ϕ λs = ϕ λs−1 · µ s . Here n λs is the refractive index of the object.
In the flowchart of the algorithm in Fig. 1 as input, we have spectral amplitudes |Ṽ (λ s )| at the sensor plane for the whole range of the HS wavelengths reconstructed at the first stage of the algorithm. For the initialization of iterations, we create the first complex-valued wavefront V 1 (λ 1 ) using for the amplitude the first component of the HS |Ṽ (λ 1 )| with zeros guess for the initial phase. After the initialization, we start the two consecutive "for" loops with running variables t and s for cube and wavelength iterations, respectively. For smooth cube iterations t, the loop on s should run consequently through all wavelengths with start and stop at the first wavelength.
Step 1 is the backward propagation from the sensor plane to the object plane, as a result, we obtain an estimate for the object wavefront A t (λ s ). P λ is a wavelength-dependent propagation operator which is defined by the angular spectrum propagation model [16] and with the superscript '−1' it indicates the backward propagation. At step 2, we perform complex-domain sparse noise suppression in A t (λ s ) by the CDBM3D filter [17], which processes amplitude and phase jointly and additionally to sparsity utilizes the correlation between amplitude and phase.
At step 3, we make phase recalculation from s to s + 1 wavelength according to (4). This spectral update is based on the following speculations. The phase delay of the object is varying according to the object thickness. A link between the phase delay and object thickness is defined by the equation: Here, h o (x, y) is a thickness of the transparent object at the point (x, y).
At step 4, we perform forward propagation of the object wavefront to the sensor plane, where at step 5 keeping the phase of V t (λ s+1 ) unchanged we make amplitude update by filtering the amplitude |V t (λ s+1 )|. This noise filtering is performed by the Sensor Noise Suppression (SNS) algorithm. This algorithm is derived based on the optimal solutions for maximum likelihood criteria for Gaussian and Poissonian noise distributions. Both solutions are published in [18].
In this paper, we use the solution for the Gaussian noise. The obtained optimal estimate of the amplitude locates between observations |Ṽ s (λ)| and the amplitude of the object wavefront propagated to the sensor plane |V s t (λ)|.
Next, in decision block, the stopping rule is applied, defined by the maximum number of iterations on t and differences between the phase estimates in successive iterations.
The phase retrieval part of the developed HSPR algorithm shown in Fig. 1 has a structure similar to the algorithms proposed earlier for multispectral phase retrieval (MPR) in the papers [19,20]. Let us discuss briefly differences between HSPR and those MPR algorithms. First of all, different optical setups considered in this and the mentioned papers. The observations in [19,20] are obtained in experiments separate for each laser bands, when in this paper all wavelengths go through the object simultaneously. Second, while the algorithms MPR and HSPR look similar due to the standard iterative structure with forward and backward wavefront propagation as it is originated from Gerchberg and Saxton [21], the efficient noise suppression in steps (2) and (5) define the originality of the developed algorithm HSPR. These filtering steps enhance the HSPR algorithm and help to obtain cleaner reconstructions with faster convergence since noise in the iterative algorithm is predisposed for corruption and ruin the reconstruction. The noise suppression becomes extremely significant in the case of HSPR because in the proposed scheme all spectrum components are observed simultaneously and, therefore, in presence of noise in one region of the spectrum it will be spread throughout the whole spectrum (Fellgett's disadvantage [22]).

Simulations
For a demonstration of the algorithm performance, we provide modeling of the HS data obtained by the system described in Section 2. We use the broadband light spectrum corresponding to a broadband supercontinuum laser source (see spectrum in Fig. 6(b)) and modeled hyperspectral data registration as interference observations according to (3) providing N = 2000 shifts in the delay line which correspond to total Z = 200 µm with ∆z = 100 nm. Since the laser spectrum is not uniform, in regions of the low laser intensity, signal-to-noise ratio (SNR) is low and these regions cannot be used for phase reconstruction. For our tests, we exploit the wavelengths from the high laser intensity region [680 : 820] nm.
As objects under investigation, we model the transparent phase objects assuming that objects' transfer function A λ (x, y) = |A λ (x, y)| · exp(jϕ λ (x, y)) with the amplitude |A λ (x, y)| = 1 and the phase images ϕ λ (x, y) given as the USAF test-target and Cameraman test-image. The distributions of the object phase delays for each λ are calculated according to (5) with maximum h o = 317 nm. Figure 2 demonstrates the USAF (a-b) and Cameraman (c-d) phase objects (depth distributions) as 2D and 3D images. The USAF test-target provides a binary object model, while the Cameraman test-image has a more complex structure with a multilevel piece-wise continuous depth map.
To check algorithm robustness, we include noise in our modeling. The noise ε is additive zero-mean i.i.d. Gaussian with the standard deviation σ noise : J noisy (z) = J(z) + ε(σ noise ). The accuracy of phase reconstruction is defined by the Relative Root-Mean-Square-Error (RRMSE) criterion: whereφ est is the reconstructed phase and ϕ true is the noiseless true phase, ||.|| stays for the Frobenius norm. RRMSE values less than 0.1 correspond to a good quality phase reconstruction.
For characterization of the noise level at the sensor plane we introduce the Peak Signal-to-Noise Ratio (PSNR) defined as PSNR = 10 log 10 max x,y,z (J(z)) σ noise , where max x,y,z (J(z)) is the maximum of the registered intensity of the beam at the sensor plane calculated by maximization on z and (x, y), where (x, y) are the spatial coordinate of 2D images; σ noise is the standard deviation of the additive noise. Descending RRMSEs for growing t indicate consecutive improvement of the reconstruction from iteration-to-iteration. It can be concluded from RRMSE surfaces a) and c) that reliable objects' reconstruction might be obtained only for low noise level with high PSNR values (> 18 dB) and iteration number t > 15. It is caused by Felgett's disadvantage with noise leakage from one spectral component to others. Surfaces b) and d) are provided for σ noise = 0.5, the RRMSE values in the region of λ = 700 nm reflect that the lower illumination intensity results in worse reconstruction quality. In general, the appearance of the RRMSE surfaces is similar for both objects under investigation, however, RRMSE Figure 4: a) PSNR dependencies from noise standard deviation σ noise for registered intensities J(z), black squares curve, and for estimated spectra |V λ S | 2 , red circles. Latter is presented as a mean value of PSNRs for the whole spectra Λ. Error bars correspond to the PSNR range for Λ. b) PSNR map of the estimated spectra for different σ noise and wavelengths λ.
values for USAF are lower than those for Cameraman. It appears that Cameraman is more complex for the considered algorithm than binary USAF.
Let us give some comments on the noise level in the intensity observations (4) and in the spectra |V (λ s )| 2 . These spectra are calculated using the Fast Fourier Transform (FFT) applied to the sequence J(z) as |V (λ s )| 2 = F F T (J(z))/N , where N is a length of this sequence, i.e. a number of shifts in z. It follows, that the noise in these estimate of |V (λ s )| 2 obtained by application of FFT to J noisy (z) is zero-mean white Gaussian, as in the observations but with the standard deviation equal to σ noise / √ N , which is the same for all spectral estimates |Ṽ (λ s )| 2 .
PSNR for observations is defined by (8), where in the nominator we have the maximum intensity calculated as the weighted sum of intensities of the spectral components. It follows that maximum spectral intensities |V (λ s )| 2 different for different λ take values lower than the observed intensity. For instance, if |V (λ s )| 2 have close values for all λ s then we can say that approximately |V (λ s )| 2 ∼ J(z)/N . (8), as well as the obtained above value for the noise standard deviation for the spectral estimate, results in the following formula for PSNR calculated in the spectral domain PSNR |V λs | 2 ∼ 10 log 10 max x,y,z (J(z)) σ noise · √ N .

Inserting this value in
Then, PSNRs calculated for the spectral estimates are smaller or much smaller than PSNR calculated for the observations. Thus, the noise level in the spectral domain is higher than that for the observations. It makes the problem of the noise removal for large N quite demanding and of the crucial importance for qualitative phase imaging.
This noise amplification phenomenon in the spectral domain is illustrated in Fig. 4 for our simulation data obtained for the Cameraman phase object.
As the spectra are varying as a function of wavelength, contrary to (9) we instead calculate PSNR according to the next formula: Here, the expression in the brackets is a ratio of the spectrum intensity for the wavelength λ s to the standard deviation corresponding to the noise in the spectrum estimates. Thus, PSNR is evaluated for each spectral component independently. Note that maximization in this definition of PSNR |V λs | 2 is restricted to the coordinates (x, y).
In Fig. 4, we show the mean value of PSNR |V λs | 2 calculated for the whole range of the spectra as well as the error bar showing variations of PSNR |V λs | 2 from minimal to maximal values. It is seen that PSNR in the spectral domain is much lower than that for the observations and this statement holds for all values of σ noise . It once more illustrates the noise amplification effect appeared in HS phase imaging. This effect is the interpretation of the well known in Fourier spectroscopy Felggett's disadvantage [22]. Figure 5 illustrates the crucial importance of iterations and filtering embedded in the developed algorithm. In Figs. 5(a-b) we show the results obtained by backward propagation of the spectral intensities obtained by FT from noisy observations. Actually, it is the initial iteration of our algorithm giving the first approximation of the objects. The comparison with the corresponding HSPR reconstructions in Figs. 5(c-d) is a clear demonstration of the performance of the algorithm and efficiency of its iterations.

Physical experiments
We have developed the HS lensless self-reference optical system shown in Fig. 6(a). The super-continuum laser, Λ = 470 − 2400 nm (YSL photonics CS-5), is used as a source of light. The light beam is split equally into two beams by beamsplitter BS1. The first beam reflected by the mirrors M1 and M2 goes undisturbed to the beamsplitter BS2. The second beam goes through the delay line M3-M4, obtaining the different optical paths with respect to the first beam. BS2 merges the beams together and after BS2 beams go through the object "O" to the registering camera "Cam" (FLIR Chameleon, 2448 × 2048, pixel size 3.45 µm).
The delay line is a piezo-based stage (Thorlabs NFL5DP20). The step-size of this stage defines the spectral resolution of the setup: the minimal step of the delay line, ∆z, should be at least twice smaller than the smallest wavelength, and the total moving distance of the delay line, Z, defines the spectral resolution for wavenumber k as ∆k = 1 2·Z = 1 2·∆z·N , here Z = N · ∆z and N is a number of steps of the stage.
In the provided experiments, the system parameters were ∆z = 59.7 nm and N = 1880, that correspond to ∆k = 44.6 cm −1 . The distance between the object and the camera was d = 16 mm. The used laser spectrum is wide, however, due to camera sensitivity we utilize only part of the spectral illumination, as it is shown in Fig. 6(b), the black dash curve demonstrates laser spectrum registered by a spectrometer (Thorlabs CCS200) and multiplied by quantum efficiency of the registration camera. The red curve in Fig. 6(b) is the spectrum |V (λ)| reconstructed by Fourier transform (2). This reconstruction is in a good agreement with the true spectral curve registered by the spectrometer. As a test object, we used PhaseFocus test target [23], which modulus and depth map are presented in Fig. 7(b,d), respectively, it is a phase object with an etched surface in fused silica glass, the depth of the etched features is 127 nm. The smallest groups from 6 to 9 have the smallest feature sizes of 10, 6.5, 3.5, 2.0 µm, respectively.
Reconstruction results by the HSPR algorithm are presented in Fig. 7(a,c), where modulus and depth map of the phase test target corresponding to λ = 687 nm are demonstrated. From these images, it is seen that HSPR provides good quality reconstructions for both modulus and phase for the whole test target. The HSPR reconstructed modulus image in Fig. 7(a) is coincide with the image provided by the manufacturer (Fig. 7(b)) with clearly resolved modulus features and almost not distinguishable features for phase regions. The difference of letters "K" and "F" in the bottom right corner is explained by the manufacturer as they provide the same image with letter "F" for all customers, however, each test target has its own distinguishing letter ("K" in our case). Other parameters are the same.
For clarity of the presented results, we provide Fig. 8 where the smallest groups 6,7,8,9 of the object presented and compared with depth and modulus reconstructed by digital holography (DH). Figure 8(a,b) shows the depth map and modulus of the object corresponding to λ = 687 nm; Fig. 8(c, d) -DH reconstruction, and Fig. 8(e) -a plot of depth cross-sections for the groups 6-9 reconstructed by HSPR (red solid curve) and DH (black dash), cross-sections are taken from lines with the same colors from Fig. 8(a) and Fig. 8(c), respectively. It is seen from the images that HSPR reconstructions correspond to the DH reconstructions and provided cross-sections show that reconstructed depth values equal to the real ones up to the elements 6 of the group 6 of the test target (region of 450 µm of cross-sections in Fig. 8(e)), which corresponds to 10 µm. Additionally, it is seen that the first four elements from group 8 appeared spatially resolved in terms of the Rayleigh criterion, but phase values are not correct. In that case, we may conclude that the given setup provides quantitative and qualitative phase estimations down to 10 µm resolution, and only qualitative results down to 4.5 µm.
Note, that the used for comparison DH system provides reference results based on a different optical setup with the use of lenses and holographic phase reconstruction algorithm [24]. As compared with the system studied in this paper, disadvantages of the reference DH system are in utilization of objective with the magnification of 40× and therefore an extremely small field of view ≈ 0.15 mm 2 (even for provided small image of groups 6-9 (Fig. 7(c)) we made stitching of 28 frames) and in lack of spectral information. a) b) Figure 6: (a) Hyperspectral phase retrieval setup. BS1-2 are beamsplitters, M1-4 are mirrors, "O" is a transparent object, "Cam" is a registering sensor, "B" is a light blocker. "Delay" is a moving delay stage. (b) Used spectrum: a black dash curve is a registered spectrum by a spectrometer and multiplied by camera quantum efficiency, a red solid curve is for the Fourier transform reconstructed spectrum.

Conclusion
A newly developed approach for hyperspectral phase retrieval in lensless self-referenced optical setup is introduced. It is based on the principles of the Fourier Transform spectroscopy and iterative phase retrieval algorithms. The essential part of the algorithm is developed for noise suppression based on sparse wavefront approximations. The provided approach is general for phase retrieval with HS illumination and can be applied for various self-reference schemes as for example [25], where self-referencing obtained by a birefringent delay line [26]. We provided a transmissive setup for investigation of transparent objects, however, this approach is applicable and for HS relief investigations in reflective setups.

Funding Information
Jane and Aatos Erkko Foundation and Finland Centennial Foundation: Computational Imaging without Lens (CIWIL) project.

Disclosures
The authors declare no conflicts of interest.