A non-iterative method for electrical property tomography based on a simple formula

TL;DR

This paper introduces a non-iterative, analytic formula-based method for reconstructing electrical properties from MRI data, simplifying the process by avoiding iterative algorithms.

Contribution

It presents a novel, exact analytic approach for electrical property tomography using harmonic electromagnetic field data at Larmor frequency.

Findings

Provides a simple, non-iterative reconstruction formula

Offers a geometric method for acquiring the full inductive field

Enhances efficiency and accuracy of electrical property imaging

Abstract

A non-iterative method of reconstruction is proposed from data of MRI system and of a harmonic electro-magnetic field at Larmor frequency. The method is based on the exact analytic formula for the contrast source function. A geometric method for acquisition of the full inductive field is discussed.

Figures

Equations

$$E=\mathbf{E}+E^{sc},H=\mathbf{H}+H^{sc},$$

$$E^{sc}=\left({\mathbf{k}}^2+\nabla \nabla \cdot \right)A,$$

$$H^{sc}={\gamma }_0\nabla \times A$$

$$A\left(p\right)={\int }_{\Omega }G\left(p-q\right)\chi \left(q\right)E\left(q\right)\mathrm{d}V\left(q\right),q,p\in {\mathbb{R}}^3$$

$$G\left(q\right)=-\frac{\exp \left(ik\left|q\right|\right)}{4\pi \left|q\right|},k=\frac{\omega }{c_0}$$

$$\chi \left(q\right)=\frac{\gamma }{{\gamma }_0}-1$$

$$E=E_x\mathrm{d}x+...,B=B_x\mathrm{d}x+...$$

$${\alpha }^{\ast }\wedge \beta$$

$${\left(a\mathrm{d}x\right)}^{\ast }$$

$$\chi =-\frac{\left(\mathbf{E}-\ast \mathrm{d}B\right)\wedge \left(\Delta +k^2\right)B^{\ast }}{\left(\mathbf{E}-\ast \mathrm{d}B\right)\wedge \left(\mathrm{d}\mathbf{E}+k^2{B}^{\ast }\right)}$$

$$\frac{\gamma }{{\gamma }_0}=\chi +1=\frac{⟨-\Delta B^{\ast }+\mathrm{d}\mathbf{E},\ast \mathrm{d}B-\mathbf{E}⟩}{⟨k^2{B}^{\ast }+\mathrm{d}\mathbf{E},\ast \mathrm{d}B-\mathbf{E}⟩}.$$

$$\gamma =\frac{1}{i\omega {\mu }_0}\frac{⟨{\nabla }^{2}H,\nabla \times H⟩}{⟨H,\nabla \times H⟩}$$

$$\mathrm{d}C=B^{\ast },{\mathrm{d}}^{\ast }C=0$$

$$C+\mathrm{d}g=A_0.$$

$$\Delta f=-{\mathrm{d}}^{\ast }{C}_0$$

$${\mathrm{d}}^{\ast }C={\mathrm{d}}^{\ast }{C}_0+{\mathrm{d}}^{\ast }\mathrm{d}f={\mathrm{d}}^{\ast }{C}_0+\Delta f=0.$$

$$B^{\ast }=\mu {\left(H^{sc}\right)}^{\ast }=\frac{i\omega }{c^2}\mathrm{d}A=\mathrm{d}{A}_0,$$

$$A=G_{\chi }E=G_{\chi }\left(E^{sc}+\mathbf{E}\right)=G_{\chi }\left(k^2+{\mathrm{dd}}^{\ast }\right)A+G_{\chi }\mathbf{E}$$

$$G_{\chi }\mathbf{E}=A-G_{\chi }{E}^{sc}=\left(I-G_{\chi }\left(k^2+{\mathrm{dd}}^{\ast }\right)\right)A_0.$$

$$G_{\chi }\mathbf{E}=\left(I-G_{\chi }{k}^2\right)C+\left(I-G_{\chi }\left(k^2+\Delta \right)\right)\mathrm{d}g$$

$$\chi \mathbf{E}$$

$$=\left(\left(1-\chi \right)k^2+\Delta \right)C+\left(1-\chi \right)\left(k^2+\Delta \right)\mathrm{d}g$$

$$\left(\Delta +k^2\right)G={\delta }_0$$

$$\theta \mathbf{E}=\left(k^2+\left(1+\theta \right)\Delta \right)C+\left(k^2+\Delta \right)\mathrm{d}g,$$

$$\theta \left(\mathbf{E}-\Delta C\right)=\left(k^2+\Delta \right)C+\left(k^2+\Delta \right)\mathrm{d}g.$$

