Real-Time Imaging With Simultaneous Use of Born and Rytov Approximations in Quantitative Microwave Holography

Microwave and millimeter-wave measurements acquire total-field responses from measurements, yet imaging algorithms instead require the data in the form of scattered-field responses. Two approaches exist for the extraction of the scattered-field data from the total-field responses, namely, the Born and the Rytov data approximations. It is well known that, depending on the target’s size, contrast, and structural complexity, one approximation can achieve an improved accuracy over the other. Yet, the Rytov approximation is rarely used in microwave and millimeter-wave imaging, likely due to phase-unwrapping problems occurring in the case of strongly heterogeneous electrically large targets. Here, we propose a method to utilize the Born and the Rytov approximations simultaneously in a single inversion step for real-time imaging with quantitative microwave holography (QMH). We show through examples with simulated and experimental data that in near-field imaging scenarios, including the imaging of a breast-tissue phantom, there are significant benefits in employing the new method.


I. INTRODUCTION
M ICROWAVE and millimeter-wave imaging methods have been extensively developed over the last halfcentury, with the more recent focus being on close-range imaging, such as concealed weapon detection and security inspection, through-the-wall imaging, nondestructive testing, and biomedical imaging [1]- [11]. The advantages include good penetration through most optically opaque nonmetallic materials, nonionizing radiation, low equipment cost, and spatial resolution on the order of millimeters.
However, many of the objects of interest have significant heterogeneity with wavelength-size structural components, which leads to high propagation complexity and significant forward-model errors corrupting the image quality. Much work is still needed in the development of advanced algorithms and application-specific technologies, especially for near-field and real-time imaging. The mathematical models of electromagnetic (EM) scattering are critically important for accurate image reconstruction. They utilize a variety of approximations in their kernels, i.e., the total electric field and/or Green's function, to reduce the computational complexity [9]. The linear models of EM scattering employ linearizing (weak scattering) approximations of the inherently nonlinear kernel, leading to real-time reconstruction speeds. Yet, image fidelity may suffer, depending on the adequacy of the approximation for the given object under test (OUT). Nonlinear models, on the other hand, enable rigorous inverse-problem solutions through alternating computations of the dielectric profile of the target and the total internal field [10]- [13]. The computational effort is significant, but the image fidelity and spatial resolution are often improved [14]. A common requirement in both the linear and nonlinear reconstruction methods is the input of the scatteredfield responses, which are extracted from the measured totalfield responses. Errors in this extraction are as detrimental to the image quality as measurement noise and uncertainties or errors in the total-field and Green's function distributions.
This has motivated significant research effort in comparing the Born and Rytov approximations as strategies for scattered data extraction [15]- [17]. To the authors' knowledge, all prior work implements both strategies separately and compares their performance. An alternative approach is offered in [18], where a unifying mathematical model is tuned through a single ad hoc parameter to operate as Born's approximation, Rytov's approximation, or in-between the two. Neither of these approaches addresses the main problem in using Rytov's approximation with coherent measurements, namely, the possible failure in unwrapping the phase of the total-field responses across the acquisition aperture [15]. Note that phase unwrapping is a prerequisite for successful image reconstruction with Rytov's approximation. Such failure may occur when the OUT is strongly heterogeneous and/or the signal path through it exceeds several wavelengths. Strong coupling and multiple scattering between the antennas and the OUT in near-field imaging may also lead to phase-unwrapping failure. This failure is the likely reason why Rytov's data approximation remains largely unused in microwave and millimeter-wave imaging with coherent data. This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ Here, we propose a method to utilize Rytov's approximation of the OUT data by completely circumventing its phase unwrapping. At the same time, the method preserves the magnitude and phase information of the kernel of the scattering model; thus, it is still capable of reconstructing the OUT complex permittivity. Furthermore, to make the best use of the complementarity in the advantages of the Born and Rytov approximations, we propose an inversion algorithm that uses both of them simultaneously. The proposed method also balances the limitations associated with each of these approximations, leading to the suppression of image artifacts. We show that the resulting image fidelity is equivalent to or better than those of the two separate algorithms (Bornand Rytov-based). The proposed method is implemented in the framework of quantitative microwave holography (QMH), which has been developed recently for the purpose of realtime quantitative imaging [17]. QMH is particularly suitable for near-field measurement scenarios since the kernel of its scattering model can be derived from the measured system point-spread function (PSF) [19]. Note that the measured PSF, too, needs scattered-field response extraction.
We note that the proposed approach of combining Born's and Rytov's data approximations is in principle applicable to any image reconstruction method where the scattered-field response is extracted from total-field measurements.
In [20], we have reported the first experimental verification of the benefits of using the Born and Rytov approximations in tandem. Higher quality images of a breast phantom with near-field measurements have been obtained while maintaining the real-time speed of the QMH reconstruction. Here, the theory underlying the combined Born-Rytov approximation approach is derived and the implementation in a straightforward QMH inversion algorithm is detailed. Common phasewrapping errors are discussed, which arises in the near-field imaging of complex structures. Examples based on simulated data highlight the situations where the Born and Rytov approximations have unique advantages over one another. The advantage of combining both in a single inversion algorithm is also demonstrated. An experimental example further demonstrates the benefits of utilizing both approximations in tandem.

II. BACKGROUND: QMH
A brief overview of QMH [9], [17], [21], [22] is provided below to aid understanding of the changes necessary to combine the Born and Rytov scattered-data extraction approaches. QMH is based on a linearized model of scattering but utilizes the measurement of a scattering probe (SP) to produce a system-specific data PSF. The measured PSF is used in the inversion process to accurately account for the field distributions in the imaged volume, providing an advantage over simulation-based or analytical approximations, especially in near-field imaging. It also enables the reconstruction of quantitative images, which depicts the estimates of the permittivity and conductivity of the imaged object. This is a significant advantage over previous linear reconstruction methods, which only produces qualitatitve images of the target reflectivity. The RO captures the incident-field portions of the measured responses, which depends on the antennas as well as the impact of the environment. The position vectors of the Tx and Rx antennas are denoted as r Tx and r Rx , respectively. The CO contains a small SP at r sp of known size and permittivity, the measurement of which produces the data PSF. A position within the OUT is denoted by r .

