**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, μ<sup>r</sup> 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 (σ<sup>i</sup> 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, μ<sup>r</sup> 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 (σ<sup>i</sup>

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 ε<sup>r</sup> 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 ε<sup>r</sup> 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 *<sup>r</sup>* 

2

r 0

0 r 0

(6a)

<sup>1</sup> x -k <sup>j</sup> *imp xE <sup>j</sup> E J* 

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

*<sup>E</sup>* . 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* <sup>0</sup> *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 *<sup>j</sup> <sup>t</sup> 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 ( ) *<sup>c</sup> J* flowing at any arbitrary point (movement of charges due to coulomb forces and the presence of electric field) as *<sup>c</sup> J E* . Thus, the current density *J* appearing in (2b) is comprised of the conduction current *<sup>c</sup> J* and the current impressed by the source as *imp <sup>J</sup>* . 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

*<sup>j</sup> <sup>t</sup> tr Er e* , Re

*<sup>j</sup> <sup>t</sup> tr H r e* , Re

*<sup>j</sup> <sup>t</sup> tr J r e* , Re

0 r xE -j <sup>H</sup>

0 r xH j E J 

 D E 

*BH H* 

> *j*

Additionally, for the conducting media assumed herein there is a conduction current ( ) *<sup>c</sup> 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 <sup>c</sup> 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 <sup>J</sup>*

 

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 *<sup>j</sup> <sup>t</sup> 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 *<sup>c</sup> J E*

complex conductivity *σc* as:

comprised of the conduction current *<sup>c</sup> 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* <sup>0</sup> *J Et j <sup>E</sup>* . 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 *<sup>r</sup>* 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 )

*<sup>D</sup> <sup>D</sup> <sup>J</sup> <sup>c</sup> <sup>t</sup> t t*

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

 *<sup>c</sup>* 

*<sup>J</sup> V V j jj <sup>V</sup> n n j*

  (3) 0 0

 ( )0 

 0 *<sup>c</sup>* 

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 *<sup>c</sup> 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*

   

*<sup>D</sup> <sup>J</sup> or J <sup>t</sup> <sup>t</sup>*

and *D E* 

 

conservation equation as:

to time harmonic fields we get:

where:

Newmann type:

Substituting the constitutive relations *J E*

 

*E EV* <sup>0</sup> (7a)

(7b)

 

(8a)

into (8a) and restricting ourselves

(10)

,

*j E* (8b)

*r V* (8c)

*r jr* (9)

on the electrode

 

*<sup>E</sup> <sup>j</sup> <sup>V</sup> <sup>j</sup>* (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 j<sup>H</sup>* produce an electric field with an identical time variation through (2a). This / *E t j<sup>E</sup>* 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ε<sup>r</sup> 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 <sup>1</sup> x -k *<sup>J</sup> 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ε<sup>r</sup> 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*

*<sup>H</sup>*

*<sup>E</sup>*

*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 <sup>r</sup> 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. Let us give a formal description of such a methodology. For each one of the possibly excited modes the field dependence in the *z*ˆ -

(,,) (,) *<sup>i</sup> <sup>j</sup> <sup>z</sup> E xyz e xye i i*

(,,) (,) *<sup>i</sup> <sup>j</sup> <sup>z</sup> 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 <sup>j</sup>*

*r r r imp*

 

 

( )

0

(15)

2

2 0

 

(13)

 

(14)

*r xy* . Due to the uniformity in the axial *z*<sup>ˆ</sup> -

(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

<sup>2</sup> *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*<sup>ˆ</sup> 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. Let us give a formal description of such a methodology. 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

<sup>1</sup> () 0 ˆ ˆ *<sup>t</sup> tz t rc z*

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

mainly occurs a product j<sup>β</sup> *te*

*t tt t z t rc t*

<sup>1</sup> ( ) ˆ ˆ *<sup>t</sup> 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

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

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 <sup>z</sup>*

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 <sup>z</sup>*

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 *<sup>c</sup> 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

<sup>0</sup> *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 ( *<sup>z</sup>*<sup>ˆ</sup> ) 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)

(.) *<sup>i</sup> <sup>j</sup> <sup>z</sup>*

(.) *<sup>i</sup> <sup>j</sup> <sup>z</sup>*

(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

<sup>1</sup> <sup>0</sup> *t tt t z t rc t*

 

*r r*

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

which for the considered inhomogeneous cross-section can only be

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

(18b)

<sup>ˆ</sup> (19)

0

(20a)

 

<sup>2</sup>

*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<sup>β</sup> *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 <sup>z</sup>* and its transverse eigenfunction is identical to a cylindrical Hankel function of the second type (2) <sup>0</sup> *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 <sup>z</sup>* and Ex=Ey=0. Besides that the two-dimensional complex dielectric profile is uniform along the z-axis or ( )/ 0 *<sup>c</sup> 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

 *<sup>E</sup>* 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

*<sup>c</sup>* (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 *<sup>c</sup> J E*

*<sup>V</sup>* <sup>1</sup> *<sup>V</sup> 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 <sup>V</sup> <sup>f</sup>*

 

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 <sup>q</sup>*

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

<sup>220</sup> () () *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 *<sup>c</sup>*

, ondrivenelectrodes

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

(26)

(25)

 

*<sup>J</sup>* (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 <sup>i</sup> <sup>j</sup>* and *<sup>j</sup> <sup>i</sup>*

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*

 

*<sup>J</sup> <sup>V</sup> <sup>j</sup> <sup>V</sup> 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 *<sup>c</sup> J E* is zero and 0 in the air. However, as frequency (ω) is increased a significant displacement current *d o J j <sup>E</sup>* 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 <sup>1</sup> () ( ) ( ) ( ) <sup>2</sup> *<sup>e</sup> 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 *<sup>e</sup>* 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 *<sup>z</sup> 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 <sup>y</sup> <sup>z</sup>* ,which states that the field is comprised of an outgoing wave with dependence *<sup>o</sup> jk r e* as *r* . For a two-dimensional case (29a) where ˆ *E Ez <sup>z</sup>* (e.g. TMz modes) or ˆ *H Hz <sup>z</sup>* (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)

*<sup>r</sup>* (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 <sup>y</sup> <sup>z</sup>* ,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*

 

<sup>1</sup> <sup>ˆ</sup> ln *V V n V*

 

*jk*

2-D: ln *V A*

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

 3-D: lim ˆ 0 *<sup>o</sup> 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 *<sup>o</sup> jk r e*

(e.g. TMz modes) or ˆ *H Hz <sup>z</sup>*

ˆ *E Ez <sup>z</sup>*

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 *<sup>e</sup>* 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 *<sup>e</sup> Ez* respectively. The scalar interpolation or basis function are given by:

*<sup>e</sup>* 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 <sup>c</sup>*

*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

<sup>1</sup> , 1,2,3

3

1 *e ee i i*

(38)

(39)

(40)

2

(41)

(42)

ln *<sup>s</sup>*

<sup>2</sup> *K jk N N d s oi <sup>j</sup> <sup>s</sup> <sup>r</sup>*

 *dl* 

*s o <sup>V</sup> K j*

(37)

*i*

2-D triangular:

*e ee e i ii i <sup>e</sup> 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: <sup>1</sup>

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**

For both the quasi static and the microwave approach the structure must be discretized utilizing a dual mesh. A fine one for the forward problem and a coarse mesh for the inverse problem, as shown in Fig.3. The coarse mesh is defined only over the object to be imaged (Fig.6), while the fine-forward mesh covers the whole solution domain up to the fictitious line or surface enclosing all constituents including antennas of electrodes. Besides that every coarse element is subdivided in a number of smaller fine mesh elements. The nodal FEM approach is employed for the MHz potential estimation as well as the 2-D scalar Helmholtz equation involving only Ez component referred bellow as scalar FEM. Conversely, the edge elements technique is adopted for the 3-D vector wave equation. Since, for the reconstruction mesh the body under consideration is split into cubic (or rectangular) elements with constant σ and εr, so a piecewise homogeneous model is constructed as:

$$
\varepsilon\_c \left( x, y \right) = \sum\_{k=1}^{E} \varepsilon\_{ck} \nu\_{k'} \qquad \qquad \nu\_k = \begin{cases} 1 & \text{within k-the element} \\ 0 & \text{elsewhere} \end{cases} \tag{36}
$$

Even though FEM is well described in numerous textbooks, e.g. [15], [16] and [17] for multigrid approaches, a short description of the basic interpolation functions and the related element matrices is given next, since they are necessary for the definition and evaluation of the Jacobian-Sensitivity matrix.

#### **3.7.5 Scalar nodal 2-D approach**

Either the potential in the functional (32) or the Ez component of (35) are discretized employing first order linear triangular elements with the unknowns (degrees of freedom) defined on their nodes as:

$$\text{2-D triangular: } \boldsymbol{\psi}^{\varepsilon} = \sum\_{i=1}^{3} N\_i^{\varepsilon} \boldsymbol{\nu}\_i^{\varepsilon} \tag{37}$$

where *e e V* or *<sup>e</sup> Ez* respectively. The scalar interpolation or basis function are given by:

$$N\_i^{\varepsilon} = \frac{1}{2\Delta^{\varepsilon}} \left( a\_i^{\varepsilon} + b\_i^{\varepsilon} \mathbf{x} + c\_i^{\varepsilon} \mathbf{y} \right), \quad i = 1, 2, 3 \tag{38}$$

*<sup>e</sup>* 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 element. The latter yields a system of equations of the form:

$$
\mathbb{E}\left[\left.Y^{\varepsilon}\left(\mathcal{E}\_{\varepsilon}^{\varepsilon}\right)\right]\right]\left[\left.Y^{\varepsilon}\right] = \left[\begin{array}{c}I^{\varepsilon} \\ \end{array}\right] \tag{39}
$$

or respectively:

168 Medical Imaging

*jk F E E E k E E dV r E r E dS*

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

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

For both the quasi static and the microwave approach the structure must be discretized utilizing a dual mesh. A fine one for the forward problem and a coarse mesh for the inverse problem, as shown in Fig.3. The coarse mesh is defined only over the object to be imaged (Fig.6), while the fine-forward mesh covers the whole solution domain up to the fictitious line or surface enclosing all constituents including antennas of electrodes. Besides that every coarse element is subdivided in a number of smaller fine mesh elements. The nodal FEM approach is employed for the MHz potential estimation as well as the 2-D scalar Helmholtz equation involving only Ez component referred bellow as scalar FEM. Conversely, the edge elements technique is adopted for the 3-D vector wave equation. Since, for the reconstruction mesh the body under consideration is split into cubic (or rectangular) elements with constant σ and εr, so a piecewise homogeneous model

1 withink-th element , , 0 elsewhere

Even though FEM is well described in numerous textbooks, e.g. [15], [16] and [17] for multigrid approaches, a short description of the basic interpolation functions and the related element matrices is given next, since they are necessary for the definition and evaluation of

Either the potential in the functional (32) or the Ez component of (35) are discretized employing first order linear triangular elements with the unknowns (degrees of freedom)

(36)

 

*<sup>o</sup> o rc*

*z z o rc z o z z z*

*S S source*

*E E F E k E dS jk E dl <sup>j</sup> J E dS x y r*

*ABCs*

*ABCs*

(35)

(34)

1 1 <sup>2</sup> 3-D: ˆ ˆ 2 2

ˆ

*s V source*

*j J EdV* 

*V S r*

<sup>2</sup> 2 2 *<sup>s</sup>*

 

2 2 <sup>1</sup> 2 2 1 1 2-D:

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

the Jacobian-Sensitivity matrix.

defined on their nodes as:

**3.7.5 Scalar nodal 2-D approach** 

1

 

 *x y* 

*E c ck k k k*

results to the functional:

condition.

is constructed as:

$$\mathbb{E}\left[\mathcal{K}^{\varepsilon}\left(\mathcal{E}\_{c}^{\varepsilon}\right)\right]\bigg|\,E\_{z}^{\varepsilon}\right] = \mathbb{E}\left[B\left(f\_{z}\right)\right] \tag{40}$$

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 triangles (as depicted in Fig.6)

The corresponding line element matrices are as follows:

$$\text{Quasi-static ABCs: } K\_s = -j\alpha \varepsilon\_o \int\_{\vec{\text{err}}\_s} \frac{V^2}{\rho \ln \rho} dl \tag{41}$$

$$\text{2-D microwave ABC.} \left[ K\_s \right] = \oint\_{\partial \Omega} \left( jk\_o + \frac{1}{2r} \right) \mathbf{N}\_i \mathbf{N}\_j d\Gamma\_s \tag{42}$$

High to Microwave Frequencies Imaging Techniques 171

6

1 *e ee i i*

1 2 12 21

The edge numbers and the associated nodes *i1* and *i2* are defined in Table 1 and in Fig.7. For

Edge *i* Node *i1* Node *i2* 1 1 2 2 1 3 3 1 4 4 2 3 5 4 2 6 3 4

(a) (b) Fig. 7. Finite elements employed for the discretization, (a) tetrahedral edge element for 3-D

Following an approach similar (but more complicated) to that for the 2-D case the resulting

*i j oc i j e*

(47)

*K N N k N N dV*

vector-tetrahedral: <sup>1</sup> <sup>2</sup> *<sup>T</sup> <sup>e</sup> e e ee e*

*V r*

a detailed definition of the quantities involved in (46) one may consult [15, 17]

*e e e e e ee N Wl L L L Ll i ii i i i i i i* (46)

(45)

is the vector interpolation

*i E NE* 

3-D:

Table 1. Edge definition for a tetrahedral element

and (b) triangular node element for 2-D.

volume element matrix reads:

or basis function given by:

*<sup>e</sup> Ei* denotes the tangential field values along the i-th edge and *<sup>e</sup> Ni*

Fig. 6. Example of coarse mesh overlapping a fine mesh, for 2-D rectangular object enclosed in a fictitious circle *(ΓS).* The contour *ΓS* is discretized using line segments which coincide with the corresponding edge of the peripheral triangles

The last term to be included here is the feeding sources contribution through the corresponding terms of equations (32) and (35). Herein the simplest possible sources are considered, i.e. infinitesimal electrodes or dipoles respectively. In both cases the current density is approximated as uniform ignoring the field/current singularity at the edges (a phenomenon which should be included in a most complete version), thus given by the ratio of the feeding current (I) to the antenna or electrode cross section (A) like *J=I/A*. The related integrals over the source yield the right hand side *[I]* or *B(Jz)* respectively. The resulting FEM system of equations takes the form:

$$\text{Quasi static FEM: } \left[ \begin{array}{c} Y(\mathcal{E}\_c) \end{array} \right] \left[ V \right] = \left[ I \right] \tag{43}$$

$$\text{2-D microowave FEM:} \left[ K\left(\varepsilon\_{c}\right) \right] \left[ E\_{z}\left(\varepsilon\_{c}\right) \right] = \left[ B\left(f\_{z}\right) \right] \tag{44}$$

Classical numerical techniques can be employed to solve systems (43) or (44), but it is useful to keep in mind that multiple solutions are necessary (for each inverse problem iteration), one for each illuminating antenna or active electrodes pair position. Hence a technique based on the inversion of the system matrix like the LU decomposition is very convenient, since the solution for each kth-right hand side is obtained by a simple multiplication as *k k* <sup>1</sup> *V YI* .

#### **3.7.6 Vector edge elements FEM**

For the volume integral of the electric field vector in the functional of (34) first order tetrahedral vector-edge elements are employed. The electric field *E* within each tetrahedral is expanded in terms of the FEM basis functions as, [15]:

Fig. 6. Example of coarse mesh overlapping a fine mesh, for 2-D rectangular object enclosed in a fictitious circle *(ΓS).* The contour *ΓS* is discretized using line segments which coincide

The last term to be included here is the feeding sources contribution through the corresponding terms of equations (32) and (35). Herein the simplest possible sources are considered, i.e. infinitesimal electrodes or dipoles respectively. In both cases the current density is approximated as uniform ignoring the field/current singularity at the edges (a phenomenon which should be included in a most complete version), thus given by the ratio of the feeding current (I) to the antenna or electrode cross section (A) like *J=I/A*. The related integrals over the source yield the right hand side *[I]* or *B(Jz)* respectively. The resulting FEM

Quasi static FEM: *Y VI*

2-D microwave FEM: *K E BJ* 

Classical numerical techniques can be employed to solve systems (43) or (44), but it is useful to keep in mind that multiple solutions are necessary (for each inverse problem iteration), one for each illuminating antenna or active electrodes pair position. Hence a technique based on the inversion of the system matrix like the LU decomposition is very convenient, since the solution for each kth-right hand side is obtained by a simple multiplication as

For the volume integral of the electric field vector in the functional of (34) first order

tetrahedral vector-edge elements are employed. The electric field *E*

is expanded in terms of the FEM basis functions as, [15]:

 

*<sup>c</sup>* (43)

*c zc z* (44)

within each tetrahedral

with the corresponding edge of the peripheral triangles

system of equations takes the form:

**3.7.6 Vector edge elements FEM** 

*k k* <sup>1</sup> *V YI* .

$$\text{3-D: } \vec{E}^{\epsilon} = \sum\_{i=1}^{6} \vec{N}\_{i}^{\epsilon} E\_{i}^{\epsilon} \tag{45}$$

*<sup>e</sup> Ei* denotes the tangential field values along the i-th edge and *<sup>e</sup> Ni* is the vector interpolation or basis function given by:

$$
\vec{N}\_i^e = \mathcal{V} \mathbf{V}\_{i\_1 i\_2} \mathbf{I}\_i^e = \left( L\_{i\_1}^e \vec{\nabla} L\_{i\_2}^e - L\_{i\_2}^e \vec{\nabla} L\_{i\_1}^e \right) \mathbf{I}\_i^e \tag{46}
$$

The edge numbers and the associated nodes *i1* and *i2* are defined in Table 1 and in Fig.7. For a detailed definition of the quantities involved in (46) one may consult [15, 17]


Table 1. Edge definition for a tetrahedral element

Fig. 7. Finite elements employed for the discretization, (a) tetrahedral edge element for 3-D and (b) triangular node element for 2-D.

Following an approach similar (but more complicated) to that for the 2-D case the resulting volume element matrix reads:

$$\text{vector-tetrahedron:}\\\left[\text{ K}^{\varepsilon}\right] = \iiint\limits\_{V} \left[\frac{1}{\mu\_{r}^{\varepsilon}} \left\{\vec{\nabla} \times \vec{N}\_{i}^{\varepsilon}\right\} \cdot \left\{\vec{\nabla} \times \vec{N}\_{j}^{\varepsilon}\right\}^{T} - k\_{o}^{2} \varepsilon\_{c}^{\varepsilon} \left\{\vec{N}\_{i}^{\varepsilon}\right\} \cdot \left\{\vec{N}\_{j}^{\varepsilon}\right\} \right] dV \tag{47}$$

High to Microwave Frequencies Imaging Techniques 173

(a) (b) Fig. 9. (a) Electric field distribution when a 2D structure is illuminated by an infinitesimal dipole operated at 800MHz and (b) electric field distribution at two planes when a 3D

The reconstruction algorithm is based on the modified perturbation method that was initially developed for the conductivity imaging in Electrical Impedance Tomography [26] and later on extended to higher frequencies (up to 10MHz),[27]. This research was carried out within the doctoral thesis of D. Drogoudis [28], which was primarily focused to extend MPM and prove its validity in the microwave regime. Indeed this proof of concept was achieved and published in [29] and [30]. For the introduction of the reader to the inverse problem approach but also to provide him with the knowledge to proceed further and ultimately achieve its practical implementation, the reasoning behind its implementation and the logical steps for its extension will be given next. Thus, let us start by reviewing the static MPM and proceed to its extension to the MHz and finally the microwave range.

First a review of the static algorithm will be given and based on that the steps towards its present time harmonic formulation will be described. According to [26] the σ-imaging

> *n n i mi j n j j M j*

1

*k*

 *<sup>j</sup>* or / *Vk*

 

*Vi*  1 1 1

*k k j* *<sup>j</sup>* constitute the elements of the Jacobian matrix.

(49)

*V*

*<sup>M</sup> mi ci <sup>i</sup>*

*V V V V*

> 1 | |

where M is the total number of linearly independent measurements, *Vmi* and Vci are the measured and calculated voltage differences at the ith port (electrodes pair) and k1 is the relaxation factor to be chosen in the range 0<k1<2 in order to provide faster convergence.

These are calculated from the voltage distributions all over the object cross section through a

structure is illuminated by an infinitesimal dipole operated at 1 GHz.

**4. Inverse problem** 

algorithm reads:

The partial derivatives /

**4.1 Review of the static MPM algorithm** 

The absorbing boundary conditions stem from (29a) and take the form of the corresponding term in (34) through an asymptotic approximation. In order to evaluate this integral, the closed fictitious surface (S) is discretized into triangles coinciding with the outer surface of the peripheral tetrahedrals as shown if Fig.8. The resulting element matrix reads:

$$\text{G-D ABCs:} \left[ K^s \right] = \bigotimes\_{S} \left[ \frac{j k\_o}{2} \left( \hat{r} \times \vec{N}\_i \right) \cdot \left( \hat{r} \times \vec{N}\_j \right) \right] dS \tag{48}$$

where *r*ˆ is the outward unit normal vector. Note that the above involved triangular edge elements shape functions are not trivial, since each triangle has a different orientation in a 3- D space. The most convenient way to extract them is through the degeneration of the corresponding peripheral tetrahedral into a surface triangle. Regarding the source modeling an infinitesimal dipole approximation is again considered as for the 2-D case.

Fig. 8. Triangles coinciding with the outer surface of the peripheral tetrahedrals

The linear system solution in the 3-D case constitutes a major difficulty due to the large number of unknowns. Thus the performance of direct inversion methods like LU decomposition becomes questionable and iterative techniques are also tried herein. Specifically, in the 3-D case the inversion of the matrix K will consume all the system memory due to the size of the matrix, so an iterative solver is preferred with the appropriate preconditioner. The Generalized Minimum Residual method (GMRES) is used for the solution of the linear system with a symmetric successive over-relaxation-vector (SSORvector) preconditioner.

Solving this system the electric fields (or the potential) on the receiving antennas (or electrodes) and at all the internal edges or nodes is calculated and stored, to be used latter within the reconstruction algorithm. An example of the electric field distribution when one of the infinitesimal dipole or a line-source antenna is activated is shown in Fig. 9.

Fig. 9. (a) Electric field distribution when a 2D structure is illuminated by an infinitesimal dipole operated at 800MHz and (b) electric field distribution at two planes when a 3D structure is illuminated by an infinitesimal dipole operated at 1 GHz.

#### **4. Inverse problem**

172 Medical Imaging

The absorbing boundary conditions stem from (29a) and take the form of the corresponding term in (34) through an asymptotic approximation. In order to evaluate this integral, the closed fictitious surface (S) is discretized into triangles coinciding with the outer surface of

the peripheral tetrahedrals as shown if Fig.8. The resulting element matrix reads:

an infinitesimal dipole approximation is again considered as for the 2-D case.

Fig. 8. Triangles coinciding with the outer surface of the peripheral tetrahedrals

vector) preconditioner.

The linear system solution in the 3-D case constitutes a major difficulty due to the large number of unknowns. Thus the performance of direct inversion methods like LU decomposition becomes questionable and iterative techniques are also tried herein. Specifically, in the 3-D case the inversion of the matrix K will consume all the system memory due to the size of the matrix, so an iterative solver is preferred with the appropriate preconditioner. The Generalized Minimum Residual method (GMRES) is used for the solution of the linear system with a symmetric successive over-relaxation-vector (SSOR-

Solving this system the electric fields (or the potential) on the receiving antennas (or electrodes) and at all the internal edges or nodes is calculated and stored, to be used latter within the reconstruction algorithm. An example of the electric field distribution when one

of the infinitesimal dipole or a line-source antenna is activated is shown in Fig. 9.

3-D ABCs: ˆ ˆ <sup>2</sup>

where *r*ˆ is the outward unit normal vector. Note that the above involved triangular edge elements shape functions are not trivial, since each triangle has a different orientation in a 3- D space. The most convenient way to extract them is through the degeneration of the corresponding peripheral tetrahedral into a surface triangle. Regarding the source modeling

*S*

*<sup>s</sup> <sup>o</sup> i j*

*jk <sup>K</sup> r N r N dS*

(48)

The reconstruction algorithm is based on the modified perturbation method that was initially developed for the conductivity imaging in Electrical Impedance Tomography [26] and later on extended to higher frequencies (up to 10MHz),[27]. This research was carried out within the doctoral thesis of D. Drogoudis [28], which was primarily focused to extend MPM and prove its validity in the microwave regime. Indeed this proof of concept was achieved and published in [29] and [30]. For the introduction of the reader to the inverse problem approach but also to provide him with the knowledge to proceed further and ultimately achieve its practical implementation, the reasoning behind its implementation and the logical steps for its extension will be given next. Thus, let us start by reviewing the static MPM and proceed to its extension to the MHz and finally the microwave range.

#### **4.1 Review of the static MPM algorithm**

First a review of the static algorithm will be given and based on that the steps towards its present time harmonic formulation will be described. According to [26] the σ-imaging algorithm reads:

$$
\sigma\_j^n = \sigma\_j^{n-1} + k\_1 \frac{\sum\_{i=1}^M \frac{V\_{mi} - V\_{ci}}{V\_{mi}} \frac{\mathcal{G} V\_i}{\mathcal{G} \sigma\_j}}{\sum\_{k=1}^M |\frac{\mathcal{G} V\_k}{\mathcal{G} \sigma\_j}|} \sigma\_j^{n-1} \tag{49}
$$

where M is the total number of linearly independent measurements, *Vmi* and Vci are the measured and calculated voltage differences at the ith port (electrodes pair) and k1 is the relaxation factor to be chosen in the range 0<k1<2 in order to provide faster convergence. The partial derivatives / *Vi <sup>j</sup>* or / *Vk <sup>j</sup>* constitute the elements of the Jacobian matrix. These are calculated from the voltage distributions all over the object cross section through a

High to Microwave Frequencies Imaging Techniques 175

(a) (b)

(c) (d)

Working on the extension of the algorithm toward complex permittivity imaging, equations

*Analyzing the constituents of MPM algorithm:* Let us clarify equation (49) starting from each term in the numerator sum which runs over all i=1, to M voltage measurement ports. The ith term corresponds to the i-th port where the difference between measured , ( ) *Vm i* and calculated , ( ) *Vc i* voltages , , ( ) *V V mi ci* is normalized by the measured value and this term is multiplied by the sensitivity of the i-th port with respect to the j-th pixel conductivity

> 

(49) it is just a term normalizing the sensitivities to unity. It is exactly the same as saying that that the infinite integral of a probability density function is equal to unity. From a different point of view, each measurement port contribute to the total updating summation by its normalized voltage deviation from the measured value weighted by the port's normalized sensitivity. The same principle can be readily applied to the complex or two variables

update (this could be even generalized to multiple variables) scheme, provided that

. As for the denominator sum in

Fig. 10. Equivalent resistor network for each jth element with locally homogenous conductivity *σ.*: (a) Triangular element defined by its coordinates *(xi,yi)*, (b) equivalent network of (a) (c) Rectangular element of dimensions a×b, (d) equivalent network of (c).

(49) and (50) should be modified for complex voltages and complex admittances.

**4.2 MPM extension to high frequencies** 

the corresponding sensitivities are available.

*<sup>j</sup>* (being currently updated), namely / / *V V <sup>i</sup> <sup>j</sup> <sup>i</sup> <sup>j</sup>*

(,) *e e j rj* 

closed form expression given by Yorkey et al. [31]. The latter is based on the electrical circuit compensation theorem as applied on a resistor network equivalent for each σ-element, as shown in Fig.10. In turn these expressions are given in [31] as:

$$J\_{ij} = \frac{\mathcal{G}V\_i}{\mathcal{G}\sigma\_j} = -\frac{1}{I\_i} \sum\_{\ell=1}^{L} S\_{\ell} V\_{ij\{\ell\}} V\_{kj\{\ell\}} \tag{50}$$

where*Vij*( ) and *Vkj*( ) are the voltages developed across the *lth* branch of the *jth* element Fig.10 when a sinusoidal current with amplitude *i k I I* is successively injected through the ports -i and -k respectively. The constants *S* are the normalized geometrical weights arising from the finite element formulation. *S* are actually given by the non-diagonal entries of the *<sup>e</sup> Y* -element matrix as: *<sup>j</sup> <sup>e</sup> <sup>e</sup> Y K* and *<sup>e</sup> S* non-diagonal entries of*Ke*.

An indicative example of these branch admittance values is through the definition of the resistance network for the rectangular element (Fig.10) and its related element matrix:

*static triangular element* 

$$\begin{bmatrix} Y \end{bmatrix}\_3 = \begin{bmatrix} G\_{11} & -G\_{12} & -G\_{13} \\ & G\_{22} & -G\_{23} \\ \text{symmetric} & & G\_{33} \end{bmatrix} = \sigma\_f \begin{bmatrix} K \end{bmatrix}\_3 \tag{50a}$$

$$\begin{bmatrix} \mathbf{K} \end{bmatrix}\_3 = \frac{1}{4\Delta} \begin{bmatrix} b\_1^2 + c\_1^2 & b\_1 b\_2 + c\_1 c\_2 & b\_1 b\_3 + c\_1 c\_3 \\ & b\_2^2 + c\_2^2 & b\_2 b\_3 + c\_2 c\_3 \\ \text{symmetric} & & b\_3^2 + c\_1^3 \end{bmatrix} \tag{50b}$$

*static rectangular element* 

$$\begin{bmatrix} Y \end{bmatrix}\_4 = \begin{bmatrix} \mathbf{G}\_{11} & -\mathbf{G}\_{12} & -\mathbf{G}\_{13} & -\mathbf{G}\_{14} \\ & \mathbf{G}\_{22} & -\mathbf{G}\_{23} & -\mathbf{G}\_{24} \\ & & \mathbf{G}\_{33} & -\mathbf{G}\_{34} \\ symmetric & & & \mathbf{G}\_{44} \end{bmatrix} \tag{50c}$$

$$\begin{bmatrix} \begin{matrix} K \end{matrix} \end{bmatrix}\_{4} = \begin{bmatrix} 2a^2 + 2b^2 & a^2 - 2b^2 & -a^2 - b^2 & -a^2 - b^2 \\\\ 2a^2 + 2b^2 & -a^2 - b^2 & -a^2 - b^2 \\\\ 2a^2 + 2b^2 & a^2 - 2b^2 \\\\ \text{symmetric} & 2a^2 + 2b^2 \end{bmatrix} \tag{50d}$$

$$\text{and}\quad \Delta = \frac{1}{2} \begin{vmatrix} 1 & x\_1 & y\_1 \\ 1 & x\_2 & y\_2 \\ 1 & x\_3 & y\_3 \end{vmatrix} = \frac{1}{2} (b\_1 c\_2 - b\_2 c\_1), \begin{vmatrix} a\_1 = x\_2 y\_3 - y\_2 x\_3, \ b\_1 = y\_2 - y\_3, c\_1 = x\_3 - x\_2 \\ a\_2 = x\_3 y\_1 - y\_3 x\_1, \ b\_2 = y\_3 - y\_1, c\_2 = x\_1 - x\_3 \\ a\_3 = x\_1 y\_2 - y\_1 x\_2, \ b\_3 = y\_1 - y\_2, c\_3 = x\_2 - x\_1 \end{vmatrix}$$

A detailed derivation of the Compensation theorem can be found in Yorkey et al. [31] as well as in Sahalos et al. [32](p.229). Guidelines for its implementation are also given therein.

closed form expression given by Yorkey et al. [31]. The latter is based on the electrical circuit compensation theorem as applied on a resistor network equivalent for each σ-element, as

1 *<sup>L</sup>*

where*Vij*( ) and *Vkj*( ) are the voltages developed across the *lth* branch of the *jth* element Fig.10 when a sinusoidal current with amplitude *i k I I* is successively injected through the ports -i and -k respectively. The constants *S* are the normalized geometrical weights arising from the finite element formulation. *S* are actually given by the non-diagonal

An indicative example of these branch admittance values is through the definition of the resistance network for the rectangular element (Fig.10) and its related element matrix:

> 11 12 13 3 3 22 23

*K b c bb cc*

*G GG Y GG K symmetric G*

 

2 2

33

1 1 12 12 13 13 2 2 2 2 23 23 3

11 12 13 14

*GGGG*

*symetric G*

4 2222

A detailed derivation of the Compensation theorem can be found in Yorkey et al. [31] as well as in Sahalos et al. [32](p.229). Guidelines for its implementation are also given therein.

2 2

2 2 2 2 22 22

*a b a b ab ab*

*symetric a b*

22 23 24

2 2 22 22

*a b ab ab*

*GGG*

33 34

22 2

*abab*

*G G*

44

2 2

2 2

1 23 23 1 2 3 1 3 2 2 31 31 2 3 1 2 1 3 3 12 12 3 1 2 3 2 1

*a xy yx b y y c x x a xy yx b y y c x x a xy yx b y y c x x*

 

, , , , , ,

*b c bb cc bb cc*

*symmetric b c*

 *j*

2 3 3 1

*ij ij kj j i <sup>V</sup> <sup>J</sup> SV V I*

 

*i*

1

() ()

(50)

(50a)

(50b)

(50c)

(50d)

and *<sup>e</sup> S* non-diagonal entries of*Ke*.

shown in Fig.10. In turn these expressions are given in [31] as:

entries of the *<sup>e</sup> Y* -element matrix as: *<sup>j</sup> <sup>e</sup> <sup>e</sup> Y K*

and 1 1

1 1 <sup>1</sup> 2 2

*x y* 

*x y*

3 3

2 2 12 21

*x y bc bc*

*K*

1

1

1 4

*Y*

4

22 2

*static triangular element* 

*static rectangular element* 

Fig. 10. Equivalent resistor network for each jth element with locally homogenous conductivity *σ.*: (a) Triangular element defined by its coordinates *(xi,yi)*, (b) equivalent network of (a) (c) Rectangular element of dimensions a×b, (d) equivalent network of (c).

#### **4.2 MPM extension to high frequencies**

Working on the extension of the algorithm toward complex permittivity imaging, equations (49) and (50) should be modified for complex voltages and complex admittances.

*Analyzing the constituents of MPM algorithm:* Let us clarify equation (49) starting from each term in the numerator sum which runs over all i=1, to M voltage measurement ports. The ith term corresponds to the i-th port where the difference between measured , ( ) *Vm i* and calculated , ( ) *Vc i* voltages , , ( ) *V V mi ci* is normalized by the measured value and this term is multiplied by the sensitivity of the i-th port with respect to the j-th pixel conductivity *<sup>j</sup>* (being currently updated), namely / / *V V <sup>i</sup> <sup>j</sup> <sup>i</sup> <sup>j</sup>* . As for the denominator sum in (49) it is just a term normalizing the sensitivities to unity. It is exactly the same as saying that that the infinite integral of a probability density function is equal to unity. From a different point of view, each measurement port contribute to the total updating summation by its normalized voltage deviation from the measured value weighted by the port's normalized sensitivity. The same principle can be readily applied to the complex or two variables (,) *e e j rj* update (this could be even generalized to multiple variables) scheme, provided that the corresponding sensitivities are available.

High to Microwave Frequencies Imaging Techniques 177

(a) (b)

Fig. 11. Equivalent admittance network for complex permittivity elements, (a) Rectangular

The complex derivative of equation (52) presents the important decoupling of its real and imaginary parts, thus actually providing separate real and imaginary sensitivities. Specifically its real part constitutes the sensitivity of the measured complex voltage with respect to the conductivity *σj.* Respectively, the imaginary part comprises the sensitivity with respect to the jth element permittivity *εrj*. Additionally each one of them can be discriminated into two terms by substituting the complex *Vi* from (51). Consequently, the complex permittivity imaging is readily provided by applying the basic MPM algorithm (49) for each one of the four sensitivity terms. Specifically, considering the conductivity

*<sup>j</sup>* update we may set-up a summation like that of (49) comprised of the contributions of the

*Vi*

(49). The real and imaginary voltage contributions are summed together and multiplied by the proper relaxation factor to form the formula for σ-imaging. Likewise, the real and imaginary voltages are weighted with the corresponding sensitivities with respect to

0 0 1 1 0 22 ( ( ) ) *nn n nn n n n cj j r j rcj j <sup>j</sup> r j r j j jj k*

> *M M re re re im im im mi ci mi ci i i re im i i mi j j mi M M re im k k k k j j*

*VV VV V V V V*

*V V*

1 1

  *Vi*

 

> 

*Wj k <sup>W</sup>* (55)

*<sup>j</sup>* . But, now the imaginary part of

 

(56)

imaging

*<sup>j</sup>* comprises another sum similar to

*<sup>j</sup>* to formulate the ε-imaging algorithm. These

11 11

 

 

> 

and (b) its equivalent networks.

permittivity / ( ) *re Vi*

> 

algorithm as:

   

> 

**4.3 Complex MPM formulation for the MHz range** 

real voltage part *re Vi* weighted by the sensitivities / *re*

 *<sup>j</sup>* and / ( ) *im Vi*

> 

1

*W*

 

two could be given either separately or combined together in a complex *<sup>c</sup>*

1 1

 

> 

the voltage *im Vi* weighted by its sensitivity / *im*

*Extending MPM to complex quantities:* Recall that when Laplace equation is solved only for a conductivity distribution the calculated voltage will be real or in-phase with the injected sinusoidal current (source), just as injecting a current in its equivalent network in Fig.10. In contrary, when permittivity is included in the model a 90 phase shifted (quadrature) component appears in the developed voltage as well. This is in turn equivalent to injecting a sinusoidal current on a complex admittance network (Fig.11) and measuring the complex voltage across any ith admittance, which reads:

$$V\_i = V\_i^{re} + jV\_i^{im} = V\_i^{re}(\sigma\_\prime \mathcal{E}\_r) + jV\_i^{im}(\sigma\_\prime \mathcal{E}\_r) \tag{51}$$

In view of the above a complex sensitivity *ij i cj J V* is required in order to establish an algorithm similar to (49). For this purpose a complex derivative can be defined as [33-36]:

$$J\_{i\dot{j}} = \frac{\partial V\_i}{\partial (\alpha \varepsilon\_{\dot{\gamma}})} = \frac{\partial V\_i}{\partial \sigma\_{\dot{\gamma}}} + j \frac{\partial V\_i}{\partial (\alpha \varepsilon\_{\dot{\gamma}})} \tag{52}$$

with:

$$
\dot{j}\rho\mathbf{e}\_{c\dot{\gamma}} = \dot{j}\rho\mathbf{e}\_{\dot{\gamma}} + \boldsymbol{\sigma}\_{\dot{\gamma}}\,\text{or}\quad \dot{j}\rho\mathbf{e}\_{0}\boldsymbol{\varepsilon}\_{c\dot{\gamma}} = \dot{j}\rho\mathbf{e}\_{0}\boldsymbol{\varepsilon}\_{r\dot{\gamma}} + \boldsymbol{\sigma}\,\text{ or}\quad \rho\mathbf{e}\_{c\dot{\gamma}} = \rho\boldsymbol{\varepsilon}\_{\dot{\gamma}} + \dot{j}\boldsymbol{\sigma}\_{\dot{\gamma}}\tag{53}
$$

This complex derivative definition has already been employed by many researchers, e.g Franchois et al. [33] and Rekanos et al. [34-36], for the definition of a complex gradient (complex Jacobian matrix) working toward the establishment of a complex permittivity microwave imaging. The definition of (52) is exploited herein in order to identify the complex Jacobian matrix *ij J* into four real submatrices as presented by Polydorides [37]

$$J\_{ij} = \begin{bmatrix} J\_{ij}^{RR} & J\_{ij}^{RI} \\ J\_{ij}^{IR} & J\_{ij}^{II} \end{bmatrix} = \begin{bmatrix} \frac{\mathcal{G}V^{ne}}{\mathcal{G}\sigma\_j} & \frac{\mathcal{G}V^{ne}}{\mathcal{G}(\alpha\sigma\_j)} \\\\ \frac{\mathcal{G}V^{im}}{\mathcal{G}\sigma\_j} & \frac{\mathcal{G}V^{im}}{\mathcal{G}(\alpha\sigma\_j)} \end{bmatrix} \tag{54}$$

An alternative approach, instead of using a complex derivative is to expand the problem into two real functions ( , ) *re V r* and ( , ) *im V r* each one depended on two real unknown variables (, ) *r* . In turn the basic MPM algorithm (49) can be employed four times for each one of the sensitivities just as those occurring in (54). As explained above this is possible not only for two but even for any arbitrary number of variables. Attention must be paid to the normalized contributing terms (summations) similar to that of (49), since all terms referring to the same variable (, ) *r* must be added together after being weighted by an appropriate relaxation factor k. Hence, measurements at each ith port along with the solutions of the complex generalized Laplace equation yields complex voltages which are separated into real ( ) *re Vi* and imaginary parts ( ) *im Vi* . These numerical solutions yield complex voltage distributions at the FEM nodes all over the model of the body to be imaged. Before proceeding to the complex MPM algorithm, the analytical expressions for the evaluation of the above four sensitivities must be extracted.

*Extending MPM to complex quantities:* Recall that when Laplace equation is solved only for a conductivity distribution the calculated voltage will be real or in-phase with the injected sinusoidal current (source), just as injecting a current in its equivalent network in Fig.10. In contrary, when permittivity is included in the model a 90 phase shifted (quadrature) component appears in the developed voltage as well. This is in turn equivalent to injecting a sinusoidal current on a complex admittance network (Fig.11) and measuring the complex

> 

( ) () *ii i*

 

 or *cj jj* 

This complex derivative definition has already been employed by many researchers, e.g Franchois et al. [33] and Rekanos et al. [34-36], for the definition of a complex gradient (complex Jacobian matrix) working toward the establishment of a complex permittivity microwave imaging. The definition of (52) is exploited herein in order to identify the complex Jacobian matrix *ij J* into four real submatrices as presented by Polydorides [37]

> *ij ij jj ij IR II im im*

*JJ V V*

An alternative approach, instead of using a complex derivative is to expand the problem

one of the sensitivities just as those occurring in (54). As explained above this is possible not only for two but even for any arbitrary number of variables. Attention must be paid to the normalized contributing terms (summations) similar to that of (49), since all terms referring

relaxation factor k. Hence, measurements at each ith port along with the solutions of the complex generalized Laplace equation yields complex voltages which are separated into real ( ) *re Vi* and imaginary parts ( ) *im Vi* . These numerical solutions yield complex voltage distributions at the FEM nodes all over the model of the body to be imaged. Before proceeding to the complex MPM algorithm, the analytical expressions for the evaluation of

 

*cjj j*

algorithm similar to (49). For this purpose a complex derivative can be defined as [33-36]:

*VV V J j*

 or 0 0 *crj rj j j* 

*RR RI*

*J J*

*J*

 *r* 

 *r* 

the above four sensitivities must be extracted.

*ij ij*

 and ( , ) *im V r* 

 

(52)

 

 

( )

*re re*

 

 

*V V*

*j j*

 

. In turn the basic MPM algorithm (49) can be employed four times for each

 

( )

must be added together after being weighted by an appropriate

(51)

is required in order to establish an

 

*j* (53)

(54)

each one depended on two real unknown

voltage across any ith admittance, which reads:

with:

(, ) (, ) *re im re im V V jV V jV ii i i r i r*

*ij*

In view of the above a complex sensitivity *ij i cj J V*

*cj jj j j*

 

into two real functions ( , ) *re V*

to the same variable (, )

variables (, ) *r* 

Fig. 11. Equivalent admittance network for complex permittivity elements, (a) Rectangular and (b) its equivalent networks.

#### **4.3 Complex MPM formulation for the MHz range**

The complex derivative of equation (52) presents the important decoupling of its real and imaginary parts, thus actually providing separate real and imaginary sensitivities. Specifically its real part constitutes the sensitivity of the measured complex voltage with respect to the conductivity *σj.* Respectively, the imaginary part comprises the sensitivity with respect to the jth element permittivity *εrj*. Additionally each one of them can be discriminated into two terms by substituting the complex *Vi* from (51). Consequently, the complex permittivity imaging is readily provided by applying the basic MPM algorithm (49) for each one of the four sensitivity terms. Specifically, considering the conductivity *<sup>j</sup>* update we may set-up a summation like that of (49) comprised of the contributions of the real voltage part *re Vi* weighted by the sensitivities / *re Vi <sup>j</sup>* . But, now the imaginary part of the voltage *im Vi* weighted by its sensitivity / *im Vi <sup>j</sup>* comprises another sum similar to (49). The real and imaginary voltage contributions are summed together and multiplied by the proper relaxation factor to form the formula for σ-imaging. Likewise, the real and imaginary voltages are weighted with the corresponding sensitivities with respect to permittivity / ( ) *re Vi <sup>j</sup>* and / ( ) *im Vi <sup>j</sup>* to formulate the ε-imaging algorithm. These two could be given either separately or combined together in a complex *<sup>c</sup>* imaging algorithm as:

$$\log(\varepsilon\_{\underline{c}\underline{j}})^n = \sigma\_{\underline{j}}^n + jo\varepsilon\_0 \varepsilon\_{r\underline{j}}^n = jo\varepsilon\_0 (\varepsilon\_{\underline{r}\underline{j}})^n = \left[\sigma\_{\underline{j}}^{n-1} + k\_1 \mathcal{W}\_1 \sigma\_{\underline{j}}^{n-1}\right] + jo\varepsilon\_0 \left[\varepsilon\_{r\underline{j}}^{n-1} + k\_2 \mathcal{W}\_2 \varepsilon\_{r\underline{j}}^{n-1}\right] \tag{55}$$

$$\mathcal{W}\_1 = \frac{\sum\_{i=1}^{M} \frac{V\_{mi}^{re} - V\_{ci}^{re}}{V\_{mi}^{re}} \frac{\mathcal{G} V\_i^{re}}{\mathcal{G} \sigma\_j} + \frac{\sum\_{i=1}^{M} \frac{V\_{mi}^{im} - V\_{ci}^{im}}{V\_{mi}^{im}} \frac{\mathcal{G} V\_i^{im}}{\mathcal{G} \sigma\_j}}{\sum\_{k=1}^{M} \left| \frac{\mathcal{G} V\_k^{im}}{\mathcal{G} \sigma\_j} \right|} \tag{56}$$

High to Microwave Frequencies Imaging Techniques 179

where M is the total number of linear independent measurements, Emi and Eci are the measured and calculated fields at the ith antenna and k1,2 are the relaxations factors that may provide faster convergence. The optimum values of k1, k2 can be obtained through a numerical investigation. Observing the MPM reconstruction equations e.g. (59) to (60) or even the original static MPM in (49) one may discuss on the expected robustness or the immunity of the method to the problem of ill-posedness. As it is well understood the ill posedness is due to very large variations (of the order of 106) of the sensitivity. For example

transmitting antenna and along the axis of transmit-receive antenna pairs. But the sensitivity away from these object regions become very small and may be even lower than the measurement or calculations inaccuracies. When exact inverse problem formulations are employed these large variations in sensitivity yields an ill posed Hessian matrix accompanied by a lot of difficulties during its required inversion. Actually, the very large variation in sensitivity corresponds to very large variation in the Hessian matrix eigenvalues which reflects to a high singularity degree. In contrary, these very low sensitivities occur as very small weighting factors in the summations involved in MPM, e.g. (60). Thus, the very small sensitivity has negligible contribution in the complex permittivity update, which from a different point of view is shadowed by the higher sensitivities. This is equivalent to discarding very small singular values of the Jacobian (or Hessian), or using relatively large regularization parameter. Hence, the resulting reconstruction algorithm is expected to be robust during the first few iterations and it was indeed proved to converge always in our previous works, e.g. [26-30, 38]. The observed monotonous convergence was independent of the electrodes or antenna locations which greatly affects the sensitivity distribution. This is in turn a clear evidence of the MPM immunity against the problem of ill-posedness. Based on these observations MPM ensures convergence without the sophistication of regularization techniques. However, the penalty paid by MPM for its robustness is the absence of any control of the implicitly involved regularization. Instead, an exact method may involve a controllable regularization parameter which is gradually reduced from one iteration to the next. The reduction of the regularization parameter allows for the small

According to Meaney et al. [39], most regularization methods require the use of a varying degree of a priori information in order to ensure convergence. In contrary MPM ensures convergence without any a priory information, but its accuracy is compromised by its inability to exploit the hidden low sensitivities. Thus, a logical hypothesis toward accuracy improvement is to formulate an exact inverse problem (e.g., a Gauss-Newton scheme) and exploit the final image provided by MPM as an initial guess (a priori information). The regularization parameter within this scheme can be gradually set to zero or to a very small value, uncovering the very low sensitivities (exploiting the very small eigenvalues) which are expected to fine tune the image within a few iterations. The proof of this hypothesis

Observing the above MPM reconstruction expression it is obvious that the Jacobiansensitivity matrix constitutes the fundamental basis of the proposed method. Its accurate and computationally efficient calculation is of primary importance in ensuring successful imaging. For this purpose the next section elaborates on an Adjoint network for the microwave imaging and its counterpart network compensation theorem for the MHz range,

both resulting in closed form expressions for the Jacobian matrix entries.

may be observed for the voxels-k around the

very large sensitivity values *Ei ck*

eigenvalues or low sensitivities to be exploited.

constitutes one of our next tasks.

$$\mathcal{W}\_2 = \frac{\sum\_{i=1}^{M} \frac{V\_{mi}^{re} - V\_{ci}^{re}}{V\_{mi}^{re}} \frac{\mathcal{G} V\_i^{re}}{\mathcal{G}(\alpha \varepsilon\_0 \varepsilon\_{\bar{\tau}\bar{\gamma}})}}{\sum\_{k=1}^{M} \left| \frac{\mathcal{G} V\_k^{re}}{\mathcal{G}(\alpha \varepsilon\_0 \varepsilon\_{\bar{\tau}\bar{\gamma}})} \right|^2 + \frac{\sum\_{i=1}^{M} \frac{V\_{mi}^{im} - V\_{ci}^{im}}{V\_{mi}^{im}} \frac{\mathcal{G} V\_i^{im}}{\mathcal{G}(\alpha \varepsilon\_0 \varepsilon\_{\bar{\tau}\bar{\gamma}})}} \tag{57}$$

M is the total number of linearly independent measurements,*Vmi* and *Vci* are the measured and calculated voltage differences at the *ith* port (electrodes pair) and 1 2 *k k*, are the relaxations factors that may provide faster convergence.

The optimum values of 1 2 *k k*, can be obtained through a numerical investigation. If these are not known then is preferable to set 1 2 *k k* 1 (no relaxation). Even though a relaxation factor could in general speed up convergence, however in all performed numerical experimentation the best performance was obtained without relaxation.

#### **4.4 MPM for microwave imaging**

Working toward the extension of MPM algorithm to microwave imaging we simply need to substitute the complex voltage by the complex electric field. Indeed, the expressions of equations (55)-(57) are readily applicable for the 2-D microwave imaging when the electric field is comprised of only one component ˆ *E Ez <sup>z</sup>* or it is a complex scalar function. However, the proposed procedure is general and applies for the arbitrary 3-D imaging where the electric field is a complex vector quantity. Restricting ourselves to the 2-D case we may repeat the definition of the complex derivative just for clarity as:

$$J\_{ij} = \frac{\partial E\_i}{\partial \mathcal{E}\_{cj}} = \frac{\partial E\_i}{\partial \mathcal{E}\_{cj}^{real}} + j \frac{\partial E\_i}{\partial \mathcal{E}\_{cj}^{imag}} \tag{58}$$

In turn, by following exactly the same rationale as for the MHz imaging the microwave reconstruction algorithm reads:

$$\boldsymbol{\varepsilon}\_{ck}^{n} = \left[\boldsymbol{\varepsilon}\_{ck}^{\text{real}(n-1)} + k\_1 \boldsymbol{\mathcal{W}}\_1 \boldsymbol{\varepsilon}\_{ck}^{\text{real}(n-1)}\right] + j \left[\boldsymbol{\varepsilon}\_{ck}^{\text{imag}(n-1)} + k\_2 \boldsymbol{\mathcal{W}}\_2 \boldsymbol{\varepsilon}\_{ck}^{\text{imag}(n-1)}\right] \tag{59}$$

$$\mathcal{M}\_{1} = \frac{\sum\_{i=1}^{M} \frac{E\_{mi}^{real} - E\_{ci}^{real}}{E\_{mi}^{real}} \frac{\partial E\_{i}^{real}}{\partial \mathcal{E}\_{ck}^{real}}}{\sum\_{i=1}^{M} \left| \frac{\partial E\_{ck}^{real}}{\partial \mathcal{E}\_{ck}^{real}} \right|} + \frac{\sum\_{i=1}^{M} \frac{E\_{mi}^{imag} - E\_{ci}^{imag}}{E\_{mi}^{imag}} \frac{\partial E\_{i}^{imag}}{\partial \mathcal{E}\_{ck}^{real}}}{\sum\_{i=1}^{M} \left| \frac{\partial E\_{k}^{imag}}{\partial \mathcal{E}\_{ck}^{real}} \right|} \tag{60a}$$

$$\mathcal{W}\_2 = \frac{\sum\_{i=1}^{M} \frac{E\_{mi}^{real} - E\_{ci}^{real}}{\mathcal{C}\_{mi}^{real}} \frac{\partial \mathcal{E}\_i^{real}}{\partial \mathcal{E}\_{ck}^{imag}}}{\sum\_{i=1}^{M} \left| \frac{\partial \mathcal{E}\_k^{real}}{\partial \mathcal{E}\_{ck}^{imag}} \right|} + \frac{\sum\_{i=1}^{M} \frac{E\_{mi}^{imag} - E\_{ci}^{imag}}{\mathcal{E}\_{mi}^{imag}} \frac{\partial \mathcal{E}\_i^{imag}}{\partial \mathcal{E}\_{ck}^{imag}}}{\sum\_{i=1}^{M} \left| \frac{\partial \mathcal{E}\_k^{imag}}{\partial \mathcal{E}\_{ck}^{imag}} \right|} \tag{60b}$$

1 1 0 0

*V V*

 

 

> 

experimentation the best performance was obtained without relaxation.

may repeat the definition of the complex derivative just for clarity as:

 

1 1

1 1

*<sup>E</sup> <sup>E</sup> <sup>W</sup>*

*<sup>E</sup> <sup>E</sup> <sup>W</sup>*

*M M re re re im im im mi ci mi ci i i re im i i mi rj mi rj M M re im k k k k rj rj*

*V V V V V V*

() ()

 

or it is a complex scalar function.

(58)

 

  (60a)

(60b)

(57)

 

 

 

1 1 0 0

M is the total number of linearly independent measurements,*Vmi* and *Vci* are the measured and calculated voltage differences at the *ith* port (electrodes pair) and 1 2 *k k*, are the

The optimum values of 1 2 *k k*, can be obtained through a numerical investigation. If these are not known then is preferable to set 1 2 *k k* 1 (no relaxation). Even though a relaxation factor could in general speed up convergence, however in all performed numerical

Working toward the extension of MPM algorithm to microwave imaging we simply need to substitute the complex voltage by the complex electric field. Indeed, the expressions of equations (55)-(57) are readily applicable for the 2-D microwave imaging when the electric

However, the proposed procedure is general and applies for the arbitrary 3-D imaging where the electric field is a complex vector quantity. Restricting ourselves to the 2-D case we

> *ii i ij real imag cj cj cj EE E J j*

( 1) ( 1) ( 1) ( 1) 1 1 2 2

 

*M real real real M imag imag imag mi ci i mi ci i real real imag real i mi ck i mi ck M real M imag k k real real i ck i ck*

*E EE E EE*

*E E*

*M M real real real imag imag imag mi ci i mi ci i real imag imag imag i i mi ck mi ck M M real imag k k imag imag i i ck ck*

*E E*

*EEE E EE*

*kW j kW* (59)

 

In turn, by following exactly the same rationale as for the MHz imaging the microwave

*n real n real n imag n imag n ck ck ck ck ck*

1 1

1 1

 

 

 

() ()

*V V*

2

relaxations factors that may provide faster convergence.

field is comprised of only one component ˆ *E Ez <sup>z</sup>*

*W*

**4.4 MPM for microwave imaging** 

reconstruction algorithm reads:

1

2

where M is the total number of linear independent measurements, Emi and Eci are the measured and calculated fields at the ith antenna and k1,2 are the relaxations factors that may provide faster convergence. The optimum values of k1, k2 can be obtained through a numerical investigation. Observing the MPM reconstruction equations e.g. (59) to (60) or even the original static MPM in (49) one may discuss on the expected robustness or the immunity of the method to the problem of ill-posedness. As it is well understood the ill posedness is due to very large variations (of the order of 106) of the sensitivity. For example very large sensitivity values *Ei ck* may be observed for the voxels-k around the transmitting antenna and along the axis of transmit-receive antenna pairs. But the sensitivity away from these object regions become very small and may be even lower than the measurement or calculations inaccuracies. When exact inverse problem formulations are employed these large variations in sensitivity yields an ill posed Hessian matrix accompanied by a lot of difficulties during its required inversion. Actually, the very large variation in sensitivity corresponds to very large variation in the Hessian matrix eigenvalues which reflects to a high singularity degree. In contrary, these very low sensitivities occur as very small weighting factors in the summations involved in MPM, e.g. (60). Thus, the very small sensitivity has negligible contribution in the complex permittivity update, which from a different point of view is shadowed by the higher sensitivities. This is equivalent to discarding very small singular values of the Jacobian (or Hessian), or using relatively large regularization parameter. Hence, the resulting reconstruction algorithm is expected to be robust during the first few iterations and it was indeed proved to converge always in our previous works, e.g. [26-30, 38]. The observed monotonous convergence was independent of the electrodes or antenna locations which greatly affects the sensitivity distribution. This is in turn a clear evidence of the MPM immunity against the problem of ill-posedness. Based on these observations MPM ensures convergence without the sophistication of regularization techniques. However, the penalty paid by MPM for its robustness is the absence of any control of the implicitly involved regularization. Instead, an exact method may involve a controllable regularization parameter which is gradually reduced from one iteration to the next. The reduction of the regularization parameter allows for the small eigenvalues or low sensitivities to be exploited.

According to Meaney et al. [39], most regularization methods require the use of a varying degree of a priori information in order to ensure convergence. In contrary MPM ensures convergence without any a priory information, but its accuracy is compromised by its inability to exploit the hidden low sensitivities. Thus, a logical hypothesis toward accuracy improvement is to formulate an exact inverse problem (e.g., a Gauss-Newton scheme) and exploit the final image provided by MPM as an initial guess (a priori information). The regularization parameter within this scheme can be gradually set to zero or to a very small value, uncovering the very low sensitivities (exploiting the very small eigenvalues) which are expected to fine tune the image within a few iterations. The proof of this hypothesis constitutes one of our next tasks.

Observing the above MPM reconstruction expression it is obvious that the Jacobiansensitivity matrix constitutes the fundamental basis of the proposed method. Its accurate and computationally efficient calculation is of primary importance in ensuring successful imaging. For this purpose the next section elaborates on an Adjoint network for the microwave imaging and its counterpart network compensation theorem for the MHz range, both resulting in closed form expressions for the Jacobian matrix entries.

High to Microwave Frequencies Imaging Techniques 181

Assuming that the *l-*th branch belongs to the equivalent network of the j-th element, then:

( ) *Y YS j j jj* 

Since the voltage is linearly dependent on the network admittances then the superposition principle can be applied in two-fold. Summing over all element branches (Fig.11) one may obtain the total variation of the i-th port voltage *ΔVi*. Additionally, in order to evaluate the four Jacobian submatrices, the above approach can be applied once when only the

*re L*

*ij Rl j i <sup>V</sup> <sup>J</sup> S V I*

 

*im L*

*ij Il j i <sup>V</sup> <sup>J</sup> S V I*

 

> <sup>1</sup> ( ) *re L*

 

> <sup>1</sup> ( ) *im L*

*ij Rl j i <sup>V</sup> <sup>J</sup> S V I*

As a point of interest showing the potential of extending this method to microwave imaging, it is proved that the employed Network Compensation Theorem constitutes a particular case of the more general ``Adjoint Network Theorem'', which is in turn a direct consequence

Regarding the microwave sensitivity a more general approach is employed based on the electromagnetic field of a particular illuminating antenna combined with an adjoint configuration field (herein that for a different transmitting antenna location) through the Reciprocity Theorem. Explicitly the components of the Jacobian are the partial derivatives

below this is in turn evaluated through closed form expressions resulting from the

of the e-th element-pixel, when the s-th antenna is activated. As explained

 

*I*

*ij Il j i <sup>V</sup> <sup>J</sup> S V*

conductivity of the jth element is varied ( ) *Y S j j*

 

and imaginary parts one may end up to the desired Jacobian submatrices as:

*RR i*

*IR i*

*RI i*

*II i*

 

 

permittivity is varied ( *Y Sj j j*

of the electromagnetism ``Reciprocity Theorem''.

(or the sensitivities) of the electric field *Er*

**5.2 Sensitivity equation based on Adjoint Network Theorem** 

where

permittivity *<sup>e</sup>*

*c* 

 

(64)

). Applying (62) for each case and separating into real

1 <sup>1</sup>

1 <sup>1</sup>

1

1

*re re im im V VV V V Rl i jl kjl ijl kjl* (66a)

*re im im re V VV V V Il i jl kjl ijl kjl* (66b)

measured at the r-th antenna to the complex

and again when only the

(65a)

(65b)

(65c)

(65d)

#### **5. Jacobian or sensitivity matrix**

As explained in the introductory section the current paths following through the object depend on the unknown σ and εr profiles. Also the current density (or field intensity) is very high around the active electrodes (or antennas). These two phenomena make the σ- and εrimaging a strongly non-linear inverse problem by creating zones with very high and others with very low sensitivity. This severe singularity presents ratios of maximum to minimum singular values σmax/σmin as high as 106. Mathematically this reflects to a strongly non-linear Jacobian function *Jij(σ,εr)*. The present method handles the non-linearity, since it calculates the Jacobian matrix at each n-th iteration, namely for each *(σ,εr)* distribution. The additional computation cost of this task is negligible, since the Jacobian is calculated from closed form expressions from already available voltage or electric field distribution data.

Let us proceed to the presentation of the methodologies resulting the closed form expressions for the Jacobian matrix.

#### **5.1 Compensation theorem on an admittance network**

The sensitivities of the MHz range imaging can be obtained by extending the resistive network compensation theorem (50) originally given by Yorkey et al. [31] and later on reviewed in detail by Sahalos et al. [32](p 229). This is equivalent to the application of the circuits compensation theorem on an admittance network. For this purpose an equivalent admittance network is considered for each triangular or rectangular element rj (,) *j* as shown in Fig.11. The admittance values are obtained from the non-diagonal entries of the corresponding element matrix as: ( ) *j j e e Y jK* where as defined *<sup>e</sup> S* nondiagonal entries of *<sup>e</sup> K* .

For a triangular element (Fig.11a-b)

$$\begin{bmatrix} Y \end{bmatrix}\_3 = \begin{bmatrix} Y\_{11} & -Y\_{12} & -Y\_{13} \\ & Y\_{22} & -Y\_{23} \\ \text{symmetric} & & Y\_{33} \end{bmatrix} = (\sigma\_j + j\alpha\varepsilon\_j) \cdot \begin{bmatrix} \mathbf{K} \end{bmatrix}\_3 \tag{61}$$

For a rectangular element (Fig.11c-d):

$$\begin{bmatrix} Y \end{bmatrix}\_{4} = \begin{bmatrix} Y\_{11} & -Y\_{12} & -Y\_{13} & -Y\_{14} \\ & Y\_{22} & -Y\_{23} & -Y\_{24} \\ & & Y\_{33} & -Y\_{34} \\ \text{symmetric} & & & Y\_{44} \end{bmatrix} = \begin{pmatrix} \sigma\_{j} + jo\sigma\_{j} \end{pmatrix} \cdot \begin{bmatrix} K \end{bmatrix}\_{4} \tag{62}$$

Where <sup>3</sup> *K* and <sup>4</sup> *K* are defined in (50) and *Y G jC ij ij ij*

According to the network compensation theorem, a differential change of the *l-*th branch admittance *Yl* by *ΔYl* causes a differential variation *ΔVi* at the *i-*th port voltage as:

$$
\Delta V\_{i\ell} = -\frac{1}{I\_i} \Delta Y\_{\ell} \cdot V\_{ij\ell} \cdot V\_{kj\ell} \tag{63}
$$

Assuming that the *l-*th branch belongs to the equivalent network of the j-th element, then:

$$
\Delta Y\_{\ell} \to \Delta Y\_{\ell j} = S\_{\ell} \left( \Delta \sigma\_{j} + j a \Delta \varepsilon\_{j} \right) \tag{64}
$$

Since the voltage is linearly dependent on the network admittances then the superposition principle can be applied in two-fold. Summing over all element branches (Fig.11) one may obtain the total variation of the i-th port voltage *ΔVi*. Additionally, in order to evaluate the four Jacobian submatrices, the above approach can be applied once when only the conductivity of the jth element is varied ( ) *Y S j j* and again when only the permittivity is varied ( *Y Sj j j* ). Applying (62) for each case and separating into real and imaginary parts one may end up to the desired Jacobian submatrices as:

$$J\_{ij}^{RR} = \frac{\mathcal{G}V\_i^{re}}{\mathcal{G}\sigma\_j} = -\frac{1}{I\_i} \sum\_{\ell=1}^{L} S\_{\ell} \cdot V\_{R\ell} \tag{65a}$$

$$J\_{i\dot{j}}^{IR} = \frac{\mathcal{B}V\_i^{im}}{\mathcal{B}\sigma\_j} = -\frac{1}{I\_i} \sum\_{\ell=1}^L \mathbf{S}\_{\ell} \cdot V\_{\text{ll}} \tag{65b}$$

$$J\_{i\dot{j}}^{RI} = \frac{\mathcal{G}V\_i^{re}}{\mathcal{G}(occ\_j)} = +\frac{1}{I\_i} \sum\_{\ell=1}^{L} S\_\ell \cdot V\_{\ell l} \tag{65c}$$

$$J\_{ij}^{II} = \frac{\mathcal{G}V\_i^{im}}{\mathcal{G}(\rho \alpha \varepsilon\_j)} = -\frac{1}{I\_i} \sum\_{\ell=1}^{L} \mathbf{S}\_{\ell} \cdot V\_{R\ell} \tag{65d}$$

where

180 Medical Imaging

As explained in the introductory section the current paths following through the object depend on the unknown σ and εr profiles. Also the current density (or field intensity) is very high around the active electrodes (or antennas). These two phenomena make the σ- and εrimaging a strongly non-linear inverse problem by creating zones with very high and others with very low sensitivity. This severe singularity presents ratios of maximum to minimum singular values σmax/σmin as high as 106. Mathematically this reflects to a strongly non-linear Jacobian function *Jij(σ,εr)*. The present method handles the non-linearity, since it calculates the Jacobian matrix at each n-th iteration, namely for each *(σ,εr)* distribution. The additional computation cost of this task is negligible, since the Jacobian is calculated from closed form

Let us proceed to the presentation of the methodologies resulting the closed form

The sensitivities of the MHz range imaging can be obtained by extending the resistive network compensation theorem (50) originally given by Yorkey et al. [31] and later on reviewed in detail by Sahalos et al. [32](p 229). This is equivalent to the application of the circuits compensation theorem on an admittance network. For this purpose an equivalent admittance network is considered for each triangular or rectangular element rj (,)

shown in Fig.11. The admittance values are obtained from the non-diagonal entries of the

 

 11 12 13

According to the network compensation theorem, a differential change of the *l-*th branch

1 *i ij kj i V YV V I*

*Y YY <sup>Y</sup> j K Y Y*

22 23 24 4 4 33 34 44

*Y YY j K*

*Y YY*

11 12 13 14

*Y YYY*

*symmetric Y*

admittance *Yl* by *ΔYl* causes a differential variation *ΔVi* at the *i-*th port voltage as:

*symmetric Y*

3 3 22 23 33

( ) *j j*

( ) *j j*

(63)

 

*ij*

   *j* as

(61)

(62)

where as defined *<sup>e</sup> S* non-

expressions from already available voltage or electric field distribution data.

**5.1 Compensation theorem on an admittance network** 

corresponding element matrix as: ( ) *j j e e Y jK*

Where <sup>3</sup> *K* and <sup>4</sup> *K* are defined in (50) and *Y G jC ij ij*

**5. Jacobian or sensitivity matrix** 

expressions for the Jacobian matrix.

diagonal entries of *<sup>e</sup> K* .

For a triangular element (Fig.11a-b)

For a rectangular element (Fig.11c-d):

$$V\_{Rl} = V\_{ijl}^{rc} V\_{kjl}^{rc} - V\_{ijl}^{im} V\_{kjl}^{im} \tag{66a}$$

$$V\_{ll} = V\_{ijl}^{re} V\_{kjl}^{im} + V\_{ijl}^{im} V\_{kjl}^{re} \tag{66b}$$

As a point of interest showing the potential of extending this method to microwave imaging, it is proved that the employed Network Compensation Theorem constitutes a particular case of the more general ``Adjoint Network Theorem'', which is in turn a direct consequence of the electromagnetism ``Reciprocity Theorem''.

#### **5.2 Sensitivity equation based on Adjoint Network Theorem**

Regarding the microwave sensitivity a more general approach is employed based on the electromagnetic field of a particular illuminating antenna combined with an adjoint configuration field (herein that for a different transmitting antenna location) through the Reciprocity Theorem. Explicitly the components of the Jacobian are the partial derivatives (or the sensitivities) of the electric field *Er* measured at the r-th antenna to the complex permittivity *<sup>e</sup> c* of the e-th element-pixel, when the s-th antenna is activated. As explained below this is in turn evaluated through closed form expressions resulting from the

High to Microwave Frequencies Imaging Techniques 183

*aa a*

 

 

*ck ck ck ck ck H H HE H E E <sup>E</sup> <sup>j</sup>*

Adding (74) and (70) the similar terms are canceled, then the use of vector identities

 (75)

Taking the volume integral over a domain enclosed within a sphere of an infinite radius and

The left hand side is identically zero, since both fields intensities tend to zero at infinity by simply considering artificially small unbounded medium losses. Hence, (76) finally gives the

*r k*

The sensitivity equation (77) can be further simplified to yield a closed form for the Jacobian

specific receiving- sensing antennas. Starting from the integrals of (77), these are in general over the whole solution domain. But, actually the left hand side is restricted over the current carrying volume (Vr) of the r-th antenna. In turn, considering the definition (36) for the complex dielectric expansion the right hand side integral of (77) is obviously restricted over the volume (Vk) of the k-th element. Further a closer look at the dot product of (77) left hand side reveals that the involved electric field is that produced by a radiating antenna at the *<sup>s</sup> r r* position (Fig. 12) illuminating the receiving antenna at *<sup>r</sup> r r* position (sensing port)

*a a <sup>a</sup>*

(76)

*<sup>a</sup>*

*ck* by first introducing the FEM basis functions and by considering the

*ck* can be identified as the sensitivity of the r-th receiving antenna (r-th

*r E EdV*

(74)

 

 

(73)

*a a a a <sup>a</sup> <sup>a</sup>*

*aaa <sup>a</sup> <sup>a</sup>*

*S V ck ck ck ck E H <sup>E</sup> <sup>E</sup> H dS <sup>j</sup>*

Likewise, dot-multiplying *<sup>a</sup> <sup>E</sup>* (68b) and subtracting *<sup>H</sup>*

results to:

*S V ck ck ck*

*V V ck <sup>E</sup> J dV <sup>j</sup>*

sensing port) with respect to the k-th element complex permittivity

with only the field component parallel to its current *rJ*

*ck ck ck ck ck H E E HE E H H E <sup>J</sup> <sup>j</sup>*

*EHE H E dS <sup>J</sup> <sup>j</sup>*

 

 

 

applying the divergence theorem ends up to:

(72)

*c r*

*rE E j*

*H*

*r k*

producing a net effect. Hence, the

*ck* when the s-th

 

*rE E*

*ck* (69a) yields:

*H j E J dV*

*r k*

(77)

*r E E dV*

*c k*

*E j*

 *V S F du F dS*

to get:

and *AB BA*

so called sensitivity equation as:

matrix entries *E*

derivative *E*

*AB BA*

reciprocity theorem and the employment of an adjoint problem. For this purpose an approach similar to that given by Oldenburg [40] and the original research cited therein is adopted. Namely, the two Maxwell Curl equations are written for the source *Js* at s-th antenna and the involved fields *E H*, are differentiated with respect to the e-th complex permittivity. For the adjoint fields , *a a E H* these two curl equations are written considering a source *Jr* at the r-th antenna (Fig.12). The four curl equations are in turn combined following the reciprocity theorem procedure. Specifically, for a dipole antenna at the s-th location the curl equations reads:

$$
\vec{\nabla} \times \vec{E} = -j\alpha\mu \vec{H} \tag{67a}
$$

$$
\vec{\nabla} \times \vec{H} = j a \kappa\_c \vec{E} + \vec{f}\_s \tag{67b}
$$

Taking the derivatives of (67a) and (67b) with respect to *ck* defined in Eq. (36):

$$
\vec{\nabla} \times \frac{\partial \vec{E}}{\partial \varepsilon\_{ck}} = -j\alpha\mu \frac{\partial \vec{H}}{\partial \varepsilon\_{ck}}\tag{68a}
$$

$$\vec{\nabla} \times \frac{\partial \vec{H}}{\partial \varepsilon\_{ck}} = jao \varepsilon\_c \frac{\partial \vec{E}}{\partial \varepsilon\_{ck}} + joob \nu\_k(\vec{r}) \vec{E} \tag{68b}$$

Consider an auxiliary (adjoint) Maxwell problem with a dipole source *Jr* located at the observation point *<sup>r</sup> r r* as shown in Fig.12.

$$
\vec{\nabla} \times \vec{E}^a = -j a \rho \mu \vec{H}^a \tag{69a}
$$

$$
\vec{\nabla} \times \vec{H}^a = j a \sigma\_c \vec{E}^a + \vec{f}\_r \tag{69b}
$$

For the application of the Reciprocity theorem procedure to (68a), (68b) and (69a), (69b) according to Balanis [24](p. 324), lets take the inner product of *<sup>a</sup> <sup>H</sup>* (68a) and *ck E* (69b) and subtract:

$$
\vec{H}^a \cdot \vec{\nabla} \times \frac{\partial \vec{E}}{\partial \varepsilon\_{ck}} - \frac{\partial \vec{E}}{\partial \varepsilon\_{ck}} \cdot \vec{\nabla} \times \vec{H}^a = -j a \rho \mu \vec{H}^a \cdot \frac{\partial \vec{H}}{\partial \varepsilon\_{ck}} - j a \varepsilon\_c \frac{\partial \vec{E}}{\partial \varepsilon\_{ck}} \cdot \vec{E}^a - \frac{\partial \vec{E}}{\partial \varepsilon\_{ck}} \cdot \vec{J}\_r \tag{70}
$$

Making use of the identity:

$$
\vec{\nabla} \cdot \left( \vec{A} \times \vec{B} \right) = \vec{B} \cdot \left( \vec{\nabla} \times \vec{A} \right) - \vec{A} \cdot \left( \vec{\nabla} \times \vec{B} \right) \tag{71}
$$

$$\text{For } \vec{B} = \vec{H}^a \text{ and } \vec{A} = \frac{\partial \vec{E}}{\partial \varepsilon\_{ck}} \text{ the left hand side of (70) can be written as } \vec{\nabla} \cdot \left( \frac{\partial \vec{E}}{\partial \varepsilon\_{ck}} \times \vec{H}^a \right).$$

In turn, one may take a volume integral of (70) over a domain enclosed by a sphere with a radius tending to infinity and apply the divergence theorem on the left side.

$$\iiint\_V \left(\vec{\nabla} \cdot \vec{F}\right) d\mu = \sflectbox{\bfthinspace}\_S \vec{F} \cdot d\vec{S} \tag{72}$$

to get:

182 Medical Imaging

reciprocity theorem and the employment of an adjoint problem. For this purpose an approach similar to that given by Oldenburg [40] and the original research cited therein is adopted. Namely, the two Maxwell Curl equations are written for the source *Js* at s-th antenna and the involved fields *E H*, are differentiated with respect to the e-th complex

a source *Jr* at the r-th antenna (Fig.12). The four curl equations are in turn combined following the reciprocity theorem procedure. Specifically, for a dipole antenna at the s-th

> *E jH*

 *H j EJ c s*

 

 

*ck ck*

*ck ck E H <sup>j</sup>*

*H E j j rE*

 

Consider an auxiliary (adjoint) Maxwell problem with a dipole source *Jr* located at the

*a a E jH* 

*a a H c r j*

For the application of the Reciprocity theorem procedure to (68a), (68b) and (69a), (69b)

*ck ck ck ck ck*

the left hand side of (70) can be written as *<sup>a</sup>*

In turn, one may take a volume integral of (70) over a domain enclosed by a sphere with a

radius tending to infinity and apply the divergence theorem on the left side.

(70)

 *H j E J*

*AB B A A B* (71)

according to Balanis [24](p. 324), lets take the inner product of *<sup>a</sup> <sup>H</sup>* (68a) and

*a a a a*

*E E H EE H H <sup>j</sup>*

 

*c k*

  these two curl equations are written considering

(67a)

(67b)

(69a)

*ck E* 

 

(69b)

defined in Eq. (36):

(68a)

*E J* (69b)

*c r*

*ck <sup>E</sup> <sup>H</sup>* 

 .

(68b)

permittivity. For the adjoint fields , *a a E H*

observation point *<sup>r</sup> r r* as shown in Fig.12.

*ck*

*<sup>E</sup> <sup>A</sup>* 

 

and subtract:

For *<sup>a</sup> B H*

and

Making use of the identity:

Taking the derivatives of (67a) and (67b) with respect to *ck*

location the curl equations reads:

$$\iiint\_{S} \left( \frac{\partial \vec{E}}{\partial \varepsilon\_{ck}} \times \vec{H}^{a} \right) \cdot d\vec{S} = \iiint\_{V} \left( -j a \mu \vec{H}^{a} \cdot \frac{\partial \vec{H}}{\partial \varepsilon\_{ck}} - j a \varepsilon\_{c} \frac{\partial \vec{E}}{\partial \varepsilon\_{ck}} \cdot \vec{E}^{a} - \frac{\partial \vec{E}}{\partial \varepsilon\_{ck}} \cdot \vec{J}\_{r} \right) dV \tag{73}$$

Likewise, dot-multiplying *<sup>a</sup> <sup>E</sup>* (68b) and subtracting *<sup>H</sup> ck* (69a) yields:

$$\vec{\nabla} \cdot \left(\frac{\partial \vec{H}}{\partial \varepsilon\_{ck}} \times \vec{E}^a \right) = \vec{E}^a \cdot \vec{\nabla} \times \frac{\partial \vec{H}}{\partial \varepsilon\_{ck}} - \frac{\partial \vec{H}}{\partial \varepsilon\_{ck}} \cdot \vec{\nabla} \times \vec{E}^a = jao\varepsilon\_c \vec{E}^a \cdot \frac{\partial \vec{E}}{\partial \varepsilon\_{ck}} + jao\nu\_k \left(\vec{r} \right) \vec{E}^a \cdot \vec{E} + jo\mu \mu \frac{\partial \vec{H}}{\partial \varepsilon\_{ck}} \cdot \vec{H}^a \tag{74}$$

Adding (74) and (70) the similar terms are canceled, then the use of vector identities *AB BA* and *AB BA* results to:

$$
\vec{\nabla} \cdot \left( \frac{\partial \vec{H}}{\partial \varepsilon\_{ck}} \times \vec{E}^a \right) + \vec{\nabla} \cdot \left( \frac{\partial \vec{E}}{\partial \varepsilon\_{ck}} \times \vec{H}^a \right) = \vec{\nabla} \cdot \left( \frac{\partial \vec{E}}{\partial \varepsilon\_{ck}} \times \vec{H}^a - \vec{E}^a \times \frac{\partial \vec{H}}{\partial \varepsilon\_{ck}} \right) = -\frac{\partial \vec{E}}{\partial \varepsilon\_{ck}} \cdot \vec{J}\_r + j\alpha \vec{\nu}\_k \left( \vec{r} \right) \vec{E}^a \cdot \vec{E} \tag{75}
$$

Taking the volume integral over a domain enclosed within a sphere of an infinite radius and applying the divergence theorem ends up to:

$$\iiint\_{S} \left( \frac{\partial \vec{E}}{\partial \mathcal{E}\_{ck}} \times \vec{H}^{a} - \vec{E}^{a} \times \frac{\partial \vec{H}}{\partial \mathcal{E}\_{ck}} \right) d\vec{S} = \iiint\_{V} \left( -\frac{\partial \vec{E}}{\partial \mathcal{E}\_{ck}} \cdot \vec{f}\_{r} + j o \nu \vec{\nu}\_{k} \left( \vec{r} \right) \vec{E}^{a} \cdot \vec{E} \right) dV \tag{76}$$

The left hand side is identically zero, since both fields intensities tend to zero at infinity by simply considering artificially small unbounded medium losses. Hence, (76) finally gives the so called sensitivity equation as:

$$\iiint\_{V} \frac{\partial \vec{E}}{\partial \mathcal{E}\_{ck}} \cdot \vec{f}\_{r}dV = \iiint\_{V} jao \vec{\nu}\_{k} \left(\vec{r}\right) \vec{E}^{a} \cdot \vec{E}dV\tag{77}$$

The sensitivity equation (77) can be further simplified to yield a closed form for the Jacobian matrix entries *E ck* by first introducing the FEM basis functions and by considering the specific receiving- sensing antennas. Starting from the integrals of (77), these are in general over the whole solution domain. But, actually the left hand side is restricted over the current carrying volume (Vr) of the r-th antenna. In turn, considering the definition (36) for the complex dielectric expansion the right hand side integral of (77) is obviously restricted over the volume (Vk) of the k-th element. Further a closer look at the dot product of (77) left hand side reveals that the involved electric field is that produced by a radiating antenna at the *<sup>s</sup> r r* position (Fig. 12) illuminating the receiving antenna at *<sup>r</sup> r r* position (sensing port) with only the field component parallel to its current *rJ* producing a net effect. Hence, the derivative *E ck* can be identified as the sensitivity of the r-th receiving antenna (r-th sensing port) with respect to the k-th element complex permittivity *ck* when the s-th

High to Microwave Frequencies Imaging Techniques 185

*k*

*sr k J EFE* 

> *k kk ij i j k*

Matrices (vectors) *[Ek], [Eak]* represent the tangential electric fields along the edges of the k-th tetrahedral element in three dimensional case or the *Ez* electric field on the nodes of the k-th triangular element for the two dimensional case. These fields are already computed during the multiple forward problem solutions performed during the setup of the calculated dataset (once for each illuminating antenna). The matrix *<sup>k</sup> Fij* is a 6 ×6 matrix in 3D or a 3 ×3 matrix in 2D and its entries can be constructed from the FEM element matrices, [15].Recall at this point that FEM is applied on a fine mesh of tetrahedral or triangular elements while the image reconstruction is carried out on a coarse rectangular (pixels) or cubical (voxels) mesh. Each reconstruction element consists of a number of forward elements and similarly a number of nodes or edges that belong to these forward elements. Equation (83) results to a partial Jacobian matrix for each triangular or tetrahedral forward element. However, the desired Jacobian is that of the rectangular-pixel or cubical-voxel reconstruction element. For each evaluation the Fk matrices are assembled together according to a classical FEM procedure to yield a m×m matrix *Me*, where *e* the global number of the reconstruction element and m the number of edges or nodes that are inside the e-th element. Note that this matrix depends only on the geometry (independent of *εrc* distribution) and its calculated only once and stored for multiple usage. Consequently, the Jacobian matrix of the e-th

*k*

*V*

where *γ=jω/I* or *jω/IΔl* is the constant term for the line source and elementary dipole

*e e ae*

(85)

*sr e J EME* 

The above z-oriented current excitations yields a primarily z-polarized electric field which is expected to interact mainly with the *εzz* component of a possibly anisotropic complex permittivity. In general, the orientation of the radiating dipoles could be exploited as an additional degree of freedom to extract information from dielectrically anisotropic

*I l* 

*sr k i i j j k ck V i j E r <sup>j</sup> <sup>J</sup> E N E N dV*

 

*k k ak*

*k k ak k*

(83)

*F N N dV* (84)

(82)

*x y* . Similarly, the sensitivity equation becomes:

The above sensitivity expressions can be written in matrix form as:

, ,

, ,

where 1 2 , ,..., *e ee e E EE Em* and 1 2 , ,..., *ae ae ae ae E EE Em* .

, ,

where 2 2 

3-D:

structures, e.g. [41, 42].

excitations and

element reads:

antenna illuminates the structure, or specifically the (*s, r, k*) entry of the Jacobian matrix, *J E srk* , , *ck* .

Fig. 12. Geometry for the application of the reciprocity theorem

Since the integration is restricted over the k-th element volume *(Vk)*, then the fields *E* and *a E* can be expanded into the FEM basis (shape) functions according to (45). In view of the above description (77) can be rewritten in the form:

$$\iiint\_{V\_k} \frac{\partial \vec{E}}{\partial \varepsilon\_{ck}} \cdot \vec{J}\_r dV\_r = j o \iiint\_{V\_k} \vec{E} \cdot \vec{E} dV\_k = j o \iiint\_{V\_k} \sum\_i E\_i^k \cdot N\_i^k \sum\_j E\_j^{ak} \cdot N\_j^k dV\_k \tag{78}$$

Working toward a closed form expression for the Jacobian matrix, the next step is to consider specific antenna types. For the two-dimensional (2-D) case when the complex permittivity is assumed homogeneous in the z-direction, along which the structure is assumed infinite, it is convenient to employ infinitely extending thin line sources as illuminating antennas with current density defined by:

$$\text{I Line source:} \qquad \vec{f}\_r = I\delta\left(\vec{r} - \vec{r}\_r\right)\hat{z} \text{ where } I \text{ is a constant current} \tag{79}$$

In this case the sensitivity entries are readily simplified as:

$$\text{2-D: } J\_{\{(s,r),k\}} = \frac{\partial \vec{E}(\vec{r})}{\partial \mathcal{E}\_{ck}} = \frac{j\alpha}{I} \oint\_{S\_k} \sum\_i E\_i^k \cdot \mathbf{N}\_i^k \sum\_j E\_j^{ak} \cdot \mathbf{N}\_j^k dS\_k \tag{80}$$

For a three dimensional structure (3-D) illuminated by thin elementary dipoles of length *(Δl)*  the current density reads:

$$\text{z-oriented elementary dipole: } \vec{J}\_r = I\delta\left(\rho - \rho^\dagger\right)\hat{z} \tag{81}$$

antenna illuminates the structure, or specifically the (*s, r, k*) entry of the Jacobian matrix,

Fig. 12. Geometry for the application of the reciprocity theorem

*k kk*

*V VV ck i j*

above description (77) can be rewritten in the form:

illuminating antennas with current density defined by:

In this case the sensitivity entries are readily simplified as:

, ,

Line source: ˆ *r r J I r rz*

2-D:

the current density reads:

Since the integration is restricted over the k-th element volume *(Vk)*, then the fields *E*

can be expanded into the FEM basis (shape) functions according to (45). In view of the

*a k k ak k*

(78)

where *I* is a constant current (79)

*k k ak k*

 ˆ *rJI z* 

(80)

(81)

'

*r r k i i j j k*

*<sup>E</sup> J dV j E EdV j E N E N dV*

Working toward a closed form expression for the Jacobian matrix, the next step is to consider specific antenna types. For the two-dimensional (2-D) case when the complex permittivity is assumed homogeneous in the z-direction, along which the structure is assumed infinite, it is convenient to employ infinitely extending thin line sources as

*k*

For a three dimensional structure (3-D) illuminated by thin elementary dipoles of length *(Δl)* 

z-oriented elementary dipole:

*sr k i i j j k ck S i j E r <sup>j</sup> <sup>J</sup> E N E N dS I* 

   and

*J E srk* , ,

*a E*

*ck* . where 2 2 *x y* . Similarly, the sensitivity equation becomes:

$$\text{3-D: } J\_{\{(s,r),k\}} = \frac{\partial \vec{E}(\vec{r})}{\partial \varepsilon\_{ck}} = \frac{jao}{I\Delta l} \iiint\_{V\_k} \left(\sum\_i E\_i^k \cdot N\_i^k\right) \cdot \left(\sum\_j E\_j^{ak} \cdot N\_j^k\right) dV\_k \tag{82}$$

The above z-oriented current excitations yields a primarily z-polarized electric field which is expected to interact mainly with the *εzz* component of a possibly anisotropic complex permittivity. In general, the orientation of the radiating dipoles could be exploited as an additional degree of freedom to extract information from dielectrically anisotropic structures, e.g. [41, 42].

The above sensitivity expressions can be written in matrix form as:

$$J\_{\{(s,r),k\}} = \mathcal{V}\left[E^k\right] \cdot \left[\boldsymbol{F}^k\right] \cdot \left[\boldsymbol{E}^{ak}\right] \tag{83}$$

where *γ=jω/I* or *jω/IΔl* is the constant term for the line source and elementary dipole excitations and

$$\mathbb{E}\left[F\_{ij}^{k}\right] = \iiint\_{V\_k} \left\{N\_i^k\right\} \cdot \left\{N\_j^k\right\} dV\_k \tag{84}$$

Matrices (vectors) *[Ek], [Eak]* represent the tangential electric fields along the edges of the k-th tetrahedral element in three dimensional case or the *Ez* electric field on the nodes of the k-th triangular element for the two dimensional case. These fields are already computed during the multiple forward problem solutions performed during the setup of the calculated dataset (once for each illuminating antenna). The matrix *<sup>k</sup> Fij* is a 6 ×6 matrix in 3D or a 3 ×3 matrix in 2D and its entries can be constructed from the FEM element matrices, [15].Recall at this point that FEM is applied on a fine mesh of tetrahedral or triangular elements while the image reconstruction is carried out on a coarse rectangular (pixels) or cubical (voxels) mesh. Each reconstruction element consists of a number of forward elements and similarly a number of nodes or edges that belong to these forward elements. Equation (83) results to a partial Jacobian matrix for each triangular or tetrahedral forward element. However, the desired Jacobian is that of the rectangular-pixel or cubical-voxel reconstruction element. For each evaluation the Fk matrices are assembled together according to a classical FEM procedure to yield a m×m matrix *Me*, where *e* the global number of the reconstruction element and m the number of edges or nodes that are inside the e-th element. Note that this matrix depends only on the geometry (independent of *εrc* distribution) and its calculated only once and stored for multiple usage. Consequently, the Jacobian matrix of the e-th element reads:

$$J\_{\left(\left(s,r\right),\varepsilon\right)} = \mathcal{V}\left[E^{\varepsilon}\right] \cdot \left[M^{\varepsilon}\right] \cdot \left[E^{ae}\right] \tag{85}$$

where 1 2 , ,..., *e ee e E EE Em* and 1 2 , ,..., *ae ae ae ae E EE Em* .

High to Microwave Frequencies Imaging Techniques 187

For comparison purposes it is more convenient to present its normalized value *(SSQN)* with respect to its maximum *(SSQmax)* occurring at the first iteration of the inverse problem as:

While SSQ is general and can be calculated in all cases, it is only an indirect indication of

safeguarded. But, this condition may be disturbed by the problem singularity degree, which in turn depends on the data collection strategy. For this purpose the well- posedeness or the singularity degree, of the sensitivity-Jacobian matrix should be preliminary examined. However a further elaboration is required. Besides this for the particular case of the "computer test" the target or true *(σ,εr)* distributions are available. Hence, the estimated

1 1

*i i*

*P P ii i r r true r calc r true r true i i*

1 1

 *error* 

elements of the reconstruction mesh. We should keep in mind that the σ-error, εr-error criteria are not applicable for practical in-vivo or even for the "laboratory test" cases where

The algorithm was applied for the circular cross-section of Fig.13a and for the rectangular cross-section of Fig.13b. A total number of 32 electrodes were used, where only two of them are active in each projection angle. An indicative example with a double anomaly is presented herein. The target σ- and εr-profiles are shown in Fig.14a. The values for conductivity are

14 / *mS cm* and for permittivity r1

current was *f* 8*MHz* . The relaxation factors 1 2 *k k*, are set equal to unity. The reconstructed

The algorithm was also tested for the rectangular model of Fig.6. An anomaly with

now *f* 9*MHz* . The relaxation factors are set equal to unity. The target model and the

For details on the algorithm convergence rate and its performance, please wait for a follow

*<sup>o</sup>* 7 / *mS cm* and ro

 

 

*<sup>r</sup>* distributions are unknown.

20 / *mS cm* and permittivity r1

*<sup>o</sup>* 7 / *mS cm* and ro

reconstructed image after 15 iterations are shown in Fig. 15.

*P P ii i true calc true true*

 *error* 

*N*

convergence. Namely, its minimization ensure *<sup>c</sup>*

 

normalized deviations of , *n n*

**6.2 Quasi-static MHz reconstruction** 

image after 15 iterations are shown in Fig.14b

up publication which is now in preparation.

 *true* and *r true* 

> ,

21 / *mS cm* and 2

conductivity 1

background of

homogeneous background of

where

1 

the objects

defined as a norm of relative error:

max

2 2

2 2

 

 

are the average values of the target profiles and P is the number of

(89)

*SSQ SSQ SSQ* (88)

*calc calc* from their true values at the n-th iteration can be

(90)

150 . The frequency of the injected current is

300 and r2

300 was introduced in a homogeneous

100 . The frequency of the injected

150 , in a

convergence when a unique solution is

#### **5.3 Sensitivities and features of the Jacobian matrix**

The Jacobian matrix is a rectangular one with dimensions M>N where M the number of linearly independent measuring points and N the number of unknown *(σi,εri)* pairs. The study of tis characteristics is quite important especially in the course toward inventing the best data collection strategy, i.e. active and sensing antenna or electrode pairs locations and their subsequent activation. Additionally, the sensitivity matrix study with the aid of the recently revisited Proper Orthogonal Decomposition [43] techniques constitutes an interesting research challenge. Even though our group is already working on this subject it will not be considered herein as it requires a separate chapter of its own. For some preliminary results regarding the SVD analysis of the Jacobian matrix please contact our previous work [30] and wait for al follow up publication.

#### **6. Numerical results & discussions**

The proposed algorithm was applied for the imaging of numerous conductivity and permittivity distributions for both circular and rectangular models and satisfactory reconstructions are obtained. The background and anomaly conductivity and permittivity values resemble these of typical human tissues.

Before proceeding to the presentation of reconstructed images the convergence criteria should be first defined.

#### **6.1 Inverse problem convergence criteria**

*Computer Test Approach:* The so called "computer test" was employed in all cases throughout this chapter. First a "target model" is considered and the forward problem is solved for each illumination position. Namely, the first antenna is activated and the forward problem is solved to calculate and store the electric field at all the remainingreceiving antennas. Each one antenna is activated in turn and the received electric fields are stored to form a complete dataset labeled as "measurements". The reconstruction algorithm starts from a homogeneous model and the desired complex permittivity profile is sought.

*Convergence Criteria:* Two convergence criteria are adopted. The more general concerns the "available information", which is determined by the difference between fields "measured" on the target model (*Em* or *Vm*) and fields calculated at the n-th iterative solution of the forward problem. As in every minimization approach, the sum of squares (SSQ) is the appropriate figure taking signs in to account. Hence the summation over all measurements ports (M) gives:

$$SSQ = \sum\_{i=1}^{M} \left(V\_{meas\_i} - V\_{calc\_i}\right)^2\tag{86}$$

$$SSQ = \sum\_{i=1}^{M} \left( E\_{\text{meas}\_i} - E\_{\text{cat}\_i} \right)^2 \tag{87}$$

The Jacobian matrix is a rectangular one with dimensions M>N where M the number of linearly independent measuring points and N the number of unknown *(σi,εri)* pairs. The study of tis characteristics is quite important especially in the course toward inventing the best data collection strategy, i.e. active and sensing antenna or electrode pairs locations and their subsequent activation. Additionally, the sensitivity matrix study with the aid of the recently revisited Proper Orthogonal Decomposition [43] techniques constitutes an interesting research challenge. Even though our group is already working on this subject it will not be considered herein as it requires a separate chapter of its own. For some preliminary results regarding the SVD analysis of the Jacobian matrix please contact our

The proposed algorithm was applied for the imaging of numerous conductivity and permittivity distributions for both circular and rectangular models and satisfactory reconstructions are obtained. The background and anomaly conductivity and permittivity

Before proceeding to the presentation of reconstructed images the convergence criteria

*Computer Test Approach:* The so called "computer test" was employed in all cases throughout this chapter. First a "target model" is considered and the forward problem is solved for each illumination position. Namely, the first antenna is activated and the forward problem is solved to calculate and store the electric field at all the remainingreceiving antennas. Each one antenna is activated in turn and the received electric fields are stored to form a complete dataset labeled as "measurements". The reconstruction algorithm starts from a homogeneous model and the desired complex permittivity profile

*Convergence Criteria:* Two convergence criteria are adopted. The more general concerns the "available information", which is determined by the difference between fields "measured" on the target model (*Em* or *Vm*) and fields calculated at the n-th iterative solution of the forward problem. As in every minimization approach, the sum of squares (SSQ) is the appropriate figure taking signs in to account. Hence the summation over all measurements

*M*

*i SSQ V V* 

*M*

*i SSQ E E* 

or <sup>2</sup>

 <sup>2</sup> <sup>1</sup> *i i*

(86)

(87)

*meas calc*

<sup>1</sup> *i i*

*meas calc*

**5.3 Sensitivities and features of the Jacobian matrix** 

previous work [30] and wait for al follow up publication.

**6. Numerical results & discussions** 

should be first defined.

is sought.

ports (M) gives:

values resemble these of typical human tissues.

**6.1 Inverse problem convergence criteria** 

For comparison purposes it is more convenient to present its normalized value *(SSQN)* with respect to its maximum *(SSQmax)* occurring at the first iteration of the inverse problem as:

$$SSQ\_N = \frac{SSQ}{SSQ\_{\text{max}}} \tag{88}$$

While SSQ is general and can be calculated in all cases, it is only an indirect indication of convergence. Namely, its minimization ensure *<sup>c</sup>* convergence when a unique solution is safeguarded. But, this condition may be disturbed by the problem singularity degree, which in turn depends on the data collection strategy. For this purpose the well- posedeness or the singularity degree, of the sensitivity-Jacobian matrix should be preliminary examined. However a further elaboration is required. Besides this for the particular case of the "computer test" the target or true *(σ,εr)* distributions are available. Hence, the estimated normalized deviations of , *n n calc calc* from their true values at the n-th iteration can be defined as a norm of relative error:

$$
\sigma - error = \left(\sum\_{i=1}^{P} \left(\sigma\_{true}^{i} - \sigma\_{calc}^{i}\right)^{2} \Bigg/ \sum\_{i=1}^{P} \left(\sigma\_{true}^{i} - \overline{\sigma}\_{true}\right)^{2}\right) \tag{89}
$$

$$
\sigma\_r - error = \left( \sum\_{i=1}^{P} \left( \boldsymbol{\varepsilon}\_{r-true}^i - \boldsymbol{\varepsilon}\_{r-value}^i \right)^2 \Bigg/ \sum\_{i=1}^{P} \left( \boldsymbol{\varepsilon}\_{r-true}^i - \overline{\boldsymbol{\varepsilon}}\_{r-true}^i \right)^2 \right) \tag{90}
$$

where *true* and *r true* are the average values of the target profiles and P is the number of elements of the reconstruction mesh. We should keep in mind that the σ-error, εr-error criteria are not applicable for practical in-vivo or even for the "laboratory test" cases where the objects ,*<sup>r</sup>* distributions are unknown.

#### **6.2 Quasi-static MHz reconstruction**

The algorithm was applied for the circular cross-section of Fig.13a and for the rectangular cross-section of Fig.13b. A total number of 32 electrodes were used, where only two of them are active in each projection angle. An indicative example with a double anomaly is presented herein. The target σ- and εr-profiles are shown in Fig.14a. The values for conductivity are 1 21 / *mS cm* and 2 14 / *mS cm* and for permittivity r1 300 and r2 150 , in a homogeneous background of *<sup>o</sup>* 7 / *mS cm* and ro 100 . The frequency of the injected current was *f* 8*MHz* . The relaxation factors 1 2 *k k*, are set equal to unity. The reconstructed image after 15 iterations are shown in Fig.14b

The algorithm was also tested for the rectangular model of Fig.6. An anomaly with conductivity 1 20 / *mS cm* and permittivity r1 300 was introduced in a homogeneous background of *<sup>o</sup>* 7 / *mS cm* and ro 150 . The frequency of the injected current is now *f* 9*MHz* . The relaxation factors are set equal to unity. The target model and the reconstructed image after 15 iterations are shown in Fig. 15.

For details on the algorithm convergence rate and its performance, please wait for a follow up publication which is now in preparation.

High to Microwave Frequencies Imaging Techniques 189

Electrical Permmitivity (ε

r )

<sup>4</sup> <sup>6</sup> <sup>8</sup>

Conductivity (σ)

cm

<sup>2</sup> <sup>4</sup> <sup>6</sup>

(a) (b) Fig. 15. (a) The target model and (b) the reconstructed σ and εr -profiles after 15 iterations for

The target model simulated as a computer phantom is presented in Fig.16. A total number of 32 antennas (line sources) were used, where only 28 are exploited as receivers for each

frequency of operation was assumed at *f* 1.1*GHz* . The image reconstructed after 9 iterations is presented in Fig.17 for the conductivity and the permittivity profile respectively. The correct location of the anomaly as well as its *σ* and *εr* peak values are

(a) (b)

Fig. 16. Target model (a) conductivity profile, (b) permittivity profile.

<sup>4</sup> <sup>6</sup> <sup>8</sup>

cm

cm <sup>2</sup> <sup>4</sup> <sup>6</sup> <sup>8</sup>

mS/cm

0.15 *S m* and permittivity

. The

*<sup>o</sup> S m* and 30 *ro*

cm

Conductivity (σ)

cm

projection angle. A single offset anomaly with conductivity 1

was introduced in a homogeneous background of 0.3

cm

<sup>4</sup> <sup>6</sup> <sup>8</sup>

the rectangular model.

cm

**6.3 2-D microwave reconstruction** 

obtained, but some artifacts are caused.

mS/cm

cm

<sup>1</sup> 15 *<sup>r</sup>* 

Electrical Permmitivity (ε

r )

Fig. 13. The FEM meshes: (a) fine forward mesh and (b) coarse reconstruction mesh.

Fig. 14. Computer phantom with two anomalies in both εr and σ. (a) The target model and (b) the reconstructed σ- and εr-profiles after 15 iterations.

(a) (b)

(a) (b)

Fig. 14. Computer phantom with two anomalies in both εr and σ. (a) The target model and

(b) the reconstructed σ- and εr-profiles after 15 iterations.

Fig. 13. The FEM meshes: (a) fine forward mesh and (b) coarse reconstruction mesh.

Fig. 15. (a) The target model and (b) the reconstructed σ and εr -profiles after 15 iterations for the rectangular model.

#### **6.3 2-D microwave reconstruction**

The target model simulated as a computer phantom is presented in Fig.16. A total number of 32 antennas (line sources) were used, where only 28 are exploited as receivers for each projection angle. A single offset anomaly with conductivity 1 0.15 *S m* and permittivity <sup>1</sup> 15 *<sup>r</sup>* was introduced in a homogeneous background of 0.3 *<sup>o</sup> S m* and 30 *ro* . The frequency of operation was assumed at *f* 1.1*GHz* . The image reconstructed after 9 iterations is presented in Fig.17 for the conductivity and the permittivity profile respectively. The correct location of the anomaly as well as its *σ* and *εr* peak values are obtained, but some artifacts are caused.

Fig. 16. Target model (a) conductivity profile, (b) permittivity profile.

High to Microwave Frequencies Imaging Techniques 191

(c) (d)

(a) (b)

(c) (d)

conductivity, (c) 2nd layer, (d) 3rd layer permittivity distributions, after 9 iterations. The object is discretized into 400 voxels and is illuminated by 48 dipole sources at a frequency of 1.4 GHz.

Fig. 19. Reconstructed profiles for the 3-D model of Fig.18 (a) 2nd layer, (b) 3rd layer

Fig. 18. Target model comprised of 4 layers to test the 3-D microwave imaging.

Fig. 17. Reconstructed profiles for the example of Fig.16, (a) conductivity and (b) permittivity distributions after 9 iterations. The object is discretized into 100 pixels and is illuminated by 32 line sources at a frequency of 1.1 GHz.

#### **6.4 3-D microwave reconstruction**

The target model comprised of 4 layers as showin in Fig.18. A total number of 48 antennas arranged in three rings of 16 antennas (elementary dipoles) were used, where only 39 are exploited as receivers for each projection angle. A single offset anomaly with conductivity *σ1= 0,15 S/m* and permittivity *εr1 = 15* was introduced at *z = -0.5 cm* (only at the second layer) in a homogeneous background of *σ = 0.6 S/m* and *εr = 60*. The frequency of operation was assumed at *f = 1.4GHz*. The image reconstructed after 9 iterations is presented in Fig.19 for the conductivity and the permittivity profiles respectively.

(a) (b)

The target model comprised of 4 layers as showin in Fig.18. A total number of 48 antennas arranged in three rings of 16 antennas (elementary dipoles) were used, where only 39 are exploited as receivers for each projection angle. A single offset anomaly with conductivity *σ1= 0,15 S/m* and permittivity *εr1 = 15* was introduced at *z = -0.5 cm* (only at the second layer) in a homogeneous background of *σ = 0.6 S/m* and *εr = 60*. The frequency of operation was assumed at *f = 1.4GHz*. The image reconstructed after 9 iterations is presented in Fig.19 for

(a) (b)

permittivity distributions after 9 iterations. The object is discretized into 100 pixels and is

Fig. 17. Reconstructed profiles for the example of Fig.16, (a) conductivity and (b)

illuminated by 32 line sources at a frequency of 1.1 GHz.

the conductivity and the permittivity profiles respectively.

**6.4 3-D microwave reconstruction** 

Fig. 18. Target model comprised of 4 layers to test the 3-D microwave imaging.

Fig. 19. Reconstructed profiles for the 3-D model of Fig.18 (a) 2nd layer, (b) 3rd layer conductivity, (c) 2nd layer, (d) 3rd layer permittivity distributions, after 9 iterations. The object is discretized into 400 voxels and is illuminated by 48 dipole sources at a frequency of 1.4 GHz.

High to Microwave Frequencies Imaging Techniques 193

[12] B. H. Brown, "Tissue impedance methods," in *Imaging with non-ionizing radiations*, Ann.

[13] C. S. Lavranos and G. A. Kyriacou, "Eigenvalue analysis of curved waveguides

[14] S. Lavdas, C. Lavranos, and G. A. Kyriacou, "A finite difference frequency domain

[16] J. L. Volakis, A. Chatterjee, and L. C. Kempel, *Finite Element Method for Electromagnetics*:

[17] Y. Zhu and A. C. Cangellaris, *Multigrid Finite Element Methods for Electromagnetic Field* 

[18] L. R. Price, "Electrical Impedance Computed Tomography (ICT): A New CT Imaging Technique", *IEEE Transactions on Nuclear Science,* vol. 26, pp. 2736-2739, 1979

[20] C. Zekios, P. Allilomes, and G. Kyriacou, "Eigenfunction expansion for the analysis of

[21] S. G. Diamantis, A. P. Orfanidis, G. A. Kyriacou, and J. N. Sahalos, "Horn Antennas

[22] P. C. Allilomes and G. A. Kyriacou, "A Nonlinear Finite-Element Leaky-Waveguide

[25] Q. Fang, "Computational methods for microwave medical imaging," Ph.D. Thesis,

[26] G. Kyriacou, C. Koukourlis, and J. Sahalos, "A reconstruction algorithm of electrical

[27] D. Drogoudis, G. Trichopoulos, and G. A. Kyriacou, "A modified pertubation method

[28] D. Drogoudis, "The inverse electromagnetic problem in high frequency tomography," Ph.D. Thesis, Democritus University of Thrace, Dec. 2009 (in Greek), http://thesis.ekt.gr/thesisBookReader/id/18153#page/4/mode/2up [29] D. Drogoudis, G. A. Kyriacou, and J. N. Sahalos, "A three dimensional microwave

http://ieeexplore.ieee.org/xpls/abs\_all.jsp?arnumber=4268422&tag=1 [23] P. Allilomes, "Electromagnetic Simulation of open-radiating structures using the Finite

[24] C. A. Balanis, *Advanced Engineering Electromagnetics*: John Wiley & Sons, 1989.

Dartmouth College Hanover, New Hampshire, 2004,

*Trans. Medical Imaging,* vol. 12, pp. 430-438, Sept. 1993

*Trans. Magnetics,* vol. 45, pp. 1686-1689, 2009,

http://ieeexplore.ieee.org/xpls/abs\_all.jsp?arnumber=5666915

closed cavities," presented at the Loughborough Antennas and Propag. Conf., 2010,

Analysis Using a Hybrid Mode Matching-Auxiliary Sources Technique " in *Proc. of* 

Solver", *IEEE Transactions on Microwave Theory and Techniques,* vol. 55, pp. 1496-

Element method," Ph.D. Thesis, Dept. of El. & Comp. Eng, Democtirus University

impedance tomography with optimal configuration of the driven electrodes", *IEEE* 

for three-dimensional time harmonic impedance tomography", *PIERS Online,* vol.

tomography exploiting an adjoint network theorem based sensitivity matrix", *IEEE* 

*IEEE Microwave Theory and Techniques,,* vol. 57, pp. 594-611, March 2009, http://ieeexplore.ieee.org/stamp/stamp.jsp?arnumber=04783093

[15] J. Jin, *The finite element method in electromagnetics*: John Wiley & Sons, 1993.

employing an orthogonal curvilinear frequency domain finite difference method",

method for the eigenanalysis of anisotropically loaded curved periodic structures," in *32nd Antenna Workshop on Antennas for Space Applications*, ESA/ESTEC,

of the NYAS, 1983,

Noordwijk, 2010,

IEEE Press, 1998.

1510, July 2007,

of Thrace, Xanthi, 2007,

1, pp. 151-155, 2005

*Modeling*: IEEE Press, 2006.

[19] J. D. Jackson, *Classical Electrodynamics*, 3rd ed.: Wiley, 1999.

*PIERS* PISA March 2004, pp. 457-460,

For details on the 3-D algorithm performance and its convergence characteristics please contact our work [30].

#### **7. Conclusions**

A modified perturbation method reconstruction algorithm for MHz and microwave tomography is successfully established at a computer test level. The key constituent of this algorithm is a close form evaluation of the sensitivity or Jacobian matrix based on an adjoint network and reciprocity theorem approach. For the MHz frequency range th Jacobian is based on an extension of the electrical networks compensation theorem but it can be obtained as a particular case fo the general Adjoint network approach. The established algorithm is proved robust and able to withstand a large amount of inverse problem ill-posedeness. Despite its simplicity this algorithm is able to successfully localize the target anomalies reaching conductivity and permittivity patterns very close to the expected global minimum at about the 6th iteration. The penalty paid for this simplicity and robustness is a compromise in the finally achieved solution mostly appearing as artefacts around the target anomalies. A further improvement, which also constitutes our next task, refers to the formulation of an exact inverse problem, which may start from the finally obtained image herein and iteratively fine tune it. This can be readily formulated exploiting the exact Jacobian matrix established herein.

#### **8. References**


For details on the 3-D algorithm performance and its convergence characteristics please

A modified perturbation method reconstruction algorithm for MHz and microwave tomography is successfully established at a computer test level. The key constituent of this algorithm is a close form evaluation of the sensitivity or Jacobian matrix based on an adjoint network and reciprocity theorem approach. For the MHz frequency range th Jacobian is based on an extension of the electrical networks compensation theorem but it can be obtained as a particular case fo the general Adjoint network approach. The established algorithm is proved robust and able to withstand a large amount of inverse problem ill-posedeness. Despite its simplicity this algorithm is able to successfully localize the target anomalies reaching conductivity and permittivity patterns very close to the expected global minimum at about the 6th iteration. The penalty paid for this simplicity and robustness is a compromise in the finally achieved solution mostly appearing as artefacts around the target anomalies. A further improvement, which also constitutes our next task, refers to the formulation of an exact inverse problem, which may start from the finally obtained image herein and iteratively fine tune it. This can be readily formulated exploiting the exact Jacobian matrix established herein.

[1] G. V. Keller, "Electrical properties of rocks and minerals", *Handbook of Physical Constants* 

[2] R. Pethig, "Dielectric properties of biological materials: Biophysical and medical

[3] C. Gabriel, S. Gabriel, and E. Corthout, "The dielectric properties of biological tissues:

[4] Q. Fang, P. M. Meaney, S. D. Geimer, A. V. Streltsov, and K. D. Paulsen, "Microwave

[5] P. M. Meaney, K. D. Paulsen, A. Hartov, and R. K. Crane, "Microwave imaging for tissue

[6] S. H. Zainud-Deen, W. M. Hassen, E. E. D. ali, and K. H. Awadalla, "Breast cancer detection

[8] W. R. B. Lionheart, J. Kaipio, and C. N. McLeond, "Generalized optimal current patterns and electrical safety in EIT", *IOP publ, Physiological Meas.,* vol. 22, pp. 85-90, 2001 [9] D. C. Barber and A. D. Seagar, "The Sheffiled data collection system", *Clin. Phys. &* 

[11] A. V. Korzhnenevskii, V. N. Kornienko, M. Y. Kultiasov, Y. S.Kultiasov, and V. A.

Applications", *Instruments for Experimental Tech.,* vol. 40, pp. 415-421, 1997

techniques", *Progress in Electromagnetics Research B,* vol. 3, pp. 35-46, 2008 [7] J. Goble and D. Isaacson, "Optial current patterns for three-dimensional computed

applications", *IEEE Transactions on Electrical Insulation,* vol. EI, pp. 453-474, Oct.1984

image reconstruction from 3-D fields coupled to 2-D parameter estimation", *IEEE* 

assessment: initial evaluation in multi-target tissue-equivalent phantoms", *IEEE* 

using a hybrid finite difference frequency domain and particle swarm optimization

tomography," in *Proc. Int. conf. IEEE Trans. Eng. in medicine & biology soc.*, 1989, pp.

Cherepenin, "Electrical Impedance Computerized Tomography for Medical

*(S.P. Clark, ed.),N.Y., Geological Society of America*, pp. 553-770, 1988

Part I, II and III", *Phys. Med. Biol.,* vol. 41, 1996

*Trans. Medical Imaging,* vol. 23, pp. 475-484, 2004

*Trans. Biomed. Eng.,* vol. 43, pp. 878-890, 1996

*Physiol. Meas. ,* vol. 8, pp. 91-97, 1987

[10] R. E. Collin, *Filed Theory of Guided Waves*, 2nd ed.: IEEE Press, 1991.

contact our work [30].

**7. Conclusions** 

**8. References** 

463-464,


http://ieeexplore.ieee.org/xpls/abs\_all.jsp?arnumber=4268422&tag=1


**9** 

*Niigata University* 

*Japan* 

**A Mutual Information-Based Image Quality** 

Information on physical image quality of medical images is important for imaging system assessment in order to promote and stimulate the development of state-of-the-art imaging systems. In this chapter, we present a method for quantifying overall image quality of digital imaging systems using mutual information (MI) metric. The MI which is a concept from information theory is used as a measure to express the amount of information that an output image contains about an input object. The MI value is considered that it can be used to express combined physical properties of image noise, resolution and contrast of an imaging system. The higher the MI value, the better the image quality. The advantages of using the MI metric are: (1) simplicity of computation, (2) simplicity of experimentation, and

The structure of this chapter is as follows. Section \*\*.2 provides a basic overview of factors that affect medical image quality. Section \*\*.3 describes the mutual information-based evaluation framework utilized in this work. An example of how to calculate MI is also given to provide a deep understanding of applying MI to the evaluation of medical imaging systems. Section \*\*.4 shows a series of computer simulations, followed by investigating the utility and superiority of MI method by evaluating the performance of two imaging-plate detectors. Section \*\*.5 presents the results that were obtained. Section \*\*.6 ends with a

In medical imaging, image quality is determined by at least five factors: contrast, resolution, noise, artifacts, and distortion. Of these factors, resolution and noise are the most commonly used physical characteristics. As is well known, they are described by the modulation transfer function (MTF) and noise power spectrum (NPS), respectively. The MTF describes the ability of an imaging system to reproduce the frequency information contained in the incident x-ray signal. The NPS describes the frequency content of the noise of an imaging system. However, one of the dilemmas in medical radiography is the extent to which theses characteristics affect image quality. In comparison of two imaging systems, for example, an imaging system may only be superior in one physical characteristic while being inferior to another in the other characteristic. To deal with this issue, the noise equivalent quanta or detective quantum efficiency (DQE), which can be calculated if the MTF, NPS, and the input

(3) combined assessment of image contrast, noise and resolution.

**1. Introduction** 

discussion and conclusions.

**2. Background** 

**Metric for Medical Imaging Systems** 

Du-Yih Tsai, Eri Matsuyama and Yongbum Lee

http://piers.org/piersonline/pdf/Vol3No8Page1217to1221.pdf

