Phaseless computational imaging with a radiating metasurface

Thomas Fromenteze; Xiaojun Liu; Michael Boyarsky; Jonah Gollub; David R. Smith

doi:10.1364/OE.24.016760

1. Introduction

Remote sensing exploits the reflection of radiated waves to localize, image, and non-destructively detect objects under study. In imaging applications, measured waveforms are usually back-propagated to the object space using techniques such as Kirchoff migration[1], and Stolt’s F-K migration [2], giving the location and reflectivity of the target of interest. This process requires the measurement of fast-time varying signals, represented by complex-valued harmonics in the Fourier domain. These phase measurements become challenging in the terahertz, visible, and X-ray ranges, and require the implementation of complex setups, often based on interferometric techniques. Thus, to overcome these hardware limitations, many phase retrieval techniques have been proposed in the last decade to allow for the reconstruction of complex field distributions solely from intensity measurements. They take inspiration from the concept of holograms defined by Gabor [3], applying alternating projection algorithms introduced by Gerchberg and Saxton [4] and Fienup [5, 6], where at each iteration a complex field source is tailored such that the absolute value of its projection in the target space matches the measured intensity. In the last decade, a resurgence of attention in the scientific community has focused on the development of efficient phase retrieval algorithms. Independent and almost simultaneous contributions by Papanicolaou et al. [7] and Candès et al. [8] proposed comparable optimization techniques, both based on semi-definite programming of a relaxed problem. The practical implementation proposed by Candès et al. is presented here as the foundation of the approach proposed in this article. This technique focuses on the reconstruction of three-dimensional objects from the measured intensity of coded diffraction patterns (Fig. 1). A reconfigurable phase plate encodes the fields diffracted from a molecule illuminated by a coherent beam, allowing the reconstruction of a 3D electron density map from multiple intensity measurements of modulated projections. The depicted X-ray imaging setup can be considered under the emerging framework of coherent computational imaging achieved on the physical layer [9, 10, 11]. The aim of this paper is thus to highligh! these similarities and demonstrate the reconstruction of near-field 3D images from phaseless measurements taken with a computational system. To this end, the mathematical derivations associated with the phaseless measurement of coded diffraction patterns are presented, followed by a concise review of recent phase retrieval techniques.

Fig. 1 The coded diffraction pattern measured from the illumination of a molecule by an X-ray source.

Download Full Size | PDF

The diffracted patterns measured on a plane made of m pixels y ∈ ℝ^m are expressed as follows:

y = {| Ax |}^{2}

where x ∈ ℂⁿ is an unknown object and A ∈ ℂ^m×n is a known transfer matrix in which each line a_i, i = 1,...,m stands for a coded diffraction. The wave diffracted from the object is thus coded by a random mask, giving an illumination pattern y^(l) of the form [12]:

y^{(l)} = {| {FD}^{(l)} x |}^{2}

where F ∈ ℂ^n,n is a discrete Fourier transform matrix and D^(l) ∈ ℂ^n,n is a diagonal matrix whose entries are the known random complex transmittance of the mask modulating the diffracted pattern. L random masks are used for the estimation of x, leading to m = nL measured points to match the initial formulation of Eq. 1. An optimization problem is formulated to find the rank-1 matrix X = xx^†. Indeed, for each line y_k:

y_{i} = {| 〈 a_{i}, x 〉 |}^{2} = Tr (x^{†} a_{i} a_{i}^{†} x) = Tr (a_{i} a_{i}^{†} x x^{†}) = Tr (a_{i} a_{i}^{†} X)

where

a_{i} a_{i}^{†}

is a rank-1 matrix. X being the outer product of two vectors, this matrix must satisfy:

X ≽ 0, rank (X) = 1, y_{i} = Tr (a_{i} a_{i}^{†} X) for i = 1, \dots, m

where X ≽ 0 means that X is positive semidefinite and † stands for the conjugate-transpose operator. Because rank minimization is computationally hard, a relaxation was proposed [8] by dropping the rank constraint and replacing it by a trace minimization, accounting for the sum of the singular values of X. The semidefinite program PhaseLift is hence defined as:

\begin{array}{l} \underset{X}{minimize} & Tr (X) \\ subject to & X ≽ 0, \\ Tr (a_{k} a_{k}^{†} X) = y_{k}, k = 1, \dots, m \end{array}

