# Permittivity profile estimation based on non-radiating equivalent source (2d case).

1. INTRODUCTION

Electromagnetic inverse scattering is one of the most promising techniques for medical imaging applications particularly where the existing medical imaging techniques are unreliable. Tooth interior imaging is an example [11]. MRI scanners are not appropriate for dental tomography because they are too expensive and not good at imaging teeth. In comparison with MRI, CT scanner is an obvious choice for dental radiology [14], but high required dose of ionizing radiation to one's head is an important concern [32]. Three-dimensional terahertz pulse imaging is another approach, which has its own challenges and difficulties [15], particularly for the large size tooth samples [11]. Three-dimensional terahertz continuous wave imaging [15] and pulse imaging [28] have been used for very low contrast objects. However, since these image reconstruction approaches are based on filter back projection (the scattered field or beam diffraction is ignored), the approaches can not be applied to the high-contrast objects or low-contrast objects, particularly with the large size one. The relative permittivities of tooth enamel and dentine have been measured by using THz spectroscopy in [2] and have been reported as 9.35 and 6.60, respectively.

Conventional dielectric profile estimation methods use Born approximation at a preliminary stage to solve the inverse scattering problem iteratively. The Born approximation initial guess has frequently been used to linearize the electromagnetic inverse scattering problem [8,15,16, 20,33-35]. This is a good initial estimate for the field inside a low-contrast scatterer as long as the scatterer size is a fraction of a wavelength [15]. This initial guess eases the formulation of the inverse scattering problem. However, the Born approximation initial guess was found in [15] to be a problematic assumption for a large size object (large in terms of wavelength).

In this paper, our focus is on the frequency domain inverse scattering algorithm. The existing methods for solving the electromagnetic inverse scattering problem in the frequency domain can be categorized under two main approaches: radiating and non-radiating.

The radiating approach takes into consideration only the radiating part of the total volume equivalent current source (VECS) and linearizes the electromagnetic inverse scattering problem. The radiating VECS is also known as the minimum energy solution [29, 30]. In a typical radiating approach, the electromagnetic inverse scattering systems are linearized by iteratively solving for the internal total electric field using the invertible part of the electromagnetic scattering Green's function. The resulting linear equation can be solved for the radiating part of VECS by means of the pseudo-inverse, mean square error [22,27,31], singular value decomposition (SVD) [25], regularization, statistical [1,5], or Fourier (holography) [37] based approaches. Initializing the internal total electric field with the incident field in the first iteration transforms the forward scattering problem into a linear equation [8, 15, 16, 33-35]. In the next iterations, the permittivity and total electric field are estimated iteratively. Iteration continues until either the scattered field estimation error or contrast factor estimation error drops below a certain threshold [8,9,12,15,16,33]. The threshold must first be set heuristically [8,9,12,23,33]. The permittivity profile of an object (scatterer) cannot be estimated with the radiating VECS alone. Signal-subspace optimization techniques are reported for permittivity profile estimation by extending the radiating objective function and minimizing the noise effects [6, 7,24].

The second approach includes the non-radiating VECS confined within the boundary of a scatterer. This approach involves the null space of the Green's function matrix of the scattering problems [4,10,17-19,26,29]. The internal scattered field inside an object is unmeasurable, and the non-radiating VECS cannot be obtained by using the invertible part of the Green's function operator in the aforementioned linearized iterative schemes. To our knowledge, no approach based on the non-radiating part of the VECS for permittivity profile reconstruction has so far been proposed.

In this paper, we propose an alternative approach based on a new non-radiating objective function for permittivity profile estimation of an unknown scatterer. We assume that an unknown scatterer consists of many homogenous regions whose boundaries are known a priori. It should be emphasized that our goal is to estimate the permittivity profile of unknown scatterers, but not the non-radiating VECS.

2. STATEMENT OF THE PROBLEM AND SUMMARY OF THE PROPOSED APPROACH

