import math
import os
from functools import lru_cache
from pathlib import Path
import numpy as np
import coacd
import igl
import Imath
import OpenEXR
import tetgen
import trimesh
from PIL import Image
import genesis as gs
from . import geom as gu
from .cache import (
OwnedCache,
SizeCappedCache,
get_cvx_cache,
get_exr_path,
get_hashkey,
get_wt_cache,
get_wth_cache,
register_cache_clear,
)
from .misc import get_assets_dir, get_src_dir
MESH_REPAIR_ERROR_THRESHOLD = 0.01
# Vertex welding quantum: co-located vertices are authored with identical coordinates, so a tight absolute quantum
# absorbs sub-nanometer parsing noise while keeping genuinely distinct vertices separate.
VERT_WELD_QUANTIZE_FACTOR = 1e-8
Y_UP_TRANSFORM = np.asarray( # translation on the bottom row
[[1.0, 0.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0], [0.0, -1.0, 0.0, 0.0], [0.0, 0.0, 0.0, 1.0]], dtype=np.float32
)
DEFAULT_PLANE_TEXTURE_PATH = "textures/checker.png" # use checkerboard texture by default
# The GPU interpolates the attributes of a triangle in low precision, so a triangle hundreds of meters wide makes the
# edges of its texture flicker and shift under a rounding change of the camera pose. Cells of this size keep them still.
PLANE_MAX_CELL_SIZE = 50.0
[docs]def color_f32_to_u8(color) -> np.ndarray:
return np.round(np.asarray(color, dtype=np.float32) * 255.0).astype(np.uint8)
[docs]def color_u8_to_f32(color) -> np.ndarray:
return np.asarray(color, dtype=np.uint8).astype(np.float32) / 255.0
[docs]def glossiness_to_roughness(glossiness: float) -> float:
return (2 / (glossiness + 2)) ** (1.0 / 4.0)
def get_signed_volume(verts: np.ndarray, faces: np.ndarray) -> float:
"""Compute the volume a closed triangle mesh encloses, positive when its faces wind outward.
'verts' holds the positions of the vertices and 'faces' the indices of the vertices of each triangle. The volume is
the sum of the signed volumes of the tetrahedra each face spans with the centroid of the vertices (divergence
theorem).
"""
# The volume is computed alone, which spares the division by the volume that the center of mass of a mesh
# enclosing no volume would require. The vertices are centered first: the volume of a closed mesh does not depend on
# its position, while the products of coordinates far from the origin cancel at the scale of those coordinates.
tris = (verts - verts.mean(axis=0))[faces]
return np.einsum("ij,ij->", tris[:, 0], np.cross(tris[:, 1], tris[:, 2])) / 6.0
def convex_hull(tmesh: trimesh.Trimesh) -> trimesh.Trimesh:
"""Compute the convex hull of a mesh, with its faces wound outward unless the vertices are coplanar.
The winding of a hull is repaired from the sign of the volume it encloses, which a flat hull does not have, so a
hull of coplanar vertices keeps the winding of each face as qhull computed it. The vertices are coplanar when the
thickness of their cloud relative to its extent is below the tolerance trimesh merges vertices with, since the
rounding of the coordinates alone already spreads the vertices of a flat mesh above machine precision.
"""
# Fewer than 4 vertices enclose no volume
is_flat = len(tmesh.vertices) < 4
if not is_flat:
singular_values = np.linalg.svd(tmesh.vertices - tmesh.vertices.mean(axis=0), compute_uv=False)
is_flat = singular_values[-1] <= trimesh.tol.merge * singular_values[0]
return trimesh.convex.convex_hull(tmesh, repair=not is_flat)
def get_wall_thickness(verts: np.ndarray, faces: np.ndarray, quantile: float = 0.25) -> float:
"""Estimate a watertight mesh's characteristic wall thickness by probing its local diameter with inward rays.
A ray is cast inward along the face normal from a stride-subsampled set of face centroids to the first opposite
surface. The `quantile` of those hit distances, weighted by probed face area so the estimate measures surface
rather than tessellation density (a thin wall spanned by a handful of large faces must not be outvoted by many
small decorative facets), is returned: a low quantile approximates the thinnest wall and the median the typical
one. The mesh must be watertight: inward rays escape through holes of an open surface, so the estimate is
meaningless there and an exception is raised (the caller is expected to skip the probe for non-watertight meshes).
The estimate is deliberately a scalar, not per-axis: the SDF grid it sizes does not only resolve walls along
their normal - it certifies contact penetrations through Lipschitz cone bounds whose slack grows with the
lattice spacing in every direction around the contact point, tangent axes included. On a closed shell each
wall's tangent directions are other walls' normals (a mug's vertical wall lies tangent to the vertical axis
even though the only wall facing that axis, the thick bottom, would justify coarse vertical cells), so relaxing
any one axis to the thickness of the walls facing it measurably degrades the certified pens of every wall
tangent to it, and no bound-side search can recover the loss (the lateral sample offset is a lattice property).
"""
cache = get_wth_cache(verts, faces, quantile)
wall_thickness = cache.load()
if wall_thickness is not None:
return wall_thickness
mesh = trimesh.Trimesh(verts, faces, process=False)
if not mesh.is_watertight:
gs.raise_exception("Wall-thickness estimation requires a watertight mesh.")
# A face of no area has no normal to cast a ray along, and a ray hitting it divides by its area, so the faces of no
# area are left out of the probe, the surface they cover being none
mesh.update_faces(mesh.nondegenerate_faces())
diag = np.linalg.norm(verts.max(axis=0) - verts.min(axis=0))
stride = max(1, len(faces) // 1000)
centers = mesh.triangles_center[::stride]
normals = mesh.face_normals[::stride]
origins = centers - 1e-3 * diag * normals
locations, ray_idx, _ = mesh.ray.intersects_location(origins, -normals, multiple_hits=False)
thickness = np.linalg.norm(locations - origins[ray_idx], axis=1)
# A hit at the ray origin's own face reads as zero thickness; drop those self-hits before taking the quantile.
is_hit_valid = thickness > 1e-4 * diag
thickness = thickness[is_hit_valid]
areas = mesh.area_faces[::stride][ray_idx][is_hit_valid]
order = np.argsort(thickness)
areas_cum = np.cumsum(areas[order])
wall_thickness = thickness[order][np.searchsorted(areas_cum, quantile * areas_cum[-1])]
cache.save(wall_thickness)
return wall_thickness
[docs]class MeshInfo:
def __init__(self):
self.surface = None
self.metadata = {}
self.verts = []
self.faces = []
self.normals = []
self.uvs = []
self.n_points = 0
def set_property(self, surface=None, metadata=None):
self.surface = surface
self.metadata = metadata
def append(self, verts, faces, normals, uvs):
faces += self.n_points
self.verts.append(verts)
self.faces.append(faces)
self.normals.append(normals)
self.uvs.append(uvs)
self.n_points += len(verts)
def export_mesh(self, scale: float, is_mesh_zup: bool) -> "gs.Mesh":
uvs = None
if self.uvs:
for i, (uvs, verts) in enumerate(zip(self.uvs, self.verts)):
if uvs is None:
self.uvs[i] = np.zeros((len(verts), 2), dtype=gs.np_float)
uvs = np.concatenate(self.uvs, axis=0)
verts = np.concatenate(self.verts, axis=0)
faces = np.concatenate(self.faces, axis=0)
normals = np.concatenate(self.normals, axis=0)
return gs.Mesh.from_attrs(
verts=verts,
faces=faces,
normals=normals,
surface=self.surface,
uvs=uvs,
scale=scale,
metadata=self.metadata,
is_mesh_zup=is_mesh_zup,
)
[docs]class MeshInfoGroup:
def __init__(self):
self.infos: dict[str, MeshInfo] = {}
def get(self, name: str):
first_created = False
mesh_info = self.infos.get(name)
if mesh_info is None:
mesh_info = self.infos.setdefault(name, MeshInfo())
first_created = True
return mesh_info, first_created
def export_meshes(self, scale, is_mesh_zup) -> "list[gs.Mesh]":
return [mesh_info.export_mesh(scale, is_mesh_zup) for mesh_info in self.infos.values()]
def get_asset_path(file):
return os.path.join(get_src_dir(), "assets", file)
[docs]def load_mesh(file):
if isinstance(file, (str, Path)):
try:
return trimesh.load_mesh(file, force="mesh", skip_texture=False)
except Exception as e:
gs.logger.warning(f"Failed to load mesh with texture: {e}")
# try loading without texture data
return trimesh.load_mesh(file, force="mesh", skip_texture=True)
return file
def compute_sdf_data(mesh, res):
"""
Convert mesh to sdf voxels and a transformation matrix from mesh frame to voxel frame.
"""
voxels_radius = 0.6
x = np.linspace(-voxels_radius, voxels_radius, res)
y = np.linspace(-voxels_radius, voxels_radius, res)
z = np.linspace(-voxels_radius, voxels_radius, res)
X, Y, Z = np.meshgrid(x, y, z, indexing="ij")
query_points = np.stack([X, Y, Z], axis=-1).reshape((-1, 3))
voxels, *_ = igl.signed_distance(query_points, mesh.vertices, mesh.faces)
voxels = voxels.reshape((res, res, res)).astype(gs.np_float, copy=False)
T_mesh_to_sdf = np.eye(4, dtype=gs.np_float)
T_mesh_to_sdf[:3, :3] *= (res - 1) / (voxels_radius * 2)
T_mesh_to_sdf[:3, 3] = (res - 1) / 2
sdf_data = {
"voxels": voxels,
"T_mesh_to_sdf": T_mesh_to_sdf,
}
return sdf_data
def surface_uvs_to_trimesh_visual(surface, uvs=None, n_verts=None):
texture = surface.get_rgba()
# 'trimesh' uses uvs starting from the top-left corner, so flip them to Genesis' convention. The flip runs at the
# precision a mesh stores its uvs at, so a mesh drawn from uvs read back from a file matches the one it was built as.
flipped_uvs = None
if uvs is not None:
flipped_uvs = uvs.astype(gs.np_float)
flipped_uvs[:, 1] = 1.0 - flipped_uvs[:, 1]
# Composite emissive additively on top of the base color, but only when the base color is the packed albedo
# (present, nonblack, and distinct from the emissive). Otherwise get_rgba already returned the emissive as the
# albedo and re-adding it would double it. Baking it into the material makes every renderer path pick it up through
# from_trimesh (rigid and deformable alike), and forces a PBRMaterial since SimpleMaterial has no emissive channel.
# An image emissive samples UVs, so it is composited only when the mesh has them (else the shader references an
# undeclared uv_0); a flat emissive color needs none.
base = surface.texture
emission = surface.emission
emissive_kwargs = {}
if emission is not None and base is not None and base is not emission and not base.is_black:
if (
isinstance(emission, gs.textures.ImageTexture)
and emission.image_array is not None
and flipped_uvs is not None
):
emissive_kwargs = dict(
emissiveTexture=Image.fromarray(emission.image_array), emissiveFactor=emission.image_color
)
elif isinstance(emission, gs.textures.ColorTexture):
emissive_kwargs = dict(emissiveFactor=emission.color)
if isinstance(texture, gs.textures.ImageTexture):
if flipped_uvs is not None:
assert texture.image_array.dtype == np.uint8
image = Image.fromarray(texture.image_array)
if emissive_kwargs:
material = trimesh.visual.material.PBRMaterial(
baseColorTexture=image, baseColorFactor=(255, 255, 255, 255), **emissive_kwargs
)
else:
material = trimesh.visual.material.SimpleMaterial(image=image, diffuse=(1.0, 1.0, 1.0, 1.0))
visual = trimesh.visual.TextureVisuals(uv=flipped_uvs, material=material)
else:
# fall back to color texture
visual = trimesh.visual.ColorVisuals(vertex_colors=np.tile(texture.mean_color, [n_verts, 1]))
elif isinstance(texture, gs.textures.ColorTexture):
if n_verts is None:
gs.raise_exception("n_verts is required for color texture.")
if emissive_kwargs:
# The flat base color applies through a factor, but an image emissive still samples the UVs, so pass them.
material = trimesh.visual.material.PBRMaterial(
baseColorFactor=color_f32_to_u8(texture.color), **emissive_kwargs
)
visual = trimesh.visual.TextureVisuals(uv=flipped_uvs, material=material)
else:
visual = trimesh.visual.ColorVisuals(vertex_colors=np.tile(np.array(texture.color), [n_verts, 1]))
assert visual.defined
else:
gs.raise_exception("Cannot get texture when generating trimesh visual.")
return visual
[docs]def convex_decompose(mesh, coacd_options):
# The decomposition of a scaled mesh is the decomposition of the mesh scaled, so the cache is keyed on vertices
# normalized by the mesh scale and the hulls are rescaled on the way out: every scaled copy of an asset shares a
# single entry.
mesh_scale = float(np.linalg.norm(mesh.extents))
assert not (np.isinf(mesh_scale) or np.isnan(mesh_scale) or mesh_scale <= 0.0)
cache = get_cvx_cache(mesh.vertices / mesh_scale, mesh.faces, coacd_options)
# loading pre-computed cache if available
is_cached_loaded = False
loaded_cache = cache.load()
if loaded_cache is not None:
gs.logger.debug("Convex decomposition file (.cvx) found in cache.")
mesh_parts = loaded_cache["mesh_parts"]
cached_mesh_scale = loaded_cache["mesh_scale"]
# rescale loaded mesh parts
if not (np.isinf(cached_mesh_scale) or np.isnan(cached_mesh_scale) or cached_mesh_scale <= 0.0):
rescale_factor = mesh_scale / cached_mesh_scale
for mesh_part in mesh_parts:
mesh_part.vertices *= rescale_factor
is_cached_loaded = True
if not is_cached_loaded:
with gs.logger.timer("Running convex decomposition."):
mesh = coacd.Mesh(mesh.vertices, mesh.faces)
args = coacd_options
result = coacd.run_coacd(
mesh,
threshold=args.threshold,
max_convex_hull=args.max_convex_hull,
preprocess_mode=args.preprocess_mode,
preprocess_resolution=args.preprocess_resolution,
resolution=args.resolution,
mcts_nodes=args.mcts_nodes,
mcts_iterations=args.mcts_iterations,
mcts_max_depth=args.mcts_max_depth,
pca=args.pca,
merge=args.merge,
decimate=args.decimate,
max_ch_vertex=args.max_ch_vertex,
extrude=args.extrude,
extrude_margin=args.extrude_margin,
apx_mode=args.apx_mode,
seed=args.seed,
)
mesh_parts = []
for vs, fs in result:
mesh_parts.append(trimesh.Trimesh(vs, fs))
cache.save({"mesh_parts": mesh_parts, "mesh_scale": mesh_scale})
return mesh_parts
# 512 MiB of processed collision geometry. Sized by the geometry footprint actually retained (vertices and faces of
# the cached meshes), so a few large assets or many small ones are both bounded without an arbitrary entry count.
_COLLISION_GEOMS_CACHE = SizeCappedCache(max_bytes=512 * 1024 * 1024)
def postprocess_collision_geoms(
g_infos,
decimate,
decimate_face_num,
decimate_aggressiveness,
convexify,
decompose_error_threshold,
coacd_options,
watertighten,
):
# Convexification / decomposition of a collision mesh (convex hull and coacd decomposition) is the dominant cost of
# adding a file-based entity, and a scene built from many copies of the same asset would otherwise redo it for every
# entity. The result is memoized on the geometry of the input collision meshes and on every option / per-geom field
# that can change the outcome. The processed meshes are immutable (their vertices are only read downstream,
# alignment is folded into the link frame), so the cached ones are shared across entities, which also collapses
# their memory footprint. Only fresh g_info dicts are handed back so the caller can re-express the per-geom pose in
# place without corrupting the template.
if not g_infos:
return []
# Keyed on the vertices themselves rather than matched within a tolerance the way the on-disk caches are: one
# entry per distinct input keeps every variant of a same topology evictable on its own, which is what bounds the
# footprint of this cache. Re-exports of one asset are what tolerant matching buys, and they cost nothing here
# since the expensive steps this memoizes hit their own tolerant on-disk cache anyway.
key_parts = [
bool(decimate),
int(decimate_face_num),
int(decimate_aggressiveness),
bool(convexify),
float(decompose_error_threshold),
-1 if watertighten is None else int(watertighten),
coacd_options.model_dump(),
]
for g_info in g_infos:
mesh = g_info["mesh"]
geom_type = g_info.get("type")
friction = g_info.get("friction")
density = g_info.get("density")
sol_params = g_info.get("sol_params")
key_parts += [
np.ascontiguousarray(mesh.verts),
np.ascontiguousarray(mesh.faces),
-1 if geom_type is None else int(geom_type),
int(g_info.get("contype", 0)),
int(g_info.get("conaffinity", 0)),
float("nan") if friction is None else float(friction),
float("nan") if density is None else float(density),
np.zeros(0) if sol_params is None else np.ascontiguousarray(sol_params, dtype=np.float64),
np.ascontiguousarray(g_info.get("pos", gu.zero_pos()), dtype=np.float64),
np.ascontiguousarray(g_info.get("quat", gu.identity_quat()), dtype=np.float64),
]
key = get_hashkey(*key_parts)
cached = _COLLISION_GEOMS_CACHE.get(key)
if cached is None:
cached = _postprocess_collision_geoms_impl(
g_infos,
decimate,
decimate_face_num,
decimate_aggressiveness,
convexify,
decompose_error_threshold,
coacd_options,
watertighten,
)
# The normals the geoms built from the cache share are computed on it, the vertex normals from the face
# normals, so both count in what it holds
n_bytes = sum(
g_info["mesh"].verts.nbytes
+ g_info["mesh"].faces.nbytes
+ g_info["mesh"].trimesh.face_normals.nbytes
+ g_info["mesh"].trimesh.vertex_normals.nbytes
for g_info in cached
)
_COLLISION_GEOMS_CACHE.put(key, cached, n_bytes)
# Hand back per-entity geoms that share the arrays of the cached geometry, and the data derived from it since they
# have its hash (see 'Mesh.tmesh_hash'), but own their trimesh and surface. This keeps appearance per-entity (each
# entity gets an independent, e.g. randomized, collision color) while the heavy geometry is computed once.
result = []
for g_info in cached:
template = g_info["mesh"]
template_tmesh = template.trimesh
mesh = gs.Mesh(
mesh=trimesh.Trimesh(
vertices=template_tmesh.vertices,
faces=template_tmesh.faces,
vertex_normals=template_tmesh.vertex_normals,
process=False,
),
surface=gs.surfaces.Collision(),
uvs=template.uvs,
metadata=template.metadata.copy(),
)
result.append({**g_info, "mesh": mesh})
return result
def _postprocess_collision_geoms_impl(
g_infos,
decimate,
decimate_face_num,
decimate_aggressiveness,
convexify,
decompose_error_threshold,
coacd_options,
watertighten,
):
# Early return if there is no geometry to process
if not g_infos:
return []
# Check whether the geometries are authored, ie they are all watertight and convex, as a cheap proxy
is_authored = all(
(tmesh := g_info["mesh"].trimesh).is_winding_consistent and tmesh.is_watertight and tmesh.is_convex
for g_info in g_infos
)
# Weld coincident vertices of non-convex collision meshes onto a separate copy. Formats that store unshared
# per-face vertices (notably STL) yield a vertex soup whose duplicates differ only by float rounding; the
# downstream exact dedup at geom build only partially fuses them, leaving a degraded (sliver) mesh that corrupts
# the SDF. The convex path skips this (the hull / decomposition replaces the surface anyway). The collision mesh
# may be the very same trimesh as the visual geom, so the weld must not be done in place: only the collision geom
# is swapped for the welded copy, leaving the visual vertices untouched. Already-shared meshes (OBJ, glTF) weld to
# no fewer vertices and keep their original geom.
if not is_authored and not convexify:
for g_info in g_infos:
if g_info["type"] != gs.GEOM_TYPE.MESH:
continue
welded = g_info["mesh"].trimesh.copy()
welded.merge_vertices()
if len(welded.vertices) < len(g_info["mesh"].trimesh.vertices):
g_info["mesh"] = gs.Mesh.from_trimesh(
mesh=welded,
surface=gs.surfaces.Collision(),
metadata=g_info["mesh"].metadata.copy(),
)
# Try the repair meshes that seems to be "broken" but not beyond repair.
# Note that this procedure is only applied if the estimated volume is significantly different before and after
# repair, to avoid altering the original mesh without actual benefit. Moreover, only duplicate faces are removed,
# which is less aggressive than `Trimesh.process(validate=True)`.
if not is_authored:
for g_info in g_infos:
mesh = g_info["mesh"]
tmesh = mesh.trimesh
if g_info["type"] != gs.GEOM_TYPE.MESH:
continue
if tmesh.is_winding_consistent and not tmesh.is_watertight:
volume_repaired = get_signed_volume(tmesh.vertices, tmesh.faces[tmesh.unique_faces()])
if abs(volume_repaired) < gs.EPS:
continue
volume = get_signed_volume(tmesh.vertices, tmesh.faces)
if abs(abs(volume / volume_repaired) - 1.0) > MESH_REPAIR_ERROR_THRESHOLD:
gs.logger.info(
"Collision mesh is not watertight and has ill-defined volume. It will be repaired by removing "
"duplicate faces."
)
tmesh.update_faces(tmesh.unique_faces())
tmesh._cache.clear()
tmesh.visual._cache.clear()
# Check which geometries can be convexified without decomposition
geoms_must_decompose = [False] * len(g_infos)
volume_err_max = 0.0
if not is_authored and convexify:
for i_g, g_info in enumerate(g_infos):
mesh = g_info["mesh"]
tmesh = mesh.trimesh
# Skip geometries that do not corresponds to mesh or have no enclosed volume
if g_info["type"] != gs.GEOM_TYPE.MESH:
continue
cmesh = convex_hull(tmesh)
cmesh_volume = get_signed_volume(cmesh.vertices, cmesh.faces)
if abs(cmesh_volume) < gs.EPS:
continue
# Fix mesh temporarily to make volume computation more reliable
if not tmesh.is_winding_consistent:
tmesh = tmesh.copy()
tmesh.process(validate=True)
# Compute volume approximation error between true geometry and its convex hull conservatively. The volume
# of a mesh wound inward is negative.
if not tmesh.is_winding_consistent:
if not math.isinf(decompose_error_threshold):
geoms_must_decompose[i_g] = True
volume_err_max = float("inf")
elif abs(volume := get_signed_volume(tmesh.vertices, tmesh.faces)) > gs.EPS:
volume_err = abs(cmesh_volume / abs(volume) - 1.0)
if volume_err > decompose_error_threshold:
geoms_must_decompose[i_g] = True
volume_err_max = max(volume_err_max, volume_err)
must_decompose = any(geoms_must_decompose)
# Fuse collision geoms that share a collision group and contact parameters into one closed surface per consistent
# sub-group, transforming each sub-mesh into its group's first geom frame so the fused geom keeps that pose. One
# surface per group exposes the whole-object topology needed to detect enclosure correctly and to estimate wall
# thickness (which sets the grid SDF cell size, hence collision accuracy), and it is cheaper - fewer vertices to
# scan and roughly one contact set per link instead of per geom. All geom types can be merged except planes and
# terrains. On the nonconvex path, sub-groups are systematically fused unless watertightening is disabled; on the
# convex path, only when one of their members requires decomposition.
geoms_is_fused = [False] * len(g_infos)
if len(g_infos) > 1 and ((not convexify and watertighten is not None) or (not is_authored and must_decompose)):
fusion_groups: list[list[int]] = []
for i, g_info in enumerate(g_infos):
# Join the first existing sub-group with a matching collision group and contact params, else open a new one.
if g_info["type"] not in (gs.GEOM_TYPE.PLANE, gs.GEOM_TYPE.TERRAIN):
for fusion_group in fusion_groups:
first_g_info = g_infos[fusion_group[0]]
if (
first_g_info["type"] not in (gs.GEOM_TYPE.PLANE, gs.GEOM_TYPE.TERRAIN)
and all(first_g_info.get(name) == g_info.get(name) for name in ("contype", "conaffinity"))
and all(
np.allclose(first_g_info.get(name, np.nan), g_info.get(name, np.nan), equal_nan=True)
# A fused mesh carries a single density, so geoms with different authored densities must
# stay separate for the density-weighted link inertial to survive fusion.
for name in ("friction", "sol_params", "density")
)
):
fusion_group.append(i)
break
else:
fusion_groups.append([i])
else:
fusion_groups.append([i])
fused_infos = []
geoms_is_fused = []
for fusion_group in fusion_groups:
if len(fusion_group) == 1 or (convexify and not any(geoms_must_decompose[i] for i in fusion_group)):
fused_infos.extend(g_infos[i] for i in fusion_group)
geoms_is_fused.extend([False] * len(fusion_group))
continue
first_g_info = g_infos[fusion_group[0]]
T_first_inv = np.linalg.inv(
gu.trans_quat_to_T(first_g_info.get("pos", gu.zero_pos()), first_g_info.get("quat", gu.identity_quat()))
)
tmeshes = []
metadata = set(first_g_info["mesh"].metadata.items())
for i in fusion_group:
g_info = g_infos[i]
mesh = g_info["mesh"]
tmesh = mesh.trimesh
T_rel = T_first_inv @ gu.trans_quat_to_T(
g_info.get("pos", gu.zero_pos()), g_info.get("quat", gu.identity_quat())
)
if not np.allclose(T_rel, np.eye(4)):
tmesh = tmesh.copy()
tmesh.apply_transform(T_rel)
metadata &= set(mesh.metadata.items())
tmeshes.append(tmesh)
tmesh = trimesh.util.concatenate(tmeshes)
mesh = gs.Mesh.from_trimesh(
mesh=tmesh, surface=gs.surfaces.Collision(), metadata=dict(metadata) | {"merged": True}
)
fused_infos.append({**first_g_info, **dict(type=gs.GEOM_TYPE.MESH, data=None, mesh=mesh)})
geoms_is_fused.append(True)
g_infos = fused_infos
# A genuinely fused mesh (a sub-group of more than one geom) re-enters the pipeline so its convex hull /
# decomposition is recomputed on the union; requiring an actual merge stops infinite recursion once every
# sub-group is either a singleton or accurate enough to skip decomposition.
if convexify and any(geoms_is_fused):
return _postprocess_collision_geoms_impl(
g_infos,
decimate,
decimate_face_num,
decimate_aggressiveness,
convexify,
decompose_error_threshold,
coacd_options,
watertighten,
)
# Nonconvex: watertighten each fused surface into a closed mesh so the grid SDF is reliable. The convex path skips
# this (its hull / decomposition replaces the surface anyway, so an alpha-wrap would be wasted work).
if not convexify and watertighten is not None:
from .watertighten import watertighten_mesh
for g_info, is_fused in zip(g_infos, geoms_is_fused):
# Fused geoms are always watertightened, as their sub-meshes may overlap while being individually
# watertight. Other geoms are skipped if they are not generic meshes or already watertight or convex.
tmesh = g_info["mesh"].trimesh
if not is_fused and (g_info["type"] != gs.GEOM_TYPE.MESH or tmesh.is_watertight or tmesh.is_convex):
continue
# On-disk cache keyed by (vertices, faces, aggressiveness): a repeated build on the same geom is a file read
# instead of a multi-second SDF + DC + QEM rebuild.
cache = get_wt_cache(tmesh.vertices, tmesh.faces, watertighten)
cached_mesh = cache.load()
if cached_mesh is None:
cached_mesh = watertighten_mesh(tmesh.vertices, tmesh.faces, aggressiveness=watertighten)
cache.save(cached_mesh)
v_out, f_out = cached_mesh
fused = trimesh.Trimesh(vertices=v_out, faces=f_out, process=False)
metadata = g_info["mesh"].metadata.copy()
metadata["watertightened"] = True
g_info["mesh"] = gs.Mesh.from_trimesh(mesh=fused, surface=gs.surfaces.Collision(), metadata=metadata)
if not is_authored and must_decompose:
if math.isinf(volume_err_max):
gs.logger.info(
"Collision mesh has inconsistent winding and 'decompose_error_threshold' != float('inf'). "
"Falling back to more expensive convex decomposition (see FileMorph options)."
)
else:
gs.logger.info(
f"Convex hull is not accurate enough for collision detection ({volume_err_max:.3f}). Falling back to "
"more expensive convex decomposition (see FileMorph options)."
)
_g_infos = []
for g_info in g_infos:
mesh = g_info["mesh"]
tmesh = mesh.trimesh
if g_info["type"] != gs.GEOM_TYPE.MESH:
volume_err = 0.0
elif not tmesh.is_winding_consistent:
volume_err = float("inf")
elif abs(volume := get_signed_volume(tmesh.vertices, tmesh.faces)) < gs.EPS:
volume_err = 0.0
else:
cmesh = convex_hull(tmesh)
volume_err = abs(get_signed_volume(cmesh.vertices, cmesh.faces) / abs(volume) - 1.0)
if volume_err > decompose_error_threshold: # Note that 'inf' is not larger than 'inf'
tmeshes = convex_decompose(tmesh, coacd_options)
meshes = [
gs.Mesh.from_trimesh(
tmesh, surface=gs.surfaces.Collision(), metadata={**mesh.metadata, "decomposed": True}
)
for tmesh in tmeshes
]
_g_infos += [{**g_info, **dict(mesh=mesh)} for mesh in meshes]
else:
_g_infos.append(g_info)
g_infos = _g_infos
# Process of meshes sequentially
_g_infos = []
for g_info in g_infos:
mesh = g_info["mesh"]
tmesh = mesh.trimesh
num_faces = len(tmesh.faces)
if not decimate and num_faces > 5000:
gs.logger.warning(
f"At least one of the meshes contain many faces ({num_faces}). Consider setting "
"'morph.decimate=True' to speed up collision detection and improve numerical stability."
)
if decimate and decimate_face_num < 100:
gs.logger.warning(
"`decimate_face_num` should be greater than 100 to ensure sufficient geometry details are preserved."
)
# Decimation is an independent step applied after watertightening (which runs its own internal QEM): when
# 'decimate' is requested it simplifies the collision mesh - including a watertighten wrap - down to
# 'decimate_face_num'. It is gated on the mesh being either above that target or watertight, since decimating
# a low-poly non-watertight mesh would be unreliable.
must_decimate = num_faces > decimate_face_num or tmesh.is_watertight
if not must_decimate:
gs.logger.debug(
"Collision mesh is not watertight. Decimate would be unreliable. Skipping as mesh is already low-poly."
)
mesh = gs.Mesh.from_trimesh(
mesh=tmesh,
convexify=convexify,
decimate=decimate and must_decimate,
decimate_face_num=decimate_face_num,
decimate_aggressiveness=decimate_aggressiveness,
surface=gs.surfaces.Collision(),
metadata=mesh.metadata.copy(),
)
_g_infos.append({**g_info, **dict(mesh=mesh)})
return _g_infos
@lru_cache(maxsize=32)
def _load_trimesh_scene_cached(path, group_by_material, mtime) -> "trimesh.Scene":
# Parsing a mesh file from disk (read + decode) is by far the dominant cost of adding a file-based entity, and a
# scene built from many copies of the same asset (e.g. a grid of identical objects) would otherwise re-parse the
# exact same bytes for every entity. The result is keyed by file modification time so an edited asset is reloaded.
# The returned scene is treated as an immutable template - callers copy each geometry out before applying scale or
# surface - so it is never mutated and is safe to share across entities.
return trimesh.load(path, force="scene", group_material=group_by_material, process=False)
register_cache_clear(_load_trimesh_scene_cached.cache_clear)
# Meshes of a file brought to the z-up frame and scaled, kept as long as any mesh built from them lives
_FILE_TMESHES_CACHE = OwnedCache()
def _load_trimesh_geometries(path, group_by_material, scale, is_mesh_zup, mtime) -> "tuple[trimesh.Trimesh, ...]":
"""Load the meshes of a file, brought to the z-up frame and scaled.
Their arrays are read-only, the meshes built from them sharing them (see 'Mesh.from_trimesh'), so that writing
into them raises rather than changing every one of those meshes.
"""
scene = _load_trimesh_scene_cached(path, group_by_material, mtime)
tmeshes = []
for tmesh in scene.geometry.values():
if not isinstance(tmesh, trimesh.Trimesh):
gs.raise_exception(f"Mesh type not supported: {path}")
tmesh = tmesh.copy(include_cache=True)
if not is_mesh_zup:
tmesh.apply_transform(Y_UP_TRANSFORM.T)
tmesh.apply_scale(np.atleast_1d(scale))
# The transforms leave the normals stored in the file strided, which every trimesh wrapping them would copy back
# to contiguous (see 'Mesh.from_trimesh' for the cache read)
face_normals, vertex_normals = tmesh._cache["face_normals"], tmesh._cache["vertex_normals"]
if face_normals is not None:
tmesh.face_normals = np.ascontiguousarray(face_normals)
if vertex_normals is not None:
tmesh.vertex_normals = np.ascontiguousarray(vertex_normals)
for array in (tmesh.vertices, tmesh.faces, tmesh._cache["face_normals"], tmesh._cache["vertex_normals"]):
if array is not None:
array.flags.writeable = False
tmeshes.append(tmesh)
return tuple(tmeshes)
def parse_mesh_trimesh(path, group_by_material, scale, is_mesh_zup, surface) -> "list[gs.Mesh]":
# The meshes of a file are processed once for all the entities instantiating it at one scale, which share them
key = (path, group_by_material, scale, is_mesh_zup, os.path.getmtime(path))
tmeshes = _FILE_TMESHES_CACHE.find(key)
if tmeshes is None:
tmeshes = _load_trimesh_geometries(*key)
meshes = [
gs.Mesh.from_trimesh(
mesh=tmesh,
surface=surface,
metadata={"mesh_path": path, "imported_as_zup": is_mesh_zup},
share_arrays=True,
)
for tmesh in tmeshes
]
for mesh in meshes:
_FILE_TMESHES_CACHE.own(mesh, key, tmeshes)
return meshes
def trimesh_to_mesh(mesh, scale, surface) -> "gs.Mesh":
return gs.Mesh.from_trimesh(mesh=mesh, scale=scale, surface=surface)
def adjust_alpha_cutoff(alpha_cutoff, alpha_mode):
if alpha_mode == 0: # OPAQUE
return 0.0
if alpha_mode == 1: # MASK
return alpha_cutoff
return None # BLEND
def PIL_to_array(image):
# A paletted image stores indices into its palette, so its texels are the colors those indices stand for
if isinstance(image, Image.Image) and image.mode == "P":
image = image.convert("RGBA" if "transparency" in image.info else "RGB")
return np.array(image)
def tonemapped(image):
exposure = 0.5
return (np.clip(np.power(image / 255 * np.power(2, exposure), 1 / 2.2), 0, 1) * 255).astype(np.uint8)
def create_texture(image, factor, encoding):
if image is not None:
return gs.textures.ImageTexture(image_array=image, image_color=factor, encoding=encoding)
if factor is not None:
return gs.textures.ColorTexture(color=factor)
return None
def apply_transform(transform, positions, normals=None):
# Note that here transform's translation is on the bottom row, different from that in Genesis and trimesh.
transformed_positions = (np.column_stack([positions, np.ones(len(positions))]) @ transform)[:, :3]
transformed_normals = normals
if normals is not None:
lin_mat = transform[:3, :3]
if not np.allclose(lin_mat, np.identity(3), atol=gs.EPS): # has rotation or scaling
# A normal is a covector, so it maps through the cofactor matrix of the linear part, whose rows are the
# cross products of the other two rows. That matrix sends the cross product of two edges to the cross
# product of the transformed edges, keeping a normal perpendicular to its own triangle under any scaling.
cofactor_mat = np.cross(np.roll(lin_mat, -1, axis=0), np.roll(lin_mat, -2, axis=0))
transformed_normals = normals @ cofactor_mat
transformed_normals /= np.linalg.norm(transformed_normals, axis=1, keepdims=True)
return transformed_positions, transformed_normals
[docs]def create_frame(
origin_radius=0.012, axis_radius=0.005, axis_length=1.0, head_radius=0.01, head_length=0.03, sections=12
):
origin = create_sphere(radius=origin_radius, subdivisions=2)
x = create_arrow(
length=axis_length,
radius=axis_radius,
l_ratio=head_length / axis_length,
r_ratio=head_radius / axis_radius,
body_color=(0.7, 0.0, 0.0, 1.0),
head_color=(0.7, 0.7, 0.7, 1.0),
sections=sections,
)
y = create_arrow(
length=axis_length,
radius=axis_radius,
l_ratio=head_length / axis_length,
r_ratio=head_radius / axis_radius,
body_color=(0.0, 0.7, 0.0, 1.0),
head_color=(0.7, 0.7, 0.7, 1.0),
sections=sections,
)
z = create_arrow(
length=axis_length,
radius=axis_radius,
l_ratio=head_length / axis_length,
r_ratio=head_radius / axis_radius,
body_color=(0.0, 0.0, 0.7, 1.0),
head_color=(0.7, 0.7, 0.7, 1.0),
sections=sections,
)
x.vertices = gu.transform_by_R(x.vertices, gu.euler_to_R((0.0, 90.0, 0.0)))
y.vertices = gu.transform_by_R(y.vertices, gu.euler_to_R((-90.0, 0.0, 0.0)))
return trimesh.util.concatenate([origin, x, y, z])
[docs]def create_camera_frustum(camera, color):
# camera
camera_mesh = trimesh.load(os.path.join(get_src_dir(), "assets", "meshes", "camera/camera.glb"), force="mesh")
camera_mesh.visual = camera_mesh.visual.to_color()
camera_mesh.apply_translation([0.0, 0.0, 1.0])
camera_mesh.apply_scale(0.05)
# frustum
near_half_height = camera.near * np.tan(np.deg2rad(camera.fov / 2))
near_half_width = near_half_height * camera.aspect_ratio
far_half_height = camera.far * np.tan(np.deg2rad(camera.fov / 2))
far_half_width = far_half_height * camera.aspect_ratio
# Define the vertices of the frustum
vertices = np.array(
[
[0, 0, 0], # apex
[-near_half_width, -near_half_height, -camera.near], # near bottom left
[near_half_width, -near_half_height, -camera.near], # near bottom right
[near_half_width, near_half_height, -camera.near], # near top right
[-near_half_width, near_half_height, -camera.near], # near top left
[-far_half_width, -far_half_height, -camera.far], # far bottom left
[far_half_width, -far_half_height, -camera.far], # far bottom right
[far_half_width, far_half_height, -camera.far], # far top right
[-far_half_width, far_half_height, -camera.far], # far top left
]
)
# Define the faces of the frustum
faces = np.array(
[
# # near face
# [1, 2, 3, 4],
# # far face
# [5, 6, 7, 8],
# side face
[2, 1, 5, 6],
[3, 2, 6, 7],
[4, 3, 7, 8],
[1, 4, 8, 5],
]
)
# Create the frustum mesh
frustum_mesh = trimesh.Trimesh(vertices=vertices, faces=faces)
frustum_mesh.visual.vertex_colors = np.asarray(color, dtype=np.float32)
return trimesh.util.concatenate([camera_mesh, frustum_mesh])
[docs]def create_tets_mesh(n_tets=1, halfsize=1.0, quats=None, randomize_halfsize=True):
"""
Create artistic tet-based mesh for rendering particles as tets.
"""
# create tet-based particles given positions
vert_per_tet = 12
face_per_tet = 20
if quats is None:
quats = np.tile(gu.random_quaternion(n_tets), [1, vert_per_tet]).reshape(-1, 4)
if randomize_halfsize:
halfsize = (
np.tile(np.random.uniform(0.3, 1.9, size=(n_tets, 1)), [1, vert_per_tet * 3]).reshape(-1, 3) * halfsize
)
halfsize = (
np.tile(np.random.uniform(0.3, 1.9, size=(n_tets * 4, 1)), [1, vert_per_tet // 4 * 3]).reshape(-1, 3)
* halfsize
)
# halfsize = np.random.uniform(0.2, 1.9, size=(n_tets * vert_per_tet, 3)) * halfsize
vertices = (
np.tile(
np.array(
[
[0.91835, 0.836701, 0.91835],
[0.91835, 0.91835, 0.836701],
[0.836701, 0.91835, 0.91835],
[-0.836701, 0.91835, -0.91835],
[-0.91835, 0.836701, -0.91835],
[-0.91835, 0.91835, -0.836701],
[-0.836701, -0.91835, 0.91835],
[-0.91835, -0.836701, 0.91835],
[-0.91835, -0.91835, 0.836701],
[0.91835, -0.836701, -0.91835],
[0.91835, -0.91835, -0.836701],
[0.836701, -0.91835, -0.91835],
]
),
[n_tets, 1],
)
* halfsize
)
vertices = gu.transform_by_quat(vertices, quats)
faces = np.tile(
np.array(
[
[0, 6, 10],
[11, 8, 4],
[2, 5, 7],
[0, 1, 2],
[3, 4, 5],
[6, 7, 8],
[9, 10, 11],
[0, 10, 9],
[9, 1, 0],
[1, 3, 5],
[5, 2, 1],
[2, 7, 6],
[6, 0, 2],
[4, 8, 7],
[7, 5, 4],
[8, 11, 10],
[10, 6, 8],
[3, 9, 11],
[11, 4, 3],
[1, 9, 3],
]
),
[n_tets, 1],
)
faces_offset = np.tile(np.arange(0, n_tets).reshape(-1, 1) * vert_per_tet, [1, face_per_tet * 3]).reshape(
n_tets * face_per_tet, 3
)
faces += faces_offset
mesh = trimesh.Trimesh(vertices=vertices, faces=faces)
return mesh
def transform_tets_mesh_verts(vertices, positions, zs=None):
vert_per_tet = 12
assert len(vertices) == len(positions) * vert_per_tet
vertices = vertices.reshape(-1, vert_per_tet, 3)
if zs is not None:
assert len(zs) == len(positions)
vertices = gu.transform_by_R(vertices, gu.z_up_to_R(zs))
return (vertices + positions[:, np.newaxis]).reshape((-1, 3))
@lru_cache(maxsize=32)
def _create_unit_sphere_impl(subdivisions):
mesh = trimesh.creation.icosphere(radius=1.0, subdivisions=subdivisions)
vertices, faces = mesh.vertices.copy(), mesh.faces.copy()
attrs = {"vertex_normals": mesh.vertex_normals.copy(), "face_normals": mesh.face_normals.copy()}
for data in (vertices, faces, *attrs.values()):
data.flags.writeable = False
return vertices, faces, attrs
[docs]def create_sphere(radius, subdivisions=3, color=(1.0, 1.0, 1.0, 1.0)):
vertices, faces, attrs = _create_unit_sphere_impl(subdivisions=subdivisions)
vertices = vertices * radius
visual = trimesh.visual.ColorVisuals()
visual._data["vertex_colors"] = np.tile(color_f32_to_u8(color), (len(vertices), 1))
mesh = trimesh.Trimesh(vertices=vertices, faces=faces, visual=visual, process=False)
mesh._cache.id_set()
mesh._cache.cache.update(attrs)
return mesh
@lru_cache(maxsize=32)
def _create_unit_cylinder_impl(sections):
mesh = trimesh.creation.cylinder(radius=1.0, height=1.0, sections=sections)
vertices, faces = mesh.vertices.copy(), mesh.faces.copy()
attrs = {"vertex_normals": mesh.vertex_normals.copy(), "face_normals": mesh.face_normals.copy()}
for data in (vertices, faces, *attrs.values()):
data.flags.writeable = False
return vertices, faces, attrs
[docs]def create_cylinder(radius, height, sections=None, color=(1.0, 1.0, 1.0, 1.0)):
vertices, faces, attrs = _create_unit_cylinder_impl(sections=sections)
vertices = vertices * (radius, radius, height)
visual = trimesh.visual.ColorVisuals()
visual._data["vertex_colors"] = np.tile(color_f32_to_u8(color), (len(vertices), 1))
mesh = trimesh.Trimesh(vertices=vertices, faces=faces, visual=visual, process=False)
mesh._cache.id_set()
mesh._cache.cache.update(attrs)
return mesh
@lru_cache(maxsize=32)
def _create_unit_cone_impl(sections):
mesh = trimesh.creation.cone(radius=1.0, height=1.0, sections=sections)
vertices, faces = mesh.vertices.copy(), mesh.faces.copy()
attrs = {"vertex_normals": mesh.vertex_normals.copy(), "face_normals": mesh.face_normals.copy()}
for data in (vertices, faces, *attrs.values()):
data.flags.writeable = False
return vertices, faces, attrs
[docs]def create_cone(radius, height, sections=None, color=(1.0, 1.0, 1.0, 1.0)):
vertices, faces, attrs = _create_unit_cone_impl(sections=sections)
vertices = vertices * (radius, radius, height)
for name, normals in attrs.items():
normals = normals / (radius, radius, height)
normals /= np.linalg.norm(normals, axis=-1, keepdims=True)
attrs[name] = normals
visual = trimesh.visual.ColorVisuals()
visual._data["vertex_colors"] = np.tile(color_f32_to_u8(color), (len(vertices), 1))
mesh = trimesh.Trimesh(vertices=vertices, faces=faces, visual=visual, process=False)
mesh._cache.id_set()
mesh._cache.cache.update(attrs)
return mesh
[docs]def create_arrow(
length=1.0,
radius=0.02,
l_ratio=0.25,
r_ratio=1.5,
body_color=(1.0, 1.0, 1.0, 1.0),
head_color=(1.0, 1.0, 1.0, 1.0),
sections=12,
):
r_head = radius * r_ratio
r_body = radius
l_head = length * l_ratio
l_body = length - l_head
head = create_cone(r_head, l_head, sections=sections, color=head_color)
body = create_cylinder(r_body, l_body, sections=sections, color=body_color)
face_normals = np.vstack((body._cache["face_normals"], head._cache["face_normals"]))
face_normals.flags.writeable = False
head._data["vertices"] += np.array([0.0, 0.0, l_body])
body._data["vertices"] += np.array([0.0, 0.0, l_body / 2])
vertices = np.vstack((body.vertices, head.vertices))
faces = np.vstack((body.faces, head.faces + len(body.vertices)))
visual = trimesh.visual.ColorVisuals()
visual._data["vertex_colors"] = np.vstack((body.visual.vertex_colors, head.visual.vertex_colors))
mesh = trimesh.Trimesh(vertices=vertices, faces=faces, visual=visual, process=False)
mesh._cache.id_set()
mesh._cache.cache["face_normals"] = face_normals
return mesh
[docs]def create_line(start, end, radius=0.002, color=(1.0, 1.0, 1.0, 1.0), sections=12):
vec = end - start
length = np.linalg.norm(vec)
mesh = create_cylinder(radius, length, sections, color) # along z-axis
mesh._data["vertices"][:, -1] += length / 2.0
mesh.vertices = gu.transform_by_trans_R(mesh._data["vertices"], start, gu.z_up_to_R(vec))
return mesh
@lru_cache(maxsize=1)
def _create_unit_box_impl():
mesh = trimesh.creation.box(extents=[1.0, 1.0, 1.0])
vertices, faces = mesh.vertices.copy(), mesh.faces.copy()
attrs = {"vertex_normals": mesh.vertex_normals.copy(), "face_normals": mesh.face_normals.copy()}
for data in (vertices, faces, *attrs.values()):
data.flags.writeable = False
return vertices, faces, attrs
[docs]def create_box(extents=None, color=(1.0, 1.0, 1.0, 1.0), bounds=None, wireframe=False, wireframe_radius=0.002):
if bounds is not None:
bounds = np.asarray(bounds)
extents = bounds[1] - bounds[0]
pos = bounds.mean(axis=0)
elif extents is not None:
extents = np.asarray(extents)
pos = np.zeros(3)
else:
gs.raise_exception("Neither `extents` nor `bounds` is provided.")
if wireframe:
box_vertices = np.asarray(
[
[-0.5, -0.5, -0.5],
[0.5, -0.5, -0.5],
[0.5, 0.5, -0.5],
[-0.5, 0.5, -0.5],
[-0.5, -0.5, 0.5],
[0.5, -0.5, 0.5],
[0.5, 0.5, 0.5],
[-0.5, 0.5, 0.5],
]
)
box_vertices = box_vertices * extents + pos
box_edges = [(0, 1), (1, 2), (2, 3), (3, 0), (4, 5), (5, 6), (6, 7), (7, 4), (0, 4), (1, 5), (2, 6), (3, 7)]
n_verts = 0
vertices, faces, attrs = [], [], {}
for v_start, v_end in box_edges:
p_start, p_end = box_vertices[v_start], box_vertices[v_end]
vec = p_end - p_start
length = np.linalg.norm(vec)
line_vertices, line_faces, line_attrs = _create_unit_cylinder_impl(sections=12)
line_vertices = line_vertices * (wireframe_radius, wireframe_radius, length)
line_vertices[:, -1] += length / 2.0
line_vertices = gu.transform_by_trans_R(line_vertices, p_start, gu.z_up_to_R(vec))
vertices.append(line_vertices)
faces.append(line_faces + n_verts)
for name, value in line_attrs.items():
attrs.setdefault(name, []).append(value)
n_verts += len(line_vertices)
for vertex in box_vertices:
sphere_vertices, sphere_faces, sphere_attrs = _create_unit_sphere_impl(subdivisions=3)
vertices.append(sphere_vertices * wireframe_radius + vertex)
faces.append(sphere_faces + n_verts)
for name, value in sphere_attrs.items():
attrs.setdefault(name, []).append(value)
n_verts += len(sphere_vertices)
vertices = np.concatenate(vertices)
faces = np.concatenate(faces)
for name, values in attrs.items():
attrs[name] = np.concatenate(values)
else:
vertices, faces, attrs = _create_unit_box_impl()
vertices = vertices * extents + pos
visual = trimesh.visual.ColorVisuals()
visual._data["vertex_colors"] = np.tile(color_f32_to_u8(color), (len(vertices), 1))
mesh = trimesh.Trimesh(vertices=vertices, faces=faces, visual=visual, process=False)
mesh._cache.id_set()
mesh._cache.cache.update(attrs)
return mesh
[docs]def create_plane(
normal=(0.0, 0.0, 1.0),
plane_size=(1e3, 1e3),
tile_size=(1, 1),
color_or_texture=DEFAULT_PLANE_TEXTURE_PATH,
double_sided=False,
):
if isinstance(color_or_texture, str):
color, texture_path = None, color_or_texture
else:
color, texture_path = color_or_texture, None
thickness = 1e-2 # for safety
mesh = trimesh.creation.box(extents=[plane_size[0], plane_size[1], thickness])
mesh.vertices[:, 2] -= thickness / 2
mesh.vertices = gu.transform_by_R(mesh.vertices, gu.z_up_to_R(np.asarray(normal, dtype=np.float32)))
# Grid of vertices spanning the plane, split in cells no larger than PLANE_MAX_CELL_SIZE (see its definition)
n_cells_x, n_cells_y = (math.ceil(size / PLANE_MAX_CELL_SIZE) for size in plane_size)
verts_x = np.linspace(-0.5 * plane_size[0], 0.5 * plane_size[0], n_cells_x + 1)
verts_y = np.linspace(-0.5 * plane_size[1], 0.5 * plane_size[1], n_cells_y + 1)
verts = np.zeros((n_cells_x + 1, n_cells_y + 1, 3), dtype=np.float32)
verts[..., 0], verts[..., 1] = np.meshgrid(verts_x, verts_y, indexing="ij")
verts = verts.reshape((-1, 3))
verts_idx = np.arange(len(verts), dtype=np.int32).reshape((n_cells_x + 1, n_cells_y + 1))
v00, v10, v11, v01 = verts_idx[:-1, :-1], verts_idx[1:, :-1], verts_idx[1:, 1:], verts_idx[:-1, 1:]
faces = np.stack((np.stack((v00, v10, v11), axis=-1), np.stack((v00, v11, v01), axis=-1)), axis=-2).reshape((-1, 3))
if double_sided:
# Add reversed faces for back-facing visibility
faces = np.vstack([faces, faces[:, ::-1]])
vmesh = trimesh.Trimesh(verts, faces, process=False)
# Align visual surface with the collision top face, both at z=0 in link-local frame, so morph.pos specifies
# the actual rendered surface position.
vmesh.vertices = gu.transform_by_R(vmesh.vertices, gu.z_up_to_R(np.asarray(normal, dtype=np.float32)))
if texture_path is not None:
# One texture tile per 'tile_size', counted from the corner of the plane. The back faces share the vertices.
uv_coords = (verts[:, :2] - verts[0, :2]) / tile_size
vmesh.visual = trimesh.visual.TextureVisuals(
uv=uv_coords,
material=trimesh.visual.material.SimpleMaterial(
image=Image.open(os.path.join(get_assets_dir(), texture_path)),
),
)
else:
vmesh.visual = trimesh.visual.ColorVisuals(
vertex_colors=np.tile(np.asarray(color, dtype=np.float32), (len(vmesh.vertices), 1))
)
return vmesh, mesh
def generate_tetgen_config_from_morph(morph):
if not isinstance(morph, gs.options.morphs.TetGenMixin):
raise TypeError(
f"Expected an instance of a class that inherits from TetGenMixin, but got an instance of {type(morph).name}."
)
return dict(
order=morph.order,
mindihedral=morph.mindihedral,
minratio=morph.minratio,
nobisect=morph.nobisect,
quality=morph.quality,
maxvolume=morph.maxvolume,
verbose=morph.verbose,
)
def make_tetgen_switches(cfg):
"""Build a TetGen switches string from a config dict."""
flags = ["p"]
if cfg.get("quality", True):
r = cfg.get("minratio", 1.1)
di = cfg.get("mindihedral", 10)
flags.append(f"q{r}/{di}")
a = cfg.get("maxvolume", -1.0)
if a > 0:
flags.append(f"a{a}")
o = cfg.get("order", 1)
if o != 1:
flags.append(f"o{o}")
if cfg.get("nobisect", False):
flags.append("Y")
v = cfg.get("verbose", 0)
if v > 0:
flags.append("V" * v)
return "".join(flags)
[docs]def merge_submeshes(verts_list, faces_list):
"""Concatenate sub-meshes into a single mesh, welding vertices that share a position.
Welding covers both duplicated vertices within a sub-mesh (e.g. texture-seam splits) and vertices shared across
sub-meshes, so that downstream consumers (e.g. tetrahedralization) see one connected surface.
Parameters
----------
verts_list : list of np.ndarray
Per-sub-mesh vertex arrays with shape (n_verts, 3).
faces_list : list of np.ndarray
Per-sub-mesh face arrays with shape (n_faces, 3), indexing into the matching vertex array.
Returns
-------
verts : np.ndarray
Welded vertices, ordered by first occurrence in the concatenated input.
faces : np.ndarray
Concatenated faces, re-indexed against the welded vertices.
verts_maps : list of np.ndarray
For each sub-mesh, the map from its local vertex indices to indices into the welded vertices.
"""
submeshes_offset = (0, *np.cumsum([len(verts) for verts in verts_list]))
stacked_verts = np.concatenate(verts_list, axis=0)
stacked_faces = np.concatenate([faces + offset for faces, offset in zip(faces_list, submeshes_offset)], axis=0)
# Weld by quantized position (see VERT_WELD_QUANTIZE_FACTOR). int64 is required: quantized metre-scale
# coordinates overflow int32.
quantized = np.round(stacked_verts / VERT_WELD_QUANTIZE_FACTOR).astype(np.int64)
_, unique_idx, remap = np.unique(quantized, axis=0, return_index=True, return_inverse=True)
# np.unique orders groups by quantized key; rank the first-occurrence indices to restore input order instead.
rank = np.empty(len(unique_idx), dtype=gs.np_int)
rank[np.argsort(unique_idx)] = np.arange(len(unique_idx))
remap = rank[remap]
verts = stacked_verts[np.sort(unique_idx)]
faces = remap[stacked_faces]
verts_maps = [remap[start:end] for start, end in zip(submeshes_offset, submeshes_offset[1:])]
return verts, faces, verts_maps
[docs]def tetrahedralize_mesh(mesh, tet_cfg):
tet = tetgen.TetGen(mesh.vertices.astype(np.float64, copy=False), mesh.faces.astype(np.int32, copy=False))
# Build and apply the switches string directly, since
# the Python wrapper sometimes ignores certain kwargs
# (e.g. maxvolume). See: https://github.com/pyvista/tetgen/issues/24
verts, elems, *_ = tet.tetrahedralize(switches=make_tetgen_switches(tet_cfg))
return verts, elems
def visualize_tet(tet, mesh, show_surface=True, plot_cell_qual=False):
grid = tet.grid
if show_surface:
grid.plot(show_edges=True)
else:
# get cell centroids
cells = grid.cells.reshape(-1, 5)[:, 1:]
cell_center = grid.points[cells].mean(axis=1)
# extract cells below the 0 xy plane
cell_ind = (cell_center[:, 2] < 0.0).nonzero(as_tuple=False)
subgrid = grid.extract_cells(cell_ind)
# advanced plotting
if plot_cell_qual:
cell_qual = subgrid.compute_cell_quality()["CellQuality"]
subgrid.plot(
scalars=cell_qual, stitle="Quality", cmap="bwr", clim=[0, 1], flip_scalars=True, show_edges=True
)
else:
# Delaying import of 'pyvista' because it is an optional dependency
import pyvista as pv
faces = np.concatenate([np.full((mesh.faces.shape[0], 1), mesh.faces.shape[1]), mesh.faces], axis=1)
pv_data = pv.PolyData(mesh.vertices, faces)
plotter = pv.Plotter()
plotter.add_mesh(subgrid, "lightgrey", lighting=True, show_edges=True)
plotter.add_mesh(pv_data, "r", "wireframe")
plotter.add_legend([[" Input Mesh ", "r"], [" Tessellated Mesh ", "black"]])
plotter.show()
def check_exr_compression(exr_path):
exr_file = OpenEXR.InputFile(exr_path)
exr_header = exr_file.header()
if exr_header["compression"].v > Imath.Compression.PIZ_COMPRESSION:
new_exr_path = get_exr_path(exr_path)
if os.path.exists(new_exr_path):
gs.logger.info(f"Assets of fixed compression detected and used: {new_exr_path}.")
else:
gs.logger.warning(
f"EXR image {exr_path}'s compression type {exr_header['compression']} is not supported. "
f"Converting to compression type ZIP_COMPRESSION and saving to {new_exr_path}."
)
channel_data = {channel: exr_file.channel(channel) for channel in exr_header["channels"]}
exr_header["compression"] = Imath.Compression(Imath.Compression.ZIP_COMPRESSION)
os.makedirs(os.path.dirname(new_exr_path), exist_ok=True)
new_exr_file = OpenEXR.OutputFile(new_exr_path, exr_header)
new_exr_file.writePixels(channel_data)
new_exr_file.close()
exr_path = new_exr_path
exr_file.close()
return exr_path