A. Measurement Setup
Holographic reconstruction (QMH included) can be applied on data collected over planar, cylindrical, or spherical surfaces [23]. Here, all examples are based on a planar scanning arrangement (along x and y), which is illustrated in Fig. 1. The acquisition planes are on both sides of the imaged object allowing for reflection and/or transmission measurements.
Two preliminary measurements are necessary before the imaging of OUTs. The first measurement is that of the reference object (RO) shown in Fig. 1(a). The RO is composed of the background medium along with the antennas and the components of the positioning mechanism. It provides the incident-field responses, which are needed in the extraction of the scattered-field responses. Note that the background medium is not necessarily air. For example, in the case of breast cancer imaging, the tissue can be embedded in a dielectric medium designed to reduce reflections at the interface between the tissue and the background [22].
The second measurement is that of the calibration object (CO) [see Fig. 1(b)]. The CO is identical to the RO with the exception of an electrically small SP embedded in the center of the background medium. This measurement provides the data PSF, which enables quantitative image reconstruction.
The selection of the size and permittivity of the SP requires careful consideration, as it impacts the quantitative accuracy of the reconstruction [19], [24]. The SP size should be less than λ min /4, where λ min is the shortest wavelength. This allows for the assumption that the electric field within the probe is uniform, simplifying the forward model. Yet, the probe cannot be too small since this results in very weak scattering and poor signal-to-noise ratio (SNR) of the acquired PSFs. The best quantitative result is achieved if the permittivity of the probe satisfies [19] sp (ω) ≈ (r ; ω) (1) or where (r ; ω) is the permittivity of the OUT at the location r = (x , y , z ) within the imaged volume V and at the angular frequency ω, sp (ω) is the permittivity of the SP, and b (ω) is the permittivity of the background medium. As per (1), the probe's permittivity should be selected to be close to that of the OUT. In the case of highly heterogeneous OUTs, the SP is selected to match particular structures of interest embedded within. For example, in breast cancer imaging, sp is selected to match closely the permittivity of cancerous tissue [17], [22], [25]. The alternative suggested by (2) requires a background permittivity that is significantly higher than that of the OUT or the SP. This is not practical in real-time image reconstruction since a large contrast between the OUT and the background leads to multiple scattering effects, which are not accounted for in the linearized scattering models. Furthermore, depending on the imaged target, it may prove challenging to find a material with sufficiently higher permittivity.

B. Forward Model of Scattering
For brevity, only a short description of the forward model of scattering employed by the QMH method is shown here. More detailed derivations are available in [9] and [17]. The QMH forward model was originally cast in terms of the scattered electric field and is derived from the vector Helmholtz equation [26], [27]. However, electric field responses do not directly translate into measured responses such as S-parameters. Beaverstone et al. [28] derived the exact forward model of scattering in terms of S-parameters as where S sc jk is the scattered portion of the S-parameter measured with the OUT, in which the receiving (Rx) antenna j is positioned at r Rx = (x Rx , y Rx , z Rx ), the transmitting (Tx) antenna k is at r Tx = (x Tx , y Tx , z Tx ), and ω is the angular frequency. Also, 0 is the permittivity of free space, α ξ , ξ = j, k, are the root-power waves exciting the respective ports, r (r ) is the permittivity contrast within the OUT at position r = (x , y , z ), E inc j is the incident field due to the Rx antenna if it were to operate in a Tx mode excited by α j , and E tot k is the total internal field due to the Tx antenna. The approximation of the total internal field as being the same as the respective incident field, E tot k ≈ E inc k , 1 leads to the linearized S-parameter model of scattering in terms of the data PSF [9], [17], [20] S sc jk (r Rx , r Tx ; ω) ≈ 1 r,sp sp · V r (r )H jk (r Rx , r Tx ; r ; ω)dr (4) where H jk is the scattered-field response of the measurement system due to an electrically small (point-like) scatterer (the SP) of volume sp and contrast r,sp = r,sp − r,b at r . This is exactly the data PSF characterizing the particular measurement system. Mathematically, the data PSF represents the linearized resolvent kernel of (3) since where the subscript r of the term in the square brackets emphasizes the position of the SP. Thus, the linearized scattering model in (4) is in essence the firstorder Born approximation of the scattered-field S-parameter response.
The model in (4) implies that the PSF H jk is needed for a probe at all positions r in the imaged volume. However, in a laterally uniform (layered) isotropic medium, the PSF can be assumed translationally invariant in (x, y) and (4) is written as · H jk r Rx xy − r xy , z Rx ; r Tx xy − r xy , z Tx ; z ; ω dx dy dz (6) where r Rx xy ≡ (x Rx , y Rx ), r Tx xy ≡ (x Tx , y Tx ), r xy ≡ (x , y ), and is the complex reflectivity function to be reconstructed. The PSF H jk (r Rx xy , z Rx ; r Tx xy , z Tx ; z ; ω) is acquired with the SP at the center of the z = const. plane and r xy = 0. Since the scan is 2-D (along x and y with fixed z Rx and z Tx ), coordinate translation cannot be employed along z unless the PSF dependence on the z distances (from the Tx and Rx antennas to the probe) is approximated analytically. Such approximation is limited to far-field behavior in a uniform background and is not considered here. Thus, in order to generate the 3-D images, N z PSFs are needed, where the SP is always centered laterally (r xy = 0) but is shifted along z to each desired range slice. This requires N z CO measurements with the SP at (0, 0, z n ), n = 1, . . . , N z .

