Commit 9f925679 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

grid_icosphere.cc: refactor.

parent 810cb336
......@@ -8,19 +8,103 @@
#include <array>
#include <map>
// Converts cartesian coordinates to geographical
static inline void
cc2gc(double xyz[], double *lon, double *lat)
class Vector3d
{
*lon = atan2(xyz[1], xyz[0]);
*lat = M_PI_2 - acos(xyz[2]);
}
private:
double X, Y, Z;
public:
// Sets all members to zero
Vector3d() : X(0), Y(0), Z(0) {};
Vector3d(const double &x, const double &y, const double &z) : X(x), Y(y), Z(z) {}
Vector3d(const double coords[3]) : X(coords[0]), Y(coords[1]), Z(coords[2]) {}
double operator[] (size_t index) const
{
switch (index)
{
case 0: return X;
case 1: return Y;
case 2: return Z;
}
return double();
}
Vector3d operator+ (const Vector3d &other) const
{
return Vector3d(X + other.X, Y + other.Y, Z + other.Z);
}
Vector3d operator- (const Vector3d &other) const
{
return Vector3d(X - other.X, Y - other.Y, Z - other.Z);
}
Vector3d operator- (void) const
{
return Vector3d(-X, -Y, -Z);
}
// Calculate the cross/outer/vector product
Vector3d operator% (const Vector3d &other) const
{
return Vector3d(Y*other.Z - Z*other.Y,
Z*other.X - X*other.Z,
X*other.Y - Y*other.X);
}
// Division by scalars
Vector3d operator/ (const double scalar) const
{
return Vector3d(X / scalar, Y / scalar, Z / scalar);
}
Vector3d operator/= (const double scalar)
{
return *this = *this / scalar;
}
// Calculate the dot/inner/scalar product
double operator* (const Vector3d &other) const
{
return (X * other.X) + (Y* other.Y) + (Z * other.Z);
}
double magnitude() const
{
return sqrt((X * X) + (Y * Y) + (Z * Z));
}
Vector3d normalised() const
{
return Vector3d(*this) / magnitude();
}
/*
Vector3d d_normalize()
{
const double dnorm = sqrt(*this * *this);
return Vector3d(*this) / dnorm;
}
*/
double longitude() const
{
return atan2(Y, X);
}
double latitude() const
{
return M_PI_2 - acos(Z);
}
};
using Index = size_t;
using vec3 = std::array<double, 3>;
using Vertex = Vector3d;
using Triangle = std::array<Index, 3>;
using TriangleList = std::vector<Triangle>;
using VertexList = std::vector<vec3>;
using VertexList = std::vector<Vertex>;
// icosahedron from ICON
namespace icosahedron
......@@ -37,8 +121,8 @@ init(void)
const double z_w = 2.0 * acos(1.0 / (2.0 * sin(pi_5)));
// set poles first - it is simple
vertices[0] = vec3{ {0.0, 0.0, 1.0} };
vertices[11] = vec3{ {0.0, 0.0, -1.0} };
vertices[0] = Vertex{ 0.0, 0.0, 1.0 };
vertices[11] = Vertex{ 0.0, 0.0, -1.0 };
// now set the vertices on the two latitude rings
int i_mdist[10];
for (int j = 1; j < 11; ++j)
......@@ -56,41 +140,17 @@ init(void)
// compute the meridian angle for the base vertex.
const double z_rlon = (1.0 + i_mdist[j - 1]) * pi_5;
// now initialize the coordinates
vertices[j] = vec3{ {sin(z_w) * cos(z_rlon), sin(z_w) * sin(z_rlon), cos(z_w) * i_msgn} };
vertices[j] = Vertex{ sin(z_w) * cos(z_rlon), sin(z_w) * sin(z_rlon), cos(z_w) * i_msgn };
}
}
// 20 triangles
static const TriangleList triangles
= { {{ 0, 1, 2 }}, {{ 0, 2, 3 }}, {{ 0, 3, 4 }}, {{ 0, 4, 5 }}, {{ 0, 5, 1 }}, {{ 6, 2, 1 }}, {{ 7, 3, 2 }},
{{ 8, 4, 3 }}, {{ 9, 5, 4 }}, {{ 10, 1, 5 }}, {{ 2, 6, 7 }}, {{ 3, 7, 8 }}, {{ 4, 8, 9 }}, {{ 5, 9, 10 }},
{{ 1, 10, 6 }}, {{ 11, 7, 6 }}, {{ 11, 8, 7 }}, {{ 11, 9, 8 }}, {{ 11, 10, 9 }}, {{ 11, 6, 10 }} };
= { { 0, 1, 2 }, { 0, 2, 3 }, { 0, 3, 4 }, { 0, 4, 5 }, { 0, 5, 1 }, { 6, 2, 1 }, { 7, 3, 2 },
{ 8, 4, 3 }, { 9, 5, 4 }, { 10, 1, 5 }, { 2, 6, 7 }, { 3, 7, 8 }, { 4, 8, 9 }, { 5, 9, 10 },
{ 1, 10, 6 }, { 11, 7, 6 }, { 11, 8, 7 }, { 11, 9, 8 }, { 11, 10, 9 }, { 11, 6, 10 } };
} // namespace icosahedron
static inline vec3
addVector(const vec3 &a, const vec3 &b)
{
vec3 c;
for (unsigned i = 0; i < 3; ++i) c[i] = a[i] + b[i];
return c;
}
static inline vec3
subVector(const vec3 &a, const vec3 &b)
{
vec3 c;
for (unsigned i = 0; i < 3; ++i) c[i] = a[i] - b[i];
return c;
}
static inline vec3
normalizeVector(const vec3 &a)
{
vec3 c;
const double magnitude = sqrt(a[0] * a[0] + a[1] * a[1] + a[2] * a[2]);
for (unsigned i = 0; i < 3; ++i) c[i] = a[i] / magnitude;
return c;
}
using Lookup = std::map<std::pair<Index, Index>, Index>;
......@@ -102,7 +162,7 @@ vertexForEdge(Lookup &lookup, VertexList &vertices, Index first, Index second)
const auto inserted = lookup.insert({ key, vertices.size() });
if (inserted.second)
vertices.push_back(normalizeVector(addVector(vertices[first], vertices[second])));
vertices.push_back((vertices[first] + vertices[second]).normalised());
return inserted.first->second;
}
......@@ -142,63 +202,42 @@ makeIcosphere(int subdivisions)
return { vertices, triangles };
}
static inline vec3
vectorProduct(const vec3 &a, const vec3 &b)
{
vec3 c;
c[0] = a[1] * b[2] - a[2] * b[1];
c[1] = a[2] * b[0] - a[0] * b[2];
c[2] = a[0] * b[1] - a[1] * b[0];
return c;
}
static inline double
dotProduct(const vec3 &a, const vec3 &b)
{
double sum = 0;
for (unsigned i = 0; i < 3; ++i) sum += a[i] * b[i];
return sum;
}
static inline void
d_normalize(vec3 &v)
d_normalize(Vertex &v)
{
const double dnorm = sqrt(dotProduct(v, v));
for (unsigned i = 0; i < 3; ++i) v[i] /= dnorm;
const double dnorm = sqrt(v * v);
v /= dnorm;
}
static vec3
circumCenterMean(const vec3 &v0, const vec3 &v1, const vec3 &v2)
static Vertex
circumCenterMean(const Vertex &v0, const Vertex &v1, const Vertex &v2)
{
/*
v0, v1, v2: the coordinates of the three triangle vertices (unit vectors) in
v0, v1, v2: the coordinates of the three triangle vertices (_dmo,nit vectors) in
counter clockwise order center: the coordinates of the circumcenter unless co-linear
*/
// cu0, cu1, cu2: vector product of center: e1 x e2
vec3 e1, e2; // edges of the underlying planar triangle: v1-v0 ands v2-v0, respectively
Vertex e1, e2; // edges of the underlying planar triangle: v1-v0 ands v2-v0, respectively
e1 = subVector(v1, v0);
e2 = subVector(v2, v0);
vec3 cu0 = vectorProduct(e1, e2);
if (dotProduct(cu0, v0) < 0.0)
for (unsigned i = 0; i < 3; ++i) cu0[i] = -cu0[i];
e1 = v1 - v0;
e2 = v2 - v0;
Vertex cu0 = e1 % e2;
if ((cu0 * v0) < 0.0) cu0 = -cu0;
d_normalize(cu0);
e1 = subVector(v2, v1);
e2 = subVector(v0, v1);
vec3 cu1 = vectorProduct(e1, e2);
if (dotProduct(cu1, v1) < 0.0)
for (unsigned i = 0; i < 3; ++i) cu1[i] = -cu1[i];
e1 = v2 - v1;
e2 = v0 - v1;
Vertex cu1 = e1 % e2;
if ((cu1 * v1) < 0.0) cu1 = -cu1;
d_normalize(cu1);
e1 = subVector(v0, v2);
e2 = subVector(v1, v2);
vec3 cu2 = vectorProduct(e1, e2);
if (dotProduct(cu2, v2) < 0.0)
for (unsigned i = 0; i < 3; ++i) cu2[i] = -cu2[i];
e1 = v0 - v2;
e2 = v1 - v2;
Vertex cu2 = e1 % e2;
if ((cu2 * v2) < 0.0) cu2 = -cu2;
d_normalize(cu2);
vec3 center = addVector(addVector(cu0, cu1), cu2);
Vertex center = cu0 + cu1 + cu2;
d_normalize(center);
return center;
......@@ -224,10 +263,16 @@ genIcosphereCoords(int subdivisions, bool lbounds, std::vector<double> &xvals, s
size_t i = 0;
for (const Triangle &t : triangles)
{
vec3 center = circumCenterMean(vertices[t[0]], vertices[t[1]], vertices[t[2]]);
cc2gc(&center[0], &xvals[i], &yvals[i]);
Vertex center = circumCenterMean(vertices[t[0]], vertices[t[1]], vertices[t[2]]);
xvals[i] = center.longitude();
yvals[i] = center.latitude();
if (lbounds)
for (size_t k = 0; k < 3; ++k) cc2gc(&vertices[t[k]][0], &xbounds[i * 3 + k], &ybounds[i * 3 + k]);
for (size_t k = 0; k < 3; ++k)
{
xbounds[i * 3 + k] = vertices[t[k]].longitude();
ybounds[i * 3 + k] = vertices[t[k]].latitude();
}
i++;
}
......@@ -247,23 +292,24 @@ main(void)
IndexedMesh mesh = makeIcosphere(0);
VertexList &vertices = mesh.first;
TriangleList &triangles = mesh.second;
if (1)
{
for (vec3 &v : vertices)
for (Vertex &v : vertices)
{
double lon, lat;
cc2gc(&v[0], &lon, &lat);
double lon = v.longitude();
double lat = v.latitude();
fprintf(stderr, "xyz:%g %g %g lon:%g lat:%g\n", v[0], v[1], v[2], rad2deg(lon), rad2deg(lat));
}
for (Triangle &t : triangles)
{
fprintf(stderr, "index: %d %d %d\n", t[0], t[1], t[2]);
fprintf(stderr, "index: %zu %zu %zu\n", t[0], t[1], t[2]);
}
for (Triangle &t : triangles)
{
double lon, lat;
vec3 center = circumCenterMean(vertices[t[0]], vertices[t[1]], vertices[t[2]]);
cc2gc(&center[0], &lon, &lat);
Vertex center = circumCenterMean(vertices[t[0]], vertices[t[1]], vertices[t[2]]);
double lon = center.longitude();
double lat = center.latitude();
fprintf(stderr, "center: %g %g\n", rad2deg(lon), rad2deg(lat));
}
for (Triangle &t : triangles)
......@@ -272,10 +318,12 @@ main(void)
printf(">\n");
for (int i = 0; i < 3; ++i)
{
cc2gc(&vertices[t[i]][0], &lon, &lat);
lon = vertices[t[i]].longitude();
lat = vertices[t[i]].latitude();
printf(" %g %g\n", rad2deg(lon), rad2deg(lat));
}
cc2gc(&vertices[t[0]][0], &lon, &lat);
lon = vertices[t[0]].longitude();
lat = vertices[t[0]].latitude();
printf(" %g %g\n", rad2deg(lon), rad2deg(lat));
}
}
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment