tetrahedron.cc
1.81 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
#include<vector>
#include<iostream>
#include<cmath>
#include<algorithm>
#include <Dense>
#include <vector>
using MatrixXd = Eigen::MatrixXd;
using Eigen::SelfAdjointEigenSolver;
using Vector = std::vector<double>;
using Matrix = std::vector<Vector>;
MatrixXd concatenate_vectors(Vector p1, Vector p2, Vector p3, Vector p4) {
MatrixXd mat(4, p1.size()); // create a MatrixXd object with 4 rows and p1.size() columns
mat.row(0) = Eigen::Map<const MatrixXd>(p1.data(), 1, p1.size());
mat.row(1) = Eigen::Map<const MatrixXd>(p2.data(), 1, p2.size());
mat.row(2) = Eigen::Map<const MatrixXd>(p3.data(), 1, p3.size());
mat.row(3) = Eigen::Map<const MatrixXd>(p4.data(), 1, p4.size());
return mat;
}
Vector getEigenvalues(const MatrixXd& matrix) {
SelfAdjointEigenSolver<MatrixXd> solver(matrix);
Vector eigenvalues(3);
eigenvalues[0] = solver.eigenvalues()[0];
eigenvalues[1] = solver.eigenvalues()[1];
eigenvalues[2] = solver.eigenvalues()[2];
return eigenvalues;
}
void substractBarycenter(Vector &v1, Vector &v2, Vector &v3, Vector &v4) {
int n = v1.size();
for (int i = 0; i < n; i++) {
double baryCenter= (v1[i] + v2[i] + v3[i] + v4[i]) / 4.0;
v1[i] -= baryCenter;
v2[i] -= baryCenter;
v3[i] -= baryCenter;
v4[i] -= baryCenter;
}
}
void compute_tetrahedron(Vector p1, Vector p2, Vector p3, Vector p4, double& elongation, double& planarity, double& characteristic) {
substractBarycenter(p1,p2,p3,p4);
// make matrix by concatenating vectors
MatrixXd mat = concatenate_vectors(p1,p2,p3,p4);
MatrixXd result = mat.transpose()*mat;
// computing eigenvalues
Vector eig = getEigenvalues(result);
// sort eigenvalues
sort(eig.begin(), eig.end());
planarity = 1. - sqrt(eig[0] / eig[1]);
elongation = 1. - sqrt(eig[1] / eig[2]);
characteristic = 2 * sqrt(eig[2]);
}