Skip to content

Geometry Processing and Visualization

Installation

All libraries (igl, meshplot, wildmeshing, polyfempy) can be installed from Conda forge with conda install -c conda-forge <lib>.

Import

import igl
import meshplot as mp
import numpy as np

Hello World

v, f = igl.read_triangle_mesh("data/bunny.obj")
mp.plot(v, f, v[:, 0])

Data Representation

Meshes are represented by a numpy array of vertex coordinates (nx3) and a numpy array of face indices (mx3) and can be loaded from 3D file formats.

v: np.array # Array of vertex coordinates (nx3)
f: np.array # Array of face indices (mx3)
v, f = igl.read_triangle_mesh("data/bunny.obj")
print(v)
print(f)

[[-0.0260146   0.112578    0.0363871 ]
 [-0.0321783   0.174119   -0.00263321]
 [-0.080718    0.152855    0.0302446 ]
 ...
 [-0.023099    0.156978   -0.00584018]
 [-0.0713101   0.15068    -0.0435721 ]
 [-0.0396435   0.152397   -0.00721968]]
[[2784 2497 2027]
 [1077  225 1060]
 [ 425  450  381]
 ...
 [3086 3203 3162]
 [3086 3162 3151]
 [3086 3151 3085]]

Visualizing Surfaces and Pointclouds

mp.plot(v, f, c=v[:, 0])

mp.plot(v, c=v[:, 0])

Scalar and Vector Field Visualization

# Calculating properties of the mesh
nf = igl.per_face_normals(v, f, np.array([0.0, 0.0, 0.0]))
nfv = np.linalg.norm(nf * 0.5 + 0.5, axis=1)
nv = igl.per_vertex_normals(v, f)
area = igl.doublearea(v, f)
avg = igl.avg_edge_length(v, f)

v1, v2, k1, k2 = igl.principal_curvature(v, f)
mean_curv = 0.5 * (k1 + k2)

# Face normals
d = mp.subplot(v, f, c=nfv, s=[2, 2, 0])

# Vertex normals as lines
mp.subplot(v, c=nv, s=[2, 2, 1], data=d)
d.rows[0][1].add_lines(v, v + nv * avg)

# Mean curvature + directions
mp.subplot(v, f, c=mean_curv, s=[2, 2, 2], data=d)
d.rows[1][0].add_lines(v + v1 * avg/3, v - v1 * avg/3, shading={"line_color": "red"})
d.rows[1][0].add_lines(v + v2 * avg/3, v - v2 * avg/3, shading={"line_color": "green"})

# Triangle area
mp.subplot(v, f, c=-area, s=[2, 2, 3], shading={"metalness": 0.0, "roughness": 1.0}, data=d)

Texture Mapping

vc, fc = igl.read_triangle_mesh("data/camelhead.off")

# Find the open boundary
bnd = igl.boundary_loop(fc)

# Map the boundary to a circle, preserving edge proportions
bnd_uv = igl.map_vertices_to_circle(vc, bnd)

# Harmonic parametrization for the internal vertices
uv = igl.harmonic_weights(vc, fc, bnd, bnd_uv, 1)

# Plotting the results
p = mp.subplot(vc, fc, uv=uv, s=[1, 2, 0])
mp.subplot(uv, fc, uv=uv, shading={"wireframe": True, "wire_color": "red", "wire_width": 1.0}, s=[1, 2, 1], data=p)

# Adding the boundary
p.rows[0][0].add_points(vc[bnd], shading={"point_size": 0.1});
lines = np.vstack([bnd, np.roll(bnd, -1)]).T
p.rows[0][0].add_edges(vc, lines, shading={"line_width": 1.0, "line_color": "red"});
p

Interactive Geometry Modification

from scipy.sparse.linalg import spsolve

v, f = igl.read_triangle_mesh("data/cow.off")
l = igl.cotmatrix(v, f)

n = igl.per_vertex_normals(v, f) * 0.5 + 0.5
c = np.linalg.norm(n, axis=1)
p = mp.plot(v, f, c, return_plot=True)

# Precalculate intermediate states
vs = [v]
cs = [c]
for i in range(10):
    m = igl.massmatrix(v, f, igl.MASSMATRIX_TYPE_BARYCENTRIC)
    s = (m - 0.001 * l)
    b = m.dot(v)
    v = spsolve(s, m.dot(v))
    n = igl.per_vertex_normals(v, f)*0.5+0.5
    c = np.linalg.norm(n, axis=1)
    vs.append(v)
    cs.append(c)

