**High to Microwave Frequencies Imaging Techniques**

George A. Kyriacou1, Ilias N. Aitidis1, Dimitrios G. Drogoudis1 and John N. Sahalos2 *1Department of Electrical and Computer Engineering, Democritus University of Thrace 2Department of Electrical & Computer Engineering, University of Nicosia 1Greece 2Cyprus*

#### **1. Introduction**

146 Medical Imaging

Laurila, J., Standertskjöld-Nordenstam, C-G., Suramo, I, Tolppanen, E-M., Tervonen, O.,

Law, A. (2007) Simulation Modeling and Analysis, Fourth Edition. New York. USA,

Monsef, Y (1997), "Modelling and simulation of complex systems," Society for computer

NAE/IOM (2005), "Building a Better Delivery System", National Academy of Engineering

Neriz, L., Ramis, F. and Sepulveda, J. (2011). "A New Approach for Healthcare Simulation",

Poole, K., Hinton, J. and Kraebber K. (2010). "The gradual leaning of health systems".

Ramis, F., Neriz, L. and Sepúlveda, J. (2008). "A Simulator to Improve Waiting Times at an

Ramis, F., Concha, P. Neriz, L. and Sepúlveda, J. (2009). "Improving Services of a Hospital

Suthummanon, S., Omachonu, V. and Akcin.M. (2005)."Applying activity-based costing to nuclear medicine unit". Health Services Magaement Research, 18, 141-150. Wideman, C. and Gallet, J. (2006) "Analog to digital workflow improvement: A quantitative

Young, T., Eatock, J., Jahangirian, M., Naseer, A. and Lilford, R. (2009)."Three critical

Emergency Unit". Industrial Engineering Research Conference, Vancouver,

Imaging Center using ABC Costing, Lean and Simulation" .Simulation Solutions,

challenges for modeling and simulation in healthcare". In Proceedings of the 2009

Proceedings Society of Health systems, SHS 2011, Orlando, USA

of the SDPS JUNE 2006, Vol. 10, No. 2, pp. 1-19

studt". Journal of Digital Imaging, 19 Suppl 1:29-34. Womack, P. and D. T. Jones.(2005). Soluciones Lean. Gestión2000, Spain.

and Institute of Medicine, USA.

Industrial Engineer. April, 50-55.

IIE Annual Conference.

Winter Simulation Conference. Flexsim user manual (2011), www.flexsim.com.

96-100.

Canada.

McGraw-Hill

Korhola, O. and Brommels, M. (2001)."The Efficacy of a continuous quality improvement (CQI) method in a radiological department". Acta Radiologica 42,

simulation international. In. Jin, X., Kagioglou M., Aouad G. 2006, Towards a dynamic healthcare process: from requirement capture to simulation. Transactions

> Microwave and high frequency tomography constitutes challenging electromagnetic inverse problems aiming at the reconstruction of its internal σ-, εr- and/or μr- distributions. The object to be imaged is embedded in a lossy homogeneous background medium. This is surrounded by a fictitious (or real) surface preferably of canonical shape (circular, rectangular, cylindrical or spherical) over which a number of antennas (electrodes at lower frequencies) are evenly distributed. The hardware implementing the modality should be able to successively activate each one of them, while setting all the other antennas to a receive state operating as sensors. Instead of requiring all antennas to operate in both transmit and receive modes, a subset of antennas can be configured in transmit mode only (activated in turn) while a preferably larger subset is configured in receive only mode (passive sensors), all of them performing simultaneous measurements. In this manner the object is illuminated each time from a different angle creating a scattered field to all possible directions, which is sampled by the receiving antennas. The locations and number of transmit antennas should be designed to cover the required spatial illuminations (projection angles), while the number of the evenly distributed receiving antennas should fulfill the spatial sampling laws. An alternative configuration mimicking the one already used in X-Ray Computer Tomography (CT) and Magnetic Resonance Imaging (MRI) seems preferable and highly recommended for the microwave imaging. Spesifically, for two-dimensional imaging a single active antenna along with an array of sensing-passive antennas could be placed on a circular-colar possibly plastic platform surrounding the object to be imaged. In turn by rotating this supporting structure the object can be illuminated from all possible projection angles. For a three-dimensinal imaging the hosting platform could be a plastic cylinder holding multiple antenna ring arrays, each ring having one active and multiple passive antennas. Each antenna is activated successively in time while the whole cylinder is roteated providing all possible illuminations and data recording-sampling by all sensing antennas "simultaneously". The information gathered by this measurement procedure constitutes the dataset to be exploited by the imaging algorithm to reconstruct the object's internal properties distributions. This constitutes a challenging mathematical-computational and engineering problem since it is proved to be a highly nonlinear and ill-posed inverse problem.

High to Microwave Frequencies Imaging Techniques 149

At this point two serious problems are introduced related to the inherent properties of the "sensitivity matrix". Firstly, its entries depend on the unknown properties (σi, εri, μri), hence the system to be solved is a non-linear one. Secondly, this sensitivity matrix is usually an illposed one or a singular matrix. This latter means that if a singular value decomposition is performed (not eigenvalue since this is a rectangular MxN matrix, where the number of measurements M should be much greater than the number (N) of unknowns as M>>N in order to confront measurement errors, noise and ill-posedeness) the maximum singular value appears a few orders of magnitude larger than the minimum one. This inherent property is closely related to the measurement (generally the data collection) setup and can be significantly improved by intuitive techniques. Returning to the nonlinearity its direct consequence is that the resulting parameters (σi, εri, μri) provided by the Minimization first iteration are not the true ones but only a better approximation, if the reconstruction process was successful. Hence, the whole procedure should be repeated again and again until the difference between the measured and the calculated dataset becomes comparable to the expected measurement error and/or the required convergence is achieved. Obviously, at

each iteration a new sensitivity matrix and a new data set is evaluated and used.

toward RF-Microwave imaging practical implementation.

for both the forward and the inverse problem solution.

constituting future research areas will be suggested.

**2. Forward problem & inverse problem characteristics**

In the following sections the reader will be introduced to the employed approaches from both the mathematical-computational as well as the electromagnetics point of view. But mostly the open research challenges will be pointed out motivating new research and paths

The forward and inverse problems general characteristics are discussed in the second section, after the introduction to high and microwave frequencies imaging modalities. The procedure to formulate and numerically solve the forward problem is given in the third section, for both the high (MHz) and microwave regimes. Within the foarth section the cost function is first setup and the Perturbation Methodology for its direct iterative solution is then introduced step-by-step from static to microwave imaging. The fifth section elaborates on the establishment of the "sensitivity Equation" based on the "Adjoint Network Theorem" for the microwave band and its equivalent "Electrical Networks Compensation Theorem" for the MHz range. Before presenting the forward and inverse problem numerical results (sixth section) the importance of the study of the sensitivity or Jacobian matrix is pointed out. This is a very promising research area especially when modern Principal Component Analysis (PCA) or its counterpart Proper Orthogonal Decomposition (POD) approaches are employed. Either PCA or POD are based on a Singular Value Decomposition of the sensitivity matrix (rectangular matrix) by an algebraic manipulation of its eigenvectors. These techniques present prospects for novel and computationally efficient methodologies

The last section is devoted to numerical results. A series of successfully reconstructed conductivity and permittivity distributions for both the MHz and the microwave regimes will be presented. The algorithm performance will be discussed and possible improvements

The first step in the course of establishing an imaging modality refers to the construction of the appropriate computer model, which should reflect the practical geometry as closely as possible. This model serves in twofold, first it should enable an approach which mimics the

Analytical methods can be employed only for simple canonical geometries but these are valuable since they may serve as exact reference methods to validate numerical techniques. Hence, for the practical arbitrary shaped and inhomogeneous bodies numerical techniques are inevitable. It is on these approaches that the following chapter is elaborating.

For the reconstruction algorithm to be implemented a realistic as far as possible computer model of the practical structure is necessary. For this purpose the geometry is discretized including the appropriate antenna (or electrodes) modeling by following an engineering compromise between accurately reflecting the structure and the required computational resources. A usual approach is to consider a "virtual body area" of canonical shape large enough to contain the unknown irregularly shaped actual object with its unknown σ, εr, μr distributions. The properties of the background media outside this virtual surface are assumed known as it is practically selected as desired. The also virtual surface carrying the transmit/receive antennas is larger than the "virtual body area" but also embedded in the known background medium. Theoretically, this background media extends to infinity, however in order to implement a numerical technique a finite solution domain should be established. Hence, another fictitious surface enclosing the whole structure constituting the solution domain "truncation surface" is considered. In turn this truncation surface must not disturb the electromagnetic field solution or to be "transparent" to the impinging waves. The numerically discretized model includes everything within the truncation surface, but the inverse problem unknowns are only the constitutive parameters within the "virtual body area". By the aid of the discretization this "virtual body" is comprised of a number of pixels for two-dimensional (2-D) or voxels for three-dimensional (3-D) geometries. Each one (ith) of them is assumed locally homogeneous with constant but unknown (σi, εri, μri) properties which are in principle different for each pixel/voxel, forming the so-called piecewise constant distribution. Now, the aim of the reconstruction is exactly the evaluation of these three vectors [σ]=[σi], [εr]=[εri] and/or [μr]=[μri], by exploiting the dataset acquired from field measurements carried out on the real object. Numerous different approaches are established for the exploitation of this information. The most usual approach which is also elaborated herein is to formulate a cost function based on the least square method, which will be in turn minimized employing some type of optimization techniques. For this purpose an initial (σi o, εrio, μrio) distribution, usually simply a homogeneous one is assumed and the measurement procedure is mimicked through computer simulations. Namely, for each active antenna the **forward problem** corresponding to the solution of a vector wave equation or in general the numerical solution of Maxwell's equations, yields the field distribution all over the receiving antennas. Exactly the field calculated on the receiving antennas, when gathered from all illuminating active antennas, it setup a calculated dataset. On the other hand the field calculated all over the structure and particularly over the "virtual body area" is exploited within the methodology elaborated herein for the evaluation of a **"Sensitivity" or "Jacobian" matrix.** This is obtained through a closed form sensitivity equation established through a combination of an "Adjoint Network Theorem" following the Electromagnetics Reciprocity Theorem approach. Its entries exactly reflect the sensitivity of the calculated field at each sensing antenna with respect to a differential change of each pixel/voxel unknown parameter. With the availability of this information an algorithm minimizing the differences between measured and calculated fields (the complete datasets) based on a least square approach is established, which herein is efficiently implemented exploiting the sensitivity matrix in its closed form expression.

Analytical methods can be employed only for simple canonical geometries but these are valuable since they may serve as exact reference methods to validate numerical techniques. Hence, for the practical arbitrary shaped and inhomogeneous bodies numerical techniques

For the reconstruction algorithm to be implemented a realistic as far as possible computer model of the practical structure is necessary. For this purpose the geometry is discretized including the appropriate antenna (or electrodes) modeling by following an engineering compromise between accurately reflecting the structure and the required computational resources. A usual approach is to consider a "virtual body area" of canonical shape large enough to contain the unknown irregularly shaped actual object with its unknown σ, εr, μr distributions. The properties of the background media outside this virtual surface are assumed known as it is practically selected as desired. The also virtual surface carrying the transmit/receive antennas is larger than the "virtual body area" but also embedded in the known background medium. Theoretically, this background media extends to infinity, however in order to implement a numerical technique a finite solution domain should be established. Hence, another fictitious surface enclosing the whole structure constituting the solution domain "truncation surface" is considered. In turn this truncation surface must not disturb the electromagnetic field solution or to be "transparent" to the impinging waves. The numerically discretized model includes everything within the truncation surface, but the inverse problem unknowns are only the constitutive parameters within the "virtual body area". By the aid of the discretization this "virtual body" is comprised of a number of pixels for two-dimensional (2-D) or voxels for three-dimensional (3-D) geometries. Each one (ith) of them is assumed locally homogeneous with constant but unknown (σi, εri, μri) properties which are in principle different for each pixel/voxel, forming the so-called piecewise constant distribution. Now, the aim of the reconstruction is exactly the evaluation of these three vectors [σ]=[σi], [εr]=[εri] and/or [μr]=[μri], by exploiting the dataset acquired from field measurements carried out on the real object. Numerous different approaches are established for the exploitation of this information. The most usual approach which is also elaborated herein is to formulate a cost function based on the least square method, which will be in turn minimized employing some type of optimization techniques. For this

o, εrio, μrio) distribution, usually simply a homogeneous one is assumed

and the measurement procedure is mimicked through computer simulations. Namely, for each active antenna the **forward problem** corresponding to the solution of a vector wave equation or in general the numerical solution of Maxwell's equations, yields the field distribution all over the receiving antennas. Exactly the field calculated on the receiving antennas, when gathered from all illuminating active antennas, it setup a calculated dataset. On the other hand the field calculated all over the structure and particularly over the "virtual body area" is exploited within the methodology elaborated herein for the evaluation of a **"Sensitivity" or "Jacobian" matrix.** This is obtained through a closed form sensitivity equation established through a combination of an "Adjoint Network Theorem" following the Electromagnetics Reciprocity Theorem approach. Its entries exactly reflect the sensitivity of the calculated field at each sensing antenna with respect to a differential change of each pixel/voxel unknown parameter. With the availability of this information an algorithm minimizing the differences between measured and calculated fields (the complete datasets) based on a least square approach is established, which herein is efficiently

implemented exploiting the sensitivity matrix in its closed form expression.

are inevitable. It is on these approaches that the following chapter is elaborating.