C. Inversion With QMH
Consider the application of the forward model (6) to the case of monostatic and/or bistatic scenarios, where the Rx and Tx antennas in each jk pair are in a fixed mutual position, i.e., r Tx and r Rx can be described by a single position vector r relative to the scanned OUT. For example, in a monostatic setup, r Tx = r Rx = r. In the bistatic setup in Fig. 1, the Rx and Tx antennas are aligned along boresight (z-axis) and scan on two opposing sides of the OUT so that r Rx ≡ r and r Tx = r −ẑ ẑ zd, whereẑ ẑ z is the unit vector along z and d is the distance between the antennas. In this case, (6) is written as · H jk (r xy − r xy ,z; z ; ω)dx dy dz (8) where r xy ≡ (x, y) andz = const. define the plane scanned by the Rx antenna. We reiterate that the scan position (r xy ,z) also defines uniquely the position of the Tx antenna since the jk antenna pair is in a mutually fixed configuration during the scan. For simpler notations, hereafter, the subscript jk is replaced with a single subscript ξ (ξ = 1, . . . , N T ), indicating the type of response.
Since (8) is a convolution in x and y, it can be efficiently solved in the 2-D Fourier space (k-space) where the convolution becomes multiplication. After discretization of the integral over z into a sum, the k-space equations are given by where tilde represents the Fourier transformed quantity,  [17]. The systems solved at all points in Fourier space are small since N T , N ω , and N z are on the order of 1-10, 10-100, and 1-10, respectively. Thus, they are solved very quickly. The number of solved systems is N x × N y , where N x and N y are the samples along x and y, respectively. With the 2-D fast Fourier transform (FFT), the number of spatial samples equals that in k-space. The QMH inversion is fast and is typically completed within seconds even in problems where the number of voxels is on the order of 10 5 . The computations are also amenable to parallel computing since each small system is solved independently.
Onceρ(κ κ κ, z m ), m = 1, . . . , N z , is found, its real-space counterpart ρ(x , y , z m ) is recovered through 2-D inverse FFT and the quantitative estimate of the complex permittivity contrast r (r ) is computed using (7). We emphasize that the quantitative output in QMH is achieved through the measured data PSF since it is quantitatively accurate with respect to the permittivity contrast of a scatterer.
At the same time, QMH has the added advantage of quantitative image output and the ability to handle equally well near-and far-field measurements. The QMH computational complexity is dominated by: 1) the 2-D FFT of the OUT data and the PSFs, z ) for a pseudoinverse solver; and 3) the 2-D slice-byslice inverse FFT of the reflectivity, O(N z · N x N y log 2 (N x N y )). On the other hand, microwave holography and the SAR RMA involve 2-D FFT of the OUT data, too, with the same computational complexity as in QMH. To cast the 3-D k-space reflectivity into real space, a 3-D inverse FFT is needed with computational complexity of O(N x N y N z log 2 (N x N y N z )). Also, Stolt's mapping (also known as Stolt's interpolation) adds a computational complexity estimated as O(N i N x N y ), where N i depends on the chosen interpolator and the length of the interpolated sequence [36]. Note that Stolt's interpolation is not needed in QMH.
Another common real-time image reconstruction method is the backprojection algorithm (BPA). In its original implementation, the BPA operates in real (x, y, z) space (see, e.g., [32]) and it has high computational complexity of O(N 2 x N 2 y N z N ω ). However, fast BPAs have been proposed, which achieves the reduced complexity of O(N 5/2 ), or even O(N 2 logN), for a 2-D reconstruction of size N × N [37]. Even the fast BPAs do not outperform the RMA and holography in terms of speed.
Finally, it is important to note that QMH is orders of magnitude faster than the traditional quantitative inversescattering methods, which employs iterative updates (e.g., nonlinear optimization) and simulation-based EM forward models. Such methods can potentially yield very accurate quantitative estimates with improved spatial resolutions, but their computational complexity is prohibitive for real-time imaging due to the use of full-wave simulators [9], [11].

III. EXTRACTING SCATTERED DATA
An important aspect of the scattering model is the approximation used to extract the scattered-field responses in the left-hand side of (4) from the measured totaland incident-field responses. Two common approaches are deduced from the Born and Rytov first-order approximations [9], [15], [18], [26]. Under the first-order Born approximation of the external field, the total measured response is a superposition of the incident-field (or RO) response and the desired scattered-field response, which, in turn, is represented by the linearized scattering model in (4). In itself, (4) is a superposition integral, which neglects higher order scattering effects such as mutual coupling and multiple scattering, leading to an effective representation of the imaged object as a collection of independent (or uncoupled) point scatterers.
The Rytov approximation, on the other hand, represents the total field as a complex phase correction on the incident field. Each approximation has distinct advantages and limitations that lead to unique image reconstructions when performed independently.
The QMH forward model (6) operates on the scatteredfield responses that are extracted from the CO measurements with the SP 2 and from the OUT measurements. The CO and RO (incident-field) measurements provide the input for the extraction of the data PSFs, i.e., H ξ ≡ S sc,CO ξ , whereas the OUT and RO measurements provide the input for the extraction of the OUT data S sc ξ . It is worth noting that the data from monostatic reflection measurements are often assumed to be purely scattered-field responses, which is not true in practice; they do include object-independent effects, such as the reflection at the antenna terminals and the reflections from the imaging-hardware components. Removing all RO signal components, measured in the absence of an OUT, is critical for maximizing the fidelity of the reconstruction.