The convexity of this formulation makes it solvable with the help of optimization software such as CVX [13]. Among the numerous applications of this approach, the excitation retrieval from phaseless far-field data of a microwave array has recently been demonstrated in numerical studies [14, 15]. This represents a promising approach for the simplification of radiation pattern characterization and far field imaging systems. However, the dimensionality “lifting” imposed by this algorithm represents the main drawback, squaring the number of unknowns to create the rank-1 matrix X.

Alternatively, novel alternating projection algorithms represent an interesting approach for solving large quadratic systems. They demonstrate the efficiency of iterative phase retrieval techniques, such as the block-Kaczmarz method [16], derived from the algebraic reconstruction technique [17], and the Wirtinger flow [18]. Recently, the truncated Wirtinger flow was proposed by Chen and Candès [12]. It has been demonstrated that these methods always converge to a solution when satisfying support constraints in the spatial domain and appropriate frequency oversampling, in contrast to the most popular methods introduced by Gerchberg, Saxton, and Fienup that can stall in local minima [8]. The truncated Wirtinger flow is the solution adopted in this article for its demonstrated efficiency, although the other mentioned methods remain compatible with quadratic formulations equivalent to Eq. 1. This algorithm computes the following equations for each iteration:

x^{(t + 1)} = x^{(t)} + \frac{μ_{t}}{m} \sum_{i \in S_{t + 1}}^{m} \nabla_{l_{i}} (x^{(t)})

\nabla_{l_{i}} (x^{(t)}) = 2 \frac{y_{i} - {| a_{i}^{*} x^{(t)} |}^{2}}{x^{(t) *} a_{i}}

For each iteration t, the value of x is thus updated by this descent gradient-like algorithm where μ_t is a step size that can be determined for example by a backtraining line search and ∇_{l_i} is a descent direction. The algorithm is computed on the adaptive index set S_t+1 as determined by Chen and Candès [19], satisfying for any i ∈ S_t:

a_{i}^{*} x^{(t)} \approx ‖ x^{(t)} ‖

\frac{y_{i} - {| a_{i}^{*} x^{(t)} |}^{2}}{a_{i}^{*} x^{(t)}} ≲ \frac{\frac{1}{m} | y_{i} - {| a_{i}^{*} x^{(t)} |}^{2} |}{‖ x^{(t)} ‖}

These constraints improved the efficiency of the initial truncated Wirtinger flow by dropping some gradient components ∇_{l_i}(x) bearing too much influence on the search direction. The efficiency of this algorithm is demonstrated in [19] considering both ideal and noisy measurements—showing that the computational effort required for solving a random quadratic system is on the same order of magnitude of that of a linear system of the same size.

2. Phase retrieval adapted to near field imaging using a radiating metasurface

We propose the adaptation of the phase retrieval framework to a microwave computational imaging system by replacing the coding phase plate presented in the depicted X-ray system (Fig. 1) by a radiating metasurface and implementing the truncated Wirtinger flow instead of the PhaseLift presented earlier. The codes made openly available by their authors, Chen and Candès, are adapted to this study. Several studies have demonstrated the possible simplification of imaging hardware exploiting the recent development of spatially resolving antennas [20]. In this context, systems radiating structured illumination patterns have been proposed, exploiting the inherent structure in the imaged targets to limit the hardware complexity [9, 10, 21]. To this end, frequency diverse [11, 22] and dynamically reconfigurable radiating apertures [23, 24] were studied, demonstrating the linear dependency between the target space and the measured signals through tailored sensing matrices. In this paper, the principle of computational imaging is thus adapted to the phase retrieval problem using a frequency diverse radiating structure, taking benefit of the hardware simplification allowed by both approaches to propose a new paradigm for imaging systems. This demonstration is proposed in the microwave range as a proof of concept to easily compare the reconstructions from complex valued and intensity measurements, paving the way for millimeter wave, terahertz, and photonic applications.

Here, we study the mathematical formalism and develop the conditions in which intensity measurement of compressed waveforms can be applied to retrieve the positions of targets in the near field. The experimental setup and the associated parameters are defined in Fig. (2).

Fig. 2 Computational imaging system used for the experimental demonstration. A metasurface radiating frequency diverse patterns is applied to the localization of field sources from the intensity measurement of a compressed frequency domain waveform.

Download Full Size | PDF

In this setup, a frequency diverse structure similar to those introduced in [9, 22, 25, 26] is considered. The large modal diversity excited by these metasurfaces allows for the radiation of structured field patterns which vary with the driving frequency. The quality factor of the metasurface is optimized to avoid modal degeneration, allowing for the radiation of a large number of pseudo-orthogonal spatial modes sensing the target space. The expression of the measured frequency signal ρ(ν) on the radiating device’s output port can be expressed as:

ρ (ν) = \int_{r_{r}} \int_{r} ϕ (r_{r}, ν) g (r_{r}, r, ν) f (r) d^{3} r d^{2} r_{r}

where ϕ(r_r, ν) stands for the near field distribution of the metasurface measured at the aperture locations r_r for each frequency ν, g(r_r, r, ν) represents the Green function modeling the propagation of field from the object space r and the metasurface’s aperture, and f(r) corresponds to a field source that is localized with this computational system. This problem can be expressed using a matrix formalism by discretizing equation Eq. 10; we represented the resulting vectors and matrices in bold notation as r = [r_i]_1≤i≤n, r_r = [r_{r_j}]_{1≤j≤n_r}, and ν = [ν_k]_1≤k≤m. A sensing matrix H ∈ ℂ^m×n is defined by the product of the Green function matrix G_{n,n_r} (ν_k) and the cavity near-field response written in the vector form ϕ_{n_r} (ν_k) for each frequency ν_k:

H_{n} (ν_{k}) = G_{n, n_{r}} (ν_{k}) ϕ_{n_{r}} (ν_{k})

The sensing matrix allows for a representation of the linear dependency between the measured frequency signal ρ ∈ ℂ^m and the object f ∈ ℂⁿ, leading to the following formulation:

ρ = Hf

Previous works demonstrated methods of reconstructing the image vector f under certain invertibility conditions of sensing matrix H, corresponding to a sufficient number of radiated orthogonal patterns interrogating the target space [9]. This work extends the frame of computational imaging by considering the measurement of the intensity of the compressed waveform ρ, described by the following quadratic equation:

{| ρ |}^{2} = {| Hf |}^{2}

Similarly to the coded diffraction problem, an object is reconstructed from phaseless data knowing the complex transfer function of the coding system. However, the frequency diversity exploited by this approach coupled to the associated derivation makes it distinguishable for its compatibility to near-field imaging, without using any actively reconfigurable radiating components. A simulation of the system depicted in Fig. 2 is studied to identify the number of independent modes required to ensure the reconstruction of the object under study.

In the numerical model of the frequency dependent and randomized radiation pattern of the metasurface, the radiating plane is set at y = 0, where y is the propagation axis. In accordance with (10), the radiation f(r) of a field source (set in the Fresnel zone of the radiating metasurface), is propagated to the aperture plane r_r, compressed by its near-field responses ϕ(r_r, ν) into a unique measurement at the port where the phaseless measurement is performed. Careful consideration of the dispersive nature of the metasurface ϕ(r_r, ν) is required to study the convergence of the computational phase retrieval problem. In the numerical simulations, the impulse response is defined in the radiating plane by a Gaussian random distribution ℛ_d(r_r, t) with mean value zero, variance one, and discretized over steps of c/2ν_max in space and 1/B in time. The quality factor Q of the metasurface is included in this model to consider the cavity’s intrinsic losses and the coupling with all of the irises, leading to a modal degeneration [20]. The random field distribution is thus multiplied by a decaying envelope or damping time τ = Q/πν₀, ν₀ being the central frequency of the studied bandwidth (Fig. 3) [20, 21].

Fig. 3 Impact of quality factor Q on the damping time of the metasurface random responses.

Download Full Size | PDF

Then, computing the Fourier transform from the time domain to the discrete set of frequency components ν, the full model of the metasurface spatial and frequency response is:

ϕ (r_{r}, ν) = 𝔉 [ℛ_{d} (r_{r}, t) exp (- t π ν_{0} / Q)]

In coherent computational imaging system, where complex-valued signals are measured, dispersive antennas presenting long lasting responses are designed to limit the correlation between each radiation pattern in the frequency domain, improving the conditioning of the sensing matrix H. In this way, an estimation of the target signature f̂ can be computed following f̂ = H⁺ρ, where the superscript + stands for the Moore-Penrose pseudo-inverse operator [27]. An ideal metasurface would presented perfectly orthogonal radiation patterns ensuring that H⁺H = I. In practice, metasurfaces are designed to have low correlation among patterns exploiting the frequency diversity, leading to a non-ideal inversion of the sensing matrix and the apparition of parasitic sidelobes. If we are limited to the measurement of intensity |ρ(ν)|², such a direct approach can not be applied. But, alternating projection algorithms may be adapted to this problem to determine of the impact of a radiating metasurface’s characteristics.

