[688] | 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 */ |
---|
[845] | 53 | |
---|
[1328] | 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 | |
---|
[688] | 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)); |
---|
[849] | 133 | // assert(t >= 0); |
---|
| 134 | if (t<1e-20) return 0. ; |
---|
[688] | 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) |
---|
[1328] | 163 | return 0; /* polygons with less than three vertices have zero area */ |
---|
[688] | 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]; |
---|
[1328] | 175 | double sc=scalarprod(crossprod(t[1] - t[0], t[2] - t[0]), t[0]) ; |
---|
[1205] | 176 | assert(sc >= -1e-10); // Error: tri a l'env (wrong orientation) |
---|
[688] | 177 | double area_gc = triarea(t[0], t[1], t[2]); |
---|
| 178 | double area_sc_gc_moon = 0; |
---|
| 179 | if (d[i]) /* handle small circle case */ |
---|
| 180 | { |
---|
| 181 | Coord m = midpoint(t[1], t[2]); |
---|
| 182 | double mext = scalarprod(m, c[i]) - d[i]; |
---|
| 183 | char sgl = (mext > 0) ? -1 : 1; |
---|
| 184 | area_sc_gc_moon = sgl * alun(arcdist(t[1], t[2]), fabs(scalarprod(t[1], pole))); |
---|
| 185 | gg_exact = gg_exact + sc_gc_moon_normalintegral(t[1], t[2], pole); |
---|
| 186 | } |
---|
| 187 | area += area_gc + area_sc_gc_moon; /* for "spherical circle segment" sum triangular part (at) and "small moon" and => account for small circle */ |
---|
| 188 | g[i] = barycentre(t, 3) * (area_gc + area_sc_gc_moon); |
---|
| 189 | } |
---|
| 190 | gg = barycentre(g, N); |
---|
| 191 | gg_exact = proj(gg_exact); |
---|
| 192 | delete[] g; |
---|
| 193 | gg = gg_exact; |
---|
| 194 | return area; |
---|
| 195 | } |
---|
| 196 | |
---|
| 197 | double polygonarea(Coord *vertices, int N) |
---|
| 198 | { |
---|
| 199 | assert(N >= 3); |
---|
| 200 | |
---|
| 201 | /* compute polygon area as sum of triangles */ |
---|
| 202 | Coord centre = barycentre(vertices, N); |
---|
| 203 | double area = 0; |
---|
| 204 | for (int i = 0; i < N; i++) |
---|
| 205 | area += triarea(centre, vertices[i], vertices[(i+1)%N]); |
---|
| 206 | return area; |
---|
| 207 | } |
---|
| 208 | |
---|
| 209 | int packedPolygonSize(const Elt& e) |
---|
| 210 | { |
---|
| 211 | return sizeof(e.id) + sizeof(e.src_id) + sizeof(e.x) + sizeof(e.val) + |
---|
| 212 | sizeof(e.n) + e.n*(sizeof(double)+sizeof(Coord)); |
---|
| 213 | } |
---|
| 214 | |
---|
| 215 | void packPolygon(const Elt& e, char *buffer, int& pos) |
---|
| 216 | { |
---|
| 217 | *((GloId *) &(buffer[pos])) = e.id; |
---|
| 218 | pos += sizeof(e.id); |
---|
| 219 | *((GloId *) &(buffer[pos])) = e.src_id; |
---|
| 220 | pos += sizeof(e.src_id); |
---|
| 221 | |
---|
| 222 | *((Coord *) &(buffer[pos])) = e.x; |
---|
| 223 | pos += sizeof(e.x); |
---|
| 224 | |
---|
| 225 | *((double*) &(buffer[pos])) = e.val; |
---|
| 226 | pos += sizeof(e.val); |
---|
| 227 | |
---|
| 228 | *((int *) & (buffer[pos])) = e.n; |
---|
| 229 | pos += sizeof(e.n); |
---|
| 230 | |
---|
| 231 | for (int i = 0; i < e.n; i++) |
---|
| 232 | { |
---|
| 233 | *((double *) & (buffer[pos])) = e.d[i]; |
---|
| 234 | pos += sizeof(e.d[i]); |
---|
| 235 | |
---|
| 236 | *((Coord *) &(buffer[pos])) = e.vertex[i]; |
---|
| 237 | pos += sizeof(e.vertex[i]); |
---|
| 238 | } |
---|
| 239 | |
---|
| 240 | } |
---|
| 241 | |
---|
| 242 | void unpackPolygon(Elt& e, const char *buffer, int& pos) |
---|
| 243 | { |
---|
| 244 | e.id = *((GloId *) &(buffer[pos])); |
---|
| 245 | pos += sizeof(e.id); |
---|
| 246 | e.src_id = *((GloId *) &(buffer[pos])); |
---|
| 247 | pos += sizeof(e.src_id); |
---|
| 248 | |
---|
| 249 | e.x = *((Coord *) & (buffer[pos]) ); |
---|
| 250 | pos += sizeof(e.x); |
---|
| 251 | |
---|
| 252 | e.val = *((double *) & (buffer[pos])); |
---|
| 253 | pos += sizeof(double); |
---|
| 254 | |
---|
| 255 | e.n = *((int *) & (buffer[pos])); |
---|
| 256 | pos += sizeof(int); |
---|
| 257 | |
---|
| 258 | for (int i = 0; i < e.n; i++) |
---|
| 259 | { |
---|
| 260 | e.d[i] = *((double *) & (buffer[pos])); |
---|
| 261 | pos += sizeof(double); |
---|
| 262 | |
---|
| 263 | e.vertex[i] = *((Coord *) & (buffer[pos])); |
---|
| 264 | pos += sizeof(Coord); |
---|
| 265 | } |
---|
| 266 | } |
---|
| 267 | |
---|
| 268 | int packIntersectionSize(const Elt& elt) |
---|
| 269 | { |
---|
[845] | 270 | return elt.is.size() * (2*sizeof(int)+ sizeof(GloId) + 5*sizeof(double)); |
---|
[688] | 271 | } |
---|
| 272 | |
---|
| 273 | void packIntersection(const Elt& e, char* buffer,int& pos) |
---|
| 274 | { |
---|
[845] | 275 | for (list<Polyg *>::const_iterator it = e.is.begin(); it != e.is.end(); ++it) |
---|
[688] | 276 | { |
---|
| 277 | *((int *) &(buffer[0])) += 1; |
---|
| 278 | |
---|
| 279 | *((int *) &(buffer[pos])) = e.id.ind; |
---|
| 280 | pos += sizeof(int); |
---|
| 281 | |
---|
[845] | 282 | *((double *) &(buffer[pos])) = e.area; |
---|
| 283 | pos += sizeof(double); |
---|
| 284 | |
---|
[688] | 285 | *((GloId *) &(buffer[pos])) = (*it)->id; |
---|
| 286 | pos += sizeof(GloId); |
---|
[845] | 287 | |
---|
[688] | 288 | *((int *) &(buffer[pos])) = (*it)->n; |
---|
| 289 | pos += sizeof(int); |
---|
| 290 | *((double *) &(buffer[pos])) = (*it)->area; |
---|
| 291 | pos += sizeof(double); |
---|
| 292 | |
---|
| 293 | *((Coord *) &(buffer[pos])) = (*it)->x; |
---|
| 294 | pos += sizeof(Coord); |
---|
| 295 | } |
---|
| 296 | } |
---|
| 297 | |
---|
| 298 | void unpackIntersection(Elt* e, const char* buffer) |
---|
| 299 | { |
---|
| 300 | int ind; |
---|
| 301 | int pos = 0; |
---|
[845] | 302 | |
---|
[688] | 303 | int n = *((int *) & (buffer[pos])); |
---|
| 304 | pos += sizeof(int); |
---|
| 305 | for (int i = 0; i < n; i++) |
---|
| 306 | { |
---|
| 307 | ind = *((int*) &(buffer[pos])); |
---|
| 308 | pos+=sizeof(int); |
---|
| 309 | |
---|
| 310 | Elt& elt= e[ind]; |
---|
[845] | 311 | |
---|
| 312 | elt.area=*((double *) & (buffer[pos])); |
---|
| 313 | pos += sizeof(double); |
---|
| 314 | |
---|
[688] | 315 | Polyg *polygon = new Polyg; |
---|
| 316 | |
---|
| 317 | polygon->id = *((GloId *) & (buffer[pos])); |
---|
| 318 | pos += sizeof(GloId); |
---|
| 319 | |
---|
| 320 | polygon->n = *((int *) & (buffer[pos])); |
---|
| 321 | pos += sizeof(int); |
---|
| 322 | |
---|
| 323 | polygon->area = *((double *) & (buffer[pos])); |
---|
| 324 | pos += sizeof(double); |
---|
| 325 | |
---|
| 326 | polygon->x = *((Coord *) & (buffer[pos])); |
---|
| 327 | pos += sizeof(Coord); |
---|
| 328 | |
---|
| 329 | elt.is.push_back(polygon); |
---|
| 330 | } |
---|
| 331 | } |
---|
| 332 | |
---|
| 333 | } |
---|