A. Born's Approximation
Born's first-order approximation of the external-field responses provides a simple method for extracting the scattered-field portion of the S-parameters from the measured incident-and total-field responses where S inc ξ is the incident-field response from the RO measurement [see Fig. 1

(a)] and S tot
ξ is the total-field response from an object's measurement [see Fig. 1(d) and (c)]. Note that the subtraction (10) is applied to the total-field responses of both the OUT, S tot ξ , and the CO, S tot,CO ξ , to obtain the OUT data and the PSFs, respectively. As implied by Born's firstorder approximation of the external field, the so extracted S sc ξ data are represented by the approximate (linearized) scattering model in (4) or (6).
The limitation of the first-order external-field Born approximation in the size and permittivity contrast of the scattering object must be kept in mind [9], [15] 2 a|k(r , ω) − k b (ω)| max < π (11) where k(r , ω) is the wavenumber inside the scatterer, k b is the background wavenumber, and a is the radius of the smallest sphere circumscribing the scatterer. In the case of the SP in the CO, k is constant in its volume. The limitation (11) is satisfied in the case of the CO as long as the probe is chosen properly (see Section II-A). However, it may not always be satisfied by the OUT, which leads to image artifacts in regions where higher order scattering effects exist due to mutual coupling and multiple scattering. Note that in near-field measurements, such effects may exist not only within the OUT but also between the OUT and the antennas.

B. Rytov's Approximation
Rytov's first-order approximation represents the totalfield response as a complex-phase correction on the 2 The selection of the SP for the CO determines the region of validity for which the forward model remains accurate. Ideally, the probe's permittivity is selected to match closely the permittivity of the target. In the case of highly heterogeneous targets, it should be close to the permittivity of structural components that are of particular interest, e.g., the cancerous tissues in breast cancer screening. Details about the selection of the CO permittivity can be found in [19]. incident-field response as [9] S tot where the scattered-field response S sc ξ in the numerator of the complex phase is represented by the linearized scattering model (4). The expression is rearranged to express the extracted data S sc ξ in terms of the measured incident and totalfield responses [17] S sc Unlike Born's approximation, Rytov's approach is limited only by the permittivity contrast of the object and not its size 3 However, Rytov's approximation is far more sensitive to errors arising when the left-hand side of (14) exceeds several percents. Also, (13) is inherently prone to errors due to phase wrapping [20]. Consider the breakdown of the logarithm in (13) into real and imaginary parts Note that the imaginary part is a difference of phases, which are wrapped between −π and π in measurements. Without unwrapping, this difference exhibits sharp discontinuities in the 3-D observation space (x, y, ω). The discontinuities in (x, y) adversely affect all inversion strategies that rely on Fourier transforms (FTs) that cast the data from the (x, y) space into the wavenumber space (k-space). Note that unwrapping along frequency ω does not ensure unwrapping in (x, y).
Unwrapping the incident-field phase S inc ξ in (x, y) is usually successful since the background is laterally uniform. However, unwrapping S tot ξ is often problematic since there is ambiguity whether a particular sharp transition in phase is due to reaching the π limit or if a true rapid phase transition has occurred. Consider Fig. 2(a), which depicts the wrapped phase of an OUT measurement of a breast-tissue phantom at a particular frequency. A rapid phase change from −π to π marks a portion of the left-side edge of the phantom, where the unwrapping is successful, as shown in Fig. 2(b). Meanwhile, other regions, including the phantom's interior, exhibit rapid phase transitions within the 2π limit as well as transitions from −π to π, which may or may not be due to the actual properties of the object. This ambiguity is the fundamental reason for the unsatisfactory unwrapped phase distribution in Fig. 2(c).
The applicability of the Rytov approximation (13) is conditional upon the successful phase unwrapping of the data in  (x, y). Unfortunately, the unwrapping may fail if significant property discontinuity exists within the OUT. This is the likely reason for the limited use of Rytov's approximation in microwave and millimeter-wave imaging with coherent data.

