PostGIS  2.5.0beta2dev-r@@SVN_REVISION@@

◆ edge_calculate_gbox()

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 contains it, taking into account the curvature of the sphere on which it is inscribed.

Any arc on the sphere defines a plane that bisects the sphere. In this plane, the arc is a portion of a unit circle. Projecting the end points of the axes (1,0,0), (-1,0,0) etc, into the plane and normalizing yields potential extrema points. Those points on the side of the plane-dividing line formed by the end points that is opposite the origin of the plane are extrema and should be added to the bounding box.

Definition at line 1377 of file lwgeodetic.c.

References dot_product(), FP_EQUALS, gbox_init_point3d(), gbox_merge_point3d(), LW_FAILURE, lw_segment_side(), LW_SUCCESS, lwerror(), normalize2d(), p3d_same(), unit_normal(), POINT2D::x, POINT3D::x, POINT2D::y, POINT3D::y, and POINT3D::z.

Referenced by ptarray_calculate_gbox_geodetic().

1378 {
1379  POINT2D R1, R2, RX, O;
1380  POINT3D AN, A3;
1381  POINT3D X[6];
1382  int i, o_side;
1383 
1384  /* Initialize the box with the edge end points */
1385  gbox_init_point3d(A1, gbox);
1386  gbox_merge_point3d(A2, gbox);
1387 
1388  /* Zero length edge, just return! */
1389  if ( p3d_same(A1, A2) )
1390  return LW_SUCCESS;
1391 
1392  /* Error out on antipodal edge */
1393  if ( FP_EQUALS(A1->x, -1*A2->x) && FP_EQUALS(A1->y, -1*A2->y) && FP_EQUALS(A1->z, -1*A2->z) )
1394  {
1395  lwerror("Antipodal (180 degrees long) edge detected!");
1396  return LW_FAILURE;
1397  }
1398 
1399  /* Create A3, a vector in the plane of A1/A2, orthogonal to A1 */
1400  unit_normal(A1, A2, &AN);
1401  unit_normal(&AN, A1, &A3);
1402 
1403  /* Project A1 and A2 into the 2-space formed by the plane A1/A3 */
1404  R1.x = 1.0;
1405  R1.y = 0.0;
1406  R2.x = dot_product(A2, A1);
1407  R2.y = dot_product(A2, &A3);
1408 
1409  /* Initialize our 3-space axis points (x+, x-, y+, y-, z+, z-) */
1410  memset(X, 0, sizeof(POINT3D) * 6);
1411  X[0].x = X[2].y = X[4].z = 1.0;
1412  X[1].x = X[3].y = X[5].z = -1.0;
1413 
1414  /* Initialize a 2-space origin point. */
1415  O.x = O.y = 0.0;
1416  /* What side of the line joining R1/R2 is O? */
1417  o_side = lw_segment_side(&R1, &R2, &O);
1418 
1419  /* Add any extrema! */
1420  for ( i = 0; i < 6; i++ )
1421  {
1422  /* Convert 3-space axis points to 2-space unit vectors */
1423  RX.x = dot_product(&(X[i]), A1);
1424  RX.y = dot_product(&(X[i]), &A3);
1425  normalize2d(&RX);
1426 
1427  /* Any axis end on the side of R1/R2 opposite the origin */
1428  /* is an extreme point in the arc, so we add the 3-space */
1429  /* version of the point on R1/R2 to the gbox */
1430  if ( lw_segment_side(&R1, &R2, &RX) != o_side )
1431  {
1432  POINT3D Xn;
1433  Xn.x = RX.x * A1->x + RX.y * A3.x;
1434  Xn.y = RX.x * A1->y + RX.y * A3.y;
1435  Xn.z = RX.x * A1->z + RX.y * A3.z;
1436 
1437  gbox_merge_point3d(&Xn, gbox);
1438  }
1439  }
1440 
1441  return LW_SUCCESS;
1442 }
double y
Definition: liblwgeom.h:342
#define LW_SUCCESS
Definition: liblwgeom.h:79
double x
Definition: liblwgeom.h:342
int gbox_init_point3d(const POINT3D *p, GBOX *gbox)
Initialize a GBOX using the values of the point.
Definition: g_box.c:246
#define LW_FAILURE
Definition: liblwgeom.h:78
double z
Definition: liblwgeom.h:342
double x
Definition: liblwgeom.h:330
void unit_normal(const POINT3D *P1, const POINT3D *P2, POINT3D *normal)
Calculates the unit normal to two vectors, trying to avoid problems with over-narrow or over-wide cas...
Definition: lwgeodetic.c:510
static double dot_product(const POINT3D *p1, const POINT3D *p2)
Convert cartesian coordinates on unit sphere to lon/lat coordinates static void cart2ll(const POINT3D...
Definition: lwgeodetic.c:415
int gbox_merge_point3d(const POINT3D *p, GBOX *gbox)
Update the GBOX to be large enough to include itself and the new point.
Definition: g_box.c:235
double y
Definition: liblwgeom.h:330
int p3d_same(const POINT3D *p1, const POINT3D *p2)
Definition: lwalgorithm.c:40
#define FP_EQUALS(A, B)
int lw_segment_side(const POINT2D *p1, const POINT2D *p2, const POINT2D *q)
lw_segment_side()
Definition: lwalgorithm.c:64
static void normalize2d(POINT2D *p)
Normalize to a unit vector.
Definition: lwgeodetic.c:493
void lwerror(const char *fmt,...)
Write a notice out to the error handler.
Definition: lwutil.c:190
Here is the call graph for this function:
Here is the caller graph for this function: