eitcom / pyeit Goto Github PK

@liubenyuan Since issue #22, I was wondering how to implement GPU support in the module. Could it be a thing to have a second module, completely separated? Like tensorflow and tensorflow-gpu? It may allow us to have gpu acceleration and have gpu specific optimizations in the program.
This would need 2 separated pipelines, but it could be of a great help to a lot of people. And the limited memory allowed on a single GPU makes it a good start to take a look at huge models support, since we would need some kind of slicing/dataflow method we could as well implement in the classic module.

Hi,
totay, i install pyeit by "conda install -c conda-forge pyeit",
but when i run "python eit_dynamic_bp.py" for test, i have some problem as follow,
i dont know why where is wrong, hope your answer.
Traceback (most recent call last):
File "E:\pyEIT\examples\fem_forward2d.py", line 38, in
protocol_obj = protocol.create(n_el, dist_exc=7, step_meas=1, parser_meas="std")
File "E:\pyEIT2\lib\site-packages\pyeit\eit\protocol.py", line 188, in create
return PyEITProtocol(ex_mat, meas_mat)
TypeError: init() takes 1 positional argument but 3 were given

BR,
xlzhang.

Hi @liubenyuan,

thanks for this great project.

I was wondering how you could work with voltage measurements taken from the real world. If you imagine that you take a voltage measurement and get a pandas time series (like with the load_csv(fstr, resample=None) method from pyeit/io/icp.py), how would you proceed in the code to generate an EIT image?

Thank you and kind regards,
Pascal

I am planning on implementing additional performance measures as described by Adler et. al. in Evaluation of EIT system performance, for the purposes of analysing measured data from EIT hardware.

  Hello everybody. I am interested in 3D reconstruction. I tried to run the example "eit_static_GN_3D"  but i got these errors:

line 24, in
mesh_obj = mesh.create(n_el, h0=0.2, bbox=bbox)
line 375, in create
p, t = build(fd, fh, pfix=p_fix, bbox=bbox, h0=h0)
line 402, in build
dm = DISTMESH(
line 102, in init
p = p[fd(p) < self.geps] # pylint: disable=E1136
line 37, in circle
return dist(pts - pc) - r
ValueError: operands could not be broadcast together with shapes (1080,3) (2,)

Why I got this error and how can i fix this? Thank you.

Originally posted by @mlucciardi in #82

Complete electrode model (CEM) could better modeling the current distributions on the electrode and the contact impedances, CEM should be incorporate into pyeit.

• block matrix calculation of the stiffness matrix.
• how to modeling the CEM into the PyEITMesh class? PyEITProtocol class could remain the same.
• Forward computing of CEM and visualization. Compare the results to the simple electrode model.

I was working with a previous version of this package, and I built myself a dataset as well, and for some reasons the new version is now unable to reconstruct properly from my saved data. The reconstruction was changed in the last 3 months maximum, and the version I was working with was still not available on pypi.

Any idea on how to fix this without reconstructing my entire dataset?

I am a student researching Electrical Impedance Tomography (EIT). While using the gn function, I encountered an inconsistency. In the paper "pyEIT: A Python-based Framework for Electrical Impedance Tomography," the update equation for x0 is stated as x0 + (J^TJ + lambR)^(-1)J^Tr0. However, in the gn function, it seems to be implemented as x0 - (J^TJ + lambR)^(-1)J^Tr0. After comparing the representations of each symbol, I found no discrepancies, yet there is a difference in the positive/negative sign. Could you kindly explain this inconsistency?
I would greatly appreciate any insights or clarification you could provide regarding this matter. Thank you for your time and assistance.

First of all thanks for the project!

I collect real world data using a 16 electrode belt wrapped around the pelvic (bladder monitoring). Measurements are taken in opposition schema resulting in a DataFrame with 448 columns (224 for each Magnitude and Phase) like this:

E01_E09_E02_E10_Magnitude | E01_E09_E02_E10_Phase | E01_E09_E03_E11_Magnitude | E01_E09_E03_E11_Phase | -- | E16_E08_E14_E06_Magnitude | E16_E08_E14_E06_Phase | E16_E08_E15_E07_Magnitude | E16_E08_E15_E07_Phase

The electrodes are not spaced equally around the pelvic, but with vacancy above the pubic bone due to hardware design:

This is the code I am using for reconstruction:

config= {
    "n_el": 16,
    "h0": 0.07,
    "dist_exc": 8,
    "step_meas": 8,
    "parser_meas": "std",
    "inverse_solver": "JAC",
    "p": 0.5,
    "lamb": 0.01,
    "method": "kotre",
    "perm": 1,
    "jac_normalized": True,
    "normalize": True
}

mesh_obj = mesh.create(n_el=config['n_el'], h0=config['h0'], fd=mesh.shape.thorax)
mesh_obj.el_pos = np.array([517, 478, 5, 302, 219, 110, 42, 17, 15, 40, 109, 274, 11, 455, 515, 528])  # Custom arrangement
protocol_obj = protocol.create(config['n_el'], dist_exc=config['dist_exc'], step_meas=config['step_meas'], parser_meas=config['parser_meas'])

eit_solver = JAC(mesh_obj, protocol_obj)
eit_solver.setup(
    p=config['p'],
    lamb=config['lamb'],
    method=config['method'],
    perm=config['perm'],
    jac_normalized=config['jac_normalized']
    )
    
ds = 224.0 * eit_solver.solve(v0, v1, normalize=config['normalize'])  # v0, v1 contain the complex representation

Reconstruction runs through, but the results are not as expected. Even when validating the hardware on an agar-agar phantom, I found the algorithm struggling with the electrode position (results are better when assuming a equally spaced electrode arrangement, but still not as expected):

Am I missing something, or is the combination of opposition schema with custom electrode arrangement just not suitable?

The latest dev branch is working on a redesign (not backward compatible) version of EITForward class. This new design isolate,

1. mesh (including pts: coordinates of nodes, tri: the connection of simplices, perm: permittivity, el_pos: electrode positions, ref: reference node),
2. and the measurement protocol (ex_mat: excitation matrix (Neumann BC only, yet), step: voltage difference electrode, parser: measurement order).

For example,

# forward simulation
fwd = EITForward(mesh, protocol)
v0 = fwd.solve_eit()
v1 = fwd.solve_eit(perm_new, init=True)
# EIT imaging
algo = JAC(mesh, protocol)
ds = algo.solve(v1, v0, normalize=True)

The mesh and protocol can be a dict, or a dataset(#44), if the code block for building these two data structures (mesh, protocol) are intuitive to use, i.e., see the construction code of mesh and protocol in examples.

This new version (maybe 2.0) of pyeit will be release based on these two milestones:

• simple, clean and easy to use API for EIT forward and inverse simulations,
• large scale Forward simulation using a fine-grained Head/Thorax phantoms.

why i don't see the value of injected current in your code. Is that important?

When the value of perm defaults to all 1, either through fwd.compute_jac() or fwd.solve_eit(), the solved v0 is clearly wrong.
As an example, when the number of electrodes is 16, using the method of adjacent excitation-adjacent measurement, the calculated v0 is not all U-shaped, but appears to plateau at some locations.

This is the v0 curve run with pyEIT, but the red circle in the figure below is clearly not correct.

This is not a coincidence! I troubleshot that the problem may be in the solve_vectorized function in pyeit/fem.py, where solving a system of linear equations using the sparse matrices kg and b appears rather obviously off at i=11.
Could you please troubleshoot and fix the problem?

After update from ver. 1.1.6 to 1.2.2 I can't run example codes.

i.e.: # eit_dynamic_bp.py

returns:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
~\AppData\Local\Temp/ipykernel_39680/3106232781.py in <module>
      8 import matplotlib.pyplot as plt
      9 import numpy as np
---> 10 import pyeit.eit.bp as bp
     11 import pyeit.eit.protocol as protocol
     12 import pyeit.mesh as mesh

C:\ProgramData\Anaconda3\lib\site-packages\pyeit\eit\__init__.py in <module>
     12   - interp2d: Spatial interpolation for EIT
     13 """
---> 14 from .bp import BP
     15 from .jac import JAC
     16 from .svd import SVD

C:\ProgramData\Anaconda3\lib\site-packages\pyeit\eit\bp.py in <module>
      8 from typing import Union
      9 import numpy as np
---> 10 from .base import EitBase
     11 
     12 

C:\ProgramData\Anaconda3\lib\site-packages\pyeit\eit\base.py in <module>
     11 from abc import ABC, abstractmethod
     12 import numpy as np
---> 13 from pyeit.eit.protocol import PyEITProtocol
     14 from pyeit.mesh import PyEITMesh
     15 

C:\ProgramData\Anaconda3\lib\site-packages\pyeit\eit\protocol.py in <module>
    143 def create(
    144     n_el: int = 16,
--> 145     dist_exc: Union[int, list[int]] = 1,
    146     step_meas: int = 1,
    147     parser_meas: Union[str, list[str]] = "std",

TypeError: 'type' object is not subscriptable

[update 25/06/2022]:

I just checked all versions one by one, and it looks like version change from 1.1.9 to 1.2.0 broke something.
Codes run until version 1.1.9 with no issues, starting from 1.2.0 one get error:

TypeError                                 Traceback (most recent call last)
~\AppData\Local\Temp/ipykernel_76780/2754546162.py in <module>
     16 """ 0. construct mesh """
     17 # Mesh shape is specified with fd parameter in the instantiation, e.g : fd=thorax , Default :fd=circle
---> 18 mesh_obj, el_pos = mesh.create(16, h0=0.1, fd=thorax)
     19 
     20 # extract node, element, alpha

TypeError: cannot unpack non-iterable PyEITMesh object

I'm trying to recreate Nir Grossman's recent Cell paper simulation using pyEIT. I've got a nice static model,
with two current sources using pyEIT, but am hoping to recreate the quasi-electrostatic model. This likely could be applicable to MFEIT also.

I realise this isn't strictly needed by pyEIT as the static model is generally considered 'sufficient', but would you have have any recommendations with regards to packages that could do quasi static modelling(ideally that even integrate with pyEIT!). What are your thoughts on incorporating time effects into EIT reconstruction?

I would like to model a larger array of varied electrode configurations, from the standard tetrapolar adjacent method, opposition method, to bipolar and other user defined electrode configurations that don't always have the same step size. Currently functions seem limited to tetrapolar adjacent/opposition only. Is there a way to allow a user defined file(example)
e_config.txt

Above is an example of a 32 electrode configuration running in random order. Ideally it should acquire the same results as running the simulator but it doesn't as there are dependencies in other files on electrode ordering.

I tried to use dynamic imaging (conductivities imaging) for back projection
ds = 192 * eit.solve(v1, v0,normalize=true)

with v1 being a 780x280 matrix and v0 a 280 row vector.

But I have a problem with the dot operation (1002x208 cant dot with 780x208 matrix) Its a common mathematical issue. Is it possible to explain a little why this error is happening? I tried to transpose the 780x208 matrix for the dot operation but I receive a 1002x780 matrix as a response for the eit.solve operation and Im not sure how to reconstruct that response on an image.
Is there something that I'm not following correctly?

thanks for your help!

Discussed in https://github.com/liubenyuan/pyEIT/discussions/39

Do you have any plans to integrate a D-Bar algorithm implementation into pyEIT? I'd love to use it if it were.

https://epubs.siam.org/doi/abs/10.1137/060656930

Hi!

I am interested in pyEIT, and my team is currently working on reconstructing images. We have successfully built our hardware and achieved good results in image reconstruction. However, we are currently facing an issue. When we place the object near the electrodes, the reconstructed image looks good. However, when we position our object right in the middle of the tank, the reconstructed image appears to be shifted:

It consistently shifts to the right like this, regardless of the number of attempts we make. We have also compared the results of pyEIT with EIDORS (Matlab), and the image from EIDORS is completely centered. There should not be an issue with our hardware:

Regarding the code we ran, we modified the eit_dynamic_JAC.py in the examples folder, we assign our voltage matrices to v0 and v1:

v0 = np.loadtxt('example_data/ref_data.txt')
v1 = np.loadtxt('example_data/diff_middle_data.txt')

and also performed some tunning. Our txt files contain a single line of 208 values, representing the voltage measured from those electrodes.

Code:

# coding: utf-8
""" demo on dynamic eit using JAC method """
# Copyright (c) Benyuan Liu. All Rights Reserved.
# Distributed under the (new) BSD License. See LICENSE.txt for more info.
from __future__ import absolute_import, division, print_function
 
import time
import matplotlib.pyplot as plt
import numpy as np
import pyeit.eit.jac as jac
import pyeit.mesh as mesh
from pyeit.eit.fem import EITForward
from pyeit.eit.interp2d import sim2pts
from pyeit.mesh.shape import thorax
import pyeit.eit.protocol as protocol
from pyeit.mesh.wrapper import PyEITAnomaly_Circle
 
""" 0. build mesh """
n_el = 16  # nb of electrodes
use_customize_shape = False
if use_customize_shape:
    # Mesh shape is specified with fd parameter in the instantiation, e.g : fd=thorax
    mesh_obj = mesh.create(n_el, h0=0.065, fd=thorax)
else:
    mesh_obj = mesh.create(n_el, h0=0.065)
 
# extract node, element, alpha
pts = mesh_obj.node
tri = mesh_obj.element
x, y = pts[:, 0], pts[:, 1]
 
""" 1. problem setup """
# mesh_obj["alpha"] = np.random.rand(tri.shape[0]) * 200 + 100 # NOT USED
anomaly = PyEITAnomaly_Circle(center=[0.5, 0.5], r=0.1, perm=1000.0)
mesh_new = mesh.set_perm(mesh_obj, anomaly=anomaly)
 
""" 2. FEM simulation """
# setup EIT scan conditions
protocol_obj = protocol.create(n_el, dist_exc=1, step_meas=1, parser_meas="std")
 
# calculate simulated data
fwd = EITForward(mesh_obj, protocol_obj)
 
 
v0 = np.loadtxt('example_data/ref_data.txt')
v1 = np.loadtxt('example_data/diff_middle_data.txt')
 
#v0 = fwd.solve_eit()
#v1 = fwd.solve_eit(perm=mesh_new.perm)
 
 
""" 3. JAC solver """
# Note: if the jac and the real-problem are generated using the same mesh,
# then, data normalization in solve are not needed.
# However, when you generate jac from a known mesh, but in real-problem
# (mostly) the shape and the electrode positions are not exactly the same
# as in mesh generating the jac, then data must be normalized.
eit = jac.JAC(mesh_obj, protocol_obj)
eit.setup(p=0.5, lamb=0.8, method="kotre", perm=1000, jac_normalized=True)
ds = eit.solve(v1, v0, normalize=True)
ds_n = sim2pts(pts, tri, np.real(ds))
print('ds=\n', ds_n)
# plot ground truth
fig, ax = plt.subplots(constrained_layout=True)
 
 
 
# plot EIT reconstruction
im = ax.tripcolor(x, y, tri, ds_n, shading="flat", cmap=plt.cm.magma)
for i, e in enumerate(mesh_obj.el_pos):
    ax.annotate(str(i + 1), xy=(x[e], y[e]), color="r")
ax.set_aspect("equal")
 
fig.colorbar(im, ax=ax)
 
# plt.savefig('../doc/images/demo_jac.png', dpi=96)
plt.show()

Sample data:

In one-line form:
Reference:   2.25 1.18 0.79 0.69 0.53 0.43 0.49 0.52 0.50 0.67 0.85 1.15 2.12 4.80 1.36 0.95 0.66 0.51 0.51 0.55 0.47 0.59 0.65 0.67 1.17 1.12 2.59 1.38 0.89 0.66 0.66 0.59 0.50 0.51 0.58 0.59 0.86 1.38 2.10 5.02 1.60 0.99 0.79 0.69 0.46 0.50 0.47 0.47 0.59 0.84 1.12 0.96 1.55 0.99 0.70 0.63 0.46 0.39 0.38 0.31 0.39 0.50 0.62 1.11 1.96 4.51 1.50 0.95 0.61 0.62 0.48 0.45 0.47 0.57 0.65 0.75 1.57 2.08 4.77 1.57 0.95 0.61 0.55 0.36 0.39 0.43 0.45 0.62 0.87 1.33 1.11 2.12 1.13 0.86 0.54 0.44 0.41 0.43 0.45 0.48 0.70 0.84 1.07 2.12 4.70 1.53 0.96 0.59 0.56 0.50 0.50 0.50 0.59 0.72 0.56 1.35 1.34 0.31 1.66 1.03 0.69 0.68 0.57 0.53 0.63 0.60 0.60 1.00 1.41 2.03 4.65 1.32 0.86 0.65 0.55 0.49 0.47 0.43 0.40 0.62 0.62 0.94 1.37 0.61 1.45 1.07 0.74 0.59 0.57 0.47 0.41 0.52 0.63 0.75 1.31 1.13 1.94 1.26 0.89 0.71 0.60 0.46 0.40 0.47 0.49 0.57 0.86 1.43 2.04 4.76 1.45 0.95 0.68 0.48 0.35 0.39 0.41 0.39 0.55 0.88 1.17 1.68 4.08 1.25 0.96 0.63 0.44 0.49 0.41 0.47 0.53 0.68 0.84 1.37 2.10 5.00 1.62 0.95 0.50 0.56 0.43 0.42 0.48 0.62 0.65 0.85 1.27 0.94
Difference:  1.52 0.95 0.70 0.85 0.76 0.68 0.77 0.78 0.66 0.69 0.72 0.91 1.78 4.19 1.14 0.99 0.82 0.76 0.83 0.88 0.69 0.73 0.69 0.59 1.03 0.82 1.90 1.20 0.82 0.84 0.97 0.96 0.80 0.75 0.70 0.66 0.77 1.08 1.76 4.63 1.40 0.99 0.98 1.00 0.78 0.78 0.71 0.64 0.69 0.71 0.98 0.59 0.71 0.92 0.72 0.78 0.66 0.64 0.58 0.51 0.58 0.59 0.65 0.98 1.71 4.13 1.46 1.35 1.15 0.81 0.78 0.72 0.73 0.82 0.80 0.50 1.48 1.22 2.12 1.34 0.94 0.80 0.71 0.58 0.66 0.66 0.70 0.88 0.97 1.05 0.73 1.61 0.98 0.93 1.10 0.60 0.63 0.68 0.68 0.78 0.87 0.91 0.91 1.94 4.06 1.24 0.88 0.67 0.70 0.75 0.77 0.83 0.77 0.87 0.67 1.52 1.56 0.96 1.96 0.96 0.75 0.79 0.81 0.86 0.97 0.94 0.81 1.11 1.23 1.76 4.13 1.21 0.76 0.63 0.66 0.69 0.74 0.70 0.60 0.78 0.47 1.15 1.73 3.29 1.25 0.89 0.77 0.70 0.82 0.79 0.66 0.82 0.81 0.77 1.16 0.86 1.40 1.16 0.84 0.68 0.75 0.74 0.60 0.75 0.66 0.66 0.62 1.08 1.81 4.12 1.27 0.78 0.73 0.63 0.52 0.66 0.63 0.57 0.66 0.84 0.90 0.73 1.50 1.03 0.76 0.73 0.59 0.76 0.66 0.67 0.72 0.79 0.76 1.15 1.86 4.53 1.37 0.87 0.59 0.78 0.70 0.71 0.75 0.81 0.70 0.63 0.99 0.74
 
In matrix form:
Reference:
2.25 1.18 0.79 0.69 0.53 0.43 0.49 0.52 0.50 0.67 0.85 1.15 2.12  
4.80 1.36 0.95 0.66 0.51 0.51 0.55 0.47 0.59 0.65 0.67 1.17 1.12  
2.59 1.38 0.89 0.66 0.66 0.59 0.50 0.51 0.58 0.59 0.86 1.38 2.10  
5.02 1.60 0.99 0.79 0.69 0.46 0.50 0.47 0.47 0.59 0.84 1.12 0.96  
1.55 0.99 0.70 0.63 0.46 0.39 0.38 0.31 0.39 0.50 0.62 1.11 1.96  
4.51 1.50 0.95 0.61 0.62 0.48 0.45 0.47 0.57 0.65 0.75 1.57 2.08  
4.77 1.57 0.95 0.61 0.55 0.36 0.39 0.43 0.45 0.62 0.87 1.33 1.11  
2.12 1.13 0.86 0.54 0.44 0.41 0.43 0.45 0.48 0.70 0.84 1.07 2.12  
4.70 1.53 0.96 0.59 0.56 0.50 0.50 0.50 0.59 0.72 0.56 1.35 1.34  
0.31 1.66 1.03 0.69 0.68 0.57 0.53 0.63 0.60 0.60 1.00 1.41 2.03  
4.65 1.32 0.86 0.65 0.55 0.49 0.47 0.43 0.40 0.62 0.62 0.94 1.37  
0.61 1.45 1.07 0.74 0.59 0.57 0.47 0.41 0.52 0.63 0.75 1.31 1.13  
1.94 1.26 0.89 0.71 0.60 0.46 0.40 0.47 0.49 0.57 0.86 1.43 2.04  
4.76 1.45 0.95 0.68 0.48 0.35 0.39 0.41 0.39 0.55 0.88 1.17 1.68  
4.08 1.25 0.96 0.63 0.44 0.49 0.41 0.47 0.53 0.68 0.84 1.37 2.10  
5.00 1.62 0.95 0.50 0.56 0.43 0.42 0.48 0.62 0.65 0.85 1.27 0.94
 
Difference:
1.52 0.95 0.70 0.85 0.76 0.68 0.77 0.78 0.66 0.69 0.72 0.91 1.78  
4.19 1.14 0.99 0.82 0.76 0.83 0.88 0.69 0.73 0.69 0.59 1.03 0.82  
1.90 1.20 0.82 0.84 0.97 0.96 0.80 0.75 0.70 0.66 0.77 1.08 1.76  
4.63 1.40 0.99 0.98 1.00 0.78 0.78 0.71 0.64 0.69 0.71 0.98 0.59  
0.71 0.92 0.72 0.78 0.66 0.64 0.58 0.51 0.58 0.59 0.65 0.98 1.71  
4.13 1.46 1.35 1.15 0.81 0.78 0.72 0.73 0.82 0.80 0.50 1.48 1.22  
2.12 1.34 0.94 0.80 0.71 0.58 0.66 0.66 0.70 0.88 0.97 1.05 0.73  
1.61 0.98 0.93 1.10 0.60 0.63 0.68 0.68 0.78 0.87 0.91 0.91 1.94  
4.06 1.24 0.88 0.67 0.70 0.75 0.77 0.83 0.77 0.87 0.67 1.52 1.56  
0.96 1.96 0.96 0.75 0.79 0.81 0.86 0.97 0.94 0.81 1.11 1.23 1.76  
4.13 1.21 0.76 0.63 0.66 0.69 0.74 0.70 0.60 0.78 0.47 1.15 1.73  
3.29 1.25 0.89 0.77 0.70 0.82 0.79 0.66 0.82 0.81 0.77 1.16 0.86  
1.40 1.16 0.84 0.68 0.75 0.74 0.60 0.75 0.66 0.66 0.62 1.08 1.81  
4.12 1.27 0.78 0.73 0.63 0.52 0.66 0.63 0.57 0.66 0.84 0.90 0.73  
1.50 1.03 0.76 0.73 0.59 0.76 0.66 0.67 0.72 0.79 0.76 1.15 1.86  
4.53 1.37 0.87 0.59 0.78 0.70 0.71 0.75 0.81 0.70 0.63 0.99 0.74

Would you mind to take a look and enlighten me on what might be wrong here? Thank you so much!

Hi @liubenyuan,
I was using pyEIT its been some time. I made some eit hardware to be able to take some real data and use it with pyEIT. Actually, I would like to use it with GREIT . Can you give me some advice about the constants, the current and the frequency to use?

Thanks you very much for your work!

Hi! Im using this code to generate the electrodes in the mesh of lungs. I'm using the 32 electrodes configuration:
# plot the mesh, visualize the electrodes locations fig, ax = plt.subplots(figsize=(6, 6)) ax.triplot(pts[:, 0], pts[:, 1], tri, linewidth=1) ax.plot(pts[mesh_obj.el_pos, 0], pts[mesh_obj.el_pos, 1], "ro") for i, el in enumerate(mesh_obj.el_pos): ax.text(pts[el, 0], pts[el, 1], str(i + 1), fontsize=15) ax.axis("equal") ax.axis([-1.2, 1.2, -1.2, 1.2])

and having this result :

is there any problem in my code to generate the good positions?
Is there any way to make a better code?

thank you

Hello,

I'm trying to use pyEIT at work for some project but I have trouble configuring it compared to EIDORS (Matlab) for an accurate size and shape of object (with dynamic jac).

When I use your example project (demo_dynamic_jac) and I change the size of the anomaly, I get this (which is perfect) :
anomaly 1 = [{'x': 0.5, 'y': 0.5, 'd': 0.1, 'perm': 1000.0}]
anomaly 2 = [{'x': 0.5, 'y': 0.5, 'd': 0.3, 'perm': 1000.0}]
reconstruction anomaly 1 << reconstruction anomaly 2 (OK)
(ie biggest anomaly = biggest activation zone)

However if I change x or y position to get closer to the center, even by a little, it's impossible to see the difference in size :
anomaly 3 = [{'x': 0.4, 'y': 0.4, 'd': 0.1, 'perm': 1000.0}]
anomaly 4 = [{'x': 0.4, 'y': 0.4, 'd': 0.3, 'perm': 1000.0}]
reconstruction anomaly 3 == reconstruction anomaly 4 (BUG)

anomaly 5 = [{'x': 0, 'y': 0, 'd': 0.1, 'perm': 1000.0}]
anomaly 6 = [{'x': 0, 'y': 0, 'd': 0.6, 'perm': 1000.0}]
reconstruction anomaly 5 == reconstruction anomaly 6 (BUG)

Thank you in advance for your explanations ! and good job anyway for pyEIT, it's awesome 👍

Hi,

During the review of the code I was wondering why you don't set a dataclass for the mesh: you could than directly add method or property and allow you to do more with the data itsel without having to unpack the dict each time etc...

for example:

@dataclass
class PyeitMesh:
    element:np.ndarray
    node:np.ndarray
    perm:np.ndarray


    @property
    def n_pts(self)->int:
        return self.node.shape[0]
    @property
    def n_dim(self)->int:
        return self.node.shape[1]
    @property
    def n_elems(self)->int:
        return self.element.shape[0]
    @property
    def n_vertices(self)->int:
        return self.element.shape[1]


    def stats(self):
        """
        print mesh or tetrahedral status

        Parameters
        ----------
        p: array_like
            coordinates of nodes (x, y) in 2D, (x, y, z) in 3D
        t: array_like
            connectives forming elements

        Notes
        -----
        a simple function for illustration purpose only.
        print the status (size) of nodes and elements
        """
        print("mesh status:")
        print(f"{self.n_pts} nodes, {self.n_elems} elements")

then you can avoid such code:

self.pts = mesh["node"]
self.tri = mesh["element"]

shape of the mesh
self.no_num, self.n_dim = self.pts.shape
self.el_num, self.n_vertices = self.tri.shape

and use that instead

self.mesh.node #(instead self.pts )
self.mesh.n_pts #(instead self.no_num ), etc.. 

also you can integrate setting of permittivity... and much more whcih is on the mesh data related!
If you are interested I could start to set it!

Hi, it's me again.
Thanks for your support, we have managed to improve our reconstructed image quality.

We are currently developing a feature to read data from the tank and plot the image along with the time. To achieve that, we modified the code as follows:

# coding: utf-8
""" demo on dynamic eit using JAC method """
# Copyright (c) Benyuan Liu. All Rights Reserved.
# Distributed under the (new) BSD License. See LICENSE.txt for more info.
from __future__ import absolute_import, division, print_function

import matplotlib.pyplot as plt
import numpy as np
import pyeit.eit.jac as jac
import pyeit.mesh as mesh
from pyeit.eit.fem import EITForward
from pyeit.eit.interp2d import sim2pts
from pyeit.mesh.shape import thorax
import pyeit.eit.protocol as protocol
from pyeit.mesh.wrapper import PyEITAnomaly_Circle

import serial
from datetime import datetime
import time
from matplotlib.animation import FuncAnimation

"""-2. Initial vars """
arduino = serial.Serial('COM8', 115200 ,timeout=4)
v0 = np.loadtxt('example_data/ref_data.txt')
fig, ax = plt.subplots(constrained_layout=True)
n_el = 16

"""-1. Functions """
'''To read data from Arduino via COM port
Each frame is 16 lines, each line has 13 values, 
in total there are 208 values representing the voltage measured from electrodes. 
Frames are separated by an enter character.
'''

def readfromArduino():
    while(True):
        try:
            data = arduino.readline().decode('ascii')
            print("data: ", data)
            break
        except UnicodeDecodeError:
            print("UnicodeDecodeError found! Retrying...")
            continue
    return data

def get_difference_img_array(n_el = n_el, NewFrameSearchFlag = 1, idx = 0):
    difference_image_array = ''
    # Read difference image f1:
    while idx < n_el:
        data = readfromArduino()
        #skip until the empty line is found to catch the whole frame
        while(NewFrameSearchFlag == 1):
            if  len(data) < 4:
                print("Searching for new frame.")
                data = readfromArduino()
                continue
            else:
                print("New frame found.")
                data = readfromArduino()
                NewFrameSearchFlag = 0
                break


        data=data.strip('\r\n')
        difference_image_array += data
        difference_image_array += ' '
        idx = idx + 1
        #print("String: {0}".format(data))
 
    return difference_image_array

#Convert data to np type
def convert_data_in(s):
    data=s
    items=[]
    for item in data.split(' '):
        item = item.strip()
        if not item:
            continue
        try:
            items.append(float(item))
        except ValueError:
            print("Value Error found! Handling...")
            items.append(float(0))
    return np.array(items)


""" 0. build mesh """
n_el = 16  # nb of electrodes
use_customize_shape = False
if use_customize_shape:
    # Mesh shape is specified with fd parameter in the instantiation, e.g : fd=thorax
    mesh_obj = mesh.create(n_el, h0=0.065, fd=thorax)
else:
    mesh_obj = mesh.create(n_el, h0=0.065)

# extract node, element, alpha
pts = mesh_obj.node
tri = mesh_obj.element
x, y = pts[:, 0], pts[:, 1]

""" 1. problem setup """
# mesh_obj["alpha"] = np.random.rand(tri.shape[0]) * 200 + 100 # NOT USED
anomaly = PyEITAnomaly_Circle(center=[0.5, 0.5], r=0.1, perm=1000.0)
mesh_new = mesh.set_perm(mesh_obj, anomaly=anomaly)

""" 2. FEM simulation """
# setup EIT scan conditions
protocol_obj = protocol.create(n_el, dist_exc=1, step_meas=1, parser_meas="fmmu")

# calculate simulated data
fwd = EITForward(mesh_obj, protocol_obj)

#v0 = fwd.solve_eit()
#v1 = fwd.solve_eit(perm=mesh_new.perm)
""" 3. JAC solver """
# Note: if the jac and the real-problem are generated using the same mesh,
# then, data normalization in solve are not needed.
# However, when you generate jac from a known mesh, but in real-problem
# (mostly) the shape and the electrode positions are not exactly the same
# as in mesh generating the jac, then data must be normalized.
eit = jac.JAC(mesh_obj, protocol_obj)
eit.setup(p=0.5, lamb=0.001, method="kotre", perm=1000, jac_normalized=True)

def animating(i):  
    while arduino.inWaiting()==0:
        print("waiting")
        time.sleep(0.5)
        pass

    s1 = get_difference_img_array()

    v1 = convert_data_in(s1)

    try:
        ds = eit.solve(v1, v0, normalize=True)
    except Exception as e:
        ani.event_source.stop()  # Stop the current animation if error occurred
        ani.event_source.start()  # Start a new animation

    ds_n = sim2pts(pts, tri, np.real(ds))
  
    # Clear the graph after each animating frame
    ax.clear()
    
    # plot EIT reconstruction
    im = ax.tripcolor(x, y, tri, ds_n, shading="flat", cmap=plt.cm.magma)
    for i, e in enumerate(mesh_obj.el_pos):
        ax.annotate(str(i + 1), xy=(x[e], y[e]), color="r")
    ax.set_aspect("equal")

    # plt.savefig('../doc/images/demo_jac.png', dpi=96)

ani = FuncAnimation(fig, animating, interval = 50, cache_frame_data= False)
plt.show()

It is noteworthy that our electrodes were messy, and there are obvious contact problems with some, which severely affect the quality of the reconstructed image. Still, it should not be a big problem since we can replace them in no time. Besides, it works. We managed to show the movement of the objects inside our tank. You can watch the video at this Google Drive link: https://drive.google.com/drive/folders/1Y_9w8cCYSJimbiGiwunk6UIfk7JaeHzk?usp=sharing

However, it's clear to be seen that our fps is a little bit modest, and some frames are bad-reconstructed (we are currently fixing our hardware so the signal will me much more stable which should fix this problem).

I would like to ask you (and all other people who are reading this issue) to take a look at my code and give me some feedback.

Is there a better way to do this real-time data plotting? Please enlighten me.

Thank you so much!

As I'm trying to get used to the package and its possibilities, the Readme and examples a couple errors, and I think the figures could be more clear as for before/after.

I'll PR with a proposition of changes soon, but as I'm on it don't hesitate if you think there are other changes I should do at the same time.

On paper, pyEIT seems stellar. In practice, pyEIT has yet to publish a stable release to either PyPI or conda-forge. That's less stellar, because it means that:

• Users without low-level CLI Git[Hub] experience cannot install pyEIT. I mean, all of us can, of course – but we're all open-source Pythonistas here. Meanwhile, the main audience is biotech researchers. They typically lack that same UNIX-oriented coder skillset.
• No other project can depend on pyEIT. This includes projects like:
  • OpenEIT, which (despite being pretty much dead now) could have benefited from just depending on pyEIT for its image reconstruction analysis.
  • py-dbar, Python's first and only D-bar EIT implementation. Did you know py-dbar even existed? I didn't. It looks awesome though, right? py-dbar even references pyEIT in its introduction: "This package implements the Dbar method for Electrical impedance tomography and follows up on the pyEIT presented in https://github.com/liubenyuan/pyEIT." Sadly, py-dbar can't even import pyEIT (e.g., to conform to the same public API), because pyEIT isn't hosted on PyPI and thus can't be referenced from the setup.py or pyproject.toml files of downstream projects that are hosted on PyPI.

Are there any significant roadblocks to publication? There probably are, because it hasn't happened yet. But if there aren't, it'd be really swell if pyEIT 2.0.0 (or whatever the current version is) could be officially released as soon as feasible. We'd all love to begin hacking around with this wonderful thing.

Thanks for all the tireless volunteerism, @liubenyuan, et al. EIT is tough stuff. Let's see if we can change that together. 😄

Cordial greetings yet again, Benyuan. Yup, it's @leycec. As noted over at #23, I'm currently packaging pyEIT for Gentoo Linux. It's happening. But...

I hit a trivial installation issue. Thankfully, it's trivial! pyEIT is currently installing the ambiguous tests package, which is kinda bad. For safety, Gentoo prohibits this from happening with a QA error resembling:

 * ERROR: dev-python/pyEIT-9999::raiagent failed (install phase):
 *   Package installs 'tests' package which is forbidden and likely a bug in the build system.

The underlying issue is this one-liner in your setup.py installer:

    packages=find_packages(exclude=["test"]),

That's... not working as intended. Firstly, there is no test package; you probably meant tests there. Secondly, even that fails to suffice. Why? Because the exclude parameter of the setuptools.find_packages() function doesn't work the way everyone wants it to work. Ideally, that parameter should implicitly exclude all subpackages of the listed packages. Right? It doesn't. You have to explicitly list an additional glob expression to exclude all subpackages.

Altogether, the above line should instead resemble:

    packages=find_packages(exclude=["tests", "tests.*"]),

Thanks for all the open-sourced excellence in EIT, Benyuan! We'll get pyEIT packaged up for you shortly and usable by at least Gentoo Linux researchers. 👋

When trying to run this example, using a fresh environment installed from the .yml file, I get this error:

Traceback (most recent call last):
  File "...\pyEIT\examples\fem_forward3d.py", line 10, in <module>
    import pyeit.mesh.plot as mplot
  File "...\pyEIT\pyeit\mesh\plot\__init__.py", line 5, in <module>
    from .tetplot import tetplot
  File "...\pyEIT\pyeit\mesh\plot\tetplot.py", line 15, in <module>
    from vispy import app, gloo, scene
  File "C:\Users\amaury\.conda\envs\pyeit-environment\lib\site-packages\vispy\app\__init__.py", line 13, in <module>
    from .canvas import Canvas, MouseEvent, KeyEvent  # noqa
  File "C:\Users\amaury\.conda\envs\pyeit-environment\lib\site-packages\vispy\app\canvas.py", line 12, in <module>
    from ..util.ptime import time
  File "C:\Users\amaury\.conda\envs\pyeit-environment\lib\site-packages\vispy\util\ptime.py", line 36, in <module>
    cstart = systime.clock()  # Required to start the clock in windows
AttributeError: module 'time' has no attribute 'clock'

Any idea why?

I've noticed that running JAC and Back Projection algorithm result in ds format which can be overlayed on the mesh as shown below:

Is there anyway to obtain ds in GREIT algorithm in the same format to make plotting it consistent across algorithms? Currently it outputs an n by n array.

hi @liubenyuan

Hope you are doing well. Just today, I installed pyEIT using pip and I am getting TypeError's while importing stuff.
I am just trying to run the fem_forward2d.py. I have attached a screenshot below for reference.

Thanks for the great work.

Hi @liubenyuan ,

I see all your examples are defined on discs (2d) or balls (3d). Do you have functions that define the problem on rectangle domain with regular mesh? Thanks!

You can open a new issue on this.

You should pass parser=fmmu if you are using rotate-meas protocol. i.e., If your electrodes are numbered as [1..8], if now (the stimulation pattern) the the source and sink of current are 4-5, then the 5 voltgaes differences are V67, V78, V81, V12, V23.

The following are the codes, you may try using interpolation for a smoothed display

# coding: utf-8
""" demo on static solving using JAC (experimental) """
# Copyright (c) Benyuan Liu. All Rights Reserved.
# Distributed under the (new) BSD License. See LICENSE.txt for more info.
from __future__ import division, absolute_import, print_function

import numpy as np
import matplotlib.pyplot as plt

# pyEIT 2D algorithms modules
from pyeit.mesh import create
from pyeit.eit.utils import eit_scan_lines
import pyeit.eit.jac as jac

""" 1. setup """
n_el = 8
mesh_obj, el_pos = create(n_el, h0=0.12)
el_dist, step = 1, 1
ex_mat = eit_scan_lines(n_el, el_dist)
pts = mesh_obj["node"]
tri = mesh_obj["element"]

fv = np.array([0.1887, 0.0921, 0.0565, 0.065, 0.1809,
               0.1602, 0.0666, 0.0598,0.1022, 0.1538,
               0.1526, 0.0867, 0.0996, 0.0932, 0.1436,
               0.1547,0.106, 0.0703, 0.0672, 0.1442,
               0.1742, 0.076, 0.0527, 0.0763,0.1893,
               0.1381, 0.0621, 0.0667, 0.1047, 0.1907,
               0.1745, 0.1185,0.1263, 0.1361, 0.183,
               0.1906, 0.124, 0.0955, 0.0823, 0.1522]) # in volt

""" 3. solve_eit using gaussian-newton (with regularization) """
# number of stimulation lines/patterns
eit = jac.JAC(mesh_obj, el_pos, ex_mat, step, perm=1.0, parser="fmmu")
eit.setup(p=0.25, lamb=1.0, method="lm")
# lamb = lamb * lamb_decay
ds = eit.gn(fv, lamb_decay=0.1, lamb_min=1e-5, maxiter=20, verbose=True)

# plot
fig, ax = plt.subplots(figsize=(6, 4))
im = ax.tripcolor(
    pts[:, 0],
    pts[:, 1],
    tri,
    np.real(ds),
    shading="flat",
    alpha=1.0,
    cmap=plt.cm.viridis,
)
fig.colorbar(im)
ax.axis("equal")
ax.set_title("Conductivities Reconstructed")
# fig.savefig('../figs/demo_static.png', dpi=96)
plt.show()

Originally posted by @liubenyuan in https://github.com/liubenyuan/pyEIT/issues/15#issuecomment-741722298

Hi I have a question.
when doing this function when transforming from mesh to voltage
f1 = fwd.solve_eit(ex_mat, step=step, perm=mesh_obj["perm"])
I received some negative voltage on the 208 array(f1.v).

I don't understand why I am receiving negative data. Using conventional data acquisition on a human, I don't receive any negative voltage. So I'm trying to have the same results from the algorithm.

Thanks for the help :)

hi benyuan, nice work! Do you have any plan to add functionality to export 3d image to png file? Your 'tetplot' function plots figures in interactive mode. But it would be nice to save them. Thanks!

hi @liubenyuan, I'm an undergraduate biomedical student and I'm new to EIT and I was interested in this pyEIT project. I looked into all the examples that are given. however, I'm curious if pyEIT can be applied in larger models like human lungs datasets or it's just coded for small prototypes. if it's possible to apply in larger models may I know it is done? thank you.

Hi.

I tested the example "eit_dynamic_bp.py".
It runs well at the default setting,
but if I change the scanning conditions,

el_dist, step = 1, 1 ==> el_dist, step = 8, 1

the result is weird like as below.

Does the BP solver not support the opposite stimulation protocol?
or Is there anything i am missing?

Thank you in advance.

Hi,

I am trying to generate a square mesh for EIT reconstruction, but not sure how to do that.
I am generating the pts and passing to the mesh.create function like this

    p1 = np.array([0, 0])  # bottom-left
    p2 = np.array([5, 5]) # top-right

    # Generate a grid of points within this rectangle
    x = np.linspace(p1[0], p2[0], num=50)  # 20 points along x
    y = np.linspace(p1[1], p2[1], num=50)  # 10 points along y
    xx, yy = np.meshgrid(x, y)
    pts_r = np.c_[xx.ravel(), yy.ravel()]
    
    rect = rectangle0(pts_r)
    mesh_obj = mesh.create(n_electrodes, fd=rect)  

however, I get the following error:
Exception has occurred: ValueError
The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

thanks

Thank you for your program, it helped me a lot. Can I ask you a question? where is the progamm, defines the input current amplitude, or the injecting current is considered/set as 1A or 1mA, automatically.

Sorry for the potentially dumb question:

I'm new to EIT but got access to lung EIT data from COVID patients, saved in .eit format. I would like to visualize the data and visually inspect it. I installed the package, but I'm lost on how to start. Do you have any standalone documenation website (e.g. made with Sphinx and hosted with github pages)?
Basically, I'm just searching for a way to open .eit files on Mac OS.

Thanks a lot!

Hello @liubenyuan,

I hope you are doing well. First of all, thank you so much for developing and sharing your pyEIT library with everyone. I am fairly new to using your library so, I have some questions for you with regards to the Forward class provided below.
https://github.com/liubenyuan/pyEIT/blob/eaa0b7981b6ab0735e25ba079c1ecc209604c438/pyeit/eit/fem.py#L17

If I understand correctly, you provide the Jacobian of the body with respect to the electric potential "u" and conductivity "sigma" using solve and solve_eit functions respectively from the above Forward class.
So, "solve" computes the forward problem by outputting the point wise potential u values and Jacobian of u.
While, "solve_eit" computes the inverse problem by outputting the Jacobian of sigma, boundary voltages at electrodes and Smear matrix.
Please correct my understanding if I am wrong here.

I am particularly interested in the Jacobian information from your library. My main questions are:

1. These Jacobians are w.r.t the simplex triangles i.e, df(body, sigma)/d sigma and df(body, u)/d u at each triangle on the mesh correct?
2. Do you know if there is a way I can compute the individual gradients du/dx, du/dy, d sigma/dx and d sigma/dy for each x, y on the mesh?

Any advice is greatly appreciated,

Thank you so much,

Hi, I've been trying to learn to use my SpectraEIT. I've installed the pyeit lib and was looking at the examples but I run into this invalid synthax error constantly.

Other similar erros I've found are in from pyeit.eit.interp2d import sim2pts with:

  File "/usr/local/lib/python3.7/dist-packages/pyeit/eit/base.py", line 13, in <module>
    from pyeit.eit.protocol import PyEITProtocol

  File "<fstring>", line 1
    (ex_mat=)

The version I am using is pyeit = 1.2.2

Is there a way to fix it? Or other version I should use? I really want to use this for my project with the SpectraEIT.

The reconstruction images of eit_static_jac.py is wrong with respect to the result of commit e6894ca

checking.

Hi,

I hope this message finds you well.

First, I want to express my appreciation for the excellent work done on pyEIT. It has been incredibly helpful for my work with electrical impedance tomography.

I have a question regarding the color scale in the images generated. When using the provided codes, I get an image that maps the conductivities in a region, along with a colourbar displaying a gradient of colors (e.g., from black to white through various shades of gray). Each color in this gradient represents a specific value. Could you please clarify what units these values on the color scale are in?

Thank you very much for your assistance and for this valuable tool.

Best regards,

Pablo Sánchez

hi benyuan, I want to ask a naive question. In your model, you do not restrict permittivity \sigma to be positive. This makes the model extremely sensitive to sigma (esp. in negative domain) and resulting log-likelihood (objective function) very rough. According to your paper in 2018, it is for analysis of multi-frequency EIT (possibly complex solutions)? Is there any realistic consideration for not requiring positive permittivity? What is it? Thanks!

While running some tests on the eit_static_jac example, I came accross an issue. It seems this method is broken, as the gtol is going crazy and the results are clearly wrong.
You can recreate it simply by putting the method flag on 'lm' for the eit.gn function.

Hi, for my master's I've been using pyEIT for some work into characterising graphene, which means we'll be injecting an excitation current of 10 uA (1e-5 A) and typically measuring voltages of the mV order. We're mainly using GREIT, which is good for identifying regions of interest, but I'm struggling to determine how the actual units of conductance match the colour bar scale. Another complication is the regularisation parameter lambda, since the actual values of the scale vary for different lambda.

When taking measurements of a resistor grid (a test stand in for graphene), we multipled all of the voltages by 1e5 to match the assumption within pyEIT that the input current is 1A. The GREIT plot scales would vary from [0.0,-0.02] for lamb=0.1 and [100, -500] for lamb=0.01. I wasn't sure how to relative the units of sieverts to the value of lambda and the scales, and hence couldn't assign units of conductivity to these scales. I could plot the L-curve for tikohonov regularisation, which would give me the 'optimal' value of lambda to use for my setup. However, this still wouldn't address the scaling and how exactly it relates to the absolute conductivity I'm measuring.

Below are the reconstructions for some oval anomaly using the opp-adj stimulation pattern for 64 electrodes around a square sample. All of the anomaly data was simulated, with a homogenous reference conductivity of 1 Sievert for the forward solution.

p = 0.2, lamb = 0.1

p= 0.2, lamb - 0.01

In the simulated case above, assuming an input current of 1 A. Would this then mean the scale on the 'true' image is in Sieverts, or Sieverts/metre, thus the blue region is -0.5 sieverts below the unit reference of 1 sievert. Hence blue = 0.5 S and yellow = 1 S. This then confuses me as to how I should relate the scale of GREIT, difference-JAC and GN to the true scale and hence obtain the absolute scaling of conductivity? Also, there must be a way of relating the scale conductivities and the lambda used to enable fixation of the units ...

Any advice on how to proceed for fixing units to this conductivity scale would be very appreciated. I'm not after a high precision of conductivity scaling, rather just the correct order of magnitudes; uS, mS, S. The only think I can think of is to use the van der Pauw method to calculate the sheet conductance of my sample and then use that value as my reference conductivity.

Many thanks,
Adam