C. Comparing the Born and Rytov Data Approximations
Consider a cube of variable edge length and permittivity when it is placed in a background medium of r,b = 10 − i 2 and is measured at 8 GHz. 4 Depending on the cube's size and permittivity, it may violate the limits of one or both of the approximations. Fig. 3(a) visualizes the limits of the Born and Rytov first-order approximations as calculated using (11) and (14) for this case. It is clear that the range of size and permittivity values where both approximations are of acceptable accuracy (area in dark blue) is limited. In this range, the two approximations yield images that are very similar to one another and the unwrapping of the phase of the OUT data for the Rytov approximation is usually successful. The complementarity of the two approximations is exemplified by the regions in yellow and light blue, where their simultaneous use can be beneficial. Another case is shown in Fig. 3(b) for a background of vacuum. In this case, the validity range of the Rytov approximation is extremely limited around the permittivity of vacuum, making it unsuitable for most dielectric targets.
The boundaries shown in Fig. 3 do not indicate strict cutoff for image reconstruction with the two linearizing approximations. The reconstruction degradation is gradual as the object's size and permittivity go beyond the limits. It manifests as image artifacts in regions where the total internal field differs significantly from its incident counterpart. Degradation in the quantitative accuracy of the reconstruction occurs as well.
The differences between the Born and Rytov data approximations are illustrated with a simple 1-D reconstruction example, where the data are acquired via simulations with the time-domain software MEFiSTo [38]. A uniform plane wave is excited as a pulse, the spectrum of which is between 10 and 25 GHz (at −3-dB level). The simulator is set to operate in a 2-D TM z mode with the electric and magnetic field vectors being along z and y, respectively. The 1-D plane-wave propagation along x is enforced by two parallel magneticwall boundaries running along x. The x-extent of the structure is 410 mm and it is terminated by absorbing boundary conditions at both ends. The excitation is at the left end of the computational domain. Two field probes in a vacuum background define ports 1 and 2, capturing reflection and transmission measurements, respectively. They are placed 20 mm away from the ends of the computational domain, i.e., they are 370 mm apart. They record the incident, scattering probe, and OUT data in three separate simulations.
The SP is a 2-mm-thick dielectric slab of r,sp = 1.1 placed midway between the field probes. The OUT is an electrically large dielectric slab of r1 = 1.05 and a length of 300 mm. This length corresponds to approximately 10 and 25 wavelengths at 10 and 25 GHz, respectively. At this length, Born's approximation is valid up to about 23 GHz, whereas Rytov's approximation is marginally valid across the frequency band.
The time-domain responses are cast in the frequency (ω) domain and the reconstructions are performed using (9), where the scattered data are extracted with either the Born or the Rytov approximation. The results are shown in Fig. 4. The two reconstructions, when using only the reflected signal at port 1, are practically identical, and both of them fail to estimate the permittivity of the slab's interior; only its edges are detected. However, when the transmission responses are added to the reflected ones, the Rytov-based reconstruction exhibits improved accuracy compared to the Born-based one. This result is important. It indicates that the benefits of using Rytov's data approximation are significant in transmission (forward-scattered) measurements. This result reinforces the extensive comparisons of the two approximations in diffraction tomography [15], [18], [39] and in transmission electron microscopy [40], [41], which considers applications with transmission measurements only. In contrast, the majority of the microwave and millimeter-wave imaging radars employ reflection responses (e.g., [2], [4], [42]), where the benefits of Rytov's approximation are insignificant. Yet, in applications such as breast cancer imaging and nondestructive testing, transmission data are common since they provide better data SNR [6], [43]- [47]. Moreover, the background permittivity in these applications is often relatively high, which expands the applicability of Rytov's approximation as already discussed in Fig. 3(a). The mathematical reason for the almost identical reconstructions when using only the reflection coefficients [see Fig. 4(a)] is a scattered response, which is much weaker than the incident one. It is well known that when the ratio of scattered-to-incident response is very small, the first-order Rytov approximation converges toward that of Born [9]. In this example, the incident response is collected by the field probe at port 1 in the absence of a scatterer and it is two orders of magnitude stronger than the scattered one. In practical reflection measurements, where the same antenna is used to transmit and receive, the incident response is the antenna reflection coefficient in the absence of a scatterer. With typical absolute values between 0.1 and 0.3 for a well-matched antenna, this response is still expected to exceed by orders of magnitude the change in the reflection coefficient due to a weak scatterer.

D. Tunable Born/Rytov Approximation
In [18], an approximation is proposed for optical diffraction tomography based on an analytical expression unifying the Born and Rytov approximations. The approximation in (13) is cast in a form that expresses the logarithmic term as a limit S sc ξ (r; ω) ≈ S inc ξ (r; ω) · lim n→∞ n ⎛ ⎝ S tot ξ (r; ω) S inc ξ (r; ω) Notice that with the setting of n = 1, (16) is identical to the Born approximation in (10). Thus, by selecting a value of 1 ≤ n < ∞ in the expression one can select the Born, the Rytov, or an in-between approximation. The study in [18] highlights the key difference between the case of n = 1 and that when n is very large: "the tendency appears to be that the Born approximation reproduces the boundary well but not the interior, while the Rytov fills the interior but removes the boundary." This conclusion motivates efforts to combine the two approximations by striking a balance between boundary and interior object reconstruction. Unfortunately, the use of (17) does not eliminate the need to unwrap S tot ξ since the complex nth root still leads to a phase ambiguity of 2π/n.

IV. COMBINED BORN/RYTOV QMH RECONSTRUCTION
It has already been demonstrated that the two approximations produce somewhat different results in near-field experiments with QMH image reconstruction [17], [20]. It was observed that the Born-based QMH recovers well small highpermittivity inclusions, whereas Rytov-based QMH recovers high-loss inclusions better. This complementarity adds to the discussions in Sections III-C and III-D.
In this work, we propose a method of utilizing the Born and Rytov approximations in tandem in a single reconstruction step where unwrapping the total-field data is not necessary. To this end, we exploit the ability of the QMH reconstruction to "stack" equations of the form (9) into a linear system of equations at each point in k-space. This allows for incorporating the equations resulting from the two approximations in a common solution.
Furthermore, to circumvent the problematic phase unwrapping in Rytov's approximation of the OUT data, only the real part of the logarithm term in (15) is used. To this end, the scattering model (8) is premultiplied as [20] S inc ξ (r xy ,z; ω) S inc ξ (r xy ,z; ω) The substitution of S sc ξ in the left-hand side of (18) with Rytov's approximation (13) results in a complex-valued expression, the real part of which is where (·) stands for (r xy ,z; ω). This represents the real part of the modified Rytov-approximated OUT data. Since the factor |S inc (18)  . (20) Note that the PSF H RA ξ is complex and it requires the unwrapping of S tot,CO ξ in (x, y). However, the probe is a weak scatterer that imparts minimal phase shift, which is reasonably handled by any 2-D phase-unwrapping algorithm [48]- [50].
Also, the multiplicative factor |S inc ξ (r; ω)|, which is present in the data and PSF terms on both sides of the scattering model, is not redundant. It preserves the magnitude of the modified Rytov-approximated OUT scattered data, which is important in the subsequent integration of the Rytov-and Born-based scattering equations.
In order to construct a system of equations containing both the Born model and the real part of the modified Rytov model (18) S sc,BA ξ,Im r xy ,z; ω ≈ V ρ Re r · H BA ξ,Im r xy − r xy ,z; z ; ω + ρ Im r · H BA ξ,Re r xy − r xy ,z; z ; ω dx dy dz (22) where the subscripts Re and Im denote the respective components of a complex quantity and superscripts BA and RA indicate the data approximation. Similar to the QMH procedure in Section II-C, a 2-D FFT is applied to all three real-valued datasets (S sc,BA ξ,Re , S sc,BA ξ,Im , S sc,RA ξ,Re ) and the four real-valued PSFs (H BA ξ,Re , H BA ξ,Im , H RA ξ,Re , H RA ξ,Im ). This allows for casting (21) and (22) Re κ κ κ, z m · H BA ξ,0,Im κ κ κ, z m ; ω n +ρ Im κ κ κ, z m · H BA ξ,0,Re κ κ κ, z m ; ω n (24) where S sc,BA ξ,Re (κ κ κ; ω n ) and S sc,BA ξ,Im (κ κ κ; ω n ) are the 2-D FTs of the real and imaginary parts of the Born-approximated OUT data S sc,BA ξ (r xy ,z; ω n ), respectively. For brevity, the position of the acquisition planez has been omitted in the argument of the Fourier-transformed data. Analogously, H BA ξ,Re (κ κ κ, z m ; ω n ) and H BA ξ,Im (κ κ κ, z m ; ω n ) are the 2-D FTs of the real and imaginary parts of the PSF H BA ξ (r xy ,z; z ; ω n ), respectively. The FTs of the real and imaginary parts of ρ(r xy , z m ) are denoted bỹ ρ Re (κ κ κ, z m ) andρ Im (κ κ κ, z m ), respectively.

