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