Visual Computing Library
Loading...
Searching...
No Matches
mesh_inertia.h
1/*****************************************************************************
2 * VCLib *
3 * Visual Computing Library *
4 * *
5 * Copyright(C) 2021-2025 *
6 * Visual Computing Lab *
7 * ISTI - Italian National Research Council *
8 * *
9 * All rights reserved. *
10 * *
11 * This program is free software; you can redistribute it and/or modify *
12 * it under the terms of the Mozilla Public License Version 2.0 as published *
13 * by the Mozilla Foundation; either version 2 of the License, or *
14 * (at your option) any later version. *
15 * *
16 * This program is distributed in the hope that it will be useful, *
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
19 * Mozilla Public License Version 2.0 *
20 * (https://www.mozilla.org/en-US/MPL/2.0/) for more details. *
21 ****************************************************************************/
22
23#ifndef VCL_SPACE_COMPLEX_MESH_INERTIA_H
24#define VCL_SPACE_COMPLEX_MESH_INERTIA_H
25
26#include <vclib/algorithms/core/polygon.h>
27#include <vclib/mesh/requirements.h>
28#include <vclib/space/core/point.h>
29
30#include <Eigen/Eigenvalues>
31
32namespace vcl {
33
54template<FaceMeshConcept MeshType>
55class MeshInertia
56{
57 using VertexType = MeshType::VertexType;
58 using FaceType = MeshType::FaceType;
59 using ScalarType = VertexType::CoordType::ScalarType;
60
61 enum { X = 0, Y = 1, Z = 2 };
62
63 int mA; // alpha
64 int mB; // beta
65 int mC; // gamma
66
67 /* projection integrals */
68 double mP1, mPa, mPb, mPaa, mPab, mPbb, mPaaa, mPaab, mPabb, mPbbb;
69
70 /* face integrals */
71 double mFa, mFb, mFc, mFaa, mFbb, mFcc;
72 double mFaaa, mFbbb, mFccc, mFaab, mFbbc, mFcca;
73
74 /* volume integrals */
75 double mT0 = 0, mT1[3] = {0, 0, 0}, mT2[3] = {0, 0, 0}, mTP[3] = {0, 0, 0};
76
77public:
78 MeshInertia(const MeshType& m)
79 {
80 double nx, ny, nz;
81
82 for (const FaceType& f : m.faces()) {
83 if (faceArea(f) > std::numeric_limits<float>::min()) {
84 Point3<ScalarType> fn = faceNormal(f);
85 fn.normalize();
86
87 nx = std::abs(fn[0]);
88 ny = std::abs(fn[1]);
89 nz = std::abs(fn[2]);
90 if (nx > ny && nx > nz)
91 mC = X;
92 else
93 mC = (ny > nz) ? Y : Z;
94 mA = (mC + 1) % 3;
95 mB = (mA + 1) % 3;
96
97 faceIntegrals(f, fn);
98
99 mT0 += fn[X] * ((mA == X) ? mFa : ((mB == X) ? mFb : mFc));
100
101 mT1[mA] += fn[mA] * mFaa;
102 mT1[mB] += fn[mB] * mFbb;
103 mT1[mC] += fn[mC] * mFcc;
104 mT2[mA] += fn[mA] * mFaaa;
105 mT2[mB] += fn[mB] * mFbbb;
106 mT2[mC] += fn[mC] * mFccc;
107 mTP[mA] += fn[mA] * mFaab;
108 mTP[mB] += fn[mB] * mFbbc;
109 mTP[mC] += fn[mC] * mFcca;
110 }
111 }
112
113 mT1[X] /= 2;
114 mT1[Y] /= 2;
115 mT1[Z] /= 2;
116 mT2[X] /= 3;
117 mT2[Y] /= 3;
118 mT2[Z] /= 3;
119 mTP[X] /= 2;
120 mTP[Y] /= 2;
121 mTP[Z] /= 2;
122 }
123
129 ScalarType volume() const { return mT0; }
130
131 Point3<ScalarType> centerOfMass() const
132 {
133 Point3<ScalarType> r;
134 r[X] = mT1[X] / mT0;
135 r[Y] = mT1[Y] / mT0;
136 r[Z] = mT1[Z] / mT0;
137 return r;
138 }
139
140 template<typename Matrix33>
141 Matrix33 inertiaTensor() const
142 {
143 Matrix33 J;
144 Point3d r;
145 r[X] = mT1[X] / mT0;
146 r[Y] = mT1[Y] / mT0;
147 r[Z] = mT1[Z] / mT0;
148 /* compute inertia tensor */
149 J(X, X) = (mT2[Y] + mT2[Z]);
150 J(Y, Y) = (mT2[Z] + mT2[X]);
151 J(Z, Z) = (mT2[X] + mT2[Y]);
152 J(X, Y) = J(Y, X) = -mTP[X];
153 J(Y, Z) = J(Z, Y) = -mTP[Y];
154 J(Z, X) = J(X, Z) = -mTP[Z];
155
156 J(X, X) -= mT0 * (r[Y] * r[Y] + r[Z] * r[Z]);
157 J(Y, Y) -= mT0 * (r[Z] * r[Z] + r[X] * r[X]);
158 J(Z, Z) -= mT0 * (r[X] * r[X] + r[Y] * r[Y]);
159 J(X, Y) = J(Y, X) += mT0 * r[X] * r[Y];
160 J(Y, Z) = J(Z, Y) += mT0 * r[Y] * r[Z];
161 J(Z, X) = J(X, Z) += mT0 * r[Z] * r[X];
162 return J;
163 }
164
173 template<typename Matrix33>
174 void inertiaTensorEigen(
175 Matrix33& eigenValues,
176 Point3<ScalarType>& eigenVector) const
177 {
178 Eigen::Matrix3d it = inertiaTensor<Eigen::Matrix3d>();
179 Eigen::SelfAdjointEigenSolver<Eigen::Matrix3d> eig(it);
180 Eigen::Vector3d c_val = eig.eigenvalues();
181 // eigenvector are stored as columns.
182 Eigen::Matrix3d c_vec = eig.eigenvectors();
183 c_vec.transposeInPlace();
184 for (uint i = 0; i < 3; i++)
185 for (uint j = 0; j < 3; j++)
186 eigenValues(i, j) = c_vec(i, j);
187 eigenVector = Point3<ScalarType>(c_val);
188 }
189
190private:
191 static ScalarType sqr(const ScalarType& x) { return x * x; }
192
193 static ScalarType cube(const ScalarType& x) { return x * x * x; }
194
195 void faceIntegrals(const FaceType& f, const Point3<ScalarType>& n)
196 {
197 ScalarType w;
198 double k1, k2, k3, k4;
199
200 projectionIntegrals(f);
201
202 w = -f.vertex(0)->coord().dot(n);
203 k1 = 1 / n[mC];
204 k2 = k1 * k1;
205 k3 = k2 * k1;
206 k4 = k3 * k1;
207
208 mFa = k1 * mPa;
209 mFb = k1 * mPb;
210 mFc = -k2 * (n[mA] * mPa + n[mB] * mPb + w * mP1);
211
212 mFaa = k1 * mPaa;
213 mFbb = k1 * mPbb;
214 mFcc = k3 * (sqr(n[mA]) * mPaa + 2 * n[mA] * n[mB] * mPab +
215 sqr(n[mB]) * mPbb +
216 w * (2 * (n[mA] * mPa + n[mB] * mPb) + w * mP1));
217
218 mFaaa = k1 * mPaaa;
219 mFbbb = k1 * mPbbb;
220 mFccc = -k4 * (cube(n[mA]) * mPaaa + 3 * sqr(n[mA]) * n[mB] * mPaab +
221 3 * n[mA] * sqr(n[mB]) * mPabb + cube(n[mB]) * mPbbb +
222 3 * w *
223 (sqr(n[mA]) * mPaa + 2 * n[mA] * n[mB] * mPab +
224 sqr(n[mB]) * mPbb) +
225 w * w * (3 * (n[mA] * mPa + n[mB] * mPb) + w * mP1));
226
227 mFaab = k1 * mPaab;
228 mFbbc = -k2 * (n[mA] * mPabb + n[mB] * mPbbb + w * mPbb);
229 mFcca = k3 * (sqr(n[mA]) * mPaaa + 2 * n[mA] * n[mB] * mPaab +
230 sqr(n[mB]) * mPabb +
231 w * (2 * (n[mA] * mPaa + n[mB] * mPab) + w * mPa));
232 }
233
238 void projectionIntegrals(const FaceType& f)
239 {
240 double a0, a1, da;
241 double b0, b1, db;
242 double a0_2, a0_3, a0_4, b0_2, b0_3, b0_4;
243 double a1_2, a1_3, b1_2, b1_3;
244 double C1, Ca, Caa, Caaa, Cb, Cbb, Cbbb;
245 double Cab, Kab, Caab, Kaab, Cabb, Kabb;
246
247 mP1 = mPa = mPb = mPaa = mPab = mPbb = mPaaa = mPaab = mPabb = mPbbb =
248 0.0;
249
250 for (uint i = 0; i < f.vertexNumber(); i++) {
251 a0 = f.vertex(i)->coord()[mA];
252 b0 = f.vertex(i)->coord()[mB];
253 a1 = f.vertexMod(i + 1)->coord()[mA];
254 b1 = f.vertexMod(i + 1)->coord()[mB];
255 da = a1 - a0;
256 db = b1 - b0;
257 a0_2 = a0 * a0;
258 a0_3 = a0_2 * a0;
259 a0_4 = a0_3 * a0;
260 b0_2 = b0 * b0;
261 b0_3 = b0_2 * b0;
262 b0_4 = b0_3 * b0;
263 a1_2 = a1 * a1;
264 a1_3 = a1_2 * a1;
265 b1_2 = b1 * b1;
266 b1_3 = b1_2 * b1;
267
268 C1 = a1 + a0;
269 Ca = a1 * C1 + a0_2;
270 Caa = a1 * Ca + a0_3;
271 Caaa = a1 * Caa + a0_4;
272 Cb = b1 * (b1 + b0) + b0_2;
273 Cbb = b1 * Cb + b0_3;
274 Cbbb = b1 * Cbb + b0_4;
275 Cab = 3 * a1_2 + 2 * a1 * a0 + a0_2;
276 Kab = a1_2 + 2 * a1 * a0 + 3 * a0_2;
277 Caab = a0 * Cab + 4 * a1_3;
278 Kaab = a1 * Kab + 4 * a0_3;
279 Cabb = 4 * b1_3 + 3 * b1_2 * b0 + 2 * b1 * b0_2 + b0_3;
280 Kabb = b1_3 + 2 * b1_2 * b0 + 3 * b1 * b0_2 + 4 * b0_3;
281
282 mP1 += db * C1;
283 mPa += db * Ca;
284 mPaa += db * Caa;
285 mPaaa += db * Caaa;
286 mPb += da * Cb;
287 mPbb += da * Cbb;
288 mPbbb += da * Cbbb;
289 mPab += db * (b1 * Cab + b0 * Kab);
290 mPaab += db * (b1 * Caab + b0 * Kaab);
291 mPabb += da * (a1 * Cabb + a0 * Kabb);
292 }
293
294 mP1 /= 2.0;
295 mPa /= 6.0;
296 mPaa /= 12.0;
297 mPaaa /= 20.0;
298 mPb /= -6.0;
299 mPbb /= -12.0;
300 mPbbb /= -20.0;
301 mPab /= 24.0;
302 mPaab /= 60.0;
303 mPabb /= -60.0;
304 }
305};
306
307} // namespace vcl
308
309#endif // VCL_SPACE_COMPLEX_MESH_INERTIA_H
vcl::MeshInertia
The MeshInertia class provides several for computing Polyhedral Mass properties (like inertia tensor,...
Definition mesh_inertia.h:56
vcl::MeshInertia::inertiaTensorEigen
void inertiaTensorEigen(Matrix33 &eigenValues, Point3< ScalarType > &eigenVector) const
Computes the inertia tensor mesh.
Definition mesh_inertia.h:174
vcl::MeshInertia::volume
ScalarType volume() const
Returns the volume (or mass) of the mesh. Meaningful only if the mesh is watertight.
Definition mesh_inertia.h:129
vcl::MeshInertia::projectionIntegrals
void projectionIntegrals(const FaceType &f)
Computes various integrations over projection of the given face.
Definition mesh_inertia.h:238
vcl::Segment
A class representing a line segment in n-dimensional space. The class is parameterized by a PointConc...
Definition segment.h:43
vcl::faceNormal
FaceType::VertexType::CoordType faceNormal(const FaceType &f)
Computes the normal of a face, without modifying the face. Works both for triangle and polygonal face...
Definition geometry.h:45
vcl::faceArea
auto faceArea(const FaceType &f)
Computes the area of a face. Works both for triangle and polygonal faces, and it is optimized in case...
Definition geometry.h:101
vcl::Point3d
Point3< double > Point3d
A convenience alias for a 3-dimensional Point with double-precision floating-point components.
Definition point.h:804