V. VALIDATION THROUGH SIMULATION
In order to demonstrate the effectiveness of the proposed reconstruction method, the measurement of an OUT composed of several objects is emulated through simulations in FEKO [52] (see Fig. 5). The objects are embedded in a uniform background medium with r,b = 10 − i 0. The first object is a large prism of dimensions 4 cm × 4 cm × 3 cm and relative permittivity r = 25 − i 5. The small cubes (1 cm on a side) are also present. Their relative permittivity is r = 40 − i 8. One of these cubes is embedded in the large prism offset from the center, and the other resides in the background medium.
The measurement is performed by six half-wavelength (at each frequency) dipole antennas. Only the two dipoles residing on the z-axis on both sides of the OUT transmit (one at a time), whereas all six antennas receive, leading to a total of 12 datasets. The antennas perform a raster scan along x and y at 5-mm increments across a 15 cm × 15 cm aperture. The frequency is swept from 3 to 7 GHz in 1-GHz intervals. The selection of the frequency bandwidth and interval is based on the desired range resolution and the maximum range associated with the imaged volume (which, in turn, ensures no aliasing along frequency) [9]. The antennas are located either at z = 60 mm or z = −60 mm.
The data PSFs are acquired in three separate simulations with the same background and antenna arrangement as in the OUT simulation, where the scattering object is a cuboidal SP of edge length 1 cm ( r,sp = 40−i 8) at the center, x = y = 0, of each of three planes z = 50, 60, and 70 mm. These are the planes where the image slices are generated.
The 3-D reconstructions of the complex permittivity using the Born, the modified Rytov, the Marks, and the proposed combined Born/Rytov QMH reconstructions are shown in Fig. 6. The Born-and modified Rytov-based reconstructions do not require the phase unwrapping of the OUT data along (x, y). We reiterate that the modified Rytov reconstruction, based on (21), employs only the absolute values of the OUT data whereas the respective PSFs are complex. Since the proposed combined Born/Rytov approach, too, utilizes (21), it does not require OUT data phase unwrapping either. Both the modified Rytov and the combined Born/Rytov reconstructions need the phase unwrapping of the CO responses in order to extract the complex Rytov-based PSFs. This phase unwrapping is performed with Itoh's method [48] and it is successful. On the other hand, the phase unwrapping of the OUT data is necessary for the reconstruction based on the Marks approximation. In this example, the OUT data unwrapping is successful.
The Born reconstruction [ Fig. 6(a) and (b)] outlines well the structures of the large prism and the external small cube. It also correctly shows the increased contrast of the external small cube relative to the larger prism. However, it fails to clearly identify the embedded cube within the prism at z = 60 mm and to estimate its higher contrast relative to that of the large prism. The reconstruction also incorrectly reconstructs a very low permittivity in the z = 60 mm slice of the large prism. The modified Rytov approximation, on the other hand, produces less structurally accurate representations of the three objects, especially that of the large prism [see Fig. 6(c) and (d)]. However, the presence of a strong scatterer in the large prism is apparent in the slice z = 60 mm. The small external cube is also reconstructed nearly identically to that of the Born-based reconstruction. Rippling artifacts are observed, which are often encountered in k-space reconstruction, especially when Rytov-based data approximations are used. These are Gibbs' effects caused by sharp transitions either in real (x, y) or in Fourier (k x , k y ) space [53].
The reconstruction using Marks' tunable approach (with n = 2) provides an in-between reconstruction between the two approximations, with slightly enhanced detection of the  I   ROOT-MEAN-SQUARE ERROR OF REAL AND IMAGINARY RELATIVE  PERMITTIVITIES RECONSTRUCTED IN THE  SIMULATION-BASED EXAMPLE embedded cube at z = 60 mm [see Fig. 6(e) and (f)]. Also, the small ad hoc value of n is low and the possibility of errors due to the 2π/n ambiguity is low. The combined Born/Rytov approach produces an image that retains the structural accuracy of the Born approximation while also detecting the embedded cube in the larger prism [see Fig. 6(g) and (h)]. The rippling artifacts are suppressed and the reconstruction is improved relative to the individual Born and modifiedRytov reconstructions.
The root-mean-square error of the real and imaginary parts of the permittivity for each reconstruction is shown in Table I. Note that while the Rytov approach has an overall lower error in the real part of the permittivity relative to the other methods, the error in the imaginary part is significantly larger than the other methods. On average, the Marks approach and the combined Born/Rytov approach show improvement versus using solely the Born or Rytov approximations. However, we must reiterate that the performance of Marks' approach is subject to the successful phase unwrapping of the OUT data, whereas the proposed combined Born/Rytov approach does not require this phase-unwrapping step, which ensures wider and more robust applicability.
Note that all the image reconstructions generate some nonphysical values, particularly in the imaginary part of the permittivity. This effect can be expected when reconstructing images near the permittivity bounds ( r = 1 and r = 0) and are due to the Gibbs effects as well as the limitation of the employed linearized scattering model.

