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