Skip to content

Commit

Permalink
generate head (symmetric) and lung meshes.
Browse files Browse the repository at this point in the history
  • Loading branch information
@liubenyuan
liubenyuan committed Apr 15, 2022
1 parent 14b17db commit e986162
Show file tree
Hide file tree
Showing 6 changed files with 198 additions and 153 deletions.
137 changes: 51 additions & 86 deletions examples/mesh_distmesh2d.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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():
Expand All @@ -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():
Expand All @@ -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():
Expand All @@ -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():
Expand All @@ -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():
Expand All @@ -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():
Expand Down Expand Up @@ -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__":
Expand All @@ -232,6 +196,7 @@ def _fd(pts):
# example3()
# example4()
# example5()
example6()
# example_thorax()
example_head_symm()
# example_voronoi_plot()
# example_intersect()
3 changes: 2 additions & 1 deletion pyeit/io/mes.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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",
)
Expand All @@ -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))
Expand Down
5 changes: 4 additions & 1 deletion pyeit/mesh/distmesh.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down
126 changes: 105 additions & 21 deletions pyeit/mesh/shape.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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),
Expand Down Expand Up @@ -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)
Loading

0 comments on commit e986162

Please sign in to comment.