A. Measurement Setup
The example discussed here is based on the transmission S-parameter measurements used previously in [20] and [25]. A compressed breast phantom is measured to evaluate the ability of QMH to detect breast-tumor simulants with planarscan data. The phantom is constructed from five 1.1-cm carbon rubber sheets, with permittivity values shown in Table II (see also Fig. 7). These permittivity values are selected to match the average permittivity of BI-RADS Type 2 breast tissue, which has scattered fibroglandular content of less than 50% of the overall tissue mass [54]. Inside two of the five sheets, sections are hollowed out and cancerous tissue simulants [circled in blue in Fig. 7(a) and (b)] as well as fibroglandular tissue simulants [white structure in Fig. 7(b)] are inserted. The remaining space is filled with a matching medium whose permittivity (see Table II) is close to that of the carbon rubber sheets. It is important to remove the air pockets as  they present high-permittivity contrast that may lead to image artifacts. All five carbon rubber sheets are stacked to form the completed phantom [see Fig. 7(c)], in which Layers 2 and 4 contain the inclusions described above. Layers 1, 3, and 5 are homogeneous. Plastic wrap is used to secure the inclusions in Layers 2 and 4 as well as to hold together the phantom. The acquisition is performed in an RF-shielded chamber [see Fig. 8(a)] with 1 TEM-horn Tx antenna [55] and a nineelement Rx array of bowtie antennas [56] [see Fig. 8(b)]. The antennas are designed to deliver maximum power when they are in direct contact with or are very close to (within 2 mm) the breast-tissue phantom. Of the nine available bowtie antennas, only the five copolarized antennas are used; the remaining are connected to 50-loads. This is because the data SNR of the cross-polarized signals is poor as determined through the quality control evaluation of the experimental setup [24]. The chamber contains a platform, which is moved laterally by stepper motors. In the current experiment, the spatial sampling step is 3 mm.
The phantom is placed in a Plexiglas tray and surrounded by an embedding medium, the permittivity of which is listed in Table II [see also Fig. 8(c)]. The embedding medium is necessary to reduce reflections at the boundary of the phantom and reduce the image artifacts. To further reduce reflections at the Plexiglas/air interface, a 4-cm-thick layer of microwave absorbing foam is placed around the sidewalls of the tray. This has been shown to significantly reduce artifacts produced by the tray itself [24]. The tray is covered with a Plexiglas lid (4 mm thick) and placed on the measurement platform.
An Advantest R3770 VNA paired with an Advantest R3970 RF switch performs a measurement at every sampled position from 3 to 8 GHz in 100-MHz increments. Since the attenuation through the embedding medium is significant (at 8 GHz, it exceeds 65 dB), a power amplifier is used between the VNA and the TEM horn, supplying 39 dBm to the TEM horn antenna. Only transmission measurements are acquired across the bowtie antennas since reflection measurements cannot be acquired because of the nonreciprocal amplification stage.
The measurement setup is calibrated through two RO (or background) measurements and one CO (SP) measurement. The CO measurement is constructed with five uncut circular carbon rubber sheets made of the same material as the sheets used to make the compressed phantom. A cylindrical SP of permittivity r,sp = 50 − i 0.05 (made of microwave ceramics), and of diameter and height of 1 cm, is placed in the central sheet (Layer 3). Constructing the CO in this manner guarantees that the SP remains fixed in the desired central position during the movement of the platform.
Two RO measurements are performed to acquire S inc for each of the CO and the OUT measurements. The first RO consists of the Plexiglas tray filled only with the matching medium. It is used to extract the scattered data from the breastphantom measurements. The second RO includes not only the tray but also the uncut carbon rubber sheets without the SP. It is used to extract the PSFs. Since the CO contains carbon rubber sheets that secure the probe in place, the highest PSF SNR is achieved with an RO that matches the CO arrangement (without the probe).
The OUT is comprised of the compressed breast phantom placed in the embedding medium in the tray.