The spatial domain sampling is determined by the dimensions of the metasurface, modeled by an array of frequency dispersive dipoles with responses ϕ(r_r, ν). According to the Rayleigh resolution limit, the transverse spatial sampling for a radiating array of dimensions D_x and D_z at a distance R is d_x = λ_cR/D_x and d_z = λ_cR/D_z. The range sampling for a wideband system is given by the width of the radiated pulses as d_y = c/(2B), where B = ν_max − ν_min is the operating bandwidth.

Assuming the measurement of the intensity of a frequency signal ρ described by (13), a numerical study is carried out to estimate the criteria required for an accurate reconstruction of a spatial field distribution f̂_I, where the subscript _I denotes a reconstruction from an intensity measurement. The performance of such a computational system can be validated against a reconstruction of complex-valued measurements, f̂, with a relative error computed for each simulation as:

ε = min_{θ \in [0, 2 π]} \frac{‖ {\hat{f}}_{I} - \hat{f} e^{j θ} ‖}{‖ \hat{f} ‖}

Because the reconstruction from phaseless techniques is unable to estimate an absolute phase, a rotation of the samples by θ = 〈f̂_I, f̂〉 is performed to align the estimations in the complex plane before the subtraction. According to [16, 18], successful estimations are considered when ε ⩽ 10⁻⁵. The field source domain is defined by a number of voxels of positions r. The ratio between the number of measured modes m and the number of reconstructed voxels n, determines the size of the sensing matrix H(r, ν). The ratio m/n is considered in this study for different values of Q to determine the relation between the number of intensity measurements and the size of the reconstructed domain. The simulation is performed on a frequency band defined between ν_min = 17.5 GHz and ν_max = 26.5 GHz (K-band), with a frequency sampling dν = (ν_max − ν_min)/m.

The metasurface delay spread τ and the equivalent quality factor Q are defined according to the frequency sampling as:

τ = \frac{α_{t}}{d ν}

Q = α_{t} π \frac{ν_{0}}{d ν}

where α_t is a frequency sampling parameter set according to the metasurface’s damping factor τ. According to [20], the optimal sampling is α_t = 1/π ≈ 0.32 considering that at most πτB orthogonal channels can coexist due to modal degeneration. The study is presented for α_t ∈ [0.1, 2] to demonstrate the impact of frequency averaging on the phase retrieval algorithm. A domain of n = 6³ = 216 voxels is considered in which a complex random field f must be retrieved. An array of 20 × 20 radiating dipoles spatially spaced at λ_min/2 in both transverse directions is considered to simulate a metasurface whose frequency response is defined by (14) (Fig. 4).

Fig. 4 Numerical simulation of the phaseless computational system applied to the localization of a field source. The radiated field is propagated to the receiving structure plane, coded by the response of this metasurface and compressed into a unique frequency domain signal.

Download Full Size | PDF

100 trials with random metasurface responses and field distributions are computed for each couple of parameters to estimate an empirical set ensuring the accurate estimation of f̂_I (Fig. 5). According to this numerical analysis, the probability of successful recovery tends to reach its maximum when m/n ⩾ 5 and $α_{t} ⩾ \frac{1}{π}$ (Fig. 6).

Fig. 5 Empirical probability of successful recovery. 100 trials are simulated for each pair of parameters m/n and α_t.

Download Full Size | PDF

Fig. 6 Empirical success rate according to the sampling m/n. α_t must be larger than 1/π to ensure that there is at least as many measured modes m as orthogonal channels available in the sensing matrix H to reconstruct n voxels.

Download Full Size | PDF

As predicted in [19], a sampling of a least m = 5n ensures a good agreement between the estimations of the field source f with and without the phase information. Furthermore, the relation between the quality factor of the designed metasurface and the frequency vector ν must satisfy:

α_{t} = \frac{Q d ν}{π ν_{0}} ⩾ \frac{1}{π}

Considering the case of a minimal sampling m = 5n and the definition of the frequency sampling dν = B/m, the maximum number of voxels imaged with this technique takes the following form:

n ⩽ \frac{Q B}{5 ν_{0}}

Having identified the critical designing and sampling parameters allowing for accurate estimation of the compressed, phaseless measurement, the impact of additive noise is studied in the next section to evaluate the performance of the algorithm in the context of near field computational imaging. The model given by Eq. 12 is modified to account for a Gaussian additive noise η:

{| ρ |}^{2} = {| H f + η |}^{2}

With a sampling m = 6n and for different values of signal-to-noise ratio (SNR), 100 trials are computed with random field distributions f and metasurface responses characterized by α_t = 2 to study the convergence of the truncated Wirtinger flow adapted to this near field computational imaging problem (Fig. 7).

Fig. 7 Convergence of the phase retrieval algorithm according to the SNR for m/n = 6 and α_t = 2. 100 trials are simulated for each SNR value.

Download Full Size | PDF

The value of α_t is deliberately chosen to be large in order to speed up the numerical convergence. According to the sampling m/n, the algorithm converges on average in less than 200 iterations. For each value of SNR, the average μ and standard deviation σ_d is computed (Fig. 8). According to the SNR value, the phase retrieval algorithm demonstrates its efficiency by leading to an average relative error μ of $1 / \sqrt{SNR}$ .

Fig. 8 Statistical distribution of the relative errors according to the SNR. In each case, the average μ of the relative error converges to the normalized noise floor $1 / \sqrt{SNR}$ . The standard deviation of each distribution is given by σ_d.

Download Full Size | PDF

Finally, a convergence study is presented considering the impact of α_t acting on the frequency oversampling. A set of numerical simulations are computed with SNR = 10⁶ and sampling m = 6n, considering 1000 trials for each value of α_t (Fig. 9).

Fig. 9 Statistical study of the convergence of the algorithm according to the factor α_t. The results are gathered in the right-hand side graphic, presenting the average μ of each distribution and the standard deviation σ_d. 1000 trials are computed for each value of α_t.

Download Full Size | PDF

This numerical study depicts the relation between the number of iterations required to reach convergence with respect to α_t, demonstrating the positive impact of a high quality factor and a fine frequency sampling on the convergence speed.

3. Pratical implementation and experimental results

The theory of phaseless computational imaging is experimentally validated using a metasurface operating in the microwave regime. To this end, a metallic leaky cavity of 28.5 × 28.5 × 15.2 cm³ was created (Fig. 10), inspired by a computational imaging prototype introduced in [22].

Fig. 10 Radiating metasurface implemented for the validation of the proposed phaseless computational technique.

Download Full Size | PDF

The front plate is perforated by a 15 × 15 cm² square array of circular holes randomly set on a 0.6 cm uniform grid. An open-ended waveguide source is set in front of the cavity and localized using the computational system. In contrast with the setup depicted in Fig. 2, two ports are connected to the back of the cavity (Fig. 11). In this way, different superpositions of modes can be measured by the two ports according to their locations in the cavity, increasing the amount of independent information in the frequency domain. A spherical reflector is also set in the cavity, providing additional mode mixing and increasing the number of uncorrelated states in the cavity [22, 28].

Fig. 11 Radiating metasurface implemented for the validation of the proposed phaseless computational technique.

Download Full Size | PDF

To match the simulations, the operating frequency range is defined in the K-band between ν_min = 17.5 GHz and ν_max = 26.5 GHz, sampled by m_ν = 3601 frequency points. In this way, m = 2m_ν points are measured for the phaseless estimation of the n voxels constituting f. The patterns radiated at each frequency and for each excitation port are measured by a near-field, single-polarized probe moved on a planar synthetic aperture by a translation stage. This field is numerically back-propagated to the metasurface plane to estimate Φ₁(r_r, ν) and Φ₂(r_r, ν), the transfer functions of the cavity when exciting ports 1 and 2, respectively. Examples of near-field distributions are presented in Fig. (12) for two consecutive frequencies of the vector ν. Different pseudo-random spatial fields distributions are thus obtained as a function of the excited port and of the frequency, due to the low level of loss of this cavity.

Fig. 12 Comparison of the near-field distributions Φ₁(r_r, ν) and Φ₂(r_r, ν) measured for the independent excitation of ports 1 and 2. The results are depicted for two consecutive frequency ν₁ = 23 GHz and ν₂ = 23.002 GHz of the frequency vector ν.

Download Full Size | PDF

The measured quality factor of the cavity is about 12000 for both ports, determined by fitting exponential functions to the measured radiation patterns expressed in the time domain using a Fourier transform and taking the root mean square over the spatial dimension r_r. A formulation matching Eq. 1 is proposed according to the measured signals ρ₁ ∈ ℂ^m_ν and ρ₂ ∈ ℂ^m_ν on ports 1 and 2 by concatenating the measured signals and the corresponding sensing matrices computed with Eq. 11:

