1 // modified code from:
2 // https://schneide.wordpress.com/2016/07/15/generating-an-icosphere-in-c
3
4 #include "vector3d.h"
5
6 #include <utility>
7 #include <limits>
8 #include <cstdio>
9 #include <iostream>
10 #include <vector>
11 #include <array>
12 #include <map>
13 #include <unordered_map>
14
15 using Index = size_t;
16 using Vertex = Vector3d;
17 using Triangle = std::array<Index, 3>;
18 using TriangleList = std::vector<Triangle>;
19 using VertexList = std::vector<Vertex>;
20
21 using Lookup = std::map<std::pair<Index, Index>, Index>;
22
23 // icosahedron from ICON
24 namespace icosahedron
25 {
26
27 // Northern hemisphere are the first 6 elements of vertices[0:5]
28 // Southern hemisphere are the other 6 elements of vertices[6:11]
29 // 12 vertices
30 static VertexList vertices(12);
31
32 void
init(void)33 init(void)
34 {
35 constexpr auto pi_5 = M_PI * 0.2;
36 // first define the vertices of the icosahedron
37 const auto z_w = 2.0 * std::acos(1.0 / (2.0 * std::sin(pi_5)));
38
39 // set poles first - it is simple
40 vertices[0] = Vertex{ 0.0, 0.0, 1.0 };
41 vertices[11] = Vertex{ 0.0, 0.0, -1.0 };
42 // now set the vertices on the two latitude rings
43 int i_mdist[10];
44 for (int j = 1; j < 11; ++j)
45 {
46 if (j % 2 == 0)
47 i_mdist[j / 2 + 4] = -1 + (j - 1) - 10 * ((j - 1) / 7);
48 else
49 i_mdist[(j + 1) / 2 - 1] = -1 + (j - 1) - 10 * ((j - 1) / 7);
50 }
51
52 for (int j = 1; j < 11; ++j)
53 {
54 // toggle the hemisphere
55 const auto i_msgn = (j >= 6) ? -1.0 : 1.0;
56 // compute the meridian angle for the base vertex.
57 const auto z_rlon = (1.0 + i_mdist[j - 1]) * pi_5;
58 // now initialize the coordinates
59 vertices[j] = Vertex{ std::sin(z_w) * std::cos(z_rlon), std::sin(z_w) * std::sin(z_rlon), std::cos(z_w) * i_msgn };
60 }
61 }
62
63 // 20 triangles
64 static const TriangleList triangles
65 = { { { 0, 1, 2 } }, { { 0, 2, 3 } }, { { 0, 3, 4 } }, { { 0, 4, 5 } }, { { 0, 5, 1 } },
66 { { 6, 2, 1 } }, { { 7, 3, 2 } }, { { 8, 4, 3 } }, { { 9, 5, 4 } }, { { 10, 1, 5 } },
67 { { 2, 6, 7 } }, { { 3, 7, 8 } }, { { 4, 8, 9 } }, { { 5, 9, 10 } }, { { 1, 10, 6 } },
68 { { 11, 7, 6 } }, { { 11, 8, 7 } }, { { 11, 9, 8 } }, { { 11, 10, 9 } }, { { 11, 6, 10 } } };
69
70 } // namespace icosahedron
71
72 static Index
vertexForEdge(Lookup & lookup,VertexList & vertices,Index first,Index second)73 vertexForEdge(Lookup &lookup, VertexList &vertices, Index first, Index second)
74 {
75 Lookup::key_type key(first, second);
76 if (key.first > key.second) std::swap(key.first, key.second);
77
78 const auto inserted = lookup.insert({ key, vertices.size() });
79 if (inserted.second) vertices.push_back((vertices[first] + vertices[second]).normalised());
80
81 return inserted.first->second;
82 }
83
84 static TriangleList
subdivide(VertexList & vertices,const TriangleList & triangles)85 subdivide(VertexList &vertices, const TriangleList &triangles)
86 {
87 Lookup lookup;
88 Triangle mid;
89
90 const auto n = triangles.size();
91 TriangleList result(4 * n);
92 for (size_t i = 0; i < n; ++i)
93 {
94 const auto &each = triangles[i];
95 for (int edge = 0; edge < 3; edge++)
96 {
97 mid[edge] = vertexForEdge(lookup, vertices, each[edge], each[(edge + 1) % 3]);
98 }
99
100 result[i * 4 + 0] = { each[0], mid[0], mid[2] };
101 result[i * 4 + 1] = { each[1], mid[1], mid[0] };
102 result[i * 4 + 2] = { each[2], mid[2], mid[1] };
103 result[i * 4 + 3] = mid;
104 }
105
106 return result;
107 }
108
109 using IndexedMesh = std::pair<VertexList, TriangleList>;
110
111 static IndexedMesh
makeIcosphere(int subdivisions)112 makeIcosphere(int subdivisions)
113 {
114 icosahedron::init();
115 auto vertices = icosahedron::vertices;
116 auto triangles = icosahedron::triangles;
117
118 for (int i = 0; i < subdivisions; ++i) triangles = subdivide(vertices, triangles);
119
120 return { vertices, triangles };
121 }
122
123 static inline void
d_normalize(Vertex & v)124 d_normalize(Vertex &v)
125 {
126 const auto dnorm = std::sqrt(v * v);
127 v /= dnorm;
128 }
129
130 static Vertex
circumCenterMean(const Vertex & v0,const Vertex & v1,const Vertex & v2)131 circumCenterMean(const Vertex &v0, const Vertex &v1, const Vertex &v2)
132 {
133 /*
134 v0, v1, v2: the coordinates of the three triangle vertices (_dmo,nit vectors) in
135 counter clockwise order center: the coordinates of the circumcenter unless co-linear
136 */
137 // cu0, cu1, cu2: vector product of center: e1 x e2
138 // e1, e2: edges of the underlying planar triangle: v1-v0 ands v2-v0, respectively
139
140 auto e1 = v1 - v0;
141 auto e2 = v2 - v0;
142 auto cu0 = e1 % e2;
143 if ((cu0 * v0) < 0.0) cu0 = -cu0;
144 d_normalize(cu0);
145
146 e1 = v2 - v1;
147 e2 = v0 - v1;
148 auto cu1 = e1 % e2;
149 if ((cu1 * v1) < 0.0) cu1 = -cu1;
150 d_normalize(cu1);
151
152 e1 = v0 - v2;
153 e2 = v1 - v2;
154 auto cu2 = e1 % e2;
155 if ((cu2 * v2) < 0.0) cu2 = -cu2;
156 d_normalize(cu2);
157
158 auto center = cu0 + cu1 + cu2;
159 d_normalize(center);
160
161 return center;
162 }
163
164 size_t
genIcosphereCoords(int subdivisions,bool lbounds,std::vector<double> & xvals,std::vector<double> & yvals,std::vector<double> & xbounds,std::vector<double> & ybounds)165 genIcosphereCoords(int subdivisions, bool lbounds, std::vector<double> &xvals, std::vector<double> &yvals,
166 std::vector<double> &xbounds, std::vector<double> &ybounds)
167 {
168 const auto mesh = makeIcosphere(subdivisions);
169 const auto &vertices = mesh.first;
170 const auto &triangles = mesh.second;
171
172 const auto ncells = triangles.size();
173 xvals.resize(ncells);
174 yvals.resize(ncells);
175 if (lbounds)
176 {
177 xbounds.resize(3 * ncells);
178 ybounds.resize(3 * ncells);
179 }
180
181 size_t i = 0;
182 for (const Triangle &t : triangles)
183 {
184 const auto center = circumCenterMean(vertices[t[0]], vertices[t[1]], vertices[t[2]]);
185 xvals[i] = center.longitude();
186 yvals[i] = center.latitude();
187 if (lbounds)
188 for (size_t k = 0; k < 3; ++k)
189 {
190 xbounds[i * 3 + k] = vertices[t[k]].longitude();
191 ybounds[i * 3 + k] = vertices[t[k]].latitude();
192 }
193
194 i++;
195 }
196
197 return ncells;
198 }
199
200 #ifdef TEST_ICO
201 #include <chrono>
202 static inline double
rad2deg(double v)203 rad2deg(double v)
204 {
205 return v * 180. / M_PI;
206 }
207
208 void
subDivdeAndVertexForEdgeTest(int subdivisions)209 subDivdeAndVertexForEdgeTest(int subdivisions)
210 {
211 icosahedron::init();
212 auto vertices = icosahedron::vertices;
213 auto triangles = icosahedron::triangles;
214
215 auto start = std::chrono::steady_clock::now();
216 for (int i = 0; i < subdivisions; ++i) triangles = subdivide(vertices, triangles);
217 auto end = std::chrono::steady_clock::now();
218 auto elapsedM = std::chrono::duration_cast<std::chrono::nanoseconds>(end - start).count();
219
220 size_t expectedNumTriangles = std::pow(4, subdivisions) * 20;
221 size_t expectedNumVerticies = std::pow(4, subdivisions) * 10 + 2;
222 if (triangles.size() != (size_t) expectedNumTriangles)
223 {
224 std::cout << "ERROR: unexpected number of triangles for T: std::map, expected " << expectedNumTriangles << " got "
225 << triangles.size() << std::endl;
226 }
227 if (vertices.size() != expectedNumVerticies)
228 {
229 std::cout << "ERROR: unexpected number of verticies for T: std::map, expected " << expectedNumVerticies << " got "
230 << vertices.size() << std::endl;
231 }
232 }
233
234 int
main(void)235 main(void)
236 {
237 int numSubdivisions = 0;
238 subDivdeAndVertexForEdgeTest(numSubdivisions);
239
240 auto mesh = makeIcosphere(numSubdivisions);
241 auto &vertices = mesh.first;
242 auto &triangles = mesh.second;
243
244 if (1)
245 {
246 for (Vertex &v : vertices)
247 {
248 auto lon = v.longitude();
249 auto lat = v.latitude();
250 fprintf(stderr, "xyz:%g %g %g lon:%g lat:%g\n", v[0], v[1], v[2], rad2deg(lon), rad2deg(lat));
251 }
252 for (Triangle &t : triangles)
253 {
254 fprintf(stderr, "index: %zu %zu %zu\n", t[0], t[1], t[2]);
255 }
256 for (Triangle &t : triangles)
257 {
258 auto center = circumCenterMean(vertices[t[0]], vertices[t[1]], vertices[t[2]]);
259 auto lon = center.longitude();
260 auto lat = center.latitude();
261 fprintf(stderr, "center: %g %g\n", rad2deg(lon), rad2deg(lat));
262 }
263 for (Triangle &t : triangles)
264 {
265 double lon, lat;
266 printf(">\n");
267 for (int i = 0; i < 3; ++i)
268 {
269 lon = vertices[t[i]].longitude();
270 lat = vertices[t[i]].latitude();
271 printf(" %g %g\n", rad2deg(lon), rad2deg(lat));
272 }
273 lon = vertices[t[0]].longitude();
274 lat = vertices[t[0]].latitude();
275 printf(" %g %g\n", rad2deg(lon), rad2deg(lat));
276 }
277 }
278 fprintf(stderr, "vertices %zu\n", vertices.size());
279 fprintf(stderr, "triangles %zu\n", triangles.size());
280 }
281 #endif
282