Reputation: 3657
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.
Upvotes: 5
Views: 4785
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
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
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
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
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:
z
can be found traversing the mesh until you find a triangle that contains x,y
coordinates 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