Introduction j 0 r 8 13 16 22 B 0 1 B V V V 0 B 0 Figure 1. Schematic of a possible coil system for MIT with 16 excitation coils and 32 receiver coils. 8 22 32 2 14 16 24 26 7 23 In contrast to electrical impedance tomography (EIT) MIT avoids the ill-defined electrode-skin interface due to its inherently contact-less operation. 1 1 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\bf y}=\varvec{\Psi}(\varvec{{\bf \kappa}})$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\kappa}$$\end{document} y y M a b a b 2 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\bf\varvec{\kappa}=\varvec{\Psi}^{-1}({\bf y}) $$\end{document} 21 3 13 33 9 15 17 27 29 Methods N i \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\kappa}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\kappa}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\kappa}$$\end{document} 3 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\bf y}=\varvec{\Psi}(\varvec{\kappa}) $$\end{document} y m \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\kappa}^{\ast}$$\end{document} 4 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \varvec{\kappa}^{\ast}=\arg\mathop{\min}\limits_{\varvec{\kappa}}\left( (\varvec{\Psi}(\varvec{\kappa})-{\bf y}_{\bf m})^T(\varvec{\Psi}(\varvec{\kappa})-{\bf y}_{\bf m})+\lambda \varvec{\kappa}^T{\bf R}^T{\bf R}\varvec{\kappa}\right). $$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\kappa}^{\ast}$$\end{document} R \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\kappa}$$\end{document} p k k \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \varvec{\kappa}_{k+1}=\varvec{\kappa}_{k}+{\bf p}_{k} $$\end{document} 5 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\bf p}_k =({\bf G}^T_k {\bf G}_k +\lambda {\bf R}^T{\bf R})^{-1}{\bf G}^T_k {\bf e}_k $$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\bf e}_k=(\varvec{\Psi}(\varvec{\kappa}_k)-{\bf y}_{\bf m})$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\bf G}_k=\frac{d\varvec{\Psi}_k}{d\varvec{\kappa}_k}$$\end{document} sensitivity matrix 4 6 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \Delta {\bf y}_{\bf m}={\bf Gp}=\frac{d\varvec{\Psi}}{d\varvec{\kappa}}{\bf p} $$\end{document} 7 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\bf p}=({\bf G}^T{\bf G}+\lambda {\bf R}^T{\bf R})^{-1}{\bf G}^T(\Delta {\bf y}_{\bf m}) $$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\bf p}=\Delta\varvec{\kappa}$$\end{document} y m Calculation of the Forward Solution and the Sensitivity Matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\Psi}(\varvec{\kappa})$$\end{document} 8 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \begin{array}{ll} \hbox{curl}{\bf H}={\bf J}& \\ \hbox{curl}{\bf E}=-j\omega{\bf B}&\hbox{in }\Omega \\ \hbox{div}{\bf B}=0&\\ {\bf B}=\mu {\bf H},{\bf J}=\kappa{\bf E},\kappa =\sigma+j\omega \varepsilon&\\ \end{array} $$\end{document} H B E J 12 , 17 A r A r A r B \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\bf G}=\frac{d\varvec{\Psi}}{d\varvec{\kappa}}$$\end{document} 31 20 dy d \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\kappa}$$\end{document} 9 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \frac{dy}{d{\varvec \kappa}}={\varvec I}_\phi \int\limits_{\Omega}{\bf L}_\phi{\bf L}_\psi {\bf d}{\varvec \Omega} $$\end{document} 10 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\bf L}_\phi =-\frac{\varvec{j}\omega {\bf A}_\phi +\nabla \varvec{V}_\phi }{{\bf I}_\phi}\quad \quad {\bf L}_\psi=\frac{\varvec{j}\omega {\bf A}_\psi +\nabla \varvec{V}_\psi}{{\bf I}_\Psi} $$\end{document} A ϕ A Ψ V ϕ V ψ I ϕ I ψ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${d}\varvec{\Psi}/{d}\varvec{\kappa}$$\end{document} 11 L ϕ L Ψ Regularization R T R I 10 1 R T R E pp T 5 5 G G U V T R T R VDV T D d i 11 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ d_i=\frac{\sigma_i}{\sqrt{c}}-\sigma_i^2 $$\end{document} i I c 19 12 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\bf p}={\bf V}_t\Sigma_t^{-1} {\bf U}_t \Delta y_m $$\end{document} t V U t R T R I R T R N N \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ N_{ij}=\left\{{\begin{array}{ll} {n}_{\rm n}&{\hbox{i}=\hbox{j}}\\ -\hbox{1}&\hbox{i,j}\quad \hbox{neighbours}\\ 0&\hbox{otherwise}\\ \end{array}}\right. $$\end{document} n n i N R T R VDV T t c c c c G T G R T R 10 Gp y m Modeling Setup x y z x y z 1 2 26 Figure 2. x y z x y z It represents a true 3-D-arrangement which delivers theoretically 512 independent measuring combinations, i.e. 512 data points for one image reconstruction. It is similar to our experimental system which employs 16 excitation sites in one plane and 14 planar gradiometers which are formed by connecting in counter-phase the coils in the upper and in the lower receiver plane. 6 The calculation of the complete sensitivity matrix required 48 forward solutions according to the Mortarelli-approach. Theoretical Limits of Resolution and Contrast/Noise Ratio 30 There is a fundamental limit for the resolution which depends on the amount of available information in the data. This information depends on the number of data points, i.e. the number of possible sensor–detector-combinations and on the degree of independence between these data points. In the case of EIT and MIT the number of independent data points is usually much lower than the number of voxels, so that the system is under-determined. Moreover the different data are correlated to a certain degree, so that the effective rank is comparatively low. In EIT, e.g. 16 electrodes provide 104 independent data points so that the information is no more than 104 ‘effective pixels’. Including some a priori-information in the form of the regularization terms leads to a defined ‘smearing’ of this information over the imaging plane and provides the typical diffuse images known from EIT. p * p 13 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \begin{aligned} &{\bf y}_{\rm m}={\bf Gp}^{\ast}\\ &{\bf p}={\bf Ay}_{\rm m}={\bf AGp}^{\ast}={\bf Mp}^{\ast} \end{aligned} $$\end{document} A A G T G R T R −1 G T A V t t −1 U t i M i 14 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\bf M}={\bf V}_t\Sigma_t^{-1}{\bf U}_t {\bf G}\quad \hbox{with}\ t=\hbox{rank} ({\bf G}^{T}{\bf G}) $$\end{document} p x 15 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox{CNR}({\bf x})=\frac{\Delta p({\bf x})}{\hbox{std}(np({\bf x}))} $$\end{document} p np p p p np x X y m np 16 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \hbox{Cov}({\bf np})={\bf AXA}^T $$\end{document} p 17 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \Delta \varvec{p}_i=\sum\limits_{\varvec{j}\in \varvec{P}} \varvec{M}_{ij} \varvec{p}_j^{\ast} $$\end{document} 18 V V V Results 0 3 1 25 Figure 3. Dependence of the theoretical normalized resolution on the noise level. Curves are depicted for noise-free data (TSVD with truncation level 512) and TSVD with truncation levels corresponding to a SNR of 44, 50 and 64 dB, respectively. TABLE 1. Tuning parameters for the regularization, chosen according to the Morozov-criterion. Regularization method TSVD Truncation level IM λ NM λ VU λ SNR = 64 dB 182 2.40E − 20 2.40E − 20 2.35E − 11 SNR = 50 dB 100 8.40E − 19 8.80E − 19 3.20E − 10 SNR = 44 dB 68 1.25E − 17 1.19E − 17 7.40E − 10 c x In analogy to EIT the PSF depends strongly on the location, showing the broadest distribution in the center of the cylinder. This is reflected by increasing resolution when moving from the center towards the periphery. Increasing the amount of regularization or decreasing the truncation level according to increasing measurement noise the PSF broadens and its center of gravity is shifted towards the borders of the cylinder. Moreover increasingly strong ringing in form of star-like patterns becomes observable close to the border (not shown explicitly in this paper). 4 Figure 4. Theoretical normalized resolution for the four regularization methods at a SNR of 50 dB. 5 5 5 Figure 5. Panel A: CNR for TSVD as a function of the normalized x-coordinate (relative to the cylinder radius) and in dependence on the noise level. Panel B: comparison of the four methods at a SNR of 50 dB. 6 7 1 Figure 6. Mean images of the Monte-Carlo study. Reconstructed Δ σ (transversal and saggittal section through the origin) for the spherical perturbations with four different regularization matrices and a SNR of 64 dB. Figure 7. Mean images of the Monte-Carlo study. Reconstructed Δ σ (transversal and saggittal section through the origin) for the spherical perturbations with four different regularization matrices and a SNR of 50 dB. In all cases the two perturbations can be recognized as diffuse bright disks. The dotted circles in the figures delineate the original position of the perturbing spheres. 6 7 2 Mean and CNR of the pixel values in the center of gravity of each reconstructed perturbation. These parameters quantify the correctness of the reconstructed values as well as their uncertainty. The theoretically expected values are listed for comparison. The center of gravity was chosen as evaluation point because the reconstructed perturbations deviate more or less from the spherical shape and show significant outward shift with increasing noise level. Radial outward shift of the spheres in the reconstructed image (fidelity of the location). This shift was determined by localizing the center of gravity for each spot within a wedge-shaped search region with a height of 2.6 times the sphere’s radius and excluding the outermost 2 mm as well as the innermost 40 mm in radial direction from the center. The restriction of the search region to this volume prevented spurious contributions from outliers and negative image values far away from the real perturbing regions. Also for this parameter we present theoretical values as expected from the PSF. TABLE 2. Summary of the performance indices defined in the text. SNR = 64 dB SNR = 50 dB SNR = 44 dB Method Performance index Sphere real Sphere theoret Background Sphere real Sphere theoret Background Sphere real Sphere theoret Background IM Mean central 0.021 0.016 0.00054 0.006 0.0043 0.0035 0.001 NM 0.022 0.016 0.00056 0.006 0.004 0.00051 0.0032 0.001 0.00035 VU 0.024 0.024 0.00055 0.01 0.0065 0.00043 0.0060 0.0036 0.00035 TSVD 0.033 0.025 0.00054 0.01 0.0074 0.00050 0.0110 0.002 0.00055 IM CNR 42.3 45.0 16.1 23.0 12 20.1 NM 41.1 44.0 16.5 22.5 10.7 19.3 VU 7.9 8.2 5.6 5.2 3.4 3.3 TSVD 23.9 25.8 13.7 18.4 9.1 13 IM Resolution 2.84 2.60 2.53 NM 2.81 2.58 2.50 VU 2.94 2.46 2.58 TSVD 2.58 2.42 2.24 IM Normalized outward shift* 0.10 0.10 0.21 0.25 0.33 0.35 NM 0.10 0.10 0.21 0.25 0.33 0.35 VU 0.10 0.06 0.11 0.15 0.21 0.15 TSVD 0.10 0.10 0.18 0.21 0.35 0.33 Performance measures for the comparison of the reconstruction methods. *Normalized to the cylinder radius. 2 2 2 0 8 Figure 8. Comparison of single-shot images for VU and IM at the three different noise levels. IM and NM perform nearly identically, also their optimal regularization parameters are almost identical. VU yields, in general, larger values in the perturbed regions but also a larger STD. Discussion 27 29 27 29 3 2 p 3 2 6 7 3 30 2 5 In the case of weak perturbation we can assume that the CNR depends approximately linearly on the conductivity difference Δσ. The dependence on the volume of the perturbation is, in general, more complicated because the PSF depends on the location and is therefore not constant throughout the whole perturbation. Only in the case of small spatial extension of the perturbation an approximately linear dependence on the volume can be assumed. The low number of significant singular values even at comparatively low noise (64 dB SNR) suggests that, similar as in EIT, a significant amount of sensor combinations does not provide enough independent information. Intuitively one would expect this finding because there exist pairs of excitation/receiving coils which nearly fulfill the reciprocity condition and hence reduce the amount of useful combinations to about half of the number of possible combinations, i.e. to 256 in our case. Further investigations should determine the maximum ‘useful’ number of sensors in one plane, i.e. that number beyond which additional sensors do not increase the resolution significantly. Adding more sensors off-plane may add more 3-D-information and hence still provide improvement. This possibility should be studied in further research. When comparing the regularization schemes after application of the Morozov-criterion, the IM and the NM approach yield the smoothest visual appearance and the highest CNR. However, they also tend to displace the perturbations towards the border of the tank. The best localization is obtained with VU, probably because the imposed variance counteracts somewhat the lower sensitivity in the center of the object. However, VU yields also the lowest CNR, i.e. the less homogeneous images and more pronounced ghosts. The failure of TSVD at a SNR of 44 dB was not expected theoretically, although, in general, it produces the poorest theoretical resolution. In terms of separability of the two perturbations VU performs best, especially when also taking into account the correct localization. In neither case, however, the single-step solution provides the correct values for Δσ. Even at a SNR of 64 dB the reconstructed differences are too low by a factor of at least 5, thus demonstrating that the method yields the correct search direction but not the correct step size. The highest mean voxel values are provided by VU and TSVD, the drop with the noise levels being lowest. However, on the other hand VU yields the highest standard deviations. Moreover VU tends to produce more pronounced ‘ghost objects’ in the homogeneous region than IM and NM. As expected from the PSF TSVD tends to produce ‘star-artifacts’ at the cylinder border, i.e. a periodic pattern with 16 peaks close to the centers of the receiving coils. This artifact gets worse at increasing noise level. First experiments with smaller models and at least 10 iterations with an iterative solver show that the solution converges towards the correct voxel values. Nevertheless the single-step method may be completely justified in cases where only qualitative changes are sought for or where proportions are to be reconstructed, e.g. in frequency differential spectroscopic imaging. Therefore the area of applicability of a single-step approach has to be analyzed carefully in future work. 28 , 26 18 28 Another open question is the influence of the mesh quality on the reconstruction results. We used a comparatively coarse non-uniform grid for the reconstruction. Therefore non-negligible numerical errors are to be expected which may explain the discrepancies between theoretically expected and the reconstructed values for CNR and radial displacement. Also the apparently somewhat worse spatial resolution in the reconstructed images than theoretically expected may be due to such numerical problems. The influence of the mesh and the optimization of mesh quality should be a major issue for further developments. The 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. A multi-frequency approach is expected to increase significantly the available information and thus the quality of the images. Possible applications may then in fact be the same as for EIT (lung function monitoring, lung edema monitoring) and hydration monitoring in the brain.