Implement EIT system performance measures. (#94)

* .ply reading
* new figures of merit example
* improvements to figures of merit
* add convenience function for rendering
* bugfix mesh plot
* add tests for figures of merit
* update render to keep image orientation
* update mesh plot to allow mesh created with create function
* figures of merit example for range of targets
* EIT system analysis
* adding drift calculation
* add open EIT format
* add option to calculate snr and detectability in decibels
* add EIT system performance example
* add colorbar function that matches height with corresponding image
* make parse oeit return numpy array
* read oeit reject non uniform lines
* add allantools to environment.yml

Author: acreegan

environment.yml

@@ -12,4 +12,5 @@ dependencies:
   - vispy
   - shapely
   - trimesh
-  - "imageio>=2.19"
+  - "imageio>=2.19"
+  - allantools

examples/eit_system_performance_simulated.py

@@ -0,0 +1,207 @@
+# coding: utf-8
+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
+import pyeit.eit.protocol as protocol
+from pyeit.mesh.wrapper import PyEITAnomaly_Circle
+from pyeit.mesh.external import place_electrodes_equal_spacing
+from numpy.random import default_rng
+from pyeit.quality.eit_system import (
+    calc_signal_to_noise_ratio,
+    calc_accuracy,
+    calc_drift,
+    calc_detectability,
+)
+from pyeit.eit.render import render_2d_mesh
+from pyeit.visual.plot import colorbar
+
+
+def main():
+    # Configuration
+    # ------------------------------------------------------------------------------------------------------------------
+    n_el = 16
+    render_resolution = (64, 64)
+    background_value = 1
+    anomaly_value = 2
+    noise_magnitude = 2e-4
+    drift_rate = 2e-7  # Per frame
+
+    n_background_measurements = 10
+    n_drift_measurements = 1800
+    measurement_frequency = 1
+
+    detectability_r = 0.5
+    distinguishability_r = 0.35
+
+    # Initialization
+    # ------------------------------------------------------------------------------------------------------------------
+    drift_period_hours = n_drift_measurements / (60 * 60 * measurement_frequency)
+    conductive_target = True if anomaly_value - background_value > 0 else False
+    det_center_range = np.arange(0, 0.667, 0.067)
+    dist_distance_range = np.arange(0.35, 1.2, 0.09)
+    rng = default_rng(0)
+
+    # Problem setup
+    # ------------------------------------------------------------------------------------------------------------------
+    sim_mesh = mesh.create(n_el, h0=0.05)
+    electrode_nodes = place_electrodes_equal_spacing(sim_mesh, n_electrodes=n_el)
+    sim_mesh.el_pos = np.array(electrode_nodes)
+    protocol_obj = protocol.create(n_el, dist_exc=1, step_meas=1, parser_meas="std")
+    fwd = EITForward(sim_mesh, protocol_obj)
+
+    recon_mesh = mesh.create(n_el, h0=0.1)
+    electrode_nodes = place_electrodes_equal_spacing(recon_mesh, n_electrodes=n_el)
+    recon_mesh.el_pos = np.array(electrode_nodes)
+    eit = jac.JAC(recon_mesh, protocol_obj)
+    eit.setup(
+        p=0.5, lamb=0.03, method="kotre", perm=background_value, jac_normalized=True
+    )
+
+    # Simulate background
+    # ------------------------------------------------------------------------------------------------------------------
+    v0 = fwd.solve_eit(perm=background_value)
+    d1 = np.array(range(n_drift_measurements)) * drift_rate  # Create drift
+    n1 = noise_magnitude * rng.standard_normal(
+        (n_drift_measurements, len(v0))
+    )  # Create noise
+
+    v0_dn = np.tile(v0, (n_drift_measurements, 1)) + np.tile(d1, (len(v0), 1)).T + n1
+
+    # Calculate background performance measures
+    # ------------------------------------------------------------------------------------------------------------------
+    snr = calc_signal_to_noise_ratio(v0_dn[:n_background_measurements], method="db")
+    accuracy = calc_accuracy(v0_dn[:n_background_measurements], v0, method="EIDORS")
+    t2, adevs = calc_drift(v0_dn, method="Allan")
+
+    drifts_delta = calc_drift(v0_dn, sample_period=10, method="Delta")
+    start = np.average(v0_dn[0:10], axis=0)
+    drifts_percent = 100 * drifts_delta / start
+
+    # Simulate detectability test
+    # ------------------------------------------------------------------------------------------------------------------
+    detectabilities = []
+    detectability_renders = []
+    for c in det_center_range:
+        anomaly = PyEITAnomaly_Circle(
+            center=[c, 0], r=detectability_r, perm=anomaly_value
+        )
+        sim_mesh_new = mesh.set_perm(
+            sim_mesh, anomaly=anomaly, background=background_value
+        )
+        v1 = fwd.solve_eit(perm=sim_mesh_new.perm)
+        n = noise_magnitude * rng.standard_normal(len(v1))
+        v1_n = v1 + n
+        ds = eit.solve(v1_n, v0_dn[0], normalize=True)
+        solution = np.real(ds)
+        image = render_2d_mesh(recon_mesh, solution, resolution=render_resolution)
+        detectability = calc_detectability(
+            image, conductive_target=conductive_target, method="db"
+        )
+        detectabilities.append(detectability)
+        detectability_renders.append(image)
+
+    # Simulate distinguishability test
+    # ------------------------------------------------------------------------------------------------------------------
+    anomaly = PyEITAnomaly_Circle(center=[0, 0], r=detectability_r, perm=anomaly_value)
+    sim_mesh_new = mesh.set_perm(sim_mesh, anomaly=anomaly, background=background_value)
+    dist_v0 = fwd.solve_eit(perm=sim_mesh_new.perm)
+    dist_v0n = dist_v0 + noise_magnitude * rng.standard_normal(len(dist_v0))
+
+    distinguishabilities = []
+    distinguishabilitiy_renders = []
+    for d in dist_distance_range:
+        a1 = PyEITAnomaly_Circle(
+            center=[d / 2, 0], r=distinguishability_r, perm=anomaly_value
+        )
+        a2 = PyEITAnomaly_Circle(
+            center=[-d / 2, 0], r=distinguishability_r, perm=anomaly_value
+        )
+        sim_mesh_new = mesh.set_perm(
+            sim_mesh, anomaly=[a1, a2], background=background_value
+        )
+        v1 = fwd.solve_eit(perm=sim_mesh_new.perm)
+        v1_n = v1 + noise_magnitude * rng.standard_normal(len(v1))
+        ds = eit.solve(v1_n, dist_v0n, normalize=True)
+        solution = np.real(ds)
+        image = render_2d_mesh(recon_mesh, solution, resolution=render_resolution)
+        # Distinguishability is detectability but with a target as the background.
+        distinguishability = calc_detectability(
+            image, conductive_target=conductive_target, method="db"
+        )
+        distinguishabilities.append(distinguishability)
+        distinguishabilitiy_renders.append(image)
+
+    # Plot results
+    # ------------------------------------------------------------------------------------------------------------------
+    fig, axs = plt.subplots(2, 2)
+
+    axs[0, 0].plot(snr)
+    axs[0, 0].set_xlabel("Channel Number")
+    axs[0, 0].set_ylabel("Signal to Noise Ratio\n(dB)")
+    axs[0, 0].title.set_text(f"Signal to Noise Ratio for {len(snr)} channels")
+
+    axs[0, 1].plot(accuracy)
+    axs[0, 1].set_xlabel("Channel Number")
+    axs[0, 1].set_ylabel("Accuracy")
+    axs[0, 1].title.set_text(f"Accuracy for {len(snr)} channels")
+
+    axs[1, 0].set_xlabel("Averaging Window (s)")
+    axs[1, 0].set_ylabel("Allan Deviation")
+    axs[1, 0].title.set_text(f"Allan Deviation for {len(snr)} channels")
+    for adev in adevs:
+        axs[1, 0].plot(t2, adev)
+
+    axs[1, 1].plot(drifts_percent)
+    axs[1, 1].title.set_text(
+        f"Drift percentage on all channels.\nDrift period (hours): {drift_period_hours}"
+    )
+    axs[1, 1].set_xlabel("Channel number")
+    axs[1, 1].set_ylabel("Drift (% of starting value)")
+
+    fig.tight_layout()
+    fig.set_size_inches((10, 6))
+
+    fig, axs = plt.subplots(1, 2)
+    axs[0].plot(det_center_range, detectabilities, ".-", label="-x axis")
+    axs[0].legend()
+    axs[0].set_xlabel("Target position (radius fraction)")
+    axs[0].set_ylabel("Detectability (dB)")
+    axs[0].title.set_text("Detectability vs radial position")
+
+    axs[1].plot(dist_distance_range, distinguishabilities)
+    axs[1].set_xlabel("Separation distance (radius fraction)")
+    axs[1].set_ylabel("Distinguishability (dB)")
+    axs[1].title.set_text("Distinguishability vs separation distance")
+
+    fig.set_size_inches((10, 4))
+    fig.tight_layout()
+
+    fig, axs = plt.subplots(1, len(det_center_range))
+    for i, c in enumerate(det_center_range):
+        axs[i].imshow(detectability_renders[i])
+        axs[i].xaxis.set_ticks([])
+        axs[i].yaxis.set_ticks([])
+
+    fig.set_size_inches((14, 2))
+    fig.suptitle("Detectability Renders")
+    fig.tight_layout()
+
+    fig, axs = plt.subplots(1, len(dist_distance_range))
+    for i, d in enumerate(dist_distance_range):
+        img = axs[i].imshow(distinguishabilitiy_renders[i])
+        axs[i].xaxis.set_ticks([])
+        axs[i].yaxis.set_ticks([])
+        colorbar(img)
+
+    fig.set_size_inches((18, 2))
+    fig.suptitle("Distinguishability Renders")
+    fig.tight_layout()
+    plt.show()
+
+
+if __name__ == "__main__":
+    main()

examples/figures_of_merit_range.py

@@ -88,15 +88,14 @@ def main():
             axs[i],
             solution,
             recon_mesh,
-            ax_kwargs={"title": f"Target pos: {plot_c[i]:.2f}/r"},
+            ax_kwargs={"title": f"Target Pos: {plot_c[i]:.2f}/r"},
         )
 