$$\mathrm{d}\left(\theta \left(\mathbf{E}-\Delta C\right)\right)\equiv \mathrm{d}\theta \wedge \left(\mathbf{E}-\Delta C\right)+\theta \mathrm{d}\left(\mathbf{E}-\Delta C\right)=\left(k^2+\Delta \right)\mathrm{d}C.$$

$$\left(\mathbf{E}-\Delta C\right)\wedge \theta \mathrm{d}\left(\mathbf{E}-\Delta C\right)=\left(\mathbf{E}-\Delta C\right)\wedge \left(k^2+\Delta \right)\mathrm{d}C.$$

$$\Delta C$$

$$\mathrm{d}\Delta C$$

$$\theta =\frac{\left(\mathbf{E}-\Delta C\right)\wedge \left(k^2+\Delta \right)\mathrm{d}C}{\left(\mathbf{E}-\Delta C\right)\wedge \mathrm{d}\left(\mathbf{E}-\Delta C\right)}=\frac{\left(\mathbf{E}-\ast \mathrm{d}B\right)\wedge \left(k^2+\Delta \right)B^{\ast }}{\left(\mathbf{E}-\ast \mathrm{d}B\right)\wedge \left(\mathrm{d}\mathbf{E}-\Delta B^{\ast }\right)}.$$

$$\chi =\frac{1}{\theta +1}-1=-\frac{\left(\mathbf{E}-\ast \mathrm{d}B\right)\wedge \left(k^2+\Delta \right)B^{\ast }}{\left(\mathbf{E}-\ast \mathrm{d}B\right)\wedge \left(\mathrm{d}\mathbf{E}+k^2{B}^{\ast }\right)}$$

$$B^{+}=\left(B_x+i{B}_y\right)/2$$

Full text

A non-iterative method for electrical property tomography based on a simple formula

Victor P. Palamodov

Tel Aviv University

**Abstract **A non-iterative method of reconstruction is proposed from data of MRI system and of a harmonic electro-magnetic field at Larmor frequency. The method is based on the exact analytic formula for the contrast source function. A geometric method for acquisition of the full inductive field is discussed.

**Key words: **contrast source function, Helmholtz operator, transmit magnetic field, acquisition geometry

MSC 2010: 35Q61 65Z99

1 Introduction

Electrical properties tomography (EPT) is a noninvasive reconstruction technique for retrieving electrical properties (conductivity and permittivity) of biological tissue from magnetic fields generated by radio frequency coils in a magnetic resonance imaging scanner. The electric properties of tissue are of great interest, since these properties can be used to aid the discrimination of cancerous tissue from benign tissue and characterize various kinds of pathological tissue of human body [4], [8]. The main benefits of EPT over other reconstruction modalities is that it uses the Larmor frequency fields of an MRI system, which can penetrate biological tissues and it does not make use of surface electrode mounting, current injection or additional hardware.

The method called now EPT was proposed by Haacke et al. [1]. The first application of this method was done by Wen [2] who used an approximation for the contrast source function $\gamma/\gamma_{0}$ (CSF). Song and Seo [6] reduced reconstruction of admittivity $\gamma$ to solution of an elliptic equation under assumption $\partial\gamma/\partial z=0$ in terms of $\mathit{B}^{+}$ field. See also [9] for an explicit construction. Ammari *et al [7] *studied the problem in terms of the quasilinear system of differential equations. In [10] a survey of local methods of reconstructions is given. C van den Berg et al [11] developed the iterative method of determination of the CSF of an object making use of a global integral approach. Arduino *et al [12] applied the iterative conjugate gradient method (1000 steps). We describe here a non-iterative algorithm *for retrieving CSF which is based on the exact global formula. This is the first algorithm of this kind so far.

2 Basic equations

We assume that the conductivity and permittivity are isotropic at the angular Larmor frequency $\omega$ and ignore spatial permeability variations $\mu=\mu_{0}$, since they are considered small for biological tissue. Let $\Omega$ be an object domain in the MRI scanner that is disjoint of the antenna (coil) generating the RF wave.  The total time harmonic electromagnetic field in $\mathbb{R}^{3}$ can be written in the form \left\{E\exp\left(i\omega t\right),H\exp\left(i\omega t\right)\right\}\where

[TABLE]

$\left\{\mathbf{E=E}^{inc},\mathbf{H=H}^{inc}\right\}$ is the incident and \left\{E^{sc},H^{sc}\right\}\the scattered fields  due to the presence of the object. Following to van der Berg et al [8] we use the integral representations

[TABLE]

[TABLE]

for the scattered fields where $\gamma_{0}=i\omega\varepsilon_{0}$,$\ \varepsilon_{0}$ is the permittivity in vacuum and

[TABLE]

is the vector potential and

[TABLE]

is the scalar Green’s function for the Helmholtz operator $\Delta+k^{2}$ in 3D background medium. The contrast source function is defined as

[TABLE]

where admittivity $\gamma=\sigma+i\varepsilon\omega$ is the per meter admittance. Conductivity $\sigma$ and permittivity $\varepsilon$ are profiles of the object, \Omega\ is a domain that contains the support of $\chi.$

3 Evaluation of the contrast function

We write the electric and magnetic fields as first order differential forms

[TABLE]

Hodge star operator $\ast$ acts linearly on differential forms $\alpha,\beta$ on $\mathbb{R}^{3}$ by the rule ([13], p.15)

[TABLE]

where \alpha,\beta\are arbitrary differential forms in \mathbb{R}^{3}\of equal  degree and \left\langle,\right\rangle\means the natural scalar product of the forms. The dual differential is defined by $\mathrm{d}^{\ast}=\ast\mathrm{d}\ast$ and Laplace operator $\Delta=\mathrm{d}^{\ast}\mathrm{d}+\mathrm{dd}^{\ast}$ commutes with $\mathrm{d}=\nabla$ and \mathrm{d}^{\ast}.\

Theorem 1

Let $\Omega$ be a domain in $\mathbb{R}^{3}$ disjoint of the RF antenna. If the fields $\mathbf{E}$ and $B$ are known then the contrast source function can be found in $\Omega$ by

[TABLE]

if the denominator does not vanish.

Corollary 2

We have

[TABLE]

Note that $\left\langle\mathrm{d}\mathbf{E,E}\right\rangle=0$ since $\mathrm{d}\mathbf{E}=-i\omega c^{-1}\mu\mathbf{H}^{\ast}$ and \mathbf{H}^{\ast}\wedge\mathbf{E}=\left\langle\mathbf{H,E}\right\rangle=0.\

**Remark. **Equation

[TABLE]

is mentioned* *by Seo [6] and attributed to Nachman et al [3].

Proof of Theorem 1**. **

Lemma 3

Let $\Omega$ be a compact set in $\mathbb{R}^{n}$ that can be contracted to a point in itself. For any k\geq 1,\and arbitrary closed differential $k$-form $a$ on $\Omega$, there exists a $k-1$-form $b$ on \Omega\such that\ \mathrm{d}b=a\on $\Omega.$

For a proof see ”Converse of the Poincaré lemma” [13], p.29.111We may assume that coefficients of $a$ and of b\are distributions since $b$ will not appear in the final formula..

Note that equation (5) does not depend on \Omega.\Therefore we may assume that$\ \Omega$ is a compact set with smooth boundary in \mathbb{R}^{3}\that fulfils the condition of Lemma 3.

Lemma 4

The system of equations for a 1-form $C$

[TABLE]

has a solution defined on \Omega.\There exists a function $g$ on $\Omega$ such that

[TABLE]

Proof. Form $B^{\ast}$ has bounded coefficients and fulfils the Gauss law $\mathrm{d}B^{\ast}=0$ on \Omega.\By Lemma 3 there exists 1-form $C_{0}$ on \Omega\that satisfies $\mathrm{d}C_{0}=B^{\mathbf{\ast}}$. Consider Dirichlet problem

[TABLE]

for a function $f$ on \Omega.\To solve it we define the function $h$ on $\mathbb{R}^{3}$ that is equal to -\mathrm{d}^{\ast}C_{0}\on $\Omega$ and $h=0$ on the complement to $\Omega.$ Set $f=G_{0}\ast h$ where $G_{0}$ is the kernel (4) with k=0.\Differential form $C=C_{0}+\mathrm{d}f$ fulfils (8) since

[TABLE]

Equation (2) implies

[TABLE]

where $A_{0}=i\omega c^{-2}A$ and $A$ is the vector potential (3). By (8) we have $\mathrm{d}\left(A_{0}-C\right)=0.$ Again, by Lemma 3 there exists a function $g$ on $\Omega$ such that $A_{0}-C=\mathrm{d}g$ and (9) follows.$\ \blacktriangleright$

By (3) and (1) we have

[TABLE]

since \nabla\nabla\cdot\ =\mathrm{dd}^{\ast}\where$\ G_{\chi}u=G\ast\left(\chi u\right)$. This yields

[TABLE]

By (9) we have

[TABLE]

since \nabla\nabla\cdot\mathrm{d}g=\mathrm{dd}^{\ast}\mathrm{d}g=\Delta\mathrm{d}g.\Apply Helmholtz operator to both parts and get

[TABLE]

since

[TABLE]

where $\delta_{0}$ is the delta-function on $\mathbb{R}^{3}$. Dividing by $1-\chi$ we obtain

[TABLE]

