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
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
|
#pragma once
#include "glGrib/Geometry.h"
#include "glGrib/Trigonometry.h"
#include "glGrib/Options.h"
#include "glGrib/Handle.h"
#include "glGrib/OpenGL.h"
#include "glGrib/Buffer.h"
#include <glm/glm.hpp>
#include <vector>
namespace glGrib
{
class GeometryGaussian : public Geometry
{
public:
explicit GeometryGaussian (HandlePtr);
explicit GeometryGaussian (int);
bool isEqual (const Geometry &) const override;
void getPointNeighbours (int, std::vector<int> *) const override;
int latlon2index (float, float) const override;
void index2latlon (int, float *, float *) const override;
int size () const override;
void applyNormScale (glGrib::BufferPtr<float> &) const override;
void applyUVangle (glGrib::BufferPtr<float> &) const override;
void sample (OpenGLBufferPtr<unsigned char> &, const unsigned char, const int) const override;
void sampleTriangle (BufferPtr<unsigned char> &, const unsigned char, const int) const override;
float resolution (int level = 0) const override
{
if (level == 0)
level = grid_gaussian.Nj;
return M_PI / level;
}
void getTriangleNeighbours (int, int [3], int [3], glm::vec3 xyz[3]) const override;
void getTriangleNeighbours (int, int [3], int [3], glm::vec2 [3]) const override;
bool triangleIsEdge (int) const override;
int getTriangle (float, float) const override;
const glm::vec2 xyz2conformal (const glm::vec3 &) const override;
const glm::vec3 conformal2xyz (const glm::vec2 &) const override;
const glm::vec2 conformal2latlon (const glm::vec2 &) const override;
void fixPeriodicity (const glm::vec2 &, glm::vec2 *, int) const override;
float getLocalMeshSize (int) const override;
void getView (View *) const override;
void setProgramParameters (Program *) const override;
private:
void setup (HandlePtr, const OptionsGeometry &) override;
const std::string md5 () const override;
class jlonlat_t
{
public:
jlonlat_t () : jlon (0), jlat (0) {}
explicit jlonlat_t (int _jlon, int _jlat) : jlon (_jlon), jlat (_jlat) {}
int jlon = 0;
int jlat = 0;
};
const jlonlat_t jlonlat (int) const;
const glm::vec2 jlonlat2merc (const jlonlat_t & jlonlat) const
{
int jlat = jlonlat.jlat;
int jlon = jlonlat.jlon;
float coordy = misc_gaussian.latgauss[jlat-1];
float coordx = 2. * M_PI * (float)(jlon-1) / (float)grid_gaussian.pl[jlat-1];
return glm::vec2 (coordx, std::log (std::tan (M_PI / 4.0f + coordy / 2.0f)));
}
const glm::vec3 jlonlat2xyz (const jlonlat_t & jlonlat) const
{
int jlat = jlonlat.jlat;
int jlon = jlonlat.jlon;
float coordy = misc_gaussian.latgauss[jlat-1];
float sincoordy = std::sin (coordy);
float lat = std::asin ((misc_gaussian.omc2 + sincoordy * misc_gaussian.opc2) / (misc_gaussian.opc2 + sincoordy * misc_gaussian.omc2));
float coslat = std::cos (lat); float sinlat = std::sin (lat);
float coordx = 2. * M_PI * (float)(jlon-1) / (float)grid_gaussian.pl[jlat-1];
float lon = coordx;
float coslon = std::cos (lon); float sinlon = std::sin (lon);
float X = coslon * coslat;
float Y = sinlon * coslat;
float Z = sinlat;
if (! misc_gaussian.rotated)
return glm::vec3 (X, Y, Z);
return misc_gaussian.rot * glm::vec3 (X, Y, Z);
}
const glm::vec2 jlonlat2lonlat (const jlonlat_t & jlonlat) const
{
if ((misc_gaussian.stretchingFactor == 1.0f) && (! misc_gaussian.rotated))
{
int jlat = jlonlat.jlat, jlon = jlonlat.jlon;
float lat = misc_gaussian.latgauss[jlat-1];
float lon = 2. * M_PI * (float)(jlon-1) / (float)grid_gaussian.pl[jlat-1];
return glm::vec2 (lon, lat);
}
glm::vec3 xyz = jlonlat2xyz (jlonlat);
float lon = atan2 (xyz.y, xyz.x);
float lat = std::asin (xyz.z);
return glm::vec2 (lon, lat);
}
int getUpperTriangle (int jglo, const jlonlat_t & jlonlat) const;
int getLowerTriangle (int jglo, const jlonlat_t & jlonlat) const;
void getTriangleVertices (int it, int jglo[3]) const;
void getTriangleNeighbours (int, int [3], int [3], jlonlat_t [3]) const;
void latlon2coordxy (float, float, float &, float &) const;
int latlon2jlatjlon (float, float, int &, int &) const;
int computeLowerTriangle (int, int) const;
int computeUpperTriangle (int, int) const;
void computeTriangleVertices (int, int [3]) const;
void checkTriangleComputation () const;
void setupSSBO ();
void setupCoordinates ();
void setupFitLatitudes ();
class latfit_t
{
public:
int kind = 0;
std::vector<float> coeff;
float error = 0;
};
void tryFitLatitudes (int, latfit_t *);
GeometryGaussian () {}
void triangulate ();
void setupSubGrid ();
private:
// Grid
struct
{
std::vector<long int> pl;
long int Nj = 0;
std::vector<int> jglooff;
BufferPtr<unsigned int> ind;
BufferPtr<int> indcnt_per_lat;
BufferPtr<int> indoff_per_lat;
BufferPtr<int> triu; // Rank of triangle above
BufferPtr<int> trid; // Rank of triangle below
} grid_gaussian;
// Coordinates
struct
{
int kind = 0; // 1=linear, 2=octahedral
std::vector<float> latfitcoeff;
double stretchingFactor = 1.0f;
double latitudeOfStretchingPoleInDegrees = 90.0f;
double longitudeOfStretchingPoleInDegrees = 0.0f;
float omc2 = 0.0f;
float opc2 = 2.0f;
glm::mat3 rot = glm::mat3 (1.0f);
bool rotated = false;
BufferPtr<double> latgauss;
OpenGLBufferPtr<int> ssbo_jglo, ssbo_jlat;
OpenGLBufferPtr<float> ssbo_glat;
} misc_gaussian;
};
}
|
