1 | /* utilities related to polygons */ |
---|
2 | |
---|
3 | #include <cassert> |
---|
4 | #include <iostream> |
---|
5 | #include "elt.hpp" |
---|
6 | #include "errhandle.hpp" |
---|
7 | |
---|
8 | #include "polyg.hpp" |
---|
9 | |
---|
10 | namespace sphereRemap { |
---|
11 | |
---|
12 | using namespace std; |
---|
13 | |
---|
14 | /* given `N` `vertex`es, `N` `edge`s and `N` `d`s (d for small circles) |
---|
15 | and `g` the barycentre, |
---|
16 | this function reverses the order of arrays of `vertex`es, `edge`s and `d`s |
---|
17 | but only if it is required! |
---|
18 | (because computing intersection area requires both polygons to have same orientation) |
---|
19 | */ |
---|
20 | void orient(int N, Coord *vertex, Coord *edge, double *d, const Coord &g) |
---|
21 | { |
---|
22 | Coord ga = vertex[0] - g; |
---|
23 | Coord gb = vertex[1] - g; |
---|
24 | Coord vertical = crossprod(ga, gb); |
---|
25 | if (N > 2 && scalarprod(g, vertical) < 0) // (GAxGB).G |
---|
26 | { |
---|
27 | for (int i = 0; i < N/2; i++) |
---|
28 | swap(vertex[i], vertex[N-1-i]); |
---|
29 | |
---|
30 | for (int i = 0; i < (N-1)/2; i++) |
---|
31 | { |
---|
32 | swap(edge[N-2-i], edge[i]); |
---|
33 | swap(d[i], d[N-2-i]); |
---|
34 | } |
---|
35 | } |
---|
36 | } |
---|
37 | |
---|
38 | |
---|
39 | void normals(Coord *x, int n, Coord *a) |
---|
40 | { |
---|
41 | for (int i = 0; i<n; i++) |
---|
42 | a[i] = crossprod(x[(i+1)%n], x[i]); |
---|
43 | } |
---|
44 | |
---|
45 | Coord barycentre(const Coord *x, int n) |
---|
46 | { |
---|
47 | if (n == 0) return ORIGIN; |
---|
48 | Coord bc = ORIGIN; |
---|
49 | for (int i = 0; i < n; i++) |
---|
50 | bc = bc + x[i]; |
---|
51 | /* both distances can be equal down to roundoff when norm(bc) < mashineepsilon |
---|
52 | which can occur when weighted with tiny area */ |
---|
53 | |
---|
54 | assert(squaredist(bc, proj(bc)) <= squaredist(bc, proj(bc * (-1.0)))); |
---|
55 | //if (squaredist(bc, proj(bc)) > squaredist(bc, proj(bc * (-1.0)))) return proj(bc * (-1.0)); |
---|
56 | |
---|
57 | return proj(bc); |
---|
58 | } |
---|
59 | |
---|
60 | /** computes the barycentre of the area which is the difference |
---|
61 | of a side ABO of the spherical tetrahedron and of the straight tetrahedron */ |
---|
62 | static Coord tetrah_side_diff_centre(Coord a, Coord b) |
---|
63 | { |
---|
64 | Coord n = crossprod(a,b); |
---|
65 | double sinc2 = n.x*n.x + n.y*n.y + n.z*n.z; |
---|
66 | assert(sinc2 < 1.0 + EPS); |
---|
67 | |
---|
68 | // exact: u = asin(sinc)/sinc - 1; asin(sinc) = geodesic length of arc ab |
---|
69 | // approx: |
---|
70 | // double u = sinc2/6. + (3./40.)*sinc2*sinc2; |
---|
71 | |
---|
72 | // exact |
---|
73 | if (sinc2 > 1.0 - EPS) /* if round-off leads to sinc > 1 asin produces NaN */ |
---|
74 | return n * (M_PI_2 - 1); |
---|
75 | double sinc = sqrt(sinc2); |
---|
76 | double u = asin(sinc)/sinc - 1; |
---|
77 | |
---|
78 | return n*u; |
---|
79 | } |
---|
80 | |
---|
81 | /* compute the barycentre as the negative sum of all the barycentres of sides of the tetrahedron */ |
---|
82 | Coord gc_normalintegral(const Coord *x, int n) |
---|
83 | { |
---|
84 | Coord m = barycentre(x, n); |
---|
85 | Coord bc = crossprod(x[n-1]-m, x[0]-m) + tetrah_side_diff_centre(x[n-1], x[0]); |
---|
86 | for (int i = 1; i < n; i++) |
---|
87 | bc = bc + crossprod(x[i-1]-m, x[i]-m) + tetrah_side_diff_centre(x[i-1], x[i]); |
---|
88 | return bc*0.5; |
---|
89 | } |
---|
90 | |
---|
91 | Coord exact_barycentre(const Coord *x, int n) |
---|
92 | { |
---|
93 | if (n >= 3) |
---|
94 | { |
---|
95 | return proj(gc_normalintegral(x, n)); |
---|
96 | } |
---|
97 | else if (n == 0) return ORIGIN; |
---|
98 | else if (n == 2) return midpoint(x[0], x[1]); |
---|
99 | else if (n == 1) return x[0]; |
---|
100 | } |
---|
101 | |
---|
102 | Coord sc_gc_moon_normalintegral(Coord a, Coord b, Coord pole) |
---|
103 | { |
---|
104 | double hemisphere = (a.z > 0) ? 1: -1; |
---|
105 | |
---|
106 | double lat = hemisphere * (M_PI_2 - acos(a.z)); |
---|
107 | double lon1 = atan2(a.y, a.x); |
---|
108 | double lon2 = atan2(b.y, b.x); |
---|
109 | double lon_diff = lon2 - lon1; |
---|
110 | |
---|
111 | // wraparound at lon=-pi=pi |
---|
112 | if (lon_diff < -M_PI) lon_diff += 2.0*M_PI; |
---|
113 | else if (lon_diff > M_PI) lon_diff -= 2.0*M_PI; |
---|
114 | |
---|
115 | // integral of the normal over the surface bound by great arcs a-pole and b-pole and small arc a-b |
---|
116 | Coord sc_normalintegral = Coord(0.5*(sin(lon2)-sin(lon1))*(M_PI_2 - lat - 0.5*sin(2.0*lat)), |
---|
117 | 0.5*(cos(lon1)-cos(lon2))*(M_PI_2 - lat - 0.5*sin(2.0*lat)), |
---|
118 | hemisphere * lon_diff * 0.25 * (cos(2.0*lat) + 1.0)); |
---|
119 | Coord p = Coord(0,0,hemisphere); // TODO assumes north pole is (0,0,1) |
---|
120 | Coord t[] = {a, b, p}; |
---|
121 | if (hemisphere < 0) swap(t[0], t[1]); |
---|
122 | return (sc_normalintegral - gc_normalintegral(t, 3)) * hemisphere; |
---|
123 | } |
---|
124 | |
---|
125 | |
---|
126 | double triarea(Coord& A, Coord& B, Coord& C) |
---|
127 | { |
---|
128 | double a = ds(B, C); |
---|
129 | double b = ds(C, A); |
---|
130 | double c = ds(A, B); |
---|
131 | double s = 0.5 * (a + b + c); |
---|
132 | double t = tan(0.5*s) * tan(0.5*(s - a)) * tan(0.5*(s - b)) * tan(0.5*(s - c)); |
---|
133 | // assert(t >= 0); |
---|
134 | if (t<1e-20) return 0. ; |
---|
135 | return 4 * atan(sqrt(t)); |
---|
136 | } |
---|
137 | |
---|
138 | /** Computes area of two two-sided polygon |
---|
139 | needs to have one small and one great circle, otherwise zero |
---|
140 | (name origin: lun is moon in french) |
---|
141 | */ |
---|
142 | double alun(double b, double d) |
---|
143 | { |
---|
144 | double a = acos(d); |
---|
145 | assert(b <= 2 * a); |
---|
146 | double s = a + 0.5 * b; |
---|
147 | double t = tan(0.5 * s) * tan(0.5 * (s - a)) * tan(0.5 * (s - a)) * tan(0.5 * (s - b)); |
---|
148 | double r = sqrt(1 - d*d); |
---|
149 | double p = 2 * asin(sin(0.5*b) / r); |
---|
150 | return p*(1 - d) - 4*atan(sqrt(t)); |
---|
151 | } |
---|
152 | |
---|
153 | /** |
---|
154 | This function returns the area of a spherical element |
---|
155 | that can be composed of great and small circle arcs. |
---|
156 | The caller must ensure this function is not called when `alun` should be called instaed. |
---|
157 | This function also sets `gg` to the barycentre of the element. |
---|
158 | "air" stands for area and "bar" for barycentre. |
---|
159 | */ |
---|
160 | double airbar(int N, const Coord *x, const Coord *c, double *d, const Coord& pole, Coord& gg) |
---|
161 | { |
---|
162 | if (N < 3) |
---|
163 | return 0; /* polygons with less then three vertices have zero area */ |
---|
164 | Coord t[3]; |
---|
165 | t[0] = barycentre(x, N); |
---|
166 | Coord *g = new Coord[N]; |
---|
167 | double area = 0; |
---|
168 | Coord gg_exact = gc_normalintegral(x, N); |
---|
169 | for (int i = 0; i < N; i++) |
---|
170 | { |
---|
171 | /* for "spherical circle segment" sum triangular part and "small moon" and => account for small circle */ |
---|
172 | int ii = (i + 1) % N; |
---|
173 | t[1] = x[i]; |
---|
174 | t[2] = x[ii]; |
---|
175 | assert(scalarprod(crossprod(t[1] - t[0], t[2] - t[0]), t[0]) >= 0); // Error: tri a l'env (wrong orientation) |
---|
176 | double area_gc = triarea(t[0], t[1], t[2]); |
---|
177 | double area_sc_gc_moon = 0; |
---|
178 | if (d[i]) /* handle small circle case */ |
---|
179 | { |
---|
180 | Coord m = midpoint(t[1], t[2]); |
---|
181 | double mext = scalarprod(m, c[i]) - d[i]; |
---|
182 | char sgl = (mext > 0) ? -1 : 1; |
---|
183 | area_sc_gc_moon = sgl * alun(arcdist(t[1], t[2]), fabs(scalarprod(t[1], pole))); |
---|
184 | gg_exact = gg_exact + sc_gc_moon_normalintegral(t[1], t[2], pole); |
---|
185 | } |
---|
186 | area += area_gc + area_sc_gc_moon; /* for "spherical circle segment" sum triangular part (at) and "small moon" and => account for small circle */ |
---|
187 | g[i] = barycentre(t, 3) * (area_gc + area_sc_gc_moon); |
---|
188 | } |
---|
189 | gg = barycentre(g, N); |
---|
190 | gg_exact = proj(gg_exact); |
---|
191 | delete[] g; |
---|
192 | gg = gg_exact; |
---|
193 | return area; |
---|
194 | } |
---|
195 | |
---|
196 | double polygonarea(Coord *vertices, int N) |
---|
197 | { |
---|
198 | assert(N >= 3); |
---|
199 | |
---|
200 | /* compute polygon area as sum of triangles */ |
---|
201 | Coord centre = barycentre(vertices, N); |
---|
202 | double area = 0; |
---|
203 | for (int i = 0; i < N; i++) |
---|
204 | area += triarea(centre, vertices[i], vertices[(i+1)%N]); |
---|
205 | return area; |
---|
206 | } |
---|
207 | |
---|
208 | int packedPolygonSize(const Elt& e) |
---|
209 | { |
---|
210 | return sizeof(e.id) + sizeof(e.src_id) + sizeof(e.x) + sizeof(e.val) + |
---|
211 | sizeof(e.n) + e.n*(sizeof(double)+sizeof(Coord)); |
---|
212 | } |
---|
213 | |
---|
214 | void packPolygon(const Elt& e, char *buffer, int& pos) |
---|
215 | { |
---|
216 | *((GloId *) &(buffer[pos])) = e.id; |
---|
217 | pos += sizeof(e.id); |
---|
218 | *((GloId *) &(buffer[pos])) = e.src_id; |
---|
219 | pos += sizeof(e.src_id); |
---|
220 | |
---|
221 | *((Coord *) &(buffer[pos])) = e.x; |
---|
222 | pos += sizeof(e.x); |
---|
223 | |
---|
224 | *((double*) &(buffer[pos])) = e.val; |
---|
225 | pos += sizeof(e.val); |
---|
226 | |
---|
227 | *((int *) & (buffer[pos])) = e.n; |
---|
228 | pos += sizeof(e.n); |
---|
229 | |
---|
230 | for (int i = 0; i < e.n; i++) |
---|
231 | { |
---|
232 | *((double *) & (buffer[pos])) = e.d[i]; |
---|
233 | pos += sizeof(e.d[i]); |
---|
234 | |
---|
235 | *((Coord *) &(buffer[pos])) = e.vertex[i]; |
---|
236 | pos += sizeof(e.vertex[i]); |
---|
237 | } |
---|
238 | |
---|
239 | } |
---|
240 | |
---|
241 | void unpackPolygon(Elt& e, const char *buffer, int& pos) |
---|
242 | { |
---|
243 | e.id = *((GloId *) &(buffer[pos])); |
---|
244 | pos += sizeof(e.id); |
---|
245 | e.src_id = *((GloId *) &(buffer[pos])); |
---|
246 | pos += sizeof(e.src_id); |
---|
247 | |
---|
248 | e.x = *((Coord *) & (buffer[pos]) ); |
---|
249 | pos += sizeof(e.x); |
---|
250 | |
---|
251 | e.val = *((double *) & (buffer[pos])); |
---|
252 | pos += sizeof(double); |
---|
253 | |
---|
254 | e.n = *((int *) & (buffer[pos])); |
---|
255 | pos += sizeof(int); |
---|
256 | |
---|
257 | for (int i = 0; i < e.n; i++) |
---|
258 | { |
---|
259 | e.d[i] = *((double *) & (buffer[pos])); |
---|
260 | pos += sizeof(double); |
---|
261 | |
---|
262 | e.vertex[i] = *((Coord *) & (buffer[pos])); |
---|
263 | pos += sizeof(Coord); |
---|
264 | } |
---|
265 | } |
---|
266 | |
---|
267 | int packIntersectionSize(const Elt& elt) |
---|
268 | { |
---|
269 | return elt.is.size() * (2*sizeof(int)+ sizeof(GloId) + 5*sizeof(double)); |
---|
270 | } |
---|
271 | |
---|
272 | void packIntersection(const Elt& e, char* buffer,int& pos) |
---|
273 | { |
---|
274 | for (list<Polyg *>::const_iterator it = e.is.begin(); it != e.is.end(); ++it) |
---|
275 | { |
---|
276 | *((int *) &(buffer[0])) += 1; |
---|
277 | |
---|
278 | *((int *) &(buffer[pos])) = e.id.ind; |
---|
279 | pos += sizeof(int); |
---|
280 | |
---|
281 | *((double *) &(buffer[pos])) = e.area; |
---|
282 | pos += sizeof(double); |
---|
283 | |
---|
284 | *((GloId *) &(buffer[pos])) = (*it)->id; |
---|
285 | pos += sizeof(GloId); |
---|
286 | |
---|
287 | *((int *) &(buffer[pos])) = (*it)->n; |
---|
288 | pos += sizeof(int); |
---|
289 | *((double *) &(buffer[pos])) = (*it)->area; |
---|
290 | pos += sizeof(double); |
---|
291 | |
---|
292 | *((Coord *) &(buffer[pos])) = (*it)->x; |
---|
293 | pos += sizeof(Coord); |
---|
294 | } |
---|
295 | } |
---|
296 | |
---|
297 | void unpackIntersection(Elt* e, const char* buffer) |
---|
298 | { |
---|
299 | int ind; |
---|
300 | int pos = 0; |
---|
301 | |
---|
302 | int n = *((int *) & (buffer[pos])); |
---|
303 | pos += sizeof(int); |
---|
304 | for (int i = 0; i < n; i++) |
---|
305 | { |
---|
306 | ind = *((int*) &(buffer[pos])); |
---|
307 | pos+=sizeof(int); |
---|
308 | |
---|
309 | Elt& elt= e[ind]; |
---|
310 | |
---|
311 | elt.area=*((double *) & (buffer[pos])); |
---|
312 | pos += sizeof(double); |
---|
313 | |
---|
314 | Polyg *polygon = new Polyg; |
---|
315 | |
---|
316 | polygon->id = *((GloId *) & (buffer[pos])); |
---|
317 | pos += sizeof(GloId); |
---|
318 | |
---|
319 | polygon->n = *((int *) & (buffer[pos])); |
---|
320 | pos += sizeof(int); |
---|
321 | |
---|
322 | polygon->area = *((double *) & (buffer[pos])); |
---|
323 | pos += sizeof(double); |
---|
324 | |
---|
325 | polygon->x = *((Coord *) & (buffer[pos])); |
---|
326 | pos += sizeof(Coord); |
---|
327 | |
---|
328 | elt.is.push_back(polygon); |
---|
329 | } |
---|
330 | } |
---|
331 | |
---|
332 | } |
---|