Add PyEITMesh, PyEITProtocol, PyEITAnomaly (#47)

* Build PyEITMesh Dataclass
- add ref_el (chekc and set (for user))
- add prop n_el
- replace internal attribute by calling directly the date from the mesh object (e.g. self.mesh.node, etc)
- now a lot of size of date can be accessed from self.mesh.... (beetr for pylance or u.a)
- use dataset of mesh directly!
- make tri_area and _hull_points for 3D-pts (x,y,z ) compatible
- use in multi_shell

* init PyEITProtocol, use of PyEITProtocl in base and fem
- remove all usued attributes
- in fem remove checks of data (less responsabiltity for fwd)
- use of the package warning to make non blocking warnings
- list all examples + selection
- use the PyeitProtocol cls to build a protocol_obj using pyeit.eit.protocol.create wrapper function
- pyeit.eit.protocol.create can accepe multiple el_dist (for stacked excitation)
- add n_el in examples..
- rename eit_scan_lines (utils.py is now oselete, ... not used, ) to build_exc_pattern_std

* Add PyEITAnomaly
- correct error in set_perm (d was not diamerter but r...)

examples/eit_dynamic_bp.py

@@ -2,57 +2,57 @@
 """ demo code for back-projection """
 # 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
+from __future__ import absolute_import, division, print_function
 
-import numpy as np
 import matplotlib.pyplot as plt
-
+import numpy as np
+import pyeit.eit.bp as bp
+import pyeit.eit.protocol as protocol
 import pyeit.mesh as mesh
 from pyeit.eit.fem import EITForward
-from pyeit.eit.utils import eit_scan_lines
 from pyeit.mesh.shape import thorax
-import pyeit.eit.bp as bp
+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(16, h0=0.1, fd=thorax)
+    mesh_obj = mesh.create(n_el, h0=0.1, fd=thorax)
 else:
-    mesh_obj = mesh.create(16, h0=0.1)
+    mesh_obj = mesh.create(n_el, h0=0.1)
 
 """ 1. problem setup """
-anomaly = [{"x": 0.5, "y": 0.5, "d": 0.1, "perm": 10.0}]
+
+anomaly = PyEITAnomaly_Circle(center=[0.5, 0.5], r=0.1, perm=10.0)
 mesh_new = mesh.set_perm(mesh_obj, anomaly=anomaly, background=1.0)
 
 """ 2. FEM forward simulations """
 # setup EIT scan conditions
-# adjacent stimulation (el_dist=1), adjacent measures (step=1)
-el_dist, step = 1, 1
-ex_mat = eit_scan_lines(16, el_dist)
-protocol = {"ex_mat": ex_mat, "step": step, "parser": "std"}
+# adjacent stimulation (dist_exc=1), adjacent measures (step_meas=1)
+protocol_obj = protocol.create(n_el, dist_exc=1, step_meas=1, parser_meas="std")
 
 # calculate simulated data
-fwd = EITForward(mesh_obj, protocol)
+fwd = EITForward(mesh_obj, protocol_obj)
 v0 = fwd.solve_eit()
-v1 = fwd.solve_eit(perm=mesh_new["perm"], init=True)
+v1 = fwd.solve_eit(perm=mesh_new.perm, init=True)
 
 """ 3. naive inverse solver using back-projection """
-eit = bp.BP(mesh_obj, protocol)
+eit = bp.BP(mesh_obj, protocol_obj)
 eit.setup(weight="none")
 ds = 192.0 * eit.solve(v1, v0, normalize=False)
 
 # extract node, element, alpha
-pts = mesh_obj["node"]
-tri = mesh_obj["element"]
+pts = mesh_obj.node
+tri = mesh_obj.element
 
 # draw
 fig, axes = plt.subplots(2, 1, constrained_layout=True, figsize=(6, 9))
 # original
 ax = axes[0]
 ax.axis("equal")
 ax.set_title(r"Input $\Delta$ Conductivities")
-delta_perm = np.real(mesh_new["perm"] - mesh_obj["perm"])
+delta_perm = np.real(mesh_new.perm - mesh_obj.perm)
 im = ax.tripcolor(pts[:, 0], pts[:, 1], tri, delta_perm, shading="flat")
 # reconstructed
 ax1 = axes[1]

examples/eit_dynamic_greit.py

@@ -9,49 +9,49 @@
 
 import pyeit.mesh as mesh
 from pyeit.eit.fem import EITForward
-from pyeit.eit.utils import eit_scan_lines
+import pyeit.eit.protocol as protocol
 from pyeit.mesh.shape import thorax
 import pyeit.eit.greit as greit
+from pyeit.mesh.wrapper import PyEITAnomaly_Circle
 
 """ 0. construct 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(16, h0=0.1, fd=thorax)
+    mesh_obj = mesh.create(n_el, h0=0.1, fd=thorax)
 else:
-    mesh_obj = mesh.create(16, h0=0.1)
+    mesh_obj = mesh.create(n_el, h0=0.1)
 
 # extract node, element, alpha
-pts = mesh_obj["node"]
-tri = mesh_obj["element"]
+pts = mesh_obj.node
+tri = mesh_obj.element
 
 """ 1. problem setup """
 # this step is not needed, actually
 # mesh_0 = mesh.set_perm(mesh_obj, background=1.0)
 
 # test function for altering the 'permittivity' in mesh
 anomaly = [
-    {"x": 0.4, "y": 0, "d": 0.1, "perm": 10},
-    {"x": -0.4, "y": 0, "d": 0.1, "perm": 10},
-    {"x": 0, "y": 0.5, "d": 0.1, "perm": 0.1},
-    {"x": 0, "y": -0.5, "d": 0.1, "perm": 0.1},
+    PyEITAnomaly_Circle(center=[0.4, 0], r=0.1, perm=10.0),
+    PyEITAnomaly_Circle(center=[-0.4, 0], r=0.1, perm=10.0),
+    PyEITAnomaly_Circle(center=[0, 0.5], r=0.1, perm=0.1),
+    PyEITAnomaly_Circle(center=[0, -0.5], r=0.1, perm=0.1),
 ]
 mesh_new = mesh.set_perm(mesh_obj, anomaly=anomaly, background=1.0)
-delta_perm = np.real(mesh_new["perm"] - mesh_obj["perm"])
+delta_perm = np.real(mesh_new.perm - mesh_obj.perm)
 
 """ 2. FEM forward simulations """
 # setup EIT scan conditions
-el_dist, step = 1, 1
-ex_mat = eit_scan_lines(16, el_dist)
-protocol = {"ex_mat": ex_mat, "step": step, "parser": "std"}
+protocol_obj = protocol.create(n_el, dist_exc=1, step_meas=1, parser_meas="std")
 
 # calculate simulated data
-fwd = EITForward(mesh_obj, protocol)
+fwd = EITForward(mesh_obj, protocol_obj)
 v0 = fwd.solve_eit()
-v1 = fwd.solve_eit(perm=mesh_new["perm"], init=True)
+v1 = fwd.solve_eit(perm=mesh_new.perm, init=True)
 
 """ 3. Construct using GREIT """
-eit = greit.GREIT(mesh_obj, protocol)
+eit = greit.GREIT(mesh_obj, protocol_obj)
 eit.setup(p=0.50, lamb=0.001)
 ds = eit.solve(v1, v0)
 x, y, ds = eit.mask_value(ds, mask_value=np.NAN)

examples/eit_dynamic_jac.py

@@ -2,54 +2,53 @@
 """ 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 division, absolute_import, print_function
+from __future__ import absolute_import, division, print_function
 
-import numpy as np
 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.utils import eit_scan_lines
-
-from pyeit.mesh.shape import thorax
-import pyeit.eit.jac as jac
 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(16, h0=0.1, fd=thorax)
+    mesh_obj = mesh.create(n_el, h0=0.1, fd=thorax)
 else:
-    mesh_obj = mesh.create(16, h0=0.1)
+    mesh_obj = mesh.create(n_el, h0=0.1)
 
 # extract node, element, alpha
-pts = mesh_obj["node"]
-tri = mesh_obj["element"]
+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
-anomaly = [{"x": 0.5, "y": 0.5, "d": 0.1, "perm": 1000.0}]
+# 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 """
-el_dist, step = 8, 1
-ex_mat = eit_scan_lines(16, el_dist)
-protocol = {"ex_mat": ex_mat, "step": step, "parser": "std"}
+# setup EIT scan conditions
+protocol_obj = protocol.create(n_el, dist_exc=8, step_meas=1, parser_meas="std")
 
 # calculate simulated data
-fwd = EITForward(mesh_obj, protocol)
+fwd = EITForward(mesh_obj, protocol_obj)
 v0 = fwd.solve_eit()
-v1 = fwd.solve_eit(perm=mesh_new["perm"], init=True)
+v1 = fwd.solve_eit(perm=mesh_new.perm, init=True)
 
 """ 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)
+eit = jac.JAC(mesh_obj, protocol_obj)
 eit.setup(p=0.5, lamb=0.01, method="kotre")
 ds = eit.solve(v1, v0, normalize=True)
 ds_n = sim2pts(pts, tri, np.real(ds))
@@ -59,14 +58,14 @@
 fig.set_size_inches(9, 4)
 
 ax = axes[0]
-delta_perm = mesh_new["perm"] - mesh_obj["perm"]
+delta_perm = mesh_new.perm - mesh_obj.perm
 im = ax.tripcolor(x, y, tri, np.real(delta_perm), shading="flat")
 ax.set_aspect("equal")
 
 # plot EIT reconstruction
 ax = axes[1]
 im = ax.tripcolor(x, y, tri, ds_n, shading="flat")
-for i, e in enumerate(mesh_obj["el_pos"]):
+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")
 

examples/eit_dynamic_jac3d.py

@@ -2,61 +2,55 @@
 """ demo on JAC 3D, extremely slow """
 # 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
+from __future__ import absolute_import, division, print_function
 
 import numpy as np
-
+import pyeit.eit.jac as jac
+import pyeit.eit.protocol as protocol
 import pyeit.mesh as mesh
-from pyeit.mesh import quality
 import pyeit.mesh.plot as mplot
 from pyeit.eit.fem import EITForward
 from pyeit.eit.interp2d import sim2pts
-from pyeit.eit.utils import eit_scan_lines
 from pyeit.mesh.shape import ball
-import pyeit.eit.jac as jac
+from pyeit.mesh.wrapper import PyEITAnomaly_Ball
 
 # build tetrahedron
 # 3D tetrahedron must have a bbox
 bbox = [[-1, -1, -1], [1, 1, 1]]
 # save calling convention as distmesh 2D
 # 3D Mesh shape is specified with fd parameter in the instantiation, e.g : fd=ball , Default in 3D :fd=ball
-mesh_obj = mesh.create(h0=0.2, bbox=bbox, fd=ball)
+n_el = 16  # nb of electrodes
+mesh_obj = mesh.create(n_el, h0=0.2, bbox=bbox, fd=ball)
 
-pts = mesh_obj["node"]
-tri = mesh_obj["element"]
+pts = mesh_obj.node
+tri = mesh_obj.element
 
-# report the status of the 2D mesh
-quality.stats(pts, tri)
+# report the status of the mesh
+mesh_obj.print_stats()
 
 """ 1. FEM forward simulations """
 # setup EIT scan conditions
-el_dist, step = 7, 1
-ex_mat = eit_scan_lines(16, el_dist)
-protocol = {"ex_mat": ex_mat, "step": step, "parser": "std"}
+protocol_obj = protocol.create(n_el, dist_exc=7, step_meas=1, parser_meas="std")
 
 # calculate simulated data
-fwd = EITForward(mesh_obj, protocol)
-
-# in python, index start from 0
-ex_line = ex_mat[2].ravel()
+fwd = EITForward(mesh_obj, protocol_obj)
 
 # change alpha
-anomaly = [{"x": 0.40, "y": 0.40, "z": 0.0, "d": 0.30, "perm": 100.0}]
+anomaly = PyEITAnomaly_Ball(center=[0.4, 0.4, 0.0], r=0.3, perm=100.0)
 mesh_new = mesh.set_perm(mesh_obj, anomaly=anomaly, background=1.0)
-tri_perm = mesh_new["perm"]
-node_perm = sim2pts(pts, tri, np.real(tri_perm))
+node_perm = sim2pts(pts, tri, np.real(mesh_new.perm))
 
 # solving once using fem
 # f, _ = fwd.solve(ex_line, tri_perm)
 # f = np.real(f)
 
 # calculate simulated data
 v0 = fwd.solve_eit()
-v1 = fwd.solve_eit(perm=mesh_new["perm"], init=True)
+v1 = fwd.solve_eit(perm=mesh_new.perm, init=True)
 
 """ 3. JAC solver """
 # number of stimulation lines/patterns
-eit = jac.JAC(mesh_obj, protocol)
+eit = jac.JAC(mesh_obj, protocol_obj)
 eit.setup(p=0.50, lamb=1e-3, method="kotre")
 ds = eit.solve(v1, v0, normalize=False)
 node_ds = sim2pts(pts, tri, np.real(ds))

examples/eit_dynamic_stack.py

@@ -2,55 +2,50 @@
 """ demo using stacked ex_mat (the devil is in the details) """
 # 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
+from __future__ import absolute_import, division, print_function
 
-import numpy as np
 import matplotlib.pyplot as plt
+import numpy as np
+import pyeit.eit.jac as jac
+import pyeit.eit.protocol as protocol
 
 # pyEIT 2D algorithm modules
 import pyeit.mesh as mesh
 from pyeit.eit.fem import EITForward
-from pyeit.eit.utils import eit_scan_lines
-
 from pyeit.mesh.shape import thorax
-import pyeit.eit.jac as jac
+from pyeit.mesh.wrapper import PyEITAnomaly_Circle
 
 """ 1. setup """
-use_thorax_model = False
-if use_thorax_model:
-    # Mesh shape is specified with fd parameter in the instantiation, e.g : fd=thorax , Default :fd=circle
-    mesh_obj = mesh.create(16, h0=0.1, fd=thorax)
+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.1, fd=thorax)
 else:
-    mesh_obj = mesh.layer_circle()
+    mesh_obj = mesh.create(n_el, h0=0.1)
 
 # test function for altering the permittivity in mesh
-anomaly = [{"x": 0.4, "y": 0.4, "d": 0.2, "perm": 100}]
+anomaly = PyEITAnomaly_Circle(center=[0.4, 0.4], r=0.2, perm=100.0)
 mesh_new = mesh.set_perm(mesh_obj, anomaly=anomaly, background=1.0)
 
 """ 2. calculate simulated data using stack ex_mat """
-el_dist, step = 7, 1
-n_el = len(mesh_obj["el_pos"])
-ex_mat1 = eit_scan_lines(n_el, el_dist)
-# TODO: a combinational el_dist of 1 and other value should also work.
-ex_mat2 = eit_scan_lines(n_el, 3)
-ex_mat = np.vstack([ex_mat1, ex_mat2])
-protocol = {"ex_mat": ex_mat, "step": step, "parser": "std"}
+protocol_obj = protocol.create(n_el, dist_exc=[7, 3], step_meas=1, parser_meas="std")
 
 # forward solver
-fwd = EITForward(mesh_obj, protocol)
+fwd = EITForward(mesh_obj, protocol_obj)
 v0 = fwd.solve_eit()
-v1 = fwd.solve_eit(perm=mesh_new["perm"], init=True)
+v1 = fwd.solve_eit(perm=mesh_new.perm, init=True)
 
 """ 3. solving using dynamic EIT """
 # number of stimulation lines/patterns
-eit = jac.JAC(mesh_obj, protocol)
+eit = jac.JAC(mesh_obj, protocol_obj)
 eit.setup(p=0.40, lamb=1e-3, method="kotre")
 ds = eit.solve(v1, v0, normalize=False)
 
 # extract node, element, alpha
-pts = mesh_obj["node"]
-tri = mesh_obj["element"]
-delta_perm = mesh_new["perm"] - mesh_obj["perm"]
+pts = mesh_obj.node
+tri = mesh_obj.element
+delta_perm = mesh_new.perm - mesh_obj.perm
 
 # show alpha
 fig, ax = plt.subplots(figsize=(6, 4))