DrBwts
DrBwts

Reputation: 3657

Remeshing 3D triangulated surface for better mesh quaility

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)

enter image description here

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.

enter image description here

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.

Upvotes: 5

Views: 4785

Answers (5)

Marc
Marc

Reputation: 83

I don't have enough reputation to post a comment so I am submitting this as an answer...

Adding to Nico's answer, it seems that the keyword args to pygalmesh.remesh_surface() have changed. The example provided in the link is:

mesh = pygalmesh.remesh_surface(
    "lion-head.off",
    max_edge_size_at_feature_edges=0.025,
    min_facet_angle=25,
    max_radius_surface_delaunay_ball=0.1,
    max_facet_distance=0.001,
    verbose=False,
)

Also, the third argument to meshio.write_points_cells defines the faces, and faces must have a shape (n_faces, 3) for triangular elements. In my case I was using a mesh in a different format from another package (PyVista) and had to reshape the arrays.

Upvotes: 4

John
John

Reputation: 71

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

Upvotes: 3

Nico Schlömer
Nico Schlömer

Reputation: 58791

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

Upvotes: 2

Darren Engwirda
Darren Engwirda

Reputation: 7015

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.

Upvotes: 2

Ripi2
Ripi2

Reputation: 7198

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.

Upvotes: 0

Related Questions