ρ = [ρ_{1}, ρ_{2}]

H = [H_{1}, H_{2}]

where ρ ∈ ℂ^m and H ∈ ℂ^m×n, dimensioned such that the target field distribution f ∈ ℂⁿ is estimated, with m = 2m_ν. The upper bound of the number of reconstructed voxels can be computed according to the quality factor and the setup presented in this experiment. Merging Eqs. 16 and 17 gives:

τ = \frac{Q}{π ν_{0}}

From the expression of the time constant τ, the maximum number of orthogonal channels is bounded by [20]:

m_{ν} < π τ B

< Q \frac{B}{ν_{0}} \approx 4900

In the considered experiment, there are at most m = 2m_ν = 9800 independent frequency points, measured on the two ports of the cavity. Under the assumption of an ideal estimation of the sensing matrix H and a sampling n = m/5, a maximum of 1960 voxels can be reconstructed. For this validation, since m_ν = 3601 frequency points are measured, a maximum of n = 1440 voxels can be estimated with a sampling of m = 5n. Considering this setup, the value of the parameter α_t is given for one port by Eq. 17:

α_{t} = \frac{Q d ν}{π ν_{0}} = \frac{Q B}{π ν_{0} m_{ν}} \approx 4.3

Under these conditions and considering that two signals are measured for the estimation of f̂_I, a fast convergence of the iterative phase retrieval algorithm is expected.

A first validation is proposed for a scene made of 10 × 10 × 10 = 1000 voxels centered around x = 0, y = 0.4 m, z = 0. In accordance with the numerical study, the spatial sampling is defined as d_x = λ_cR/D_x ≈ 3.6 cm, d_z = λ_cR/D_z ≈ 3.6 cm and d_y = c/(2B) ≈ 1.7 cm, with R = 0.4 m, D_x = D_z = 0.15 cm, and B = 9 GHz. A probe is first set in front of the radiating metasurface in the middle of the pre-defined region of interest. A comparison between the retrieved field distributions f̂ and f̂_I is presented in Fig. 13.

Fig. 13 Localization of a field source on a domain of 10 × 10 × 10 voxels, with and without the phase information. The blue square represents the array of equivalent dipoles constituting the radiating metasurface.

Download Full Size | PDF

The field distribution reconstructed from intensity-only measurements matches the estimation from the complex measurements, taken as a reference in this study. A higher level of noise is observed in the phaseless case, most likely due to a non-ideal estimation of the sensing matrix H considering that a supplementary mounting structure (source of diffraction) was added after the near-field characterization of the leaky cavity leading to a relative error ε = 0.57. A more precise comparison of the two estimated fields f̂ and f̂_I is proposed, extracting the x, y, and z-cuts from the maximum values (Fig. 14).

Fig. 14 Comparison of the x, y, and z-cuts extracted at the maximum value of the reconstructed fields f̂ and f̂_I. The orange solid lines correspond to the phaseless results f̂_I, and are compared to the dashed blue lines standing for the reconstructions from complex measurements f̂.

Download Full Size | PDF

Comparable locations and resolutions are observed in both cases, validating the fidelity of the truncated Wirtinger flow applied to this phaseless computational system. A larger domain is considered for the last part of this experimental validation, studying the impact of the sampling m/n presented earlier in a practical situation. A domain of 20 × 20 × 10 = 4000 voxels is thus considered this time, conserving the same spatial sampling and centered at the same location. We now consider a sampling of m/n = 9800/4000 = 2.45 (with the approximation of two independent measurements ρ₁ and ρ₂). In contrast with the numerical simulations where spatial random field distributions were considered, the experimental cases are focused on the reconstruction of a punctual point like object. Even if a relative error of ε < 10⁻⁵ may not be reachable, we are interested in determining whether a localization of the field source is possible in the given conditions. A comparison of the reconstructions achieved with and without the phase information is once again presented from the same measurements as in the previous case (Fig. 15).

Fig. 15 Localization of a field source set in front of the radiating metasurface in a domain of 20 × 20 × 10 voxels, with and without the phase information. The blue square represents the array of equivalent dipoles constituting the radiating metasurface.

Download Full Size | PDF

A relative error of ε = 0.82 is obtained for this phaseless reconstruction, which is consistent with expectations of a larger ε than in the previous case due to the lower sampling m/n in this case. Despite this larger error, the localization of the field source remains possible, as depicted in the comparison between the field cuts presented in Fig. 16.