purpose an initial (σi

At this point two serious problems are introduced related to the inherent properties of the "sensitivity matrix". Firstly, its entries depend on the unknown properties (σi, εri, μri), hence the system to be solved is a non-linear one. Secondly, this sensitivity matrix is usually an illposed one or a singular matrix. This latter means that if a singular value decomposition is performed (not eigenvalue since this is a rectangular MxN matrix, where the number of measurements M should be much greater than the number (N) of unknowns as M>>N in order to confront measurement errors, noise and ill-posedeness) the maximum singular value appears a few orders of magnitude larger than the minimum one. This inherent property is closely related to the measurement (generally the data collection) setup and can be significantly improved by intuitive techniques. Returning to the nonlinearity its direct consequence is that the resulting parameters (σi, εri, μri) provided by the Minimization first iteration are not the true ones but only a better approximation, if the reconstruction process was successful. Hence, the whole procedure should be repeated again and again until the difference between the measured and the calculated dataset becomes comparable to the expected measurement error and/or the required convergence is achieved. Obviously, at each iteration a new sensitivity matrix and a new data set is evaluated and used.

In the following sections the reader will be introduced to the employed approaches from both the mathematical-computational as well as the electromagnetics point of view. But mostly the open research challenges will be pointed out motivating new research and paths toward RF-Microwave imaging practical implementation.

The forward and inverse problems general characteristics are discussed in the second section, after the introduction to high and microwave frequencies imaging modalities. The procedure to formulate and numerically solve the forward problem is given in the third section, for both the high (MHz) and microwave regimes. Within the foarth section the cost function is first setup and the Perturbation Methodology for its direct iterative solution is then introduced step-by-step from static to microwave imaging. The fifth section elaborates on the establishment of the "sensitivity Equation" based on the "Adjoint Network Theorem" for the microwave band and its equivalent "Electrical Networks Compensation Theorem" for the MHz range. Before presenting the forward and inverse problem numerical results (sixth section) the importance of the study of the sensitivity or Jacobian matrix is pointed out. This is a very promising research area especially when modern Principal Component Analysis (PCA) or its counterpart Proper Orthogonal Decomposition (POD) approaches are employed. Either PCA or POD are based on a Singular Value Decomposition of the sensitivity matrix (rectangular matrix) by an algebraic manipulation of its eigenvectors. These techniques present prospects for novel and computationally efficient methodologies for both the forward and the inverse problem solution.

The last section is devoted to numerical results. A series of successfully reconstructed conductivity and permittivity distributions for both the MHz and the microwave regimes will be presented. The algorithm performance will be discussed and possible improvements constituting future research areas will be suggested.

## **2. Forward problem & inverse problem characteristics**

The first step in the course of establishing an imaging modality refers to the construction of the appropriate computer model, which should reflect the practical geometry as closely as possible. This model serves in twofold, first it should enable an approach which mimics the

High to Microwave Frequencies Imaging Techniques 151

matrix (as [J][J] T). As will be explained latter the severity of ill-conditioning is defined by the ratio of the larger to the smaller singular value of the Jacobian matrix (singular instead of eigenvalue since the Jacobian matrix is rectangular due to the requirement of a number of

Summarizing the above, the imaging modality elaborated herein constitutes a highly nonlinear and ill-conditioned (singular) inverse problem. The next issue to be discussed refers to

Theoretically an object can be represented by a two-dimensional model (2-D) if it is homogeneous and uniformly extends to "infinity" along a direction perpendicular to the modeled cross-section, which can be arbitrary shaped and inhomogeneous. The question is now, under what assumptions the current density, the electric potential or electric field intensity of a practical three-dimensional (3-D) object can be replicated over one of its crosssections (2-D model)? Once again recall that this approach is very well suited for x-rays CT mainly because x-rays travel along straight lines and hence their measured intensity involves information relevant to the inhomogeneities (varying tissue) along the straight line from the transmiter to the receiving sensor. Thus by attaching an array of sensors in the same plane with the source as shown in Fig.1a, all measured intensities depend on the specific enclosed cross-section. Besides that the specific source location yields an illumination from a corresponding angle of view, called a "projection angle" from physical optics. It is, thus, found very convenient to locate the source and the receivers-sensors equidistantly along a circular bar (a plastic collar). Rotating this circular structure results to an illumination of the imaged cross-section from any possible angle of view, to be defined at

A similar configuration could be in principle setup for the electrical impedance (EIT) or microwave tomography by planning an array of electrodes or antennas on a single plane

Focusing on EIT or MHz tomography the current density injected through the driven-active electrode (source) is curved around the objects of lower conductivity (Fig.1), but the same phenomenon occurs along the third dimension perpendicular to the cross-section under study. Consequently, the injected current density which offers the means to extract the information regarding the conductivity distribution to be imaged, cannot be restricted to flow only across the electrodes plane. Instead this current is spread in all directions around the active electrode and respectively it is collimated from all directions toward the current-sink active electrode. Thus, the coplanar voltage-sensing electrodes yield measurements infected by conductivity inhomogeneities above or below the studied cross-section. Additionally, it is impossible to estimate the total current flowing across the electrodes cross section. To explain the difficulty, assume that a 5mA current (I=5mA) is injected during the actual EIT measurements, the question is what is the current value to be applied to the active electrode in the 2-D model? This should be only the fraction of the 5mA actually flown through the electrode cross-section, otherwise by assuming a 5mA value yields overestimated calculated voltages at the electrodes. Besides this problem, by restricting the current to flow in a single plane in the 2-D model, its density becomes higher within the highly conducting areas. This phenomenon results to higher sensitivity related to these areas and causes the Jacobian matrix to be fictitiously more

measurements higher than the number of unknown σ, εr degrees of freedom).

the assumptions involved in a two-dimensional approximation.

**2.2 Two-versus three-dimensional modeling**

the desired number enabling the image reconstruction.

around the desired cross-section as shown in Fig.1b and 1c respectively.

measurement procedure and secondly it should offer the ability to represent the unknown distribution (σ, εr or both) as a set of successively improved-updated variables. Mimicking the measurement procedure refers to the solution of the governing differential equations, identified as a generalized Laplace or Poisson equation in MHz range, while the full wave Maxwell equations must be considered in the microwave regime. The practical bodies of interest are of arbitrary shape and present a complicated inhomogeneous internal structure in both σ and εr.

In general, these objects should be represented as three dimensional (3-D) models, however in this case both the forward and the inverse problems become very complicated with high computational demands. Besides these, the data collection or measurement strategy is also complex. The possibility of a two-dimensional modeling even an approximate one could be very convenient, since it could simplify the forward problem solution, restricting the data collection approaches to just a few and overally resulting to a straightforward reconstruction algorithm. The question is when the involved approximations are acceptable and whether these introduce any fictitious mathematical or numerical complications? Let as elaborate next on the issue of a 3-D versus a 2-D modeling.

#### **2.1 Non-linearity and Ill-conditioning of the inverse problem**

The issue of two-dimensional modeling is not rhetorical, since it has been extensively exploited in most of the established imaging modalities like x-rays computerized tomography (CT). However, there is a significant difference between x-rays and electromagnetic tomography in that x-rays propagate along straight lines, but in contrary electric current and electromagnetic waves in general flow (stream-) lines are curved whenever they intercept an interface where the conductivity or permittivity changes [e.g. from (εr1, σ1) to (εrz, σz)], as it usually happens in biological structures. This causes severe problems not only in worsening the option of a 2-D approximation but also in rendering a non-linear inverse electromagnetic problem. A fundamental reasoning for this difficulty is explained by the fact that current streamlines or wave propagation directions curvature strongly depends on the conductivity (σz/σ1) and permittivity (εz/ε1) contrast at these interfaces. But these contrasts are indeed the unknown quantities to be sought by the electromagnetic imaging algorithm. Namely, when the inverse problem is finally formulated into a discrete system of equations the stiffness matrix elements will depend on the unknown (σ, εr) distributions hence comprising a non-linear system. In general, this phenomenon can be identified as an internal scattering and/or diffraction which is equally present in Acoustical or Ultrasonic imaging.

Besides the non-linearity, the streamlines curvature appearing in bodies with very high σ or εr contrast (e.g. between blood and fat or bone tissues) results to regions with very high current densities (e.g. blood) while at others this is very low (in fat or bone). Similar high variations occur in the electric field intensities mainly due to the high *εr* constrast. This is exactly what causes the high variation in the sensitivity of the measured quantity (voltage or electric field) with respect to the unknown *σ* or *εr* values. To understand this phenomenon assume a tissue (area) which is isolated from the current flowing through the body, this will in turn present a very low sensitivity as it cannot affect the injected current, which will not be able to "see" it. Hence, its σ or εr values cannot be reliably reconstructed. Mathematically, this high variation in the sensitivity over the body to be imaged is depicted as a high singularity or high ill-conditioning in the "Jacobian" matrix [J], which comprises the inverse problem

measurement procedure and secondly it should offer the ability to represent the unknown distribution (σ, εr or both) as a set of successively improved-updated variables. Mimicking the measurement procedure refers to the solution of the governing differential equations, identified as a generalized Laplace or Poisson equation in MHz range, while the full wave Maxwell equations must be considered in the microwave regime. The practical bodies of interest are of arbitrary shape and present a complicated inhomogeneous internal structure

In general, these objects should be represented as three dimensional (3-D) models, however in this case both the forward and the inverse problems become very complicated with high computational demands. Besides these, the data collection or measurement strategy is also complex. The possibility of a two-dimensional modeling even an approximate one could be very convenient, since it could simplify the forward problem solution, restricting the data collection approaches to just a few and overally resulting to a straightforward reconstruction algorithm. The question is when the involved approximations are acceptable and whether these introduce any fictitious mathematical or numerical complications? Let as elaborate

The issue of two-dimensional modeling is not rhetorical, since it has been extensively exploited in most of the established imaging modalities like x-rays computerized tomography (CT). However, there is a significant difference between x-rays and electromagnetic tomography in that x-rays propagate along straight lines, but in contrary electric current and electromagnetic waves in general flow (stream-) lines are curved whenever they intercept an interface where the conductivity or permittivity changes [e.g. from (εr1, σ1) to (εrz, σz)], as it usually happens in biological structures. This causes severe problems not only in worsening the option of a 2-D approximation but also in rendering a non-linear inverse electromagnetic problem. A fundamental reasoning for this difficulty is explained by the fact that current streamlines or wave propagation directions curvature strongly depends on the conductivity (σz/σ1) and permittivity (εz/ε1) contrast at these interfaces. But these contrasts are indeed the unknown quantities to be sought by the electromagnetic imaging algorithm. Namely, when the inverse problem is finally formulated into a discrete system of equations the stiffness matrix elements will depend on the unknown (σ, εr) distributions hence comprising a non-linear system. In general, this phenomenon can be identified as an internal scattering and/or diffraction which

Besides the non-linearity, the streamlines curvature appearing in bodies with very high σ or εr contrast (e.g. between blood and fat or bone tissues) results to regions with very high current densities (e.g. blood) while at others this is very low (in fat or bone). Similar high variations occur in the electric field intensities mainly due to the high *εr* constrast. This is exactly what causes the high variation in the sensitivity of the measured quantity (voltage or electric field) with respect to the unknown *σ* or *εr* values. To understand this phenomenon assume a tissue (area) which is isolated from the current flowing through the body, this will in turn present a very low sensitivity as it cannot affect the injected current, which will not be able to "see" it. Hence, its σ or εr values cannot be reliably reconstructed. Mathematically, this high variation in the sensitivity over the body to be imaged is depicted as a high singularity or high ill-conditioning in the "Jacobian" matrix [J], which comprises the inverse problem

in both σ and εr.

next on the issue of a 3-D versus a 2-D modeling.

is equally present in Acoustical or Ultrasonic imaging.

**2.1 Non-linearity and Ill-conditioning of the inverse problem**

matrix (as [J][J] T). As will be explained latter the severity of ill-conditioning is defined by the ratio of the larger to the smaller singular value of the Jacobian matrix (singular instead of eigenvalue since the Jacobian matrix is rectangular due to the requirement of a number of measurements higher than the number of unknown σ, εr degrees of freedom).

Summarizing the above, the imaging modality elaborated herein constitutes a highly nonlinear and ill-conditioned (singular) inverse problem. The next issue to be discussed refers to the assumptions involved in a two-dimensional approximation.

## **2.2 Two-versus three-dimensional modeling**

Theoretically an object can be represented by a two-dimensional model (2-D) if it is homogeneous and uniformly extends to "infinity" along a direction perpendicular to the modeled cross-section, which can be arbitrary shaped and inhomogeneous. The question is now, under what assumptions the current density, the electric potential or electric field intensity of a practical three-dimensional (3-D) object can be replicated over one of its crosssections (2-D model)? Once again recall that this approach is very well suited for x-rays CT mainly because x-rays travel along straight lines and hence their measured intensity involves information relevant to the inhomogeneities (varying tissue) along the straight line from the transmiter to the receiving sensor. Thus by attaching an array of sensors in the same plane with the source as shown in Fig.1a, all measured intensities depend on the specific enclosed cross-section. Besides that the specific source location yields an illumination from a corresponding angle of view, called a "projection angle" from physical optics. It is, thus, found very convenient to locate the source and the receivers-sensors equidistantly along a circular bar (a plastic collar). Rotating this circular structure results to an illumination of the imaged cross-section from any possible angle of view, to be defined at the desired number enabling the image reconstruction.

A similar configuration could be in principle setup for the electrical impedance (EIT) or microwave tomography by planning an array of electrodes or antennas on a single plane around the desired cross-section as shown in Fig.1b and 1c respectively.

Focusing on EIT or MHz tomography the current density injected through the driven-active electrode (source) is curved around the objects of lower conductivity (Fig.1), but the same phenomenon occurs along the third dimension perpendicular to the cross-section under study. Consequently, the injected current density which offers the means to extract the information regarding the conductivity distribution to be imaged, cannot be restricted to flow only across the electrodes plane. Instead this current is spread in all directions around the active electrode and respectively it is collimated from all directions toward the current-sink active electrode. Thus, the coplanar voltage-sensing electrodes yield measurements infected by conductivity inhomogeneities above or below the studied cross-section. Additionally, it is impossible to estimate the total current flowing across the electrodes cross section. To explain the difficulty, assume that a 5mA current (I=5mA) is injected during the actual EIT measurements, the question is what is the current value to be applied to the active electrode in the 2-D model? This should be only the fraction of the 5mA actually flown through the electrode cross-section, otherwise by assuming a 5mA value yields overestimated calculated voltages at the electrodes. Besides this problem, by restricting the current to flow in a single plane in the 2-D model, its density becomes higher within the highly conducting areas. This phenomenon results to higher sensitivity related to these areas and causes the Jacobian matrix to be fictitiously more

High to Microwave Frequencies Imaging Techniques 153

A more sophisticated trend toward singularity reduction or as sometimes assumed "toward optimal imaging setup" aims at the establishment of the optimum current density distribution for EIT or in the MHz range and respectively "optimal electric field intensity" distribution. As explained above all the involved complications stem from the current or field streamline curvature which depends on the unknown (σ, εr) discontinuities. But it is intuitively expected for EIT that these curvature may become smoother (or streamlines tending toward straight lines) when the active electrodes (both injecting and current sink) is increased. This phenomenon can be diaesthetically understood by realizing that all current (or field) streamlines are emanating from the source or converging toward the sink electrode, traveling toward all possible directions. Thus, their curvature becomes maximum for point electrodes and smooths down as its size is increased. This approach was forced to its limits in the EIT case by the groups of Isaacson [7] and Lionheart [8], where the active electrodes were increased at the limit of almost covering the entire object's surface, retaining only small gaps between them for isolation. This, voltage or electric field sensing electrode or antennas do not presume any large surface and they could retain a small size even almost infinitesimal (point sensing electrodes). Recall that the primary aim of this trend is to reduce the sensitivity or Jacobian matrix singularity, explicitly by reducing the ratio of its larger to the smaller singular value. Conversely means smoothing out the sensitivity over the body surface or volume by forcing more current to flow through low conductivity regions and lowering the current density within high conductivity areas. However, the Authors experience shows that one should not exaggerate in increasing the electrode or antenna surface, but the current should be allowed to follow its natural paths through the object, since this is indeed the mechanism of extracting information from the object interior. As a rule of thumb the electrode or antenna size could be increased only up to the point that the singularity becomes manageable by the reconstruction algorithm, which should withstand relatively high sensitivity variations. In turn the rotation of the active antenna illuminating the structure from all possible angles would ensure a higher overall-total information extraction.

Returning back to the different forward (3-D) and inverse 2-D model, another characteristic should be taken into account, which is related to the calculation accuracy and the desired imaging spatial resolution. Starting from the latter it depends on the available number of linearly independent measurements, which is in turn defined by the number of active and sensing electrodes/antennas, as it will be explained later on. Also, keep in mind that measured values are corrupted by noise and the measurement inaccuracies, including quantization noise introduced during the analog to digital conversion. Thus a reliable inversion scheme requires a quite higher number of measurements (M equal to the number of equations) than the number of unknowns (N) in order to cancel out these inaccuracies through minimization. As a rule of thumb the number of equations should be almost double than that of the unknowns (M2N). Conversely the number of electrodes/antennas corresponds to the achievable spatial voltage or filed intensity sampling, but it is actually defined by the practical setup as the electrode/antenna size and the hosting object dimensions on which these will be attached. Concluding the above parameters define an upper limit in the achievable number of unknowns (N) and hence the offered inverse problem spatial resolution-discretization. Thus the in principle continuous (σ, εr)

**2.3 Dual mesh discretization**

singular than its actual 3-D form. Exactly similar behavior is observed in MHz and microwave imaging where the role of high conductivity regions is also played by the high permittivity areas, especially as the source frequency is increased.

Fig. 1. (a) Current density distribution of an inhomogeneous 2-D model at low frequencies (b) Electric field distribution of an inhomogenous dielectric *(εr)* at microwave frequency.

The above observation regarding the higher singularity of the 2-D model as compared to a 3-D model is very important and should be utilized as a guide toward the establishment of more robust imaging techniques. In particular, most of these complexities stems from the attempt to solve (simulate) the forward problem as a 2-D cross-sectional model, rather than caused from the 2-D inverse problem formulation. Hence, whenever the simplification of the 2-D setup is sought, it is a very good idea to adopt the corresponding electrode/antennas coplanar setup (e.g. Fig.2b, 2c) along with a uniform (σ, εr) distribution in the axial direction, which yields a 2-D inverse problem formulation. However, a finite axial length is considered and a 3-D numerical technique is employed to solve the forward problem. The important benefit of this approach refers to the removal of the deviations between measured and calculated voltages or electric field intensities. Besides that, the fictitious higher sensitivity matrix singularity is reduced to the inherently existing as well.

b) An electrode array for cross-sectional impedance tomography.

