git.haldean.org d / a9824e3
start working on delaunay in so3 Haldean Brown 3 years ago
3 changed file(s) with 124 addition(s) and 0 deletion(s). Raw diff Collapse all Expand all
0 EIGEN_DIR=/usr/include/eigen3
1
2 libal: main.cpp
3 g++ -O0 -Wall -Werror -ggdb -fsanitize=address -o $@ -isystem $(EIGEN_DIR) $<
0 #include <iostream>
1 #include <fstream>
2 #include <vector>
3 #include <Eigen/Dense>
4
5 typedef Eigen::Matrix3Xd SO3Set;
6 typedef Eigen::Vector3d SO3;
7 typedef size_t PointId;
8 typedef std::tuple<PointId, PointId, PointId> Triangle;
9
10 struct SO3Mesh
11 {
12 SO3Set points;
13 std::vector<Triangle> triangles;
14 };
15
16 bool InsideCircumcircle(
17 SO3Mesh mesh, PointId t0, PointId t1, PointId t2, PointId p)
18 {
19 Eigen::Matrix3d circ;
20 Eigen::Vector3d A = mesh.points.row(t0);
21 Eigen::Vector3d B = mesh.points.row(t1);
22 Eigen::Vector3d C = mesh.points.row(t2);
23 Eigen::Vector3d D = mesh.points.row(p);
24 circ(0, 0) = A[0] - D[0];
25 circ(0, 1) = A[1] - D[1];
26 circ(0, 2) = A[0] * A[0] - D[0] * D[0] + A[1] * A[1] - D[1] * D[1];
27 circ(1, 0) = B[0] - D[0];
28 circ(1, 1) = B[1] - D[1];
29 circ(1, 2) = B[0] * B[0] - D[0] * D[0] + B[1] * B[1] - D[1] * D[1];
30 circ(2, 0) = C[0] - D[0];
31 circ(2, 1) = C[1] - D[1];
32 circ(2, 2) = C[0] * C[0] - D[0] * D[0] + C[1] * C[1] - D[1] * D[1];
33 return circ.determinant() > 0;
34 }
35
36 void Triangulate(SO3Mesh &mesh)
37 {
38 double theta_min, theta_max, phi_min, phi_max;
39 double theta_range, phi_range;
40 double theta_center, phi_center;
41 double r_max;
42
43 theta_min = mesh.points.col(0).minCoeff();
44 theta_max = mesh.points.col(0).maxCoeff();
45 phi_min = mesh.points.col(1).minCoeff();
46 phi_max = mesh.points.col(1).maxCoeff();
47 r_max = mesh.points.col(2).maxCoeff();
48
49 theta_range = theta_max - theta_min;
50 phi_range = phi_max - phi_min;
51
52 theta_center = (theta_max + theta_min) / 2.;
53 phi_center = (phi_max + phi_min) / 2.;
54
55 /* add 3 fake points to form the supertriangle */
56 PointId fpbase = mesh.points.cols();
57 mesh.points.resize(Eigen::NoChange, fpbase + 3);
58
59 mesh.points.col(fpbase + 0) = SO3(
60 theta_min - 10 * theta_range,
61 phi_min - phi_range,
62 r_max);
63 mesh.points.col(fpbase + 1) = SO3(
64 theta_center,
65 phi_max + 10 * phi_range,
66 r_max);
67 mesh.points.col(fpbase + 2) = SO3(
68 theta_max + 10 * theta_range,
69 phi_min - phi_range,
70 r_max);
71
72 std::vector<bool> ok;
73 triangles.push_back(Triangle(fpbase, fpbase + 1, fpbase + 2));
74 ok.push_back(false);
75
76 for (PointId p = 0; p < fpbase; p++)
77 {
78
79 }
80 }
81
82 void ReadFile(SO3Set &res, std::istream &input)
83 {
84 size_t r = 0;
85 size_t n_cols = std::count(
86 std::istreambuf_iterator<char>(input),
87 std::istreambuf_iterator<char>(),
88 '\n');
89 input.seekg(0);
90 res.resize(Eigen::NoChange, n_cols);
91 std::string line;
92 while (std::getline(input, line))
93 {
94 double theta, phi, r;
95 sscanf(line.c_str(), "%lf %lf %lf", &theta, &phi, &r);
96 res(0, r) = theta; res(1, r) = phi; res(2, r) = r;
97 r++;
98 }
99 std::cout << "loaded points:" << std::endl << res << std::endl;
100 }
101
102 int main() {
103 SO3Mesh mesh;
104 std::ifstream infile("simple-points.lpf");
105 ReadFile(mesh.points, infile);
106 return 0;
107 }
0 0 0 20
1 10 0 20
2 20 0 20
3 30 0 20
4 0 10 20
5 10 10 20
6 20 10 20
7 30 10 20
8 0 -10 20
9 10 -10 20
10 20 -10 20
11 30 -10 20