Hyperspectral phase retrieval: spectral–spatial data processing with sparsity-based complex domain cube filter

Vladimir Y. Katkovnik; Igor A. Shevkunov; Karen O. Eguiazarian

doi:10.1117/1.OE.60.1.013108

29 January 2021 Hyperspectral phase retrieval: spectral–spatial data processing with sparsity-based complex domain cube filter

Vladimir Y. Katkovnik, Igor A. Shevkunov, Karen O. Eguiazarian

Author Affiliations +

Optical Engineering, Vol. 60, Issue 1, 013108 (January 2021). https://doi.org/10.1117/1.OE.60.1.013108

Abstract

Hyperspectral (HS) imaging retrieves information from data obtained across broadband spectral channels. Information to retrieve is a 3D cube, where two coordinates are spatial and the third one is spectral. This cube is complex-valued with varying amplitude and phase. We consider shearography optical setup, in which two phase-shifted broadband copies of the object projections are interfering at a sensor. Registered observations are intensities summarized over spectral channels. For phase reconstruction, the variational setting of the phase retrieval problem is used to derive the iterative algorithm, which includes the original proximity spectral analysis operator and the sparsity modeling of the complex-valued object 3D cube. We resolve the HS phase retrieval problem without random phase coding of wavefronts typical for the most conventional phase retrieval techniques. We show the performance of the algorithm for object phase and thickness imaging in simulation and experimental tests.

1. Introduction

Hyperspectral (HS) imaging retrieves information from data collected across hundreds to thousands of spectral channels. Conventionally, these data are two-dimensional (2D) images stacked together in 3D cubes, where the first two coordinates are spatial $(x, y)$ and the third one is the spectral channel, which is usually represented by a wavelength $λ$ .

Over the last few years, HS coherent diffractive imaging (HSCDI) shows a strong progress, for instance, in high-resolution microscopy. Interference between reference and object wavefronts registered as intensity diffraction patterns is used conventionally in spectrally resolved interferometry,¹^,² Fourier transform holography,³^,⁴ and HS digital holography.²^,⁵^,⁶ Along with the traditional intensity imaging, HSCDI provides phase imaging, which brings additional information about an object under investigation, e.g., label-free multispectral refractive index estimation.⁷

In this paper, we consider a class of HSCDI in a shearography optical setup,⁸ where two coherent beams, the main (basic) and its time-delay copy, go through a transparent object simultaneously. Thus two identical but phase-shifted broadband patterns of the object are superimposed on the sensor. This scenario, typical for Fourier-transform spectrometry,⁹ leads to the phase loss problem, where a complex-valued 3D object should be reconstructed from indirect intensity observations as solutions of the ill-posed inverse problem. It results in serious amplification of measurement and processing noise. In addition, the high spectral resolution separates the energy obtained by an imaging sensor between many narrow wavebands, which results in small values of signal-to-noise ratio (SNR).

HS phase retrieval in the shear geometry has not been studied broadly. We have found only two papers dedicated to such a problem.¹⁰^,¹¹ Kalenkov¹⁰ realized phase retrieval by special optical system design, in which they produce spatial zero-order filtering of one beam and thanks to it the phase reconstruction comes down to a conventional HS holographic technique of direct phase reconstruction by Fourier transform. Reference 11 is our first solution to the HS phase retrieval problem in the lensless self-reference HS setup. There we perform phase reconstruction by a subsequent forward–backward propagation of HS cube slices by analogy to classical multiwavelength phase retrieval,¹² with the hypothesis that the object wavefronts produced by neighboring wavelengths are similar. We used this hypothesis for producing phase recalculation from one wavelength to the subsequent one.

In this paper, we propagate HS cube slices all together, without phase recalculation from one slice to another, as it is in Ref. 11. Such processing frees our method from the adjacent slice similarity hypothesis and by that extends the field of operation for objects with a sharp spectral response, where the subsequent slices might be unique.

The contribution of this paper can be summarized as follows. We formalize the HS phase retrieval as a variational optimization problem for noisy Gaussian observations. One of the key components of the derived iterative algorithm is an original proximal operator enabling both the spectral analysis of intensity observations and their denoising. Another important component of the algorithm is the original complex domain filtering of 3D complex-valued object arrays following from the proposed sparsity modeling for complex-valued 3D arrays. It is demonstrated by simulation experiments and processing of experimental data that the HS phase retrieval in the proposed setup can be resolved.

We have published the preliminary results about the proposed algorithm in the conference proceeding.¹³ In this paper, we present more stimulation and as well as experimental results with detailed discussion.

2. Problem Formulation

Let $U_{o} (x, y, k) \in C^{n \times m}$ be a 2D slice of the complex-valued 3D HS cube with spatial coordinates $(x, y)$ , a spectral component $k$ , and $Q_{K} (x, y) = {U_{0} (x, y, k), k \in K}$ , $Q_{K} \in C^{n \times m \times l_{K}}$ , be the total spectral cube composed of $l_{K}$ spectral 2D slices. The size of the cube is $n \times m \times l_{K}$ .

The lines of $Q_{K} (x, y)$ contain $l_{K}$ spectral observations corresponding to the coordinate $(x, y)$ . This is a HS model of the object $U_{0}$ .

The squared magnitude (intensity) of observations may be written as

Eq. (1)

Y_{t} = \sum_{k \in K} {| U_{t, k} |}^{2}, U_{t, k} = A_{t, k} U_{o, k}, t \in T .

Here we use vector representations for the slices $U_{0} (x, y, k)$ , $U_{o, k} \in C^{N}$ , $N = n m$ , and $A_{t, k} \in C^{M \times N}$ are linear operators of image formation modeling the propagation of 2D object wavefronts from the object plane to the sensor plane.

Thus $U_{t, k} = A_{t, k} U_{o, k} \in C^{M}$ are the complex-valued signals, which intensities are registered by the sensor as $Y_{t} \in R_{+}^{M}$ . The summation on $k$ in Eq. (1) says that the measurements are intensities calculated over the entire spectral components. The squared absolute value ${| \cdot |}^{2}$ in Eq. (1) is element-wise for components of the vectors $U_{t, k}$ . For the additive noise case, the intensity observations become

Eq. (2)

Z_{t} = Y_{t} + ε_{t}, t \in T,

where

ε_{t}

stays for a Gaussian noise ∼

N (0, σ)