c) An antenna array for cross-sectional microwave imaging, that could be setup on a plastic "transparent to electromagnetic waves" collar just like the x-rays CT.

singular than its actual 3-D form. Exactly similar behavior is observed in MHz and microwave imaging where the role of high conductivity regions is also played by the high permittivity

(a) (b)

Fig. 1. (a) Current density distribution of an inhomogeneous 2-D model at low frequencies (b) Electric field distribution of an inhomogenous dielectric *(εr)* at microwave frequency.

The above observation regarding the higher singularity of the 2-D model as compared to a 3-D model is very important and should be utilized as a guide toward the establishment of more robust imaging techniques. In particular, most of these complexities stems from the attempt to solve (simulate) the forward problem as a 2-D cross-sectional model, rather than caused from the 2-D inverse problem formulation. Hence, whenever the simplification of the 2-D setup is sought, it is a very good idea to adopt the corresponding electrode/antennas coplanar setup (e.g. Fig.2b, 2c) along with a uniform (σ, εr) distribution in the axial direction, which yields a 2-D inverse problem formulation. However, a finite axial length is considered and a 3-D numerical technique is employed to solve the forward problem. The important benefit of this approach refers to the removal of the deviations between measured and calculated voltages or electric field intensities. Besides that, the fictitious higher sensitivity

Fig. 2. a) An x-ray setup for cross-sectional imaging, different illumination by rotating the x-

c) An antenna array for cross-sectional microwave imaging, that could be setup on a plastic

matrix singularity is reduced to the inherently existing as well.

b) An electrode array for cross-sectional impedance tomography.

"transparent to electromagnetic waves" collar just like the x-rays CT.

ray platform.

areas, especially as the source frequency is increased.

A more sophisticated trend toward singularity reduction or as sometimes assumed "toward optimal imaging setup" aims at the establishment of the optimum current density distribution for EIT or in the MHz range and respectively "optimal electric field intensity" distribution. As explained above all the involved complications stem from the current or field streamline curvature which depends on the unknown (σ, εr) discontinuities. But it is intuitively expected for EIT that these curvature may become smoother (or streamlines tending toward straight lines) when the active electrodes (both injecting and current sink) is increased. This phenomenon can be diaesthetically understood by realizing that all current (or field) streamlines are emanating from the source or converging toward the sink electrode, traveling toward all possible directions. Thus, their curvature becomes maximum for point electrodes and smooths down as its size is increased. This approach was forced to its limits in the EIT case by the groups of Isaacson [7] and Lionheart [8], where the active electrodes were increased at the limit of almost covering the entire object's surface, retaining only small gaps between them for isolation. This, voltage or electric field sensing electrode or antennas do not presume any large surface and they could retain a small size even almost infinitesimal (point sensing electrodes). Recall that the primary aim of this trend is to reduce the sensitivity or Jacobian matrix singularity, explicitly by reducing the ratio of its larger to the smaller singular value. Conversely means smoothing out the sensitivity over the body surface or volume by forcing more current to flow through low conductivity regions and lowering the current density within high conductivity areas. However, the Authors experience shows that one should not exaggerate in increasing the electrode or antenna surface, but the current should be allowed to follow its natural paths through the object, since this is indeed the mechanism of extracting information from the object interior. As a rule of thumb the electrode or antenna size could be increased only up to the point that the singularity becomes manageable by the reconstruction algorithm, which should withstand relatively high sensitivity variations. In turn the rotation of the active antenna illuminating the structure from all possible angles would ensure a higher overall-total information extraction.

#### **2.3 Dual mesh discretization**

Returning back to the different forward (3-D) and inverse 2-D model, another characteristic should be taken into account, which is related to the calculation accuracy and the desired imaging spatial resolution. Starting from the latter it depends on the available number of linearly independent measurements, which is in turn defined by the number of active and sensing electrodes/antennas, as it will be explained later on. Also, keep in mind that measured values are corrupted by noise and the measurement inaccuracies, including quantization noise introduced during the analog to digital conversion. Thus a reliable inversion scheme requires a quite higher number of measurements (M equal to the number of equations) than the number of unknowns (N) in order to cancel out these inaccuracies through minimization. As a rule of thumb the number of equations should be almost double than that of the unknowns (M2N). Conversely the number of electrodes/antennas corresponds to the achievable spatial voltage or filed intensity sampling, but it is actually defined by the practical setup as the electrode/antenna size and the hosting object dimensions on which these will be attached. Concluding the above parameters define an upper limit in the achievable number of unknowns (N) and hence the offered inverse problem spatial resolution-discretization. Thus the in principle continuous (σ, εr)

High to Microwave Frequencies Imaging Techniques 155

though the forward mesh elements comprising each reconstruction cell have the same constitutive (σi, εri) parameters, this subdivision in necessary in order to fullfill the voltage/field interpolation requirements. Explicitly, in most cases linear interpolation functions are considered within the established numerical techniques (MoM, FEM, FDTD). In turn the field variation especially around fine structures and around conductors are very high, thus asking either for highly non-linear interpolation functions or conversely very fine

Overally, up to now the basic configuration and limiting characteristics for the forward and inverse problems are explained. We are, thus, ready to proceed to the description of forward

Electromagnetic problems are obviously governed by Maxwell equations in their full vector differential form, which in general include a temporal variation [10]. Even well-known, their solute, ion for complicated extremely inhomogeneous structures like the Human body constitutes a formidable task, especially when the conductivity anisotropy of skeletal and myocardial muscle is to be considered. Since, the inverse problem asks for multiple forward solutions (at least one for each illumination-projection angle) at each of its iterations, then a computationally efficient technique is inevitable. For this purpose a series of simplifications

The temporal dependence could be employed in imaging modalities, since radar type techniques could be employed, but these are efficient for electrically large bodies (dimensions of multiple wavelengths) located at distances of at least a few wavelengths away from the illuminating source. Besides that pulsed-waveform excitation could be exploited for imaging in general, where the received pulse delay *(Δτ)* in respect with the transmitted one provides the useful phase difference *(Δφ)*. In turn based on *Δφ* the ratio of imaginary to real part of the voltage or electric field is readily calculated, while the in principle required phase sensitive sensors are replaced by "scalar" measuring only their amplitude. As will be explained next, this information is directly proportional to *(εr/σ)* or inversely proportional to the dielectric loss tangent *[tanδ=σ/(ωε0εr)]* of the media profile to be reconstructed. This imaging approach could be classified as a "multistatic phase" radar and it is very challenging as well as a promising modality, especially regarding the involved relative simplification of the sensing-measuring instrumentation. However, it has only received a very limited attention and to the Authors knowledge, there was only one Russian group working on that, [11]. Conversely, most of the research groups activated in the field in MHz and microwave imaging (including our own) assume time harmonic (sinusoidal) excitation due to its simplicity in the forward problem solution, thus the following analysis will be restricted to this case. Note that this excitation eliminates the temporal dependence but also avoids (dues not account for) the materials frequency dispersion *σ(ω), εr(ω)* which introduce more difficulties. In contrary this dispersion can be exploited as an additional degree of freedom enabling "dynamic imaging" or calibration purposes and even a spectral imaging by performing measurements and reconstructing the media profile at multiple

mesh of linearly interpolated elements (λ/10 up to λ/32).

problem governing equations and their numerical solutions.

are necessary in order to reduce the involved complexity.

frequencies (σ(ωi), εr(ωi)) selected at appropriate steps (Δωi).

**3. Forward problem**

**3.1 Excitation rationale**

distribution should be discretized into N locally homogeneous elements comprising the so called piece-wise constant distribution (σi, εri, i=1-N). This is actually formed by the "reconstruction-mesh", which is in general a 3-D one (Fig.3a), but it could also restricted to a 2-D assembly of bars with in general arbitrary cross-section (Fig.3b with brick bars).

Fig. 3. Examples of reconstruction coarse meshes and forward fine meshes as: reconstruction (a) 3-D and (b) 2-D cases, forward: (c) 3-D and (d) 2-D cases corresponding to the structures (a) and (b) respectively.

The reasoning described above yields usually a coarse reconstruction mesh, which is always inappropriate for the forward problem solution. This one purpose is to ensure calculations providing an accuracy of the same order as the available measurements. In the low frequency EIT case it was proved by Barber and Seagar [9] that an accuracy of the order of 0.1% is necessary for a reliable reconstruction. On the other end, in the microwave regime it is widely understood that an acceptable simulation asks for a discretization with elements size smaller than one tenth of the wavelength (λ/10) for Moment Method (MoM) or Finite Element techniques (FEM), while the finest mesh of λ/32 is asked by Finite Difference Time Domain (FDTD) to ensure its stability. In any case trying to fulfill these requirements a very fine mesh is necessary for the forward problem solution, which also yields piece-wise constant (σ, εr) distributions. However, the inverse problem will be iteratively linearized first around the 0 0 (,) *r* initial guess and subsequently around the most recently updated (,) *k k r* distribution at its k-th iteration. Multiple forward problem solutions are then required on this (,) *k k r* distribution, hence the "forward" and "reconstruction" meshes should be compatible. A very straight forward approach is to subdivide each of the reconstruction mesh into smaller forward mesh elements, as shown in Fig.3(c)-(d). Even though the forward mesh elements comprising each reconstruction cell have the same constitutive (σi, εri) parameters, this subdivision in necessary in order to fullfill the voltage/field interpolation requirements. Explicitly, in most cases linear interpolation functions are considered within the established numerical techniques (MoM, FEM, FDTD). In turn the field variation especially around fine structures and around conductors are very high, thus asking either for highly non-linear interpolation functions or conversely very fine mesh of linearly interpolated elements (λ/10 up to λ/32).

