Remeshing 3D triangulated surface for better mesh quaility

Question

I am using marching cubes to extract a 2D surface from a volume. In this case a Gyroid.

import numpy as np
from numpy import sin, cos, pi
from skimage import measure
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

def gyroid(x, y, z, t):
    return cos(x)*sin(y) + cos(y)*sin(z) + cos(z)*sin(x) - t

lattice_param = 1.0
strut_param = 0.0
resolution = 31j

x, y, z = pi*np.mgrid[-1:1:resolution, -1:1:resolution, -1:1:resolution] * lattice_param
vol = gyroid(x, y, z, strut_param)

verts, faces = measure.marching_cubes(vol, 0, spacing=(0.1, 0.1, 0.1)) # , normals, values

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_trisurf(verts[:, 0], verts[:, 1], faces, verts[:, 2], cmap='ocean', lw=1)

This all works fine but the mesh quailty is appalling in a lot of places. I cant run any FEA on the meshes as many of the elements/faces have near zero area or are highly distorted.

Is there a way of either remeshing given the vertices & ensuring particular elementface/facet metrics (such as aspect ratio) or forcing marching cubes to do such a thing?

I am not bothered about moving vertices as long as the mesh is a fair approximation.

Answer 1

One option could be to 're-mesh' the marching cubes output using a surface meshing package. Essentially this means the marching cubes triangulation would serve as an initial surface definition to be re-triangulated.

There are many techniques one could employ to do this. A few options that may be useful (all C++ / C implementations):

• JIGSAW: a restricted, frontal-delaunay algorithm^^ that will generally build very high-quality surface Delaunay triangulations. For the type of object shown I'd expect it to work well. In the demos included (provided in MATLAB) several examples address the re-meshing of marching cubes output.
• CGAL: a restricted, delaunay-refinement approach that can also build surface Delaunay triangulations, but uses a slightly different algorithm to JIGSAW and also includes CVT-type mesh optimisation schemes.
• MMG: a collection of re-meshing/optimisation strategies (as I understand it) that can be used to transform (and hence improve) an initial mesh via the iterative application of local modifications.

^^ I'm the author of JIGSAW, so, shameless promotion here basically.

Answer 2

I suggest using 3d adaptive remeshing with smoothing. The Geogram programming library (http://alice.loria.fr/software/geogram/doc/html/index.html) provides a nice implementation. Please see: https://twitter.com/brunolevy01/status/1132343120690122752?lang=en

Answer 3

Use can use pygalmesh for surface remeshing, too. (It's a easy-to-use interface to CGAL.)

import pygalmesh

# create verts, faces

meshio.write_points_cells("in.vtu", verts, [("triangle", faces)])

mesh = pygalmesh.remesh_surface(
    "in.vtu",
    edge_size=0.025,
    facet_angle=25,
    facet_size=0.1,
    facet_distance=0.001,
    verbose=False,
)
# mesh.points, mesh.cells

Answer 4

Weird triangles come from weird data, not from the method used for triangulation.

I can say that a Delaunay triangulation achieves the best triangle-area/triangle-perimeter ratio (best theorethical ratio is with equilateral triangles). But you can't use if for your mesh because Delaunay triangulation outputs a convex mesh.

You have a tight task ahead. Some thoughts:

• Remove cuasi-duplicated points before meshing. This will avoid many near-zero-area triangles.
• Detect thin triangles. Subdivide them into more 'regular' triangles.
• You may create a dense regular grid. Each z can be found traversing the mesh until you find a triangle that contains x,ycoordinates of the point in the grid and then interpolating the z. If you have some extra info, like neigbours or hierachy, the search can be done faster than inspecting every triangle. Triangulating the grid is quite simple.