.

The HS phase retrieval is formulated as a reconstruction of the HS object cube $Q_{K} (x, y)$ from the intensity observations $Y_{t}$ or $Z_{t}$ , $t \in T$ .

A dependence of the object $U_{o} (x, y, k)$ on the spectral $k$ is a specific feature of the considered phase retrieval setup as compared with the phase retrieval for a monochromatic object $U_{o}$ with observations $Z_{t} = {| A_{t} U_{o} |}^{2} + ε_{t}$ , or with the recent multispectral phase retrieval¹⁴ from observations $Y_{t} = \sum_{k \in K} {| A_{t, k} U_{o} |}^{2}$ , where the object $U_{o}$ is again invariant in the spectral domain and only the propagation operator is varying on $k$ .

In this paper, we restrict a class of operators $A_{t, k}$ to the form appeared in Fourier transform spectrometry¹⁵ with the intensities $Y_{t}$ in the form:

Eq. (3)

Y_{t} = \sum_{k \in K} {| (1 + e^{- j 2 π k t}) A_{k} U_{o, k} |}^{2},

Eq. (4)

A_{t, k} = (1 + e^{- j 2 π k t}) A_{k}, t \in T .

Due to this $A_{t, k}$ , two identical but phase-shifted copies of the object wavefront $U_{o, k}$ are superimposed on the sensor plane: main $A_{k} U_{o, k}$ and phase-shifted $e^{- j 2 π k t} A_{k} U_{o, k}$ .

3. Algorithm Development

3.1.

Approach

We assume the following for sampling on $t$ and $k$ in the above observation model, $K = 0,1, \dots, N / 2 - 1$ , $T = 0,1, \dots, N - 1$ , then

Eq. (5)

Y_{t} = 2 \sum_{k \in K} [1 + \cos (\frac{2 π}{N} k t)] {| B_{k} |}^{2}, t \in T,

where

B_{k} = A_{k} U_{o, k}

,

N

is a number of observations and the integer discrete frequency

k

covers the low-frequency interval

{0,1, \dots, N / 2 - 1}

.

The restriction of this frequency interval to length $N / 2$ follows from the periodicity on $t$ of the observation function Eq. (5). Provided that ${| B_{k} |}^{2} = 0$ for $k = N / 2, \dots, N - 1$ , the observations ${Y_{t}}$ uniquely define the intensity spectrum ${| B_{k} |}^{2} \in K$ (see Ref. 16).

Let us replace the noisy observations $Z_{t}$ by the differences $Δ Z_{t} = Z_{t} - \frac{1}{N} \sum_{t = 0}^{N - 1} Z_{t}$ . For large $N$ , the noise level in the mean value $\frac{1}{N} \sum_{t = 0}^{N - 1} Z_{t}$ is of the order smaller than that in $Z_{t}$ , and we can assume that

Eq. (6)

Δ Z_{t} = Z_{t} - \frac{1}{N} \sum_{t = 0}^{N - 1} Y_{t} = 2 \sum_{k = 1}^{N / 2 - 1} \cos (\frac{2 π}{N} k t) {| B_{k} |}^{2} + ε_{t},

since

\frac{1}{N} \sum_{t = 0}^{N - 1} \cos (\frac{2 π}{N} k t) = δ_{k, 0}

, where

δ_{k, 0}

is the Kronecker symbol.

Then the criterion proposed for the algorithm design can be given in the form:

Eq. (7)

J = \frac{1}{σ^{2}} \sum_{t = 0}^{N - 1} {‖ Δ Z_{t} - 2 \sum_{k = 1}^{N / 2 - 1} {| B_{k} |}^{2} \cos (\frac{2 π}{N} k t) ‖}_{2}^{2} + \frac{1}{γ} \sum_{k = 1}^{N / 2 - 1} {‖ B_{k} - A_{k} U_{o, k} ‖}_{2}^{2} + f_{reg} ({U_{o, k}}_{1}^{N / 2 - 1}) .

The first summand stays for the minus log-likelihood provided a Gaussian noise. The second summand $J$ in the criterion is penalizing the difference between $B_{k}$ and $A_{k} U_{o, k}$ with the weight $1 / γ$ , $γ > 0$ . The last third item is a penalty using the sparsity hypothesis on the cube of the object images ${U_{o, k}}_{1}^{N / 2 - 1}$ . The object reconstruction is obtained by minimizing $J$ on ${U_{o, k}}_{1}^{N / 2 - 1}$ . The wavefronts at the sensor plane $B_{k}$ serve as splitting variables in this optimization, separating the observations from the object images.

Note that the zero-order DC term of the object $U_{o, 0}$ is dropped from the consideration of Eq. (7) as we are interested only in higher-order components of the spectrum.

The developed algorithm iterates $\min_{{B_{k}}} J$ on $B_{k}$ provided given ${U_{o, k}}_{1}^{N / 2 - 1}$ and $\min_{{U_{o, k}}} J$ provided given ${B_{k}}$ . Minimization on $B_{k}$ concerns the two first summands in Eq. (7). We show that this minimization is implemented by the derived proximity operator explicitly separating the spectral wavefronts $B_{k}$ at the sensor plane. Minimization on $U_{o, k}$ concerns the last two summands in Eq. (7). In this paper, we prefer to do not to formalize the penalty function $f_{reg} ({U_{o, k}}_{1}^{N / 2 - 1})$ but to use the special high-quality complex domain filter. This approach is in line with the recent tendency in inverse imaging to use high-quality filters instead of regularizers formally obtained by solving optimization problems.

3.2.

Minimization on $B_{l}$

For minimization $\min_{{B_{l}}} J$ , we solve the equations $\partial J / \partial B_{l}^{*} = 0$ , $l = 1, \dots, N / 2 - 1$ , what leads to the set of the complex-valued equations for $B_{l}$ :

Eq. (8)

[\frac{- 4}{σ^{2}} \sum_{t = 0}^{N - 1} Δ Z_{t} \cos (\frac{2 π}{N} l t) + \frac{4 N}{σ^{2}} {| B_{l} |}^{2} + \frac{1}{γ}] B_{l} = \frac{1}{γ} A_{l} U_{o, l} .

Derivation of Eq. (8) can be seen in Appendix A. Inserting $B_{l} = | B_{l} | e^{j φ_{B_{l}}}$ in Eq. (8), we may conclude that