# Add interactive visulization
@mp.interact(level=(0, 9))
def mcf(level=0):
    p.update_object(vertices=vs[level], colors=cs[level])
p

Saving Results and Offline Visualization

# Save the previous result
p.save("test1.html")

# Load a new mesh
v, f = igl.read_triangle_mesh("data/bumpy.off")

# Switch to offline plotting
mp.offline()
mp.plot(v, f, c=np.random.rand(*f.shape))

# Switch to jupyter plotting
mp.jupyter()
p = mp.plot(v, f, c=np.random.rand(*f.shape), shading={"width":900, "height": 1000}, return_plot=True)
p.add_mesh(v + 5, f, c=v[:,1]);
p.add_points(v - 5, c=v[:,2], shading={"point_size": 1.0})
p.save("test2.html")
p

Plot saved to file test1.html.
Plot saved to file a96186d9-a92c-45fb-8266-c10648ee498d.html.



Renderer(camera=PerspectiveCamera(aspect=0.9, children=(DirectionalLight(color='white', intensity=0.6, positio…


Plot saved to file test2.html.

Advanced Examples

import scipy as sp

v, f = igl.read_triangle_mesh("data/beetle.off")
l = -igl.cotmatrix(v, f)
m = igl.massmatrix(v, f, igl.MASSMATRIX_TYPE_VORONOI)

d, u = sp.sparse.linalg.eigsh(l, 10, m, sigma=0, which="LM")

u = (u - np.min(u)) / (np.max(u) - np.min(u))
bbd = 0.5 * np.linalg.norm(np.max(v, axis=0) - np.min(v, axis=0))

p = mp.plot(v, f, bbd * u[:, 0], shading={"wireframe":False, "flat": False}, return_plot=True)

@mp.interact(ev=[("EV %i"%i, i) for i in range(10)])
def sf(ev):
    p.update_object(colors=u[:, ev])
p

v, f = igl.read_triangle_mesh("data/bunny.obj")

# Select a vertex from which the distances should be calculated
vs = np.array([0])
# All vertices are the targets
vt = np.arange(v.shape[0])

d = igl.exact_geodesic(v, f, vs, vt)

# The function should be 1 on each integer coordinate
c = np.abs(np.sin((d / 0.04 * np.pi)))
mp.plot(v, f, c, shading={"wireframe": False}, return_plot=True)

Renderer(camera=PerspectiveCamera(children=(DirectionalLight(color='white', intensity=0.6, position=(-0.016860…

v, f = igl.read_triangle_mesh("data/decimated-max.obj")
v[:,[0, 2]] = v[:,[2, 0]] # Swap X and Z axes
u = v.copy()

s = igl.read_dmat("data/decimated-max-selection.dmat")
b = np.array([[t[0] for t in [(i, s[i]) for i in range(0, v.shape[0])] if t[1] >= 0]]).T

## Boundary conditions directly on deformed positions
u_bc = np.zeros((b.shape[0], v.shape[1]))
v_bc = np.zeros((b.shape[0], v.shape[1]))

for bi in range(b.shape[0]):
    v_bc[bi] = v[b[bi]]

    if s[b[bi]] == 0: # Don't move handle 0
        u_bc[bi] = v[b[bi]]
    elif s[b[bi]] == 1: # Move handle 1 down
        u_bc[bi] = v[b[bi]] + np.array([[0, -50, 0]])
    else: # Move other handles forward
        u_bc[bi] = v[b[bi]] + np.array([[-25, 0, 0]])

p = mp.plot(v, f, s, shading={"colormap": "tab10"}, return_plot=True)

@mp.interact(deformation_field=True, step=(0.0, 2.0))
def update(deformation_field, step=0.0):
    # Determine boundary conditions
    u_bc_anim = v_bc + step * (u_bc - v_bc)

    if deformation_field:
        d_bc = u_bc_anim - v_bc
        d = igl.harmonic_weights(v, f, b, d_bc, 2)
        u = v + d
    else:
        u = igl.harmonic_weights(v, f, b, u_bc_anim, 2)
    p.update_object(vertices=u)
p