Figure 1 shows the electromagnetic inverse scattering system (EMISS). EMISS includes scatterers, a transmitter antenna, and multiple observation points (antenna array, R's). We assume that the unknown scatterer, [v.sub.2], is located inside the EMISS region, [v.sub.1]. Generated by the impressed or known source, [[??].sub.im], the total electric fields are measured at the observation points within [v.sub.1] in the presence of the scatterer. The unknown scatterer consists of many homogenous regions surrounded with a background medium.

The proposed approach estimates the permittivity profile using the data collected at observation points within [v.sub.1]. The proposed method is summarized as follows:

(i) Measuring the incident electric fields at observation points in EMISS in the absence of any scatterers,

(ii) Illuminating the scatterer by the incident field,

(iii) Measuring the total electric fields at the observation points,

(iv) Estimating the permittivity profile by minimizing an objective function including both radiating and non-radiating parts of the equivalent source respectively.

The focus of this paper is the detailed formulation of the fourth step of the above procedure, which is explained in the next two sections. Scattered field generated by the radiating part of the VECS is explained in Section 3. The new minimization (optimization) problem including non-radiating source, which is the essential part of the proposed method is presented in Section 4. The effectiveness of the proposed method is tested against both low-contrast and high contrast objects by conducting simulations. Section 5 presents the simulation results, and Section 6 concludes the paper.

3. RADIATING VECS AND "RADIATING" CONTRAST FACTOR

In this section, the radiating VECS and the radiating contrast factor are obtained by solving the forward scattering equation for the total electric fields measured at observation points. For a medium with a homogenous magnetic permeability profile, the total electric field in EMISS satisfies the complex vector wave equation:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (1)

where

[[??].sub.tot] = [[??].sub.inc] + [[??].sub.scat], (2)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (3)

and [[??].sub.tot], [[??].sub.scat], [[??].sub.inc], [[??].sub.tot], [[??].sub.im], [[??].sub.eq], [omega], [[epsilon].sub.0], and [[mu].sub.0] are the total electric field, scattered electric field, incident electric field, total electric current density (total current source), impressed current source, VECS, the angular frequency, free space permittivity, and free space permeability, respectively. The incident field is generated by [[??].sub.im] in the absence of any scattering object. The scattered field is generated by the VECS in a homogeneous medium. As shown in (3), the VECS is proportional to the contrast factor, [[kappa].sub.r],

[[kappa].sub.r]([??]) = [[epsilon].sub.r]([??]) - 1, (4)

and the total electric field, where [[epsilon].sub.r] is the relative permittivity. The scattered electric field can be obtained [21] as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (5)

where k is the wave number ([omega] [square root of ([micro][epsilon])]), and [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is the magnetic vector potential Green's function,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (6)

where [??] and [??]' are the position vectors that locate the observation point and the source point, respectively, as shown in Figure 1. Equation (4) indicates that the contrast factor is non-zero within the scatterer and is zero for a homogenous background. The scattered electric field Equation (5) can be discretized by dividing [v.sub.2] into the q number of elements. Using the Method of Moments (MOM), the scattered electric field Equation (5) can be written for the p number of observation points in a matrix form as follows:

[E.sub.scat] = [G.sub.e][J.sub.eq], (7)

where [E.sub.scat] is the p x 1 single column matrix, and [J.sub.eq] is the q x 1 total VECS single column matrix, and is defined as,

[J.sub.eq] = j[omega][[epsilon].sub.0][[kappa].sub.r] [E.sub.tot], (8)

[E.sub.tot] is the q x 1 single column matrix whose nth element is the average total electric field at the nth discretized element of [v.sub.2]; [[kappa].sub.r] is a q x q diagonal matrix whose nth diagonal element is the average contrast factor at the nth element of [v.sub.2]; [G.sub.e] is the p x q Green's function matrix wherein the mth row and nth column element of the electric field Green's function matrix, [G.sub.e mn], for two-dimensional TMz case, is:

[G.sub.e mn] = -j[omega][mu] [integral] [G.sub.a]{[??], [??]'} [f.sub.n]([??]') dv', (9)

[f.sub.n]([??]') is the nth basis function in the expansion of [J.sub.eq]. The pulse basis function is defined as,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (10)

Equation (7) is referred to as the forward scattering equation. The Green's function matrix is ill-conditioned [3,13] and rank deficient. The singular values of the Green's function matrix are first calculated by employing singular value decomposition. The singular values below a certain threshold are eliminated to improve the condition number of the Green's function matrix. The threshold can be calculated once for an EMISS at the calibration stage. The Green's function matrix in its now improved condition is used to solve (7) for [J.sub.eq]. The VECS is written as,

[J.sub.eq] = [J.sup.RA.sub.eq] + [J.sup.NR.sub.eq], (11)

where [J.sup.RA.sub.eq] (radiating VECS) is the radiating part of the equivalent source, whereas [J.sup.NR.sub.eq] (non-radiating VECS) is the remaining part of the equivalent source, which does not generate any field outside the scatterer. The inverse solution of Equation (7) yields only the radiating VECS.

Here, we assume that the relative contrast factor can be represented as,

[[kappa].sub.r] = [[kappa].sup.RA.sub.r] + [[kappa].sup.NR.sub.r], (12)

where [k.sup.RA.sub.r] is the relative radiating contrast factor, and [[kappa].sup.NR.sub.r] is the relative non-radiating contrast factor. Subdivided into m' number of homogenous subregions, the contrast factor of an inhomogeneous region is defined as a q x q block diagonal matrix as follows,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (13)

[I.sub.t] is the [q.sub.t] x [q.sub.t] identity matrix; m' is the number of homogenous clusters within a scatterer; [q.sub.t] is the number of elements in the tth subregion of [v.sub.2]; the sum of all the qt's, (t = 1, 2, ..., m') is the total elements, q, considered for the contrast factor estimation; [[kappa.sup.t.sub.r] is a scalar that represents the tth subregion's contrast factor.

The radiating contrast factor is entirely based on the radiating portion

of the VECS, and for each point, is defined as,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (14)

where [J.sup.RA h.sub.eq] and [E.sup.int RA h.sub.tot] are the corresponding radiating VECS and total electric field at a point denoted by h. The radiating internal total electric field, [E.sup.inta RA.sub.tot] is now defined by

[E.sup.int RA.sub.tot] = [E.sub.inc] + [G.sup.int.sub.e][J.sup.RA.sub.eq]. (15)

The Green's function matrix, [G.sup.int.sub.e], is evaluated for the observation points inside the scatterer. Note that the radiating contrast factor estimated above is one of the two parts of the contrast factor (12). To find the contrast factor or permittivity profile, we need to obtain the non-radiating part, which will be described in the next section.

4. NON-RADIATING VECS, NON-RADIATING CONTRAST FACTOR, AND THE OBJECTIVE FUNCTION

The non-radiating part of VECS cannot be obtained by solving the forward scattering equation directly, as the non-radiating part of VECS generates zero electric field outside a scatterer. The radiating VECS rigorously reproduces the external scattering field but fails to provide the correct internal scattered field via the forward scattering equations inside an object, particularly when the scatterer has a high contrast factor. Hence,

[E.sup.ext.sub.scat] = [G.sup.ext.sub.e] [J.sup.RA.sub.eq], (16)

[E.sup.int.sub.scat] [not equal to] [G.sup.int.sub.e] [J.sup.RA.sub.eq], (17)

where [E.sup.ext.sub.scat], [E.sup.int.sub.scat], [G.sup.ext.sub.e], and [G.sup.int.sub.e] are the external scattered field, internal scattered field, external electric field Green's function matrix, and internal electric field Green's function matrix, respectively. The non-radiating part of the VECS does not generate any fields outside the scatterer,

0 = [G.sup.ext.sub.e] [J.sup.NR.sub.eq]. (18)

The solutions to Equation (18) form the null space of the electromagnetic Green's function operator. Therefore, the VECS from (11) is non-unique.

The internal scattered field can be expressed in terms of the radiating and non-radiating parts of the total VECS within the scatterer:

[E.sup.int.sub.scat] = [G.sup.int.sub.e] ([J.sup.RA.sub.eq] + [J.sup.NR.sub.eq]). (19)

The total internal scattered field, [E.sup.int.sub.scat], can be decomposed into two parts, namely, the radiating internal scattered field, [E.sup.int RA.sub.scat], and the non-radiating internal scattered field, [E.sup.int NR.sub.scat],

[E.sup.int.sub.scat] = [E.sup.int RA.sub.scat] + [E.sup.int NR.sub.scat], (20)

where

[E.sup.int NR.sub.scat] = [G.sup.int.sub.e] [J.sup.NR.sub.eq]. (21)

Equations (2) and (3) can then be rewritten as follows by considering (11), (12), (20), and boundary information:

[E.sub.int.sub.tot] = [E.sub.inc] + [E.sup.int RA.sub.scat] + [E.sup.int NR.sub.scat], [J.sub.eq] = j[omega][[epsilon].sub./0]([[kappa].sup.RA.sub.r] + [[kappa].sup.NR.sub.r]) ([E.sub.inc] + [E.sup.int RA.sub.scat] + [E.sup.int NR.sub.scat]). (22)

The radiating VECS formulation can be written in a matrix form based on (14),

[J.sup.RA.sub.eq] = j[omega][[epsilon].sub.0][[kappa].sup.RA.sub.r] ([E.sup.int RA.sub.tot]). (23)

The non-radiating VECS can be obtained by replacing the VECS and the radiating VECS from (22) and (23), respectively, into (11):

[J.sup.NR.sub.eq] = j[omega][[epsilon].sub.0] (([[kappa].sup.RA.sub.r] + [[kappa].sup.NR.sub.r]) [E.sup.int NR.sub.scat] + [[kappa].sup.NR.sub.r] [E.sup.int RA.sub.tot]), (24)

where based on (15),

[E.sup.int RA.sub.tot] = [E.sub.inc] + [E.sup.int RA.sub.scat]. (25)

The non-radiating VECS given by (24) contains two unknowns, namely the non-radiating contrast factor, [[kappa].sup.NR.sub.r], and non-radiating internal scattered field, [E.sup.int NR.sub.tot]. Using (21) and (24), the non-radiating internal scattered field can be expressed in terms of the non-radiating contrast factor:

[E.sup.int NR.sub.scat] = j[omega][[epsilon].sub.0]Q[G.sup.int.sub.e][[kappa].sup.NR.sub.r][E.sup.int RA.sub.tot], (26)

where Q is

Q = [(I - j[omega][[epsilon].sub.0][G.sup.int.sub.e] ([[kappa].sup.RA.sub.r] + [[kappa].sup.NR.sub.r])).sup.-1], (27)

and Equation (24) can now be rewritten, as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (28)

Equations (18), (12), and (28) should be solved simultaneously to determine the non-radiating contrast factor and the contrast factor. For this purpose, rather than trying to solve (18) directly, we seek for the contrast factors that minimize the non-radiating objective function, [R.sub.NR](.), which is defined as the [l.sup.2]-Norm of the external scattered field due to the non-radiating VECS. Thus, the optimum answer to the permittivity profile estimation is the contrast factor minimizing the non-radiating objective function and is expressed as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (29)

where

[J.sup.NR i.sub.eq] = j[omega][[epsilon].sub.0][[kappa].sup.i.sub.r] (j[omega][[epsilon].sub.0][Q.sup.i][G.sup.int.sub.e] + j[omega][[epsilon].sub.0]) ([[kappa].sup.i.sub.r] - [[kappa].sup.RA.sub.r]) [E.sup.int RA.sub.tot], (30)

[Q.sup.i] = [(I - j[omega][[epsilon].sub.0][G.sup.int.sub.e][[kappa.sup.i.sub.r]]).sup.-1], (31)

[[kappa].sup.i.sub.r] and [[kappa].sup.RA i.sub.r] are the diagonal contrast factor matrix and the diagonal radiating contrast factor matrix for the ith test permittivity set, respectively. The contrast factor of an inhomogeneous region is defined as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (32)

where [I.sub.t] is the [q.sub.t] x [q.sub.t] identity matrix; [q.sub.t] is the number of elements in the tth subregion of [v.sub.2]; [[kappa].sup.i t.sub.r] is a scalar that represents tth subregion's contrast factor from the ith test permittivity set. Equation (29) in conjunction with (30) can be considered as the accurate formulation for estimating the contrast factor. The simulation results confirm that a unique contrast factor can be obtained by using Equation (29) in conjunction with (30).

The proposed objective function based on the non-radiating VECS includes a single unknown, the total contrast factor, while the radiating objective function linearized by applying the Born approximation initial guess includes two unknowns, the total contrast factor and the internal total electric field. To perform permittivity profile estimation, the search space dimension for the proposed approach is (n'), whereas the search space dimension for the Born iterative approach is (2 x n'). Therefore, the search complexity for the proposed approach is half that of the Born iterative approach. An interesting aspect of the above proposed approach is that it can provide a unique contrast factor that is not affected by the non-uniqueness of the non-radiating VECS problem.

5. SIMULATION RESULTS

In this section, we have evaluated the proposed approach for permittivity profile estimation by conducting several simulations. The permittivity profile of an inhomogeneous object can be represented with a number of homogenous subregions. We have considered a 2-D problem and a region bounded by a circle with a radius of 10[lambda]/3. We have chosen 255 observation points, distributed uniformly on the perimeter of the circle to ensure enough samples to accurately capture the full scattered electric field. To prevent the inverse crime scenario [36], the electromagnetic forward scattering problem is independently simulated using a finite element method (FEM). The FEM is used to illuminate the object under test with a plane wave and to collect data (the total electric field) at the observation points. The collected data generated independently by FEM is used and considered as the input for permittivity profile estimation. The proposed approach is verified by its application to two different case studies as explained below.

In the first case study, we have considered the EMISS geometrical configuration shown in Figure 2(a). The object under test, [v.sub.2], includes two homogenous subregions, [v.sup.1.sub.2] and [v.sup.2.sub.2], in a free space background. The first subregion is a circular cylindrical dielectric one wavelength in diameter, and its center lies at (0.6[lambda], 1.8[lambda]). The second subregion is a rectangular cylindrical dielectric each side one wavelength wide, and its center lies at (-0.6[lambda], -1.8[lambda]). Without loss of generality, the dielectric scatterers are considered to be lossless and divided into 5810 elements.

To verify the performance of the proposed approach for estimating the permittivity profile of a low-contrast object, first, the scatterer whose discretized permittivity profile is shown in Figure 3(a) is illuminated by a plane wave. The forward scattering is solved independently once by FEM, and data is collected at the observation points. The search space can always be restricted to the prior known range of the dielectric constant of a structure. As an example, if the approach is applied to a tooth structure, and its relative permittivity range is known for tooth enamel and dentine, the search space is bounded within this range (between 1.0 and 10.0). Searching over the permittivity range enables us to distinguish between the tooth enamel and dentine subregions in this particular example. Figure 4(a) presents the proposed objective function based on the non-radiating VECS for the object with two subregions with the test permittivities ranging from 1.25 to 9.00 (n' = 961). The top (x-y) view of the objective functions is shown in Figure 4(b). The Born iterative method's objective function is converted to a linear one by initializing one of the unknowns, the internal total electric field. However, there is no such simplification in the proposed objective function, and the objective function remains non-linear. The non-radiating objective function is minimal when the estimated permittivity of the two regions is 1.25, which is the correct value.

Secondly, we investigated the performance of the proposed objective function for estimating the permittivity profile of a high-contrast object. As shown in Figure 3(b), the high-contrast scatterer is illuminated by a plane wave. The non-radiating objective function for the object with two subregions with the permittivities ranging from 1.25 to 9.00 (n' = 961) is presented in Figure 5. As shown in Figures 5(a) and (b), the non-radiating objective function and the non-radiating approach enable us to estimate the relative permittivity of the [v.sup.1.sub.2] and [v.sup.2.sub.2] regions as 6.00, which is the correct value.

Next, we investigated the performance of the proposed nonradiating objective function for estimating the permittivity profile of an object including two subregions with different permittivity values. The permittivity profile of the scatterer under test is shown in Figure 3(c). As illustrated in Figures 6(a) and (b), the non-radiating objective function has distinctive minima when the relative permittivities of the [v.sup.1.sub.2] and [v.sup.2.sub.2] regions are, respectively, 3.00 and 6.00, which are the correct values.

In the second case study, the performance of the proposed approach is evaluated for the object including more than two homogenous subregions. In this case, we considered the EMISS geometrical configuration of Figure 2(b). The object under test, [v.sub.2], includes three homogenous subregions: [v.sup.1.sub.2], [v.sup.2.sub.2] and [v.sup.3.sub.2] in a free space background. The first subregion is a circular cylindrical dielectric with a diameter of one wavelength centered at (0.6[lambda], 1.8[lambda]). The second subregion is a rectangular cylindrical dielectric with one wavelength side centered at (1.5[lambda], -1.5[lambda]). The third subregion is a rectangular cylindrical dielectric with a half wavelength width and a wavelength length centered at (-1.8[lambda], 0). The dielectric scatterers are considered to be lossless and are subdivided into 5046 elements.

We investigated the performance of the proposed non-radiating objective function for estimating the permittivity profile of a high-contrast object. The high-contrast scatterer's permittivity profile under test is shown in Figure 7(a). The non-radiating objective function for the object with three subregions, with permittivities ranging from 1.25 to 9.00 (n' = 729) are presented in Figure 7(b). The non-radiating objective function is minimal when the relative permittivity of regions [v.sup.1.sub.2], [v.sup.2.sub.2] and [v.sup.3.sub.2] is 6.00, which is the correct value.

Then, we investigated the performance of the proposed approach for estimating the permittivity profile of an object including three subregions with different permittivity values as shown in Figure 8(a). The non-radiating objective function for the object with three subregions is presented in Figure 8(b). The non-radiating objective function is at minimum when the permittivity profile amplitudes of regions [v.sup.1.sub.2], [v.sup.2.sub.2] and [v.sup.3.sub.2] are, respectively, 3.00, 5.00 and 8.00, which are the correct values.

So far, the proposed approach has been verified for the noise free data. To evaluate the proposed approach performance with the noisy data, the high-contrast object with the permittivity profile shown in Figure 3(b) is illuminated with a plane wave, and a white Gaussian noise is added to the FEM simulation results. The non-radiating objective function is evaluated by sweeping the relative permittivity with the 0.25 step size between 1 and 9 for the following signal to noise ratio (SNR): 60db, 40db, 20db, 10db, 5db, 3db, 2db, and 1db. The results indicate that the minimum values for the non-radiating objective function occur at the correct relative permittivity ([[epsilon].sub.r] = 6) for the 60 db, 40 db, 20 db, and 10 db SNRs, but it deviates from the true relative permittivity when the SNR drops below 10 db, as shown in Figure 9(a). Thus, the permittivity profile estimation based on the proposed approach is also valid for noisy measured data. The permittivity profile estimation error for the 5db, 3db, and 2db SNRs is %8, and for the 1db SNR is %14, as shown in Figure 9(b).

The simulation results indicate that the proposed approach can be successfully applied to both low-contrast and high-contrast permittivity profiles of a large object (in terms of wavelength). The simulation results illustrate that the proposed approach can correctly estimate the permittivity profile of the object under test in a noisy measurement environment.

6. CONCLUSION

A new solution to the permittivity profile estimation problem is presented for a scatterer whose permittivity profile can be divided into a number of homogenous clusters. To address the Born approximation initial guess issue pointed out in [15] for the electromagnetic inverse scattering problem, we have introduced the non-radiating contrast factor through which an objective function for permittivity profile estimation is formulated. The solutions to the inverse source problem are non-unique. However, our method yields the correct and unique permittivity profile of an unknown object by minimizing the non-radiating objective function and applying the boundary information. The conventional Born's approximation works well for a low-contrast object as long as the object size is a fraction of a wavelength. In contrast, the proposed non-radiating objective function can be used for both low-contrast and high-contrast permittivity profiles with even large size objects. The method has been verified by extensive simulations. The proposed approach has a much smaller number of unknowns and is therefore computationally more efficient than the permittivity profile estimation approaches based on Born approximation. The simulation results also depict that the proposed approach can accurately estimate the permittivity profile of an object under test in a noisy environment.

ACKNOWLEDGMENT

This work has been supported by NSERC (Natural Sciences and Engineering Research Council) of Canada and BlackBerry (former Research In Motion).

The author would also like to thank Prof. J. Orchard, Dr. D. Saeedkia, Dr. M. Neshat, Dr. D. Hailu, and Ms. N. Ranjkesh for the valuable technical discussions, and Ms. M. McPherson for her invaluable help in editing this paper.

REFERENCES

[1.] Autierin, R., M. Durso, C. Eyraud, A. Litman, V. Pascazio, and T. Isernia, "3d inversion of fresneldatabase with 3d Markov random field," Session 3P6 Inverse Scattering Probelms: Open Problems and New Challenges, 559, 2010.

[2.] Berry, E., A. J. Fitzgerald, N. N. Zinov'ev, G. C. Walker, S. Homer-Vanniasinkam, C. D. Sudworth, R. E. Miles, J. M. Chamberlain, and M. A. Smith, "Optical properties of tissue measured using terahertz pulsed imaging," Proceedings of SPIE, Vol. 5030, 459-470, 2003.

[3.] Bertero, M., T. A. Poggio, and V. Torre, "Ill-posed problems in early vision," Proceedings of the IEEE, Vol. 76, No. 8, 869-889, 1988.

[4.] Bojarski, N. N., "Inverse scattering," Technical report, DTIC Document, 1973.

[5.] Caorsi, S., G. Gragnani, M. Pastorino, G. Zunino, et al, "Microwave imaging based on a Markov random field model," IEEE Transactions on Antennas and Propagation, Vol, 42, No. 3, 293-303, 1994.

[6.] Chen, X., "Application of signal-subspace and optimization method in reconstructing extended scatterers," JOSA A, Vol. 26, 1022-1026, 2009.

[7.] Chen, X., "Subspace-based optimization method for solving inverse scattering problems," IEEE Trans. Geosci. Remote Sens, Vol. 48, No. 1, 42-49, 2010.

[8.] Chew, W. C. and Y. M. Wang, "Reconstruction of two-dimensional permittivity distribution using the distorted born iterative method," IEEE Transactions on Medical Imaging, Vol. 9, No. 2, 218-225, 1990.

[9.] Chew, W. C., Waves and Fields in Inhomogeneous Media, IEEE Press, 1995.

[10.] Chiappe, M. and G. L. Gragnani, "Analytical solution to inverse electromagnetic scattering: Shape and position reconstruction of dielectric objects," IEEE International Workshop on Imaging Systems and Techniques, 1-6, 2007.

[11.] Crawley, D., C. Longbottomm, V. P. Wallace, B. Cole, D. Arnone, and M. Pepper, "Three-dimensional terahertz pulse imaging of dental tissue," Journal of Biomedical Optics, Vol. 8, No. 2, 303-307, 2003.

[12.] Cui, T. J., W. C. Chew, A. A. Aydiner, and S. Chen, "Inverse scattering of two-dimensional dielectric objects buried in a lossy earth using the distorted born iterative method," IEEE Trans. Geosci. Remote Sens., Vol. 39, No. 2, 339-346, 2001.

[13.] Fieguth, P., Statistical Image Processing and Multidimensional Modeling, Vol. 155, Spring, 2010

[14.] Scaf, C. E. S. G., J. Morais, and L. C. M. Loffredo, "A survey of radiographic measurement estimation in assessment of dental implant length," Oral Surgery, Oral Medicine, Oral Pathology, Oral, Radiology and Endodontology, Vol. 103, No. 53, 2007.

[15.] Gilmore, C., P. Mojabi, A. Zakaria, M. Ostadrahimi, C. Kaye, S. Noghanian, L. Shafai, S. Pistorius, and J. LoVetri, "A wideband microwave tomography system with a novel frequency selection procedure," IEEE Transactions on Biomedical Engineering, Vol. 57, No. 4, 894-903, 2010.

[16.] Gilmore, C., P. Mojabi, A. Zakaria, S. Pistorius, and J. LoVetri, "On super-resolution with an experimental microwave tomography system," IEEE Antennas and Wireless Propagation Letters, Vol. 9, 393-396, 2010.

[17.] Gragnani, G., "Electromagnetic imaging using closed-form radiating and non-radiating currents," Mediterranean Microwave Symposium, 1-6, November 2009.

[18.] Gragnani, G., "Two-dimensional nonradiating currents for imaging systems: Theoretical development and preliminary assessment," IET Microwaves, Antennas Propagation, Vol. 3, No. 8, 1164-1171, 2009.

[19.] Gragnani, G. and M. D. Mendez, "An improved electromagnetic imaging procedure using non-radiating sources," Mediterranean Microwave Symposium, 188-191, 2010.

[20.] Hajihashemi, M. R. and M. El-Shenawee, "Inverse scattering of three-dimensional PEC objects using the level-set method," Progress In Electromagnetics Research, Vol. 116, 23-47, 2011.

[21.] Harrington, R. F., Time-Harmonic Electromagnetic Fields, The IEEE Press Series, Pscataway, NJ, 2001.

[22.] Johnson, S., Y. Tae-Hoon, and R. Jung-Woong, "Inverse scattering solutions of scalar Helmholtz wave equation by a multiple source moment method," Electron. Lett., Vol. 19, 130-132, 1983.

[23.] Lavarello, R. J. and M. L. Oelze, "Tomographic reconstruction of three-dimensional volumes using the distorted born iterative method," IEEE Transactions on Medical Imaging, Vol. 28, No. 10, 1643-1653, 2009.

[24.] Lee, K. C., J.-S. Ou, and M. C. Fang, "Application of SVD noise-reduction technique to PCA based radar target recognition," Progress In Electromagnetics Research, Vol. 81, 447-459, 2008.

[25.] Markel, V. A., V. Mital, and J. C. Schotland, "Inverse problem in optical diffusion tomography iii inversion formulas and singular value decomposition," J. Opt. Soc. Am. A, Vol. 20, No. 5, 890-902, 2003.

[26.] Muller, C., "Electromagnetic radiation patterns and sources," IRE Transactions on Antennas and Propagation, Vol. 4, No. 3, 224-232, 1956.

[27.] Ney, M. M., A. M. Smith, and S. S. Stuchly, "A solution of electromagnetic imaging using pseudoinverse transformation," IEEE Transactions on Medical Imaging, Vol. 3, No. 4, 1984.

[28.] Nguyen, K. L., M. L. Johns, and L. F. Gladden, "Three dimensional imaging with a terahertz quantum cascade laser," Optics Express, Vol. 14, No. 6, 2006.

[29.] Porter, R. and A. Devaney, "Holography and the inverse source problem," Journal of Optical Society of America, Vol. 72, 327-330, 1982.

[30.] Porter, R. P., "Diffraction-limited, scalar image formation with holograms of arbitrary shape," Journal of the Optical Society of America, Vol. 60, No. 8, 1051-1059, 1970.

[31.] Richmond, J. H., "Scattering by a dielectric cylinder of arbitrary cross section shape," IEEE Antennas and Wireless Propagation, Vol. 13, No. 3, 334-341, 1965.

[32.] Tachibana, H. and K. Matsumoto, "Applicability of x-ray computerized tomography in endodontics," Dental Traumatology, Vol. 6, No. 1, 16-20, 1990.

[33.] Wang, Y. and W. Chew, "An iterative solution of the two-dimensional electromagnetic inverse scattering problem," International Journal of Imaging Systems and Technology, Vol. 1, No. 1, 100-108, 1989.

[34.] Winters, D. W., B. D. Van Veen, and S. C. Hagness, "Estimation of the frequency-dependent average dielectric properties of breast tissue using a time-domain inverse scattering technique," IEEE Transactions on Antennas and Propagation, Vol. 54, 3517-3528, 2006.

[35.] Winters, D. W., B. D. Van Veen, and S. C. Hagness, "A sparsity regularization approach to the electromagnetic inverse scattering problem. IEEE Transactions on Antennas Propagation, Vol. 58, 145-154, 2010.

[36.] Wirgin, A., "The inverse crime," ArXiv preprint mathph/0401050, 2004.

[37.] Wolf, E., "Three-dimensional structure determination of semitransparent objects from holographic data," Optical Communication, Vol. 1, No. 4, 153-156, 1969.

Shahed Shahir *, Mehrbod Mohajer, Arash Rohani, and Safieddin Safavi-Naeini

Department of Electrical and Computer Engineering, University of Waterloo, Waterloo, Ontario N2L 3G1, Canada

Received 23 January 2013, Accepted 28 March 2013, Scheduled 2 April 2013

* Corresponding author: Shahed Shahir (sshahir@maxwell.uwaterloo.ca).

Electromagnetic inverse scattering is one of the most promising techniques for medical imaging applications particularly where the existing medical imaging techniques are unreliable. Tooth interior imaging is an example [11]. MRI scanners are not appropriate for dental tomography because they are too expensive and not good at imaging teeth. In comparison with MRI, CT scanner is an obvious choice for dental radiology [14], but high required dose of ionizing radiation to one's head is an important concern [32]. Three-dimensional terahertz pulse imaging is another approach, which has its own challenges and difficulties [15], particularly for the large size tooth samples [11]. Three-dimensional terahertz continuous wave imaging [15] and pulse imaging [28] have been used for very low contrast objects. However, since these image reconstruction approaches are based on filter back projection (the scattered field or beam diffraction is ignored), the approaches can not be applied to the high-contrast objects or low-contrast objects, particularly with the large size one. The relative permittivities of tooth enamel and dentine have been measured by using THz spectroscopy in [2] and have been reported as 9.35 and 6.60, respectively.

Conventional dielectric profile estimation methods use Born approximation at a preliminary stage to solve the inverse scattering problem iteratively. The Born approximation initial guess has frequently been used to linearize the electromagnetic inverse scattering problem [8,15,16, 20,33-35]. This is a good initial estimate for the field inside a low-contrast scatterer as long as the scatterer size is a fraction of a wavelength [15]. This initial guess eases the formulation of the inverse scattering problem. However, the Born approximation initial guess was found in [15] to be a problematic assumption for a large size object (large in terms of wavelength).

In this paper, our focus is on the frequency domain inverse scattering algorithm. The existing methods for solving the electromagnetic inverse scattering problem in the frequency domain can be categorized under two main approaches: radiating and non-radiating.

The radiating approach takes into consideration only the radiating part of the total volume equivalent current source (VECS) and linearizes the electromagnetic inverse scattering problem. The radiating VECS is also known as the minimum energy solution [29, 30]. In a typical radiating approach, the electromagnetic inverse scattering systems are linearized by iteratively solving for the internal total electric field using the invertible part of the electromagnetic scattering Green's function. The resulting linear equation can be solved for the radiating part of VECS by means of the pseudo-inverse, mean square error [22,27,31], singular value decomposition (SVD) [25], regularization, statistical [1,5], or Fourier (holography) [37] based approaches. Initializing the internal total electric field with the incident field in the first iteration transforms the forward scattering problem into a linear equation [8, 15, 16, 33-35]. In the next iterations, the permittivity and total electric field are estimated iteratively. Iteration continues until either the scattered field estimation error or contrast factor estimation error drops below a certain threshold [8,9,12,15,16,33]. The threshold must first be set heuristically [8,9,12,23,33]. The permittivity profile of an object (scatterer) cannot be estimated with the radiating VECS alone. Signal-subspace optimization techniques are reported for permittivity profile estimation by extending the radiating objective function and minimizing the noise effects [6, 7,24].

The second approach includes the non-radiating VECS confined within the boundary of a scatterer. This approach involves the null space of the Green's function matrix of the scattering problems [4,10,17-19,26,29]. The internal scattered field inside an object is unmeasurable, and the non-radiating VECS cannot be obtained by using the invertible part of the Green's function operator in the aforementioned linearized iterative schemes. To our knowledge, no approach based on the non-radiating part of the VECS for permittivity profile reconstruction has so far been proposed.

In this paper, we propose an alternative approach based on a new non-radiating objective function for permittivity profile estimation of an unknown scatterer. We assume that an unknown scatterer consists of many homogenous regions whose boundaries are known a priori. It should be emphasized that our goal is to estimate the permittivity profile of unknown scatterers, but not the non-radiating VECS.

2. STATEMENT OF THE PROBLEM AND SUMMARY OF THE PROPOSED APPROACH

Figure 1 shows the electromagnetic inverse scattering system (EMISS). EMISS includes scatterers, a transmitter antenna, and multiple observation points (antenna array, R's). We assume that the unknown scatterer, [v.sub.2], is located inside the EMISS region, [v.sub.1]. Generated by the impressed or known source, [[??].sub.im], the total electric fields are measured at the observation points within [v.sub.1] in the presence of the scatterer. The unknown scatterer consists of many homogenous regions surrounded with a background medium.

The proposed approach estimates the permittivity profile using the data collected at observation points within [v.sub.1]. The proposed method is summarized as follows:

(i) Measuring the incident electric fields at observation points in EMISS in the absence of any scatterers,

(ii) Illuminating the scatterer by the incident field,

(iii) Measuring the total electric fields at the observation points,

(iv) Estimating the permittivity profile by minimizing an objective function including both radiating and non-radiating parts of the equivalent source respectively.

The focus of this paper is the detailed formulation of the fourth step of the above procedure, which is explained in the next two sections. Scattered field generated by the radiating part of the VECS is explained in Section 3. The new minimization (optimization) problem including non-radiating source, which is the essential part of the proposed method is presented in Section 4. The effectiveness of the proposed method is tested against both low-contrast and high contrast objects by conducting simulations. Section 5 presents the simulation results, and Section 6 concludes the paper.

3. RADIATING VECS AND "RADIATING" CONTRAST FACTOR

In this section, the radiating VECS and the radiating contrast factor are obtained by solving the forward scattering equation for the total electric fields measured at observation points. For a medium with a homogenous magnetic permeability profile, the total electric field in EMISS satisfies the complex vector wave equation:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (1)

where

[[??].sub.tot] = [[??].sub.inc] + [[??].sub.scat], (2)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (3)

and [[??].sub.tot], [[??].sub.scat], [[??].sub.inc], [[??].sub.tot], [[??].sub.im], [[??].sub.eq], [omega], [[epsilon].sub.0], and [[mu].sub.0] are the total electric field, scattered electric field, incident electric field, total electric current density (total current source), impressed current source, VECS, the angular frequency, free space permittivity, and free space permeability, respectively. The incident field is generated by [[??].sub.im] in the absence of any scattering object. The scattered field is generated by the VECS in a homogeneous medium. As shown in (3), the VECS is proportional to the contrast factor, [[kappa].sub.r],

[[kappa].sub.r]([??]) = [[epsilon].sub.r]([??]) - 1, (4)

and the total electric field, where [[epsilon].sub.r] is the relative permittivity. The scattered electric field can be obtained [21] as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (5)

where k is the wave number ([omega] [square root of ([micro][epsilon])]), and [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is the magnetic vector potential Green's function,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (6)

where [??] and [??]' are the position vectors that locate the observation point and the source point, respectively, as shown in Figure 1. Equation (4) indicates that the contrast factor is non-zero within the scatterer and is zero for a homogenous background. The scattered electric field Equation (5) can be discretized by dividing [v.sub.2] into the q number of elements. Using the Method of Moments (MOM), the scattered electric field Equation (5) can be written for the p number of observation points in a matrix form as follows:

[E.sub.scat] = [G.sub.e][J.sub.eq], (7)

where [E.sub.scat] is the p x 1 single column matrix, and [J.sub.eq] is the q x 1 total VECS single column matrix, and is defined as,

[J.sub.eq] = j[omega][[epsilon].sub.0][[kappa].sub.r] [E.sub.tot], (8)

[E.sub.tot] is the q x 1 single column matrix whose nth element is the average total electric field at the nth discretized element of [v.sub.2]; [[kappa].sub.r] is a q x q diagonal matrix whose nth diagonal element is the average contrast factor at the nth element of [v.sub.2]; [G.sub.e] is the p x q Green's function matrix wherein the mth row and nth column element of the electric field Green's function matrix, [G.sub.e mn], for two-dimensional TMz case, is:

[G.sub.e mn] = -j[omega][mu] [integral] [G.sub.a]{[??], [??]'} [f.sub.n]([??]') dv', (9)

[f.sub.n]([??]') is the nth basis function in the expansion of [J.sub.eq]. The pulse basis function is defined as,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (10)

Equation (7) is referred to as the forward scattering equation. The Green's function matrix is ill-conditioned [3,13] and rank deficient. The singular values of the Green's function matrix are first calculated by employing singular value decomposition. The singular values below a certain threshold are eliminated to improve the condition number of the Green's function matrix. The threshold can be calculated once for an EMISS at the calibration stage. The Green's function matrix in its now improved condition is used to solve (7) for [J.sub.eq]. The VECS is written as,

[J.sub.eq] = [J.sup.RA.sub.eq] + [J.sup.NR.sub.eq], (11)

where [J.sup.RA.sub.eq] (radiating VECS) is the radiating part of the equivalent source, whereas [J.sup.NR.sub.eq] (non-radiating VECS) is the remaining part of the equivalent source, which does not generate any field outside the scatterer. The inverse solution of Equation (7) yields only the radiating VECS.

Here, we assume that the relative contrast factor can be represented as,

[[kappa].sub.r] = [[kappa].sup.RA.sub.r] + [[kappa].sup.NR.sub.r], (12)

where [k.sup.RA.sub.r] is the relative radiating contrast factor, and [[kappa].sup.NR.sub.r] is the relative non-radiating contrast factor. Subdivided into m' number of homogenous subregions, the contrast factor of an inhomogeneous region is defined as a q x q block diagonal matrix as follows,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (13)

[I.sub.t] is the [q.sub.t] x [q.sub.t] identity matrix; m' is the number of homogenous clusters within a scatterer; [q.sub.t] is the number of elements in the tth subregion of [v.sub.2]; the sum of all the qt's, (t = 1, 2, ..., m') is the total elements, q, considered for the contrast factor estimation; [[kappa.sup.t.sub.r] is a scalar that represents the tth subregion's contrast factor.

The radiating contrast factor is entirely based on the radiating portion

of the VECS, and for each point, is defined as,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (14)

where [J.sup.RA h.sub.eq] and [E.sup.int RA h.sub.tot] are the corresponding radiating VECS and total electric field at a point denoted by h. The radiating internal total electric field, [E.sup.inta RA.sub.tot] is now defined by

[E.sup.int RA.sub.tot] = [E.sub.inc] + [G.sup.int.sub.e][J.sup.RA.sub.eq]. (15)

The Green's function matrix, [G.sup.int.sub.e], is evaluated for the observation points inside the scatterer. Note that the radiating contrast factor estimated above is one of the two parts of the contrast factor (12). To find the contrast factor or permittivity profile, we need to obtain the non-radiating part, which will be described in the next section.

4. NON-RADIATING VECS, NON-RADIATING CONTRAST FACTOR, AND THE OBJECTIVE FUNCTION

The non-radiating part of VECS cannot be obtained by solving the forward scattering equation directly, as the non-radiating part of VECS generates zero electric field outside a scatterer. The radiating VECS rigorously reproduces the external scattering field but fails to provide the correct internal scattered field via the forward scattering equations inside an object, particularly when the scatterer has a high contrast factor. Hence,

[E.sup.ext.sub.scat] = [G.sup.ext.sub.e] [J.sup.RA.sub.eq], (16)

[E.sup.int.sub.scat] [not equal to] [G.sup.int.sub.e] [J.sup.RA.sub.eq], (17)

where [E.sup.ext.sub.scat], [E.sup.int.sub.scat], [G.sup.ext.sub.e], and [G.sup.int.sub.e] are the external scattered field, internal scattered field, external electric field Green's function matrix, and internal electric field Green's function matrix, respectively. The non-radiating part of the VECS does not generate any fields outside the scatterer,

0 = [G.sup.ext.sub.e] [J.sup.NR.sub.eq]. (18)

The solutions to Equation (18) form the null space of the electromagnetic Green's function operator. Therefore, the VECS from (11) is non-unique.

The internal scattered field can be expressed in terms of the radiating and non-radiating parts of the total VECS within the scatterer:

[E.sup.int.sub.scat] = [G.sup.int.sub.e] ([J.sup.RA.sub.eq] + [J.sup.NR.sub.eq]). (19)

The total internal scattered field, [E.sup.int.sub.scat], can be decomposed into two parts, namely, the radiating internal scattered field, [E.sup.int RA.sub.scat], and the non-radiating internal scattered field, [E.sup.int NR.sub.scat],

[E.sup.int.sub.scat] = [E.sup.int RA.sub.scat] + [E.sup.int NR.sub.scat], (20)

where

[E.sup.int NR.sub.scat] = [G.sup.int.sub.e] [J.sup.NR.sub.eq]. (21)

Equations (2) and (3) can then be rewritten as follows by considering (11), (12), (20), and boundary information:

[E.sub.int.sub.tot] = [E.sub.inc] + [E.sup.int RA.sub.scat] + [E.sup.int NR.sub.scat], [J.sub.eq] = j[omega][[epsilon].sub./0]([[kappa].sup.RA.sub.r] + [[kappa].sup.NR.sub.r]) ([E.sub.inc] + [E.sup.int RA.sub.scat] + [E.sup.int NR.sub.scat]). (22)

The radiating VECS formulation can be written in a matrix form based on (14),

[J.sup.RA.sub.eq] = j[omega][[epsilon].sub.0][[kappa].sup.RA.sub.r] ([E.sup.int RA.sub.tot]). (23)

The non-radiating VECS can be obtained by replacing the VECS and the radiating VECS from (22) and (23), respectively, into (11):

[J.sup.NR.sub.eq] = j[omega][[epsilon].sub.0] (([[kappa].sup.RA.sub.r] + [[kappa].sup.NR.sub.r]) [E.sup.int NR.sub.scat] + [[kappa].sup.NR.sub.r] [E.sup.int RA.sub.tot]), (24)

where based on (15),

[E.sup.int RA.sub.tot] = [E.sub.inc] + [E.sup.int RA.sub.scat]. (25)

The non-radiating VECS given by (24) contains two unknowns, namely the non-radiating contrast factor, [[kappa].sup.NR.sub.r], and non-radiating internal scattered field, [E.sup.int NR.sub.tot]. Using (21) and (24), the non-radiating internal scattered field can be expressed in terms of the non-radiating contrast factor:

[E.sup.int NR.sub.scat] = j[omega][[epsilon].sub.0]Q[G.sup.int.sub.e][[kappa].sup.NR.sub.r][E.sup.int RA.sub.tot], (26)

where Q is

Q = [(I - j[omega][[epsilon].sub.0][G.sup.int.sub.e] ([[kappa].sup.RA.sub.r] + [[kappa].sup.NR.sub.r])).sup.-1], (27)

and Equation (24) can now be rewritten, as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (28)

Equations (18), (12), and (28) should be solved simultaneously to determine the non-radiating contrast factor and the contrast factor. For this purpose, rather than trying to solve (18) directly, we seek for the contrast factors that minimize the non-radiating objective function, [R.sub.NR](.), which is defined as the [l.sup.2]-Norm of the external scattered field due to the non-radiating VECS. Thus, the optimum answer to the permittivity profile estimation is the contrast factor minimizing the non-radiating objective function and is expressed as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (29)

where

[J.sup.NR i.sub.eq] = j[omega][[epsilon].sub.0][[kappa].sup.i.sub.r] (j[omega][[epsilon].sub.0][Q.sup.i][G.sup.int.sub.e] + j[omega][[epsilon].sub.0]) ([[kappa].sup.i.sub.r] - [[kappa].sup.RA.sub.r]) [E.sup.int RA.sub.tot], (30)

[Q.sup.i] = [(I - j[omega][[epsilon].sub.0][G.sup.int.sub.e][[kappa.sup.i.sub.r]]).sup.-1], (31)

[[kappa].sup.i.sub.r] and [[kappa].sup.RA i.sub.r] are the diagonal contrast factor matrix and the diagonal radiating contrast factor matrix for the ith test permittivity set, respectively. The contrast factor of an inhomogeneous region is defined as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (32)

where [I.sub.t] is the [q.sub.t] x [q.sub.t] identity matrix; [q.sub.t] is the number of elements in the tth subregion of [v.sub.2]; [[kappa].sup.i t.sub.r] is a scalar that represents tth subregion's contrast factor from the ith test permittivity set. Equation (29) in conjunction with (30) can be considered as the accurate formulation for estimating the contrast factor. The simulation results confirm that a unique contrast factor can be obtained by using Equation (29) in conjunction with (30).

The proposed objective function based on the non-radiating VECS includes a single unknown, the total contrast factor, while the radiating objective function linearized by applying the Born approximation initial guess includes two unknowns, the total contrast factor and the internal total electric field. To perform permittivity profile estimation, the search space dimension for the proposed approach is (n'), whereas the search space dimension for the Born iterative approach is (2 x n'). Therefore, the search complexity for the proposed approach is half that of the Born iterative approach. An interesting aspect of the above proposed approach is that it can provide a unique contrast factor that is not affected by the non-uniqueness of the non-radiating VECS problem.

5. SIMULATION RESULTS

In this section, we have evaluated the proposed approach for permittivity profile estimation by conducting several simulations. The permittivity profile of an inhomogeneous object can be represented with a number of homogenous subregions. We have considered a 2-D problem and a region bounded by a circle with a radius of 10[lambda]/3. We have chosen 255 observation points, distributed uniformly on the perimeter of the circle to ensure enough samples to accurately capture the full scattered electric field. To prevent the inverse crime scenario [36], the electromagnetic forward scattering problem is independently simulated using a finite element method (FEM). The FEM is used to illuminate the object under test with a plane wave and to collect data (the total electric field) at the observation points. The collected data generated independently by FEM is used and considered as the input for permittivity profile estimation. The proposed approach is verified by its application to two different case studies as explained below.

In the first case study, we have considered the EMISS geometrical configuration shown in Figure 2(a). The object under test, [v.sub.2], includes two homogenous subregions, [v.sup.1.sub.2] and [v.sup.2.sub.2], in a free space background. The first subregion is a circular cylindrical dielectric one wavelength in diameter, and its center lies at (0.6[lambda], 1.8[lambda]). The second subregion is a rectangular cylindrical dielectric each side one wavelength wide, and its center lies at (-0.6[lambda], -1.8[lambda]). Without loss of generality, the dielectric scatterers are considered to be lossless and divided into 5810 elements.

To verify the performance of the proposed approach for estimating the permittivity profile of a low-contrast object, first, the scatterer whose discretized permittivity profile is shown in Figure 3(a) is illuminated by a plane wave. The forward scattering is solved independently once by FEM, and data is collected at the observation points. The search space can always be restricted to the prior known range of the dielectric constant of a structure. As an example, if the approach is applied to a tooth structure, and its relative permittivity range is known for tooth enamel and dentine, the search space is bounded within this range (between 1.0 and 10.0). Searching over the permittivity range enables us to distinguish between the tooth enamel and dentine subregions in this particular example. Figure 4(a) presents the proposed objective function based on the non-radiating VECS for the object with two subregions with the test permittivities ranging from 1.25 to 9.00 (n' = 961). The top (x-y) view of the objective functions is shown in Figure 4(b). The Born iterative method's objective function is converted to a linear one by initializing one of the unknowns, the internal total electric field. However, there is no such simplification in the proposed objective function, and the objective function remains non-linear. The non-radiating objective function is minimal when the estimated permittivity of the two regions is 1.25, which is the correct value.

Secondly, we investigated the performance of the proposed objective function for estimating the permittivity profile of a high-contrast object. As shown in Figure 3(b), the high-contrast scatterer is illuminated by a plane wave. The non-radiating objective function for the object with two subregions with the permittivities ranging from 1.25 to 9.00 (n' = 961) is presented in Figure 5. As shown in Figures 5(a) and (b), the non-radiating objective function and the non-radiating approach enable us to estimate the relative permittivity of the [v.sup.1.sub.2] and [v.sup.2.sub.2] regions as 6.00, which is the correct value.

Next, we investigated the performance of the proposed nonradiating objective function for estimating the permittivity profile of an object including two subregions with different permittivity values. The permittivity profile of the scatterer under test is shown in Figure 3(c). As illustrated in Figures 6(a) and (b), the non-radiating objective function has distinctive minima when the relative permittivities of the [v.sup.1.sub.2] and [v.sup.2.sub.2] regions are, respectively, 3.00 and 6.00, which are the correct values.

In the second case study, the performance of the proposed approach is evaluated for the object including more than two homogenous subregions. In this case, we considered the EMISS geometrical configuration of Figure 2(b). The object under test, [v.sub.2], includes three homogenous subregions: [v.sup.1.sub.2], [v.sup.2.sub.2] and [v.sup.3.sub.2] in a free space background. The first subregion is a circular cylindrical dielectric with a diameter of one wavelength centered at (0.6[lambda], 1.8[lambda]). The second subregion is a rectangular cylindrical dielectric with one wavelength side centered at (1.5[lambda], -1.5[lambda]). The third subregion is a rectangular cylindrical dielectric with a half wavelength width and a wavelength length centered at (-1.8[lambda], 0). The dielectric scatterers are considered to be lossless and are subdivided into 5046 elements.

We investigated the performance of the proposed non-radiating objective function for estimating the permittivity profile of a high-contrast object. The high-contrast scatterer's permittivity profile under test is shown in Figure 7(a). The non-radiating objective function for the object with three subregions, with permittivities ranging from 1.25 to 9.00 (n' = 729) are presented in Figure 7(b). The non-radiating objective function is minimal when the relative permittivity of regions [v.sup.1.sub.2], [v.sup.2.sub.2] and [v.sup.3.sub.2] is 6.00, which is the correct value.

Then, we investigated the performance of the proposed approach for estimating the permittivity profile of an object including three subregions with different permittivity values as shown in Figure 8(a). The non-radiating objective function for the object with three subregions is presented in Figure 8(b). The non-radiating objective function is at minimum when the permittivity profile amplitudes of regions [v.sup.1.sub.2], [v.sup.2.sub.2] and [v.sup.3.sub.2] are, respectively, 3.00, 5.00 and 8.00, which are the correct values.

So far, the proposed approach has been verified for the noise free data. To evaluate the proposed approach performance with the noisy data, the high-contrast object with the permittivity profile shown in Figure 3(b) is illuminated with a plane wave, and a white Gaussian noise is added to the FEM simulation results. The non-radiating objective function is evaluated by sweeping the relative permittivity with the 0.25 step size between 1 and 9 for the following signal to noise ratio (SNR): 60db, 40db, 20db, 10db, 5db, 3db, 2db, and 1db. The results indicate that the minimum values for the non-radiating objective function occur at the correct relative permittivity ([[epsilon].sub.r] = 6) for the 60 db, 40 db, 20 db, and 10 db SNRs, but it deviates from the true relative permittivity when the SNR drops below 10 db, as shown in Figure 9(a). Thus, the permittivity profile estimation based on the proposed approach is also valid for noisy measured data. The permittivity profile estimation error for the 5db, 3db, and 2db SNRs is %8, and for the 1db SNR is %14, as shown in Figure 9(b).

The simulation results indicate that the proposed approach can be successfully applied to both low-contrast and high-contrast permittivity profiles of a large object (in terms of wavelength). The simulation results illustrate that the proposed approach can correctly estimate the permittivity profile of the object under test in a noisy measurement environment.

6. CONCLUSION

A new solution to the permittivity profile estimation problem is presented for a scatterer whose permittivity profile can be divided into a number of homogenous clusters. To address the Born approximation initial guess issue pointed out in [15] for the electromagnetic inverse scattering problem, we have introduced the non-radiating contrast factor through which an objective function for permittivity profile estimation is formulated. The solutions to the inverse source problem are non-unique. However, our method yields the correct and unique permittivity profile of an unknown object by minimizing the non-radiating objective function and applying the boundary information. The conventional Born's approximation works well for a low-contrast object as long as the object size is a fraction of a wavelength. In contrast, the proposed non-radiating objective function can be used for both low-contrast and high-contrast permittivity profiles with even large size objects. The method has been verified by extensive simulations. The proposed approach has a much smaller number of unknowns and is therefore computationally more efficient than the permittivity profile estimation approaches based on Born approximation. The simulation results also depict that the proposed approach can accurately estimate the permittivity profile of an object under test in a noisy environment.

ACKNOWLEDGMENT

This work has been supported by NSERC (Natural Sciences and Engineering Research Council) of Canada and BlackBerry (former Research In Motion).

The author would also like to thank Prof. J. Orchard, Dr. D. Saeedkia, Dr. M. Neshat, Dr. D. Hailu, and Ms. N. Ranjkesh for the valuable technical discussions, and Ms. M. McPherson for her invaluable help in editing this paper.

REFERENCES

[1.] Autierin, R., M. Durso, C. Eyraud, A. Litman, V. Pascazio, and T. Isernia, "3d inversion of fresneldatabase with 3d Markov random field," Session 3P6 Inverse Scattering Probelms: Open Problems and New Challenges, 559, 2010.

[2.] Berry, E., A. J. Fitzgerald, N. N. Zinov'ev, G. C. Walker, S. Homer-Vanniasinkam, C. D. Sudworth, R. E. Miles, J. M. Chamberlain, and M. A. Smith, "Optical properties of tissue measured using terahertz pulsed imaging," Proceedings of SPIE, Vol. 5030, 459-470, 2003.

[3.] Bertero, M., T. A. Poggio, and V. Torre, "Ill-posed problems in early vision," Proceedings of the IEEE, Vol. 76, No. 8, 869-889, 1988.

[4.] Bojarski, N. N., "Inverse scattering," Technical report, DTIC Document, 1973.

[5.] Caorsi, S., G. Gragnani, M. Pastorino, G. Zunino, et al, "Microwave imaging based on a Markov random field model," IEEE Transactions on Antennas and Propagation, Vol, 42, No. 3, 293-303, 1994.

[6.] Chen, X., "Application of signal-subspace and optimization method in reconstructing extended scatterers," JOSA A, Vol. 26, 1022-1026, 2009.

[7.] Chen, X., "Subspace-based optimization method for solving inverse scattering problems," IEEE Trans. Geosci. Remote Sens, Vol. 48, No. 1, 42-49, 2010.

[8.] Chew, W. C. and Y. M. Wang, "Reconstruction of two-dimensional permittivity distribution using the distorted born iterative method," IEEE Transactions on Medical Imaging, Vol. 9, No. 2, 218-225, 1990.

[9.] Chew, W. C., Waves and Fields in Inhomogeneous Media, IEEE Press, 1995.

[10.] Chiappe, M. and G. L. Gragnani, "Analytical solution to inverse electromagnetic scattering: Shape and position reconstruction of dielectric objects," IEEE International Workshop on Imaging Systems and Techniques, 1-6, 2007.

[11.] Crawley, D., C. Longbottomm, V. P. Wallace, B. Cole, D. Arnone, and M. Pepper, "Three-dimensional terahertz pulse imaging of dental tissue," Journal of Biomedical Optics, Vol. 8, No. 2, 303-307, 2003.

[12.] Cui, T. J., W. C. Chew, A. A. Aydiner, and S. Chen, "Inverse scattering of two-dimensional dielectric objects buried in a lossy earth using the distorted born iterative method," IEEE Trans. Geosci. Remote Sens., Vol. 39, No. 2, 339-346, 2001.

[13.] Fieguth, P., Statistical Image Processing and Multidimensional Modeling, Vol. 155, Spring, 2010

[14.] Scaf, C. E. S. G., J. Morais, and L. C. M. Loffredo, "A survey of radiographic measurement estimation in assessment of dental implant length," Oral Surgery, Oral Medicine, Oral Pathology, Oral, Radiology and Endodontology, Vol. 103, No. 53, 2007.

[15.] Gilmore, C., P. Mojabi, A. Zakaria, M. Ostadrahimi, C. Kaye, S. Noghanian, L. Shafai, S. Pistorius, and J. LoVetri, "A wideband microwave tomography system with a novel frequency selection procedure," IEEE Transactions on Biomedical Engineering, Vol. 57, No. 4, 894-903, 2010.

[16.] Gilmore, C., P. Mojabi, A. Zakaria, S. Pistorius, and J. LoVetri, "On super-resolution with an experimental microwave tomography system," IEEE Antennas and Wireless Propagation Letters, Vol. 9, 393-396, 2010.

[17.] Gragnani, G., "Electromagnetic imaging using closed-form radiating and non-radiating currents," Mediterranean Microwave Symposium, 1-6, November 2009.

[18.] Gragnani, G., "Two-dimensional nonradiating currents for imaging systems: Theoretical development and preliminary assessment," IET Microwaves, Antennas Propagation, Vol. 3, No. 8, 1164-1171, 2009.

[19.] Gragnani, G. and M. D. Mendez, "An improved electromagnetic imaging procedure using non-radiating sources," Mediterranean Microwave Symposium, 188-191, 2010.

[20.] Hajihashemi, M. R. and M. El-Shenawee, "Inverse scattering of three-dimensional PEC objects using the level-set method," Progress In Electromagnetics Research, Vol. 116, 23-47, 2011.

[21.] Harrington, R. F., Time-Harmonic Electromagnetic Fields, The IEEE Press Series, Pscataway, NJ, 2001.

[22.] Johnson, S., Y. Tae-Hoon, and R. Jung-Woong, "Inverse scattering solutions of scalar Helmholtz wave equation by a multiple source moment method," Electron. Lett., Vol. 19, 130-132, 1983.

[23.] Lavarello, R. J. and M. L. Oelze, "Tomographic reconstruction of three-dimensional volumes using the distorted born iterative method," IEEE Transactions on Medical Imaging, Vol. 28, No. 10, 1643-1653, 2009.

[24.] Lee, K. C., J.-S. Ou, and M. C. Fang, "Application of SVD noise-reduction technique to PCA based radar target recognition," Progress In Electromagnetics Research, Vol. 81, 447-459, 2008.

[25.] Markel, V. A., V. Mital, and J. C. Schotland, "Inverse problem in optical diffusion tomography iii inversion formulas and singular value decomposition," J. Opt. Soc. Am. A, Vol. 20, No. 5, 890-902, 2003.

[26.] Muller, C., "Electromagnetic radiation patterns and sources," IRE Transactions on Antennas and Propagation, Vol. 4, No. 3, 224-232, 1956.

[27.] Ney, M. M., A. M. Smith, and S. S. Stuchly, "A solution of electromagnetic imaging using pseudoinverse transformation," IEEE Transactions on Medical Imaging, Vol. 3, No. 4, 1984.

[28.] Nguyen, K. L., M. L. Johns, and L. F. Gladden, "Three dimensional imaging with a terahertz quantum cascade laser," Optics Express, Vol. 14, No. 6, 2006.

[29.] Porter, R. and A. Devaney, "Holography and the inverse source problem," Journal of Optical Society of America, Vol. 72, 327-330, 1982.

[30.] Porter, R. P., "Diffraction-limited, scalar image formation with holograms of arbitrary shape," Journal of the Optical Society of America, Vol. 60, No. 8, 1051-1059, 1970.

[31.] Richmond, J. H., "Scattering by a dielectric cylinder of arbitrary cross section shape," IEEE Antennas and Wireless Propagation, Vol. 13, No. 3, 334-341, 1965.

[32.] Tachibana, H. and K. Matsumoto, "Applicability of x-ray computerized tomography in endodontics," Dental Traumatology, Vol. 6, No. 1, 16-20, 1990.

[33.] Wang, Y. and W. Chew, "An iterative solution of the two-dimensional electromagnetic inverse scattering problem," International Journal of Imaging Systems and Technology, Vol. 1, No. 1, 100-108, 1989.

[34.] Winters, D. W., B. D. Van Veen, and S. C. Hagness, "Estimation of the frequency-dependent average dielectric properties of breast tissue using a time-domain inverse scattering technique," IEEE Transactions on Antennas and Propagation, Vol. 54, 3517-3528, 2006.

[35.] Winters, D. W., B. D. Van Veen, and S. C. Hagness, "A sparsity regularization approach to the electromagnetic inverse scattering problem. IEEE Transactions on Antennas Propagation, Vol. 58, 145-154, 2010.

[36.] Wirgin, A., "The inverse crime," ArXiv preprint mathph/0401050, 2004.

[37.] Wolf, E., "Three-dimensional structure determination of semitransparent objects from holographic data," Optical Communication, Vol. 1, No. 4, 153-156, 1969.

Shahed Shahir *, Mehrbod Mohajer, Arash Rohani, and Safieddin Safavi-Naeini

Department of Electrical and Computer Engineering, University of Waterloo, Waterloo, Ontario N2L 3G1, Canada

Received 23 January 2013, Accepted 28 March 2013, Scheduled 2 April 2013

* Corresponding author: Shahed Shahir (sshahir@maxwell.uwaterloo.ca).

Printer friendly Cite/link Email Feedback | |

Author: | Shahir, Shahed; Mohajer, Mehrbod; Rohani, Arash; Safavi-Naeini, Safieddin |
---|---|

Publication: | Progress In Electromagnetics Research B |

Article Type: | Report |

Geographic Code: | 1CANA |

Date: | May 1, 2013 |

Words: | 5568 |

Previous Article: | Detailed study of millimeter wave EBG guide: broadbanding techniques, modal structure, and crosstalk behavior. |

Next Article: | On performance of high-efficiency ferrite meander antenna (HEMA) for MIMO communications. |

Topics: |