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