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