Skip to content

Commit d4e8556

Browse files
committed
Merge branch 'main' into 440-add-inertia-verification-correction
2 parents c979014 + 67a17ff commit d4e8556

11 files changed

Lines changed: 555 additions & 122 deletions

File tree

docs/generated/core/newton.ModelBuilder.rst

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -49,6 +49,7 @@
4949
~ModelBuilder.add_tetrahedron
5050
~ModelBuilder.add_triangle
5151
~ModelBuilder.add_triangles
52+
~ModelBuilder.approximate_meshes
5253
~ModelBuilder.collapse_fixed_joints
5354
~ModelBuilder.color
5455
~ModelBuilder.finalize

newton/examples/example_g1.py

Lines changed: 19 additions & 42 deletions
Original file line numberDiff line numberDiff line change
@@ -14,19 +14,21 @@
1414
# limitations under the License.
1515

1616
###########################################################################
17-
# Example Sim G1
17+
# Example G1
1818
#
1919
# Shows how to set up a simulation of a rigid-body humanoid articulation
20-
# from a xml using the newton.ModelBuilder().
20+
# from a XML using the newton.ModelBuilder().
2121
# Note this example does not include a trained policy.
22+
#
23+
# Users can pick bodies by right-clicking and dragging with the mouse.
24+
#
2225
###########################################################################
2326

2427
import warp as wp
2528

2629
wp.config.enable_backward = False
2730

2831
import newton
29-
import newton.examples
3032
import newton.utils
3133

3234

@@ -45,61 +47,27 @@ def __init__(self, stage_path="example_g1.usd", num_envs=8, use_cuda_graph=True)
4547
up_axis="Z",
4648
enable_self_collisions=False,
4749
)
48-
simplified_meshes = {}
49-
try:
50-
import tqdm # noqa: PLC0415
51-
52-
meshes = tqdm.tqdm(articulation_builder.shape_geo_src, desc="Simplifying meshes")
53-
except ImportError:
54-
meshes = articulation_builder.shape_geo_src
55-
for i, m in enumerate(meshes):
56-
if m is None:
57-
continue
58-
hash_m = hash(m)
59-
if hash_m in simplified_meshes:
60-
articulation_builder.shape_geo_src[i] = simplified_meshes[hash_m]
61-
else:
62-
simplified = newton.geometry.utils.remesh_mesh(
63-
m, visualize=False, method="convex_hull", recompute_inertia=False
64-
)
65-
# simplified = newton.geometry.utils.remesh_mesh(
66-
# simplified, visualize=False, target_reduction=None, target_count=32, recompute_inertia=False
67-
# )
68-
# simplified = newton.geometry.utils.remesh_mesh(
69-
# simplified, visualize=False, method="convex_hull", recompute_inertia=False
70-
# )
71-
# simplified = newton.geometry.utils.remesh_mesh(
72-
# simplified, visualize=False, method="convex_hull", alpha=0.01, recompute_inertia=False
73-
# )
74-
# simplified = newton.geometry.utils.remesh_mesh(
75-
# m, visualize=False, method="ftetwild", edge_length_fac=0.5, optimize=True, recompute_inertia=False
76-
# )
77-
articulation_builder.shape_geo_src[i] = simplified
78-
simplified_meshes[hash_m] = simplified
50+
articulation_builder.approximate_meshes("bounding_box")
7951

8052
spacing = 3.0
8153
sqn = int(wp.ceil(wp.sqrt(float(self.num_envs))))
8254