Fig. 16 Comparison of the x, y, and z-cuts extracted at the maximum value of the reconstructed fields f̂ and f̂_I. The orange solid lines correspond to the phaseless results f̂_I, and are compared to the dashed blue lines standing for the reconstructions from complex measurements f̂.

Download Full Size | PDF

Taking benefit of this larger domain, the source is translated at an off-center location to be imaged (Fig. 17). With a relative error of ε = 1.03, the three-dimensional reconstructions with and without the phase information remain comparable. The x, y, and z-cuts are extracted from the maximum value of both reconstructions for a finer analysis (Fig. 18).

Fig. 17 Localization of a field source shifted from the center in a domain of 20 × 20 × 10 voxels, with and without the phase information. The blue square represents the array of equivalent dipoles constituting the radiating metasurface.

Download Full Size | PDF

Fig. 18 Comparison of the x, y, and z-cuts extracted at the maximum value of the reconstructed fields f̂ and f̂_I for a source field shifted from the center of the imaging domain. The orange solid lines correspond to the phaseless results f̂_I, and are compared to the dashed blue lines standing for the reconstructions from complex measurements f̂.

Download Full Size | PDF

While the fields extracted from the transverse axis are almost identical, a discrepancy is noted in the range axis. Indeed, the range information is coded by the time of arrival of propagating waves, mainly represented by phase ramps in the frequency domain that cannot be exploited in intensity measurements. Despite the lack of phase data and the under-sampling of the considered case, it remains possible to estimate the location of the field source with a degraded resolution compared to the computational imaging case based on complex valued measurements.

4. Conclusion

An application of a phase retrieval algorithm to a computational imaging system has been presented, allowing for the spatial reconstruction of field distributions from phaseless measurements. The truncated Wirtinger flow has been adapted in this study to determine the position of field sources from the measurements of a metasurface radiating pseudo-orthogonal patterns in the frequency domain. In contrast with the coded X-ray diffraction experiments simulated by Candès et al., there is no need of a reconfigurable random lens since the information is coded in the frequency domain by a static and passive device. While this application has been presented in the microwave range to where it is possible to compare the results with and without the phase information, the most useful applications stand at higher frequencies where the burdensome measurement of phase is problematic and limits the realization of 3D imaging systems. Future studies will thus focus on extending the implementation of such a technique to the terahertz, visible and infrared domains.

Acknowledgments

This work was supported by the Air Force Office of Scientific Research (AFOSR, Grant No. FA9550-12-1-0491).

References and links

1. X. Zhuge, A. G. Yarovoy, T. Savelyev, and L. Ligthart, “Modified kirchhoff migration for uwb mimo array-based radar imaging,” Geoscience and Remote Sensing, IEEE Transactions on 48, 2692–2703 (2010). [CrossRef]

2. J. M. Lopez-Sanchez and J. Fortuny-Guasch, “3-d radar imaging using range migration techniques,” Antennas and Propagation, IEEE Transactions on 48, 728–737 (2000). [CrossRef]

3. D. Gabor, “A new microscopic principle,” Nature 161, 777–778 (1948). [CrossRef] [PubMed]

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

5. J. R. Fienup, “Reconstruction of an object from the modulus of its fourier transform,” Opt. Lett. 3, 27–29 (1978). [CrossRef] [PubMed]

6. J. R. Fienup, “Phase retrieval algorithms: a comparison,” Appl. Opt. 21, 2758–2769 (1982). [CrossRef] [PubMed]

7. A. Chai, M. Moscoso, and G. Papanicolaou, “Array imaging using intensity-only measurements,” Inverse Problems 27, 015005 (2010). [CrossRef]

8. E. J. Candès, Y. C. Eldar, T. Strohmer, and V. Voroninski, “Phase retrieval via matrix completion,” SIAM Review 57, 225–251 (2015). [CrossRef]

9. J. Hunt, T. Driscoll, A. Mrozack, G. Lipworth, M. Reynolds, D. Brady, and D. R. Smith, “Metamaterial apertures for computational imaging,” Science 339, 310–313 (2013). [CrossRef] [PubMed]

10. B. Sun, M. P. Edgar, R. Bowman, L. E. Vittert, S. Welsh, A. Bowman, and M. Padgett, “3d computational imaging with single-pixel detectors,” Science 340, 844–847 (2013). [CrossRef] [PubMed]