Eq. (9)

φ_{B_{l}} = φ_{A_{l} U_{o, l}},

where

φ_{A_{l} U_{o, l}}

is the phase of

A_{l} U_{o, l}

.

This result follows from Eq. (8) because the expression in square brackets is real valued. Further, the magnitude $| B_{l} |$ is defined as a non-negative solution of the polynomial equation, cubic with respect to $| B_{l} |$ :

Eq. (10)

[\frac{- 4}{σ^{2}} \sum_{t = 0}^{N - 1} Δ Z_{t} \cos (\frac{2 π}{N} l t) + \frac{4 N}{σ^{2}} {| B_{l} |}^{2} + \frac{1}{γ}] | B_{l} | - \frac{1}{γ} | A_{l} U_{o, l} | = 0 .

There is no second power in this equation and the free term is negative, then there is only one positive solution $| B_{l} |$ ,¹⁷ which is computed by Cardano’s formula.

The calculations in Eqs. (8)–(10) are produced in the pixel-wise manner for each pixel separately: the amplitudes for $B_{l}$ are given by Eq. (10) and the phases by Eq. (9).

This solution of $\min_{{B_{l}}} J$ can be interpreted as the proximity operator¹⁴^,¹⁸ with a compact notation:

Eq. (11)

B_{l} = {prox}_{f γ} (A_{l} U_{o, l}), l = 1, \dots, N / 2 - 1,

where

γ > 0

is a parameter of the quadratic item in Eq. (7) and

f

stays for the minus log-likelihood part of the criterion Eq. (7),

f ≜ \frac{1}{σ^{2}} \sum_{t = 0}^{N - 1} {‖ Δ Z_{t} - 2 \sum_{k = 1}^{N / 2 - 1} {| B_{k} |}^{2} \cos (\frac{2 π}{N} k t) ‖}_{2}^{2}

.

The proximity solution ${B_{l}}$ resolves two problems. First, complex domain spectral components $B_{l}$ are extracted from the intensity observations. Thus we yield the spectral analysis of the observed total intensities. Second, the noisy observations are filtered with the power controlled by the parameter $γ$ compromising the noisy observations $Z_{t}$ and the power of the predicted signal $A_{l} U_{o, l}$ at the sensor plane.

The calculation of the cosine transform in Eqs. (8) and (10) can be produced using FFT (Proposition 2 in Ref. 16). In MATLAB notation, it looks as follows:

Eq. (12)

\frac{1}{N} \sum_{t = 0}^{N - 1} Δ Z_{t} \cos (\frac{2 π}{N} k t) = real [ifft 2 (Δ Z_{t \in T})], for k = 1,2, \dots, N / 2 - 1,

where

Δ Z_{t \in T}

is a sequence

Δ Z_{t}

,

t = 0, \dots, N - 1

and

ifft 2

stays for the 2D inverse FFT.

This equation is valid for the nonsymmetric sampling interval $t \in T = {0,1, \dots, N - 1}$ . For the symmetric $T = {- N / 2, - N / 2 + 1, \dots, N / 2 - 1}$ , the formula with FFT is different (Proposition 4 in Ref. 16):

Eq. (13)

\frac{1}{N} \sum_{t = - N / 2}^{N / 2 - 1} Δ Z_{t} \cos (\frac{2 π}{N} k t) = real {ifft 2 [fftshift (Δ Z_{t \in T})]}, for k = 1,2, \dots, N / 2 - 1,

where “fftshift” denotes the MATLAB shift operator in frequency domain. This symmetric sampling is often used in applications. In the noiseless scenario, Eqs. (12) and (13) give precise values of

{| B_{k} |}^{2}

, the intensity spectrum of observations.

3.3.

Regularization by Sparsity-Based Filters

Let $U_{o} (x, y, k)$ be a transfer function of a transparent thin object. The phase of this object can be written in the form:¹⁹

Eq. (14)

φ_{U_{o}} (x, y, λ) = \frac{2 π}{λ} h (x, y) (n_{λ} - 1),

where

h (x, y)

is the object thickness invariant on wavelength

λ

, and the phase of

φ_{U_{o}}

is smooth with respect to wavelength

λ

in the denominator of this equation for a visual interval of wavelengths. The refractive index

n_{λ}

is a slowly varying function of wavelength for many materials. Thus the phase of

φ_{U_{o}}

as defined by Eq. (14) is a slowly varying function of wavelength and the object slices

U_{o} (x, y, k) = b_{U_{o}} (x, y) \exp (j φ_{U_{o}})

close to each other for nearby wavelengths. Here

b_{U_{o}}

is an amplitude of the slice.

This simple example, being of a general nature, allows to assume that on many occasions the object slices $U_{o} (x, y, k)$ are slowly varying functions of wavelength, which means that these slices are similar for nearby $k$ . Then the spectral lines of $Q_{k} (x, y)$ live in $l_{k}$ -dimensional subspaces with $l_{k} ≪ l_{K}$ . Therefore, the concept of sparsity in the spectral domain can be applied for the modeling of this phenomenon.

In this paper, we exploit the spectral sparsity hypothesis in the special filter designed for joint filtering of spectral slices of the complex-valued object cube. This filter is used as a regularizer of the inverse problem instead of formalizing the last regularization item in the criterion $J$ .

This approach is in-line with the recent progress in computational imaging where efficient plug-and-play filters have been recognized as powerful tools for prior and regularization to resolve inverse imaging problems.²⁰^,²¹

The following algorithm is exploited for the joint sparsity-based processing of 3D HS cubes:

Eq. (15)

{{\hat{U}}_{o, k}, k \in K} = CCF {U_{o, k}, k \in K} .

Complex domain cube filter (CCF) processes 3D cube data ${U_{o, k}, k \in K}$ jointly and provides 3D estimates ${{\hat{U}}_{o, k}, k \in K}$ for all $k \in K$ . The CCF algorithm starts from the singular value decomposition (SVD) analysis of the HS cube. As a next step, preliminary estimates of complex-valued signal and noise correlation matrices are used to select a small-sized subset (eigenimages) of the SVD eigenvectors that best approximate the signal subspace in the least square error sense. Further, 2D complex domain filtering is applied to the denoising of this small number of eigenimages. The complex domain block-matching 3D algorithm²²^,²³ is used for this filtering of eigenimages. Finally, the filtered eigenimages are used to reconstruct the estimates for the entire 3D cube, i.e., for $k \in K$ . Overall, CCF can be classified as the adaptive SVD algorithm with optimal selection of a small number of eigenimages.

This CCF algorithm is a complex domain modification of the fast algorithms developed for real-valued HS observations in Refs. 24 and 25. The modification of this algorithm for complex domain HS data is produced recently in Refs. 26 and 27. A sliding window version of CCF was developed for objects with discontinuous and fast varying spectral characteristics.²⁷

3.4.

HS Phase Retrieval Algorithm

Table 1 presents a block scheme of the developed algorithm. The calculation of the initial spectral guess $U_{o, k}^{(0)}$ uses the complex-valued 2D object model with the wavelength-varying phase according to Eq. (14) obtained from an uniform initial guess for $h (x, y)$ . Stage 1 is the forward propagation from the object plane to the sensor plane according to the operator $A_{k}$ . In stage 2, the proximity operator produces the spectral analysis of the intensity observations and provides the updated version of the wavefront at the sensor plane defined as $B_{k}^{(s)}$ . In stage 3, the backward propagation from the sensor plane to the object plane gives the update for the object spectral image $U_{o, k}^{(s)}$ . In stage 4, the CCF algorithm produces the joint filtering of the object spectral estimates using optimal SVD eigenimages and updates the object cube estimate as $U_{o, k}^{(s + 1)}$ . These iterations are repeated $\max_{iter}$ times.

Table 1

HS phase retrieval algorithm.

	Initialization: Uo,k(0), k∈K;
	Main iterations: For $s = 1,2, \dots, \max_{iter}$ do
