Home Search Collections Journals About Contact us My IOPscience

Solution of the inverse problem of magnetic induction tomography (MIT)

This article has been downloaded from IOPscience. Please scroll down to see the full text article.

2005 Physiol. Meas. 26 S241

(http://iopscience.iop.org/0967-3334/26/2/023)

View the table of contents for this issue, or go to the journal homepage for more

Download details:
IP Address: 134.99.128.41
The article was downloaded on 08/04/2013 at 16:32

Please note that terms and conditions apply.

INSTITUTE OF PHYSICS PUBLISHING PHYSIOLOGICAL MEASUREMENT
Physiol. Meas. 26 (2005) S241–S250 doi:10.1088/0967-3334/26/2/023

Solution of the inverse problem of magnetic induction

tomography (MIT)

Robert Merwa1, Karl Hollaus2, Patricia Brunner1

and Hermann Scharfetter1
1 Institute for Medical Engineering, Graz University of Technology, Krenngasse 37,

8010 Graz, Austria

2 Institute for Fundamentals and Theory in Electrical Engineering, Graz University of

Technology, Kopernikusgasse 24, 8010 Graz, Austria

E-mail: [email protected]

Received 13 October 2004, accepted for publication 29 November 2004

Published 29 March 2005
Online at stacks.iop.org/PM/26/S241

Abstract
Magnetic induction tomography (MIT) of biological tissue is used to
reconstruct the changes in the complex conductivity distribution inside an
object under investigation. The measurement principle is based on determining
the perturbation
B of a primary alternating magnetic field B0 , which is
coupled from an array of excitation coils to the object under investigation.
The corresponding voltages
V and V0 induced in a receiver coil carry
the information about the passive electrical properties (i.e. conductivity,
permittivity and permeability). The reconstruction of the conductivity
distribution requires the solution of a 3D inverse eddy current problem. As in
EIT the inverse problem is ill-posed and on this account some regularization
scheme has to be applied. We developed an inverse solver based on the Gauss–
Newton-one-step method for differential imaging, and we implemented and
tested four different regularization schemes: the first and second approaches
employ a classical smoothness criterion using the unit matrix and a differential
matrix of first order as the regularization matrix. The third method is based
on variance uniformization, and the fourth method is based on the truncated
singular value decomposition. Reconstructions were carried out with synthetic
measurement data generated with a spherical perturbation at different locations
within a conducting cylinder. Data were generated on a different mesh and
1% random noise was added. The model contained 16 excitation coils and
32 receiver coils which could be combined pairwise to give 16 planar
gradiometers. With 32 receiver coils all regularization methods yield fairly
good 3D-images of the modelled changes of the conductivity distribution,
and prove the feasibility of difference imaging with MIT. The reconstructed
perturbations appear at the right location, and their size is in the expected
range. With 16 planar gradiometers an additional spurious feature appears
mirrored with respect to the median plane with negative sign. This demonstrates
0967-3334/05/020241+10$30.00 © 2005 IOP Publishing Ltd Printed in the UK S241
S242 R Merwa et al

that a symmetrical arrangement with one ring of planar gradiometers cannot

distinguish between a positive conductivity change at the true location and a
negative conductivity change at the mirrored location.

Keywords: magnetic induction tomography, sensitivity, inverse problem,

regularization, gradiometer

1. Introduction

Magnetic induction tomography (MIT) is a noninvasive and contactless imaging method for
reconstructing the changes of the complex conductivity distribution
κ in a target object
(Griffiths et al 1999, Griffiths 2001, Korzhenevsky and Cherepenin 1997, 1999, Korjenevsky
et al 2000, Peyton 1996). This technique when applied at multiple frequencies appears
especially attractive for the monitoring of pathologies in the brain, which are correlated with
local fluid shifts (Merwa et al 2004), e.g. oedema, haemorrhages or epileptic events. MIT
requires an array of excitation and receiving coils. Each excitation coil couples an alternating
magnetic field B0 to the object under investigation. Changes
κ of the complex conductivity
κ = σ + j ω in a target region cause a field perturbation
B due to the induction of eddy
currents and magnetic dipoles in the object under investigation. The perturbation induces
voltage changes
V in the receiver coils. A difficulty in MIT is that the signal due to the eddy
currents in the sample is superimposed on the primary signal from the applied magnetic field.
Subtraction of the primary signal in order to improve the sensitivity of measurement has been
achieved using a planar gradiometer (Casanas et al 2004, Rosell et al 2001, Scharfetter et al
2004).
The reconstruction of the absolute conductivity in a target region requires the solution of
a complex inverse eddy current problem. Let
V = Ψ(κ) (1)
be the discretized nonlinear forward mapping (Merwa et al 2003, Hollaus et al 2002) of the
conductivity vector κ to the vector of induced voltages V. V contains a • b entries, a being the
number of excitation coils and b that of receiver coils. The corresponding inverse problem
κ = Ψ−1 (V) (2)
is ill-posed and usually underdetermined. The generic approach for the solution of this type
of nonlinear problem is the application of an iterative scheme such as the regularized Gauss–
Newton approach, including appropriate regularization methods. This paper is dedicated to
demonstrating the feasibility of the 3D reconstruction of a spherical perturbation within a
cylindrical conducting body in consideration on the one hand of normal coils and on the other
hand of planar gradiometers.

2. Methods

The solution of the inverse problem in (2) requires the discretization of the target region
into N voxels. Within each voxel i the assigned component κi of the vector of conductivity
κ is assumed to be constant. In the presented work a grid of tetrahedral finite elements
was employed. In general κ is found with an iterative scheme but this procedure is very
time-consuming, and a complete iteration requires significant computer power. However, in
Solution of the inverse problem in MIT S243

EIT it could be shown that in practice most features of the image can already be recognized
satisfactorily after the first iteration step. This fact led to the development of the Newton-one-
step reconstructor (Cheney et al 1990). It is especially appropriate for dynamical imaging
where only the change in the conductivity between two different states of the object under
investigation is of interest. In this case the first Newton step corresponds to the solution of the
linearized forward problem
dΨ

V =
κ = S
κ (3)
dκ
hence

κ = (ST S)−1 ST
V (4)
and if a regularization matrix R is implemented, (4) becomes

κ = (ST S + λRT R)−1 ST
V (5)
where
κ is the change of the conductivity between two states of the observed object,
V
is the corresponding change of the measured data, and λ is a regularization parameter. The
inversion requires an efficient calculation of the sensitivity matrix S thus an extension of the
Geselowitz relationship (Geselowitz 1971) to the general electromagnetic case (Mortarelli
1980) is exploited for this purpose. The matrices S and R contain m × n entries, m being the
number of independent measurements, and n is the number of reconstructed conductivities
(finite elements). A detailed description for the calculation of the sensitivity matrix is given
by Hollaus et al (2004).
In the case of comparatively small changes of the conductivity the inversion of (3)
according to (5) yields a fairly correct localization of the perturbed regions.
In order to demonstrate the feasibility of this kind of reconstruction in MIT we
implemented a Newton-one-step reconstructor approach according to (5) under consideration
of four different regularization methods. An alternative way for dynamical imaging by means
of Kalman filter approach has been already suggested by Vauhkonen et al (1998).

2.1. Regularization
In EIT the regularization matrix RT R most frequently is either the unit matrix I or a discrete
spatial derivative operator of first or second order. The latter two can be regarded as smoothness
criteria for the solution, and such approaches have been discussed extensively in the literature.
Even better results can be achieved with variance uniformization (Cohen-Bacrie et al 1997),
choosing R according to statistical a priori information. The objective is to uniformize
the expected variance of the reconstructed conductivity changes over the target region, thus
providing approximately equal image noise at the centre as in the periphery. The algorithm
has been described in detail in Cohen-Bacrie et al (1997).
Alternatively also truncated singular value decomposition (TSVD) is very common in EIT
reconstruction, hence this approach was also implemented for MIT. In this case the inverse
solution becomes

κ = Wt Σ−1 T
t Ut
V (6)
where t denotes the truncation level of the original matrices U, Σ and W, thus removing the
contributions of singular values with index >t.
In this paper the results obtained with four different regularization schemes were
compared:
(a) RT R = E. Using the unit matrix is the most simple Tikhonov-regularization method,
penalizing high values of the reconstructed conductivity changes.
S244 R Merwa et al

(b) RT R = N. N is the neighbouring matrix. N is defined as



nn for i = j
Nij = −1 for i, j neighbours (7)


0 otherwise
nn is the number of neighbouring elements for element i, whereas only elements with
common faces are considered as neighbours. N is an approximation of the spatial
derivative operator of first order.
(c) λRT R = WDWT . According to the variance uniformization approach which implies the
singular value decomposition of S = UΣWT and a diagonal matrix D which is defined as
Σ
D = √ − Σ2 (8)
c
where c is a ‘free’ scalar parameter.
(d) TSVD. While choosing the truncation level t so as to remove all singular values below
the noise level of the measurement data.
Methods 1 and 2 require the regularization parameter λ to be chosen appropriately while
method 3 implies the choice of the tuning parameter c. In practice it turns out that the choice
of this parameter is uncritical over a very wide range of values, and c was chosen to 0.1.
It is not possible to calculate the sensitivity matrix frequently on a single PC due to the
time-consuming calculation thus the NOSER algorithm (Cheney et al 1990) was implemented.
Hence it is not possible to calculate the regularization parameter λ with the method of L-curves.
The regularization parameter λ accounts for the degree of smoothness of the reconstructed
image and determines the condition number of the term (ST S + λRT R) in (5). Selecting the
‘optimal’ λ manually by visual inspection of the solution is very time-consuming because
this parameter varies in a wide range (10−11 to 10−22 ). It is much more efficient to select
the ‘optimal’ condition number manually by visual inspection and to calculate the associated
λ with the aid of a numerical optimization procedure. The optimal condition number can
be found with a small number of test simulations because it is within a range of 3 decimal
powers. Moreover the optimal condition number is independent of the frequency of the
alternating magnetic field B0 in contrast to the regularization parameter λ. Another method
to determine the regularization parameter λ for large problems is to use the generalized cross
validation approach described by Hanke and Hansen (1993).

2.2. Modelling set-up

The inverse solver was tested with a simple 3D model depicted in figure 1 comprising a
cylindrical conductor with a spherical perturbation and an array of 16 excitation coils and
32 receiving coils (16 planar gradiometers). The outer cylinder had a radius of 96 mm and
a height of 140 mm, and the perturbation had a radius of 20 mm. The solenoid excitation
coils had a diameter of 155 mm, a height of 21.5 mm and a thickness of 34 mm. They
were placed with their centres onto a ring with a radius of 152 mm in the median transversal
plane of the cylinder. The rectangular receiver coils with an edge length of 38.7 and 20 mm
were placed with their centres on two symmetrically arranged parallel rings with a radius
of 117 mm, each comprising 16 evenly spaced coils. The measured data were simulated
when changing the conductivity of the spherical perturbation from 1 S m−1 (=background
conductivity, homogeneous cylinder, initial state) to 2 S m−1. The permittivity was kept
constant with 80 in the cylinder and in the sphere. The excitation frequency was 100 kHz.
Solution of the inverse problem in MIT S245

Figure 1. Simulation arrangement for solving the inverse problem and generating artificial
measurement data: (a) front view, (b) top view (all measurements in mm).

Table 1. Tuning parameters and corresponding regularization parameters for the used
regularization methods at a frequency of 100 kHz using 32 receiver coils.

Regularization methd Condition number λ Cutting level

Unit matrix 10−4 1.77 × 10−20 –

Neigbouring matrix 10−5 2.14 × 10−21 –
Variance uniformization 10−2 4.25 × 10−11 –
TSVD – – 160

Two different meshes A and B were used for the generation of an artificial dataset and for
the reconstruction (calculation of the sensitivity matrix). Mesh A comprised 16 000 elements
within the cylinder and 400 elements within the sphere while mesh B comprised 14 000
elements for the homogeneous cylinder without perturbation. In both cases the diameter of
the surrounding spherical surface which approximated the far boundary was chosen with 1 m,
requiring to a number of approximately 54 000 elements. One per cent uncorrelated Gaussian
noise was added to the voltage data in order to simulate the noise of the receiver channels.
The gradiometers are defined as the coils spaced one above the other (see figure 1(a)), and the
resulting voltage difference
VGrad in one gradiometer is calculated as follows:

V =
VGrad =
Vcoil up −
Vcoil down (9)
where
Vcoil up is the induced voltage difference in the corresponding upper coil, and
Vcoil down
is the induced voltage difference in the corresponding lower coil.

3. Results

Figures 2–5 show the reconstructed images using 32 receiver coils for the four different
positions of the perturbing sphere including the four different regularization methods. The
tuning parameters were chosen as listed in table 1. Figures 6–9 show the reconstructed
images using 16 planar gradiometers and the corresponding tuning parameters are listed in
table 2.
S246 R Merwa et al

Figure 2. Reconstruction of
κ for the position of the sphere at [x, y, z] = [30, 0, 23.75] with
(a) unit matrix, (b) neighbouring matrix, (c) variance uniformization, (d) TSVD. Left: top view,
right: front view, with 32 coils.

Figure 3. Reconstruction of
κ for the position of the sphere at [x, y, z] = [50, 0, 23.75] with
(a) unit matrix, (b) neighbouring matrix, (c) variance uniformization, (d) TSVD. Left: top view,
right: front view, with 32 coils.

Table 2. Tuning parameters and corresponding regularization parameters for the used
regularization methods at a frequency of 100 kHz using 16 planar gradiometers.

Regularization methd Condition number λ Cutting level

Unit matrix 10−4 2.50 × 10−21 –

Neigbouring matrix 10−5 3.10 × 10−22 –
Variance uniformization 10−2 1.60 × 10−11 –
TSVD – – 100

Because of the non-iterative nature of the algorithm, the absolute values of the results
cannot be interpreted quantitatively, because the first iteration step in (5) yields the correct
search direction but usually a wrong step width, due to the lack of a line search. On this
account the colourmaps are not quantified only the location and the visual appearance of the
Solution of the inverse problem in MIT S247

Figure 4. Reconstruction of
κ for the position of the sphere at [x, y, z] = [30, 0, 10] with
(a) unit matrix, (b) neighbouring matrix, (c) variance uniformization, (d) TSVD. Left: top view,
right: front view, with 32 coils.

Figure 5. Reconstruction of
κ for the position of the sphere at [x, y, z] = [50, 0, 10] with
(a) unit matrix, (b) neighbouring matrix, (c) variance uniformization, (d) TSVD. Left: top view,
right: front view, with 32 coils.

Figure 6. Reconstruction of
κ for the position of the sphere at [x, y, z] = [30, 0, 23.75] with
(a) unit matrix, (b) neighbouring matrix, (c) variance uniformization, (d) TSVD. Left: top view,
right: front view, with 16 gradiometers.
S248 R Merwa et al

Figure 7. Reconstruction of
κ for the position of the sphere at [x, y, z] = [50, 0, 23.75] with
(a) unit matrix, (b) neighbouring matrix, (c) variance uniformization, (d) TSVD. Left: top view,
right: front view, with 16 gradiometers.

Figure 8. Reconstruction of
κ for the position of the sphere at [x, y, z] = [30, 0, 10] with
(a) unit matrix, (b) neighbouring matrix, (c) variance uniformization, (d) TSVD. Left: top view,
right: front view, with 16 gradiometers.

Figure 9. Reconstruction of
κ for the position of the sphere at [x, y, z] = [50, 0, 10] with
(a) unit matrix, (b) neighbouring matrix, (c) variance uniformization, (d) TSVD. Left: top view,
right: front view, with 16 gradiometers.
Solution of the inverse problem in MIT S249

reconstructed perturbation were assessed. The black-dotted lines in the figures demonstrate
the original position of the disturbed sphere.

4. Discussion

The shown results demonstrate the feasibility of image reconstruction from MIT data with
the same NOSER method as used for EIT. The main difference is in the underlying partial
differential equation thus leading to significantly different sensitivity distribution in EIT and
MIT (Scharfetter et al 2002, 2004).
The reconstruction if 32 receiver coils are used is depicted in figures 2–5, and the images
show a good localization of the perturbed sphere. All perturbations are slightly displaced
towards the nearest border of the cylinder whereas this effect is clearly visible in figures 3(c)
and 5(c). The first reason may be as shown by Scharfetter et al (2002, 2004) that in MIT
the sensitivity is not concentrated within a field tube between excitation and receiver coil but
increases with the distance from the tube axis, according to the increase of the eddy current
density. The second reason may be that the inverse problem is solved by means of a linear
approximation which means that the sensitivity matrix for the first iteration step corresponds
to the unperturbed distribution. But the sensitivity matrix itself is affected by the perturbation,
so the correct position could be achieved with a nonlinear solver until the next iterations with
the corrected sensitivity matrix.
To improve the dynamic range and to obtain maximum sensitivity to changes in the
conductivity, it is essential to decrease the induced voltage V0 in the receiver coils. One way
is to subtract the signals of two coils called a gradiometer. The reconstruction with 16 planar
gradiometers is shown in figures 6–9. In contrast to the use of 32 receiver coils, the localization
of the perturbed sphere is better in this case but the disadvantage is that the reconstructed
perturbation appears at two places. One part at the right location (black-dotted line), and the
other (white-dotted line) down mirror-inverted with negative conductivity changes.
In fact the planar gradiometer arrangement cannot distinguish between the upper and
lower objects, because it is always possible to find two corresponding conductivities so that
the induced voltages
V in the planar gradiometers are the same. Obviously, in this ambiguous
situation, the algorithm favours the splitted solution according to the sensitivity distribution.
This splitting disappears immediately when placing two parallel rings of planar gradiometers
around the cylinder (Scharfetter et al 2005).
When comparing the regularization schemes, the unit matrix and the variance
uniformization approach yield the smoothest solutions, as expected. The most homogeneous
background is also obtained with variance uniformization, according to the variance constraint
based on this regularization scheme. The condition number in order to optimize the
regularization parameter λ is the same for the reconstructions with 32 receiver coils and
with 16 planar gradiometers. Differences are only in the case of the TSVD. The cutting level
is 160 if 32 receiver coils are used and this level decreases in the case of planar gradiometers
because the number of independent equations in (6) decreases.
The shown results were obtained at a single frequency only. Future work should
concentrate on the exploitation of the frequency dependence of the tissue conductivity and
measurements at frequencies up to several MHz (Scharfetter et al 2001, 2003). Multifrequency
measurements can decrease the degree of ill-conditioning of the inverse problem if additional
a priori information is included in the form of spectral tissue models. The feasibility of
this approach has been recently demonstrated for multifrequency EIT in combination with
Cole-models in 2D in (Brandstätter et al 2003).
S250 R Merwa et al

Acknowledgments

This work was supported by the Austrian Science Foundation FWF, project P14990.

References

Brandstätter B, Hollaus K, Hutten H, Mayer M, Merwa R and Scharfetter H 2003 Direct estimation of Cole parameters
in multifrequency EIT using a regularized Gauss–Newton method Physiol. Meas. 24 437–48
Casanas R, Scharfetter H, Altes A, Remacha A, Sarda P, Sierra J, Merwa R, Hollaus K and Rosell J 2004 Measurement
of liver iron overload by magnetic induction using a planar gradiometer: preliminary human results Physiol.
Meas. 25 315–23
Cheney M, Isaacson D, Newell J C, Simske S and Goble J 1990 NOSER: an algorithm for solving the inverse
conductivity problem Int. J. Imaging Syst. Technol. 2 66–75
Cohen-Bacrie C, Goussard Y and Guardo R 1997 Regularized reconstruction in electrical impedance tomography
using a variance uniformization constraint IEEE Trans. Med. Imaging 16/5 562–71
Geselowitz D B 1971 An application of electrocdiographic lead theory to impedance plethysmography IEEE Trans.
Biomed. Eng. 18 38-41
Griffiths H 2001 Magnetic induction tomography Meas. Sci. Technol. 12 1126–31
Griffiths H, Steward W R and Gough W 1999 Magnetic induction tomography: a measuring system for biological
tissues Ann. New York Acad. Sci. 873 335–45
Hanke M and Hansen P C 1993 Regularization methods for large-scale problems Surv. Math. Ind. 3 253–315
Hollaus K, Magele C, Brandstätter B, Merwa R and Scharfetter H 2002 Numerical simulation of the forward problem
in magnetic induction tomography of biological tissue Proc. 10th IGTE Symp. 25 159–68
Hollaus K, Magele C, Merwa R and Scharfetter H 2004 Fast calculation of the sensitivity matrix in magnetic induction
tomography by tetrahedral edge finite elements and the reciprocity theorem Physiol. Meas. 25 159–68
Korjenevsky A and Cherepenin V 1999 Progress in realization of magnetic induction tomography Ann. NY Acad. Sci.
873 346–52
Korjenevsky A, Cherepenin V and Sapetsky S 2000 Magnetic induction tomography: experimental realization Physiol.
Meas. 21 89–94
Korzhenevskii A and Cherepenin V 1997 Magnetic induction tomography J. Commun. Technol. Electron. 42 469–74
Merwa R, Hollaus K, Biro O and Scharfetter H 2004 Detection of brain oedema using magnetic induction tomography:
a feasibility study of the likely sensitivity and detectability Physiol. Meas. 25 347–54
Merwa R, Hollaus K, Brandstätter B and Scharfetter H 2003 Numerical solution of the general 3D eddy current
problem for magnetic induction tomography (spectroscopy) Physiol. Meas. 24 545–54
Mortarelli J R 1980 A generalization of the Geselowitz relationship useful in impedance plethysmographic field
calculations IEEE Trans. Biomed. Eng. 27 665–7
Peyton A 1996 An overview of electromagnetic inductance tomography: description of three different systems Meas.
Sci. Technol. 7 261–71
Rosell J, Casanas R and Scharfetter H 2001 Sensitivity maps and system requirements in magnetic induction
tomography using a planar gradiometer Physiol. Meas. 22 121–30
Scharfetter H, Casanas R and Rosell J 2003 Biological tissue characterization by magnetic induction spectroscopy
(MIS): requirements and limitations IEEE. Trans. Biomed. Eng. 50 870–80
Scharfetter H, Lackner H and Rosell J 2001 Magnetic induction tomography: hardware for multi-frequency
measurements in biological tissues Physiol. Meas. 22 131–46
Scharfetter H, Merwa R and Pilz K 2005 A new type of gradiometer for the receiving circuit of magnetic induction
tomography (MIT) Physiol. Meas. 26 S307–18
Scharfetter H, Rauchenzauner S, Merwa R, Biro O and Hollaus K 2004 Planar gradiometer for magnetic induction
tomography (MIT): theoretical and experimental sensitivity maps for a low-contrast phantom Physiol. Meas. 25
325–33
Scharfetter H, Riu P, Populo M and Rosell J 2002 Sensitivity maps for low-contrast perturbations within conducting
background in magnetic induction tomography Physiol. Meas. 23 195–202
Vauhkonen M, Karjalainen P A and Kaipio J P 1998 A Kalman filter approach to track fast impedance changes in
electrical impedance tomography IEEE. Trans. Biomed. Eng. 45 486–93