Overally, up to now the basic configuration and limiting characteristics for the forward and inverse problems are explained. We are, thus, ready to proceed to the description of forward problem governing equations and their numerical solutions.

## **3. Forward problem**

154 Medical Imaging

distribution should be discretized into N locally homogeneous elements comprising the so called piece-wise constant distribution (σi, εri, i=1-N). This is actually formed by the "reconstruction-mesh", which is in general a 3-D one (Fig.3a), but it could also restricted to a

(a) (b)

(c) (d) Fig. 3. Examples of reconstruction coarse meshes and forward fine meshes as: reconstruction (a) 3-D and (b) 2-D cases, forward: (c) 3-D and (d) 2-D cases corresponding to the structures

The reasoning described above yields usually a coarse reconstruction mesh, which is always inappropriate for the forward problem solution. This one purpose is to ensure calculations providing an accuracy of the same order as the available measurements. In the low frequency EIT case it was proved by Barber and Seagar [9] that an accuracy of the order of 0.1% is necessary for a reliable reconstruction. On the other end, in the microwave regime it is widely understood that an acceptable simulation asks for a discretization with elements size smaller than one tenth of the wavelength (λ/10) for Moment Method (MoM) or Finite Element techniques (FEM), while the finest mesh of λ/32 is asked by Finite Difference Time Domain (FDTD) to ensure its stability. In any case trying to fulfill these requirements a very fine mesh is necessary for the forward problem solution, which also yields piece-wise constant (σ, εr) distributions. However, the inverse problem will be iteratively linearized first

initial guess and subsequently around the most recently updated

distribution, hence the "forward" and "reconstruction" meshes

distribution at its k-th iteration. Multiple forward problem solutions are then

should be compatible. A very straight forward approach is to subdivide each of the reconstruction mesh into smaller forward mesh elements, as shown in Fig.3(c)-(d). Even

(a) and (b) respectively.

around the 0 0 (,) *r*

required on this (,) *k k*

 *r*

(,) *k k r*

2-D assembly of bars with in general arbitrary cross-section (Fig.3b with brick bars).

Electromagnetic problems are obviously governed by Maxwell equations in their full vector differential form, which in general include a temporal variation [10]. Even well-known, their solute, ion for complicated extremely inhomogeneous structures like the Human body constitutes a formidable task, especially when the conductivity anisotropy of skeletal and myocardial muscle is to be considered. Since, the inverse problem asks for multiple forward solutions (at least one for each illumination-projection angle) at each of its iterations, then a computationally efficient technique is inevitable. For this purpose a series of simplifications are necessary in order to reduce the involved complexity.

## **3.1 Excitation rationale**

The temporal dependence could be employed in imaging modalities, since radar type techniques could be employed, but these are efficient for electrically large bodies (dimensions of multiple wavelengths) located at distances of at least a few wavelengths away from the illuminating source. Besides that pulsed-waveform excitation could be exploited for imaging in general, where the received pulse delay *(Δτ)* in respect with the transmitted one provides the useful phase difference *(Δφ)*. In turn based on *Δφ* the ratio of imaginary to real part of the voltage or electric field is readily calculated, while the in principle required phase sensitive sensors are replaced by "scalar" measuring only their amplitude. As will be explained next, this information is directly proportional to *(εr/σ)* or inversely proportional to the dielectric loss tangent *[tanδ=σ/(ωε0εr)]* of the media profile to be reconstructed. This imaging approach could be classified as a "multistatic phase" radar and it is very challenging as well as a promising modality, especially regarding the involved relative simplification of the sensing-measuring instrumentation. However, it has only received a very limited attention and to the Authors knowledge, there was only one Russian group working on that, [11]. Conversely, most of the research groups activated in the field in MHz and microwave imaging (including our own) assume time harmonic (sinusoidal) excitation due to its simplicity in the forward problem solution, thus the following analysis will be restricted to this case. Note that this excitation eliminates the temporal dependence but also avoids (dues not account for) the materials frequency dispersion *σ(ω), εr(ω)* which introduce more difficulties. In contrary this dispersion can be exploited as an additional degree of freedom enabling "dynamic imaging" or calibration purposes and even a spectral imaging by performing measurements and reconstructing the media profile at multiple frequencies (σ(ωi), εr(ωi)) selected at appropriate steps (Δωi).

High to Microwave Frequencies Imaging Techniques 157

Note that the term related to the actual permittivity *(ε)* is known as displacement

sources has actually no place herein since their spectral content is restricted to less than about 10KHz. Thus the impressed current in equation (4) is solely due to the source injecting current or illuminating the body. For safety reasons the injected current density should be less than 10mA/cm2 according to recommendation included in [12]. Besides these, recall the possible anisotropy issue which only concerns the conductivity of skeletal muscular tissues and the myocardium. Both of them exhibit a fibrous structure, where each fiber has a high conductivity interior surrounded by a thin low conductivity membrane. It is this structure that causes the conductivity anisotropy where parallel (along) the fibers it is σ//=7.3mS/cm while transverse to them it is only σ=0.49mS/cm. The question is how to include this anisotropy into the computer model. Mathematically the anisotropy is accounted for by a

be oriented parallel to the muscle fibers. However, these fibers are curved and even twisted,

computational complexity of full tensor parameters is tοo high to be considered in most

Returning back to Maxwell equations (2), it is well understood that the divergence equations (2c) and (2d) can be derived from the curl equations (2a) and (2b) by using continuity equation (3). Thus one may easily conclude that it is only necessary to solve the implicit pair of curl equations (2a), (2b) in conjunction with the boundary conditions for the tangential to the various media-tissue interfaces electric and magnetic field components. Recall that the latter are indeed obtained from the corresponding integral form laws resulting from the same curl equations, [10]. However, one should keep in mind that this conclusion is exactly valid only when an exact analytical solution is sought. In contrary when numerical solutions are adopted along with the internal structure discretization accompanied by approximate interpolation functions within each piece-wise homogeneous element, then the aforementioned "exact" conditions can be violated, resulting to field distributions not obeying the divergence conditions. Even though this could be kept in mind as a possible source of inaccuracies, usually their effects can be negligible. Thus there is indeed a very powerful method to discretized and solve the interdependent pair of vector curl differential equations (2a) and (2b), such as the Finite Difference Frequency Domain (FDFD) method established in our previous work, Lavranos [13] in two-dimensional curvilinear coordinates and currently been extended to 3-D structures as in Lavdas [14]. This approach can easily account for material anisotropy as well, within the usual FD limitations regarding the interfaces spatial resolution. Even though this is a promising approach we have not yet employed it within inverse problem solutions, but instead the Finite Element Method (FEM) is mostly employed for this purpose. Since FEM is formulated in an integral form, it is more convenient to decouple the two curl equations in order to formulate two separate vector wave equations for the electric and the magnetic field

 or *r*

2

r 0

0 r 0

(6a)

1 x -k j *imp xE j E J*

cases. Thus, herein the inherent muscles anisotropy will not be considered.

*E* . A usual question regarding the body inherent bioelectric

, but this could be only convenient if one of the coordinate axis could

. It is, indeed, very easy to

parameters, but the mathematical and

current *d r* 0 *J Et j*

second rank tensor

model a diagonal tensor of either

respectively as, e.g. Jin [15] or Volakis [16](p.5) :

hence they could be only represented by a full 3×3 tensorial

#### **3.2 Time harmonic fields & currents**

Maxwell equations for the time harmonic excitation of the form *j t e* are significantly simplified since the temporal variation is substituted as / *t j* , while the actual field *tr tr* ,, , and source quantities *t r*, are replaced by the corresponding complex phasors *Er Hr J r* , , as:

$$\vec{\mathcal{E}}\left(t,\vec{r}\right) = \text{Re}\left(\vec{E}\left(\vec{r}\right)e^{j\phi t}\right) \tag{1a}$$

$$\vec{\mathcal{H}}(t,\vec{r}) = \text{Re}\left(\vec{H}(\vec{r})e^{j\alpha t}\right) \tag{1b}$$

$$\vec{\mathcal{J}}\left(t,\vec{r}\right) = \text{Re}\left(\vec{J}\left(\vec{r}\right)e^{j\alpha t}\right) \tag{1c}$$

In turn Maxwell equations for time harmonic fields read:

$$
\vec{\nabla} \mathbf{x} \vec{\mathbf{E}} = \text{-j}\rho \mu\_0 \mu\_\mathbf{r} \vec{\mathbf{H}} \tag{2a}
$$

$$
\vec{\nabla} \mathbf{x} \vec{\mathbf{H}} = \mathbf{j} \rho \mathbf{z}\_0 \mathbf{z}\_r \vec{\mathbf{E}} + \vec{\mathbf{J}} \tag{2b}
$$

$$
\vec{\nabla} \cdot \vec{\mathbf{D}} = \rho \Rightarrow \vec{\nabla} \cdot \boldsymbol{\varepsilon} \vec{\mathbf{E}} = \rho \tag{2c}
$$

$$
\vec{\nabla} \cdot \vec{B} = 0 \implies \vec{\nabla} \cdot \mu \vec{H} = 0 \Big|\_{\mu=\text{scalar}} \to \vec{\nabla} \cdot \vec{H} = 0 \tag{2d}
$$

Where the source quantities current ( )*J* and charge *(ρ)* densities should obey the continuity equation:

$$
\vec{\nabla} \cdot \vec{\J} = -\hat{\mathbf{c}} \rho \;/\,\hat{\mathbf{c}}
\\
\mathbf{t} = -j\alpha\rho\tag{3}
$$

Additionally, for the conducting media assumed herein there is a conduction current ( ) *c J* flowing at any arbitrary point (movement of charges due to coulomb forces and the presence of electric field) as *c J E* . Thus, the current density *J* appearing in (2b) is comprised of the conduction current *c J* and the current impressed by the source as *imp J* . Substituting in (2a) we may define an equivalent complex permittivity *εc* or an equivalent complex conductivity *σc* as:

$$\vec{\nabla} \mathbf{x} \vec{\mathbf{H}} = \mathbf{j} \rho \boldsymbol{\varepsilon}\_0 \boldsymbol{\varepsilon}\_\mathbf{r} \vec{E} + \sigma \vec{E} + \vec{J}\_{imp} = \mathbf{j} \rho \boldsymbol{\varepsilon}\_0 \boldsymbol{\varepsilon}\_\mathbf{r} \vec{E} + \vec{J}\_{imp} = \sigma\_\mathbf{c} \vec{E} + \vec{J}\_{imp} \tag{4}$$

where

$$
\varepsilon\_{\rm c} = \varepsilon\_0 \varepsilon\_{\rm rc} = \varepsilon\_0 \varepsilon\_{\rm r} (1 - j \tan \delta) = \varepsilon\_0 \varepsilon\_{\rm r} \left( 1 - j \frac{\sigma}{\alpha \varepsilon\_0 \varepsilon\_{\rm r}} \right) \tag{5a}
$$

$$
\sigma\_c = \sigma\_{eff} = \sigma + j\alpha \varepsilon\_0 \varepsilon\_\mathbf{r} \tag{5b}
$$

*tr tr* ,, , and source quantities *t r*, are replaced by the corresponding complex

*j t tr Er e* , Re

*j t tr H r e* , Re

*j t tr J r e* , Re

0 r xE -j H

0 r xH j E J

 D E

*BH H*

> *j*

Additionally, for the conducting media assumed herein there is a conduction current ( ) *c J*

flowing at any arbitrary point (movement of charges due to coulomb forces and the

Substituting in (2a) we may define an equivalent complex permittivity *εc* or an equivalent

*c eff* 0 r

0 r 0 rc c xH j j

(1 tan ) 1 *j j*

c 0 rc 0 r 0 r

>

>

00 0 *scalar*

(2a)

(2b)

.

appearing in (2b) is

(5a)

(1a)

(1b)

(1c)

(2c)

and charge *(ρ)* densities should obey the continuity

(3)

and the current impressed by the source as *imp J*

0 r

*j* (5b)

(2d)

. Thus, the current density *J*

 *E EJ EJ EJ imp imp imp* (4)

are significantly

, while the actual field

Maxwell equations for the time harmonic excitation of the form *j t e*

simplified since the temporal variation is substituted as / *t j*

In turn Maxwell equations for time harmonic fields read:

**3.2 Time harmonic fields & currents**

Where the source quantities current ( )*J*

presence of electric field) as *c J E*

complex conductivity *σc* as:

comprised of the conduction current *c J*

/ *J t*

equation:

where

phasors *Er Hr J r* , , as:

Note that the term related to the actual permittivity *(ε)* is known as displacement current *d r* 0 *J Et j E* . A usual question regarding the body inherent bioelectric sources has actually no place herein since their spectral content is restricted to less than about 10KHz. Thus the impressed current in equation (4) is solely due to the source injecting current or illuminating the body. For safety reasons the injected current density should be less than 10mA/cm2 according to recommendation included in [12]. Besides these, recall the possible anisotropy issue which only concerns the conductivity of skeletal muscular tissues and the myocardium. Both of them exhibit a fibrous structure, where each fiber has a high conductivity interior surrounded by a thin low conductivity membrane. It is this structure that causes the conductivity anisotropy where parallel (along) the fibers it is σ//=7.3mS/cm while transverse to them it is only σ=0.49mS/cm. The question is how to include this anisotropy into the computer model. Mathematically the anisotropy is accounted for by a second rank tensor , but this could be only convenient if one of the coordinate axis could be oriented parallel to the muscle fibers. However, these fibers are curved and even twisted, hence they could be only represented by a full 3×3 tensorial . It is, indeed, very easy to model a diagonal tensor of either or *r* parameters, but the mathematical and computational complexity of full tensor parameters is tοo high to be considered in most cases. Thus, herein the inherent muscles anisotropy will not be considered.