1.	Forward propagation:
1.	$U_{t, k}^{(s)} = A_{k} U_{o, k}^{(s)}$ , $k \in K$ ;
2.	Proximity operation:
2.	$B_{k}^{(s)} = {prox}_{f γ} [U_{t, k}^{(s)}]$ , $k \in K$ ;
3.	Backward propagation:
3.	$U_{o, k}^{(s)} = A_{k}^{#} B_{k}^{(s)}$ , $k \in K$ ;
4.	Spectral filtering by CCF:
4.	${U_{o, k}^{(s + 1)}, k \in K} = CCF {U_{o, k}^{(s)}, k \in K}$ ;
	Output:
	$U_{o, k}^{(\max_{iter} + 1)}$ , $k \in K$ .

The backward propagation operator $A_{k}^{#}$ at stage 3 is an inverse operator for the forward propagation $A_{k}$ .

The observation noise is filtered at stage 2 by the proximity operator, whereas the noisy components in the object reconstructions are filtered by CCF.

3.5.

On Interpretation of Hyperspectral Complex Domain Data

Intensity HS imaging in remote sensing technology mainly works with reflected signals and enables a lot of applications requiring fine identification of materials and physical parameters. HS complex domain sensing is of different nature and of different interpretations. Let us discuss this issue assuming that the specimen of interest is thin and transparent. The proposed and discussed algorithm estimates both the amplitude and phase of the specimen and both estimates are spectral. The amplitude gives information on specimen transparency and material density. This data can be interpreted in a manner similar to that used for HS intensity data.

The reconstructed (measured, imaged) phase, in general, differs essentially from Eq. (14) as this object phase is replaced by the corresponding wrapped phase $wrap [φ_{U_{o}} (x, y, λ)]$ , where $wrap (φ) = \mod (φ + p i, 2 p i) - p i$ is the operator wrapping the phase $φ$ to the interval $[- π, π)$ .

In what follows, we assume that the phase $φ_{U_{o}} (x, y, λ)$ belongs to the basic phase interval $[- π, π)$ , then the phase measurements are indeed equal to the object phase. This assumption simplifies the discussion. The reconstructed HS phase $φ_{U_{o}} (x, y, λ)$ according to Eq. (14) defines a product of two specimen parameters: thickness $h$ and refractive index $n_{λ}$ :

Eq. (16)

φ_{U_{o}} (x, y, λ) \frac{λ}{2 π} = h (x, y) (n_{λ} - 1) .

Calculating in Eq. (16) the averages of spatial variables $(x, y)$ and spectral $λ$ , we may separate estimates of thickness and refractive index in the following way:

Eq. (17)

{mean}_{λ} [φ_{U_{o}} (x, y, λ) \frac{λ}{2 π}] = h (x, y) c_{n}, {mean}_{(x, y)} [φ_{U_{o}} (x, y, λ) \frac{λ}{2 π}] = n_{λ} c_{h},

where

c_{n} = {mean}_{λ} (n_{λ} - 1)

and

c_{h} = {mean}_{h} [h (x, y)]

are the constants.

These equations give estimates for variations of $h (x, y)$ and $n_{λ}$ within invariant factors $c_{n}$ , $c_{h}$ . These reconstructions can be useful for many applications, e.g., for phase plate calibration to find variations in a surface and in the optical properties of materials.

If $h (x, y)$ or $n_{λ}$ are given, the estimates for $h (x, y)$ and $n_{λ}$ , respectively, take the form:

Eq. (18)

h (x, y) = \frac{1}{2 π} {mean}_{λ} [φ_{U_{o}} (x, y, λ) λ / (n_{λ} - 1)], n_{λ} - 1 = \frac{1}{2 π} {mean}_{(x, y)} [φ_{U_{o}} (x, y, λ) λ / h (x, y)] .

A much more complex case occurs when we need to measure the refractive index varying spatially. Then $φ_{U_{o}} (x, y, λ) = \frac{2 π}{λ} h (x, y) [n_{λ} (x, y) - 1]$ , and for estimation of refractive index $n_{λ} (x, y)$ we need measurements or extra information on object thickness.

Provided a fixed specimen thickness, we have a refractive index estimate in the explicit form $n_{λ} (x, y) = 1 + \frac{1}{2 π} φ_{U_{o}} (x, y, λ) λ / h (x, y)$ .

If the thickness is slowly varying on $(x, y)$ , then the windowed on spatial coordinates HS phase estimates can be used to image fast variations in $n_{λ} (x, y)$ .

All estimates discussed in this section are implemented as a postprocessing of HS cube phase $φ_{U_{o}} (x, y, λ)$ obtained by the proposed algorithm.

If the variations of the object’s phase go beyond the basic interval $[- π, π)$ , the wrapping effect becomes essential and the above equations for the thickness and refractive index do not work. The phase unwrapping²⁸ is one instrument to resolve the problem, otherwise, a proper interpretation of phase imaging, even visual, without attempts to separate thickness and refractive index, may become problematic. As far as we may conclude from publications, modern phase imaging, at least in biomedicine, is restricted to phase imaging in the basic interval $[- π, π)$ .²⁹

4. Results

Simulation and experimental studies of the algorithm performance have been produced for various system parameters and various objects to be reconstructed. Some of these results are shown and discussed in what follows.

4.1.

Simulation Tests

We model a transparent phase specimen assuming an invariant amplitude and phase defined by Eq. (14), where $h (x, y)$ is varying according to the USAF test image [see Figs. 1(a) and 1(b)]. The refractive index $n_{λ}$ is calculated for each $λ$ according to Cauchy’s equation ³⁰ with coefficients taken for the glass BK7.³¹ As an illumination source, we model a broadband light beam in the range of $Λ = [450 : 900] nm$ with a uniform intensity distribution. The number of wavelengths is 250. The light beam goes through the specimen and after freely propagates to the sensor where the intensity observations are registered as a 3D cube $Y_{t} (x, y)$ of length $N = 2000$ , which corresponds to $N$ steps of the phase-shift modulation Eq. (5). We implement this modulation according to the simulated motion of the phase delay stage of the Fourier spectrometer. Each step of this simulated motion was equal to 100 nm to resolve the smallest wavelengths of the spectral range of the light beam.¹⁵ The wavelength-dependent image formation operator $A_{k}$ is calculated according to the angular spectrum propagation technique.¹⁹ A distance from the specimen to the sensor was equal to 10 mm and the camera pixel size was $3.45 μ m$ .

Fig. 1

The true thickness of the transparent phase object: (a) 2D and (b) 3D images; (c) 3D noisy intensity observations: diffractive data cube $Z_{t}$ , $t \in T$ ; and (d) the intensity distribution for the central pixel of the 3D data cube $Z_{t}$ as a function of the experiment number $t$ (slice number in the data cube).

The proximity HS analysis is produced for 250 wavelengths. However, we use only 50 of these wavelengths (range [680:820] nm) corresponding to the forthcoming physical experiments, where only spectra of higher signal-to-noise ratio are exploited for phase retrieval processing. Noisy observations are created according to Eq. (2).

The registered noisy observation cube $Z_{t} (x, y)$ is demonstrated in Fig. 1(c) and a single-pixel intensity distribution along $t$ coordinate is presented in Fig. 1(d). We reconstruct a thickness map of the object for 50 wavelengths and the results are presented as HS cube $64 \times 64 \times 50$ . The accuracy of the reconstruction is characterized by the relative-root-mean-square error (RRMSE) criterion calculated for each wavelength:

Eq. (19)

RRMSE = \frac{\sqrt{{‖ {\hat{h}}_{est} - h_{true} ‖}_{2}^{2}}}{\sqrt{{‖ h_{true} ‖}_{2}^{2}}},

where

{\hat{h}}_{est}

and

h_{true}

are the reconstructed and true thickness maps, respectively. RRMSE values

< 0.1

correspond to visually high-quality 2D imaging.

We perform experiments with an additive Gaussian noise of different $σ$ to test the algorithm’s robustness to noise. The noise level with respect to the spectral signals at the sensor plane is characterized by the peak-signal-to-noise ratio (PSNR) calculated as

Eq. (20)

PSNR = 10 \log_{10} [\frac{\max_{x, y, k} (| B_{k} |)}{σ}],

where

\max_{x, y, k} (| B_{k} |)

is the maximum of the intensity spectra at the sensor plane.

The estimates of the object thickness obtained from the phases of the entire HS cube according to Eq. (18) are presented in Fig. 2 as 2D and 3D images. A nearly perfect correspondence of these estimates to the true thickness images in Fig. 1 is clear. In these tests: iteration number $s = 30$ , $σ = 0$ , and $σ = 15$ . The low values of RRMSEs prove the high-accuracy performance of the algorithm achieved for a small number of iterations.

Fig. 2

Reconstruction of the object thickness obtained by the HS phase retrieval, iteration number $s = 30$ , (a), (b) HSPR, $PSNR = 29.8 dB$ , $RRMSE = 0.031$ , $σ = 0$ and (c), (d) HSPR, $PSNR = 12.2 dB$ , $RRMSE = 0.11$ , $σ = 15$ .

Let us compare two ways of thickness reconstruction: from the phase estimate obtained for a single wavelength and from the reconstructed entire complex-valued cube where the data for all wavelengths are used simultaneously. In the first case, the thickness is calculated by Eq. (16) and, in the second case, according to Eq. (18).

A red dashed curve in Fig. 3(a), the low-noise case ( $PSNR = 29.8 dB$ ) is obtained for separate estimates of thickness from single-wavelength phase estimates. The minimum value of RRMSE is achieved near the middle point of the wavelength interval. The solid-black curve shows RRMSE achieved by the estimate obtained for thickness of the entire reconstructed phase cube. Naturally, this RRMSE value is invariant on wavelength. Note that this RRMSE value is quite near to the smallest value achieved by the red-dashed curve.

Fig. 3

RRMSE curves for reconstructed thicknesses from observations: (a) $σ = 0$ ( $PSNR = 29.8 dB$ ) and (b) $σ = 15$ ( $PSNR = 12.2 dB$ ). Dash red curves are for reconstructed objects’ thickness from separate wavelengths and solid black curves are for reconstructed objects’ thickness averaged along the whole range of wavelengths.

The curves in Fig. 3(b) are obtained for much noisier observations ( $PSNR = 12.2 dB$ ). The reconstruction by Eq. (18) (solid-black line) is wavelength invariant and shows higher RRMSE as compared with the corresponding curve in Fig. 3(a). A red-dashed curve is varying with a lot of up and down peaks. This behavior of the estimates differs from what we had in Fig. 3(a). We may conclude that for the noisier data, the separate thickness estimates are very noisy and unstable with wavelength. At the same time, the estimate based on the joint use of all wavelength estimates gives the results that can be treated as much better than those for most wavelengths. As we do not know in advance the best wavelength for thickness estimation, the estimate [Eq. (18)] guarantees the high accuracy for thickness reconstruction for both noisy and noiseless cases.

The fast convergence rate of the algorithm is demonstrated in Fig. 4 showing RRMSE as a function of iteration number $s$ for different noise levels and different wavelengths. Figure 4(a) presents RRMSE values averaged over $λ$ for the entire HS cube depending on the given noise level, measured by PSNR, Eq. (20). Figure 4(b) shows the RRMSE values for $PSNR = 12.2 dB$ and for each wavelength. These images prove that the algorithm is robust to noise and provides eligible low RRMSE values even for very noisy observations of $PSNR = 12.2 dB$ . Due to the spectral uniformity of the modeled illumination source and the joint HS cube processing with CCF, RRMSEs for different wavelengths achieve similar values. For both RRMSE maps, the dark blue regions correspond to high-quality reconstructions.

Fig. 4

RRMSE maps for the reconstructions in different noise conditions. (a) Map of RRMSE values averaged over $λ$ as function of iteration number $s$ and PSNR and (b) RRMSE for each $λ$ depending on iteration number $s$ for the case of $PSNR = 12.2 dB$ .

4.2.

Physical Tests

For experimental validation, we consider a biological specimen, a fly wing that has both amplitude and phase-varying features. The experimental data were obtained by the lensless optical system implementing the principles of Fourier spectrometry with a supercontinuum laser source $Λ = [650 : 900] nm$ and a piezo-driven stage in the delay line.¹¹ The optical setup is presented in Fig. 5(a), where the interferometric scheme is realized. Figure 5(b) shows the spectrum of the laser source registered by the commercial spectrometer (Thorlabs CCS200), dashed blue curve, and reconstructed by a Fourier transform, solid orange curve. For the commercial spectrometer curve, we applied the quantum efficiency of the camera. These two curves are in close agreement with each other.

Fig. 5

(a) HS phase retrieval setup. BS1 and BS2 are beamsplitters, M1 to M4 are mirrors, “O” is a transparent specimen, “Cam” is a registering sensor, “B” is a light blocker, and “Delay” is a moving delay stage. (b) Used light spectrum: a blue-dashed curve is a spectrum registered by a spectrometer and multiplied by the camera quantum efficiency. An orange solid curve is a spectrum reconstructed by our algorithm.

We recorded $l_{K} = 1880$ observations $Z_{t} (x, y)$ with the phase-delaying step of 59.7 nm. The total distance of the delay line defines the spectral resolution of the wavenumber as $Δ k = 44.6 {cm}^{- 1}$ . Due to the non-uniformity of the laser spectrum, the object reconstruction is performed for only 50 wavelengths with relatively high-intensity spectra in the range $Λ = [681 : 802] nm$ . The results of the amplitude and phase reconstruction for slice $λ = 736 nm$ are shown in Fig. 6. It is clearly seen that both amplitude and phase features are well distinguishable in the provided images. The wrapped phase regions appear in the phase image (sharp black lines), however, they do not corrupt the reconstructions since the algorithm CCF provides filtering for complex-domain wavefronts, where no distinction exists between absolute and wrapped phases. The HS cube of phase/amplitude reconstruction for each slice (wavelength) is presented in Video 1.

Fig. 6

The amplitude and phase reconstructions corresponding to the $λ = 736 nm$ slice of the 3D cube. The amplitudes and phases for the entire reconstructed HS cube can be seen in Video 1 (Video 1, MP4, 1.193 MB [URL: https://doi.org/10.1117/1.OE.60.1.013108.1]).

For demonstration of the spatial differences in the spectral responses of the specimen, we show phase and amplitude spectra obtained for three different points of the specimen [see Figs. 7(a) and 7(b)]. Each of the spectral curves corresponds to one of the numbered colored points in the amplitude image [Fig. 7(c)] with the same corresponding colors.

Fig. 7

(a) Amplitude and (b) phase spectra for different points of the object. The curves in the left and middle plots correspond to the points numbered in (c) the amplitude image, $λ = 745 nm$ with the same colors. Black circles curve in plot (b) is a phase difference between curves from points 2 and 3.

The point “1” is located out of the specimen area and, therefore, the amplitude spectrum in Fig. 7(a) represents the laser spectrum multiplied on the registration camera quantum efficiency. The corresponding phase spectrum curve in Fig. 7(b) is nearly invariant on wavelength.

The points “2” and “3” belong to the specimen and reflect different spectrum absorbance and disturbances in different areas of the wing. These spatial variations in phase and amplitude cannot be interpreted unambiguously as variations in both thickness and refractive index can cause them.

The black circle curve in Fig. 7(b) is a difference of the phase curves for points “2” and “3.” It is growing approximately linearly with a growing wavelength. If we assume that the wavelength-dependent refractive index is the same for these two points, such a behavior of this curve indicates that the object in point 2 is thinner than in point 3. In this way, sometimes the ambiguity of thickness/refractive index can be overcome.

For the experiments, we used MATLAB R2020a on a computer with 32 GB of RAM and 1.9 GHz Intel^® Core™ i7-08665U CPU. A single iteration for HS cube with dimensions $1500 \times 1500 \times 50$ took around 150 s.

5. Conclusion

An HS phase retrieval algorithm has been developed and presented, which allows the reconstruction of the complex amplitude of a transparent specimen in the lensless optical setup based on Fourier spectrometry principles in a shearographic geometry. This self-referencing lensless system enables a large field of view and robustness of the results with respect to vibration. A specimen phase lost in intensity measurements is reconstructed by the developed iterative algorithm utilizing HS complex domain sparsity of wavefronts. The advantage of HS data processing compared to the results available for monochromatic experiments is demonstrated. With the growing interest in HS measurements, we expect that the developed algorithm finds applications in various bio- and medical tasks concerning noninvasive quantitative phase imaging.

6. Appendix A: Minimization on $B_{l}$

For minimization $\min_{{B_{l}}} J$ , we calculate the derivative $\partial J / \partial B_{l}^{*}$ and solve the equation $\partial J / \partial B_{l}^{*} = 0$ what leads to the complex valued solution of the form $B_{l} = | B_{l} | e^{j φ_{B_{l}}}$ .

The following manipulations define the derivatives $\partial J / \partial B_{l}^{*}$ , $l = 1,2, \dots, N / 2 - 1$ :

Eq. (21)

\partial J / \partial B_{l}^{*} = \frac{- 2}{σ^{2}} \sum_{t = 0}^{N - 1} [Δ Z_{t} - 2 \sum_{k = 1}^{N / 2 - 1} {| B_{k} |}^{2} \cos (\frac{2 π}{N} k t)] [2 \cos (\frac{2 π}{N} l t) B_{l}] + \frac{1}{γ} (B_{l} - A_{l} U_{o, l}) = \frac{- 4}{σ^{2}} \sum_{t = 0}^{N - 1} Δ Z_{t} \cos (\frac{2 π}{N} l t) B_{l} + \frac{8}{σ^{2}} \sum_{k = 1}^{N / 2 - 1} {| B_{k} |}^{2} \frac{N}{2} δ_{k, l} B_{l} + \frac{1}{γ} (B_{l} - A_{l} U_{o, l}) = \frac{- 4}{σ^{2}} \sum_{t = 0}^{N - 1} Δ Z_{t} \cos (\frac{2 π}{N} l t) B_{l} + \frac{4 N}{σ^{2}} {| B_{l} |}^{2} B_{l} + \frac{1}{γ} (B_{l} - A_{l} U_{o, l}) .

It is used in the second line of these equations that

\sum_{t = 0}^{N - 1} \cos (\frac{2 π}{N} k t) \cos (\frac{2 π}{N} l t) = \frac{N}{2} δ_{k, l}, for k, l = 1, \dots, N / 2 - 1,

where

δ_{k, l}

is the Kronecker’s symbol,

δ_{k, l} = 1

for

k = l

and otherwise

δ_{k, l} = 0

.

It follows from the third line in Eq. (22) that the complex-valued $B_{l}$ is a solution of Eq. (8).

Acknowledgments

This research was supported by Jane and Aatos Erkko Foundation, Finland, and Finland Centennial Foundation: Computational Imaging without Lens (CIWIL) project. The authors declare no conflict of interest. The founding sponsors had no role in the design of the study; in the collection, analysis, or interpretation of data; in the writing of the manuscript, and in the decision to publish the results.

References

1.

C. Dorrer et al., “Spectral resolution and sampling issues in Fourier-transform spectral interferometry,” J. Opt. Soc. Am. B, 17 (10), 1795 –1802 (2000). https://doi.org/10.1364/JOSAB.17.001795 JOBPDE 0740-3224 Google Scholar

2.

S. Kalenkov, G. Kalenkov and A. Shtanko, “Spectrally-spatial Fourier-holography,” Opt. Express, 21 (21), 24985 –24990 (2013). https://doi.org/10.1364/OE.21.024985 OPEXFF 1094-4087 Google Scholar

3.

V. T. Tenner, K. S. Eikema and S. Witte, “Fourier transform holography with extended references using a coherent ultra-broadband light source,” Opt. Express, 22 (21), 25397 –25409 (2014). https://doi.org/10.1364/OE.22.025397 OPEXFF 1094-4087 Google Scholar

4.

S. Eisebitt et al., “Lensless imaging of magnetic nanostructures by x-ray spectro-holography,” Nature, 432 (7019), 885 (2004). https://doi.org/10.1038/nature03139 Google Scholar

5.

D. N. Naik et al., “Spectrally resolved incoherent holography: 3D spatial and spectral imaging using a Mach–Zehnder radial-shearing interferometer,” Opt. Lett., 39 1857 –1860 (2014). https://doi.org/10.1364/OL.39.001857 Google Scholar

6.

D. Claus et al., “Accuracy enhanced and synthetic wavelength adjustable optical metrology via spectrally resolved digital holography,” J. Opt. Soc. Am. A: Opt. Image Sci. Vision, 35 (4), 546 –552 (2018). https://doi.org/10.1364/JOSAA.35.000546 Google Scholar

7.

Á. B. Peña et al., “Refractive index properties of the retina accessed by multi-wavelength digital holographic microscopy,” Proc. SPIE, 10883 108830X (2019). https://doi.org/10.1117/12.2508976 Google Scholar

8.

Y. Hung, “Shearography: a new optical method for strain measurement and nondestructive testing,” Opt. Eng., 21 (3), 213391 (1982). https://doi.org/10.1117/12.7972920 Google Scholar

9.

R. J. Bell, Introductory Fourier Transform Spectroscopy, Academic Press(1972). Google Scholar

10.

S. G. Kalenkov, G. S. Kalenkov and A. E. Shtanko, “Self-reference hyperspectral holographic microscopy,” J. Opt. Soc. Am. A, 36 A34 (2019). https://doi.org/10.1364/JOSAA.36.000A34 JOAOD6 0740-3232 Google Scholar

11.

I. Shevkunov, V. Katkovnik and K. Egiazarian, “Lensless hyperspectral phase imaging in a self-reference setup based on Fourier transform spectroscopy and noise suppression,” Opt. Express, 28 (12), 17944 (2020). https://doi.org/10.1364/OE.393009 OPEXFF 1094-4087 Google Scholar

12.

P. Bao et al., “Lensless phase microscopy using phase retrieval with multiple illumination wavelengths,” Appl. Opt., 51 (22), 5486 –5494 (2012). https://doi.org/10.1364/AO.51.005486 APOPAI 0003-6935 Google Scholar

13.

V. Y. Katkovnik, I. A. Shevkunov and K. Eguiazarian, “Hyperspectral phase retrieval,” Proc. SPIE, 11351 113511G (2020). https://doi.org/10.1117/12.2553713 PSISDG 0277-786X Google Scholar

14.

B. Roig-Solvas, L. Makowski and D. H. Brooks, “A proximal operator for multispectral phase retrieval problems,” SIAM J. Optim., 29 (4), 2594 –2607 (2019). https://doi.org/10.1137/18M120227X SJOPE8 1095-7189 Google Scholar

15.

R. Bell, Introductory Fourier Transform Spectroscopy, Elsevier(2012). Google Scholar

16.

V. Katkovnik, I. Shevkunov and K. Egiazarian, “Hyperspectral holography and spectroscopy: computational features of inverse discrete cosine transform,” (2019). Google Scholar

17.

F. Soulez et al., “Proximity operators for phase retrieval,” Appl. Opt., 55 (26), 7412 –7421 (2016). https://doi.org/10.1364/AO.55.007412 APOPAI 0003-6935 Google Scholar

18.

N. Parikh et al., “Proximal algorithms,” Found. Trends® Optim., 1 (3), 127 –239 (2014). https://doi.org/10.1561/2400000003 Google Scholar

19.

J. W. Goodman, Introduction to Fourier Optics, Roberts and Company Publishers(2005). Google Scholar

20.

S. V. Venkatakrishnan, C. A. Bouman and B. Wohlberg, “Plug-and-play priors for model based reconstruction,” in IEEE Global Conf. Signal Inf. Process., 945 –948 (2013). https://doi.org/10.1109/GlobalSIP.2013.6737048 Google Scholar

21.

C. A. Metzler, A. Maleki and R. G. Baraniuk, “Bm3d-prgamp: compressive phase retrieval based on bm3d denoising,” in IEEE Int. Conf. Image Process. (ICIP), 2504 –2508 (2016). https://doi.org/10.1109/ICIP.2016.7532810 Google Scholar

22.

V. Katkovnik, M. Ponomarenko and K. Egiazarian, “Sparse approximations in complex domain based on BM3D modeling,” Signal Process., 141 96 –108 (2017). https://doi.org/10.1016/j.sigpro.2017.05.032 Google Scholar

23.

V. Katkovnik and K. Egiazarian, “Sparse phase imaging based on complex domain nonlocal BM3D techniques,” Digital Signal Process.: Rev. J., 63 72 –85 (2017). https://doi.org/10.1016/j.dsp.2017.01.002 Google Scholar

24.

J. M. Bioucas-Dias and J. M. Nascimento, “Hyperspectral subspace identification,” IEEE Trans. Geosci. Remote Sens., 46 (8), 2435 –2445 (2008). https://doi.org/10.1109/TGRS.2008.918089 IGRSD2 0196-2892 Google Scholar

25.

L. Zhuang and J. M. Bioucas-Dias, “Fast hyperspectral image denoising and inpainting based on low-rank and sparse representations,” IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens., 11 730 –742 (2018). https://doi.org/10.1109/JSTARS.2018.2796570 Google Scholar

26.

V. Katkovnik et al., “Complex-domain joint broadband hyperspectral image denoising,” Sens. Transducers, 233 (5), 33 –39 (2019). Google Scholar

27.

I. Shevkunov et al., “Hyperspectral phase imaging based on denoising in complex-valued eigensubspace,” Opt. Lasers Eng., 127 105973 (2020). https://doi.org/10.1016/j.optlaseng.2019.105973 Google Scholar

28.

J. Bioucas-Dias and G. Valadão, “Multifrequency absolute phase estimation via graph cuts,” in Eur. Signal Process. Conf., 1389 –1393 (2009). Google Scholar

29.

Y. Park, C. Depeursinge and G. Popescu, “Quantitative phase imaging in biomedicine,” Nat. Photonics, 12 (10), 578 –589 (2018). https://doi.org/10.1038/s41566-018-0253-x NPAHBY 1749-4885 Google Scholar

30.

M. Born and E. Wolf, Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light, Elsevier(2013). Google Scholar

31.

P. Hartmann, Optical Glass, SPIE Press, Bellingham, Washington (2014). Google Scholar

Biography

Vladimir Katkovnik received his PhD and DSc degrees in technical cybernetics from Leningrad Polytechnic Institute (LPI) in 1964 and 1974, respectively. From 1964 to 1991, he was an associate professor and then a professor in the Department of Mechanics and Control Processes at LPI. Since 2003, he has been in the Department of Signal Processing, Tampere University of Technology, Finland. He published 6 books and more than 350 refereed journal and conference papers. His research interests include stochastic image/signal processing, nonparametric estimation, computational imaging, and computational phase imaging.

Igor Shevkunov received his PhD in optics from St. Petersburg State University, Saint Petersburg, Russia, in 2013. He has been a postdoctoral researcher at Tampere University since 2017. He is the author of more than 45 refereed papers. His main research interests are digital holography, phase retrieval, and interferometry.

Karen Egiazarian received his MSc degree from Yerevan State University in 1981, his PhD from Moscow State University, Moscow, Russia, in 1986, and his DTech degree from Tampere University of Technology (TUT), Tampere, Finland, in 1994. He is a professor leading the Computational Imaging Group of the ICT Faculty at Tampere University. He has authored about 650 refereed journal and conference papers. His research interests include computational imaging, sparse coding, and image and video restoration. He serves as an associate editor for the IEEE Transactions on Image Processing and was the editor-in-chief of the Journal of Electronic Imaging from 2016-2020.

CC BY: © The Authors. Published by SPIE under a Creative Commons Attribution 4.0 Unported License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.

Citation Download Citation

Vladimir Y. Katkovnik, Igor A. Shevkunov, and Karen O. Eguiazarian "Hyperspectral phase retrieval: spectral–spatial data processing with sparsity-based complex domain cube filter," Optical Engineering 60(1), 013108 (29 January 2021). https://doi.org/10.1117/1.OE.60.1.013108

Received: 7 May 2020; Accepted: 11 January 2021; Published: 29 January 2021

Access the abstract

JOURNAL ARTICLE
14 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY

GET CITATION

KEYWORDS

Phase retrieval

Algorithm development

Data processing

Sensors

Reconstruction algorithms

Refractive index

Optical engineering

1.

Introduction

2.

Problem Formulation

Eq. (1)

Eq. (2)

Eq. (3)

Eq. (4)

3.

Algorithm Development

3.1.

Approach

Eq. (5)

Eq. (6)

Eq. (7)

3.2.

Minimization on Bl

Eq. (8)

Eq. (9)

Eq. (10)

Eq. (11)

Eq. (12)

Eq. (13)

3.3.

Regularization by Sparsity-Based Filters

Eq. (14)

Eq. (15)

3.4.

HS Phase Retrieval Algorithm

Table 1

3.5.

On Interpretation of Hyperspectral Complex Domain Data

Eq. (16)

Eq. (17)

Eq. (18)

4.

Results

4.1.

Simulation Tests

Fig. 1

Eq. (19)

Eq. (20)

Fig. 2

Fig. 3

Fig. 4

4.2.

Physical Tests

Fig. 5

Fig. 6

Fig. 7

5.

Conclusion

6.

Appendix A: Minimization on Bl

Eq. (21)

Acknowledgments

References

Biography

Show All Keywords

Keywords/Phrases

Search In:

Publication Years

Minimization on $B_{l}$

Appendix A: Minimization on $B_{l}$