11. T. Fromenteze, C. Decroze, and D. Carsenat, “Waveform coding for passive multiplexing: Application to microwave imaging,” Antennas and Propagation, IEEE Transactions on 63, 593–600 (2015). [CrossRef]

12. Y. Chen and E. Candes, “Solving random quadratic systems of equations is nearly as easy as solving linear systems,” Advances in Neural Information Processing Systems pp. 739–747 (2015).

13. M. Grant, S. Boyd, and Y. Ye, “Cvx: Matlab software for disciplined convex programming,” (2008).

14. B. Fuchs and L. Le Coq, “Phase retrieval procedure for microwave linear arrays,” Antennas and Propagation (EuCAP), 2015 9th European Conference on pp. 1–4 (2015).

15. B. Fuchs and L. Le Coq, “Excitation retrieval of microwave linear arrays from phaseless far-field data,” Antennas and Propagation, IEEE Transactions on 63, 748–754 (2015). [CrossRef]

16. K. Wei, “Solving systems of phaseless equations via kaczmarz methods: a proof of concept study,” Inverse Problems 31, 125008 (2015). [CrossRef]

17. S. Kaczmarz, “Angenäherte auflösung von systemen linearer gleichungen,” Bulletin International de l’Académie Polonaise des Sciences et des Lettres 35, 355–357 (1937).

18. E. J. Candès, X. Li, and M. Soltanolkotabi, “Phase retrieval via wirtinger flow: Theory and algorithms,” Information Theory, IEEE Transactions on 61, 1985–2007 (2015). [CrossRef]

19. E. J. Candès and X. Li, “Solving quadratic equations via phaselift when there are about as many equations as unknowns,” Foundations of Computational Mathematics 14, 1017–1026 (2014). [CrossRef]

20. D. L. Marks, J. Gollub, and D. R. Smith, “Spatially resolving antenna arrays using frequency diversity,” J. Opt. Soc. Am. A 33, 899–912 (2016). [CrossRef]

21. T. Fromenteze, E. Kpré, D. Carsenat, C. Decroze, and T. Sakamoto, “Single-shot compressive multiple-inputs multiple-outputs radar imaging using a two-port passive device,” IEEE Access 4, 1050–1060 (2016). [CrossRef]

22. T. Fromenteze, O. Yurduseven, M. F. Imani, J. Gollub, C. Decroze, D. Carsenat, and D. R. Smith, “Computational imaging using a mode-mixing cavity at microwave frequencies,” Applied Physics Letters 106, 194104 (2015). [CrossRef]

23. T. Sleasman, M. F. Imani, J. N. Gollub, and D. R. Smith, “Dynamic metamaterial aperture for microwave imaging,” Applied Physics Letters 107, 204104 (2015). [CrossRef]

24. T. Sleasman, M. Boyarsky, M. F. Imani, J. N. Gollub, and D. R. Smith, “Design considerations for a dynamic metamaterial aperture for computational imaging at microwave frequencies,” JOSA B 33, 1098–1111 (2016). [CrossRef]

25. G. Lipworth, A. Rose, O. Yurduseven, V. R. Gowda, M. F. Imani, H. Odabasi, P. Trofatter, J. Gollub, and D. R. Smith, “Comprehensive simulation platform for a metamaterial imaging system,” Applied optics 54, 9343–9353 (2015). [CrossRef] [PubMed]

26. O. Yurduseven, V. R. Gowda, J. N. Gollub, and D. R. Smith, “Printed aperiodic cavity for computational and microwave imaging,” IEEE Microwave and Wireless Components Letters 26, 367–369 (2016). [CrossRef]

27. T. Fromenteze, E. L. Kpre, C. Decroze, D. Carsenat, O. Yurduseven, M. Imani, J. Gollub, and D. R. Smith, “Unification of compressed imaging techniques in the microwave range and deconvolution strategy,” European Microwave Week 2015-Eurad (2015).

28. G. Montaldo, D. Palacio, M. Tanter, and M. Fink, “Building three-dimensional images using a time-reversal chaotic cavity,” Ultrasonics, Ferroelectrics, and Frequency Control, IEEE Transactions on 52, 1489–1497 (2005). [CrossRef]

Phaseless computational imaging with a radiating metasurface

Abstract

1. Introduction

2. Phase retrieval adapted to near field imaging using a radiating metasurface

3. Pratical implementation and experimental results

4. Conclusion

Acknowledgments

References and links

Cited By

Figures (18)

Equations (26)

Optics Express