Returning back to Maxwell equations (2), it is well understood that the divergence equations (2c) and (2d) can be derived from the curl equations (2a) and (2b) by using continuity equation (3). Thus one may easily conclude that it is only necessary to solve the implicit pair of curl equations (2a), (2b) in conjunction with the boundary conditions for the tangential to the various media-tissue interfaces electric and magnetic field components. Recall that the latter are indeed obtained from the corresponding integral form laws resulting from the same curl equations, [10]. However, one should keep in mind that this conclusion is exactly valid only when an exact analytical solution is sought. In contrary when numerical solutions are adopted along with the internal structure discretization accompanied by approximate interpolation functions within each piece-wise homogeneous element, then the aforementioned "exact" conditions can be violated, resulting to field distributions not obeying the divergence conditions. Even though this could be kept in mind as a possible source of inaccuracies, usually their effects can be negligible. Thus there is indeed a very powerful method to discretized and solve the interdependent pair of vector curl differential equations (2a) and (2b), such as the Finite Difference Frequency Domain (FDFD) method established in our previous work, Lavranos [13] in two-dimensional curvilinear coordinates and currently been extended to 3-D structures as in Lavdas [14]. This approach can easily account for material anisotropy as well, within the usual FD limitations regarding the interfaces spatial resolution. Even though this is a promising approach we have not yet employed it within inverse problem solutions, but instead the Finite Element Method (FEM) is mostly employed for this purpose. Since FEM is formulated in an integral form, it is more convenient to decouple the two curl equations in order to formulate two separate vector wave equations for the electric and the magnetic field respectively as, e.g. Jin [15] or Volakis [16](p.5) :

$$\vec{\nabla} \mathbf{x} \left(\frac{1}{\mu\_{\rm r}} \vec{\nabla} \mathbf{x} \vec{E}\right) \cdot \mathbf{k}\_0^2 \left(\boldsymbol{\varepsilon}\_{\rm r} - j \frac{\boldsymbol{\sigma}}{\alpha \varepsilon\_0}\right) \vec{E} = -j \alpha \mu\_0 \vec{J}\_{imp} \tag{6a}$$

High to Microwave Frequencies Imaging Techniques 159

*J j*

*E V*

Equations (7b) and (7c) seems to be in contradiction since, they call for different voltage definition or making V to voltage uniqueness. Here comes the original Maxwell observation, e.g. see Jackson [19]( p.238), that one could differentiate (2c) and substitute in charge

0 (2 )

*D D J c t t t*

In turn utilizing (7a) the so called **generalized complex Laplace equation** is obtained:

 *c*

*J V V j jj V n n j*

  (3) 0 0

 ( )0

 0 *c*

However, besides this classical approach, (8c) can be obtained by dividing (7c) with *jω* and adding the result to (7b). Charge density at the operating frequency indeed exist over the metallic electrodes obeying (7c), but these can be taken into account through the **boundary conditions** to be enforced on the voltage distribution during the FEM formulation of the

> on the body air interface *c o*

where V/n denotes the normal derivative at the interface. At this point is important to note that (10) applies to both active-driven as well as passive (sensing) metallic electrodes, since current and charge density is induced on them. Most important is that the field and current singularity at metallic edges is a local phenomenon and thus it occurs at any frequency from DC to microwaves. Referring to a classical electromagnetics text book, e.g. Collin [10](p.25), field components normal to the edge as well as current density components parallel to the edge exhibit a singular behavior tending to infinity as 1 /

where ρ the distance from the edge. Conversely, field components parallel to the edge and

*n*

*D J or J t t*

and *D E*

conservation equation as:

to time harmonic fields we get:

where:

Newmann type:

Substituting the constitutive relations *J E*

*E EV* 0 (7a)

(7b)

(8a)

into (8a) and restricting ourselves

(10)

,

*j E* (8b)

*r V* (8c)

*r jr* (9)

on the electrode

*E j V j* (7c)

$$\bar{\nabla}\mathbf{x}\left(\frac{1}{\varepsilon\_{\rm rc}}\vec{\nabla}\mathbf{x}\vec{H}\right)\cdot\mathbf{k}\_{0}^{2}\mu\_{\rm r}\vec{H}=\nabla\mathbf{x}\left(\frac{\vec{J}}{\varepsilon\_{\rm rc}}\right)\tag{6b}$$

Either wave equation could be solved employing FEM after constructing a functional (or weak form) to be minimized. This task is achieved utilizing either a variational approach or the Galerkin approach, e.g. Cangellaris [17]. The choice among the two wave equations is mainly undermined by the more convenient enforcement of boundary conditions. For the biomedical applications the media-tissues are non-magnetic (μ=μ0) and the inhomogeneity is approximated by step changes in conductivity and permittivity. Thus the related boundary conditions are directly enforced to the electric field as Dirichlet type. Conversely, if ones solves for the magnetic field the same boundary condition would be translated into a Newmann type and a differentiation after the solution (prone to numerical errors) would be performed to evaluate the desired electric field. Hence, it is preferable to choose a straightforward solution of the electric field wave equation (6a).

Observing equation (6a) and following its solution in the next section one would realize that this is anything else but a trivial task, involving huge computational resources to solve over a Human body or just a Human torso. For this reason two simplifications will be presented next, a quasi-static approach valid for the frequency range up to a few MHz and a second approach with unidirectional dipole sources (aligned along the z-axis) and even a 2-D scalar Helmholtz wave equations excited by "line sources" which significantly reduce the involved complexity.

#### **3.3 Quasi-static approach**

Observing equations (2a) and (2b) one may recall the basic understanding that the electromagnetic waves are created and propagated due to the temporal variation of impressed current source which is inherited to the generated field. Explicitly, a time depended impressed current generates a temporally changing magnetic field through (2b). In turn the varying magnetic field / *H t jH* produce an electric field with an identical time variation through (2a). This / *E t jE* regenerates a changing magnetic field through (2b) and this cycle is infinitely repeated producing a propagating wave. However, this "chain" becomes very weak or breaks when either jωμ=jωμ0μr or jωε=jωε0εr quantities become very small (where ε0=8.854x10-12F/m and μ0=4π10-7Η/m). Biological tissues are non-magnetic μr=1 while due to their high water content they present a very high dielectric constant which starts from about εr1000 (or more) at 100KHz range and reduced to about εr80 in the microwave (GHz) range. Based on this reasoning it was long ago realized, e.g. Price 1979 [18], that magnetic effects can be ignored for frequencies lower than about 10MHz. Namely, the electric filed generated by magnetic field temporal variations becomes negligible. From the wave propagation point of view the wavelength in a media εr100 at f=10MHz is 0 3 *g r m* (and *λ0=c/f*), which means that at a distance d=λg/4=75cm the field or voltages changes due to wave propagation from its maximum value to zero and vice-versa. Hence for a thorax with larger dimension of about 40cm the 10MHz frequency constitutes indeed an upper limit.

In view of the above, the ignorant question asks then "from where comes the varying electric field and the desired electric potential"? Now we have to reconsider the divergence equation (2c) and the charge conservation law (3) as:

$$
\vec{\nabla} \times \vec{E} \approx 0 \quad \Rightarrow \quad \vec{E} = -\vec{\nabla} V \tag{7a}
$$

$$
\vec{\nabla} \cdot \varepsilon \vec{E} = \rho \quad \Rightarrow \quad \vec{\nabla} \cdot \varepsilon \vec{\nabla} V = -\rho \tag{7b}
$$

$$
\vec{\nabla} \cdot \vec{J} = -j\alpha\rho\sigma \implies \vec{\nabla} \cdot \sigma \vec{E} = -j\alpha\rho\sigma \implies \vec{\nabla} \cdot \sigma \vec{\nabla}V = j\alpha\rho\sigma \tag{7c}
$$

Equations (7b) and (7c) seems to be in contradiction since, they call for different voltage definition or making V to voltage uniqueness. Here comes the original Maxwell observation, e.g. see Jackson [19]( p.238), that one could differentiate (2c) and substitute in charge conservation equation as:

$$\begin{aligned} \text{(2c)} \quad &\Rightarrow \quad \vec{\nabla} \cdot \frac{\partial \vec{D}}{\partial t} = \frac{\partial \rho}{\partial t} \Bigg|\_{} \vec{\nabla} \cdot \vec{J} + \vec{\nabla} \cdot \frac{\partial \vec{D}}{\partial t} = 0 \\ \text{(3)} \quad &\Rightarrow \quad \vec{\nabla} \cdot \vec{J} + \frac{\partial \rho}{\partial t} = 0 \Bigg|\_{} \text{or} \quad \vec{\nabla} \cdot \left( \vec{J} + \frac{\partial \vec{D}}{\partial t} \right) = 0 \end{aligned} \tag{8a}$$

Substituting the constitutive relations *J E* and *D E* into (8a) and restricting ourselves to time harmonic fields we get:

$$
\overline{\nabla} \cdot (\sigma + jo\nu\varepsilon) \overline{E} = 0 \tag{8b}
$$

In turn utilizing (7a) the so called **generalized complex Laplace equation** is obtained:

$$
\nabla \cdot \varepsilon\_c \left( \vec{r} \right) \vec{\nabla} V = 0 \tag{8c}
$$

where:

158 Medical Imaging

2 0 r rc rc 1 x -k *J xH H x*

Either wave equation could be solved employing FEM after constructing a functional (or weak form) to be minimized. This task is achieved utilizing either a variational approach or the Galerkin approach, e.g. Cangellaris [17]. The choice among the two wave equations is mainly undermined by the more convenient enforcement of boundary conditions. For the biomedical applications the media-tissues are non-magnetic (μ=μ0) and the inhomogeneity is approximated by step changes in conductivity and permittivity. Thus the related boundary conditions are directly enforced to the electric field as Dirichlet type. Conversely, if ones solves for the magnetic field the same boundary condition would be translated into a Newmann type and a differentiation after the solution (prone to numerical errors) would be performed to evaluate the desired electric field. Hence, it is preferable to choose a

Observing equation (6a) and following its solution in the next section one would realize that this is anything else but a trivial task, involving huge computational resources to solve over a Human body or just a Human torso. For this reason two simplifications will be presented next, a quasi-static approach valid for the frequency range up to a few MHz and a second approach with unidirectional dipole sources (aligned along the z-axis) and even a 2-D scalar Helmholtz wave equations excited by "line sources" which significantly reduce the involved complexity.

Observing equations (2a) and (2b) one may recall the basic understanding that the electromagnetic waves are created and propagated due to the temporal variation of impressed current source which is inherited to the generated field. Explicitly, a time depended impressed current generates a temporally changing magnetic field through (2b).

*H t j*

field through (2b) and this cycle is infinitely repeated producing a propagating wave. However, this "chain" becomes very weak or breaks when either jωμ=jωμ0μr or jωε=jωε0εr quantities become very small (where ε0=8.854x10-12F/m and μ0=4π10-7Η/m). Biological tissues are non-magnetic μr=1 while due to their high water content they present a very high dielectric constant which starts from about εr1000 (or more) at 100KHz range and reduced to about εr80 in the microwave (GHz) range. Based on this reasoning it was long ago realized, e.g. Price 1979 [18], that magnetic effects can be ignored for frequencies lower than about 10MHz. Namely, the electric filed generated by magnetic field temporal variations becomes negligible. From the wave propagation point of view the wavelength in a media

d=λg/4=75cm the field or voltages changes due to wave propagation from its maximum value to zero and vice-versa. Hence for a thorax with larger dimension of about 40cm the

In view of the above, the ignorant question asks then "from where comes the varying electric field and the desired electric potential"? Now we have to reconsider the divergence

*E t j*

*H*

*E*

*g r m* (and *λ0=c/f*), which means that at a distance

produce an electric field with an

regenerates a changing magnetic

(6b)

straightforward solution of the electric field wave equation (6a).

**3.3 Quasi-static approach**

In turn the varying magnetic field /

identical time variation through (2a). This /

εr100 at f=10MHz is 0 3

10MHz frequency constitutes indeed an upper limit.

equation (2c) and the charge conservation law (3) as:

$$
\varepsilon\_c \left( \rho \right) = \varepsilon \left( \vec{r} \right) - j \left. \sigma \left( \vec{r} \right) \right| \rho \tag{9}
$$

However, besides this classical approach, (8c) can be obtained by dividing (7c) with *jω* and adding the result to (7b). Charge density at the operating frequency indeed exist over the metallic electrodes obeying (7c), but these can be taken into account through the **boundary conditions** to be enforced on the voltage distribution during the FEM formulation of the Newmann type:

$$-j\alpha \varepsilon\_c \frac{\partial V}{\partial \mathfrak{n}} = -j\alpha \left(\varepsilon - j\frac{\sigma}{\alpha}\right) \frac{\partial V}{\partial \mathfrak{n}} = \begin{cases} \mathcal{J} & \text{on the electrode} \\ -j\alpha \varepsilon\_o \frac{\partial V}{\partial \mathfrak{n}} & \text{on the body air interface} \end{cases} \tag{10}$$

where V/n denotes the normal derivative at the interface. At this point is important to note that (10) applies to both active-driven as well as passive (sensing) metallic electrodes, since current and charge density is induced on them. Most important is that the field and current singularity at metallic edges is a local phenomenon and thus it occurs at any frequency from DC to microwaves. Referring to a classical electromagnetics text book, e.g. Collin [10](p.25), field components normal to the edge as well as current density components parallel to the edge exhibit a singular behavior tending to infinity as 1 / , where ρ the distance from the edge. Conversely, field components parallel to the edge and

High to Microwave Frequencies Imaging Techniques 161

Substituting equation (12b) into (11) and the resulting expression back into (6a) the desired

0 0

( ) () () *r r r Kr k r j*

() () () ()() ( ) ( ) *r imp Er K r Er K rEr j Jr K r*

The convenience offered by a two dimensional imaging are already explained in section II.2 along with the necessity of solving a three-dimensional forward problem. The geometry considered in this case is a cylindrical structure with an arbitrary shaped cross-section but uniform along its axial *z*ˆ -direction. Likewise the material complex permittivity to be sought is inhomogeneous in the (xy) cross section but uniform along the z-direction, like an

>