8355
builder = newton.ModelBuilder()
8456
for i in range(self.num_envs):
85-
pos = wp.vec3((i % sqn) * spacing, (i // sqn) * spacing, 2)
57+
pos = wp.vec3((i % sqn) * spacing, (i // sqn) * spacing, 0.0)
8658
builder.add_builder(articulation_builder, xform=wp.transform(pos, wp.quat_identity()))
8759
builder.add_ground_plane()
8860

8961
self.sim_time = 0.0
9062
fps = 600
9163
self.frame_dt = 1.0 / fps
9264

93-
self.sim_substeps = 5
94-
self.sim_dt = self.frame_dt / self.sim_substeps
95-
9665
# finalize model
9766
self.model = builder.finalize()
9867

9968
self.control = self.model.control()
100-
# self.solver = newton.solvers.FeatherstoneSolver(self.model)
101-
# self.solver = newton.solvers.SemiImplicitSolver(self.model, joint_attach_kd=100, joint_attach_ke= 1000)
10269
if self.use_mujoco:
70+
self.sim_substeps = 4
10371
self.solver = newton.solvers.MuJoCoSolver(
10472
self.model,
10573
use_mujoco=False,
@@ -111,7 +79,15 @@ def __init__(self, stage_path="example_g1.usd", num_envs=8, use_cuda_graph=True)
11179
ncon_per_env=150,
11280
)
11381
else:
114-
self.solver = newton.solvers.XPBDSolver(self.model, iterations=20)
82+
self.sim_substeps = 10
83+
self.solver = newton.solvers.XPBDSolver(
84+
self.model,
85+
iterations=20,
86+
angular_damping=0.01,
87+
joint_angular_compliance=1e-3,
88+
)
89+
90+
self.sim_dt = self.frame_dt / self.sim_substeps
11591

11692
self.renderer = None
11793

@@ -120,7 +96,6 @@ def __init__(self, stage_path="example_g1.usd", num_envs=8, use_cuda_graph=True)
12096
path=stage_path,
12197
model=self.model,
12298
scaling=1.0,
123-
up_axis=str(newton.Axis.Z),
12499
screen_width=1280,
125100
screen_height=720,
126101
camera_pos=(0, 1, 4),
@@ -146,6 +121,8 @@ def simulate(self):
146121
self.contacts = self.model.collide(self.state_0)
147122
for _ in range(self.sim_substeps):
148123
self.state_0.clear_forces()
124+
if self.renderer and hasattr(self.renderer, "apply_picking_force"):
125+
self.renderer.apply_picking_force(self.state_0)
149126
self.solver.step(self.state_0, self.state_1, self.control, self.contacts, self.sim_dt)
150127
self.state_0, self.state_1 = self.state_1, self.state_0
151128

newton/geometry/kernels.py

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -766,7 +766,8 @@ def create_soft_contacts(
766766
# region Rigid body collision detection
767767

768768

769-
@wp.kernel(enable_backward=False)
769+
# NOTE: Kernel is in a unique module to speed up cold-start ModelBuilder.finalize() time
770+
@wp.kernel(enable_backward=False, module="unique")
770771
def count_contact_points(
771772
contact_pairs: wp.array(dtype=wp.vec2i),
772773
shape_type: wp.array(dtype=int),

newton/geometry/types.py

Lines changed: 20 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -133,6 +133,26 @@ def __init__(
133133
self.mass = 1.0
134134
self.com = wp.vec3()
135135

136+
def copy(
137+
self,
138+
vertices: Sequence[Vec3] | nparray | None = None,
139+
indices: Sequence[int] | nparray | None = None,
140+
recompute_inertia: bool = False,
141+
):
142+
if vertices is None:
143+
vertices = self.vertices
144+
if indices is None:
145+
indices = self.indices
146+
m = Mesh(
147+
vertices, indices, compute_inertia=recompute_inertia, is_solid=self.is_solid, maxhullvert=self.maxhullvert
148+
)
149+
if not recompute_inertia:
150+
m.I = self.I
151+
m.mass = self.mass
152+
m.com = self.com
153+
m.has_inertia = self.has_inertia
154+
return m
155+
136156
@property
137157
def vertices(self):
138158
return self._vertices

newton/geometry/utils.py

Lines changed: 171 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -59,6 +59,83 @@ def compute_shape_radius(geo_type: int, scale: Vec3, src: Mesh | SDF | None) ->
5959
return 10.0
6060

6161

62+
def compute_aabb(vertices: nparray) -> tuple[Vec3, Vec3]:
63+
"""Compute the axis-aligned bounding box of a set of vertices."""
64+
min_coords = np.min(vertices, axis=0)
65+
max_coords = np.max(vertices, axis=0)
66+
return min_coords, max_coords
67+
68+
69+
def compute_obb(vertices: nparray) -> tuple[wp.transform, wp.vec3]:
70+
"""Compute the oriented bounding box of a set of vertices.
71+
72+
Args:
73+
vertices: A numpy array of shape (N, 3) containing the vertex positions.
74+
75+
Returns:
76+
A tuple containing:
77+
- transform: The transform of the oriented bounding box
78+
- extents: The half-extents of the box along its principal axes
79+
"""
80+
if len(vertices) == 0:
81+
return wp.transform_identity(), wp.vec3(0.0, 0.0, 0.0)
82+
if len(vertices) == 1:
83+
return wp.transform(wp.vec3(vertices[0]), wp.quat_identity()), wp.vec3(0.0, 0.0, 0.0)
84+
85+
# Center the vertices
86+
center = np.mean(vertices, axis=0)
87+
centered_vertices = vertices - center
88+
89+
# Compute covariance matrix with handling for degenerate cases
90+
if len(vertices) < 3:
91+
# For 2 points, create a line-aligned OBB
92+
direction = centered_vertices[1] if len(vertices) > 1 else np.array([1, 0, 0])
93+
direction = direction / np.linalg.norm(direction) if np.linalg.norm(direction) > 1e-6 else np.array([1, 0, 0])
94+
# Create orthogonal basis
95+
if abs(direction[0]) < 0.9:
96+
perpendicular = np.cross(direction, [1, 0, 0])
97+
else:
98+
perpendicular = np.cross(direction, [0, 1, 0])
99+
perpendicular = perpendicular / np.linalg.norm(perpendicular)
100+
third = np.cross(direction, perpendicular)
101+
eigenvectors = np.column_stack([direction, perpendicular, third])
102+
else:
103+
cov_matrix = np.cov(centered_vertices.T)
104+
eigenvalues, eigenvectors = np.linalg.eigh(cov_matrix)
105+
# Sort by eigenvalues in descending order
106+
sorted_indices = np.argsort(eigenvalues)[::-1]
107+
eigenvalues = eigenvalues[sorted_indices]
108+
eigenvectors = eigenvectors[:, sorted_indices]
109+
110+
# Ensure right-handed coordinate system
111+
if np.linalg.det(eigenvectors) < 0:
112+
eigenvectors[:, 2] *= -1
113+
114+
# Project vertices onto principal axes
115+
projected = centered_vertices @ eigenvectors
116+
117+
# Compute extents
118+
min_coords = np.min(projected, axis=0)
119+
max_coords = np.max(projected, axis=0)
120+
extents = (max_coords - min_coords) / 2.0
121+
122+
# Calculate the center in the projected coordinate system
123+
center_offset = (max_coords + min_coords) / 2.0
124+
# Transform the center offset back to the original coordinate system
125+
center = center + center_offset @ eigenvectors.T
126+
127+
# Convert rotation matrix to quaternion
128+
# The rotation matrix should transform from the original coordinate system to the principal axes
129+
# eigenvectors is the rotation matrix from original to principal axes
130+
rotation_matrix = eigenvectors
131+
132+
# Convert to quaternion using Warp's quat_from_matrix function
133+
# First convert numpy array to Warp matrix
134+
orientation = wp.quat_from_matrix(wp.mat33(rotation_matrix))
135+
136+
return wp.transform(wp.vec3(center), orientation), wp.vec3(extents)
137+
138+
62139
def load_mesh(filename: str, method: str | None = None):
63140
"""
64141
Loads a 3D triangular surface mesh from a file.
@@ -339,6 +416,13 @@ def remesh_convex_hull(vertices, maxhullvert: int = 0):
339416
if np.dot(normal, a - centre) < 0:
340417
faces[i] = tri[[0, 2, 1]]
341418

419+
# trim vertices to only those that are used in the faces
420+
unique_verts = np.unique(faces.flatten())
421+
verts = verts[unique_verts]
422+
# update face indices to use the new vertex indices
423+
mapping = {v: i for i, v in enumerate(unique_verts)}
424+
faces = np.array([mapping[v] for v in faces.flatten()], dtype=np.int32).reshape(faces.shape)
425+
342426
return verts, faces
343427

344428

@@ -384,23 +468,105 @@ def remesh(
384468
return new_vertices, new_faces
385469

386470

387-
def remesh_mesh(mesh: Mesh, method: RemeshingMethod = "quadratic", recompute_inertia=False, **remeshing_kwargs) -> Mesh:
471+
def remesh_mesh(
472+
mesh: Mesh,
473+
method: RemeshingMethod = "quadratic",
474+
recompute_inertia: bool = False,
475+
inplace: bool = False,
476+
**remeshing_kwargs,
477+
) -> Mesh:
388478
"""Remesh a mesh using the specified method.
389479
Args:
390480
mesh: The mesh to remesh.
391481
method: The remeshing method to use. One of "ftetwild", "quadratic", "convex_hull", or "alphashape".
392482
recompute_inertia: Whether to recompute the inertia of the mesh.
483+
inplace: Whether to modify the mesh in place or return a new mesh.
393484
**remeshing_kwargs: Additional keyword arguments passed to the remeshing function.
394485
Returns:
395486
The remeshed mesh.
396487
"""
397488
if method == "convex_hull":
398489
remeshing_kwargs["maxhullvert"] = mesh.maxhullvert
399-
mesh.vertices, mesh.indices = remesh(mesh.vertices, mesh.indices.reshape(-1, 3), method=method, **remeshing_kwargs)
400-
mesh.indices = mesh.indices.flatten()
401-
if recompute_inertia:
402-
mesh.mass, mesh.com, mesh.I, _ = compute_mesh_inertia(1.0, mesh.vertices, mesh.indices, is_solid=mesh.is_solid)
490+
vertices, indices = remesh(mesh.vertices, mesh.indices.reshape(-1, 3), method=method, **remeshing_kwargs)
491+
if inplace:
492+
mesh.vertices = vertices
493+
mesh.indices = indices.flatten()
494+
if recompute_inertia:
495+
mesh.mass, mesh.com, mesh.I, _ = compute_mesh_inertia(1.0, vertices, indices, is_solid=mesh.is_solid)
496+
else:
497+
return mesh.copy(vertices=vertices, indices=indices, recompute_inertia=recompute_inertia)
403498
return mesh
404499

405500

501+
def create_box_mesh(half_extents: Vec3) -> tuple[nparray, nparray]:
502+
x_extent, y_extent, z_extent = half_extents
503+
vertices = np.array(
504+
[
505+
[-x_extent, -y_extent, -z_extent],
506+
[x_extent, -y_extent, -z_extent],
507+
[x_extent, y_extent, -z_extent],
508+
[-x_extent, y_extent, -z_extent],
509+
[-x_extent, -y_extent, z_extent],
510+
[x_extent, -y_extent, z_extent],
511+
[x_extent, y_extent, z_extent],
512+
[-x_extent, y_extent, z_extent],
513+
],
514+
dtype=np.float32,
515+
)
516+
indices = np.array(
517+
[
518+
# Bottom face (z = -z_extent)
519+
0,
520+
2,
521+
1,
522+
0,
523+
3,
524+
2,
525+
# Top face (z = z_extent)
526+
4,
527+
5,
528+
6,
529+
4,
530+
6,
531+
7,
532+
# Front face (y = -y_extent)
533+
0,
534+
1,
535+
5,
536+
0,
537+
5,
538+
4,
539+
# Back face (y = y_extent)
540+
2,
541+
3,
542+
7,
543+
2,
544+
7,
545+
6,
546+
# Left face (x = -x_extent)
547+
0,
548+
4,
549+
7,
550+
0,
551+
7,
552+
3,
553+
# Right face (x = x_extent)
554+
1,
555+
2,
556+
6,
557+
1,
558+
6,
559+
5,
560+
],
561+
dtype=np.int32,
562+
)
563+
return vertices, indices
564+
565+
566+
def transform_points(points: nparray, transform: wp.transform, scale: Vec3 | None = None) -> nparray:
567+
if scale is not None:
568+
points = points * np.array(scale, dtype=np.float32)
569+
return points @ np.array(wp.quat_to_matrix(transform.q)).reshape(3, 3) + transform.p
570+
571+
406572
__all__ = ["compute_shape_radius", "load_mesh", "visualize_meshes"]

0 commit comments

Comments
 (0)