PostGIS  2.1.10dev-r@@SVN_REVISION@@
lwgeodetic.c
Go to the documentation of this file.
1 /**********************************************************************
2  * $Id: lwgeodetic.c 14487 2015-12-14 13:10:24Z dbaston $
3  *
4  * PostGIS - Spatial Types for PostgreSQL
5  * Copyright 2009 Paul Ramsey <pramsey@cleverelephant.ca>
6  * Copyright 2009 David Skea <David.Skea@gov.bc.ca>
7  *
8  * This is free software; you can redistribute and/or modify it under
9  * the terms of the GNU General Public Licence. See the COPYING file.
10  *
11  **********************************************************************/
12 
13 #include "liblwgeom_internal.h"
14 #include "lwgeodetic.h"
15 #include "lwgeom_log.h"
16 
23 
27 double longitude_radians_normalize(double lon)
28 {
29  if ( lon == -1.0 * M_PI )
30  return M_PI;
31  if ( lon == -2.0 * M_PI )
32  return 0.0;
33 
34  if ( lon > 2.0 * M_PI )
35  lon = remainder(lon, 2.0 * M_PI);
36 
37  if ( lon < -2.0 * M_PI )
38  lon = remainder(lon, -2.0 * M_PI);
39 
40  if ( lon > M_PI )
41  lon = -2.0 * M_PI + lon;
42 
43  if ( lon < -1.0 * M_PI )
44  lon = 2.0 * M_PI + lon;
45 
46  if ( lon == -2.0 * M_PI )
47  lon *= -1.0;
48 
49  return lon;
50 }
51 
55 double latitude_radians_normalize(double lat)
56 {
57 
58  if ( lat > 2.0 * M_PI )
59  lat = remainder(lat, 2.0 * M_PI);
60 
61  if ( lat < -2.0 * M_PI )
62  lat = remainder(lat, -2.0 * M_PI);
63 
64  if ( lat > M_PI )
65  lat = M_PI - lat;
66 
67  if ( lat < -1.0 * M_PI )
68  lat = -1.0 * M_PI - lat;
69 
70  if ( lat > M_PI_2 )
71  lat = M_PI - lat;
72 
73  if ( lat < -1.0 * M_PI_2 )
74  lat = -1.0 * M_PI - lat;
75 
76  return lat;
77 }
78 
83 double longitude_degrees_normalize(double lon)
84 {
85  if ( lon > 360.0 )
86  lon = remainder(lon, 360.0);
87 
88  if ( lon < -360.0 )
89  lon = remainder(lon, -360.0);
90 
91  if ( lon > 180.0 )
92  lon = -360.0 + lon;
93 
94  if ( lon < -180.0 )
95  lon = 360 + lon;
96 
97  if ( lon == -180.0 )
98  return 180.0;
99 
100  if ( lon == -360.0 )
101  return 0.0;
102 
103  return lon;
104 }
105 
110 double latitude_degrees_normalize(double lat)
111 {
112 
113  if ( lat > 360.0 )
114  lat = remainder(lat, 360.0);
115 
116  if ( lat < -360.0 )
117  lat = remainder(lat, -360.0);
118 
119  if ( lat > 180.0 )
120  lat = 180.0 - lat;
121 
122  if ( lat < -180.0 )
123  lat = -180.0 - lat;
124 
125  if ( lat > 90.0 )
126  lat = 180.0 - lat;
127 
128  if ( lat < -90.0 )
129  lat = -180.0 - lat;
130 
131  return lat;
132 }
133 
137 void point_shift(GEOGRAPHIC_POINT *p, double shift)
138 {
139  double lon = p->lon + shift;
140  if ( lon > M_PI )
141  p->lon = -1.0 * M_PI + (lon - M_PI);
142  else
143  p->lon = lon;
144  return;
145 }
146 
148 {
149  return FP_EQUALS(g1->lat, g2->lat) && FP_EQUALS(g1->lon, g2->lon);
150 }
151 
157 void geographic_point_init(double lon, double lat, GEOGRAPHIC_POINT *g)
158 {
161 }
162 
164 double
166 {
167  double d[6];
168  int i;
169  double zmin = MAXFLOAT;
170  double zmax = -1 * MAXFLOAT;
171  POINT3D pt;
172 
173  /* Take a copy of the box corners so we can treat them as a list */
174  /* Elements are xmin, xmax, ymin, ymax, zmin, zmax */
175  memcpy(d, &(gbox->xmin), 6*sizeof(double));
176 
177  /* Generate all 8 corner vectors of the box */
178  for ( i = 0; i < 8; i++ )
179  {
180  pt.x = d[i / 4];
181  pt.y = d[2 + (i % 4) / 2];
182  pt.z = d[4 + (i % 2)];
183  normalize(&pt);
184  if ( pt.z < zmin ) zmin = pt.z;
185  if ( pt.z > zmax ) zmax = pt.z;
186  }
187  return asin(zmax) - asin(zmin);
188 }
189 
191 double
193 {
194  double d[6];
195  int i, j;
196  POINT3D pt[3];
197  double maxangle;
198  double magnitude;
199 
200  /* Take a copy of the box corners so we can treat them as a list */
201  /* Elements are xmin, xmax, ymin, ymax, zmin, zmax */
202  memcpy(d, &(gbox->xmin), 6*sizeof(double));
203 
204  /* Start with the bottom corner */
205  pt[0].x = gbox->xmin;
206  pt[0].y = gbox->ymin;
207  magnitude = sqrt(pt[0].x*pt[0].x + pt[0].y*pt[0].y);
208  pt[0].x /= magnitude;
209  pt[0].y /= magnitude;
210 
211  /* Generate all 8 corner vectors of the box */
212  /* Find the vector furthest from our seed vector */
213  for ( j = 0; j < 2; j++ )
214  {
215  maxangle = -1 * MAXFLOAT;
216  for ( i = 0; i < 4; i++ )
217  {
218  double angle, dotprod;
219  POINT3D pt_n;
220 
221  pt_n.x = d[i / 2];
222  pt_n.y = d[2 + (i % 2)];
223  magnitude = sqrt(pt_n.x*pt_n.x + pt_n.y*pt_n.y);
224  pt_n.x /= magnitude;
225  pt_n.y /= magnitude;
226  pt_n.z = 0.0;
227 
228  dotprod = pt_n.x*pt[j].x + pt_n.y*pt[j].y;
229  angle = acos(dotprod > 1.0 ? 1.0 : dotprod);
230  if ( angle > maxangle )
231  {
232  pt[j+1] = pt_n;
233  maxangle = angle;
234  }
235  }
236  }
237 
238  /* Return the distance between the two furthest vectors */
239  return maxangle;
240 }
241 
243 int
244 gbox_centroid(const GBOX* gbox, POINT2D* out)
245 {
246  double d[6];
248  POINT3D pt;
249  int i;
250 
251  /* Take a copy of the box corners so we can treat them as a list */
252  /* Elements are xmin, xmax, ymin, ymax, zmin, zmax */
253  memcpy(d, &(gbox->xmin), 6*sizeof(double));
254 
255  /* Zero out our return vector */
256  pt.x = pt.y = pt.z = 0.0;
257 
258  for ( i = 0; i < 8; i++ )
259  {
260  POINT3D pt_n;
261 
262  pt_n.x = d[i / 4];
263  pt_n.y = d[2 + ((i % 4) / 2)];
264  pt_n.z = d[4 + (i % 2)];
265  normalize(&pt_n);
266 
267  pt.x += pt_n.x;
268  pt.y += pt_n.y;
269  pt.z += pt_n.z;
270  }
271 
272  pt.x /= 8.0;
273  pt.y /= 8.0;
274  pt.z /= 8.0;
275  normalize(&pt);
276 
277  cart2geog(&pt, &g);
280 
281  return LW_SUCCESS;
282 }
283 
293 static int gbox_check_poles(GBOX *gbox)
294 {
295  int rv = LW_FALSE;
296  LWDEBUG(4, "checking poles");
297  LWDEBUGF(4, "gbox %s", gbox_to_string(gbox));
298  /* Z axis */
299  if ( gbox->xmin < 0.0 && gbox->xmax > 0.0 &&
300  gbox->ymin < 0.0 && gbox->ymax > 0.0 )
301  {
302  if ( (gbox->zmin + gbox->zmax) > 0.0 )
303  {
304  LWDEBUG(4, "enclosed positive z axis");
305  gbox->zmax = 1.0;
306  }
307  else
308  {
309  LWDEBUG(4, "enclosed negative z axis");
310  gbox->zmin = -1.0;
311  }
312  rv = LW_TRUE;
313  }
314 
315  /* Y axis */
316  if ( gbox->xmin < 0.0 && gbox->xmax > 0.0 &&
317  gbox->zmin < 0.0 && gbox->zmax > 0.0 )
318  {
319  if ( gbox->ymin + gbox->ymax > 0.0 )
320  {
321  LWDEBUG(4, "enclosed positive y axis");
322  gbox->ymax = 1.0;
323  }
324  else
325  {
326  LWDEBUG(4, "enclosed negative y axis");
327  gbox->ymin = -1.0;
328  }
329  rv = LW_TRUE;
330  }
331 
332  /* X axis */
333  if ( gbox->ymin < 0.0 && gbox->ymax > 0.0 &&
334  gbox->zmin < 0.0 && gbox->zmax > 0.0 )
335  {
336  if ( gbox->xmin + gbox->xmax > 0.0 )
337  {
338  LWDEBUG(4, "enclosed positive x axis");
339  gbox->xmax = 1.0;
340  }
341  else
342  {
343  LWDEBUG(4, "enclosed negative x axis");
344  gbox->xmin = -1.0;
345  }
346  rv = LW_TRUE;
347  }
348 
349  return rv;
350 }
351 
356 {
357  p->x = cos(g->lat) * cos(g->lon);
358  p->y = cos(g->lat) * sin(g->lon);
359  p->z = sin(g->lat);
360 }
361 
366 {
367  g->lon = atan2(p->y, p->x);
368  g->lat = asin(p->z);
369 }
370 
374 void ll2cart(const POINT2D *g, POINT3D *p)
375 {
376  double x_rad = M_PI * g->x / 180.0;
377  double y_rad = M_PI * g->y / 180.0;
378  double cos_y_rad = cos(y_rad);
379  p->x = cos_y_rad * cos(x_rad);
380  p->y = cos_y_rad * sin(x_rad);
381  p->z = sin(y_rad);
382 }
383 
397 static double dot_product(const POINT3D *p1, const POINT3D *p2)
398 {
399  return (p1->x*p2->x) + (p1->y*p2->y) + (p1->z*p2->z);
400 }
401 
405 static void cross_product(const POINT3D *a, const POINT3D *b, POINT3D *n)
406 {
407  n->x = a->y * b->z - a->z * b->y;
408  n->y = a->z * b->x - a->x * b->z;
409  n->z = a->x * b->y - a->y * b->x;
410  return;
411 }
412 
416 void vector_sum(const POINT3D *a, const POINT3D *b, POINT3D *n)
417 {
418  n->x = a->x + b->x;
419  n->y = a->y + b->y;
420  n->z = a->z + b->z;
421  return;
422 }
423 
427 static void vector_difference(const POINT3D *a, const POINT3D *b, POINT3D *n)
428 {
429  n->x = a->x - b->x;
430  n->y = a->y - b->y;
431  n->z = a->z - b->z;
432  return;
433 }
434 
438 static void vector_scale(POINT3D *n, double scale)
439 {
440  n->x *= scale;
441  n->y *= scale;
442  n->z *= scale;
443  return;
444 }
445 
446 // static inline double vector_magnitude(const POINT3D* v)
447 // {
448 // return sqrt(v->x*v->x + v->y*v->y + v->z*v->z);
449 // }
450 
454 double vector_angle(const POINT3D* v1, const POINT3D* v2)
455 {
456  POINT3D v3, normal;
457  double angle, x, y;
458 
459  cross_product(v1, v2, &normal);
460  normalize(&normal);
461  cross_product(&normal, v1, &v3);
462 
463  x = dot_product(v1, v2);
464  y = dot_product(v2, &v3);
465 
466  angle = atan2(y, x);
467  return angle;
468 }
469 
473 static void normalize2d(POINT2D *p)
474 {
475  double d = sqrt(p->x*p->x + p->y*p->y);
476  if (FP_IS_ZERO(d))
477  {
478  p->x = p->y = 0.0;
479  return;
480  }
481  p->x = p->x / d;
482  p->y = p->y / d;
483  return;
484 }
485 
490 void unit_normal(const POINT3D *P1, const POINT3D *P2, POINT3D *normal)
491 {
492  double p_dot = dot_product(P1, P2);
493  POINT3D P3;
494 
495  /* If edge is really large, calculate a narrower equivalent angle A1/A3. */
496  if ( p_dot < 0 )
497  {
498  vector_sum(P1, P2, &P3);
499  normalize(&P3);
500  }
501  /* If edge is narrow, calculate a wider equivalent angle A1/A3. */
502  else if ( p_dot > 0.95 )
503  {
504  vector_difference(P2, P1, &P3);
505  normalize(&P3);
506  }
507  /* Just keep the current angle in A1/A3. */
508  else
509  {
510  P3 = *P2;
511  }
512 
513  /* Normals to the A-plane and B-plane */
514  cross_product(P1, &P3, normal);
515  normalize(normal);
516 }
517 
522 void vector_rotate(const POINT3D* v1, const POINT3D* v2, double angle, POINT3D* n)
523 {
524  POINT3D u;
525  double cos_a = cos(angle);
526  double sin_a = sin(angle);
527  double uxuy, uyuz, uxuz;
528  double ux2, uy2, uz2;
529  double rxx, rxy, rxz, ryx, ryy, ryz, rzx, rzy, rzz;
530 
531  /* Need a unit vector normal to rotate around */
532  unit_normal(v1, v2, &u);
533 
534  uxuy = u.x * u.y;
535  uxuz = u.x * u.z;
536  uyuz = u.y * u.z;
537 
538  ux2 = u.x * u.x;
539  uy2 = u.y * u.y;
540  uz2 = u.z * u.z;
541 
542  rxx = cos_a + ux2 * (1 - cos_a);
543  rxy = uxuy * (1 - cos_a) - u.z * sin_a;
544  rxz = uxuz * (1 - cos_a) + u.y * sin_a;
545 
546  ryx = uxuy * (1 - cos_a) + u.z * sin_a;
547  ryy = cos_a + uy2 * (1 - cos_a);
548  ryz = uyuz * (1 - cos_a) - u.x * sin_a;
549 
550  rzx = uxuz * (1 - cos_a) - u.y * sin_a;
551  rzy = uyuz * (1 - cos_a) + u.x * sin_a;
552  rzz = cos_a + uz2 * (1 - cos_a);
553 
554  n->x = rxx * v1->x + rxy * v1->y + rxz * v1->z;
555  n->y = ryx * v1->x + ryy * v1->y + ryz * v1->z;
556  n->z = rzx * v1->x + rzy * v1->y + rzz * v1->z;
557 
558  normalize(n);
559 }
560 
565 {
566  double d = sqrt(p->x*p->x + p->y*p->y + p->z*p->z);
567  if (FP_IS_ZERO(d))
568  {
569  p->x = p->y = p->z = 0.0;
570  return;
571  }
572  p->x = p->x / d;
573  p->y = p->y / d;
574  p->z = p->z / d;
575  return;
576 }
577 
578 
584 {
585  double lon_qpp = (q->lon + p->lon) / -2.0;
586  double lon_qmp = (q->lon - p->lon) / 2.0;
587  double sin_p_lat_minus_q_lat = sin(p->lat-q->lat);
588  double sin_p_lat_plus_q_lat = sin(p->lat+q->lat);
589  double sin_lon_qpp = sin(lon_qpp);
590  double sin_lon_qmp = sin(lon_qmp);
591  double cos_lon_qpp = cos(lon_qpp);
592  double cos_lon_qmp = cos(lon_qmp);
593  a->x = sin_p_lat_minus_q_lat * sin_lon_qpp * cos_lon_qmp -
594  sin_p_lat_plus_q_lat * cos_lon_qpp * sin_lon_qmp;
595  a->y = sin_p_lat_minus_q_lat * cos_lon_qpp * cos_lon_qmp +
596  sin_p_lat_plus_q_lat * sin_lon_qpp * sin_lon_qmp;
597  a->z = cos(p->lat) * cos(q->lat) * sin(q->lon-p->lon);
598 }
599 
600 void x_to_z(POINT3D *p)
601 {
602  double tmp = p->z;
603  p->z = p->x;
604  p->x = tmp;
605 }
606 
607 void y_to_z(POINT3D *p)
608 {
609  double tmp = p->z;
610  p->z = p->y;
611  p->y = tmp;
612 }
613 
614 
616 {
617  double sign_s = signum(s->lon);
618  double sign_e = signum(e->lon);
619  double ss = fabs(s->lon);
620  double ee = fabs(e->lon);
621  if ( sign_s == sign_e )
622  {
623  return LW_FALSE;
624  }
625  else
626  {
627  double dl = ss + ee;
628  if ( dl < M_PI )
629  return LW_FALSE;
630  else if ( FP_EQUALS(dl, M_PI) )
631  return LW_FALSE;
632  else
633  return LW_TRUE;
634  }
635 }
636 
642 static int
644 {
645  POINT3D normal, pt;
646  double w;
647  /* Normal to the plane defined by e */
648  robust_cross_product(&(e->start), &(e->end), &normal);
649  normalize(&normal);
650  geog2cart(p, &pt);
651  /* We expect the dot product of with normal with any vector in the plane to be zero */
652  w = dot_product(&normal, &pt);
653  LWDEBUGF(4,"dot product %.9g",w);
654  if ( FP_IS_ZERO(w) )
655  {
656  LWDEBUG(4, "point is on plane (dot product is zero)");
657  return 0;
658  }
659 
660  if ( w < 0 )
661  return -1;
662  else
663  return 1;
664 }
665 
669 static double
671 {
672  POINT3D normal1, normal2;
673  robust_cross_product(b, a, &normal1);
674  robust_cross_product(b, c, &normal2);
675  normalize(&normal1);
676  normalize(&normal2);
677  return sphere_distance_cartesian(&normal1, &normal2);
678 }
679 
689 static double
691 {
692  double angle_a, angle_b, angle_c;
693  double area_radians = 0.0;
694  int side;
695  GEOGRAPHIC_EDGE e;
696 
697  angle_a = sphere_angle(b,a,c);
698  angle_b = sphere_angle(a,b,c);
699  angle_c = sphere_angle(b,c,a);
700 
701  area_radians = angle_a + angle_b + angle_c - M_PI;
702 
703  /* What's the direction of the B/C edge? */
704  e.start = *a;
705  e.end = *b;
706  side = edge_point_side(&e, c);
707 
708  /* Co-linear points implies no area */
709  if ( side == 0 )
710  return 0.0;
711 
712  /* Add the sign to the area */
713  return side * area_radians;
714 }
715 
716 
717 
725 {
726  int side = edge_point_side(e, p);
727  if ( side == 0 )
728  return LW_TRUE;
729 
730  return LW_FALSE;
731 }
732 
738 {
739  POINT3D vcp, vs, ve, vp;
740  double vs_dot_vcp, vp_dot_vcp;
741  geog2cart(&(e->start), &vs);
742  geog2cart(&(e->end), &ve);
743  /* Antipodal case, everything is inside. */
744  if ( vs.x == -1.0 * ve.x && vs.y == -1.0 * ve.y && vs.z == -1.0 * ve.z )
745  return LW_TRUE;
746  geog2cart(p, &vp);
747  /* The normalized sum bisects the angle between start and end. */
748  vector_sum(&vs, &ve, &vcp);
749  normalize(&vcp);
750  /* The projection of start onto the center defines the minimum similarity */
751  vs_dot_vcp = dot_product(&vs, &vcp);
752  LWDEBUGF(4,"vs_dot_vcp %.19g",vs_dot_vcp);
753  /* The projection of candidate p onto the center */
754  vp_dot_vcp = dot_product(&vp, &vcp);
755  LWDEBUGF(4,"vp_dot_vcp %.19g",vp_dot_vcp);
756  /* If p is more similar than start then p is inside the cone */
757  LWDEBUGF(4,"fabs(vp_dot_vcp - vs_dot_vcp) %.39g",fabs(vp_dot_vcp - vs_dot_vcp));
758 
759  /*
760  ** We want to test that vp_dot_vcp is >= vs_dot_vcp but there are
761  ** numerical stability issues for values that are very very nearly
762  ** equal. Unfortunately there are also values of vp_dot_vcp that are legitimately
763  ** very close to but still less than vs_dot_vcp which we also need to catch.
764  ** The tolerance of 10-17 seems to do the trick on 32-bit and 64-bit architectures,
765  ** for the test cases here.
766  ** However, tuning the tolerance value feels like a dangerous hack.
767  ** Fundamentally, the problem is that this test is so sensitive.
768  */
769 
770  /* 1.1102230246251565404236316680908203125e-16 */
771 
772  if ( vp_dot_vcp > vs_dot_vcp || fabs(vp_dot_vcp - vs_dot_vcp) < 2e-16 )
773  {
774  LWDEBUG(4, "point is in cone");
775  return LW_TRUE;
776  }
777  LWDEBUG(4, "point is not in cone");
778  return LW_FALSE;
779 }
780 
785 {
786  GEOGRAPHIC_EDGE g;
788  double slon = fabs((e->start).lon) + fabs((e->end).lon);
789  double dlon = fabs(fabs((e->start).lon) - fabs((e->end).lon));
790  double slat = (e->start).lat + (e->end).lat;
791 
792  LWDEBUGF(4, "e.start == GPOINT(%.6g %.6g) ", (e->start).lat, (e->start).lon);
793  LWDEBUGF(4, "e.end == GPOINT(%.6g %.6g) ", (e->end).lat, (e->end).lon);
794  LWDEBUGF(4, "p == GPOINT(%.6g %.6g) ", p->lat, p->lon);
795 
796  /* Copy values into working registers */
797  g = *e;
798  q = *p;
799 
800  /* Vertical plane, we need to do this calculation in latitude */
801  if ( FP_EQUALS( g.start.lon, g.end.lon ) )
802  {
803  LWDEBUG(4, "vertical plane, we need to do this calculation in latitude");
804  /* Supposed to be co-planar... */
805  if ( ! FP_EQUALS( q.lon, g.start.lon ) )
806  return LW_FALSE;
807 
808  if ( ( g.start.lat <= q.lat && q.lat <= g.end.lat ) ||
809  ( g.end.lat <= q.lat && q.lat <= g.start.lat ) )
810  {
811  return LW_TRUE;
812  }
813  else
814  {
815  return LW_FALSE;
816  }
817  }
818 
819  /* Over the pole, we need normalize latitude and do this calculation in latitude */
820  if ( FP_EQUALS( slon, M_PI ) && ( signum(g.start.lon) != signum(g.end.lon) || FP_EQUALS(dlon, M_PI) ) )
821  {
822  LWDEBUG(4, "over the pole...");
823  /* Antipodal, everything (or nothing?) is inside */
824  if ( FP_EQUALS( slat, 0.0 ) )
825  return LW_TRUE;
826 
827  /* Point *is* the north pole */
828  if ( slat > 0.0 && FP_EQUALS(q.lat, M_PI_2 ) )
829  return LW_TRUE;
830 
831  /* Point *is* the south pole */
832  if ( slat < 0.0 && FP_EQUALS(q.lat, -1.0 * M_PI_2) )
833  return LW_TRUE;
834 
835  LWDEBUG(4, "coplanar?...");
836 
837  /* Supposed to be co-planar... */
838  if ( ! FP_EQUALS( q.lon, g.start.lon ) )
839  return LW_FALSE;
840 
841  LWDEBUG(4, "north or south?...");
842 
843  /* Over north pole, test based on south pole */
844  if ( slat > 0.0 )
845  {
846  LWDEBUG(4, "over the north pole...");
847  if ( q.lat > FP_MIN(g.start.lat, g.end.lat) )
848  return LW_TRUE;
849  else
850  return LW_FALSE;
851  }
852  else
853  /* Over south pole, test based on north pole */
854  {
855  LWDEBUG(4, "over the south pole...");
856  if ( q.lat < FP_MAX(g.start.lat, g.end.lat) )
857  return LW_TRUE;
858  else
859  return LW_FALSE;
860  }
861  }
862 
863  /* Dateline crossing, flip everything to the opposite hemisphere */
864  else if ( slon > M_PI && ( signum(g.start.lon) != signum(g.end.lon) ) )
865  {
866  LWDEBUG(4, "crosses dateline, flip longitudes...");
867  if ( g.start.lon > 0.0 )
868  g.start.lon -= M_PI;
869  else
870  g.start.lon += M_PI;
871  if ( g.end.lon > 0.0 )
872  g.end.lon -= M_PI;
873  else
874  g.end.lon += M_PI;
875 
876  if ( q.lon > 0.0 )
877  q.lon -= M_PI;
878  else
879  q.lon += M_PI;
880  }
881 
882  if ( ( g.start.lon <= q.lon && q.lon <= g.end.lon ) ||
883  ( g.end.lon <= q.lon && q.lon <= g.start.lon ) )
884  {
885  LWDEBUG(4, "true, this edge contains point");
886  return LW_TRUE;
887  }
888 
889  LWDEBUG(4, "false, this edge does not contain point");
890  return LW_FALSE;
891 }
892 
893 
898 {
899  double d_lon = e->lon - s->lon;
900  double cos_d_lon = cos(d_lon);
901  double cos_lat_e = cos(e->lat);
902  double sin_lat_e = sin(e->lat);
903  double cos_lat_s = cos(s->lat);
904  double sin_lat_s = sin(s->lat);
905 
906  double a1 = POW2(cos_lat_e * sin(d_lon));
907  double a2 = POW2(cos_lat_s * sin_lat_e - sin_lat_s * cos_lat_e * cos_d_lon);
908  double a = sqrt(a1 + a2);
909  double b = sin_lat_s * sin_lat_e + cos_lat_s * cos_lat_e * cos_d_lon;
910  return atan2(a, b);
911 }
912 
916 double sphere_distance_cartesian(const POINT3D *s, const POINT3D *e)
917 {
918  return acos(FP_MIN(1.0, dot_product(s, e)));
919 }
920 
924 double sphere_direction(const GEOGRAPHIC_POINT *s, const GEOGRAPHIC_POINT *e, double d)
925 {
926  double heading = 0.0;
927  double f;
928 
929  /* Starting from the poles? Special case. */
930  if ( FP_IS_ZERO(cos(s->lat)) )
931  return (s->lat > 0.0) ? M_PI : 0.0;
932 
933  f = (sin(e->lat) - sin(s->lat) * cos(d)) / (sin(d) * cos(s->lat));
934  if ( FP_EQUALS(f, 1.0) )
935  heading = 0.0;
936  else if ( fabs(f) > 1.0 )
937  {
938  LWDEBUGF(4, "f = %g", f);
939  heading = acos(f);
940  }
941  else
942  heading = acos(f);
943 
944  if ( sin(e->lon - s->lon) < 0.0 )
945  heading = -1 * heading;
946 
947  return heading;
948 }
949 
950 #if 0 /* unused */
951 
962 static double sphere_excess(const GEOGRAPHIC_POINT *a, const GEOGRAPHIC_POINT *b, const GEOGRAPHIC_POINT *c)
963 {
964  double a_dist = sphere_distance(b, c);
965  double b_dist = sphere_distance(c, a);
966  double c_dist = sphere_distance(a, b);
967  double hca = sphere_direction(c, a, b_dist);
968  double hcb = sphere_direction(c, b, a_dist);
969  double sign = signum(hcb-hca);
970  double ss = (a_dist + b_dist + c_dist) / 2.0;
971  double E = tan(ss/2.0)*tan((ss-a_dist)/2.0)*tan((ss-b_dist)/2.0)*tan((ss-c_dist)/2.0);
972  return 4.0 * atan(sqrt(fabs(E))) * sign;
973 }
974 #endif
975 
976 
982 {
983  if ( edge_point_in_cone(e, p) && edge_point_on_plane(e, p) )
984  /* if ( edge_contains_coplanar_point(e, p) && edge_point_on_plane(e, p) ) */
985  {
986  LWDEBUG(4, "point is on edge");
987  return LW_TRUE;
988  }
989  LWDEBUG(4, "point is not on edge");
990  return LW_FALSE;
991 }
992 
996 double z_to_latitude(double z, int top)
997 {
998  double sign = signum(z);
999  double tlat = acos(z);
1000  LWDEBUGF(4, "inputs: z(%.8g) sign(%.8g) tlat(%.8g)", z, sign, tlat);
1001  if (FP_IS_ZERO(z))
1002  {
1003  if (top) return M_PI_2;
1004  else return -1.0 * M_PI_2;
1005  }
1006  if (fabs(tlat) > M_PI_2 )
1007  {
1008  tlat = sign * (M_PI - fabs(tlat));
1009  }
1010  else
1011  {
1012  tlat = sign * tlat;
1013  }
1014  LWDEBUGF(4, "output: tlat(%.8g)", tlat);
1015  return tlat;
1016 }
1017 
1023 int clairaut_cartesian(const POINT3D *start, const POINT3D *end, GEOGRAPHIC_POINT *g_top, GEOGRAPHIC_POINT *g_bottom)
1024 {
1025  POINT3D t1, t2;
1026  GEOGRAPHIC_POINT vN1, vN2;
1027  LWDEBUG(4,"entering function");
1028  unit_normal(start, end, &t1);
1029  unit_normal(end, start, &t2);
1030  LWDEBUGF(4, "unit normal t1 == POINT(%.8g %.8g %.8g)", t1.x, t1.y, t1.z);
1031  LWDEBUGF(4, "unit normal t2 == POINT(%.8g %.8g %.8g)", t2.x, t2.y, t2.z);
1032  cart2geog(&t1, &vN1);
1033  cart2geog(&t2, &vN2);
1034  g_top->lat = z_to_latitude(t1.z,LW_TRUE);
1035  g_top->lon = vN2.lon;
1036  g_bottom->lat = z_to_latitude(t2.z,LW_FALSE);
1037  g_bottom->lon = vN1.lon;
1038  LWDEBUGF(4, "clairaut top == GPOINT(%.6g %.6g)", g_top->lat, g_top->lon);
1039  LWDEBUGF(4, "clairaut bottom == GPOINT(%.6g %.6g)", g_bottom->lat, g_bottom->lon);
1040  return LW_SUCCESS;
1041 }
1042 
1049 {
1050  POINT3D t1, t2;
1051  GEOGRAPHIC_POINT vN1, vN2;
1052  LWDEBUG(4,"entering function");
1053  robust_cross_product(start, end, &t1);
1054  normalize(&t1);
1055  robust_cross_product(end, start, &t2);
1056  normalize(&t2);
1057  LWDEBUGF(4, "unit normal t1 == POINT(%.8g %.8g %.8g)", t1.x, t1.y, t1.z);
1058  LWDEBUGF(4, "unit normal t2 == POINT(%.8g %.8g %.8g)", t2.x, t2.y, t2.z);
1059  cart2geog(&t1, &vN1);
1060  cart2geog(&t2, &vN2);
1061  g_top->lat = z_to_latitude(t1.z,LW_TRUE);
1062  g_top->lon = vN2.lon;
1063  g_bottom->lat = z_to_latitude(t2.z,LW_FALSE);
1064  g_bottom->lon = vN1.lon;
1065  LWDEBUGF(4, "clairaut top == GPOINT(%.6g %.6g)", g_top->lat, g_top->lon);
1066  LWDEBUGF(4, "clairaut bottom == GPOINT(%.6g %.6g)", g_bottom->lat, g_bottom->lon);
1067  return LW_SUCCESS;
1068 }
1069 
1075 {
1076  POINT3D ea, eb, v;
1077  LWDEBUGF(4, "e1 start(%.20g %.20g) end(%.20g %.20g)", e1->start.lat, e1->start.lon, e1->end.lat, e1->end.lon);
1078  LWDEBUGF(4, "e2 start(%.20g %.20g) end(%.20g %.20g)", e2->start.lat, e2->start.lon, e2->end.lat, e2->end.lon);
1079 
1080  LWDEBUGF(4, "e1 start(%.20g %.20g) end(%.20g %.20g)", rad2deg(e1->start.lon), rad2deg(e1->start.lat), rad2deg(e1->end.lon), rad2deg(e1->end.lat));
1081  LWDEBUGF(4, "e2 start(%.20g %.20g) end(%.20g %.20g)", rad2deg(e2->start.lon), rad2deg(e2->start.lat), rad2deg(e2->end.lon), rad2deg(e2->end.lat));
1082 
1083  if ( geographic_point_equals(&(e1->start), &(e2->start)) )
1084  {
1085  *g = e1->start;
1086  return LW_TRUE;
1087  }
1088  if ( geographic_point_equals(&(e1->end), &(e2->end)) )
1089  {
1090  *g = e1->end;
1091  return LW_TRUE;
1092  }
1093  if ( geographic_point_equals(&(e1->end), &(e2->start)) )
1094  {
1095  *g = e1->end;
1096  return LW_TRUE;
1097  }
1098  if ( geographic_point_equals(&(e1->start), &(e2->end)) )
1099  {
1100  *g = e1->start;
1101  return LW_TRUE;
1102  }
1103 
1104  robust_cross_product(&(e1->start), &(e1->end), &ea);
1105  normalize(&ea);
1106  robust_cross_product(&(e2->start), &(e2->end), &eb);
1107  normalize(&eb);
1108  LWDEBUGF(4, "e1 cross product == POINT(%.12g %.12g %.12g)", ea.x, ea.y, ea.z);
1109  LWDEBUGF(4, "e2 cross product == POINT(%.12g %.12g %.12g)", eb.x, eb.y, eb.z);
1110  LWDEBUGF(4, "fabs(dot_product(ea, eb)) == %.14g", fabs(dot_product(&ea, &eb)));
1111  if ( FP_EQUALS(fabs(dot_product(&ea, &eb)), 1.0) )
1112  {
1113  LWDEBUGF(4, "parallel edges found! dot_product = %.12g", dot_product(&ea, &eb));
1114  /* Parallel (maybe equal) edges! */
1115  /* Hack alert, only returning ONE end of the edge right now, most do better later. */
1116  /* Hack alart #2, returning a value of 2 to indicate a co-linear crossing event. */
1117  if ( edge_contains_point(e1, &(e2->start)) )
1118  {
1119  *g = e2->start;
1120  return 2;
1121  }
1122  if ( edge_contains_point(e1, &(e2->end)) )
1123  {
1124  *g = e2->end;
1125  return 2;
1126  }
1127  if ( edge_contains_point(e2, &(e1->start)) )
1128  {
1129  *g = e1->start;
1130  return 2;
1131  }
1132  if ( edge_contains_point(e2, &(e1->end)) )
1133  {
1134  *g = e1->end;
1135  return 2;
1136  }
1137  }
1138  unit_normal(&ea, &eb, &v);
1139  LWDEBUGF(4, "v == POINT(%.12g %.12g %.12g)", v.x, v.y, v.z);
1140  g->lat = atan2(v.z, sqrt(v.x * v.x + v.y * v.y));
1141  g->lon = atan2(v.y, v.x);
1142  LWDEBUGF(4, "g == GPOINT(%.12g %.12g)", g->lat, g->lon);
1143  LWDEBUGF(4, "g == POINT(%.12g %.12g)", rad2deg(g->lon), rad2deg(g->lat));
1144  if ( edge_contains_point(e1, g) && edge_contains_point(e2, g) )
1145  {
1146  return LW_TRUE;
1147  }
1148  else
1149  {
1150  LWDEBUG(4, "flipping point to other side of sphere");
1151  g->lat = -1.0 * g->lat;
1152  g->lon = g->lon + M_PI;
1153  if ( g->lon > M_PI )
1154  {
1155  g->lon = -1.0 * (2.0 * M_PI - g->lon);
1156  }
1157  if ( edge_contains_point(e1, g) && edge_contains_point(e2, g) )
1158  {
1159  return LW_TRUE;
1160  }
1161  }
1162  return LW_FALSE;
1163 }
1164 
1166 {
1167  double d1 = 1000000000.0, d2, d3, d_nearest;
1168  POINT3D n, p, k;
1169  GEOGRAPHIC_POINT gk, g_nearest;
1170 
1171  /* Zero length edge, */
1172  if ( geographic_point_equals(&(e->start), &(e->end)) )
1173  {
1174  *closest = e->start;
1175  return sphere_distance(&(e->start), gp);
1176  }
1177 
1178  robust_cross_product(&(e->start), &(e->end), &n);
1179  normalize(&n);
1180  geog2cart(gp, &p);
1181  vector_scale(&n, dot_product(&p, &n));
1182  vector_difference(&p, &n, &k);
1183  normalize(&k);
1184  cart2geog(&k, &gk);
1185  if ( edge_contains_point(e, &gk) )
1186  {
1187  d1 = sphere_distance(gp, &gk);
1188  }
1189  d2 = sphere_distance(gp, &(e->start));
1190  d3 = sphere_distance(gp, &(e->end));
1191 
1192  d_nearest = d1;
1193  g_nearest = gk;
1194 
1195  if ( d2 < d_nearest )
1196  {
1197  d_nearest = d2;
1198  g_nearest = e->start;
1199  }
1200  if ( d3 < d_nearest )
1201  {
1202  d_nearest = d3;
1203  g_nearest = e->end;
1204  }
1205  if (closest)
1206  *closest = g_nearest;
1207 
1208  return d_nearest;
1209 }
1210 
1217 {
1218  double d;
1219  GEOGRAPHIC_POINT gcp1s, gcp1e, gcp2s, gcp2e, c1, c2;
1220  double d1s = edge_distance_to_point(e1, &(e2->start), &gcp1s);
1221  double d1e = edge_distance_to_point(e1, &(e2->end), &gcp1e);
1222  double d2s = edge_distance_to_point(e2, &(e1->start), &gcp2s);
1223  double d2e = edge_distance_to_point(e2, &(e1->end), &gcp2e);
1224 
1225  d = d1s;
1226  c1 = gcp1s;
1227  c2 = e2->start;
1228 
1229  if ( d1e < d )
1230  {
1231  d = d1e;
1232  c1 = gcp1e;
1233  c2 = e2->end;
1234  }
1235 
1236  if ( d2s < d )
1237  {
1238  d = d2s;
1239  c1 = e1->start;
1240  c2 = gcp2s;
1241  }
1242 
1243  if ( d2e < d )
1244  {
1245  d = d2e;
1246  c1 = e1->end;
1247  c2 = gcp2e;
1248  }
1249 
1250  if ( closest1 ) *closest1 = c1;
1251  if ( closest2 ) *closest2 = c2;
1252 
1253  return d;
1254 }
1255 
1256 
1261 int sphere_project(const GEOGRAPHIC_POINT *r, double distance, double azimuth, GEOGRAPHIC_POINT *n)
1262 {
1263  double d = distance;
1264  double lat1 = r->lat;
1265  double lon1 = r->lon;
1266  double lat2, lon2;
1267 
1268  lat2 = asin(sin(lat1)*cos(d) + cos(lat1)*sin(d)*cos(azimuth));
1269 
1270  /* If we're going straight up or straight down, we don't need to calculate the longitude */
1271  /* TODO: this isn't quite true, what if we're going over the pole? */
1272  if ( FP_EQUALS(azimuth, M_PI) || FP_EQUALS(azimuth, 0.0) )
1273  {
1274  lon2 = r->lon;
1275  }
1276  else
1277  {
1278  lon2 = lon1 + atan2(sin(azimuth)*sin(d)*cos(lat1), cos(d)-sin(lat1)*sin(lat2));
1279  }
1280 
1281  if ( isnan(lat2) || isnan(lon2) )
1282  return LW_FAILURE;
1283 
1284  n->lat = lat2;
1285  n->lon = lon2;
1286 
1287  return LW_SUCCESS;
1288 }
1289 
1290 
1292 {
1293  int steps = 1000000;
1294  int i;
1295  double dx, dy, dz;
1296  double distance = sphere_distance(&(e->start), &(e->end));
1297  POINT3D pn, p, start, end;
1298 
1299  /* Edge is zero length, just return the naive box */
1300  if ( FP_IS_ZERO(distance) )
1301  {
1302  LWDEBUG(4, "edge is zero length. returning");
1303  geog2cart(&(e->start), &start);
1304  geog2cart(&(e->end), &end);
1305  gbox_init_point3d(&start, gbox);
1306  gbox_merge_point3d(&end, gbox);
1307  return LW_SUCCESS;
1308  }
1309 
1310  /* Edge is antipodal (one point on each side of the globe),
1311  set the box to contain the whole world and return */
1312  if ( FP_EQUALS(distance, M_PI) )
1313  {
1314  LWDEBUG(4, "edge is antipodal. setting to maximum size box, and returning");
1315  gbox->xmin = gbox->ymin = gbox->zmin = -1.0;
1316  gbox->xmax = gbox->ymax = gbox->zmax = 1.0;
1317  return LW_SUCCESS;
1318  }
1319 
1320  /* Walk along the chord between start and end incrementally,
1321  normalizing at each step. */
1322  geog2cart(&(e->start), &start);
1323  geog2cart(&(e->end), &end);
1324  dx = (end.x - start.x)/steps;
1325  dy = (end.y - start.y)/steps;
1326  dz = (end.z - start.z)/steps;
1327  p = start;
1328  gbox->xmin = gbox->xmax = p.x;
1329  gbox->ymin = gbox->ymax = p.y;
1330  gbox->zmin = gbox->zmax = p.z;
1331  for ( i = 0; i < steps; i++ )
1332  {
1333  p.x += dx;
1334  p.y += dy;
1335  p.z += dz;
1336  pn = p;
1337  normalize(&pn);
1338  gbox_merge_point3d(&pn, gbox);
1339  }
1340  return LW_SUCCESS;
1341 }
1342 
1355 int edge_calculate_gbox(const POINT3D *A1, const POINT3D *A2, GBOX *gbox)
1356 {
1357  POINT2D R1, R2, RX, O;
1358  POINT3D AN, A3;
1359  POINT3D X[6];
1360  int i, o_side;
1361 
1362  /* Initialize the box with the edge end points */
1363  gbox_init_point3d(A1, gbox);
1364  gbox_merge_point3d(A2, gbox);
1365 
1366  /* Zero length edge, just return! */
1367  if ( p3d_same(A1, A2) )
1368  return LW_SUCCESS;
1369 
1370  /* Error out on antipodal edge */
1371  if ( FP_EQUALS(A1->x, -1*A2->x) && FP_EQUALS(A1->y, -1*A2->y) && FP_EQUALS(A1->z, -1*A2->z) )
1372  {
1373  lwerror("Antipodal (180 degrees long) edge detected!");
1374  return LW_FAILURE;
1375  }
1376 
1377  /* Create A3, a vector in the plane of A1/A2, orthogonal to A1 */
1378  unit_normal(A1, A2, &AN);
1379  unit_normal(&AN, A1, &A3);
1380 
1381  /* Project A1 and A2 into the 2-space formed by the plane A1/A3 */
1382  R1.x = 1.0;
1383  R1.y = 0.0;
1384  R2.x = dot_product(A2, A1);
1385  R2.y = dot_product(A2, &A3);
1386 
1387  /* Initialize our 3-space axis points (x+, x-, y+, y-, z+, z-) */
1388  memset(X, 0, sizeof(POINT3D) * 6);
1389  X[0].x = X[2].y = X[4].z = 1.0;
1390  X[1].x = X[3].y = X[5].z = -1.0;
1391 
1392  /* Initialize a 2-space origin point. */
1393  O.x = O.y = 0.0;
1394  /* What side of the line joining R1/R2 is O? */
1395  o_side = lw_segment_side(&R1, &R2, &O);
1396 
1397  /* Add any extrema! */
1398  for ( i = 0; i < 6; i++ )
1399  {
1400  /* Convert 3-space axis points to 2-space unit vectors */
1401  RX.x = dot_product(&(X[i]), A1);
1402  RX.y = dot_product(&(X[i]), &A3);
1403  normalize2d(&RX);
1404 
1405  /* Any axis end on the side of R1/R2 opposite the origin */
1406  /* is an extreme point in the arc, so we add the 3-space */
1407  /* version of the point on R1/R2 to the gbox */
1408  if ( lw_segment_side(&R1, &R2, &RX) != o_side )
1409  {
1410  POINT3D Xn;
1411  Xn.x = RX.x * A1->x + RX.y * A3.x;
1412  Xn.y = RX.x * A1->y + RX.y * A3.y;
1413  Xn.z = RX.x * A1->z + RX.y * A3.z;
1414 
1415  gbox_merge_point3d(&Xn, gbox);
1416  }
1417  }
1418 
1419  return LW_SUCCESS;
1420 }
1421 
1422 void lwpoly_pt_outside(const LWPOLY *poly, POINT2D *pt_outside)
1423 {
1424  /* Make sure we have boxes */
1425  if ( poly->bbox )
1426  {
1427  gbox_pt_outside(poly->bbox, pt_outside);
1428  return;
1429  }
1430  else
1431  {
1432  GBOX gbox;
1433  lwgeom_calculate_gbox_geodetic((LWGEOM*)poly, &gbox);
1434  gbox_pt_outside(&gbox, pt_outside);
1435  return;
1436  }
1437 }
1438 
1443 void gbox_pt_outside(const GBOX *gbox, POINT2D *pt_outside)
1444 {
1445  double grow = M_PI / 180.0 / 60.0; /* one arc-minute */
1446  int i;
1447  GBOX ge;
1448  POINT3D corners[8];
1449  POINT3D pt;
1450  GEOGRAPHIC_POINT g;
1451 
1452  while ( grow < M_PI )
1453  {
1454  /* Assign our box and expand it slightly. */
1455  ge = *gbox;
1456  if ( ge.xmin > -1 ) ge.xmin -= grow;
1457  if ( ge.ymin > -1 ) ge.ymin -= grow;
1458  if ( ge.zmin > -1 ) ge.zmin -= grow;
1459  if ( ge.xmax < 1 ) ge.xmax += grow;
1460  if ( ge.ymax < 1 ) ge.ymax += grow;
1461  if ( ge.zmax < 1 ) ge.zmax += grow;
1462 
1463  /* Build our eight corner points */
1464  corners[0].x = ge.xmin;
1465  corners[0].y = ge.ymin;
1466  corners[0].z = ge.zmin;
1467 
1468  corners[1].x = ge.xmin;
1469  corners[1].y = ge.ymax;
1470  corners[1].z = ge.zmin;
1471 
1472  corners[2].x = ge.xmin;
1473  corners[2].y = ge.ymin;
1474  corners[2].z = ge.zmax;
1475 
1476  corners[3].x = ge.xmax;
1477  corners[3].y = ge.ymin;
1478  corners[3].z = ge.zmin;
1479 
1480  corners[4].x = ge.xmax;
1481  corners[4].y = ge.ymax;
1482  corners[4].z = ge.zmin;
1483 
1484  corners[5].x = ge.xmax;
1485  corners[5].y = ge.ymin;
1486  corners[5].z = ge.zmax;
1487 
1488  corners[6].x = ge.xmin;
1489  corners[6].y = ge.ymax;
1490  corners[6].z = ge.zmax;
1491 
1492  corners[7].x = ge.xmax;
1493  corners[7].y = ge.ymax;
1494  corners[7].z = ge.zmax;
1495 
1496  LWDEBUG(4, "trying to use a box corner point...");
1497  for ( i = 0; i < 8; i++ )
1498  {
1499  normalize(&(corners[i]));
1500  LWDEBUGF(4, "testing corner %d: POINT(%.8g %.8g %.8g)", i, corners[i].x, corners[i].y, corners[i].z);
1501  if ( ! gbox_contains_point3d(gbox, &(corners[i])) )
1502  {
1503  LWDEBUGF(4, "corner %d is outside our gbox", i);
1504  pt = corners[i];
1505  normalize(&pt);
1506  cart2geog(&pt, &g);
1507  pt_outside->x = rad2deg(g.lon);
1508  pt_outside->y = rad2deg(g.lat);
1509  LWDEBUGF(4, "returning POINT(%.8g %.8g) as outside point", pt_outside->x, pt_outside->y);
1510  return;
1511  }
1512  }
1513 
1514  /* Try a wider growth to push the corners outside the original box. */
1515  grow *= 2.0;
1516  }
1517 
1518  /* This should never happen! */
1519  lwerror("BOOM! Could not generate outside point!");
1520  return;
1521 }
1522 
1523 
1529 static POINTARRAY*
1530 ptarray_segmentize_sphere(const POINTARRAY *pa_in, double max_seg_length)
1531 {
1532  POINTARRAY *pa_out;
1533  int hasz = ptarray_has_z(pa_in);
1534  int hasm = ptarray_has_m(pa_in);
1535  int pa_in_offset = 0; /* input point offset */
1536  POINT4D p1, p2, p;
1537  POINT3D q1, q2, q, qn;
1538  GEOGRAPHIC_POINT g1, g2, g;
1539  double d;
1540 
1541  /* Just crap out on crazy input */
1542  if ( ! pa_in )
1543  lwerror("ptarray_segmentize_sphere: null input pointarray");
1544  if ( max_seg_length <= 0.0 )
1545  lwerror("ptarray_segmentize_sphere: maximum segment length must be positive");
1546 
1547  /* Empty starting array */
1548  pa_out = ptarray_construct_empty(hasz, hasm, pa_in->npoints);
1549 
1550  /* Add first point */
1551  getPoint4d_p(pa_in, pa_in_offset, &p1);
1552  ptarray_append_point(pa_out, &p1, LW_FALSE);
1553  geographic_point_init(p1.x, p1.y, &g1);
1554  pa_in_offset++;
1555 
1556  while ( pa_in_offset < pa_in->npoints )
1557  {
1558  getPoint4d_p(pa_in, pa_in_offset, &p2);
1559  geographic_point_init(p2.x, p2.y, &g2);
1560 
1561  /* Skip duplicate points (except in case of 2-point lines!) */
1562  if ( (pa_in->npoints > 2) && p4d_same(&p1, &p2) )
1563  {
1564  /* Move one offset forward */
1565  p1 = p2;
1566  g1 = g2;
1567  pa_in_offset++;
1568  continue;
1569  }
1570 
1571  /* How long is this edge? */
1572  d = sphere_distance(&g1, &g2);
1573 
1574  /* We need to segmentize this edge */
1575  if ( d > max_seg_length )
1576  {
1577  int nsegs = 1 + d / max_seg_length;
1578  int i;
1579  double dx, dy, dz, dzz = 0, dmm = 0;
1580 
1581  geog2cart(&g1, &q1);
1582  geog2cart(&g2, &q2);
1583 
1584  dx = (q2.x - q1.x) / nsegs;
1585  dy = (q2.y - q1.y) / nsegs;
1586  dz = (q2.z - q1.z) / nsegs;
1587 
1588  /* The independent Z/M values on the ptarray */
1589  if ( hasz ) dzz = (p2.z - p1.z) / nsegs;
1590  if ( hasm ) dmm = (p2.m - p1.m) / nsegs;
1591 
1592  q = q1;
1593  p = p1;
1594 
1595  for ( i = 0; i < nsegs - 1; i++ )
1596  {
1597  /* Move one increment forwards */
1598  q.x += dx; q.y += dy; q.z += dz;
1599  qn = q;
1600  normalize(&qn);
1601 
1602  /* Back to spherical coordinates */
1603  cart2geog(&qn, &g);
1604  /* Back to lon/lat */
1605  p.x = rad2deg(g.lon);
1606  p.y = rad2deg(g.lat);
1607  if ( hasz )
1608  p.z += dzz;
1609  if ( hasm )
1610  p.m += dmm;
1611  ptarray_append_point(pa_out, &p, LW_FALSE);
1612  }
1613 
1614  ptarray_append_point(pa_out, &p2, LW_FALSE);
1615  }
1616  /* This edge is already short enough */
1617  else
1618  {
1619  ptarray_append_point(pa_out, &p2, (pa_in->npoints==2)?LW_TRUE:LW_FALSE);
1620  }
1621 
1622  /* Move one offset forward */
1623  p1 = p2;
1624  g1 = g2;
1625  pa_in_offset++;
1626  }
1627 
1628  return pa_out;
1629 }
1630 
1637 LWGEOM*
1638 lwgeom_segmentize_sphere(const LWGEOM *lwg_in, double max_seg_length)
1639 {
1640  POINTARRAY *pa_out;
1641  LWLINE *lwline;
1642  LWPOLY *lwpoly_in, *lwpoly_out;
1643  LWCOLLECTION *lwcol_in, *lwcol_out;
1644  int i;
1645 
1646  /* Reflect NULL */
1647  if ( ! lwg_in )
1648  return NULL;
1649 
1650  /* Clone empty */
1651  if ( lwgeom_is_empty(lwg_in) )
1652  return lwgeom_clone(lwg_in);
1653 
1654  switch (lwg_in->type)
1655  {
1656  case MULTIPOINTTYPE:
1657  case POINTTYPE:
1658  return lwgeom_clone_deep(lwg_in);
1659  break;
1660  case LINETYPE:
1661  lwline = lwgeom_as_lwline(lwg_in);
1662  pa_out = ptarray_segmentize_sphere(lwline->points, max_seg_length);
1663  return lwline_as_lwgeom(lwline_construct(lwg_in->srid, NULL, pa_out));
1664  break;
1665  case POLYGONTYPE:
1666  lwpoly_in = lwgeom_as_lwpoly(lwg_in);
1667  lwpoly_out = lwpoly_construct_empty(lwg_in->srid, lwgeom_has_z(lwg_in), lwgeom_has_m(lwg_in));
1668  for ( i = 0; i < lwpoly_in->nrings; i++ )
1669  {
1670  pa_out = ptarray_segmentize_sphere(lwpoly_in->rings[i], max_seg_length);
1671  lwpoly_add_ring(lwpoly_out, pa_out);
1672  }
1673  return lwpoly_as_lwgeom(lwpoly_out);
1674  break;
1675  case MULTILINETYPE:
1676  case MULTIPOLYGONTYPE:
1677  case COLLECTIONTYPE:
1678  lwcol_in = lwgeom_as_lwcollection(lwg_in);
1679  lwcol_out = lwcollection_construct_empty(lwg_in->type, lwg_in->srid, lwgeom_has_z(lwg_in), lwgeom_has_m(lwg_in));
1680  for ( i = 0; i < lwcol_in->ngeoms; i++ )
1681  {
1682  lwcollection_add_lwgeom(lwcol_out, lwgeom_segmentize_sphere(lwcol_in->geoms[i], max_seg_length));
1683  }
1684  return lwcollection_as_lwgeom(lwcol_out);
1685  break;
1686  default:
1687  lwerror("lwgeom_segmentize_sphere: unsupported input geometry type: %d - %s",
1688  lwg_in->type, lwtype_name(lwg_in->type));
1689  break;
1690  }
1691 
1692  lwerror("lwgeom_segmentize_sphere got to the end of the function, should not happen");
1693  return NULL;
1694 }
1695 
1696 
1701 double
1703 {
1704  int i;
1705  const POINT2D *p;
1706  GEOGRAPHIC_POINT a, b, c;
1707  double area = 0.0;
1708 
1709  /* Return zero on nonsensical inputs */
1710  if ( ! pa || pa->npoints < 4 )
1711  return 0.0;
1712 
1713  p = getPoint2d_cp(pa, 0);
1714  geographic_point_init(p->x, p->y, &a);
1715  p = getPoint2d_cp(pa, 1);
1716  geographic_point_init(p->x, p->y, &b);
1717 
1718  for ( i = 2; i < pa->npoints-1; i++ )
1719  {
1720  p = getPoint2d_cp(pa, i);
1721  geographic_point_init(p->x, p->y, &c);
1722  area += sphere_signed_area(&a, &b, &c);
1723  b = c;
1724  }
1725 
1726  return fabs(area);
1727 }
1728 
1729 
1730 static double ptarray_distance_spheroid(const POINTARRAY *pa1, const POINTARRAY *pa2, const SPHEROID *s, double tolerance, int check_intersection)
1731 {
1732  GEOGRAPHIC_EDGE e1, e2;
1733  GEOGRAPHIC_POINT g1, g2;
1734  GEOGRAPHIC_POINT nearest1, nearest2;
1735  POINT3D A1, A2, B1, B2;
1736  POINT2D p;
1737  double distance;
1738  int i, j;
1739  int use_sphere = (s->a == s->b ? 1 : 0);
1740 
1741  /* Make result really big, so that everything will be smaller than it */
1742  distance = MAXFLOAT;
1743 
1744  /* Empty point arrays? Return negative */
1745  if ( pa1->npoints == 0 || pa2->npoints == 0 )
1746  return -1.0;
1747 
1748  /* Handle point/point case here */
1749  if ( pa1->npoints == 1 && pa2->npoints == 1 )
1750  {
1751  getPoint2d_p(pa1, 0, &p);
1752  geographic_point_init(p.x, p.y, &g1);
1753  getPoint2d_p(pa2, 0, &p);
1754  geographic_point_init(p.x, p.y, &g2);
1755  /* Sphere special case, axes equal */
1756  distance = s->radius * sphere_distance(&g1, &g2);
1757  if ( use_sphere )
1758  return distance;
1759  /* Below tolerance, actual distance isn't of interest */
1760  else if ( distance < 0.95 * tolerance )
1761  return distance;
1762  /* Close or greater than tolerance, get the real answer to be sure */
1763  else
1764  return spheroid_distance(&g1, &g2, s);
1765  }
1766 
1767  /* Handle point/line case here */
1768  if ( pa1->npoints == 1 || pa2->npoints == 1 )
1769  {
1770  /* Handle one/many case here */
1771  int i;
1772  const POINTARRAY *pa_one;
1773  const POINTARRAY *pa_many;
1774 
1775  if ( pa1->npoints == 1 )
1776  {
1777  pa_one = pa1;
1778  pa_many = pa2;
1779  }
1780  else
1781  {
1782  pa_one = pa2;
1783  pa_many = pa1;
1784  }
1785 
1786  /* Initialize our point */
1787  getPoint2d_p(pa_one, 0, &p);
1788  geographic_point_init(p.x, p.y, &g1);
1789 
1790  /* Initialize start of line */
1791  getPoint2d_p(pa_many, 0, &p);
1792  geographic_point_init(p.x, p.y, &(e1.start));
1793 
1794  /* Iterate through the edges in our line */
1795  for ( i = 1; i < pa_many->npoints; i++ )
1796  {
1797  double d;
1798  getPoint2d_p(pa_many, i, &p);
1799  geographic_point_init(p.x, p.y, &(e1.end));
1800  /* Get the spherical distance between point and edge */
1801  d = s->radius * edge_distance_to_point(&e1, &g1, &g2);
1802  /* New shortest distance! Record this distance / location */
1803  if ( d < distance )
1804  {
1805  distance = d;
1806  nearest2 = g2;
1807  }
1808  /* We've gotten closer than the tolerance... */
1809  if ( d < tolerance )
1810  {
1811  /* Working on a sphere? The answer is correct, return */
1812  if ( use_sphere )
1813  {
1814  return d;
1815  }
1816  /* Far enough past the tolerance that the spheroid calculation won't change things */
1817  else if ( d < tolerance * 0.95 )
1818  {
1819  return d;
1820  }
1821  /* On a spheroid and near the tolerance? Confirm that we are *actually* closer than tolerance */
1822  else
1823  {
1824  d = spheroid_distance(&g1, &nearest2, s);
1825  /* Yes, closer than tolerance, return! */
1826  if ( d < tolerance )
1827  return d;
1828  }
1829  }
1830  e1.start = e1.end;
1831  }
1832 
1833  /* On sphere, return answer */
1834  if ( use_sphere )
1835  return distance;
1836  /* On spheroid, calculate final answer based on closest approach */
1837  else
1838  return spheroid_distance(&g1, &nearest2, s);
1839  }
1840 
1841  /* Initialize start of line 1 */
1842  getPoint2d_p(pa1, 0, &p);
1843  geographic_point_init(p.x, p.y, &(e1.start));
1844  geog2cart(&(e1.start), &A1);
1845 
1846 
1847  /* Handle line/line case */
1848  for ( i = 1; i < pa1->npoints; i++ )
1849  {
1850  getPoint2d_p(pa1, i, &p);
1851  geographic_point_init(p.x, p.y, &(e1.end));
1852  geog2cart(&(e1.end), &A2);
1853 
1854  /* Initialize start of line 2 */
1855  getPoint2d_p(pa2, 0, &p);
1856  geographic_point_init(p.x, p.y, &(e2.start));
1857  geog2cart(&(e2.start), &B1);
1858 
1859  for ( j = 1; j < pa2->npoints; j++ )
1860  {
1861  double d;
1862 
1863  getPoint2d_p(pa2, j, &p);
1864  geographic_point_init(p.x, p.y, &(e2.end));
1865  geog2cart(&(e2.end), &B2);
1866 
1867  LWDEBUGF(4, "e1.start == GPOINT(%.6g %.6g) ", e1.start.lat, e1.start.lon);
1868  LWDEBUGF(4, "e1.end == GPOINT(%.6g %.6g) ", e1.end.lat, e1.end.lon);
1869  LWDEBUGF(4, "e2.start == GPOINT(%.6g %.6g) ", e2.start.lat, e2.start.lon);
1870  LWDEBUGF(4, "e2.end == GPOINT(%.6g %.6g) ", e2.end.lat, e2.end.lon);
1871 
1872  if ( check_intersection && edge_intersects(&A1, &A2, &B1, &B2) )
1873  {
1874  LWDEBUG(4,"edge intersection! returning 0.0");
1875  return 0.0;
1876  }
1877  d = s->radius * edge_distance_to_edge(&e1, &e2, &g1, &g2);
1878  LWDEBUGF(4,"got edge_distance_to_edge %.8g", d);
1879 
1880  if ( d < distance )
1881  {
1882  distance = d;
1883  nearest1 = g1;
1884  nearest2 = g2;
1885  }
1886  if ( d < tolerance )
1887  {
1888  if ( use_sphere )
1889  {
1890  return d;
1891  }
1892  else
1893  {
1894  d = spheroid_distance(&nearest1, &nearest2, s);
1895  if ( d < tolerance )
1896  return d;
1897  }
1898  }
1899 
1900  /* Copy end to start to allow a new end value in next iteration */
1901  e2.start = e2.end;
1902  B1 = B2;
1903  }
1904 
1905  /* Copy end to start to allow a new end value in next iteration */
1906  e1.start = e1.end;
1907  A1 = A2;
1908  }
1909  LWDEBUGF(4,"finished all loops, returning %.8g", distance);
1910 
1911  if ( use_sphere )
1912  return distance;
1913  else
1914  return spheroid_distance(&nearest1, &nearest2, s);
1915 }
1916 
1917 
1924 double lwgeom_area_sphere(const LWGEOM *lwgeom, const SPHEROID *spheroid)
1925 {
1926  int type;
1927  double radius2 = spheroid->radius * spheroid->radius;
1928 
1929  assert(lwgeom);
1930 
1931  /* No area in nothing */
1932  if ( lwgeom_is_empty(lwgeom) )
1933  return 0.0;
1934 
1935  /* Read the geometry type number */
1936  type = lwgeom->type;
1937 
1938  /* Anything but polygons and collections returns zero */
1939  if ( ! ( type == POLYGONTYPE || type == MULTIPOLYGONTYPE || type == COLLECTIONTYPE ) )
1940  return 0.0;
1941 
1942  /* Actually calculate area */
1943  if ( type == POLYGONTYPE )
1944  {
1945  LWPOLY *poly = (LWPOLY*)lwgeom;
1946  int i;
1947  double area = 0.0;
1948 
1949  /* Just in case there's no rings */
1950  if ( poly->nrings < 1 )
1951  return 0.0;
1952 
1953  /* First, the area of the outer ring */
1954  area += radius2 * ptarray_area_sphere(poly->rings[0]);
1955 
1956  /* Subtract areas of inner rings */
1957  for ( i = 1; i < poly->nrings; i++ )
1958  {
1959  area -= radius2 * ptarray_area_sphere(poly->rings[i]);
1960  }
1961  return area;
1962  }
1963 
1964  /* Recurse into sub-geometries to get area */
1965  if ( type == MULTIPOLYGONTYPE || type == COLLECTIONTYPE )
1966  {
1967  LWCOLLECTION *col = (LWCOLLECTION*)lwgeom;
1968  int i;
1969  double area = 0.0;
1970 
1971  for ( i = 0; i < col->ngeoms; i++ )
1972  {
1973  area += lwgeom_area_sphere(col->geoms[i], spheroid);
1974  }
1975  return area;
1976  }
1977 
1978  /* Shouldn't get here. */
1979  return 0.0;
1980 }
1981 
1982 
1992 LWPOINT* lwgeom_project_spheroid(const LWPOINT *r, const SPHEROID *spheroid, double distance, double azimuth)
1993 {
1994  GEOGRAPHIC_POINT geo_source, geo_dest;
1995  POINT4D pt_dest;
1996  double x, y;
1997  POINTARRAY *pa;
1998  LWPOINT *lwp;
1999 
2000  /* Check the azimuth validity, convert to radians */
2001  if ( azimuth < -2.0 * M_PI || azimuth > 2.0 * M_PI )
2002  {
2003  lwerror("Azimuth must be between -2PI and 2PI");
2004  return NULL;
2005  }
2006 
2007  /* Check the distance validity */
2008  if ( distance < 0.0 || distance > (M_PI * spheroid->radius) )
2009  {
2010  lwerror("Distance must be between 0 and %g", M_PI * spheroid->radius);
2011  return NULL;
2012  }
2013 
2014  /* Convert to ta geodetic point */
2015  x = lwpoint_get_x(r);
2016  y = lwpoint_get_y(r);
2017  geographic_point_init(x, y, &geo_source);
2018 
2019  /* Try the projection */
2020  if( spheroid_project(&geo_source, spheroid, distance, azimuth, &geo_dest) == LW_FAILURE )
2021  {
2022  LWDEBUGF(3, "Unable to project from (%g %g) with azimuth %g and distance %g", x, y, azimuth, distance);
2023  lwerror("Unable to project from (%g %g) with azimuth %g and distance %g", x, y, azimuth, distance);
2024  return NULL;
2025  }
2026 
2027  /* Build the output LWPOINT */
2028  pa = ptarray_construct(0, 0, 1);
2029  pt_dest.x = rad2deg(longitude_radians_normalize(geo_dest.lon));
2030  pt_dest.y = rad2deg(latitude_radians_normalize(geo_dest.lat));
2031  pt_dest.z = pt_dest.m = 0.0;
2032  ptarray_set_point4d(pa, 0, &pt_dest);
2033  lwp = lwpoint_construct(r->srid, NULL, pa);
2035  return lwp;
2036 }
2037 
2038 
2047 double lwgeom_azumith_spheroid(const LWPOINT *r, const LWPOINT *s, const SPHEROID *spheroid)
2048 {
2049  GEOGRAPHIC_POINT g1, g2;
2050  double x1, y1, x2, y2;
2051 
2052  /* Convert r to a geodetic point */
2053  x1 = lwpoint_get_x(r);
2054  y1 = lwpoint_get_y(r);
2055  geographic_point_init(x1, y1, &g1);
2056 
2057  /* Convert s to a geodetic point */
2058  x2 = lwpoint_get_x(s);
2059  y2 = lwpoint_get_y(s);
2060  geographic_point_init(x2, y2, &g2);
2061 
2062  /* Same point, return NaN */
2063  if ( FP_EQUALS(x1, x2) && FP_EQUALS(y1, y2) )
2064  {
2065  return NAN;
2066  }
2067 
2068  /* Do the direction calculation */
2069  return spheroid_direction(&g1, &g2, spheroid);
2070 }
2071 
2078 double lwgeom_distance_spheroid(const LWGEOM *lwgeom1, const LWGEOM *lwgeom2, const SPHEROID *spheroid, double tolerance)
2079 {
2080  uint8_t type1, type2;
2081  int check_intersection = LW_FALSE;
2082  GBOX gbox1, gbox2;
2083 
2084  gbox_init(&gbox1);
2085  gbox_init(&gbox2);
2086 
2087  assert(lwgeom1);
2088  assert(lwgeom2);
2089 
2090  LWDEBUGF(4, "entered function, tolerance %.8g", tolerance);
2091 
2092  /* What's the distance to an empty geometry? We don't know.
2093  Return a negative number so the caller can catch this case. */
2094  if ( lwgeom_is_empty(lwgeom1) || lwgeom_is_empty(lwgeom2) )
2095  {
2096  return -1.0;
2097  }
2098 
2099  type1 = lwgeom1->type;
2100  type2 = lwgeom2->type;
2101 
2102  /* Make sure we have boxes */
2103  if ( lwgeom1->bbox )
2104  gbox1 = *(lwgeom1->bbox);
2105  else
2106  lwgeom_calculate_gbox_geodetic(lwgeom1, &gbox1);
2107 
2108  /* Make sure we have boxes */
2109  if ( lwgeom2->bbox )
2110  gbox2 = *(lwgeom2->bbox);
2111  else
2112  lwgeom_calculate_gbox_geodetic(lwgeom2, &gbox2);
2113 
2114  /* If the boxes aren't disjoint, we have to check for edge intersections */
2115  if ( gbox_overlaps(&gbox1, &gbox2) )
2116  check_intersection = LW_TRUE;
2117 
2118  /* Point/line combinations can all be handled with simple point array iterations */
2119  if ( ( type1 == POINTTYPE || type1 == LINETYPE ) &&
2120  ( type2 == POINTTYPE || type2 == LINETYPE ) )
2121  {
2122  POINTARRAY *pa1, *pa2;
2123 
2124  if ( type1 == POINTTYPE )
2125  pa1 = ((LWPOINT*)lwgeom1)->point;
2126  else
2127  pa1 = ((LWLINE*)lwgeom1)->points;
2128 
2129  if ( type2 == POINTTYPE )
2130  pa2 = ((LWPOINT*)lwgeom2)->point;
2131  else
2132  pa2 = ((LWLINE*)lwgeom2)->points;
2133 
2134  return ptarray_distance_spheroid(pa1, pa2, spheroid, tolerance, check_intersection);
2135  }
2136 
2137  /* Point/Polygon cases, if point-in-poly, return zero, else return distance. */
2138  if ( ( type1 == POLYGONTYPE && type2 == POINTTYPE ) ||
2139  ( type2 == POLYGONTYPE && type1 == POINTTYPE ) )
2140  {
2141  POINT2D p;
2142  LWPOLY *lwpoly;
2143  LWPOINT *lwpt;
2144  double distance = MAXFLOAT;
2145  int i;
2146 
2147  if ( type1 == POINTTYPE )
2148  {
2149  lwpt = (LWPOINT*)lwgeom1;
2150  lwpoly = (LWPOLY*)lwgeom2;
2151  }
2152  else
2153  {
2154  lwpt = (LWPOINT*)lwgeom2;
2155  lwpoly = (LWPOLY*)lwgeom1;
2156  }
2157  getPoint2d_p(lwpt->point, 0, &p);
2158 
2159  /* Point in polygon implies zero distance */
2160  if ( lwpoly_covers_point2d(lwpoly, &p) )
2161  {
2162  return 0.0;
2163  }
2164 
2165  /* Not inside, so what's the actual distance? */
2166  for ( i = 0; i < lwpoly->nrings; i++ )
2167  {
2168  double ring_distance = ptarray_distance_spheroid(lwpoly->rings[i], lwpt->point, spheroid, tolerance, check_intersection);
2169  if ( ring_distance < distance )
2170  distance = ring_distance;
2171  if ( distance < tolerance )
2172  return distance;
2173  }
2174  return distance;
2175  }
2176 
2177  /* Line/polygon case, if start point-in-poly, return zero, else return distance. */
2178  if ( ( type1 == POLYGONTYPE && type2 == LINETYPE ) ||
2179  ( type2 == POLYGONTYPE && type1 == LINETYPE ) )
2180  {
2181  POINT2D p;
2182  LWPOLY *lwpoly;
2183  LWLINE *lwline;
2184  double distance = MAXFLOAT;
2185  int i;
2186 
2187  if ( type1 == LINETYPE )
2188  {
2189  lwline = (LWLINE*)lwgeom1;
2190  lwpoly = (LWPOLY*)lwgeom2;
2191  }
2192  else
2193  {
2194  lwline = (LWLINE*)lwgeom2;
2195  lwpoly = (LWPOLY*)lwgeom1;
2196  }
2197  getPoint2d_p(lwline->points, 0, &p);
2198 
2199  LWDEBUG(4, "checking if a point of line is in polygon");
2200 
2201  /* Point in polygon implies zero distance */
2202  if ( lwpoly_covers_point2d(lwpoly, &p) )
2203  return 0.0;
2204 
2205  LWDEBUG(4, "checking ring distances");
2206 
2207  /* Not contained, so what's the actual distance? */
2208  for ( i = 0; i < lwpoly->nrings; i++ )
2209  {
2210  double ring_distance = ptarray_distance_spheroid(lwpoly->rings[i], lwline->points, spheroid, tolerance, check_intersection);
2211  LWDEBUGF(4, "ring[%d] ring_distance = %.8g", i, ring_distance);
2212  if ( ring_distance < distance )
2213  distance = ring_distance;
2214  if ( distance < tolerance )
2215  return distance;
2216  }
2217  LWDEBUGF(4, "all rings checked, returning distance = %.8g", distance);
2218  return distance;
2219 
2220  }
2221 
2222  /* Polygon/polygon case, if start point-in-poly, return zero, else return distance. */
2223  if ( ( type1 == POLYGONTYPE && type2 == POLYGONTYPE ) ||
2224  ( type2 == POLYGONTYPE && type1 == POLYGONTYPE ) )
2225  {
2226  POINT2D p;
2227  LWPOLY *lwpoly1 = (LWPOLY*)lwgeom1;
2228  LWPOLY *lwpoly2 = (LWPOLY*)lwgeom2;
2229  double distance = MAXFLOAT;
2230  int i, j;
2231 
2232  /* Point of 2 in polygon 1 implies zero distance */
2233  getPoint2d_p(lwpoly1->rings[0], 0, &p);
2234  if ( lwpoly_covers_point2d(lwpoly2, &p) )
2235  return 0.0;
2236 
2237  /* Point of 1 in polygon 2 implies zero distance */
2238  getPoint2d_p(lwpoly2->rings[0], 0, &p);
2239  if ( lwpoly_covers_point2d(lwpoly1, &p) )
2240  return 0.0;
2241 
2242  /* Not contained, so what's the actual distance? */
2243  for ( i = 0; i < lwpoly1->nrings; i++ )
2244  {
2245  for ( j = 0; j < lwpoly2->nrings; j++ )
2246  {
2247  double ring_distance = ptarray_distance_spheroid(lwpoly1->rings[i], lwpoly2->rings[j], spheroid, tolerance, check_intersection);
2248  if ( ring_distance < distance )
2249  distance = ring_distance;
2250  if ( distance < tolerance )
2251  return distance;
2252  }
2253  }
2254  return distance;
2255  }
2256 
2257  /* Recurse into collections */
2258  if ( lwtype_is_collection(type1) )
2259  {
2260  int i;
2261  double distance = MAXFLOAT;
2262  LWCOLLECTION *col = (LWCOLLECTION*)lwgeom1;
2263 
2264  for ( i = 0; i < col->ngeoms; i++ )
2265  {
2266  double geom_distance = lwgeom_distance_spheroid(col->geoms[i], lwgeom2, spheroid, tolerance);
2267  if ( geom_distance < distance )
2268  distance = geom_distance;
2269  if ( distance < tolerance )
2270  return distance;
2271  }
2272  return distance;
2273  }
2274 
2275  /* Recurse into collections */
2276  if ( lwtype_is_collection(type2) )
2277  {
2278  int i;
2279  double distance = MAXFLOAT;
2280  LWCOLLECTION *col = (LWCOLLECTION*)lwgeom2;
2281 
2282  for ( i = 0; i < col->ngeoms; i++ )
2283  {
2284  double geom_distance = lwgeom_distance_spheroid(lwgeom1, col->geoms[i], spheroid, tolerance);
2285  if ( geom_distance < distance )
2286  distance = geom_distance;
2287  if ( distance < tolerance )
2288  return distance;
2289  }
2290  return distance;
2291  }
2292 
2293 
2294  lwerror("arguments include unsupported geometry type (%s, %s)", lwtype_name(type1), lwtype_name(type1));
2295  return -1.0;
2296 
2297 }
2298 
2299 
2300 int lwgeom_covers_lwgeom_sphere(const LWGEOM *lwgeom1, const LWGEOM *lwgeom2)
2301 {
2302  int type1, type2;
2303  GBOX gbox1, gbox2;
2304  gbox1.flags = gbox2.flags = 0;
2305 
2306  assert(lwgeom1);
2307  assert(lwgeom2);
2308 
2309  type1 = lwgeom1->type;
2310  type2 = lwgeom2->type;
2311 
2312  /* Currently a restricted implementation */
2313  if ( ! ( (type1 == POLYGONTYPE || type1 == MULTIPOLYGONTYPE || type1 == COLLECTIONTYPE) &&
2314  (type2 == POINTTYPE || type2 == MULTIPOINTTYPE || type2 == COLLECTIONTYPE) ) )
2315  {
2316  lwerror("lwgeom_covers_lwgeom_sphere: only POLYGON covers POINT tests are currently supported");
2317  return LW_FALSE;
2318  }
2319 
2320  /* Make sure we have boxes */
2321  if ( lwgeom1->bbox )
2322  gbox1 = *(lwgeom1->bbox);
2323  else
2324  lwgeom_calculate_gbox_geodetic(lwgeom1, &gbox1);
2325 
2326  /* Make sure we have boxes */
2327  if ( lwgeom2->bbox )
2328  gbox2 = *(lwgeom2->bbox);
2329  else
2330  lwgeom_calculate_gbox_geodetic(lwgeom2, &gbox2);
2331 
2332 
2333  /* Handle the polygon/point case */
2334  if ( type1 == POLYGONTYPE && type2 == POINTTYPE )
2335  {
2336  POINT2D pt_to_test;
2337  getPoint2d_p(((LWPOINT*)lwgeom2)->point, 0, &pt_to_test);
2338  return lwpoly_covers_point2d((LWPOLY*)lwgeom1, &pt_to_test);
2339  }
2340 
2341  /* If any of the first argument parts covers the second argument, it's true */
2342  if ( lwtype_is_collection( type1 ) )
2343  {
2344  int i;
2345  LWCOLLECTION *col = (LWCOLLECTION*)lwgeom1;
2346 
2347  for ( i = 0; i < col->ngeoms; i++ )
2348  {
2349  if ( lwgeom_covers_lwgeom_sphere(col->geoms[i], lwgeom2) )
2350  {
2351  return LW_TRUE;
2352  }
2353  }
2354  return LW_FALSE;
2355  }
2356 
2357  /* Only if all of the second arguments are covered by the first argument is the condition true */
2358  if ( lwtype_is_collection( type2 ) )
2359  {
2360  int i;
2361  LWCOLLECTION *col = (LWCOLLECTION*)lwgeom2;
2362 
2363  for ( i = 0; i < col->ngeoms; i++ )
2364  {
2365  if ( ! lwgeom_covers_lwgeom_sphere(lwgeom1, col->geoms[i]) )
2366  {
2367  return LW_FALSE;
2368  }
2369  }
2370  return LW_TRUE;
2371  }
2372 
2373  /* Don't get here */
2374  lwerror("lwgeom_covers_lwgeom_sphere: reached end of function without resolution");
2375  return LW_FALSE;
2376 
2377 }
2378 
2384 int lwpoly_covers_point2d(const LWPOLY *poly, const POINT2D *pt_to_test)
2385 {
2386  int i;
2387  int in_hole_count = 0;
2388  POINT3D p;
2389  GEOGRAPHIC_POINT gpt_to_test;
2390  POINT2D pt_outside;
2391  GBOX gbox;
2392  gbox.flags = 0;
2393 
2394  /* Nulls and empties don't contain anything! */
2395  if ( ! poly || lwgeom_is_empty((LWGEOM*)poly) )
2396  {
2397  LWDEBUG(4,"returning false, geometry is empty or null");
2398  return LW_FALSE;
2399  }
2400 
2401  /* Make sure we have boxes */
2402  if ( poly->bbox )
2403  gbox = *(poly->bbox);
2404  else
2405  lwgeom_calculate_gbox_geodetic((LWGEOM*)poly, &gbox);
2406 
2407  /* Point not in box? Done! */
2408  geographic_point_init(pt_to_test->x, pt_to_test->y, &gpt_to_test);
2409  geog2cart(&gpt_to_test, &p);
2410  if ( ! gbox_contains_point3d(&gbox, &p) )
2411  {
2412  LWDEBUG(4, "the point is not in the box!");
2413  return LW_FALSE;
2414  }
2415 
2416  /* Calculate our outside point from the gbox */
2417  gbox_pt_outside(&gbox, &pt_outside);
2418 
2419  LWDEBUGF(4, "pt_outside POINT(%.18g %.18g)", pt_outside.x, pt_outside.y);
2420  LWDEBUGF(4, "pt_to_test POINT(%.18g %.18g)", pt_to_test->x, pt_to_test->y);
2421  LWDEBUGF(4, "polygon %s", lwgeom_to_ewkt((LWGEOM*)poly));
2422  LWDEBUGF(4, "gbox %s", gbox_to_string(&gbox));
2423 
2424  /* Not in outer ring? We're done! */
2425  if ( ! ptarray_contains_point_sphere(poly->rings[0], &pt_outside, pt_to_test) )
2426  {
2427  LWDEBUG(4,"returning false, point is outside ring");
2428  return LW_FALSE;
2429  }
2430 
2431  LWDEBUGF(4, "testing %d rings", poly->nrings);
2432 
2433  /* But maybe point is in a hole... */
2434  for ( i = 1; i < poly->nrings; i++ )
2435  {
2436  LWDEBUGF(4, "ring test loop %d", i);
2437  /* Count up hole containment. Odd => outside boundary. */
2438  if ( ptarray_contains_point_sphere(poly->rings[i], &pt_outside, pt_to_test) )
2439  in_hole_count++;
2440  }
2441 
2442  LWDEBUGF(4, "in_hole_count == %d", in_hole_count);
2443 
2444  if ( in_hole_count % 2 )
2445  {
2446  LWDEBUG(4,"returning false, inner ring containment count is odd");
2447  return LW_FALSE;
2448  }
2449 
2450  LWDEBUG(4,"returning true, inner ring containment count is even");
2451  return LW_TRUE;
2452 }
2453 
2454 
2459 int getPoint2d_p_ro(const POINTARRAY *pa, int n, POINT2D **point)
2460 {
2461  uint8_t *pa_ptr = NULL;
2462  assert(pa);
2463  assert(n >= 0);
2464  assert(n < pa->npoints);
2465 
2466  pa_ptr = getPoint_internal(pa, n);
2467  /* printf( "pa_ptr[0]: %g\n", *((double*)pa_ptr)); */
2468  *point = (POINT2D*)pa_ptr;
2469 
2470  return LW_SUCCESS;
2471 }
2472 
2474 {
2475  int i;
2476  int first = LW_TRUE;
2477  const POINT2D *p;
2478  POINT3D A1, A2;
2479  GBOX edge_gbox;
2480 
2481  assert(gbox);
2482  assert(pa);
2483 
2484  gbox_init(&edge_gbox);
2485  edge_gbox.flags = gbox->flags;
2486 
2487  if ( pa->npoints == 0 ) return LW_FAILURE;
2488 
2489  if ( pa->npoints == 1 )
2490  {
2491  p = getPoint2d_cp(pa, 0);
2492  ll2cart(p, &A1);
2493  gbox->xmin = gbox->xmax = A1.x;
2494  gbox->ymin = gbox->ymax = A1.y;
2495  gbox->zmin = gbox->zmax = A1.z;
2496  return LW_SUCCESS;
2497  }
2498 
2499  p = getPoint2d_cp(pa, 0);
2500  ll2cart(p, &A1);
2501 
2502  for ( i = 1; i < pa->npoints; i++ )
2503  {
2504 
2505  p = getPoint2d_cp(pa, i);
2506  ll2cart(p, &A2);
2507 
2508  edge_calculate_gbox(&A1, &A2, &edge_gbox);
2509 
2510  /* Initialize the box */
2511  if ( first )
2512  {
2513  gbox_duplicate(&edge_gbox, gbox);
2514  first = LW_FALSE;
2515  }
2516  /* Expand the box where necessary */
2517  else
2518  {
2519  gbox_merge(&edge_gbox, gbox);
2520  }
2521 
2522  A1 = A2;
2523  }
2524 
2525  return LW_SUCCESS;
2526 }
2527 
2528 static int lwpoint_calculate_gbox_geodetic(const LWPOINT *point, GBOX *gbox)
2529 {
2530  assert(point);
2531  return ptarray_calculate_gbox_geodetic(point->point, gbox);
2532 }
2533 
2534 static int lwline_calculate_gbox_geodetic(const LWLINE *line, GBOX *gbox)
2535 {
2536  assert(line);
2537  return ptarray_calculate_gbox_geodetic(line->points, gbox);
2538 }
2539 
2540 static int lwpolygon_calculate_gbox_geodetic(const LWPOLY *poly, GBOX *gbox)
2541 {
2542  GBOX ringbox;
2543  int i;
2544  int first = LW_TRUE;
2545  assert(poly);
2546  if ( poly->nrings == 0 )
2547  return LW_FAILURE;
2548  ringbox.flags = gbox->flags;
2549  for ( i = 0; i < poly->nrings; i++ )
2550  {
2551  if ( ptarray_calculate_gbox_geodetic(poly->rings[i], &ringbox) == LW_FAILURE )
2552  return LW_FAILURE;
2553  if ( first )
2554  {
2555  gbox_duplicate(&ringbox, gbox);
2556  first = LW_FALSE;
2557  }
2558  else
2559  {
2560  gbox_merge(&ringbox, gbox);
2561  }
2562  }
2563 
2564  /* If the box wraps a poly, push that axis to the absolute min/max as appropriate */
2565  gbox_check_poles(gbox);
2566 
2567  return LW_SUCCESS;
2568 }
2569 
2570 static int lwtriangle_calculate_gbox_geodetic(const LWTRIANGLE *triangle, GBOX *gbox)
2571 {
2572  assert(triangle);
2573  return ptarray_calculate_gbox_geodetic(triangle->points, gbox);
2574 }
2575 
2576 
2578 {
2579  GBOX subbox;
2580  int i;
2581  int result = LW_FAILURE;
2582  int first = LW_TRUE;
2583  assert(coll);
2584  if ( coll->ngeoms == 0 )
2585  return LW_FAILURE;
2586 
2587  subbox.flags = gbox->flags;
2588 
2589  for ( i = 0; i < coll->ngeoms; i++ )
2590  {
2591  if ( lwgeom_calculate_gbox_geodetic((LWGEOM*)(coll->geoms[i]), &subbox) == LW_SUCCESS )
2592  {
2593  /* Keep a copy of the sub-bounding box for later */
2594  if ( coll->geoms[i]->bbox )
2595  lwfree(coll->geoms[i]->bbox);
2596  coll->geoms[i]->bbox = gbox_copy(&subbox);
2597  if ( first )
2598  {
2599  gbox_duplicate(&subbox, gbox);
2600  first = LW_FALSE;
2601  }
2602  else
2603  {
2604  gbox_merge(&subbox, gbox);
2605  }
2606  result = LW_SUCCESS;
2607  }
2608  }
2609  return result;
2610 }
2611 
2613 {
2614  int result = LW_FAILURE;
2615  LWDEBUGF(4, "got type %d", geom->type);
2616 
2617  /* Add a geodetic flag to the incoming gbox */
2618  gbox->flags = gflags(FLAGS_GET_Z(geom->flags),FLAGS_GET_M(geom->flags),1);
2619 
2620  switch (geom->type)
2621  {
2622  case POINTTYPE:
2623  result = lwpoint_calculate_gbox_geodetic((LWPOINT*)geom, gbox);
2624  break;
2625  case LINETYPE:
2626  result = lwline_calculate_gbox_geodetic((LWLINE *)geom, gbox);
2627  break;
2628  case POLYGONTYPE:
2629  result = lwpolygon_calculate_gbox_geodetic((LWPOLY *)geom, gbox);
2630  break;
2631  case TRIANGLETYPE:
2632  result = lwtriangle_calculate_gbox_geodetic((LWTRIANGLE *)geom, gbox);
2633  break;
2634  case MULTIPOINTTYPE:
2635  case MULTILINETYPE:
2636  case MULTIPOLYGONTYPE:
2637  case POLYHEDRALSURFACETYPE:
2638  case TINTYPE:
2639  case COLLECTIONTYPE:
2640  result = lwcollection_calculate_gbox_geodetic((LWCOLLECTION *)geom, gbox);
2641  break;
2642  default:
2643  lwerror("lwgeom_calculate_gbox_geodetic: unsupported input geometry type: %d - %s",
2644  geom->type, lwtype_name(geom->type));
2645  break;
2646  }
2647  return result;
2648 }
2649 
2650 
2651 
2652 static int ptarray_check_geodetic(const POINTARRAY *pa)
2653 {
2654  int t;
2655  POINT2D pt;
2656 
2657  assert(pa);
2658 
2659  for (t=0; t<pa->npoints; t++)
2660  {
2661  getPoint2d_p(pa, t, &pt);
2662  /* printf( "%d (%g, %g)\n", t, pt.x, pt.y); */
2663  if ( pt.x < -180.0 || pt.y < -90.0 || pt.x > 180.0 || pt.y > 90.0 )
2664  return LW_FALSE;
2665  }
2666 
2667  return LW_TRUE;
2668 }
2669 
2670 static int lwpoint_check_geodetic(const LWPOINT *point)
2671 {
2672  assert(point);
2673  return ptarray_check_geodetic(point->point);
2674 }
2675 
2676 static int lwline_check_geodetic(const LWLINE *line)
2677 {
2678  assert(line);
2679  return ptarray_check_geodetic(line->points);
2680 }
2681 
2682 static int lwpoly_check_geodetic(const LWPOLY *poly)
2683 {
2684  int i = 0;
2685  assert(poly);
2686 
2687  for ( i = 0; i < poly->nrings; i++ )
2688  {
2689  if ( ptarray_check_geodetic(poly->rings[i]) == LW_FALSE )
2690  return LW_FALSE;
2691  }
2692  return LW_TRUE;
2693 }
2694 
2695 static int lwtriangle_check_geodetic(const LWTRIANGLE *triangle)
2696 {
2697  assert(triangle);
2698  return ptarray_check_geodetic(triangle->points);
2699 }
2700 
2701 
2703 {
2704  int i = 0;
2705  assert(col);
2706 
2707  for ( i = 0; i < col->ngeoms; i++ )
2708  {
2709  if ( lwgeom_check_geodetic(col->geoms[i]) == LW_FALSE )
2710  return LW_FALSE;
2711  }
2712  return LW_TRUE;
2713 }
2714 
2716 {
2717  if ( lwgeom_is_empty(geom) )
2718  return LW_TRUE;
2719 
2720  switch (geom->type)
2721  {
2722  case POINTTYPE:
2723  return lwpoint_check_geodetic((LWPOINT *)geom);
2724  case LINETYPE:
2725  return lwline_check_geodetic((LWLINE *)geom);
2726  case POLYGONTYPE:
2727  return lwpoly_check_geodetic((LWPOLY *)geom);
2728  case TRIANGLETYPE:
2729  return lwtriangle_check_geodetic((LWTRIANGLE *)geom);
2730  case MULTIPOINTTYPE:
2731  case MULTILINETYPE:
2732  case MULTIPOLYGONTYPE:
2733  case POLYHEDRALSURFACETYPE:
2734  case TINTYPE:
2735  case COLLECTIONTYPE:
2736  return lwcollection_check_geodetic((LWCOLLECTION *)geom);
2737  default:
2738  lwerror("lwgeom_check_geodetic: unsupported input geometry type: %d - %s",
2739  geom->type, lwtype_name(geom->type));
2740  }
2741  return LW_FALSE;
2742 }
2743 
2745 {
2746  int t;
2747  int changed = LW_FALSE;
2748  POINT4D pt;
2749 
2750  assert(pa);
2751 
2752  for ( t=0; t < pa->npoints; t++ )
2753  {
2754  getPoint4d_p(pa, t, &pt);
2755  if ( pt.x < -180.0 || pt.x > 180.0 || pt.y < -90.0 || pt.y > 90.0 )
2756  {
2757  pt.x = longitude_degrees_normalize(pt.x);
2758  pt.y = latitude_degrees_normalize(pt.y);
2759  ptarray_set_point4d(pa, t, &pt);
2760  changed = LW_TRUE;
2761  }
2762  }
2763  return changed;
2764 }
2765 
2767 {
2768  assert(point);
2769  return ptarray_force_geodetic(point->point);
2770 }
2771 
2773 {
2774  assert(line);
2775  return ptarray_force_geodetic(line->points);
2776 }
2777 
2779 {
2780  int i = 0;
2781  int changed = LW_FALSE;
2782  assert(poly);
2783 
2784  for ( i = 0; i < poly->nrings; i++ )
2785  {
2786  if ( ptarray_force_geodetic(poly->rings[i]) == LW_TRUE )
2787  changed = LW_TRUE;
2788  }
2789  return changed;
2790 }
2791 
2793 {
2794  int i = 0;
2795  int changed = LW_FALSE;
2796  assert(col);
2797 
2798  for ( i = 0; i < col->ngeoms; i++ )
2799  {
2800  if ( lwgeom_force_geodetic(col->geoms[i]) == LW_TRUE )
2801  changed = LW_TRUE;
2802  }
2803  return changed;
2804 }
2805 
2807 {
2808  switch ( lwgeom_get_type(geom) )
2809  {
2810  case POINTTYPE:
2811  return lwpoint_force_geodetic((LWPOINT *)geom);
2812  case LINETYPE:
2813  return lwline_force_geodetic((LWLINE *)geom);
2814  case POLYGONTYPE:
2815  return lwpoly_force_geodetic((LWPOLY *)geom);
2816  case MULTIPOINTTYPE:
2817  case MULTILINETYPE:
2818  case MULTIPOLYGONTYPE:
2819  case COLLECTIONTYPE:
2820  return lwcollection_force_geodetic((LWCOLLECTION *)geom);
2821  default:
2822  lwerror("unsupported input geometry type: %d", lwgeom_get_type(geom));
2823  }
2824  return LW_FALSE;
2825 }
2826 
2827 
2829 {
2830  GEOGRAPHIC_POINT a, b;
2831  double za = 0.0, zb = 0.0;
2832  POINT4D p;
2833  int i;
2834  int hasz = LW_FALSE;
2835  double length = 0.0;
2836  double seglength = 0.0;
2837 
2838  /* Return zero on non-sensical inputs */
2839  if ( ! pa || pa->npoints < 2 )
2840  return 0.0;
2841 
2842  /* See if we have a third dimension */
2843  hasz = FLAGS_GET_Z(pa->flags);
2844 
2845  /* Initialize first point */
2846  getPoint4d_p(pa, 0, &p);
2847  geographic_point_init(p.x, p.y, &a);
2848  if ( hasz )
2849  za = p.z;
2850 
2851  /* Loop and sum the length for each segment */
2852  for ( i = 1; i < pa->npoints; i++ )
2853  {
2854  seglength = 0.0;
2855  getPoint4d_p(pa, i, &p);
2856  geographic_point_init(p.x, p.y, &b);
2857  if ( hasz )
2858  zb = p.z;
2859 
2860  /* Special sphere case */
2861  if ( s->a == s->b )
2862  seglength = s->radius * sphere_distance(&a, &b);
2863  /* Spheroid case */
2864  else
2865  seglength = spheroid_distance(&a, &b, s);
2866 
2867  /* Add in the vertical displacement if we're in 3D */
2868  if ( hasz )
2869  seglength = sqrt( (zb-za)*(zb-za) + seglength*seglength );
2870 
2871  /* Add this segment length to the total */
2872  length += seglength;
2873 
2874  /* B gets incremented in the next loop, so we save the value here */
2875  a = b;
2876  za = zb;
2877  }
2878  return length;
2879 }
2880 
2881 double lwgeom_length_spheroid(const LWGEOM *geom, const SPHEROID *s)
2882 {
2883  int type;
2884  int i = 0;
2885  double length = 0.0;
2886 
2887  assert(geom);
2888 
2889  /* No area in nothing */
2890  if ( lwgeom_is_empty(geom) )
2891  return 0.0;
2892 
2893  type = geom->type;
2894 
2895  if ( type == POINTTYPE || type == MULTIPOINTTYPE )
2896  return 0.0;
2897 
2898  if ( type == LINETYPE )
2899  return ptarray_length_spheroid(((LWLINE*)geom)->points, s);
2900 
2901  if ( type == POLYGONTYPE )
2902  {
2903  LWPOLY *poly = (LWPOLY*)geom;
2904  for ( i = 0; i < poly->nrings; i++ )
2905  {
2906  length += ptarray_length_spheroid(poly->rings[i], s);
2907  }
2908  return length;
2909  }
2910 
2911  if ( type == TRIANGLETYPE )
2912  return ptarray_length_spheroid(((LWTRIANGLE*)geom)->points, s);
2913 
2914  if ( lwtype_is_collection( type ) )
2915  {
2916  LWCOLLECTION *col = (LWCOLLECTION*)geom;
2917 
2918  for ( i = 0; i < col->ngeoms; i++ )
2919  {
2920  length += lwgeom_length_spheroid(col->geoms[i], s);
2921  }
2922  return length;
2923  }
2924 
2925  lwerror("unsupported type passed to lwgeom_length_sphere");
2926  return 0.0;
2927 }
2928 
2935 static int
2937 {
2938 
2939  int i;
2940  POINT4D p;
2941  int altered = LW_FALSE;
2942  int rv = LW_FALSE;
2943  static double tolerance = 1e-10;
2944 
2945  if ( ! pa )
2946  lwerror("ptarray_nudge_geodetic called with null input");
2947 
2948  for(i = 0; i < pa->npoints; i++ )
2949  {
2950  getPoint4d_p(pa, i, &p);
2951  if ( p.x < -180.0 && (-180.0 - p.x < tolerance) )
2952  {
2953  p.x = -180.0;
2954  altered = LW_TRUE;
2955  }
2956  if ( p.x > 180.0 && (p.x - 180.0 < tolerance) )
2957  {
2958  p.x = 180.0;
2959  altered = LW_TRUE;
2960  }
2961  if ( p.y < -90.0 && (-90.0 - p.y < tolerance) )
2962  {
2963  p.y = -90.0;
2964  altered = LW_TRUE;
2965  }
2966  if ( p.y > 90.0 && (p.y - 90.0 < tolerance) )
2967  {
2968  p.y = 90.0;
2969  altered = LW_TRUE;
2970  }
2971  if ( altered == LW_TRUE )
2972  {
2973  ptarray_set_point4d(pa, i, &p);
2974  altered = LW_FALSE;
2975  rv = LW_TRUE;
2976  }
2977  }
2978  return rv;
2979 }
2980 
2987 int
2989 {
2990  int type;
2991  int i = 0;
2992  int rv = LW_FALSE;
2993 
2994  assert(geom);
2995 
2996  /* No points in nothing */
2997  if ( lwgeom_is_empty(geom) )
2998  return LW_FALSE;
2999 
3000  type = geom->type;
3001 
3002  if ( type == POINTTYPE )
3003  return ptarray_nudge_geodetic(((LWPOINT*)geom)->point);
3004 
3005  if ( type == LINETYPE )
3006  return ptarray_nudge_geodetic(((LWLINE*)geom)->points);
3007 
3008  if ( type == POLYGONTYPE )
3009  {
3010  LWPOLY *poly = (LWPOLY*)geom;
3011  for ( i = 0; i < poly->nrings; i++ )
3012  {
3013  int n = ptarray_nudge_geodetic(poly->rings[i]);
3014  rv = (rv == LW_TRUE ? rv : n);
3015  }
3016  return rv;
3017  }
3018 
3019  if ( type == TRIANGLETYPE )
3020  return ptarray_nudge_geodetic(((LWTRIANGLE*)geom)->points);
3021 
3022  if ( lwtype_is_collection( type ) )
3023  {
3024  LWCOLLECTION *col = (LWCOLLECTION*)geom;
3025 
3026  for ( i = 0; i < col->ngeoms; i++ )
3027  {
3028  int n = lwgeom_nudge_geodetic(col->geoms[i]);
3029  rv = (rv == LW_TRUE ? rv : n);
3030  }
3031  return rv;
3032  }
3033 
3034  lwerror("unsupported type (%s) passed to lwgeom_nudge_geodetic", lwtype_name(type));
3035  return rv;
3036 }
3037 
3038 
3042 static int
3043 point_in_cone(const POINT3D *A1, const POINT3D *A2, const POINT3D *P)
3044 {
3045  POINT3D AC; /* Center point of A1/A2 */
3046  double min_similarity, similarity;
3047 
3048  /* The normalized sum bisects the angle between start and end. */
3049  vector_sum(A1, A2, &AC);
3050  normalize(&AC);
3051 
3052  /* The projection of start onto the center defines the minimum similarity */
3053  min_similarity = dot_product(A1, &AC);
3054 
3055  /* The projection of candidate p onto the center */
3056  similarity = dot_product(P, &AC);
3057 
3058  /* If the point is more similar than the end, the point is in the cone */
3059  if ( similarity > min_similarity || fabs(similarity - min_similarity) < 2e-16 )
3060  {
3061  return LW_TRUE;
3062  }
3063  return LW_FALSE;
3064 }
3065 
3066 
3070 static int
3071 point3d_equals(const POINT3D *p1, const POINT3D *p2)
3072 {
3073  return FP_EQUALS(p1->x, p2->x) && FP_EQUALS(p1->y, p2->y) && FP_EQUALS(p1->z, p2->z);
3074 }
3075 
3080 static int
3081 dot_product_side(const POINT3D *p, const POINT3D *q)
3082 {
3083  double dp = dot_product(p, q);
3084 
3085  if ( FP_IS_ZERO(dp) )
3086  return 0;
3087 
3088  return dp < 0.0 ? -1 : 1;
3089 }
3090 
3095 int
3096 edge_intersects(const POINT3D *A1, const POINT3D *A2, const POINT3D *B1, const POINT3D *B2)
3097 {
3098  POINT3D AN, BN, VN; /* Normals to plane A and plane B */
3099  double ab_dot;
3100  int a1_side, a2_side, b1_side, b2_side;
3101  int rv = PIR_NO_INTERACT;
3102 
3103  /* Normals to the A-plane and B-plane */
3104  unit_normal(A1, A2, &AN);
3105  unit_normal(B1, B2, &BN);
3106 
3107  /* Are A-plane and B-plane basically the same? */
3108  ab_dot = dot_product(&AN, &BN);
3109  if ( FP_EQUALS(fabs(ab_dot), 1.0) )
3110  {
3111  /* Co-linear case */
3112  if ( point_in_cone(A1, A2, B1) || point_in_cone(A1, A2, B2) ||
3113  point_in_cone(B1, B2, A1) || point_in_cone(B1, B2, A2) )
3114  {
3115  rv |= PIR_INTERSECTS;
3116  rv |= PIR_COLINEAR;
3117  }
3118  return rv;
3119  }
3120 
3121  /* What side of plane-A and plane-B do the end points */
3122  /* of A and B fall? */
3123  a1_side = dot_product_side(&BN, A1);
3124  a2_side = dot_product_side(&BN, A2);
3125  b1_side = dot_product_side(&AN, B1);
3126  b2_side = dot_product_side(&AN, B2);
3127 
3128  /* Both ends of A on the same side of plane B. */
3129  if ( a1_side == a2_side && a1_side != 0 )
3130  {
3131  /* No intersection. */
3132  return PIR_NO_INTERACT;
3133  }
3134 
3135  /* Both ends of B on the same side of plane A. */
3136  if ( b1_side == b2_side && b1_side != 0 )
3137  {
3138  /* No intersection. */
3139  return PIR_NO_INTERACT;
3140  }
3141 
3142  /* A straddles B and B straddles A, so... */
3143  if ( a1_side != a2_side && (a1_side + a2_side) == 0 &&
3144  b1_side != b2_side && (b1_side + b2_side) == 0 )
3145  {
3146  /* Have to check if intersection point is inside both arcs */
3147  unit_normal(&AN, &BN, &VN);
3148  if ( point_in_cone(A1, A2, &VN) && point_in_cone(B1, B2, &VN) )
3149  {
3150  return PIR_INTERSECTS;
3151  }
3152 
3153  /* Have to check if intersection point is inside both arcs */
3154  vector_scale(&VN, -1);
3155  if ( point_in_cone(A1, A2, &VN) && point_in_cone(B1, B2, &VN) )
3156  {
3157  return PIR_INTERSECTS;
3158  }
3159 
3160  return PIR_NO_INTERACT;
3161  }
3162 
3163  /* The rest are all intersects variants... */
3164  rv |= PIR_INTERSECTS;
3165 
3166  /* A touches B */
3167  if ( a1_side == 0 )
3168  {
3169  /* Touches at A1, A2 is on what side? */
3170  rv |= (a2_side < 0 ? PIR_A_TOUCH_RIGHT : PIR_A_TOUCH_LEFT);
3171  }
3172  else if ( a2_side == 0 )
3173  {
3174  /* Touches at A2, A1 is on what side? */
3175  rv |= (a1_side < 0 ? PIR_A_TOUCH_RIGHT : PIR_A_TOUCH_LEFT);
3176  }
3177 
3178  /* B touches A */
3179  if ( b1_side == 0 )
3180  {
3181  /* Touches at B1, B2 is on what side? */
3182  rv |= (b2_side < 0 ? PIR_B_TOUCH_RIGHT : PIR_B_TOUCH_LEFT);
3183  }
3184  else if ( b2_side == 0 )
3185  {
3186  /* Touches at B2, B1 is on what side? */
3187  rv |= (b1_side < 0 ? PIR_B_TOUCH_RIGHT : PIR_B_TOUCH_LEFT);
3188  }
3189 
3190  return rv;
3191 }
3192 
3201 int ptarray_contains_point_sphere(const POINTARRAY *pa, const POINT2D *pt_outside, const POINT2D *pt_to_test)
3202 {
3203  POINT3D S1, S2; /* Stab line end points */
3204  POINT3D E1, E2; /* Edge end points (3-space) */
3205  POINT2D p; /* Edge end points (lon/lat) */
3206  int count = 0, i, inter;
3207 
3208  /* Null input, not enough points for a ring? You ain't closed! */
3209  if ( ! pa || pa->npoints < 4 )
3210  return LW_FALSE;
3211 
3212  /* Set up our stab line */
3213  ll2cart(pt_to_test, &S1);
3214  ll2cart(pt_outside, &S2);
3215 
3216  /* Initialize first point */
3217  getPoint2d_p(pa, 0, &p);
3218  ll2cart(&p, &E1);
3219 
3220  /* Walk every edge and see if the stab line hits it */
3221  for ( i = 1; i < pa->npoints; i++ )
3222  {
3223  LWDEBUGF(4, "testing edge (%d)", i);
3224  LWDEBUGF(4, " start point == POINT(%.12g %.12g)", p.x, p.y);
3225 
3226  /* Read next point. */
3227  getPoint2d_p(pa, i, &p);
3228  ll2cart(&p, &E2);
3229 
3230  /* Skip over too-short edges. */
3231  if ( point3d_equals(&E1, &E2) )
3232  {
3233  continue;
3234  }
3235 
3236  /* Our test point is on an edge end! Point is "in ring" by our definition */
3237  if ( point3d_equals(&S1, &E1) )
3238  {
3239  return LW_TRUE;
3240  }
3241 
3242  /* Calculate relationship between stab line and edge */
3243  inter = edge_intersects(&S1, &S2, &E1, &E2);
3244 
3245  /* We have some kind of interaction... */
3246  if ( inter & PIR_INTERSECTS )
3247  {
3248  /* If the stabline is touching the edge, that implies the test point */
3249  /* is on the edge, so we're done, the point is in (on) the ring. */
3250  if ( (inter & PIR_A_TOUCH_RIGHT) || (inter & PIR_A_TOUCH_LEFT) )
3251  {
3252  return LW_TRUE;
3253  }
3254 
3255  /* It's a touching interaction, disregard all the left-side ones. */
3256  /* It's a co-linear intersection, ignore those. */
3257  if ( inter & PIR_B_TOUCH_RIGHT || inter & PIR_COLINEAR )
3258  {
3259  /* Do nothing, to avoid double counts. */
3260  LWDEBUGF(4," edge (%d) crossed, disregarding to avoid double count", i, count);
3261  }
3262  else
3263  {
3264  /* Increment crossingn count. */
3265  count++;
3266  LWDEBUGF(4," edge (%d) crossed, count == %d", i, count);
3267  }
3268  }
3269  else
3270  {
3271  LWDEBUGF(4," edge (%d) did not cross", i);
3272  }
3273 
3274  /* Increment to next edge */
3275  E1 = E2;
3276  }
3277 
3278  LWDEBUGF(4,"final count == %d", count);
3279 
3280  /* An odd number of crossings implies containment! */
3281  if ( count % 2 )
3282  {
3283  return LW_TRUE;
3284  }
3285 
3286  return LW_FALSE;
3287 }
void ptarray_set_point4d(POINTARRAY *pa, int n, const POINT4D *p4d)
Definition: lwgeom_api.c:501
double x
Definition: liblwgeom.h:308
#define LINETYPE
Definition: liblwgeom.h:61
GBOX * gbox_copy(const GBOX *box)
Return a copy of the GBOX, based on dimensionality of flags.
Definition: g_box.c:362
int gbox_centroid(const GBOX *gbox, POINT2D *out)
Computes the average(ish) center of the box and returns success.
Definition: lwgeodetic.c:244
void robust_cross_product(const GEOGRAPHIC_POINT *p, const GEOGRAPHIC_POINT *q, POINT3D *a)
Computes the cross product of two vectors using their lat, lng representations.
Definition: lwgeodetic.c:583
static int lwtriangle_calculate_gbox_geodetic(const LWTRIANGLE *triangle, GBOX *gbox)
Definition: lwgeodetic.c:2570
static int ptarray_nudge_geodetic(POINTARRAY *pa)
When features are snapped or sometimes they are just this way, they are very close to the geodetic bo...
Definition: lwgeodetic.c:2936
double sphere_distance(const GEOGRAPHIC_POINT *s, const GEOGRAPHIC_POINT *e)
Given two points on a unit sphere, calculate their distance apart in radians.
Definition: lwgeodetic.c:897
static int lwline_force_geodetic(LWLINE *line)
Definition: lwgeodetic.c:2772
POINTARRAY * ptarray_construct(char hasz, char hasm, uint32_t npoints)
Construct an empty pointarray, allocating storage and setting the npoints, but not filling in any inf...
Definition: ptarray.c:49
double longitude_radians_normalize(double lon)
Convert a longitude to the range of -PI,PI.
Definition: lwgeodetic.c:27
GBOX * bbox
Definition: liblwgeom.h:354
POINTARRAY * points
Definition: liblwgeom.h:389
int sphere_project(const GEOGRAPHIC_POINT *r, double distance, double azimuth, GEOGRAPHIC_POINT *n)
Given a starting location r, a distance and an azimuth to the new point, compute the location of the ...
Definition: lwgeodetic.c:1261
double m
Definition: liblwgeom.h:308
double edge_distance_to_edge(const GEOGRAPHIC_EDGE *e1, const GEOGRAPHIC_EDGE *e2, GEOGRAPHIC_POINT *closest1, GEOGRAPHIC_POINT *closest2)
Calculate the distance between two edges.
Definition: lwgeodetic.c:1216
Two-point great circle segment from a to b.
Definition: lwgeodetic.h:42
double gbox_angular_height(const GBOX *gbox)
Returns the angular height (latitudinal span) of the box in radians.
Definition: lwgeodetic.c:165
#define PIR_B_TOUCH_RIGHT
Definition: lwgeodetic.h:77
double spheroid_direction(const GEOGRAPHIC_POINT *r, const GEOGRAPHIC_POINT *s, const SPHEROID *spheroid)
Computes the direction of the geodesic joining two points on the spheroid.
Definition: lwspheroid.c:155
void normalize(POINT3D *p)
Normalize to a unit vector.
Definition: lwgeodetic.c:564
static double ptarray_distance_spheroid(const POINTARRAY *pa1, const POINTARRAY *pa2, const SPHEROID *s, double tolerance, int check_intersection)
Definition: lwgeodetic.c:1730
double lwgeom_distance_spheroid(const LWGEOM *lwgeom1, const LWGEOM *lwgeom2, const SPHEROID *spheroid, double tolerance)
Calculate the distance between two LWGEOMs, using the coordinates are longitude and latitude...
Definition: lwgeodetic.c:2078
char * r
Definition: cu_in_wkt.c:25
int getPoint2d_p_ro(const POINTARRAY *pa, int n, POINT2D **point)
This function can only be used on LWGEOM that is built on top of GSERIALIZED, otherwise alignment err...
Definition: lwgeodetic.c:2459
void gbox_duplicate(const GBOX *original, GBOX *duplicate)
Copy the values of original GBOX into duplicate.
Definition: g_box.c:369
void lwfree(void *mem)
Definition: lwutil.c:190
int edge_contains_point(const GEOGRAPHIC_EDGE *e, const GEOGRAPHIC_POINT *p)
Returns true if the point p is on the minor edge defined by the end points of e.
Definition: lwgeodetic.c:981
double y
Definition: liblwgeom.h:296
char * gbox_to_string(const GBOX *gbox)
Allocate a string representation of the GBOX, based on dimensionality of flags.
Definition: g_box.c:328
int npoints
Definition: liblwgeom.h:327
static int point3d_equals(const POINT3D *p1, const POINT3D *p2)
Utility function for ptarray_contains_point_sphere()
Definition: lwgeodetic.c:3071
int crosses_dateline(const GEOGRAPHIC_POINT *s, const GEOGRAPHIC_POINT *e)
Definition: lwgeodetic.c:615
#define PIR_A_TOUCH_RIGHT
Definition: lwgeodetic.h:75
uint32_t lwgeom_get_type(const LWGEOM *geom)
Return LWTYPE number.
Definition: lwgeom.c:785
int lwgeom_covers_lwgeom_sphere(const LWGEOM *lwgeom1, const LWGEOM *lwgeom2)
Calculate covers predicate for two lwgeoms on the sphere.
Definition: lwgeodetic.c:2300
double z_to_latitude(double z, int top)
Used in great circle to compute the pole of the great circle.
Definition: lwgeodetic.c:996
#define POLYGONTYPE
Definition: liblwgeom.h:62
Datum area(PG_FUNCTION_ARGS)
char * lwgeom_to_ewkt(const LWGEOM *lwgeom)
Return an alloced string.
Definition: lwgeom.c:425
double b
Definition: liblwgeom.h:270
uint8_t flags
Definition: liblwgeom.h:353
double xmax
Definition: liblwgeom.h:249
#define NAN
Definition: lwgeodetic.h:23
POINTARRAY * ptarray_construct_empty(char hasz, char hasm, uint32_t maxpoints)
Create a new POINTARRAY with no points.
Definition: ptarray.c:57
#define MULTIPOINTTYPE
Definition: liblwgeom.h:63
int gbox_contains_point3d(const GBOX *gbox, const POINT3D *pt)
Return true if the point is inside the gbox.
Definition: g_box.c:205
#define LW_SUCCESS
Definition: liblwgeom.h:55
void cart2geog(const POINT3D *p, GEOGRAPHIC_POINT *g)
Convert cartesion coordinates on unit sphere to spherical coordinates.
Definition: lwgeodetic.c:365
#define signum(a)
Ape a java function.
Definition: lwgeodetic.h:66
double radius
Definition: liblwgeom.h:274
int edge_calculate_gbox_slow(const GEOGRAPHIC_EDGE *e, GBOX *gbox)
Definition: lwgeodetic.c:1291
#define LWDEBUG(level, msg)
Definition: lwgeom_log.h:50
#define TRIANGLETYPE
Definition: liblwgeom.h:73
static int dot_product_side(const POINT3D *p, const POINT3D *q)
Utility function for edge_intersects(), signum with a tolerance in determining if the value is zero...
Definition: lwgeodetic.c:3081
double x
Definition: liblwgeom.h:296
#define POLYHEDRALSURFACETYPE
Definition: liblwgeom.h:72
void point_shift(GEOGRAPHIC_POINT *p, double shift)
Shift a point around by a number of radians.
Definition: lwgeodetic.c:137
void y_to_z(POINT3D *p)
Definition: lwgeodetic.c:607
double lwgeom_length_spheroid(const LWGEOM *geom, const SPHEROID *s)
Calculate the geodetic length of a lwgeom on the unit sphere.
Definition: lwgeodetic.c:2881
static void vector_scale(POINT3D *n, double scale)
Scale a vector out by a factor.
Definition: lwgeodetic.c:438
int gbox_init_point3d(const POINT3D *p, GBOX *gbox)
Initialize a GBOX using the values of the point.
Definition: g_box.c:197
int lwgeom_force_geodetic(LWGEOM *geom)
Force coordinates of LWGEOM into geodetic range (-180, -90, 180, 90)
Definition: lwgeodetic.c:2806
static int lwline_calculate_gbox_geodetic(const LWLINE *line, GBOX *gbox)
Definition: lwgeodetic.c:2534
LWGEOM * lwgeom_clone_deep(const LWGEOM *lwgeom)
Deep clone an LWGEOM, everything is copied.
Definition: lwgeom.c:389
GBOX * bbox
Definition: liblwgeom.h:409
LWPOLY * lwgeom_as_lwpoly(const LWGEOM *lwgeom)
Definition: lwgeom.c:125
double longitude_degrees_normalize(double lon)
Convert a longitude to the range of -180,180.
Definition: lwgeodetic.c:83
#define FP_MIN(A, B)
Point in spherical coordinates on the world.
Definition: lwgeodetic.h:33
POINTARRAY * point
Definition: liblwgeom.h:367
int lwgeom_has_z(const LWGEOM *geom)
Return LW_TRUE if geometry has Z ordinates.
Definition: lwgeom.c:792
char ** result
Definition: liblwgeom.h:218
int32_t srid
Definition: liblwgeom.h:355
LWGEOM * lwpoly_as_lwgeom(const LWPOLY *obj)
Definition: lwgeom.c:239
#define FP_IS_ZERO(A)
static double sphere_angle(const GEOGRAPHIC_POINT *a, const GEOGRAPHIC_POINT *b, const GEOGRAPHIC_POINT *c)
Returns the angle in radians at point B of the triangle formed by A-B-C.
Definition: lwgeodetic.c:670
int spheroid_project(const GEOGRAPHIC_POINT *r, const SPHEROID *spheroid, double distance, double azimuth, GEOGRAPHIC_POINT *g)
Given a location, an azimuth and a distance, computes the location of the projected point...
Definition: lwspheroid.c:232
int gbox_merge(const GBOX *new_box, GBOX *merge_box)
Update the merged GBOX to be large enough to include itself and the new box.
Definition: g_box.c:215
#define PIR_A_TOUCH_LEFT
Definition: lwgeodetic.h:76
#define LW_FAILURE
Definition: liblwgeom.h:54
double latitude_degrees_normalize(double lat)
Convert a latitude to the range of -90,90.
Definition: lwgeodetic.c:110
double z
Definition: liblwgeom.h:296
void lwerror(const char *fmt,...)
Write a notice out to the error handler.
Definition: lwutil.c:67
double x
Definition: liblwgeom.h:284
int clairaut_geographic(const GEOGRAPHIC_POINT *start, const GEOGRAPHIC_POINT *end, GEOGRAPHIC_POINT *g_top, GEOGRAPHIC_POINT *g_bottom)
Computes the pole of the great circle disk which is the intersection of the great circle with the lin...
Definition: lwgeodetic.c:1048
void unit_normal(const POINT3D *P1, const POINT3D *P2, POINT3D *normal)
Calculates the unit normal to two vectors, trying to avoid problems with over-narrow or over-wide cas...
Definition: lwgeodetic.c:490
double lwpoint_get_x(const LWPOINT *point)
Definition: lwpoint.c:48
int ptarray_contains_point_sphere(const POINTARRAY *pa, const POINT2D *pt_outside, const POINT2D *pt_to_test)
This routine returns LW_TRUE if the stabline joining the pt_outside and pt_to_test crosses the ring a...
Definition: lwgeodetic.c:3201
static int lwpoint_force_geodetic(LWPOINT *point)
Definition: lwgeodetic.c:2766
double zmax
Definition: liblwgeom.h:253
double ymin
Definition: liblwgeom.h:250
const char * lwtype_name(uint8_t type)
Return the type name string associated with a type number (e.g.
Definition: lwutil.c:164
static double dot_product(const POINT3D *p1, const POINT3D *p2)
Convert cartesion coordinates on unit sphere to lon/lat coordinates static void cart2ll(const POINT3D...
Definition: lwgeodetic.c:397
LWGEOM * lwline_as_lwgeom(const LWLINE *obj)
Definition: lwgeom.c:249
void gbox_init(GBOX *gbox)
Zero out all the entries in the GBOX.
Definition: g_box.c:34
double xmin
Definition: liblwgeom.h:248
int ptarray_append_point(POINTARRAY *pa, const POINT4D *pt, int allow_duplicates)
Append a point to the end of an existing POINTARRAY If allow_duplicate is LW_TRUE, then a duplicate point will not be added.
Definition: ptarray.c:141
#define LW_FALSE
Definition: liblwgeom.h:52
const POINT2D * getPoint2d_cp(const POINTARRAY *pa, int n)
Returns a POINT2D pointer into the POINTARRAY serialized_ptlist, suitable for reading from...
Definition: lwgeom_api.c:458
uint8_t flags
Definition: liblwgeom.h:325
#define rad2deg(r)
Definition: lwgeodetic.h:61
GEOGRAPHIC_POINT start
Definition: lwgeodetic.h:44
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:51
LWLINE * lwline_construct(int srid, GBOX *bbox, POINTARRAY *points)
Definition: lwline.c:29
double edge_distance_to_point(const GEOGRAPHIC_EDGE *e, const GEOGRAPHIC_POINT *gp, GEOGRAPHIC_POINT *closest)
Definition: lwgeodetic.c:1165
int edge_point_in_cone(const GEOGRAPHIC_EDGE *e, const GEOGRAPHIC_POINT *p)
Returns true if the point p is inside the cone defined by the two ends of the edge e...
Definition: lwgeodetic.c:737
double lwgeom_area_sphere(const LWGEOM *lwgeom, const SPHEROID *spheroid)
Calculate the area of an LWGEOM.
Definition: lwgeodetic.c:1924
LWGEOM ** geoms
Definition: liblwgeom.h:465
double ptarray_length_spheroid(const POINTARRAY *pa, const SPHEROID *s)
Definition: lwgeodetic.c:2828
int edge_contains_coplanar_point(const GEOGRAPHIC_EDGE *e, const GEOGRAPHIC_POINT *p)
True if the longitude of p is within the range of the longitude of the ends of e. ...
Definition: lwgeodetic.c:784
void ll2cart(const POINT2D *g, POINT3D *p)
Convert lon/lat coordinates to cartesion coordinates on unit sphere.
Definition: lwgeodetic.c:374
int clairaut_cartesian(const POINT3D *start, const POINT3D *end, GEOGRAPHIC_POINT *g_top, GEOGRAPHIC_POINT *g_bottom)
Computes the pole of the great circle disk which is the intersection of the great circle with the lin...
Definition: lwgeodetic.c:1023
#define TINTYPE
Definition: liblwgeom.h:74
void vector_rotate(const POINT3D *v1, const POINT3D *v2, double angle, POINT3D *n)
Rotates v1 through an angle (in radians) within the plane defined by v1/v2, returns the rotated vecto...
Definition: lwgeodetic.c:522
void x_to_z(POINT3D *p)
Definition: lwgeodetic.c:600
int lwgeom_calculate_gbox_geodetic(const LWGEOM *geom, GBOX *gbox)
Calculate the geodetic bounding box for an LWGEOM.
Definition: lwgeodetic.c:2612
int gbox_merge_point3d(const POINT3D *p, GBOX *gbox)
Update the GBOX to be large enough to include itself and the new point.
Definition: g_box.c:186
static void vector_difference(const POINT3D *a, const POINT3D *b, POINT3D *n)
Calculate the difference of two vectors.
Definition: lwgeodetic.c:427
uint8_t * getPoint_internal(const POINTARRAY *pa, int n)
Definition: ptarray.c:1645
static int lwpoint_calculate_gbox_geodetic(const LWPOINT *point, GBOX *gbox)
Definition: lwgeodetic.c:2528
static int lwcollection_force_geodetic(LWCOLLECTION *col)
Definition: lwgeodetic.c:2792
POINTARRAY ** rings
Definition: liblwgeom.h:413
int count
Definition: genraster.py:57
int lwtype_is_collection(uint8_t type)
Determine whether a type number is a collection or not.
Definition: lwgeom.c:955
static POINTARRAY * ptarray_segmentize_sphere(const POINTARRAY *pa_in, double max_seg_length)
Create a new point array with no segment longer than the input segment length (expressed in radians!)...
Definition: lwgeodetic.c:1530
void lwpoly_pt_outside(const LWPOLY *poly, POINT2D *pt_outside)
Definition: lwgeodetic.c:1422
int nrings
Definition: liblwgeom.h:411
double ymax
Definition: liblwgeom.h:251
GEOGRAPHIC_POINT end
Definition: lwgeodetic.h:45
int lwgeom_nudge_geodetic(LWGEOM *geom)
When features are snapped or sometimes they are just this way, they are very close to the geodetic bo...
Definition: lwgeodetic.c:2988
char * s
Definition: cu_in_wkt.c:24
double y
Definition: liblwgeom.h:284
static int point_in_cone(const POINT3D *A1, const POINT3D *A2, const POINT3D *P)
Utility function for checking if P is within the cone defined by A1/A2.
Definition: lwgeodetic.c:3043
int ptarray_calculate_gbox_geodetic(const POINTARRAY *pa, GBOX *gbox)
Calculate geodetic (x/y/z) box and add values to gbox.
Definition: lwgeodetic.c:2473
int getPoint2d_p(const POINTARRAY *pa, int n, POINT2D *point)
Definition: lwgeom_api.c:434
#define FLAGS_GET_Z(flags)
Macros for manipulating the 'flags' byte.
Definition: liblwgeom.h:106
double gbox_angular_width(const GBOX *gbox)
Returns the angular width (longitudinal span) of the box in radians.
Definition: lwgeodetic.c:192
double z
Definition: liblwgeom.h:308
tuple x
Definition: pixval.py:53
static int lwpolygon_calculate_gbox_geodetic(const LWPOLY *poly, GBOX *gbox)
Definition: lwgeodetic.c:2540
#define POW2(x)
Definition: lwgeodetic.h:28
LWGEOM * lwgeom_clone(const LWGEOM *lwgeom)
Clone LWGEOM object.
Definition: lwgeom.c:351
Datum distance(PG_FUNCTION_ARGS)
static int lwline_check_geodetic(const LWLINE *line)
Definition: lwgeodetic.c:2676
void lwgeom_set_geodetic(LWGEOM *geom, int value)
Set the FLAGS geodetic bit on geometry an all sub-geometries and pointlists.
Definition: lwgeom.c:814
LWLINE * lwgeom_as_lwline(const LWGEOM *lwgeom)
Definition: lwgeom.c:89
int32_t srid
Definition: liblwgeom.h:366
uint8_t flags
Definition: liblwgeom.h:247
uint8_t gflags(int hasz, int hasm, int geodetic)
Construct a new flags char.
Definition: g_util.c:131
void geog2cart(const GEOGRAPHIC_POINT *g, POINT3D *p)
Convert spherical coordinates to cartesion coordinates on unit sphere.
Definition: lwgeodetic.c:355
void gbox_pt_outside(const GBOX *gbox, POINT2D *pt_outside)
Given a unit geocentric gbox, return a lon/lat (degrees) coordinate point point that is guaranteed to...
Definition: lwgeodetic.c:1443
int gbox_geocentric_slow
For testing geodetic bounding box, we have a magic global variable.
Definition: lwgeodetic.c:22
#define MULTIPOLYGONTYPE
Definition: liblwgeom.h:65
void vector_sum(const POINT3D *a, const POINT3D *b, POINT3D *n)
Calculate the sum of two vectors.
Definition: lwgeodetic.c:416
double latitude_radians_normalize(double lat)
Convert a latitude to the range of -PI/2,PI/2.
Definition: lwgeodetic.c:55
int edge_intersects(const POINT3D *A1, const POINT3D *A2, const POINT3D *B1, const POINT3D *B2)
Returns non-zero if edges A and B interact.
Definition: lwgeodetic.c:3096
double sphere_distance_cartesian(const POINT3D *s, const POINT3D *e)
Given two unit vectors, calculate their distance apart in radians.
Definition: lwgeodetic.c:916
int ptarray_has_m(const POINTARRAY *pa)
Definition: ptarray.c:30
int geographic_point_equals(const GEOGRAPHIC_POINT *g1, const GEOGRAPHIC_POINT *g2)
Definition: lwgeodetic.c:147
LWGEOM * lwgeom_segmentize_sphere(const LWGEOM *lwg_in, double max_seg_length)
Create a new, densified geometry where no segment is longer than max_seg_length.
Definition: lwgeodetic.c:1638
double a
Definition: liblwgeom.h:269
LWCOLLECTION * lwgeom_as_lwcollection(const LWGEOM *lwgeom)
Definition: lwgeom.c:143
static int lwcollection_check_geodetic(const LWCOLLECTION *col)
Definition: lwgeodetic.c:2702
double zmin
Definition: liblwgeom.h:252
void geographic_point_init(double lon, double lat, GEOGRAPHIC_POINT *g)
Initialize a geographic point.
Definition: lwgeodetic.c:157
double ptarray_area_sphere(const POINTARRAY *pa)
Returns the area of the ring (ring must be closed) in square radians (surface of the sphere is 4*PI)...
Definition: lwgeodetic.c:1702
#define POINTTYPE
LWTYPE numbers, used internally by PostGIS.
Definition: liblwgeom.h:60
static int lwtriangle_check_geodetic(const LWTRIANGLE *triangle)
Definition: lwgeodetic.c:2695
double lwgeom_azumith_spheroid(const LWPOINT *r, const LWPOINT *s, const SPHEROID *spheroid)
Calculate a bearing (azimuth) given a source and destination point.
Definition: lwgeodetic.c:2047
LWPOLY * lwpoly_construct_empty(int srid, char hasz, char hasm)
Definition: lwpoly.c:66
#define FLAGS_GET_M(flags)
Definition: liblwgeom.h:107
int p3d_same(const POINT3D *p1, const POINT3D *p2)
Definition: lwalgorithm.c:38
#define deg2rad(d)
Conversion functions.
Definition: lwgeodetic.h:60
#define MAXFLOAT
Largest float value.
int gbox_overlaps(const GBOX *g1, const GBOX *g2)
Return LW_TRUE if the GBOX overlaps, LW_FALSE otherwise.
Definition: g_box.c:241
uint8_t type
Definition: liblwgeom.h:352
static double sphere_signed_area(const GEOGRAPHIC_POINT *a, const GEOGRAPHIC_POINT *b, const GEOGRAPHIC_POINT *c)
Computes the spherical area of a triangle.
Definition: lwgeodetic.c:690
#define FP_EQUALS(A, B)
double lwpoint_get_y(const LWPOINT *point)
Definition: lwpoint.c:58
LWGEOM * lwpoint_as_lwgeom(const LWPOINT *obj)
Definition: lwgeom.c:254
LWPOINT * lwgeom_project_spheroid(const LWPOINT *r, const SPHEROID *spheroid, double distance, double azimuth)
Calculate a projected point given a source point, a distance and a bearing.
Definition: lwgeodetic.c:1992
static int lwpoly_check_geodetic(const LWPOLY *poly)
Definition: lwgeodetic.c:2682
double spheroid_distance(const GEOGRAPHIC_POINT *a, const GEOGRAPHIC_POINT *b, const SPHEROID *spheroid)
Computes the shortest distance along the surface of the spheroid between two points.
Definition: lwspheroid.c:59
int lwpoly_add_ring(LWPOLY *poly, POINTARRAY *pa)
Add a ring, allocating extra space if necessary.
Definition: lwpoly.c:154
static int lwpoly_force_geodetic(LWPOLY *poly)
Definition: lwgeodetic.c:2778
int ptarray_has_z(const POINTARRAY *pa)
Definition: ptarray.c:23
#define PIR_NO_INTERACT
Bitmask elements for edge_intersects() return value.
Definition: lwgeodetic.h:72
int lwpoly_covers_point2d(const LWPOLY *poly, const POINT2D *pt_to_test)
Given a polygon (lon/lat decimal degrees) and point (lon/lat decimal degrees) and a guaranteed outsid...
Definition: lwgeodetic.c:2384
LWPOINT * lwpoint_construct(int srid, GBOX *bbox, POINTARRAY *point)
Definition: lwpoint.c:96
int lwgeom_is_empty(const LWGEOM *geom)
Return true or false depending on whether a geometry is an "empty" geometry (no vertices members) ...
Definition: lwgeom.c:1229
#define PIR_B_TOUCH_LEFT
Definition: lwgeodetic.h:78
static int lwpoint_check_geodetic(const LWPOINT *point)
Definition: lwgeodetic.c:2670
int lw_segment_side(const POINT2D *p1, const POINT2D *p2, const POINT2D *q)
lw_segment_side()
Definition: lwalgorithm.c:62
double y
Definition: liblwgeom.h:308
#define MULTILINETYPE
Definition: liblwgeom.h:64
LWCOLLECTION * lwcollection_construct_empty(uint8_t type, int srid, char hasz, char hasm)
Definition: lwcollection.c:81
int edge_point_on_plane(const GEOGRAPHIC_EDGE *e, const GEOGRAPHIC_POINT *p)
Returns true if the point p is on the great circle plane.
Definition: lwgeodetic.c:724
double vector_angle(const POINT3D *v1, const POINT3D *v2)
Angle between two unit vectors.
Definition: lwgeodetic.c:454
static int gbox_check_poles(GBOX *gbox)
Check to see if this geocentric gbox is wrapped around a pole.
Definition: lwgeodetic.c:293
static int edge_point_side(const GEOGRAPHIC_EDGE *e, const GEOGRAPHIC_POINT *p)
Returns -1 if the point is to the left of the plane formed by the edge, 1 if the point is to the righ...
Definition: lwgeodetic.c:643
static void cross_product(const POINT3D *a, const POINT3D *b, POINT3D *n)
Calculate the cross product of two vectors.
Definition: lwgeodetic.c:405
#define LWDEBUGF(level, msg,...)
Definition: lwgeom_log.h:55
static void normalize2d(POINT2D *p)
Normalize to a unit vector.
Definition: lwgeodetic.c:473
static int ptarray_force_geodetic(POINTARRAY *pa)
Definition: lwgeodetic.c:2744
#define PIR_COLINEAR
Definition: lwgeodetic.h:74
int edge_calculate_gbox(const POINT3D *A1, const POINT3D *A2, GBOX *gbox)
The magic function, given an edge in spherical coordinates, calculate a 3D bounding box that fully co...
Definition: lwgeodetic.c:1355
int lwgeom_has_m(const LWGEOM *geom)
Return LW_TRUE if geometry has M ordinates.
Definition: lwgeom.c:799
LWCOLLECTION * lwcollection_add_lwgeom(LWCOLLECTION *col, const LWGEOM *geom)
Appends geom to the collection managed by col.
Definition: lwcollection.c:174
int lwgeom_check_geodetic(const LWGEOM *geom)
Check that coordinates of LWGEOM are all within the geodetic range (-180, -90, 180, 90)
Definition: lwgeodetic.c:2715
int edge_intersection(const GEOGRAPHIC_EDGE *e1, const GEOGRAPHIC_EDGE *e2, GEOGRAPHIC_POINT *g)
Returns true if an intersection can be calculated, and places it in *g.
Definition: lwgeodetic.c:1074
static int lwcollection_calculate_gbox_geodetic(const LWCOLLECTION *coll, GBOX *gbox)
Definition: lwgeodetic.c:2577
tuple y
Definition: pixval.py:54
int getPoint4d_p(const POINTARRAY *pa, int n, POINT4D *point)
Definition: lwgeom_api.c:217
#define COLLECTIONTYPE
Definition: liblwgeom.h:66
int p4d_same(const POINT4D *p1, const POINT4D *p2)
Definition: lwalgorithm.c:29
#define FP_MAX(A, B)
#define PIR_INTERSECTS
Definition: lwgeodetic.h:73
POINTARRAY * points
Definition: liblwgeom.h:378
LWGEOM * lwcollection_as_lwgeom(const LWCOLLECTION *obj)
Definition: lwgeom.c:219
double sphere_direction(const GEOGRAPHIC_POINT *s, const GEOGRAPHIC_POINT *e, double d)
Given two points on a unit sphere, calculate the direction from s to e.
Definition: lwgeodetic.c:924
static int ptarray_check_geodetic(const POINTARRAY *pa)
Definition: lwgeodetic.c:2652