direction the structure can be considered as an inhomogeneously loaded lossy dielectric open waveguide of arbitrary cross-section which is excited-illuminated by different type of sources. In view of this uniformity the electromagnetic field within and outside this open waveguide could be expressed as a superposition (expansion) of an in general infinite number of modes either propagating or evanescent just as classically done in the mode-Matching technique. Even though this is a very promising approach to our knowledge it has not yet being exploited for imaging purposes.  For each one of the possibly excited modes the field dependence in the *z*ˆ -

(,,) (,) *i j z E xyz e xye i i*

(,,) (,) *i j z H xyz h xye i i*

where *βi* is the ith mode complex propagation constant. When β is real or exhibits a small imaginary part due to losses, then it represents a propagating wave. Conversely, when β is purely imaginary then it represents an evanescent or non-propagating wave which can be excited by the source but is exponentially attenuated away from that. Namely, it exists only

>

*E j Ek j E j J j*

*r r r imp*

( )

0

(15)

2

2 0

(13)

(14)

*r xy* . Due to the uniformity in the axial *z*ˆ -

(16a)

(16b)

in the usual

0

2 2 0

It is now convenient to define a complex inhomogeneous wavenumber *K r*( )

electric field wave equation for biological media is obtained:

2 2

form as:

where ( ) *r*

With the aid of (14) the electric wave equation reads

2 2

denotes the source spatial vector.

assembly of different dielectric bars, as () (, ) *c c*

direction can be denoted as:

**3.5 Two dimensional structure-inhomogeneous cross-section**

current density normal to that vanish as ρ0. This is an important behavior and it should be taken into account when an accurate modeling is sought. To clarify this behavior the current density on a patch electrode driven by a current 1 is shown in Fig.4a and the induced current density on a passive electrode is sketched in Fig.4b. Note that the corresponding integrals denote the total current I on the driven and zero on the passive electrode. Obeying equation (3) the charge density on the active or passive electrodes exhibits identical singular behavior.

Fig. 4. Current density on metallic strip electrodes exhibiting the inherent singularity as a) emanating from an active electrode and b) induced on a passive electrode, e.g. [19](p.165)

#### **3.4 Wave equation for microwave sources in biological media**

The forward problem for the microwave imaging is governed by the wave equations (6a) or (6b). It was explained before that the electric field wave equation (6a) offers a straightforward enforcement of boundary condition at the different lossy dielectric interfaces as a rigid-Dirichlet type. Besides that, observing the first term of equation (6a) this could be significantly simplified for the non-magnetic biological media as *μ=μ0=const* or *μr=1*. In contrary the corresponding first term of (6b) is complicated as it contains the spatially dependent complex ( ) *rc rc r* . In view of this clarification the electric field wave equation (6a) can be simplified by utilizing the identity:

$$
\nabla \times \overline{\nabla} \times \overline{E} = \vec{\nabla} (\vec{\nabla} \cdot \vec{E}) - \nabla^2 \vec{E} \tag{11}
$$

For the divergence of *E* we may again exploit the original Maxwell observation [19](p.238) as given in equation (8c) which is generally valid and using the symbolism (5a) *εc=σ+jωε*, then (8c) reads:

$$\vec{\nabla} \cdot \left( \varepsilon\_c(\vec{r}) \vec{E}(\vec{r}) \right) = 0 \quad \rightarrow \quad \varepsilon\_c(\vec{r}) \vec{\nabla} \cdot \vec{E}(\vec{r}) + \vec{E}(\vec{r}) \cdot \vec{\nabla} \varepsilon\_c(\vec{r}) = 0 \tag{12a}$$

or

$$
\vec{\nabla} \cdot \vec{E}(\vec{r}) = -\frac{\vec{E}(\vec{r}) \cdot \vec{\nabla} \varepsilon\_c(\vec{r})}{\varepsilon\_c(\vec{r})} \tag{12b}
$$

current density normal to that vanish as ρ0. This is an important behavior and it should be taken into account when an accurate modeling is sought. To clarify this behavior the current density on a patch electrode driven by a current 1 is shown in Fig.4a and the induced current density on a passive electrode is sketched in Fig.4b. Note that the corresponding integrals denote the total current I on the driven and zero on the passive electrode. Obeying equation (3) the charge density on the active or passive electrodes exhibits identical singular behavior.

(a) (b) Fig. 4. Current density on metallic strip electrodes exhibiting the inherent singularity as a) emanating from an active electrode and b) induced on a passive electrode, e.g. [19](p.165)

The forward problem for the microwave imaging is governed by the wave equations (6a) or (6b). It was explained before that the electric field wave equation (6a) offers a straightforward enforcement of boundary condition at the different lossy dielectric interfaces as a rigid-Dirichlet type. Besides that, observing the first term of equation (6a) this could be significantly simplified for the non-magnetic biological media as *μ=μ0=const* or *μr=1*. In contrary the corresponding first term of (6b) is complicated as it contains the

2 *E EE* ( )

as given in equation (8c) which is generally valid and using the symbolism (5a) *εc=σ+jωε*,

() () ( ) ( )

*Er r E r*

*c*

*r*

we may again exploit the original Maxwell observation [19](p.238)

*cc c* ()() 0 () () () () 0 *rEr r Er Er r* (12a)

*c*

(12b)

. In view of this clarification the electric field wave

(11)

**3.4 Wave equation for microwave sources in biological media**

 *r*

equation (6a) can be simplified by utilizing the identity:

spatially dependent complex ( ) *rc rc*

For the divergence of *E*

then (8c) reads:

or

Substituting equation (12b) into (11) and the resulting expression back into (6a) the desired electric field wave equation for biological media is obtained:

$$\left(\vec{\nabla}^2 \vec{E} + k\_0^2 \mu\_r \left(\varepsilon\_r - j\frac{\sigma}{\alpha \varepsilon\_0}\right) \vec{E} + \vec{\nabla} \left(\frac{\vec{E} \cdot \vec{\nabla} (\sigma + j\alpha \varepsilon)}{\sigma + j\alpha \varepsilon}\right)\right) = j\alpha \mu\_0 \mu\_r \vec{J}\_{imp} \tag{13}$$

It is now convenient to define a complex inhomogeneous wavenumber *K r*( ) in the usual form as:

$$K^2(\vec{r}) = k\_0^2 \mu\_r \left( \varepsilon\_r(\vec{r}) - j \frac{\sigma(\vec{r})}{\alpha \varepsilon\_0} \right) \tag{14}$$

With the aid of (14) the electric wave equation reads

$$\nabla^2 \vec{E}(\vec{r}) + K^2(\vec{r}) \vec{E}(\vec{r}) + \vec{\nabla} \left( \frac{\vec{E}(\vec{r}) \cdot \vec{\nabla} K^2(\vec{r})}{K^2(\vec{r})} \right) = j\rho\mu\_0\mu\_r \vec{J}\_{imp}(\vec{r}') \tag{15}$$

where ( ) *r* denotes the source spatial vector.

#### **3.5 Two dimensional structure-inhomogeneous cross-section**

The convenience offered by a two dimensional imaging are already explained in section II.2 along with the necessity of solving a three-dimensional forward problem. The geometry considered in this case is a cylindrical structure with an arbitrary shaped cross-section but uniform along its axial *z*ˆ -direction. Likewise the material complex permittivity to be sought is inhomogeneous in the (xy) cross section but uniform along the z-direction, like an assembly of different dielectric bars, as () (, ) *c c r xy* . Due to the uniformity in the axial *z*ˆ direction the structure can be considered as an inhomogeneously loaded lossy dielectric open waveguide of arbitrary cross-section which is excited-illuminated by different type of sources. In view of this uniformity the electromagnetic field within and outside this open waveguide could be expressed as a superposition (expansion) of an in general infinite number of modes either propagating or evanescent just as classically done in the mode-Matching technique. Even though this is a very promising approach to our knowledge it has not yet being exploited for imaging purposes.  For each one of the possibly excited modes the field dependence in the *z*ˆ direction can be denoted as:

$$
\vec{E}\_i(\mathbf{x}, y, z) = \vec{e}\_i(\mathbf{x}, y)e^{-j\beta\_i z} \tag{16a}
$$

$$
\vec{H}\_i(\mathbf{x}, y, z) = \vec{h}\_i(\mathbf{x}, y)e^{-j\beta\_i z} \tag{16b}
$$

where *βi* is the ith mode complex propagation constant. When β is real or exhibits a small imaginary part due to losses, then it represents a propagating wave. Conversely, when β is purely imaginary then it represents an evanescent or non-propagating wave which can be excited by the source but is exponentially attenuated away from that. Namely, it exists only

High to Microwave Frequencies Imaging Techniques 163

1 () 0 ˆ ˆ *t tz t rc z*

The pair of equations (20) constitutes the eigenproblem to be solved numerically. Observing

mainly occurs a product jβ *te*

*t tt t z t rc t*

1 ( ) ˆ ˆ *t tz t rc z*

*e e z k ez*

Equations (21) constitute the desired eigenproblem for the *β* eigenvalue. The interested reader may contact our previous work [22] on open waveguides which was carried out

Even though the above 2-D analysis offered some simplification the forward problem still remains complicated and computationally demanding. A truly two-dimensional approach presumes the sources to be two dimensional, i.e. tending to infinity along the z-axis, just as

The field of an infinite electric line source embedded in a homogeneous media is well studied and can be found in any classical electromagnetics textbook, e.g. Balanis [24]( p.571). Since the line source is infinite its current density must be constant and if it is considered

> 2 () ( ) ˆ ˆ *L I Jr J z z A*

where I is a sinusoidal current (I=I0ejωt) and ρ the polar radial coordinate transverse to *z*ˆ , as

For a homogeneous medium the field generated by an infinite line source is a Transverse

assuming that the line source is at ρ=0, for details in [24]( p.571). For the case of the inhomogeneous cross-section the TM nature of the wave is exactly preserved, thus ˆ *E Ez z*

and Ex=Ey=0. Besides that the two-dimensional complex dielectric profile is uniform along the

above simplifications the third term of equation (15) vanishes as also explained by Fang [25] :

() () 0 ( )/ 0 *E E x y Er K r*

0

and consequently from equation (14) it is *Kr z* ( )/ 0 . In view of the

2 2

(23)

eigenfunction is identical to a cylindrical Hankel function of the second type (2)

and multiplying (20a) by jβ yields:

2 2 2

1 1

*r r*

*r*

within Dr. Allilomes doctoral thesis [23] also available online (in Greek).

*e j e z k ez*

*r*

**3.6 Two dimensional object excited by infinite line sources**

A let the line source cross section assumed very small.

Magnetic (TMz), Hz=0 with only and axial electric field ˆ *E Ez z*

2 2

*Kr z*

that the transverse component *te*

in the case of an "infinite line source".

aligned along the z-axis, then:

z-axis or ( )/ 0 *c d r dz*

Line Source:

simplified by letting *t t je e*

2 0

2 2

*e ee k e*

(20b)

0

(21b)

(22)

and its transverse

0 *H k*( )

(21a)

0

, the eigenproblem in further

around the source, it does not transfer energy but it only stores electric and magnetic energy in its neighborhood, thus contributing only to the reactive part of the source (antenna) input impedance. In order to utilize the eigenmode expansion approach the characteristics of every possible mode, its eigenfunctions ( , ), ( , ) *i i e xy h xy* and the complex propagation (βi) constant should be estimated first. The corresponding expansion formally reads:

$$\vec{E} = \sum\_{i=1}^{\infty} \boldsymbol{w}\_{ei} \vec{E}\_i = \sum\_{i=1}^{\infty} \boldsymbol{w}\_{ei} \vec{e}\_i(\mathbf{x}, \boldsymbol{y}) e^{-j\beta\_i z} \tag{17a}$$

$$
\vec{H} = \sum\_{i=1}^{\infty} w\_{mi} \vec{H}\_i = \sum\_{i=1}^{\infty} w\_{ei} \vec{h}\_i(\mathbf{x}, \mathbf{y}) e^{-j\beta\_i z} \tag{17b}
$$

Where wei and wmi represent the modal amplitudes or weighting factors which depend on the source type and its location, see for example our work [20]. These can be estimated through a power conservation low which is the basis of a mode matching technique. Explicitly, each specific source can be enclosed inside a fictitious box inside which the field is described with the aid of a numerical technique on a fine mesh exploiting the knowledge of source modeling, e.g. Volakis [16](p.238). Outside the fictitious surface the eigenfunctions (17) is considered and the weighting factors wei and wmi are evaluated through field continuity conditions across this psudo surface. It is important to keep in mind that the power conservation low should be enforced, which is preferably formulated in conjunction with the modes orthogonality properties, e.g. [21].

The above mode matching approach pressumes the knowledge of the eigenvalues *(βi)* and the eigenfunctions i *ei* , h which for the considered inhomogeneous cross-section can only be acquired numerically. For this purpose the wave equation (6a) shall be formulated as an eigenvalue problem and the usual approach is based on the separation of both the electric field *e* and the nabla differential operator into axial ( *z*ˆ ) and transverse *(t)* components as:

$$
\vec{e} = \vec{e}\_t + \hat{z}e\_z \quad \text{and} \quad \vec{e}\_t = e\_x\hat{\mathbf{x}} + e\_y\hat{\mathbf{y}}, \qquad \vec{h} = \vec{h}\_t + \hat{z}h\_z \tag{18a}
$$

$$
\vec{\nabla} = \vec{\nabla}\_t + \hat{z}\frac{\vec{\mathcal{O}}}{\partial \mathbf{z}} = \vec{\nabla}\_t - j\beta\hat{z}, \qquad \vec{\nabla}\_t = \hat{x}\frac{\vec{\mathcal{O}}}{\partial \mathbf{x}} + \hat{y}\frac{\vec{\mathcal{O}}}{\partial y} \tag{18b}
$$

The detailed procedure given in a lot of textbooks, e.g. Volakis et al [16](p.98) is summarized as follows:

$$
\vec{\nabla} \times \vec{\varepsilon} = \vec{\nabla}\_t \times \vec{e}\_t + (\vec{\nabla}\_t e\_z + j\mathcal{J}\vec{e}\_t) \times \hat{z} \tag{19}
$$

Which can be substitute into (6a) to yield a pair of differential equations for the axial (ez) and transverse ( *te* ) field components as:

$$\vec{\nabla}\_t \times \left(\frac{1}{\mu\_r} \vec{\nabla}\_t \times \vec{e}\_t\right) - j \frac{\mathcal{J}}{\mu\_r} \left(\vec{\nabla}\_t e\_z + j\beta \vec{e}\_t\right) - k\_0^2 \varepsilon\_{re} \vec{e}\_t = 0 \tag{20a}$$

around the source, it does not transfer energy but it only stores electric and magnetic energy in its neighborhood, thus contributing only to the reactive part of the source (antenna) input impedance. In order to utilize the eigenmode expansion approach the characteristics of every possible mode, its eigenfunctions ( , ), ( , ) *i i e xy h xy* and the complex propagation (βi)

(.) *i j z*

(.) *i j z*

(17b)

(17a)

constant should be estimated first. The corresponding expansion formally reads:

1 1

1 1

*i i*

with the modes orthogonality properties, e.g. [21].

) field components as:

eigenfunctions i *ei* , h

*e*

as follows:

transverse ( *te*

*i i*

*ei i ei i*

*mi i ei i*

Where wei and wmi represent the modal amplitudes or weighting factors which depend on the source type and its location, see for example our work [20]. These can be estimated through a power conservation low which is the basis of a mode matching technique. Explicitly, each specific source can be enclosed inside a fictitious box inside which the field is described with the aid of a numerical technique on a fine mesh exploiting the knowledge of source modeling, e.g. Volakis [16](p.238). Outside the fictitious surface the eigenfunctions (17) is considered and the weighting factors wei and wmi are evaluated through field continuity conditions across this psudo surface. It is important to keep in mind that the power conservation low should be enforced, which is preferably formulated in conjunction

The above mode matching approach pressumes the knowledge of the eigenvalues *(βi)* and the

acquired numerically. For this purpose the wave equation (6a) shall be formulated as an eigenvalue problem and the usual approach is based on the separation of both the electric field

> ˆ ˆ ˆˆ , *tt t z jz x y z x y*

The detailed procedure given in a lot of textbooks, e.g. Volakis et al [16](p.98) is summarized

*t t tz t e e e je z* ( )

Which can be substitute into (6a) to yield a pair of differential equations for the axial (ez) and

1 0 *t tt t z t rc t*

*r r*

and the nabla differential operator into axial ( *z*ˆ ) and transverse *(t)* components as:

which for the considered inhomogeneous cross-section can only be

ˆ ˆˆ ˆ , *t z tx y t z e e ze and e e x e y h h zh* (18a)

(18b)

ˆ (19)

0

(20a)

2

*e j e je k e*

*H w H w h xy e*

*E w E w e xy e*

$$\vec{\nabla}\_t \times \left(\frac{1}{\mu\_r} (\vec{\nabla}\_t e\_z + j\beta \vec{e}\_t) \times \hat{z}\right) - k\_0^2 \varepsilon\_{rc} e\_z \hat{z} = 0 \tag{20b}$$

The pair of equations (20) constitutes the eigenproblem to be solved numerically. Observing that the transverse component *te* mainly occurs a product jβ *te* , the eigenproblem in further simplified by letting *t t je e* and multiplying (20a) by jβ yields:

$$\vec{\nabla}\_t \times \left(\frac{1}{\mu\_r} \vec{\nabla}\_t \times \vec{e}\_t'\right) + \beta^2 \frac{1}{\mu\_r} \left(\vec{\nabla}\_t e\_z + \vec{e}\_t'\right) = k\_0^2 \varepsilon\_{rc} \vec{e}\_t' \tag{21a}$$

$$
\beta \,\beta^2 \vec{\nabla}\_t \times \left(\frac{1}{\mu\_r} (\vec{\nabla}\_t e\_z + \vec{e}\_t') \times \hat{z}\right) = \beta^2 k\_0^2 \varepsilon\_{rc} \cdot e\_z \hat{z} \tag{21b}
$$

Equations (21) constitute the desired eigenproblem for the *β* eigenvalue. The interested reader may contact our previous work [22] on open waveguides which was carried out within Dr. Allilomes doctoral thesis [23] also available online (in Greek).

Even though the above 2-D analysis offered some simplification the forward problem still remains complicated and computationally demanding. A truly two-dimensional approach presumes the sources to be two dimensional, i.e. tending to infinity along the z-axis, just as in the case of an "infinite line source".

#### **3.6 Two dimensional object excited by infinite line sources**

The field of an infinite electric line source embedded in a homogeneous media is well studied and can be found in any classical electromagnetics textbook, e.g. Balanis [24]( p.571). Since the line source is infinite its current density must be constant and if it is considered aligned along the z-axis, then:

Line Source:

$$\overline{J}\left(\vec{r}'\right) = I\_{2L}\hat{z} = \frac{I}{A}\mathcal{S}(\vec{\rho} - \vec{\rho}')\hat{z} \tag{22}$$

where I is a sinusoidal current (I=I0ejωt) and ρ the polar radial coordinate transverse to *z*ˆ , as A let the line source cross section assumed very small.

For a homogeneous medium the field generated by an infinite line source is a Transverse Magnetic (TMz), Hz=0 with only and axial electric field ˆ *E Ez z* and its transverse eigenfunction is identical to a cylindrical Hankel function of the second type (2) 0 *H k*( ) assuming that the line source is at ρ=0, for details in [24]( p.571). For the case of the inhomogeneous cross-section the TM nature of the wave is exactly preserved, thus ˆ *E Ez z* and Ex=Ey=0. Besides that the two-dimensional complex dielectric profile is uniform along the z-axis or ( )/ 0 *c d r dz* and consequently from equation (14) it is *Kr z* ( )/ 0 . In view of the above simplifications the third term of equation (15) vanishes as also explained by Fang [25] :

$$\begin{aligned} E\_x &= E\_y = 0\\ \left(\hat{\sigma}^2 K^2(\vec{r}) / \hat{\sigma} z = 0 \right) \end{aligned} \qquad \vec{E}(\vec{r}) \cdot \nabla^2 K^2(\vec{r}) = 0 \tag{23}$$

High to Microwave Frequencies Imaging Techniques 165

At this point recall that the Generalized Laplace equation solution domain is basically

related scalar potential (V) extends in the air to "infinity" obeying in general the Sommerfeld radiation condition. Explicitly, for very low frequencies the static potential in

0 in the air. However, as frequency (ω) is increased a significant displacement current

 *E* flows in the air around the structure and this phenomenon must be taken into account. But, for a numerical solution the analysis domain must be finite and restricted if possible to the actual body of interest. For this purpose the solution domain is enclosed within a fictitious closed surface of canonical shape on which "transparent Absorbing Boundary Conditions (ABC)" are enforced. Transparent ABCs means that they should not disturb the field distribution but conversely to behave like perfect absorbers (if possible), absorbing any wave incident on them without any reflection, just like a termination load. For the two dimensional forward mesh of Fig.5a enclosed within a fictitious circular

and the

is zero and

*c* (27)

(28a)

infinite, since for high frequency sources (i.e. in the MHz range) the electric field *E*

the air surrounding the structure vanishes since the conduction current *c J E*

*V* 1 *V ln*

on circle

The three dimensional structure of Fig.6b is enclosed within a fictitious cylindrical surface on which the appropriate ABCs should be imposed. In order to understand the ABCs in both 2-D and 3-D structures it worths to recall that these are extracted from a generalized type of boundary conditions, [15](p137). For this purpose let us write the generalized

*V VV V f*

on S1 (28b)

on S2 (28c)

*cx cy cz*

Considering the solution domain enclosed within a closed surface *SS S* 1 2 on which two

*VVV n V q*

where *n*ˆ the outward unit normal vector. Besides these the field and consequently the scalar electric potential should obey the Sommerfeld radiation condition at infinity which reads,

*x xy yz z*

ˆ *cx cy cz*

*xyz*

contour-C, the quasi-static potential ABCs read [15](p.97):

*x y* the radial polar coordinate.

Laplace (8c) or Poisson equation in the general form of:

General (mixed Dirichlet & Newmann) Boundary Condition

type of boundary conditions can be imposed as:

*V*

Dirichlet Boundary Condition:

[15](p.8) or [16]( p.8):

**Absorbing boundary conditions**

*d o J j*

where 2 2

Consequently, the wave equation (15) is reduced to a scalar Helmholtz equation:

2-D & Line Source:

$$\dot{\nabla}^2 E\_2(\vec{r}) + K^2 E\_2(\vec{r}) = j\alpha\mu\_0\mu\_r l\_{Z\perp} \tag{24}$$

A more practical structure to be considered herein is comprised by the same lossy dielectric profile *εc(x,y)* but with a finite extend in the z-direction and illuminated by a circular array of infinitesimal *z*ˆ -directed dipoles located αt the mid z-plane and encircle the object to be imaged. Even though this can again be classified as a 2-D inverse problem and solved accordingly, the corresponding forward problem is governed by the 3-D wave equation (6a) or (15) and it will be solved with the aid of a 3-D finite element method.

#### **3.7 Variational formulation**

The finite element will be employed for the numerical solution of either the generalized Laplace equation (8c), the scalar wave equation (24) or the general vector wave equation (6a). For each case the structure will be discretized with the aid of the appropriate forward fine mesh. The first step of this procedure refers to the formulation of the differential equation into an appropriate integral functional or weak formulation which will, in turn, be minimized to obtain the numerical solution. Let us start from the simpler case.

#### **3.7.1 Functional for the generalized laplace equation**

For these relatively simple scalar problems the variational formulation is usually adopted. According to the approach, e.g. Jin [15]( p.72), the differential equation (8c) along with its boundary conditions yields the functional to be minimized. These are as follows:

1. Neumann boundary conditions at the active current electrodes and along the bodysurrounding air interface:

$$-j\alpha \varepsilon\_c \frac{\mathcal{g}V}{\mathcal{g}n} = \begin{cases} \mathcal{J}\_{\prime} \text{ ondrivenelectrodes} \\ -j\alpha \varepsilon\_0 \frac{\mathcal{g}V}{\mathcal{g}n} \text{,on body-air interface} \end{cases} \tag{25}$$

where *V n* / is the potential derivative in the direction normal to the unknown object surface.

2. Mixed Dirichlet and Neumann boundary conditions between different media regionsobject elements (i,j) as:

$$V\_i = V\_j \quad \text{and} \quad \varepsilon\_{ri}^\* \frac{\mathcal{G}V\_i}{\mathcal{G}\mathfrak{n}} = \varepsilon\_{rj}^\* \frac{\mathcal{G}V\_j}{\mathcal{G}\mathfrak{n}} \tag{26}$$

These are natural boundary conditions which are an inherent property of the differential equation (8c). Namely, these are automatically imposed through the FEM solution of (8c). Note that the above does not include any Dirichlet condition, namely only the derivatives are defined (delta change), hence the resulting solution will be floating and thus nonunique. To avoid this problem at least one point must be grounded (V=0)

#### **Absorbing boundary conditions**

164 Medical Imaging

220 () () *E r KE r r ZL j*

A more practical structure to be considered herein is comprised by the same lossy dielectric profile *εc(x,y)* but with a finite extend in the z-direction and illuminated by a circular array of infinitesimal *z*ˆ -directed dipoles located αt the mid z-plane and encircle the object to be imaged. Even though this can again be classified as a 2-D inverse problem and solved accordingly, the corresponding forward problem is governed by the 3-D wave equation (6a)

The finite element will be employed for the numerical solution of either the generalized Laplace equation (8c), the scalar wave equation (24) or the general vector wave equation (6a). For each case the structure will be discretized with the aid of the appropriate forward fine mesh. The first step of this procedure refers to the formulation of the differential equation into an appropriate integral functional or weak formulation which will, in turn, be

For these relatively simple scalar problems the variational formulation is usually adopted. According to the approach, e.g. Jin [15]( p.72), the differential equation (8c) along with its

1. Neumann boundary conditions at the active current electrodes and along the body-

,on body-airinterface *c*

, ondrivenelectrodes

is the potential derivative in the direction normal to the unknown object

(26)

(25)

*J* (24)

Consequently, the wave equation (15) is reduced to a scalar Helmholtz equation:

2 2

or (15) and it will be solved with the aid of a 3-D finite element method.

minimized to obtain the numerical solution. Let us start from the simpler case.

boundary conditions yields the functional to be minimized. These are as follows:

0

*n*

2. Mixed Dirichlet and Neumann boundary conditions between different media regions-

*V V i j* and *j i*

These are natural boundary conditions which are an inherent property of the differential equation (8c). Namely, these are automatically imposed through the FEM solution of (8c). Note that the above does not include any Dirichlet condition, namely only the derivatives are defined (delta change), hence the resulting solution will be floating and thus non-

*ri rj V V n n*

*J V j V n j*

unique. To avoid this problem at least one point must be grounded (V=0)

**3.7.1 Functional for the generalized laplace equation**

2-D & Line Source:

**3.7 Variational formulation**

surrounding air interface:

object elements (i,j) as:

where *V n* /

surface.

At this point recall that the Generalized Laplace equation solution domain is basically infinite, since for high frequency sources (i.e. in the MHz range) the electric field *E* and the related scalar potential (V) extends in the air to "infinity" obeying in general the Sommerfeld radiation condition. Explicitly, for very low frequencies the static potential in the air surrounding the structure vanishes since the conduction current *c J E* is zero and 0 in the air. However, as frequency (ω) is increased a significant displacement current *d o J j E* flows in the air around the structure and this phenomenon must be taken into account. But, for a numerical solution the analysis domain must be finite and restricted if possible to the actual body of interest. For this purpose the solution domain is enclosed within a fictitious closed surface of canonical shape on which "transparent Absorbing Boundary Conditions (ABC)" are enforced. Transparent ABCs means that they should not disturb the field distribution but conversely to behave like perfect absorbers (if possible), absorbing any wave incident on them without any reflection, just like a termination load. For the two dimensional forward mesh of Fig.5a enclosed within a fictitious circular contour-C, the quasi-static potential ABCs read [15](p.97):

$$\frac{\mathcal{G}V}{\mathcal{G}\rho} \approx \frac{1}{\rho \ln \rho} V \text{ on circle } \rho = \rho\_c \tag{27}$$

where 2 2 *x y* the radial polar coordinate.

The three dimensional structure of Fig.6b is enclosed within a fictitious cylindrical surface on which the appropriate ABCs should be imposed. In order to understand the ABCs in both 2-D and 3-D structures it worths to recall that these are extracted from a generalized type of boundary conditions, [15](p137). For this purpose let us write the generalized Laplace (8c) or Poisson equation in the general form of:

$$-\frac{\partial}{\partial \mathbf{x}} \left( \varepsilon\_{cx} \frac{\partial V}{\partial \mathbf{x}} \right) - \frac{\partial}{\partial y} \left( \varepsilon\_{cy} \frac{\partial V}{\partial y} \right) - \frac{\partial}{\partial z} \left( \varepsilon\_{cz} \frac{\partial V}{\partial z} \right) + \mathcal{J}V = f \tag{28a}$$

Considering the solution domain enclosed within a closed surface *SS S* 1 2 on which two type of boundary conditions can be imposed as:

Dirichlet Boundary Condition:

$$V = \rho \text{ on } \mathbf{S}\_{\text{I}} \tag{28b}$$

General (mixed Dirichlet & Newmann) Boundary Condition

$$\left(\varepsilon\_{cx}\frac{\partial V}{\partial x} + \varepsilon\_{cy}\frac{\partial V}{\partial y} + \varepsilon\_{cz}\frac{\partial V}{\partial z}\right)\cdot\hat{n} + \gamma V = q \quad \text{on } \mathbf{S}\_2 \tag{28c}$$

where *n*ˆ the outward unit normal vector. Besides these the field and consequently the scalar electric potential should obey the Sommerfeld radiation condition at infinity which reads, [15](p.8) or [16]( p.8):

High to Microwave Frequencies Imaging Techniques 167

(a) (b) Fig. 5. The object to be imaged is surrounded by the illumintating or sensing antennas array and is enclosed within a fictitious surface "transparent" to electromagnetic waves on which

For the scalar microwave case where *ψ=Εz* or *Hz* equation (29a) yields the absorbing

With the aid of the variational technique the generalized Laplace equation (8c) along with the boundary conditions (25), (26) or their general form (28) and the corresponding ABCs in

2 2

1 1 1 () ( ) ( ) ( ) 2 *e M M e e*

*e e V V F V dxdy JVd F V x y*

where M is the total number of elements and *e* is the area of each e-element. Note that the variational formula (32) is valid for either real or complex quantities according to [15](p.76) and the references therein. For *F(V)* to be minimized, its partial derivatives with respect to

( ) 0 1,2,...,

A variational technique can be employed to formulate the microwave functional for the vector wave equation (6a) along with the appropriate boundary conditions and the

corresponding Absorbing Boundary Conditions. The resulting functional reads [15]:

*F V for i n*

(27), are found to be equivalent to the minimization of the following functional:

*i*

*V*

**3.7.3 Functional for the vector and scalar wave equations**

1 2 *z o z jk*

  (32)

(33)

(31c)

ABCs are imposed: (a) 2-D and (b) 3-D configuration.

boundary condition along a cruved fictitious surface as,

**3.7.2 Functional for generalized Laplace equation**

*c*

the elements nodal voltages must be zero, namely:

where n is the total number of nodes.

$$\text{3-D: } \lim\_{r \to \ast} r \left\{ \nabla \times \left( \vec{\Psi} \right) + jk\_o \hat{r} \times \left( \vec{\Psi} \right) \right\} = 0 \tag{29a}$$

where or *E H* and 2 22 *r x y z* ,which states that the field is comprised of an outgoing wave with dependence *o jk r e* as *r* . For a two-dimensional case (29a) where ˆ *E Ez z* (e.g. TMz modes) or ˆ *H Hz z* (for TEz modes) equation (29a) reduces to:

$$\lim\_{\varphi\_{r\to\infty}} \rho \left\{ \frac{\partial}{\partial \rho} (\Psi\_z) + j k\_o \times \left(\Psi\_z\right) \right\} = 0 \tag{29b}$$

where or *zz z E H* . The question is now how to express (29b) in the general form of (28c) to be applied on the fictitious boundary truncating the solution domain, by keeping in mind that *E V* . Additionally it has been proved that for a 2-D case the potential (the static or the quasi-static approximation) asymptotically behaves as [15](p.97):

$$\text{2-D: } V = A \ln \rho \text{ as } \; \rho \to \infty \tag{30a}$$

where A is a constant.

For a three-dimensional domain the potential in the homogenous background (i.e. air) and at a large distance from its sources behaves as [15](p150):

$$a \approx \frac{A}{r} \tag{30b}$$

where 2 22 *r x y z* .

For the two dimensional case let *V z* 0 and the 2-D ABC (27) is obtained from (28c) by letting *q=0* and *γ=1/(ρlnρ)*. Note, that inside the fictitious surface-S there is always an air layer (*εrx= εry= εrz=1*), hence (28c) reads:

$$
\left(\frac{\partial V}{\partial \mathbf{x}} + \frac{\partial V}{\partial y}\right)\hat{\mathbf{n}} = -\frac{1}{\rho \ln \rho}V\tag{31a}
$$

For a circular fictitious surface it is *n*ˆ ˆ and for a large enough radius the left hand side of (29) is reduced to *V* to yield the ABCs of (27), (the related vector identities see e.g. Balanis [24] (appendix II)

Likewise, for the three dimensional domain eq.(28c) yields the absorbing boundary conditions.

For a spherical fictitious surface with radius *rc* large enough the absorbing boundary conditions take the form:

$$\frac{\partial V}{\partial r} \approx -\frac{1}{r}V \tag{31b}$$

which, in turn, can be expressed through (28c) by letting *γ=1/r*, *q=0* and *n*ˆ ˆ .

(29a)

as *r* . For a two-dimensional case (29a) where

(29b)

(30a)

*r* (30b)

and for a large enough radius the left hand side of

(31b)

to yield the ABCs of (27), (the related vector identities see e.g.

(31a)

.

*r jk r*

where or *E H* and 2 22 *r x y z* ,which states that the field is comprised of an

lim 0 *zoz*

where or *zz z E H* . The question is now how to express (29b) in the general form of (28c) to be applied on the fictitious boundary truncating the solution domain, by keeping in mind

. Additionally it has been proved that for a 2-D case the potential (the static or

 as

For a three-dimensional domain the potential in the homogenous background (i.e. air) and

*A a*

For the two dimensional case let *V z* 0 and the 2-D ABC (27) is obtained from (28c) by letting *q=0* and *γ=1/(ρlnρ)*. Note, that inside the fictitious surface-S there is always an air

Likewise, for the three dimensional domain eq.(28c) yields the absorbing boundary

For a spherical fictitious surface with radius *rc* large enough the absorbing boundary

*V* 1

*r r*

which, in turn, can be expressed through (28c) by letting *γ=1/r*, *q=0* and *n*ˆ ˆ

*V*

*x y*

1 ˆ ln *V V n V*

*jk*

2-D: ln *V A*

(for TEz modes) equation (29a) reduces to:

 3-D: lim ˆ 0 *o r*

*r*

at a large distance from its sources behaves as [15](p150):

the quasi-static approximation) asymptotically behaves as [15](p.97):

outgoing wave with dependence *o jk r e*

(e.g. TMz modes) or ˆ *H Hz z*

ˆ *E Ez z*

that *E V*

where A is a constant.

where 2 22 *r x y z* .

(29) is reduced to *V*

Balanis [24] (appendix II)

conditions take the form:

conditions.

layer (*εrx= εry= εrz=1*), hence (28c) reads:

For a circular fictitious surface it is *n*ˆ ˆ

Fig. 5. The object to be imaged is surrounded by the illumintating or sensing antennas array and is enclosed within a fictitious surface "transparent" to electromagnetic waves on which ABCs are imposed: (a) 2-D and (b) 3-D configuration.

For the scalar microwave case where *ψ=Εz* or *Hz* equation (29a) yields the absorbing boundary condition along a cruved fictitious surface as,

$$\frac{\partial \mathcal{W}\_z}{\partial \rho} = \left(-jk\_o - \frac{1}{2\rho}\right)\mathcal{W}\_z\tag{31c}$$

#### **3.7.2 Functional for generalized Laplace equation**

With the aid of the variational technique the generalized Laplace equation (8c) along with the boundary conditions (25), (26) or their general form (28) and the corresponding ABCs in (27), are found to be equivalent to the minimization of the following functional:

$$F(V) = \frac{1}{2} \sum\_{c=1}^{M} \varepsilon\_c \int \left[ \left( \frac{\mathcal{G}V}{\mathcal{G}\mathbf{x}} \right)^2 + \left( \frac{\mathcal{G}V}{\mathcal{G}y} \right)^2 \right] d\mathbf{x} dy + \oint fV d\ell = \sum\_{c=1}^{M} F^\varepsilon(V^\varepsilon) \tag{32}$$

where M is the total number of elements and *e* is the area of each e-element. Note that the variational formula (32) is valid for either real or complex quantities according to [15](p.76) and the references therein. For *F(V)* to be minimized, its partial derivatives with respect to the elements nodal voltages must be zero, namely:

$$\frac{\partial \mathcal{F}(V)}{\partial V\_i} = 0 \quad \text{for} \quad i = 1, 2, \dots, n \tag{33}$$

where n is the total number of nodes.

#### **3.7.3 Functional for the vector and scalar wave equations**

A variational technique can be employed to formulate the microwave functional for the vector wave equation (6a) along with the appropriate boundary conditions and the corresponding Absorbing Boundary Conditions. The resulting functional reads [15]:

High to Microwave Frequencies Imaging Techniques 169

*N*

*V* or *e Ez* respectively. The scalar interpolation or basis function are given by:

*e* area of the e-th element, 1 23 23 1 2 3 1 3 2 , , *e ee ee e e e e e e a x y y x b y y cxx* , 232323 ,,,,, *eeeeee aabbcc* can be found by cyclic interchange. The above interpolation functions include unknown coefficients related to the element geometry and its material constitutive parameters. These are estimated by applying the interpolation functions at the points where the field or potential is sampled, in this case the element nodes. The resulting system of equations is analytically solved for these unknowns in terms of the potential or field values at the nodes. These are, then, substituted back into the interpolation function. In turn, this can be readily exploited in the conditions minimizing the FEM functional like that of equation (33) through an analytical evaluation of the derivatives with respect to the unknown potentials (*Vi*) or field components *Ezi*. Observe at this point that the involved surface integrals over the solution domain (e.g. (32) or (35)) are discriminated into a sum of integrals over each

> *ee e e Y VI c*

*ee e K E BJ cz z*

Within the classical FEM the above element equations are assembled to from a unified system matrix representing the whole solution domain. Hence, this matrix should include all boundary conditions along with the illuminating sources and passive substructures like sensing electrodes or antenna models as well. The latter is useful and it should be accounted at least when the methodology is matured, but it is for now neglected. Regarding the important absorbing boundary conditions these are usually imposed through the closed contour integral over ΓS in equations (32) and (35). For this purpose the contour ΓS is discretized into line segments coinciding with the corresponding edge of the peripheral

1 , 1,2,3

3

1 *e ee i i*

(38)

(39)

(40)

2

(41)

(42)

ln *s*

2 *K jk N N d s oi j s r*

 *dl*

*s o V K j*

(37)

*i*

2-D triangular:

*e ee e i ii i e N a bx cy i*

2

element. The latter yields a system of equations of the form:

The corresponding line element matrices are as follows:

Quasi-static ABCs:

2-D microwave ABCs: 1

where *e e*

or respectively:

triangles (as depicted in Fig.6)

$$\begin{aligned} \text{(3-D: } F\left(\vec{E}\right) &= \frac{1}{2} \iiint\limits\_{V} \left[\frac{1}{\mu\_{r}} (\vec{\nabla} \times \vec{E}) \cdot (\vec{\nabla} \times \vec{E}) - k\_{o}^{2} \boldsymbol{\varepsilon}\_{rc} \vec{E} \cdot \vec{E}\right] dV + \iiint\limits\_{S} \left[\frac{j k\_{o}}{2} (\hat{\boldsymbol{\tau}} \times \vec{E}) \cdot (\hat{\boldsymbol{\tau}} \times \vec{E})\right] dS\\ &+ \iiint\limits\_{V} \underbrace{j o \mu\_{o} \hat{\boldsymbol{\delta}}\_{s} \cdot \vec{E}}\_{\text{source}} dV \end{aligned} \tag{34}$$

Likewise, the scalar wave equation (24) when a 2-D structure is illuminated by a line source results to the functional:

$$\text{2-D: } F(\vec{E}) = \frac{1}{2} \iiint\_{S} \left[ \left( \frac{\partial E\_z}{\partial x} \right)^2 + \left( \frac{\partial E\_z}{\partial y} \right)^2 - k\_o^2 \varepsilon\_{rc} E\_z \right] dS + \frac{1}{2} \underbrace{\oint\_{\Gamma\_\*} \underbrace{\left(jk\_o + \frac{1}{2r}\right)}\_{A \text{BCs}} E\_z^2 dl}\_{A \text{BCs}} + \iint\_S \underbrace{\left[ \alpha \mu\_o I\_z E\_z \right]}\_{source} dS \quad \text{(35)}$$

The surface integral in eq.(34) is defined over the fictitious cylinder truncating the solution domain. Likewise the line integral in eq. (35) is defined along a fictitious circle (ΓS) terminating the 2-D mesh. Both these integrals express the absorbing boundary condition.

#### **3.7.4 Finite element solution of the forward problem**