[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 | { |
---|
[1579] | 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]); |
---|
[688] | 29 | |
---|
[1579] | 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 | } |
---|
[688] | 36 | } |
---|
| 37 | |
---|
| 38 | |
---|
| 39 | void normals(Coord *x, int n, Coord *a) |
---|
| 40 | { |
---|
[1579] | 41 | for (int i = 0; i<n; i++) |
---|
| 42 | a[i] = crossprod(x[(i+1)%n], x[i]); |
---|
[688] | 43 | } |
---|
| 44 | |
---|
| 45 | Coord barycentre(const Coord *x, int n) |
---|
| 46 | { |
---|
[1579] | 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 | |
---|
[1579] | 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)); |
---|
[688] | 56 | |
---|
[1579] | 57 | return proj(bc); |
---|
[688] | 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 | { |
---|
[1579] | 64 | Coord n = crossprod(a,b); |
---|
[688] | 65 | double sinc2 = n.x*n.x + n.y*n.y + n.z*n.z; |
---|
[1579] | 66 | assert(sinc2 < 1.0 + EPS); |
---|
[688] | 67 | |
---|
| 68 | // exact: u = asin(sinc)/sinc - 1; asin(sinc) = geodesic length of arc ab |
---|
[1579] | 69 | // approx: |
---|
[688] | 70 | // double u = sinc2/6. + (3./40.)*sinc2*sinc2; |
---|
| 71 | |
---|
[1579] | 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); |
---|
[688] | 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 | { |
---|
[1579] | 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; |
---|
[688] | 89 | } |
---|
| 90 | |
---|
| 91 | Coord exact_barycentre(const Coord *x, int n) |
---|
| 92 | { |
---|
[1579] | 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]; |
---|
[688] | 100 | } |
---|
| 101 | |
---|
| 102 | Coord sc_gc_moon_normalintegral(Coord a, Coord b, Coord pole) |
---|
| 103 | { |
---|
[1579] | 104 | double hemisphere = (a.z > 0) ? 1: -1; |
---|
[688] | 105 | |
---|
[1579] | 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; |
---|
[688] | 110 | |
---|
[1579] | 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; |
---|
[688] | 114 | |
---|
[1579] | 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; |
---|
[688] | 123 | } |
---|
| 124 | |
---|
| 125 | |
---|
[1579] | 126 | double triarea(const Coord& A, const Coord& B, const Coord& C) |
---|
[688] | 127 | { |
---|
[1579] | 128 | double a = ds(B, C); |
---|
| 129 | double b = ds(C, A); |
---|
| 130 | double c = ds(A, B); |
---|
| 131 | double tmp ; |
---|
| 132 | |
---|
| 133 | if (a<b) { tmp=a ; a=b ; b=tmp ; } |
---|
| 134 | if (c > a) { tmp=a ; a=c ; c=b, b=tmp; } |
---|
| 135 | else if (c > b) { tmp=c ; c=b ; b=tmp ; } |
---|
| 136 | |
---|
| 137 | double s = 0.5 * (a + b + c); |
---|
| 138 | 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))); |
---|
| 139 | if (t>0) return 4 * atan(sqrt(t)); |
---|
| 140 | else |
---|
| 141 | { |
---|
| 142 | std::cout<<"double triarea(const Coord& A, const Coord& B, const Coord& C) : t < 0 !!! t="<<t<<endl ; |
---|
| 143 | return 0 ; |
---|
| 144 | } |
---|
[688] | 145 | } |
---|
| 146 | |
---|
| 147 | /** Computes area of two two-sided polygon |
---|
| 148 | needs to have one small and one great circle, otherwise zero |
---|
| 149 | (name origin: lun is moon in french) |
---|
| 150 | */ |
---|
| 151 | double alun(double b, double d) |
---|
| 152 | { |
---|
[1579] | 153 | double a = acos(d); |
---|
| 154 | assert(b <= 2 * a); |
---|
| 155 | double s = a + 0.5 * b; |
---|
| 156 | double t = tan(0.5 * s) * tan(0.5 * (s - a)) * tan(0.5 * (s - a)) * tan(0.5 * (s - b)); |
---|
| 157 | double r = sqrt(1 - d*d); |
---|
| 158 | double p = 2 * asin(sin(0.5*b) / r); |
---|
| 159 | return p*(1 - d) - 4*atan(sqrt(t)); |
---|
[688] | 160 | } |
---|
| 161 | |
---|
| 162 | /** |
---|
| 163 | This function returns the area of a spherical element |
---|
| 164 | that can be composed of great and small circle arcs. |
---|
| 165 | The caller must ensure this function is not called when `alun` should be called instaed. |
---|
| 166 | This function also sets `gg` to the barycentre of the element. |
---|
| 167 | "air" stands for area and "bar" for barycentre. |
---|
| 168 | */ |
---|
| 169 | double airbar(int N, const Coord *x, const Coord *c, double *d, const Coord& pole, Coord& gg) |
---|
| 170 | { |
---|
[1579] | 171 | if (N < 3) |
---|
| 172 | return 0; /* polygons with less then three vertices have zero area */ |
---|
| 173 | Coord t[3]; |
---|
| 174 | t[0] = barycentre(x, N); |
---|
| 175 | Coord *g = new Coord[N]; |
---|
| 176 | double area = 0; |
---|
| 177 | Coord gg_exact = gc_normalintegral(x, N); |
---|
| 178 | for (int i = 0; i < N; i++) |
---|
| 179 | { |
---|
| 180 | /* for "spherical circle segment" sum triangular part and "small moon" and => account for small circle */ |
---|
| 181 | int ii = (i + 1) % N; |
---|
| 182 | t[1] = x[i]; |
---|
| 183 | t[2] = x[ii]; |
---|
[950] | 184 | double sc=scalarprod(crossprod(t[1] - t[0], t[2] - t[0]), t[0]) ; |
---|
[1579] | 185 | assert(sc >= -1e-10); // Error: tri a l'env (wrong orientation) |
---|
| 186 | double area_gc = triarea(t[0], t[1], t[2]); |
---|
| 187 | double area_sc_gc_moon = 0; |
---|
| 188 | if (d[i]) /* handle small circle case */ |
---|
| 189 | { |
---|
| 190 | Coord m = midpoint(t[1], t[2]); |
---|
| 191 | double mext = scalarprod(m, c[i]) - d[i]; |
---|
| 192 | char sgl = (mext > 0) ? -1 : 1; |
---|
| 193 | area_sc_gc_moon = sgl * alun(arcdist(t[1], t[2]), fabs(scalarprod(t[1], pole))); |
---|
| 194 | gg_exact = gg_exact + sc_gc_moon_normalintegral(t[1], t[2], pole); |
---|
| 195 | } |
---|
| 196 | area += area_gc + area_sc_gc_moon; /* for "spherical circle segment" sum triangular part (at) and "small moon" and => account for small circle */ |
---|
| 197 | g[i] = barycentre(t, 3) * (area_gc + area_sc_gc_moon); |
---|
| 198 | } |
---|
| 199 | gg = barycentre(g, N); |
---|
| 200 | gg_exact = proj(gg_exact); |
---|
| 201 | delete[] g; |
---|
| 202 | gg = gg_exact; |
---|
| 203 | return area; |
---|
[688] | 204 | } |
---|
| 205 | |
---|
| 206 | double polygonarea(Coord *vertices, int N) |
---|
| 207 | { |
---|
[1579] | 208 | assert(N >= 3); |
---|
[688] | 209 | |
---|
[1579] | 210 | /* compute polygon area as sum of triangles */ |
---|
| 211 | Coord centre = barycentre(vertices, N); |
---|
| 212 | double area = 0; |
---|
| 213 | for (int i = 0; i < N; i++) |
---|
| 214 | area += triarea(centre, vertices[i], vertices[(i+1)%N]); |
---|
| 215 | return area; |
---|
[688] | 216 | } |
---|
| 217 | |
---|
| 218 | int packedPolygonSize(const Elt& e) |
---|
| 219 | { |
---|
[1614] | 220 | return sizeof(e.id) + sizeof(e.src_id) + sizeof(e.x) + sizeof(e.val) + sizeof(e.given_area)+ |
---|
[1579] | 221 | sizeof(e.n) + e.n*(sizeof(double)+sizeof(Coord)); |
---|
[688] | 222 | } |
---|
| 223 | |
---|
| 224 | void packPolygon(const Elt& e, char *buffer, int& pos) |
---|
| 225 | { |
---|
[1579] | 226 | *((GloId *) &(buffer[pos])) = e.id; |
---|
| 227 | pos += sizeof(e.id); |
---|
| 228 | *((GloId *) &(buffer[pos])) = e.src_id; |
---|
| 229 | pos += sizeof(e.src_id); |
---|
[688] | 230 | |
---|
[1579] | 231 | *((Coord *) &(buffer[pos])) = e.x; |
---|
| 232 | pos += sizeof(e.x); |
---|
[688] | 233 | |
---|
[1579] | 234 | *((double*) &(buffer[pos])) = e.val; |
---|
| 235 | pos += sizeof(e.val); |
---|
[688] | 236 | |
---|
[1614] | 237 | *((double*) &(buffer[pos])) = e.given_area; |
---|
| 238 | pos += sizeof(e.val); |
---|
| 239 | |
---|
[1579] | 240 | *((int *) & (buffer[pos])) = e.n; |
---|
| 241 | pos += sizeof(e.n); |
---|
[688] | 242 | |
---|
[1579] | 243 | for (int i = 0; i < e.n; i++) |
---|
| 244 | { |
---|
| 245 | *((double *) & (buffer[pos])) = e.d[i]; |
---|
| 246 | pos += sizeof(e.d[i]); |
---|
[688] | 247 | |
---|
[1579] | 248 | *((Coord *) &(buffer[pos])) = e.vertex[i]; |
---|
| 249 | pos += sizeof(e.vertex[i]); |
---|
| 250 | } |
---|
[688] | 251 | |
---|
| 252 | } |
---|
| 253 | |
---|
| 254 | void unpackPolygon(Elt& e, const char *buffer, int& pos) |
---|
| 255 | { |
---|
[1579] | 256 | e.id = *((GloId *) &(buffer[pos])); |
---|
| 257 | pos += sizeof(e.id); |
---|
| 258 | e.src_id = *((GloId *) &(buffer[pos])); |
---|
| 259 | pos += sizeof(e.src_id); |
---|
[688] | 260 | |
---|
[1579] | 261 | e.x = *((Coord *) & (buffer[pos]) ); |
---|
| 262 | pos += sizeof(e.x); |
---|
[688] | 263 | |
---|
[1579] | 264 | e.val = *((double *) & (buffer[pos])); |
---|
| 265 | pos += sizeof(double); |
---|
[688] | 266 | |
---|
[1614] | 267 | e.given_area = *((double *) & (buffer[pos])); |
---|
| 268 | pos += sizeof(double); |
---|
| 269 | |
---|
[1579] | 270 | e.n = *((int *) & (buffer[pos])); |
---|
| 271 | pos += sizeof(int); |
---|
[688] | 272 | |
---|
[1579] | 273 | for (int i = 0; i < e.n; i++) |
---|
| 274 | { |
---|
| 275 | e.d[i] = *((double *) & (buffer[pos])); |
---|
| 276 | pos += sizeof(double); |
---|
[688] | 277 | |
---|
[1579] | 278 | e.vertex[i] = *((Coord *) & (buffer[pos])); |
---|
| 279 | pos += sizeof(Coord); |
---|
| 280 | } |
---|
[688] | 281 | } |
---|
| 282 | |
---|
| 283 | int packIntersectionSize(const Elt& elt) |
---|
| 284 | { |
---|
[1579] | 285 | return elt.is.size() * (2*sizeof(int)+ sizeof(GloId) + 5*sizeof(double)); |
---|
[688] | 286 | } |
---|
| 287 | |
---|
| 288 | void packIntersection(const Elt& e, char* buffer,int& pos) |
---|
| 289 | { |
---|
[845] | 290 | for (list<Polyg *>::const_iterator it = e.is.begin(); it != e.is.end(); ++it) |
---|
[1579] | 291 | { |
---|
| 292 | *((int *) &(buffer[0])) += 1; |
---|
[688] | 293 | |
---|
[1579] | 294 | *((int *) &(buffer[pos])) = e.id.ind; |
---|
| 295 | pos += sizeof(int); |
---|
[688] | 296 | |
---|
[845] | 297 | *((double *) &(buffer[pos])) = e.area; |
---|
| 298 | pos += sizeof(double); |
---|
| 299 | |
---|
[1614] | 300 | |
---|
[1579] | 301 | *((GloId *) &(buffer[pos])) = (*it)->id; |
---|
| 302 | pos += sizeof(GloId); |
---|
[845] | 303 | |
---|
[1579] | 304 | *((int *) &(buffer[pos])) = (*it)->n; |
---|
| 305 | pos += sizeof(int); |
---|
| 306 | *((double *) &(buffer[pos])) = (*it)->area; |
---|
| 307 | pos += sizeof(double); |
---|
[688] | 308 | |
---|
[1579] | 309 | *((Coord *) &(buffer[pos])) = (*it)->x; |
---|
| 310 | pos += sizeof(Coord); |
---|
| 311 | } |
---|
[688] | 312 | } |
---|
| 313 | |
---|
| 314 | void unpackIntersection(Elt* e, const char* buffer) |
---|
| 315 | { |
---|
[1579] | 316 | int ind; |
---|
| 317 | int pos = 0; |
---|
[845] | 318 | |
---|
[1579] | 319 | int n = *((int *) & (buffer[pos])); |
---|
| 320 | pos += sizeof(int); |
---|
| 321 | for (int i = 0; i < n; i++) |
---|
| 322 | { |
---|
| 323 | ind = *((int*) &(buffer[pos])); |
---|
| 324 | pos+=sizeof(int); |
---|
[688] | 325 | |
---|
[1579] | 326 | Elt& elt= e[ind]; |
---|
[845] | 327 | |
---|
| 328 | elt.area=*((double *) & (buffer[pos])); |
---|
[1579] | 329 | pos += sizeof(double); |
---|
[845] | 330 | |
---|
[1614] | 331 | |
---|
[1579] | 332 | Polyg *polygon = new Polyg; |
---|
[688] | 333 | |
---|
[1579] | 334 | polygon->id = *((GloId *) & (buffer[pos])); |
---|
| 335 | pos += sizeof(GloId); |
---|
[688] | 336 | |
---|
[1579] | 337 | polygon->n = *((int *) & (buffer[pos])); |
---|
| 338 | pos += sizeof(int); |
---|
[688] | 339 | |
---|
[1579] | 340 | polygon->area = *((double *) & (buffer[pos])); |
---|
| 341 | pos += sizeof(double); |
---|
[688] | 342 | |
---|
[1579] | 343 | polygon->x = *((Coord *) & (buffer[pos])); |
---|
| 344 | pos += sizeof(Coord); |
---|
[688] | 345 | |
---|
[1579] | 346 | elt.is.push_back(polygon); |
---|
| 347 | } |
---|
[688] | 348 | } |
---|
| 349 | |
---|
| 350 | } |
---|