where \theta\doteqdot\chi/\left(1-\chi\right)\hence

[TABLE]

Differential operators $\mathrm{d}$ and $k^{2}+\Delta$ commute, therefore

[TABLE]

Multiplying by form $\mathbf{E}-\Delta C$ we kill the term with $\mathrm{d}\theta:$

[TABLE]

By (9)\ \mathrm{d}C=\mathrm{d}A_{0}=B^{\ast},\hence

[TABLE]

which yields

[TABLE]

Finally

[TABLE]

and (5) follows.$\ \blacktriangleright$

4 Acquisition geometries

Any acquisition geometry should be rich enough to guarantee reconstruction of field $B$ on $\Omega.$ Let $x,y,z$ be an euclidean positively oriented coordinate system on the physical space such that the static magnetic field has the form $B_{0}=\left(0,0,b\right),$ $b>0$ in this system. The field

[TABLE]

is called positively rotating part of B\or transmit field. Several methods (sequences) are known that provide positively rotating part of a magnetic field for ex. [5]. According to [4] the ”negatively rotating” part $\mathit{B}^{-}=\left(B_{x}-iB_{y}\right)/2$ can not be determined in this way.** **The asymmetric role of coordinates $x$ and $y$ in (10) is defined with respect to $z$ coordinate which means that the coordinate system $x,y,z$ is positively oriented. Note that the orientation is indispensable feature of the Maxwell system which is illustrated by Maxwell’s right hand rule. It follows that such measurements of the transmit field can be made for any positively oriented coordinate system that is in the system obtained by rotation of the system $x,y,z.$

**1. **The simple acquisition geometry for magnetic field $B$ is to fix the body on the bed and rotate both around the central axis $y$ that is parallel to the bed. Let $x,y,z$ be the laboratory euclidean positively oriented (left-handed) system of coordinates,  and let $\xi,\eta,\zeta$ be the positively oriented system of coordinates that are constant on the bed and the body when rotating that is

[TABLE]

where $\varphi$ is the rotation angle see Fig.1.

Write magnetic field B\left(\varphi\right)=\beta\by means of both coordinate systems.

[TABLE]

and note that the body components of $\beta_{\xi},\beta_{\eta},\beta_{\zeta}$ of this field do not depend on \varphi.\Keeping in view that

[TABLE]

we obtain

[TABLE]

Suppose that measurements of the transmit field  are available for \varphi=-\psi,0,\psi\and consider system of linear equations

[TABLE]

The solution exists and is unique for any $\psi\neq 0,\pi:$

[TABLE]

hence magnetic form $\beta$ is uniquely reconstructed by (11). Now we can apply (5) to the fields $B=B\left(0\right)=\beta$ and $\mathbf{E.}$

**2. **Another geometry for scanning a body is as follows: the bed with the body rotating in a tilted plane around the center $O$ in $B_{0}$ field. Let $\alpha,\ 0<\alpha<\pi/2$ be the constant tilting angle and e_{3}=\left(0,-\sin\alpha,\cos\alpha\right).\For any $\varphi,\ 0\leq\varphi<2\pi,$the vectors

[TABLE]

and e_{3}\form the orthogonal frame since $e_{1}$ and $e_{2}$ belong to the plane P\through the origin orthogonal to $e_{3}$. This frame is a positively oriented. Functions

[TABLE]

are euclidean coordinates that are constant on $P$. To express the transmit field (10) by means of the body coordinates we write

[TABLE]

and substitute in (11). This gives

[TABLE]

and

[TABLE]

Evaluating the transmit field for $\varphi=\varphi_{1},\varphi_{2},\varphi_{3}$ we obtain the system like (12) for unknown functions $\beta_{\xi},\beta_{\eta},\beta_{\zeta}$ that do not depend on $\varphi.$The determinant of this system equals

[TABLE]

where $\varphi_{ij}=\left(\varphi_{i}-\varphi_{j}\right)/2.$ It does vanish if $\alpha\neq 0,\pi$ and \varphi_{ij}\neq 0.\It follows that for arbitrary $\alpha\neq 0,\pi$, arbitrary different angles \varphi_{i},\ ,\functions $\beta_{\xi},\beta_{\eta},\beta_{\zeta}$ are uniquely determined form data of fields$\ \mathit{B}^{+}\left(\varphi_{i}\right).$

Algorithm for computation of Contrast source function

1. Make three measurements of the transmit field according to geometry 1 or 2.

2. Calculate  the field $\beta$ in the coordinate system of the body. and set $B=B\left(0\right)=\beta.$

3. Apply (5) to the fields $\beta$ and $\mathbf{E}$ and calculate the quotient as the function \chi\of the body coordinates\ \xi,\eta,\zeta.\