Coverage for trimesh/inertia.py: 89%
83 statements
« prev ^ index » next coverage.py v7.14.1, created at 2026-06-24 04:40 +0000
1"""
2inertia.py
3-------------
5Functions for dealing with inertia tensors.
7Results validated against known geometries and checked for
8internal consistency.
9"""
11import numpy as np
13from .typed import ArrayLike, NDArray, Number, float64
14from .util import multi_dot
17def cylinder_inertia(
18 mass: Number, radius: Number, height: Number, transform: ArrayLike | None = None
19) -> NDArray[float64]:
20 """
21 Return the inertia tensor of a cylinder.
23 Parameters
24 ------------
25 mass : float
26 Mass of cylinder
27 radius : float
28 Radius of cylinder
29 height : float
30 Height of cylinder
31 transform : (4, 4) float
32 Transformation of cylinder
34 Returns
35 ------------
36 inertia : (3, 3) float
37 Inertia tensor
38 """
39 h2, r2 = height**2, radius**2
40 diagonal = np.array(
41 [
42 ((mass * h2) / 12) + ((mass * r2) / 4),
43 ((mass * h2) / 12) + ((mass * r2) / 4),
44 (mass * r2) / 2,
45 ]
46 )
47 inertia = diagonal * np.eye(3)
49 if transform is not None:
50 inertia = transform_inertia(transform, inertia)
52 return inertia
55def sphere_inertia(mass: Number, radius: Number) -> NDArray[float64]:
56 """
57 Return the inertia tensor of a sphere.
59 Parameters
60 ------------
61 mass : float
62 Mass of sphere
63 radius : float
64 Radius of sphere
66 Returns
67 ------------
68 inertia : (3, 3) float
69 Inertia tensor
70 """
71 return (2.0 / 5.0) * (radius**2) * mass * np.eye(3)
74def points_inertia(
75 points: ArrayLike,
76 weights: None | ArrayLike | Number = None,
77 at_center_mass: bool = True,
78) -> NDArray[float64]:
79 """
80 Calculate an inertia tensor for an array of point masses
81 at the center of mass.
83 Parameters
84 ----------
85 points : (n, 3)
86 Points in space.
87 weights : (n,) or number
88 Per-point weight to use.
89 at_center_mass
90 Calculate at the center of mass of the points, or if False
91 at the original origin.
93 Returns
94 -----------
95 tensor : (3, 3)
96 Inertia tensor for point masses.
97 """
98 if weights is None:
99 # by default make the total weight 1.0 to match
100 # the default mass in other functions, and so that
101 # if a user didn't specify anything it doesn't blow
102 # up the scale depending on the number of points
103 weights = np.full(len(points), 1.0 / float(len(points)), dtype=np.float64)
104 elif isinstance(weights, (float, np.integer, int)):
105 # "is it a number" check
106 weights = np.full(len(points), float(weights), dtype=np.float64)
107 else:
108 weights = np.array(weights)
109 if len(weights) != len(points):
110 raise ValueError(
111 f"Weights must correspond to points! {len(weights)} != {len(points)}"
112 )
114 # make sure the points are an array of correct shape
115 points = np.asanyarray(points, dtype=np.float64)
116 if len(points.shape) != 2 or points.shape[1] != 3:
117 raise ValueError(f"Points must be `(n, 3)` not {points.shape}")
119 if at_center_mass:
120 # get the center of mass of the points
121 center_mass = np.average(points, weights=weights, axis=0)
122 # get the points with the origin at their center of mass
123 points_com = points - center_mass
124 else:
125 # calculate at original origin
126 points_com = points
128 # expand into shorthand for the expressions
129 x, y, z = points_com.T
130 x2, y2, z2 = (points_com**2).T
132 # calculate tensors per-point in a flattened (9, n) array
133 # from physics.stackexchange.com/questions/614094
134 tensors = np.array(
135 [y2 + z2, -x * y, -x * z, -x * y, x2 + z2, -y * z, -x * z, -y * z, x2 + y2],
136 dtype=np.float64,
137 )
139 # combine the weighted tensors and reshape
140 tensor = (tensors * weights).sum(axis=1).reshape((3, 3))
142 return tensor
145def principal_axis(inertia: ArrayLike):
146 """
147 Find the principal components and principal axis
148 of inertia from the inertia tensor.
150 Parameters
151 ------------
152 inertia : (3, 3) float
153 Inertia tensor
155 Returns
156 ------------
157 components : (3,) float
158 Principal components of inertia
159 vectors : (3, 3) float
160 Row vectors pointing along the
161 principal axes of inertia
162 """
163 inertia = np.asanyarray(inertia, dtype=np.float64)
164 if inertia.shape != (3, 3):
165 raise ValueError("inertia tensor must be (3, 3)!")
167 # you could any of the following to calculate this:
168 # np.linalg.svd, np.linalg.eig, np.linalg.eigh
169 # moment of inertia is square symmetric matrix
170 # eigh has the best precision in tests
171 components, vectors = np.linalg.eigh(inertia)
173 # eigh returns them as column vectors, change them to row vectors
174 vectors = vectors.T
176 return components, vectors
179def transform_inertia(
180 transform: ArrayLike,
181 inertia_tensor: ArrayLike,
182 parallel_axis: bool = False,
183 mass: Number | None = None,
184):
185 """
186 Transform an inertia tensor to a new frame.
188 Note that in trimesh `mesh.moment_inertia` is *axis aligned*
189 and at `mesh.center_mass`.
191 So to transform to a new frame and get the moment of inertia at
192 the center of mass the translation should be ignored and only
193 rotation applied.
195 If parallel axis is enabled it will compute the inertia
196 about a new location.
198 More details in the MIT OpenCourseWare PDF:
199 ` MIT16_07F09_Lec26.pdf`
202 Parameters
203 ------------
204 transform : (3, 3) or (4, 4) float
205 Transformation matrix
206 inertia_tensor : (3, 3) float
207 Inertia tensor.
208 parallel_axis : bool
209 Apply the parallel axis theorum or not.
210 If the passed inertia tensor is at the center of mass
211 and you want the new post-transform tensor also at the
212 center of mass you DON'T want this enabled as you *only*
213 want to apply the rotation. Use this to get moment of
214 inertia at an arbitrary frame that isn't the center of mass.
216 Returns
217 ------------
218 transformed : (3, 3) float
219 Inertia tensor in new frame.
220 """
221 # check inputs and extract rotation
222 transform = np.asanyarray(transform, dtype=np.float64)
223 if transform.shape == (4, 4):
224 rotation = transform[:3, :3]
225 elif transform.shape == (3, 3):
226 rotation = transform
227 else:
228 raise ValueError("transform must be (3, 3) or (4, 4)!")
230 inertia_tensor = np.asanyarray(inertia_tensor, dtype=np.float64)
231 if inertia_tensor.shape != (3, 3):
232 raise ValueError("inertia_tensor must be (3, 3)!")
234 if parallel_axis:
235 if transform.shape == (3, 3):
236 # shorthand for "translation"
237 a = np.zeros(3, dtype=np.float64)
238 else:
239 # get the translation
240 a = transform[:3, 3]
241 # First the changed origin of the new transform is taken into
242 # account. To calculate the inertia tensor
243 # the parallel axis theorem is used
244 M = np.array(
245 [
246 [a[1] ** 2 + a[2] ** 2, -a[0] * a[1], -a[0] * a[2]],
247 [-a[0] * a[1], a[0] ** 2 + a[2] ** 2, -a[1] * a[2]],
248 [-a[0] * a[2], -a[1] * a[2], a[0] ** 2 + a[1] ** 2],
249 ]
250 )
251 aligned_inertia = inertia_tensor + mass * M
253 return multi_dot([rotation.T, aligned_inertia, rotation])
255 return multi_dot([rotation, inertia_tensor, rotation.T])
258def radial_symmetry(mesh):
259 """
260 Check whether a mesh has radial symmetry.
262 Returns
263 -----------
264 symmetry : None or str
265 None No rotational symmetry
266 'radial' Symmetric around an axis
267 'spherical' Symmetric around a point
268 axis : None or (3,) float
269 Rotation axis or point
270 section : None or (3, 2) float
271 If radial symmetry provide vectors
272 to get cross section
273 """
275 # shortcuts to avoid typing and hitting cache
276 scalar = mesh.principal_inertia_components.copy()
278 # exit early if inertia components are all zero
279 if (scalar < 1e-30).any():
280 return None, None, None
282 # normalize the PCI so we can compare them
283 scalar = scalar / np.linalg.norm(scalar)
284 vector = mesh.principal_inertia_vectors
285 # the sorted order of the principal components
286 order = scalar.argsort()
288 # we are checking if a geometry has radial symmetry
289 # if 2 of the PCI are equal, it is a revolved 2D profile
290 # if 3 of the PCI (all of them) are equal it is a sphere
291 diff = np.abs(np.diff(scalar[order]))
292 # diffs that are within tol of zero
293 diff_zero = diff < 1e-4
295 if diff_zero.all():
296 # this is the case where all 3 PCI are identical
297 # this means that the geometry is symmetric about a point
298 # examples of this are a sphere, icosahedron, etc
299 axis = vector[0]
300 section = vector[1:]
302 return "spherical", axis, section
304 elif diff_zero.any():
305 # this is the case for 2/3 PCI are identical
306 # this means the geometry is symmetric about an axis
307 # probably a revolved 2D profile
309 # we know that only 1/2 of the diff values are True
310 # if the first diff is 0, it means if we take the first element
311 # in the ordered PCI we will have one of the non- revolve axis
312 # if the second diff is 0, we take the last element of
313 # the ordered PCI for the section axis
314 # if we wanted the revolve axis we would just switch [0,-1] to
315 # [-1,0]
317 # since two vectors are the same, we know the middle
318 # one is one of those two
319 section_index = order[np.array([[0, 1], [1, -1]])[diff_zero]].flatten()
320 section = vector[section_index]
322 # we know the rotation axis is the sole unique value
323 # and is either first or last of the sorted values
324 axis_index = order[np.array([-1, 0])[diff_zero]][0]
325 axis = vector[axis_index]
326 return "radial", axis, section
328 return None, None, None
331def scene_inertia(scene, transform: ArrayLike | None = None) -> NDArray[float64]:
332 """
333 Calculate the inertia of a scene about a specific frame.
335 Parameters
336 ------------
337 scene : trimesh.Scene
338 Scene with geometry.
339 transform : None or (4, 4) float
340 Homogeneous transform to compute inertia at.
342 Returns
343 ----------
344 moment : (3, 3)
345 Inertia tensor about requested frame
346 """
347 # shortcuts for tight loop
348 graph = scene.graph
349 geoms = scene.geometry
351 # get the matrix ang geometry name for
352 nodes = [graph[n] for n in graph.nodes_geometry]
353 # get the moment of inertia with the mesh moved to a location
354 moments = np.array(
355 [
356 geoms[g].moment_inertia_frame(np.dot(np.linalg.inv(mat), transform))
357 for mat, g in nodes
358 if hasattr(geoms[g], "moment_inertia_frame")
359 ],
360 dtype=np.float64,
361 )
363 return moments.sum(axis=0)