From e98616297d553faec56ecd644f5304976d17219f Mon Sep 17 00:00:00 2001 From: byliu Date: Fri, 15 Apr 2022 20:28:40 +0800 Subject: [PATCH] generate head (symmetric) and lung meshes. --- examples/mesh_distmesh2d.py | 137 ++++++++++++++---------------------- pyeit/io/mes.py | 3 +- pyeit/mesh/distmesh.py | 5 +- pyeit/mesh/shape.py | 126 +++++++++++++++++++++++++++------ pyeit/mesh/utils.py | 21 ++++++ pyeit/mesh/wrapper.py | 59 ++++------------ 6 files changed, 198 insertions(+), 153 deletions(-) diff --git a/examples/mesh_distmesh2d.py b/examples/mesh_distmesh2d.py index 3596ee9..a6a2a98 100644 --- a/examples/mesh_distmesh2d.py +++ b/examples/mesh_distmesh2d.py @@ -11,7 +11,36 @@ from pyeit.mesh import shape from pyeit.mesh import distmesh from pyeit.mesh.plot import voronoi_plot -from pyeit.mesh.shape import thorax, area_uniform + + +def mesh_plot(p, t, el_pos=[]): + """helper function to plot mesh""" + fig, ax = plt.subplots() + ax.triplot(p[:, 0], p[:, 1], t) + mesh_center = np.array([np.median(p[:, 0]), np.median(p[:, 1])]) + if len(el_pos) > 0: + ax.plot(p[el_pos, 0], p[el_pos, 1], "ro") # ro : red circle + for i, el in enumerate(el_pos): + xy = np.array([p[el, 0], p[el, 1]]) + text_offset = (xy - mesh_center) * [1, -1] * 0.05 + ax.annotate( + str(i + 1), + xy=xy, + xytext=text_offset, + textcoords="offset points", + color="k", + fontsize=15, + ha="center", + va="center", + ) + xmax, ymax = np.max(p, axis=0) + xmin, ymin = np.min(p, axis=0) + ax.set_xlim([1.2 * xmin, 1.2 * xmax]) + ax.set_ylim([1.2 * ymin, 1.2 * ymax]) + ax.set_aspect("equal") + plt.show() + + return fig def example1(): @@ -29,21 +58,12 @@ def _fh(pts): # build fix points, may be used as the position for electrodes num = 16 p_fix = shape.fix_points_circle(ppl=num) - # firs num nodes are the positions for electrodes el_pos = np.arange(num) # build triangle p, t = distmesh.build(_fd, _fh, pfix=p_fix, h0=0.05) - - # plot - fig, ax = plt.subplots() - ax.triplot(p[:, 0], p[:, 1], t) - ax.plot(p[el_pos, 0], p[el_pos, 1], "ro") - ax.set_aspect("equal") - ax.set_xlim([-1.5, 1.5]) - ax.set_ylim([-1.1, 1.1]) - plt.show() + mesh_plot(p, t, el_pos) def example2(): @@ -54,12 +74,7 @@ def _fd(pts): # build triangle p, t = distmesh.build(_fd, shape.area_uniform, h0=0.1) - - # plot - fig, ax = plt.subplots() - ax.triplot(p[:, 0], p[:, 1], t) - ax.set_aspect("equal") - plt.show() + mesh_plot(p, t) def example3(): @@ -77,14 +92,7 @@ def _fh(pts): # build triangle p, t = distmesh.build(_fd, _fh, h0=0.025) - - # plot - fig, ax = plt.subplots() - ax.triplot(p[:, 0], p[:, 1], t) - ax.set_aspect("equal") - ax.set_xlim([-1.2, 1.2]) - ax.set_ylim([-1, 1]) - plt.show() + mesh_plot(p, t) def example4(): @@ -98,12 +106,7 @@ def _fd(pts): # build triangle p, t = distmesh.build(_fd, shape.area_uniform, bbox=[[-2, -1], [2, 1]], h0=0.15) - - # plot - fig, ax = plt.subplots() - ax.triplot(p[:, 0], p[:, 1], t) - ax.set_aspect("equal") - plt.show() + mesh_plot(p, t) def example5(): @@ -127,58 +130,28 @@ def _fd(pts): # build p, t = distmesh.build(_fd, shape.area_uniform, pfix=p_fix, h0=0.15) - - # plot - fig, ax = plt.subplots() - ax.triplot(p[:, 0], p[:, 1], t) - ax.plot(p_fix[:, 0], p_fix[:, 1], "ro") - ax.set_aspect("equal") - ax.set_xlim([-1.2, 1.2]) - ax.set_ylim([-1.2, 1.2]) - plt.show() + mesh_plot(p, t) -def example6(): +def example_thorax(): """Thorax mesh""" - # build fix points, may be used as the position for electrodes - num = 16 - - el_pos = np.arange(num) - p_fix = [ - (0.1564, 0.6571), - (0.5814, 0.6353), - (0.8298, 0.433), - (0.9698, 0.1431), - (0.9914, -0.1767), - (0.8359, -0.449), - (0.5419, -0.5833), - (0.2243, -0.6456), - (-0.098, -0.6463), - (-0.4181, -0.6074), - (-0.7207, -0.4946), - (-0.933, -0.2647), - (-0.9147, 0.0543), - (-0.8022, 0.3565), - (-0.5791, 0.5864), - (-0.1653, 0.6819), - ] + el_pos = np.arange(16) # build triangles - p, t = distmesh.build(thorax, fh=area_uniform, pfix=p_fix, h0=0.1) - # plot - fig, ax = plt.subplots() - ax.triplot(p[:, 0], p[:, 1], t) - ax.plot(p[el_pos, 0], p[el_pos, 1], "ro") # ro : red circle - ax.set_aspect("equal") - ax.set_xlim([-1.5, 1.5]) - ax.set_ylim([-1.1, 1.1]) - ax.set_title("Thorax mesh") - plt.show() + p, t = distmesh.build( + fd=shape.thorax, fh=shape.area_uniform, pfix=shape.thorax_pfix, h0=0.1 + ) + mesh_plot(p, t, el_pos) def example_head_symm(): """head phantom (symmetric)""" - pass + el_pos = np.arange(16) + # build triangles + p, t = distmesh.build( + fd=shape.head_symm, fh=shape.area_uniform, pfix=shape.head_symm_pfix, h0=0.1 + ) + mesh_plot(p, t, el_pos) def example_voronoi_plot(): @@ -210,20 +183,11 @@ def _fd(pts): # create equal-distributed electrodes p_fix = shape.fix_points_fd(_fd) - + el_pos = np.arange(len(p_fix)) # generate mesh bbox = [[-2, -2], [2, 2]] p, t = distmesh.build(_fd, shape.area_uniform, pfix=p_fix, bbox=bbox, h0=0.1) - - # plot - fig, ax = plt.subplots() - # Draw a unstructured triangular grid as lines and/or markers. - ax.triplot(p[:, 0], p[:, 1], t) - ax.plot(p_fix[:, 0], p_fix[:, 1], "ro") - ax.set_aspect("equal") - ax.set_xlim([-1.5, 1.5]) - ax.set_ylim([-1.2, 1.2]) - plt.show() + mesh_plot(p, t, el_pos) if __name__ == "__main__": @@ -232,6 +196,7 @@ def _fd(pts): # example3() # example4() # example5() - example6() + # example_thorax() + example_head_symm() # example_voronoi_plot() # example_intersect() diff --git a/pyeit/io/mes.py b/pyeit/io/mes.py index ae61238..75567ef 100644 --- a/pyeit/io/mes.py +++ b/pyeit/io/mes.py @@ -162,6 +162,7 @@ def extract_el(fh): def mesh_plot(ax, mesh_obj, el_pos, imstr="", title=None): """plot and annotate mesh""" p, e, perm = mesh_obj["node"], mesh_obj["element"], mesh_obj["perm"] + mesh_center = np.array([np.median(p[:, 0]), np.median(p[:, 1])]) annotate_color = "k" if os.path.exists(imstr): im = plt.imread(imstr) @@ -179,6 +180,7 @@ def mesh_plot(ax, mesh_obj, el_pos, imstr="", title=None): xytext=text_offset, textcoords="offset points", color=annotate_color, + fontsize=15, ha="center", va="center", ) @@ -197,7 +199,6 @@ def mesh_plot(ax, mesh_obj, el_pos, imstr="", title=None): # print the size e, pts, perm = mesh_obj["element"], mesh_obj["node"], mesh_obj["perm"] - mesh_center = np.array([np.median(pts[:, 0]), np.median(pts[:, 1])]) # print('tri size = (%d, %d)' % e.shape) # print('pts size = (%d, %d)' % pts.shape) fig, ax = plt.subplots(1, figsize=(6, 6)) diff --git a/pyeit/mesh/distmesh.py b/pyeit/mesh/distmesh.py index 81f7c89..e03a5a5 100644 --- a/pyeit/mesh/distmesh.py +++ b/pyeit/mesh/distmesh.py @@ -350,6 +350,10 @@ def build( Parameters ---------- + fd, fh : function + signed distance and quality control function + pfix : ndarray + fixed points on the boundary, usually used for electrodes positions maxiter : int, optional maximum iteration numbers, default=1000 @@ -372,7 +376,6 @@ def build( Also, the user should be aware that, equal-edged tetrahedron cannot fill space without gaps. So, in 3D, you can lower dptol, or limit the maximum iteration steps. - """ # parsing arguments # make sure : g_Fscale < 1.5 diff --git a/pyeit/mesh/shape.py b/pyeit/mesh/shape.py index 8e50905..5b8ad4c 100644 --- a/pyeit/mesh/shape.py +++ b/pyeit/mesh/shape.py @@ -340,20 +340,45 @@ def area_uniform(p): return np.ones(p.shape[0]) -def thorax(pts): +def lshape(pts): + """L_shaped mesh (for testing)""" + return dist_diff( + rectangle(pts, p1=[-1, -1], p2=[1, 1]), rectangle(pts, p1=[0, 0], p2=[1, 1]) + ) - """ - p : Nx2 ndarray - returns boolean ndarray specifiying whether the point is inside thorax or not - e.g : - array([[False, False], - [ True, True], - [ True, True]]) +lshape_pfix = np.array( + [ + [1, 0], + [1, -1], + [0, -1], + [-1, -1], + [-1, 0], + [-1, 1], + [0, 1], + [0, 0], + ] +) + +def fd_polygon(poly, pts): + """return signed distance of polygon""" + pts_ = [Point(p) for p in pts] + # calculate signed distance + dist = [poly.exterior.distance(p) for p in pts_] + sign = np.sign([-int(poly.contains(p)) + 0.5 for p in pts_]) + + return sign * dist + + +def thorax(pts): + """ + thorax polygon signed distance function + + Thorax contour points coordinates are taken from + a thorax simulation based on EIDORS """ - # Thorax contour points coordinates are taken from a thorax simulation based on EIDORS - thrx_pts = [ + poly = [ (0.0487, 0.6543), (0.1564, 0.6571), (0.2636, 0.6697), @@ -405,17 +430,76 @@ def thorax(pts): (-0.1653, 0.6819), (-0.0581, 0.6699), ] - thrx = Polygon(thrx_pts) - pts_ = [Point(p) for p in pts] - # calculate signed distance - dist = [thrx.exterior.distance(p) for p in pts_] - sign = np.sign([-int(thrx.contains(p)) + 0.5 for p in pts_]) + poly_obj = Polygon(poly) + return fd_polygon(poly_obj, pts) - return sign * dist - -# L_shaped mesh (for testing) -def L_shaped(pts): - return dist_diff( - rectangle(pts, p1=[-1, -1], p2=[1, 1]), rectangle(pts, p1=[0, 0], p2=[1, 1]) +thorax_pfix = np.array( + [ + (-0.098, -0.6463), + (-0.4181, -0.6074), + (-0.7207, -0.4946), + (-0.933, -0.2647), + (-0.9147, 0.0543), + (-0.8022, 0.3565), + (-0.5791, 0.5864), + (-0.1653, 0.6819), + (0.1564, 0.6571), + (0.5814, 0.6353), + (0.8298, 0.433), + (0.9698, 0.1431), + (0.9914, -0.1767), + (0.8359, -0.449), + (0.5419, -0.5833), + (0.2243, -0.6456), + ] +) + + +head_symm_poly = ( + np.array( + [ + [197, 0], + [188, 43], + [174, 83], + [150, 120], + [118, 148], + [81, 168], + [41, 182], + [0, 186], + [-41, 182], + [-81, 168], + [-118, 148], + [-150, 120], + [-174, 83], + [-188, 43], + [-197, 0], + [-201, -35], + [-194, -70], + [-185, -106], + [-169, -141], + [-148, -177], + [-123, -213], + [-88, -241], + [-45, -259], + [0, -263], + [45, -259], + [88, -241], + [123, -213], + [148, -177], + [169, -141], + [185, -106], + [194, -70], + [201, -35], + ] ) + / 255.0 +) + +head_symm_pfix = np.array(head_symm_poly[::-2]) + + +def head_symm(pts): + """symmetric head polygon""" + poly_obj = Polygon(head_symm_poly) + return fd_polygon(poly_obj, pts) diff --git a/pyeit/mesh/utils.py b/pyeit/mesh/utils.py index 1295622..17b070b 100644 --- a/pyeit/mesh/utils.py +++ b/pyeit/mesh/utils.py @@ -194,6 +194,27 @@ def tet_volume(xyz): return v_tot +def to_polar(xy, shift=True, sort=True): + vec = xy + if shift: + pc = np.median(xy, axis=0) + print(pc) + vec = vec - pc + dist = np.sqrt(np.sum(vec**2, axis=1)) + deg = np.rad2deg(np.arctan2(vec[:, 1], vec[:, 0])) + deg = deg % 360 + if sort: + ind = np.argsort(deg) + dist, deg = dist[ind], deg[ind] + return dist, deg + + +def to_xy(dist, deg): + x = dist * np.cos(np.deg2rad(deg)) + y = dist * np.sin(np.deg2rad(deg)) + return x, y + + if __name__ == "__main__": # test 'edge_project' def fd_test(p): diff --git a/pyeit/mesh/wrapper.py b/pyeit/mesh/wrapper.py index 402fa8e..f123504 100644 --- a/pyeit/mesh/wrapper.py +++ b/pyeit/mesh/wrapper.py @@ -10,11 +10,10 @@ from .distmesh import build from .mesh_circle import MeshCircle from .utils import check_order -from .shape import circle, area_uniform, ball, thorax, L_shaped -from .shape import fix_points_fd, fix_points_ball +from . import shape -def create(n_el=16, fd=None, fh=area_uniform, h0=0.1, p_fix=None, bbox=None): +def create(n_el=16, fd=None, fh=shape.area_uniform, h0=0.1, p_fix=None, bbox=None): """ Generating 2D/3D meshes using distmesh (pyEIT built-in) @@ -42,22 +41,21 @@ def create(n_el=16, fd=None, fh=area_uniform, h0=0.1, p_fix=None, bbox=None): # test conditions if fd or/and bbox are none if bbox is None: - if fd != ball: + if fd != shape.ball: bbox = np.array([[-1, -1], [1, 1]]) else: bbox = [[-1.2, -1.2, -1.2], [1.2, 1.2, 1.2]] - bbox = np.array( - bbox - ) # list is converted to Numpy array so we can use it then (calling shape method..) + # list is converted to Numpy array so we can use it then (calling shape method..) + bbox = np.array(bbox) n_dim = bbox.shape[1] # bring dimension # infer dim if fd is None: if n_dim == 2: - fd = circle + fd = shape.circle elif n_dim == 3: - fd = ball + fd = shape.ball if n_dim not in [2, 3]: raise TypeError("distmesh only supports 2D or 3D") @@ -66,44 +64,17 @@ def create(n_el=16, fd=None, fh=area_uniform, h0=0.1, p_fix=None, bbox=None): if p_fix is None: if n_dim == 2: - if fd == thorax: - # thorax shape is generated so far without fixed points (to be updated later) - p_fix = [ - (-0.098, -0.6463), - (-0.4181, -0.6074), - (-0.7207, -0.4946), - (-0.933, -0.2647), - (-0.9147, 0.0543), - (-0.8022, 0.3565), - (-0.5791, 0.5864), - (-0.1653, 0.6819), - (0.1564, 0.6571), - (0.5814, 0.6353), - (0.8298, 0.433), - (0.9698, 0.1431), - (0.9914, -0.1767), - (0.8359, -0.449), - (0.5419, -0.5833), - (0.2243, -0.6456), - ] - p_fix = np.array(p_fix) - elif fd == L_shaped: - p_fix = [ - [1, 0], - [1, -1], - [0, -1], - [-1, -1], - [-1, 0], - [-1, 1], - [0, 1], - [0, 0], - ] # values brought from distmesh2D L shaped mesh example - p_fix = np.array(p_fix) + if fd == shape.thorax: + p_fix = shape.thorax_pfix + elif fd == shape.head_symm: + p_fix = shape.head_symm_pfix + elif fd == shape.lshape: + p_fix = shape.lshape_pfix h0 = 0.15 else: - p_fix = fix_points_fd(fd, n_el=n_el) + p_fix = shape.fix_points_fd(fd, n_el=n_el) elif n_dim == 3: - p_fix = fix_points_ball(n_el=n_el) + p_fix = shape.fix_points_ball(n_el=n_el) # 1. build mesh p, t = build(fd, fh, pfix=p_fix, bbox=bbox, h0=h0)