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