From 8c4e22880650041e7d30f4ee29afbefde17fac11 Mon Sep 17 00:00:00 2001 From: DavidMetzIMT Date: Fri, 6 May 2022 02:27:18 +0200 Subject: [PATCH] Separate solve, solve_eit, compute_b, and compute_jac (#43) - Docstrings style back to numpy-style, NDArray>> np.ndarray - Move the computation of H in setup() - transform EitBase to an absclass with absmethod "setup" and "_compute_h" - add "_check_solver_is_ready()" which raise SolverNotReadyError if setup() not called before solving for example. - modify msg of _check_solver_is_ready and fix missing call in jac.gn - bp, jac, greit and svd modified to fit ne EitBase - "map" method only in eitBase (same impletation in fac, bp, greit) - add/enhance doctring (numpy-format) in all 4 files and add typehints - spliting computation of potential dist, jac, b (see new methods compute_*) - using those directly in EitSolvers - solve and solve_eit (assuring bck compatibility) with new implementation - the method still return Fwd results (to minimize back compatility ) but contain only the voltages, can be extended with futher computed data e.g. diff_op or others info - fix gn solver jac and v0 were not actualizated, now passing x0 - put old test in unit test frame! now bay calling the script all test are run and a protocols is printed in term - to set the ref electrode a new set method has been added to Forward - solve >> handle ex_line (return f[0, :].ravel()) - jac=np.vstack(jac) also in compute_jac returns jac and not jac_i like the older solve - b=np.vstack(b_matrix) - Suppression of K^-1 for nodes potential - the inversion is still used but only for building the jac - correction of the output of f (ex_mat has not been cheke in the scope of of solve) - add intern _get_boundary _volt method to avoid multple code reapting and add note for v0 after calculation of jac! --- examples/eit_sensitivity2d.py | 8 +- examples/fem_forward2d.py | 2 +- examples/fem_forward3d.py | 2 +- examples/mesh_multi_shell.py | 2 +- examples/softx/figure02.py | 2 +- examples/softx/figure02b.py | 9 +- pyeit/eit/__init__.py | 5 + pyeit/eit/base.py | 253 +++++++++++++------ pyeit/eit/bp.py | 118 ++++++--- pyeit/eit/fem.py | 454 +++++++++++++++++++++++++--------- pyeit/eit/greit.py | 184 +++++++++----- pyeit/eit/interp2d.py | 271 ++++++++++++-------- pyeit/eit/jac.py | 296 ++++++++++++++++------ pyeit/eit/svd.py | 36 ++- pyeit/eit/utils.py | 26 +- pyeit/mesh/__init__.py | 10 +- pyeit/mesh/utils.py | 9 + tests/run_examples.py | 20 ++ tests/test_eit.py | 28 +++ tests/test_fem.py | 234 ++++++++++++------ tests/test_mesh.py | 54 ++++ 21 files changed, 1454 insertions(+), 569 deletions(-) create mode 100644 tests/run_examples.py diff --git a/examples/eit_sensitivity2d.py b/examples/eit_sensitivity2d.py index 592df84..352a8fc 100644 --- a/examples/eit_sensitivity2d.py +++ b/examples/eit_sensitivity2d.py @@ -28,17 +28,17 @@ quality.stats(pts, tri) -def calc_sens(fwd, ex_mat): +def calc_sens(fwd:Forward, ex_mat): """ see Adler2017 on IEEE TBME, pp 5, figure 6, Electrical Impedance Tomography: Tissue Properties to Image Measures """ # solving EIT problem - p = fwd.solve_eit(ex_mat=ex_mat, parser="fmmu") - v0 = p.v + jac = fwd.compute_jac(ex_mat=ex_mat, parser="fmmu") + v0 = fwd.v0 # normalized jacobian (note: normalize affect sensitivity) v0 = v0[:, np.newaxis] - jac = p.jac # / v0 # (normalize or not) + jac = jac # / v0 # (normalize or not) # calculate sensitivity matrix s = np.linalg.norm(jac, axis=0) ae = tri_area(pts, tri) diff --git a/examples/fem_forward2d.py b/examples/fem_forward2d.py index 098a4fa..a0553df 100644 --- a/examples/fem_forward2d.py +++ b/examples/fem_forward2d.py @@ -37,7 +37,7 @@ # calculate simulated data using FEM fwd = Forward(mesh_obj, el_pos) -f, _ = fwd.solve(ex_line, perm=perm) +f = fwd.solve(ex_line, perm=perm) f = np.real(f) """ 2. plot """ diff --git a/examples/fem_forward3d.py b/examples/fem_forward3d.py index 4f8a167..186f0f6 100644 --- a/examples/fem_forward3d.py +++ b/examples/fem_forward3d.py @@ -40,7 +40,7 @@ node_perm = sim2pts(pts, tri, np.real(tri_perm)) # solving once using fem -f, _ = fwd.solve(ex_line, perm=tri_perm) +f = fwd.solve(ex_line, perm=tri_perm) f = np.real(f) # mplot.tetplot(p, t, edge_color=(0.2, 0.2, 1.0, 1.0), alpha=0.01) diff --git a/examples/mesh_multi_shell.py b/examples/mesh_multi_shell.py index 294decd..84c8553 100644 --- a/examples/mesh_multi_shell.py +++ b/examples/mesh_multi_shell.py @@ -70,7 +70,7 @@ ex_line = ex_mat[0].ravel() # solving once using fem -f, _ = fwd.solve(ex_line, perm=tri_perm) +f = fwd.solve(ex_line, perm=tri_perm) f = np.real(f) vf = np.linspace(min(f), max(f), 32) diff --git a/examples/softx/figure02.py b/examples/softx/figure02.py index 80622e8..a5604f7 100644 --- a/examples/softx/figure02.py +++ b/examples/softx/figure02.py @@ -31,7 +31,7 @@ # calculate simulated data using FEM fwd = Forward(mesh_obj, el_pos) -f, _ = fwd.solve(ex_line, perm=perm) +f = fwd.solve(ex_line, perm=perm) f = np.real(f) """ 2. plot """ diff --git a/examples/softx/figure02b.py b/examples/softx/figure02b.py index 79adb34..1376b9a 100644 --- a/examples/softx/figure02b.py +++ b/examples/softx/figure02b.py @@ -27,17 +27,18 @@ quality.stats(pts, tri) -def calc_sens(fwd, ex_mat): +def calc_sens(fwd:Forward, ex_mat): """ see Adler2017 on IEEE TBME, pp 5, figure 6, Electrical Impedance Tomography: Tissue Properties to Image Measures """ # solving EIT problem - p = fwd.solve_eit(ex_mat=ex_mat, parser="fmmu") - v0 = p.v + # p = fwd.solve_eit(ex_mat=ex_mat, parser="fmmu") + jac= fwd.compute_jac(ex_mat=ex_mat, parser="fmmu") + v0 = fwd.v0 # normalized jacobian (note: normalize affect sensitivity) v0 = v0[:, np.newaxis] - jac = p.jac # / v0 + jac = jac # / v0 # calculate sensitivity matrix s = np.linalg.norm(jac, axis=0) ae = tri_area(pts, tri) diff --git a/pyeit/eit/__init__.py b/pyeit/eit/__init__.py index bbea564..0faad82 100644 --- a/pyeit/eit/__init__.py +++ b/pyeit/eit/__init__.py @@ -11,3 +11,8 @@ - utils: EIT related helper function - interp2d: Spatial interpolation for EIT """ +from .bp import BP + +__all__ = [ + "BP", +] diff --git a/pyeit/eit/base.py b/pyeit/eit/base.py index b3db2ba..6d66348 100644 --- a/pyeit/eit/base.py +++ b/pyeit/eit/base.py @@ -8,6 +8,8 @@ # 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 abc import ABC, abstractmethod +from typing import Tuple, Union import numpy as np @@ -15,39 +17,49 @@ from .utils import eit_scan_lines -class EitBase: +class SolverNotReadyError(BaseException): + """Is raised if solver.setup() not called before using solver""" + + +class EitBase(ABC): """ - A base EIT solver. + Base EIT solver. """ def __init__( self, - mesh, - el_pos, - ex_mat=None, - step=1, - perm=None, - jac_normalized=False, - parser="std", + mesh: dict, + el_pos: np.ndarray, + ex_mat: np.ndarray = None, + step: int = 1, + perm: np.ndarray = None, + jac_normalized: bool = False, + parser: Tuple[str, list[str]] = "std", **kwargs, - ): + ) -> None: """ + An EIT solver. + + WARNING: Before using it run solver.setup() get set the solver ready! + Parameters ---------- - mesh: dict + mesh : dict mesh structure - el_pos: array_like + el_pos : np.ndarray position (numbering) of electrodes - ex_mat: array_like, optional (default: opposition) - 2D array, each row is one stimulation pattern/line - step: int, optional - measurement method - perm: array_like, optional - initial permittivity in generating Jacobian - jac_normalized: Boolean (default is False) - normalize the jacobian using f0 computed from input perm - parser: str, optional, default is 'std' - parsing the format of each frame in measurement/file + ex_mat : np.ndarray, optional + 2D array, each row is one stimulation pattern/line, by default None + step : int, optional + measurement method, by default 1 + perm : np.ndarray, optional + initial permittivity in generating Jacobian, by default None + jac_normalized : bool, optional + normalize the jacobian using f0 computed from input perm, by + default False + parser : Tuple[str, list[str]], optional + parsing the format of each frame in measurement/file, by + default "std" Notes ----- @@ -60,8 +72,7 @@ def __init__( perm = mesh["perm"] # build forward solver - fwd = Forward(mesh, el_pos) - self.fwd = fwd + self.fwd = Forward(mesh, el_pos) # solving mesh structure self.mesh = mesh @@ -75,74 +86,173 @@ def __init__( self.parser = parser # user may specify a scalar for uniform permittivity - if np.size(perm) == 1: - self.perm = perm * np.ones(self.el_num) - else: - self.perm = perm - + self.perm = perm * np.ones(self.el_num) if np.size(perm) == 1 else perm # solving configurations self.ex_mat = ex_mat self.step = step + self.jac_normalized = jac_normalized + + # initialize other parameters + self.params = None + self.xg = None + self.yg = None + self.mask = None + # user must run solver.setup() manually to get correct H + self.H = None + self.is_ready = False - # solving Jacobian using uniform sigma distribution - res = fwd.solve_eit(ex_mat, step=step, perm=self.perm, parser=self.parser) - self.J, self.v0, self.B = res.jac, res.v, res.b_matrix + @abstractmethod + def setup(self) -> None: + """ + Setup EIT solver - # Jacobian normalization: divide each row of J (J[i]) by abs(v0[i]) - if jac_normalized: - self.J = self.J / np.abs(self.v0[:, None]) + 1. memory parameters in self.params + 2. compute some other stuff needed for 3. + 3. compute self.H used for solving inv problem by using + >> self.H=self._compute_h() + 4. set flag self.is_ready to `True` + """ - # initialize other parameters - self.params = {} - self.xg = [] - self.yg = [] - self.mask = [] + @abstractmethod + def _compute_h(self) -> np.ndarray: + """ + Compute H matrix for solving inv problem - # mapping matrix - self.H = self.B.T - # warning: user must run solver.setup() manually to get correct H - # self.setup() + To be used in self.setup() + >> self.H=self._compute_h() - def setup(self): - """setup EIT solver""" - raise NotImplementedError + Returns + ------- + np.ndarray + H matrix + """ - def solve(self, v1, v0, normalize=False, log_scale=False): + def solve( + self, + v1: np.ndarray, + v0: np.ndarray, + normalize: bool = False, + log_scale: bool = False, + ) -> np.ndarray: """ - dynamic imaging (conductivities imaging) + Dynamic imaging (conductivities imaging) Parameters ---------- - v1: NDArray + v1: np.ndarray current frame - v0: NDArray + v0: np.ndarray referenced frame, d = H(v1 - v0) normalize: Bool, optional - true for conducting normalization + true for conducting normalization, by default False log_scale: Bool, optional - remap reconstructions in log scale + remap reconstructions in log scale,by default False + + Raises + ------ + SolverNotReadyError: raised if solver not ready + (see self._check_solver_is_ready()) Returns ------- - ds: NDArray - complex-valued NDArray, changes of conductivities + ds: np.ndarray + complex-valued np.ndarray, changes of conductivities """ - if normalize: - dv = self.normalize(v1, v0) - else: - dv = v1 - v0 - + self._check_solver_is_ready() + dv = self._normalize(v1, v0) if normalize else v1 - v0 ds = -np.dot(self.H, dv.transpose()) # s = -Hv if log_scale: ds = np.exp(ds) - 1.0 - return ds - def map(self): - """simple mat using projection matrix""" - raise NotImplementedError + def map(self, dv: np.ndarray) -> np.ndarray: + """ + (NOT USED, Deprecated?) simple mat using projection matrix + + return -H*dv, dv should be normalized. + + Parameters + ---------- + dv : np.ndarray + voltage measurement frame difference (reference frame - current frame) + + Raises + ------ + SolverNotReadyError + raised if solver not ready (see self._check_solver_is_ready()) + + Returns + ------- + np.ndarray + -H*dv + """ + self._check_solver_is_ready() + return -np.dot(self.H, dv.transpose()) + + def _compute_jac_matrix( + self, perm: Union[int, float, np.ndarray] = None, allow_jac_norm: bool = True + ) -> np.ndarray: + """ + Return Jacobian matrix correspoding to the fwd + + Parameters + ---------- + perm : Union[int, float, np.ndarray], optional + permittivity, by default None + (see Foward._get_perm for more details, in fem.py) + allow_jac_norm : bool, optional + flag allowing the Jacobian to be normalized according to + `self.jac_normalized` intern flag, by default True + (e.g. for `jac.gn` or `greit` no normalization is needed!) + + Returns + ------- + np.ndarray + Jacobian matrix + + Notes + ----- + - initial boundary voltage meas. extimation v0 can be accessed + after computation through call self.fwd.v0 + """ + + return self.fwd.compute_jac( + ex_mat=self.ex_mat, + step=self.step, + perm=perm if perm is not None else self.perm, + parser=self.parser, + normalize=self.jac_normalized and allow_jac_norm, + ) + + def _compute_b_matrix(self) -> np.ndarray: + """ + Return BP matrix correspoding to the fwd + + Returns + ------- + np.ndarray + BP matrix + """ + return self.fwd.compute_b_matrix( + ex_mat=self.ex_mat, step=self.step, perm=self.perm, parser=self.parser + ) + + def _check_solver_is_ready(self) -> None: + """ + Check if solver is ready for solving - def normalize(self, v1, v0): + Addtionaly test also if self.H not `None` + + Raises + ------ + SolverNotReadyError + raised if solver not ready + """ + if not self.is_ready or self.H is None: + msg = "User must first run solver.setup() before using solver for solving purpose" + raise SolverNotReadyError(msg) + + def _normalize(self, v1: np.ndarray, v0: np.ndarray) -> np.ndarray: """ Normalize current frame using the amplitude of the reference frame. Boundary measurements v are complex-valued, we can use the real part of v, @@ -150,11 +260,14 @@ def normalize(self, v1, v0): Parameters ---------- - v1: NDArray + v1 : np.ndarray current frame, can be a Nx192 matrix where N is the number of frames - v0: NDArray + v0 : np.ndarray referenced frame, which is a row vector - """ - dv = (v1 - v0) / np.abs(v0) - return dv + Returns + ------- + np.ndarray + Normalized current frame difference dv + """ + return (v1 - v0) / np.abs(v0) diff --git a/pyeit/eit/bp.py b/pyeit/eit/bp.py index 3e0f198..03f19d7 100644 --- a/pyeit/eit/bp.py +++ b/pyeit/eit/bp.py @@ -12,62 +12,118 @@ class BP(EitBase): """A naive inversion of (Euclidean) back projection.""" - def setup(self, weight="none"): - """setup BP""" + def setup(self, weight: str = "none") -> None: + """ + Setup BP solver + + Parameters + ---------- + weight : str, optional + BP parameter, by default "none" + """ self.params = {"weight": weight} # build the weighting matrix # BP: in node imaging, H is the smear matrix (transpose of B) - if weight == "simple": - weights = self.simple_weight(self.B.shape[0]) - wb_mat = weights * self.B - self.H = wb_mat.T + self.B = self._compute_b_matrix() + self.H = self._compute_h(b_matrix=self.B) + self.is_ready = True - def normalize(self, v1, v0): + def _compute_h(self, b_matrix: np.ndarray) -> np.ndarray: """ - redefine normalize for BP (without amplitude normalization) using - only the sign of v0.real. [experimental] + Compute H matrix for BP solver + + Parameters + ---------- + b_matrix : np.ndarray + BP matrix + + Returns + ------- + np.ndarray + H matrix + """ + if self.params["weight"] == "simple": + weights = self._simple_weight(b_matrix.shape[0]) + b_matrix = weights * b_matrix + return b_matrix.T + + # -------------------------------------------------------------------------- + # Special methods for BP + # -------------------------------------------------------------------------- + + def solve_gs(self, v1: np.ndarray, v0: np.ndarray) -> np.ndarray: """ - dv = (v1 - v0) / np.sign(v0.real) - return dv + Solving using gram-schmidt orthogonalization - def map(self, dv): - """return -H*dv, dv should be normalized.""" - x = -dv - return np.dot(self.H, x.transpose()) + Parameters + ---------- + v1: np.ndarray + current frame + v0: np.ndarray + referenced frame + + Raises + ------ + SolverNotReadyError + raised if solver not ready (see self._check_solver_is_ready()) - def solve_gs(self, v1, v0): - """solving using gram-schmidt orthogonalization""" + Returns + ------- + np.ndarray + complex-valued np.ndarray, changes of conductivities + """ + self._check_solver_is_ready() a = np.dot(v1, v0) / np.dot(v0, v0) vn = -(v1 - a * v0) / np.sign(v0.real) - ds = np.dot(self.H, vn.transpose()) - return ds + return np.dot(self.H, vn.transpose()) - def simple_weight(self, num_voltages): + def _normalize(self, v1: np.ndarray, v0: np.ndarray) -> np.ndarray: """ - building weighting matrix : simple, normalize by radius. + redefine normalize for BP (without amplitude normalization) using + only the sign of v0.real. [experimental] - Note - ---- - as in fem.py, we could either smear at, + Normalize current frame using the amplitude of the reference frame. + Boundary measurements v are complex-valued - (1) elements, using the center co-ordinates (x,y) of each element - >> center_e = np.mean(self.pts[self.tri], axis=1) - (2) nodes. + Parameters + ---------- + v1: np.ndarray + current frame + v0: np.ndarray + referenced frame + + Returns + ------- + np.ndarray + Normalized current frame difference dv + """ + return (v1 - v0) / np.sign(v0.real) + + def _simple_weight(self, num_voltages: int) -> np.ndarray: + """ + Build weighting matrix : simple, normalize by radius. Parameters ---------- - num_voltages: int + num_voltages : int number of equal-potential lines Returns ------- - w: NDArray + np.ndarray weighting matrix + + Notes + ----- + as in fem.py, we could either smear at, + + (1) elements, using the center co-ordinates (x,y) of each element + >> center_e = np.mean(self.pts[self.tri], axis=1) + (2) nodes. """ d = np.sqrt(np.sum(self.pts**2, axis=1)) r = np.max(d) w = (1.01 * r - d) / (1.01 * r) # weighting by element-wise multiplication W with B - weights = np.dot(np.ones((num_voltages, 1)), w.reshape(1, -1)) - return weights + return np.dot(np.ones((num_voltages, 1)), w.reshape(1, -1)) diff --git a/pyeit/eit/fem.py b/pyeit/eit/fem.py index c4b5fb7..4fa49d5 100644 --- a/pyeit/eit/fem.py +++ b/pyeit/eit/fem.py @@ -11,27 +11,22 @@ import numpy as np import numpy.linalg as la from scipy import sparse +import scipy.linalg from .utils import eit_scan_lines @dataclass class FwdResult: - """Summarize the PDE results from solving the fwd problem + """Summarize the results from solving the eit fwd problem Attributes ---------- - jac: np.ndarray - number of measures x n_E complex array, the Jacobian; shape(n_exc, n_el, n_tri) v: np.ndarray number of measures x 1 array, simulated boundary measures; shape(n_exc, n_el) - b_matrix: np.ndarray - back-projection mappings (smear matrix); shape(n_exc, n_pts, 1), dtype= bool """ v: np.ndarray # number of measures x 1 array, simulated boundary measures - jac: np.ndarray # number of measures x n_E complex array, the Jacobian - b_matrix: np.ndarray # back-projection mappings (smear matrix) class Forward: @@ -39,6 +34,8 @@ class Forward: def __init__(self, mesh: dict[str, np.ndarray], el_pos: np.ndarray) -> None: """ + FEM forward solver + A good FEM forward solver should only depend on mesh structure and the position of electrodes @@ -64,97 +61,253 @@ def __init__(self, mesh: dict[str, np.ndarray], el_pos: np.ndarray) -> None: self.el_pos = el_pos # reference electrodes [ref node should not be on electrodes] - ref_el = 0 - while ref_el in self.el_pos: - ref_el += 1 - self.ref = ref_el + self.set_ref_el() # infer dimensions from mesh self.n_pts, self.n_dim = self.pts.shape self.n_tri, self.n_vertices = self.tri.shape self.n_el = el_pos.size + # temporary memory attributes for computation (e.g. jac) + self._r_matrix = None + self._ke = None + + def solve( + self, + ex_mat: np.ndarray = None, + perm: np.ndarray = None, + ) -> np.ndarray: + """ + Calculate and compute the potential distribution (complex-valued) + corresponding to the permittivity distribution `perm ` for all + excitations contained in the excitation pattern `ex_mat` + + Currently, only simple electrode model is supported, + CEM (complete electrode model) is under development. + + Parameters + ---------- + ex_mat : np.ndarray, optional + stimulation/excitation matrix, of shape (n_exc, 2), by default `None`. + (see _get_ex_mat for more details) + perm : Union[int, float, np.ndarray], optional + permittivity on elements ; shape (n_tri,), by default `None`. + Must be the same size with self.tri_perm + If `None`, `self.tri_perm` will be used + If perm is int or float, uniform permittivity on elements will be used + (see _get_perm for more details) + + Returns + ------- + np.ndarray + potential on nodes ; shape (n_exc, n_pts) + + Notes + ------- + For compatibility with some scripts in /examples a single excitation + line can be passed instead of the whole excitation pattern `ex_mat` + (e.g. [0,7] or np.array([0,7]) or ex_mat[0].ravel). In that case a + simplified version of `f` with shape (n_pts,) + """ + f = self._compute_potential_distribution(ex_mat=ex_mat, perm=perm) + # case ex_line has been passed instead of ex_mat + # we return simplified version of f with shape (n_pts,) + if f.shape[0] == 1: + return f[0, :].ravel() + return f + def solve_eit( self, ex_mat: np.ndarray = None, step: int = 1, - perm: np.ndarray = None, + perm: Union[int, float, np.ndarray] = None, parser: Union[str, list[str]] = None, - **kwargs, ) -> FwdResult: """ - EIT simulation, generate perturbation matrix and forward v + EIT simulation, generate forward v measurement Parameters ---------- - ex_mat: np.ndarray - numLines x n_el array, stimulation matrix - step: int - the configuration of measurement electrodes (default: apposition) - perm: np.ndarray - Mx1 array, initial x0. must be the same size with self.tri_perm - parser: str - see voltage_meter for more details. + ex_mat : np.ndarray, optional + stimulation/excitation matrix, of shape (n_exc, 2), by default `None`. + (see _get_ex_mat for more details) + step: int, optional + the configuration of measurement electrodes, by default 1 (adjacent). + perm : Union[int, float, np.ndarray], optional + permittivity on elements ; shape (n_tri,), by default `None`. + Must be the same size with self.tri_perm + If `None`, `self.tri_perm` will be used + If perm is int or float, uniform permittivity on elements will be used + (see _get_perm for more details) + parser: Union[str, list[str]], optional + see voltage_meter for more details, by default `None`. Returns ------- - fwd_results: FwdResult - PDE results from solving the fwd problem (see `FwdResult`) + FwdResult + Foward results comprising + v: np.ndarray + simulated boundary voltage measurements; shape(n_exc, n_el) + """ + f = self._compute_potential_distribution(ex_mat=ex_mat, perm=perm) + # boundary measurements, subtract_row-voltages on electrodes + diff_op = voltage_meter(ex_mat, n_el=self.n_el, step=step, parser=parser) + return FwdResult(v=self._get_boundary_voltages(f, diff_op)) - Note - ---- - Single PDE results are: - - fwd_results.jac: np.ndarray - number of measures x n_E complex array, the Jacobian - - fwd_results.v: np.ndarray - number of measures x 1 array, simulated boundary measures - - fwd_results.b_matrix: np.ndarray - back-projection mappings (smear matrix) + def _get_boundary_voltages(self, f: np.ndarray, diff_op: np.ndarray) -> np.ndarray: + """ + Compute boundary voltages from potential distribution - (see `FwdResult` for more details) + Parameters + ---------- + f : np.ndarray + potential on nodes ; shape (n_exc, n_pts) + diff_op : np.ndarray + measurements pattern / subtract_row pairs [N, M]; shape (n_exc, n_meas_per_exc, 2) + Returns + ------- + np.ndarray + simulated boundary voltage measurements; shape(n_exc, n_el) + """ + f_el = f[:, self.el_pos] + v = subtract_row(f_el, diff_op) + return np.hstack(v) + def compute_jac( + self, + ex_mat: np.ndarray = None, + step: int = 1, + perm: Union[int, float, np.ndarray] = None, + parser: Union[str, list[str]] = None, + normalize: bool = False, + ) -> np.ndarray: """ - # initialize/extract the scan lines (default: apposition) - if ex_mat is None: - ex_mat = eit_scan_lines(16, 8) + Compute the Jacobian matrix + + Parameters + ---------- + ex_mat : np.ndarray, optional + stimulation/excitation matrix, of shape (n_exc, 2), by default `None`. + (see _get_ex_mat for more details) + step: int, optional + the configuration of measurement electrodes, by default 1 (adjacent). + perm : Union[int, float, np.ndarray], optional + permittivity on elements ; shape (n_tri,), by default `None`. + Must be the same size with self.tri_perm + If `None`, `self.tri_perm` will be used + If perm is int or float, uniform permittivity on elements will be used + (see _get_perm for more details) + parser: Union[str, list[str]], optional + see voltage_meter for more details, by default `None`. + normalize : bool, optional + flag for Jacobian normalization, by default False. + If True the Jacobian is normalized + + Returns + ------- + np.ndarray + Jacobian matrix + + Notes + ----- + - initial boundary voltage meas. extimation v0 can be accessed + after computation through call fwd.v0 + + """ + f = self._compute_potential_distribution( + ex_mat=ex_mat, perm=perm, memory_4_jac=True + ) + + # Build Jacobian matrix column wise (element wise) + # Je = Re*Ke*Ve = (nex3) * (3x3) * (3x1) + jac_i = np.zeros((ex_mat.shape[0], self.n_el, self.n_tri), dtype=perm.dtype) + + r_el = self._r_matrix[self.el_pos] + + def jac_init(jac, k): + for (i, e) in enumerate(self.tri): + jac[:, i] = np.dot(np.dot(r_el[:, e], self._ke[i]), f[k, e]) + return jac + + jac_i = np.array(list(map(jac_init, jac_i, np.arange(ex_mat.shape[0])))) + + self._r_matrix = None # clear memory + self._ke = None # clear memory - # initialize the permittivity on element - if perm is None: - perm0 = self.tri_perm - elif np.isscalar(perm): - perm0 = np.ones(self.n_tri, dtype=np.float) - else: - assert perm.shape == (self.n_tri,) - perm0 = perm - - f, jac_i = self.solve(ex_mat, perm0) - f_el = f[:, self.el_pos] - # boundary measurements, subtract_row-voltages on electrodes diff_op = voltage_meter(ex_mat, n_el=self.n_el, step=step, parser=parser) - v = subtract_row(f_el, diff_op) jac = subtract_row(jac_i, diff_op) + self.v0 = self._get_boundary_voltages(f, diff_op) + jac = np.vstack(jac) + + # Jacobian normalization: divide each row of J (J[i]) by abs(v0[i]) + + return jac / np.abs(self.v0[:, None]) if normalize else jac + + def compute_b_matrix( + self, + ex_mat: np.ndarray = None, + step: int = 1, + perm: Union[int, float, np.ndarray] = None, + parser: Union[str, list[str]] = None, + ) -> np.ndarray: + """ + Compute back-projection mappings (smear matrix) + + Parameters + ---------- + ex_mat : np.ndarray, optional + stimulation/excitation matrix, of shape (n_exc, 2), by default `None`. + (see _get_ex_mat for more details) + step: int, optional + the configuration of measurement electrodes, by default 1 (adjacent). + perm : Union[int, float, np.ndarray], optional + permittivity on elements ; shape (n_tri,), by default `None`. + Must be the same size with self.tri_perm + If `None`, `self.tri_perm` will be used + If perm is int or float, uniform permittivity on elements will be used + (see _get_perm for more details) + parser: Union[str, list[str]], optional + see voltage_meter for more details, by default `None`. + + Returns + ------- + np.ndarray + back-projection mappings (smear matrix); shape(n_exc, n_pts, 1), dtype= bool + """ + f = self._compute_potential_distribution(ex_mat=ex_mat, perm=perm) + f_el = f[:, self.el_pos] # build bp projection matrix # 1. we can either smear at the center of elements, using # >> fe = np.mean(f[:, self.tri], axis=1) # 2. or, simply smear at the nodes using f - b_matrix = smear( - f, f_el, diff_op, new=True - ) # set new to `False` to get computation from ChabaneAmaury - # update output, now you can call p.jac, p.v, p.b_matrix - return FwdResult( - jac=np.vstack(jac), v=np.hstack(v), b_matrix=np.vstack(b_matrix) + diff_op = voltage_meter(ex_mat, n_el=self.n_el, step=step, parser=parser) + # set new to `False` to get smear-computation from ChabaneAmaury + b_matrix = smear(f, f_el, diff_op, new=True) + return np.vstack(b_matrix) + + def set_ref_el(self, val: int = None) -> None: + """ + Set reference electrode node + + Parameters + ---------- + val : int, optional + node number of reference electrode, by default None + + """ + self.ref_el = ( + val if val is not None and val not in self.el_pos else max(self.el_pos) + 1 ) - def solve( - self, ex_mat: np.ndarray, perm: np.ndarray - ) -> tuple[np.ndarray, np.ndarray]: + def _compute_potential_distribution( + self, ex_mat: np.ndarray, perm: np.ndarray, memory_4_jac: bool = False + ) -> np.ndarray: """ Calculate and compute the potential distribution (complex-valued) corresponding to the permittivity distribution `perm ` for all excitations contained in the excitation pattern `ex_mat` - The calculation of Jacobian can be skipped. Currently, only simple electrode model is supported, CEM (complete electrode model) is under development. @@ -163,69 +316,118 @@ def solve( ex_mat: np.ndarray stimulation/excitation matrix ; shape (n_exc, 2) perm: np.ndarray - permittivity on elements (initial) ; shape (n_tri,) + permittivity on elements ; shape (n_tri,) + memory_4_jac : bool, optional + flag to memory r_matrix to self._r_matrix and ke to self._ke, + by default False. Returns ------- - f: np.ndarray + np.ndarray potential on nodes ; shape (n_exc, n_pts) - jac: np.ndarray - Jacobian ; shape (n_exc, n_el, n_tri) - Notes - ------- - For compatibility with some scripts in /examples a single excitation - line can be passed instead of the whole excitation pattern `ex_mat` - (e.g. [0,7] or np.array([0,7]) or ex_mat[0].ravel). In that case a - simplified version of `f` with shape (n_pts,) and of - `jac` with shape (n_el, n_tri) are returned """ - # case ex_line has been passed instead of ex_mat - if isinstance(ex_mat, list) and len(ex_mat) == 2: - ex_mat = np.array([ex_mat]).reshape((1, 2)) # build a 2D array - elif isinstance(ex_mat, np.ndarray) and ex_mat.ndim == 1: - ex_mat = ex_mat.reshape((-1, 2)) - elif ( - isinstance(ex_mat, np.ndarray) and ex_mat.ndim == 2 and ex_mat.shape[1] == 2 - ): - pass # correct value - else: - raise TypeError( - f"Wrong value of {ex_mat=} expected an ndarray ; shape (n_exc, 2)" - ) + ex_mat = self._get_ex_mat(ex_mat) # check/init stimulation + perm = self._get_perm(perm) # check/init permitivity # 1. calculate local stiffness matrix (on each element) ke = calculate_ke(self.pts, self.tri) - # 2. assemble to global K - kg = assemble(ke, self.tri, perm, self.n_pts, ref=self.ref) + kg = assemble(ke, self.tri, perm, self.n_pts, ref=self.ref_el) - # 3. calculate electrode impedance matrix R = K^{-1} - r_matrix = la.inv(kg) - r_el = r_matrix[self.el_pos] + if memory_4_jac: + # save + # 3. calculate electrode impedance matrix R = K^{-1} + self._r_matrix = la.inv(kg) + self._ke = ke # 4. solving nodes potential using boundary conditions b = self._natural_boundary(ex_mat) - f = np.dot(r_matrix, b[:, None]).T.reshape(b.shape[:-1]) + return ( + scipy.linalg.solve(kg, b.swapaxes(0, 1)) + .swapaxes(0, 1) + .reshape(b.shape[0:2]) + ) - # 5. build Jacobian matrix column wise (element wise) - # Je = Re*Ke*Ve = (nex3) * (3x3) * (3x1) - jac = np.zeros((ex_mat.shape[0], self.n_el, self.n_tri), dtype=perm.dtype) + def _get_perm(self, perm: Union[int, float, np.ndarray] = None) -> np.ndarray: + """ + Check/init the permittivity on element - def jac_init(jac, k): - for (i, e) in enumerate(self.tri): - jac[:, i] = np.dot(np.dot(r_el[:, e], ke[i]), f[k, e]) - return jac + Parameters + ---------- + perm : Union[int, float, np.ndarray], optional + permittivity on elements ; shape (n_tri,), by default `None`. + Must be the same size with self.tri_perm + If `None`, `self.tri_perm` will be used + If perm is int or float, uniform permittivity on elements will be used - jac = np.array(list(map(jac_init, jac, np.arange(ex_mat.shape[0])))) + Returns + ------- + np.ndarray + permittivity on elements ; shape (n_tri,) - # case ex_line has been passed instead of ex_mat - # we return simplified version of f with shape (n_pts,) and jac with shape (ne, n_tri) - if ex_mat.shape[0] == 1: - return f[0, :].ravel(), jac[0, :, :] + Raises + ------ + TypeError + raised if perm is not ndarray and of shape (n_tri,) + """ - return f, jac + if perm is None: + return self.tri_perm + elif isinstance(perm, (int, float)): + return np.ones(self.n_tri, dtype=float) * perm + + if not isinstance(perm, np.ndarray) or perm.shape != (self.n_tri,): + raise TypeError( + f"Wrong type/shape of {perm=}, expected an ndarray; shape (n_tri, )" + ) + return perm + + def _get_ex_mat(self, ex_mat: np.ndarray = None) -> np.ndarray: + """ + Check/init stimulation + + Parameters + ---------- + ex_mat : np.ndarray, optional + stimulation/excitation matrix, of shape (n_exc, 2), by default `None`. + If `None` initialize stimulation matrix for 16 electrode and + apposition mode (see function `eit_scan_lines(16, 8)`) + If single stimulation (ex_line) is passed only a list of length 2 + and np.ndarray of size 2 will be treated. + + Returns + ------- + np.ndarray + stimulation matrix + + Raises + ------ + TypeError + Only accept, `None`, list of length 2, np.ndarray of size 2, + or np.ndarray of shape (n_exc, 2) + """ + if ex_mat is None: + # initialize the scan lines for 16 electrodes (default: apposition) + ex_mat = eit_scan_lines(16, 8) + elif isinstance(ex_mat, list) and len(ex_mat) == 2: + # case ex_line has been passed instead of ex_mat + ex_mat = np.array([ex_mat]).reshape((1, 2)) # build a 2D array + elif isinstance(ex_mat, np.ndarray) and ex_mat.size == 2: + # case ex_line np.ndarray has been passed instead of ex_mat + ex_mat = ex_mat.reshape((-1, 2)) + + if ( + not isinstance(ex_mat, np.ndarray) + or ex_mat.ndim != 2 + or ex_mat.shape[1] != 2 + ): + raise TypeError( + f"Wrong shape of {ex_mat=} expected an ndarray ; shape (n_exc, 2)" + ) + + return ex_mat def _natural_boundary(self, ex_mat: np.ndarray) -> np.ndarray: """ @@ -241,7 +443,7 @@ def _natural_boundary(self, ex_mat: np.ndarray) -> np.ndarray: Returns ---------- - b: np.ndarray + np.ndarray Global boundary condition on pts ; shape (n_exc, n_pts, 1) """ drv_a_global = self.el_pos[ex_mat[:, 0]] @@ -259,6 +461,8 @@ def _smear(f: np.ndarray, fb: np.ndarray, pairs: np.ndarray) -> np.ndarray: """ Build smear matrix B for bp for one exitation + used for the smear matrix computation by @ChabaneAmaury + Parameters ---------- f: np.ndarray @@ -293,11 +497,14 @@ def smear( potential on adjacent electrodes; shape (n_exc, n_el) meas_pattern: np.ndarray electrodes numbering pairs; shape (n_exc, n_meas, 2) + new : bool, optional + flag to use new matrices based computation, by default False. + If `False` to smear-computation from ChabaneAmaury is used Returns ------- - B: np.ndarray - back-projection matrix; shape (n_exc, n_meas, n_pts), dtype= bool + np.ndarray + back-projection (smear) matrix; shape (n_exc, n_meas, n_pts), dtype= bool """ if new: # new implementation not much faster! :( @@ -341,8 +548,8 @@ def subtract_row(v: np.ndarray, meas_pattern: np.ndarray) -> np.ndarray: Returns ------- - v_diff: np.ndarray - difference measurements + np.ndarray + difference measurements v_diff """ if v.shape[:1] != meas_pattern.shape[:1]: @@ -381,13 +588,14 @@ def voltage_meter( Parameters ---------- - ex_mat: np.ndarray + ex_mat : np.ndarray Nx2 array, [positive electrode, negative electrode]. ; shape (n_exc, 2) - n_el: int - number of total electrodes. - step: int - measurement method (two adjacent electrodes are used for measuring). - parser: str | list[str] + n_el : int, optional + number of total electrodes, by default 16 + step : int, optional + measurement method (two adjacent electrodes are used for measuring), by default 1 (adjacent) + parser : Union[str, list[str]], optional + parsing the format of each frame in measurement/file, by default None if parser contains 'fmmu', or 'rotate_meas' then data are trimmed, boundary voltage measurements are re-indexed and rotated, start from the positive stimulus electrode start index 'A'. @@ -399,8 +607,8 @@ def voltage_meter( Returns ------- - meas_pattern : np.ndarray - subtract_row pairs [N, M]; shape (n_exc, n_meas_per_exc, 2) + np.ndarray + measurements pattern / subtract_row pairs [N, M]; shape (n_exc, n_meas_per_exc, 2) """ # local node ex_mat = ex_mat.astype(int) @@ -445,12 +653,12 @@ def assemble( n_tri x 1 conductivities on elements n_pts: int number of nodes - ref: int - reference electrode + ref: int, optional + reference electrode, by default 0 Returns ------- - K: np.ndarray + np.ndarray NxN array of complex stiffness matrix Notes @@ -506,7 +714,7 @@ def calculate_ke(pts: np.ndarray, tri: np.ndarray) -> np.ndarray: Returns ------- - ke_array: np.ndarray + np.ndarray n_tri x (n_dim x n_dim) 3d matrix """ n_tri, n_vertices = tri.shape @@ -536,7 +744,7 @@ def calculate_ke(pts: np.ndarray, tri: np.ndarray) -> np.ndarray: def _k_triangle(xy: np.ndarray) -> np.ndarray: """ - given a point-matrix of an element, solving for Kij analytically + Given a point-matrix of an element, solving for Kij analytically using barycentric coordinates (simplex coordinates) Parameters @@ -546,7 +754,7 @@ def _k_triangle(xy: np.ndarray) -> np.ndarray: Returns ------- - ke_matrix: np.ndarray + np.ndarray local stiffness matrix """ # edges (vector) of triangles @@ -579,7 +787,7 @@ def _k_tetrahedron(xy: np.ndarray) -> np.ndarray: Returns ------- - ke_matrix: np.ndarray + np.ndarray local stiffness matrix Notes diff --git a/pyeit/eit/greit.py b/pyeit/eit/greit.py index 27ea51a..b1cccb3 100644 --- a/pyeit/eit/greit.py +++ b/pyeit/eit/greit.py @@ -12,6 +12,7 @@ # 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 typing import Tuple import numpy as np import scipy.linalg as la @@ -23,26 +24,40 @@ class GREIT(EitBase): """The GREIT algorithm""" - def setup(self, method="dist", w=None, p=0.20, lamb=1e-2, n=32, s=20.0, ratio=0.1): + def setup( + self, + method: str = "dist", + w: np.ndarray = None, + p: float = 0.20, + lamb: float = 1e-2, + n: int = 32, + s: float = 20.0, + ratio: float = 0.1, + ) -> None: """ - setup GREIT + Setup GREIT solver Parameters ---------- - method: str, optional - 'set' or 'dist' - w: NDArray, optional - weight on each element - p: float, optional - noise covariance - lamb: float - regularization parameters - n: int, optional - grid size - s: float, optional - control the blur - ratio: float, optional - desired ratio + method : str, optional + only 'dist' accepted, by default "dist" + w : np.ndarray, optional + weight on each element, by default None + p : float, optional + noise covariance, by default 0.20 + lamb : float, optional + regularization parameters, by default 1e-2 + n : int, optional + grid size, by default 32 + s : float, optional + control the blur, by default 20.0 + ratio : float, optional + desired ratio, by default 0.1 + + Raises + ------ + ValueError + raised if method != "dist" References ---------- @@ -52,68 +67,119 @@ def setup(self, method="dist", w=None, p=0.20, lamb=1e-2, n=32, s=20.0, ratio=0. "GREIT: a unified approach to 2D linear EIT reconstruction of lung images." Physiological measurement 30.6 (2009): S35. """ + + if method != "dist": + raise ValueError(f"method {method} not supported yet") + # parameters for GREIT projection if w is None: w = np.ones_like(self.mesh["perm"]) self.params = {"w": w, "p": p, "lamb": lamb, "n": n, "s": s, "ratio": ratio} - # action (currently only support 'dist') - if method == "dist": - w_mat, self.xg, self.yg, self.mask = self._build_grid() - self.H = self._build_dist(w_mat) - else: - raise ValueError("method " + method + " not supported yet") - - def map(self, dv): - """return H*v""" - return -np.dot(self.H, dv.transpose()) - - def _build_dist(self, w_mat): - """generate R using distribution method.""" - lamb, p = self.params["lamb"], self.params["p"] - f = self.fwd.solve_eit( - self.ex_mat, - step=self.step, - perm=self.perm, - parser=self.parser, - ) - jac = f.jac + # Build grids and mask + self.xg, self.yg, self.mask = meshgrid(self.pts, n=n) + + w_mat = self._compute_grid_weights(self.xg, self.yg) + self.J = self._compute_jac_matrix() + self.H = self._compute_h(jac=self.J, w_mat=w_mat) + self.is_ready = True + + def _compute_h(self, jac: np.ndarray, w_mat: np.ndarray) -> np.ndarray: + """ + Generate H (or R) using distribution method for GREIT solver + + Args: + jac (np.ndarray): Jacobian matrix + w_mat (np.ndarray): meights matrix + + Returns: + np.ndarray: H + """ + lamb, p = self.params["lamb"], self.params["p"] # E[yy^T], it is more efficient to use left pinv than right pinv j_j_w = np.dot(jac, jac.T) r_mat = np.diag(np.diag(j_j_w) ** p) jac_inv = la.inv(j_j_w + lamb * r_mat) # RM = E[xx^T] / E[yy^T] - h_mat = np.dot(np.dot(w_mat.T, jac.T), jac_inv) - - return h_mat + return np.dot(np.dot(w_mat.T, jac.T), jac_inv) - def _build_grid(self): - """build grids and mask""" - # initialize grids - n = self.params["n"] - xg, yg, mask = meshgrid(self.pts, n=n) - # mapping from values on triangles to values on grids - xy = np.mean(self.pts[self.tri], axis=1) - xyi = np.vstack((xg.flatten(), yg.flatten())).T - # GREIT is using sigmod as weighting function (global) - ratio, s = self.params["ratio"], self.params["s"] - w_mat = weight_sigmod(xy, xyi, ratio=ratio, s=s) - return w_mat, xg, yg, mask + # -------------------------------------------------------------------------- + # Special methods for GREIT + # -------------------------------------------------------------------------- - def get_grid(self): - """get grids and mask""" + def get_grid(self) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: + """ + Return masking grid data + + Raises + ------ + SolverNotReadyError + raised if solver not ready (see self._check_solver_is_ready()) + + Returns + ------- + Tuple[np.ndarray, np.ndarray, np.ndarray] + x grid, y grid and masking data, which denotes nodes outside + 2D mesh + """ + self._check_solver_is_ready() return self.xg, self.yg, self.mask - def mask_value(self, ds, mask_value=0): - """(plot only) mask values on nodes outside 2D mesh.""" + def mask_value( + self, ds: np.ndarray, mask_value: float = 0 + ) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: + """ + Set mask values on nodes outside 2D mesh. (for plot only) + + Parameters + ---------- + ds : np.ndarray + conductivity data on nodes + mask_value : float, optional + mask conductivity value to set on nodes outside 2D mesh, by + default 0 + + Raises + ------ + SolverNotReadyError + raised if solver not ready (see self._check_solver_is_ready()) + + Returns + ------- + Tuple[np.ndarray, np.ndarray, np.ndarray] + x grid, y grid and "masked" conductivity data on nodes + """ + self._check_solver_is_ready() ds[self.mask] = mask_value ds = ds.reshape(self.xg.shape) return self.xg, self.yg, ds + def _compute_grid_weights(self, xg: np.ndarray, yg: np.ndarray) -> np.ndarray: + """ + Compute weights for given grid (xg,yg) + + Parameters + ---------- + xg : np.ndarray + x grid + yg : np.ndarray + y grid + + Returns + ------- + np.ndarray + weights + """ + # mapping from values on triangles to values on grids + xy = np.mean(self.pts[self.tri], axis=1) + xyi = np.vstack((xg.flatten(), yg.flatten())).T + # GREIT is using sigmod as weighting function (global) + ratio, s = self.params["ratio"], self.params["s"] + return weight_sigmod(xy, xyi, ratio=ratio, s=s) + @staticmethod - def build_set(x, y): + def build_set(x: np.ndarray, y: np.ndarray) -> np.ndarray: """generate R from a set of training sets (deprecate).""" # E_w[yy^T] y_y_t = la.inv(np.dot(y, y.transpose())) - h_matrix = np.dot(np.dot(x, y), y_y_t) - return h_matrix + return np.dot(np.dot(x, y), y_y_t) diff --git a/pyeit/eit/interp2d.py b/pyeit/eit/interp2d.py index 2a26a5d..0ad760f 100644 --- a/pyeit/eit/interp2d.py +++ b/pyeit/eit/interp2d.py @@ -4,6 +4,7 @@ # 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 typing import Tuple import numpy as np import scipy.linalg as la @@ -16,37 +17,63 @@ from pyeit.mesh import layer_circle, set_perm -def meshgrid(pts, n=32, ext_ratio=0, gc=False): +def meshgrid( + pts: np.ndarray, n: int = 32, ext_ratio: float = 0.0, gc: bool = False +) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: """ build xg, yg, mask grids from triangles point cloud function for interpolating regular grids Parameters ---------- - pts: NDArray + pts: np.ndarray nx2 array of points (x, y) - el_pos: NDArray (optional) - the location of electrodes (for extract the convex hull of pts) n: int - the number of meshgrid per dimension + the number of meshgrid per dimension, by default 32 ext_ratio: float - extend the boundary of meshgrid by ext_ratio*d + extend the boundary of meshgrid by ext_ratio*d, by default 0.0 gc: bool - grid_correction, offset xgrid and ygrid by half step size + grid_correction, offset xgrid and ygrid by half step size , by default + False + + Returns + ------- + Tuple[np.ndarray, np.ndarray, np.ndarray] + x grid, y grid, mask Notes ----- mask denotes points outside mesh. """ xg, yg = _build_grid(pts, n=n, ext_ratio=ext_ratio, gc=gc) - # pts_edges = pts[el_pos] pts_edges = _hull_points(pts) mask = _build_mask(pts_edges, xg, yg) return xg, yg, mask -def _build_grid(pts, n=32, ext_ratio=0, gc=False): - """generating mesh grids""" +def _build_grid( + pts: np.ndarray, n: int = 32, ext_ratio: float = 0.0, gc: bool = False +) -> Tuple[np.ndarray, np.ndarray]: + """ + Generating mesh grid from triangles point cloud + + Parameters + ---------- + pts: np.ndarray + nx2 array of points (x, y) + n: int + the number of meshgrid per dimension, by default 32 + ext_ratio: float + extend the boundary of meshgrid by ext_ratio*d, by default 0.0 + gc: bool + grid_correction, offset xgrid and ygrid by half step size , by default + False + + Returns + ------- + Tuple[np.ndarray, np.ndarray] + x grid, y grid + """ x, y = pts[:, 0], pts[:, 1] x_min, x_max = min(x), max(x) y_min, y_max = min(y), max(y) @@ -66,29 +93,75 @@ def _build_grid(pts, n=32, ext_ratio=0, gc=False): return xg, yg -def _build_mask(pts_edges, xg, yg): - """find whether meshgrids is interior of mesh""" +def _build_mask(pts_edges: np.ndarray, xg: np.ndarray, yg: np.ndarray) -> np.ndarray: + """ + find whether meshgrids is interior of mesh + + Parameters + ---------- + pts_edges : np.ndarray + points on the edges of the mesh + xg : np.ndarray + x grid + yg : np.ndarray + x grid + + Returns + ------- + np.ndarray + mask (denotes points outside mesh.) + """ # 1. create mask based on meshes points = np.vstack((xg.flatten(), yg.flatten())).T # 2. extract edge points using el_pos path = Path(pts_edges, closed=False) mask = path.contains_points(points) - return ~mask -def _hull_points(pts): - """return the convex hull points""" +def _hull_points(pts: np.ndarray) -> np.ndarray: + """ + return the convex hull points from a point cloud + + Parameters + ---------- + pts: np.ndarray + nx2 array of points (x, y) + + Returns + ------- + np.ndarray + convex hull points (edge points) + """ cv = ConvexHull(pts) hull_nodes = cv.vertices return pts[hull_nodes, :] -def _distance2d(x, y, center="mean"): +def _distance2d( + x: np.ndarray, y: np.ndarray, center: Tuple[str, list] = "mean" +) -> np.ndarray: """ Calculate radius given center. This function can be OPTIMIZED using numba or cython. + + Parameters + ---------- + x : np.ndarray + nx1 array of x coordiante + y : np.ndarray + nx1 array of y coordiante + center : Tuple[str, list], optional + center definition, by default "mean". + If center is `None`, [0,0] will be used. + If center is "mean", [np.mean(x), np.mean(y)] will be used. + If center is list, [center[0], center[1]] will be used. + + Returns + ------- + np.ndarray + nx1 array of distance from center to points (x,y) """ if center is None: xc, yc = 0, 0 @@ -96,16 +169,25 @@ def _distance2d(x, y, center="mean"): xc, yc = np.mean(x), np.mean(y) else: xc, yc = center[0], center[1] + # distance 2d + return np.sqrt((x - xc) ** 2 + (y - yc) ** 2).ravel() - d = np.sqrt((x - xc) ** 2 + (y - yc) ** 2).ravel() - return d - -def _distance_matrix2d(xy, xyi): +def _distance_matrix2d(xy: np.ndarray, xyi: np.ndarray) -> np.ndarray: """ - Description - ----------- (2D only) return element-wise distance matrix (pair-wise) + + Parameters + ---------- + xy : np.ndarray + nx2 array of points (x, y) + xyi : np.ndarray + points pairs + + Returns + ------- + np.ndarray + distance matrix between pairwise observations """ # Make a distance matrix between pairwise observations # Note: from @@ -117,27 +199,27 @@ def _distance_matrix2d(xy, xyi): return np.hypot(d0, d1) -def weight_sigmod(xy, xyi, ratio=0.05, s=20.0): +def weight_sigmod( + xy: np.ndarray, xyi: np.ndarray, ratio: float = 0.05, s: float = 20.0 +) -> np.ndarray: """ - Description - ----------- (2D only) local weight/interpolate by sigmod function (GREIT3D) Parameters ---------- - xy: NDArray + xy: np.ndarray (x, y) of values - xyi: NDArray + xyi: np.ndarray (xi, yi) of interpolated locations ratio: float - R0 = d_max * ratio + R0 = d_max * ratio, by default 0.05. s: float - control the decay ratio + control the decay ratio, by default 20.0. Returns ------- - w_mat: NDArray + w_mat: np.ndarray weighting matrix mapping from xy to xyi (xy meshgrid) """ d_mat = _distance_matrix2d(xy, xyi) @@ -148,33 +230,31 @@ def weight_sigmod(xy, xyi, ratio=0.05, s=20.0): r0 = 5.0 * ratio # weights is the sigmod function weight = 1.0 / (1 + np.exp(s * (d_mat - r0))) - # normalized - w_mat = weight / weight.sum(axis=0) - - return w_mat + # weighting matrix normalized + return weight / weight.sum(axis=0) -def weight_idw(xy, xyi, k=6, p=1.0): +def weight_idw( + xy: np.ndarray, xyi: np.ndarray, k: int = 6, p: float = 1.0 +) -> np.ndarray: """ - Description - ----------- (2D only) local weight/interpolate by inverse distance Parameters ---------- - xy: NDArray + xy: np.ndarray (x, y) of values - xyi: NDArray + xyi: np.ndarray (xi, yi) of interpolated locations k: int - number of nearest neighbores + number of nearest neighbores, by default 6. p: float - scaling distance + scaling distance, by default 1.0. Returns ------- - w_mat: NDArray + w_mat: np.ndarray weighting matrix mapping from xy to xy_mesh """ d_mat = _distance_matrix2d(xy, xyi) @@ -184,67 +264,58 @@ def weight_idw(xy, xyi, k=6, p=1.0): for w in weight.T: sort_indices = np.argsort(w) np.put(w, sort_indices[:-k], 0) - # normalized - w_mat = weight / weight.sum(axis=0) - + # weighting matrix normalized # xy times xyi size, use w_mat.T to multiply - return w_mat + return weight / weight.sum(axis=0) -def weight_linear_rbf(xy, xyi, z): +def weight_linear_rbf(xy: np.ndarray, xyi: np.ndarray, z: np.ndarray) -> np.ndarray: """ - Description - ----------- (2D only) local weight/interpolate by linear rbf function (z value required) Parameters ---------- - xy: NDArray + xy: np.ndarray (x, y) of values - xyi: NDArray + xyi: np.ndarray (xi, yi) of interpolated locations + z: np.ndarray + z values Returns ------- - w_mat: NDArray + w_mat: np.ndarray weighting matrix mapping from xy to xy_mesh """ internal_dist = _distance_matrix2d(xy, xy) weights = la.solve(internal_dist, z) - interp_dist = _distance_matrix2d(xy, xyi) - zi = np.dot(interp_dist.T, weights) - - return zi + return np.dot(interp_dist.T, weights) def weight_barycentric_gradient(): """ - Description - ----------- (2D only) local weight/interpolate by barycentric gradient Parameters ---------- - xy: NDArray + xy: np.ndarray (x, y) of values - xyi: NDArray + xyi: np.ndarray (xi, yi) of interpolated locations Returns ------- - w_mat: NDArray + w_mat: np.ndarray weighting matrix mapping from xy to xy_mesh """ raise NotImplementedError() -def sim2pts(pts, sim, sim_values): +def sim2pts(pts: np.ndarray, sim: np.ndarray, sim_values: np.ndarray) -> np.ndarray: """ - Description - ----------- (2D/3D) compatible. Interp values on points using values on simplex, @@ -256,6 +327,16 @@ def sim2pts(pts, sim, sim_values): where r_e is the value on triangles who share the node n, S_e is the area of triangle e. + Parameters + ---------- + pts_values: np.ndarray + Nx1 array, real/complex valued + sim: np.ndarray + Mx3, Mx4 array, elements or simplex + triangles denote connectivity [[i, j, k]] + tetrahedrons denote connectivity [[i, j, m, n]] + sim_value: np.ndarray + Notes ----- This function is similar to pdeprtni of MATLAB pde. @@ -282,10 +363,8 @@ def sim2pts(pts, sim, sim_values): return f / w -def pts2sim(sim, pts_values): +def pts2sim(sim: np.ndarray, pts_values: np.ndarray) -> np.ndarray: """ - Description - ----------- (2D/3D) compatible. Given values on nodes, calculate interpolated values on simplex, @@ -294,16 +373,16 @@ def pts2sim(sim, pts_values): Parameters ---------- - sim: NDArray + sim: np.ndarray Mx3, Mx4 array, elements or simplex triangles denote connectivity [[i, j, k]] tetrahedrons denote connectivity [[i, j, m, n]] - pts_values: NDArray + pts_values: np.ndarray Nx1 array, real/complex valued Returns ------- - el_value: NDArray + el_value: np.ndarray Mx1 array, real/complex valued Notes @@ -311,32 +390,28 @@ def pts2sim(sim, pts_values): This function is similar to pdfinterp of MATLAB pde. """ # averaged over 3 nodes of a triangle - el_value = np.mean(pts_values[sim], axis=1) - return el_value + return np.mean(pts_values[sim], axis=1) -def tri_area(pts, sim): +def tri_area(pts: np.ndarray, sim: np.ndarray) -> np.ndarray: """ calculate the area of each triangle Parameters ---------- - pts: NDArray + pts: np.ndarray Nx2 array, (x,y) locations for points - sim: NDArray + sim: np.ndarray Mx3 array, elements (triangles) connectivity Returns ------- - a: NDArray + a: np.ndarray Areas of triangles """ a = np.zeros(np.shape(sim)[0]) for i, e in enumerate(sim): xy = pts[e] - # s1 = xy[2, :] - xy[1, :] - # s2 = xy[0, :] - xy[2, :] - # s3 = xy[1, :] - xy[0, :] # which can be simplified to # s = xy[[2, 0, 1]] - xy[[1, 2, 0]] s = xy[[2, 0]] - xy[[1, 2]] @@ -347,20 +422,20 @@ def tri_area(pts, sim): return a * 0.5 -def tet_volume(pts, sim): +def tet_volume(pts: np.ndarray, sim: np.ndarray) -> np.ndarray: """ calculate the area of each triangle Parameters ---------- - pts: NDArray + pts: np.ndarray Nx3 array, (x,y, z) locations for points - sim: NDArray + sim: np.ndarray Mx4 array, elements (tetrahedrons) connectivity Returns ------- - v: NDArray + v: np.ndarray Volumes of tetrahedrons """ v = np.zeros(np.shape(sim)[0]) @@ -374,10 +449,8 @@ def tet_volume(pts, sim): return v / 6.0 -def pdetrg(pts, tri): +def pdetrg(pts: np.ndarray, tri: np.ndarray) -> np.ndarray: """ - Description - ----------- (Deprecated) analytical calculate the Area and grad(phi_i) using barycentric coordinates (simplex coordinates) @@ -390,18 +463,18 @@ def pdetrg(pts, tri): Parameters ---------- - pts: NDArray + pts: np.ndarray Nx2 array, (x,y) locations for points - tri: NDArray + tri: np.ndarray Mx3 array, elements (triangles) connectivity Returns ------- - a: NDArray + a: np.ndarray Mx1 array, areas of elements - grad_phi_x: NDArray + grad_phi_x: np.ndarray Mx3 array, x-gradient on elements' local coordinate - grad_phi_y: NDArray + grad_phi_y: np.ndarray Mx3 array, y-gradient on elements' local coordinate """ m = np.size(tri, 0) @@ -429,10 +502,10 @@ def pdetrg(pts, tri): return a, grad_phi_x, grad_phi_y -def pdegrad(pts, tri, node_value): +def pdegrad( + pts: np.ndarray, tri: np.ndarray, node_value: np.ndarray +) -> Tuple[np.ndarray, np.ndarray]: """ - Description - ----------- (Deprecated) given values on nodes, calculate the averaged-grad on elements this function was tested and equivalent to MATLAB 'pdegrad' @@ -440,16 +513,16 @@ def pdegrad(pts, tri, node_value): Parameters ---------- - pts: NDArray + pts: np.ndarray Nx2 array, (x,y) locations for points - tri: NDArray + tri: np.ndarray Mx3 array, elements (triangles) connectivity - node_value: NDArray + node_value: np.ndarray Nx1 array, real/complex valued Returns ------- - el_grad: NDArray + el_grad: np.ndarray el_grad, Mx2 array, real/complex valued """ m = np.size(tri, 0) @@ -460,7 +533,7 @@ def pdegrad(pts, tri, node_value): return grad_el_x, grad_el_y -def demo(): +def demo() -> None: """demo shows how to interpolate on regular/irregular grids""" # 1. create mesh mesh_obj, _ = layer_circle(n_layer=8, n_fan=6) diff --git a/pyeit/eit/jac.py b/pyeit/eit/jac.py index edb157c..a1666bd 100644 --- a/pyeit/eit/jac.py +++ b/pyeit/eit/jac.py @@ -7,6 +7,8 @@ # Distributed under the (new) BSD License. See LICENSE.txt for more info. from __future__ import division, absolute_import, print_function +from typing import Tuple, Union + import numpy as np import scipy.linalg as la @@ -16,48 +18,142 @@ class JAC(EitBase): """A sensitivity-based EIT imaging class""" - def setup(self, p=0.20, lamb=0.001, method="kotre"): + def setup( + self, p: float = 0.20, lamb: float = 0.001, method: str = "kotre" + ) -> None: """ - JAC, Jacobian matrix based reconstruction. + Setup JAC solver + + Jacobian matrix based reconstruction. Parameters ---------- - p, lamb: float - JAC parameters - method: str - regularization methods + p : float, optional + JAC parameters, by default 0.20 + lamb : float, optional + JAC parameters, by default 0.001 + method : str, optional + regularization methods ("kotre", "lm", "dgn" ), by default "kotre" """ # passing imaging parameters self.params = {"p": p, "lamb": lamb, "method": method} # pre-compute H0 for dynamical imaging # H = (J.T*J + R)^(-1) * J.T - self.H = h_matrix(self.J, p, lamb, method) + self.J = self._compute_jac_matrix() + self.H = self._compute_h(self.J, p, lamb, method) + self.is_ready = True - def map(self, dv): - """return Hv""" - return -np.dot(self.H, dv.transpose()) + def _compute_h( + self, jac: np.ndarray, p: float, lamb: float, method: str = "kotre" + ) -> np.ndarray: + """ + Compute self.H matrix for JAC solver + + JAC method of dynamic EIT solver: + H = (J.T*J + lamb*R)^(-1) * J.T + + Parameters + ---------- + jac : np.ndarray + Jacobian + p : float + regularization parameter + lamb : float + regularization parameter + method : str, optional + regularization method, ("kotre", "lm", "dgn" ), by default "kotre" + + Returns + ------- + np.ndarray + H matrix, pseudo-inverse matrix of JAC + """ + j_w_j = np.dot(jac.transpose(), jac) + if method == "kotre": + # see adler-dai-lionheart-2007 + # p=0 : noise distribute on the boundary ('dgn') + # p=0.5 : noise distribute on the middle + # p=1 : noise distribute on the center ('lm') + r_mat = np.diag(np.diag(j_w_j)) ** p + elif method == "lm": + # Marquardt–Levenberg, 'lm' for short + # or can be called NOSER, DLS + r_mat = np.diag(np.diag(j_w_j)) + else: + # Damped Gauss Newton, 'dgn' for short + r_mat = np.eye(jac.shape[1]) - def solve_gs(self, v1, v0): - """solving by weighted frequency""" + # build H + return np.dot(la.inv(j_w_j + lamb * r_mat), jac.transpose()) + + # -------------------------------------------------------------------------- + # Special method for JAC + # -------------------------------------------------------------------------- + + def solve_gs(self, v1: np.ndarray, v0: np.ndarray) -> np.ndarray: + """ + Solving by weighted frequency + + Parameters + ---------- + v1: np.ndarray + current frame + v0: np.ndarray + referenced frame + + Raises + ------ + SolverNotReadyError + raised if solver not ready (see self._check_solver_is_ready()) + + Returns + ------- + np.ndarray + complex-valued np.ndarray, changes of conductivities + """ + self._check_solver_is_ready() a = np.dot(v1, v0) / np.dot(v0, v0) dv = v1 - a * v0 - ds = -np.dot(self.H, dv.transpose()) - # return average epsilon on element - return ds + # return ds average epsilon on element + return -np.dot(self.H, dv.transpose()) - def jt_solve(self, v1, v0, normalize=True): + def jt_solve( + self, v1: np.ndarray, v0: np.ndarray, normalize: bool = True + ) -> np.ndarray: """ a 'naive' back projection using the transpose of Jac. This scheme is the one published by kotre (1989): - [1] Kotre, C. J. (1989). - A sensitivity coefficient method for the reconstruction of - electrical impedance tomograms. - Clinical Physics and Physiological Measurement, - 10(3), 275–281. doi:10.1088/0143-0815/10/3/008 + Parameters + ---------- + v1: np.ndarray + current frame + v0: np.ndarray + referenced frame + normalize : bool, optional + flag to log-normalize the current frame difference dv, by default + True. The input (dv) and output (ds) is log-normalized. + + Raises + ------ + SolverNotReadyError + raised if solver not ready (see self._check_solver_is_ready()) + + Returns + ------- + np.ndarray + complex-valued np.ndarray, changes of conductivities + + Notes + ----- + [1] Kotre, C. J. (1989). + A sensitivity coefficient method for the reconstruction of + electrical impedance tomograms. + Clinical Physics and Physiological Measurement, + 10(3), 275–281. doi:10.1088/0143-0815/10/3/008 - The input (dv) and output (ds) is log-normalized. """ + self._check_solver_is_ready() if normalize: dv = np.log(np.abs(v1) / np.abs(v0)) * np.sign(v0.real) else: @@ -68,45 +164,56 @@ def jt_solve(self, v1, v0, normalize=True): def gn( self, - v, - x0=None, - maxiter=1, - gtol=1e-4, - p=None, - lamb=None, - lamb_decay=1.0, - lamb_min=0, - method="kotre", - verbose=False, + v: np.ndarray, + x0: Union[int, float, np.ndarray] = None, + maxiter: int = 1, + gtol: float = 1e-4, + p: float = None, + lamb: float = None, + lamb_decay: float = 1.0, + lamb_min: float = 0.0, + method: str = "kotre", + verbose: bool = False, **kwargs, - ): + ) -> np.ndarray: """ Gaussian Newton Static Solver You can use a different p, lamb other than the default ones in setup Parameters ---------- - v: NDArray + v : np.ndarray boundary measurement - x0: NDArray, optional - initial guess - maxiter: int, optional - number of maximum iterations - p, lamb: float - JAC parameters (can be overridden) - lamb_decay: float - decay of lamb0, i.e., lamb0 = lamb0 * lamb_delay of each iteration - lamb_min: float - minimal value of lamb - method: str, optional - 'kotre' or 'lm' - verbose: bool, optional - print debug information + x0 : Union[int, float, np.ndarray], optional + initial permittivity guess, by default None + (see Foward._get_perm for more details, in fem.py) + maxiter : int, optional + number of maximum iterations, by default 1 + gtol : float, optional + convergence threshold, by default 1e-4 + p : float, optional + JAC parameters (can be overridden), by default None + lamb : float, optional + JAC parameters (can be overridden), by default None + lamb_decay : float, optional + decay of lamb0, i.e., lamb0 = lamb0 * lamb_delay of each iteration, + by default 1.0 + lamb_min : float, optional + minimal value of lamb, by default 0.0 + method : str, optional + regularization methods ("kotre", "lm", "dgn" ), by default "kotre" + verbose : bool, optional + verbose flag, by default False + + Raises + ------ + SolverNotReadyError + raised if solver not ready (see self._check_solver_is_ready()) Returns ------- - sigma: NDArray - Complex-valued conductivities + np.ndarray + Complex-valued conductivities, sigma Note ---- @@ -116,6 +223,7 @@ def gn( R = diag(J^TJ)**p r0 (residual) = real_measure - forward_v """ + self._check_solver_is_ready() if x0 is None: x0 = self.perm if p is None: @@ -130,16 +238,13 @@ def gn( for i in range(maxiter): - # forward solver - fs = self.fwd.solve_eit( - self.ex_mat, step=self.step, perm=x0, parser=self.parser - ) + # forward solver, + jac, v0 = self._gn_single_jac_v0(x0) # Residual - r0 = v - fs.v - jac = fs.jac + r0 = v - v0 # Damped Gaussian-Newton - h_mat = h_matrix(jac, p, lamb, method) + h_mat = self._compute_h(jac, p, lamb, method) # update d_k = np.dot(h_mat, r0) @@ -156,45 +261,83 @@ def gn( # update regularization parameter # lambda can be given in user defined decreasing lists lamb *= lamb_decay - if lamb < lamb_min: - lamb = lamb_min - + lamb = max(lamb, lamb_min) return x0 - def project(self, ds): + def _gn_single_jac_v0( + self, perm0: Union[int, float, np.ndarray] + ) -> Tuple[np.ndarray, np.ndarray]: + """ + Compute Jacobian and first meas. extimation for gn single step + + Parameters + ---------- + x0 : Union[int, float, np.ndarray] + previous permittivity guess, by default None + (see Foward._get_perm for more details, in fem.py) + + Returns + ------- + Tuple[np.ndarray, np.ndarray] + _description_ + """ + + jac = self._compute_jac_matrix(perm=perm0, allow_jac_norm=False) + v0 = self.fwd.v0 + return jac, v0 + + def project(self, ds: np.ndarray) -> np.ndarray: + """ + Project ds using spatial difference filter (deprecated) + + Parameters + ---------- + ds : np.ndarray + delta sigma (conductivities) + + Returns + ------- + np.ndarray + _description_ + """ """project ds using spatial difference filter (deprecated) Parameters ---------- - ds: NDArray + ds: np.ndarray delta sigma (conductivities) Returns ------- - NDArray + np.ndarray """ d_mat = sar(self.tri) return np.dot(d_mat, ds) -def h_matrix(jac, p, lamb, method="kotre"): +def h_matrix( + jac: np.ndarray, p: float, lamb: float, method: str = "kotre" +) -> np.ndarray: """ + (NOT USED in JAC solver) JAC method of dynamic EIT solver: H = (J.T*J + lamb*R)^(-1) * J.T Parameters ---------- - jac: NDArray + jac : np.ndarray Jacobian - p, lamb: float - regularization parameters - method: str, optional - regularization method + p : float + regularization parameter + lamb : float + regularization parameter + method : str, optional + regularization method, ("kotre", "lm", "dgn" ), by default "kotre" Returns ------- - H: NDArray - pseudo-inverse matrix of JAC + np.ndarray + H matrix, pseudo-inverse matrix of JAC """ j_w_j = np.dot(jac.transpose(), jac) if method == "kotre": @@ -212,23 +355,22 @@ def h_matrix(jac, p, lamb, method="kotre"): r_mat = np.eye(jac.shape[1]) # build H - h_mat = np.dot(la.inv(j_w_j + lamb * r_mat), jac.transpose()) - return h_mat + return np.dot(la.inv(j_w_j + lamb * r_mat), jac.transpose()) -def sar(el2no): +def sar(el2no: np.ndarray) -> np.ndarray: """ - extract spatial difference matrix on the neighbors of each element + Extract spatial difference matrix on the neighbors of each element in 2D fem using triangular mesh. Parameters ---------- - el2no: NDArray + el2no : np.ndarray triangle structures Returns ------- - D: NDArray + np.ndarray SAR matrix """ ne = el2no.shape[0] diff --git a/pyeit/eit/svd.py b/pyeit/eit/svd.py index 1c3388a..b735cb7 100644 --- a/pyeit/eit/svd.py +++ b/pyeit/eit/svd.py @@ -13,21 +13,23 @@ class SVD(JAC): """implementing a sensitivity-based EIT imaging class""" - def setup(self, n=25, rcond=1e-2, method="svd"): + def setup(self, n: int = 25, rcond: float = 1e-2, method: str = "svd") -> None: """ - SVD, singular value decomposition based reconstruction. + Setup of SVD solver, singular value decomposition based reconstruction. Parameters ---------- - n: int - largest n eigenvalues to be kept - rcond: double - r-condition number of pinv - method: string + n : int, optional + largest n eigenvalues to be kept, by default 25 + rcond : float, optional + r-condition number of pinv, by default 1e-2 + method : str, optional + reconstruction method, by default "svd" 'svd': SVD truncation, 'pinv': pseudo inverse """ # correct n_ord + self.J = self._compute_jac_matrix() nm, ne = self.J.shape n_ord = np.min([nm, ne, n]) @@ -35,7 +37,10 @@ def setup(self, n=25, rcond=1e-2, method="svd"): self.params = {"n": n_ord, "rcond": rcond, "method": method} # pre-compute H0 for dynamical imaging - if method == "svd": + if method == "pinv": + self.H = np.linalg.pinv(self.J, rcond=rcond) + + elif method == "svd": JtJ = np.dot(self.J.T, self.J) # using svd @@ -52,5 +57,16 @@ def setup(self, n=25, rcond=1e-2, method="svd"): # pseudo inverse JtJ_inv = np.dot(U, np.dot(np.diag(s**-1), U.T)) self.H = np.dot(JtJ_inv, self.J.T) - elif method == "pinv": - self.H = np.linalg.pinv(self.J, rcond=rcond) + self.is_ready = True + + def gn(self): + """deactivate gn""" + raise NotImplementedError() + + def solve_gs(self): + """deactivate solve_gs""" + raise NotImplementedError() + + def jt_solve(self): + """deactivate jt_solve""" + raise NotImplementedError() diff --git a/pyeit/eit/utils.py b/pyeit/eit/utils.py index cc346d7..7d1d85f 100644 --- a/pyeit/eit/utils.py +++ b/pyeit/eit/utils.py @@ -15,30 +15,38 @@ def eit_scan_lines(n_el: int = 16, dist: int = 1) -> np.ndarray: """ Generate scan matrix, `ex_mat` ( or excitation pattern), see notes - Args: - n_el (int, optional): number of electrodes. Defaults to `16`. - dist (int, optional): distance (number of electrodes) of A to B. Defaults to `1`. + Parameters + ---------- + n_el : int, optional + number of electrodes, by default 16 + dist : int, optional + distance (number of electrodes) of A to B, by default 1 For 'adjacent'- or 'neighbore'-mode (default) use `1` , and - for 'apposition'-mode use `n_el/2`, (see Examples) + for 'apposition'-mode use `n_el/2` (see Examples). - Returns: - np.ndarray: stimulation matrix; shape (n_exc, 2) + Returns + ------- + np.ndarray + stimulation matrix; shape (n_exc, 2) - Notes: + Notes + ----- - in the scan of EIT (or stimulation matrix), we use 4-electrodes mode, where A, B are used as positive and negative stimulation electrodes and M, N are used as voltage measurements. - `1` (A) for positive current injection, `-1` (B) for negative current sink - Examples: + Examples + -------- n_el=16 if mode=='neighbore': ex_mat = eit_scan_lines(n_el=n_el) elif mode=='apposition': ex_mat = eit_scan_lines(dist=n_el/2) - WARNING: + WARNING + ------- `ex_mat` is a local index, where it is ranged from 0...15, within the range of the number of electrodes. In FEM applications, you should convert `ex_mat` to global index using the (global) `el_pos` parameters. diff --git a/pyeit/mesh/__init__.py b/pyeit/mesh/__init__.py index 84ada25..55ca332 100644 --- a/pyeit/mesh/__init__.py +++ b/pyeit/mesh/__init__.py @@ -1,5 +1,13 @@ """ main module for 2D/3D mesh """ from .wrapper import create, set_perm, layer_circle from .shell import multi_shell, multi_circle +from .utils import check_ccw -__all__ = ["create", "set_perm", "layer_circle", "multi_shell", "multi_circle"] +__all__ = [ + "create", + "set_perm", + "layer_circle", + "multi_shell", + "multi_circle", + "check_ccw", +] diff --git a/pyeit/mesh/utils.py b/pyeit/mesh/utils.py index 17b070b..4792c99 100644 --- a/pyeit/mesh/utils.py +++ b/pyeit/mesh/utils.py @@ -127,6 +127,15 @@ def edge_list(tri): return bars[np.array(ix)].view("i") +def check_ccw(no2xy, el2no): + """ + check whether the simplices are CCW ordered, triangles only + """ + xys = no2xy[el2no] + a = [tri_area(xy) > 0 for xy in xys] + return np.all(a) + + def check_order(no2xy, el2no): """ loop over all elements, calculate the Area of Elements (aoe) diff --git a/tests/run_examples.py b/tests/run_examples.py new file mode 100644 index 0000000..1123de0 --- /dev/null +++ b/tests/run_examples.py @@ -0,0 +1,20 @@ + + +import os +import subprocess + +folder = ".\examples" +example_files=[ + "eit_dynamic_bp.py", + "eit_dynamic_jac.py", + "eit_static_jac.py", + "eit_dynamic_greit.py", + "fem_forward2d.py", + "fem_forward3d.py", +] + +for ex in example_files: + path= os.path.join(folder,ex) + cmd= f"python {path}" + print(f"runs >> {cmd}") + subprocess.call(cmd, shell=True) diff --git a/tests/test_eit.py b/tests/test_eit.py index e69de29..9cd1fa6 100644 --- a/tests/test_eit.py +++ b/tests/test_eit.py @@ -0,0 +1,28 @@ +# test for eit +import numpy as np +import pyeit.eit + + +def _mesh_obj(): + """build a simple mesh, which is used in FMMU.CEM""" + node = np.array([[0, 0], [1, 1], [1, 2], [0, 1]]) + element = np.array([[0, 1, 3], [1, 2, 3]]) + perm = np.array([3.0, 1.0]) # assemble should not use perm.dtype + mesh = {"node": node, "element": element, "perm": perm, "ref": 3} + el_pos = np.array([1, 2]) + + return mesh, el_pos + + +def test_bp(): + """test back projection""" + mesh, el_pos = _mesh_obj() + ex_mat = np.array([[0, 1], [1, 0]]) + solver = pyeit.eit.BP(mesh, el_pos, ex_mat, parser="meas_current") + solver.setup() + + assert solver.B.shape[0] > 0 + + +if __name__ == "__main__": + pass diff --git a/tests/test_fem.py b/tests/test_fem.py index fb1cf15..d292a89 100644 --- a/tests/test_fem.py +++ b/tests/test_fem.py @@ -1,4 +1,5 @@ # test for fem.py +import unittest import numpy as np import pyeit.eit.fem @@ -19,89 +20,166 @@ def _assemble(ke, tri, perm, n): return k - def _voltage_meter(ex_line, n_el, dist, parser): """a simple voltage meter""" meas_current = parser == "meas_current" rel_electrode = parser in ["fmmu", "rotate_meas"] - if rel_electrode: - i0 = ex_line[0] - else: - i0 = 0 + i0 = ex_line[0] if rel_electrode else 0 m = np.arange(i0, i0 + n_el) % n_el n = np.arange(i0 + dist, i0 + dist + n_el) % n_el v = np.array([[ni, mi] for ni, mi in zip(n, m)]) keep = [~np.any(np.isin(vi, ex_line), axis=0) for vi in v] - if meas_current: - return v - else: - return v[keep] - - -def test_ke_triangle(): - """test ke calculation using triangle (2D)""" - pts = np.array([[0, 1], [0, 0], [1, 0]]) - tri = np.array([[0, 1, 2]]) - k_truth = np.array([[1, -1, 0], [-1, 2, -1], [0, -1, 1]]) - area = 0.5 - ke = pyeit.eit.fem.calculate_ke(pts, tri) - assert ke.shape == (1, 3, 3) - assert np.allclose(ke[0], k_truth * area) - - -def test_ke_tetrahedron(): - """test ke calculation using tetrahedron (3D)""" - pts = np.array([[0, 0, 0], [0, 0, 1], [0, 1, 0], [1, 0, 0]]) - tri = np.array([[0, 1, 2, 3]]) - k_truth = np.array([[3, -1, -1, -1], [-1, 1, 0, 0], [-1, 0, 1, 0], [-1, 0, 0, 1]]) - volumn = 1 / 6.0 - ke = pyeit.eit.fem.calculate_ke(pts, tri) - assert ke.shape == (1, 4, 4) - assert np.allclose(ke[0], k_truth * volumn) - - -def test_assemble(): - """test assembling coefficients matrix, K""" - np.random.seed(0) - n, ne = 10, 42 - nodes = np.arange(n) - ke = np.random.randn(ne, 3, 3) - tri = np.array([nodes[np.random.permutation(n)[:3]] for _ in range(ne)]) - perm = np.random.randn(ne) - k_truth = _assemble(ke, tri, perm, n) - k = pyeit.eit.fem.assemble(ke, tri, perm, n) - assert np.allclose(k, k_truth) - - -def test_voltage_meter(): - """test voltage meter""" - n_el = 16 - np.random.seed(42) - for parser in ["meas_current", "fmmu", "rotate_meas"]: - ex_lines = [np.random.permutation(n_el)[:2] for _ in range(10)] - for ex_line in ex_lines: - ex_mat = np.array([ex_line]) # two dim - diff_truth = _voltage_meter(ex_line, n_el, 1, parser) - diff = pyeit.eit.fem.voltage_meter(ex_mat, n_el, 1, parser) - assert np.allclose(diff, diff_truth) - - -def test_subtract_row(): - """ - subtract the last dimension - v is [n_exe, n_el, 1] where n_exe is the number of execution (stimulations) - and n_el is the number of voltage sensing electrode, - meas_pattern is [n_exe, n_meas, 2], where n_meas is the effective number - of voltage differences. - - for simplification, we let n_meas=1 - """ - n_exe = 10 - n_el = 16 - v = np.random.randn(n_exe, n_el, 1) - diff_pairs = np.array([np.random.permutation(n_el)[:2] for _ in range(n_exe)]) - vd_truth = np.array([v[i, d[0]] - v[i, d[1]] for i, d in enumerate(diff_pairs)]) - meas_pattern = diff_pairs.reshape(n_exe, 1, 2) - vd = pyeit.eit.fem.subtract_row(v, meas_pattern) - assert vd_truth.size == vd.size - assert np.allclose(vd.ravel(), vd_truth.ravel()) + return v if meas_current else v[keep] + +def _mesh_obj(): + """build a simple mesh, which is used in FMMU.CEM""" + node = np.array([[0.13, 0.15], [0.2, 0.2], [0.1, 0.1], [0.18, 0.12]]) + element = np.array([[0, 2, 3], [0, 3, 1]]) + # assemble uses perm.dtype, perm MUST not be np.int + perm = np.array([3.0, 1.0]) + mesh = {"node": node, "element": element, "perm": perm, "ref": 3} + el_pos = np.array([1, 2]) + + return mesh, el_pos + +class TestFem(unittest.TestCase): + + def test_ke_triangle(self): + """test ke calculation using triangle (2D)""" + pts = np.array([[0, 1], [0, 0], [1, 0]]) + tri = np.array([[0, 1, 2]]) + k_truth = np.array([[1, -1, 0], [-1, 2, -1], [0, -1, 1]]) + area = 0.5 + ke = pyeit.eit.fem.calculate_ke(pts, tri) + + self.assertTrue(ke.shape == (1, 3, 3)) + self.assertTrue(np.allclose(ke[0], k_truth * area)) + + def test_ke_tetrahedron(self): + """test ke calculation using tetrahedron (3D)""" + pts = np.array([[0, 0, 0], [0, 0, 1], [0, 1, 0], [1, 0, 0]]) + tri = np.array([[0, 1, 2, 3]]) + k_truth = np.array([[3, -1, -1, -1], [-1, 1, 0, 0], [-1, 0, 1, 0], [-1, 0, 0, 1]]) + volumn = 1 / 6.0 + ke = pyeit.eit.fem.calculate_ke(pts, tri) + + self.assertTrue(ke.shape == (1, 4, 4)) + self.assertTrue(np.allclose(ke[0], k_truth * volumn)) + + + def test_assemble(self): + """test assembling coefficients matrix, K""" + np.random.seed(0) + n, ne = 10, 42 + nodes = np.arange(n) + ke = np.random.randn(ne, 3, 3) + tri = np.array([nodes[np.random.permutation(n)[:3]] for _ in range(ne)]) + perm = np.random.randn(ne) + k_truth = _assemble(ke, tri, perm, n) + k = pyeit.eit.fem.assemble(ke, tri, perm, n) + + self.assertTrue(np.allclose(k, k_truth)) + + def test_voltage_meter(self): + """test voltage meter""" + n_el = 16 + np.random.seed(42) + for parser in ["meas_current", "fmmu", "rotate_meas"]: + ex_lines = [np.random.permutation(n_el)[:2] for _ in range(10)] + for ex_line in ex_lines: + ex_mat = np.array([ex_line]) # two dim + diff_truth = _voltage_meter(ex_line, n_el, 1, parser) + diff = pyeit.eit.fem.voltage_meter(ex_mat, n_el, 1, parser) + + self.assertTrue(np.allclose(diff, diff_truth)) + + def test_subtract_row(self): + """ + subtract the last dimension + v is [n_exe, n_el, 1] where n_exe is the number of execution (stimulations) + and n_el is the number of voltage sensing electrode, + meas_pattern is [n_exe, n_meas, 2], where n_meas is the effective number + of voltage differences. + + for simplification, we let n_meas=1 + """ + n_exe = 10 + n_el = 16 + v = np.random.randn(n_exe, n_el, 1) + diff_pairs = np.array([np.random.permutation(n_el)[:2] for _ in range(n_exe)]) + vd_truth = np.array([v[i, d[0]] - v[i, d[1]] for i, d in enumerate(diff_pairs)]) + meas_pattern = diff_pairs.reshape(n_exe, 1, 2) + vd = pyeit.eit.fem.subtract_row(v, meas_pattern) + + self.assertTrue(vd_truth.size == vd.size) + self.assertTrue(np.allclose(vd.ravel(), vd_truth.ravel())) + + def test_k(self): + """test K using a simple mesh structure""" + k_truth = np.array( + [ + [3.7391, -0.1521, -1.5, 0.0], + [-0.1521, 0.3695, 0.0, 0.0], + [-1.5, 0.0, 1.5, 0.0], + [0.0, 0.0, 0.0, 1.0], + ] + ) + mesh, _ = _mesh_obj() + n_pts = mesh["node"].shape[0] + ke = pyeit.eit.fem.calculate_ke(mesh["node"], mesh["element"]) + # fix ref to be exactly the one in mesh + k = pyeit.eit.fem.assemble( + ke, mesh["element"], mesh["perm"], n_pts, ref=mesh["ref"] + ) + + self.assertTrue(np.allclose(k, k_truth, rtol=0.01)) + + def test_solve(self): + """test solve using a simple mesh structure""" + mesh, el_pos = _mesh_obj() + f_truth = np.array([-0.27027027, 2.59459459, -0.93693694, 0.0]) + fwd = pyeit.eit.fem.Forward(mesh, el_pos) + # fix ref to be exactly the one in mesh + ex_mat = np.array([[0, 1]]) + fwd.set_ref_el(mesh["ref"]) + f= fwd.solve(ex_mat, perm=mesh["perm"]) + + self.assertTrue(np.allclose(f, f_truth)) + + def test_solve_eit(self): + """test solve using a simple mesh structure""" + mesh, el_pos = _mesh_obj() + f_truth = np.array([-0.27027027, 2.59459459, -0.93693694, 0.0]) + + fwd = pyeit.eit.fem.Forward(mesh, el_pos) + # fix ref to be exactly the one in mesh + fwd.set_ref_el(mesh["ref"]) + ex_mat = np.array([[0, 1], [1, 0]]) + # include voltage differences on driving electrodes + fwd = fwd.solve_eit(ex_mat, parser="meas_current") + vdiff_truth = f_truth[el_pos[1]] - f_truth[el_pos[0]] + v_truth = vdiff_truth * np.array([1, -1, -1, 1]) + + self.assertTrue(np.allclose(fwd.v, v_truth)) + + def test_compute_jac(self): + """test solve using a simple mesh structure""" + #TODO @ liubenyuan please checkt this test + # compute_jac return jac with the "subtract_row part" here you wanted to test only the jac_i get from old solve method + # the jac_i_truth correspond to the jac_truth! + mesh, el_pos = _mesh_obj() + jac_i_truth = np.array([[-0.02556611, 2.67129291], [-0.28431134, -0.08400292]]) + jac_truth = np.array([[-0.25874523, -2.75529584], [ 0.25874523, 2.75529584]]) + + # testing solve + ex_mat = np.array([[0, 1]]) + fwd = pyeit.eit.fem.Forward(mesh, el_pos) + # fix ref to be exactly the one in mesh + fwd.set_ref_el(mesh["ref"]) + jac = fwd.compute_jac(ex_mat, perm=mesh["perm"], parser="meas_current" ) + + self.assertTrue(np.allclose(jac, jac_truth)) + + +if __name__ == "__main__": + unittest.main() \ No newline at end of file diff --git a/tests/test_mesh.py b/tests/test_mesh.py index e69de29..a7bb8a9 100644 --- a/tests/test_mesh.py +++ b/tests/test_mesh.py @@ -0,0 +1,54 @@ +# test for mesh.py +import numpy as np +import pyeit.mesh as pm + + +def test_shape_circle(): + """unit circle mesh""" + + def _fd(pts): + """shape function""" + return pm.shape.circle(pts, pc=[0, 0], r=1.0) + + def _fh(pts): + """distance function""" + r2 = np.sum(pts**2, axis=1) + return 0.2 * (2.0 - r2) + + # build fix points, may be used as the position for electrodes + num = 16 + p_fix = pm.shape.fix_points_circle(ppl=num) + p, t = pm.distmesh.build(_fd, _fh, pfix=p_fix, h0=0.12) + assert p.shape[0] > 0 + assert t.shape[0] > 0 + assert pm.check_ccw(p, t) + + +def test_thorax(): + """Thorax mesh""" + p, t = pm.distmesh.build( + fd=pm.shape.thorax, + fh=pm.shape.area_uniform, + pfix=pm.shape.thorax_pfix, + h0=0.15, + ) + assert p.shape[0] > 0 + assert t.shape[0] > 0 + assert pm.check_ccw(p, t) + + +def test_head_symm(): + """Head mesh""" + p, t = pm.distmesh.build( + fd=pm.shape.head_symm, + fh=pm.shape.area_uniform, + pfix=pm.shape.head_symm_pfix, + h0=0.15, + ) + assert p.shape[0] > 0 + assert t.shape[0] > 0 + assert pm.check_ccw(p, t) + + +if __name__ == "__main__": + pass