Examples¶
The examples/ folder contains one Python script per topic. Run them from the repository root; each one opens a polyscope viewer if it is installed (pip install polyscope), or pass --no-gui:
Read a tet mesh and extract its boundary surface¶
examples/01_read_and_extract_surface.py
"""Read a tet mesh and extract its boundary surface.
python examples/01_read_and_extract_surface.py
"""
import libmat
from _common import data, parse_args, polyscope_or_none
args = parse_args(__doc__)
V, T = libmat.read_tet_mesh(data("lblock.tet"))
SV, SF, sf2tet = libmat.extract_surface(V, T)
print(f"tet mesh: {len(V)} vertices, {len(T)} tets")
print(f"surface: {len(SV)} vertices, {len(SF)} triangles")
ps = polyscope_or_none(args)
if ps:
tm = ps.register_volume_mesh("tet mesh", V, tets=T)
tm.set_cull_whole_elements(True)
ps.register_surface_mesh("surface", SV, SF, enabled=False)
ps.show()
Detect sharp edges, concave edges and corners of a CAD-like model¶
examples/02_detect_sharp_features.py
"""Detect sharp edges, concave edges and corners of a CAD-like model.
python examples/02_detect_sharp_features.py
"""
import numpy as np
import libmat
from _common import data, parse_args, polyscope_or_none
args = parse_args(__doc__)
V, F = libmat.read_surface_mesh(data("mesh_69058.obj"))
E, E_type, E_line, corners, corner_type = libmat.detect_sharp_features(
V, F, thres_concave=0.2, thres_convex=40.0)
print(f"{np.sum(E_type == 1)} sharp edges, {np.sum(E_type == 2)} concave edges")
print(f"{len(np.unique(E_line))} feature lines")
print(f"{np.sum(corner_type == 1)} sharp corners, "
f"{np.sum(corner_type == 2)} concave corners")
ps = polyscope_or_none(args)
if ps:
ps.register_surface_mesh("surface", V, F, color=(0.9, 0.9, 0.9))
sharp = ps.register_curve_network("sharp edges", V, E[E_type == 1],
color=(0.1, 0.3, 0.9), radius=0.004)
concave = ps.register_curve_network("concave edges", V, E[E_type == 2],
color=(0.9, 0.2, 0.1), radius=0.006)
lines = ps.register_curve_network("feature lines", V, E, enabled=False)
lines.add_scalar_quantity("line id", E_line, defined_on="edges",
enabled=True, cmap="rainbow")
ps.register_point_cloud("corners", V[corners], radius=0.01)
ps.show()
Sample medial spheres with the ball-shrinking algorithm¶
examples/03_shrink_spheres.py
"""Sample medial spheres with the ball-shrinking algorithm.
python examples/03_shrink_spheres.py [--num-spheres 2000]
"""
import numpy as np
import libmat
from _common import data, output, parse_args, polyscope_or_none
args = parse_args(__doc__, num_spheres=dict(type=int, default=2000),
mesh=dict(default=data("mesh_69058.obj")))
V, F = libmat.read_surface_mesh(args.mesh)
C, R, types = libmat.shrink_spheres(V, F, num_spheres=args.num_spheres,
preserve_features=True)
print(f"{len(C)} medial spheres, radius in [{R.min():.4f}, {R.max():.4f}]")
for t, n in zip(*np.unique(types, return_counts=True)):
print(f" {libmat.SPHERE_TYPES.get(t, t)}: {n}")
libmat.write_ma(output("shrink_spheres.ma"), C, R)
print("saved", output("shrink_spheres.ma"))
ps = polyscope_or_none(args)
if ps:
ps.register_surface_mesh("surface", V, F, transparency=0.3)
pc = ps.register_point_cloud("medial spheres", C)
pc.add_scalar_quantity("radius", R, enabled=True)
pc.set_point_radius_quantity("radius", autoscale=False)
ps.register_point_cloud("sphere centers", C, radius=0.002, enabled=False)
ps.show()
Neighborhood of medial spheres: regular triangulation (power diagram dual)¶
examples/04_regular_triangulation.py
"""Neighborhood of medial spheres: regular triangulation (power diagram dual).
python examples/04_regular_triangulation.py
"""
import numpy as np
import libmat
from _common import data, parse_args, polyscope_or_none
args = parse_args(__doc__)
V, F = libmat.read_surface_mesh(data("mesh_69058.obj"))
C, R, _ = libmat.shrink_spheres(V, F, num_spheres=1000)
T, E, valid = libmat.regular_triangulation(C, R)
print(f"{len(C)} spheres, {(~valid).sum()} hidden")
print(f"regular triangulation: {len(E)} edges, {len(T)} tets")
degree = np.bincount(E.ravel(), minlength=len(C))
print(f"neighbors per sphere: mean {degree.mean():.1f}, max {degree.max()}")
ps = polyscope_or_none(args)
if ps:
ps.register_surface_mesh("surface", V, F, transparency=0.2)
ps.register_point_cloud("centers", C, radius=0.003)
cn = ps.register_curve_network("power diagram neighbors", C, E, radius=0.0007)
cn.add_scalar_quantity("degree", degree, defined_on="nodes")
ps.show()
Visualize a medial mesh (.ma) and the shape it represents¶
examples/05_medial_mesh_visualization.py
"""Visualize a medial mesh (.ma) and the shape it represents.
The shape is the union of medial spheres, cones and slabs. Its boundary is
made of the triangles tangent to the three spheres of every medial face
(plus pieces of cones, not shown here).
python examples/05_medial_mesh_visualization.py [--ma data/bug.ma]
"""
import numpy as np
import libmat
from _common import data, parse_args, polyscope_or_none
args = parse_args(__doc__, ma=dict(default=data("bug.ma")),
surface=dict(default=data("bug.obj")))
V, R, E, F = libmat.read_ma(args.ma)
print(f"medial mesh: {len(V)} spheres, {len(E)} edges, {len(F)} faces")
print(f"Euler characteristic V - E + F = {len(V) - len(E) + len(F)}")
T0, T1, valid = libmat.slab_tangent_triangles(V, R, F)
print(f"{valid.sum()}/{len(F)} slabs have tangent triangles")
# outward oriented triangle soup
tris = np.concatenate([T0[valid], T1[valid][:, [0, 2, 1]]])
env_V = tris.reshape(-1, 3)
env_F = np.arange(len(env_V)).reshape(-1, 3)
ps = polyscope_or_none(args)
if ps:
mm = ps.register_surface_mesh("medial mesh", V, F)
mm.set_back_face_policy("identical")
mm.add_scalar_quantity("radius", R, enabled=True)
ps.register_curve_network("medial edges", V, E, radius=0.0005, enabled=False)
ps.register_surface_mesh("slab envelope", env_V, env_F, enabled=False,
transparency=0.5)
SV, SF = libmat.read_surface_mesh(args.surface)
ps.register_surface_mesh("input surface", SV, SF, enabled=False)
ps.show()
Measure how well a medial mesh approximates a surface¶
examples/06_distance_to_medial_mesh.py
"""Measure how well a medial mesh approximates a surface.
Computes the signed distance from every surface vertex to the volume of the
medial mesh (union of spheres, cones and slabs), on the CPU and, if libmat
was built with CUDA, on the GPU.
python examples/06_distance_to_medial_mesh.py
"""
import time
import numpy as np
import libmat
from _common import data, parse_args, polyscope_or_none
args = parse_args(__doc__)
SV, SF = libmat.read_surface_mesh(data("bug.obj"))
V, R, E, F = libmat.read_ma(data("bug.ma"))
diag = np.linalg.norm(SV.max(0) - SV.min(0))
devices = ["cpu"] + (["cuda"] if libmat.has_cuda() else [])
for device in devices:
t = time.time()
D, prim_type, prim_id = libmat.distance_to_medial_mesh(SV, V, R, E, F,
device=device)
print(f"[{device}] {len(SV)} points x {len(V) + len(E) + len(F)} "
f"primitives in {time.time() - t:.2f}s")
print(f"[{device}] one-sided Hausdorff distance surface -> MAT: "
f"{np.abs(D).max():.4f} ({100 * np.abs(D).max() / diag:.3f}% of bbox)")
print(f"[{device}] closest primitives: "
f"{np.bincount(prim_type, minlength=3)} (spheres, cones, slabs)")
ps = polyscope_or_none(args)
if ps:
sm = ps.register_surface_mesh("surface", SV, SF)
sm.add_scalar_quantity("distance to MAT", np.abs(D), enabled=True,
cmap="reds")
sm.add_scalar_quantity("closest primitive type", prim_type)
mm = ps.register_surface_mesh("medial mesh", V, F, enabled=False)
mm.set_back_face_policy("identical")
ps.show()
Topology-preserving thinning of a medial mesh¶
examples/07_thin_medial_mesh.py
"""Topology-preserving thinning of a medial mesh.
A small medial mesh with a medial tet (4 spheres whose 4 triangles are all
medial faces) attached to a sheet. Thinning removes the tet and the faces of
low importance, without changing the topology.
python examples/07_thin_medial_mesh.py
"""
import numpy as np
import libmat
from _common import parse_args, polyscope_or_none
args = parse_args(__doc__)
# a 4x4 grid sheet in the plane z = 0 ...
n = 4
xs, ys = np.meshgrid(np.arange(n), np.arange(n))
C = np.stack([xs.ravel(), ys.ravel(), np.zeros(n * n)], axis=1).astype(float)
F = []
for i in range(n - 1):
for j in range(n - 1):
a, b, c, d = i * n + j, i * n + j + 1, (i + 1) * n + j, (i + 1) * n + j + 1
F += [[a, b, d], [a, d, c]]
# ... plus a medial tet on top of the triangle (5, 6, 10)
C = np.vstack([C, [1.3, 1.6, 1.0]])
apex = len(C) - 1
F += [[5, 6, apex], [6, 10, apex], [5, 10, apex]]
F = np.array(F)
R = np.full(len(C), 0.2)
# importance: 1 for the sheet, small for the tet faces above it
importance = np.ones(len(F))
importance[-3:] = [0.1, 0.2, 0.3]
E_out, F_out = libmat.thin_medial_mesh(C, R, None, F, importance, threshold=0.5)
def euler(E, F):
"""V - E + F - T, counting medial tets as (filled) 3-cells."""
faces = {tuple(sorted(f)) for f in F}
verts = np.unique(np.concatenate([E.ravel(), F.ravel()]))
tets = {tuple(sorted(f + (d,))) for f in faces for d in verts
if d not in f and all(tuple(sorted(g)) in faces
for g in [(f[0], f[1], d), (f[0], f[2], d), (f[1], f[2], d)])}
return len(verts) - len(E) + len(F) - len(tets), len(tets)
E_in = np.unique(np.sort(np.concatenate([F[:, [0, 1]], F[:, [1, 2]], F[:, [0, 2]]]), axis=1), axis=0)
for name, (E_, F_) in [("before", (E_in, F)), ("after", (E_out, F_out))]:
chi, num_tets = euler(E_, F_)
print(f"{name:6s}: {len(F_)} faces, {len(E_)} edges, {num_tets} medial tets, "
f"Euler characteristic {chi}")
ps = polyscope_or_none(args)
if ps:
ps.register_surface_mesh("before", C, F, transparency=0.4).set_back_face_policy("identical")
after = ps.register_surface_mesh("after", C, F_out)
after.set_back_face_policy("identical")
after.translate((n + 1.0, 0, 0))
ps.show()
Full medial mesh pipeline on the GPU (requires libmat built with CUDA)¶
examples/08_medial_mesh_from_spheres_cuda.py
"""Full medial mesh pipeline on the GPU (requires libmat built with CUDA).
tet mesh -> boundary surface -> ball-shrinking spheres -> restricted power
diagram (CUDA) -> medial mesh -> thinning -> .ma file.
python examples/08_medial_mesh_from_spheres_cuda.py
"""
import sys
import numpy as np
import libmat
from _common import data, output, parse_args, polyscope_or_none
args = parse_args(__doc__, num_spheres=dict(type=int, default=1500),
threshold=dict(type=float, default=0.3))
if not libmat.has_cuda():
print("libmat was built without CUDA; rebuild with LIBMAT_WITH_CUDA=ON:")
print(" LIBMAT_WITH_CUDA=ON pip install .")
sys.exit(0)
V, T = libmat.read_tet_mesh(data("lblock.tet"))
SV, SF, _ = libmat.extract_surface(V, T)
C, R, types = libmat.shrink_spheres(SV, SF, num_spheres=args.num_spheres)
print(f"{len(C)} medial spheres")
E, F, importance = libmat.medial_mesh_from_spheres(V, T, C, R)
print(f"medial mesh: {len(E)} edges, {len(F)} faces")
E, F = libmat.thin_medial_mesh(C, R, E, F, importance, args.threshold)
print(f"after thinning: {len(E)} edges, {len(F)} faces")
libmat.write_ma(output("lblock_mat.ma"), C, R, E, F)
print("saved", output("lblock_mat.ma"))
ps = polyscope_or_none(args)
if ps:
ps.register_surface_mesh("surface", SV, SF, transparency=0.3)
mm = ps.register_surface_mesh("medial mesh", C, F)
mm.set_back_face_policy("identical")
mm.add_scalar_quantity("radius", R, enabled=True)
ps.show()