B. Results
QMH reconstructions of the phantom using the Born, the modified Rytov, the Marks, and the proposed combined Born/Rytov approach are shown in Fig. 9. The phase unwrapping of the OUT data, which is necessary for the reconstruction based on the Marks approximation, fails at several frequencies, as shown in Fig. 2. Nonetheless, the Marks deembedding formula (17) is applied with n = 2 for comparison.
The 2-D images in Fig. 9 are projections of the 3-D complex permittivity distribution of the phantom onto the z = 0 plane. The 3-D imaging is not feasible with this setup since it provides transmission measurements only, where the Rx antennas are either aligned with the Tx-antenna boresight or are only slightly off-boresight. This results in loss of range resolution [9], [57]. The range resolution can be improved by increasing the lateral (x y) distance between the Tx and Rx antennas. Unfortunately, this increases the signal path and thus the attenuation, leading to poor data SNR. When operating in 2-D mode, the QMH algorithm uses only one PSF, namely, the PSF corresponding to the position z = 0 of the SP. In this case, the system matrix in (28) has only one column, i.e., N z = 1, and the k-space solution vectorρ ρ ρ(κ κ κ) in (32) has a single value at each κ κ κ ≡ (k x , k y ). In effect, the permittivity values shown in Fig. 9 are averaged over the thickness (5.5 cm) of the phantom. This is why the locations of the tumor simulants yield a relatively small permittivity increase over that of the carbon rubber material.
We note that the 2-D reconstructions in Fig. 9 demonstrate a marked improvement of the quantitative accuracy compared to those reported in [20]. The improvement is due to a new k-space low-pass filtering strategy. Low-pass filtering is common in all methods performing the inversion in k-space in order to remove the high-frequency artifacts that are generated by the 2-D FFT of the data (see [17] for the application of such filtering with QMH). Here, an energy-conserving lowpass filter is proposed. In [17], a fourth-order Butterworth lowpass filter with a maximum value of 1 is used. This setting reduces significantly the aggregate complex permittivity image energy in k-space, resulting in underestimating the permittivity values in the (x, y) space. To conserve the image energy, the filter's maximum value M is determined so that its overall energy equals that of the unity square filter, leading to the expression whereB(k x , k y ) is the Butterworth filter of unity maximum and B(k x , k y ) = M ·B(k x , k y ) is its energy-conserving counterpart. The number of samples along k x and k y is N x and N y , respectively. A k-space low-pass filter is employed to suppress spatial frequencies beyond κ max > 2π/λ min along with an apodization filter in (x, y) [17].
The Born reconstruction [see Fig. 9(a)] outlines the fibroglandular structure and the cancerous inclusions marginally well in the real part of the permittivity. The imaginary permittivity image identifies only the tumor inclusions inside the fibroglandular tissue simulant. Gibb's (ripple) artifacts at the phantombackground interface are strong. They are caused by the abrupt change in the imaginary part of the permittivity at this interface (see Table I for the imaginary permittivity values of the carbon rubber and matching materials) as well as the gaps introduced by the plastic wrap and the air gaps trapped in it. Also, the size of the phantom is large, violating the first-order Born approximation limit.
In the modified Rytov reconstruction [see Fig. 9(b)], the two tumor-simulant groups, particularly in the imaginary part of the permittivity, are visible very well. Indeed, the response in the imaginary part of the permittivity is so significant that Gibbs artifacts are generated at their boundaries. This is due to the loss affecting significantly the magnitude of the measured OUT transmission coefficients. Importantly, the images are much less affected by ripple artifacts at the phantom edges. This is explained by the fact that only the magnitude of the OUT data is used, which is less sensitive to the presence of the very small gaps formed by the plastic wrap.
Utilizing the Marks approximation (with n = 2) leads to an increase in the overall noise of the reconstruction [see Fig. 9(c)] due to phase-wrapping problems. The artifacts clutter the image and reduce its overall fidelity, in particular in the imaginary response. Nonetheless, the real permittivity image localizes the tumor simulants.
The proposed combined Born/Rytov approach produces images, which contains structural features from both the Born-and Rytov-based reconstructions, which is expected. The significant rippling artifacts from the Born approximation are reduced. This results in a much clearer estimate of the fibroglandular structure while at the same time preserving the image features due to the embedded tumors. Overall, the image quality exceeds that of the images generated with the previous three approaches.
We note that the individual tumor simulants (two at each location) cannot be differentiated in neither of the images. The tumor simulants are spherical with 1 cm diameter and they are touching each other at both locations, i.e., the center-to-center distance is 1 cm. This is expected bearing in mind that the shortest wavelength in the carbon rubber material is 11.3 mm and the target viewing angle is severely limited by attenuation.

VII. CONCLUSION
In this work, the combined Born/Rytov QMH approach is proposed and explored via simulation-based and experimental imaging examples. The Born and Rytov approximations, when applied separately, yield markedly different images in the cases where their respective limits are violated. Depending on the specific imaging scenario, one approximation outperforms the other. The Born approximation performs better in recovering small scatterers and it tolerates a much wider range of object permittivity values in a low-permittivity background. However, when electrically large objects in a high-permittivity background are imaged, the Rytov-based reconstruction is advantageous and it exhibits an ability to better identify smaller scatterers embedded in larger objects, especially with transmission measurements. A major drawback of Rytov's approximation is its susceptibility to phase-unwrapping errors. This drawback is eliminated in the proposed approach, which removes the need to phase-unwrap the data. It only requires phase unwrapping of the system PSF, which is generally successful because the PSF is the response to a weak point-like scatterer. The integration of the two approximations into one approach allows to leverage their advantages without increasing the computational complexity of the reconstruction and preserving its real-time performance. We note that the experimental-and simulation-based reconstructions presented here are computed in MATLAB [51] within 2-3 s on a personal laptop without any parallel computing or code optimization.
One of the main applications of QMH is near-field biomedical imaging. In this application, the real-time reconstruction ability is critical since it dramatically increases the diagnostic capability. Within a minute, various image sets can be generated and image enhancement can be performed by: 1) tuning the ad hoc parameters of apodization and Fourier-space filters; 2) changing the data approximation strategy (e.g., Born or the proposed combined approach); and 3) employing system PSFs acquired in different backgrounds and with different SPs. Furthermore, it can potentially enable imaging during dynamic processes (e.g., breathing and circulation) in real time.
Future work aims to explore the use of multiple PSFs acquired with varying permittivities of the background and the SP to improve the quantitative accuracy of the images. The best quantitative accuracy is achieved in regions of the imaged objects, the permittivity of which is close to that of the SP [19]. Since the permittivity of the SP is known, it is in principle possible to determine which regions in the images have high quantitative accuracy and which need further enhancement. Also, further explorations into acquiring 3-D images without measuring SPs at each z slice could yield significant benefits by reducing system calibration effort.