[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 | */ |
---|
[2500] | 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 |
---|
[2500] | 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 | { |
---|
[2500] | 37 | swap(edge[N-2-i], edge[i]); |
---|
| 38 | swap(d[i], d[N-2-i]); |
---|
[1579] | 39 | } |
---|
[2500] | 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)); |
---|
[2500] | 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]; |
---|
[688] | 105 | } |
---|
| 106 | |
---|
[2500] | 107 | /* other methode to compute barycenter of spherical polygon |
---|
| 108 | for a spherical polygon, the moment is half the sum of (a x b) / ||a x b|| * (angle between a and b) |
---|
| 109 | for each pair of consecutive vertices a,b. |
---|
| 110 | */ |
---|
| 111 | |
---|
| 112 | Coord new_barycentre(const Coord *x, int n) |
---|
| 113 | { |
---|
| 114 | if (n >= 3) |
---|
| 115 | { |
---|
| 116 | Coord sum=ORIGIN ; |
---|
| 117 | for (int i = 0; i < n; i++) |
---|
| 118 | { |
---|
| 119 | 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]))) ; |
---|
| 120 | } |
---|
| 121 | return proj(sum) ; |
---|
| 122 | } |
---|
| 123 | else if (n == 0) return ORIGIN; |
---|
| 124 | else if (n == 2) return midpoint(x[0], x[1]); |
---|
| 125 | else if (n == 1) return x[0]; |
---|
| 126 | } |
---|
| 127 | |
---|
| 128 | |
---|
| 129 | |
---|
[688] | 130 | Coord sc_gc_moon_normalintegral(Coord a, Coord b, Coord pole) |
---|
| 131 | { |
---|
[1579] | 132 | double hemisphere = (a.z > 0) ? 1: -1; |
---|
[688] | 133 | |
---|
[1579] | 134 | double lat = hemisphere * (M_PI_2 - acos(a.z)); |
---|
| 135 | double lon1 = atan2(a.y, a.x); |
---|
| 136 | double lon2 = atan2(b.y, b.x); |
---|
| 137 | double lon_diff = lon2 - lon1; |
---|
[688] | 138 | |
---|
[1579] | 139 | // wraparound at lon=-pi=pi |
---|
| 140 | if (lon_diff < -M_PI) lon_diff += 2.0*M_PI; |
---|
| 141 | else if (lon_diff > M_PI) lon_diff -= 2.0*M_PI; |
---|
[688] | 142 | |
---|
[1579] | 143 | // integral of the normal over the surface bound by great arcs a-pole and b-pole and small arc a-b |
---|
| 144 | Coord sc_normalintegral = Coord(0.5*(sin(lon2)-sin(lon1))*(M_PI_2 - lat - 0.5*sin(2.0*lat)), |
---|
| 145 | 0.5*(cos(lon1)-cos(lon2))*(M_PI_2 - lat - 0.5*sin(2.0*lat)), |
---|
| 146 | hemisphere * lon_diff * 0.25 * (cos(2.0*lat) + 1.0)); |
---|
| 147 | Coord p = Coord(0,0,hemisphere); // TODO assumes north pole is (0,0,1) |
---|
| 148 | Coord t[] = {a, b, p}; |
---|
| 149 | if (hemisphere < 0) swap(t[0], t[1]); |
---|
| 150 | return (sc_normalintegral - gc_normalintegral(t, 3)) * hemisphere; |
---|
[688] | 151 | } |
---|
| 152 | |
---|
| 153 | |
---|
[1579] | 154 | double triarea(const Coord& A, const Coord& B, const Coord& C) |
---|
[688] | 155 | { |
---|
[1579] | 156 | double a = ds(B, C); |
---|
| 157 | double b = ds(C, A); |
---|
| 158 | double c = ds(A, B); |
---|
| 159 | double tmp ; |
---|
| 160 | |
---|
| 161 | if (a<b) { tmp=a ; a=b ; b=tmp ; } |
---|
| 162 | if (c > a) { tmp=a ; a=c ; c=b, b=tmp; } |
---|
| 163 | else if (c > b) { tmp=c ; c=b ; b=tmp ; } |
---|
| 164 | |
---|
| 165 | double s = 0.5 * (a + b + c); |
---|
| 166 | 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))); |
---|
| 167 | if (t>0) return 4 * atan(sqrt(t)); |
---|
| 168 | else |
---|
| 169 | { |
---|
| 170 | std::cout<<"double triarea(const Coord& A, const Coord& B, const Coord& C) : t < 0 !!! t="<<t<<endl ; |
---|
| 171 | return 0 ; |
---|
| 172 | } |
---|
[688] | 173 | } |
---|
| 174 | |
---|
| 175 | /** Computes area of two two-sided polygon |
---|
| 176 | needs to have one small and one great circle, otherwise zero |
---|
| 177 | (name origin: lun is moon in french) |
---|
| 178 | */ |
---|
| 179 | double alun(double b, double d) |
---|
| 180 | { |
---|
[1579] | 181 | double a = acos(d); |
---|
| 182 | assert(b <= 2 * a); |
---|
| 183 | double s = a + 0.5 * b; |
---|
| 184 | double t = tan(0.5 * s) * tan(0.5 * (s - a)) * tan(0.5 * (s - a)) * tan(0.5 * (s - b)); |
---|
| 185 | double r = sqrt(1 - d*d); |
---|
| 186 | double p = 2 * asin(sin(0.5*b) / r); |
---|
| 187 | return p*(1 - d) - 4*atan(sqrt(t)); |
---|
[688] | 188 | } |
---|
| 189 | |
---|
| 190 | /** |
---|
| 191 | This function returns the area of a spherical element |
---|
| 192 | that can be composed of great and small circle arcs. |
---|
| 193 | The caller must ensure this function is not called when `alun` should be called instaed. |
---|
| 194 | This function also sets `gg` to the barycentre of the element. |
---|
| 195 | "air" stands for area and "bar" for barycentre. |
---|
| 196 | */ |
---|
| 197 | double airbar(int N, const Coord *x, const Coord *c, double *d, const Coord& pole, Coord& gg) |
---|
| 198 | { |
---|
[1579] | 199 | if (N < 3) |
---|
| 200 | return 0; /* polygons with less then three vertices have zero area */ |
---|
| 201 | Coord t[3]; |
---|
| 202 | t[0] = barycentre(x, N); |
---|
| 203 | Coord *g = new Coord[N]; |
---|
| 204 | double area = 0; |
---|
| 205 | Coord gg_exact = gc_normalintegral(x, N); |
---|
| 206 | for (int i = 0; i < N; i++) |
---|
| 207 | { |
---|
| 208 | /* for "spherical circle segment" sum triangular part and "small moon" and => account for small circle */ |
---|
| 209 | int ii = (i + 1) % N; |
---|
| 210 | t[1] = x[i]; |
---|
| 211 | t[2] = x[ii]; |
---|
[950] | 212 | double sc=scalarprod(crossprod(t[1] - t[0], t[2] - t[0]), t[0]) ; |
---|
[1579] | 213 | assert(sc >= -1e-10); // Error: tri a l'env (wrong orientation) |
---|
| 214 | double area_gc = triarea(t[0], t[1], t[2]); |
---|
| 215 | double area_sc_gc_moon = 0; |
---|
| 216 | if (d[i]) /* handle small circle case */ |
---|
| 217 | { |
---|
| 218 | Coord m = midpoint(t[1], t[2]); |
---|
| 219 | double mext = scalarprod(m, c[i]) - d[i]; |
---|
| 220 | char sgl = (mext > 0) ? -1 : 1; |
---|
| 221 | area_sc_gc_moon = sgl * alun(arcdist(t[1], t[2]), fabs(scalarprod(t[1], pole))); |
---|
| 222 | gg_exact = gg_exact + sc_gc_moon_normalintegral(t[1], t[2], pole); |
---|
| 223 | } |
---|
| 224 | area += area_gc + area_sc_gc_moon; /* for "spherical circle segment" sum triangular part (at) and "small moon" and => account for small circle */ |
---|
| 225 | g[i] = barycentre(t, 3) * (area_gc + area_sc_gc_moon); |
---|
| 226 | } |
---|
| 227 | gg = barycentre(g, N); |
---|
| 228 | gg_exact = proj(gg_exact); |
---|
| 229 | delete[] g; |
---|
| 230 | gg = gg_exact; |
---|
| 231 | return area; |
---|
[688] | 232 | } |
---|
| 233 | |
---|
| 234 | double polygonarea(Coord *vertices, int N) |
---|
| 235 | { |
---|
[1579] | 236 | assert(N >= 3); |
---|
[688] | 237 | |
---|
[1579] | 238 | /* compute polygon area as sum of triangles */ |
---|
| 239 | Coord centre = barycentre(vertices, N); |
---|
| 240 | double area = 0; |
---|
| 241 | for (int i = 0; i < N; i++) |
---|
| 242 | area += triarea(centre, vertices[i], vertices[(i+1)%N]); |
---|
| 243 | return area; |
---|
[688] | 244 | } |
---|
| 245 | |
---|
| 246 | int packedPolygonSize(const Elt& e) |
---|
| 247 | { |
---|
[1614] | 248 | return sizeof(e.id) + sizeof(e.src_id) + sizeof(e.x) + sizeof(e.val) + sizeof(e.given_area)+ |
---|
[1579] | 249 | sizeof(e.n) + e.n*(sizeof(double)+sizeof(Coord)); |
---|
[688] | 250 | } |
---|
| 251 | |
---|
| 252 | void packPolygon(const Elt& e, char *buffer, int& pos) |
---|
| 253 | { |
---|
[1579] | 254 | *((GloId *) &(buffer[pos])) = e.id; |
---|
| 255 | pos += sizeof(e.id); |
---|
| 256 | *((GloId *) &(buffer[pos])) = e.src_id; |
---|
| 257 | pos += sizeof(e.src_id); |
---|
[688] | 258 | |
---|
[1579] | 259 | *((Coord *) &(buffer[pos])) = e.x; |
---|
| 260 | pos += sizeof(e.x); |
---|
[688] | 261 | |
---|
[1579] | 262 | *((double*) &(buffer[pos])) = e.val; |
---|
| 263 | pos += sizeof(e.val); |
---|
[688] | 264 | |
---|
[1614] | 265 | *((double*) &(buffer[pos])) = e.given_area; |
---|
| 266 | pos += sizeof(e.val); |
---|
| 267 | |
---|
[1579] | 268 | *((int *) & (buffer[pos])) = e.n; |
---|
| 269 | pos += sizeof(e.n); |
---|
[688] | 270 | |
---|
[1579] | 271 | for (int i = 0; i < e.n; i++) |
---|
| 272 | { |
---|
| 273 | *((double *) & (buffer[pos])) = e.d[i]; |
---|
| 274 | pos += sizeof(e.d[i]); |
---|
[688] | 275 | |
---|
[1579] | 276 | *((Coord *) &(buffer[pos])) = e.vertex[i]; |
---|
| 277 | pos += sizeof(e.vertex[i]); |
---|
| 278 | } |
---|
[688] | 279 | |
---|
| 280 | } |
---|
| 281 | |
---|
| 282 | void unpackPolygon(Elt& e, const char *buffer, int& pos) |
---|
| 283 | { |
---|
[1579] | 284 | e.id = *((GloId *) &(buffer[pos])); |
---|
| 285 | pos += sizeof(e.id); |
---|
| 286 | e.src_id = *((GloId *) &(buffer[pos])); |
---|
| 287 | pos += sizeof(e.src_id); |
---|
[688] | 288 | |
---|
[1579] | 289 | e.x = *((Coord *) & (buffer[pos]) ); |
---|
| 290 | pos += sizeof(e.x); |
---|
[688] | 291 | |
---|
[1579] | 292 | e.val = *((double *) & (buffer[pos])); |
---|
| 293 | pos += sizeof(double); |
---|
[688] | 294 | |
---|
[1614] | 295 | e.given_area = *((double *) & (buffer[pos])); |
---|
| 296 | pos += sizeof(double); |
---|
| 297 | |
---|
[1579] | 298 | e.n = *((int *) & (buffer[pos])); |
---|
| 299 | pos += sizeof(int); |
---|
[688] | 300 | |
---|
[1579] | 301 | for (int i = 0; i < e.n; i++) |
---|
| 302 | { |
---|
| 303 | e.d[i] = *((double *) & (buffer[pos])); |
---|
| 304 | pos += sizeof(double); |
---|
[688] | 305 | |
---|
[1579] | 306 | e.vertex[i] = *((Coord *) & (buffer[pos])); |
---|
| 307 | pos += sizeof(Coord); |
---|
| 308 | } |
---|
[688] | 309 | } |
---|
| 310 | |
---|
| 311 | int packIntersectionSize(const Elt& elt) |
---|
| 312 | { |
---|
[1579] | 313 | return elt.is.size() * (2*sizeof(int)+ sizeof(GloId) + 5*sizeof(double)); |
---|
[688] | 314 | } |
---|
| 315 | |
---|
| 316 | void packIntersection(const Elt& e, char* buffer,int& pos) |
---|
| 317 | { |
---|
[845] | 318 | for (list<Polyg *>::const_iterator it = e.is.begin(); it != e.is.end(); ++it) |
---|
[1579] | 319 | { |
---|
| 320 | *((int *) &(buffer[0])) += 1; |
---|
[688] | 321 | |
---|
[1579] | 322 | *((int *) &(buffer[pos])) = e.id.ind; |
---|
| 323 | pos += sizeof(int); |
---|
[688] | 324 | |
---|
[845] | 325 | *((double *) &(buffer[pos])) = e.area; |
---|
| 326 | pos += sizeof(double); |
---|
| 327 | |
---|
[1614] | 328 | |
---|
[1579] | 329 | *((GloId *) &(buffer[pos])) = (*it)->id; |
---|
| 330 | pos += sizeof(GloId); |
---|
[845] | 331 | |
---|
[1579] | 332 | *((int *) &(buffer[pos])) = (*it)->n; |
---|
| 333 | pos += sizeof(int); |
---|
| 334 | *((double *) &(buffer[pos])) = (*it)->area; |
---|
| 335 | pos += sizeof(double); |
---|
[688] | 336 | |
---|
[1579] | 337 | *((Coord *) &(buffer[pos])) = (*it)->x; |
---|
| 338 | pos += sizeof(Coord); |
---|
| 339 | } |
---|
[688] | 340 | } |
---|
| 341 | |
---|
| 342 | void unpackIntersection(Elt* e, const char* buffer) |
---|
| 343 | { |
---|
[1579] | 344 | int ind; |
---|
| 345 | int pos = 0; |
---|
[845] | 346 | |
---|
[1579] | 347 | int n = *((int *) & (buffer[pos])); |
---|
| 348 | pos += sizeof(int); |
---|
| 349 | for (int i = 0; i < n; i++) |
---|
| 350 | { |
---|
| 351 | ind = *((int*) &(buffer[pos])); |
---|
| 352 | pos+=sizeof(int); |
---|
[688] | 353 | |
---|
[1579] | 354 | Elt& elt= e[ind]; |
---|
[845] | 355 | |
---|
| 356 | elt.area=*((double *) & (buffer[pos])); |
---|
[1579] | 357 | pos += sizeof(double); |
---|
[845] | 358 | |
---|
[1614] | 359 | |
---|
[1579] | 360 | Polyg *polygon = new Polyg; |
---|
[688] | 361 | |
---|
[1579] | 362 | polygon->id = *((GloId *) & (buffer[pos])); |
---|
| 363 | pos += sizeof(GloId); |
---|
[688] | 364 | |
---|
[1579] | 365 | polygon->n = *((int *) & (buffer[pos])); |
---|
| 366 | pos += sizeof(int); |
---|
[688] | 367 | |
---|
[1579] | 368 | polygon->area = *((double *) & (buffer[pos])); |
---|
| 369 | pos += sizeof(double); |
---|
[688] | 370 | |
---|
[1579] | 371 | polygon->x = *((Coord *) & (buffer[pos])); |
---|
| 372 | pos += sizeof(Coord); |
---|
[688] | 373 | |
---|
[1579] | 374 | elt.is.push_back(polygon); |
---|
| 375 | } |
---|
[688] | 376 | } |
---|
| 377 | |
---|
| 378 | } |
---|