-    fig.set_size_inches(15, 2)
+    fig.set_size_inches(12, 2)
     fig.tight_layout()
 
     figs_list = np.array(figs_list)
-    fig, axs = plt.subplots(5, 1, sharex=True)
-    axs[4].set_xlabel("Target pos/r")
+    fig, axs = plt.subplots(1, 5)
     titles = [
         "Average Amplitude",
         "Position Error",
@@ -106,9 +105,12 @@ def main():
     ]
     for i in range(5):
         axs[i].plot(c_range, figs_list[:, i])
-        axs[i].set_title(titles[i], size="small")
+        axs[i].set_title(f"{titles[i]}\nvs Target Pos")
+        axs[i].set_xlabel("Target Pos/r")
+        axs[i].set_ylabel(titles[i])
 
-    plt.tight_layout()
+    fig.set_size_inches(15, 3)
+    fig.tight_layout()
     plt.show()
 
 

examples/figures_of_merit_single.py

@@ -69,20 +69,20 @@ def main():
     # Print figures of merit
     print("")
     print(f"Amplitude: Average pixel value in reconstruction image is {figs[0]:.4f}")
-    print(f"Position Error: {100*figs[1]:.2f}% of widest axis")
+    print(f"Position Error: {100 * figs[1]:.2f}% of widest axis")
     print(
-        f"Resolution: Reconstructed point radius {100*figs[2]:.2f}% of image equivalent radius"
+        f"Resolution: Reconstructed point radius {100 * figs[2]:.2f}% of image equivalent radius"
     )
     print(
-        f"Shape Deformation: {100*figs[3]:.2f}% of pixels in the thresholded image are outside the equivalent circle"
+        f"Shape Deformation: {100 * figs[3]:.2f}% of pixels in the thresholded image are outside the equivalent circle"
     )
     print(
-        f"Ringing: Ringing pixel amplitude is  {100*figs[4]:.2f}% of image amplitude in thresholded region"
+        f"Ringing: Ringing pixel amplitude is  {100 * figs[4]:.2f}% of image amplitude in thresholded region"
     )
 
     # Create mesh plots
     fig, axs = plt.subplots(1, 2)
-    create_mesh_plot(axs[0], sim_mesh, ax_kwargs={"title": "Sim mesh"})
+    create_mesh_plot(axs[0], sim_mesh_new, ax_kwargs={"title": "Sim mesh"})
     create_mesh_plot(axs[1], recon_mesh, ax_kwargs={"title": "Recon mesh"})
     fig.set_size_inches(10, 4)
 

pyeit/io/oeit.py

@@ -0,0 +1,35 @@
+import numpy as np
+from scipy import stats
+
+
+def load_oeit_data(file_name):
+    with open(file_name, "r") as f:
+        lines = f.readlines()
+
+    data = []
+    for line in lines:
+        eit = parse_oeit_line(line)
+        if eit is not None:
+            data.append(eit)
+
+    mode_len = stats.mode([len(item) for item in data], keepdims=False)
+    data = [item for item in data if len(item) == mode_len.mode]
+
+    return np.array(data)
+
+
+def parse_oeit_line(line):
+    try:
+        _, data = line.split(":", 1)
+    except (ValueError, AttributeError):
+        return None
+    items = []
+    for item in data.split(","):
+        item = item.strip()
+        if not item:
+            continue
+        try:
+            items.append(float(item))
+        except ValueError:
+            return None
+    return np.array(items)