This code is intended primarily for reproducibility of results of J. Darde, N. Hyvönen, T. Kuutela and T. Valkonen - Electrodeless Electrode Model for Electrical Impedance Tomography.

The code has been tested with Matlab 2020b on Linux Fedora 33 and Matlab 2020b on Windows 10 Professional.

DEPENDENCIES

• C++ compiler compatible with Matlab. Tested with GCC 10.2.1 on Linux Fedora 33 and MSVC 2015 on Windows 10.
• MTIMESX - Fast matrix multiply with multidimensional support
  • Mex compilation of MTIMESX may require some tweaking depending on your platform. Especially, on Linux you must include a BLAS library to the compilation command. For example, on Fedora 33 the library can be compiled by calling in Matlab mex('-DDEFINEUNIX', '<path_to>/mtimesx.c', '/usr/lib64/libblas.so.3.9.0').
  • See the comments section at mathworks.com for discussion and further instructions.
• Triangle - A Two-Dimensional Quality Mesh Generator and Delaunay Triangulator (Modified)
  • See the original Triangle by J. R. Shewchuk
  • Triangle can be compiled by simply running make. Then, the resulting executable triangle must be copied to a location in your $PATH.
  • On Linux: Triangle must be installed so that it can be executed with the shell command triangle. For example, put the compiled executable in /usr/bin/triangle.
  • On Windows 10: The command bash -c "triangle" must be executable. For example, use Windows Subsystem for Linux and install Triangle as you would on Linux.
• Torso tank dataset
  • Extract the dataset as a folder under the same path as this repository.
• pdfcrop and convert (ImageMagick) are used if they are available on Linux.

INSTRUCTIONS

To reproduce our figures and tables, run generate_figures.m. Several figure windows are left open once the script is finished, and some image files are saved in the pics folder. The resulting figures should be equivalent to the figures in our publication. However, there are slight differences in floating point computations between MSVC and GCC compiled C++ components. Furthermore, minor deviations are expected as reproducible plotting on different hardware is remarkably difficult with Matlab.

After running generate_figures.m, you may also run useful_prints.m which prints Latex arrays for the measures used in the publication.

The Algorithm 1 in our article is mostly implemented in paper/EIT_minimization.m and paper/Generate_data_*_tank.m. These are largely undocumented beyond the descriptions in our article, but some notes are included.

ACKNOWLEDGEMENTS

Significant parts of the finite element system matrix construction and the mesh structure handling are based on work by Antti Hannukainen (Aalto University).

NOTES & KNOWN WARTS

This is research code and therefore provided AS IS for verification of our numerical results. There are several messy parts which will not be fixed. The computations use rather significant global state as this has been found useful for debugging and injecting modifications during experimentation. However, it certainly doesn’t ease the understanding of the code. Also note that this repository is a pruned version of a larger, unreleased, code base. Therefore there may be references to parts which are not available.

The following notes may help in understanding what is going on:

• Variables named “fm” refer to forward model. The FEM for electrical conduction within the domain is implemented directly in these objects, while the contact conduction/current injection is in the member fm.electrode_model.
  • The zero-mean condition for current injections is implemented using “canonical basis” in the electrode models. FEM computations are only done using a (fixed) zero-mean current-basis. The actual forward solutions are then computed using projections as the model is linear. Same applies to derivatives.
• The boundary conductivity is always implemented as a function, including for the standard CEM. This allows us to define the contacts at arbitrary positions on the boundary, instead of per edge. See Functional_electrodes.m.
• Files named Peak* correspond to the Parametric Hat (PH) contact conductivity parametrization in our article, while P1CEM* corresponds to Piecewise Linear (PL) parametrizations.
• Resistance_CEM_electrode_parametrization is the standard CEM parametrization used as reference. Note that the computations are actually computed in terms of contact resistivity, not contact resistance in this case!
  • Resistivity is preferred over conductivity for the numerical stability of the total Jacobian.
• The printed finalerrs-values can be interpreted as follows:
  1. Domain mean (3D) conductivity. For empty domain reconstructions this is the conductivity of water. For the domains with inclusions, this is mean-of-nodewise-conductivity, which is biased because of uneven density of elements in the mesh.
  2. l2-Residual, i.e. norm(U_r - U_m), where U_r are the voltages of the forward model at the end of reconstruction and U_m are the real life measured voltages. Note that U_m are in fact normalized to correspond to current injection of 1 mA, which can be done due to linearity of the model. See also discussion in the beginning of the numerical experiments section in our article.
  3. Domain prior cost, i.e. squares of the second column in Table 4.4. in our article.
  4. Electrode prior cost, i.e. squares of the third column in Table 4.4.
  5. The remaining values are used as prior parameters. These were used while searching for reasonable values. Units of these values have been adjusted for the article. Furthermore, the tank height is accounted for in the values of the article.
• Similar inconsistency between conductivity and conductance can be found in some sections of the code.
• The printed log_contact_conductances_mean and log_contact_conductances_sd correspond to the Table 4.3.
• The printed midpoint_errors correspond to the Table 4.5.
• The ring artifact in the reconstruction of inclusions in the 0 mm extension PL case seems to be caused by the domain conductivity regularization and its Gaussian smoothness prior.
• See the notes in torso tank dataset regarding the precission of the domain boundary and electrode positions. These are known to be imperfect.
  • In particular, we are aware that the CEM produces terrible results on torso6 dataset, but surprisingly good results if the electrodes are extended by even just a single edge element.
  • We believe this is caused by the imperfection of electrode positions.
• In our paper we claimed that the results generalize to the vast majority of our dataset. In fact, the above torso6 is the only measurement in which the standard CEM seemingly outperforms PL and/or PH relatively often (w.r.t. different electrode extensions). We also tested our method with the Finnish Inverse Problems Society Electrical Impedance Tomography (EIT) Dataset and the results were qualitatively equivalent to our dataset.
• The folder paper contains most of the top-level scripts used to generate the figures and tables in our article. Read “hard” as the domain with inclusions and “empty” as the empty tank test cases. “constant” scripts are used for the empty tank and “gauss” for the hard tank.
• We are aware that the backtracking implementation could be improved by e.g. actually searching for Wolfe-condition filling step length. Especially CEM seems to have some trouble converging.
• The stopping criteria of step length being less than 10^-4 is not really justified beyond seemingly being good enough. Decreasing this value, and allowing more iterations, improves the final residual of CEM. However, the domain reconstructions remain qualitatively the same.

LICENSE

The Matlab code (including C++ components) is licensed under CC-BY-4.0.

Copyright © 2020-2021 Topi Kuutela, Antti Hannukainen